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

はじめに
特異値分解
特異値分解と最適化問題
リーマン多様体上での特異値分解
$\mathrm{St}(p,m)\times\mathrm{St}(p,n)$ 上の最適化
リーマン多様体上での共役勾配法
Pymanopt による求解
おわりに
参考文献
おまけ

はじめに

以前（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 から Stiefel と Product をインポートする．また，今回は最適化手法として共役勾配法を用いるため，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

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

Euclidean Space: Euclidean
Sphere Manifold: Sphere
Grassmann Manifold: Grassmann
Complex Circle: ComplexCircle
Special Orthogonal Group: SpecialOrthogonalGroup
Oblique Manifold: Oblique
Symmetric Positive Definite Matrices: SymmetricPositiveDefinite
Positive Semidefinite Matrices: PSDFixedRank
Fixed-Rank Matrices: FixedRankEmbedded
Positive Matrices: Positive

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

Nelder-Mead Algorithm: NelderMead
Particle Swarms: ParticleSwarm
Steepest Descent: SteepestDescent
Riemannian Trust Regions Algorithm: TrustRegions

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

問題を定義する方法について述べる．今回は行列 $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})$ を最急降下法と共役勾配法で解いたときの勾配ノルムの推移を描いた結果を載せておく．

冷めたコーヒー

Weniger, aber besser