如何实现 Matlab 的 mldivide（又名反斜杠运算符“\”）

Question

我目前正在尝试开发一个小型的面向矩阵的数学库（我使用Eigen 3 进行矩阵数据结构和运算）并且我想实现一些方便的 Matlab 函数，例如广泛使用的反斜杠运算符 (等价于mldivide)，用于计算线性系统的解（以矩阵形式表示）。

有没有关于如何实现这一点的详细解释？ （我已经用经典的 SVD 分解实现了 Moore-Penrose 伪逆 pinv 函数，但我在某处读到 A\b 并不总是 pinv(A)*b ，至少 Matalb 不会简单地这样做）

谢谢

Answer 1

对于x = A\b，backslash 运算符包含多个algorithms 来处理不同类型的输入矩阵。因此矩阵A被诊断出来，并根据其特征选择执行路径。

以下page 以伪代码描述A 为全矩阵时：

if size(A,1) == size(A,2)         % A is square
    if isequal(A,tril(A))         % A is lower triangular
        x = A \ b;                % This is a simple forward substitution on b
    elseif isequal(A,triu(A))     % A is upper triangular
        x = A \ b;                % This is a simple backward substitution on b
    else
        if isequal(A,A')          % A is symmetric
            [R,p] = chol(A);
            if (p == 0)           % A is symmetric positive definite
                x = R \ (R' \ b); % a forward and a backward substitution
                return
            end
        end
        [L,U,P] = lu(A);          % general, square A
        x = U \ (L \ (P*b));      % a forward and a backward substitution
    end
else                              % A is rectangular
    [Q,R] = qr(A);
    x = R \ (Q' * b);
end

对于非方阵，使用QR decomposition。对于方三角矩阵，它执行简单的forward/backward substitution。对于平方对称正定矩阵，使用Cholesky decomposition。否则LU decomposition 用于一般方阵。

更新： MathWorks 更新了mldivide 文档页面中的algorithm section 并添加了一些漂亮的流程图。请参阅 here 和 here（完整和稀疏案例）。

所有这些算法在@​​987654331@ 中都有对应的方法，实际上这可能是 MATLAB 正在做的事情（请注意，最新版本的 MATLAB 附带了优化的 Intel MKL 实现）。

使用不同方法的原因是它试图使用最具体的算法来求解利用系数矩阵的所有特征的方程组（因为它会更快或更数值稳定）。所以你当然可以使用通用求解器，但它不会是最有效的。

事实上，如果您事先知道A 是什么样的，您可以通过调用linsolve 并直接指定选项来跳过额外的测试过程。

如果A 是矩形或奇异的，您还可以使用PINV 来找到最小范数最小二乘解（使用SVD decomposition 实现）：

x = pinv(A)*b

---

以上所有内容都适用于稠密矩阵，稀疏矩阵则完全不同。通常在这种情况下使用iterative solvers。我相信 MATLAB 将 UMFPACK 和 SuiteSpase 包中的其他相关库用于直接求解器。

使用稀疏矩阵时，您可以打开诊断信息并查看使用spparms 执行的测试和选择的算法：

spparms('spumoni',2)
x = A\b;

---

此外，反斜杠运算符也适用于 gpuArray，在这种情况下，它依赖于 cuBLAS 和 MAGMA 在 GPU 上执行。

它也为distributed arrays 实现，它在分布式计算环境中工作（工作在一组计算机之间分配，每个工作人员只有数组的一部分，可能整个矩阵不能一次全部存储在内存中） .底层实现使用ScaLAPACK。

如果你想自己实现所有这些，这是一个相当高的要求:)