Ejemplo práctico de por qué no es bueno invertir una matriz

16

Soy consciente de que invertir una matriz para resolver un sistema lineal no es una buena idea, ya que no es tan preciso y eficiente como resolver directamente el sistema o usar la descomposición LU, Cholesky o QR.

Sin embargo, no he podido comprobar esto con un ejemplo práctico. He intentado este código (en 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))

y los residuos son siempre del mismo orden (10 ^ -13).

¿Podría alguien proporcionar un ejemplo práctico en el que inv (A) * b sea mucho menos impreciso que A \ b?

------ Actualización de la pregunta ------

Gracias por sus respuestas. Sin embargo, supongamos que tenemos que resolver veces un sistema A x = b , donde A es siempre la misma matriz. Considere esonorteUNX=siUN

- es completa, y por lo tanto A - 1 requiere el mismo almacenamiento de memoria de A .UNUN-1UN

-El número de condición de es pequeño, por lo tanto, A - 1 se puede calcular con precisión.UNUN-1

En ese caso, ¿no sería más eficiente calcular lugar de usar una descomposición LU? Por ejemplo, he probado este código de Matlab:UN-1

%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))

Para una matriz con la condición número de 450, los residuos son en ambos casos, pero se necesita 19 segundos para resolver el sistema de n veces usando la descomposición LU, mientras que usando la inversa de A sólo toma 9 segundos.O(10-11)

Manu
fuente
8
la página de ayuda de MATLAB para inv da un buen ejemplo. Mire debajo de la sección titulada Solve Linear System .
GoHokies
1
por cierto, ¿cuál es el número de condición de su matriz ? No tengo MATLAB en mi PC en el trabajo, así que no puedo comprobarlo, pero supongo que es lo suficientemente pequeño como para que puedas obtener un inverso preciso ...UN
GoHokies
2
Eché un vistazo a Trefethen y Bau (ejercicio 21.4), y lo describen únicamente en términos de costo de cómputo, flops vs. 22norte3flops. Entonces, incluso si encuentra residuales similares (¿ha intentado verificar matrices más mal acondicionadas, como en el comentario de GoHokies?), El costo de cómputo innecesario solo probablemente valga la pena el consejo. 23norte3
Kirill
3
El tamaño de su matriz es demasiado pequeño y está bien acondicionado para esta comparación. No es que no haya problemas relevantes donde tenga tales matrices, pero la opinión recibida de que no debe invertir está destinada a un entorno diferente (por ejemplo, el mencionado por Chris Rackauckas en su respuesta). De hecho, para matrices pequeñas y, certificablemente, bien acondicionadas, calcular la inversa puede ser la mejor opción. Un caso extremo sería matrices de rotación 3x3 (o, más realista, transformación afín).
Christian Clason
1
Si necesita resolver repetidamente Ax=bcon lo mismo Ay es lo suficientemente pequeño como para tomar el inverso, puede guardar la factorización LU y reutilizarla.
Chris Rackauckas

Respuestas:

11

Normalmente hay algunas razones principales para preferir resolver un sistema lineal con respecto al uso del inverso. Brevemente:

  • problema con el número condicional (comentario de @GoHokies)
  • problema en el caso escaso (respuesta @ChrisRackauckas)
  • eficiencia (comentario de @Kirill)

De todos modos, como comentó @ChristianClason en los comentarios, puede haber algunos casos en los que el uso del inverso sea una buena opción.

En la nota / artículo de Alex Druinsky, Sivan Toledo, ¿Qué precisión tiene inv (A) * b? Hay alguna consideración sobre este problema.

X

inversoEl |El |XV-XEl |El |O(κ2(UN)ϵmetrounChyonortemi) estable hacia atrás (LU, QR, ...)El |El |XsiunCkwunrre-stunsilmi-XEl |El |O(κ(UN)ϵmetrounChyonortemi)

XV

V

V

VEl |El |XVEl |El |El |El |XEl |El |

siUN

Entonces, la oportunidad de usar o no el inverso depende de la aplicación, puede consultar el artículo y ver si su caso cumple la condición para obtener la estabilidad hacia atrás o si no la necesita.

En general, en mi opinión, es más seguro resolver el sistema lineal.

Mauro Vanzetto
fuente
12

Δtu

tut=Δtu+F(t,tu).

UN

tut=UNtu+F(t,tu)

UNyo-γUNSpecialMatrices.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

norteO(3norte)O(1)

Sin embargo, digamos que queremos invertir la matriz.

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(norte2)

\IterativeSolvers.jlUNX=siUN-1UN

Como otros han mencionado, el número de condición y el error numérico es otra razón, pero el hecho de que el inverso de una matriz dispersa sea denso da una clara "esta es una mala idea".

Chris Rackauckas
fuente