20

私は現在、小さな行列指向の数学ライブラリを開発しようとしています (行列のデータ構造と操作にEigen 3を使用しています)。広く使用されているバックスラッシュ演算子などのいくつかの便利な Matlab 関数を実装したいと考えていました (これはmldivide ) を使用して、線形システム (行列形式で表現) の解を計算します。

これをどのように達成できるかについて、詳細な説明はありますか? (ムーア・ペンローズの疑似逆pinvA\b関数を古典的な SVD 分解で既に実装していますが、常にそうであるとは限らない場所を読んだことがあります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 分解が使用されます。正三角形行列の場合、単純な前方/後方置換を実行します。正定正定正定行列の場合、コレスキー分解が使用されます。それ以外の場合、一般的な正方行列にはLU 分解が使用されます。

更新: MathWorks はのドキュメント ページのアルゴリズム セクションmldivideをいくつかの優れたフローチャートで更新しました。hereとhereを参照してください（完全な場合と疎な場合）。

mldivide_full

これらのアルゴリズムはすべてLAPACKに対応するメソッドがあり、実際、おそらく MATLAB が行っていることです (MATLAB の最近のバージョンには、最適化されたIntel MKL実装が付属していることに注意してください)。

異なる方法を使用する理由は、最も特殊なアルゴリズムを使用して、係数行列のすべての特性を利用する連立方程式を解こうとするためです (その方が高速であるか、数値的に安定しているからです)。したがって、確かに一般的なソルバーを使用できますが、最も効率的ではありません。

実際、事前にどのようなものかがわかっている場合は、オプションを直接A呼び出して指定することで、追加のテスト プロセスをスキップできます。linsolve

Aが長方形または特異な場合、PINV を使用して最小ノルム最小二乗解を見つけることもできます ( SVD 分解を使用して実装):

x = pinv(A)*b

上記はすべて密行列に適用されますが、疎行列はまったく別の話です。通常、このような場合には反復ソルバーが使用されます。MATLAB はUMFPACKと SuiteSpase パッケージのその他の関連ライブラリを直接ソルバーに使用していると思います。

疎行列を操作する場合、診断情報をオンにして、実行されたテストと選択されたアルゴリズムを次のように表示できspparmsます。

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

さらに、バックスラッシュ演算子はgpuArrayでも機能します。この場合、GPU での実行はcuBLASとMAGMAに依存します。

また、分散コンピューティング環境で動作する分散配列にも実装されます (各ワーカーが配列の一部しか持たず、おそらく行列全体を一度にメモリに格納できないコンピューターのクラスター間で分割された作業)。基礎となる実装はScaLAPACKを使用しています。

そのすべてを自分で実装したい場合、それはかなり難しい注文です:)

于 2013-08-31T23:30:22.287 に答える