计算科学 - 为什么反转矩阵不好的实际示例 - 吾爱随笔录

为什么反转矩阵不好的实际示例

计算科学矩阵准确性逆

2021-12-18 22:15:17

我知道对矩阵求逆来求解线性系统并不是一个好主意，因为它不如直接求解系统或使用 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}$ 可以准确计算。

在那种情况下，计算不是更有效吗 $A^{-1}$ 而不是使用LU分解？例如，我试过这个 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))

对于条件数约为 450 的矩阵，残差为 $O(10^{-11})$ 在这两种情况下，但使用 LU 分解求解系统 n 次需要 19 秒，而使用 A 的逆分解只需要 9 秒。

2个回答

这是一个与 PDE 中的内存使用相关的非常实用的快速示例。当离散拉普拉斯算子时， $\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)$

典型的一维示例是 Strang 矩阵。隐式方法需要反转 $A$ ，或某种形式的 $I-\gamma A$ . 对于空间中有 5 个点的离散化，让我们看看这个算子会发生什么。我们可以使用 Julia 在 Julia 中轻松生成它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

这个特殊的矩阵是三对角矩阵（许多其他离散化是带状的，因此仍然是稀疏的）。这意味着您可以通过仅存储 3 个数组来存储它，或者在这种情况下，它可以是“惰性的”（即，不需要数组，您可以拥有一个按需生成值的“伪数组”类型）。这意味着即使对于非常大的空间离散化 $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

请注意，这个稀疏矩阵的逆矩阵是稠密的。因此，如果我们不知道逆的解析解（对于 PDE 离散化产生的大多数稀疏矩阵都是如此），则逆所需的内存量为 $\mathcal{O}(n^2)$ . 这意味着逆矩阵会占用大量额外内存，因此您的内存限制将由密集逆矩阵的大小决定。

相反，直接求解器方法和迭代求解\器（如IterativeSolvers.jl $Ax=b$ 无需计算 $A^{-1}$ ，因此只有必须存储的内存要求 $A$ . 这可以极大地扩展您可以求解的 PDE 的大小。

正如其他人所提到的，条件数和数值误差是另一个原因，但是稀疏矩阵的逆是密集的这一事实给出了一个非常明确的“这是一个坏主意”。

通常有一些主要原因倾向于解决线性系统而不是使用逆。简要地：

条件数问题（@GoHokies 评论）
稀疏情况下的问题（@ChrisRackauckas 回答）
效率（@Kirill 评论）

无论如何，正如@ChristianClason 在评论中所说，在某些情况下使用倒数是一个不错的选择。

在 Alex Druinsky、Sivan Toledo 的注释/文章中，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})$

现在可以改进对逆的估计，在逆的某些条件下，参见论文中的定理 1，但是 $x_V$ 可以是有条件的准确的，而不是向后稳定的。

该论文显示了发生这种情况的情况（ $V$ 是相反的）

(1) $V$ 不是一个好的右逆，或者

(2) $V$ 是这样一个可怜的左逆 $||x_V||$ 远小于 $||x||$ ，或者

(3) 投影 $b$ 在左侧奇异向量 $A$ 与小的奇异值相关联的值很小。

所以使用与否的机会取决于应用程序，您可以查看文章，看看您的案例是否满足获得后向稳定性的条件，或者您是否不需要它。

一般来说，在我看来，解决线性系统更安全。

其它你可能感兴趣的问题

上一篇哪些方法可以确保物理量在整个 PDE 模拟中保持为正？下一篇小波如何应用于 PDE？