LoginSignup
9
9

More than 5 years have passed since last update.

@infhyroyage(Takeru KIDO)

2次錐計画問題(SOCP)から半正定値計画問題(SDP)への帰着

Last updated at Posted at 2018-02-24

任意の2次錐計画問題(Second Order Cone Programming、SOCP)が半正定値計画問題(SemiDefinite Programming、SDP)に属しているという性質、言い換えるとSOCPはSDPの特殊形であるという性質の証明みたいな説明を述べていきます。線形代数と数理最適化の基礎知識を習得している人は理解できると思います。日本語で書き残しているエビデンスが見当たらなかったので、備忘録を兼ねて書き残しました。

扱う最適化問題のモデル

2次錐計画問題(SOCP)

\begin{aligned}
& \text{minimize} & & \mathbf{p}^\top \mathbf{x} \\
& \text{subject to} & & \left\| \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right\| \leq {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \ \ \ (^\forall i = 1 , \ldots , \alpha)
\end{aligned}

各記号、表記の意味は以下の通りとします。

  • $\mathbf{x} \in \mathbb{R}^n$ : 変数
  • $\mathbf{p} \in \mathbb{R}^n$、$\mathbf{A}^{(i)} \in \mathbb{R}^{n_i \times n}$、$\mathbf{b}^{(i)} \in \mathbb{R}^{n_i}$、$\mathbf{c}^{(i)} \in \mathbb{R}^n$、$d^{(i)} \in \mathbb{R}$ : 定数($^\forall i = 1 , \ldots , \alpha$)
  • $|| \cdot ||$ : $n_i(i = 1 , \ldots , \alpha)$次元実数空間上に定義されるユークリッドノルム

半正定値計画問題(SDP)

\begin{aligned}
& \text{minimize} & & \mathbf{Q} \bullet \mathbf{Y} \\
& \text{subject to} & & \mathbf{F}^{(i)} \bullet \mathbf{Y} = g^{(i)} \ \ \ (^\forall i = 1 , \ldots , \beta) \\
& & & \mathbf{Y} \succeq \mathbf{0}_m
\end{aligned}

各記号、表記の意味は以下の通りとします。

  • $\mathbb{S}^m$ : $m$次対称実数行列全体の集合
  • $\mathbf{0}_m$ : $m$次零行列
  • $\mathbf{Y} \in \mathbb{S}^m$ : 変数
  • $\mathbf{Q} \in \mathbb{R}^{m \times m}$、$\mathbf{F}^{(i)} \in \mathbb{R}^{m \times m}$、$g^{(i)} \in \mathbb{R}$ : 定数($^\forall i = 1 , \ldots , \beta$)
  • $\mathbf{F}^{(i)} \bullet \mathbf{Y}$ : $\mathbf{F}^{(i)} , \mathbf{Y}$の$j$行$k$列の要素をそれぞれ$F_{jk}^{(i)} , Y_{jk}$としたとき、$\mathbf{F}^{(i)} \bullet \mathbf{Y} := \sum_{j = 1}^m \sum_{k = 1}^m F_{jk}^{(i)} Y_{jk}$
  • $\mathbf{Y} \succeq \mathbf{0}_m$ : 行列$\mathbf{Y}$は半正定値、すなわち任意の$\mathbf{w} \in \mathbb{R}^m$に対して$\mathbf{w}^\top \mathbf{Yw} \geq 0$が成り立つ

また、上記の双対問題は以下のように表わされます。

\begin{aligned}
& \text{maximize} & & \sum_{i = 1}^\beta g^{(i)} z_i \\
& \text{subject to} & & \mathbf{Q} - \sum_{i = 1}^\beta \mathbf{F}^{(i)} z_i \succeq \mathbf{0}_m
\end{aligned}

ここで、$z_1 , \ldots , z_\beta \in \mathbb{R}$は変数です。

Schur Complementとその性質

以下の$(m + n)$次対称行列$\mathbf{M} \in \mathbb{S}^{m + n}$を考えます。

\mathbf{M} =
\begin{pmatrix}
\mathbf{A} & \mathbf{B} \\
\mathbf{B}^\top & \mathbf{C}
\end{pmatrix}

ここで、$\mathbf{A} \in \mathbb{S}^m , \mathbf{B} \in \mathbb{R}^{m \times n} , \mathbf{C} \in \mathbb{S}^n$であり、$\det A \neq 0$を満たすとします。このとき、行列$\mathbf{C} - \mathbf{B}^\top \mathbf{A}^{-1} \mathbf{B}$を、$\mathbf{M}$における$\mathbf{A}$のSchur Complementと呼びます。実は、$\mathbf{A}$が正定値行列である場合、$\mathbf{M} \succeq \mathbf{0}_{m + n} \Leftrightarrow \mathbf{C} - \mathbf{B}^\top \mathbf{A}^{-1} \mathbf{B} \succeq \mathbf{0}_n$が成り立ちます。この証明は、Boyd&VandenbergheのComvex OptimizationのA.5.5など参考にしてください。上記の性質は凸最適化論の中で非常に重要であり、SDPの半正定値条件式を導き出す際によく利用されています。本記事でも、SOCPのモデルからSchur Complementを利用してSDPのモデルへと帰着していきます。

SOCPの条件式の変形によるSDPへの帰着

さきほど紹介したSOCPのモデルでの条件式である

\left\| \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right\| \leq {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)}

について変形していきます。もし、${\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} = 0$を満たす場合、明らかにこの不等式は両辺とも0となるので、任意の$\mathbf{x}$に対して常に成立してしまい、SOCPの条件式としての意味がありません。そのため、${\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} > 0$の場合のみ考えれば十分です。
SOCPの条件式の両辺を2乗し、以下のように変形します。

\begin{aligned}
\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right)^2 & - \left\| \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right\|^2 & \geq 0 \\
\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right)^2 & - \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right)^\top \mathbf{I}_{n_i} \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right) & \geq 0 \\
\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right) & - \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right)^\top \frac{1}{{\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)}} \mathbf{I}_{n_i} \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right) & \geq 0 \\
\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right) & - \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right)^\top \left( \left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right) \mathbf{I}_{n_i} \right)^{-1} \left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right) & \geq 0 \\
\end{aligned}

ここで、$\mathbf{I}_{n_i}$は$n_i$次単位行列です。そのため、$\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right) \mathbf{I}_{n_i}$は正定値行列です。したがって、Schur Complementの性質からSOCPのモデルでの条件式は

\begin{pmatrix}
\left( {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)} \right) \mathbf{I}_{n_i} & \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \\
\left( \mathbf{A}^{(i)} \mathbf{x} + \mathbf{b}^{(i)} \right)^\top & {\mathbf{c}^{(i)}}^\top \mathbf{x} + d^{(i)}
\end{pmatrix}
\succeq \mathbf{0}_{n_i + 1}

となります。さらに、各係数ベクトルや行列$\mathbf{A}^{(i)}$の第$j$列の要素を、ベクトル$\mathbf{a}_j^{(i)} \in \mathbb{R}^{n_i}$を用いて

\mathbf{A}^{(i)} :=
\begin{pmatrix}
\mathbf{a}_1^{(i)} & \cdots & \mathbf{a}_n^{(i)}
\end{pmatrix}
, \mathbf{c}^{(i)} :=
\begin{pmatrix}
c_1^{(i)} \\
\vdots \\
c_n^{(i)}
\end{pmatrix}
, \mathbf{p} :=
\begin{pmatrix}
p_1 \\
\vdots \\
p_n
\end{pmatrix}
, \mathbf{x} :=
\begin{pmatrix}
x_1 \\
\vdots \\
x_n
\end{pmatrix}

のように各列ごとに表わすと、上記の関係式は以下のように変形できます。

\begin{aligned}
\begin{pmatrix}
\left( \sum_{j = 1}^n c_j^{(i)} x_j + d^{(i)} \right) \mathbf{I}_{n_i} & \sum_{j = 1}^n \mathbf{a}_j^{(i)} x_j + \mathbf{b}^{(i)} \\
\sum_{j = 1}^n \left( \mathbf{a}_j^{(i)}  \right)^\top x_j + \left( \mathbf{b}^{(i)} \right)^\top & \sum_{j = 1}^n c_j^{(i)} x_j + d^{(i)}
\end{pmatrix}
& \succeq \mathbf{0}_{n_i + 1} \\
\begin{pmatrix}
d^{(i)} \mathbf{I}_{n_i} & \mathbf{b}^{(i)} \\
\left( \mathbf{b}^{(i)} \right)^\top & d^{(i)}
\end{pmatrix}
- \sum_{j = 1}^n
\begin{pmatrix}
- c_j^{(i)} \mathbf{I}_{n_i} & - \mathbf{a}_j^{(i)} \\
- \left( \mathbf{a}_j^{(i)} \right)^\top & - c_j^{(i)}
\end{pmatrix}
x_j
& \succeq \mathbf{0}_{n_i + 1}
\end{aligned}

以上より、下側の半正定値条件式がSDPの双対問題の条件式の形に一致している…ように見えますが、この半正定値条件式は任意の$i = 1 , \ldots , \alpha$について成立していることに注意してください。要するに、下側の半正定値条件式は$\alpha$個も存在しています。一方、SDPの双対問題の条件式は1個のみです。そのため、下側の半正定値条件式を1個に圧縮する必要があります。ではどうすれば良いのかと言いますと、各$i = 1 , \ldots , \alpha$に対して以下のようにBlockwiseな対角行列を考えれば良いのです。

\begin{pmatrix}
d^{(1)} \mathbf{I}_{n_1} & \mathbf{b}^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\left( \mathbf{b}^{(1)} \right)^\top & d^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\vdots & \vdots & \ddots \vdots & \vdots & \\
\mathbf{0} & \mathbf{0} & \cdots & d^{(\alpha)} \mathbf{I}_{n_\alpha} & \mathbf{b}^{(\alpha)} \\
\mathbf{0} & \mathbf{0} & \cdots & \left( \mathbf{b}^{(\alpha)} \right)^\top & d^{(\alpha)}
\end{pmatrix}
- \sum_{j = 1}^n
\begin{pmatrix}
- c_j^{(1)} \mathbf{I}_{n_1} & - \mathbf{a}_j^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
- \left( \mathbf{a}_j^{(1)} \right)^\top & - c_j^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & - c_j^{(\alpha)} \mathbf{I}_{n_\alpha} & - \mathbf{a}_j^{(\alpha)} \\
\mathbf{0} & \mathbf{0} & \cdots & - \left( \mathbf{a}_j^{(\alpha)} \right)^\top & - c_j^{(\alpha)}
\end{pmatrix}
x_j
\succeq \mathbf{0}_\gamma

ここで、$\gamma := \sum_{i = 1}^\alpha n_i + \alpha$です。なぜならば、

\mathbf{U}_1 \succeq \mathbf{0}_{n_1} , \ldots , \mathbf{U}_\alpha \succeq \mathbf{0}_{n_\alpha}
\Leftrightarrow
\begin{pmatrix}
\mathbf{U}_1 & \cdots & \mathbf{0} \\
\vdots & \ddots & \vdots \\
\mathbf{0} & \cdots & \mathbf{U}_\alpha
\end{pmatrix}
\succeq \mathbf{0}_{\sum_{i = 1}^\alpha n_i}

の性質が、半正定値行列の定義より成り立つからです。興味のある人は証明してみてください。以上より、この半正定値条件式とSDPの双対問題の条件式を照らし合わせることによって、

\begin{aligned}
\beta & \gets n \\
g^{(i)} & \gets p_i \ \ \ (^\forall i = 1 , \ldots , \beta) \\
\mathbf{Q} & \gets
\begin{pmatrix}
d^{(1)} \mathbf{I}_{n_1} & \mathbf{b}^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\left( \mathbf{b}^{(1)} \right)^\top & d^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\vdots & \vdots & \ddots \vdots & \vdots & \\
\mathbf{0} & \mathbf{0} & \cdots & d^{(\alpha)} \mathbf{I}_{n_\alpha} & \mathbf{b}^{(\alpha)} \\
\mathbf{0} & \mathbf{0} & \cdots & \left( \mathbf{b}^{(\alpha)} \right)^\top & d^{(\alpha)}
\end{pmatrix} \\
\mathbf{F}^{(i)} & \gets
\begin{pmatrix}
- c_i^{(1)} \mathbf{I}_{n_1} & - \mathbf{a}_i^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
- \left( \mathbf{a}_i^{(1)} \right)^\top & - c_i^{(1)} & \cdots & \mathbf{0} & \mathbf{0} \\
\vdots & \vdots & \ddots & \vdots & \vdots \\
\mathbf{0} & \mathbf{0} & \cdots & - c_i^{(\alpha)} \mathbf{I}_{n_\alpha} & - \mathbf{a}_i^{(\alpha)} \\
\mathbf{0} & \mathbf{0} & \cdots & - \left( \mathbf{a}_i^{(\alpha)} \right)^\top & - c_i^{(\alpha)}
\end{pmatrix}
 \ \ \ (^\forall i = 1 , \ldots , \beta)
\end{aligned}

と係数を用意することで、任意のSOCPがSDPに帰着できることがわかりました。これによって、SDPAといったSDP専用のソルバーを使ってSOCPを解くことができます！

例題

1つのSOCPの例題を考え、それをSDPに帰着します。また、本当に正しく変形できているのかSOCP、SDPのソルバーを用いて検証してみます。正しく変形できているならば、両者の最適値は一致しているはずです。本記事ではCVXOPTのソルバーを扱い、SOCPの例題はそのUser's Guideに記載されているSecond-Order Cone Programmingから抜粋します。

SOCPのまま解く場合

以下に、その例題を再掲します。

\begin{aligned}
& \text{minimize} & & -2 x_1 + x_2 +5 x_3 \\
& \text{subject to} & & \left\|
\begin{pmatrix}
-13 x_1 +3 x_2 +5 x_3 -3 \\
-12 x_1 +12 x_2 -6 x_3 -2
\end{pmatrix}
\right\| \leq -12 x_1 -6 x_2 +5 x_3 -12 \\
& & & \left\|
\begin{pmatrix}
-3 x_1 +6 x_2 +2 x_3 \\
x_1 +9 x_2 +2 x_3 +3 \\
- x_1 -19 x_2 +3 x_3 -42
\end{pmatrix}
\right\| \leq -3 x_1 +6 x_2 -10 x_3 +27
\end{aligned}

これより、係数は以下のようになります。

\begin{aligned}
\mathbf{p} =
\begin{pmatrix}
-2 \\
1 \\
5
\end{pmatrix}
, \\
\mathbf{a}_1^{(1)} =
\begin{pmatrix}
-13 \\
-12
\end{pmatrix}
, \mathbf{a}_2^{(1)} =
\begin{pmatrix}
3 \\
12
\end{pmatrix}
, \mathbf{a}_3^{(1)} =
\begin{pmatrix}
5 \\
-6
\end{pmatrix}
, \\
\mathbf{a}_1^{(2)} =
\begin{pmatrix}
-3 \\
1 \\
-1
\end{pmatrix}
, \mathbf{a}_2^{(2)} =
\begin{pmatrix}
6 \\
9 \\
-19
\end{pmatrix}
, \mathbf{a}_3^{(2)} =
\begin{pmatrix}
2 \\
2 \\
3
\end{pmatrix}
, \\
\mathbf{b}^{(1)} =
\begin{pmatrix}
-3 \\
-2
\end{pmatrix}
, \mathbf{b}^{(2)} =
\begin{pmatrix}
0 \\
3 \\
-42
\end{pmatrix}
, \\
\mathbf{c}^{(1)} =
\begin{pmatrix}
-12 \\
-6 \\
5
\end{pmatrix}
, \mathbf{c}^{(2)} =
\begin{pmatrix}
-3 \\
6 \\
-10
\end{pmatrix}
, \\
d^{(1)} = -12 , d^{(2)} = 27
\end{aligned}

これらの係数を用いて、CVXOPTのcvxopt.solvers.socp関数でSOCPを解いてみます。User's Guideに記載されている実装とほぼ同じです。

>>> from cvxopt import matrix, solvers
>>> p = matrix([-2., 1., 5.])
>>> Gq = [ matrix( [[12., 13., 12.], [6., -3., -12.], [-5., -5., 6.]] ) ]
>>> Gq += [ matrix( [[3., 3., -1., 1.], [-6., -6., -9., 19.], [10., -2., -2., -3.]] ) ]
>>> hq = [ matrix( [-12., -3., -2.] ),  matrix( [27., 0., 3., -42.] ) ]
>>> sol = solvers.socp(p, Gq=Gq, hq=hq)
     pcost       dcost       gap    pres   dres   k/t
 0:  4.9969e+00 -1.7285e+01  6e+01  3e-01  4e+00  1e+00
 1: -1.6732e+00 -7.0431e+00  1e+01  7e-02  1e+00  6e-01
 2: -1.6221e+01 -3.5417e+01  2e+02  3e-01  5e+00  7e+00
 3: -2.1832e+01 -2.2849e+01  3e+01  4e-02  6e-01  2e+00
 4: -3.5265e+01 -3.5594e+01  1e+01  1e-02  2e-01  9e-01
 5: -3.8303e+01 -3.8314e+01  3e-01  4e-04  6e-03  2e-02
 6: -3.8342e+01 -3.8342e+01  1e-02  1e-05  2e-04  7e-04
 7: -3.8346e+01 -3.8346e+01  9e-04  1e-06  2e-05  7e-05
 8: -3.8346e+01 -3.8346e+01  4e-05  6e-08  9e-07  4e-06
 9: -3.8346e+01 -3.8346e+01  2e-06  3e-09  4e-08  2e-07
Optimal solution found.
>>> print(sol['x'])
[-5.01e+00]
[-5.77e+00]
[-8.52e+00]

>>> print(p.trans() * sol['x'])
[-3.83e+01]

以上より、最適値は-3.83e+01と求まりました。

SDPに帰着して解く場合

まず上式で述べたように、SOCPでの係数からSDP用の係数$\mathbf{Q} , \mathbf{F}^{(i)} , g^{(i)} ( ^\forall i = 1 , \ldots , \beta )$を用意します。

\begin{aligned}
\mathbf{Q} =
\begin{pmatrix}
-12 & 0 & -3 & 0 & 0 & 0 & 0 \\
0 & -12 & -2 & 0 & 0 & 0 & 0 \\
-3 & -2 & -12 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 27 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 27 & 0 & 3 \\
0 & 0 & 0 & 0 & 0 & 27 & -42 \\
0 & 0 & 0 & 0 & 3 & -42 & 27
\end{pmatrix}
, \\
\mathbf{F}^{(1)} = 
\begin{pmatrix}
12 & 0 & 13 & 0 & 0 & 0 & 0 \\
0 & 12 & 12 & 0 & 0 & 0 & 0 \\
13 & 12 & 12 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 3 & 0 & 0 & 3 \\
0 & 0 & 0 & 0 & 3 & 0 & -1 \\
0 & 0 & 0 & 0 & 0 & 3 & 1 \\
0 & 0 & 0 & 3 & -1 & 1 & 3
\end{pmatrix}
, \\
\mathbf{F}^{(2)} =
\begin{pmatrix}
6 & 0 & -3 & 0 & 0 & 0 & 0 \\
0 & 6 & -12 & 0 & 0 & 0 & 0 \\
-3 & -12 & 6 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & -6 & 0 & 0 & -6 \\
0 & 0 & 0 & 0 & -6 & 0 & -9 \\
0 & 0 & 0 & 0 & 0 & -6 & 19 \\
0 & 0 & 0 & -6 & -9 & 19 & -6
\end{pmatrix}
, \\
\mathbf{F}^{(3)} =
\begin{pmatrix}
-5 & 0 & -5 & 0 & 0 & 0 & 0 \\
0 & -5 & 6 & 0 & 0 & 0 & 0 \\
-5 & 6 & -5 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 10 & 0 & 0 & -2 \\
0 & 0 & 0 & 0 & 10 & 0 & -2 \\
0 & 0 & 0 & 0 & 0 & 10 & -3 \\
0 & 0 & 0 & -2 & -2 & -3 & 10
\end{pmatrix}
, \\
g^{(1)} = -2 , g^{(1)} = 1 , g^{(1)} = 5
\end{aligned}

これらの係数を用いて、CVXOPTのcvxopt.solvers.sdp関数でSDPを解いてみます。行列のコーディングは疎性を利用してcvxopt.spmatrixを使うと、SDPの求解性能がより良くなるかもしれませんが、本記事では変数の数が少なく面倒くさいので目を瞑ります。

>>> from cvxopt import matrix, solvers
>>> from itertools import chain
>>> Q = [[-12,   0,  -3,  0,  0,   0,   0],
...      [  0, -12,  -2,  0,  0,   0,   0],
...      [ -3,  -2, -12,  0,  0,   0,   0],
...      [  0,   0,   0, 27,  0,   0,   0],
...      [  0,   0,   0,  0, 27,   0,   3],
...      [  0,   0,   0,  0,  0,  27, -42],
...      [  0,   0,   0,  0,  3, -42,  27]]
>>> F1 = [[12,  0, 13,  0,  0,  0,  0],
...       [ 0, 12, 12,  0,  0,  0,  0],
...       [13, 12, 12,  0,  0,  0,  0],
...       [ 0,  0,  0,  3,  0,  0,  3],
...       [ 0,  0,  0,  0,  3,  0, -1],
...       [ 0,  0,  0,  0,  0,  3,  1],
...       [ 0,  0,  0,  3, -1,  1,  3]]
>>> F2 = [[ 6,   0,  -3,  0,  0,  0,  0],
...       [ 0,   6, -12,  0,  0,  0,  0],
...       [-3, -12,   6,  0,  0,  0,  0],
...       [ 0,   0,   0, -6,  0,  0, -6],
...       [ 0,   0,   0,  0, -6,  0, -9],
...       [ 0,   0,   0,  0,  0, -6, 19],
...       [ 0,   0,   0, -6, -9, 19, -6]]
>>> F3 = [[-5,  0, -5,  0,  0,  0,  0],
...       [ 0, -5,  6,  0,  0,  0,  0],
...       [-5,  6, -5,  0,  0,  0,  0],
...       [ 0,  0,  0, 10,  0,  0, -2],
...       [ 0,  0,  0,  0, 10,  0, -2],
...       [ 0,  0,  0,  0,  0, 10, -3],
...       [ 0,  0,  0, -2, -2, -3, 10]]
>>> p = matrix([-2., 1., 5.])
>>> Gs = [matrix([list(chain.from_iterable(F1)), list(chain.from_iterable(F2)), list(chain.from_iterable(F3))], tc='d')]
>>> hs = [matrix(Q, tc='d')]
>>> sol = solvers.sdp(p, Gs=Gs, hs=hs)
     pcost       dcost       gap    pres   dres   k/t
 0:  6.5618e+00 -7.8254e+01  2e+02  4e-01  1e+01  1e+00
 1:  6.4655e-01 -2.3215e+01  5e+01  1e-01  3e+00  1e+00
 2: -5.8497e+00 -1.7747e+01  3e+01  7e-02  2e+00  1e+00
 3: -2.7469e+01 -4.1282e+01  1e+02  1e-01  3e+00  5e+00
 4: -3.7717e+01 -3.7743e+01  3e+00  3e-03  7e-02  5e-01
 5: -3.8338e+01 -3.8342e+01  8e-02  6e-05  2e-03  8e-03
 6: -3.8346e+01 -3.8346e+01  3e-03  2e-06  6e-05  3e-04
 7: -3.8346e+01 -3.8346e+01  1e-04  9e-08  2e-06  1e-05
 8: -3.8346e+01 -3.8346e+01  2e-06  2e-09  4e-08  2e-07
Optimal solution found.
>>> print(sol['x'])
[-5.01e+00]
[-5.77e+00]
[-8.52e+00]

>>> print(p.trans() * sol['x'])
[-3.83e+01]

以上より、最適値は-3.83e+01と求まり、SOCPの場合での最適値と一致しました！SOCPの問題がSDPに帰着されていることが検証できましたね。

注意点

任意のSOCPがSDPの形に帰着できることがわかったので、今度からSOCP用のソルバーを捨ててSDP用のソルバーのみ使っていけばいいじゃん…残念ながら、そんなうまい話は通用しません。
というのも、SOCPのモデルを無理矢理SDPのモデルに落としこんでいるため、変数の数が膨大になってしまうからです。最初で紹介したモデルだと、SOCPの変数は$\mathbf{x} \in \mathbb{R}^n$であり$n$個のみでしたが、SOCPの変数は$\mathbf{Y} \in \mathbb{S}^\gamma$のため、$\gamma ( \gamma + 1 ) / 2 = O \left( \left( \sum_{i = 1}^\alpha n_i + \alpha \right)^2 \right)$個も用意する必要があります。そのため。$n_i , \alpha$の値によってはSDPの変数の個数が非常に多くなるため、SOCPを変形してSDPにしてからSDP用のソルバーを用いると、計算時間がかえって悪化する場合があります。
一概には言えませんが、SOCPを解きたい場合において、使いたいライブラリにSOCP用のソルバーがあれば、それをそのまま使った方が無難だと思います。

9
9
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
9
9