equilibrate
条件を改善するための行列のスケーリング
説明
例
条件を改善するために行列を平衡化
大きな条件数をもつ行列を平衡化して、反復ソルバー gmres
による線形系の解の効率と安定性を改善します。
行列 west0479
を読み込みます。これは 479 行 479 列の実数値のスパース行列です。condest
を使用して、行列の推定条件数を計算します。
load west0479
A = west0479;
c1 = condest(A)
c1 = 1.4244e+12
反復 450 回、許容誤差 1e-11
を指定して gmres
を使用し、線形系 を解いてみます。gmres
が反復ごとに解の残差ノルムを返すように、5 つの出力を指定します (~
を使用して不要な出力を抑制)。残差ノルムを片対数プロットでプロットします。このプロットは、gmres
が妥当な残差ノルムに到達できず、 について計算された解が信頼できないことを示します。
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
により返された出力を使用して、問題 を式 に書き換えることができます。ここで、、 です。この形式では、元の系の解を で復元できます。ただし、復元された解は、元の系の望ましい残差誤差をもたない可能性があります。詳細については、再スケーリングによる線形系の求解を参照してください。
gmres
を使用して を について解き、反復ごとの残差ノルムを再プロットします。このプロットは、行列の平衡化によって問題の安定性が高まり、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
出力形式の指定
6 行 6 列の魔方陣行列を作成し、equilibrate
を使用して行列を因数分解します。既定では、equilibrate
は置換情報およびスケーリング係数を行列として返します。
A = magic(6)
A = 6×6
35 1 6 26 19 24
3 32 7 21 23 25
31 9 2 22 27 20
8 28 33 17 10 15
30 5 34 12 14 16
4 36 29 13 18 11
[P,R,C] = equilibrate(A)
P = 6×6
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 1 0 0
1 0 0 0 0 0
0 0 1 0 0 0
0 1 0 0 0 0
R = 6×6
0.1852 0 0 0 0 0
0 0.1749 0 0 0 0
0 0 0.1909 0 0 0
0 0 0 0.1588 0 0
0 0 0 0 0.1793 0
0 0 0 0 0 0.1966
C = 6×6
0.1799 0 0 0 0 0
0 0.1588 0 0 0 0
0 0 0.1588 0 0 0
0 0 0 0.2422 0 0
0 0 0 0 0.2066 0
0 0 0 0 0 0.2035
行列の因子 P
、R
、および C
を使用して平衡行列 Bmatrix
を作成します。
Bmatrix = R*P*A*C
Bmatrix = 6×6
1.0000 0.1471 1.0000 0.5385 0.5358 0.6031
0.1259 1.0000 0.8056 0.5509 0.6506 0.3916
0.2747 0.8485 1.0000 0.7859 0.3943 0.5825
1.0000 0.0252 0.1513 1.0000 0.6233 0.7754
1.0000 0.2562 0.0569 0.9553 1.0000 0.7295
0.1061 0.9988 0.2185 1.0000 0.9341 1.0000
次に、"vector"
オプションを指定して出力をベクトルとして返します。入力行列が大きい場合、出力をベクトルとして返すとメモリ消費量が少なくなり、効率性が向上します。
[p,r,c] = equilibrate(A,"vector")
p = 6×1
5
6
4
1
3
2
r = 6×1
0.1852
0.1749
0.1909
0.1588
0.1793
0.1966
c = 6×1
0.1799
0.1588
0.1588
0.2422
0.2066
0.2035
ベクトル出力 p
、r
、および c
を使用して平衡行列 Bvector
を作成します。
Bvector = r.*A(p,:).*c'
Bvector = 6×6
1.0000 0.1471 1.0000 0.5385 0.5358 0.6031
0.1259 1.0000 0.8056 0.5509 0.6506 0.3916
0.2747 0.8485 1.0000 0.7859 0.3943 0.5825
1.0000 0.0252 0.1513 1.0000 0.6233 0.7754
1.0000 0.2562 0.0569 0.9553 1.0000 0.7295
0.1061 0.9988 0.2185 1.0000 0.9341 1.0000
Bvector
と Bmatrix
を比較します。結果で等しいことが示されます。
norm(Bmatrix - Bvector)
ans = 0
入力引数
A
— 入力行列
正方行列
入力行列。正方行列として指定します。A
は密でもスパースでも構いませんが、sprank
の判定で構造的に非特異でなければなりません。
A
がスパースである場合、出力 P
、R
、および C
もスパースになります。
データ型: single
| double
複素数のサポート: あり
outputForm
— 出力形式
"matrix"
(既定値) | "vector"
出力引数
P
— 置換情報
行列 (既定値) | ベクトル
置換情報。行列またはベクトルとして返されます。P*A
(または "vector"
オプションの場合は A(P,:)
) は、対角要素の積の絶対値が最大になるように A
を並べ替えたものです。
R
— 行のスケーリング
対角行列 (既定値) | ベクトル
行のスケーリング。対角行列またはベクトルとして返されます。R
と C
の非ゼロのエントリは正の実数です。
C
— 列のスケーリング
対角行列 (既定値) | ベクトル
列のスケーリング。対角行列またはベクトルとして返されます。R
と C
の非ゼロのエントリは正の実数です。
詳細
再スケーリングによる線形系の求解
線形系の解 x = A\b
では、A
の条件数が計算の精度と効率性のために重要になります。equilibrate
では、基底ベクトルを再スケーリングすることで A
の条件数を改善させることができます。これによって、b
と x
の両方を表現できる新しい座標系が効果的に形成されます。
equilibrate
は、元の系 x = A\b
のベクトル b
および x
のスケールが無関係である場合に最も役立ちます。ただし、b
と x
のスケールが関連 "している" 場合は、equilibrate
を使用して、線形系の求解に対して A
のみを再スケーリングすることはお勧めできません。得られた解では通常、元の基底ベクトルを使用して表された場合でも、元の系の小さい残差が得られません。
参照
[1] Duff, I. S., and J. Koster. “On Algorithms For Permuting Large Entries to the Diagonal of a Sparse Matrix.” SIAM Journal on Matrix Analysis and Applications 22, no. 4 (January 2001): 973–96.
拡張機能
スレッドベースの環境
MATLAB® の backgroundPool
を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool
を使用してコードを高速化します。
この関数はスレッドベースの環境を完全にサポートしています。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。
バージョン履歴
R2019a で導入R2022a: 出力形式の指定
outputForm
を "vector"
または "matrix"
に指定することで、equilibrate
が出力引数をベクトルと行列のどちらとして返すかを制御します。因数分解が大規模な場合、出力をベクトルとして返すとメモリ消費量が少なくなり、効率性が向上します。
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)