Fundamentos de álgebra lineal en

Acerca de Fundamentos de álgebra lineal

¿Qué es el «álgebra lineal»?

Existen muchas definiciones técnicas, en Wikipedia, en libros de texto y en sitios excelentes como 3Blue1Brown, pero para nuestros propósitos 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 advierto: esto es, inevitablemente, un concepto bastante matemático.

¿Qué es un vector?

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

• En Julia, Vector{T} es solo un alias de Array{T, 1}, y el eltype T puede ser cualquier cosa.
• En (gran 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á fija en el espacio, con su cola en el origen. Esto hace que el cuerpo 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 de caracteres. En este concepto trabajaremos con tipos numéricos: Int, Float o Complex.

¿Qué es una matriz?

De nuevo, encontrarás una variedad de respuestas.

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

Algunos tipos de matrices cuadradas 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

Triangular superior: Valores distintos de cero sobre la diagonal y por encima de ella, ceros por debajo. Triangular inferior ya te lo imaginas.

Transponer una matriz

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

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

Para los propósitos del álgebra lineal, la función transpose() es más útil, ya que crea rápidamente un envoltorio perezoso alrededor de la matriz original.

Aún más útil es la función adjoint(), que además cambia 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 tan común que basta 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 lo repetiremos en los ejemplos (demasiado ruido visual).

Vectores, producto elemento a elemento

Esto 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 punto

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

Un producto punto 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: lo que 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 tabulador) como azúcar sintáctico para la función dot(). No hace falta especificar la adjunta, porque este detalle se maneja 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 debe ser el conjugado (con el signo de la parte imaginaria cambiado). 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 cruz

¡Un poco de nostalgia para quien haya tomado un curso de Electricidad y Magnetismo en el pasado! También, para ingenieros familiarizados con el cálculo de vectores de torque o de momento angular.

Mientras que el producto punto convierte dos vectores en un scalar, el producto cruz 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 nos ha dejado a muchos 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 que esta operación está restringida a vectores de longitud 3 (equivalente al espacio euclidiano 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 corriendo o cambiamos de tema rápidamente en ese punto!

Normas

¿Qué tan «grande» es un vector?

La norm es un intento de capturar esto, reduciendo el vector a un escalar apropiado.

Se define toda una familia de normas, pero de lejos 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 desde el origen (en un espacio de N dimensiones). Si visualizamos el vector como una flecha, con su cola en el origen, la norma 2 es la longitud de la flecha.

Cualquier norma p se puede calcular proporcionando p como segundo argumento. La norma 1 a veces es ú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 gran parte de los últimos 80 años convirtiendo cada cálculo en una serie de multiplicaciones de matrices.

Esta es una operación en la que las computadoras son muy buenas:

• Es altamente repetitiva.
• Se puede paralelizar de manera eficiente.
• Se ha desarrollado mucho hardware especializado para hacerla más rápida, desde el Cray-1 de 80 millones de dólares en la década de 1970 hasta la GPU que probablemente está integrada en tu computadora 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 como resultado 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 punto 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 las columnas de la matriz de la derecha.

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

De manera equivalente, podríamos decir que la primera fila de A se multiplica punto 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 como esta, y quienes estén interesados pueden ver toda una clase del MIT 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 simple entre A*B y B*A.

Probablemente sea difícil de visualizar para quien recién se acerca a la multiplicación de matrices, solo leyendo las palabras. YouTube tiene muchos videos que la 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 punto de dos vectores requiere que 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 punto. 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 punto. Las dimensiones son (1, 3) * (3, 1), y Julia simplifica la salida (1, 1) a un escalar (a diferencia de, por ejemplo, R).

Alternativamente, podríamos usar u * v', con dimensiones (3, 1) * (1, 3) -> (3, 3), para multiplicar todos los pares posibles de elementos y expandir los vectores hasta formar 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 producto interno.

El (menos común) u * v' es, correspondientemente, el producto externo. Esto se relaciona con el producto tensorial, aunque eso está mucho más allá de nuestro alcance.

Rotaciones

Un tipo de multiplicación de matrices particularmente común involucra matrices de rotación.

En 2D, hay una matriz relativamente simple 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 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 elementos (o, de manera equivalente, una columna de la 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 la inclinación, la perspectiva, etc.

Rendimiento

En los departamentos de ciencias de la computación 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 para su propio beneficio.

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

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

El sistema usado fue una computadora pequeña de 430 dólares (USD): 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 tardó unos pocos milisegundos, y las matrices de 100 millones de elementos tardaron unos pocos segundos. Añadir varios órdenes de magnitud más necesitará hardware mejor...