PICOSによる錐計画

PICOSとは?

PICOSは,錐計画のモデリングのためのインターフェースです.錐計画は線形計画問題を拡張したものと考えることができ,二次錐計画問題(Second-Order Cone Programming, SOCP)や,半正定値計画問題(SemiDefinite Programming,SDP)が錐計画の例です.

LPをPICOSで書く

制約式を1つずつ書く

PICOSのtutorialにある,次のLPを解いてみましょう.

In [66]:
from IPython.display import Image
fig=Image(url='http://mathopt.sakura.ne.jp/LP1.png')
fig
Out[66]:
In [67]:
import picos as pic
P=pic.Problem()
x=P.add_variable('x',2)
P.add_constraint(x[0]>x[1])
P.add_constraint(x[0]<=3)
P.add_constraint(x[0]+x[1]<=4)
objective=0.5*x[0]+x[1]
P.set_objective('max',objective)
print P
P.solve(solver='cvxopt',verbose=False);
print P.obj_value()
---------------------
optimization problem  (LP):
2 variables, 3 affine constraints

x 	: (2, 1), continuous

	maximize 0.5*x[0] + x[1]
such that
  x[0] > x[1]
  x[0] < 3.0
  x[0] + x[1] < 4.0
---------------------
3.0000000002

制約式を行列で書く

上の定義では,制約式をひとつずつ定義しましたが,おなじことを行列とベクトルを使って行うことができます.行列とベクトルを用いて書くとこのLPは次のようになります.

In [68]:
fig2=Image(url='http://mathopt.sakura.ne.jp/LP2.png')
fig2
Out[68]:

このためには,制約行列A, ベクトルb, そして目的関数を定めるベクトルcを定義します.Aの定義にはcvxoptの行列生成機能を使います.cvx.matrix()の引数には,行列の列を順に指定します.

In [69]:
import cvxopt as cvx
A=pic.new_param('A',cvx.matrix([[-1,1,1],[1,0,1]]))
b=[0,3,4]
c=[0.5,1]

P=pic.Problem()
x=P.add_variable('x',2)
objective=(c|x)
P.set_objective('max',objective)
P.add_constraint(A*x<=b)
print P
P.solve(solver='cvxopt',verbose=False);
print P.obj_value()
---------------------
optimization problem  (LP):
2 variables, 3 affine constraints

x 	: (2, 1), continuous

	maximize 〈 [ 2 x 1 MAT ] | x 〉
such that
  A*x < [ 3 x 1 MAT ]
---------------------
3.0000000002

ここでは,ベクトルcとベクトルxの内積を定義するのに,(c|x)という表現を用いています.picosでは,適宜オーバーロードされており,直感的な記述が可能です.制約式の定義では,A*xという表現を使っています.これは,行列Aとベクトルxの積を定義しています.

SDPをPICOSで書く

さて,これをもとに,簡単な半正定値計画問題を定義してみます.半正定値計画問題は,数式では次のように表わせます.

In [70]:
SDP1=Image(url='http://mathopt.sakura.ne.jp/SDP1.png')
SDP1
Out[70]:

ここで,2番目の制約式は,行列Yが半正定値であることを表わします.LPは半正定値計画問題として書くことができます.上の半正定値計画問題にあわせて,先のLPを書き直してみます.

In [71]:
SDP2=Image(url='http://mathopt.sakura.ne.jp/SDP2.png')
SDP2
Out[71]:

さらに書き換えます.

In [72]:
SDP3=Image(url='http://mathopt.sakura.ne.jp/SDP3.png')
SDP3
Out[72]:

この各ベクトルを(3,3)行列として書き換えます.

In [73]:
SDP4=Image(url='http://mathopt.sakura.ne.jp/SDP4.png')
SDP4
Out[73]:

ということで,

In [74]:
SDP5=Image(url='http://mathopt.sakura.ne.jp/SDP5.png')
SDP5
Out[74]:

とすれば,上記の形式の半正定値計画問題として表すことができます.

さて,この半正定値計画問題をpicosで定義してみましょう.

In [75]:
A={}
A[0]=pic.new_param('A0',cvx.matrix([[0,0,0],[0,3,0],[0,0,4]]))
A[1]=pic.new_param('A1',cvx.matrix([[-1,0,0],[0,1,0],[0,0,1]]))
A[2]=pic.new_param('A2',cvx.matrix([[1,0,0],[0,0,0],[0,0,1]]))
b=[0.5,1.0]
In [76]:
sdp=pic.Problem()
z=sdp.add_variable('z',2)
objective=(b|z)
sdp.set_objective('max',objective)
sdp.add_constraint(A[0]-z[0]*A[1]-z[1]*A[2]>>0)
print sdp
---------------------
optimization problem  (SDP):
2 variables, 0 affine constraints, 6 vars in 1 SD cones

z 	: (2, 1), continuous

	maximize 〈 [ 2 x 1 MAT ] | z 〉
such that
  A0 -z[0]*A1 -z[1]*A2 ≽ |0|
---------------------

ここでは,>>という二項演算子を用いています.これは,X>>0は行列Xが半正定値であることを表わすようにオーバーロードされています.したがって,上のようにadd_constraint()を用いて半正定値計画問題を定義することができます.

これを解くと,LPとして定義した時と同じ最適値が得られます.

In [77]:
sdp.solve(solver='cvxopt',verbose=False);
print sdp.obj_value()
3.0000000002

SOCPをPICOSで書く

今度は,二次錐計画問題(socp)を定義してみます.二次錐計画問題は,数式では次のように表わせます.

In [78]:
SOCP1=Image(url='http://mathopt.sakura.ne.jp/SOCP1.png')
SOCP1
Out[78]:

ここで,a_iとyは適当な次元のベクトルとします.また,2番目の制約式は,ベクトルy=(y_0,y_1,y_2,...,y_n)が2次錐に入っている,すなわち,y_0 >= ¥sqrt{ y_1^2 + y_2^2 + ... + y_n^2}を満たすことを意味します.

In [79]:
import math
import random
a={}
a[0]=[random.randint(1,5)*10 for i in range(1,4)]
a[1]=[random.randint(1,5)*10 for i in range(1,4)]
a[2]=[random.randint(1,5)*10 for i in range(1,4)]
b=[0.5,1]
In [80]:
socp=pic.Problem()
z=socp.add_variable('z',2)
objective=(b|z)
socp.set_objective('max',objective)
socp.add_constraint(abs((a[0]-z[0]*a[1]-z[1]*a[2])[1:] )< (a[0]-z[0]*a[1]-z[1]*a[2])[0])
print socp
---------------------
optimization problem  (SOCP):
2 variables, 0 affine constraints, 3 vars in 1 SO cones

z 	: (2, 1), continuous

	maximize 〈 [ 2 x 1 MAT ] | z 〉
such that
  ||( -z[0]*[ 3 x 1 MAT ] + [ 3 x 1 MAT ] -z[1]*[ 3 x 1 MAT ] )[1:]|| < ( -z[0]*[ 3 x 1 MAT ] + [ 3 x 1 MAT ] -z[1]*[ 3 x 1 MAT ] )[0]
---------------------
In [81]:
socp.solve(solver='cvxopt');
print socp.obj_value()
--------------------------
  cvxopt CONELP solver
--------------------------
     pcost       dcost       gap    pres   dres   k/t
 0: -5.8955e-01 -5.8955e-01  7e-01  4e-01  3e-16  1e+00
 1: -1.0115e-01 -1.0507e-01  7e-02  4e-02  4e-16  1e-01
 2: -1.0089e-01 -1.0120e-01  2e-03  1e-03  2e-16  3e-03
 3: -1.0031e-01 -1.0031e-01  3e-05  1e-05  2e-16  4e-05
 4: -1.0030e-01 -1.0031e-01  7e-07  4e-07  4e-15  9e-07
 5: -1.0030e-01 -1.0030e-01  2e-08  9e-09  5e-13  2e-08
 6: -1.0030e-01 -1.0030e-01  2e-10  1e-10  4e-12  3e-10
Optimal solution found.
cvxopt status: optimal
0.100304937556