Resolución de ecuaciones lineales en

Acerca de Resolución de ecuaciones lineales

¿Te acuerdas de cuando aprendiste álgebra por primera vez, en el instituto?

Normalmente nos daban un sistema de ecuaciones y nos pedían que despejáramos x e y.

4x - 3y = -2
2x + 7y = 16

Dos ecuaciones, dos incógnitas: un poco de sustitución y enseguida ves que x = 1, y = 2.

Fíjate de nuevo en el lado izquierdo del ejemplo anterior. Se parece a una matriz [4 -3; 2 7] multiplicando un vector desconocido [x, y] para dar el vector [-2, 16].

En notación de álgebra lineal: A x = b, siguiendo la convención habitual de letras mayúsculas para las matrices y minúsculas para los vectores.

A x = 0 y el espacio nulo

Un caso especial importante es cuando todas nuestras ecuaciones tienen cero en el lado derecho.

4x - 3y = 0
2x + 7y = 0

En el ejemplo anterior, la única solución es cuando x = y = 0, lo cual no suele ser interesante.

Las soluciones no triviales de A x = 0 solo existen cuando A tiene un determinante nulo.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Considera en cambio estas ecuaciones:

4x - 3y = 0
8x + 6y = 0

La segunda ecuación es simplemente el doble de la primera, así que no aporta información nueva. Cualquier par de valores en los que x = 0.75y será una solución.

Estos valores (x, y) se encuentran a lo largo de una recta en el espacio bidimensional. Hay un número infinito de soluciones, y forman el espacio nulo de la matriz.

En este caso, las filas de la matriz no son linealmente independientes, y el rango de la matriz es menor que el número de filas o de columnas.

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

La función nullspace() devuelve un único punto, elegido de modo que esté normalizado a un vector unitario, pero cualquier múltiplo escalar de este vector también pertenece al espacio nulo.

Resolver A x = b

Ahora supongamos que tienes 1000 ecuaciones con 1000 incógnitas. Si eso te parece una tontería, recuerda que los transductores digitales inalámbricos ahora son baratos y versátiles (miden deformaciones, velocidad del viento, aceleración en 3 ejes, lo que sea...). Tener 1000 de ellos monitorizando una estructura moderna como un puente colgante es perfectamente razonable. Hace falta un sistema de control que interprete el flujo de datos y esté listo para lanzar una alerta si la cosa se vuelve alarmante.

De nuevo, A x = b, donde tenemos A (del diseño de ingeniería) y b (los valores medidos por los transductores), pero necesitamos encontrar x.

Una primera idea, para quien no esté familiarizado con esto, es «dividir entre A» de algún modo para pasarlo al lado derecho.

De hecho, A tiene una inversa (la mayoría de las matrices cuadradas la tienen, aunque no todas: el determinante debe ser distinto de cero), y el cálculo funciona (¡aunque despacio!).

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

Por desgracia, calcular la inversa es lento para matrices pequeñas y glacial para las grandes. No es lo ideal si tu puente se derrumba durante el cálculo.

Por suerte, otros algoritmos son muchísimo más rápidos (en este caso, la eliminación gaussiana).

Julia (copiando a Matlab) simplemente usa una barra invertida para resolver (técnicamente, «división por la izquierda»).

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

Este es un ejemplo trivial, pero hay algunos detalles a los que conviene prestar atención.

• Los coeficientes deben quedar alineados en columnas, así que inserta ceros según sea necesario para los términos que falten.
• Las ecuaciones (y por tanto las filas de la matriz) deben ser linealmente independientes. Si la fila 3 es simplemente la suma de las filas 1 y 2, no aporta información nueva y el problema está subdeterminado.

El comando rank() te dará el número de filas o columnas linealmente independientes, así que procura que sea igual al número de variables desconocidas.

julia> rank(A)
2

Matrices rectangulares

En algún momento, siendo un preadolescente, probablemente te enseñaron que necesitas un sistema de N ecuaciones para resolver N incógnitas.

Como ocurre con casi todo en matemáticas, la realidad es un poco más matizada (y no solo por el detalle de la independencia lineal que mencionamos en la sección anterior).

«Pocas» filas

Con N-1 ecuaciones (filas de la matriz A), el problema está subdeterminado. ¿Nos rendimos desesperados?

¡Depende!

Las N-1 ecuaciones siguen conteniendo mucha información. En términos geométricos, necesitamos N para determinar la solución como un punto en el espacio N-dimensional, pero N-1 nos dirá que la solución debe estar en algún lugar de una recta.

Esto nos deja con un número infinito de soluciones, pero también con una cantidad infinita de espacio de parámetros sin solución.

¿Atraviesa tu recta de soluciones una región del espacio de parámetros por la que un ingeniero de buena reputación debería preocuparse? ¿O que podría justificar cerrar la estructura al público (los puentes colgantes pueden ponerse muy movidos con vientos laterales fuertes)? Quizá merezca la pena hacer algunos cálculos más para determinarlo.

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

¡El solucionador con barra invertida nos da una solución! De manera no oficial, al parecer esta es la solución con la norma más pequeña (la más cercana al origen, geométricamente).

Para obtener las demás soluciones, necesitamos el nullspace de A2.

Sumar cualquier múltiplo escalar del espacio nulo a la solución original dará otra solución válida.

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

De forma similar, N-2 ecuaciones restringirán las soluciones (infinitas) a un plano, y se aplican los mismos principios.

«Demasiadas» filas

La situación opuesta se da cuando hay más ecuaciones que incógnitas y, aun así, siguen siendo linealmente independientes.

En la ingeniería del mundo real, esto es totalmente normal, ¡y se considera algo bueno!

Los transductores tienen precisión limitada, el vector b tiene precisión limitada y hay ruido en el cálculo. Ahora la solución es un ajuste por mínimos cuadrados a los datos ruidosos.

La técnica más sencilla usa la pseudoinverse de la matriz, que Julia implementa como la función pinv().

Se puede usar de forma muy parecida a la inversa de una matriz cuadrada (no singular), y da una estimación por mínimos cuadrados de las variables en lugar de una solución exacta.

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

Valores y vectores propios

La sección anterior describía las soluciones de A x = b, equivalentes al álgebra que ya conoces.

Esta sección trata de las soluciones de A x = λ x, donde λ es un escalar.

Este es un concepto característico del álgebra lineal, pero podemos interpretarlo geométricamente:

Para valores adecuados de x y λ, la matriz cuadrada A escala la longitud del vector x distinto de cero por un factor λ, sin cambiar su dirección (salvo que un λ negativo la invierte).

Esto suena muy especializado, pero resulta ser ridículamente útil.

Terminología: Los valores válidos de λ son los valores propios de A, y los valores correspondientes de x son los vectores propios de A. Por desgracia, no nos queda otra que convivir con palabras que cambian de idioma a mitad de camino.

Normalmente se enseña a los estudiantes a calcular los valores y vectores propios de matrices de 2×2 a mano, pero usar un ordenador es mucho más fácil (aunque sigue siendo bastante lento, y escala como O(n^3) para una matriz de n×n en el caso general).

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

En general, una matriz de n×n tendrá n valores propios, aunque no siempre distintos. Piensa en ellos como las raíces de un polinomio de orden n (llamado polinomio característico), que pueden repetirse y a menudo son complejas incluso para una matriz con valores reales.

Cada vector propio representa una dirección, y cualquier múltiplo escalar también es un vector propio válido. Por comodidad en cálculos posteriores, Julia devuelve vectores unitarios con una norma de 1.

Aplicaciones

Explicar por qué los vectores propios son importantes es una tarea más propia de un libro de texto de 500 páginas que de unos pocos párrafos. Este tema está presente en muchísimas áreas de las matemáticas aplicadas modernas.

A grandes rasgos, los vectores propios representan los ejes (direcciones) «más importantes» de un conjunto de datos. Los valores propios indican la «importancia relativa» de cada eje (sujeto a algunas suposiciones sobre la normalización de las entradas).

Lo que esto significa en realidad depende de la aplicación.

Análisis de componentes principales

En un problema típico de ciencia de datos, podemos tener 100 o más «características», almacenadas como columnas de datos. Casi inevitablemente, habrá ruido, redundancia y correlaciones no deseadas.

Necesitamos reducir la dimensionalidad, y el ACP es una forma de poner orden en el caos:

• Calcula la matriz de covarianza de los datos.
• Ahora tenemos una matriz cuadrada, así que a continuación calcula los valores y vectores propios.
• Ordena los vectores propios en orden descendente según sus valores propios (más concretamente, según el valor absoluto |λ|).
• Los k primeros vectores propios son ahora tus principal components, donde k es bastante menor que el número original de características.
• Proyecta el conjunto de datos completo sobre estos k ejes y empieza a buscar patrones interesantes.

En esta forma, el ACP lo han usado grupos tan diversos como científicos médicos que intentan mejorar las propiedades moleculares en el diseño de fármacos, y activistas políticos que intentan entender las preferencias de los votantes a partir de datos de encuestas. Probablemente también por grupos de marketing que intentan venderte cosas, pero cualquier tecnología puede usarse para bien o para mal.

El ACP no forma parte de una instalación mínima de Julia, pero (fuera de Exercism) el paquete MultivariateStats contiene lo que necesitas.

Procesamiento de imágenes

Llevando la idea del ACP un paso más allá, las imágenes digitales no son más que matrices de valores de píxeles, y podemos calcular sus componentes principales.

Esto se ha usado durante décadas en la compresión de imágenes, donde el ACP sirve de guía sobre qué es más importante conservar al reducir el tamaño del archivo.

Cada vez más, el ACP es una parte esencial de los clasificadores de imágenes, como el reconocimiento facial. La próxima vez que pases por un aeropuerto, el Gran Hermano no solo te observa, sino que también usa álgebra lineal para entender lo que ve.

Ingeniería mecánica

De forma menos polémica, los ejes principales de un componente mecánico son los vectores propios de su tensor de momento de inercia I (una matriz con un nombre ligeramente distinto).

Cada rueda de tu coche probablemente tiene uno o más contrapesos, ajustados para anular los elementos fuera de la diagonal de este tensor I. El mecánico del taller sin duda evita hacer los cálculos (a diferencia de los diseñadores de aviones y los ingenieros de cohetes), pero «vibración» no es más que una palabra cotidiana para describir los términos cruzados del tensor, que en este caso harían tu viaje menos cómodo y aumentarían el desgaste mecánico de los rodamientos.