以前に Python で使える数理最適化モデリングツール CVXPY に関する記事を書きました:
そこでは，CVXPY の導入方法と簡単な例を解説しました．導入方法に関しては内容が古くなっている可能性が高いので，公式ドキュメントを参照されることをお勧めします．自分も最近 Ubuntu にインストールしましたが，基本的に公式ドキュメントに沿って進めていけば O.K. です．ただし，Ubuntu 環境下で Jupyter Notebook を用いて import cvxpy をすると，パスの読み込みが上手くいかない（異なるパスを参照しようとする）ことがあるようなので注意してください．
本記事の目的
本記事では，等式制約を含んだ最適化問題を CVXPY を用いて解くことを考えます．解くべき最適化問題は以下のように記述される問題とします:
\begin{align} \min_{x\in\mathbb{R}^n} \quad &f(x) \\ \mathrm{s.t.} \quad &h(x) = 0. \end{align}
ここで，$f\colon\mathbb{R}^n\to\mathbb{R}$，$h\colon\mathbb{R}^n\to\mathbb{R}$ であり，それぞれ連続であると仮定する（収束解析をするためには凸性などの仮定を設ける必要がありそうですが，今回は数値実験をするだけなのでその辺は大らかにやっていきます）．
解く具体的な問題
等式 $Ax = b$ で記述される制約条件の下で，$\|x\|_0$ を最小にする問題を扱います．この $\|x\|_0$ は $\ell_0$ ノルムと呼ばれ，ベクトル $x\in\mathbb{R}^n$ の非零成分の個数を表します．解くべき最適化問題は次のように記述されます:
\begin{align} \min_{x\in\mathbb{R}^n} \quad &\|x\|_0 \\ \mathrm{s.t.} \quad &Ax - b = 0. \end{align}
ここで，$A\in\mathbb{R}^{m\times n}$，$b\in\mathbb{R}^{m}$ とします．この問題は線形方程式 $Ax=b$ の解で非零の成分がなるべく多いような解を求める問題とも見ることができます．
CVXPY を用いた実装
素直に CVXPY を用いて解く場合は次のようにすれば O.K. です．
x = cp.Variable(shape=n)
f = cp.norm(x, 1)
# Solve with CVXPY.
cp.Problem(cp.Minimize(f), [A@x == b]).solve()
print("Optimal value from CVXPY: {}".format(f.value))
実行結果は次のようになります．本当に興味があるのは最適値ではなく最適解ではあるのですが，人工データなので良しとしましょう．
Optimal value from CVXPY: 2.811766447134533
ペナルティ関数法による実装
これで終わりだと詰まらないので，等式制約付きの最適化問題に対する汎用性の高い手法である「ペナルティ関数法」と「拡張ラグランジュ関数法」で同じ問題を解いてみましょう．なお，それぞれの手法の詳細な説明は，opt さん の以下の記事を参照してください（丸投げ；すごく分かりやすいのでみんな読みましょう）:
ペナルティ関数法では，元の制約付き問題の代わりに以下で記述される制約なしの最適化問題を解く方法です:
\begin{align} \min_{x\in\mathbb{R}^n} \quad &f(x) + \frac{\mu}{2} \| h(x) \|_2^2. \end{align}
ここで，$\mu > 0$ はパラメータであり，元の問題の最適解を得るためにはその値を適切に定める必要があります．ペナルティ関数法の実装は次のようになります．本来であれば，最急降下法などの手法を用いて実装する必要があるのですが，CVXPY を使えばこの程度の実装でお手軽に最適化することができます．
mu = 10.0
resid = A@x - b
# penalty function method
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
print("Optimal value from method of penalty function method: {}".format(f.value))
実行結果は次のようになります．
Optimal value from method of penalty function method: 2.693372134825692
この結果は元の問題に対して CVXPY を適用した最適値とかなり異なることが分かります．これは，ペナルティ関数の定数 $\mu$ の値が小さいことに起因すると考えられるので，より大きな値を用いて実行してみようと思います．
mu = 200.0 # changed
resid = A@x - b
# penalty function method
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
print("Optimal value from method of penalty function method: {}".format(f.value))
このときの実行結果は次のようになります．
Optimal value from method of penalty function method: 2.8040009637834995
解が改善されたことが分かります．最後に，様々な定数 $\mu > 0$ を用いて実験した結果を提示しておきます．ここでは，$\mu = 50, 100, \dots, 1000$ の計 20 通りで試しています．
f_opts = []
mus = [50.0 * i for i in range(1, 21)]
for mu in mus:
x = cp.Variable(shape=n)
f = cp.norm(x, 1)
resid = A@x - b
p = f + (mu/2) * cp.sum_squares(resid)
cp.Problem(cp.Minimize(p)).solve()
f_opts.append(f.value)
print("Optimal value from method of penalty function method: {}".format(f.value))
結果は省略し，得られたそれぞれのパラメータに対する解における目的関数値を図示しておきます．赤い点線は CVXPY で解くべき問題をそのまま求めたときの最適値を，青い点は各パラメータ $\mu > 0$ でペナルティ化した関数を最適化したときの最適値をそれぞれ表しています．
この結果を見ると，パラメータ $\mu > 0$ の値を $1000$ まで上げても CVXPY で求めた最適値から僅かながら離れていることが確認されます．なお，可視化のためのコードは次のように書くことができます．
x_range = np.array([50 * i for i in range(1, 21)])
plt.scatter(x_range, f_opts)
plt.hlines([2.811766447134533], 50, 1000, "red", linestyles='dashed')
plt.xlabel(r"$\mu$")
plt.ylabel("the optimal value")
plt.grid(linestyle='dotted', linewidth=1)
# plt.savefig("1.jpg") # 保存する場合はコメントを外す
拡張ラグランジュ関数法による実装
次に拡張ラグランジュ関数法の実装を行います．拡張ラグランジュ関数法では，上で述べたペナルティ関数法によって得られた最適解 $\hat{x}$ の情報を用いて，次の最適化問題を解くことを考えます:
\begin{align} \min_{x\in\mathbb{R}^n} \quad f(x) + \left< \hat{\lambda}, h(x) \right> + \frac{\mu}{2} \| h(x) \|_2^2. \end{align}
ここで，$\hat{\lambda} := \mu h(\hat{x})$ で定義されます．ペナルティ関数法と異なり，拡張ラグランジュ関数法では，問題を一度解いて得られる解を採用するのではなく，以下のような反復を行うことを考えます:
- (step 0): 適当な $\mu > 0$ を選び，$\hat{\lambda} \leftarrow 0$ とする．
- (step 1): 以下を交互に繰り返す:
- 問題を適当な方法（e.g., 最急降下法）で解き，得られた解を $\hat{x}$ とする．
- $\hat{\lambda} \leftarrow \hat{\lambda} + \mu h(\hat{x})$ とする．
上記の拡張ラグランジュ関数法を実装したコードは次のようになります．今回はパラメータ $\mu > 0$ の値を $50.0$ としました．
mu = 50.0
lamhat = np.zeros(m)
MAX_ITERS = 10
f_opts = []
x = cp.Variable(shape=n)
f = cp.norm(x, 1)
resid = A@x - b
lagr = f + lamhat.T@resid + (mu/2) * cp.sum_squares(resid)
# repeat
for t in range(MAX_ITERS):
cp.Problem(cp.Minimize(lagr)).solve()
lamhat += mu * resid.value
f_opts.append(f.value)
各反復で得られる暫定解 $\hat{x}$ での目的関数値 $f(\hat{x})$ をプロットした図を載せておきます．
ペナルティ関数法と比べて，しっかり赤い線（CVXPY で求めた最適値）にほぼ一致していることが確認されます．拡張ラグランジュ関数法では，ペナルティ関数法で適切に定める必要があった定数 $\mu > 0$ の値を今回のように大きくとらなくても解が求まるところが嬉しいポイントですね．
なお，可視化用のコードは次のように書くことができます．
x_range = np.arange(1, MAX_ITERS + 1)
plt.scatter(x_range, f_opts)
plt.hlines([2.811766447134533], 1, MAX_ITERS + 1, "red", linestyles='dashed')
plt.xlabel("the number of trials of augmented Lagrangian method")
plt.ylabel(r"the objective value when $x = \hat{x}$")
plt.grid(linestyle='dotted', linewidth=1)
# plt.savefig("2.jpg") # 保存する場合はコメントを外す
おまけ
最適解がどうなっているのか，一応記述しておきましょう．元の問題を純粋に CVXPY で解いたときの最適解は次のようなベクトルとなっていました．
print("Optimal solution from CVXPY: {}".format(x.value))
Optimal vector from CVXPY: [ 2.58880942e-12 1.68008414e-12 1.65765397e-11 6.99201791e-02
7.25909486e-12 -3.87874516e-01 1.43492888e-11 -2.89455826e-11
-1.52916435e-01 -3.99941194e-01 4.63658087e-01 -1.03224670e-12
-4.80387390e-12 -2.95786459e-12 4.57547694e-01 -7.12834725e-12
5.16633759e-02 -1.28747489e-10 6.96156269e-12 -8.57859234e-11
1.59492733e-12 1.03410048e-12 3.10912764e-01 1.42863849e-11
-2.03657886e-11 -6.78293885e-13 -1.49928277e-11 -1.21466280e-11
-4.31004230e-03 -5.13022159e-01]
もちろん，厳密に $0$ とはなっていませんが，ほとんどの要素が $0$ になっている解が得られていますね．
