lsqr (MATLAB Function 参考文献)

正規方程式に対して、共役勾配法による LSQR 実現を行います。

表示

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

詳細

x = lsqr(A,b) は、A に矛盾がない場合、A*x=b の線形方程式を、x について解こうとします。他の場合は、norm(b-A*x) を最小にする最小二乗解 x を求めようとします。m 行 n 列の係数行列 A は、正方である必要はありませんが、列ベクトル b は、長さが、m でなければなりません。 A は、A*x を出力するような関数 afun(x) です。すなわち、afun(x,'transp') は、A'*x を戻します。

lsqr が収束する場合、その影響に関するメッセージが表示されます。lsqr が、設定された最大繰り返し回数に達しても、また、何らかの理由で収束しない場合、ワーニングメッセージが、相対残差 norm(b-A*x)/norm(b) と、計算が停止した繰り返し回数と共に印刷表示されます。

lsqr(A,b,tol) は、手法に関するトレランスを指定します。tol が、[] の場合、lsqr は、デフォルトの 1e-6 を使います。

lsqr(A,b,tol,maxit) は、繰り返し最大回数を設定します。maxit が、[] の場合、lsqr は、デフォルト値、すなわち、min([m,n,20]) を使います。

lsqr(A,b,tol,maxit,M1) と lsqr(A,b,tol,maxit,M1,M2) は、n 行 n 列の前置条件子 M 、または、M = M1*M2 を使い、システム A*inv(M)*y = b を y について解きます。ここで、x = M*y の関係があります。M が、[] の場合、lsqr は、前置条件子を適用しません。M は、mfun(x) が、M\x を戻し、mfun(x,'transp') が、M'\x を戻します。

lsqr(A,b,tol,maxit,M1,M2,x0) は、code>n 行 1 列の初期推定を指定します。 x0 が、[] の場合、lsqr は、デフォルトのすべてゼロベクトルを使います。

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

[x,flag] = lsqr(A,b,tol,maxit,M1,M2,x0) は、収束フラグも戻します。

フラグ
収束

0
lsqr は、maxit の繰り返しの回数の中で、希望するトレランス tol に収束します。

1
lsqr は、maxit 回の繰り返しに対して、収束していないことを示しています。

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

3
lsqr の結果に変化が見られません(連続する2回の計算で、結果が全く同じです)。

4

lsqr の実行中に計算されるスカラ量の一つが、計算を続けるために、大き過ぎるか、小さ過ぎます。

フラグ	収束
`0`	`lsqr` は、`maxit` の繰り返しの回数の中で、希望するトレランス `tol` に収束します。
`1`	`lsqr` は、`maxit` 回の繰り返しに対して、収束していないことを示しています。
`2`	前置条件子 `M` は、条件数が良くありません。
`3`	`lsqr` の結果に変化が見られません(連続する2回の計算で、結果が全く同じです)。
`4`	`lsqr` の実行中に計算されるスカラ量の一つが、計算を続けるために、大き過ぎるか、小さ過ぎます。

flag が、0 でない場合は、常に、出力される解 x は、繰り返し全体に渡り、計算される最小のノルム残差内になります。flag 出力が設定されている場合、メッセージは表示されません。

[x,flag,relres] = lsqr(A,b,tol,maxit,M1,M2,x0) は、相対残差 norm(b-A*x)/norm(b) も出力します。flag が、0 の場合、relres <= tol となります。

[x,flag,relres,iter] = lsqr(A,b,tol,maxit,M1,M2,x0) は、x が計算されたときの繰り返し回数も出力します。ここで、0 <= iter <= maxit です。

[x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,x0) は、各繰り返しで、norm(b-A*x0) に含まれる残差ノルムの推定のベクトルも出力します。

例題

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 = lsqr(A,b,tol,maxit,M1,M2,[]);
lsqr converged at iteration 12 to a solution with relative 
residual 3.5e-009

また、つぎの行列-ベクトル積関数を使います。

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

lsqr への入力と同じです。

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

参考

bicg, bicgstab, cgs, gmres, minres, pcg, qmr, 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.

Paige, C. C. and M. A. Saunders, "LSQR: An アルゴリズム for Sparse Linear Equations And Sparse Least Squares," ACM Trans. Math. Soft., Vol.8, 1982, pp. 43-71.

lsqnonneg lu