MATLAB反斜杠运算符如何求解A x = bAx=b对于方阵？

计算科学 线性代数 表现 matlab
2021-12-11 19:50:34

我正在将我的一些代码与“库存”MATLAB 代码进行比较。我对结果感到惊讶。

我运行了一个示例代码（稀疏矩阵）

n = 5000;
a = diag(rand(n,1));
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('Inv(A)*B');
tic;inv(a)*b;toc;

结果 ：

    For a\b
    Elapsed time is 0.052838 seconds.

    For LU
    Elapsed time is 7.441331 seconds.

    For Conj Grad
    Elapsed time is 3.819182 seconds.

    Inv(A)*B
    Elapsed time is 38.511110 seconds.

对于密集矩阵：

n = 2000;
a = rand(n,n);
b = rand(n,1);
disp('For a\b');
tic;a\b;toc;
disp('For LU');
tic;LULU;toc;
disp('For Conj Grad');
tic;conjgrad(a,b,1e-8);toc;
disp('For INV(A)*B');
tic;inv(a)*b;toc;

结果：

For a\b
Elapsed time is 0.575926 seconds.

For LU
Elapsed time is 0.654287 seconds.

For Conj Grad
Elapsed time is 9.875896 seconds.

Inv(A)*B
Elapsed time is 1.648074 seconds.

a\b 怎么这么厉害？

4个回答

在 Matlab 中，'\' 命令调用一种算法，该算法取决于矩阵 A 的结构并包括对 A 属性的检查（小开销）。

  1. 如果 A 是稀疏且带状的，则使用带状求解器。
  2. 如果 A 是上三角矩阵或下三角矩阵，则采用后向替换算法。
  3. 如果 A 是对称的并且具有实数正对角元素，请尝试 Cholesky 分解。如果 A 是稀疏的，则首先使用重新排序以最小化填充。
  4. 如果上述条件均不满足，则使用带有部分旋转的高斯消元进行一般三角分解。
  5. 如果 A 是稀疏的，则使用 UMFPACK 库。
  6. 如果 A 不是正方形，则对未定系统采用基于 QR 分解的算法。

为了减少开销，可以使用 Matlab 中的 linsolve 命令并自己在这些选项中选择合适的求解器。

如果您想查看a\b特定矩阵的作用，您可以设置spparms('spumoni',1)并准确计算出您印象深刻的算法。例如：

spparms('spumoni',1);
A = delsq(numgrid('B',256));
b = rand(size(A,2),1);
mldivide(A,b);  % another way to write A\b

将输出

sp\: bandwidth = 254+1+254.
sp\: is A diagonal? no.
sp\: is band density (0.01) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? yes.
sp\: is CHOLMOD's symbolic Cholesky factorization (with automatic reordering) successful? yes.
sp\: is CHOLMOD's numeric Cholesky factorization successful? yes.
sp\: is CHOLMOD's triangular solve successful? yes.

所以我可以看到在这种情况下“\”最终使用了“CHOLMOD”。

对于稀疏矩阵，Matlab 使用 UMFPACK 进行“ \”操作，在您的示例中，该操作基本上遍历 的值a，反转它们，并将它们与 的值相乘b。但是，对于此示例，您应该使用b./diag(a)，这要快得多。

对于密集系统，反斜杠运算符有点复杂。这里给出了什么时候做的简要描述。根据该描述，在您的示例中，Matlab 将a\b使用反向替换来解决。对于一般方阵，使用 LU 分解。

根据经验，如果您有一个具有合理复杂性的稀疏矩阵（即，它不必是 5 点模板，但实际上可以是每行非零数为的 Stokes 方程的离散化远大于 5），那么如果问题不超过大约 100,000 个未知数，那么像 UMFPACK 这样的稀疏直接求解器通常会击败迭代 Krylov 求解器。

换句话说，对于大多数由 2d 离散化产生的稀疏矩阵，直接求解器是最快的替代方案。只有对于快速获得超过 100,000 个未知数的 3d 问题，才必须使用迭代求解器。

其它你可能感兴趣的问题
上一篇并行 ODE 方法的最新技术是什么？ 下一篇发表强调软件的论文的场所