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

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 :

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 :

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 :

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 :

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.

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

x =
      10
     -12
       5

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

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

P = pinv(A)*b

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

rank(A)

ans =

     2

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

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

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

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 :

À 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.

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 :

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é.