如何实现Matlab的mldivide(又名反斜杠运算符"\")

Question

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

有关如何实现这一目标的详细解释吗？(我已经用经典的SVD分解实现了Moore-Penrose伪反向pinv函数,但我读过的地方A\b并不总是如此pinv(A)*b,至少Matalb并不是那么简单)

谢谢

Answer 1

因为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已使用一些不错的流程图更新了doc页面中的算法部分mldivide.看这里和这里(完整和稀疏的情况).

所有这些算法在LAPACK中都有相应的方法,事实上它可能就是MATLAB正在做的事情(注意最近版本的MATLAB附带优化的英特尔MKL实现).

使用不同方法的原因在于它尝试使用最具体的算法来求解利用系数矩阵的所有特征的方程组(因为它将更快或更数值稳定).所以你当然可以使用一般解算器,但它不是最有效的.

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

如果A是矩形或单数,您还可以使用PINV查找最小范数最小二乘解(使用SVD分解实现):

x = pinv(A)*b

以上所有都适用于密集矩阵,稀疏矩阵是完全不同的故事.通常在这种情况下使用迭代求解器.我相信MATLAB使用来自SuiteSpase包的UMFPACK和其他相关库来进行直接求解.

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

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

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

它也适用于在分布式计算环境中工作的分布式阵列(工作在一组计算机之间划分,其中每个工作者只有部分数组,可能整个矩阵不能同时存储在内存中).底层实现使用ScaLAPACK.

如果你想自己实现所有这些,这是一个非常高的顺序:)