20

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

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

谢谢

4

1 回答 1

40

对于x = A\b，反斜杠运算符包含许多算法来处理不同类型的输入矩阵。因此对矩阵A进行诊断，并根据其特征选择执行路径。

以下页面以伪代码描述了何时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 分解。对于方三角矩阵，它执行简单的前向/后向替换。对于平方对称正定矩阵，使用Cholesky 分解。否则LU 分解用于一般方阵。

更新： MathWorks用一些漂亮的流程图更新了文档页面中的算法部分。mldivide请参见此处和此处（完整和稀疏案例）。

mldivide_full

所有这些算法在LAPACK中都有相应的方法，实际上这可能是 MATLAB 正在做的事情（请注意，最新版本的 MATLAB 附带优化的英特尔 MKL实现）。

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

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

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

x = pinv(A)*b

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

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

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

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

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

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

于 2013-08-31T23:30:22.287 回答