冷めたコーヒー

Weniger, aber besser

リーマン多様体上の最適化の初歩と Pymanopt による数値実験

数学 Python 最適化

新しい記事を書きました．Pymanopt による実装は以下の記事を参照してください．

mirucacule.hatenablog.com

はじめに

この記事では，リーマン多様体上の最適化問題について扱う．まず，リーマン多様体に関する簡単な説明を行い，通常の最適化問題と比較する形でリーマン多様体上の最適化問題を定式化する．そして，Pymanoptと呼ばれるソルバーを用いて数値実験を行う．

予備知識

リーマン多様体

「リーマン多様体」とタイトルで使っているので，多様体の定義くらいは述べておきたいと思う．しかしながら，厳密に定義しようとすると，位相空間論の知識が必要になってしまうので，厳密性にこだわらない場合は，曲がった空間を抽象化したものくらいの認識でいいと思われる．さて，定義を述べる．

定義(リーマン多様体) $C ^{\infty}$ 級多様体 $\mathcal{M}$ 上の $2$ 次の対称テンソル場 $\omega$ が，$\mathcal{M}$ の各点 $p$ において正定値であるとき，$\omega$ を $\mathcal{M}$ 上のリーマン計量(Riemannian metric)という．また，リーマン計量 $g$ がひとつ与えられた多様体 $(\mathcal{M},g)$ のことをリーマン多様体(Riemannian manifold)という．

正直，この定義だけでリーマン多様体を理解することは不可能だと思う．詳しくは『多様体入門』(松島)や『多様体の基礎』(松本)などを参照されたい．

リーマン多様体の例

リーマン多様体の最も基本的な例として，$(n-1)$ 次元球面 $S ^{n-1} := \{ x\in\R ^n \mid x ^\top x = 1 \}$ がある．他の例としては次のような空間が挙げられる．

  • $n$ 次直交群 $\mathcal{O}(n) = \{ X \in \R ^{n \times n} \mid X ^{\top} X = I_n \} \subset \R ^{n \times n}$
  • シュティーフェル多様体 $\mathrm{St}(p,t) = \{ Y \in \R ^{n \times p} \mid Y ^{\top} Y = I_p \} \subset \R ^{n \times p}$
  • グラスマン多様体 $\mathrm{Grass}(p,n)$：$\R ^n$ の $p$ 次元部分空間全体

定式化

ここでは，一般の( $\R ^n$ 上の )最適化問題とリーマン多様体上の最適化問題の定式化を行う．

一般の最適化問題

一般の( $\R ^n$ 上の )最適化問題は次のように定式化される． $$ \begin{align} \mathrm{Minimize} \quad &f(x) \\ \mathrm{subject\ to} \quad &x \in S \end{align} $$ ここに，実数値関数 $f:\R ^n \to \R$ は目的関数を表す．また，集合 $S$ は $n$ 次元ユークリッド空間 $\R ^n$ の部分空間であり，特に $S=\R ^n$ の場合，最適化問題は無制約最適化問題と呼ばれる．

リーマン多様体上の最適化問題

リーマン多様体 $\mathcal{M}$ 上の最適化問題は次のように定式化される． $$ \begin{align} \mathrm{Minimize} \quad &f(x) \\ \mathrm{subject\ to} \quad &x \in \mathcal{M} \end{align} $$ ここに，実数値関数 $f:\mathcal{M} \to \R$ は目的関数を表す．一見，制約付き最適化問題にように見えるかもしれないが， $$ \begin{align} \underset{x\in \mathcal{M}}{\mathrm{Minimize}} \quad &f(x) \end{align} $$ と同じ意味であり，リーマン多様体上で考えれば制約なしと見なすことができる．

リーマン多様体上の最適化問題の例

ここでは，リーマン多様体上の最適化問題として定式化される問題の例を挙げる．

レイリー商最小化問題

レイリー商(Rayleigh quotient)最小化問題は次のように定式化される問題である． $$ \begin{align} \mathrm{Minimize} \quad &\tilde{f} (x) = \frac{x ^\top A x}{x ^\top x} \\ \mathrm{subject\ to} \quad &x \in \R ^n - \{ 0 \} \end{align} $$ ここに，$A$ は所与の $n$ 次元実対称行列を表す．制約条件で原点 $\{0\}$ を除いているため完全には無制約ではないが，原点から十分に離れたところで考えれば無制約と見なすことができる．関数 $\tilde{f}$ は任意の非零な実数 $\lambda$ に対して，次の性質を満たす． $$ \begin{align} \tilde{f}( \lambda x) &= \frac{(\lambda x) ^\top A (\lambda x)}{(\lambda x) ^\top (\lambda x)} \\ &= \frac{\lambda ^2 \cdot x ^\top A x}{\lambda ^2 \cdot x ^\top x} \\ &= \frac{x ^\top A x}{x ^\top x} \\ &= \tilde{f}(x) \end{align} $$ これより，関数 $\tilde{f}$ の値は引数 $x$ の方向のみに依存しており，ユークリッドノルムには依存していないことが分かる． したがって，$x$ のノルムを $1$ に固定したとしても，関数 $\tilde{f}$ を最小にする $x$ には影響がない．そこで，探索領域を $x ^\top x = 1$ 上に制限すると，上記最適化問題は次のように書き換えることができる． $$ \begin{align} \underset{x\in \R ^n}{\mathrm{Minimize}} \quad &f (x) = x ^\top A x \\ \mathrm{subject\ to} \quad &x ^\top x = 1 \end{align} $$ この問題は，通常の最適化の文脈では制約付き最適化問題となる． さて，$x ^\top x = 1$ を満たすようなベクトル $x$ の集合は，$(n-1)$ 次元球面 $S ^{n-1}$ で表されるので，上記の問題は $$ \begin{align} \mathrm{Minimize} \quad &f (x) = x ^\top A x \\ \mathrm{subject\ to} \quad &x \in S ^{n-1} \end{align} $$ と定式化することもできる．これは，$(n-1)$ 次元球面上の無制約最小化問題と見なすことができる．

所与の正方行列を近似する低ランク半正定値行列の可解問題

与えられた $n$ 次元実対称行列 $A$ に対して，$A$ を最も良く近似するようなランク $k<n$ の正定値行列 $S$ を求める問題を考える．ここで，行列 $A$ と半正定値行列 $S$ との近似度は次式で測られるとする． $$ \begin{align} L_{\delta} (S, A) := \sum_{i=1} ^{n} \sum_{j=1} ^{n} H_{\delta} (s_{i,j} - a_{i,j}) \end{align} $$ ここに，$\delta > 0$ は正の定数であり，$H_{\delta}(x) := \sqrt{x ^2 + \delta ^2} - \delta$ て定義される関数であり，"pseudo-Huber loss function" と呼ばれる． この問題を，$\R ^n$ 上で定式化すると次のようになる． $$ \begin{align} \underset{S\in \R ^{n\times n}}{\mathrm{Minimize}} \quad &L_{\delta} (S,A) \\ \mathrm{subject\ to} \quad &S \succcurlyeq O , \, \mathrm{rank}(S) = k \end{align} $$ さて，ランク $k$ の $n$ 次元半正定値行列全体の集合を $\mathcal{PSD}_{k} ^{n}$ で表す．すなわち，$\mathcal{PSD}_{k} ^{n} := \{ X \in \R ^{n \times n} \mid X \succcurlyeq O ,\, \mathrm{rank}(X) = k \}$ である．$\mathcal{PSD}_{k} ^{n}$ はリーマン多様体であるため，$\mathcal{PSD}_{k} ^{n}$ 上での最適化として定式化し直すと次のようになる． $$ \begin{align} \underset{S\in \mathcal{PSD}_{k} ^{n}}{\mathrm{Minimize}} \quad &L_{\delta} (S,A) \end{align} $$

その他の例

上記で挙げた例の他にも数多くの応用先があることが知られている．特に近年では，機械学習の分野における諸問題がリーマン多様体上の最適化問題として定式化することによって効率的に最適解を求解できることもあり，理論・応用の両面から盛んに研究されているようである．

アルゴリズム

$\R ^n$ における最適化アルゴリズム(直線探索)

$\R ^n$ 上の無制約最適化問題に対するアルゴリズムを述べる．ここでは，最も基本的な反復法について述べる．

  • Step0 : 初期点 $x_0$ を設定し，$k=0$ とする
  • Step1 : 終了判定を行い，条件を満たすならば $x_k$ を最適解として出力する
  • Step2 : 探索方向 $\xi_k \in \R ^n$ とステップ幅 $\alpha_k > 0$ を決定する
  • Step3 : 暫定解を $x_{k+1} = x_k + \alpha_k \xi_k$ で更新し，Step1 へ戻る

多様体上の最適化アルゴリズム

$\R ^n$ 上の無制約最適化問題に対するアルゴリズムを拡張する形でリーマン多様体上の最適化アルゴリズムを構築する．ここで注意しなければならない点として，多様体上では Step3 における "$+$" が定義されていないということである．つまり，通常の意味での加算の結果，考えている多様体からはみ出してしまう場合があるということである．そこで，はみ出したところから多様体に引き戻す操作を行う必要があるわけだが，そのような概念を一般化したレトラクションと呼ばれる写像 $R$ を定義することができる．なお，本記事では具体的な探索方向ベクトルの求め方やレクトラクションの構成方法については触れない．

  • Step0 : 初期点 $x_0$ を設定し，$k=0$ とする
  • Step1: 終了判定を行い，条件を満たすならば $x_k$ を最適解として出力する
  • Step2 : 探索方向 $\xi_k \in \R ^n$ とステップ幅 $\alpha_k > 0$ を決定する
  • Step3 : 暫定解を $x_{k+1} = R_{x_k} (\alpha_k \xi_k)$ で更新し，Step1 へ戻る

実装

リーマン多様体上の最適化問題を Python で解くことを考える．本記事では自前で実装することはせず，Pymanoptという Python の toolbox を使って求解することを考える．

Pymanopt の概要

Pymanoptは多様体上での最適化を行うための Python toolbox であり，勾配やヘッセ行列を自動で計算することができる．ユーザは以下の3点を定義するだけで最適化を行うことができる．

  1. 最適化を行う多様体 $\mathcal{M}$ を定義する
  2. 目的関数 $f:\mathcal{M} \to \R$ を定義する
  3. Pymanoptをインスタンス化する

より詳しい解説は公式サイトや文献を参照されたい．

Pymanopt の導入

Pymanoptはpipを使うことによって簡単にインストールできる．

  • pip install --user pymanopt

また，自動微分を行うために，目的関数をAutograd Theano TensorFlowのいずれかで定義する必要がある． これらのいずれにも馴染みのない場合，Autogradを使用することが推奨されている．AutogradはNumPyをラップしているため，NumPyに馴染みのある場合，簡単に実装できると思われる．なお，Autogradはpipによってインストールできる．

  • pip install autograd

実装

多様体の定義

Pymanoptでサポートされている多様体は以下の通りである．

  • Euclidean manifold：ユークリッド空間
  • Symmetric matrices：対称行列全体の集合
  • Sphere manifold：球面
  • Complex circle：複素平面上の円 $\{x \in \mathbb{C}^n \mid |x_1| = |x_2| =\dots = |x_n| = 1\}$
  • Rotations manifold：$\mathrm{SO}(n) := \{X\in\mathbb{R}^{n\times n} \mid X ^\top X=I_n, \mathrm{det}(X)=1\}$
  • Stiefel manifold：シュティーフェル多様体
  • Grassmann manifold：グラスマン多様体
  • Oblique manifold：Manifold of $m \times n$ matrices with unit-norm columns
  • Positive-definite matrices：正定値行列全体の集合
  • Elliptope manifold：Manifold of $n \times n$ positive-semidefinite matrices with rank $k$ and unit diagonal elements.
  • Fixed rank positive-semidefinite matrices：ランク $k$ の半正定値行列全体の集合
  • Fixed rank matrices：ランク $k$ の実対称行列全体の集合
  • Product manifold：積多様体

例えば，Sphere manifold を使用したい場合，以下のようにモジュールをインポートする．

  • from pymanopt.manifolds import Sphere

その上で，多様体を次のように定義する．例えば，$3$ 次元球面を定義したい場合には次のようにする．

  • manifold = Sphere(3)

他の多様体を用いる際の具体的な定義の仕方については公式サイトを参照されたい．

目的関数の定義

Autogradを使う場合，目的関数は次のように定義すればよい．

def cost(X):
    return np.exp(-np.sum(X**2))
多様体と目的関数の統合

多様体と目的関数は Pymanopt Problem オブジェクトによって結合し，その上でPymanoptソルバーに渡す必要がある．具体的には次のようにすればよい．

from pymanopt import Problem
problem = Problem(manifold=manifold, cost=cost)
ソルバー

以上で多様体の定義と目的関数の定義を終えた．後はソルバーに渡して解かせるだけである．ソルバーのパラメータは，ソルバーをインスタンス化したときするときに指定される．以下にデフォルト値を載せておく．

  • maxtime=1000：ソルバーが実行される最大時間(秒)
  • maxiter=1000：反復回数の最大回数
  • mingradnorm=1e-6：終了判定に用いられる勾配のノルム
  • minstepsize=1e-10：終了判定に用いられるステップサイズ
  • maxcostvals=5000：目的関数を評価する最大回数
  • logverbosity=0：ソルバーの動作中にログに記録される情報のレベル．0を指定すると最終的な最適解の情報のみを返す．2を指定すると最終的な最適解の情報とログ情報を保持する辞書を返す．

具体的な実装方法は次のように記述する．

xoptimal = solver.solve(problem, logverbosity=0)
xoptimal, optlog = solver.solve(problem, logverbosity=2)

Pymanoptに実装されているソルバーを以下に示す．これらのソルバーは Manopt をベースにしている．なお，ソルバーを自前で実装する場合，Solver abstract base class から継承する必要がある．

  • Trust regions()
  • SteepestDescent()
  • ConjugateGradient()
  • NelderMead()
  • ParticleSwarm()

例えば，Trust regions を用いる場合，次のように記述する．

from pymanopt.solvers import TrustRegions
solver = TrustRegions()
Xopt = solver.solve(problem)

数値実験

Pymanoptの説明が終わったので，Pymanoptを使って数値実験を行う．ここでは，上で例として挙げた問題を実験に用いる． 以下では，必要なモジュールは次のようにインポート済みであるとする．

import autograd.numpy as np
from autograd import grad
from pymanopt import Problem
from pymanopt.solvers import TrustRegions

レイリー商最小化問題

$2$ 次対称行列 $A = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ におけるレイリー商最小化問題を扱う．実は，レイリー商最小化問題の最適化 $x ^{\star}$ は行列 $A$ の最小固有値に属する固有ベクトルであることが知られている．そのため，この問題に対する最適解は固有値および固有ベクトルを計算することによって解析的に計算することができる．そのため，アルゴリズムによって求まった最適化と比較することができる．なお，行列 $A = \begin{pmatrix} 3 & 1 \\ 1 & 3 \end{pmatrix}$ に対する固有値は，$\lambda_1 = 4 ,\, \lambda_2 = 2$ であり，$\lambda_2 = 2$ に属する固有ベクトルは，$t \left( -\frac{1}{\sqrt{2}} , \, \frac{1}{\sqrt{2}} \right) ^{\top} \; (\forall t \in \R)$ となる．

実装
from pymanopt.manifolds import Sphere
A = np.array([[3,1],[1,2]])
manifold = Sphere(2)
def cost(x):
    return np.dot(x, np.dot(A, x))
problem = Problem(manifold=manifold, cost=cost)
solver = TrustRegions(logverbosity=2)
結果

実行結果は以下のようになる．

xopt, optlog = solver.solve(problem)
Compiling cost function...
Computing gradient of cost function...
Computing Hessian of cost function...
Optimizing...
                                            f: +2.637939e+00   |grad|: 1.864309e+00
acc TR+   k:     1     num_inner:     0     f: +2.100392e+00   |grad|: 8.733960e-01   exceeded trust region
acc       k:     2     num_inner:     0     f: +2.000295e+00   |grad|: 4.859006e-02   maximum inner iterations
acc       k:     3     num_inner:     0     f: +2.000000e+00   |grad|: 7.173227e-06   maximum inner iterations
acc       k:     4     num_inner:     0     f: +2.000000e+00   |grad|: 4.440892e-16   maximum inner iterations
Terminated - min grad norm reached after 4 iterations, 0.02 seconds.

optlog の中身は次のように確認できる．

print(optlog)
{'solver': 'TrustRegions',
 'stoppingcriteria': {'maxtime': 1000,
  'maxiter': 1000,
  'mingradnorm': 1e-06,
  'minstepsize': 1e-10,
  'maxcostevals': 5000},
 'solverparams': None,
 'stoppingreason': 'Terminated - min grad norm reached after 4 iterations, 0.02 seconds.',
 'final_values': {'x': array([-0.70710678,  0.70710678]),
  'f(x)': 2.0,
  'time': 0.024547100067138672,
  'stepsize': inf,
  'gradnorm': 4.440892098500626e-16,
  'iterations': 4}}

final_values が最適解に対応しており，$x ^{\star} = \begin{pmatrix} -0.70710678 \\ 0.70710678 \end{pmatrix}$ であることが分かる．また，f(x) は最適値に対応しており，$f(x ^{\star}) = 2.0$ であることが分かる．これらは，確かに行列 $A$ の最小固有値とそれに属する固有ベクトルであることが確認できる．

低ランク半正定値行列の可解問題

例として，$3$ 次対称行列 $A = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \\ 3 & 5 & 6 \end{pmatrix}$ を最も良く近似するような正定値行列 $S$ を求める問題を扱う．なお，$S$ のランクは $2$ となるようにする．

実装
from pymanopt.manifolds import PSDFixedRank
A = np.array([[1,2,3],[2,4,5],[3,5,6]])
manifold = PSDFixedRank(A.shape[0], A.shape[0]-1)
def cost(Y):
    S = np.dot(Y, Y.T)
    delta = 0.5
    return np.sum(np.sqrt((S-A)**2 + delta**2) - delta)
problem = Problem(manifold=manifold, cost=cost)
solver = TrustRegions(logverbosity=2)
結果

実行結果は以下のようになる．

Yopt, optlog = solver.solve(problem)
Compiling cost function...
Computing gradient of cost function...
Computing Hessian of cost function...
Optimizing...
                                            f: +2.401402e+01   |grad|: 5.877665e+00
acc TR+   k:     1     num_inner:     0     f: +9.994776e+00   |grad|: 9.898422e+00   negative curvature
REJ TR-   k:     2     num_inner:     0     f: +9.994776e+00   |grad|: 9.898422e+00   exceeded trust region
acc TR+   k:     3     num_inner:     0     f: +3.883474e+00   |grad|: 1.009627e+01   exceeded trust region
REJ TR-   k:     4     num_inner:     1     f: +3.883474e+00   |grad|: 1.009627e+01   exceeded trust region
acc TR+   k:     5     num_inner:     0     f: +1.049904e+00   |grad|: 7.140704e+00   exceeded trust region
REJ TR-   k:     6     num_inner:     3     f: +1.049904e+00   |grad|: 7.140704e+00   reached target residual-kappa (linear)
acc       k:     7     num_inner:     0     f: +4.655802e-01   |grad|: 4.194963e+00   exceeded trust region
acc TR+   k:     8     num_inner:     3     f: +2.900976e-01   |grad|: 1.679291e+00   exceeded trust region
acc       k:     9     num_inner:     3     f: +2.542600e-01   |grad|: 6.877936e-02   reached target residual-kappa (linear)
acc       k:    10     num_inner:     1     f: +2.542176e-01   |grad|: 4.442510e-03   reached target residual-theta (superlinear)
acc       k:    11     num_inner:     4     f: +2.542130e-01   |grad|: 1.874113e-05   reached target residual-theta (superlinear)
acc       k:    12     num_inner:     4     f: +2.542130e-01   |grad|: 6.544584e-10   reached target residual-theta (superlinear)
Terminated - min grad norm reached after 12 iterations, 0.43 seconds.

optlog の中身は次のように確認できる．

print(optlog)
{'solver': 'TrustRegions',
 'stoppingcriteria': {'maxtime': 1000,
  'maxiter': 1000,
  'mingradnorm': 1e-06,
  'minstepsize': 1e-10,
  'maxcostevals': 5000},
 'solverparams': None,
 'stoppingreason': 'Terminated - min grad norm reached after 12 iterations, 0.13 seconds.',
 'final_values': {'x': array([[0.77928005, 0.82330077],
         [0.74491914, 1.87118565],
         [1.37193605, 2.07299432]]),
  'f(x)': 0.25421300993652163,
  'time': 0.13328790664672852,
  'stepsize': inf,
  'gradnorm': 2.2973372294669235e-11,
  'iterations': 12}}

最適解はfinal_valuesを参照すると， $$ \begin{align} Y_{\star} = \begin{pmatrix} 0.77928005 & 0.82330077 \\ 0.74491914 & 1.87118565 \\ 1.37193605 & 2.07299432 \end{pmatrix} \end{align} $$ であることが確認できる．それでは，最適解における $Y_{\star} Y_{\star} ^{\top}$ を見てみよう．

print(np.dot(Yopt, Yopt.T))
array([[1.28510159, 2.12104918, 2.77582023],
       [2.12104918, 4.05624029, 4.90093862],
       [2.77582023, 4.90093862, 6.17951397]])

いかがだろうか．与えられた行列 $A = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 4 & 5 \\ 3 & 5 & 6 \end{pmatrix}$ と見比べると確かに似通っていることが見て取れる．では，実際にこの行列のランクが $2$ になっているのか確認してみよう．

S = np.dot(Yopt, Yopt.T)
print(np.linalg.matrix_rank(S))
2

numpy.linalg.matrix_rankは特異値分解をベースに与えられた行列のランクを計算してくれる関数である．この結果より，確かに得られた行列のランクは $2$ になっていることが分かる．

おわりに

今回はリーマン多様体上の最適化を簡単に説明した上で，Pymanoptを導入し，ソルバーを用いて数値実験を行なった． ほとんどソルバーに丸投げしたため，多様体上での最適化アルゴリズムを構築する理論的な解説を行うことができなかったが，自分自身も理解が不十分であるためご容赦いただければ幸いである．自前で実装する例に関しては，ksknw 様のリーマン多様体の最適化(シュティーフェル多様体上の最適化まで)を参照されたい．本記事では詳しく触れなかったシュティーフェル多様体上の最適化やリーマン多様体上での勾配について説明されている．また，本格的にこの分野を勉強したい場合は，[2] を参照されたい．リーマン多様体上の最適化に関して体系的にまとめられている数少ない書籍である．

References

[1] Townsend, James, Niklas Koep, and Sebastian Weichwald. "Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation." The Journal of Machine Learning Research 17.1 (2016): 4755-4759.
[2] Absil, P-A., Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009.