線形方程式の反復法

数値線形代数の最も重要で一般的な用途の 1 つは、A*x = b の形で表せる線形方程式を解くことです。A が大規模なスパース行列である場合、線形方程式を反復法で解くことができます。この方法では、計算の実行時間と解の精度をトレードオフできます。このトピックでは、MATLAB^® で方程式 A*x = b を解くために使用できる反復法について説明します。

直接法と反復法

線形方程式 A*x = b を解く方法は 2 種類あります。

• 直接法は、ガウスの消去法の変型です。この方法では、LU 分解、QR 分解、コレスキー分解などの行列演算によって、個々の行列要素を直接に使用します。直接法を使用すると線形方程式を高い精度で解くことができますが、大規模なスパース行列では演算が遅くなる場合があります。直接法を使用して線形方程式を解く際の速度は、係数行列の密度と非ゼロのパターンに大きく依存します。

  たとえば、次のコードは小規模な線形方程式を解きます。

  A = magic(5);
  b = sum(A,2);
  x = A\b;
  norm(A*x-b)

  ans =
  
     1.4211e-14

  MATLAB は、行列の除算演算子 / および \ や、decomposition、lsqminnorm、linsolve などの関数を通して直接法を実装します。

• 反復法では、有限回の反復後に、線形方程式の近似解を得ます。この方法は、精度を犠牲にして実行時間を短くすることに合理性がある大規模な方程式の場合に有用です。反復法は行列とベクトルの積または抽象的な線形作用素により、係数行列を間接的にのみ使用します。反復法はどのような行列にも使用できますが、通常は、直接法で時間がかかるような大規模なスパース行列に適用されます。間接的な方法で線形方程式を解く速度は、直接法ほど係数行列の非ゼロのパターンに大きく依存しません。ただし、反復法を使用する場合、通常はそれぞれ特定の問題に合わせてパラメーターを調整する必要があります。

  たとえば、次のコードは、対称正定値の係数行列をもつ大規模でスパースな線形方程式を解きます。

  A = delsq(numgrid('L',400));
  b = ones(size(A,1),1);
  x = pcg(A,b,[],1000);
  norm(b-A*x)

  pcg converged at iteration 796 to a solution with relative residual 9.9e-07.
  
  ans =
  
     3.4285e-04

  MATLAB は、係数行列 A のプロパティによって異なる強みと弱みをもつ、さまざまな反復法を実装しています。

直接法は通常、その実行に使用できる十分なストレージがある場合、間接的な方法に比べてより速く、より一般的に適用可能なアプローチです。一般的には、最初に x = A\b の使用を試すべきです。直接法が遅すぎる場合に、反復法を試します。

汎用的な反復アルゴリズム

線形方程式を解くほとんどの反復アルゴリズムは、類似するプロセスに従います。

1. まず、解ベクトル x0 の初期推定を行います (より良い推定を指定しない限り、一般的にこれはゼロ ベクトルです)。

2. 残差ノルム res = norm(b-A*x0) を計算します。

3. 残差を、指定した許容誤差と比較します。res <= tol の場合、計算を終了して、x0 の計算解を返します。

4. A*x0 を適用し、残差などの計算された量の値に基づいて、ベクトル x0 の大きさと方向を更新します。この手順で、ほとんどの計算は完了します。

5. x0 の値が許容誤差を十分に満たすまで、手順 2 から 4 を繰り返します。

手順 4 で x0 の大きさと方向を更新する方法は各反復法で異なり、手順 2 と 3 の収束条件にもわずかな違いがある場合がありますが、上記はすべての反復ソルバーが従う基本的なプロセスを捉えています。

各反復法の概要

MATLAB には、線形方程式系に対する反復法を実装する複数の関数があります。これらの手法は、Ax = b を解いたり、ノルム ||b – Ax|| を最小化するために設計されています。これらの手法のいくつかは類似点をもち、同じ複数のアルゴリズムに基づいていますが、それぞれのアルゴリズムは特定の状況でメリットがあります ([1]、[2])。

説明	メモ
pcg (前処理付き共役勾配法)	係数行列は対称正定値でなければなりません。 格納する必要のあるベクトルの数が少ないため、対称正定値システムに最も効果的なソルバーです。
lsqr (最小二乗法)	方形システムに使用できる唯一のソルバーです。 正規方程式 (A'*A)*x = A'*b に適用された共役勾配法 (PCG) と解析的に同等です。
minres (最小残差法)	係数行列は対称でなければなりませんが、正定値である必要はありません。 このアルゴリズムは、各反復で 2 ノルムでの残差誤差が最小化されるため、ステップごとに進行があることが保証されています。 ブレークダウン (アルゴリズムが解に向かって進行することができなくなり、停止すること) が発生しません。
symmlq (対称 LQ 分解法)	係数行列は対称でなければなりませんが、正定値である必要はありません。 射影システムを解き、残差をそれ以前のすべての残差と直交状態に保ちます。 ブレークダウン (アルゴリズムが解に向かって進行することができなくなり、停止すること) が発生しません。
bicg (双共役傾斜法)	係数行列は正方でなければなりません。 bicg では計算量が少なくすみますが、収束が不規則で信頼できません。 bicg は、他の多くの反復アルゴリズムがこれを改善して開発されているため、歴史的に重要です。
bicgstab (双共役傾斜安定化法)	係数行列は正方でなければなりません。 BiCG ステップを GMRES(1) ステップと交互に使用して、安定性を向上させています。
bicgstabl (双共役傾斜安定化法 (l))	係数行列は正方でなければなりません。 BiCG ステップを GMRES(2) ステップと交互に使用して、安定性を向上させています。
cgs (共役傾斜二乗法)	係数行列は正方でなければなりません。 反復ごとに必要な演算数は bicg と同じですが、二乗残差を使用することで転置を避けています。
gmres (一般化最小残差法)	係数行列は正方でなければなりません。 残差ノルムが反復ごとに最小化されるため、最も信頼できるアルゴリズムの 1 つです。 必要な計算とストレージは、反復数に応じて線形的に上昇します。 不必要な計算とストレージを避けるためには、適切な restart 値を選択することがきわめて重要です。
qmr (準最小残差法)	係数行列は正方でなければなりません。 bicg よりも反復ごとのオーバーヘッドがわずかに多くなっていますが、安定性が向上しています。
tfqmr (転置なし準最小残差法)	係数行列は正方でなければなりません。 メモリが限られている状況で対称不定システムを解こうとする場合に最適なソルバーです。

反復ソルバーの選択

次に示す MATLAB の反復ソルバーのフローチャートで、各ソルバーについて有用となる状況の概要がわかります。一般的に、非対称正方のほとんどすべての問題には gmres を使用できます。双共役傾斜法系のアルゴリズム (bicg、bicgstab、cgs など) が gmres よりも効率的な場合もありますが、収束動作が予測できないため、最初の選択肢としては gmres が優れていることが多くなります。

前処理行列

反復法の収束率は、係数行列のスペクトル (固有値) に依存します。そのため、より好ましいスペクトル (分類された固有値、または 1 に近い条件数) をもつように線形システムを変換することで、ほとんどの反復法の収束率と安定性を改善できます。この変換を行うには、システムに "前処理行列" と呼ばれる 2 番目の行列を適用します。このプロセスにより、次の線形方程式

$$Ax=b$$

が

次の等価の方程式に変換されます。

$$\tilde{A}\tilde{x}=\tilde{b}.$$

理想的な前処理行列は、係数行列 A を単位行列に変換します。そのような前処理行列があれば、あらゆる反復法が 1 回の反復で収束するためです。実際には、優れた前処理行列を見つけるにはトレードオフが必要です。変換は、左前処理、右前処理、分離前処理の 3 つの方法のいずれかで実行します。

1 番目のケースは、前処理行列 M が A の左側に現れるため、"左前処理" と呼ばれます。

$$\left(M^{-1}A\right)x=\left(M^{-1}b\right).$$

次の反復ソルバーで、左前処理が使用されています。

• bicg

• gmres

• qmr

"右前処理" では、M が A の右側に現れます。

$$\left(A{M}^{-1}\right)\left(Mx\right)=b.$$

次の反復ソルバーで、右前処理が使用されています。

• lsqr

• bicgstab

• bicgstabl

• cgs

• tfqmr

最後に、対称係数行列 A では、"分離前処理" を行うことで変換後のシステムも確実に対称行列のままとなります。前処理行列 $$M=H{H}^T$$ は分割され、因子が A の両側に表れます。

$$\left(H^{-1}A{H}^{-T}\right)H^{T}x=\left(H^{-1}b\right)$$

分離前処理を適用したシステムのソルバー アルゴリズムは上記の方程式に基づきますが、実際には H を計算する必要はありません。ソルバー アルゴリズムは、直接 M を使用して乗算および求解します。

次の反復ソルバーで、分離前処理が使用されています。

• pcg

• minres

• symmlq

いずれの場合も、反復法の収束を加速するために前処理行列 M が選択されています。反復法の残差誤差に変化が見られなくなったり、反復間の進行が小さくなった場合、前処理行列を生成し、問題に組み込む必要があることを意味している場合がよくあります。

MATLAB の反復ソルバーでは、単一の前処理行列 M、または M = M₁M₂ となる 2 つの前処理行列因子を指定できます。これにより、M = LU のように、前処理行列を分解した形式で簡単に指定できるようになります。分離前処理において M = HH^T も成立する場合、M1 および M2 の入力と H の因子との間に関係はありません。

前処理行列は与えられた問題の数学的なモデルの中に自然に存在することもあります。自然な前処理行列がない場合、次の表に示す不完全分解のいずれかを使用して前処理行列を生成できます。不完全分解とは実質的に、短時間で計算できる不完全な直接法です。

関数	分解	説明
ilu	$$\text{A}\approx \text{LU}$$	正方行列または方形行列用の不完全 LU 分解です。
ichol	$$\text{A}\approx \text{L}{\text{L}}^{*}$$	対称正定値行列用の不完全コレスキー分解です。

ilu と ichol の詳細については、不完全分解を参照してください。

前処理行列の例

正方、2 次元領域上のラプラス方程式の 5 点の有限差分近似を考えます。以下のコマンドでは、前処理を使用した共役勾配 (PCG) 法を前処理行列 M = L*L' と共に使用します。ここで L は、A のゼロで埋めた不完全コレスキー因子です。この方程式系では、前処理行列を指定しなければ pcg は解を求めることができません。

A = delsq(numgrid('S',250));
b = ones(size(A,1),1);
tol = 1e-3;
maxit = 100;
L = ichol(A);
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 92 to a solution with relative residual 0.00076.

pcg は、指定した許容誤差に達するために 92 回の反復が必要です。しかし、別の前処理行列を使用することで、結果が改善する場合があります。たとえば、ichol を使用して修正不完全コレスキー分解を作成すると、pcg は指定した許容誤差にわずか 39 回でたどり着きます。

L = ichol(A,struct('type','nofill','michol','on'));
x = pcg(A,b,tol,maxit,L,L');

pcg converged at iteration 39 to a solution with relative residual 0.00098.

平衡と並べ替え

計算の難しい問題では、ilu または ichol が直接生成したものよりも優れた前処理行列が必要になることがあります。たとえば、より高品質の前処理行列を生成するか、実行する計算の量を最小化することが望ましい場合があります。これらの場合では、"平衡" を使用して係数行列をより対角優位にし (これにより前処理行列の品質が改善されることがある)、さらに "並べ替え" によって行列因子内の非ゼロの数を最小化する (これによりメモリ要件が緩和され、後続の計算の効率が向上することがある) ことができます。

前処理行列の生成に平衡と並べ替えの両方を使用する場合、プロセスは次のようになります。

1. 係数行列に equilibrate を使用します。

2. dissect や symrcm などのスパース行列の並べ替え関数を使用して、平衡行列を並べ替えます。

3. ilu または ichol を使用して、最終的な前処理行列を生成します。

平衡と並べ替えを使用してスパース係数行列の前処理行列を生成する例を次に示します。

1. 線形方程式の右辺について、係数行列 A と 1 のベクトル b を作成します。A の条件数の推定を計算します。

   load west0479;
   A = west0479;
   b = ones(size(A,1),1);
   condest(A)

   ans =
   
      1.4244e+12

   equilibrate を使用して、係数行列の条件数を改善します。

   [P,R,C] = equilibrate(A);
   Anew = R*P*A*C;
   bnew = R*P*b;
   condest(Anew)

   ans =
   
      5.1042e+04

2. dissect を使用して、平衡行列を並べ替えます。

   q = dissect(Anew);
   Anew = Anew(q,q);
   bnew = bnew(q);

3. 不完全 LU 分解を使用して、前処理行列を生成します。

   [L,U] = ilu(Anew);

4. 前処理行列、許容誤差 1e-10、外側の最大反復回数 50、内側の反復回数 30 を使用して、gmres で線形方程式を解きます。

   tol = 1e-10;
   maxit = 50;
   restart = 30;
   [xnew, flag, relres] = gmres(Anew,bnew,restart,tol,maxit,L,U);
   x(q) = xnew;
   x = C*x(:);

   次に、gmres (前処理行列を含む) によって返された相対残差 relres を、前処理行列なしの相対残差 resnew および平衡なしの相対残差 res と比較します。結果は、すべての線形方程式が等価であるにもかかわらず、各因子に適用される重みは手法ごとに異なり、これにより残差の値が大きく影響を受けることを示しています。

   relres
   resnew = norm(Anew*xnew - bnew) / norm(bnew)
   res = norm(A*x - b) / norm(b)

   relres =
      8.7537e-11
   resnew =
      3.6805e-08
   res =
      5.1415e-04

行列の代わりとしての線形演算子の使用

MATLAB の反復ソルバーでは、A に数値行列を指定することは "必須" ではありません。ソルバーによって実行される計算では行列-ベクトル乗算 A*x または A'*x の結果を使用しますが、それに代えて、これらの線形演算の結果を計算する関数を指定することができます。これらの量を計算する関数は、よく "線形演算子" と呼ばれます。

係数行列 A の代わりに線形演算子を使用するだけでなく、前処理行列 M の代わりにも線形演算子を使用できます。その場合、ソルバーのリファレンス ページに示されているように、関数は M\x または M'\x を計算する必要があります。

線形演算子を使用すると、A または M のパターンを利用して、ソルバーが行列を明示的に使用して非スパース行列とベクトルとの乗算を実行した場合よりも効率的に線形演算の値を計算できるようになります。また、線形演算子は通常、行列をまったく形成せずに行列-ベクトル乗算の結果を計算するため、係数行列や前処理行列を格納するためのメモリが不要になります。

たとえば、次の係数行列を考えます。

A = [2 -1  0  0  0  0;
    -1  2 -1  0  0  0;
     0 -1  2 -1  0  0;
     0  0 -1  2 -1  0;
     0  0  0 -1  2 -1;
     0  0  0  0 -1  2];

A がベクトルを乗算する場合、結果のベクトルのほとんどの要素はゼロとなります。結果の非ゼロ要素は、A の非ゼロの三重対角要素に対応します。したがって、与えられたベクトル x に対し、線形演算子関数が A*x の値を計算するために必要なのは、単純に 3 つのベクトルを加算することだけです。

function y = linearOperatorA(x)
y = -1*[0; x(1:end-1)] ...
  + 2*x ...
  + -1*[x(2:end); 0];
end

ほとんどの反復ソルバーは、A*x の値を返すために A の線形演算子関数が必要です。同様に、前処理行列 M に対しても、関数は一般的に M\x を計算しなければなりません。ソルバー lsqr、qmr、bicg では、要求に応じて、線形演算子関数は A'*x または M'\x の値も返さなければなりません。線形演算子関数の例と説明については、反復ソルバーのリファレンス ページを参照してください。

参照

[1] Barrett, R., M. Berry, T. F. Chan, et al., Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM, Philadelphia, 1994.

[2] Saad, Yousef, Iterative Methods for Sparse Linear Equations. PWS Publishing Company, 1996.

関連するトピック

• 線形方程式