6.13. Álgebra lineal¶

El álgebra lineal en una cámara es trabajo con matrices pequeñas: una rotación 3x3 que fusiona una muestra de la IMU en el marco del mundo, una matriz de calibración que rectifica una lente, la actualización del estado y la covarianza de un filtro de Kalman, un ajuste polinómico cuyas ecuaciones normales se resuelven como un pequeño sistema lineal. Los submódulos numpy.linalg y scipy.linalg cubren exactamente esa escala.

Las funciones residen en dos módulos. Las operaciones de producto de matrices están en el nivel superior de numpy; las descomposiciones y la inversa de matriz están bajo numpy.linalg; los solucionadores dedicados de sistemas lineales están bajo scipy.linalg. Una multiplicación de matrices 2x2, por ejemplo:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. Qué está disponible¶

dot() – producto de matrices o de vectores.
cross() – producto vectorial 3-D.
trace() – suma de la diagonal.
inv() – inversa de matriz.
det() – determinante.
cholesky() – descomposición de Cholesky (entrada simétrica definida positiva).
eig() – valores y vectores propios de una matriz real simétrica.
norm() – norma 2 de un vector o una matriz.
qr() – descomposición QR con mode='reduced' (por defecto) o mode='complete'.
solve_triangular() – resuelve A @ x = b cuando A es triangular.
cho_solve() – resuelve A @ x = b dado un factor de Cholesky de A.

6.13.2. dot, cross, trace¶

dot() es como se escribe la multiplicación de matrices en la cámara. El operador @ que el numpy de escritorio proporciona para la misma tarea no está implementado en ndarray, por lo que cada producto matriz-vector, matriz-matriz y producto escalar de vectores pasa por la llamada a dot():

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

El resultado de dot() siempre es de dtype float. cross() es el producto vectorial de dos vectores de 3 componentes y trace() la suma de la diagonal principal de una matriz cuadrada.

6.13.3. inv y det¶

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

La inversa se calcula mediante eliminación de Gauss-Jordan, por lo que inv() lanza ValueError cuando la matriz es singular (una entrada diagonal se vuelve cero durante la eliminación). El coste de RAM es aproximadamente el doble del tamaño de la entrada.

El determinante reutiliza la misma eliminación: el tiempo de ejecución es esencialmente el mismo que el de la inversa.

Cuando el objetivo de la aplicación es resolver un sistema lineal, no inviertas y multipliques: prefiere los solucionadores dedicados que se muestran a continuación. Ambos son más rápidos y se comportan mejor numéricamente.

6.13.4. cholesky¶

Para una matriz simétrica definida positiva A, cholesky() devuelve una L triangular inferior tal que A = L @ L.T:

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

Si la entrada no es definida positiva o no es simétrica, se lanza ValueError.

El factor de Cholesky requiere la mitad del trabajo de una factorización LU y es el punto de partida adecuado para cualquier problema cuya matriz se sabe que es simétrica definida positiva (actualizaciones de covarianza, ecuaciones normales de un ajuste por mínimos cuadrados).

6.13.5. eig¶

eig() funciona solo con matrices reales simétricas. Las matrices no simétricas lanzan ValueError. Devuelve una tupla de 2 elementos (eigenvalues, eigenvectors):

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Notas:

Los valores propios se devuelven sin un orden particular. Aplica sort() (y la misma permutación a los vectores propios mediante argsort()) cuando importe un orden ordenado.
Un vector propio es único solo salvo un escalar distinto de cero, por lo que el signo de cada vector propio no está definido de forma única. Dos ejecuciones correctas pueden producir vectores de signo opuesto; esto es inofensivo.

6.13.6. norm¶

La norma euclídea (de Frobenius) de un vector o una matriz:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

La palabra clave opcional axis= calcula la norma a lo largo de un único eje en lugar de sobre todo el array.

6.13.7. qr¶

qr() factoriza una matriz rectangular A (forma (M, N)) en una Q ortonormal y una R triangular superior tales que A == Q @ R:

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

La descomposición se implementa mediante rotaciones de Givens sucesivas. Es la opción adecuada para problemas de mínimos cuadrados donde la matriz no es simétrica.

6.13.8. Resolución de sistemas¶

Los dos solucionadores dedicados bajo ulab.scipy.linalg son a la vez más rápidos y más precisos que np.dot(np.linalg.inv(A), b):

solve_triangular(a, b, lower=False)() – resuelve a @ x = b suponiendo que a es triangular:

A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

cho_solve(L, b)() – dado un factor de Cholesky L, resuelve A @ x = b donde A = L @ L.T:
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

Recurre a estos en lugar de invertir siempre que la estructura de A lo permita: ahorran trabajo de eliminación y la inversa explícita.

6.13.9. Resolución de un sistema lineal pequeño¶

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

El mismo problema puede expresarse factorizando A y llamando al solucionador adecuado: más rápido y más preciso cuando A tiene la estructura adecuada.

Para la referencia completa a nivel de argumentos, consulta numpy.linalg — Rutinas de álgebra lineal y scipy.linalg — Rutinas de álgebra lineal.