Exemple pratique de pourquoi il n'est pas bon d'inverser une matrice

16

Je suis conscient du fait qu'inverser une matrice pour résoudre un système linéaire n'est pas une bonne idée, car il n'est pas aussi précis et aussi efficace que de résoudre directement le système ou d'utiliser la décomposition LU, Cholesky ou QR.

Cependant, je n'ai pas pu vérifier cela avec un exemple pratique. J'ai essayé ce code (dans MATLAB)

M   = 500;    
A   = rand(M,M);
A   = real(expm(1i*(A+A.')));
b   = rand(M,1);

x1  = A\b;
x2  = inv(A)*b;

disp(norm(b-A*x1))
disp(norm(b-A*x2))

et les résidus sont toujours du même ordre (10 ^ -13).

Quelqu'un pourrait-il fournir un exemple pratique dans lequel inv (A) * b est beaucoup moins inexact que A \ b?

------ Mise à jour de la question ------

Merci pour vos réponses. Supposons cependant que nous devions résoudre fois un système , où est toujours la même matrice. Considérez queA x = b AnAx=bA

- A est plein, et donc A1 nécessite la même mémoire de stockage que A .

-Le nombre de conditions de A est petit, donc A1 peut être calculé avec précision.

Dans ce cas, ne serait-il pas plus efficace de calculer A1 plutôt que d'utiliser une décomposition LU? Par exemple, j'ai essayé ce code Matlab:

%Set A and b:
M           = 1000; 
A           = rand(M,M);
A           = real(expm(1i*(A+A.')));
b           = rand(M,1);

%Times we solve the system:
n           = 3000;

%Performing LU decomposition:
disp('Performing LU decomposition')
tic
[L,U,P]     = lu(A);
toc
fprintf('\n')

%Solving the system n times with LU decomposition:
optsL.LT    = true;   %Options for linsolve
optsU.UT    = true;
disp('Solving the system n times using LU decomposition')
tic
for ii=1:n
    x1      = linsolve(U, linsolve(L,P*b,optsL) , optsU);
end
toc
fprintf('\n')

%Computing inverse of A:
disp('Computing inverse of A')
tic
Ainv        = inv(A);
toc
fprintf('\n')

%Solving the system n times with Ainv:
disp('Solving the system n times with A inv')
tic
for ii=1:n
    x2  = Ainv*b;
end
toc
fprintf('\n')

disp('Residuals')
disp(norm(b-A*x1))
disp(norm(b-A*x2))

disp('Condition number of A')
disp(cond(A))

Pour une matrice dont le nombre de conditions est d'environ 450, les résidus sont O(1011) dans les deux cas, mais il faut 19 secondes pour résoudre le système n fois en utilisant la décomposition LU, alors qu'en utilisant l'inverse de A, cela prend seulement 9 secondes.

Manu
la source
8
la page d'aide MATLAB pour inv en donne un bon exemple. Regardez sous la section intitulée Solve Linear System .
GoHokies
1
btw, quel est le numéro de condition de votre matrice ? Je n'ai pas MATLAB sur mon PC au travail donc je ne peux pas vérifier, mais je suppose qu'il est assez petit pour que vous obteniez une inverse précise ...A
GoHokies
2
J'ai jeté un œil à Trefethen et Bau (exercice 21.4), et ils le décrivent purement en termes de coût de calcul, flops contre flops. Donc, même si vous trouvez des résidus similaires (avez-vous essayé de vérifier des matrices plus mal conditionnées, comme dans le commentaire de GoHokies?), Le coût de calcul inutile à lui seul vaut probablement le conseil. 22n323n3
Kirill
3
La taille de votre matrice est beaucoup trop petite et bien conditionnée pour cette comparaison. Non pas qu'il n'y ait pas de problèmes pertinents lorsque vous avez de telles matrices, mais l'opinion reçue que vous ne devriez pas inverser est destinée à un paramètre différent (par exemple, celui mentionné par Chris Rackauckas dans sa réponse). En fait, pour des matrices petites et - certes - bien conditionnées, le calcul de l'inverse peut en effet être la meilleure option. Un cas extrême serait les matrices de rotation 3x3 (ou, plus réaliste, de transformation affine).
Christian Clason
1
Si vous devez résoudre plusieurs fois Ax=ble même problème Aet qu'il est suffisamment petit pour prendre l'inverse, vous pouvez enregistrer la factorisation LU et la réutiliser à la place.
Chris Rackauckas

Réponses:

11

Normalement, il y a quelques raisons principales de préférer résoudre un respect de système linéaire pour utiliser l'inverse. Brièvement:

  • problème avec le numéro conditionnel (commentaire @GoHokies)
  • problème dans le cas rare (réponse @ChrisRackauckas)
  • efficacité (commentaire @Kirill)

Quoi qu'il en soit, comme l'a fait remarquer @ChristianClason dans les commentaires, il peut y avoir des cas où l'utilisation de l'inverse est une bonne option.

Dans la note / l'article d'Alex Druinsky, Sivan Toledo, How Accurate is inv (A) * b? il y a une certaine réflexion sur ce problème.

Selon l'article, la principale raison de la préférence générale pour résoudre le système linéaire est à l'intérieur de ces deux estimations ( est la vraie solution): inverseX

inverse||XV-X||O(κ2(UNE)ϵmunechjene) arrière stable (LU, QR, ...)||Xbuneckwuner-stuneble-X||O(κ(UNE)ϵmunechjene)

Maintenant, l'estimation de l'inverse peut être améliorée, sous certaines conditions par rapport à l'inverse, voir le théorème 1 dans l'article, mais peut être conditionnellement précis et non stable en arrière.XV

Le papier montre les cas où cela se produit ( est l'inverse)V

(1) n'est pas une bonne inverse droite, ouV

(2) est un mauvais inverse gauche tel queest beaucoup plus petit que, ouV||XV||||X||

(3) la projection de sur les vecteurs singuliers gauches de associée à de petites valeurs singulières est petite.bUNE

Ainsi, la possibilité d'utiliser ou non l'inverse dépend de l'application, vous pouvez vérifier l'article et voir si votre cas remplit la condition pour obtenir la stabilité en arrière ou si vous n'en avez pas besoin.

En général, à mon avis, il est plus sûr de résoudre le système linéaire.

Mauro Vanzetto
la source
12

Voici un exemple rapide très pratique lié à l'utilisation de la mémoire dans les PDE. Quand on discrétise l'opérateur de Laplace, , par exemple, dans l'équation de chaleurΔu

ut=Δu+F(t,u).

UNE

ut=UNEu+F(t,u)

UNEje-γUNESpecialMatrices.jl

julia> using SpecialMatrices
julia> Strang(5)
5×5 SpecialMatrices.Strang{Float64}:
 2.0  -1.0   0.0   0.0   0.0
-1.0   2.0  -1.0   0.0   0.0
 0.0  -1.0   2.0  -1.0   0.0
 0.0   0.0  -1.0   2.0  -1.0
 0.0   0.0   0.0  -1.0   2.0

nO(3n)O(1)

Cependant, disons que nous voulons inverser la matrice.

julia> inv(collect(Strang(5)))
5×5 Array{Float64,2}:
 0.833333  0.666667  0.5  0.333333  0.166667
 0.666667  1.33333   1.0  0.666667  0.333333
 0.5       1.0       1.5  1.0       0.5
 0.333333  0.666667  1.0  1.33333   0.666667
 0.166667  0.333333  0.5  0.666667  0.833333

O(n2)

\IterativeSolvers.jlUNEX=bUNE-1UNE

Comme d'autres l'ont mentionné, le numéro de condition et l'erreur numérique sont une autre raison, mais le fait que l'inverse d'une matrice clairsemée soit dense donne une idée très claire "c'est une mauvaise idée".

Chris Rackauckas
la source