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

我目前正在开发一个小型的面向matrix的math库（我使用Eigen 3进行matrix数据结构和操作），我想实现一些方便的Matlab函数，比如广泛使用的反斜杠运算符（相当于mldivide ）来计算线性系统的解（用matrixforms表示）。

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

谢谢

对于x = A\b ， 反斜杠运算符包含许多algorithm来处理不同types的inputmatrix。 因此，matrixA被诊断并且根据其特征select执行path。

当A是一个完整的matrix时，下面的页面用伪代码描述： 

 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

对于非正方形matrix，使用QR分解 。 对于正方形三angular形matrix，它执行简单的向前/向后replace 。 对于平方对称正定matrix，使用Cholesky分解 。 否则LU分解被用于一般的正方形matrix。

更新： MathWorks更新了mldivide doc页面中的algorithm部分 ， mldivide带有一些不错的stream程图。 看到这里和这里 （完整和稀疏的情况下）。

mldivide_full

所有这些algorithm在LAPACK中都有相应的方法，实际上它可能是MATLAB正在做的（请注意，MATLAB的最新版本随优化的英特尔MKL实现一起提供）。

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

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

如果A是矩形或单数，也可以使用PINV来find一个最小范数最小二乘解（用SVD分解来实现）： 

 x = pinv(A)*b

以上所有适用于密集matrix，稀疏matrix是一个完全不同的故事。 通常在这种情况下使用迭代求解器 。 我相信MATLAB使用UMFPACK和SuiteSpase包中的其他相关库来直接求解。

使用稀疏matrix时，可以打开诊断信息，并使用spparms查看执行的testing和select的algorithm： 

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

而且，反斜线运算符也可以在gpuArray上运行 ，在这种情况下，它依靠cuBLAS和MAGMA在GPU上执行。

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

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

Interesting Posts