qmr (MATLAB Function 参考文献)

準最小残差法

表示

x = qmr(A,b)
qmr(A,b,tol)
qmr(A,b,tol,maxit)
qmr(A,b,tol,maxit,M)
qmr(A,b,tol,maxit,M1,M2)
qmr(A,b,tol,maxit,M1,M2,x0)
qmr(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...)
[x,flag] = qmr(A,b,...)
[x,flag,relres] = qmr(A,b,...)
[x,flag,relres,iter] = qmr(A,b,...)
[x,flag,relres,iter,resvec] = qmr(A,b,...)

詳細

x = qmr(A,b) は、線形方程式 A*x=b を x について、解きます。n 行 n 列の係数行列 A は正方で、列ベクトル b は長さ n である必要があります。A は、afun(x) が A*x を出力し、afun(x,'transp') が A'*x を出力するような afun 関数です。

qmr が収束する場合、そのことを示すメッセージが表示されます。qmr が、最大繰り返し設定回数に達した後でも収束できなかったり、何らかの理由で、繰り返しが停止された場合、相対残差 norm(b-A*x)/norm(b) や停止したときや収束しなかったときに使用した繰り返し回数を含む情報がワーニングメッセージに表れます。

qmr(A,b,tol) は、許容誤差を指定します。tol が [] の場合、qmr は、デフォルトの 1e-6 を使います。

qmr(A,b,tol,maxit) は、さらに最大の繰り返し回数を指定します。maxit が[] の場合、qmr は、デフォルト min(n,20) を使います。

qmr(A,b,tol,maxit,M) と qmr(A,b,tol,maxit,M1,M2) は、前提条件子 M、または、M = M1*M2 を使い、システム inv(M1)*A*inv(M2)*y = inv(M1)*b を y について、効率良く解きます。ここで、x = inv(M2)*y の関係が成り立ちます。M が、[] の場合、qmr は、前提条件子を適用しません。M は、mfun(x) が M\x を戻したり、mfun(x,'transp') が M'\x 戻したりする関数 mfun です。

qmr(A,b,tol,maxit,M1,M2,x0) は、初期推定を指定します。x0 が、[] の場合、qmr は、デフォルトのゼロベクトルを使います。

qmr(afun,b,tol,maxit,m1fun,m2fun,x0,p1,p2,...) は、パラメータ p1,p2,... を関数afun(x,p1,p2,...) と afun(x,p1,p2,...,'transp') に渡し、同時に、前提条件子関数 m1fun と m2fun にも、それぞれ渡します。

[x,flag] = qmr(A,b,...) は、収束フラグも出力します。

フラグ
収束

0
qmr は、繰り返し回数maxit内で、希望する許容誤差tolに収束しました。

1
qmr は、maxit回の繰り返しを行いましたが、収束しませんでした。

2
前提条件子 M は、条件数が良くありません。

3
qmr は計算を進めていません(連続する2回の繰り返しの結果が同じでした)。

4
qmr 中で計算されたスカラ量が、繰り返しを続行するのに小さする、または大きすぎるため、この手法は停止されました。

フラグ	収束
`0`	`qmr` は、繰り返し回数`maxit`内で、希望する許容誤差`tol`に収束しました。
`1`	`qmr` は、`maxit`回の繰り返しを行いましたが、収束しませんでした。
`2`	前提条件子 `M` は、条件数が良くありません。
`3`	`qmr` は計算を進めていません(連続する2回の繰り返しの結果が同じでした)。
`4`	`qmr` 中で計算されたスカラ量が、繰り返しを続行するのに小さする、または大きすぎるため、この手法は停止されました。

flag が 0 でないとき、出力された解 x は、すべての繰り返しで計算された中で、最小のノルム残差をもつ解です。flag の出力が指定されている場合は、メッセージは表示されません。

[x,flag,relres] = qmr(A,b,...) は、相対誤差norm(b-A*x)/norm(b)も出力します。flagが0 の場合、これは relres <= tolです。

[x,flag,relres,iter] = qmr(A,b,...) は、xが計算された繰り返し回数を出力します。これは、常に0 <= iter <= maxit です。

[x,flag,relres,iter,resvec] = qmr(A,b,...) は、norm(b-A*x0) を含んで、各繰り返しで、残差ノルムのベクトルも計算します。

例題

例題 1.

n = 100;
on = ones(n,1);
A = spdiags([-2*on 4*on -on],-1:1,n,n);
b = sum(A,2);
tol = 1e-8; maxit = 15;
M1 = spdiags([on/(-2) on],-1:0,n,n);
M2 = spdiags([4*on -on],0:1,n,n);
x = qmr(A,b,tol,maxit,M1,M2,[]);

一方、つぎの行列-ベクトル乗算関数を使います。

function y = afun(x,n,transp_flag)
if (nargin > 2) & strcmp(transp_flag,'transp')
   y = 4 * x;
   y(1:n-1) = y(1:n-1) - 2 * x(2:n);
   y(2:n) = y(2:n) - x(1:n-1);
else
   y = 4 * x;
   y(2:n) = y(2:n) - 2 * x(1:n-1);
   y(1:n-1) = y(1:n-1) - x(2:n);
end

を、qmr へ入力とします。

x1 = qmr(@afun,b,tol,maxit,M1,M2,[],n);

例題 2.

load west0479;
A = west0479;
b = sum(A,2);
[x,flag] = qmr(A,b)

qmr は、デフォルトの20回の繰り返し回数以内で、デフォルトの許容範囲 1e-6 に収束しなかったので、flag は 1 になります。

[L1,U1] = luinc(A,1e-5);
[x1,flag1] = qmr(A,b,1e-6,20,L1,U1)

上三角 U1 が対角要素にゼロをもつので、flag1 は、2 になります。そして、U1*y = r for y をバックシュラッシュを使って解く場合、qmr は、最初の繰り返しで停止します。

[L2,U2] = luinc(A,1e-6);
[x2,flag2,relres2,iter2,resvec2] = qmr(A,b,1e-15,10,L2,U2)

qmr は、8回の繰り返しで、1.6571e-016(relres2 の値) の許容範囲で収束するので、1e-6 のドロップトレランスを使った不完全 LU 分解による前置条件子の場合、flag2 は、0 になります。resvec2(1) = norm(b) と resvec2(9) = norm(b-A*x2) です。初期推定(繰り返し回数 0)から始めて、各繰り返しで、相対残差をプロットすることで、qmr の進行を見れます。

semilogy(0:iter2,resvec2/norm(b),'-o')
xlabel('iteration number')
ylabel('relative residual')

参考

bicg, bicgstab, cgs, gmres, lsqr, luinc, minres, pcg, symmlq

@ (関数ハンドル), \ (バックスラッシュ)

参考文献

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

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.

pwd qr