SOCPのソルバーでQPを解く話

PythonRoboticsの一例でQPのソルバーとしてECOSが使われていました。ECOSは二次錐計画問題 (SOCP)のソルバーです。 QP ⊂ QCQP(Quadratically constrained quadratic program)⊂ SOCP であって、
制約なしのQP形式の問題

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ \end{array} \]

は、変数tを追加して

\[ \begin{array}{ll} \min_{x,t} & t \\ {\rm s.t.} & \frac{1}{2} x^T Q x + q^T x + r \leq t \\ \end{array} \]

と変換すればSOCPとして解ける(制約ありも同様)とのことでしたので、ちょっと実装してみます。

まず改めて、$ x \in \mathbb{R}^n $ のとき、

\[ \begin{aligned} y &\equiv \left[ x, t\right] \\ c &\equiv [0_n,1] \\ Q &= U^T U \\ Q_s &\equiv \begin{bmatrix} U & 0 \end{bmatrix} \\ q_s &\equiv [q,0] \\ \end{aligned} \]

とすると表現がこう変わります

\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ {\rm s.t.} & \frac{1}{2} y^T Q_s^T Q_s y + q_s^T y + r \leq c^T \cdot y \\ \end{array} \]

制約式をノルムの形式に変換すると

\[ \left\| \begin{matrix} \left( 1 + (q_s+c)^{T}y+r \right) / \sqrt{2} \\ Q_s y \end{matrix} \right\|_{2} \leq \frac{1-(q_s+c)^{T}y-r}{\sqrt{2}} \]

最後にCVXOPTなどでも使われるスラック変数を用いた表現に変形します。

\[ \begin{bmatrix} - (q_s + c)/\sqrt{2} \\ (q_s + c)/\sqrt{2} \\ Q_s \end{bmatrix} \cdot y + \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} = \left[ \begin{array}{c} s_0 \\ s_1 \\s_2 \end{array} \right], s_0 \geq \left\Vert \begin{array}{c} s_1 \\ s_2 \end{array} \right\Vert_2, s_0, s_1 \in \mathbb{R}, s_2 \in \mathbb{R}^{n} \]

続けて $ G \cdot x \preceq_K h $ という形式にあわせて、

\[ \begin{bmatrix} (q_s + c)/\sqrt{2} \\ - (q_s + c)/\sqrt{2} \\ - Q_s \end{bmatrix} \cdot y + \begin{bmatrix} s_0 \\ s_1 \\ s_2 \end{bmatrix} = \begin{bmatrix} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0 \end{bmatrix} \]

QPをSOCPとして解く

くどいですが、以上をまとめると

\[ \begin{array}{ll} \min_{x} & \frac{1}{2} x^T Q x + q^T x + r \\ where & x \in \mathbb{R}^n \end{array} \]
というQPは、
\[ \begin{array}{ll} \min_{y} & c^T \cdot y \\ \text{s.t.} & G \cdot y \preceq_{\mathcal{K}} h \quad \left( \Leftrightarrow G \cdot y + s = h, s \in \mathcal{K} \right) \\ \text{where} & y \equiv [x\ t], t \in \mathbb{R} \\ & c = [0_{1 \times n}\ 1] \\ & Q = U^T U \\ & \mathcal{K} \equiv \{ (u_0, u_1) \in \mathbb{R} \times \mathbb{R}^{n+1} \ |\ u_0 \geq \Vert u_1 \Vert_2 \} \\ & G = \left[ \begin{array}{cc} q/\sqrt{2} & 1/\sqrt{2} \\ -q/\sqrt{2} & -1/\sqrt{2} \\ U & 0 \end{array} \right], h = \left[ \begin{array}{cc} (1-r)/\sqrt{2} \\ (1+r)/\sqrt{2} \\ 0_{n \times 1} \end{array} \right] \\ \end{array} \]
というSOCPとして解くことができます。

余談

ECOSに投げるときは、Usage from C · embotech/ecos Wikiを参照しつつ

// Nは元となるQPのxの次元
idxint l = 0; // もし制約付QPを解く場合は l>0, 配列Gに対してはconeの定義にあたる係数よりも前方に配置する
idxint ncones = 1;
idxint q[] = {N+2};

とします。

補足

$ Q = U^T U $ は、コレスキー分解です。

コレスキー分解でベンチマークを取った先日のポストも広告:Pythonと関数単位での高速化(cython, numba, swig) - 夜間飛行

参考

No comments:

Post a Comment