リーマン多様体上の最適化 (シュティーフェル多様体上の最適化まで)

概要

  • リーマン多様体上での最適化のお勉強をしている.
  • 多様体上での最適化では接空間上での勾配を求めたり,多様体上へレトラクションしたりする.
  • 例として球面上で,自明な最適化とレイリー商(対称行列の固有ベクトルを求める)をやる.
  • レイリー商から複数の固有値を同時に求めようとするとシュティーフェル多様体上での最適化になる.
  • シュティーフェル多様体はわかった気がするけど,グラスマン多様体がわからんのでできてない.
  • 実装は色んなライブラリがあるけど,pytorchを使う.autogradがあればなんでもいい.

はじめに

年末からお正月で多様体上の最適化について勉強している. 今回は自分の理解のために勉強したことをまとめるが,多様体がそもそもよくわかってないレベルなので色々許してください.

やりたいのはユークリッド空間上での制約つき最適化. 制約付きの最適化はラグランジュの未定乗数法とかで解くイメージだが,ニューラルネットとかがくっついているとどうやって解くのかわからない(あるのかもしれないが). 簡単な方法としてペナルティ項として制約条件を表し,制約なし最適化として解く方法が考えられるが,制約条件を満たす保証はなかったり,パラメータ調整が難しかったりする.

ユークリッド空間上での制約付きの最適化は,制約条件を満たす全ての点があるリーマン多様体上にあるとき,リーマン多様体上での制約なし問題として解くことができる. 例えば,以下は後で使うレイリー商(を変形した式)である.この制約付きの最適化は半径1の球面上の制約なし最適化問題として解釈できる.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

ただし{\bf x} \in \mathbb{R}^n ,\ \ A \in \mathbb{R}^{n\times n}

普通の勾配法に加わる操作としては以下の2つである.

  1. ユークリッド空間で求めた勾配 \nabla _x\  f(x)多様体の接空間に射影する( \mathrm{grad}\ f(x)).
  2. 更新後の点(x_t + \xi多様体からはみ出す)を多様体上に戻す(レトラクション  R_x (x_t+\xi)).

図にすると以下のような感じ.

f:id:ksknw:20180103100907p:plain

まだあんまり理解できてないので数学的な議論を避けて,具体例をやっていく. 最適化にもニュートン法とか共役勾配法とか色々あるけど,今回は単純な最急降下法だけを考える.(学習率も固定) pytorchを使ってfloatで計算しているので,数値誤差がそれなりにある.

こちらなどを参考にして勉強した(している). 数学的な話を含めてわかりやすいように思うのでおすすめ.

以下の内容は主に,これらに書かれていることを自分の理解した範疇で書きなおして,理解のためにpython (pytorch) で実装したというもの.

またこれが有名らしい(?)けど,あまり読めてない (英語苦手...)

3次元空間中の球面の例

まずは簡単な例として球面上の最適化を考える. 最適化する関数は以下のようなもの

\mathrm{argmax}_{\bf x \in \mathbb{R}^3} f({\bf x})={\bf x}[2] \ \ s.t.\ || {\bf x} ||_2=1

解は自明で{\bf x}=[0,0,1]^\topだが,これを多様体上の最急降下法で解くことにする.

ペナルティ項としての実装

リーマン多様体上での最適化を行う前に比較手法として,制約をペナルティ項として表すことを考える. つまり以下の関数を制約なし最適化問題として解く.

\mathrm{argmin}_{\bf x \in \mathbb{R}^3} - x[2] + \gamma |(1-||x||_2)|

コードは以下.

求まった解は

([  5.31967416e-05,   3.33191572e-05,   1.02169939e+00])

で0.02ぐらいずれていて,制約を満たせていない.

収束途中の遷移は以下. 青い球が制約条件で,オレンジの点が上に向かって進んでいる.

多様体からはみだしたり戻ったりを繰り返しながらガタガタしつつ進んでいく. 最急降下法で最適解に落とすには学習率を減衰させる必要があるが,今回はやってないので,こんな感じになる.

f:id:ksknw:20180102211837p:plain

これでもこの問題に対してはある程度の解は得られるが,今回はリーマン多様体上によって制約条件を常に(数値的には)厳密に満たしつつ最急降下法を行う.

多様体上の最適化

以下の最適化問題をリーマン多様体上での最適化として解く. つまり球面をS^{n-1}とするとき,目的関数を以下のように変形する.

\mathrm{argmax}_{\bf x \in S^{n-1}} f({\bf x})={\bf x}[2]

そもそもリーマン多様体が何かあまり理解できてないが, 多様体上の任意の点の接空間でリーマン計量\lt \cdot, \cdot> (内積的なもの)が定義できていればいい(たぶん). 内積が定義できれば接ベクトルの長さが定義できるので,測地線の長さがわかるということ(なのかな?)

ということで,まずは球面S^{n-1}とするときの内積を以下のように定義する.

\lt \xi,\eta> =(\xi^\top \eta)\ \ \xi,\eta\in T_xS^{n-1}\ \ x\in S^{n-1}

(これはいわゆる普通の内積だが,後で行列多様体を考える練習として順にやる.)

接空間上での勾配を考える. ユークリッド空間上での勾配は[0,0,1]^\topであるが,これは多様体の接空間上にはない.

リーマン多様体上での点pにおける勾配(\mathrm{grad}\ f)_pは以下を満たすものとして定義される.  \xiは任意の接ベクトル.

\lt (\mathrm{grad}\ f)_p, \xi>_p = \xi^\top \nabla f

なので,多様体上での勾配は以下のようになる.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

\xi^\top x=0なのでこの勾配は上の定義を満たす.

次にレトラクションを考える.今回は球面上に点を戻せばいいので,レトラクションは単純にノルムでわるだけにする.

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

以上から,勾配法で解くプログラムは以下のようになる.

得られた解は以下. float32なのでこれぐらいの誤差は許してほしいところ.

array([  1.29403008e-04,  -1.01325495e-04,   1.00000000e+00], dtype=float32)

制約条件は途中の点でも全て同様に満たされている.

f:id:ksknw:20180102223542p:plain

レイリー商の例

レイリー商(を変形したもの)である以下の最適化問題を考える.

 \mathrm{min}_{\bf x} \ {\bf x}^\top A {\bf x}\ \ s.t. \ \ ||x||_2=1

これも上と同様にリーマン多様体(球面)上の最適化問題として解ける. ちなみにこの問題の解は Aの最小固有値に対応する固有ベクトルになる.

このとき勾配,レトラクションは上と同様.

\mathrm{grad}f\left({\bf x}\right)= \left(\mathbb{I}-{\bf xx}^\top  \right) \nabla f

 R_x(\eta)=\frac{x+\eta}{||x+\eta||_2}

プログラムは以下.

図中青がペナルティ法での点の軌跡.オレンジが多様体上で最適化した時の点の軌跡. 紫は解(固有ベクトル)を表している.

オレンジは多様体上を移動しているのに対して,青は制約条件を満たさない点を経由して最適解へと近づいていく. また,(学習率などの問題もあるが)今回の問題では多様体上での最適化のほうが圧倒的に収束が速かった.

f:id:ksknw:20180102230356p:plain

レイリー商(シュティーフェル多様体上での最急降下法)

ここまではベクトル{\bf x}の最適化について考えてきた. 次に行列 X\in \mathbb{R}^{n\times m}の最適化について考える.

解く問題は先ほどと同様に対称行列A固有ベクトルを求めるもので以下の最適化問題である.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X\ \ s.t. \ \ X^\top X=\mathbb{I}_n

ここでA固有値n個あるがここではそのうちm個を求めるときを考えている。つまりXは正方行列に限らない。

これを解きたいがトレースは行列を並び替えても値が変わらないのでこのままでは解けない。 そこで以下のBrockett cost functionを導入する.

 \mathrm{min}_{X\in \mathbb{R}^{n\times m}} \ \mathrm{trace} \ X^\top A X N \ \ s.t. \ \ X^\top X=\mathbb{I}_n

N=\mathrm{diag}(\mu_1,...,\mu_n),\ \ 0\lt\mu_1\lt...\lt \mu_n

Nによって重み付けされるので,固有値の大きさによって解が順序付けられた状態のみが最適解になる.

制約条件 X^\top X=\mathbb{I}_nを満たす行列からなる多様体をシュティーフェル多様体\mathrm{St(m,n)}という。 まずは 接空間上のリーマン計量を以下の標準内積で定義する。 \lt \xi,\eta>_X = \mathrm{trace} (\xi^\top \eta) \ \ \xi,\eta \in T_X\mathrm{St}(m,n)

(あんまりわかってないので中略) このとき,多様体上の勾配は以下(たぶん). \mathrm{grad}\ f(X) = \nabla f(X) - X \mathrm{symm}(X^\top  \nabla f(X) ), \ \ \mathrm{symm}(D) = 0.5 (D+D^\top)

次にレトラクションを考える. シュティーフェル多様体のレトラクションとして

  • 極分解
  • qr分解を使う方法

などがあるが今回はqr分解を使う方法を用いる.

以上よりプログラムは以下.

これにより2つの固有値を同時に求めることができた.

f:id:ksknw:20180102231107p:plain

おわりに

リーマン多様体上での最急降下をやった. ペナルティで実装するよりは良さそう. 接空間を求めてからどうやって接空間上の勾配を求めれば良いのかがまだわかっていない. あとはQR分解とかの計算コストがどれぐらいなのかを考えたほうが良いかもしれない.

今回はpytorchで実装したが,pymanopt (python manifold optimization) というプロジェクトがあり,このライブラリを使うとシュティーフェル多様体上での最適化などが簡単にできるっぽいので良さそう.(Tensorflowなどとも連携できる(?))

そのうち中略をちゃんと埋める. また,Brockett cost functionではなく,グラスマン多様体として最適化する方法もやっていきたい.

参考