Практический пример того, почему не хорошо инвертировать матрицу

Мне известно о том, что инвертировать матрицу для решения линейной системы не очень хорошая идея, поскольку она не так точна и эффективна, как непосредственное решение системы или использование разложения LU, Cholesky или QR.

Однако я не смог проверить это на практическом примере. Я пробовал этот код (в 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))

и остатки всегда одного и того же порядка (10 ^ -13).

Может ли кто-нибудь привести практический пример, в котором inv (A) * b гораздо менее неточно, чем A \ b?

------ Обновление вопроса ------

Спасибо за ваши ответы. Однако предположим, что мы должны решить раз систему , где всегда одна и та же матрица. Считают, что $n$ $Ax = b$ $A$

- полно, и , таким образом требует такого же память для хранения , чем . $A$ $A^{-1}$ $A$

-Кондиционное число мало, поэтому может быть вычислено с точностью. $A$ $A^{-1}$

В таком случае, не будет ли эффективнее вычислить , чем использовать разложение LU? Например, я пробовал этот код 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))

Для матрицы с номером условия около 450 невязки в обоих случаях составляют , но для решения системы n раз с использованием разложения LU требуется 19 секунд, тогда как при использовании обратной величины A это занимает всего 9 секунд. $O(10^{-11})$

matrix accuracy inverse Manu
источник

страница справки MATLAB для inv дает хороший пример. Посмотрите в разделе « Решить линейную систему» .

GoHokies

Кстати, каков номер вашей матрицы

? У меня нет MATLAB на моем ПК на работе, поэтому я не могу проверить, но я предполагаю, что он достаточно мал, чтобы вы могли получить точное обратное значение ...

A

$A$

GoHokies

Я взглянул на Trefethen и Bau (упражнение 21.4), и они описывают это чисто с точки зрения стоимости вычислений,

флопов против

2 n^{3}

$2n^3$

флопа. Таким образом, даже если вы находите схожие остатки (пытались ли вы проверить более плохо обусловленные матрицы, как в комментарии GoHokies?), Одна только ненужная стоимость вычислений, вероятно, стоит совета.

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

$\frac23 n^3$

Кирилл

Ваш размер матрицы слишком мал и хорошо подготовлен для этого сравнения. Не то чтобы не было соответствующих проблем, когда у вас есть такие матрицы, но полученное мнение о том, что вы не должны инвертировать, предназначено для другой настройки (например, той, которая была упомянута Крисом Ракауцкасом в его ответе). На самом деле, для небольших и - конечно же - хорошо обусловленных матриц, вычисление обратного может быть лучшим вариантом. В крайнем случае это матрицы вращения 3х3 (или, более реалистично, аффинного преобразования).

Кристиан

Если вам нужно многократно решать Ax=bс одним Aи тем же, и он достаточно мал, чтобы принять обратное, вы можете вместо этого сохранить LU-факторизацию и использовать ее повторно.

Крис

Ответы:

Обычно есть несколько основных причин, чтобы предпочесть решить линейную систему относительно использования обратного. Кратко:

проблема с условным числом (комментарий @GoHokies)
проблема в редком случае (ответ @ChrisRackauckas)
эффективность (комментарий @Kirill)

В любом случае, как заметил @ChristianClason в комментариях, в некоторых случаях использование обратного является хорошим вариантом.

В заметке / статье Алекса Друинского, Сиван Толедо, Насколько точен inv (A) * b? Есть некоторые соображения по поводу этой проблемы.

$x$

inverse | | x_{V} - x | | \leq O (κ^{2} (A) ϵ_{m a c h i n e}) backward stable (LU, QR,...) | | x_{b a c k w a r d - s t a b l e} - x | | \leq O (κ (A) ϵ_{m a c h i 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})$

$x_V$

$V$

$V$

$V$ $||x_V||$ $||x||$

$b$ $A$

Таким образом, возможность использовать или не использовать обратное зависит от приложения. Можете ли вы проверить статью и посмотреть, удовлетворяет ли ваш случай условию получения обратной стабильности или вам она не нужна.

В общем, на мой взгляд, безопаснее решить линейную систему.

Мауро Ванзетто
источник

$\Delta u$

u_{t} = Δ u + f (t, u) .

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

$A$

u_{t} = A u + f (t, u)

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

Однако, допустим, мы хотим инвертировать матрицу.

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$

Как уже упоминали другие, условное число и числовая ошибка - это еще одна причина, но тот факт, что инверсия разреженной матрицы плотная, дает очень четкое «это плохая идея».

Крис Ракауцкас
источник