冷めたコーヒー

Weniger, aber besser

混合整数計画問題のソルバー scipy.optimize.milp と Python-MIP でのデフォルトのソルバー CBC の簡単な比較

前回の記事で,SciPy==1.9.0 で新たに導入された scipy.optimize.milp に関する簡単な紹介をしました.

mirucacule.hatenablog.com

本記事では,上記の記事で扱った問題と同様の問題を Pyhton-MIP で求解する方法を紹介します.Python-MIP のより詳細な使用方法については公式ドキュメント [1] を参照してください.

解くべき混合整数計画問題(再掲)

次の混合整数計画問題を考えます:

\begin{align} \min_{\mathbf{x}} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{A}\mathbf{x} \leqslant \mathbf{b},\\ &\mathbf{x} \in \mathbb{Z}^n,\,\mathbf{x} \geqslant \mathbf{0}. \end{align}

ここで,$\mathbf{c}\in\mathbb{R}^n, \, \mathbf{A}\in\mathbb{R}^{m\times n}, \, \mathbf{b}\in\mathbb{R}^{m}$ とします.

scipy.optimize.milp による求解(再掲)

必要なパッケージを import しておきます.

import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds

解くべき問題の係数行列などを定義します.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

制約条件に関わる部分を定義します.

l = np.zeros_like(c)
u = np.full_like(c, np.inf)
b_u = b
b_l = np.full_like(b_u, -np.inf)

bounds = Bounds(l, u)
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)

求解します.

%%time
res = milp(
    c=c,
    bounds=bounds,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下の通りです.

            fun: 7.441844504063867
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 7.44111388697805
        mip_gap: 9.8176881473139e-05
 mip_node_count: 49538
         status: 0
        success: True
              x: array([1., 0., 1., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0., 4., 0., 0., 0.,
       0., 0., 0., 2., 1., 0., 0., 0., 0., 3., 1., 1., 0., 1., 3., 0., 0.,
       1., 0., 0., 1., 0., 0.])
CPU times: user 17.8 s, sys: 665 ms, total: 18.4 s
Wall time: 17.6 s

Python-MIP の CBC による求解

必要なパッケージを import しておきます.

import numpy as np
import mip

同様の問題を扱うため,係数行列などの定義は同様です.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

次に,決定変数,目的関数,制約条件をそれぞれ定義して求解を行います.

model = mip.Model()
# define decision variable
x = {}
for i in range(n):
    x[i] = model.add_var(lb=0, ub=mip.INF, var_type=mip.INTEGER)
# define objective function
model.objective = mip.minimize(mip.xsum(c[i] * x[i] for i in range(n)))
# define constraints
for i in range(m):
    model.add_constr(
        mip.xsum(A[i,:][j] * x[j] for j in range(n)) <= b[i]
    )
# solve
model.optimize()

model.optimize() を実行すると以下のような出力が得られます.実際には,最適化計算の各反復における情報が出力されますが,ここでは省略します.

Cbc0032I Strong branching done 29094 times (142226 iterations), fathomed 305 nodes and fixed 796 variables
Cbc0041I Maximum depth 33, 3997 variables fixed on reduced cost (complete fathoming 1136 times, 157339 nodes taking 667811 iterations)
Total time (CPU seconds):       8.99   (Wallclock seconds):       9.63

<OptimizationStatus.OPTIMAL: 0>

得られた最適値と最適解を確認してみましょう.

print('Objective value: {model.objective_value:.15}'.format(**locals()))
print('Solution: ', end='')
for v in model.vars:
    if v.x > 1e-5:
        print('{v.name} = {v.x}'.format(**locals()))
        print('          ', end='')

実行すると以下のような結果が得られます.

Objective value: 7.44184450406387
Solution: var(0) = 1.0
          var(2) = 1.0
          var(5) = 1.0
          var(6) = 1.0
          var(10) = 1.0
          var(13) = 4.0
          var(20) = 2.0
          var(21) = 1.0
          var(26) = 3.0
          var(27) = 1.0
          var(28) = 1.0
          var(30) = 1.0
          var(31) = 3.0
          var(34) = 1.0
          var(37) = 1.0

最適値,最適解ともに scipy.optimize.milp で計算した結果と整合していることが分かります.なお,求解にかかった時間は CPU 時間だけで比較すると,milp は 17.8 秒,Python-MIP (cbc) は 8.99 秒でした.

参考

[1] Python MIP Documentation

SciPy==1.9.0 で使用可能になった混合整数計画問題のソルバー scipy.optimize.milp に関する備忘録

はじめに

この記事の目的:

  • SciPy==1.9.0 で導入された scipy.optimize.milp の簡単な説明
  • scipy.optimize.milp の使い方の解説

この記事で扱わないこと

  • scipy.optimize.milp で採用されている最適化手法
  • 他のソルバーとの比較

[追記]

  • 公式ドキュメントにナップザック問題を解くチュートリアルが掲載されたので参考 [3] にリンクを追記しました (2022/09/08)

SciPy==1.9.0 のリリース

Atsushi さんによる以下の記事を参照してください.

myenigma.hatenablog.com

先行記事

書籍 [2] で挙げられている問題を scipy.optimize.milp を用いて解く方法について解説されています.

github.com

scipy.optimize.milp で解ける最適化問題

scipy.optimize.milp では以下の形式の混合整数計画問題を解くことができる:

\begin{align} \min_{\mathbf{x} \in \mathbb{R}^n} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{b}_l \leqslant\mathbf{A}\mathbf{x} \leqslant \mathbf{b}_u,\\ &\mathbf{l} \leqslant \mathbf{x} \leqslant \mathbf{u},\\ &x_i \in \mathbb{Z},\, i \in \mathcal{I} \subseteq [n]. \end{align}

ここで,$\mathcal{I}$ は添字集合 $[n] := \{1,2,\dots,n\}$ の部分集合 (e.g., $\mathcal{I} = \{2, 3, 5 \}$) を表す.

残念ながら,scipy.optimize.milp を用いて上記の形式でない混合整数計画問題を解くことは今のところできません.頑張って上記の形式に書き直すか,他の最適化モデリング言語などを用いましょう.

scipy.optimize.milp の使い方

導入

以下で scipy==1.9.0 のインストールをしましょう:

python -m pip install scipy==1.9.0

以下では,scipy.optimize に含まれるパッケージを使用します.ついでに numpy も使用するので import しておきましょう.

import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds

例 1

まず,簡単な問題として Mathworksintlinprog のドキュメントで扱われている以下のような混合整数計画問題を考えます:

\begin{align} \min_{x_1,x_2} \quad &8x_1 + x_2 \\ \text{subject to} \quad&x_1 + 2x_2 \geqslant -14,\\ &-4x_1 - x_2 \leqslant -33,\\ &2x_1 + x_2 \leqslant 20,\\ &x_2 \in \mathbb{Z}. \end{align}

この問題はベクトルと行列を用いて次のように表すことができます:

\begin{align} \min_{x_1,x_2} \quad &\begin{bmatrix}8 \\ 1\end{bmatrix}^\top \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\ \text{subject to} \quad&\begin{bmatrix}-1 & -2 \\ -4 & -1 \\2 & 1\end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix} \leqslant \begin{bmatrix}14 \\ -33 \\ 20\end{bmatrix},\\ &x_2 \in \mathbb{Z}. \end{align}

この問題を scipy.minimize.milp を用いて解けるように係数ベクトルや行列を以下のように定めます.なお,この問題は $(x_1,x_2)^\top$ に関する上下限制約はないので特別指定する必要はありません.また,$\mathbf{Ax}$ に関する下限制約もないので,以下のように b_l の各要素に -np.inf を設定します.整数制約は $x_2$ に対してのみ課されるため,以下のように integrality = np.array([0, 1]) と設定します.なお,0 は連続変数,1 は整数制約を表します.他に 23 を設定することもできますが,それらの説明は 公式ドキュメント を参照してください.また,integrality を指定しなかった場合,デフォルトでは 0 が指定されます.

c = np.array([8, 1])
A = np.array([[-1, -2], [-4, -1], [2, 1]])
b_u = np.array([14, -33, 20])
b_l = np.full_like(b_u, -np.inf)
constraints = LinearConstraint(A, b_l, b_u) $b_l <= Ax <= b_u
integrality = np.array([0, 1]) # x_2 is integer

求解は以下のように行います.

%%time
res = milp(
    c=c,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下のようになります.

            fun: 59.0
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 59.0
        mip_gap: 0.0
 mip_node_count: 1
         status: 0
        success: True
              x: array([6.5, 7. ])
CPU times: user 5.49 ms, sys: 0 ns, total: 5.49 ms
Wall time: 3.24 ms

Mathworksドキュメント に掲載されている結果と同じ結果が得られていることを確認できます.

例 2

次の混合整数計画問題を考えます:

\begin{align} \min_{\mathbf{x}} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{A}\mathbf{x} \leqslant \mathbf{b},\\ &\mathbf{x} \in \mathbb{Z}^n,\,\mathbf{x} \geqslant \mathbf{0}. \end{align}

ここで,$\mathbf{c}\in\mathbb{R}^n, \, \mathbf{A}\in\mathbb{R}^{m\times n}, \, \mathbf{b}\in\mathbb{R}^{m}$ とします.この問題の実行可能領域 $\mathcal{F}=\left\{ \mathbf{x} \in \mathbb{Z}^n \mid \mathbf{A}\mathbf{x} \leqslant \mathbf{b},\, \mathbf{x} \geqslant \mathbf{0} \right\}$ に実行可能解が存在することを保証するために,$\overline{\mathbf{x}} \in \mathbf{Z}^n$ をランダムに生成し,その $\overline{\mathbf{x}} \in \mathbf{Z}^n$ が実行可能領域 $\mathcal{F}$ に含まれるように $\mathbf{A}$ と $\mathbf{b}$ を設定します.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

決定変数 $\mathbf{x}$ に対する下限制約が存在するため,scipy.optimize.Bounds を使用します.決定変数 $\mathbf{x}$ に対する上限制約は存在しないので,以下のように u の各要素に np.inf を設定します.また,今回の問題では決定変数のすべての要素は整数であることを課しているため,integrality のすべての要素を 1 とします.

l = np.zeros_like(c)
u = np.full_like(c, np.inf)
b_u = b
b_l = np.full_like(b_u, -np.inf)

bounds = Bounds(l, u)
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)

求解します.

%%time
res = milp(
    c=c,
    bounds=bounds,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下の通りです.

            fun: 7.441844504063867
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 7.44111388697805
        mip_gap: 9.8176881473139e-05
 mip_node_count: 49538
         status: 0
        success: True
              x: array([1., 0., 1., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0., 4., 0., 0., 0.,
       0., 0., 0., 2., 1., 0., 0., 0., 0., 3., 1., 1., 0., 1., 3., 0., 0.,
       1., 0., 0., 1., 0., 0.])
CPU times: user 17.8 s, sys: 665 ms, total: 18.4 s
Wall time: 17.6 s

決定変数の次元 $n=40$ の問題を解いてみました.およそ 20 秒弱の時間で最適解を計算することができました.

おわりに

この記事では SciPy==1.9.0 で新たに導入された scipy.optimize.milp を用いて混合整数計画問題を解く方法について解説をしました.

参考

[1] scipy.optimize.milp
[2] 『Pythonではじめる数理最適化: ケーススタディモデリングのスキルを身につけよう』

[3] Knapsack problem example

リーマン多様体上の最適化―特異値分解の例を通して―

はじめに

以前(2019 年 11 月)に「リーマン多様体上の最適化の初歩と Pymanopt による数値実験」 という記事を投稿した.

mirucacule.hatenablog.com

記事内で用いた最適化 Toolbox である Pymanopt のバージョンアップに伴い,実行方法に変更が生じたため改めて書き直そうと思ったのが本記事を執筆するに至った経緯である.また,リーマン多様体上での共役勾配法を簡単に勉強したので,それを用いて行列の特異値分解を行う問題について扱いたい.

なお,本記事を執筆するに際して参考にした文献は 参考文献 にまとめてある.特に,本記事の内容は [1] の Chapter 4 に書かれている内容の重要と思われる箇所を抜粋してまとめたものであり,より深く知りたい方は [1] を参照されたい.本記事で扱っていない Newton 法や,共役勾配法と Newton 法を組み合わせたアルゴリズムに関する記述もあり楽しめると思う.

また,Pymanopt の使い方についてのみ知りたい読者は Pymanopt による求解 の節を確認されたい.より詳細な使い方や例については以下を参照されたい.

特異値分解

一般に,行列 $A \in \mathbb{R}^{m \times n},\ (m \ge n)$ は直交行列 $U \in \mathbb{R}^{m \times p},\, V \in \mathbb{R}^{n \times p}$ を用いて $$ A = U D V^\top \tag{1} \label{eq:1} $$

のように分解することができる.ここで,$D \in \mathbb{R}^{n\times n}$ は対角行列であり,その対角成分を特異値と呼ぶ.また,分解 (\ref{eq:1}) を行列 $A$ の特異値分解と呼ぶ.具体的な例としては [6] などを参照されたい.

特異値分解最適化問題

特異値分解は以下の最適化問題と密接な関係がある.

$$ \begin{align} (\mathrm{P}_1)\qquad \min_{U,V}\quad& -\mathrm{tr}(U^\top AVN)\\ \text{subject to}\quad& U\in \mathbb{R}^{m\times p},\, V \in \mathbb{R}^{n\times p},\, U^\top U = V^\top V = I_p. \end{align} $$

ここで,$A \in \mathbb{R}^{m\times n}$ は所与の行列,$N \in \mathbb{R}^{p\times p}$ は所与の対角行列をそれぞれ表す. 特に,$N$ の対角成分は $\mu_1 > \mu_2 > \dots > \mu_p > 0$ を満たす.また,$1 \le p \le n$ である.このとき,以下の命題が成り立つ.

命題 [4, Proposition 2.1] $m \ge n$ とする.最適化問題 ($\mathrm{P}_1$) の最適解を $(U_\star, V_\star)$ とする.このとき,$U_\star$ と $V_\star$ の列ベクトルたちは,それぞれ $A$ の小さい方から $p$ 個の特異値に属する左特異ベクトル,右特異ベクトルとなる.

この命題より,$U_\star = (\mathbf{u}_1,\dots,\mathbf{u}_p),\, V_\star = (\mathbf{v}_1,\dots,\mathbf{v}_p)$ としたとき,行列 $A$ の特異値分解

$$ A = \sum_{i=1}^{p} d_i \mathbf{u}_i \mathbf{v}_i^\top $$

と表される.

リーマン多様体上での特異値分解

ユークリッド空間上の最適化問題 ($\mathrm{P}_1$) は次のリーマン多様体上の最適化問題と等価である.

$$ \begin{align} (\mathrm{P})\qquad \min_{U,V}\quad& F(U, V) = -\mathrm{tr}(U^\top AVN)\\ \text{subject to}\quad& (U,V) \in \mathrm{St}(p,m) \times \mathrm{St}(p,n). \end{align} $$

ここで,制約条件に現れる $\mathrm{St}(p,n)$ は Stiefel 多様体といい,

$$\mathrm{St}(p,n) := \left\{Y \in \mathbb{R}^{n \times p} \mid Y^\top Y = I_p \right\}$$

で定義される.また,$\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ は $\mathrm{St}(p,m)$ と $\mathrm{St}(p,n)$ の積多様体を表す.

$\mathrm{St}(p,m)\times\mathrm{St}(p,n)$ 上の最適化

リーマン多様体 $\mathcal{M}$ 上での最適化を考える上で必要となるレトラクション,勾配について述べる.一般的な定義や導出方法については [1] や [2] を参照されたい.ここでは単に Stiefel 多様体の積多様体におけるレトラクションと勾配について,その結果についてのみ提示するに留める.以下の図は,リーマン多様体上の最適化手法で生成される点列のイメージである.リーマン多様体 $\mathcal{M}$ 上の点 $\mathbf{x}_k$ を探索方向 $\mathbf{\xi}_k \in T_{\mathbf{x}_k} \mathcal{M}$ 方向に $t_k>0$ でスケールして動かし,レトラクション $R_{\mathbf{x}_k}$ という操作によってリーマン多様体 $\mathcal{M}$ 上に写すことを繰り返し行いながら探索を行う.

リーマン多様体上のレトラクション

接空間

多様体 $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ 上の点 $(U, V)$ における接空間 $T_{(U,V)}(\mathrm{St}(p,m) \times \mathrm{St}(p,n))$ は次で与えられる.

\begin{align} T_{(U,V)}(\mathrm{St}(p,m) \times \mathrm{St}(p,n)) &\simeq T_{U}(\mathrm{St}(p,m)) \times T_{V}(\mathrm{St}(p,n)) \\ &= \left\{ (\xi, \zeta) \in \mathbb{R}^{m \times p} \times \mathbb{R}^{n \times p} \mid \xi^\top U + U^\top\xi = \zeta^\top V + V^\top\zeta = O_p \right\}. \end{align}

ここで,$O_p$ は $p$ 次零行列を表す.

レトラクション $R_{(U,V)}$

レトラクションは接空間 $T_{(U, V)} (\mathrm{St}(p,m) \times \mathrm{St}(p,n))$ からリーマン多様体 $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ への写像である.図の緑色の破線がレトラクションに対応する.点 $(U,V)\in\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ における $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ のレトラクションは次のように定められる.

$$ R_{(U, V)} (\xi, \zeta) = \left( \mathrm{qf}( U + \xi ), \mathrm{qf}( V + \zeta ) \right), \quad (\xi, \zeta) \in T_{(U, V)} (\mathrm{St}(p,m) \times \mathrm{St}(p,n)). $$

ここで,$\mathrm{qf}:\mathbb{R}^{m \times n} \to \mathrm{St}(p, n)$ であり,具体的にはフルランク行列 $X \in \mathbb{R}^{n \times p}$ に対して QR 分解を行ったときの $Q \in \mathrm{St}(p, n)$ を返す写像として定義される.ここで,$R$ は $p$ 次上三角行列であり,その対角成分はすべて正とする.

勾配 $\mathrm{grad} F(U,V)$

目的関数 $F(U,V) = -\mathrm{tr}(U^\top AVN)$ の点 $(U, V) \in \mathrm{St}(p,m)\times\mathrm{St}(p,n)$ における勾配は以下で定義される.

$$ \mathrm{grad} F(U, V) = \left( U \mathrm{sym}(U^\top AVN) - AVN, V \mathrm{sym}(V^\top A^\top UN) - A^\top UN \right). $$

ここで,$\mathrm{sym}(X)$ は正方行列 $X \in \mathbb{R}^{k\times k}$ に対して $\mathrm{sym}(X) = \frac{X + X^\top}{2}$ で定義される.

リーマン多様体上での共役勾配法

ユークリッド空間上での最適化アルゴリズムとして共役勾配法(Conjugate gradient descent method; CG 法)がある.この手法は,暫定解における探索方向を上手く選ぶことによって最急降下法のようなジグザグとした挙動を抑えようというアルゴリズムである.もともとは線形方程式系 $Ax = b$ を解く手法として考案されたアルゴリズムだったが,関数を最小化する問題に対しても適用できるということで様々な拡張が考えられてきた.この辺の詳しい話は,例えば [3] の Chapter 5 に詳しく書かれているので興味がある方は参照されたい. また,リーマン多様体上での共役勾配法アルゴリズムに関する和文記事は例えば [5] を参照されたい.

Pymanopt による求解

以下では,最適化問題 (P) を Pymanopt で解く手順について説明する.Pymanopt はリーマン多様体上での最適化を行うフレームワークである.本記事における各モジュールのバージョンは以下の通りである.

  • pymanopt == 2.0.0
  • autograd == 1.4
  • numpy == 1.23.0

モジュールのインポート

まずは必要なモジュールをインポートする.今回,リーマン多様体として Stiefel 多様体と Stiefel 多様体の積多様体を考えるため,pymanopt.manifolds から StiefelProduct をインポートする.また,今回は最適化手法として共役勾配法を用いるため,pymanopt.optimizers から ConjugateGradient をインポートする.

from pymanopt.manifolds import Stiefel, Product
from pymanopt.optimizers import ConjugateGradient
from pymanopt import Problem
import pymanopt
import autograd.numpy as np

なお,リーマン多様体に関するその他の選択肢は以下の通りである.

また,最適化手法に関するその他の選択肢は以下の通りである.

解くべき最適化問題の定義

問題を定義する方法について述べる.今回は行列 $A$ の左特異ベクトル,右特異ベクトルが既知となるように直交行列 $U_r \in \mathbb{R}^{m\times n}, \, V_r \in \mathbb{R}^{n\times n}$ をランダムに生成し,特異値を対角成分に持つ対角行列 $S$ は $S = \mathrm{diag} [n, n-1, \dots, 1]$ となるように定義する.これによって,アルゴリズムによって得られた解を最適解と比較することができる.

np.random.seed(2022)

# setting problem
m, n, p = 5, 4, 3
N = np.diag(np.arange(p, 0, -1))
Ur, _ = np.linalg.qr(np.random.normal(size=(m, n)))
Vr, _ = np.linalg.qr(np.random.normal(size=(n, n)))
S = np.diag(np.arange(n, 0, -1))
A = Ur @ S @ Vr.T

次に,解くべき最適化問題を Pymanopt の記法に従って記述する.まず,多様体は次のように定める.今回は,Stiefel 多様体の積多様体として記述されるリーマン多様体上での最適化問題を扱うため,次のように書くことができる.

# setting manifold
stiefel_1 = Stiefel(m, p)
stiefel_2 = Stiefel(n, p)
manifold = Product([stiefel_1, stiefel_2])

目的関数は次のように記述することができる.以前のバージョンと異なり,@pymanopt.function.autograd(manifold) といったデコレータが必要である.このデコレータを付けないと ValueError: Function 'cost' must be decorated with a backend decorator というエラーが返ってくる.

# define cost function
@pymanopt.function.autograd(manifold)
def cost(U, V):
    return -np.trace(U.T @ A @ V @ N)

以上の下で解くべき最適化問題Problem クラスのインスタンスとして次のように定める.引数にリーマン多様体と目的関数を与えればよい.

# define problem
problem = Problem(manifold=manifold, cost=cost)

最適化手法の定義

最適化問題を解くための最適化手法を定める.最適化アルゴリズム共役勾配法を扱うため ConjugateGradient クラスのインスタンスとして生成する.引数は以下の通りで,verbosity は出力されるログレベルを表し,0 が出力なしを表し,2 が最大限出力することを課す.詳しくは コード を確認されたい.

# define solver
solver = ConjugateGradient(
    max_time=100,
    max_iterations=100,
    min_gradient_norm=1e-6,
    min_step_size=1e-10,
    max_cost_evaluations=5000,
    verbosity=2,
    log_verbosity=0,
    beta_rule="PolakRibiere",
)

初期値は以下のように与える.初期値を自前で定義して与えることもできるが,pymanopt.manifolds から import したクラス(e.g., Stiefel, Product)は random_point というメソッドを持っており,これを使うと対象とするリーマン多様体上の変数をランダムに生成できる.

# initial guess
U0, _ = np.linalg.qr(np.random.rand(m, p))
V0, _ = np.linalg.qr(np.random.rand(n, p))
x0 = np.array([U0, V0])
# x0 = manifold.random_point()

最後に実行する.引数には,Problem クラスのインスタンスと初期値を与えればよい.

# run
res = solver.run(problem, initial_point=x0)
print(res)

出力内容

出力内容は以下の通り.反復回数 60 回で終了条件を満たし停止していることが確認できる.

Optimizing...
Iteration    Cost                       Gradient norm
---------    -----------------------    --------------
  1          +4.7770270621793562e+00    9.86333970e+00
  2          -4.1494305543650158e+00    9.87110012e+00
  3          -1.1198367922895978e+01    6.04524074e+00
  4          -1.3925445689443663e+01    3.76110809e+00
  5          -1.5989745840394253e+01    2.93949033e+00
(中略)
 56          -1.9999999999970512e+01    8.87200818e-06
 57          -1.9999999999987086e+01    7.52846533e-06
 58          -1.9999999999997979e+01    3.67079919e-06
 59          -1.9999999999999311e+01    1.58361854e-06
 60          -1.9999999999999545e+01    9.85344173e-07
Terminated - min grad norm reached after 60 iterations, 0.09 seconds.

OptimizerResult(point=[array([[-3.27751202e-04, -2.00772692e-01,  6.05821509e-02], 
       [ 1.75156668e-01,  7.33257570e-01, -5.71602349e-04],
       [ 2.31897948e-01,  5.11838715e-01,  4.82436752e-01],
       [ 7.09519651e-01, -3.96506778e-01,  4.94409705e-01],
       [ 6.41969813e-01,  5.31711091e-02, -7.20515996e-01]]), array([[-0.13814535, 
-0.50517496,  0.67497762],
       [-0.88146048, -0.32414775, -0.21981553],
       [ 0.14490487, -0.1186366 ,  0.56149162],
       [-0.4277217 ,  0.79098157,  0.42522171]])], cost=-19.999999999999545, iterations=60, stopping_criterion='Terminated - min grad norm reached after 60 iterations, 0.09 seconds.', time=0.09300017356872559, cost_evaluations=60, step_size=2.3126464564821429e-07, gradient_norm=9.853441734587333e-07, log={'optimizer': 'ConjugateGradient', 'stopping_criteria': {'max_time': 100, 'max_iterations': 100, 'min_gradient_norm': 1e-06, 'min_step_size': 1e-10, 'max_cost_evaluations': 5000}, 'optimizer_parameters': {'beta_rule': 'PolakRibiere', 'orth_value': inf, 'line_searcher': <pymanopt.optimizers.line_search.AdaptiveLineSearcher object at 0x0000022E9CA75C70>}, 'iterations': None})

今回は,行列 $A$ の構成の仕方から Pymanopt で得られた(近似)最適解の質を評価することができる.

def F_opt(S, N):
    p = N.shape[0]
    S = S[:p, :p]
    return -np.trace(S @ N)

U_star, V_star = res.point[0], res.point[1]

print('Optimal value obtained by conjugate descent: {}'.format(cost(U_star, V_star)))
print('Optimal value: {}'.format(F_opt(S, N)))

結果は以下の通り.

Optimal value obtained by conjugate descent: -19.999999999999545
Optimal value: -20

おわりに

本記事では,特異値分解を Stiefel 多様体の積集合上で最適化する問題について扱った.また,最適化問題 $(\mathrm{P})$ をリーマン多様体上の最適化を行うパッケージである Pymanopt を用いて求解する方法について述べた.

参考文献

[1] Sato, H. (2013). Riemannian Optimization Algorithms and Their Applications to Numerical Linear Algebra.
[2] Absil, P. A., Mahony, R., & Sepulchre, R. (2009). Optimization algorithms on matrix manifolds. In Optimization Algorithms on Matrix Manifolds. Princeton University Press.
[3] Nocedal, J., & Wright, S. J. (Eds.). (1999). Numerical optimization. New York, NY: Springer New York.
[4] Sato, H., & Iwai, T. (2013). A Riemannian optimization approach to the matrix singular value decomposition. SIAM Journal on Optimization, 23(1), 188-212.
[5] 【Python】リーマン多様体上の共役勾配法
[6] https://manabitimes.jp/math/1280
[7] https://pymanopt.org/docs/stable/

おまけ

最適化問題 $(\mathrm{P})$ を最急降下法共役勾配法で解いたときの勾配ノルムの推移を描いた結果を載せておく.

勾配ノルムの推移

双対問題の作り方―LP の例―

双対問題の作り方

以下の最小化問題を考える:

\begin{align} (\text{P})\qquad\min_{x} &\quad f(x) \\ \text{subject to} &\quad g_i(x) \le 0 \quad (i=1,2,\dots,m). \end{align}

関数 $f,\,g_i$ に対する連続性や微分可能性は適当に性質の良いものを考えることにしておく.また,以下で min や max の操作を多用するが,それぞれの存在性を仮定しておく.

ラグランジュ関数 $L:\mathbb{R}^n\times\mathbb{R}^m\to\mathbb{R}$ を以下で定義する:

$$ L(x,\lambda) := f(x) + \sum_{i=1}^{m}\lambda_i g_i(x), \quad (\lambda \geq 0). $$

このとき,適当に $\lambda \ge 0$ を固定して元の問題 P に対する制約条件を "緩和" した次の最適化問題を考える:

\begin{align} (\text{P}(\lambda))\qquad\min_{x} &\quad L(x,\lambda) := f(x) + \sum_{i=1}^{m} \lambda_i g_i(x)\\ \text{subject to} &\quad x \in \mathbb{R}^n. \end{align}

この問題を問題 P に対する "ラグランジュ緩和問題" といい,P($\lambda$) で表すこととする.また,そのときの最適値関数を $L(\lambda)$ で表すこととする.すなわち,

$$ L(\lambda) := \min_x L(x,\lambda) $$

である.

問題 P に対する最適解を $x^\star$ とすると,$x^\star$ は問題 P の制約条件を満たすので

$$ f(x^\star) \ge f(x^\star) + \sum_{i=1}^{m} \lambda_i g_i(x^\star) = L(x^\star,\lambda) $$

を満たす.すなわち,最適化問題 P の最適値 $f(x^\star)$ は,ラグランジュ緩和問題 P($\lambda$) を解くことによってその下界値(ここでいうところの $L(x^\star,\lambda)$)を得ることができる.このとき,緩和に用いたラグランジュ乗数 $\lambda \ge 0$ によって得られる下界値が変わってくる.

そこで,最適なラグランジュ乗数 $\lambda^\star$ が何なのか気になってくる.どうすればよいかというと,ラグランジュ緩和問題の最適値関数 $L(\lambda)$ を $\lambda$ に関して最大化した問題を考えればよい.すなわち,以下の最大化問題:

\begin{align} (\text{D}(\lambda))\qquad\max_{\lambda} &\quad L(\lambda)\\ \text{subject to} &\quad \lambda \geq 0 \end{align}

を解いて得られる $\lambda$ が気になるわけであり,この D($\lambda$) が最適化問題 P に対する双対問題というわけである.

おまけ: LP の場合

標準的な LP に対するラグランジュ双対問題がどのように導出されるのか見てみよう.以下の LP を考える:

\begin{align} (\text{P})\qquad\min_{x} &\quad c^\top x \\ \text{subject to} &\quad Ax \ge b,\,x\ge 0. \end{align}

ラグランジュ関数 $L$ を以下で定義する:

$$ L(x,\lambda) := c^\top x + \lambda^\top(b-Ax), \quad (\lambda \geq 0). $$

ラグランジュ緩和問題は次のように書ける:

\begin{align} (\text{P}(\lambda))\qquad\min_{x} &\quad L(x,\lambda) := c^\top x + \lambda^\top(b-Ax) \\ \text{subject to} &\quad x \ge 0. \end{align}

ここで,$x$ に関する非負制約もラグランジュ関数に取り込んでラグランジュ緩和問題を $x$ に関する無制約問題としてもいいのだが,$x\ge 0$ を残しても同様の議論ができるのでこのまま進めることにする.さて,ラグランジュ緩和問題 P($\lambda$) の目的関数は $x$ の項と $\lambda$ の項について整理してあげると

$$ L(x,\lambda) = (c-A^\top \lambda)^\top x + b^\top\lambda $$

となり,$b^\top\lambda$ の部分は $x$ に依存しない.一方,$c-A^\top\lambda$ について,$(A^\top\lambda)_i < 0$ なる $i\in[m]$が存在すると,対応する $x_i$ の値を大きくすることで幾らでも目的関数値を改善することができてしまうので,$c-A^\top\lambda \ge 0$ であることが必要である.また,ラグランジュ緩和問題に対する最適値関数 $L(\lambda)$ は $L(\lambda):=\min\{L(x,\lambda)\mid x\ge 0\}$ で定義されるが,$\lambda$ に関する部分は $b^\top\lambda$ のみであることに注意すると,ラグランジュ双対問題は以下で定式化される:

\begin{align} (\text{D}(\lambda))\qquad\max_{\lambda} &\quad b^\top\lambda\\ \text{subject to} &\quad A^\top\lambda\le c,\,\lambda \geq 0. \end{align}

Overleaf 環境を Docker でローカルに構築する手順

github.com

  1. overleaf/toolkit を clone する
$ cd ~
$ git clone https://github.com/overleaf/toolkit.git ./overleaf
  1. 設定ファイルの作成用コマンドを実行する
$ cd overleaf
$ bin/init

以下が表示されれば OK:

Copying config files to 'config/'
  1. Overleaf を起動するコマンドを実行する
$ bin/up

※ 内部で docker-compose を呼んでいます.初回は 10 分程度の時間がかかります.

  1. localhost:80/launchpad を開き,アカウントを作成する.
HAPPY OVERLEAF LIFE !