Fundamentos de álgebra lineal en

Acerca de Fundamentos de álgebra lineal

¿Qué es el «álgebra lineal»?

Existen muchísimas definiciones técnicas, en Wikipedia, en los libros de texto y en sitios web excelentes como 3Blue1Brown, pero para lo que nos ocupa podemos ser más informales.

A los mantenedores de Julia (y de Python) les encanta este tipo de cosas. A la dirección de Exercism, no tanto.

Este concepto se limitará a los aspectos más sencillos de un tema enorme, pero te aviso: esto es, inevitablemente, un concepto bastante matemático.

¿Qué es un vector?

Depende de a quién preguntes y de cómo quieras visualizar algo bastante abstracto.

• En Julia, Vector{T} no es más que un alias de Array{T, 1}, y el eltype T puede ser cualquier cosa.
• En (buena parte de) la física, un vector es una flecha con una length y una direction en un espacio de N dimensiones, pero sin una posición fija.
• En álgebra lineal, la flecha de los físicos está fijada en el espacio, con su cola en el origen. Esto hace que el asta de la flecha sea redundante, así que podemos representar el vector como un point en un espacio de N dimensiones, donde los N elementos representan la distancia al origen a lo largo de cada uno de los N ejes.

Por ahora, ignora los vectores de strings o caracteres. En este concepto vamos a trabajar con tipos numéricos: Int, Float o Complex.

¿Qué es una matriz?

De nuevo, encontrarás respuestas muy variadas.

• Un array rectangular bidimensional de números (sin bordes irregulares). Una matriz cuadrada es un caso especial muy común.
• Una linear combination de vectores columna, colocados uno junto a otro.
• Una transformation que se puede aplicar a un vector, igual que una function se aplica a otros tipos de entrada.

Algunos tipos de matriz cuadrada son lo bastante comunes como para tener nombres especiales.

Matriz diagonal: todas las entradas distintas de cero están en la main diagonal (de arriba a la izquierda a abajo a la derecha).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Matriz identidad: una matriz diagonal con solo unos en la diagonal (por razones que quedarán más claras después). A menudo se abrevia como I.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Matriz triangular superior: valores distintos de cero en la diagonal y por encima de ella, y ceros por debajo. La matriz triangular inferior ya te la imaginas.

Trasponer una matriz

Si de verdad quieres intercambiar las filas por las columnas, la función permutedims() lo hace y te da una matriz nueva.

Esto implica copiar, lo cual es lento y consume mucha memoria cuando se trabaja con matrices grandes.

Para los fines del álgebra lineal, la función transpose() es más útil, ya que crea rápidamente una envoltura perezosa en torno a la matriz original.

Aún más útil es la función adjoint(), que además invierte el signo de la parte imaginaria de los números complejos (si te preguntas por qué es tan importante, una respuesta rápida es la mecánica cuántica). Esta operación es lo bastante común como para que baste con añadir un apóstrofo ' al nombre de la variable para crear la adjunta.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Multiplicación

El resto de este documento incluye funcionalidad del módulo LinearAlgebra. Una línea using LinearAlgebra lo traerá al espacio de nombres, pero no vamos a repetirlo en los ejemplos (demasiado ruido visual).

Vectores, producto elemento a elemento

Esto ya se trató en el concepto Operaciones con vectores. El operador es .*, que opera sobre los vectores de entrada por pares y da una salida del mismo tamaño y tipo que las entradas.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vectores, producto escalar

Esta operación tan común se escribe en los libros de texto como u ⋅ v, y se llama producto «punto».

Un producto escalar equivale a la suma del producto elemento a elemento.

Dos vectores se pueden multiplicar con el operador * habitual, pero solo si el vector de la izquierda es una adjoint: se escribe cómodamente como u' * v. Esto quedará más claro en la sección posterior sobre multiplicación de matrices.

El punto elevado (centrado) está disponible en Julia (se escribe \cdot y tabulación) como azúcar sintáctico de la función dot(). No hace falta especificar la adjunta, ya que este detalle se gestiona automáticamente.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Para vectores con valores complejos, el vector de la izquierda tiene que ser el conjugado (con el signo de la parte imaginaria invertido). La función dot lo hace automáticamente, la sintaxis u' lo hace de forma explícita, pero sum(u .* v) fallará si u y v son complejos.

Vectores, producto vectorial

Un poco de nostalgia para quien haya cursado alguna vez una asignatura de electricidad y magnetismo. También para los ingenieros familiarizados con el cálculo del par o del momento angular.

Mientras que el producto escalar convierte dos vectores en un scalar, el producto vectorial convierte dos vectores de longitud 3 en un tercer vector de longitud 3.

En la representación geométrica del espacio vectorial, el nuevo vector es perpendicular al plano que contiene los dos vectores de entrada. Si las dos entradas son paralelas (salvo por el signo), no definen un plano, así que la salida será [0, 0, 0].

El orden importa: u × v == -(v × u). Esta es la famosa regla de la mano derecha, que a muchos nos ha dejado mirándonos el pulgar y dos dedos mientras los girábamos en el espacio (a menudo con cara de desconcierto).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

Fíjate en que esta operación está restringida a vectores de longitud 3 (equivalentes al espacio euclídeo con ejes x, y, z ortogonales). Los matemáticos se complacen en definir la operación para otras dimensiones, con una ensalada de palabras incomprensible (tensores de orden superior, producto cuña, producto exterior, multivectores...) como resultado. ¡La mayoría salimos huyendo o cambiamos de tema rápidamente en ese momento!

Normas

¿Cómo de «grande» es un vector?

La norm es un intento de captar esto, reduciendo el vector a un escalar adecuado.

Se define toda una familia de normas, pero con diferencia la más común es la norma 2, que es √(v ⋅ v).

Esta operación de raíz cuadrática media es la distancia pitagórica al origen (en un espacio de N dimensiones). Si visualizamos el vector como una flecha con la cola en el origen, la norma 2 es la longitud de la flecha.

Cualquier norma p se puede calcular pasando p como segundo argumento. La norma 1 a veces resulta útil: es simplemente la suma de los valores absolutos de las entradas (así que es muy rápida y fácil de calcular).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Multiplicación de matrices

A veces parece que todos los matemáticos aplicados del mundo han pasado buena parte de los últimos 80 años convirtiendo cada cálculo en una serie de multiplicaciones de matrices.

Es una operación en la que los ordenadores son muy buenos:

• Es muy repetitiva.
• Se puede paralelizar de forma eficiente.
• Se ha desarrollado mucho hardware especializado para hacerla más rápida, desde el Cray-1 de 80 millones de dólares en los años setenta hasta la GPU que probablemente está integrada en tu portátil.

Los detalles son bastante sencillos, aunque no muy intuitivos a primera vista.

Considera una matriz A que multiplica a un vector v y da una salida w (por convención del álgebra lineal, usaremos letras mayúsculas para las matrices y minúsculas para los vectores).

La primera fila de A se multiplica escalarmente con v para obtener el primer elemento de w, la segunda fila da el segundo elemento, y así sucesivamente hacia abajo.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

Una multiplicación matriz * matriz extiende esto a lo largo de las columnas de la matriz de la derecha.

Para C = A * B, podemos considerar que A multiplica cada columna de B y da la columna correspondiente de C: una serie de multiplicaciones matriz-vector.

De forma equivalente, podríamos decir que la primera fila de A se multiplica escalarmente con cada columna de B para dar la primera fila de C, la segunda fila da la segunda fila, y así sucesivamente hacia abajo.

Hay varias representaciones mentales de este tipo, y quienes tengan interés pueden ver una clase del MIT entera que las analiza.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Como se muestra en el ejemplo anterior, la multiplicación de matrices no es conmutativa: no hay una relación sencilla entre A*B y B*A.

Probablemente sea difícil visualizar la multiplicación de matrices para quien se acerca a ella por primera vez solo leyendo las palabras. YouTube tiene muchos vídeos que lo demuestran gráficamente, así que busca «multiplicación de matrices» y elige uno con tu estilo, nivel de detalle e idioma preferidos.

Dimensiones

El producto escalar de dos vectores depende de que ambos tengan la misma longitud.

Por extensión, en la multiplicación de matrices el número de columnas de la matriz de la izquierda debe coincidir con el número de filas de la matriz de la derecha.

Si expresamos los tamaños como tuplas (nrows, ncols), tal como las devuelve size(A), tenemos (a, b) * (b, c) -> (a, c). Las dimensiones «internas», aquí b y b, son compatibles para calcular productos escalares. Las dimensiones «externas», aquí a y c, determinan las dimensiones de la salida.

Un ejemplo con matrices rectangulares:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Multiplicar pares de vectores al estilo matricial tiene dos posibilidades.

Por convención, usaríamos u' * v como equivalente al producto escalar. Las dimensiones son (1, 3) * (3, 1), y Julia simplifica la salida (1, 1) a un escalar (a diferencia, por ejemplo, de R).

Como alternativa, podríamos usar u * v', con dimensiones (3, 1) * (1, 3) -> (3, 3), para multiplicar todos los pares de elementos posibles y expandir los vectores a una matriz.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

A u' * v a veces se le llama el producto interno.

El (menos común) u * v' es, correspondientemente, el producto externo. Esto está relacionado con el producto tensorial, aunque eso queda muy lejos de nuestro alcance.

Rotaciones

Un tipo de multiplicación de matrices especialmente común implica matrices de rotación.

En 2D, hay una matriz relativamente sencilla para rotar un vector en sentido antihorario θ radianes.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

Las rotaciones en 3D necesitan una matriz más compleja, con dos ángulos. Es difícil de mostrar en un documento Markdown sin incrustar mucho LaTeX, así que consulta la Wikipedia para ver la fórmula.

Para los gráficos en 3D, como OpenGL y sus sucesores, la convención es trabajar con coordenadas homogéneas.

Cada vértice se representa con un vector de 4 componentes (o, de forma equivalente, una columna de matriz): [x, y, z, 1.0] para un punto en (x, y, z).

Esto permite matrices de transformación más elaboradas. La rotación sigue estando en A[1:3, 1:3], las traslaciones son [Δx, Δy, Δz] en A[1:3, 4], el escalado está en la diagonal, y existen otras posibilidades para el sesgo, la perspectiva, etc.

Rendimiento

En los departamentos de informática de todo el mundo hay muchas tesis doctorales con uno o más capítulos sobre pequeñas mejoras incrementales en los algoritmos de multiplicación de matrices. Esto es realmente importante, y varias organizaciones están dispuestas a financiar la investigación en su propio beneficio.

El rendimiento es un tema enorme que no podemos tratar en profundidad, pero a modo de ilustración podemos probar a multiplicar matrices aleatorias de varios tamaños.

El código de abajo es una estimación rápida y poco rigurosa (usa BenchmarkTools.jl para un enfoque mejor).

El sistema utilizado era un PC pequeño de 430 dólares estadounidenses: procesador Ryzen 9, 32 GB de RAM, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

A grandes rasgos, un par de matrices de un millón de elementos tardaron unos milisegundos, y las de 100 millones de elementos, unos segundos. ¡Añadir varios órdenes de magnitud más requerirá un hardware mejor...