pcg

x = pcg(A,b) は、前処理付き共役勾配法を使用して、線形方程式系 A*x = b を x について解きます。試行が正常に完了した場合、pcg は収束を確認するメッセージを表示します。pcg が最大反復回数に達しても収束しないか、何らかの理由で実行を停止した場合は、相対残差 norm(b-A*x)/norm(b) と停止時の反復回数を含む診断メッセージが表示されます。

x = pcg(A,b,tol) は、このメソッドの許容誤差を指定します。既定の許容誤差は 1e-6 です。

x = pcg(A,b,tol,maxit) は、使用する最大反復回数を指定します。pcg は、maxit 以内の反復で収束しない場合、診断メッセージを表示します。

x = pcg(A,b,tol,maxit,M) は前処理行列 M を指定し、y について系 $H^{- 1} A H^{- T} y = H^{- 1} b$ を実質的に解くことにより x を計算します。ここで、 $y = H^{T} x$ および $H = M^{1 / 2} = {(M_{1} M_{2})}^{1 / 2}$ です。アルゴリズムは H を明示的に形成しません。前処理行列を使用すると、問題の数値的なプロパティと計算の効率を向上させることができます。

x = pcg(A,b,tol,maxit,M1,M2) は M = M1*M2 となるような前処理行列 M の因子を指定します。

x = pcg(A,b,tol,maxit,M1,M2,x0) は解のベクトル x の初期推定を指定します。既定値はゼロのベクトルです。

[x,flag] = pcg(___) は、アルゴリズムが正常に収束したかどうかを示すフラグを返します。flag = 0 の場合、収束は正常に実行されています。この出力構文は、前記のすべての入力引数の組み合わせで使用できます。flag 出力を指定する場合、pcg は診断メッセージを表示しません。

[x,flag,relres] = pcg(___) は、相対残差 norm(b-A*x)/norm(b) も返します。flag が 0 の場合は、relres <= tol になります。

[x,flag,relres,iter] = pcg(___) は、x が計算されたときの反復回数 iter も返します。

[x,flag,relres,iter,resvec] = pcg(___) は、各反復で、最初の残差 norm(b-A*x0) を含む残差ノルムのベクトルも返します。

例

線形システムの反復法

既定の設定で pcg を使用して正方線形システムを解いてから、解法プロセスで使用される許容誤差と反復回数を調整します。

乱数対称スパース行列 A を作成します。また、真の解 $x$ が 1 のベクトルになるように、 $Ax = b$ の右辺の A の行の和のベクトル b を作成します。

rng default
A = sprand(400,400,.5);
A = A'*A;
b = sum(A,2);

pcg を使用して $Ax = b$ を解きます。出力の表示には相対残差誤差 $\frac{‖ b - Ax ‖}{‖ b ‖}$ の値が含まれます。

x = pcg(A,b);

pcg stopped at iteration 20 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 20) has relative residual 3.6e-06.

既定では pcg は 20 回の反復と 1e-6 の許容誤差を使用し、アルゴリズムはこの行列について、その 20 回の反復内で収束できません。ただし、残差は許容誤差に近いため、アルゴリズムは収束するまでにより多くの反復が必要になる可能性が高くなります。

1e-7 の許容誤差と 150 回の反復を使用して、再度方程式系を解きます。

x = pcg(A,b,1e-7,150);

pcg converged at iteration 129 to a solution with relative residual 1e-07.

前処理行列を使用した `pcg` の使用

線形システムを解くために前処理行列を pcg と使用する効果を調べます。

対称正定値の、バンド係数行列を作成します。

A = delsq(numgrid('S',102));

線形方程式 $Ax = b$ の右辺の b を定義します。

b = ones(size(A,1),1);

許容誤差と最大反復回数を設定します。

tol = 1e-8;
maxit = 100;

pcg を使用して、要求された許容誤差と反復回数で解を求めます。解法プロセスに関する情報を返す出力を 5 つ指定します。

x は A*x = b の計算された解です。
fl0 はアルゴリズムが収束したかどうかを示すフラグです。
rr0 は計算解 x の相対残差です。
it0 は x が計算されたときの反復回数です。
rv0 は $‖ b - Ax ‖$ の残差履歴のベクトルです。

[x,fl0,rr0,it0,rv0] = pcg(A,b,tol,maxit);
fl0

fl0 = 1

rr0

rr0 = 0.0131

it0

it0 = 100

pcg は、要求された反復回数 100 回以内で、要求された許容誤差 1e-8 に収束しなかったため、fl0 は、1 です。

遅い収束への対応として、前処理行列を指定できます。A は対称であるため、ichol を使用して前処理行列 $M = L L^{T}$ を生成します。pcg への入力として L および L' を指定して、前処理された方程式系を解きます。

L = ichol(A);
[x1,fl1,rr1,it1,rv1] = pcg(A,b,tol,maxit,L,L');
fl1

fl1 = 0

rr1

rr1 = 8.0992e-09

it1

it1 = 79

ichol 前処理行列を使用すると、79 回目の反復で 1e-8 の指定の許容誤差より少ない相対残差が生成されます。出力 rv1(1) は norm(b) であり、rv1(end) は norm(b-A*x1) です。

次に、michol オプションを使用して修正不完全コレスキー前処理行列を作成します。

L = ichol(A,struct('michol','on'));
[x2,fl2,rr2,it2,rv2] = pcg(A,b,tol,maxit,L,L');
fl2

fl2 = 0

rr2

rr2 = 9.9565e-09

it2

it2 = 47

この前処理行列は、この例の係数行列に対するゼロ充てんの不完全コレスキー分解によって作成される前処理行列より優れているため、pcg はさらに迅速に収束させることができます。

初期推定 (反復数 0) から始めて各残差履歴をプロットすることで、前提条件が pcg の収束レートにどのように影響するかを確認できます。指定された許容誤差の線を追加します。

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(b),'-o')
semilogy(0:length(rv2)-1,rv2/norm(b),'-o')
yline(tol,'r--');
legend('No Preconditioner','Default ICHOL','Modified ICHOL','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Figure contains an axes object. The axes object with xlabel Iteration number, ylabel Relative residual contains 4 objects of type line, constantline. These objects represent No Preconditioner, Default ICHOL, Modified ICHOL, Tolerance.

初期推定の指定

pcg に解の初期推定を指定する効果を調べます。

三重対角スパース行列を作成します。 $x$ の想定される解が 1 のベクトルとなるよう、 $Ax = b$ の右辺のベクトルとして各行の合計を使用します。

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

pcg を使用して $Ax = b$ を 2 回解きます。1 回は既定の初期推定、もう 1 回は解の適切な初期推定を使用します。両方の解に対して 200 回の反復と既定の許容誤差を使用します。すべての要素が 0.99 と等価のベクトルとして初期推定を 2 番目の解に指定します。

maxit = 200;
x1 = pcg(A,b,[],maxit);

pcg converged at iteration 35 to a solution with relative residual 9.5e-07.

x0 = 0.99*e;
x2 = pcg(A,b,[],maxit,[],[],x0);

pcg converged at iteration 7 to a solution with relative residual 8.7e-07.

この場合、初期推定を指定すると pcg をより迅速に収束させることができます。

中間結果を返す

for ループで pcg を呼び出して、初期推定を使用して中間結果を取得することもできます。ソルバーを呼び出すたびに、数回の反復が行われ、計算された解が格納されます。その後、その解を次の反復のバッチに対する初期ベクトルとして使用します。

たとえば、次のコードは 100 回の反復を 4 回実行し、for ループを通過するたびに、解のベクトルを格納します。

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = pcg(A,b,tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

X(:,k) は、for ループの反復 k で計算された解のベクトルであり、R(k) はその解の相対残差です。

数値行列に代わる関数ハンドルの使用

pcg に、係数行列 A の代わりに A*x を計算する関数ハンドルを与えて線形システムを解きます。

gallery を使用して、20 行 20 列の正定値三重対角行列を生成します。上対角と下対角には 1 があり、主対角要素の数は 20 から 1 までです。行列をプレビューします。

n = 20;
A = gallery('tridiag',ones(n-1,1),n:-1:1,ones(n-1,1));
full(A)

ans = 20×20

    20     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1    19     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1    18     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1    17     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1    16     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1    15     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1    14     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1    13     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1    12     1     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1    11     1     0     0     0     0     0     0     0     0     0
      ⋮

この三重対角行列は特殊な構造であるため、演算 A*x を関数ハンドルで表すことができます。A がベクトルを乗算する場合、結果のベクトルのほとんどの要素はゼロとなります。結果の非ゼロ要素は、A の非ゼロの三重対角要素に対応します。さらに、主対角のみに 1 と等しくない非ゼロ要素があります。

式 $Ax$ は次のようになります。

$[\begin{array}{c} 20 & 1 & 0 & \dots & \dots & \dots & 0 & 0 \\ 1 & 19 & 1 & 0 & 0 \\ 0 & 1 & 18 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 17 & 1 & 0 \\ 0 & 1 & 16 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 15 & 1 & 0 \\ 0 & 1 & 14 & 1 & 0 & ⋮ \\ ⋮ & 0 & 1 & 13 & ⋱ & 0 \\ 0 & 0 & ⋱ & ⋱ & 1 \\ 0 & 0 & \dots & \dots & \dots & 0 & 1 & 1 \end{array}] [\begin{array}{c} x_{1} \\ x_{2} \\ x_{3} \\ x_{4} \\ x_{5} \\ ⋮ \\ ⋮ \\ x_{20} \end{array}] = [\begin{array}{c} 2 {0 x}_{1} + x_{2} \\ x_{1} + 19 x_{2} + x_{3} \\ x_{2} + 18 x_{3} + x_{4} \\ ⋮ \\ x_{18} + 2 x_{19} + x_{20} \\ x_{19} + x_{20} \end{array}]$ .

結果のベクトルは、3 つのベクトルの合計として記述できます。

$[\begin{array}{c} 2 {0 x}_{1} + x_{2} \\ x_{1} + 19 x_{2} + x_{3} \\ x_{2} + 18 x_{3} + x_{4} \\ ⋮ \\ x_{18} + 2 x_{19} + x_{20} \\ x_{19} + x_{20} \end{array}]$ = $[\begin{array}{c} 0 \\ x_{1} \\ ⋮ \\ x_{19} \end{array}] + [\begin{array}{c} 20 x_{1} \\ 19 x_{2} \\ ⋮ \\ x_{20} \end{array}] + [\begin{array}{c} x_{2} \\ ⋮ \\ x_{20} \\ 0 \end{array}]$

MATLAB® で、これらのベクトルを作成および合算することにより A*x の値を与える関数を記述します。

function y = afun(x)
y = [0; x(1:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 0];
end

(この関数は、ローカル関数として例の最後に保存されています)

ここで、pcg に A*x を計算する関数ハンドルを与えて、線形システム $Ax = b$ を解きます。1e-12 の許容誤差と 50 回の反復を使用します。

b = ones(20,1);
tol = 1e-12;  
maxit = 50;
x1 = pcg(@afun,b,tol,maxit)

pcg converged at iteration 20 to a solution with relative residual 4.4e-16.

afun(x1) が 1 のベクトルを生成していることを確認します。

afun(x1)

ローカル関数

function y = afun(x)
y = [0; x(1:19)] + ...
    [(20:-1:1)'].*x + ...
    [x(2:20); 0];
end

入力引数

`A` — 係数行列
行列 | 関数ハンドル

係数行列。対称正定値行列または関数ハンドルとして指定します。この行列は線形システム A*x = b の係数行列です。一般に、A は大規模なスパース行列、または大規模なスパース行列と列ベクトルの積を返す関数ハンドルです。A が対称正定値であることを確認する方法の詳細については、行列が対称正定値かどうかの判別を参照してください。

関数ハンドルとしての `A` の指定

オプションで、係数行列を、行列ではなく関数ハンドルとして指定できます。この関数ハンドルは、係数行列全体を作成するのではなく、行列とベクトルの積を返し、計算をより効率的にします。

関数ハンドルを使用するには、関数シグネチャ function y = afun(x) を使用します。関数のパラメーター化では、必要な場合に関数 afun に追加のパラメーターを指定する方法を説明しています。関数呼び出し afun(x) は A*x の値を返さなければなりません。

データ型: double | function_handle
複素数のサポート: あり

`b` — 線形方程式の右辺
列ベクトル

線形方程式の右辺。列ベクトルとして指定します。b の長さは size(A,1) と等しくなければなりません。

データ型: double
複素数のサポート: あり

`tol` — メソッドの許容誤差
`[]` または `1e-6` (既定値) | 正のスカラー

メソッドの許容誤差。正のスカラーとして指定します。この入力を使用して計算の精度と実行時間とをトレードオフします。成功するには、pcg は、許容される反復回数内に許容誤差を満たさなければなりません。tol を小さい値にすると、解がより正確でなければ計算が成功しないことを意味します。

データ型: double

`maxit` — 最大反復回数
`[]` または `min(size(A,1),20)` (既定値) | 正のスカラー整数

最大反復回数。正のスカラー整数として指定します。maxit の値を増やして、pcg が許容誤差 tol を満たすためにより多くの反復を行えるようにします。一般に、tol の値が小さいほど、計算を正常に完了するための反復がより多く必要であるということを意味します。

`M`, `M1`, `M2` — 前処理行列 (個別の引数として)
`eye(size(A))` (既定値) | 行列 | 関数ハンドル

前処理行列。行列または関数ハンドルの個別の引数として指定します。前処理行列 M またはその行列の因子 M = M1*M2 を指定して線形システムの数値的側面を改善し、pcg での迅速な収束を容易にします。不完全行列因数分解関数 ilu および ichol を使用して前処理行列を生成できます。また、因数分解の前に equilibrate を使用して、係数行列の条件数を改善できます。前処理行列の詳細については、線形方程式の反復法を参照してください。

pcg は指定なし前処理行列を単位行列として扱います。

関数ハンドルとしての `M` の指定

オプションで、M、M1 または M2 のいずれかを、行列ではなく関数ハンドルとして指定できます。関数ハンドルは、前処理行列全体を作成するのではなく行列ベクトル演算を行い、計算をより効率化します。

関数ハンドルを使用するには、関数シグネチャ function y = mfun(x) を使用します。関数のパラメーター化では、必要な場合に関数 mfun に追加のパラメーターを指定する方法を説明しています。関数呼び出し mfun(x) は M\x または M2\(M1\x) の値を返さなければなりません。

データ型: double | function_handle
複素数のサポート: あり

`x0` — 初期推定
`[]` またはゼロの列ベクトル (既定値) | 列ベクトル

初期推定。size(A,2) に等しい長さの列ベクトルとして指定します。pcg に、既定値のゼロのベクトルではなく、より妥当な初期推定 x0 を与えることができれば、計算時間が節約でき、アルゴリズムがより速く収束するようにできます。

データ型: double
複素数のサポート: あり

出力引数

`x` — 線形システムの解
列ベクトル

線形システムの解。列ベクトルとして返されます。この出力は線形システム A*x = b の近似解を与えます。計算が正常に完了した場合 (flag = 0)、relres は tol 以下になります。

計算が正常に完了しない場合 (flag ~= 0) は常に、pcg によって返される解 x は、すべての反復にわたり計算された最小残差ノルムをもちます。

`flag` — 収束フラグ
スカラー

収束フラグ。次の表のいずれかのスカラー値として返されます。収束フラグは、計算が正常に実行されたかどうかを示し、さまざまな失敗の形式を区別します。

フラグの値	収束
`0`	成功 — `pcg` は希望する許容誤差 `tol` に反復回数 `maxit` 以内で収束しました。
`1`	失敗 — `pcg` は `maxit` 回の反復を行いましたが収束しませんでした。
`2`	失敗 — 前処理行列 `M` または `M = M1*M2` が悪条件になっています。
`3`	失敗 — 連続する 2 回の反復が同じであり、`pcg` は計算を進めていません。
`4`	失敗 — `pcg` アルゴリズムで計算されたスカラー量の 1 つが、計算を続けるには大きすぎるかまたは小さすぎます。

`relres` — 相対残差誤差
スカラー

相対残差誤差。スカラーとして返されます。相対残差誤差 relres = norm(b-A*x)/norm(b) は解の正確さを示しています。計算が許容誤差 tol に反復回数 maxit 以内で収束する場合、relres <= tol になります。

データ型: double

`iter` — 反復回数
スカラー

反復回数。スカラーとして返されます。この出力は、x に対する計算解が計算された時点での反復回数を示します。

データ型: double

`resvec` — 残差誤差
ベクトル

残差誤差。ベクトルとして返されます。残差誤差 norm(b-A*x) によって、アルゴリズムが x の特定の値に対する収束にどの程度近いかがわかります。resvec の要素数は反復の回数と同じです。resvec の内容を調べて、tol または maxit の値を変更するかどうかを決定することができます。

データ型: double

詳細