Bases de l’algèbre linéaire en

À propos de Bases de l’algèbre linéaire

Qu'est-ce que « l'algèbre linéaire » ?

Les définitions techniques ne manquent pas : sur Wikipédia, dans les manuels et sur d'excellents sites comme 3Blue1Brown. Mais pour ce qui nous intéresse, on peut se montrer plus décontracté.

Les mainteneurs de Julia (et de Python) adorent ce genre de choses. La direction d'Exercism, beaucoup moins.

Ce concept se limitera aux aspects les plus simples d'un immense sujet, mais attention : c'est inévitablement un concept assez mathématique.

Qu'est-ce qu'un vecteur ?

Tout dépend à qui on pose la question, et de la façon dont on souhaite visualiser une notion assez abstraite.

• En Julia, Vector{T} n'est qu'un alias de Array{T, 1}, et l'eltype T peut être n'importe quoi.
• En physique (du moins en grande partie), un vecteur est une flèche avec une length et une direction dans un espace à N dimensions, mais sans position fixe.
• En algèbre linéaire, la flèche des physiciens est fixée dans l'espace, sa queue à l'origine. La hampe de la flèche devient donc superflue, ce qui permet de représenter le vecteur comme un point dans un espace à N dimensions, les N éléments représentant la distance à l'origine le long de chacun des N axes.

Pour ce qui nous intéresse ici, on ignore les vecteurs de strings ou de caractères. Dans ce concept, on travaillera avec des types numériques : Int, Float ou Complex.

Qu'est-ce qu'une matrice ?

Là encore, les réponses varient.

• Un tableau rectangulaire à deux dimensions de nombres (sans bords irréguliers). Une matrice carrée en est un cas particulier courant.
• Une linear combination de vecteurs colonnes, placés côte à côte.
• Une transformation que l'on peut appliquer à un vecteur, tout comme une function s'applique à d'autres types d'entrées.

Certains types de matrices carrées sont assez courants pour avoir des noms spéciaux.

Matrice diagonale : Tous les éléments non nuls se trouvent sur la main diagonal (d'en haut à gauche jusqu'en bas à droite).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Matrice identité : Une matrice diagonale dont la diagonale ne contient que des uns (pour des raisons qui deviendront plus claires par la suite). Souvent abrégée en I.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Triangulaire supérieure : Des valeurs non nulles sur la diagonale et au-dessus, des zéros en dessous. Triangulaire inférieure, tu devines.

Transposer une matrice

Si tu veux vraiment échanger les lignes et les colonnes, la fonction permutedims() s'en charge et te renvoie une nouvelle matrice.

Cela implique une copie, ce qui est lent et gourmand en mémoire lorsqu'on travaille avec de grandes matrices.

Pour l'algèbre linéaire, la fonction transpose() est plus utile, car elle crée rapidement une enveloppe paresseuse autour de la matrice d'origine.

Encore plus utile, la fonction adjoint(), qui inverse aussi le signe de la partie imaginaire des nombres complexes (si tu te demandes pourquoi c'est si important, une réponse rapide est : la mécanique quantique). Cette opération est assez courante pour qu'on puisse simplement ajouter une apostrophe ' au nom de la variable afin de créer l'adjointe.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Multiplication

Le reste de ce document utilise des fonctionnalités du module LinearAlgebra. Une ligne using LinearAlgebra les fera entrer dans l'espace de noms, mais on ne le répétera pas dans les exemples (trop de bruit visuel).

Vecteurs, produit élément par élément

Cette opération a été présentée dans le concept Opérations sur les vecteurs. L'opérateur est .*, qui agit sur les vecteurs d'entrée deux à deux, et produit un résultat de même taille et de même type que les entrées.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vecteurs, produit scalaire

Cette opération extrêmement courante s'écrit u ⋅ v dans les manuels, où on l'appelle le produit scalaire.

Un produit scalaire équivaut à la somme du produit élément par élément.

On peut multiplier deux vecteurs avec l'opérateur * habituel, mais seulement si le vecteur de gauche est un adjoint : ce que l'on écrit commodément u' * v. Cela deviendra plus clair dans la section sur la multiplication matricielle.

Le point surélevé (centré) est disponible en Julia (saisi avec \cdot puis tabulation) comme sucre syntaxique de la fonction dot(). Nul besoin de préciser l'adjointe, car ce détail est géré automatiquement.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Pour des vecteurs à valeurs complexes, le vecteur de gauche doit être le conjugué (signe inversé sur la partie imaginaire). La fonction dot s'en charge automatiquement, la syntaxe u' le fait explicitement, mais sum(u .* v) échouera si u et v sont complexes.

Vecteurs, produit vectoriel

Un brin de nostalgie pour quiconque a déjà suivi un cours d'électricité et de magnétisme ! De même pour les ingénieurs habitués à calculer des vecteurs de couple ou de moment cinétique.

Alors que le produit scalaire transforme deux vecteurs en un scalar, le produit vectoriel transforme deux vecteurs à 3 composantes en un troisième vecteur à 3 composantes.

Dans la représentation géométrique de l'espace vectoriel, le nouveau vecteur est perpendiculaire au plan contenant les deux vecteurs d'entrée. Si les deux entrées sont parallèles (au signe près), elles ne définissent pas de plan, et le résultat sera donc [0, 0, 0].

L'ordre compte : u × v == -(v × u). C'est la fameuse règle de la main droite, qui nous a tous laissés, un jour ou l'autre, à fixer notre pouce et deux doigts en les tordant dans l'espace (souvent avec une expression perplexe).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

On remarque que cette opération se limite aux vecteurs de longueur 3 (équivalent à l'espace euclidien avec des axes orthogonaux x, y, z). Les mathématiciens définissent volontiers cette opération pour d'autres dimensions, avec pour résultat un salmigondis incompréhensible (tenseurs d'ordre supérieur, produit cunéaire, produit extérieur, multivecteurs...). La plupart d'entre nous prennent leurs jambes à leur cou ou changent vite de sujet à ce moment-là !

Normes

Quelle est la « taille » d'un vecteur ?

La norm est une tentative de capturer cela, en réduisant le vecteur à un scalaire approprié.

Toute une famille de normes est définie, mais de loin la plus courante est la norme 2, qui vaut √(v ⋅ v).

Cette opération de moyenne quadratique correspond à la distance de Pythagore depuis l'origine (dans un espace à N dimensions). Si on visualise le vecteur comme une flèche dont la queue est à l'origine, la norme 2 est la longueur de la flèche.

On peut calculer n'importe quelle norme p en passant p comme deuxième argument. La norme 1 est parfois utile : c'est simplement la somme des valeurs absolues des éléments (donc très rapide et facile à calculer).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Multiplication matricielle

Il semble parfois que tous les mathématiciens appliqués du monde aient passé une bonne partie des 80 dernières années à transformer chaque calcul en une série de multiplications matricielles.

C'est une opération dans laquelle les ordinateurs excellent :

• Elle est très répétitive.
• Elle peut être parallélisée efficacement.
• Beaucoup de matériel spécialisé a été développé pour l'accélérer, du Cray-1 à 80 millions de dollars des années 1970 au GPU probablement intégré à ton ordinateur portable.

Les détails sont assez simples, même s'ils ne sont pas très intuitifs au premier abord.

Considère une matrice A qui multiplie un vecteur v pour donner un résultat w (par convention en algèbre linéaire, on utilise des majuscules pour les matrices et des minuscules pour les vecteurs).

La première ligne de A est multipliée scalairement par v pour obtenir le premier élément de w, la deuxième ligne donne le deuxième élément, et ainsi de suite.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

La multiplication matrice * matrice étend cela à toutes les colonnes de la matrice de droite.

Pour C = A * B, on peut considérer que A multiplie chaque colonne de B pour donner la colonne correspondante de C : une série de multiplications matrice-vecteur.

De façon équivalente, on peut dire que la première ligne de A est multipliée scalairement par chaque colonne de B pour donner la première ligne de C, la deuxième ligne donne la deuxième, et ainsi de suite.

Il existe plusieurs représentations mentales de ce genre, et les apprenants intéressés peuvent regarder un cours du MIT entier qui les explore.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Comme le montre l'exemple ci-dessus, la multiplication matricielle ne commute pas : il n'existe pas de relation simple entre A*B et B*A.

La multiplication matricielle est probablement difficile à visualiser pour qui la découvre, à la simple lecture des mots. YouTube regorge de vidéos qui la démontrent graphiquement : cherche « matrix multiplication » et choisis celle dont le style, le niveau de détail et la langue te conviennent.

Dimensions

Le produit scalaire de deux vecteurs suppose qu'ils aient la même longueur.

Par extension, pour la multiplication matricielle, le nombre de colonnes de la matrice de gauche doit correspondre au nombre de lignes de la matrice de droite.

En exprimant les tailles sous forme de tuples (nrows, ncols), comme les renvoie size(A), on a (a, b) * (b, c) -> (a, c). Les dimensions « internes », ici b et b, sont compatibles pour le produit scalaire. Les dimensions « externes », ici a et c, déterminent les dimensions du résultat.

Un exemple avec des matrices rectangulaires :

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Multiplier des paires de vecteurs, à la manière des matrices, offre deux possibilités.

Par convention, on utiliserait u' * v comme équivalent du produit scalaire. Les dimensions sont (1, 3) * (3, 1), et Julia simplifie le résultat (1, 1) en un scalaire (contrairement, par exemple, à R).

Autrement, on pourrait utiliser u * v', de dimensions (3, 1) * (1, 3) -> (3, 3), pour multiplier toutes les paires d'éléments possibles et transformer les vecteurs en matrice.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v est parfois appelé le produit intérieur.

Le (moins courant) u * v' est quant à lui le produit extérieur. Cela est lié au produit tensoriel, bien que cela dépasse largement le cadre de ce document.

Rotations

Un type particulièrement courant de multiplication matricielle fait intervenir des matrices de rotation.

En 2D, il existe une matrice relativement simple pour faire tourner un vecteur de θ radians dans le sens antihoraire.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

Les rotations en 3D nécessitent une matrice plus complexe, avec deux angles. C'est difficile à représenter dans un document Markdown sans y intégrer beaucoup de LaTeX, alors consulte Wikipédia pour la formule.

Pour la 3D graphique, comme OpenGL et ses successeurs, la convention est de travailler avec des coordonnées homogènes.

Chaque sommet est représenté par un vecteur à 4 composantes (ou, de façon équivalente, une colonne de matrice) : [x, y, z, 1.0] pour un point en (x, y, z).

Cela permet des matrices de transformation plus élaborées. La rotation reste en A[1:3, 1:3], les translations sont [Δx, Δy, Δz] en A[1:3, 4], la mise à l'échelle est sur la diagonale, et d'autres possibilités existent pour l'inclinaison, la perspective, etc.

Performance

Dans les départements d'informatique du monde entier, de nombreuses thèses de doctorat consacrent un ou plusieurs chapitres à de petites améliorations progressives des algorithmes de multiplication matricielle. C'est vraiment important, et diverses organisations sont prêtes à financer ces recherches pour leur propre bénéfice.

La performance est un vaste sujet que nous ne pouvons pas approfondir ici, mais pour illustrer, on peut essayer de multiplier des matrices aléatoires de tailles diverses.

Le code ci-dessous donne une estimation à la va-vite (utilise BenchmarkTools.jl pour une approche plus rigoureuse).

Le système utilisé était un petit PC à 430 $ (US) : processeur Ryzen 9, 32 Go de RAM, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

En gros, une paire de matrices d'un million d'éléments a pris quelques millisecondes, des matrices de 100 millions d'éléments ont pris quelques secondes. Ajouter plusieurs ordres de grandeur supplémentaires demandera un meilleur matériel...