6.13. Algèbre linéaire¶

L’algèbre linéaire sur une caméra concerne les petites matrices : une rotation 3x3 qui intègre un échantillon d’IMU dans le repère du monde, une matrice de calibration qui rectifie un objectif, la mise à jour de la covariance d’état d’un filtre de Kalman, un ajustement polynomial dont les équations normales se résolvent comme un minuscule système linéaire. Les sous-modules numpy.linalg et scipy.linalg couvrent exactement cette échelle.

Les fonctions résident dans deux modules. Les opérations de produit matriciel se trouvent au niveau supérieur de numpy ; les décompositions et l’inverse de matrice sont sous numpy.linalg ; les solveurs de systèmes linéaires dédiés sont sous scipy.linalg. Une multiplication de matrices 2 par 2, par exemple

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. Ce qui est disponible¶

dot() – produit matriciel ou vectoriel.
cross() – produit vectoriel 3D.
trace() – somme de la diagonale.
inv() – inverse de matrice.
det() – déterminant.
cholesky() – décomposition de Cholesky (entrée symétrique définie positive).
eig() – valeurs propres et vecteurs propres d’une matrice réelle symétrique.
norm() – norme 2 d’un vecteur ou d’une matrice.
qr() – décomposition QR avec mode='reduced' (par défaut) ou mode='complete'.
solve_triangular() – résout A @ x = b lorsque A est triangulaire.
cho_solve() – résout A @ x = b à partir d’un facteur de Cholesky de A.

6.13.2. dot, cross, trace¶

dot() est la manière d’écrire la multiplication de matrices sur la caméra. L’opérateur @ que numpy de bureau fournit pour la même tâche n’est pas implémenté sur ndarray, de sorte que tout produit matrice-vecteur, matrice-matrice et produit scalaire de vecteurs passe par l’appel 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

Le résultat de dot() est toujours de dtype float. cross() est le produit vectoriel de deux vecteurs à 3 composantes, et trace() la somme de la diagonale principale d’une matrice carrée.

6.13.3. inv et 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))

L’inverse est calculé par élimination de Gauss-Jordan, donc inv() lève ValueError lorsque la matrice est singulière (une entrée diagonale devient nulle au cours de l’élimination). Le coût en RAM est environ le double de la taille de l’entrée.

Le déterminant réutilise la même élimination – le temps d’exécution est essentiellement le même que pour l’inverse.

Lorsque l’objectif de l’application est de résoudre un système linéaire, n’inversez pas pour ensuite multiplier – préférez les solveurs dédiés ci-dessous. Les deux sont plus rapides et numériquement mieux conditionnés.

6.13.4. cholesky¶

Pour une matrice symétrique définie positive A, cholesky() retourne une matrice triangulaire inférieure L telle que A = L @ L.T

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

Si l’entrée n’est pas définie positive ou n’est pas symétrique, ValueError est levée.

Le facteur de Cholesky représente la moitié du travail d’une factorisation LU et constitue le bon point de départ pour tout problème dont la matrice est connue comme symétrique définie positive (mises à jour de covariance, équations normales d’un ajustement par moindres carrés).

6.13.5. eig¶

eig() ne fonctionne que sur les matrices réelles symétriques. Les matrices non symétriques lèvent ValueError. Elle retourne un 2-uplet (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)

Remarques :

Les valeurs propres reviennent dans un ordre quelconque. Appliquez sort() (et la même permutation aux vecteurs propres via argsort()) lorsqu’un ordre trié est important.
Un vecteur propre n’est unique qu’à un scalaire non nul près, donc le signe des vecteurs propres individuels n’est pas défini de manière unique. Deux exécutions correctes peuvent produire des vecteurs de signe opposé ; cela est sans conséquence.

6.13.6. norm¶

La norme euclidienne (de Frobenius) d’un vecteur ou d’une matrice

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...

Le mot-clé optionnel axis= calcule la norme le long d’un seul axe au lieu de l’ensemble du tableau.

6.13.7. qr¶

qr() factorise une matrice rectangulaire A (de forme (M, N)) en une matrice orthonormale Q et une matrice triangulaire supérieure R telles 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 décomposition est implémentée par rotations de Givens successives. C’est le bon choix pour les problèmes de moindres carrés où la matrice n’est pas symétrique.

6.13.8. Résolution de systèmes¶

Les deux solveurs dédiés sous ulab.scipy.linalg sont à la fois plus rapides et plus précis que np.dot(np.linalg.inv(A), b) :

solve_triangular(a, b, lower=False)() – résout a @ x = b en supposant que a est triangulaire

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)() – à partir d’un facteur de Cholesky L, résout A @ x = b où A = L @ L.T
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

Privilégiez ceux-ci plutôt que l’inversion dès que la structure de A le permet – ils économisent le travail d’élimination et l’inverse explicite.

6.13.9. Résoudre un petit système linéaire¶

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

Le même problème peut s’exprimer en factorisant A et en appelant le solveur approprié – plus rapide et plus précis lorsque A a la bonne structure.

Pour la référence complète au niveau des arguments, consultez numpy.linalg — Routines d’algèbre linéaire et scipy.linalg — Routines d’algèbre linéaire.