Matrices en Python : produit, puissance, inverse et système

En Python, une matrice se représente par une liste de listes : une liste par ligne, chaque ligne contenant ses coefficients. À partir de là, quelques fonctions de dix lignes suffisent pour additionner, multiplier, élever à une puissance, inverser exactement une matrice et résoudre un système AX = B. Ce chapitre construit ces fonctions avec Python seul (sans bibliothèque à installer), dans l'esprit de l'E3, où l'on te demande d'écrire, de tracer ou de compléter un algorithme (en CCF, le sujet te rappelle alors la méthode mathématique), et montre comment les tracer et les tester.

Représenter une matrice

Q = [[2, 3, 1],
     [4, 2, 5]]

print(len(Q))      # nombre de lignes : 2
print(len(Q[0]))   # nombre de colonnes : 3
print(Q[1][2])     # ligne 2, colonne 3 en maths : 5

Une matrice colonne est une liste de listes à un élément : [[149], [25], [59]], et non [149, 25, 59], qui serait une simple liste (un vecteur au sens de Python, sans dimensions de matrice). Garder la forme liste de listes permet d'utiliser la même fonction produit partout.

Les bases du langage (variables, boucles, fonctions, listes) sont dans le cours de Python. Parcourir un tableau à deux dimensions et construire une matrice avec des boucles imbriquées sont expliqués dans le cours d'algorithmique : tableaux à deux dimensions, matrices et dictionnaires.

Un module de fonctions

Voici les fonctions de base, à enregistrer dans un fichier matrices.py. Chacune correspond à une définition des chapitres précédents.

from fractions import Fraction


def dimensions(A):
    """Renvoie (nombre de lignes, nombre de colonnes)."""
    return len(A), len(A[0])


def afficher(A):
    """Affiche une matrice ligne par ligne, colonnes alignées."""
    largeur = max(len(str(x)) for ligne in A for x in ligne)
    for ligne in A:
        print("[ " + "  ".join(str(x).rjust(largeur) for x in ligne) + " ]")


def nulle(m, n):
    """Matrice nulle de dimensions m × n."""
    return [[0] * n for _ in range(m)]


def identite(n):
    """Matrice identité d'ordre n."""
    return [[1 if i == j else 0 for j in range(n)] for i in range(n)]


def somme(A, B):
    """A + B, case par case (mêmes dimensions)."""
    if dimensions(A) != dimensions(B):
        raise ValueError("somme impossible : dimensions différentes")
    m, n = dimensions(A)
    return [[A[i][j] + B[i][j] for j in range(n)] for i in range(m)]


def fois_reel(k, A):
    """k × A : chaque coefficient multiplié par k."""
    return [[k * x for x in ligne] for ligne in A]

Le produit : trois boucles imbriquées

Le coefficient cᵢⱼ est la somme des aᵢₖ × bₖⱼ pour k allant de 1 à p. En Python, k va de 0 à p − 1, et il faut une troisième boucle pour faire la somme :

def produit(A, B):
    """A × B, ligne par colonne."""
    m, p = dimensions(A)
    p2, n = dimensions(B)
    if p != p2:
        raise ValueError("produit impossible : colonnes de A != lignes de B")
    C = nulle(m, n)
    for i in range(m):
        for j in range(n):
            for k in range(p):
                C[i][j] = C[i][j] + A[i][k] * B[k][j]
    return C

Les deux boucles extérieures parcourent les cases du résultat ; la boucle intérieure calcule la somme de la case (i, j). La vérification des dimensions reproduit la condition du chapitre 3 : avec produit(Q, Q), le programme s'arrête sur le message « produit impossible ».

P2 = [[149, 139],
      [25, 29],
      [59, 49]]
afficher(produit(Q, P2))

[ 432  414 ]
[ 941  859 ]

On retrouve les coûts par service et par fournisseur du chapitre 3.

La puissance

def puissance(A, n):
    """A à la puissance n (A carrée, n >= 0), par produits successifs."""
    R = identite(len(A))
    for _ in range(n):
        R = produit(R, A)
    return R


afficher(puissance([[1, 2], [0, 3]], 3))

[  1  26 ]
[  0  27 ]

On part de l'identité (A⁰ = I) et l'on multiplie n fois par A. Pour n = 0, la boucle ne tourne pas et la fonction renvoie bien I.

Avec la matrice des mots binaires du chapitre 4, puissance([[1, 1], [1, 0]], 30)[0][1] donne 832040, le 30e terme de la suite de Fibonacci, en une fraction de seconde.

L'inverse exact avec des fractions

Pour inverser une matrice, Python a besoin d'un algorithme. Celui des calculatrices est la méthode de Gauss-Jordan : on écrit A et l'identité côte à côte, [A | I], puis on combine les lignes jusqu'à transformer A en I ; la partie de droite est alors devenue A⁻¹. La méthode est hors programme à la main (chapitre 6), mais tu peux utiliser la fonction suivante comme un outil, au même titre que la touche x⁻¹ de ta calculatrice.

Le module standard fractions fournit le type Fraction, qui calcule avec des fractions exactes : Fraction(1, 3) + Fraction(1, 6) vaut exactement Fraction(1, 2), sans les arrondis des nombres à virgule.

def inverse(A):
    """Inverse exact de A (fractions), ou None si A n'est pas inversible.
    Méthode de Gauss-Jordan : on transforme [A | I] en [I | A⁻¹]."""
    n = len(A)
    M = [[Fraction(x) for x in A[i]] + [Fraction(int(i == j)) for j in range(n)]
         for i in range(n)]
    for c in range(n):
        pivot = next((r for r in range(c, n) if M[r][c] != 0), None)
        if pivot is None:
            return None
        M[c], M[pivot] = M[pivot], M[c]
        p = M[c][c]
        M[c] = [x / p for x in M[c]]
        for r in range(n):
            if r != c:
                f = M[r][c]
                M[r] = [a - f * b for a, b in zip(M[r], M[c])]
    return [ligne[n:] for ligne in M]

Pour chaque colonne c, la fonction cherche une ligne dont le coefficient dans cette colonne n'est pas nul (le pivot), la place en position c, la divise pour obtenir un 1, puis retranche le bon multiple de cette ligne à toutes les autres pour obtenir des 0 dans la colonne. Si aucune ligne ne convient, la matrice n'est pas inversible et la fonction renvoie None.

A = [[1, 2, 0],
     [0, 1, 3],
     [1, 0, 2]]
A_inv = inverse(A)
afficher(A_inv)
print(produit(A, A_inv) == identite(3))
print(inverse([[2, 4], [1, 2]]))

[  1/4  -1/2   3/4 ]
[  3/8   1/4  -3/8 ]
[ -1/8   1/4   1/8 ]
True
None

On retrouve l'inverse calculé au chapitre 5, et la vérification A × A⁻¹ = I se fait par une simple comparaison de listes : Fraction(1) est égale à l'entier 1. La matrice [[2, 4], [1, 2]], qui n'a pas d'inverse, donne None.

Résoudre un système AX = B

def resoudre(A, B):
    """Solution X de AX = B, soit X = A⁻¹ × B (A inversible)."""
    A_inv = inverse(A)
    if A_inv is None:
        raise ValueError("A n'est pas inversible")
    return produit(A_inv, B)


commandes = [[5, 10, 3],
             [2, 4, 6],
             [8, 0, 2]]
totaux = [[1172], [752], [1310]]
afficher(resoudre(commandes, totaux))

[ 149 ]
[  25 ]
[  59 ]

C'est le problème des prix unitaires du chapitre 6 : un écran à 149 €, un clavier à 25 €, un casque à 59 €. L'ordre produit(A_inv, B) est celui de la formule X = A⁻¹B.

Tester ses fonctions

Une fonction de calcul se teste sur des cas dont on connaît la réponse. Quelques assertions suffisent à attraper la plupart des erreurs d'indices :

assert produit([[1, 3], [2, 4]], [[2, 0], [1, 5]]) == [[5, 15], [8, 20]]
assert puissance([[1, 2], [0, 3]], 0) == identite(2)
assert puissance([[2, 0], [0, 5]], 3) == [[8, 0], [0, 125]]
assert produit([[2, 1], [5, 3]], inverse([[2, 1], [5, 3]])) == identite(2)
print("Tous les tests passent.")

Si une assertion échoue, Python s'arrête sur une AssertionError qui indique la ligne : c'est le cas à reprendre à la main.

Exercices

def produit_faux(A, B):
    m, p = dimensions(A)
    n = len(B[0])
    C = nulle(m, n)
    for i in range(m):
        for j in range(n):
            for k in range(p):
                C[i][j] = A[i][k] * B[k][j]
    return C

def est_inverse(A, B):
    """Renvoie True si B est l'inverse de A."""
    return produit(A, B) == identite(len(A))


print(est_inverse([[2, 1], [5, 3]], [[3, -1], [-5, 2]]))   # True
print(est_inverse([[2, 1], [5, 3]], [[3, 1], [5, 2]]))     # False

X = produit(puissance([[1, 1], [1, 0]], 9), [[1], [1]])
print(X)                    # [[89], [55]]
print(X[0][0] + X[1][0])    # 144