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

lsqr

線形方程式系の求解 — 最小二乗法

説明

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

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

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

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

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

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

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

また、[x,flag,relres] = lsqr(___) は、計算された解 x の残差誤差を返します。flag0 の場合、xnorm(b-A*x) を最小化する最小二乗解です。relres が小さい場合、x も定数の解です。これは、relresnorm(b-A*x)/norm(b) を表すためです。

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

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

[x,flag,relres,iter,resvec,lsvec] = lsqr(___) はまた lsvec を返します。これは、各反復でのスケーリングされた正規方程式の誤差の推定です。

すべて折りたたむ

lsqr を既定の設定で使用して方形線形システムを解き、次にその求解のプロセスで使用される許容誤差と反復の回数を調整します。

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

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

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

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

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

1e-4 の許容誤差と 70 回の反復を使用して、再度方程式系を解きます。計算された解の相対残差 relres を返す 6 つの出力を、残差履歴 resvec および最小二乗残差履歴 lsvec とともに指定します。

[x,flag,relres,iter,resvec,lsvec] = lsqr(A,b,1e-4,70);
flag
flag = 0

flag が 0 であるため、アルゴリズムは指定された反復回数内に、想定される許容誤差を満たすことができました。一般的に、このようにして、許容誤差と反復回数を共に調整して、速度と精度の間のトレードオフを図ることができます。

計算された解の相対残差と最小二乗残差を調べます。

relres
relres = 0.2625
lsres = lsvec(end)
lsres = 2.7640e-04

これらの残差ノルムは、x が最小二乗解であることを示しています。relres が指定された 1e-4 の許容誤差よりも小さいためです。線形システムに定数解が存在しないため、ソルバーとしての最善は、最小二乗残差が許容誤差を満足するようにすることです。

残差履歴をプロットします。相対残差 resvec は最小値にすばやく到達し、それ以上進むことはできません。一方、最小二乗残差 lsvec は後続の反復で継続して最小化し続けます。

N = length(resvec);
semilogy(1:N,lsvec,'--o',1:N,resvec,'-o')
legend("Least-squares residual","Relative residual")

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

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

load west0479
A = west0479;

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

b = sum(A,2);

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

tol = 1e-12;
maxit = 20;

lsqr を使用して、要求された許容誤差と反復回数で解を求めます。この求解プロセスについての情報を返す 6 つの出力を指定します。

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

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

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

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

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

  • lsrv0 は最小二乗残差履歴のベクトルです。

[x0,fl0,rr0,it0,rv0,lsrv0] = lsqr(A,b,tol,maxit);
fl0
fl0 = 1
rr0
rr0 = 0.0017
it0
it0 = 20

fl0 = 1 であるため、アルゴリズムは最大反復回数内に指定された許容誤差に収束しませんでした。

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

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

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

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

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')

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

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

A = sprand(700,900,0.1);
b = sum(A,2);

lsqr を使用して Ax=b を 2 回解きます。1 回は既定のオプション値を使用し、1 回は解の初期推定を使用します。両方の求解で 75 回の反復を使用し、すべての要素が 0.99 に等しいベクトルとして初期推定を指定します。

maxit = 75;
x1 = lsqr(A,b,[],maxit);
lsqr converged at iteration 64 to a solution with relative residual 8.7e-07.
x0 = 0.99*ones(size(A,2),1);
x2 = lsqr(A,b,[],maxit,[],[],x0);
lsqr converged at iteration 26 to a solution with relative residual 9.6e-07.

初期推定が想定される解に近いため、lsqr はより少ない反復で収束できます。

中間結果を返す

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

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

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

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

lsqr に、係数行列 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

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

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

b = full(sum(A,2));
tol = 1e-6;  
maxit = 25;
x1 = lsqr(@afun,b,tol,maxit)
lsqr converged at iteration 21 to a solution with relative residual 5.4e-13.
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
複素数のサポート: あり

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

データ型: double

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

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

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

関数ハンドルとしての 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) に等しい長さの列ベクトルとして指定します。lsqr に、既定値のゼロのベクトルではなく、より妥当な初期推定 x0 を与えることができれば、計算時間が節約でき、アルゴリズムがより速く収束するようにできます。

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

出力引数

すべて折りたたむ

線形システムの解。ベクトルとして返されます。この出力は線形システム A*x = b の近似解を与えます。

  • flag0 であり、かつ relres <= tol である場合、xA*x = b の定数解です。

  • flag0 であるが relres > tol である場合、xnorm(b-A*x) を最小化する最小二乗解です。

計算が成功しない (flag ~= 0) 場合はいつでも、lsqr によって返される解 x は、すべての反復で計算された最小ノルム残差をもつ解です。

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

フラグの値

収束

0

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

1

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

2

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

3

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

4

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

相対残差誤差。スカラーとして返されます。相対残差誤差は、返された答え x がどれだけ正確かを示すものです。lsqr は、求解のプロセスの各反復で相対残差および最小二乗残差を追跡し、"いずれかの" 残差が指定された許容誤差 tol を満たしたときにアルゴリズムは収束します。relres の出力には、収束した残差の値 (相対残差または最小二乗残差) が含まれています。

  • 相対残差誤差は norm(b-A*x)/norm(b) に等しく、一般には、lsqr が収束するときの許容誤差 tol を満たす残差です。resvec の出力は、すべての反復にわたってこの残差の履歴を追跡します。

  • 最小二乗残差誤差は norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro') に等しくなります。この残差は lsqr が収束する回数が相対残差よりも少なくなる原因となります。lsvec の出力は、すべての反復にわたってこの残差の履歴を追跡します。

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

データ型: double

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

データ型: double

スケーリングされた正規方程式の誤差。ベクトルとして返されます。各反復について、lsvec にはスケーリングされた正規方程式の残差 norm((A*inv(M))'*(B-A*X))/norm(A*inv(M),'fro') の推定が含まれています。lsvec の要素数は反復の回数と同じです。

詳細

すべて折りたたむ

最小二乗法

最小二乗 (LSQR) アルゴリズムは、方形行列に対する共役勾配 (CG) メソッドの適応です。解析的には、A*x = b の LSQR は、正規方程式 A'*A*x = A'*b の CG と同じ残差を生じますが、LSQR はより適した数値特性をもち、そのため一般的にはより信頼性があります[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] 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 より前に導入