非線形最適化のアルゴリズムとソフトウエア

数理計画の実務と

半正定値対称行列

田辺隆人

[email protected]

[email protected]

（株）NTTデータ数理システム

数理計画法アルゴリズム

• 線形計画法（単体法、内点法）

• 二次計画法（有効制約法、内点法）

• 非線形計画法（逐次二次計画法，内点法）

• 非線形半正定値計画（内点法）

• 線形混合整数計画法（分枝限定法）

• 混合整数二次計画法（分枝限定法）

• 非線形整数計画法（分枝限定法）

• 重み付き制約充足問題（タブ・サーチ）

数理計画法アルゴリズム

• 線形計画法（単体法、内点法）

• 二次計画法（有効制約法、内点法）

• 非線形計画法（逐次二次計画法，内点法）

• 非線形半正定値計画（内点法）

• 線形混合整数計画法（分枝限定法）

• 混合整数二次計画法（分枝限定法）

• 非線形整数計画法（分枝限定法）

• 重み付き制約充足問題（タブ・サーチ）

数理計画法アルゴリズム

• 線形計画法（単体法、

内点法

_）

• 二次計画法（有効制約法、

内点法

_）

• 非線形計画法（逐次二次計画法，

内点法

_）

• 非線形半正定値計画（

内点法

_）

• 線形混合整数計画法（分枝限定法）

• 混合整数二次計画法（分枝限定法）

• 非線形整数計画法（分枝限定法）

• 重み付き制約充足問題（タブ・サーチ）

線形計画問題（標準形）

0

x

b,

Ax

,

x

x,

c

t

n







条　件　

R

最小化

,

n

m

A



R



変数の分割..

A

x

c

b



m

n

変数の分割

B

A

x

_B

N

c

b



m

n

N

A

N

x

B

c

適当な分割：

とした解（基底解）の中に

最小化問題の解は存在する

線形計画問題の性質

基底

非基底

)

0

|

(

)

|

(

x

_B

x

_N

x

_B

x





最適解の全体

基底解

単体法アルゴリズム

• 分割を探す⇔（基底）解を探す

)

0

|

(

)

|

(

x

_B

x

_N

x

_B

x





b

x

A

x

A

_B

_B



_N

_N



b

A

x

_B



_B



1

分割があれば

解の取得は

容易

)

|

(

x

_B

x

_N

x



単体法アルゴリズム

変数一つを選択し

て基底変数に

対応する一つを

非基底変数に

できるだけ目的関数に実質

貢献するものを選ぼう

（プライシング）

行列

効率よく更新しよう．．

1

T

B

B

A



と

A



線形計画法のＫＫＴ条件

（最適解の必要十分条件）

,

0,

0,

0,

0

t

Ax

b

c

A y

z

Xz

x

z





 







( )

X



diag x

,

0,

0,

0,

0

t

Ax

b

c

A y

z

Xz

x

z





 







T

B

B

y



A





c

(

_B

| 0)

x



x

(0 |

_N

_N

t

)

z



c



A y

最適解の全体

最適な基底解

1

B

B

x



A b



最適な基底解は

KKT条件を満たす

最適な基底解は

KKT条件を満たす

,

0,

0,

0,

0

t

Ax

b

c

A y

z

Xz

x

z





 







T

B

B

y



A





c

(

_B

| 0)

x



x

(0 |

_N

_N

t

)

z



c



A y

1

B

B

x



A b



単体法

内点法

問題の性質に特化

より汎用的

最適解の全体

最適な基底解

• 次を「

非

_{線形方程式」として解こう！}

内点法アルゴリズム

（線形計画法）

,

0,

0,

0,

0

t

Ax

b

c

A y

z

Xz

x

z





 







( )

X



diag x

• 次を「

非

_{線形方程式」として解こう！}

,

0,

,

0 ,

0

t

Ax

b

c

A y

z

Xz

x

z





 



内点法アルゴリズム

（線形計画法）







解法テクニック

Newton法

(

)

( )

( )

0

f x

 

x

f x



f x

 

x

1

(

( ))

( )

x

f x



f x

   

x

( )

f x

ステップサイズの取得（一次方程式解法）

に計算負荷が集中

• 一次方程式の形

内点法のステップ方向取得

t

A



A



D

O

x



y





_b

• Normal Equation Form (Adler,Karmarkar’89)

正定値対称行列への帰着

t

A



A



D

O

AD A



1

t

ピボット選択不要

コレスキー分解の存在保障

B



不定値行列

正定値対称行列

for(k)｛

for(i,j)｛

B[i,j] +=

A[k,i]・A[k,j]

・D[k,k]

｝

｝

反復中不変

特殊なデータ構造

係数行列の更新

• Dense column 問題

Normal Equation Formの欠点

1

t

AD A



_

_B

疎行列

密行列

t

A



A



D

O

• Column 分割(Gonzio’94)

• Rank-One-Update の利用 (Andersen’94)

• Normal Equation Form の放棄

_{Augumented System}

_{Approach(Maros,Meszaros’95)}

• Augumented System Approach

エレガントな解決：

脱・正定値対称行列

t

A



A



D

O



A

s

D

s

A

s

t

疎行列

疎行列

列の並べ替え・変形

不定値だが

ピボット選択

は不要

Meszaros’95

密

な

列

• 最適性条件に依拠

⇒汎用的な解法の枠組み

• 汎用かつ効率的な実装？

正定値対称行列への帰着は必須ではない

• NLP用Augumented System Approach

非線形最適化問題の実装

（脱・正定値対称行列）

t

A



A



D

G



O



A

D

A

疎な不定値

_行列

Bunch-Palett

分解

Bunch,

Palett

’71

ならば分解可能

並べ替え・変形

密

な

列

（後付けの）まとめ

• 正定値対称性は

直接法による一次方程式解法には

「オーバースペック」

金融工学アプリケーションと

正定値対称性との出会い

• 相関行列 に従った乱数列が欲しい

• を感覚的に手修正

• のコレスキー分解が失敗

⇒ソフトベンダーに電話：

「本件の修正には

どのくらいの工数が必要でしょうか」

G

G

G

,

,

(

)

t

t

t t

G



L L R R



I G



R L LR

エレガントな解決：

半正定値計画問題（ＳＤＰ）を内点法で解く

F

G



X

:

X

 



I

半正定値

最小化

制約

最も近い相関行列の生成：

Matrix Calibration



：最小固有値

ポートフォリオ

( )

_j

_j

t

j Asset

R x

r x

r x









1 ,

j

j Asset

x







x

_j



0

ポートフォリオ最適化



( )



0

t

j

j

j Asset

E R x

r x

r x

r















,

( )

_ij

_i

_j

t

i j Asset

V R x

Q x x

x Qx









最小化



 



1

(

)

(

)

k

k

k

k

ij

i

i

j

j

k

Q

r

E r

r

E r

K











均等分配

分散最小化

均等分配

定式化の選択

• Full Covariance

• コンパクト分解

,

ij

i

j

i j Asset

Q x x





最小化

2

k

k Sample

s





最小化

(

)

k

kj

j

j

j Asset

s

R

R x









制約

定式化の選択

• Full Covariance

• コンパクト分解

,

ij

i

j

i j Asset

Q x x





最小化

2

k

k Sample

s





最小化

(

)

k

kj

j

j

j Asset

s

R

R x









制約

3000銘柄でも1～20秒程度

半正定値

銘柄数制約付き

平均分散モデル

最小化

制約

t

x Qx

1 ,

j

j Asset

x







0

[ ,

]

j

j

j

j

x



または　

x



l u

( )

supp x



K

( )

{ |

_i

0}

supp x



i x



とりあえず

分枝限定法（

MIQP）

,

{0,1}

j

j

j

j

j

j

l

   



x

u

 



t

x Qx

1 ,

j

j Asset

x







j

j

K







最小化

制約

「場合分け」＋限定

限界あり．二次の目的関数は特に難しい

凸緩和してSDP/SOCPに帰着

目的関数の書きかえとラグランジュ緩和

x



z

Ax



b

1

j

j Asset

x







0

[ ,

]

j

j

j

j

z



または　

z



l u

( )

supp z



K

Az



b

1

j

j Asset

z







(

( ))

( )

t

t

t

x Qx



x Q diag d x





z diag d z

正定値になるよ

うに

d >0を選ぶ

Xiaojin Zheng, Xiaoling Sun, Duan Li,

Convex Relaxations and Mixed-Integer Quadratic

(

( ))

( )

t

t

t

x Qx



x Q diag d x





z diag d z

目的関数の書きかえとラグランジュ緩和

x



z

1

j

j Asset

x







0

[ ,

]

j

j

j

j

z



または　

z



l u

( )

supp z



K

1

j

j Asset

z







緩和（乗数：

_λ）

凸緩和してSDP/SOCPに帰着

min{ ( ; , )} min{ ( ; , )}

_x

_z

x

f x



d



z

f z



d

下界値の最大化問題

最大化

を解く問題⇒SDP， を解く問題⇒SOCP

_{( , )}



_d

( )



0

[ ,

]

j

j

j

j

z



または　

z



l u

( )

supp z



K

1

j

j Asset

z







1

j

j Asset

x







凸緩和してSDP/SOCPに帰着

• SDP/SOCPの定式化

⇒下界値の導出

⇒ MIQPの補強による高速化

上下界値の改善

凸緩和してSDP/SOCPに帰着

エレガントな解決：

StQPに関する事実の応用

最小化

制約

t

x Qx

1 ,

0

j

j

j Asset

x

x









StQP

StQP の最適解は

が強凸であるようなフェース

_の

relative interior に 存在する．

( )

I



t

x Qx

Francesco Cesarone,Andrea Scozzari, and Fabio

Tardella, A new method for mean-variance portfolio

optimization with cardinality constraints

半正定値に

限らない

最小化

制約

t

x Qx

1 ,

0

j

j

j I

Asset

x

x

 







これも

_StQP

が強凸であるようなフェース を

について

調べ尽くせば銘柄数制約付き分散モデルの

大域的最適解を見逃すことはない

t

x Qx



( )

I

,

I

Asset

I

K

 



I



K

エレガントな解決：

StQPに関する事実の応用

正定値な

_{の部分行列の}

インデックスの集合の中に はある

⇒

_{インデックス について}

の部分行列が正定値でなければ，

を含むあらゆる集合は に成り得ない

_{⇒ Increasing Set Algorithm}

Q

0

I

*

I

0

I

*

I

最適解を与える集合

複数のベンチマーク問題に厳密解を

与えることに成功

Q

エレガントな解決：

StQPに関する事実の応用

StQPに関する事実の応用

そういえば．．

• ポートフォリオ最適化問題で

選択される銘柄は一般に少ない

（⇒有効制約法が有利）

ヘッセ行列が

低ランク

なことが理由では？

long/short ポジションを許容する

リスク最小化

最小化

制約

1 ,

0

L

L

j

j

j Asset

x

x









0

0

L

S

j

j

x



または　

x



(

x

L



x

S t

)

Q x

(

L



x

S

)

1 ,

0

S

S

j

j

j Asset

x

x









とりあえず

分枝限定法（

MIQP）

最小化

制約

1 ,

0

L

L

j

j

j Asset

x

x









(

x

L



x

S t

)

Q x

(

L



x

S

)

1 ,

0

S

S

j

j

j Asset

x

x









, (1

)

,

{0,1}

L

S

j

x

j

j

x

j

j















• 連続緩和問題の補強

とりあえず

分枝限定法（

MIQP）

PROBLEM_NAME iqp NUMBER_OF_VARIABLES 600 (#INTEGER/DISCRETE) 200 NUMBER_OF_FUNCTIONS 603 PROBLEM_TYPE MINIMIZATION METHOD ACTIVE_SET_QP <preprocess begin>...<preprocess end>

<iteration begin> .1.2B

up: 1e+050 lo:7.2290069e-020 time: 0.5s:mem(Mb)=77/57:avail(Mb)=4039/1872 llen:0 #prob:0 #piv:173

#1 up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time: 33.4s:mem(Mb)=111/91:avail(Mb)=4004/1837 llen:1286 #prob:2702 #piv:181822

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time: 61.8s:mem(Mb)=139/119:avail(Mb)=3977/1808 llen:2364 #prob:5329 #piv:345892

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:124.7s:mem(Mb)=153/133:avail(Mb)=3963/1793 llen:2904 #prob:7589 #piv:720381

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:184.0s:mem(Mb)=177/157:avail(Mb)=3938/1768 llen:3808 #prob:9569 #piv:1046687

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:187.2s:mem(Mb)=178/157:avail(Mb)=3939/1767 llen:3856 #prob:9689 #piv:1066187

• 制約なしの場合．．

PROBLEM_NAME iqp NUMBER_OF_VARIABLES 600 (#INTEGER/DISCRETE) 200 NUMBER_OF_FUNCTIONS 603 PROBLEM_TYPE MINIMIZATION METHOD ACTIVE_SET_QP <preprocess begin>...<preprocess end>

<iteration begin> .1.2B

up: 1e+050 lo:7.2290069e-020 time: 0.5s:mem(Mb)=77/57:avail(Mb)=4039/1872 llen:0 #prob:0 #piv:173

#1 up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time: 33.4s:mem(Mb)=111/91:avail(Mb)=4004/1837 llen:1286 #prob:2702 #piv:181822

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time: 61.8s:mem(Mb)=139/119:avail(Mb)=3977/1808 llen:2364 #prob:5329 #piv:345892

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:124.7s:mem(Mb)=153/133:avail(Mb)=3963/1793 llen:2904 #prob:7589 #piv:720381

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:184.0s:mem(Mb)=177/157:avail(Mb)=3938/1768 llen:3808 #prob:9569 #piv:1046687

up: 211.13232 lo:7.2290069e-020 gap: 211.13232 time:187.2s:mem(Mb)=178/157:avail(Mb)=3939/1767 llen:3856 #prob:9689 #piv:1066187

• 制約なしの場合．．

分枝限定法は全探索と同じ：

緩和解は相補性を満たさず

目的関数は零

• 相補性条件を緩和した問題

( )

t

L

L

S

S

x

x

M

x



x

 

 

 

 

 

 

M

( )

Q

Q

Q

Q







 





 

_{ }







最小化

制約

L

_j

1

j Asset

x







S

_j

1

j Asset

x







,

0,

L

S

j

j

x

x



j



Asset

が半正定値ならば解ける！

( )

L



( )

M



緩和⇒ＳＤＰを利用した補強の導出

• の最小値を最大化する問題 （ＳＤＰ）

最小化

制約

( )

L





,( )

,( )

,( )

,( )

,( )

,( )

T

L i

L i

L i

S i

j

j

j

S i

S i

j S

Q

Q

x

x

x

x

Q

Q

x

x























_

_

_

_

_

_



_



















1,..,

i



I

( )

0

M





（半正定値）

過去の反復で得られた解

_変数

緩和⇒ＳＤＰを利用した補強の導出

• n=200での実験結果

• λを用いて補強されたMIQP

最小化

制約

L

_{1 ,}

L

₀

j

j

j Asset

x

x









1 ,

0

S

S

j

j

j Asset

x

x









( )

t

L

L

S

S

x

x

M

x



x

 

 

 

 

 

 

, (1

)

,

{0,1}

L

S

j

x

j

j

x

j

j















緩和⇒ＳＤＰを利用した補強の導出

PROBLEM_NAME qp_int NUMBER_OF_VARIABLES 600 (#INTEGER/DISCRETE) 200 NUMBER_OF_FUNCTIONS 603 PROBLEM_TYPE MINIMIZATION METHOD ACTIVE_SET_QP <preprocess begin>...<preprocess end>

<iteration begin> .1...2B

up: 1e+050 lo: 65.069268 time: 1.4s:mem(Mb)=51/32:avail(Mb)=4064/1897 llen:0 #prob:0 #piv:797

up: 1e+050 lo: 65.202836 time: 65.1s:mem(Mb)=59/39:avail(Mb)=4058/1890 llen:1 #prob:349 #piv:136291

up: 1e+050 lo: 65.507298 time:127.3s:mem(Mb)=59/40:avail(Mb)=4055/1888 llen:360 #prob:1119 #piv:357707

up: 1e+050 lo: 65.507298 time:189.1s:mem(Mb)=75/56:avail(Mb)=4041/1874 llen:985 #prob:2369 #piv:649083

#1 up: 78.939007 lo: 65.507298 gap: 13.431709 time:192.7s:mem(Mb)=77/57:avail(Mb)=4039/1872 llen:1034 #prob:2472 #piv:670000

#2 up: 77.947918 lo: 65.507298 gap: 12.44062 time:192.8s:mem(Mb)=77/58:avail(Mb)=4038/1871 llen:1007 #prob:2473 #piv:670194

#3 up: 77.930633 lo: 65.507298 gap: 12.423335 time:193.5s:mem(Mb)=76/57:avail(Mb)=4040/1873 llen:1006 #prob:2498 #piv:674716

#4 up: 77.928705 lo: 65.507298 gap: 12.421407 time:193.8s:mem(Mb)=76/57:avail(Mb)=4039/1872 llen:1006 #prob:2508 #piv:676891

up: 77.928705 lo: 65.507298 gap: 12.421407 time:249.1s:mem(Mb)=88/68:avail(Mb)=4028/1860 llen:1455 #prob:3488 #piv:919435

50銘柄程度ならば

厳密解が得られる場合も．．

Copositive Programmingに帰着

• 等価な定式化

（Copositive Programming の 双対形）

completely positive

matrix

の

cone

の集合

エレガントな解決：

？

偏微分方程式制約付き最適化

PDE constrained optimization

2

(

_i

_i

)

i

 







( , )

0

F

 



最小化

変数

 

_,

_i

_i



制約

離散点での物理量

パラメータ

観測点での値

離散化された

物理方程式

パラメータと物理量を同時に解く

偏微分方程式制約付き最適化

動機づけ

Gassan S Abdoulaev, Kui Ren and Andreas H Hielscher,

Optical tomography as a PDE-constrained optimization problem,

Inverse Problems 21 (2005) 1507–1530

偏微分方程式制約付き最適化

「最適化」の位置づけ

偏微分方程式制約付き最適化

最適化問題としての側面

○ 不等式制約はほぼ効かない

○ ＫＫＴ条件は線形性が強い

× 超大規模（⇒反復法必須）

× ヘッセ行列の正定値性は保障できない

（⇒正規化項の導入）

偏微分方程式制約付き最適化

解くべき一次方程式

t

x

y

b

x

H

A

b

y

A

W



 



_{  }

  



_

_{  }

_

 





_{ }

またお会いしましたね．．

“

_{Saddle Point System”}

正規化項を加えて優対角にはできる

Augumented Lagrangean

のペナルティの逆数

偏微分方程式制約付き最適化

内点法の場合との違い

61

t

x

y

b

x

H

A

b

y

A

W



 



  

  



_

  

_

 





 

超大規模（⇒反復法必須

）

正定値性への目配りが必要に

前処理なし反復法での限界

• 正規化項を手厚くする

• 次の形に変形してＣＧ法

1

t

H



AW A



  

x

b









似たような方に以前お目にかか

ったような．．

計算の安定化（ＣＧ反復減少）

正規化項増大で非物理的解

の構成が鍵

（

_{approximate Schur complement）}

エレガントな解決？

Constraint Preconditioner

1

1

t

t

t

H

A

H

A

A

AH A

S

A

W







 





_

 

_





 



1

,

t

H

H

S

AH A





W

あなたは．．!

S

エレガントな解決？

Approximate

Schur complement

を左辺とする一次方程式解法

•

を対角行列に置き換えてコレスキー分解

• multigrid 法による前処理行列でCG法

• ノルムの定義を変更してCG法

• 対角スケーリングしてCG法を数回適用したものを前処理

行列としてGMRES法

• ピボット選択付きの分解を数ステップ行ったものを前処理

行列としてGMRES法

1

,

t

H

H

S

AH A





W

今度は近似で良いので自由度が高い

まとめ

• 半正定値性は結果を保障の「拠り所」

ー convex QP

ー SDP緩和

• 「距離感」の調整も必要

ー 内点法の実装（線形計画法）

ー cardinality constrained portfolio