1. No category

数理計画の実務と
半正定値対称行列
田辺隆人
[email protected]
[email protected]
（株）NTTデータ数理システム
数理計画法アルゴリズム
•
•
•
•
•
•
•
•
線形計画法（単体法、内点法）
二次計画法（有効制約法、内点法）
非線形計画法（逐次二次計画法，内点法）
非線形半正定値計画（内点法）
線形混合整数計画法（分枝限定法）
混合整数二次計画法（分枝限定法）
非線形整数計画法（分枝限定法）
重み付き制約充足問題（タブ・サーチ）
数理計画法アルゴリズム
•
•
•
•
線形計画法（単体法、内点法）
二次計画法（有効制約法、内点法）
非線形計画法（逐次二次計画法，内点法）
非線形半正定値計画（内点法）
•
•
•
•
線形混合整数計画法（分枝限定法）
混合整数二次計画法（分枝限定法）
非線形整数計画法（分枝限定法）
重み付き制約充足問題（タブ・サーチ）
数理計画法アルゴリズム
•
•
•
•
線形計画法（単体法、内点法）
二次計画法（有効制約法、内点法）
非線形計画法（逐次二次計画法，内点法）
非線形半正定値計画（内点法）
•
•
•
•
線形混合整数計画法（分枝限定法）
混合整数二次計画法（分枝限定法）
非線形整数計画法（分枝限定法）
重み付き制約充足問題（タブ・サーチ）
線形計画問題（標準形）
最小化
条　件　xR ,
mn
Ax  b, A  R ,
t
c x,
x0
n
変数の分割..
n
m
b
A
x
c
変数の分割
n
m
AB
AN
xB  b
cB
xN
cN
線形計画問題の性質
適当な分割：
非基底
x  ( xB | xN )  ( xB | 0)
基底
とした解（基底解）の中に
最小化問題の解は存在する
最適解の全体
基底解
単体法アルゴリズム
• 分割を探す⇔（基底）解を探す
AB xB  AN xN  b
x  ( xB | xN )  ( xB | 0)
1
B
xB  A b
分割があれば
解の取得は容易
単体法アルゴリズム
変数一つを選択し
て基底変数に
できるだけ目的関数に実質
貢献するものを選ぼう
（プライシング）
x  ( xB | x N )
行列
対応する一つを
非基底変数に
AB1とABT
効率よく更新しよう．．
線形計画法のＫＫＴ条件
（最適解の必要十分条件）
Ax  b,
c  A y  z  0,
t
Xz  0,
x  0, z  0
X  diag ( x)
最適な基底解は
KKT条件を満たす
Ax  b,
c  A y  z  0,
x  ( xB | 0)
t
z  (0 | cN  A y)
t
N
T
B
B
1
B
y  A c
xB  A b
Xz  0,
x  0, z  0
最適解の全体
最適な基底解
最適な基底解は
KKT条件を満たす
問題の性質に特化
単体法
より汎用的
内点法
Ax  b,
c  A y  z  0,
x  ( xB | 0)
t
z  (0 | cN  A y)
t
N
T
B
B
1
B
y  A c
xB  A b
Xz  0,
x  0, z  0
最適解の全体
最適な基底解
内点法アルゴリズム
（線形計画法）
• 次を「非線形方程式」として解こう！
Ax  b,
c  A y  z  0,
t
Xz  0,
x  0, z  0
X  diag ( x)
内点法アルゴリズム
（線形計画法）
• 次を「非線形方程式」として解こう！
Ax  b,
c  A y  z  0,
解法テクニック
Xz   ,
t
x 0 , z  0
Newton法
f ( x  x)
f ( x)  f ( x)x  0
f ( x)
1
x  (f ( x)) f ( x)
x
ステップサイズの取得（一次方程式解法）
に計算負荷が集中
内点法のステップ方向取得
• 一次方程式の形
D
A
x
A
O
y
t
 b
正定値対称行列への帰着
• Normal Equation Form (Adler,Karmarkar’89)
D
A
A
t
O
不定値行列
ピボット選択不要
コレスキー分解の存在保障
1
AD A  B
t
正定値対称行列
係数行列の更新
for(k)｛
for(i,j)｛
B[i,j] += A[k,i]・A[k,j]・D[k,k]
｝
反復中不変
｝
特殊なデータ構造
Normal Equation Formの欠点
• Dense column 問題
D
A
疎行列
A
t
1
O
t
AD A  B
密行列
Dense Column 問題の解決
• Column 分割(Gonzio’94)
• Rank-One-Update の利用 (Andersen’94)
• Normal Equation Form の放棄
Augumented System
Approach(Maros,Meszaros’95)
エレガントな解決：
脱・正定値対称行列
• Augumented System Approach
列の並べ替え・変形
D
密
な
列
A
疎行列
A
t
O
 As Ds As
疎行列
t
不定値だが
ピボット選択
は不要
Meszaros’95
内点法の特徴
• 最適性条件に依拠
⇒汎用的な解法の枠組み
• 汎用かつ効率的な実装？
正定値対称行列への帰着は必須ではない
非線形最適化問題の実装
（脱・正定値対称行列）
• NLP用Augumented System Approach
並べ替え・変形
GD
A
密
な
列
Bunch-Palett
分解
Bunch,
Palett
’71
ならば分解可能
A
t
O
 As Ds As
t
疎な不定値行列
（後付けの）まとめ
• 正定値対称性は
直接法による一次方程式解法には
「オーバースペック」
• 単純な実装
≠ よい実装
金融工学アプリケーションと
正定値対称性との出会い
G  L L, R R  I , G  R L (LR)
t
t
t t
• 相関行列 G に従った乱数列が欲しい
• G を感覚的に手修正
• G のコレスキー分解が失敗
⇒ソフトベンダーに電話：
「本件の修正には
どのくらいの工数が必要でしょうか」
エレガントな解決：
半正定値計画問題（ＳＤＰ）を内点法で解く
最小化
GX
制約
X    I : 半正定値
F
 ：最小固有値
最も近い相関行列の生成：
Matrix Calibration
実行
ポートフォリオ
R( x) 

rj x j  r x

xj 1 , xj  0
jAsset
jAsset
t
ポートフォリオ最適化
E  R( x)  
V  R( x)  

rj x j  r x  r0

Qij xi x j  x Qx
jAsset
i , jAsset
t
t
最小化
1
k
k
k
k
Qij  ri  E (ri ) rj  E (rj )
K k
ポートフォリオ最適化
分散最小化
均等分配
定式化の選択
• Full Covariance
最小化

i , jAsset
Qij xi x j
• コンパクト分解
最小化
制約

s
sk 

kSample
2
k
jAsset
( Rkj  R j )x j
定式化の選択
• Full Covariance
最小化
半正定値

i , jAsset
• コンパクト分解
最小化
制約
Qij xi x j
3000銘柄でも1～20秒程度

s
sk 

kSample
2
k
jAsset
( Rkj  R j )x j
銘柄数制約付き
平均分散モデル
最小化
制約
t
x Qx

jAsset
xj 1 ,
x j  0 または　x j [l j , u j ]
supp( x)  K
supp( x)  {i | xi  0}
とりあえず分枝限定法（MIQP）
「場合分け」＋限定
最小化
制約
t
x Qx

xj 1 ,

j
jAsset
K
j
l j   j  x j  u j   j ,  j {0,1}
限界あり．二次の目的関数は特に難しい
凸緩和してSDP/SOCPに帰着
目的関数の書きかえとラグランジュ緩和
x Qx  x (Q  diag (d )) x  z diag (d ) z
t
正定値になるよ
うに d >0を選ぶ
t
t
xz
Ax  b

jAsset
xj 1
Az  b

jAsset
zj 1
z j  0 または　z j [l j , u j ]
Xiaojin Zheng, Xiaoling Sun, Duan Li,
Convex Relaxations and Mixed-Integer Quadratic
Reformulations for Cardinality Constrained Quadratic Programs
supp( z )  K
凸緩和してSDP/SOCPに帰着
目的関数の書きかえとラグランジュ緩和
x Qx  x (Q  diag (d )) x  z diag (d ) z
t
t
t
xz

jAsset
xj 1

jAsset
緩和（乗数：λ）
zj 1
z j  0 または　z j [l j , u j ]
supp( z )  K
凸緩和してSDP/SOCPに帰着
下界値の最大化問題
最大化
min{ f x ( x;  , d )}  min{ f z ( z; , d )}
x
z


jAsset
xj 1
jAsset
zj 1
z j  0 または　z j [l j , u j ]
supp( z )  K
( , d ) を解く問題⇒SDP，( ) を解く問題⇒SOCP
凸緩和してSDP/SOCPに帰着
• SDP/SOCPの定式化
⇒下界値の導出
⇒ MIQPの補強による高速化
上下界値の改善
エレガントな解決：
StQP
StQPに関する事実の応用
最小化
制約
t
x Qx

jAsset
xj 1 , xj  0
半正定値に
限らない
StQP の最適解は
t
x Qx が強凸であるようなフェース ( I ) の
relative interior に 存在する．
Francesco Cesarone,Andrea Scozzari, and Fabio
Tardella, A new method for mean-variance portfolio
optimization with cardinality constraints
エレガントな解決：
StQPに関する事実の応用
最小化
制約
これも StQP
t
x Qx

jI  Asset
xj 1 , xj  0
I K
x Qx が強凸であるようなフェース ( I ) を
t
I  Asset , I  K について
調べ尽くせば銘柄数制約付き分散モデルの
大域的最適解を見逃すことはない
エレガントな解決：
StQPに関する事実の応用
最適解を与える集合
正定値な Q の部分行列の
*
インデックスの集合の中に I はある
⇒ インデックス I 0 について
Q の部分行列が正定値でなければ，I 0
を含むあらゆる集合は I * に成り得ない
⇒ Increasing Set Algorithm
複数のベンチマーク問題に厳密解を
与えることに成功
StQPに関する事実の応用
そういえば．．
ヘッセ行列が
低ランク
なことが理由では？
• ポートフォリオ最適化問題で
選択される銘柄は一般に少ない
（⇒有効制約法が有利）
long/short ポジションを許容する
リスク最小化
最小化 ( x  x ) Q( x  x )
L
L
制約
x

1
,
x
 j
j 0
L
S t
L
jAsset

jAsset
x Sj  1 , x Sj  0
x Lj  0 または　x Sj  0
S
とりあえず分枝限定法（MIQP）
最小化 ( x  x ) Q( x  x )
L
L
制約
x

1
,
x
 j
j 0
L
S t
L
S
jAsset

jAsset
x Sj  1 , x Sj  0
 j  x , (1   j )  x ,  j {0,1}
L
j
S
j
とりあえず分枝限定法（MIQP）
• 連続緩和問題の補強
収益率制約や複雑な投資制約が助けになっていた．．
とりあえず分枝限定法（MIQP）
• 制約なしの場合．．
PROBLEM_NAME
NUMBER_OF_VARIABLES
(#INTEGER/DISCRETE)
NUMBER_OF_FUNCTIONS
PROBLEM_TYPE
METHOD
<preprocess begin>..........<preprocess end>
<iteration begin>
.1.2B
up:
1e+050 lo:7.2290069e-020
#1
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
iqp
600
200
603
MINIMIZATION
ACTIVE_SET_QP
time: 0.5s:mem(Mb)=77/57:avail(Mb)=4039/1872
llen:0 #prob:0 #piv:173
211.13232 time: 33.4s:mem(Mb)=111/91:avail(Mb)=4004/1837
llen:1286 #prob:2702 #piv:181822
211.13232 time: 61.8s:mem(Mb)=139/119:avail(Mb)=3977/1808
llen:2364 #prob:5329 #piv:345892
211.13232 time:124.7s:mem(Mb)=153/133:avail(Mb)=3963/1793
llen:2904 #prob:7589 #piv:720381
211.13232 time:184.0s:mem(Mb)=177/157:avail(Mb)=3938/1768
llen:3808 #prob:9569 #piv:1046687
211.13232 time:187.2s:mem(Mb)=178/157:avail(Mb)=3939/1767
llen:3856 #prob:9689 #piv:1066187
とりあえず分枝限定法（MIQP）
• 制約なしの場合．．
PROBLEM_NAME
NUMBER_OF_VARIABLES
(#INTEGER/DISCRETE)
NUMBER_OF_FUNCTIONS
PROBLEM_TYPE
METHOD
<preprocess begin>..........<preprocess end>
<iteration begin>
.1.2B
up:
1e+050 lo:7.2290069e-020
iqp
600
200
603
MINIMIZATION
ACTIVE_SET_QP
分枝限定法は全探索と同じ：
#1
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
up:
211.13232 lo:7.2290069e-020 gap:
time: 0.5s:mem(Mb)=77/57:avail(Mb)=4039/1872
llen:0 #prob:0 #piv:173
211.13232 time: 33.4s:mem(Mb)=111/91:avail(Mb)=4004/1837
llen:1286 #prob:2702 #piv:181822
211.13232 time: 61.8s:mem(Mb)=139/119:avail(Mb)=3977/1808
llen:2364 #prob:5329 #piv:345892
211.13232 time:124.7s:mem(Mb)=153/133:avail(Mb)=3963/1793
llen:2904 #prob:7589 #piv:720381
211.13232 time:184.0s:mem(Mb)=177/157:avail(Mb)=3938/1768
llen:3808 #prob:9569 #piv:1046687
211.13232 time:187.2s:mem(Mb)=178/157:avail(Mb)=3939/1767
llen:3856 #prob:9689 #piv:1066187
緩和解は相補性を満たさず
目的関数は零
緩和⇒ＳＤＰを利用した補強の導出
• 相補性条件を緩和した問題 L( )
t
L
x
x




最小化
 S  M ( )  S 
x 
x 
L
制約

jAsset
x 1
L
j

jAsset
 Q
M ( )  
 Q  
x Sj  1
x , x  0, j  Asset
L
j
S
j
M ( ) が半正定値ならば解ける！
Q   
Q 
緩和⇒ＳＤＰを利用した補強の導出
• L( ) の最小値を最大化する問題 （ＳＤＰ）
最小化

制約
過去の反復で得られた解
変数
L ,( i )
Q

Q
x
x

 


L ,( i )
S ,( i )
   S ,( i )  
  j  xj  xj


S ,( i ) 
 x   Q Q   x  jS
L ,( i )
T
M ( )  0 （半正定値）
i  1,.., I
緩和⇒ＳＤＰを利用した補強の導出
•
n=200での実験結果
緩和⇒ＳＤＰを利用した補強の導出
•
λを用いて補強されたMIQP
t
最小化
x 
x 
 S  M ( )  S 
x 
x 
制約

x 1 , x  0

x Sj  1 , x Sj  0
L
jAsset
jAsset
L
L
j
L
j
 j  x , (1   j )  x ,  j {0,1}
L
j
S
j
緩和⇒ＳＤＰを利用した補強の導出
PROBLEM_NAME
NUMBER_OF_VARIABLES
(#INTEGER/DISCRETE)
NUMBER_OF_FUNCTIONS
PROBLEM_TYPE
METHOD
<preprocess begin>..........<preprocess end>
<iteration begin>
.1.......2B
up:
1e+050 lo:
65.069268
qp_int
600
200
603
MINIMIZATION
ACTIVE_SET_QP
up:
1e+050 lo:
65.202836
up:
1e+050 lo:
65.507298
up:
1e+050 lo:
65.507298
#1
up:
78.939007 lo:
65.507298 gap:
13.431709
#2
up:
77.947918 lo:
65.507298 gap:
12.44062
#3
up:
77.930633 lo:
65.507298 gap:
12.423335
#4
up:
77.928705 lo:
65.507298 gap:
12.421407
up:
77.928705 lo:
65.507298 gap:
12.421407
50銘柄程度ならば
厳密解が得られる場合も．．
time: 1.4s:mem(Mb)=51/32:avail(Mb)=4064/1897
llen:0 #prob:0 #piv:797
time: 65.1s:mem(Mb)=59/39:avail(Mb)=4058/1890
llen:1 #prob:349 #piv:136291
time:127.3s:mem(Mb)=59/40:avail(Mb)=4055/1888
llen:360 #prob:1119 #piv:357707
time:189.1s:mem(Mb)=75/56:avail(Mb)=4041/1874
llen:985 #prob:2369 #piv:649083
time:192.7s:mem(Mb)=77/57:avail(Mb)=4039/1872
llen:1034 #prob:2472 #piv:670000
time:192.8s:mem(Mb)=77/58:avail(Mb)=4038/1871
llen:1007 #prob:2473 #piv:670194
time:193.5s:mem(Mb)=76/57:avail(Mb)=4040/1873
llen:1006 #prob:2498 #piv:674716
time:193.8s:mem(Mb)=76/57:avail(Mb)=4039/1872
llen:1006 #prob:2508 #piv:676891
time:249.1s:mem(Mb)=88/68:avail(Mb)=4028/1860
llen:1455 #prob:3488 #piv:919435
Copositive Programmingに帰着
• 等価な定式化（Copositive Programming の 双対形）
completely positive
matrixのconeの集合
エレガントな解決：
？
偏微分方程式制約付き最適化
PDE constrained optimization
離散点での物理量
パラメータ
変数
 , i
i 
最小化  (i i )
2
i
制約
F ( , )  0
パラメータと物理量を同時に解く
観測点での値
離散化された
物理方程式
偏微分方程式制約付き最適化
動機づけ
Gassan S Abdoulaev, Kui Ren and Andreas H Hielscher,
Optical tomography as a PDE-constrained optimization problem,
Inverse Problems 21 (2005) 1507–1530
偏微分方程式制約付き最適化
「最適化」の位置づけ
協働の絶好の素材ではないでしょうか．．．
偏微分方程式制約付き最適化
最適化問題としての側面
○ 不等式制約はほぼ効かない
○ ＫＫＴ条件は線形性が強い
× 超大規模（⇒反復法必須）
× ヘッセ行列の正定値性は保障できない
（⇒正規化項の導入）
偏微分方程式制約付き最適化
解くべき一次方程式
正規化項を加えて優対角にはできる
H

A
A   x  bx 

  y  b 
W     y 
t
またお会いしましたね．．
“Saddle Point System”
Augumented Lagrangean
のペナルティの逆数
偏微分方程式制約付き最適化
内点法の場合との違い
超大規模（⇒反復法必須）
H

A
A   x  bx 
  y   b 
W     y 
t
正定値性への目配りが必要に
61
前処理なし反復法での限界
• 正規化項を手厚くする
• 次の形に変形してＣＧ法
似たような方に以前お目にかか
ったような．．
 H  AW A  x  b
t
1
計算の安定化（ＣＧ反復減少）
正規化項増大で非物理的解
エレガントな解決？
Constraint Preconditioner
H

A
1
A
 H
 
1 t
AH A  S   A
t
H
H, S
A 

W 
AH A  W
1
t
あなたは．．!
S の構成が鍵
（approximate Schur complement）
t
エレガントな解決？
Approximate Schur complement
を左辺とする一次方程式解法
•
•
•
•
•
を対角行列に置き換えてコレスキー分解
multigrid 法による前処理行列でCG法
ノルムの定義を変更してCG法
対角スケーリングしてCG法を数回適用したものを前処理
行列としてGMRES法
ピボット選択付きの分解を数ステップ行ったものを前処理
行列としてGMRES法
H
H, S
AH A  W
1
t
今度は近似で良いので自由度が高い
まとめ
• 半正定値性は結果を保障の「拠り所」
ー convex QP
ー SDP緩和
• 「距離感」の調整も必要
ー 内点法の実装（線形計画法）
ー cardinality constrained portfolio