ドキュメンテーション

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

equilibrate

条件を改善するための行列のスケーリング

説明

[P,R,C] = equilibrate(A) は、新しい行列 B = R*P*A*C の対角要素が大きさ 1 であり、非対角要素の大きさが 1 未満になるように、行列 A を並べ替えて再スケーリングします。

すべて折りたたむ

大きな条件数をもつ行列を平衡化して、反復ソルバー gmres による線形系の解の効率と安定性を改善します。

行列 west0479 を読み込みます。これは 479 行 479 列の実数値のスパース行列です。condest を使用して、行列の推定条件数を計算します。

load west0479
A = west0479;
c1 = condest(A)
c1 = 1.4244e+12

反復 450 回、許容誤差 1e-11 を指定して gmres を使用し、線形系 Ax=b を解いてみます。gmres が反復ごとに解の残差ノルムを返すように、5 つの出力を指定します (~ を使用して不要な出力を抑制)。残差ノルムを片対数プロットでプロットします。このプロットは、gmres が妥当な残差ノルムに到達できず、x について計算された解が信頼できないことを示します。

b = ones(size(A,1),1);
tol = 1e-11;
maxit = 450;
[x,flx,~,~,rvx] = gmres(A,b,[],tol,maxit);
semilogy(rvx)
title('Residual Norm at Each Iteration')

equilibrate を使用して、A の並べ替えと再スケーリングをします。より良い条件数をもち、かつ対角要素が 1 と -1 のみの新しい行列 B = R*P*A*C を作成します。

[P,R,C] = equilibrate(A);
B = R*P*A*C;
c2 = condest(B)
c2 = 5.1036e+04

equilibrate により返された出力を使用して、問題 Ax=b を式 By=d に書き換えることができます。ここで、B=RPACd=RPb です。この形式では、元の系の解を x=Cy で回復できます。

gmres を使用して By=dy について解き、反復ごとの残差ノルムを再プロットします。このプロットは、行列の平衡化によって問題の安定性が高まり、gmres が反復 200 回未満で目的の許容誤差 1e-11 に収束することを示します。

d = R*P*b;
[y,fly,~,~,rvy] = gmres(B,d,[],tol,maxit);
hold on
semilogy(rvy)
legend('Original', 'Equilibrated', 'Location', 'southeast')
title('Relative Residual Norms (No Preconditioner)')
hold off

前提条件子を使用した解の改善

行列 B を得た後に、gmres と共に使用する前提条件子を計算して、さらに問題の安定性を改善することができます。B の数値特性は元の行列 A の数値特性よりも優れているため、平衡化した行列を使用して前提条件子を計算することが望まれます。

ilu を使用して 2 つの異なる前提条件子を計算し、それらを gmres の入力として使用して問題を再び解きます。比較のために、平衡化したノルムと同じプロットに反復ごとの残差ノルムをプロットします。このプロットは、平衡化した行列を使用した前提条件子の計算によって問題の安定性が大幅に高まり、gmres が反復 30 回未満で目的の許容誤差に達することを示します。

semilogy(rvy)
hold on

[L1,U1] = ilu(B,struct('type','ilutp','droptol',1e-1,'thresh',0));
[yp1,flyp1,~,~,rvyp1] = gmres(B,d,[],tol,maxit,L1,U1);
semilogy(rvyp1)

[L2,U2] = ilu(B,struct('type','ilutp','droptol',1e-2,'thresh',0));
[yp2,flyp2,~,~,rvyp2] = gmres(B,d,[],tol,maxit,L2,U2);
semilogy(rvyp2)

legend('No preconditioner', 'ILUTP(1e-1)', 'ILUTP(1e-2)')
title('Relative Residual Norms with ILU Preconditioner (Equilibrated)')
hold off

入力引数

すべて折りたたむ

入力行列。正方行列として指定します。A は密でもスパースでも構いませんが、sprank の判定で構造的に非特異でなければなりません。

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

出力引数

すべて折りたたむ

置換行列。スパース行列として返されます。P*A は、A の対角要素の積の絶対値を最大にする A の並べ替えです。

行のスケーリング。スパース対角行列として返されます。RC の対角要素は正の実数です。

列のスケーリング。スパース対角行列として返されます。RC の対角要素は正の実数です。

R2019a で導入