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

16

Мне известно о том, что инвертировать матрицу для решения линейной системы не очень хорошая идея, поскольку она не так точна и эффективна, как непосредственное решение системы или использование разложения 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?

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

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

- полно, и , таким образом требует такого же память для хранения , чем .AA1A

-Кондиционное число мало, поэтому может быть вычислено с точностью.AA1

В таком случае, не будет ли эффективнее вычислить , чем использовать разложение LU? Например, я пробовал этот код Matlab:A1

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

Manu
источник
8
страница справки MATLAB для inv дает хороший пример. Посмотрите в разделе « Решить линейную систему» .
GoHokies
1
Кстати, каков номер вашей матрицы ? У меня нет MATLAB на моем ПК на работе, поэтому я не могу проверить, но я предполагаю, что он достаточно мал, чтобы вы могли получить точное обратное значение ...A
GoHokies
2
Я взглянул на Trefethen и Bau (упражнение 21.4), и они описывают это чисто с точки зрения стоимости вычислений, флопов против 22n3флопа. Таким образом, даже если вы находите схожие остатки (пытались ли вы проверить более плохо обусловленные матрицы, как в комментарии GoHokies?), Одна только ненужная стоимость вычислений, вероятно, стоит совета. 23n3
Кирилл
3
Ваш размер матрицы слишком мал и хорошо подготовлен для этого сравнения. Не то чтобы не было соответствующих проблем, когда у вас есть такие матрицы, но полученное мнение о том, что вы не должны инвертировать, предназначено для другой настройки (например, той, которая была упомянута Крисом Ракауцкасом в его ответе). На самом деле, для небольших и - конечно же - хорошо обусловленных матриц, вычисление обратного может быть лучшим вариантом. В крайнем случае это матрицы вращения 3х3 (или, более реалистично, аффинного преобразования).
Кристиан
1
Если вам нужно многократно решать Ax=bс одним Aи тем же, и он достаточно мал, чтобы принять обратное, вы можете вместо этого сохранить LU-факторизацию и использовать ее повторно.
Крис

Ответы:

11

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

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

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

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

x

inverse||xVx||O(κ2(A)ϵmachine) backward stable (LU, QR,...)||xbackwardstablex||O(κ(A)ϵmachine)

xV

V

V

V||xV||||x||

bA

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

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

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

Δu

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

A

ut=Au+f(t,u)

AIγASpecialMatrices.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)

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

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.jlAx=bA1A

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

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