Self-Dual Embeddingを使った線形計画法に対する内点法

この記事はJij Advent Calendar 2024の7日目の記事です．

この記事では，線形計画問題を自己双対線形計画問題へとマップして内点法で解く手法を紹介する．
オリジナルな仕事はYe, Todd and Mizuno[1]であり，水野先生が公開しているテキスト[2]でもその内容を読むことができる．私は[3]でこの手法を知ったので，[3]をベースに説明したいと思う

この手法は，人工的に性質の良い問題(自己双対線形計画問題)を作り，それを内点法で解く手法である．まずは，この性質の良い人工的な問題を作るところを説明し，つづいて，それらの性質を利用して，内点法を実装する部分を説明することにする．

自己双対線形計画問題の導出

$n$個の変数を持ち，制約条件が$m$本ある線形計画問題を考える．

$$(P)\quad \min_{\bm{x}} \{ \bm{c} \cdot \bm{x} : A\bm{x} \geq \bm{b},~\bm{x} \geq 0 \}$$

ここで$A \in \R^{m \times n}$，$b \in \R^{m}$，$c \in \R^{n}$である．どんな線形計画問題でもこの形に変形することができる．

この問題に対する双対問題は，

$$(D)\quad \max_{\bm{y}}\{\bm{c}\cdot \bm{y} : A^{\top} \bm{y} \leq \bm{c},~\bm{y} \geq 0 \}$$

と書くことができる．

(P)-(D)の最適解は，双対ギャップ$\bm{c} \cdot\bm{x} - \bm{b}\cdot \bm{y}$を0にするので，以下の不等式を満たす解が存在するときにのみ最適解が存在する．

$$\tag{1} \begin{cases} A\bm{x} &\geq \bm{b}, \quad \bm{x} \geq 0, \\ -A^{\top} \bm{y} &\geq -\bm{c}, \quad \bm{y} \geq 0, \\ \bm{b} \cdot \bm{y} - \bm{c} \cdot \bm{x} &\geq 0, \end{cases}$$

上の式で，最初の2行は(P),(D)の実行可能解であることを示している．
(P),(D)の実行可能解に対して，弱双対性から $\bm{b} \cdot \bm{y} \leq \bm{c} \cdot \bm{x}$ が成立するので，3番目の制約を満たためには $\bm{b} \cdot \bm{y} = \bm{c} \cdot \bm{x}$でなければいけないことがわかり，双対ギャップが閉じた解の存在を要求している．

ここで，homogenizing variable $\beta$を導入する．$\beta = 1$の時，問題(1)は以下の問題と等価である．

$$\tag{2} \bar{M} \coloneqq \begin{bmatrix} 0_{m \times m} & A & - \bm{b} \\ -A^{\top} & 0_{n \times n} & \bm{c} \\ \bm{b}^{\top} & - \bm{c}^{\top} & 0 \end{bmatrix}, \quad \bar{\bm{z}}\coloneqq \begin{bmatrix} \bm{y} \\ \bm{x} \\ \beta \end{bmatrix},\quad \bar{M}\bar{\bm{z}} \geq 0, \quad \bar{\bm{z}} \geq 0$$

ここで，$\bar{M}$は歪対称行列であり，$\bar{M}^\top = - \bar{M}$である.
ここで，この式は$(\bm{x},\bm{y},\beta)$が解となる時，正の定数$\lambda >0$を用いた，$\lambda(\bm{x},\bm{y},\beta)$も解になる．
つまり，$\beta > 0$である問題(2)の解$(\bm{x},\bm{y},\beta)$が得られた場合には，$(\bm{x}/\beta,\bm{y}/\beta,1)$は，問題(1)の解になっている．一方で，$\beta = 0$の解しか存在しなければ，問題(1)は解を持たない．よって，問題(1)の代わりに問題(2)を解くことができる．

さらに，問題(2)に対して，normalizing variableと呼ばれる新たな変数$\theta$を追加した問題を考える．
以下のように行列$M$と変数ベクトル$\bm{z}$を定義する．

$$M \coloneqq \begin{bmatrix} \bar{M} & \bm{r} \\ -\bm{r}^{\top} & 0 \end{bmatrix} = \begin{bmatrix} \begin{bmatrix} 0_{m \times m} & A & - \bm{b} \\ -A^{\top} & 0_{n \times n} & c \\ \bm{b}^{\top} & - \bm{c}^{\top} & 0 \end{bmatrix} & \bm{r} \\ -\bm{r}^{\top} & 0 \end{bmatrix}, \quad \bm{z} \coloneqq \begin{bmatrix} \bar{\bm{z}} \\ \theta \end{bmatrix} = \begin{bmatrix} \bm{y} \\ \bm{x} \\ \beta \\ \theta \end{bmatrix},\quad \bm{q} = \begin{bmatrix} 0_{\bar{n} - 1} \\ \bar{n} \end{bmatrix}$$

$M \in \R^{\bar{n} \times \bar{n}}$は再び歪対称行列である($M^\top = - M$)．
ここで，$\bar{n} = m + n + 2$とすると，$\bm{r}$は以下のように定義している．

$$\bm{r} = e_{\bar{n}-1} - \bar{M} e_{\bar{n}-1}$$

$e_\ell$は長さが$\ell$ですべての要素が1のベクトルである．

これらを用いて，以下の問題を考える．

$$\tag{3} M\bm{z} \geq - \bm{q}, \quad \bm{z} \geq 0,$$

明らかに，問題(3)の解が$\bm{z} = (\bar{\bm{z}}, 0)$であるとき，これは問題(2)の解となる．
よって，$\theta = 0$かつ$\beta > 0$の解$\bm{z}$は問題(2)の解となり，これは(P)と(D)の双対ギャップを0にするような解を与える．

さらに問題(3)を変形して，自己双対になるように目的関数を定義した問題

$$\tag{HLP} \min \{\bm{q}\cdot \bm{z} : M\bm{z} \geq -\bm{q}, \bm{z} \geq 0\}$$

を考える．この問題の双対問題もまたこの形になる．

問題(3)と$(HLP)$の解の関係を考える．
$(HLP)$の目的関数について$\bm{q}\geq 0,\ \bm{z} \geq 0$であるので，$\bm{q}\cdot \bm{z} \geq 0$である．また，$\bm{q} \geq 0$であるから，$\bm{z} = 0$も実行可能解なので，(HLP)の最適値は$0$である．また，$\bm{q} \cdot \bm{z}$はその定義から$\bm{q} \cdot \bm{z} = \bar{n}\theta$である．
よって，(3)が$\theta=0,\ \beta >0$の解を持つと仮定すると，$\bm{z} = (\bar{\bm{z}}, 0)$であり，$\bar{n}\theta = 0$となるので，これは$(HLP)$の最適解であることがわかる．
逆に，$(HLP)$が$\beta > 0$の最適解を持つとすると，最適解の目的関数値は0なので，$\theta = 0$になり，このことから(3)が$\theta=0,\ \beta >0$という解を持つことがわかる．

よって，式(3)が$\theta=0,\ \beta >0$の解を持つのは$(HLP)$が$\beta > 0$である最適解を持つときに限られる．

Self-dual embeddingの実装

この手法では，インフィージブル内点法の場合のように，行列の中に各ステップでの値などが含まれないので，問題の情報だけで$M$を構築できる．なので，この時点で，$M$を実装してしまおう．
ここで，A,b,cは全て(P)の形式で定義された線形計画問題の定数である．

import numpy as np

def generate_embedding_model(
    A: np.ndarray,
    b: np.ndarray,
    c: np.ndarray,
) -> np.ndarray:
    """Generate self-dual embedding model"""
    m, n = A.shape

    M_bar = np.block(
        [
            [np.zeros((m, m)), A, -b.reshape(-1, 1)],
            [-A.T, np.zeros((n, n)), c.reshape(-1, 1)],
            [b.reshape(1, -1), -c.reshape(1, -1), np.zeros((1, 1))],
        ]
    )

    e = np.ones(m + n + 1)
    r = e - M_bar @ e

    M = np.block([[M_bar, r.reshape(-1, 1)], [-r.reshape(1, -1), np.zeros((1, 1))]])

    return M

Central Path

slack変数$\bm{s}(\bm{z})$を以下のように定義する．

$$\bm{s}(\bm{z}) = M\bm{z} + \bm{q}$$

すると，$(HLP)$の実行可能領域は$\bm{z}\geq 0 ,\bm{s}\geq 0$で記述できる．

さて，ここで$(HLP)$のいくつかの性質について議論する．
内点法では，初期解として内点を用意する必要がある．興味深いことに$(HLP)$は，すべての要素が1のベクトル$e_{\bar{n}}$を自明な内点として持っている．これから，この点を示そう．
$e_{\bar{n}} > 0$は明らかなので，$\bm{s}(e_{\bar{n}}) > 0$を示そう．直接代入すると，

$$\bm{s}(e_{\bar{n}}) = Me_{\bar{n}} + \bm{q} = \begin{bmatrix} \bar{M} & \bm{r} \\ -\bm{r}^{\top} & 0 \end{bmatrix} \begin{bmatrix} e_{\bar{n} - 1} \\ 1 \end{bmatrix} + \begin{bmatrix} 0 \\ \bar{n} \end{bmatrix} = \begin{bmatrix} \bar{M} e_{\bar{n} - 1} + \bm{r} \\ -\bm{r}^{\top}\cdot e_{\bar{n} - 1} + \bar{n} \end{bmatrix} = \begin{bmatrix} e_{\bar{n} - 1} \\ 1 \end{bmatrix}> 0$$

さらに，(証明はしないが)正の定数$\mu > 0$に対して

$$\bm{z} \odot \bm{s}(\bm{z}) = \mu e, \bm{z} \geq 0, \bm{s}(\bm{z}) \geq 0$$

を満たすような非負の$\bm{z}$が高々一つ存在する．ここで$\odot$は要素積を表す．
$\mu > 0$に対して上の式を満たす解$z(\mu)$を(HLP)の$\mu$-centerと呼び，その集合を

$$\{z(\mu): \mu >0\}$$

と書く．これは$(HLP)$の実行可能空間上の曲線となり，これをcentral pathと呼ぶ．

このcentral path上の目的関数値を考えると

$$\bm{q}\cdot \bm{z}(\mu) =(\bm{s}(\bm{z}(\mu)) - M\bm{z}(\mu)) \cdot \bm{z}(\mu) = \bm{s}(\bm{z}(\mu)) \cdot \bm{z}(\mu) = \mu \bar{n}$$

となる．ここで$M$が歪対称行列であることを用いた．
central path上の目的関数は$\mu$について線形なので，$\mu$を0に近づけていくと単調に目的関数値は0に近づく．つまり，最適解に$\bm{z}(\mu)$は近づいていく．

更新式

Central Path上で徐々に$\mu$を小さくしていくことで，最適解に近づいていくことがわかった．
また，既に$\bm{z} = e_{\bar{n}}$の時には，$s(e_{\bar{n}}) = e_{\bar{n}}$となることはわかっており，これは$\mu=1$の時の$\mu$-centerである．よって，この$\mu=1$の点から初めて，central pathを辿るようにNewton法を行えば，最適解に辿り着くことができる．
なので，ここでは更新式について考えることにしよう．

$\bm{z}$が$\bm{s}(\bm{z}) > 0$を満たす$(HLP)$の解であるとする．
この時，$\bm{z}$からわずかに移動した点$\bm{z}^+$も$\mu$-centerに含まれるかを考える．
$\bm{z}$から$\Delta \bm{z}$動いた先に$\bm{z}^+= \bm{z}+ \Delta \bm{z}$があるとする．この時，移動した先での$\bm{s}(\bm{z}^+)$は

$$\bm{s}(\bm{z}^+) = M (\bm{z} + \Delta \bm{z} ) = \bm{s}(\bm{z}) + M\Delta \bm{z}$$

であるので，$\bm{s}$の変化量は

$$\Delta \bm{s} = \bm{s}(\bm{z}^+) - \bm{s}(\bm{z}) = M\Delta \bm{z}$$

である．これは$\Delta \bm{z}\cdot \Delta \bm{s} =\Delta \bm{z}^\top M\Delta \bm{z} = 0$なので，$\Delta \bm{z}$と$\Delta \bm{s}$が互いに直交することがわかる．

$\bm{z}^+$が$\mu$-centerとなるような$\Delta \bm{z}$方向が存在するかが知りたいので，

$$\bm{z}^+\odot \bm{s}(\bm{z}^+) = \bm{z}\odot \bm{s}(\bm{z}) + \Delta \bm{z} \odot \bm{s}(\bm{z}) + \bm{z}\odot\Delta \bm{s} + \Delta \bm{z} \odot\Delta \bm{s} = \mu e$$

と，先ほど得た$\Delta \bm{s} = M\Delta \bm{z}$という式を組み合わせて$\Delta \bm{s} ,\Delta \bm{z}$を求めることにする．2次の項$\Delta \bm{z} \odot\Delta \bm{s}$を落とすと，計算すべき連立方程式は

$$\left\{ \begin{aligned} &M\Delta \bm{z} - \Delta \bm{s} = 0\\ &\Delta \bm{z} \odot \bm{s}(\bm{z}) + \bm{z}\odot\Delta \bm{s} = \mu e - \bm{z}\odot \bm{s}(\bm{z}) \end{aligned} \right.$$

となる．上式から$\Delta \bm{s}=M\Delta \bm{z}$がわかっており，これを下の式に代入して，$\Delta \bm{z}$について解けば，

$$\begin{cases} \Delta \bm{z} &= (S + ZM)^{-1}(\mu e - \bm{z}\odot \bm{s})\\ \Delta \bm{s} &= M\Delta \bm{z} \end{cases}$$

となる．ここで，$Z = \mathrm{diag}(z), S = \mathrm{diag}(s)$である．

この差分は以下のように実装できる．

def compute_newton_direction(
    M : np.ndarray, z: np.ndarray, s: np.ndarray, mu: float
) -> tuple[np.ndarray, np.ndarray]:
    """Compute Newton direction for the central path"""

    S = np.diag(s)
    Z = np.diag(z)
    e = np.ones(len(z))

    dz = np.linalg.inv(S + Z @ M) @ (mu * e - Z @ s)
    ds = M @ dz
    return dz, ds

この$\Delta \bm{z}, \Delta \bm{s}$を用いて，$\bm{z},\bm{s}$を以下のように更新すれば良い．

$$\begin{cases} \bm{z}^{(k+1)} &= \bm{z}^{(k)} + \Delta\bm{z}\\ \bm{s}^{(k+1)} &= \bm{s}^{(k)} + \Delta\bm{s} \end{cases}$$

最適解を得るためには，$\mu$も小さくしていかなければいけない．しかし，$\mu$の変化が大きすぎてはcentral pathから大きく外れてはいけないので，ある程度の変化に止める必要がある．ここでは，$\mu \leftarrow (1-\gamma)\mu$と更新し，$\gamma = 1/\sqrt{2\bar{n}}$と設定する．
また，$(HLP)$は最適解で目的関数値が0になることがわかっているので，$\varepsilon > 0$を十分に小さい定数として，$\mu\bar{n} < \varepsilon$を収束条件として使うことができる．

また，$(HLP)$で得られる最適解は$\beta=1$とは限らないので，(P),(D)の最適解を得るためには，$(HLP)$の最適解での$\beta$の値でスケールしなおす必要がある点に注意する必要がある，

Self-Dual Embeddingを使った内点法の実装

以上をまとめると，Self-Dual Embeddingを使った内点法としては以下のように実装できる．

def solve_lp(
    A: np.ndarray,
    b: np.ndarray,
    c: np.ndarray,
    max_iter: int = 100,
    tol: float = 1e-8,
) -> tuple[np.ndarray, np.ndarray]:
    """Solve linear programming problem using self-dual embedding"""
    
    m, n = A.shape
    n_bar = n + m + 2
    # initialize
    z = np.ones(n_bar)
    M, q = generate_embedding_model(A, b, c)

    mu = 1.0
    gamma = 1.0 / np.sqrt(2.0 * n_bar)

    for _ in range(max_iter):
        if n_bar * mu < tol:
            break

        mu = (1.0 - gamma) * mu
        s = M @ z + q

        dz, ds = compute_newton_direction(M, z, s, mu)
        z = z + dz
        s = s + ds

    if n_bar * mu < tol:
        x = z[m : n + m]
        y = z[: m]
        beta = z[m + n]
        return (
            x / beta,
            y / beta,
        )

    raise ValueError("Failed to converge")

最後に

Jijでは数理最適化・量子コンピュータに関連するソフトウェア開発を行っています。
OSSも複数開発、提供しています。Discordのコミュニティへぜひ参加してください！
https://discord.gg/Km5dKF9JjG

また，Jijでは各ポジションを絶賛採用中です！
ご興味ございましたらカジュアル面談からでもぜひ！ご連絡お待ちしております。
数理最適化エンジニア（採用情報）
量子コンピューティングソフトウェアエンジニア（採用情報）
Rustエンジニア(アルゴリズムエンジニア)（採用情報）
研究チームPMO（採用情報）
オープンポジション等その他の採用情報はこちら

参考文献

[1] Y. Ye, M. J. Todd, and S. Mizuno, "An O(√nL)-iteration homogeneous and self-dual linear programming algorithm", Mathematics of Operations Research, 19,53–67, 1994.
[2] 水野 眞治, 「自己双対線形計画問題と内点法」,学習用テキスト 内点法
[3] C. Roos, T. Terlaky and J. P. Vial, "Interior Point Methods for Linear Optimization", Springer 2005