制約条件付問題より生じる線形方程式反復解法の理論的な諸問題について鷲尾巧（ＪＳＴ・東京大学）久田俊明、鈴木健二（東京大学）制約条件付き変分問題で生じる行列エネルギー: 1 T E( x) : x Ax  xT f 2 A: n n x, f : n の制約条件: Bx  g  0 B : m  n (m  n) g : m rank(B)  m のもとでの変分問題を考える。ラグランジュ未定乗数法以下の汎関数の停留点を求める。 L( x, y) : E( x)  yT (Bx  g ) L  xT ( Au  BT y  f )  yT (Bx  g)  0, x, y A  B BT   x   f       0   y  g  •ラグランジュ未定乗数法でも変分問題の解を求めるわけであるが、後にみるように係数行列は不定値になる。 •保存則が含まれる方程式系では、ラグランジュ未定乗数法が陽に適用されていないにせよ、上のような構造の行列が現れる。例えば回路問題の場合第二式は電流の保存則を示すキルヒホフの法則 v3 ik Rk  v  vk , k  1,2,3 i1  i2  i3  0 i3 i1 v1 v i2 v2  R1 0  0   1 0 0 R2 0 0 R3 1 1 1  i1  v1  1  i2  v2   1  i3  v3      0  v  0  中心点での電圧ｖがラグランジュ未定乗数に対応する。制約条件付き問題の固有値は？対称行列での場合  A BT   x  x     y    y   B  C     A BT   A 0  A1 0   A BT       1   B S  0  S  0 S  B  C  S  C  BA1BT シルベスターの慣性律より正の固有値数はdim(A), 負の固有値数はrank(B) C=0 の場合は Bx = λyとなるので以下のように固有値は評価できる。   x, Ax  2Bx, y   x, Ax  1    x, Ax  2 2  Bx, Bx x, Ax2  8Bx, Bx  制約条件付き問題の固有値は？歪対称化した場合  A BT   x  x  x , x    y , y   1    y    y    B C    T  x  A B T T x ,y    y    B C      x, Ax  y, Cy  2   1 Im x, BT y   固有値は右半面に入るが、虚軸方向の分布が大きくなる可能性がある。安定な離散化とは？小さな残差ｒが大きな誤差eを生まない!! 楕円型方程式の場合はメッシュサイズｈに依存しない定数ｃが存在して Ah eh  rh  eh A  c rh 2 ラグランジュ未定乗数に関しては   Sh eh  Ch  Bh Ah1BhT eh  rh  eh 1  A1  A1BT S 1BA1 A B      1 1 S BA B  C T 2  c rh A1BT S 1   1 S  2 安定性のためにＡとＢが満たすべき条件は？ C=0 の場合を考える。 Shの最小固有値をλｈその固有ベクトルをehとすると rh  Bh Ah 1BhT eh  heh eh 2  c rh 2  c h eh 2  h  1/ c λｈがhに依存せずに下から抑えられなければならない。   e, BA B e  A 1  1/ 2  sup b, A b 2 1 1/ 2 T   1/ 2 T B e, A  2  2 b 2 1   sup x, B e  sup x, B e A1 / 2 x 1 2 T x A 1 Be 1/ 2 B e  sup A T T T  T b, B e  2  2 inf-sup条件 Shの最小固有値がｈに依存せず下から抑えられる。ｈに依存しないδ>0が存在して   inf sup xh , BhT yh  inf sup Bh xh , yh    yh 2 1 x h Ah 1 yh 2 1 x h Ah 1 inf-sup条件は簡単には満たされない（例１）  vx vy  dxdy  0  p B v  0   x y   p [-1,1] [1,1] δp [-1,-1] [1,-1] 上の非圧縮性制約条件式を四角形要素上で圧力ｐに対しては定数で、流速ｖに対しては双一次関数で補間し、離散化したとする。 BBTを計算してみると縦横方向がdecoupleした行列が生成される。 0 0 2 Cavity flow問題での圧力振動例 Tol = 1.e-6 Tol = 1.e-10 Tol = 1.e-8 ペナルティ法 vs ラグランジュ未定乗数法  A BT   x   f     y   g   B  I      y 1  Bx  g  yは簡単に消去できて、以下の正定値問題に帰着できる。 1 T  1 T  A  B B x  f  B g       これはペナルティ法の考え方に他ならないように見える。しかし、BTBはランク落ちの行列であるため、上記係数行列の条件数は小さな εに対して極端に悪くなる。逆に、ペナルティ法で与えられた問題を逆の変形をたどってラグランジュ型にすることにより固有値が分離され解き易くなるのでは？例：   W  C e  1 / 2  CcomprJ  ln J  J  1 Q  体積変化に関わるポテンシャル項     Q  b ff E ff 2  bss Ess 2  bnn Enn2  b fs E fs 2  Esf 2  b fn E fn 2  Enf 2  bns Ens2  Esn 2  二つの前処理行列 ILU factorization with special fill-ins ~T 0  DA1 0  DA  U A DA  LA B P1   ~   1    B DB  LB   0  DB   0 DB  U B  Inexact block LU factorization QA 0  QA1 0  QA B T  P2      1   B Q B    0  QB   0 QB  P1におけるフィルインの決定法 for i  1,,n for k  1,,i 1 (i, k )   for j  k  1,, n (k , j)   , (i, j)   fij : fij  fik f kk1 f kj k k i j (bi,bk,bj)=(block(i),block(k),block(j)) (bi,bk,bj) Allowed fill-ins (2,1,1) , (1,1,2) Level 0 or 1 (2,1,2) All fill-ins (1,1,1), (2,2,2) No fill-ins P1における分解安定性について ~T 0  DA1 0  DA  U A DA  LA B P1   ~   1    B DB  LB   0  DB   0 DB  U B  DB  LB DB DB UB  1 ~ 1 ~T BDA B のILU（０）分解より与えられる。 ~ ~ BDA1BT のILU分解安定性は理論的に保証できない。以下のような場合ILU分解でDBにマイナス成分が現れる - inf-sup条件が満たされていないか十分な安定化が施されていない。 - メッシュが破綻している。この場合は、 DAですでに逆符号が現れる。 ILU分解の安定性とM-matrix 理論的にILU分解の安定性が保証されている代表的な行列クラスとしてM行列が挙げられる。 AがM行列であるための条件 1. aii  0 2. aij  0, i  j 3. A : non - singular 4. A1  0 ILU分解のアルゴリズム for i  1,,n for k  1,,i 1 (i, k )   for j  k  1,, n (k , j)   , (i, j)   fij : fij  fik f kk1 f kj しかし、有限要素離散化の場合は、Poisson方程式においてさえ２番目の条件は満足されない。一方で、高次の要素を用いない限りは、複雑なポテンシャルで定義される連続体でさえ、ほとんどの場合ILU分解は安定となっている。より本質的なILU分解安定性の条件はないか？ P1-1Kの固有値分布の評価 ILU分解が成功した場合には、いつもP1-1Kの固有値は右半面に分布しているが、どうしてか？ 1 1 I 0  X 0   A BT   X T Y T  1 P1 K ~   Y Z  B  C   T  0  I 0 Z        0  A BT  I  X T Y T   X T 0   X 1 0  I     1 1   T  YX  I B  C 0 I 0 Z 0 Z          X 1 0   A   X T 0   AX T Y T  BT  1   1 T T 1 T 1 1 T   T  0 Z YX A  B BX Y  YX B  C  YX AX Y 0 Z       この部分の正定値性が保証できない。 P2 におけるQA とQB の作成法について  BT QA 0  QA1 0  QA B T  QA   B Q   1   1 T B    0  QB   0 QB  B BQA B  QB   QA ~ A , QB ~ BQA1 B T or QB ~ BA1 B T BA-1BTは密行列、しかしinf-sup条件のことを考えると条件数はO(1) 例えば： QB   DB  LB DB1DB UB  in P1 P2-1Kの固有値分布の評価  M 1 P2 K ~ H    N(I  M) (I  M)NT  T N( 2I  M)N  M  QA1/ 2 AQA1/ 2  T  N  QA1/ 2 BT QB1/ 2  z1 x  2 2 X    , x  y  1, z1  z2  1  z2 y  X*HX  z1 x*Mx  z2 y* N( 2I  M)NT y 2  2 1 Im(z1* z2 x*(I  M)NT y) | x (I  M)N y | * T x*Mx 2 y* N( 2I  M)NT y 前処理に対する二つの考え方 P-1Kを単位行列に近づけようとする考え方 QA 0  QA1 0  QA BT  P   1   B Q 0 Q B   B 0  QB    M (I  M)NT  P K~ T  N(I  M) N( 2I  M)N  -1 前処理後の行列の対角ブロックは正定値、非対角ブロックは歪対称となる。前処理後も対称性を保存しようとする考え方 QA 0  QA1 0  QA BT  P   1   B Q 0 Q 0 Q   B   B B   M -1 P K~  N(I  M)  T  N( 2I  M)N  (I  M)NT 前処理後の行列は対称、しかし第２対角ブロックは負値行列となる。収束性の数値実験結果２次元非圧縮性超弾性体問題に対する収束までの反復回数（tol=1.e-8) 4/3c （MINI)要素 h 1/8 1/16 1/32 1/64 非対称前処理 41 58 113 258 対称前処理 70 127 258 541 •どちらの場合も反復数は、ほぼ1/hに比例している。 •対称前処理は、非対称の２倍の反復数を要する。対称問題の場合（CG, MINRES法） CGの場合 ek e0 A A  min pk , pk ( 0) 1 max i 1,...,n i  [ a, b] ( a  0)  a 0 pk (i ) ek e0 MINRESの場合 rk  min max pk (i ) r0 pk , pk (0)1 i 1,...,n b A A  b / a 1    2    b / a 1  a k a=O(h2), b=O(1)ならば＃Iterations=O(1/h) b c d 0 c=O(h2), a,b,d=O(1)ならば #Iterations=O(1/h)  ad / bc 1  rk  i  [ a, b] [ c, d ] ( b  0  c)   2  ad / bc  1  r0   k /2 前処理された行列の固有値分布対称前処理の場合前処理された行列の固有値分布非対称前処理の場合非対称の場合は、固有値分布の評価より妥当な反復数評価（O(1/h))を得ることは困難。局所的な緩和法は安定化か？ rloc  (b  Axold)loc Alocxloc  rloc x new  x  xloc old E( x new  x)  ( x new  x)T A( x new  x)  E( xold  x)  xlocT Alocxloc 係数行列Ａが対称である限りは、上の等式は成立するが、Ａが不定値の場合は、安定性が不明。局所的な緩和法はMultigrid法にとって重要 Element by Element緩和法の可能性久田研究室鈴木健二博士論文より C=εh2Δhが適当な大きさならば安定な緩和法が実現できる。 P2の安定性解析 E  i 1    Ei  ノルムを定義する行列 1 with   1  A B 2  QA  BT  1 T   B QB  BQA B  3  0   AQA  A  (1   A ) QA 2  0   BQB  BA1BT  (2   B )QB  あらゆる問題に登場する制約条件 • • • • 連続体非圧縮性電磁気電荷や磁荷の保存回路問題キルヒホフの法則量子力学パウリの排他律固有値問題も E( A)  x, Ax 1  x, x  0 L A,    x, Ax  1  x, x  L  2x, Ax  1  x, x  2 x, x  2x, Ax  x  1  x, x