Transcript Lagrange未定乗数法

制約条件付問題より生じる
線形方程式反復解法 の
理論的な諸問題について
鷲尾 巧 （ＪＳＴ・東京大学）
久田 俊明、鈴木 健二（東京大学）
制約条件付き変分問題
で生じる行列
エネルギー:
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