Main Content

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

線形方程式の反復法

数値線形代数の最も重要で一般的な用途の 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 は、行列の除算演算子 /\ や、lsqminnormdecompositionlinsolve などの関数を通じて、直接法を実装します。

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

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

    A = delsq(numgrid('L',400));
    b = ones(size(A,1),1);
    x = pcg(A,b,[],1000);
    norm(A*x-b)
    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. 残差などの計算された量の値に基づいて、ベクトル x0 の大きさと方向を更新します。

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

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

各反復法の概要

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

説明

メモ

lsqr (最小二乗法)

  • 係数行列は方形にできます。

  • 方形システムに使用できる唯一のソルバーです。

  • 正規方程式 (A'*A)*x = A'*b に適用された共役勾配法 (PCG) と解析的には同等ですが、より好ましい数値特性をもっています。

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 を使用できます。双共役傾斜法系のアルゴリズム (bicgbicgstabcgs など) が gmres よりも効率的な場合もありますが、収束動作が予測できないため、最初の選択肢としては gmres が優れていることが多くなります。

前処理行列

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

Ax=b

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

A˜x˜=b˜.

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

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

(M1A)x=(M1b).

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

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

(AM1)(Mx)=b.

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

最後に、対称係数行列では、"分離前処理" を行うことで新しいシステムも確実に対称行列のままとなります。前処理行列 M=HHT は A の両側に現れます。

(H1AHT)HTx=(H1b)

分離前処理を適用したシステムのソルバー アルゴリズムは上記の方程式に基づきますが、実際には H を計算する必要はありません。ソルバー アルゴリズムは M または inv(M) を直接適用するだけです。

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

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

MATLAB の反復ソルバーでは、単一の前処理行列 M、または M = M1M2 となる 2 つの前処理行列因子を指定できます。これにより、M = LU のように因数分解から計算される前処理行列を簡単に指定できます。分離前処理で M = HHT も成立する場合、H 因子をもつ M1 入力と M2 入力の間には直接接続はありません。

多くの場合、前処理行列は与えられた問題の数学的なモデルの中に自然に存在します。たとえば変数係数をもつ偏微分方程式は、定数係数をもつ偏微分方程式で近似されます。自然な前処理行列がない場合、次の表に示す不完全因数分解のいずれかを使用できます。

関数因数分解説明
ilu

ALU

正方行列または方形行列用の不完全 LU 分解です。
ichol

ALL*

対称正定行列用の不完全コレスキー分解です。

iluichol の詳細については、不完全因数分解を参照してください。

前処理行列の例

正方、2 次元領域上のラプラス方程式の 5 点の有限差分近似を考えます。次のコマンドは、前処理を使用した共役勾配法 (PCG) の前処理行列 M = L*L' = *' を使用します。ここで L は、A のゼロで埋めた不完全コレスキー因子です。

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. dissectsymrcm などのスパース行列の並べ替え関数を使用して、平衡行列を並べ替えます。

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

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

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

    load west0479
    A = west0479;
    n = size(A,1);
    b = ones(n,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 の三重対角因子を乗算する因子のみです。したがって、与えられたベクトル 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 を計算しなければなりません。ソルバー lsqrqmrbicg については、線形演算子関数は要求されたときに 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.

関連するトピック