22.3 疎行列に適用される反復手法
前のセクションで説明した左除算演算子\と右除算演算子は、直接ソルバを使用して、または の形式の線形方程式を解きます。Octave には、反復手法を使用して疎線形方程式を解くための関数もいくつか含まれています。 /x = A \ bx = b / A
: x = pcg (A, b, tol, maxit, m1, m2, x0, …)
: x = pcg (A, b, tol, maxit, M, [], x0, …)
: [x, flag, relres, iter, resvec, eigest] = pcg (A, b, …)
前処理付き共役勾配反復法を使用して 線形方程式系を解きます 。A * x = b
入力引数は次のとおりです。
Aは線形システムの行列であり、正方行列でなければなりません。 A はAfcn、 のような行列、関数ハンドル、またはインライン関数として渡すことができます 。 x0Afcn(x) = A * xの後に、追加のパラメータ Afcnを渡すこともできます。
A はエルミートかつ正定値 (HPD) である必要があります。A が 正pcg定値ではないことが検出された場合、警告が印刷され、フラグ出力が設定されます。
bは右側のベクトルです。
tol は残差誤差に必要な相対許容値です 。 ≤ の場合、反復は停止します 。tolが省略されるか空の場合、許容値 1e-6 が使用されます。 b - A * xnorm (b - A * x)tol * norm (b)
maxit は許容される最大反復回数です。maxit が省略されるか空の場合、 値 20 が使用されます。
mは HPD 前処理行列です。 が HPD であるような 任意の分解に対して 、共役勾配法は形式的には線型システム に適用されます (分割前処理)。実際には、共役勾配法の各反復で、行列mを持つ線型システムは で解かれます。特定の因数分解 が使用できる場合 (たとえば、 の不完全コレスキー分解)、2 つの行列 m1とm2を渡して、相対的な線型システムを 演算子で解くことができます。前処理行列を適切に選択すると、この方法の全体的なパフォーマンスが大幅に向上する可能性があることに注意してください。行列m1と m2の代わりに、ユーザーは、ベクトルにm1とm2の逆関数を適用した結果を返す 2 つの関数を渡すことができます。 m1が省略されるか空の場合、前処理は適用されません。 mの因数分解が使用できない場合は、m2 を 省略するか [] のままにして、入力変数m1 を使用して前処理行列mを渡すことができます。 m = p1 * p2inv (p1) * A * inv (p2)inv (p1) * A * inv (p2) * y = inv (p1) * bx = inv (p2) * ymldividem = m1 * m2mldivide[]
x0は初期推定値です。x0 が省略または空の場合、関数はデフォルトでx0 をゼロ ベクトルに設定します。
x0に続く引数はパラメータとして扱われ、に与えられた関数 ( Aまたはm1または m2 )のいずれかに適切な方法で渡されますpcg。詳細については、以下の例を参照してください。
出力引数は次のとおりです。
x はの解の計算された近似値です 。アルゴリズムが収束しなかった場合、xは残差が最小となる反復です。 A * x = b
収束に関する フラグレポート:
0: アルゴリズムは規定の許容範囲内に収束しました。
1: アルゴリズムが収束せず、反復の最大回数に達しました。
2: 前処理行列は特異である。
3: アルゴリズムが停滞しました。つまり、現在の反復xと前の反復の差の絶対値が 未満です。 eps * norm (x,2)
4: アルゴリズムは、入力（前処理済み）行列が HPD ではないことを検出します。
relres は、ユークリッドノルムで測定された最終残差と初期値の比です。
iter はxが計算された反復回数を示します。出力x は最小残差解に対応するため、メソッドが実行した反復回数の合計は で与えられますlength(resvec) - 1。
resvec は、この方法の収束履歴を表します。 は残差のユークリッドノルムで、 は ( i -1) 回目の反復後の前処理された残差ノルムです。 前処理された残差ノルムはとして 定義されます 。mの説明も参照してください。eigestが不要な場合は、 のみが 返されます。 resvec (i, 1)resvec (i, 2)i = 1, 2, …, iter+1r' * (m \ r)r = b - A * xresvec (:, 1)
eigest は、前処理された行列 の 最小eigest(1) および最大の固有値の推定値を返します。特に、前処理が使用されない場合は、 Aの極限固有値の推定値が返されます。 は過大評価であり、 は過小評価であるため、 は の下限値になりますが 、それでも限界では理論的には条件数の実際の値に等しくなるはずです。 eigest(2)P = m \ Aeigest(1)eigest(2)eigest(2) / eigest(1)cond (P, 2)
三角行列に関する簡単な問題を考えてみましょう
n = 10; A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n)); b = A * ones (n, 1); M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)' M2 = M1'; M = M1 * M2; Afcn = @(x) A * x; Mfcn = @(x) M \ x; M1fcn = @(x) M1 \ x; M2fcn = @(x) M2 \ x;
EXAMPLE 1: 最も簡単なpcg使用法 x = pcg (A, b)
EXAMPLE 2:pcg計算する関数 A * x x = pcg (Afcn, b)
EXAMPLE 3: pcg前処理行列M x = pcg (A, b, 1e-06, 100, M)
EXAMPLE 4: pcg関数を前処理として用いる x = pcg (Afcn, b, 1e-6, 100, Mfcn)
EXAMPLE 5: pcg前処理行列M1 とM2 x = pcg (A, b, 1e-6, 100, M1, M2)
EXAMPLE 6: pcg関数を前処理として用いる x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn)
EXAMPLE 7: pcg引数を必要とする関数を入力として
function y = Ap (A, x, p) # compute A^p * x
y = x;
for i = 1:p
y = A * y;
endfor
endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = pcg (Apfcn, b, [], [], [], [], [], 2);
EXAMPLE 8:pcg分割前処理を使用する ことを示す明示的な例 M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed M2 = M1'; M = M1 * M2; ## reference solution computed by pcg after two iterations [x_ref, fl] = pcg (A, b, [], 2, M) ## split preconditioning [y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2) x = M2 \ y # compare x and x_ref
References:
C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995. (the base PCG algorithm)
Y. Saad, Iterative Methods for Sparse Linear Systems, PWS 1996. (condition number estimate from PCG) Revised version of this book is available online at https://www-users.cs.umn.edu/~saad/books.html
See also: sparse, pcr, gmres, bicg, bicgstab, cgs.
: x = pcr (A, b, tol, maxit, m, x0, …)
: [x, flag, relres, iter, resvec] = pcr (…)
前処理付き共役残差反復法を使用して 線形方程式系を解きます。A * x = b
入力引数は
A は、正方行列（できれば疎行列）か、関数ハンドル、インライン関数、または を計算する関数の名前を含む文字列のいずれかです。原則として、A は対称かつ非特異である必要があります。Aが数値的に特異であると判断された場合、警告メッセージが表示され、フラグ出力パラメータが設定されます。 A * xpcr
bは右側のベクトルです。
tol は残差誤差に必要な相対許容値です 。 の場合、反復は停止します 。tolが空または省略されている場合、関数は デフォルトで設定します。 b - A * xnorm (b - A * x) <= tol * norm (b - A * x0)tol = 1e-6
maxit は許容される反復の最大回数です。maxit[]に が指定されている場合、または引数が少ない場合は、デフォルト値の 20 が使用されます。 pcr
mは（左）前処理行列であるため、この反復は（理論的には） を用いて を解くことと等価です。前処理行列を適切に選択すると、この方法の全体的なパフォーマンスが大幅に向上することに注意してください。行列mの代わりに、 mの逆関数をベクトルに適用した結果を返す関数を渡すこともできます（通常、これは前処理行列を使用する好ましい方法です）。 mに が指定されている場合、またはmが省略されている場合、前処理は適用されません。 pcr P * x = m \ bP = m \ A[]
x0は初期推定値です。x0 が空または省略されている場合、関数はデフォルトでx0 をゼロ ベクトルに設定します。
x0に続く引数はパラメータとして扱われ、に渡される関数 ( Aまたはm )のいずれかに適切な方法で渡されますpcr。詳細については、以下の例を参照してください。
出力引数は
x はの解の計算された近似値です 。 A * x = b
フラグは収束を報告します。 flag = 0解が収束し、tolで指定された許容基準が満たされたことを意味します。 反復回数の上限に達したことを意味します。 内訳flag = 1を報告します 。詳細については[1]を参照してください。 flag = 3pcr
relres は、ユークリッドノルムで測定された最終残差と初期値の比です。
iter は実際に実行された反復回数です。
resvecは方法の収束履歴を記述し、 ( i -1)回目の反復resvec (i)後の残差のユークリッドノルムが含まれます。 i = 1,2, …, iter+1
対角行列に関する簡単な問題を考えてみましょう（Aのスパース性を利用します）
n = 10; A = sparse (diag (1:n)); b = rand (N, 1);
EXAMPLE 1: 最も簡単な使用法pcr x = pcr (A, b)
EXAMPLE 2: pcrを計算する関数を使用します 。 A * x
function y = apply_a (x)
y = [1:10]' .* x;
endfunction
x = pcr ("apply_a", b)
EXAMPLE 3: 完全な診断を伴う前処理反復。前処理子(元の行列 Aでさえ自明なので、かなり奇妙である)は関数として定義される。
function y = apply_m (x)
k = floor (length (x) - 2);
y = x;
y(1:k) = x(1:k) ./ [1:k]';
endfunction
[x, flag, relres, iter, resvec] = ...
pcr (A, b, [], [], "apply_m")
semilogy ([1:iter+1], resvec);
EXAMPLE 4: 最後に、パラメータkに依存する前処理子。
function y = apply_m (x, varargin)
k = varargin{1};
y = x;
y(1:k) = x(1:k) ./ [1:k]';
endfunction
[x, flag, relres, iter, resvec] = ...
pcr (A, b, [], [], "apply_m"', [], 3)
Reference:
W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, section 9.5.4; Springer, 1994
See also: sparse, pcg.
反復ソルバーが解に収束する速度は、前処理行列Mを使用することで加速できます。この場合、代わりに線形方程式が解かれます。一般的な前処理行列は、元の行列の部分因数分解です。 M^-1 * x = M^-1 * A \ b
: L = ichol (A)
: L = ichol (A, opts)
疎正方行列 Aの不完全コレスキー分解を計算します。
デフォルトでは、Aicholの下三角のみを使用し、 A を 近似するような下三角因子Lを生成します。 L*L'
このルーチンによって与えられた係数は、PCG (前処理付き共役勾配法) などの反復法によって解かれる線形方程式系の前処理子として役立つ場合があります。
分解は、構造体 optsでオプションを渡すことによって変更できます。オプション名は構造体のフィールドであり、設定はフィールドの値です。名前と指定子は大文字と小文字が区別されます。
type 因数分解の種類。
"nofill" (default) 埋め込みのない不完全なコレスキー分解（IC(0)）.
"ict" 閾値低下を伴う不完全コレスキー分解 (ICT).
diagcomp Aの代わりに不完全コレスキー分解に使用する 非負のスカラーアルファ。これは、 Aが正定値でない場合に便利です。デフォルト値は 0 です。 A + alpha * diag (diag (A))
droptol ICT を実行する場合の因数分解のドロップ許容値を指定する非負のスカラー。デフォルト値は 0 で、完全なコレスキー分解を生成します。 Lの非対角要素は、次の場合を除き0に設定されます。 abs (L(i,j)) >= droptol * norm (A(j:end, j), 1).
michol 修正不完全コレスキー分解: "off" (default) 行と列の合計は必ずしも保持されるわけではありません。 "on" Lの対角線は、因数分解中に要素が削除された場合でも行 (および列) の合計が保存されるように変更されます。保存される関係は次のとおりです。ここで、e は 1 のベクトルです。 A * e = L * L' * e
shape "lower" (default) Aの下三角のみを使用し、 Aに近似するような下三角因子 Lを返します。 L*L' "upper" Aの上三角のみを使用し、 A に近似するような上三角因子 Uを返します。 U'*U
例
次の問題は、完全なコレスキー分解と不完全なコレスキー分解を使用して、サンプル対称正定値行列を因数分解する方法を示しています。
A = [ 0.37, -0.05, -0.05, -0.07;
-0.05, 0.116, 0.0, -0.05;
-0.05, 0.0, 0.116, -0.05;
-0.07, -0.05, -0.05, 0.202];
A = sparse (A);
nnz (tril (A))
ans = 9
L = chol (A, "lower");
nnz (L)
ans = 10
norm (A - L * L', "fro") / norm (A, "fro")
ans = 1.1993e-16
opts.type = "nofill";
L = ichol (A, opts);
nnz (L)
ans = 9
norm (A - L * L', "fro") / norm (A, "fro")
ans = 0.019736
分解のもう 1 つの例は、単位正方形上の境界値問題を解くために使用される有限差分行列です。
nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
[-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
[-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans = 239400
opts.type = "nofill";
L = ichol (A, opts);
nnz (tril (A))
ans = 239400
norm (A - L * L', "fro") / norm (A, "fro")
ans = 0.062327
References for implemented algorithms:
[1] Y. Saad. "Preconditioning Techniques." Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996.
[2] M. Jones, P. Plassmann: An Improved Incomplete Cholesky Factorization, 1992.
See also: chol, ilu, pcg.
: LUA = ilu (A)
: LUA = ilu (A, opts)
: [L, U] = ilu (…)
: [L, U, P] = ilu (…)
疎正方行列Aの不完全 LU 分解を計算します。
iluは、 を近似する単位下三角行列L、上三角行列U、およびオプションで置換行列Pを返します。 L*UP*A
このルーチンによって与えられた係数は、BICG (BiConjugate Gradients) や GMRES (Generalized Minimum Residual Method) などの反復法によって解かれる線形方程式系の前提条件として役立つ場合があります。
分解は、構造体 optsでオプションを渡すことによって変更できます。オプション名は構造体のフィールドであり、設定はフィールドの値です。名前と指定子は大文字と小文字が区別されます。
type 因数分解の種類 "nofill" (default) 埋め込みなしのILU分解（ILU(0)）. サポートされている追加オプション: milu. "crout" ILU 因数分解 (ILUC) の Crout バージョン. サポートされている追加オプション: milu、droptol "ilutp" 因数分解のドロップ許容値を指定する非負のスカラー。デフォルト値は 0 で、完全な LU 因数分解を生成します。
droptol 因数分解のドロップ許容値を指定する非負のスカラー。デフォルト値は 0 で、完全な LU 因数分解を生成します。
Uの非対角要素は、次の場合を除き0に設定されます。 abs (U(i,j)) >= droptol * norm (A(:,j)).
Lの非対角要素は、次の場合を除き0に設定されます abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j).
milu 修正不完全LU分解: "row" 行和を修正した不完全 LU 分解。この分解では行和が保存されます: 、ここで e は 1 のベクトルです。 A * e = L * U * e "col" 列和を修正した不完全 LU 分解。この分解では列の合計が保持されます 。 e' * A = e' * L * U "off" (default) 行と列の合計は必ずしも保持されるわけではありません。
udiag true の場合、上三角因子の対角線上のゼロは、ローカルドロップ許容値に置き換えられます 。デフォルトは false です。 droptol * norm (A(:,j))/U(j,j)
thresh 因数分解のピボットしきい値。範囲は 0 (対角ピボット) から 1 (デフォルト) までで、列の最大絶対値のエントリがピボットとして選択されます。
が 1 つの出力のみで呼び出された場合ilu、返される行列は です 。ここで、Lは単位下三角行列、Uは上三角行列です。 L + U - speye (size (A))
2 つの出力で、ilu単位下三角行列L と上三角行列Uを返します。opts .type == の場合、係数の 1 つはopts"ilutp" .miluの値に基づいて並べ替えられます 。opts .milu ==の場合、U は列が並べ替えられた上三角係数です。それ以外の場合、Lは行が並べ替えられた単位下三角係数です。 "row"
名前付き出力が 3 つあり、opts .milu != の場合"row"、 LとU がの不完全因数となるPが返されます。opts .milu ==の場合、LとU がの不完全因数 となるP が返されます。 P*A"row"A*P
例
A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840
nnz (lu (A))
ans = 126478
nnz (ilu (A, opts))
ans = 7840
これは、Aには 7,840 個の非ゼロ要素があり、完全な LU 分解には 126,478 個の非ゼロ要素があり、0 レベルの埋め込みによる不完全な LU 分解にはAと同じ数の 7,840 個の非ゼロ要素があることを示しています。出典: https://www.mathworks.com/help/matlab/ref/ilu.html
A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)
この例では、大きな条件数を持つランダム FEM マトリックスの前提条件として ILU を使用します。LとUがないと、BICG は収束しません。
See also: lu, ichol, bicg, gmres.