math - GSL / BLAS：行列と逆行列を乗算します

Question

GNUGSLを使用して行列計算を行っています。行列Bに行列Aの逆行列を乗算しようとしています。

ここで、GSLのBLAS部分にこれを行う機能があることに気付きましたが、これはAが三角形の場合に限られます。これには特別な理由がありますか？また、この計算を行うための最速の方法は何でしょうか？LU分解を使用してAを反転する必要がありますか、それともより良い方法がありますか？

FWIW、Aの形式はP'G Pです。ここで、Pは正規行列、P'はその逆行列、Gは対角行列です。

本当にありがとう ：）

Answer 1

BLAS が正方行列の逆関数を持たない理由について、Adrien は正しいと思います。逆行列の計算を最適化するために使用している行列によって異なります。

一般に、任意の正方行列に有効な LU 分解を使用できます。つまり、次のようなものです。

gsl_linalg_LU_decomp(A, p, signum);
gsl_linalg_LU_invert(A, p, invA);

ここで、A は逆行列を求める正方行列、p はgsl_permutationオブジェクト (順列行列がエンコードされる順列オブジェクト)、signum は順列の符号、invA は A の逆行列です。

正規で対角線であると述べてA = P' G Pいるので、おそらく正規行列です。私はしばらくそれらを使用していませんが、それらの因数分解定理がGSL リファレンス マニュアルの第 14 章または線形代数の本で見つかるはずです。PGA

ほんの一例ですが、対称正定行列があり、その逆行列を見つけたい場合は、その種類の行列用に最適化されたコレスキー分解を使用できます。次に、GSL で関数gsl_linalg_cholesky_decomp()とgsl_linalg_cholesky_invert()関数を使用して、GSL を効率的にすることができます。

それが役立つことを願っています！

Answer 2

手短に：

三角行列のみがサポートされているという事実は、この操作が三角行列に対して実行するのが非常に簡単であるため (逆置換と呼ばれます)、BLAS は低レベル関数のルーチンのみを提供するためです。より高いレベルの関数は通常、ブロック単位の操作に BLAS を内部的に使用する LAPACK にあります。

あなたが扱っている特定のケースでは: A:=P'*G*P、Pが通常の行列の場合、 の逆をA書くことができますinv(A) := P'*inv(G)*P。次に、2 つのオプションがあります。

1. ビルドinv(A)して掛けBます。
2. B を のさまざまな係数で順次乗算します。(左側) をinv(A)乗算し、結果の各行を の対角要素の逆数で再スケーリングし、(左側) を再度乗算します。BPGP'

特定のケースに応じて、方法を選択できます。とにかく、倍精度を仮定するdgemmと、(行列-行列乗算、Lvl3 BLAS) とdscal(行のスケーリング、Lvl 1 BLAS) が必要になります。

お役に立てれば。

A.