二次関数の凸性
原点を通る二次関数 $f(\mathbf{x}) = \left< \mathbf{x}, \mathbf{Ax} \right>$ は行列 $\mathbf{A}$が半正定値対称であるとき凸であることについて以下の記事で述べた.
証明の技法としてはシンプルで,以下の事実を基に式変形を行うだけである.
証明ではスカラー $\lambda$ として $0<\lambda < 1$ として取っているが,$\lambda = 0,\, 1$ のときは自明に $f(\lambda \mathbf{x} + (1-\lambda )\mathbf{y}) = \lambda f(\mathbf{x}) + (1 - \lambda ) f(\mathbf{y}) $ が成立するため $0<\lambda < 1$ の場合についてのみ検討すれば十分である.
さて,一般の二次関数は行列 $\mathbf{A}\in\mathbb{R}^{n\times n}$,ベクトル $\mathbf{b}\in\mathbb{R}^n$,スカラー $c\in\mathbb{R}$ によって
$$ f(\mathbf{x}) = \left< \mathbf{x}, \mathbf{Ax} \right> + 2 \left< \mathbf{b}, \mathbf{x} \right> + c $$
と表すことができる.これは,次のように変形することができる:
$$ f(\mathbf{x}) = \begin{pmatrix} 1 & \mathbf{x}^\top \end{pmatrix} \begin{pmatrix} c & \mathbf{b}^\top \\ \mathbf{b} & \mathbf{A} \end{pmatrix} \begin{pmatrix} 1 & \mathbf{x} \end{pmatrix}. $$
したがって,行列 $\begin{pmatrix} c & \mathbf{b}^\top \\ \mathbf{b} & \mathbf{A} \end{pmatrix}$ が半正定値行列であれば,この関数は凸であることが言える. → 行列 $\mathbf{A}$ が半正定値であれば $ \left< \mathbf{x}, \mathbf{Ax} \right>$ は凸であり,$\mathbf{x}$ に関する線型項 $ 2 \left< \mathbf{b}, \mathbf{x} \right> + c$ も凸であるので,実数値関数 $f$ は $\mathbf{A}$ が半正定値でありさえすれば凸であることが言える.
非線型一変数関数の最小化(SciPy 使用)
これは何?
- 一変数関数 $f\colon\mathbb{R}\to\mathbb{R}$ を与えられた閉区間 $\Omega \subseteq \mathbb{R}$ 上で最小化する $x^\star$ を求める方法について扱う
- 最適化問題として記述すると以下のように表される: \begin{align} \min \quad &f(x) \\ \mathrm{s.t.} \quad &x \in \Omega. \end{align}
- アルゴリズムの詳細には一切触れず,
scipy.optimize
モジュールに含まれる以下の関数を用いる:minimize_scalar
shgo
minimize
differential_evolution
dual_annealing
- 詳しい使い方は公式ドキュメントを参照されたい
実装
適当な関数 $f$ が必要なので,ここでは以下の非線型関数を扱う.
$$ f(x) = \frac{3}{7} \cos\left(\frac{3}{2} x\right) - \sin\left(\frac{1}{3} x\right) + \frac{1}{3} \sin(x^2) - 1 $$
制約条件は,$\Omega = [0, 2\pi]$ とする.グラフの概形は以下の通り.
モジュールの import
import numpy as np from scipy.optimize import minimize_scalar, shgo, minimize import matplotlib.pyplot as plt from scipy.optimize import differential_evolution, dual_annealing
minimize_scalar
res = minimize_scalar(f, bounds=(0, 2*np.pi), method="bounded") print(res.x) print(res.fun)
- 結果は
res.x
で最適解,res.fun
で最適値がそれぞれ得られる - 制約付きのため
method="bounded"
を指定 - 公式ドキュメント
shgo
res = shgo(f, bounds=[(0, 2*np.pi)]) print(res.x) print(res.fun)
minimize_scalar
と制約の与え方が異なるので注意- 公式ドキュメント
minimize
div = 100 x = np.linspace(0, 2*np.pi, div) y = f(x) x0 = (2 * np.pi / div) * (np.argmin(y) + 1) # initial point bnds = ((0, 2 * np.pi),) res = minimize(f, x0, method='Nelder-Mead', bounds=bnds, tol=1e-6) print(res.x) print(res.fun)
minimize
には初期点が必要- 初期点は閉区間 $\Omega=[0, 2\pi]$ を 100 分割し,その中で最も関数値が小さいところを採用
- 制約を
((0, 2 * np.pi),)
のように与えているが最後の,
は誤字ではない - 公式ドキュメント
differential_evolution
res = differential_evolution(f, bounds=[(0, 2*np.pi)], disp=False) print(res.x) print(res.fun)
dual_annealing
lw = [0] up = [2 * np.pi] ret = dual_annealing(f, bounds=list(zip(lw, up))) print(res.x) print(res.fun)
- Dual Annealing optimization
- 公式ドキュメント
結果の比較
結果は以下の通り.
- いずれの場合でも(局所的)最適解が得られている
minimize
を使った場合のみ大域的最適解が得られているminimize_scalar
およびshgo
は局所的最適解に収束しているminimize
は初期点を大域的最適解の近傍で生成しているので上手く大域的最適解に収束しているdifferential_evolution
およびdual_annealing
も大域的最適解が得られているが,あくまでヒューリスティック解法なので収束先が異なる場合があり得る
追記記録
differential_evolution
とdual_annealing
の記述を追加
フルランク行列をランダムに生成する手法(Python による実装付き)
フルランクな行列をランダムに生成する手法を教えて頂きました:
https://t.co/mp7HaadpwZ
— とみや🐉これならわかる機械学習 (@TomiyaAkio) August 8, 2021
ランダムなゼロでない固有値を対角にならべてから、ランダムなユニタリー行列で相似変換すれば良さそうです。ランダムなユニタリー行列はランダムなエルミート行列から作れます。
実装
手順:
- ランダムに対角行列 $D$ を生成する
- ランダムな直行行列 $Q$ を生成する
- 行列 $D$ を $Q$ によって相似変換する(i.e., $A = Q^\top D Q$)
from scipy.stats import ortho_group
def make_full_rank_matrix(d, seed=0):
"""
>>> make_full_rank_matrix(3)
array([[ 0.59537814, 0.03454805, -0.05188665],
[ 0.03454805, 0.6643958 , -0.04420196],
[-0.05188665, -0.04420196, 0.60699231]])
"""
np.random.seed(seed)
D = np.diag(np.random.rand(d))
Q = ortho_group.rvs(dim=d)
A = Q.T@D@Q
return A
# test for whether the matrix made by the above method is full rank or not
num_iter = 10
success = 0
for seed in range(num_iter):
d = 500
S = make_full_rank_matrix(d, seed)
if np.linalg.matrix_rank(S) == d:
success += 1
print(success)
別の手法
次の手法で生成される行列を考える:
- $n\times k$ 次元の行列 $W$ をランダムに生成する
- 対角成分が正であるような対角行列 $D$ をランダムに生成する
- 行列 $B$ を $B = WW^\top + D$ で生成する
- 行列 $B$ を相似変換して行列 $A$ を得る,i.e., $A = \text{diag}[1/\sqrt{B_{ii}}] B \text{diag}[1/\sqrt{B_{ii}}]$.
実装:
def make_full_rank_matrix(d, k, seed=0):
"""
>>> make_full_rank_matrix(3, 2)
array([[1. , 0.62140315, 0.74064102],
[0.62140315, 1. , 0.70247127],
[0.74064102, 0.70247127, 1. ]])
"""
np.random.seed(seed)
W = np.random.rand(d, k)
B = W@W.T + np.diag(np.random.rand(d))
A = np.diag(1/np.sqrt(np.diag(B))) @ B @ np.diag(1/np.sqrt(np.diag(B)))
return A
# test for whether the matrix made by the above method is full rank or not
num_iter = 10
success = 0
for seed in range(num_iter):
d = 500
S = make_full_rank_matrix(d, seed)
if np.linalg.matrix_rank(S) == d:
success += 1
print(success)
この手法で生成される行列 $A$ の固有値は,$k$ の値によって特徴付けられる.図を見た方が早いと思うので,以下に $k = 100, 50, 20, 10, 5, 2$ のそれぞれの場合における固有値を大きい順に並べて図示した 図を載せておく(行列のサイズは $100\times 100$ である):
図のように,上記の手法で生成した $n\times n$ の行列 $A$ の固有値を大きい順に並べたとき,$k$ 番目に大きい固有値までは大きい値をとり,それ以降の固有値は微小の値をとる($k=2$ のときはそうでもないように見えるが...).
Reference
『妄想する頭 思考する手 想像を超えるアイデアのつくり方』を読んだ感想
思いかけずツイートが伸びてしまいました.
「M2 病」の話が面白かった。即ち、大学院に入学したての頃は無邪気で柔軟に発想することができるが、修士 2 回くらいになると先行文献などを読み漁って知識が付いた分、自分に残されたフロンティアなど存在しないと思ってしまい、ある種のスランプ状態に陥ってしまう病。わたしは万年これ…。 pic.twitter.com/JiR1Tm271Z
— みるか (@mirucaaura) July 25, 2021
折角なので,読んだ感想を認めておきたいと思います.こういうのはシュッと読んでシュッと記録しておくに限ります.
感想
本書を読んで,以下の 3 点が印象に残りました:
- 「天使度」と「悪魔度」
- 妄想の言語化: クレーム
- 見る前に跳ぶことの効用
「天使度」と「悪魔度」
本書を通して「天使度」と「悪魔度」という指標が多様される.これは研究テーマの良し悪しを評価するための尺度で暦本研で使われているらしい.「天使度」は発想の大胆さを表す尺度であり,「悪魔度」は技術レベルの高さを表す尺度として使われている.すなわち,人をポカンとさせるようなアイディアであれば天使度は高いし,実現するのに必要な技術レベルが高ければ悪魔度が高いことになる.この 2 軸のバランスがちょうどよい研究テーマは誰にでもその価値が明らかであり,かつ,破壊力のあるイノベーションとなり得る.
とかく技術競争の激しい現代では,悪魔度を高めることに走りがちであるように感じる.自分の研究の方向性としてもそういう風になっている傾向がある.確かに技術度は高くてパッと見すごいんだけど,「結局,何が嬉しいの?」と言われる典型で,研究それ自体の面白さが伝わらない.もちろん,分野にも依拠するので,それが悪いとは言い切れないが,天使度を高められるような柔軟な発想もできるように視野広く活動していきたい.
妄想の言語化: クレーム
クレームとは「私はこの研究ではここを主張します」という言明のこと.例えば,「DNA は二重螺旋構造をしている」などが例として挙げられている.クレームは必ずしも正しくある必要はなく(その時点で正しいかどうかなんて分からないので),検証可能な仮説となっていればよい1.その際に,ダラダラと長い文章で表現するのではなく,できるだけ短く一行で言い切ったものであるべきと著者は説いている.思考に耽っていると,様々なアイディアが頭の中でモヤモヤと無限に広がっていくことが度々ある.それを放置せずに,モヤモヤした思考の中からクレームとして切り出せることは何なのかと考え,アイディアを洗練させていくことを習慣にしていきたい.
見る前に跳ぶことの効用
「良いアイディアを思い付いたらさっさと手を動かせ」ということ.熟考してばかりでは事態は何も前進しないが,手を動かせばたとえ失敗したとしても熟考の何倍もの発見が期待できる.本書では「GAN2」を考案したグッドフェローの話を例として取り上げている.GAN は彼が仲間と一緒に夕飯をとりながら話をしているときに閃いたらしい.そして,食事を終えるとすぐにアイディアを検証し,その翌日に論文を書いて発表している3.
「このアイディアは面白そうだけど,本当にうまくいくだろうか」などと考えて,重い腰が上がらずに検証に時間がかかってしまうことがあったりするが,必要なのは熟考ではなくシュッと手を動かして実験してみることというのがよく分かる話だと思った.
Refference
CVXPY を用いた等式制約付き最適化問題の求解
以前に 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$ になっている解が得られていますね.