Main Content

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

qmr

線形方程式系の求解 — 準最小残差法

説明

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

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

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

x = qmr(A,b,tol,maxit,M) は前処理行列 M を指定し、システム M1Ax=M1b を x について効果的に解くことにより x を計算します。前処理行列を使用すると、問題の数値的なプロパティと計算の効率を向上させることができます。

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

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

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

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

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

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

すべて折りたたむ

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

密度が 50% の乱数スパース行列 A を作成します。また、Ax=b の右辺に乱数のベクトル b を作成します。

rng default
A = sprand(400,400,.5);
A = A'*A;
b = rand(400,1);

qmr を使用して Ax=b を解きます。出力の表示には相対残差誤差 Ax-bb の値が含まれます。

x = qmr(A,b);
qmr 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 0.12.

既定では qmr は 20 回の反復と 1e-6 の許容誤差を使用し、アルゴリズムではこの行列について、その 40 回の反復で収束できません。残差がまだ大きいため、より多くの反復 (または前処理行列) が必要であることを示しています。許容誤差を減らして、アルゴリズムの収束をより簡単にすることもできます。

1e-4 の許容誤差と 100 回の反復を使用して、再度システムを解きます。

x = qmr(A,b,1e-4,100);
qmr stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.063.

許容誤差が緩く、反復回数が多い場合でも、残差誤差はあまり向上しません。反復アルゴリズムがこの方法で停滞する場合、前処理行列が必要であることを示しています。

A の不完全コレスキー分解を計算し、qmr への前処理行列入力として L' 因子を使用します。

L = ichol(A);
x = qmr(A,b,1e-4,100,L');
qmr converged at iteration 57 to a solution with relative residual 6.1e-05.

前処理行列を使用すると、qmr が収束できるまで問題の数値特性が向上します。

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

479 行 479 列の非対称実スパース行列 west0479 を読み込みます。

load west0479
A = west0479;

真の解がすべて 1 のベクトルになるように、b を定義します。

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

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

  • x0A*x0 = b の計算された解です。

  • fl0 はアルゴリズムが収束したかどうかを示すフラグです。

  • rr0 は計算解 x0 の残差です。

  • it0x0 が計算されたときの反復回数です。

  • rv0Ax-b の残差履歴のベクトルです。

[x0,fl0,rr0,it0,rv0] = qmr(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.7984
it0
it0 = 17

fl0 は、qmr が要求した反復回数 20 回以内に要求した許容誤差 1e-12 に収束しなかったため、1 となります。17 番目の反復が最良の近似解であり、it0 = 17 の指定によって返されます。

遅い収束への対応として、前処理行列を指定できます。A は非対称であるため、ilu を使用して前処理行列 M=LU を生成します。棄却許容誤差を指定して、1e-6 よりも小さい値をもつ非対角エントリを無視します。qmr への入力として L および U を指定して、前処理されたシステム M-1Ax=M-1b を解きます。

setup = struct('type','ilutp','droptol',1e-6);
[L,U] = ilu(A,setup);
[x1,fl1,rr1,it1,rv1] = qmr(A,b,tol,maxit,L,U);
fl1
fl1 = 0
rr1
rr1 = 4.1557e-14
it1
it1 = 6

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

各反復での相対残差をプロットして、qmr の進行状況を確認できます。指定された許容誤差のラインと共に、それぞれの解の残差履歴をプロットします。

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

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

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

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

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

maxit = 200;
x1 = qmr(A,b,[],maxit);
qmr converged at iteration 27 to a solution with relative residual 9.5e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = qmr(A,b,[],maxit,[],[],x0);
qmr converged at iteration 7 to a solution with relative residual 6.7e-07.

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

中間結果を返す

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

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

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

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

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

非対称の三重対角行列を作成します。行列をプレビューします。

A = gallery('wilk',21) + diag(ones(20,1),1)
A = 21×21

    10     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     2     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     2     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     2     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     2     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     2     0     0     0     0     0     0     0     0     0     0
      ⋮

この三重対角行列は特殊な構造であるため、演算 A*x を関数ハンドルで表すことができます。A の各行は x の要素を乗算するため、いくつかの結果のみが非ゼロです (三重対角行列の非ゼロに対応)。

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

Ax=[1020019200120010010200110][x1x2x3x21]=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21]

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

Ax=[10x1+2x2x1+9x2+2x3x19+9x20+2x21x20+10x21]=[0x1x2x20]+[10x19x29x2010x21]+2[x2x3x210]

同様に、ATx の式は次のようになります。

ATx=[1010029100210020010100210][x1x2x3x21]=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21]

ATx=[10x1+x22x1+9x2+x32x19+9x20+x212x20+10x21]=2[0x1x2x20]+[10x19x29x2010x21]+[x2x3x210]

MATLAB® では、これらのベクトルを作成する関数を記述し、それらを合計して、フラグの入力に基づいて A*x または A'*x の値を与えるようにします。

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

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

ここで、qmrA*x および A'*x を計算する関数ハンドルを与えて、線形システム Ax=b を解きます。1e-6 の許容誤差と 25 回の反復を使用します。A の行の合計として b を指定し、x の真の解が 1 のベクトルとなるようにします。

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = qmr(@afun,b,tol,maxit)
qmr converged at iteration 19 to a solution with relative residual 4.7e-07.
x1 = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

ローカル関数

function y = afun(x,flag)
if strcmp(flag,'notransp') % Compute A*x
    y = [0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + 2*[x(2:end); 0];
elseif strcmp(flag,'transp') % Compute A'*x
    y = 2*[0; x(1:20)] ...
        + [(10:-1:0)'; (1:10)'].*x ...
        + [x(2:end); 0];
end
end

入力引数

すべて折りたたむ

係数行列。正方行列または関数ハンドルとして指定します。この行列は線形システム A*x = b の係数行列です。一般に、A は大規模なスパース行列、または大規模なスパース行列と列ベクトルの積を返す関数ハンドルです。

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

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

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

  • afun(x,'notransp') は、A*x の積を返します。

  • afun(x,'transp') は、A'*x の積を返します。

許容される関数の例を次に示します。

function y = afun(x,opt,B,C,n)
if strcmp(opt,'notransp')
    y = [B*x(n+1:end); C*x(1:n)];
else
    y = [C'*x(n+1:end); B'*x(1:n)];
end
関数 afunB および C を使用して、実際にスパース行列 A = [zeros(n) B; C zeros(n)] 全体を作成することなく、指定されたフラグに応じて A*x または A'*x を計算します。これは行列のスパース パターンを利用して、A*x および A'*x を計算するときにメモリを節約します。

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

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

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

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

データ型: double

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

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

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

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

MM1、または M2 のいずれも、行列ではなく関数ハンドルとして指定して、計算でメモリを節約できます。関数ハンドルは、前処理行列全体を作成するのではなく行列ベクトル演算を行い、計算をより効率化します。

関数ハンドルを使用するには、まずシグネチャ function y = mfun(x,opt) をもつ関数を作成します。関数のパラメーター化では、必要な場合に関数 mfun に追加のパラメーターを指定する方法を説明しています。関数 mfun は、次の条件を満たさなければなりません。

  • mfun(x,'notransp')M\x または M2\(M1\x) の値を返します。

  • mfun(x,'transp')M'\x または M1'\(M2'\x) の値を返します。

許容される関数の例を次に示します。

function y = mfun(x,opt,a,b)  
if strcmp(opt,'notransp')
    y = x.*a;
else
    y = x.*b;
end
end
この例では、関数 mfuna および b を使用して、実際にスパース行列 M 全体を作成することなく、指定されたフラグに応じて M\x = x*a または M'\x = x*b を計算します。

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

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

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

出力引数

すべて折りたたむ

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

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

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

フラグの値

収束

0

成功 — qmr は希望する許容誤差 tol に反復回数 maxit 以内で収束しました。

1

失敗 — qmrmaxit 回の反復を行いましたが収束しませんでした。

2

失敗 — 前処理行列 M または M = M1*M2 が悪条件になっています。

3

失敗 — 連続する 2 回の反復が同じであり、qmr は計算を進めていません。

4

失敗 — qmr アルゴリズムで計算されたスカラー量の 1 つが、計算を続けるには大きすぎるかまたは小さすぎます。

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

データ型: double

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

データ型: double

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

データ型: double

詳細

すべて折りたたむ

準最小残差法

QMR アルゴリズムは、BiCG を改良したものとして開発されました。GMRES は Krylov 部分空間の直交基底を使用して最小残差解を計算しますが、QMR は双直交基底を使用するため、準最小残差解を計算します。

QMR は通常、BiCG よりもさらに滑らかに収束します。また、先読みの手法を使用してほぼすべての場合にブレークダウンを回避します。QMR の計算コストは、BiCG より若干高くなるだけです[1]

ヒント

  • 最大反復のメソッドの収束は、係数行列の条件数 cond(A) に依存します。equilibrate を使用して A の条件数を改善することができ、それ自体で最大反復のソルバーが収束しやすくなります。ただし、equilibrate を使用することでも、平衡化した行列 B = R*P*A*C を次に因子分解する際に前処理行列の品質が向上します。

  • dissectsymrcm などの行列の並べ替え関数を使用して係数行列の行と列を並べ替え、係数行列が因子分解されて前処理行列が生成される際に非ゼロの数を最小化できます。これによって、後で前処理を使用した線形システムを解くのに必要なメモリと時間を削減できます。

参照

[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] Freund, Roland W. and Nöel M. Nachtigal, “QMR: A quasi-minimal residual method for non-Hermitian linear systems,” SIAM Journal: Numer. Math. 60, 1991, pp. 315–339.

拡張機能

R2006a より前に導入