22.2 疎行列上の線形代数
Octave には、疎行列用の多態的ソルバーが含まれています。行列を因数分解するために使用される正確なソルバーは、疎行列自体の特性に依存します。一般に、行列タイプを決定するコストは、行列自体を因数分解するコストに比べて小さいですが、いずれの場合も、行列タイプは計算されるとキャッシュされるため、線形方程式で使用されるたびに再決定されることはありません。
線形方程式を解く方法の選択ツリーは
行列が対角行列の場合は直接解いて8に進みます。
行列が順列対角行列である場合、順列を考慮して直接解きます。Goto 8
行列が正方で、バンド化されており、バンド密度が指定された値より小さい場合はspparms ("bandden")続行し、そうでない場合は 4 に進みます。
行列が三角行列で、右側がスパースでない場合は続行し、そうでない場合は 3b に進みます。
行列がエルミート行列で、対角成分が正の実数である場合は、LAPACK xPTSV を使用してコレスキー分解を試みます。
上記が失敗した場合、または行列が正の実対角を持つエルミート行列でない場合は、 LAPACK xGTSV を使用してピボットによるガウス消去法を使用し、8 に進みます。
行列が正の実対角を持つエルミート行列の場合、LAPACK xPBTRF を使用してコレスキー分解を試みます。
上記が失敗した場合、または行列が正の実対角を持つエルミート行列でない場合は、 LAPACK xGBTRF を使用してピボットによるガウス消去法を使用し、8 に進みます。
行列が上三角または下三角の場合は、疎な前方または後方置換を実行し、8に進みます。
行列が列置換のある上三角行列または行置換のある下三角行列である場合、疎な前方置換または後方置換を実行し、8に進みます。
行列が正方エルミート行列で、対角要素が正の実数である場合、CHOLMODを使用してスパース コレスキー分解を試みます。
スパース コレスキー分解が失敗した場合、または行列が実数の正の対角を持つエルミート行列ではなく、行列が正方行列である場合は、UMFPACKを使用して因数分解し、解き、1 回の改良反復を実行します。
行列が正方でない場合、または以前のソルバーのいずれかが特異行列または特異行列に近い行列としてフラグを立てた場合は、 CXSPARSE 10を使用して最小ノルム解を見つけます。
バンド密度は、バンド内の非ゼロ値の数をバンド全体の値の総数で割ったものとして定義されます。 バンド行列ソルバーは、spparmsbanddenを使用して1 (つまりspparms ("bandden", 1)) に設定する ことで完全に無効にすることができます。
QR ソルバーは、問題を Dulmage-Mendelsohn 分解で因数分解し、過剰決定として扱えるブロック、複数の適切に決定されたブロック、および最終的な過剰決定ブロックに問題を分割します。強く接続されたノードのブロックを持つ行列の場合、LU 分解を多くのブロックに使用できるため、これは大きな利点です。また、単に NaNのベクトルを返すのではなく、過剰決定の問題に対するソリューションを見つける可能性が大幅に向上します。
上記のソルバーはすべて、条件数の推定値を計算できます。これを使用して、解の数値安定性の問題を検出し、最小ノルム解を強制的に使用できます。ただし、狭帯域、三角、または対角行列の場合、条件数の計算コストは大きく、実際には行列の因数分解のコストを超える可能性があります。したがって、これらのケースでは条件数は計算されず、Octave は特異行列を検出するためにより単純な手法を使用するか、帯域行列の場合は 基礎となるLAPACKコードを使用します。
ユーザーは、関数を使用して行列の型を強制することができますmatrix_type 。これにより、行列の型を検出するコストを回避できます。ただし、行列の型を誤って識別すると予期しない結果につながるため、matrix_type注意して使用する必要があります。
: nest = normest (A)
: nest = normest (A, tol)
: [nest, iter] = normest (…)
べき級数解析を使用して 行列Aの 2 ノルムを推定します。
これは通常、計算コストが 法外に高く、2 ノルムへの近似値が許容される大規模な行列に使用されます。 norm (A)
tolは 2 ノルムを計算する許容値です。デフォルトでは tolは 1e-6 です。
オプションの出力iter は、収束するのに必要な反復回数を返しますnormest。
See also: normest1, norm, cond, condest.
: nest = normest1 (A)
: nest = normest1 (A, t)
: nest = normest1 (A, t, x0)
: nest = normest1 (Afcn, t, x0, p1, p2, …)
: [nest, v] = normest1 (A, …)
: [nest, v, w] = normest1 (A, …)
: [nest, v, w, iter] = normest1 (A, …)
ブロックアルゴリズムを使用して 行列Aの 1 ノルムを推定します。
normest1は、ノルムの推定値のみが必要な大規模な疎行列に最適です。小規模から中規模の行列の場合は、 の使用を検討してください 。さらに、行列ベクトル積とが簡単に計算できる場合、 は線形演算子Aの 1 ノルムの推定に使用できます。この場合、行列Aの代わりに関数 が使用されます。この関数は次を返す必要があります。 norm (A, 1)normest1A * xA' * xAfcn (flag, x)
Aの次元n（フラグが"dim"
Aが実数演算子の場合はtrue 、フラグが"real"
フラグがA * x"notransp"
フラグがA' * x"transp"
典型的なケースは、によって定義されるAであり、次のように明示的に形成しなくても結果を計算できます。 b ^ mA * xb ^ m
y = x; for i = 1:m y = b * y; endfor
パラメータp1、p2、… は の引数です 。 Afcn (flag, x, p1, p2, …)
tのデフォルト値は 2 です。アルゴリズムには、サイズがn x nおよびn x tの行列間積が必要です。
初期行列x0 は単位 1 ノルムの列を持つ必要があります。デフォルトの初期行列x0には最初の列があり 、t > 1 の場合は、残りの列にはnで割ったランダム要素 、が含まれます。 ones (n, 1) / n-1 / n1 / n
出力では、nestは目的の推定値、vおよびwは、 =と なるベクトルです。iter に は、反復回数 (最大値は 5 にハードコードされています) と、アルゴリズムによって実行された 積 または合計数が含まれます。w = A * vnorm (w, 1)c * norm (v, 1)iter(1)iter(2)A * xA' * x
アルゴリズムに関する注意:normest1は評価中に乱数を使用します。したがって、一貫した結果が必要な場合は、"state"を呼び出す前に乱数ジェネレータの を固定する必要がありますnormest1。
参考文献:N. J. Higham and F. Tisseur, A block algorithm for matrix 1-norm estimation, with and application to 1-norm pseudospectra, SIAM J. Matrix Anal. Appl., pp. 1185–1201, Vol 21, No. 4, 2000.
See also: normest, norm, cond, condest.
: cest = condest (A)
: cest = condest (A, t)
: cest = condest (A, Ainvfcn)
: cest = condest (A, Ainvfcn, t)
: cest = condest (A, Ainvfcn, t, p1, p2, …)
: cest = condest (Afcn, Ainvfcn)
: cest = condest (Afcn, Ainvfcn, t)
: cest = condest (Afcn, Ainvfcn, t, p1, p2, …)
: [cest, v] = condest (…)
t検定ベクトルとランダム化 1 ノルム推定量 を使用して、 正方行列Aの 1 ノルム条件数を推定します。
オプションの入力tはテスト ベクトルの数を指定します (デフォルトは 5)。
入力は行列Aになります(このアルゴリズムは、特に大規模な疎行列に適しています)。あるいは、行列の動作を関数によって暗黙的に定義することもできます。暗黙的な定義を使用する場合は、condest次の関数が必要です。
Afcn (flag, x) which must return Aの次元n（フラグが"dim" Aが実数演算子の場合はtrue 、フラグが"real" フラグが「notransp」の 場合の結果A * x フラグが「transp」の 場合の結果A' * x
Ainvfcn (flag, x) which must return の次元n、フラグが inv (A)"dim" が実数演算子の場合はtrue 、フラグが inv (A)"real" フラグが「notransp」の 場合の結果inv (A) * x フラグが「transp」の 場合の結果inv (A)' * x
任意のパラメータp1、p2 、… は、 およびの追加引数です 。 Afcn (flag, x, p1, p2, …)Ainvfcn (flag, x, p1, p2, …)
主な出力は1ノルムの条件数推定値cestです。
オプションの 2 番目の出力vは近似ヌル ベクトルであり、次の式を満たします。 norm (A*v, 1) == norm (A, 1) * norm (v, 1) / cest
アルゴリズムに関する注意:condestはランダム化アルゴリズムを使用して 1 ノルムを近似します。したがって、一貫した結果が必要な場合は、 "state"を呼び出す前に乱数ジェネレータの を固定する必要があります condest。
References:
N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. SIMAX vol 21, no 4, pp 1185–1201. https://dx.doi.org/10.1137/S0895479899356080
N.J. Higham and F. Tisseur, A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra. https://citeseer.ist.psu.edu/223007.html
See also: cond, rcond, norm, normest1, normest.
: spparms ()
: vals = spparms ()
: [keys, vals] = spparms ()
: val = spparms (key)
: spparms (vals)
: spparms ("default")
: spparms ("tight")
: spparms (key, val)
スパース ソルバーと因数分解関数で使用されるパラメーターを照会または設定します。
上記の最初の 4 つの呼び出しは現在の設定に関する情報を取得し、その他の呼び出しは現在の設定を変更します。パラメータはキーと値のペアとして保存されます。値はすべて浮動小数点数で、キーは次の文字列のいずれかです。
‘spumoni’ ソルバーのデバッグ情報の印刷レベル（デフォルト 0）
‘ths_rel’ 互換性のために含まれています。使用されません。(デフォルト 1)
‘ths_abs’ 互換性のために含まれています。使用されません。(デフォルト 1)
‘exact_d’ 互換性のために含まれています。使用されません。(デフォルト 0)
‘supernd’ 互換性のために含まれています。使用されません。(デフォルト 3)
‘rreduce’ 互換性のために含まれています。使用されません。(デフォルト 3)
‘wh_frac’ 互換性のために含まれています。使用されません。(デフォルト 0.5)
‘autommd’ LU/QR および '\' および '/' 演算子がスパース性を保持する mmd 関数を自動的に使用するかどうかのフラグ (デフォルト 1)
‘autoamd’ LU および '\' および '/' 演算子がスパース性を保持する amd 関数を自動的に使用するかどうかのフラグ (デフォルト 1)
‘piv_tol’ UMFPACKソルバーのピボット許容値（デフォルト0.1）
‘sym_tol’ UMFPACK対称ソルバーのピボット許容値（デフォルト0.001）
‘bandden’ LAPACKバンド ソルバーで処理される前のバンド マトリックス内の非ゼロ要素の密度(デフォルト 0.5)
‘umfpack’ LU、'\'、'/' 演算に UMFPACKまたは mmd ソルバーを使用するかどうかを示すフラグ(デフォルト 1)
個々のキーの値は で設定できます 。 デフォルト値は、特殊キーワード で復元できます 。 特殊キーワードを使用すると、実行時間が長くなる可能性はありますが、よりスパースなソリューションを試行するように mmd ソルバーを設定できます。 spparms (key, val)"default""tight"
See also: chol, colamd, lu, qr, symamd.
: p = sprank (S)
疎行列Sの構造ランクを計算します。
この計算では、三角形をブロックするためのダルメージ-メンデルゾーン置換に基づいて、行列の構造のみが使用されることに注意してください。そのため、行列Sの数値ランクは で制限されます 。浮動小数点エラーは無視します 。 sprank (S) >= rank (S)sprank (S) == rank (S)
See also: dmperm.
: [count, h, parent, post, R] = symbfact (S)
: […] = symbfact (S, typ)
: […] = symbfact (S, typ, mode)
疎行列Sの記号分解解析を実行します。
入力変数は
S Sは実数または複素数のスパース行列です。
typ
因数分解の型であり、次のいずれかになります。
"sym" (default)
Sを因数分解します。S が対称であると仮定し、行列の上三角部分を使用します。
"col"
因数分解します。 S' * S
"row"
因数分解します。 S * S'
"lo"
因数分解しますS'。Sが対称であると仮定し、行列の下三角部分を使用します。
mode
modeが指定されていない場合は、 Rのコレスキー分解を返します 。mode が または の場合は、"lower"下三角因子である"L"共役転置を返します。共役転置バージョンはより高速でメモリ使用量も少なくなりますが、他のすべての出力 ( count、h、 parent、post )R'に対しては同じ値を返します。
出力変数は次のとおりです。
count typによって決定されるコレスキー分解の行数 。 を使用して真の分解を実行する計算上の困難さはcholです。 sum (count .^ 2)
h 除去ツリーの高さ。
parent 除去ツリー自体。
post typによって決定されるコレスキー分解の構造を持つ疎ブール行列。
See also: chol, etree, treelayout.
非正方行列の場合、ユーザーはspaugment 関数を利用して線形方程式の最小二乗解を見つけることもできます。
: s = spaugment (A, c)
Aの拡張行列を作成します。
これは次のように与えられる。
[c * eye(m, m), A;
A', zeros(n, n)]
これはA \ bの最小二乗解と関係があり 、
s * [ r / c; x] = [ b, zeros(n, columns(b)) ]
ここでrは残差誤差である
r = b - A * x
行列sは対称不定値なので、 で因数分解することができ lu、したがって因数分解を必要とせずに最小ノルム解を求めることができるqr。残差誤差は 劣決定問題では となり、例は zeros (m, m)
m = 11; n = 10; mn = max (m, n);
A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],
[-1, 0, 1], m, n);
x0 = A \ ones (m,1);
s = spaugment (A);
[L, U, P, Q] = lu (s);
x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)])));
x1 = x1(end - n + 1 : end);
過剰決定問題の解を見つけるには残差誤差rの推定値が必要なので、関数を使用して最小ノルム解を定式化するのはより複雑ですspaugment。
一般に、左除算演算子はspaugment関数を使用するよりも安定しており、高速です。
See also: mldivide.
最後に、この関数はeigs、選択基準に基づいて限られた数の固有値と固有ベクトルを計算するために使用でき、同様にsvds限られた数の特異値と特異ベクトルを計算するために使用できます。
: d = eigs (A)
: d = eigs (A, k)
: d = eigs (A, k, sigma)
: d = eigs (A, k, sigma, opts)
: d = eigs (A, B)
: d = eigs (A, B, k)
: d = eigs (A, B, k, sigma)
: d = eigs (A, B, k, sigma, opts)
: d = eigs (Af, n)
: d = eigs (Af, n, k)
: d = eigs (Af, n, k, sigma)
: d = eigs (Af, n, k, sigma, opts)
: d = eigs (Af, n, B)
: d = eigs (Af, n, B, k)
: d = eigs (Af, n, B, k, sigma)
: d = eigs (Af, n, B, k, sigma, opts)
: [V, D] = eigs (…)
: [V, D, flag] = eigs (…)
選択基準に基づいて、限られた数の固有値と固有ベクトルを計算します。
デフォルトでは、eigs対応する固有ベクトルである方程式を解きます。正定値行列 Bが与えられた場合、eigs一般固有値方程式を解きます。
入力A は、 n行n列の次元の正方行列です。通常、Aも大きく、疎です。
一般化固有値問題の入力B は、 Aと同じサイズ( n行n列)の正方行列です。通常、Bも大きく、疎です。
計算する固有値と固有ベクトルの数はkで指定され 、デフォルトは 6 です。
引数sigma は、返される固有値を決定します。sigma はスカラーまたは文字列のいずれかになります。sigma がスカラーの場合、 sigmaに最も近いk 個の固有値が返されます。sigmaが文字列の場合 、次のいずれかの値である必要があります。
"lm" 最大の大きさ（デフォルト）
"sm" 最小のマグニチュード
"la" 最大代数（実対称問題にのみ有効）
"sa" 最小代数（実対称問題にのみ有効）
"be" 両端。kが奇数の場合は上限から 1 つ増えます(実対称問題にのみ有効)
"lr" 最大の実数部（複雑な問題または非対称の問題にのみ有効）
"sr" 最小の実数部 (複雑な問題または非対称の問題にのみ有効)
"li" 最大虚数部（複雑な問題または非対称の問題にのみ有効）
"si" 最小の虚数部（複雑な問題または非対称の問題にのみ有効）
optsが指定されている場合は、使用すべきオプションを定義する構造体です 。opts構造eigs体のフィールドは次のとおりです。
issym Afが指定されている場合、このフラグ (true/false) は関数Af が対称問題を定義するかどうかを決定します。行列 Aが指定されている場合は無視されます。デフォルトは false です。
isreal Afが指定されている場合、このフラグ (true/false) は関数Afが実際の問題を定義するかどうかを決定します。行列A が指定されている場合は無視されます。デフォルトは true です。
tol 必要な収束許容値を定義します。これは として計算されます tol * norm (A)。デフォルトは ですeps。
maxit 反復の最大回数。デフォルトは 300 です。
p 使用する Lanczos 基底ベクトルの数。ベクトルの数が多いほど収束は速くなりますが、メモリの使用量は多くなります。 の最適値はp問題に依存し、 k + 1～nの範囲でなければなりません。デフォルト値は です。 2 * k
v0 アルゴリズムの開始ベクトル。最終ベクトルに近い初期ベクトルは、収束を高速化します。デフォルトでは、ARPACKは開始ベクトルをランダムに生成します。指定する場合は、n行 1 列のベクトルでv0ある必要があります。 n = rows (A)
disp 診断出力のレベル (0|1|2)。disp0 の場合、診断は無効になります。デフォルト値は 0 です。
cholB 一般化固有値問題が計算されている場合、このフラグ (true/false) は、B入力が行列 B を表す か、単に行列B を表すかを指定します。デフォルトは false です。 chol (B)
permB が真の場合 のBのコレスキー分解の順列ベクトルcholB。これは によって得られます 。デフォルトは です 。 [R, ~, permB] = chol (B, "vector")1:n
A をAfで示される関数で表すこともできます。 Af の後には、 Afが受け入れるベクトル引数の長さを定義するスカラー引数nを続ける必要があります。 Af は、関数ハンドル、インライン関数、または文字列にすることができます。Afが文字列の場合、使用する関数の名前が保持されます。
Af は、 Afy = Af (x)の必要な戻り値がsigmaの値によって決定される形式の関数です。可能な形式は4つあります。
A * x sigmaが指定されていないか、または「sm」以外の文字列である 場合
A \ x sigmaが 0 または "sm" の 場合
(A - sigma * I) \ x sigmaが 0 以外のスカラーの場合、 AIと同じサイズの単位行列になります
(A - sigma * B) \ x 一般的な固有値問題の場合
戻り引数とその形式は、要求された戻り引数の数によって異なります。戻り引数が 1 つの場合、見つかったk 個の固有値を含む長さ kの列ベクトルd が返されます。戻り引数が 2 つの場合、 Vはn行k列の行列で、その列は返された固有値に対応するk 個の固有ベクトルです。固有値自体は、対角要素が固有値である k行k列の行列の形式でDに返されます。
3 番目の戻り引数flag は収束の状態を返します。flag が 0 の場合、すべての固有値は収束しています。その他の値は収束に失敗したことを示します。
プログラミング ノート: n < 500 の小さな問題の場合は、 の使用を検討してください 。 eig (full (A))
ARPACK が収束しない場合は、Lanczos ベクトルの数を増やす ( opt .p )、反復回数を増やす ( opt .maxiter )、または許容値を減らす ( opt .tol ) ことを検討してください。
参考: この関数は、R. Lehoucq、K. Maschhoff、D. Sorensen、および C. Yang によって作成されたARPACKパッケージに基づいています。詳細については、http://www.caam.rice.edu/software/ARPACK/を参照してください。
See also: eig, svds.
: s = svds (A)
: s = svds (A, k)
: s = svds (A, k, sigma)
: s = svds (A, k, sigma, opts)
: [u, s, v] = svds (…)
: [u, s, v, flag] = svds (…)
行列Aの特異値をいくつか見つけます。
特異値は次のように計算される
[m, n] = size (A);
s = eigs ([sparse(m, m), A;
A', sparse(n, n)])
によって返される固有値は、Aeigsの特異値に対応します。計算する特異値の数はkによって指定され 、デフォルトは 6 です。
引数sigma は、どの特異値を検索するかを指定します。 sigmaが文字列'L'（デフォルト）の場合、Aの最大の特異値が検索されます。 それ以外の場合、sigma は実数スカラーでなければならず、 sigmaに最も近い特異値が検索されます。 結果として、 は最小の特異値を検索します。 sigmasigma = 0の値が比較的小さい場合、要求された数の特異値が見つからない可能性があることに注意してください。 その場合は、 sigmaを増やす必要があります。
opts はsvds、 に渡されるオプションを定義する構造体ですeigs。この構造体の可能なフィールドは に記載されています eigs。デフォルトでは、svdsは次の 3 つのフィールドを設定します。
tol 特異値に必要な収束許容値。デフォルト値は 1e-10 です。 eigsが渡されますtol / sqrt (2)
maxit 反復の最大回数。デフォルトは 300 です。
disp 診断出力のレベル (0|1|2)。disp0 の場合、診断は無効になります。デフォルト値は 0 です。
複数の出力が要求された場合、Asvdsの特異値分解の近似値を返します。
A_approx = u*s*v'
ここで、A _approx はサイズがAだがランクがk の行列です。
アルゴリズムが正常に収束した場合はフラグは0を返し、そうでない場合は1を返す。収束のテストは
norm (A*v - u*s, 1) <= tol * norm (A, 1)
svds大きな疎行列から少数の特異値のみを見つけるのに最適です。それ以外の場合は、の方が効率的である可能性があります。 svd (full (A))
See also: svd, eigs.
Footnotes
(10)
CHOLMOD 、UMFPACK、CXSPARSE パッケージは Tim Davis によって作成され、 http://faculty.cse.tamu.edu/davis/suitesparse.htmlから入手できます。