ドキュメンテーション

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

lsqr

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)
[x,flag] = lsqr(A,b,tol,maxit,M1,M2,x0)
[x,flag,relres] = lsqr(A,b,tol,maxit,M1,M2,x0)
[x,flag,relres,iter] = lsqr(A,b,tol,maxit,M1,M2,x0)
[x,flag,relres,iter,resvec] = lsqr(A,b,tol,maxit,M1,M2,x0)
[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,tol,maxit,M1,M2,x0)

説明

x = lsqr(A,b) は、A に矛盾がない場合、A*x=b の線形方程式系を x について解きます。それ以外の場合は、norm(b-A*x) を最小にする最小二乗解 x を求めようとします。mn 列の係数行列 A は、正方である必要はありませんが、大規模なスパース行列である必要があります。列ベクトル b の長さは m でなければなりません。A は、関数ハンドル afun として指定することもできます (たとえば、afun(x,'notransp')A*x を返す、afun(x,'transp')A'*x を返すなど)。

関数のパラメーター化では、必要に応じて関数 afun および以下で説明する前提条件関数 mfun に追加のパラメーターを設定する方法が説明されています。

関数 lsqr が収束する場合、それを伝えるメッセージが表示されます。関数 lsqr が設定された最大反復回数に達しても収束しないか、または何らかの理由で停止した場合、警告メッセージが、相対残差 norm(b-A*x)/norm(b) と、計算が停止した反復回数と共に表示されます。

lsqr(A,b,tol) は、このメソッドの許容誤差を指定します。tol[] の場合、関数 lsqr は既定の設定の 1e-6 を使用します。

lsqr(A,b,tol,maxit) は、最大反復回数を指定します。

lsqr(A,b,tol,maxit,M) および lsqr(A,b,tol,maxit,M1,M2) は、nn 列の前提条件子 M あるいは M = M1*M2 を使用し、方程式 A*inv(M)*y = by について効率的に解きます。ここで、y = M*x です。M[] の場合、関数 lsqr は前提条件を適用しません。M は、mfun(x,'notransp')M\x を返し、mfun(x,'transp')M'\x を返すような関数 mfun です。

lsqr(A,b,tol,maxit,M1,M2,x0) は、n1 列の初期推定を指定します。x0[] の場合、関数 lsqr は既定の設定のゼロ ベクトルを使用します。

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

フラグ

収束

0

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

1

関数 lsqr は、maxit 回の反復に対して、収束しませんでした。

2

前提条件 M は悪条件です。

3

関数 lsqr の出力結果に変化が見られません(連続する 2 回の計算で、結果がまったく同じでした)。

4

関数 lsqr で計算されたスカラー量の 1 つが計算を続けるには、大きすぎるかまたは小さすぎます。

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

[x,flag,relres] = lsqr(A,b,tol,maxit,M1,M2,x0) は、相対残差の推定 norm(b-A*x)/norm(b) も返します。flag0 の場合は、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) を含めて、残差ノルムの推定のベクトルも返します。

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,tol,maxit,M1,M2,x0) は、各反復で、スケーリングされた正規方程式の残差の推定値を要素としたベクトルも返します。すなわち、norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro') です。norm(A*inv(M),'fro') の推定は変化し、各反復で希望する改良を加えることができます。

例 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 = lsqr(A,b,tol,maxit,M1,M2);

以下のメッセージが表示されます。

lsqr converged at iteration 11 to a solution with relative 
residual 3.5e-009

例 2

この例では、例 1 の行列 A を行列ベクトル作成関数 afun のハンドルで置き換えます。この例は、以下を行う関数 run_lsqr に含まれています。

  • 第 1 引数として関数ハンドル @afun を指定し、関数 lsqr を呼び出す。

  • run_lsqr のすべての変数が afun に適用できるように、入れ子関数として afun を含みます。

以下に run_lsqr のコードを示します。

function x1 = run_lsqr
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);
x1 = lsqr(@afun,b,tol,maxit,M1,M2);

    function y = afun(x,transp_flag)
       if strcmp(transp_flag,'transp')      % y = A'*x
          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);
       elseif strcmp(transp_flag,'notransp') % y = A*x
          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
    end
end

次のように入力した場合、

x1=run_lsqr;

MATLAB® ソフトウェアが以下のメッセージを表示します。

lsqr converged at iteration 11 to a solution with relative 
residual 3.5e-009

参考文献

[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] Paige, C. C. and M. A. Saunders, "LSQR: An Algorithm for Sparse Linear Equations And Sparse Least Squares," ACM Trans. Math. Soft., Vol.8, 1982, pp. 43-71.

拡張機能

R2006a より前に導入