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) - 夜間飛行
参考
- ECOS: An SOCP Solver for Embedded Systems
- Decison Tree for Optimization Software
- Second Order Cone Programming with CVXOPT | pwnetics
No comments:
Post a Comment