La computación numérica es el territorio donde Julia nació y donde todavía brilla más. No hay capas de Python encima de C: el código Julia que escribe un científico ya es código nativo. Cuando multiplicas dos matrices con A * B, Julia llama directamente a DGEMM de BLAS sin intermediarios. El resultado es que los bucles numéricos en Julia compiten con Fortran.
Arrays: la estructura central
Los arrays en Julia son multidimensionales por diseño. Un Array{Float64, 2} es una matriz, y Array{Float64, 1} es un vector. Las operaciones básicas de creación son directas:
# Crear arrays v = zeros(5) # vector de 5 ceros Float64 m = ones(3, 4) # matriz 3x4 de unos r = rand(100) # 100 valores aleatorios uniformes en [0,1] n = randn(50) # 50 valores de distribución normal # Rango con paso x = range(0, 2*pi, length=100) # 100 puntos de 0 a 2? y = range(1, 10, step=0.5) # de 1 a 10 con paso 0.5 # Array literal A = [1 2 3; 4 5 6; 7 8 9] # matriz 3x3 # Dimensiones e información size(A) # (3, 3) length(A) # 9 eltype(A) # Int64 ndims(A) # 2
La indexación en Julia empieza en 1. Puede parecer extraño para quien viene de Python, pero es consistente con Matlab y Fortran, los lenguajes de referencia en computación científica. La notación de rangos es A[2:4, :] para filas 2 a 4, todas las columnas.
Broadcasting
El operador punto (.) aplica cualquier función o operador elemento a elemento. Es uniforme para toda la librería, no solo para funciones de NumPy:
v = [1.0, 2.0, 3.0, 4.0, 5.0] # Funciones matemáticas elemento a elemento sin.(v) # sin de cada elemento log.(v) # log de cada elemento sqrt.(v) # raíz cuadrada de cada elemento # Operaciones aritméticas 2 .* v .+ 1 # 2x + 1 para cada x v .^ 2 # cuadrado de cada elemento v ./ sum(v) # normalizar # Broadcasting entre arrays de diferente forma A = rand(3, 3) b = rand(3) A .+ b' # suma broadcasting: b' se extiende por filas
Lo mejor del broadcasting con punto es que aplica a cualquier función, incluidas las que escribas tú. Si tienes una función f(x), f.(v) la aplica a cada elemento de v sin que tengas que hacer nada especial.
LinearAlgebra: BLAS y LAPACK nativos
La librería estándar LinearAlgebra proporciona acceso directo a las rutinas de álgebra lineal de alto rendimiento. OpenBLAS viene por defecto, y se puede cambiar a MKL de Intel para hardware específico.
using LinearAlgebra A = rand(4, 4) b = rand(4) # Multiplicación de matrices (llama a DGEMM) C = A * A' # A por su transpuesta # Factorizaciones F_lu = lu(A) # factorización LU F_qr = qr(A) # factorización QR F_svd = svd(A) # SVD: valores singulares F_eig = eigen(A) # valores y vectores propios # Resolver sistema lineal Ax = b x = A b # equivalente a solve(A, b), usa LU por defecto # Normas norm(b) # norma L2 norm(A, Inf) # norma infinito opnorm(A) # norma de operador # Determinante, traza, rango det(A) tr(A) rank(A)
El operador para resolver sistemas lineales es uno de los aciertos de diseño de Julia. Internamente elige el algoritmo más apropiado según la estructura de la matriz: LU general, Cholesky para matrices definidas positivas, o métodos especiales para matrices dispersas.
Matrices dispersas
Para matrices con muchos ceros, SparseArrays de la librería estándar es eficiente en memoria y en operaciones:
using SparseArrays # Crear matriz dispersa S = sprand(1000, 1000, 0.01) # 1000x1000 con 1% de densidad # Desde coordenadas i = [1, 3, 5]; j = [2, 4, 6]; v = [1.0, 2.0, 3.0] S2 = sparse(i, j, v, 6, 6) # Las operaciones estándar funcionan igual b = rand(1000) x = S b # resolver sistema disperso # Convertir dense = Matrix(S) # dispersa a densa (cuidado con matrices grandes)
Vistas sin copia
Una de las causas más frecuentes de código Julia más lento de lo esperado es crear copias innecesarias al hacer slicing. @view evita eso:
A = rand(1000, 1000)
# Sin @view: crea una copia de la submatriz
sub = A[1:500, 1:500] # 500*500*8 bytes copiados
# Con @view: referencia a la memoria original, sin copia
sub_view = @view A[1:500, 1:500]
# Para múltiples vistas en una expresión
@views begin
fila = A[3, :]
col = A[:, 7]
bloque = A[1:10, 1:10]
end
Números de punto flotante y precisión
Julia tiene soporte nativo para diferentes precisiones:
# Tipos de punto flotante x32 = Float32(3.14) # 32 bits (más rápido en GPU) x64 = Float64(3.14) # 64 bits (defecto, máxima compatibilidad) xbig = BigFloat(3.14) # precisión arbitraria (más lento) # Aritmética de intervalos (con IntervalArithmetic.jl) # Propagación de errores de redondeo controlada # Números complejos nativos z = 3.0 + 4.0im abs(z) # módulo: 5.0 angle(z) # argumento en radianes conj(z) # conjugado
Esta variedad de tipos numéricos, combinada con el multiple dispatch que ya vimos en el artículo sobre el sistema de tipos, permite escribir algoritmos numéricos genéricos que funcionan igual con Float32, Float64 o BigFloat. En el contexto de machine learning con Flux.jl, esto significa que el mismo código de entrenamiento puede correr en CPU con Float64 o en GPU con Float32 sin cambios.
Para quien viene de Python y NumPy, la diferencia más notable no está en la sintaxis sino en que estos bucles se ejecutan a velocidad nativa. No hay el problema del «bucle lento de Python» que obliga a vectorizar todo: en Julia un bucle explícito sobre millones de elementos es tan rápido como la versión vectorizada.
Imagen: Pexels / Pranjall Kumar