计算科学 - 为什么 LU 分解的各个部分的速度如此不同？ - 吾爱随笔录

为什么 LU 分解的各个部分的速度如此不同？

计算科学线性求解器稀疏矩阵矩阵表现

2021-12-22 00:06:08

我知道解决矩阵问题的简单方法

A \cdot x = b

$A\cdot x=b$ 是 LU 分解

\begin{aligned} L, U & = lu (A) \\ y & = solve (L, b) \\ x & = solve (U, y) \end{aligned}

$\begin{split} L,\,U&=\text{lu}(A)\\ y&=\text{solve}(L,\,b)\\ x&=\text{solve}(U,\,y) \end{split}$ 因此我在matlab中实现了以下代码：

r_num = 10000;
r_rounds = 10000;

A = -30*sparse(diag(ones(r_num,1),0))+16*sparse(diag(ones(r_num-1,1),1))+16*sparse(diag(ones(r_num-1,1),-1))-sparse(diag(ones(r_num-2,1),-2))-sparse(diag(ones(r_num-2,1),2));

b = rand(r_num,1);

tic;
for i = 1:r_rounds
    c1 = A\b;
end
toc;

tic;
[L, U] = lu(A);
for i = 1:r_rounds
    y = L\b;

end
toc;
tic;
for i = 1:r_rounds
    c2 = U\y;
end
toc;

我已经知道了 $A\setminus b$ 正在后台进行 LU 分解，但我仍然对结果感兴趣。现在我的脚本返回了计算 $L\setminus b$ 几乎和直接计算一样昂贵 $A\setminus b$ ，而计算 $U\setminus y$ 非常快。但差异从何而来？还是我的测量方法有缺陷？

当分别计时 LU 分解时，我得到

A\b
Elapsed time is 13.380155 seconds.
lu(A)
Elapsed time is 0.002662 seconds.
y=L\b
Elapsed time is 10.557989 seconds.
c2=U\y
Elapsed time is 1.229101 seconds.

单独计时 LU 分解的原因（在循环之外）是因为对于我以后的应用程序，我将不得不处理一个常量 $A$ , 但变化量很大 $b$ -价值观。因此，我尝试将循环中的元素数量保持在尽可能低的水平。

1个回答

首先，不要忘记在循环中对 LU 分解进行计时！否则，这不是一个公平的比较。如果我这样做，我会得到以下时间：

A\b:   Elapsed time is 10.294272 seconds.
lu(A): Elapsed time is 26.675160 seconds.
L\b:   Elapsed time is 9.934748 seconds.
U\y:   Elapsed time is 1.469325 seconds.

您可以看到两者A\b和lu(A)的数量级相同，这是意料之中的。由于本例中的矩阵实际上是对称的，因此您还希望 Matlab不会进行 LU 分解。替换lu为- 非常接近 Matlab 的反斜杠chol。10.067633 seconds

（其实你可以通过设置询问Matlab它在做什么spparms('spumoni',2)，它会告诉你A\b：

sp\: bandwidth = 2+1+2.
sp\: is A diagonal? no.
sp\: is band density (1) > bandden (0.5) to try banded solver? yes.
sp\: is LAPACK's banded solver successful? yes.

所以它实际上是使用专门的求解器。）

要了解为什么L\b比慢得多U\y，看看两者的稀疏模式是有启发性的 via spy：

可以看到，它的稀疏模式比（上三对角线，对应的稀疏模式）L复杂得多，因此相应的替换步骤会慢得多。事实上，再次询问 Matlab 它在做什么：UAL\b

sp\: bandwidth = 2+1+707.
sp\: is A diagonal? no.
sp\: is band density (0.0043793) > bandden (0.5) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? yes.
sp\: permute and solve.

所以你看到这里正在做额外的工作。

其它你可能感兴趣的问题

上一篇浮点算术中的关联性因两个值而失败下一篇从 FEM 软件中提取系统矩阵