PICOSは,錐計画のモデリングのためのインターフェースです.錐計画は線形計画問題を拡張したものと考えることができ,二次錐計画問題(Second-Order Cone Programming, SOCP)や,半正定値計画問題(SemiDefinite Programming,SDP)が錐計画の例です.
PICOSのtutorialにある,次のLPを解いてみましょう.
from IPython.display import Image
fig=Image(url='http://mathopt.sakura.ne.jp/LP1.png')
fig
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は次のようになります.
fig2=Image(url='http://mathopt.sakura.ne.jp/LP2.png')
fig2
このためには,制約行列A, ベクトルb, そして目的関数を定めるベクトルcを定義します.Aの定義にはcvxoptの行列生成機能を使います.cvx.matrix()の引数には,行列の列を順に指定します.
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の積を定義しています.
さて,これをもとに,簡単な半正定値計画問題を定義してみます.半正定値計画問題は,数式では次のように表わせます.
SDP1=Image(url='http://mathopt.sakura.ne.jp/SDP1.png')
SDP1
ここで,2番目の制約式は,行列Yが半正定値であることを表わします.LPは半正定値計画問題として書くことができます.上の半正定値計画問題にあわせて,先のLPを書き直してみます.
SDP2=Image(url='http://mathopt.sakura.ne.jp/SDP2.png')
SDP2
さらに書き換えます.
SDP3=Image(url='http://mathopt.sakura.ne.jp/SDP3.png')
SDP3
この各ベクトルを(3,3)行列として書き換えます.
SDP4=Image(url='http://mathopt.sakura.ne.jp/SDP4.png')
SDP4
ということで,
SDP5=Image(url='http://mathopt.sakura.ne.jp/SDP5.png')
SDP5
とすれば,上記の形式の半正定値計画問題として表すことができます.
さて,この半正定値計画問題をpicosで定義してみましょう.
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]
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として定義した時と同じ最適値が得られます.
sdp.solve(solver='cvxopt',verbose=False);
print sdp.obj_value()
3.0000000002
今度は,二次錐計画問題(socp)を定義してみます.二次錐計画問題は,数式では次のように表わせます.
SOCP1=Image(url='http://mathopt.sakura.ne.jp/SOCP1.png')
SOCP1
ここで,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}を満たすことを意味します.
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]
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] ---------------------
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