La traduction de cette page n'est pas à jour. Cliquez ici pour voir la dernière version en anglais.

Algèbre linéaire

Matrices dans l’environnement MATLAB

Ce thème comprend une introduction à la création de matrices et à la réalisation de calculs matriciels de base dans MATLAB^®.

L’environnement MATLAB utilise le terme matrice pour indiquer une variable contenant des nombres réels ou complexes arrangés selon une grille à deux dimensions. Un tableau est, plus généralement, un vecteur, une matrice ou une grille de nombres comportant plus de dimensions. Tous les tableaux dans MATLAB sont rectangulaires, dans le sens où les vecteurs composants, quelle que soit la dimension du tableau, sont tous de la même longueur. Les opérations mathématiques définies sur des matrices sont le sujet de l’algèbre linéaire.

Création de matrices

MATLAB comporte de nombreuses fonctions qui créent différents types de matrices. Par exemple, vous pouvez créer une matrice symétrique avec des entrées basées sur le triangle de Pascal :

A = pascal(3)

A =
       1     1     1
       1     2     3
       1     3     6

Ou vous pouvez créer une matrice de carré magique asymétrique, dont les sommes des lignes et des colonnes sont égales :

B = magic(3)

B =
       8     1     6
       3     5     7
       4     9     2

Un autre exemple est une matrice rectangulaire 3 x 2 composée de nombres entiers aléatoires. Dans ce cas, la première entrée de randi décrit la plage des valeurs possibles pour les nombres entiers, et les deux entrées suivantes décrivent le nombre de lignes et de colonnes.

C = randi(10,3,2)

C =

     9    10
    10     7
     2     1

Un vecteur colonne est une matrice m x 1, un vecteur ligne est une matrice de 1 x n et un scalaire est une matrice de 1 x 1. Pour définir une matrice manuellement, utilisez des crochets [ ] pour indiquer le début et la fin de la matrice. À l’intérieur des crochets, utilisez un point-virgule ; pour indiquer la fin d’une ligne. Dans le cas d’un scalaire (matrice de 1 x 1), les crochets ne sont pas requis. Par exemple, les instructions suivantes produisent un vecteur colonne, un vecteur ligne et un scalaire :

u = [3; 1; 4]

v = [2 0 -1]

s = 7

u =
       3
       1
       4

v =
       2     0    -1

s =
       7

Pour plus d’informations sur la création et l’utilisation des matrices, consultez Creating, Concatenating, and Expanding Matrices.

Addition et soustraction de matrices

L’addition et la soustraction de matrices et de tableaux sont effectuées éléments par éléments. Par exemple, ajouter A à B, puis soustraire A du résultat restitue B :

X = A + B

X =
       9     2     7
       4     7    10
       5    12     8

Y = X - A

Y =
       8     1     6
       3     5     7
       4     9     2

L’addition et la soustraction nécessitent que les deux matrices aient des dimensions compatibles. Si les dimensions sont incompatibles, cela génère une erreur :

X = A + C

Error using  + 
Matrix dimensions must agree.

Pour plus d’informations, consultez Array vs. Matrix Operations.

Produits et transposition de vecteurs

Il est possible de multiplier un vecteur ligne et un vecteur colonne de la même longueur dans n’importe quel ordre. Le résultat est soit un scalaire, appelé le produit intérieur, soit une matrice, appelée le produit extérieur :

u = [3; 1; 4];
v = [2 0 -1];
x = v*u

x =

     2

X = u*v

X =

     6     0    -3
     2     0    -1
     8     0    -4

Pour les matrices réelles, latransposée échange a_ij et a_ji. Pour les matrices complexes, un autre élément à prendre en compte est d'utiliser le conjugué complexe des entrées complexes de la matrice pour former la transposée conjuguée complexe. MATLAB utilise l’opérateur apostrophe (') pour effectuer une transposée conjuguée complexe et l’opérateur point-apostrophe (.') pour transposer sans conjugaison. Pour les matrices ne contenant que des éléments réels, les deux opérateurs donnent le même résultat.

La matrice proposée en exemple A = pascal(3) est symétrique, par conséquent, A' est égal à A. Cependant, B = magic(3) n’est pas symétrique, donc les éléments de B' se trouvent le long de la diagonale principale :

B = magic(3)

B =

     8     1     6
     3     5     7
     4     9     2

X = B'

X =

     8     3     4
     1     5     9
     6     7     2

Pour les vecteurs, la transposée change un vecteur ligne en un vecteur colonne (et inversement) :

x = v'

x =
       2
       0
      -1

Si x et y sont tous deux des vecteurs colonnes réels, alors le produit x*y n’est pas défini, mais les deux produits

x'*y

et

y'*x

produisent le même résultat scalaire. Cette quantité est utilisée tellement souvent qu’elle porte trois noms différents : produit intérieur, produit scalaire ou produit point. Il existe même une fonction dédiée pour les produits points appelée dot.

Pour un vecteur ou une matrice complexe, z, la quantité z' transpose non seulement le vecteur ou la matrice, mais convertit également chaque élément complexe en son conjugué complexe. Autrement dit, le signe de la partie imaginaire de chaque élément complexe change. Par exemple, soit la matrice complexe suivante

z = [1+2i 7-3i 3+4i; 6-2i 9i 4+7i]

z =

   1.0000 + 2.0000i   7.0000 - 3.0000i   3.0000 + 4.0000i
   6.0000 - 2.0000i   0.0000 + 9.0000i   4.0000 + 7.0000i

La transposée conjuguée complexe de z est :

z'

ans =

   1.0000 - 2.0000i   6.0000 + 2.0000i
   7.0000 + 3.0000i   0.0000 - 9.0000i
   3.0000 - 4.0000i   4.0000 - 7.0000i

La transposée complexe non conjuguée, où la partie complexe de chaque élément conserve son signe, est indiquée par z.' :

z.'

ans =

   1.0000 + 2.0000i   6.0000 - 2.0000i
   7.0000 - 3.0000i   0.0000 + 9.0000i
   3.0000 + 4.0000i   4.0000 + 7.0000i

Pour les vecteurs complexes, les deux produits scalaires x'*y et y'*x sont des conjugués complexes l’un de l’autre, et le produit scalaire x'*x d’un vecteur complexe avec lui-même est un réel.

Multiplication de matrices

La multiplication de matrices est définie d’une manière qui reflète la composition des transformations linéaires sous-jacentes et permet la représentation compacte de systèmes d’équations linéaires simultanées. Le produit matriciel C = AB est défini lorsque la dimension des colonnes de A est égale à la dimension des lignes de B, ou lorsque l’une des deux est scalaire. Si A est de dimension m x p et si B est de dimension p x n, leur produit C est de dimension m x n. Le produit peut être défini à l’aide de boucles for, de la notation colon et de produits points vectoriels de MATLAB :

A = pascal(3);
B = magic(3);
m = 3; 
n = 3;
for i = 1:m
     for j = 1:n
        C(i,j) = A(i,:)*B(:,j);
     end
end

MATLAB utilise un astérisque pour indiquer la multiplication de matrices, comme dans C = A*B. La multiplication de matrices n’est pas commutative, autrement dit, A*B n’est généralement pas égal à B*A :

X = A*B

X =
      15    15    15
      26    38    26
      41    70    39

Y = B*A

Y =
      15    28    47
      15    34    60
      15    28    43

Une matrice peut être multipliée à droite par un vecteur colonne et à gauche par un vecteur ligne :

u = [3; 1; 4];
x = A*u

x =

     8
    17
    30

v = [2 0 -1];
y = v*B

y =

    12    -7    10

Les multiplications de matrices rectangulaires doivent satisfaire aux conditions de compatibilité des dimensions. Étant donné que A est de dimension 3 x 3 et C de dimension 3 x 2, vous pouvez les multiplier pour obtenir un résultat de dimension 3 x 2 (la dimension intérieure commune s’annule) :

X = A*C

X =

    24    17
    47    42
    79    77

Cependant, la multiplication ne fonctionne pas dans l’ordre inverse :

Y = C*A

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns 
in the first matrix matches the number of rows in the second matrix. To perform 
elementwise multiplication, use '.*'.

Vous pouvez multiplier n’importe quoi par un scalaire :

s = 10;
w = s*y

w =

   120   -70   100

Lorsque vous multipliez un tableau par un scalaire, le scalaire s’étend implicitement pour avoir la même taille que l’autre entrée. On appelle souvent cela l’expansion scalaire.

Matrice identité

La notation mathématique généralement acceptée utilise la lettre majuscule I pour indiquer les matrices identité. Ce sont des matrices de taille différente avec des 1 sur la diagonale principale et des 0 ailleurs. Ces matrices ont pour propriété que AI = A et IA = A lorsque les dimensions sont compatibles.

La version d’origine de MATLAB ne pouvait pas utiliser I à cette fin, car elle ne faisait pas la distinction entre les lettres majuscules et minuscules, et i était déjà utilisé pour les indices et les unités complexe. Un jeu de mots en anglais a donc été introduit. La fonction

eye(m,n)

renvoie une matrice identité rectangulaire de dimension m x n et eye(n) renvoie une matrice identité carrée n x n.

Inverse de matrice

Si une matrice A est carrée et non singulière (déterminant non nul), alors les équations AX = I et XA = I ont la même solution X. La solution est appelée l’inverse de A et est notée A^-1. La fonction inv et l’expression A^-1 calculent toutes deux l’inverse de la matrice.

A = pascal(3)

A =
       1     1     1
       1     2     3
       1     3     6

X = inv(A)

X =

    3.0000   -3.0000    1.0000
   -3.0000    5.0000   -2.0000
    1.0000   -2.0000    1.0000

A*X

ans =

    1.0000         0         0
    0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000

Le déterminant calculé par det est une mesure du facteur d’échelle de la transformation linéaire décrite par la matrice. Lorsque le déterminant est exactement égal à zéro, la matrice est singulière et aucun inverse n’existe.

d = det(A)

d =

     1

Certaines matrices sont presque singulières et, malgré le fait qu’une matrice inverse existe, le calcul s’expose à des erreurs numériques. La fonction cond calcule la mesure de conditionnement pour l’inversion, qui donne une indication de la validité des résultats de l’inversion de la matrice. La mesure de conditionnement est comprise entre 1 pour une matrice numériquement stable et Inf pour une matrice singulière.

c = cond(A)

c =

   61.9839

Il est rarement nécessaire de former l’inverse explicite d’une matrice. Une mauvaise utilisation de inv se produit fréquemment lors de la résolution d’un système d’équations linéaires Ax = b. La meilleure manière de résoudre cette équation, d'un point de vue temps d’exécution et exactitude numérique, est d’utiliser l’opérateur matriciel antislash x = A\b. Consultez mldivide pour plus d’informations.

Produit tensoriel de Kronecker

Le produit de Kronecker, kron(X,Y), de deux matrices est la matrice plus grande formée à partir de tous les produits possibles des éléments de X avec ceux de Y. Si X est de dimension m x n et si Y est de dimension p x q, alors kron(X,Y) est de dimension mp x nq. Ces éléments sont arrangés de manière à ce que chaque élément de X soit multiplié par la matrice entière Y :

[X(1,1)*Y  X(1,2)*Y  . . .  X(1,n)*Y
                     . . .
 X(m,1)*Y  X(m,2)*Y  . . .  X(m,n)*Y]

Le produit de Kronecker est souvent utilisé avec des matrices de 0 et de 1 pour constituer des copies répétées de petites matrices. Par exemple, si X est une matrice de dimension 2 x 2

X = [1   2
     3   4]

et I = eye(2,2) est la matrice identité de dimension 2 x 2, alors :

kron(X,I)

ans =

     1     0     2     0
     0     1     0     2
     3     0     4     0
     0     3     0     4

et

kron(I,X)

ans =

     1     2     0     0
     3     4     0     0
     0     0     1     2
     0     0     3     4

Outre kron, d’autres fonctions utiles pour répliquer des matrices sont repmat, repelem et blkdiag.

Normes de vecteurs et de matrices

La norme p d’un vecteur x,

$${\Vert x\Vert }_p={\left({\displaystyle \sum {\left|x_i\right|}^p}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$p$}\right.},$$

est calculée par norm(x,p). Cette opération est définie pour toute valeur de p > 1, mais les valeurs les plus courantes de p sont 1, 2 et ∞. La valeur par défaut est p = 2, ce qui correspond à la distance euclidienne ou à l’amplitude du vecteur :

v = [2 0 -1];
[norm(v,1) norm(v) norm(v,inf)]

ans =

    3.0000    2.2361    2.0000

La norme p d’une matrice A,

$${\Vert A\Vert }_p=\underset{x}{\max }\frac{{\Vert Ax\Vert }_p}{{\Vert x\Vert }_p},$$

peut être calculée pour p = 1, 2 et ∞ par norm(A,p). Là encore, la valeur par défaut est p = 2 :

A = pascal(3);
[norm(A,1) norm(A) norm(A,inf)]

ans =

   10.0000    7.8730   10.0000

Dans les cas où vous souhaitez calculer la norme de chaque ligne ou colonne d’une matrice, vous pouvez utiliser vecnorm :

vecnorm(A)

ans =

    1.7321    3.7417    6.7823

Utilisation d’un calcul multithread avec des fonctions d’algèbre linéaire

MATLAB supporte le calcul multithread pour plusieurs fonctions d’algèbre linéaire et fonctions numériques calculées éléments par éléments. Ces fonctions s’exécutent automatiquement sur plusieurs threads. Pour qu’une fonction ou une expression s’exécute plus rapidement sur plusieurs processeurs, plusieurs conditions doivent être remplies :

1. La fonction effectue des opérations qui peuvent facilement être partitionnées en sections qui s’exécutent simultanément. Ces sections doivent pouvoir s’exécuter avec une communication limitée entre les processus. Elles doivent nécessiter peu d’opérations séquentielles.

2. La taille des données est suffisamment importante pour que les avantages d’une exécution simultanée fassent plus que compenser le temps requis pour partitionner les données et gérer des threads d’exécution séparés. Par exemple, la plupart des fonctions s’accélèrent uniquement lorsque la matrice contient au minimum plusieurs milliers d’éléments.

3. L’opération n’est pas liée à la mémoire ; la durée de traitement ne dépend pas de la durée d’accès à la mémoire. En règle générale, les fonctions plus complexes sont davantage accélérées que les fonctions simples.

Les opérateurs de multiplication (X*Y) et de puissance (X^p) de matrices démontrent une vitesse d'exécution significativement accrue sur des matrices de grandes tailles en double précision (de l’ordre de 10 000 éléments). Les fonctions d’analyse de matrices det, rcond, hess et expm démontrent également une vitesse sensiblement accrue sur des matrices de grandes tailles en double précision.

Sujets associés

• Systèmes d’équations linéaires
• Factorisations
• Valeurs propres
• Valeurs singulières

Sites web externes

• Méthodes matricielles de l’algèbre linéaire (Ressources pédagogiques MathWorks)

Systèmes d’équations linéaires

• Considérations liées au calcul

• Solution générale

• Systèmes carrés

• Systèmes surdéterminés

• Systèmes sous-déterminés

• Résolution pour plusieurs côtés droits

• Méthodes itératives

• Calcul multithread

Considérations liées au calcul

L’un des problèmes les plus importants en calcul scientifique est la résolution de systèmes d’équations linéaires simultanées.

En notation matricielle, le problème général prend la forme suivante : soit deux matrices A et b, existe-t-il une matrice unique x telle que Ax = b ou xA = b ?

Il est instructif d’étudier un exemple de dimension 1 x 1. Par exemple, l’équation

7x = 21

a-t-elle une solution unique ?

La réponse, bien sûr, est oui. L’équation a pour solution unique x = 3. La solution est facilement obtenue à l’aide de la division :

x = 21/7 = 3.

La solution n’est pas obtenue habituellement en calculant l’inverse de 7, soit 7^-1 = 0,142857…, puis en multipliant 7^-1 par 21. Cela représenterait plus de travail et, si 7^-1 est représenté avec un nombre fini de chiffres, cela réduirait la précision. Des considérations similaires s’appliquent aux ensembles d’équations linéaires avec plus d’une inconnue, MATLAB résout ces équations sans calculer l’inverse de la matrice.

Bien que cela ne constitue pas la notation mathématique standard, MATLAB utilise la terminologie de division habituelle dans le cas d’un scalaire pour décrire la résolution d’un système général d’équations simultanées. Les deux symboles de la division, le slash, /, et l'antislash, \, correspondent aux deux fonctions MATLAB mrdivide et mldivide. Ces opérateurs sont utilisés pour les deux situations où la matrice inconnue apparaît à gauche ou à droite de la matrice de coefficients :

x = b/A	Indique la solution de l’équation matricielle xA = b, obtenue en utilisant mrdivide.
x = A\b	Indique la solution de l’équation matricielle Ax = b, obtenue en utilisant mldivide.

Considérez que l’on « divise » les deux côtés de l’équation Ax = b ou xA = b par A. La matrice de coefficients A est toujours le « dénominateur ».

Les conditions de compatibilité de dimension pour x = A\b nécessitent que les deux matrices A et b comportent le même nombre de lignes. La solution x comporte alors le même nombre de colonnes que b, et son nombre de lignes est égale au nombre de colonnes de A. Pour x = b/A, les rôles des lignes et des colonnes sont échangés.

En pratique, les équations linéaires de la forme Ax = b sont plus fréquentes que celles de la forme xA = b. Par conséquent, l'antislash est utilisé beaucoup plus fréquemment que le slash. Le reste de cette section se concentre sur l’opérateur antislash ; les propriétés correspondantes de l’opérateur de slash peuvent être déduites de l'égalité :

(b/A)' = (A'\b').

La matrice de coefficients A ne doit pas être carrée. Si A a une taille de m x n, trois cas sont possibles :

m = n	Système carré. Cherchez une solution exacte.
m > n	Système surdéterminé, avec plus d’équations que d’inconnues. Trouvez une solution au sens des moindres carrés.
m < n	Système sous-déterminé, avec moins d’équations que d’inconnues. Trouvez une solution basique avec au maximum m composants non nuls.

L’algorithme mldivide.  L’opérateur mldivide utilise des solveurs différents pour traiter différents types de matrices de coefficients. Les différents cas sont déterminés automatiquement en examinant la matrice de coefficients. Pour plus d’informations, consultez la section « Algorithmes » de la page de référence mldivide.

Solution générale

La solution générale d’un système d’équations linéaires Ax = b décrit toutes les solutions possibles. Vous pouvez trouver la solution générale :

1. En résolvant le système homogène correspondant Ax = 0. Pour cela, utilisez la commande null en tapant null(A). Celle-ci renvoie une base pour l’espace de solution de Ax = 0. Toute solution est une combinaison linéaire de vecteurs de base.

2. En trouvant une solution particulière au système non homogène Ax = b.

Vous pouvez alors écrire toute solution à Ax = b comme la somme de la solution particulière à Ax = b de l’étape 2, plus une combinaison linéaire des vecteurs de base de l’étape 1.

Le reste de cette section décrit comment utiliser MATLAB pour trouver une solution particulière à Ax = b, comme exposé dans l’étape 2.

Systèmes carrés

La situation la plus fréquente implique une matrice carrée de coefficients A et un vecteur colonne unique du côté droit de l'équation b.

Matrice non singulière de coefficients.  Si la matrice A est non singulière, alors la solution x = A\b est de la même taille que b. Par exemple :

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

Il est possible de confirmer que A*x est exactement égal à u.

Si A et b sont carrées et de la même taille, x= A\b est également de cette taille :

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

Il est possible de confirmer que A*x est exactement égal à b.

Ces deux exemples ont des solutions exactes et entières. Cela est dû au fait que la matrice de coefficients choisie est pascal(3), une matrice de rang plein (non singulière).

Matrice de coefficients singulière.  Une matrice carrée A est singulière si elle ne comporte pas de colonnes linéairement indépendantes. Si A est singulière, la solution à Ax = b soit n’existe pas, soit n’est pas unique. L’opérateur antislash, A\b, affiche un avertissement si A est presque singulière ou s’il détecte une singularité exacte.

Si A est singulière et si Ax = b a une solution, vous pouvez trouver une solution particulière qui n’est pas unique, en tapant

P = pinv(A)*b

pinv(A) est un pseudo-inverse de A. Si Ax = b n’a pas de solution exacte, alors pinv(A) renvoie une solution des moindres carrés.

Par exemple :

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

est singulière, comme vous pouvez le vérifier en tapant

rank(A)

ans =

     2

Étant donné que A n’est pas de rang plein, elle comporte certaines valeurs singulières égales à zéro.

Solutions exactes. Pour b =[5;2;12], l’équation Ax = b a une solution exacte, donnée par

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

Vérifiez que pinv(A)*b est une solution exacte en tapant

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

Solutions au sens des moindres carrés. Toutefois, si b = [3;6;0], Ax = b n’a pas de solution exacte. Dans ce cas, pinv(A)*b renvoie une solution au sens des moindres carrés. Si vous tapez

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

vous n’obtenez pas le vecteur d’origine b.

Vous pouvez déterminer si Ax = b a une solution exacte en trouvant la forme échelonnée réduite de la matrice augmentée [A b]. Pour ce faire, dans cet exemple, tapez

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Étant donné que la ligne du bas contient uniquement des zéros à l’exception de la dernière entrée, l’équation n’a pas de solution. Dans ce cas, pinv(A) renvoie une solution au sens des moindres carrés.

Systèmes surdéterminés

Cet exemple montre comment des systèmes surdéterminés sont souvent rencontrés dans divers types de courbes correspondant à des données expérimentales.

Une quantité y est mesurée à différentes valeurs de temps t pour produire les observations suivantes. Vous pouvez saisir les données et les afficher dans une table avec les instructions suivantes.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)

B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Essayez de modéliser les données avec une fonction exponentielle décroissante

$y(t)=c_1+c_2{e}^{-t}$.

L’équation précédente indique que le vecteur y doit être estimé à l’aide d’une combinaison linéaire de deux autres vecteurs. L’un est un vecteur constant contenant uniquement des 1 et l’autre est le vecteur avec les composants exp(-t). Les coefficients inconnus, $$c_1$$ et $$c_2$$, peuvent être calculés en effectuant un ajustement par la méthode des moindres carrés, qui minimise la somme des carrés des écarts entre les données et le modèle. Il y a six équations à deux inconnues, représentées par une matrice de 6 x 2.

E = [ones(size(t)) exp(-t)]

E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

Utilisez l’opérateur antislash pour obtenir la solution au sens des moindres carrés.

c = E\y

c = 2×1

    0.4760
    0.3413

En d’autres termes, l’ajustement par la méthode des moindres carrés aux données est

$$y(t)=0.4760+0.3413e^{-t}.$$

Les instructions suivantes évaluent le modèle à des incréments régulièrement espacés dans t, puis établissent un graphique du résultat avec les données d’origine :

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*c n’est pas exactement égal à y, mais la différence peut être tout à fait inférieure aux erreurs de mesure dans les données d’origine.

Une matrice rectangulaire A présente une déficience de rang si elle ne comporte pas des colonnes linéairement indépendantes. Si A présente une déficience de rang, alors la solution au sens des moindres carrés à AX = B n’est pas unique. A\B affiche un avertissement si A présente une déficience de rang et produit une solution au sens des moindres carrés. Vous pouvez utiliser lsqminnorm pour trouver la solution X qui comporte la norme minimum parmi toutes les solutions.

Systèmes sous-déterminés

Cet exemple montre comment la solution à des systèmes sous-déterminés n’est pas unique. Les systèmes linéaires sous-déterminés comportent plus d’inconnues que d’équations. L’opération matricielle de division à gauche dans MATLAB trouve une solution basique au sens des moindres carrés, qui comporte au maximum m composants non nuls pour une matrice de coefficients de dimension m par n.

Voici un petit exemple au hasard :

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)

R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

Le système linéaire Rp = b comprend deux équations à quatre inconnues. Étant donné que la matrice de coefficients contient de petits entiers, il est approprié d’utiliser la commande format pour afficher la solution au format rationnel. La solution particulière est obtenue avec

format rat
p = R\b

p =

       0       
      17/7     
       0       
     -29/7    

L’un des composants non nuls est p(2), car R(:,2) est la colonne de R avec la norme la plus grande. L’autre composant non nul est p(4), parce que R(:,4) domine après l’élimination de R(:,2).

La solution générale complète au système sous-déterminé peut être caractérisée en ajoutant p à une combinaison linéaire arbitraire des vecteurs spatiaux nuls, qui peuvent être trouvés à l’aide de la fonction null avec une option demandant une base rationnelle.

Z = null(R,'r')

Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

Il est possible de confirmer que R*Z vaut zéro et que le résidu R*x - b est petit pour tout vecteur x, où

x = p + Z*q

Étant donné que les colonnes de Z sont des vecteurs spatiaux nuls, le produit Z*q est une combinaison linéaire de ces vecteurs :

$$Zq=\left(\begin{array}{cc}{\stackrel{\rightharpoonup }{x}}_1& {\stackrel{\rightharpoonup }{x}}_2\end{array}\right)\left(\begin{array}{c}u\\ w\end{array}\right)=u{\stackrel{\rightharpoonup }{x}}_1+w{\stackrel{\rightharpoonup }{x}}_2.$$

À titre d’illustration, choisissez un q arbitraire et construisez x.

q = [-2; 1];
x = p + Z*q;

Calculez la norme du résidu.

format short
norm(R*x - b)

ans =

   2.6645e-15

Lorsqu’une infinité de solutions sont possibles, la solution avec la norme minimum présente un intérêt particulier. Vous pouvez utiliser lsqminnorm pour calculer la solution de norme minimum au sens des moindres carrés. Cette solution présente la plus petite valeur possible pour norm(p).

p = lsqminnorm(R,b)

p =

    -207/137   
     365/137   
      79/137   
    -424/137  

Résolution pour plusieurs côtés droits

Certains problèmes portent sur la résolution de systèmes linéaires qui comportent la même matrice de coefficients A, mais différents côtés droits b. Lorsque les différentes valeurs de b sont disponibles simultanément, vous pouvez construire b comme une matrice avec plusieurs colonnes et résoudre tous les systèmes d’équations simultanément à l’aide d’une seule commande de antislash : X = A\[b1 b2 b3 …].

Toutefois, il arrive que différentes valeurs de b ne soient pas toutes disponibles en même temps, ce qui signifie que vous devez résoudre plusieurs systèmes d’équations consécutivement. Lorsque vous résolvez l’un de ces systèmes d’équations à l’aide du slash (/) ou de l'antislash (\), l’opérateur factorise la matrice de coefficients A et utilise cette décomposition de la matrice pour calculer la solution. Cependant, à chaque fois que vous résoudrez un système d’équations similaire avec un b différent, l’opérateur calculera la même décomposition de A, ce qui est un calcul redondant.

La solution à ce problème est de précalculer la décomposition de A, puis de réutiliser les facteurs à résoudre pour les différentes valeurs de b. Cependant, en pratique, précalculer la décomposition de cette manière peut être difficile, car vous devez savoir quelle décomposition calculer (LU, LDL, Cholesky, etc.) et comment multiplier les facteurs pour résoudre le problème. Par exemple, avec la décomposition LU, vous devez résoudre deux systèmes linéaires pour résoudre le système d’origine Ax = b :

[L,U] = lu(A);
x = U \ (L \ b);

Pour résoudre des systèmes linéaires avec plusieurs côtés droit consécutifs, il est plutôt recommandé d’utiliser des objets decomposition. Ces objets vous permettent de tirer profit de l'avantage en terme de performance du précalcul de la décomposition matricielle, mais ils ne nécessitent pas de savoir comment utiliser les facteurs matriciels. Vous pouvez remplacer la décomposition LU précédente par :

dA = decomposition(A,'lu');
x = dA\b;

Si vous n’êtes pas certain de la décomposition à utiliser, decomposition(A) choisit le type correct en fonction des propriétés de A, comme le fait l'antislash.

Voici un test simple des avantages possibles de cette approche en terme de performance. Ce test résout le même système linéaire creux 100 fois avec l'antislash (\) et la fonction decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc

Elapsed time is 9.006156 seconds.

% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc

Elapsed time is 0.374347 seconds.

Pour ce problème, la solution de decomposition est beaucoup plus rapide que l’utilisation de l'antislash seul, mais la syntaxe reste simple.

Méthodes itératives

Si la matrice de coefficients A est grande et creuse, les méthodes de factorisation ne sont généralement pas efficaces. Les méthodes itératives génèrent une série de solutions approximatives. MATLAB propose plusieurs méthodes itératives pour traiter de grandes matrices creuses.

Fonction	Description
pcg	Méthode du gradient conjugué avec préconditionnement. Cette méthode est appropriée pour la matrice hermitienne définie positive A.
bicg	Méthode du gradient biconjugué
bicgstab	Méthode stabilisée du gradient biconjugué
bicgstabl	Méthode BiCGStab(l)
cgs	Méthode carrée du gradient conjugué
gmres	Généralisation de la méthode de minimisation du résidu
lsqr	Méthode LSQR
minres	Méthode de minimisation du résidu. Cette méthode est appropriée pour la matrice hermitienne A.
qmr	Méthode de quasi-minimisation du résidu
symmlq	Méthode LQ symétrique
tfqmr	Méthode QMR sans transposition

Calcul multithread

MATLAB supporte le calcul multithread pour plusieurs fonctions d’algèbre linéaire et des fonctions numériques calculées éléments par éléments. Ces fonctions s’exécutent automatiquement sur plusieurs threads. Pour qu’une fonction ou une expression s’exécute plus rapidement sur plusieurs processeurs, plusieurs conditions doivent être remplies :

1. La fonction effectue des opérations qui peuvent facilement être partitionnées en sections qui s’exécutent simultanément. Ces sections doivent pouvoir s’exécuter avec une communication limitée entre les processus. Elles doivent nécessiter peu d’opérations séquentielles.

2. La taille des données est suffisamment importante pour que les avantages d’une exécution simultanée fassent plus que compenser le temps requis pour partitionner les données et gérer des threads d’exécution séparés. Par exemple, la plupart des fonctions s’accélèrent uniquement lorsque la matrice contient au minimum plusieurs milliers d’éléments.

3. L’opération n’est pas liée à la mémoire ; la durée de traitement ne dépend pas de la durée d’accès à la mémoire. En règle générale, les fonctions plus complexes sont davantage accélérées que les fonctions simples.

inv, lscov, linsolve et mldivide démontrent une vitesse d'exécution significativement accrue sur des matrices de grandes tailles en double précision (de l’ordre de 10 000 éléments ou plus) lorsque le multithreading est activé.

Voir aussi

mldivide | mrdivide | pinv | decomposition | lsqminnorm

Sujets associés

• Matrices dans l’environnement MATLAB
• Factorisations
• Valeurs propres
• Valeurs singulières

Sites web externes

• Méthodes matricielles de l’algèbre linéaire (Ressources pédagogiques MathWorks)

Factorisations

• Introduction

• Factorisation de Cholesky

• Factorisation LU

• Factorisation QR

• Utilisation du calcul multithread pour la factorisation

Introduction

Les trois factorisations matricielles abordées dans cette section utilisent des matrices triangulaires, où tous les éléments au-dessus ou au-dessous de la diagonale sont égaux à zéro. Les systèmes d’équations linéaires impliquant des matrices triangulaires sont faciles et rapides à résoudre à l’aide d’une substitution avant ou arrière.

Factorisation de Cholesky

La factorisation de Cholesky exprime une matrice symétrique comme le produit d’une matrice triangulaire et de sa transposée

A = R′R,

où R est une matrice triangulaire supérieure.

Toutes les matrices symétriques ne peuvent pas être factorisées de cette manière ; les matrices qui présentent une factorisation de ce type sont dites définies positives. Cela implique que tous les éléments de la diagonale de A sont positifs et que les éléments hors diagonale ne sont « pas trop grands ». Les matrices de Pascal offrent un exemple intéressant. Tout au long de ce chapitre, la matrice A utilisée comme exemple était la matrice de Pascal de dimension 3 x 3. Passez temporairement à la matrice de dimension 6 par 6 :

A = pascal(6)

A =
       1     1     1     1     1     1
       1     2     3     4     5     6
       1     3     6    10    15    21
       1     4    10    20    35    56
       1     5    15    35    70   126
       1     6    21    56   126   252

Les éléments de A sont des coefficients binomiaux. Chaque élément est la somme de ses voisins au nord et à l’ouest. La factorisation de Cholesky est

R = chol(A)

R =
     1     1     1     1     1     1
     0     1     2     3     4     5
     0     0     1     3     6    10
     0     0     0     1     4    10
     0     0     0     0     1     5
     0     0     0     0     0     1

Les éléments sont à nouveau des coefficients binomiaux. Le fait que R'*R est égal à A démontre une identité impliquant les sommes de produits des coefficients binomiaux.

Remarque

La factorisation de Cholesky s’applique également aux matrices complexes. Toute matrice complexe ayant une factorisation de Cholesky satisfait

A′ = A

et est dite hermitienne définie positive.

La factorisation de Cholesky permet de remplacer le système linéaire

Ax = b

par

R′Rx = b.

Étant donné que l’opérateur antislash reconnaît les systèmes triangulaires, ceci peut être résolu rapidement dans l’environnement MATLAB avec

x = R\(R'\b)

Si A a une taille de n x n, la complexité algorithmique de chol(A) est O(n³), mais la complexité de la solution précédente utilisant l'operateur antislash est seulement de O(n²).

Factorisation LU

La factorisation LU, ou élimination de Gauss-Jordan, exprime toute matrice carrée A comme le produit d’une permutation d’une matrice triangulaire inférieure et d’une matrice triangulaire supérieure

A = LU,

où L est une permutation d’une matrice triangulaire inférieure avec des 1 sur sa diagonale et U est une matrice triangulaire supérieure.

Les permutations sont nécessaires pour des raisons à la fois théoriques et liées au calcul. La matrice

$$\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$$

ne peut pas être exprimée comme le produit de matrices triangulaires sans échanger ses deux lignes. Bien que la matrice

$$\left[\begin{array}{cc}\epsilon & 1\\ 1& 0\end{array}\right]$$

puisse être exprimée comme le produit de matrices triangulaires, lorsque ε est petite, les éléments dans les facteurs sont grands et amplifient les erreurs, c’est pourquoi même si les permutations ne sont pas strictement nécessaires, elles sont souhaitables. Le pivotement partiel garantit que les éléments de L sont bornés à 1 en amplitude et que les éléments de U ne sont pas beaucoup plus grands que ceux de A.

Par exemple :

[L,U] = lu(B)

L =
    1.0000         0         0
    0.3750    0.5441    1.0000
    0.5000    1.0000         0

U =
    8.0000    1.0000    6.0000
         0    8.5000   -1.0000
         0         0    5.2941

La factorisation LU de A permet de résoudre rapidement le système linéaire

A*x = b

avec

x = U\(L\b)

Les déterminants et les inverses sont calculés à partir de la factorisation LU avec

det(A) = det(L)*det(U)

et

inv(A) = inv(U)*inv(L)

Vous pouvez également calculer les déterminants à l’aide de det(A) = prod(diag(U)), bien que les signes des déterminants puissent être inversés.

Factorisation QR

Une matrice orthogonale, ou une matrice avec des colonnes orthonormales, est une matrice réelle dont les colonnes sont perpendiculaires entre elles et de norme 1. Si Q est orthogonale, alors

Q^TQ = I,

où I est la matrice identité.

Les matrices orthogonales les plus simples sont des rotations de coordonnées à deux dimensions :

$$\left[\begin{array}{cc}\cos (\theta )& \sin (\theta )\\ -\sin (\theta )& \cos (\theta )\end{array}\right].$$

Pour les matrices complexes, le terme correspondant est unitaire. Les matrices orthogonales et unitaires sont intéressantes pour le calcul numérique, car elles conservent les longueurs ainsi que les angles et n’amplifient pas les erreurs.

La factorisation orthogonale, ou QR, exprime toute matrice rectangulaire comme le produit d’une matrice orthogonale ou unitaire et d’une matrice triangulaire supérieure. Elle peut également impliquer une permutation de colonnes :

A = QR

ou

AP = QR,

où Q est orthogonale ou unitaire, R est triangulaire supérieure et P est une permutation.

Il existe quatre variantes de factorisation QR – complète ou réduite, et avec ou sans permutation de colonnes.

Les systèmes linéaires surdéterminés impliquent une matrice rectangulaire avec plus de lignes que de colonnes, c’est-à-dire de dimension m x n avec m > n. La factorisation QR complete produit une matrice orthogonale carrée Q de dimension m x m et une matrice triangulaire supérieure rectangulaire R de dimension m x n :

C=gallery('uniformdata',[5 4], 0);
[Q,R] = qr(C)

Q =

    0.6191    0.1406   -0.1899   -0.5058    0.5522
    0.1506    0.4084    0.5034    0.5974    0.4475
    0.3954   -0.5564    0.6869   -0.1478   -0.2008
    0.3167    0.6676    0.1351   -0.1729   -0.6370
    0.5808   -0.2410   -0.4695    0.5792   -0.2207


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648
         0         0         0         0

Souvent, les m – n dernières colonnes de Q ne sont pas nécessaires, car elles sont multipliées par les 0 de la portion basse de R. La factorisation QR réduite produit donc une matrice rectangulaire Q de dimensionm x n avec des colonnes orthonormales et une matrice triangulaire supérieure carrée R de dimension n x n. Pour l’exemple de dimension 5 x 4, le gain est limité, mais pour des matrices plus grandes et très rectangulaires, les gains en temps et en mémoire peuvent être assez significatifs :

[Q,R] = qr(C,0)	
Q =

    0.6191    0.1406   -0.1899   -0.5058
    0.1506    0.4084    0.5034    0.5974
    0.3954   -0.5564    0.6869   -0.1478
    0.3167    0.6676    0.1351   -0.1729
    0.5808   -0.2410   -0.4695    0.5792


R =

    1.5346    1.0663    1.2010    1.4036
         0    0.7245    0.3474   -0.0126
         0         0    0.9320    0.6596
         0         0         0    0.6648

Contrairement à la factorisation LU, la factorisation QR ne nécessite pas de pivotement ou de permutations. Mais une permutation facultative de colonne, déclenchée par la présence d’un troisième argument de sortie, est utile pour détecter une singularité ou une déficience de rang. À chaque étape de la factorisation, la colonne de la matrice restante non factorisée, dont la norme est la plus grande, est utilisée comme base pour cette étape. Ceci permet de garantir que les éléments de la diagonale de R apparaissent par ordre décroissant et que toute dépendance linéaire au niveau des colonnes est presque certainement révélée par l’examen de ces éléments. Pour le petit exemple donné ici, la deuxième colonne de C a une norme plus grande que la première, les deux colonnes sont donc échangées :

[Q,R,P] = qr(C)

Q =
   -0.3522    0.8398   -0.4131
   -0.7044   -0.5285   -0.4739
   -0.6163    0.1241    0.7777

R =
  -11.3578   -8.2762
         0    7.2460
         0         0

P =
     0     1
     1     0

Lorsque les méthodes de factorisation réduite et de permutations de colonnes sont combinées, le troisième argument de sortie est un vecteur de permutation, plutôt qu’une matrice de permutation :

[Q,R,p] = qr(C,0)

Q =
   -0.3522    0.8398
   -0.7044   -0.5285
   -0.6163    0.1241

R =
  -11.3578   -8.2762
         0    7.2460


p =
     2     1

La factorisation QR transforme un système linéaire surdéterminé en un système triangulaire équivalent. L’expression

norm(A*x - b)

équivaut à

norm(Q*R*x - b)

La multiplication par les matrices orthogonales préserve la norme euclidienne, ce qui fait que cette expression est aussi égale à

norm(R*x - y)

où y = Q'*b. Étant donné que les dernières m-n lignes de R sont nulles, cette expression se scinde en deux parties :

norm(R(1:n,1:n)*x - y(1:n))

et

norm(y(n+1:m))

Lorsque A est de rang plein, il est possible de la résoudre pour x de manière à ce que la première de ces expressions soit nulle. La seconde expression donne alors la norme du résidu. Lorsque A n’est pas de rang plein, la structure triangulaire de R offre la possibilité de trouver une solution basique au problème des moindres carrés.

Utilisation du calcul multithread pour la factorisation

Le logiciel MATLAB supporte le calcul multithread pour plusieurs fonctions d’algèbre linéaire et des fonctions numériques calculées éléments par éléments. Ces fonctions s’exécutent automatiquement sur plusieurs threads. Pour qu’une fonction ou une expression s’exécute plus rapidement sur plusieurs processeurs, plusieurs conditions doivent être remplies :

1. La fonction effectue des opérations qui peuvent facilement être partitionnées en sections qui s’exécutent simultanément. Ces sections doivent pouvoir s’exécuter avec une communication limitée entre les processus. Elles doivent nécessiter peu d’opérations séquentielles.

2. La taille des données est suffisamment importante pour que les avantages d’une exécution simultanée fassent plus que compenser le temps requis pour partitionner les données et gérer des threads d’exécution séparés. Par exemple, la plupart des fonctions s’accélèrent uniquement lorsque la matrice contient au minimum plusieurs milliers d’éléments.

3. L’opération n’est pas liée à la mémoire ; la durée de traitement ne dépend pas de la durée d’accès à la mémoire. En règle générale, les fonctions plus complexes sont davantage accélérées que les fonctions simples.

lu et qr démontrent une vitesse d'exécution significativement accrue sur des matrices de grandes tailles en double précision (de l’ordre de 10 000 éléments).

Voir aussi

lu | qr | chol

Sujets associés

• Matrices dans l’environnement MATLAB

Puissances et exponentielles

Ce thème montre comment calculer des puissances et des exponentielles de matrices à l’aide de diverses méthodes.

Puissances d’entiers positifs

Si A est une matrice carrée et p un entier positif, alors A^p multiplie A par elle-même p-1 fois. Par exemple :

A = [1 1 1
     1 2 3
     1 3 6];
A^2

ans = 3×3

     3     6    10
     6    14    25
    10    25    46

Puissances inverses et fractionnelles

Si A est carrée et non singulière, alors A^(-p) multiplie inv(A) par elle-même p-1 fois.

A^(-3)

ans = 3×3

  145.0000 -207.0000   81.0000
 -207.0000  298.0000 -117.0000
   81.0000 -117.0000   46.0000

MATLAB® calcule inv(A) et A^(-1) avec le même algorithme, de sorte que les résultats sont exactement les mêmes. inv(A) et A^(-1) génèrent toutes deux des avertissements si la matrice est proche d’être singulière.

isequal(inv(A),A^(-1))

ans = logical
   1

Les puissances fractionnelles, telles que A^(2/3), sont également permises. Les résultats obtenus avec des puissances fractionnelles dépendent de la distribution des valeurs propres de la matrice.

A^(2/3)

ans = 3×3

    0.8901    0.5882    0.3684
    0.5882    1.2035    1.3799
    0.3684    1.3799    3.1167

Puissances élément par élément

L’opérateur .^ calcule les puissances élément par élément. Par exemple, pour élever au carré chaque élément d’une matrice, vous pouvez utiliser A.^2.

A.^2

ans = 3×3

     1     1     1
     1     4     9
     1     9    36

Racines carrées

La fonction sqrt est une manière pratique de calculer la racine carrée de chaque élément d’une matrice. Une autre manière consiste à utiliser A.^(1/2).

sqrt(A)

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    1.4142    1.7321
    1.0000    1.7321    2.4495

Pour d’autres racines, vous pouvez utiliser nthroot. Par exemple, calculez A.^(1/3).

nthroot(A,3)

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    1.2599    1.4422
    1.0000    1.4422    1.8171

Ces racines éléments par éléments diffèrent de la racine carrée matricielle, qui calcule une seconde matrice $\mathit{B}$ telle que $\mathit{A}=\mathrm{BB}$. La fonction sqrtm(A) calcule A^(1/2) à l’aide d’un algorithme plus précis. Le m de sqrtm distingue cette fonction de sqrt(A), qui, comme A.^(1/2), procède élément par élément.

B = sqrtm(A)

B = 3×3

    0.8775    0.4387    0.1937
    0.4387    1.0099    0.8874
    0.1937    0.8874    2.2749

B^2

ans = 3×3

    1.0000    1.0000    1.0000
    1.0000    2.0000    3.0000
    1.0000    3.0000    6.0000

Bases scalaires

En plus d’élever une matrice à une puissance, vous pouvez également élever un scalaire à la puissance d’une matrice.

2^A

ans = 3×3

   10.4630   21.6602   38.5862
   21.6602   53.2807   94.6010
   38.5862   94.6010  173.7734

Lorsque vous élevez un scalaire à la puissance d’une matrice, MATLAB utilise les valeurs propres et les vecteurs propres d’une matrice pour calculer la puissance de la matrice. Si [V,D] = eig(A), alors $2^{\mathit{A}}={\mathit{V}\mathrm{2}}^{\mathit{D}}{\mathit{V}}^{-1}$.

[V,D] = eig(A);
V*2^D*V^(-1)

ans = 3×3

   10.4630   21.6602   38.5862
   21.6602   53.2807   94.6010
   38.5862   94.6010  173.7734

Exponentielles de matrices

L’exponentielle de matrice est un cas particulier d’élévation d’un scalaire à une puissance de matrice. La base d’une exponentielle de matrice est le nombre d’Euler e = exp(1).

e = exp(1);
e^A

ans = 3×3
103 ×

    0.1008    0.2407    0.4368
    0.2407    0.5867    1.0654
    0.4368    1.0654    1.9418

La fonction expm est une manière plus pratique de calculer les exponentielles d’une matrice.

expm(A)

ans = 3×3
103 ×

    0.1008    0.2407    0.4368
    0.2407    0.5867    1.0654
    0.4368    1.0654    1.9418

L’exponentielle de la matrice peut être calculée de plusieurs manières. Consultez Matrix Exponentials pour plus d’informations.

Traitement des petits nombres

Les fonctions MATLAB log1p et expm1 calculent $\log \left(1+\mathit{x}\right)$ et ${\mathit{e}}^{\mathit{x}}-1$ de manière précise pour de très petites valeurs de $\mathit{x}$. Par exemple, si vous essayez d’additionner un nombre inférieur à la précision de la machine à 1, le résultat est arrondi à 1.

log(1+eps/2)

ans = 0

Cependant, log1p peut renvoyer une réponse plus précise.

log1p(eps/2)

ans = 1.1102e-16

De même pour ${\mathit{e}}^{\mathit{x}}-1$, si $\mathit{x}$ est très petit, il est arrondi à zéro.

exp(eps/2)-1

ans = 0

Là encore, expm1 est capable de renvoyer une réponse plus précise.

expm1(eps/2)

ans = 1.1102e-16

Voir aussi

exp | expm | expm1 | power | mpower | sqrt | sqrtm | nthroot

Sujets associés

• Matrices dans l’environnement MATLAB
• Matrix Exponentials

Valeurs propres

• Décomposition d'une matrice en éléments propres

• Valeurs propres multiples

• Décomposition de Schur

Décomposition d'une matrice en éléments propres

Une valeur propre et un vecteur propre d’une matrice carrée A sont, respectivement, un scalaire λ et un vecteur non nul υ qui satisfont à

Aυ = λυ.

Avec les valeurs propres sur la diagonale d’une matrice diagonale Λ et les vecteurs propres correspondants qui forment les colonnes d’une matrice V, vous avez

AV = VΛ.

Si V est non singulière, cette expression devient la décomposition en éléments propres

A = VΛV^–1.

La matrice de coefficients de l’équation différentielle dx/dt = Ax est un bon exemple :

A =
     0    -6    -1
     6     2   -16
    -5    20   -10

La solution à cette équation est exprimée en terme d’exponentielle de matrice x(t) = e^tAx(0). L’instruction

lambda = eig(A)

donne un vecteur colonne contenant les valeurs propres de A. Pour cette matrice, les valeurs propres sont complexes :

lambda =
     -3.0710         
     -2.4645+17.6008i
     -2.4645-17.6008i

La partie réelle de chacune des valeurs propres est négative, donc e^λt s’approche de zéro lorsque t augmente. La partie imaginaire non nulle de deux des valeurs propres, ±ω, fait contribuer le composant oscillatoire, sin(ωt), à la solution de l’équation différentielle.

Avec deux arguments de sortie, eig calcule les vecteurs propres et stocke les valeurs propres dans une matrice diagonale :

[V,D] = eig(A)

V =
  -0.8326         0.2003 - 0.1394i   0.2003 + 0.1394i
  -0.3553        -0.2110 - 0.6447i  -0.2110 + 0.6447i
  -0.4248        -0.6930            -0.6930          

D =
  -3.0710                 0                 0         
        0           -2.4645+17.6008i        0         
        0                 0           -2.4645-17.6008i

Le premier vecteur propre est réel et les deux autres vecteurs sont des conjugués complexes l’un de l’autre. Les trois vecteurs sont normalisés de manière à avoir une longueur euclidienne, norm(v,2), égale à 1.

La matrice V*D*inv(V), qui peut être écrite de manière plus brève V*D/V, est dans la marge d'erreur d'arrondi de A. Et inv(V)*A*V, ou V\A*V, est dans la marge d'erreur d'arrondi de D.

Valeurs propres multiples

Certaines matrices n’ont pas de décomposition en éléments propres. Ces matrices ne sont pas diagonalisables. Par exemple :

A = [ 1    -2    1 
      0     1    4 
      0     0    3 ]

Pour cette matrice

[V,D] = eig(A)

produit

V =

    1.0000    1.0000   -0.5571
         0    0.0000    0.7428
         0         0    0.3714


D =

     1     0     0
     0     1     0
     0     0     3

Il existe une double valeur propre à λ = 1. La première et la deuxième colonne de V sont identiques. Pour cette matrice, il n’existe pas d’ensemble complet de vecteurs propres linéairement indépendants.

Décomposition de Schur

De nombreux calculs matriciels avancés ne nécessitent pas de décompositions en éléments propres. Ils s’appuient plutôt sur la décomposition de Schur

A = USU ′ ,

où U est une matrice orthogonale et S est une matrice triangulaire supérieure par blocs avec des blocs de 1 x 1 et de 2 x 2 sur la diagonale. Les valeurs propres sont révélées par les éléments de la diagonale et les blocs de S, tandis que les colonnes de U fournissent une base orthogonale, qui a des propriétés numériques bien meilleures qu’un ensemble de vecteurs propres.

Par exemple, comparez les décompositions de Schur et en éléments propres de cette matrice défectueuse :

A = [ 6    12    19 
     -9   -20   -33 
      4     9    15 ];

[V,D] = eig(A)

V =

  -0.4741 + 0.0000i  -0.4082 - 0.0000i  -0.4082 + 0.0000i
   0.8127 + 0.0000i   0.8165 + 0.0000i   0.8165 + 0.0000i
  -0.3386 + 0.0000i  -0.4082 + 0.0000i  -0.4082 - 0.0000i


D =

  -1.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   1.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.0000 - 0.0000i

[U,S] = schur(A)

U =

   -0.4741    0.6648    0.5774
    0.8127    0.0782    0.5774
   -0.3386   -0.7430    0.5774


S =

   -1.0000   20.7846  -44.6948
         0    1.0000   -0.6096
         0    0.0000    1.0000

La matrice A est défectueuse, car elle ne comporte pas d’ensemble complet de vecteurs propres linéairement indépendants (les deuxième et troisième colonnes de V sont identiques). Étant donné que toutes les colonnes de V ne sont pas linéairement indépendantes, elle a un nombre de conditions élevé environ égal à 1e8. Cependant, schur peut calculer trois vecteurs de base différents dans U. Étant donné que U est orthogonale, cond(U) = 1.

La matrice S comporte la valeur propre réelle comme première entrée sur la diagonale et la valeur propre répétée représentée par le bloc inférieur droit de dimension 2 x 2. Les valeurs propres du bloc de 2 x 2 sont également les valeurs propres de A :

eig(S(2:3,2:3))

ans =

   1.0000 + 0.0000i
   1.0000 - 0.0000i

Voir aussi

eig | schur

Sujets associés

• Matrices dans l’environnement MATLAB
• Factorisations
• Valeurs singulières

Sites web externes

• Méthodes matricielles de l’algèbre linéaire (Ressources pédagogiques MathWorks)

Valeurs singulières

Une valeur singulière et les vecteurs singuliers correspondants d’une matrice rectangulaire A sont, respectivement, un scalaire σ et une paire de vecteurs u et v qui satisfont à

$$\begin{array}{l}Av=\sigma u\\ A^{H}u=\sigma v,\end{array}$$

où $$A^H$$ est la transposée hermitienne de A. Les vecteurs singuliers u et v sont généralement mis à l’échelle pour avoir une norme de 1. En outre, si u et v sont des vecteurs singuliers de A, alors -u et -v sont des vecteurs singuliers de A également.

Les valeurs singulières σ sont toujours réelles et non négatives, même si A est complexe. Avec les valeurs singulières dans une matrice diagonale Σ et les vecteurs singuliers correspondants formant les colonnes de deux matrices orthogonales U et V, vous obtenez les équations

$$\begin{array}{l}AV=U\Sigma \\ A^{H}U=V\Sigma .\end{array}$$

Étant donné que U et V sont des matrices unitaires, multiplier la première équation par $$V^H$$ à droite donne l’équation de décomposition en valeurs singulières

$$A=U\Sigma V^H.$$

La décomposition en valeurs singulières complète d’une matrice de dimension m x n implique :

• La matrice U de dimension m x m

• La matrice Σ de dimension m x n

• La matrice V de dimension n x n

En d’autres termes, U et V sont toutes deux carrées, et Σ est de la même taille que A. Si A a beaucoup plus de lignes que de colonnes (m > n), alors la matrice résultante U de dimension m x m est grande. Cependant, la plupart des colonnes de U sont multipliées par des 0 dans Σ. Dans cette situation, la décomposition économique permet d’économiser du temps et du stockage en produisant une matrice U de dimension m x n, une matrice Σ de dimension n x n et la même matrice V :

La décomposition en éléments propres est l’outil approprié pour analyser une matrice lorsqu’elle représente un mapping d’un espace vectoriel dans lui-même, comme c’est le cas pour une équation différentielle ordinaire. Cependant, la décomposition en valeurs singulières est l’outil approprié pour analyser un mapping d’un espace vectoriel dans un autre espace vectoriel, avec potentiellement une dimension différente. La plupart des systèmes d’équations linéaires simultanées appartiennent à cette seconde catégorie.

Si A est carrée, symétrique et définie positive, alors ses décompositions en éléments propres et en valeurs singulières sont identiques. Mais, à mesure que la symétrie et le caractère défini positif de A diminuent, la différence entre les deux décompositions s’accroît. En particulier, la décomposition en valeurs singulières d’une matrice réelle est toujours réelle, mais la décomposition en éléments propres d’une matrice réelle non symétrique peut être complexe.

Pour la matrice donnée comme exemple

A = [9     4
     6     8
     2     7];

la décomposition en valeurs singulières complète est

[U,S,V] = svd(A)

U =

   -0.6105    0.7174    0.3355
   -0.6646   -0.2336   -0.7098
   -0.4308   -0.6563    0.6194


S =

   14.9359         0
         0    5.1883
         0         0


V =

   -0.6925    0.7214
   -0.7214   -0.6925

Vous pouvez vérifier que U*S*V' est égal à A dans la marge d'erreur d'arrondi. Pour ce petit problème, la décomposition réduite est seulement légèrement plus petite.

[U,S,V] = svd(A,"econ")

U =

   -0.6105    0.7174
   -0.6646   -0.2336
   -0.4308   -0.6563


S =

   14.9359         0
         0    5.1883


V =

   -0.6925    0.7214
   -0.7214   -0.6925

Là encore, U*S*V' est égal à A dans la marge d'erreur d'arrondi.

Calcul de la SVD en batch

Si vous devez décomposer un ensemble important de matrices de la même dimension, il ne sera pas efficace d’effectuer toutes les décompositions en boucle avec svd. Mieux vaut concaténer toutes les matrices au sein d’un tableau multidimensionnel et utiliser pagesvd pour effectuer les décompositions en valeurs singulières sur l'ensemble des couches (pages) du tableau avec un seul appel de fonction.

Fonction	Utilisation
pagesvd	Utilisez pagesvd pour effectuer les décompositions en valeurs singulières sur les couches d’un tableau multidimensionnel. C’est une manière efficace d’effectuer une SVD sur un ensemble important de matrices ayant toutes la même dimension.

Prenons l’exemple d’un ensemble de trois matrices de dimension 2 x 2. Concaténez les matrices dans un tableau 2 x 2 x 3 avec la fonction cat.

A = [0 -1; 1 0];
B = [-1 0; 0 -1];
C = [0 1; -1 0];
X = cat(3,A,B,C);

À présent, utilisez pagesvd pour effectuer simultanément les trois décompositions.

[U,S,V] = pagesvd(X);

À chaque couche de X correspondent des couches dans les sorties U, S, et V. Par exemple, la matrice A se trouve sur la première couche de X, et sa décomposition est donnée par U(:,:,1)*S(:,:,1)*V(:,:,1)'.

Approximation SVD de rang faible

Pour les grandes matrices creuses, l’utilisation de svd pour calculer toutes les valeurs singulières et tous les vecteurs singuliers n’est pas toujours pratique. Par exemple, si vous ne devez connaître que quelques-unes des valeurs singulières les plus grandes, le fait de calculer toutes les valeurs singulières d’une matrice creuse de 5 000 x 5 000 représente un travail supplémentaire.

Dans les cas où seul un sous-ensemble des valeurs singulières et des vecteurs singuliers est requis, les fonctions svds et svdsketch sont préférables à svd.

Fonction	Utilisation
svds	Utilisez svds pour calculer une approximation de rang k de la SVD. Vous pouvez spécifier si le sous-ensemble de valeurs singulières doit être le plus grand, le plus petit ou le plus proche d’un nombre spécifique. svds calcule généralement la meilleure approximation de rang k possible.
svdsketch	Utilisez svdsketch pour calculer une SVD partielle de la matrice d’entrée et satisfaisant une tolérance spécifiée. Alors que svds nécessite que vous spécifiiez le rang, svdsketch détermine de manière adaptative le rang de l’esquisse de matrice en fonction de la tolérance spécifiée. L’approximation de rang k que svdsketch utilise satisfait à la tolérance, mais, contrairement à svds, il n’est pas garanti que ce soit la meilleure possible.

Par exemple, supposez une matrice creuse aléatoire de 1 000 x 1 000 avec une densité d’environ 30 %.

n = 1000;
A = sprand(n,n,0.3);

Les six valeurs singulières les plus grandes sont

S = svds(A)

S =

  130.2184
   16.4358
   16.4119
   16.3688
   16.3242
   16.2838

Et les six valeurs singulières les plus petites sont

S = svds(A,6,"smallest")

S =

    0.0740
    0.0574
    0.0388
    0.0282
    0.0131
    0.0066

Pour des matrices plus petites capables de tenir en mémoire sous la forme d’une matrice entière, full(A), l’utilisation de svd(full(A)) peut être plus rapide que svds ou svdsketch. Cependant pour des matrices vraiment grandes et creuses, il devient nécessaire d’utiliser svds ou svdsketch.

Voir aussi

svd | svds | svdsketch | gsvd | pagesvd

Sujets associés

• Matrices dans l’environnement MATLAB
• Factorisations
• Valeurs propres