Résolution d’équations linéaires en

À propos de Résolution d’équations linéaires

Tu te souviens de tes premiers cours d'algèbre, au lycée ?

En général, on nous donnait un système d'équations et on nous demandait de trouver x et y.

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

Deux équations, deux inconnues : un peu de substitution, et on voit vite que x = 1, y = 2.

Regarde à nouveau le membre de gauche de l'exemple précédent. On dirait une matrice [4 -3; 2 7] qui multiplie un vecteur inconnu [x, y] pour donner le vecteur [-2, 16].

En notation d'algèbre linéaire : A x = b, selon la convention habituelle des majuscules pour les matrices et des minuscules pour les vecteurs.

A x = 0 et le noyau

Un cas particulier important est celui où toutes nos équations ont zéro dans le membre de droite.

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

Dans l'exemple ci-dessus, la seule solution est x = y = 0, ce qui n'est généralement pas très intéressant.

Les solutions non triviales de A x = 0 n'existent que lorsque le déterminant de A est nul.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Prenons plutôt ces équations :

4x - 3y = 0
8x + 6y = 0

La deuxième équation est simplement le double de la première, elle n'apporte aucune information nouvelle. Toutes les valeurs telles que x = 0.75y sont solutions.

Ces valeurs (x, y) se situent sur une droite dans l'espace à 2 dimensions. Il existe une infinité de solutions, et elles forment le noyau de la matrice.

Dans ce cas, les lignes de la matrice ne sont pas linéairement indépendantes, et le rang de la matrice est inférieur au nombre de lignes ou de colonnes.

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 fonction nullspace() renvoie un seul point, choisi pour être normalisé en un vecteur unitaire, mais tout multiple scalaire de ce vecteur appartient lui aussi au noyau.

Résous A x = b

Suppose maintenant que tu as 1000 équations à 1000 inconnues ? Si ça te semble absurde, rappelle-toi que les capteurs numériques sans fil sont aujourd'hui bon marché et polyvalents (mesure de déformation, vitesse du vent, accélération sur 3 axes, et ainsi de suite). En avoir 1000 qui surveillent une structure moderne comme un pont suspendu est parfaitement raisonnable. Il faut un système de contrôle qui interprète le flux de données et se tienne prêt à déclencher une alerte si la situation devient inquiétante.

À nouveau, A x = b, où l'on connaît A (issue de la conception technique) et b (les valeurs mesurées par les capteurs), mais où l'on cherche x.

Une première idée, quand on ne connaît pas le sujet, consiste à « diviser les deux membres par A » d'une manière ou d'une autre pour le déplacer dans le membre de droite.

En fait, A admet un inverse (la plupart des matrices carrées en ont un, mais pas toutes : le déterminant doit être non nul), et le calcul fonctionne (lentement !).

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

Malheureusement, calculer l'inverse est lent pour les petites matrices et carrément glacial pour les grandes. Pas idéal, si ton pont s'effondre pendant le calcul !

Heureusement, d'autres algorithmes sont bien plus rapides (dans ce cas, le pivot de Gauss).

Julia (comme Matlab) utilise simplement une barre oblique inverse pour le solveur (techniquement, la « division à gauche »).

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

C'est un exemple trivial, mais quelques détails méritent l'attention.

• Les coefficients doivent être alignés en colonnes, alors ajoute des zéros au besoin pour les termes manquants.
• Les équations (et donc les lignes de la matrice) doivent être linéairement indépendantes. Si la ligne 3 n'est que la somme des lignes 1 et 2, elle n'apporte aucune information nouvelle et le problème est sous-déterminé.

La commande rank() donne le nombre de lignes ou de colonnes linéairement indépendantes ; cherche donc à ce qu'il soit égal au nombre d'inconnues.

julia> rank(A)
2

Matrices rectangulaires

Quand tu étais enfant, on t'a probablement appris qu'il faut un système de N équations pour résoudre N inconnues.

Comme pour la plupart des choses en mathématiques, la réalité est un peu plus nuancée (et pas seulement à cause du détail d'indépendance linéaire mentionné dans la section précédente).

« Trop peu » de lignes

Avec N-1 équations (lignes de la matrice A), le problème est sous-déterminé. Faut-il abandonner de désespoir ?

Ça dépend !

Les N-1 équations contiennent encore beaucoup d'informations. En termes géométriques, il faut N équations pour déterminer la solution comme un point dans un espace à N dimensions, mais N-1 nous indiquent que la solution se trouve quelque part sur une droite.

Il nous reste donc une infinité de solutions, mais aussi une infinité d'espace de paramètres sans aucune solution.

Ta droite de solutions traverse-t-elle une région de l'espace des paramètres qui devrait inquiéter un ingénieur sérieux ? Ou qui justifierait de fermer la structure au public (les ponts suspendus peuvent s'agiter sous l'effet de vents latéraux violents) ? Ça vaut peut-être la peine de faire quelques calculs supplémentaires pour le déterminer !

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

Le solveur à barre oblique inverse nous donne une solution ! Non officiellement, il s'agit apparemment de la solution de plus petite norme (la plus proche de l'origine, géométriquement).

Pour obtenir les autres solutions, il faut le nullspace de A2.

Ajouter à la solution d'origine n'importe quel multiple scalaire du noyau donne une autre solution valide.

# 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 même, N-2 équations contraignent les solutions (infinies) à un plan, et les mêmes principes s'appliquent.

« Trop » de lignes

La situation inverse se produit lorsqu'il y a plus d'équations que d'inconnues, et qu'elles restent pourtant linéairement indépendantes.

Dans l'ingénierie réelle, c'est parfaitement normal, et même considéré comme une bonne chose !

Les capteurs ont une précision limitée, le vecteur b a une précision limitée, et le calcul comporte du bruit. La solution est alors un ajustement par moindres carrés aux données bruitées.

La technique la plus simple utilise la pseudoinverse de la matrice, que Julia implémente sous forme de fonction pinv().

On peut l'utiliser un peu comme l'inverse d'une matrice carrée (non singulière), ce qui donne une estimation des variables au sens des moindres carrés plutôt qu'une solution exacte.

# 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

Valeurs propres et vecteurs propres

La section précédente décrivait les solutions de A x = b, ce qui correspond à de l'algèbre familière.

Cette section traite des solutions de A x = λ x, où λ est un scalaire.

C'est un concept caractéristique de l'algèbre linéaire, mais on peut l'interpréter géométriquement :

Pour des valeurs appropriées de x et de λ, la matrice carrée A multiplie la longueur du vecteur non nul x par un facteur λ, sans changer sa direction (sauf pour un λ négatif, qui l'inverse).

Ça peut sembler anecdotique, mais c'est en réalité ridiculement utile !

Terminologie : les valeurs valides de λ sont les valeurs propres de A, et les valeurs correspondantes de x sont les vecteurs propres de A. Malheureusement, il faut bien faire avec des mots qui changent de langue en cours de route.

On apprend généralement aux élèves à calculer à la main les valeurs propres et les vecteurs propres d'une matrice 2×2, mais un ordinateur simplifie grandement les choses (même si le calcul reste assez lent, avec une complexité en O(n^3) pour une matrice n×n dans le cas général).

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 général, une matrice n×n possède n valeurs propres, qui ne sont pas toujours distinctes. Vois-les comme les racines d'un polynôme de degré n (appelé polynôme caractéristique), qui peuvent être multiples et sont souvent complexes, même pour une matrice à coefficients réels.

Chaque vecteur propre représente une direction, et tout multiple scalaire est aussi un vecteur propre valide. Pour faciliter les calculs ultérieurs, Julia renvoie des vecteurs unitaires de norme 1.

Applications

Expliquer pourquoi les vecteurs propres sont importants est idéalement le travail d'un manuel de 500 pages, plutôt que de quelques courts paragraphes. Ce sujet imprègne une grande partie des mathématiques appliquées modernes.

Au premier niveau, les vecteurs propres représentent les axes (directions) « les plus importants » d'un jeu de données. Les valeurs propres indiquent l'« importance relative » de chaque axe (sous réserve de certaines hypothèses sur la normalisation des entrées).

Ce que cela signifie réellement dépend de l'application.

Analyse en composantes principales

Dans un problème classique de science des données, on peut avoir 100 « caractéristiques » ou plus, stockées sous forme de colonnes de données. Presque inévitablement, il y aura du bruit, de la redondance et des corrélations indésirables.

Il faut réduire la dimensionnalité, et l'ACP est un moyen de mettre de l'ordre dans le chaos :

• Calcule la matrice de covariance des données.
• On dispose alors d'une matrice carrée : calcule ensuite les valeurs propres et les vecteurs propres.
• Trie les vecteurs propres par ordre décroissant de leurs valeurs propres (plus précisément, de leur valeur absolue |λ|).
• Les k premiers vecteurs propres sont désormais tes principal components, où k est nettement plus petit que le nombre initial de caractéristiques.
• Projette le jeu de données complet sur ces k axes, puis commence à chercher des motifs intéressants.

Sous cette forme, l'ACP a été utilisée par des groupes aussi divers que des scientifiques cherchant à améliorer les propriétés moléculaires dans la conception de médicaments, ou des militants politiques cherchant à comprendre les préférences des électeurs à partir de données de sondage. Probablement aussi par des équipes marketing qui cherchent à te vendre des choses, mais toute technologie peut servir le bien comme le mal.

L'ACP ne fait pas partie d'une installation minimale de Julia, mais (en dehors d'Exercism) le paquet MultivariateStats contient ce qu'il te faut.

Traitement d'images

En poussant l'idée de l'ACP un peu plus loin, une image numérique n'est qu'une matrice de valeurs de pixels, dont on peut calculer les composantes principales.

Cette approche est utilisée depuis des décennies pour la compression d'images, où l'ACP indique ce qu'il est le plus important de préserver lorsqu'on réduit la taille du fichier.

De plus en plus, l'ACP joue un rôle essentiel dans les classificateurs d'images, comme la reconnaissance faciale. La prochaine fois que tu traverseras un aéroport, Big Brother ne se contentera pas de te regarder : il utilise aussi l'algèbre linéaire pour comprendre ce qu'il voit !

Génie mécanique

De manière moins polémique, les axes principaux d'une pièce mécanique sont les vecteurs propres de son tenseur d'inertie I (une matrice sous un nom légèrement différent).

Chaque roue de ta voiture possède probablement une ou plusieurs masses d'équilibrage, réglées pour annuler les éléments hors diagonale de ce tenseur I. Le mécanicien du garage évite sans aucun doute de faire les calculs (contrairement aux concepteurs d'avions et aux ingénieurs en fusées), mais le « ballottement » n'est qu'un mot du quotidien pour décrire les termes croisés du tenseur, qui, dans ce cas, rendraient ton trajet moins confortable et augmenteraient l'usure mécanique des roulements.