为什么使用迭代方法:AMG 预处理 PCG 比 Matlab 直接方法 'A\b' 慢?

计算科学 matlab 线性求解器 稀疏矩阵 迭代法 预处理
2021-12-12 23:50:52

最近遇到一个问题

对于大型线性系统,有句俗话说:由于直接方法的内存问题,需要迭代方法。

但是当我在Matlab中进行一些实验时,我发现了一个奇怪的现象:我在2D中使用泊松方程并使用Q1元素,网格输入为10,我得到了一个系统: where is,大而稀疏。

Ax=f
A1050625×1050625

原则上,我们应该使用迭代方法,例如 PCG 或 minres 内置 Matlab,与 AMG preconditioner。但是当我A\f在命令窗口输入时,Matlab 直接法只需要秒,速度很快。4.588234

然后,我想用 AMG 预处理器测试 PCG。我发现设置 AMG 预调节器的时间非常长。一开始以为是系统规模不够大,然后我用grid input=11,内存坏了,无法生成 ,所以在我的PC上,我无法得到结果对于大型稀疏系统,迭代方法优于直接方法。为什么?是不是尺寸不够的原因?但更大的系统无法在我们的个人电脑中生成。A

我也认为对于大型稀疏系统,迭代方法是必要的,但数值结果却相反:Matlab 的A\b速度如此之快。

我们应该如何理解“大型稀疏系统迭代方法优于直接方法”的说法?你能给我一些建议吗?

2个回答

在这个实验中有几件事需要考虑:

为什么 Matlab 稀疏直接可能“这么快”:

(针对您的特定测试)

  • 在二维(当然,取决于问题)中,您的矩阵在 FEM 离散化后产生,经过一些重新排序后可能看起来是“接近带状”结构。的带宽越小,稀疏直接方法处理它的效率就越高。如果带宽足够小,Matlab 可能会选择专门的带状求解器AA
  • 我不知道您使用的是什么特定的代数多重网格 (AMG) 预处理器。虽然它内部可能存在内部性能问题,但 AMG 本身对于您的问题可能是过度杀伤力。
  • 稀疏线性求解器现在是一回事。过去几年取得了很大进展,您所看到的至少部分效果应该归功于它。

为什么人们可能想要使用迭代求解器,无论如何:

  • 一般来说,即使是分解过程中最好的稀疏直接求解器,仍然会生成填充。因此,与原始矩阵相比,分解后的矩阵将占用[显着]更多空间,这对于大型问题来说肯定是一个问题。
  • 迭代求解器的加速是一个比较成熟的领域。在不同的应用领域,有许多快速算法(例如,基于快速傅里叶变换 - FFT、快速多极子方法 - FMM 等)可以加速矩阵向量乘积,自然适合迭代线性求解器路线。

要尝试的事情:

  • 尝试使用带有香草对角线预处理器的迭代求解器来解决您的问题。看,迭代次数是多少,并决定是否值得考虑更好的预处理器(如 AMG)。
  • 如果这符合您的应用领域的兴趣,请考虑进行 3-D 测试。
  • 检查由二维泊松产生的带宽和稀疏模式。可能会发生,您正在解决一个非常特殊的情况,而不是一般的稀疏矩阵。A

如何理解“说法”:

对于大型线性系统:由于直接方法的内存问题,需要迭代方法。

你可能想批判性地看待它。几乎与任何一般性的一揽子声明一样。

虽然迭代方法有其优势,但某些问题仍然需要直接方法。此外,还有许多快速直接方法,其中也可以加速分解(例如,应用于 FEM 的分层矩阵)。我想说,在过去的 15 年里,非加速稀疏线性求解器领域变得更有前途。所以,这句话在 20 年前可能是教条,现在至少要弱得多。

谢谢大家的关注。以下是一位教授的回复:

MATLAB 稀疏求解器是求解与二维拉普拉斯算子相关的线性系统的一种非常有效的方法。原因之一是 CHOLMOD 求解器是非常有效的多线程,因此它可以在求解过程中使用所有可用的处理器。例如,我的 Apple 笔记本电脑是 I9 六核架构,当我解决您在下面讨论的问题时,我可以看到所有六核都得到了充分利用。相比之下,AMG 网格设置是解释代码,正如您所观察到的,在 MATLAB 环境中非常慢。但是,它的内存效率很高。

我已经尝试了 3D 数值实验,使用 5 点差来离散泊松方程:

{Δu=f,(x,y,z)G=(1,1)3,u=g,(x,y,z)G.

当系统大小变为 1,000,000 X 1,000,000 时,matlab 命令 A\b 内存不足。Matlab代码如下:

%%  poisson in 2D and 3D 5 points difference matrix 
clc;clear;
n=10;
e=ones(n,1);
B = [-1 2 -1].*e;
d = [-1 0 1];
Tn = spdiags(B,d,n,n);
e=ones(n-1,1);
I = speye(n);
%  2D
Tn_I = kron(Tn,I);
I_Tn = kron(I,Tn);
A = Tn_I+I_Tn;
%  3D
A = kron(Tn_I,I)+kron(I,Tn_I)+kron(I,I_Tn);
b = sum(A,2);
tic;
A\b;toc