Chapitre 10
Algèbre linéaire appliquée

Nous entamons la Partie B : les applications. Jusqu'ici on a bâti des outils (types, ownership, modules). À partir de maintenant, on résout de vrais problèmes de robotique, et le premier langage de la robotique, c'est l'algèbre linéaire. Position d'un bras articulé, orientation d'un drone, fusion de capteurs, filtre de Kalman : tout se ramène à des vecteurs et des matrices. Ce chapitre pose des bases mathématiquement rigoureuses puis les traduit en Rust correct.

1. Rappels : vecteurs, matrices, notations

Un vecteur de dimension n est une liste ordonnée de n réels. On le note en colonne, en gras minuscule, avec ses composantes indicées à partir de 1 (convention mathématique ; en Rust on indexera à partir de 0) :

Une matrice m×n (m lignes, n colonnes) est un tableau rectangulaire de réels, noté en gras majuscule. L'élément à la ligne i et la colonne j est A[i][j], ou en notation mathématique a_ij :

Un vecteur colonne de dimension n est simplement une matrice n×1. La dimension est capitale : la plupart des erreurs (et des panics) viennent de dimensions incompatibles. En robotique on rencontre surtout :

• Vecteurs 3D (ℝ³) : position, vitesse, force, axe de rotation.
• Matrices 3×3 : matrices de rotation, tenseur d'inertie, matrice de covariance.
• Matrices 4×4 : transformations homogènes (rotation + translation), voir chapitre 11.

Pourquoi l'algèbre linéaire est-elle partout ?

Deux grandes raisons. D'abord les transformations géométriques : faire tourner un repère, composer les articulations d'un bras, projeter un point de caméra — ce sont des produits matrice-vecteur. Ensuite l'estimation d'état : un filtre (Kalman, moindres carrés) manipule des vecteurs d'état et des matrices de covariance à chaque cycle. Une boucle de contrôle à 1 kHz peut donc effectuer des milliers de produits matriciels par seconde : il faut que ce soit à la fois juste et rapide.

2. Opérations de base

L'addition de deux vecteurs (ou matrices) de mêmes dimensions se fait composante par composante. Le produit par un scalaire λ multiplie chaque composante :

Le produit scalaire (dot product) de deux vecteurs de même dimension renvoie un réel. C'est la brique de la projection et de la mesure d'angle :

La norme euclidienne (longueur) d'un vecteur se déduit du produit scalaire de v avec lui-même :

3. Produit matrice-vecteur et matrice-matrice

Le produit matrice-vecteur A·v n'est défini que si le nombre de colonnes de A égale la dimension de v. Si A est m×n et v ∈ ℝⁿ, alors A·v ∈ ℝᵐ. Chaque composante du résultat est le produit scalaire d'une ligne de A avec v :

Le produit matrice-matrice C = A·B n'est défini que si le nombre de colonnes de A égale le nombre de lignes de B. Si A est m×p et B est p×n, alors C est m×n. L'élément c_ij est le produit scalaire de la ligne i de A par la colonne j de B :

Complexité. Un produit matrice-vecteur m×n coûte O(m·n) multiplications. Un produit matrice-matrice m×p par p×n coûte O(m·p·n) : pour des matrices n×n, c'est O(n³) avec l'algorithme naïf. Pour une 3×3, c'est 27 multiplications — négligeable ; mais dans une boucle temps réel, ces constantes s'additionnent.

4. Transposée, identité, déterminant, inverse

La transposée Aᵀ échange lignes et colonnes : l'élément (i, j) de Aᵀ est l'élément (j, i) de A. Une matrice m×n devient n×m.

La matrice identité I_n est carrée, avec des 1 sur la diagonale et des 0 ailleurs. Elle est neutre pour le produit : A·I = I·A = A.

Le déterminant det(A) est un scalaire défini pour les matrices carrées. Il mesure le facteur d'échelle des volumes sous la transformation. Fait clé : A est inversible si et seulement si det(A) ≠ 0.

La matrice inverse A⁻¹ (quand elle existe) vérifie A·A⁻¹ = A⁻¹·A = I. Elle « annule » la transformation. Formule 2×2 explicite (l'inverse existe ssi det ≠ 0) :

Pour la 3×3, l'inverse vaut (1/det) fois la transposée de la matrice des cofacteurs (comatrice). C'est calculable à la main mais fastidieux et numériquement fragile — d'où la section suivante.

5. Résoudre plutôt qu'inverser

On veut souvent trouver x tel que A·x = b (par exemple : quelles vitesses articulaires donnent telle vitesse d'effecteur ?). La tentation est d'écrire x = A⁻¹·b. Mauvaise idée en pratique.

• Coût. Calculer A⁻¹ puis A⁻¹·b fait plus d'opérations que résoudre directement.
• Précision. L'inversion explicite accumule des erreurs d'arrondi et se comporte très mal si A est mal conditionnée (presque singulière). Une résolution par décomposition (LU, QR, Cholesky) est plus stable.
• Robustesse. Les bibliothèques exposent des solveurs (lu().solve(&b)) conçus exactement pour ça.

Règle : n'inversez une matrice que si vous avez réellement besoin de la matrice inverse elle-même (rare). Pour résoudre un système, utilisez un solveur. On y revient au chapitre sur la cinématique inverse.

6. Approche 1 — une Mat3 « maison »

Pour comprendre ce qui se passe sous le capot, implémentons une matrice 3×3 avec un tableau [[f64; 3]; 3] — trois lignes de trois f64. On choisit f64 et pas f32 : la précision double limite l'accumulation d'erreurs dans les produits et les déterminants.

/// Matrice 3×3 stockée en lignes : data[ligne][colonne].
#[derive(Debug, Clone, Copy, PartialEq)]
struct Mat3 {
    data: [[f64; 3]; 3],
}

impl Mat3 {
    /// Constructeur direct à partir d'un tableau 3×3.
    fn new(data: [[f64; 3]; 3]) -> Self {
        Mat3 { data }
    }

    /// Matrice identité I₃.
    fn identity() -> Self {
        Mat3 {
            data: [
                [1.0, 0.0, 0.0],
                [0.0, 1.0, 0.0],
                [0.0, 0.0, 1.0],
            ],
        }
    }

    /// Produit matrice-vecteur : A·v, avec v ∈ ℝ³.
    /// Chaque composante = produit scalaire d'une ligne par v.
    fn mul_vec(&self, v: [f64; 3]) -> [f64; 3] {
        let mut out = [0.0; 3];
        for i in 0..3 {
            for j in 0..3 {
                out[i] += self.data[i][j] * v[j];
            }
        }
        out
    }

    /// Produit matrice-matrice : C = A·B.
    /// c[i][j] = Σ_k a[i][k] · b[k][j]
    fn mul_mat(&self, other: &Mat3) -> Mat3 {
        let mut out = [[0.0; 3]; 3];
        for i in 0..3 {
            for j in 0..3 {
                for k in 0..3 {
                    out[i][j] += self.data[i][k] * other.data[k][j];
                }
            }
        }
        Mat3 { data: out }
    }

    /// Transposée : (Aᵀ)[i][j] = A[j][i].
    fn transpose(&self) -> Mat3 {
        let mut out = [[0.0; 3]; 3];
        for i in 0..3 {
            for j in 0..3 {
                out[i][j] = self.data[j][i];
            }
        }
        Mat3 { data: out }
    }

    /// Déterminant 3×3 (règle de Sarrus / cofacteurs 1re ligne).
    fn det(&self) -> f64 {
        let m = &self.data;
        m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])
    }

    /// Inverse : renvoie None si la matrice est (quasi) singulière.
    /// A⁻¹ = (1/det) · comatriceᵀ.
    fn inverse(&self) -> Option {
        let d = self.det();
        // Seuil pour éviter les divisions par un déterminant ~0.
        if d.abs() < 1e-12 {
            return None;
        }
        let inv_d = 1.0 / d;
        let m = &self.data;

        // Cofacteurs, déjà transposés (adjugate).
        let mut out = [[0.0; 3]; 3];
        out[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]) * inv_d;
        out[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]) * inv_d;
        out[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]) * inv_d;
        out[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]) * inv_d;
        out[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]) * inv_d;
        out[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]) * inv_d;
        out[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]) * inv_d;
        out[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]) * inv_d;
        out[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]) * inv_d;
        Some(Mat3 { data: out })
    }
}

fn main() {
    let a = Mat3::new([
        [2.0, 0.0, 1.0],
        [1.0, 3.0, 2.0],
        [1.0, 0.0, 1.0],
    ]);

    println!("det(A) = {}", a.det());        // det(A) = 3
    let v = a.mul_vec([1.0, 2.0, 3.0]);
    println!("A·v = {:?}", v);               // A·v = [5.0, 13.0, 4.0]

    match a.inverse() {
        Some(inv) => {
            // Vérification : A·A⁻¹ doit valoir (numériquement) l'identité.
            let prod = a.mul_mat(&inv);
            println!("A·A⁻¹ ≈ {:?}", prod.data);
        }
        None => println!("Matrice singulière, pas d'inverse."),
    }
}

Notez le retour Option<Mat3> : c'est la façon idiomatique en Rust d'exprimer « peut ne pas exister ». Le C renverrait un code d'erreur ou écraserait un pointeur ; ici, l'absence d'inverse est dans le type, et le compilateur force l'appelant à traiter le cas None.

7. Approche 2 — nalgebra (celle qu'on utilisera)

La Mat3 maison est pédagogique, mais en production on utilise nalgebra. Elle offre des types de taille fixe (Matrix3, Vector3, Matrix4…), des opérateurs surchargés, des décompositions robustes et une inversion numériquement fiable via try_inverse (qui renvoie Option, exactement notre convention).

use nalgebra::{Matrix3, Vector3};

fn main() {
    // Attention : Matrix3::new remplit LIGNE par ligne (row-major à la saisie),
    // même si le stockage interne de nalgebra est column-major.
    let a = Matrix3::new(
        2.0, 0.0, 1.0,
        1.0, 3.0, 2.0,
        1.0, 0.0, 1.0,
    );
    let v = Vector3::new(1.0, 2.0, 3.0);

    // Opérateurs naturels : produit matrice-vecteur, matrice-matrice.
    let av: Vector3<f64> = a * v;          // [5, 13, 4]
    let a2: Matrix3<f64> = a * a;           // produit matriciel
    let at = a.transpose();
    let d = a.determinant();               // 3.0

    println!("A·v = {}", av);
    println!("det(A) = {}", d);

    // Inversion sûre : Option.
    match a.try_inverse() {
        Some(inv) => {
            let check = a * inv;           // ≈ identité
            println!("A·A⁻¹ ≈\n{}", check);
        }
        None => println!("A n'est pas inversible."),
    }

    // Pour RÉSOUDRE A·x = b, on ne fait pas x = a.try_inverse() * b :
    // on utilise une décomposition, plus stable.
    let b = Vector3::new(5.0, 13.0, 4.0);
    if let Some(x) = a.lu().solve(&b) {
        println!("x = {}", x);             // ≈ [1, 2, 3]
    }
}

En quelques lignes on obtient produit, transposée, déterminant, inverse et un solveur LU. C'est plus court, plus rapide (code optimisé, éventuellement vectorisé) et plus sûr numériquement que notre version maison. La Mat3 reste utile pour comprendre ; nalgebra pour produire.

Exercices

fn dot(u: [f64; 3], v: [f64; 3]) -> f64 {
    let mut s = 0.0;
    for i in 0..3 {
        s += u[i] * v[i];
    }
    s
}

/// ‖v‖ = √(v · v)
fn norm(v: [f64; 3]) -> f64 {
    dot(v, v).sqrt()
}

fn main() {
    let u = [1.0, 0.0, 0.0];
    let v = [0.0, 2.0, 0.0];
    println!("u·v = {}", dot(u, v));      // 0  (orthogonaux)
    println!("‖(3,4,0)‖ = {}", norm([3.0, 4.0, 0.0])); // 5
}

#[derive(Debug, Clone, Copy)]
struct Mat2 {
    data: [[f64; 2]; 2],
}

impl Mat2 {
    fn mul_vec(&self, v: [f64; 2]) -> [f64; 2] {
        let mut out = [0.0; 2];
        for i in 0..2 {
            for j in 0..2 {
                out[i] += self.data[i][j] * v[j];
            }
        }
        out
    }
}

fn main() {
    // Rotation de +90° dans le plan.
    let rot = Mat2 { data: [[0.0, -1.0], [1.0, 0.0]] };
    let r = rot.mul_vec([1.0, 0.0]);
    println!("{:?}", r); // [0.0, 1.0]  → l'axe x devient l'axe y
}

impl Mat2 {
    fn det(&self) -> f64 {
        // det = a·d − b·c
        self.data[0][0] * self.data[1][1]
            - self.data[0][1] * self.data[1][0]
    }

    fn inverse(&self) -> Option {
        let d = self.det();
        if d.abs() < 1e-12 {
            return None; // singulière
        }
        let inv_d = 1.0 / d;
        let a = self.data[0][0];
        let b = self.data[0][1];
        let c = self.data[1][0];
        let dd = self.data[1][1];
        // A⁻¹ = (1/det) · [[ d, −b], [−c, a]]
        Some(Mat2 {
            data: [
                [ dd * inv_d, -b * inv_d],
                [ -c * inv_d,  a * inv_d],
            ],
        })
    }
}

fn main() {
    let m = Mat2 { data: [[4.0, 7.0], [2.0, 6.0]] };
    println!("det = {}", m.det());        // 10
    println!("{:?}", m.inverse());        // Some(...)

    let sing = Mat2 { data: [[1.0, 2.0], [2.0, 4.0]] };
    println!("{:?}", sing.inverse());     // None (det = 0)
}

use nalgebra::Matrix2;

fn main() {
    // nalgebra : saisie ligne par ligne.
    let m = Matrix2::new(
        4.0, 7.0,
        2.0, 6.0,
    );

    let inv = m.try_inverse().expect("matrice inversible");
    println!("nalgebra inverse =\n{}", inv);

    // Inverse attendue (calcul maison) :
    //   det = 4·6 − 7·2 = 10
    //   A⁻¹ = 1/10 · [[6, −7], [−2, 4]]
    let attendu = [[0.6, -0.7], [-0.2, 0.4]];

    for i in 0..2 {
        for j in 0..2 {
            let ecart = (inv[(i, j)] - attendu[i][j]).abs();
            assert!(ecart < 1e-12, "écart trop grand en ({i},{j})");
        }
    }
    println!("Les deux inverses coïncident à 1e-12 près.");

    // Contrôle final : A·A⁻¹ ≈ identité.
    println!("A·A⁻¹ =\n{}", m * inv);
}

Récapitulatif