Resolución de ecuaciones lineales en

Acerca de Resolución de ecuaciones lineales

¿Recuerdas cuando aprendiste álgebra por primera vez, en la secundaria?

Por lo general, nos daban un sistema de ecuaciones y nos pedían resolver para x y y.

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

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

Mira de nuevo 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 del 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 igual a cero.

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, y no aporta información nueva. Cualquier valor donde x = 0.75y será una solución.

Estos valores (x, y) se encuentran a lo largo de una recta en el espacio 2-D. Hay un número infinito de soluciones, y estas 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 solo punto, elegido de modo que esté normalizado a un vector unitario, pero cualquier múltiplo escalar de este vector también está en el espacio nulo.

Resolver A x = b

¿Y si ahora tienes 1000 ecuaciones con 1000 incógnitas? Si eso te suena absurdo, recuerda que los transductores digitales inalámbricos ahora son baratos y versátiles (miden deformación, velocidad del viento, aceleración en 3 ejes, lo que sea...). Tener 1000 de ellos monitoreando una estructura moderna como un puente colgante es perfectamente razonable. Se necesita un sistema de control que interprete el flujo de datos, listo para emitir una alerta si las cosas se vuelven alarmantes.

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.

Lo primero que se le ocurre a cualquiera que no esté familiarizado con esto es «dividir todo entre A» 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 (¡lentamente!).

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 el solucionador (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 prestar atención.

• Los coeficientes deben alinearse en columnas, así que inserta ceros según sea necesario para los términos que falten.
• Las ecuaciones (y por lo 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/columnas linealmente independientes, así que procura que este sea igual al número de variables desconocidas.

julia> rank(A)
2

Matrices rectangulares

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

Como con la mayoría de las cosas 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).

Cuando hay «muy pocas» filas

Con N-1 ecuaciones (filas de la matriz A), el problema está subdeterminado. ¿Nos rendimos y caemos en la desesperación?

¡Depende!

Las N-1 ecuaciones todavía contienen mucha información. En términos geométricos, necesitamos N ecuaciones para determinar la solución como un punto en un espacio de N dimensiones, pero N-1 nos dirá que la solución debe estar en algún lugar sobre 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 soluciones.

¿Tu recta de soluciones pasa por una región del espacio de parámetros que debería preocupar a un ingeniero respetable? ¿O que podría justificar cerrar la estructura al acceso público (los puentes colgantes se ponen animados con vientos cruzados fuertes)? ¡Quizás vale 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! Extraoficialmente, 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 otras 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 manera similar, N-2 ecuaciones restringirán las soluciones (infinitas) a un plano, y se aplican los mismos principios.

Cuando hay «demasiadas» filas

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

En la ingeniería del mundo real, esto es completamente 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 simple 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), lo que 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 soluciones para A x = b, equivalente 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 específico, 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, tenemos que vivir con palabras que cambian de idioma a mitad de camino.

Por lo general, a los estudiantes se les enseña a calcular valores y vectores propios de matrices 2×2 a mano, pero usar una computadora es mucho más fácil (aunque sigue siendo bastante lento, y escala como O(n^3) para una matriz 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 n×n tendrá n valores propios, aunque no siempre distintos. Piensa en ellos como las raíces de un polinomio de grado n (llamado polinomio característico), que pueden repetirse y a menudo son complejas incluso para una matriz de 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, idealmente, tarea de un libro de texto de 500 páginas, y no de unos pocos párrafos breves. Este tema está presente en tantísimas partes de la matemática aplicada moderna.

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, podríamos tener 100 «características» o más, almacenadas como columnas de datos. Casi inevitablemente, habrá ruido, redundancia y correlaciones no deseadas.

Necesitamos reducir la dimensión, y el PCA es una forma de poner orden en el caos:

• Calcular la matriz de covarianza de los datos.
• Ahora tenemos una matriz cuadrada, así que a continuación hay que calcular los valores y vectores propios.
• Ordenar los vectores propios en orden descendente según sus valores propios (más precisamente, el valor absoluto |λ|).
• Los k vectores propios principales son ahora tus principal components, donde k es mucho menor que el número original de características.
• Proyectar el conjunto de datos completo sobre estos k ejes y empezar a buscar patrones interesantes.

En esta forma, el PCA ha sido usado por 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 PCA 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 PCA un paso más allá, las imágenes digitales no son más que matrices de valores de píxeles, y podemos calcular componentes principales para ellas.

Esto se ha usado durante décadas en la compresión de imágenes, donde el PCA es una guía de lo que es más importante preservar al reducir el tamaño del archivo.

Cada vez más, el PCA es una parte vital de los clasificadores de imágenes, como el reconocimiento facial. La próxima vez que camines por un aeropuerto, el Gran Hermano no solo te está observando: 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 un poco distinto).

Cada rueda de tu auto probablemente tenga uno o más contrapesos de balanceo, ajustados para anular los elementos fuera de la diagonal de este tensor I. El mecánico del taller sin duda evita hacer las cuentas (a diferencia de los diseñadores de aeronaves y los ingenieros de cohetes), pero «bamboleo» es solo 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.