問題描述

我目前正在嘗試開發(fā)一個(gè)面向矩陣的小型數(shù)學(xué)庫(我正在使用 中都有對應(yīng)的方法，實(shí)際上這可能是MATLAB正在做的(請注意，最新版本的 MATLAB 附帶優(yōu)化的 Intel MKL 實(shí)現(xiàn)).

采用不同方法的原因是它試圖使用最具體的算法來求解利用系數(shù)矩陣所有特征的方程組(或者因?yàn)樗鼤?huì)更快或更數(shù)值穩(wěn)定).所以你當(dāng)然可以使用通用求解器，但它不會(huì)是最有效的.

事實(shí)上，如果您事先知道 A 是什么樣的，您可以通過調(diào)用 linsolve 并直接指定選項(xiàng).

如果 A 是矩形或單數(shù)，你也可以使用 PINV 找到最小范數(shù)最小二乘解(使用 SVD 分解):

x = pinv(A)*b

以上所有內(nèi)容都適用于稠密矩陣，稀疏矩陣則完全不同.通常迭代求解器用于此類案件.我相信 MATLAB 使用 UMFPACK 和 SuiteSpase 包中的其他相關(guān)庫用于直接求解器.

在處理稀疏矩陣時(shí)，您可以打開診斷信息并查看使用 sparms:

sparms('spumoni',2)x = A;

此外，反斜杠運(yùn)算符也適用于 gpuArray，在在這種情況下，它依賴于 cuBLAS 和 MAGMA 在 GPU 上執(zhí)行.

它也適用于 分布式數(shù)組，適用于分布式計(jì)算環(huán)境(工作在一組計(jì)算機(jī)中分配，其中每個(gè)工作人員只有數(shù)組的一部分，可能整個(gè)矩陣不能一次全部存儲(chǔ)在內(nèi)存中).底層實(shí)現(xiàn)使用 ScaLAPACK.

如果你想自己實(shí)現(xiàn)所有這些，那是一個(gè)非常艱巨的任務(wù):)

I'm currently trying to develop a small matrix-oriented math library (I'm using Eigen 3 for matrix data structures and operations) and I wanted to implement some handy Matlab functions, such as the widely used backslash operator (which is equivalent to mldivide ) in order to compute the solution of linear systems (expressed in matrix form).

Is there any good detailed explanation on how this could be achieved ? (I've already implemented the Moore-Penrose pseudoinverse pinv function with a classical SVD decomposition, but I've read somewhere that A isn't always pinv(A)*b , at least Matalb doesn't simply do that)

Thanks

For x = A, the backslash operator encompasses a number of algorithms to handle different kinds of input matrices. So the matrix A is diagnosed and an execution path is selected according to its characteristics.

The following page describes in pseudo-code when A is a full 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

For non-square matrices, QR decomposition is used. For square triangular matrices, it performs a simple forward/backward substitution. For square symmetric positive-definite matrices, Cholesky decomposition is used. Otherwise LU decomposition is used for general square matrices.

Update: MathWorks has updated the algorithm section in the doc page of mldivide with some nice flow charts. See here and here (full and sparse cases).

All of these algorithms have corresponding methods in LAPACK, and in fact it's probably what MATLAB is doing (note that recent versions of MATLAB ship with the optimized Intel MKL implementation).

The reason for having different methods is that it tries to use the most specific algorithm to solve the system of equations that takes advantage of all the characteristics of the coefficient matrix (either because it would be faster or more numerically stable). So you could certainly use a general solver, but it wont be the most efficient.

In fact if you know what A is like beforehand, you could skip the extra testing process by calling linsolve and specifying the options directly.

if A is rectangular or singular, you could also use PINV to find a minimal norm least-squares solution (implemented using SVD decomposition):

x = pinv(A)*b

All of the above applies to dense matrices, sparse matrices are a whole different story. Usually iterative solvers are used in such cases. I believe MATLAB uses UMFPACK and other related libraries from the SuiteSpase package for direct solvers.

When working with sparse matrices, you can turn on diagnostic information and see the tests performed and algorithms chosen using spparms:

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

What's more, the backslash operator also works on gpuArray's, in which case it relies on cuBLAS and MAGMA to execute on the GPU.

It is also implemented for distributed arrays which works in a distributed computing environment (work divided among a cluster of computers where each worker has only part of the array, possibly where the entire matrix cannot be stored in memory all at once). The underlying implementation is using ScaLAPACK.

That's a pretty tall order if you want to implement all of that yourself :)

這篇關(guān)于如何實(shí)現(xiàn) Matlab 的 mldivide(又名反斜杠運(yùn)算符“")的文章就介紹到這了，希望我們推薦的答案對大家有所幫助，也希望大家多多支持html5模板網(wǎng)！