Exemplo prático de por que não é bom inverter uma matriz

Estou ciente de que inverter uma matriz para resolver um sistema linear não é uma boa ideia, pois não é tão preciso e eficiente como solucionar diretamente o sistema ou usar a decomposição LU, Cholesky ou QR.

No entanto, não pude verificar isso com um exemplo prático. Eu tentei este código (em 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))

e os resíduos são sempre da mesma ordem (10 ^ -13).

Alguém poderia fornecer um exemplo prático em que inv (A) * b é muito menos impreciso que A \ b?

------ Atualização da pergunta ------

Obrigado por suas respostas. No entanto, suponha que tenhamos que resolver vezes um sistema , onde é sempre a mesma matriz. Considere isso $n$ $Ax = b$ $A$

- está cheio, e, assim, requer o mesmo armazenamento de memória do que . $A$ $A^{-1}$ $A$

-O número da condição de é pequeno, portanto pode ser calculado com precisão. $A$ $A^{-1}$

Nesse caso, não seria mais eficiente calcular do que usar uma decomposição de LU? Por exemplo, eu tentei este código do Matlab: $A^{-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 uma matriz com o número de condição em torno de 450, os resíduos são nos dois casos, mas são necessários 19 segundos para resolver o sistema n vezes usando a decomposição da LU, enquanto que o inverso de A leva apenas 9 segundos. $O(10^{-11})$

matrix accuracy inverse Manu
fonte

a página de ajuda do MATLAB para inv fornece um bom exemplo. Procure na seção intitulada Solve Linear System .

GoHokies

btw, qual é o número da condição da sua matriz ? Eu não tenho MATLAB no meu PC no trabalho, então eu não posso verificar, mas eu presumo que é pequeno o suficiente para você obter o inverso preciso ...

A

$A$

GoHokies

Dei uma olhada em Trefethen e Bau (exercício 21.4), e eles descrevem apenas os termos do custo de computação, flops vs. flops. Portanto, mesmo que você encontre resíduos semelhantes (você tentou verificar matrizes mais mal condicionadas, como no comentário da GoHokies?), Apenas o custo desnecessário da computação provavelmente vale o conselho.

2 n^{3}

$2n^3$

\frac{2}{3} n^{3}

$\frac23 n^3$

21417 Kirill

O tamanho da sua matriz é muito pequeno e está bem condicionado para essa comparação. Não que não haja problemas relevantes em que você tenha essas matrizes, mas a opinião recebida de que você não deve inverter se destina a uma configuração diferente (por exemplo, a mencionada por Chris Rackauckas em sua resposta). De fato, para matrizes pequenas e - certificadamente - bem condicionadas, calcular o inverso pode de fato ser a melhor opção. Um caso extremo seria matrizes de rotação 3x3 (ou, mais realista, transformação afim).

Christian Clason

Se você precisar resolver repetidamente Ax=bcom o mesmo Ae for pequeno o suficiente para suportar o inverso, poderá salvar a fatoração da LU e reutilizá-la.

Chris Rackauckas 19/03

Respostas:

Normalmente, existem algumas razões principais para preferir resolver um sistema linear em relação ao uso do inverso. Resumidamente:

problema com o número condicional (comentário do @GoHokies)
problema no caso escasso (@ChrisRackauckas answer)
eficiência (@Kirill comment)

De qualquer forma, como o @ChristianClason observou nos comentários, pode haver alguns casos em que o uso do inverso é uma boa opção.

Na nota / artigo de Alex Druinsky, Sivan Toledo, Quão preciso é inv (A) * b? Há alguma consideração sobre esse problema.

De acordo com o artigo, a principal razão para a preferência geral usar resolver o sistema linear está dentro dessas duas estimativas ( é a solução verdadeira): $x$

inverso | | x_{V} - x | | \leq O (κ^{2} (UMA) ϵ_{m uma c h Eu n e}) estável para trás (LU, QR, ...) | | x_{b uma c k W uma r d - s t uma b eu e} - x | | \leq O (κ (UMA) ϵ_{m uma c h Eu n e})

$\text{inverse} \quad || x_V - x || \leq O(\kappa^2(A) \epsilon_{machine})\\ \text{ backward stable (LU, QR,...)} \quad || x_{backward-stable} - x || \leq O(\kappa(A) \epsilon_{machine})$

Agora, a estimativa para o inverso pode ser melhorada, sob alguma condição sobre o inverso, veja o teorema 1 no artigo, mas pode ser condicionalmente preciso e não ser estável para trás. $x_V$

O artigo mostra os casos em que isso acontece ( é o inverso) $V$

(1) não é um bom inverso à direita, ou $V$

(2) é uma inversa à esquerda tão fraca queé muito menor queou $V$ $||x_V||$ $||x||$

(3) a projeção de nos vetores singulares esquerdos de associados a pequenos valores singulares é pequena. $b$ $A$

Portanto, a oportunidade de usar ou não o inverso depende do aplicativo, você pode verificar o artigo e ver se o seu caso satisfaz a condição de obter a estabilidade para trás ou se não precisa.

Em geral, na minha opinião, é mais seguro resolver o sistema linear.

Mauro Vanzetto
fonte

$\Delta u$

{você}_{t} = Δ você + f (t, você) .

$u_t = \Delta u + f(t,u) .$

$A$

{você}_{t} = UMA você + f (t, você)

$u_t = Au + f(t,u)$

$A$ $I-\gamma A$ SpecialMatrices.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

$n$ $\mathcal{O}(3n)$ $\mathcal{O}(1)$

No entanto, digamos que queremos inverter a 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

$\mathcal{O}(n^2)$

\IterativeSolvers.jl $Ax=b$ $A^{-1}$ $A$

Como outros já mencionaram, o número da condição e o erro numérico é outra razão, mas o fato de a inversa de uma matriz esparsa ser densa dá uma ideia muito clara "essa é uma má idéia".

Chris Rackauckas
fonte