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

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

关于如何实现这一点，有没有很好的详细解释？（我已经实现了摩尔定律 - Penrose逆pinv功能与经典的SVD分解，但我在什么地方读到A\b并不总是pinv(A)*b，至少做的Matalb不是简单地做到这一点）

感谢

来源

2013-08-31 865719

这不是简单的，我会强烈建议不要试图。只需将它与LAPACK的DGEMV连接即可。很多聪明人花了很多时间精细调整'mldivide'。 – 2013-08-31 22:22:15

类似的问题：[MATLAB的反斜杠运算符如何求解方形矩阵的Ax = b？]（http://scicomp.stackexchange.com/questions/1001/how-does-the-matlab-backslash-operator-solve-ax -b-for-square-matrices） – Amro

@Amro，woohoo我的问题。 LOL – 2013-09-01 00:01:32

对于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（完整和稀疏的情况）。

所有这些算法在LAPACK相应的方法，实际上它可能是什么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。

，如果你想实现这一切这是一个相当艰巨的任务自己:)

来源

2013-08-31 23:30:22 Amro

谢谢Amro，这是一个非常有建设性的答案。 “A”最可能是一个非方矩阵，所以，如果我理解正确，最合适的方法是使用QR分解。是的，我愿意（至少尝试）实现所需的东西来创建一个自包含的库，所以我并不想使用其他数学包（算法的效率是，目前，不是我的主要目标;））。这是一个个人项目，现在没有什么严重的... – 865719

@ M4XMX：是的，您可以使用QR分解：'Ax = b - > QRx = b - > Rx = Q'b - > x = R \（Q'b）'（最后一步是一个简单的后向替换过程）。请参阅[此处]（https://inst.eecs.berkeley.edu/~ee127a/book/login/l_lineqs_solving.html）以获取解释。这里还有一个相关的[问题]（http://stackoverflow.com/q/7949229/97160），展示了如何使用各种库来解决'Ax = b'：GSL，Armadillo，Eigen。无需重新发明轮子:) – Amro

我认为MATLAB也可能使用[pivoting]（https://en.wikipedia.org/wiki/QR_decomposition#Column_pivoting）来提高数值的准确性。这是分解是：'AP = QR'其中'P'是一个置换矩阵 – Amro

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

回答

相关问题