数理解析研究所講究録
第 1829 巻 2013 年 39-53

39

リーマン多様体上の共役勾配法

およびその特異値分解問題への応用
京都大学大学院情報学研究科数理工学専攻
佐藤寛之 (Hiroyuki Sato), 岩井敏洋 (Toshihiro Iwai)

Department of Applied Mathematics and Physics,
Kyoto University
概要

行列の特異値分解は，異なるサイズの 2 つのシュティーフェル多様体の積からなる多様体の上の最適化問題に帰着される．本稿では，その問題に対する最適化アルゴリズムとして最急降下法，共役勾配法，ニュートン法を導出し，それらの特徴につぃて詳しく述べる．また，共役勾配法とニュートン法を組み合わせることによって，新しい特異値分解手法を提案する．

1

序論
$m,$

$n$

$m\geq n$ とする．
$m\cross n$ 行列
を整数とし，

$A$

は

(1.1)

$A=U_{0}\Sigma_{0}V_{0}^{T}, U_{0}\in O(m), V_{0}\in O(n), \Sigma_{0}=(\begin{array}{l}\sum_{1}0\end{array})$diagと分解される．ここで，$\Sigma_{1}=$$(\sigma_{1}, \ldots, \sigma_{n}),$$\sigma_{1}\geq\cdots\geq\sigma_{n}\geq 0$である．この分解(1.1) を $A$ の特異値分解といい，$i=1,$との列ベクトルとする．つまり，$\sigma_{i},$$v_{1},$$\ldots,$$v_{n}$

をそれぞれ$U_{0}$である．
う．これらの列ベクトル$V_{0}=(v_{1}, \ldots, v_{n})$$u_{i},$$u_{i},$$\ldots,$

$n$

を

$A$

の特異値という [4,7].

$u_{1},$$\ldots,$$u_{m}$$V_{0}$$v_{i}$

をそれぞれ $A$ の左特異ベクトル，右特異ベクトルとい$v_{i}$$U_{0}=(u_{1}, \ldots, u_{m})$を使えば $A$ の特異値分解は(1.2)$A= \sum_{i=1}^{n}\sigma_{i}u_{i}v_{i}^{T}$と書くこともできる．

特異値分解は，次の最適化問題と密接な関係がある．

問題 1.1.

maximizesubject to,と

tr $(U^{T}AVN)$ ,$U\in \mathbb{R}^{m\cross p},$$V\in \mathbb{R}^{n\cross p},$

(1.3)
$U^{T}U=V^{T}V=I_{p}$

.

(1.4)

40
$1\leq p\leq n$ であり，また $N=$ diag
ここで，$(\mu_{1}, \ldots, \mu_{p}),$$\mu_{1}>\cdots>\mu_{p}>0$である．実際，次の命題が成り立つ．
$m\cross n$ 行列 $A$ の特異値分解は (1.2) の形をしており，
命題 1.1. $m\geq n$ とする．

$\sigma_{1}\geq$...

のであるとする．
を問題 1.1 の最適解とする．このとき， と列ベクトルたちは，それぞれ の小さい方から 個の特異値に属する左特異ベクトル，$\geq\sigma_{n}\geq 0$$(U_{*}, V_{*})$$U_{*}$$A$$V_{*}$$p$

右特異ベクトルとなる．つまり，(1.5)$U_{*}=(u_{1}, \ldots, u_{p}), V_{*}=(v_{1}, \ldots, v_{p})$が成り立つ．

この命題の証明は，ラグランジュの未定乗数法による．詳しくは [9] を参照されたい．
ここで，問題 1.1 の制約条件について注目する．一つの行列変数$U$を見ると，これは集合 $\{U\in \mathbb{R}^{n\cross p}|U^{T}U=I_{p}\}$ の上を動く．ところで，この集合には多様体の構造が入ることが知られており，シュティーフェル多様体と呼ばれている．つまり，シュティーフェル多様体 $St(p, n)$ は

St $(p, n)=\{Y\in \mathbb{R}^{n\cross p}|Y^{T}Y=I_{p}\}$(1.6)$U^{T}U=V^{T}V=I_{p}$ は，と定義される．すると，問題 1.1 の制約条件$(U, V)\in St(p, m)\cross St(p, n)$ と書き換えることができる．さらに，最適化の分野では最小化問題を扱うのが普通であるから，目的関数に $-1$ を乗じて，問題 1.1 を次の最適化問題$U\in \mathbb{R}^{m\cross p},$$V\in \mathbb{R}^{n\cross p},$に変換する．

問題 1.2.

minimize $F(U, V)=-$ tr $(U^{T}AVN)$ ,(1.7)

St $(p, m)\cross$ St $(p, n)$ .

(1.8)

subject to$(U, V)\in$

以降の節では問題 1.1 の解法の導出およびその適用について議論する．その準備として，
一般のリーマン多様体上の最適化手法についても簡単に紹介する．

2

リーマン多様体上の最適化手法
本節では，[1] に従ってリーマン多様体上の最適化手法の一般論を展開する．

2.1

ユークリッド空間からリーマン多様体へ

最初に，次のユークリッド空間$\mathbb{R}^{N}$における制約条件なしの最適化問題を考える．

41

問題 2.1.

minimizesubject to$f(x)$ ,$x\in \mathbb{R}^{N}.$

制約条件なしの問題に対しては，最急降下法や共役勾配法，ニュートン法などのアルゴ

リズムが知られているが，それらのアルゴリズムは次のように共通する枠組みを持つ．
Algorithm 2.1
1: 初期点$\mathbb{R}^{N}$$x_{0}\in \mathbb{R}^{N}$

2: for $k=0,1,2,$3:4:上の最適化アルゴリズム
を選ぶ．

$\ldots$do

探索方向とステップサイズ $t_{k}>0$ を計算する．
次の点 $x_{k+1}$ を $x_{k+1}:=x_{k}+t_{k}\eta_{k}$ によって計算する．
$\eta_{k}\in \mathbb{R}^{N}$

5: end for

この枠組みにおいて，探索方向の決め方が各最適化アルゴリズムを特徴付ける．探索方向$\eta_{k}\in \mathbb{R}^{N}$は，たとえば最急降下法では$x_{k}$

における $f$ の逆勾配として，ニュートン法ではニュートン方程式の解として定める．共役勾配法については後述する．

これらの概念をリーマン多様体上に拡張する．つまり，次のリーマン多様体 $(M, g)$ 上の制約条件なしの最適化問題の解法を考えたい．

問題 2.2.

まず，反復目における探索方向

minimize$f(x)$ ,subject to$x\in M.$$k$り，$\eta_{k}\in T_{x_{k}}M$$\eta_{k}$は点 $x_{k}\in M$ における接ベクトルとして選ぶ．つまとなるようにする．

また，探索方向$\eta_{k}\in T_{x_{k}}M$とステップサイズ $t_{k}>0$ を何らかの方法で求めたとして，ユークリッド空間上のアルゴリズムにおける更新の式(2.1)$x_{k+1}:=x_{k}+t_{k}\eta_{k}$は多様体上では一般には意味をなさない．実際，多様体はユークリッド空間に埋め込まれているとは限らないので加法が一般には定義されないし，仮にユークリッド空間
埋め込まれていて$x_{k}+t_{k}\eta_{k}\in \mathbb{R}^{N}$はない．そこで，が定まるとしても，これは一般には多様体$\gamma(0)=x_{k},\dot{\gamma}(0)=\eta_{k}$なる $M$ 上の曲線$\gamma$$\mathbb{R}^{N}$に

上の点でに沿って次の点 $x_{k+1}$ を探索す$M$42
る．そのために，「探索をする上で妥当な」 曲線を定める写像 $R:TMarrow M$ が見つかれ$R_{x}:=R|_{T_{x}M}$ として，
ば，
$x_{k+1}:=R_{x_{k}}(t_{k}\eta_{k})$$R$ に課すべき条件は，なる更新式を得られる．

$R_{x_{k}}(0)=x_{k}$(2.2)$\frac{d}{dt}R_{x_{k}}(t\eta_{k})|_{t=0}=\eta_{k}$

(2.3)

かつ

であることが曲線の初期点と初速度についての考察から分かる．このことを踏まえて，レ
トラクション $R:TMarrow M$ を次のように定義する．
$R$$R_{x}:=R|_{T_{x}M}$
定義 2.1. 写像 $R:TMarrow M$ が次を満たすとき，をレトラクションという．

とする．

.
$\bullet$$R_{x}(0_{x})=x$, ここで， は $T_{x}M$ の零ベクトルである．
$0_{x}$$T_{0_{x}}T_{x}M\simeq T_{x}M$なる同一視の下で，$DR_{x}(0_{x})=id_{T_{x}M}.$

2.2

リーマン多様体上の最適化アルゴリズム$M$ 上の最適化アルゴリズム
さて，リーマン多様体 $M$ 上のレトラクション $R$ を用いて，は次のように記述される:

Algorithm 2.2 リ $-$ マン多様体 $M$ 上の最適化アルゴリズム
初期点 $x_{0}\in M$ を選ぶ．

for $k=0,1,2,$
探索方向

次の点

$\ldots$do$\eta_{k}\in T_{x_{k}}M$$x_{k+1}$とステップサイズ $t_{k}>0$ を求める．

を，$x_{k+1}:=R_{x_{k}}(t_{k}\eta_{k})$によって計算する．

end for
ステップサイズの決め方としては，アルミホの方法が有名であるが，これもリーマン多様体上に拡張することができる．
$R$
定義 2.2. $f$ をリーマン多様体 $M$ 上の目的関数とし，をレトラクションとする．スカ$\sigma\in(0,1)$ は与えられた定数とする．
ラーに与えられた点 $x\in M$ と接ベクトル $\eta\in T_{x}M$ に対して，アルミホポイントを$m$ はよって定義する．ここで，$\overline{\alpha}>0,$$\beta,$$\eta^{A}=\beta^{m}\overline{\alpha}\eta$$f(x)-f(R_{x}(\beta^{m}\overline{\alpha}\eta))\geq-\sigma\langle gradf(x), \beta^{m}\overline{\alpha}\eta\rangle_{x}.$を満たす最小の非負の整数である．
をアルミホステップサイズという．
また，$t^{A}=\beta^{m}\overline{\alpha}$

43
これらを用いて，最急降下法やニュートン法，共役勾配法などをリーマン多様体上に拡張することができるが，それらについては実際に問題 1.2 に対するアルゴリズムを導出する際に述べることにする．

3積多様体 St $(p, m)\cross$ St $(p, n)$ の幾何$St(p, m)\cross St(p, n)$$St(p, m)\cross St(p, n)$ と上の問題 1.2 に対する解法を導出するために，

目的関数 $F(U, V)=tr(U^{T}AVN)$ の幾何について議論する．St $(p, n)$ そのものの幾何については [1,3] を参照されたい．

3.1$Y\in$

接空間
St $(p, n)$ における接空間$T_{Y}$$T_{Y}$

St $(p, n)$ は

St $(p, n)=\{\xi\in \mathbb{R}^{n\cross p}|\xi^{T}Y+Y^{T}\xi=0\}$(3.1)$St(p, m)\cross St(p, n)$ 上の点 $(U, V)$ における接空間はと書けるから，$T_{(U,V)}(St(p, m)\cross St(p, n))\simeq T_{U}St(p, m)\cross T_{V}St(p, n)$(3.2)$=\{(\xi, \eta)\in \mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}|\xi^{T}U+U^{T}\xi=\eta^{T}V+V^{T}\eta=0\}$と書かれる．St $(p, n)$ はユークリッド空間$\mathbb{R}^{n\cross p}$の部分多様体であるから，$\mathbb{R}^{n\cross p}$における通常の内積からの誘導計量$\langle\xi_{1},$$\xi_{2}\rangle_{Y}$

$:=$

tr

$(\xi_{1}^{T}\xi_{2})$,$\xi_{1},$(3.3)$\xi_{2}\in T_{Y}St(p, n)$$St(p, m)\cross St(p, n)$を入れることができる．このリーマン計量を自然に拡張することで，をリーマン多様体と見なす．つまり，$\langle(\xi_{1}, \eta_{1}),$$:=\langle\xi_{1},$$(\xi_{2}, \eta_{2})\rangle_{(U_{)}V)}$$\xi_{2}\rangle_{U}+\langle\eta_{1},$ $\eta_{2}\rangle_{V}=$tr$(\xi_{1}^{T}\xi_{2})+$tr$(\eta_{1}^{T}\eta_{2})$,$(\xi_{1}, \eta_{1}), (\xi_{2}, \eta_{2})\in T_{(U,V)}(St(p, m)\cross St(p, n))$(3.4)というリーマン計量を入れる．
これを用いて接空間 $T_{(U,V)}(St(p, m)\cross St(p, n))$ への射影を次のように記述できる．

こ

こで，$T_{(U,V)}(St(p, m)\cross St(p, n))\simeq T_{U}St(p, m)\cross T_{V}St(p, n)$(3.5)に注意する．

命題 3.1. 任意の$(U, V)\in St(p, m)\cross St(p, n)$ におけるに対して，接空間 $T_{(U,V)}(St(p, m)\cross St(p, n))$ への直交射影 $P_{(U,V)}$ は，$(B, C)\in \mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$$P_{(U,V)}(B, C)=(P_{U}(B), P_{V}(C))$(3.6)$P_{U}(B)=B-U$ sym $(U^{T}B),$ $P_{V}(C)=C-V$ sym $(V^{T}C)$(3.7)

で与えられる．ここで，はそれぞれ $T_{U}St(p, m),$ $T_{V}St(p, n)$ への直交射影である．なお，$sym(B):=(B+B^{T})/2$
は行列 $B$ の対称部分を表す．

44
3.2

レトラクション

$St(p, m)\cross St(p, n)$る．

$n\cross p$$QR$ 分解に基づくものを紹介す上のレトラクションの 1 つとして，のフルランク行列 $B$ の $QR$ 分解とは

(3.8)

$B=QR, Q\in St(p, n), R\in S_{upp}^{+}(p)$上三角行列全体の集合では，対角成分がすべて正のが (3.8) のようにある [4, 7]. $qf()$ を $QR$ 分解の 成分を返す写像とする．つまり，行列$St(p,m)\cross St(p, n)$$QR$ 分解を用いて$QR$ 分解されるなら，$qf(B)=Q$ である．次のように，なるものである．ここで，$S_{upp}^{+}(p)$$p\cross p$

$B$

$Q$

上のレトラクションを定義することができる．
$R_{(U,V)}$ を $T_{(U,V)}(St(p, m)\cross St(p, n))$
命題 3.2. 任意の点 $(U, V)\in St(p, m)\cross St(p, n)$ において，から St $(p, m)\cross$ St $(p, n)$ への写像で$R_{(U,V)}(\xi, \eta)=(qf(U+\xi), qf(V+\eta))$ ,となるものと定義する．すると，$R_{(U,V)}$$(\xi, \eta)\in T_{(U,V)}(St(p, m)\cross St(p, n))$(3.9)

をすべての $(U, V)\in St(p, m)\cross St(p, n)$ にわたって集め合わせてできる $R:T(St(p, m)\cross St(p, n))arrow St(p, m)\cross St(p, n)$ は $St(p, m)\cross St(p, n)$
上のレトラクションとなる．

Proof. $(qf(U+\xi), qf(V+\eta))\in$ St$(p, m)\cross$

St $(p, n)$ は $qf$ の定義から明らか．また，(3.10)$R_{(U,V)}(0,0)=(qf(U), qf(V))=(U, V)$ .

さらに，$Dqf(Y)[\triangle]=\triangle$

が任意の $\triangle\in T_{Y}St(p, n)$ について成り立つ [1] から，(3.11)$DR_{(U,V)}(0,0)[(\xi, \eta)]=(Dqf(U)[\xi], Dqf(V)[\eta])=(\xi, \eta)$が成り立つ．したがって $R$ はレトラクションである口

この $R$ を $QR$ レトラクションと呼ぶことにする．

3.3

目的関数の勾配とヘシアン

目的関数 $F$ の $(U, V)\in St(p, m)\cross St(p, n)$ における勾配 $gradF(U, V)$ は$\langle gradF(U, V),$ $(\xi, \eta)\rangle_{(U,V)}=DF(U, V)[(\xi, \eta)],$$(\xi, \eta)\in T_{(U,V)}(St(p, m)\cross St(p, n))$(3.12)

を満たす接ベクトルとして一意に定義される．また，ヘシアン $HessF(U, V)$ は，$T_{(U,V)}(St(p, m)\cross St(p, n))$上の線形変換であり，共変微分$\nabla_{(\xi,\eta)g}radF$の $(U, V)$ における値として定義される．つまり，$HessF(U, V)[(\xi, \eta)]$$:=\nabla_{(\xi,\eta)}gradF,$$(\xi, \eta)\in T_{(U,V)}(St(p, m)\cross St(p, n))$

ここで，共変微分は St $(p, m)\cross$ St $(p, n)$ 上のレヴイチビタ接続$\nabla$.

(3.13)

によって定義される．

45
さて，ここでは $F$ は問題 1.2 の目的関数(3.14)$F(U, V)=-$ tr $(U^{T}AVN)$とする．

まずは勾配についてである．
命題 3.3. 目的関数 (3.14) の $(U, V)\in St(p, m)\cross St(p, n)$ における勾配は$gradF(U, V)=(U$ sym $(U^{T}AVN)-AVN,$ $V$ sym $(V^{T}A^{T}UN)-A^{T}UN)$(3.15)と書かれる．

Proof. $gradF(U, V)$ は$F$の $(U, V)$ におけるユークリッド勾配に直交射影すれば得られるが，直交射影$P_{(U,V)}$$F_{(U,V)}^{\neg}$

を $T_{(u,v)}(St(p, m)\cross St(p, n))$は式 (3.6) と式 (3.7) で分かっているから，$gradF(U, V)=P_{(U,V)}(F_{(U,V)})=P_{(U,V)}(-AVN, -A^{T}UN)$$=(-P_{U}(AVN), -P_{V}(A^{T}UN))$
$=$

$(U$

sym $(U^{T}AVN)-AVN,$$V$

sym $(V^{T}A^{T}UN)-A^{T}UN)$ .

よって命題は示された(3.16)
口

続いてヘシアンについてである．
命題 3.4.
$S_{1}=$$(\xi, \eta)$

を $(U, V)\in$ St $(p, m)\cross St(p, n)$ における接ベクトルとする．

sym $(U^{T}AVN),$ $S_{2}=$ sym $(V^{T}A^{T}UN)$ とする．目的関数 (3.14) の $(U, V)$ におけるヘシアンは $T_{(U,V)}$ $(St(p, m)\cross St(p, n))$ 上の線形変換であり，$HessF(U, V)[(\xi, \eta)]=(\xi S_{1}-A\eta N-U$ sym $(U^{T}(\xi S_{1}-A\eta N))$ ,$\eta S_{2}-A^{T}\xi N-V$

sym $(V^{T}(\eta S_{2}-A^{T}\xi N)))$(3.17)のように作用する．

Proof. $(U(t), V(t))$ を $(U(0), V(0))=(U, V),$
St$(p, n)$上の測地線とする．

$U(t)$と$V(t)$$\ddot{U}(0)=-U\xi^{T}\xi,$を満たす [1, 9].

$(\dot{U}(0),\dot{V}(0))=(\xi, \eta)$を満たす St $(p, m)\cross$は(3.18)$\ddot{V}(0)=-V\eta^{T}\eta$$\langle HessF(U, V)[(\xi, \eta)],$ $(\xi, \eta)\rangle_{(u,v)}$は$\frac{d^{2}}{dt^{2}}F(U, V)$の $t=0$ における値であるから，$\langle HessF(U, V)[(\xi, \eta)], (\xi, \eta)\rangle_{(U,V)}=\frac{d^{2}}{dt^{2}}F(U(t), V(t))|_{t=0}$$=$

tr $(\xi^{T}\xi U^{T}AVN+U^{T}AV\eta^{T}\eta N-2\xi^{T}A\eta N)$ .

(3.19)

46
ヘシアンは対称かつ線形だから，接ベクトル$(\xi, \eta)$と$(\zeta, \chi)$に対しては，$\langle HessF(U, V)[(\xi, \eta)], (\zeta, \chi)\rangle_{(U,V)}$$= \frac{1}{2}(\langle HessF(U, V)[(\xi, \eta)+(\zeta, \chi)], (\xi, \eta)+(\zeta, \chi)\rangle_{(U,V)}$$-\langle HessF(U, V)[(\xi, \eta)], (\xi, \eta)\rangle_{(U,V)}-\langle HessF(U, V)[(\zeta, \chi)], (\zeta, \chi)\rangle_{(U,V)})$(3.20)$=\langle P_{(U,V)}((\xi S_{1}-A\eta N), (\eta S_{2}-A^{T}\xi N)), (\zeta, \chi)\rangle_{(U,V)}$となる．よって，$HessF(U, V)[(\xi, \eta)]=P_{(U,V)}((\xi S_{1}-A\eta N), (\eta S_{2}-A^{T}\xi N))$$=$$(\xi S_{1}-A\eta N-U$

sym $(U^{T}(\xi S_{1}-A\eta N)),$ $\eta S_{2}-A^{T}\xi N-V$ sym $(V^{T}(\eta S_{2}-A^{T}\xi N)))$ .

(3.21)
これで示せた4$St(p, m)\cross St(p, n)$

口上の最適化アルゴリズム

前の節で得られた結果を用いて，問題 1.2 の解法を考える．具体的には，最急降下法．
共役勾配法ニュートン法を導出する．

4.1$St(p, m)\cross St(p, n)$ 上の最急降下法

問題 2.2 $())$ 一マン多様体における制約なし問題) では， における探索方向$x_{k}$$\triangle_{k}\in T_{x_{k}}M$

を $f$ の $x_{k}\in M$ における逆勾配として選ぶ．つまり，$\Delta_{k}=-gradf(x_{k})$ とする．すると，レトラクション $R$ を用いて更新の式は$x_{k+1}=R_{x_{k}}(t_{k}\triangle_{k})$(4.1)と書ける．なお，砺はアルミホのステップサイズである．
したがって最急降下法を問題 1.2 に対して導出すると，$-gradF(U_{k}, V_{k})=(AV_{k}N-U_{k}$ sym $(U_{k}^{T}AV_{k}N),$ $A^{T}U_{k}N-V_{k}$ sym $(V_{k}^{T}A^{T}U_{k}N))$に注意して，次のようなアルゴリズムが得られる．

(4.2)

47
Algorithm 4.1 問題 1.2 に対する最急降下法
1: 初期点 $(U_{0}, V_{0})\in St(p, m)\cross St(p, n)$ を選ぶ．

2: for $k=0,1,2,$
探索方向

3:

$\ldots$do$(\xi_{k}, \eta_{k})\in\tau_{(U,V)}$$\xi_{k}=AV_{k}N-U_{k}$$(St(p, m)\cross St(p, n))$を

sym $(U_{k}^{T}AV_{k}N),$ $\eta_{k}=A^{T}U_{k}N-V_{k}$ sym$(V_{k}^{T}A^{T}U_{k}N)$(4.3)

によって計算する．
4:

アルミホのステップサイズ $t_{k}>0$ を計算する．

5:

次の点を

$(U_{k+1}, V_{k+1})=R_{(u_{k},v_{k})}(t_{k}(\xi_{k,\eta_{k}}))$

によって計算する．ここで， は$R$

St $(p, m)\cross$ St $(p, n)$ 上のレトラクションである．
6: end for
リーマン多様体 $M$ がコンパクトで目的関数 $f$ が滑らかなら，最急降下法で得られる点列$\{x_{k}\}$

は

(4.4)

$\lim_{karrow\infty}\Vert gradf(x_{k})\Vert_{x_{k}}=0$を満たすことが知られている [1]. $St(p, m)\cross St(p, n)$ はコンパクトだから，アルゴリズム
4.1 によって得られる点列は $F$ の臨界点に収束する．ところが，実際に数値実験してみると[9], その収束は非常に遅く，問題 1.2 に対しては実用的とはいえないことが分かる．そこで，共役勾配法を考えてみる．

4.2

St $(p, m)\cross$ St $(p, n)$ 上の共役勾配法

まずユークリッド空間

習しておく．まず，$\mathbb{R}^{N}$$x_{k}\in \mathbb{R}^{N}$における問題 2.1 に対する (非線形) 共役勾配法を簡単に復における探索方向は，$\triangle_{k}\in \mathbb{R}^{N}$(4.5)$\triangle_{k}=-gradf(x_{k})+\beta_{k}\triangle_{k-1}, k\geq 1$$\triangle_{0}=-gradf(x_{0})$ とする．ここで，として選ぶ．ただし，れているが，本稿では Polak-Ribiere による$\beta_{k}\in \mathbb{R}$の選び方は多数提案さ$\beta_{k}=\frac{gradf(x_{k})^{T}(gradf(x_{k})-gradf(x_{k-1}))}{gradf(x_{k-1})^{T}gradf(x_{k-1})}$(4.6)

を採用する．

共役勾配法をリーマン多様体 $M$ 上に拡張する際に問題となる部分を見ていこう．まず，$gradf(x_{k})-gradf(x_{k-1})$ という部
魚についてであるが，(4.6) の右辺の分子において，$M$ 上では gradf $(x_{k})\in T_{x_{k}}M,$ $gradf(x_{k-1})\in T_{x_{k-1}}M$ は異な分がある．リーマン多様体
る接空間の元であるから，この 2 つの差をとることはできない．同様に，(4.5) の右辺に
$M$ 上では $gradf(x_{k})\in T_{x_{k}}M,$
おいても，となるから，やはり (4.5) の右辺は計算できない (定義されない). そこで，Smith [10] は平行移動を用いてこの問題点を解決することを提案した．ところが，平行移動は必ずしも計算できるとは限らないし，$\triangle_{k-1}\in T_{x_{k-1}}M$

48
計算できたとしてもさほど効率的ではない．そこで，Absil [1] は vector transport という概念を導入した．
定義 4.1. 写像(4.7)$\mathcal{T}:TM\oplus TMarrow TM:(\eta_{x}, \xi_{x})\mapsto \mathcal{T}_{\eta_{x}}(\xi_{x})\in TM$

が $M$ 上の vector tmnsport であるとは，任意の $x\in M$ に対して次の条件が満たされるときをいう．

1. レトラクション $R$ が存在して，(4.8)$\pi(\mathcal{T}_{\eta_{x}}(\xi_{x}))=R(\eta_{x})$が成り立つ．ここで，$\pi(\mathcal{T}_{\eta_{x}}(\xi_{x}))$

2. 任意の
8

$\xi$

。

に対して

$\in T_{x}M$は接ベクトル$\mathcal{T}_{0_{x}}(\xi_{x})=\xi_{x}$$\mathcal{T}_{\eta_{x}}(a\xi_{x}+b\zeta_{x})=a\mathcal{T}_{\eta_{x}}(\xi_{x})+b\mathcal{T}_{\eta_{x}}(\zeta_{x})$この vector transport$\mathcal{T}_{\eta_{x}}(\xi_{x})$の足を表す．

が成り立つ．

が成り立つ．

を用いて，リーマン多様体 $M$ 上の共役勾配法の探索方向を(4.9)$\eta_{k+1}=-gradf(x_{k+1})+\beta_{k+1}\mathcal{T}_{\alpha_{k}\eta_{k}}(\eta_{k})$のように決める．また，魚は(4.10)$\beta_{k}=\frac{\langle gradf(x_{k}),gradf(x_{k})-\mathcal{T}_{\alpha_{k-1}\eta_{k-1}}(gradf(x_{k-1}))\rangle_{x_{k}}}{\langle gradf(x_{k-1}),gradf(x_{k-l})\rangle_{x_{k-1}}}$とする．

さて，St $(p, m)\cross$ St $(p, n)$ 上の vector transport$\mathcal{T}$の一つは次のようにして計算できる．

$\mathcal{T}:TM\oplus TMarrow TM$ を
命題 4.1. $M=$ St $(p, m)\cross$ St $(p, n)$ に対して，$\mathcal{T}_{(\xi,\eta)}(\zeta, \chi)$$:=(\zeta-$

qf $(U+\xi)$ sym $(qf(U+\xi)^{T}\zeta),$ $\chi-$ qf $(V+\eta)$ sym $(qf(V+\eta)^{T}\chi))$(4.11)

で定義する．ここで，$(\xi, \eta),$である．このとき， は$(\zeta, \chi)\in\tau_{(U,V)}(St(p, m)\cross St(p, n))$$\mathcal{T}$St $(p,m)\cross$ St $(p, n)$ 上の vector transpo 弼である．

Proof. 一般に， がユークリッド空間の部分多様体であり，レトラクション
$M$

$R$

を持つな

らば，
$\mathcal{T}_{\eta_{x}}(\xi_{x})=P_{R_{x}(\eta_{x})}(\xi_{x}) , \eta_{x}, \xi_{x}\in T_{x}M$によって定義される$\mathcal{T}$は $M$ 上の vector transport となる[1].

(4.12),

ここで，$P_{R_{x}(\eta_{x})}$

は $M$ 上の $R_{x}(\eta_{x})$ における接空間への射影である．よって $M=St(p, m)\cross St(p, n)$ に $QR$

レト

ラクション $R((3.9))$ を付随させて考えれば，既知である直交射影の式から (4.11) は従う．

したがって，(4. 11) によって定義される$\mathcal{T}$は vector transport である．
$\square$49
また，ユークリッド空間における通常の共役勾配法ではウルフのステップサイズが用いられることが多い．そこで，リーマン多様体上にもウルフのステップサイズを拡張し，それを用いて共役勾配法を行うことが望ましいと考えられる．このことについても，詳細は [9] を参照されたい．

これで問題 1.2 に対しても共役勾配法を適用することができる．ただし紙面の都合上，詳細なアルゴリズムは割愛する．数値実験によれば，最急降下法よりかなり速く臨界点に収束することが分かる [9].

4.

$3$

$St(p, m)\cross$

St $(p, n)$ 上のニュートン法

次にニュートン法について考える．一般のリーマン多様体 $M$ 上の最適化問題 2.2 に対$x_{k}\in M$ における探索方向
するニュートン法では，$\triangle_{k}\in T_{x_{k}}M$はニュートン方程式(4. 13)$Hessf(x_{k})[\triangle_{k}]=-gradf(x_{k})$の解として決める．なお，ステップサイズは $t_{k}=1$ とすることにする．問題 1.2 の目的
関数の勾配やヘシアンは既知であるから，それらを用いてニュートン法のアルゴリズムは次のようになる．

Algorithm 4.2 問題 1.2 に対するニュ $-$ トン法
1: 初期点 $(U_{0}, V_{0})\in St(p, m)\cross St(p, n)$ を選ぶ．
2: for $k=0,1,2,$3:$\ldots$do

ニュートン方程式$\{\begin{array}{l}\xi_{k}S_{1,k}-A\eta_{k}N-U_{k} sym (U_{k}^{T}(\xi_{k}S_{1,k}-A\eta_{k}N))=AV_{k}N-U_{k}S_{1,k},(4.14)\eta_{k}S_{2,k}-A^{T}\xi_{k}N-V_{k} sym (V_{k}^{T}(\eta_{k}S_{2,k}-A^{T}\xi_{k}N))=A^{T}U_{k}N-V_{k}S_{2,k}\end{array}$を$(\xi_{k}, \eta_{k})$$\in$$T_{(U_{k},V_{k})}(St(p, m)\cross St(p, n))$sym $(U_{k}^{T}AV_{k}N)$ , $S_{2,k}=$ sym
4:$(V_{k}^{T}A^{T}U_{k}N)$について解く．ここで，$S_{1,k}$

$=$

である．

$R$ は $St(p, m)\cross$
次の点を $(U_{k+1}, V_{k+1})=R_{(U_{k},V_{k})}((\xi_{k}, \eta_{k}))$ によって計算する．ここで，St $(p, n)$ 上のレトラクションである．
5: end for

ただし，このアルゴリズムにおけるニュートン方程式 (4.14) を解くのは困難である．そ$p=1$ の場合を考える．この場合，ニュートン方程式 (4.14) はこで，,(4.15)$S_{k}\eta_{k}-(I_{n}-V_{k}V_{k}^{T})A^{T}\xi_{k}=(I_{n}-V_{k}V_{k}^{T})A^{T}U_{k}$(4.16)$S_{k}\xi_{k}-(I_{m}-U_{k}U_{k}^{T})A\eta_{k}=(I_{m}-U_{k}U_{k}^{T})AV_{k}$$S_{1,k}=S_{2,k}=$ : 亀としており，さらにという簡単な形になる．ここで，ある．

$S_{k}$

はスヵラーで

50
かつ $S_{k}^{2}I_{n}-(I_{n}-V_{k}V_{k}^{T})A^{T}(I_{m}-U_{k}U_{k}^{T})A$ が正則であれば，方程式 (4.15),
もし臨
newtoneqp12 は$\neq 0$$\eta_{k}=(S_{k}^{2}I_{n}-(I_{n}-V_{k}V_{k}^{T})A^{T}(I_{m}-U_{k}U_{k}^{T})A)^{-1}(I_{n}-V_{k}V_{k}^{T})A^{T}AV_{k}$(4.17),(4.18)$\xi_{k}=S_{k}^{-1}(I_{m}-U_{k}U_{k}^{T})A(\eta_{k}+V_{k})$と解くことができる．
こうして $p=1$ のときのニュートン法は次のように書くことができる．

Algorithm 4.3 $p=1$ のときの問題 1.2 に対するニュ $-$ トン法
1: 初期点 $(u_{0}, v_{0})\in St(1, m)\cross St(1, n)=S^{m-1}\cross S^{n-1}$ を選ぶ．

2: for $k=0,1,2,$3:$\eta_{k}$と$\xi_{k}$$\ldots$doを$\eta_{k}=(s_{k}^{2}I_{n}-(I_{n}-v_{k}v_{k}^{T})A^{T}(I_{m}-u_{k}u_{k}^{T})A)^{-1}(I_{n}-v_{k}v_{k}^{T})A^{T}Av_{k}$$\xi_{k}=s_{k}^{-1}(I_{m}-u_{k}u_{k}^{T})A(\eta_{k}+v_{k})$によって計算する．ここで，$s_{k}=u_{k}^{T}Av_{k}$4:,(4.19)(4.20)

である．

次の点を $(u_{k+1}, v_{k+1})=R_{(u_{k},v_{k})}((\xi_{k,\eta_{k}}))$ によって計算する．

5: end forニュートン法は，ヘシアンが退化していない臨界点の近傍に初期点を取れば，その臨界点に二次収束する点列を生成することが知られている．しかし，必ずしも最適解に収束するという保証はない．そこで，共役勾配法とニュートン法を組み合わせた方法を次節で提案する．

5$St(p, m)\cross St(p, n)$上の最適化手法に基づく新しい特異値分

解手法$p=1$ のときの問題 1.2 の解法から，行列
まず，の最大特異値とそれに属する特異ベクトルを求めるアルゴリズムを提案する．それは，前節で議論したように共役勾配法とニュートン法を組み合わせたものである．

$A$

51
Algorithm 5.1 $p=1$ のときの問題 1.2 に対する
1: 初期点 $(U_{0}, V_{0})\in S^{m-1}\cross S^{n-1}$

イブリを選び，パラメ $-$ タ

ド法$\epsilon>0$$k:=0$ とおく．
を決める．

2:
$(\xi_{0}, \eta_{0})=-gradF(U_{0}, V_{0})$$=(AV_{0}N-U_{0}$ sym $(U_{0}^{T}AV_{0}N),$ $A^{T}U_{0}N-V_{0}$ sym $(V_{0}^{T}A^{T}U_{0}N))$(5.1)$=(AV_{0}-U_{0}U_{0}^{T}AV_{0}, A^{T}U_{0}-V_{0}V_{0}^{T}A^{T}U_{0})$とし，$(\overline{\xi}_{0}, \eta_{0})=(\xi_{0}, \eta 0)$とする．

3: while $\Vert gradF(U_{k}, V_{k})\Vert>\epsilon$ do4:共役勾配法の反復部分を実行する．
$k:=k+1.$
6: end while
7: $(u_{0}, v_{0}):=(U_{k}, V_{k})$ とおき，
$k:=0$ とおく．
8: アルゴリズム 4.3 を実行する．
5:

しかし，実際にはもちろん一般の に対して問題の解法アルゴリズムを得たい．ただし，先に述べたように，ニュートン法をそのまま実行するのは $p=1$ の場合を除いて困難
だから，問題を分割してアルゴリズム 4.3 ($p=1$ のニュートン法) を 回反復するとい
$p$

$p$

う方法を提案する．具体的には，まず最初に共役勾配法で最適解に十分近い点$St(p, m)\cross St(p, n)$を得ておく．

$u_{1},$$\ldots,$ $u_{p}$と$v_{1},$$\ldots,$とする．つまり，$\tilde{U}=(u_{1}, \ldots, u_{p}),\tilde{V}=(v_{1}, \ldots, v_{p})$$v_{p}$

をそれぞれ$\tilde{U}$とである．すると， や$u_{i}$$(\tilde{U},\tilde{V})\in$

V の列ベクトル$v_{i}$

は，大きい方番目の特異値に属する特異ベクトルやに十分近い．すると，関数 trの値は特異値と十分近くなる．そこで，を一つ固定するたびに先ほどの $p=1$ の場
合のニュートン法を点価， ) に適用すると，得られる点列はに二次収束する．こから$i$$u_{i}^{*}$$i$$\sigma_{i}$$(u_{i}^{*}, v_{i}^{*})$$v_{i}$

うしてニュートン法でべて$U_{*}$と$V_{*}$$(u_{i}^{T}Av_{i})$$v_{i}^{*}$$u_{i}^{*}$と$v_{i}^{*}$

を $i=1,$ $\ldots,p$ について得ることができれば，それらを並を作る．つまり，$U_{*}=(u_{1}^{*}, \ldots, u_{p}^{*}),$$V_{*}=(v_{1}^{*}, \ldots, v_{p}^{*})$とする．特異値が互$U_{*}\in St(p, m)$ といに異なるとき，特異ベクトルは互いに直交するから，言える．アルゴリズムは次の通り．

$V.$$\in St(p, n)$が52
Algorithm 5.2 問題 1.2 に対するハイブリッド法$k:=0$ と
1: 初期点 $(U_{0}, V_{0})\in St(p, m)\cross St(p, n)$ を選び，パラメ $-$ タ $\epsilon>0$ を決める．
おく．

2:$(\xi_{0}, \eta_{0})=-gradF(U_{0}, V_{0})$$=$$(AV_{0}N-U_{0}$とし，$(\overline{\xi}_{0,\overline{\eta}0})=(\xi_{0,\eta 0})$sym $(U_{0}^{T}AV_{0}N),$ $A^{T}U_{0}N-V_{0}$ sym$(V_{0}^{T}A^{T}U_{0}N))$

(5.2)

とする．

3: whiledo
共役勾配法の反復部分を実行する．
4:$\Vert gradF(U_{k}, V_{k})\Vert>\epsilon$$k:=k+1.$

5:

6: end while$:=(U_{k}, V_{k})$ とおく．
7:$(\tilde{U},\tilde{V})$s: for $i=1,2,$9:$\ldots$

,

$p$

do

$(u_{0}, v_{0});=(\tilde{U}_{i},\tilde{V}_{i})$$k:=0$とおき，

とおく．ここで， と は，それぞれ$\tilde{U}_{i}$$\tilde{V}_{i}$$\tilde{U}$と番目の列ベクトルである．
アルゴリズム 4.3 を実行する．

$\tilde{V}$

10:

の

$i$

11: end for6終わりに行列の特異値分解を，積多様体 St $(p, m)\cross$ St $(p, n)$ 上の最適化問題に帰着させ，その解法アルゴリズムを導出することによって，新たな特異値分解手法を得た．本稿では数値実験については扱わなかったが，詳細は [9] にある．その実験結果も踏まえて提案手法は次のような利点を持つ．まず一つには，通常の特異値分解手法に必要な前処理を必要としない．あえて言えば，アルゴリズム 5.2 において，共役勾配法で近似解を得るところが前処理と見なせるかもしれない．見方を変えれば，収束性の要はあくまでニュートン法であり，近似解を求めるところは共役勾配法でやらずとも良い．たとえば MATLAB の svd 関数
で求めた特異値分解を初期点としてニュートン法を適用することで，より良い解を得るということもできる．これが第二の利点であり，精度の高い特異値分解を実現することがで
きる．

また，行列

$A$

の特異値が縮退している場合についてのアルゴリズムの振る舞いについても詳しく調べる必要があるが，そちらについては [9] を参照されたい．そのような場合
でも，提案手法で特異値分解は達成されることが分かっている．
なお，最後に，本稿で扱ったリーマン多様体上の共役勾配法の収束性についてはさほど一般的な研究がなされておらず，詳しく解析されるべき事柄である．

53
参考文献
[1]

.-A. ABSIL, R. MAHONY, AND R. SEPULCHRE, optimization Algorithms on Matrix
Manifolds, Princeton University Press, Princeton, $NJ$ , 2008.
$P$[2] R. L. ADLER, J. P. DEDIEU, J. Y. MARGULIES, M. MARTENS, AND M. SHUB,
Newton’s method on Riemannian manifolds and a geometric model for the humanspine, $IMA$ J. Numer. Anal., 22(3) (2002), pp. 359-390.
[3] A. EDELMAN, T. A. ARIAS, AND S. T. SMITH, The geometry of algorithms withorthogonality constraints, SIAM J. Matrix Anal. Appl., 20(2) (1998), pp. 303-353.

[4] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computations, Third ed., The Johns
Hopkins University Press, Baltimore, $MD$ , 1996.
[5] U. HELMKE, J. B. MOORE, optimization and Dynamical Systems, in: Communications and Control Engineering, Springer, London, 1994.
[6] M. R. HESTENES AND E. STIEFEL, Methods of conjugate gmdients for solving linear
systems, J. Res. Nat. Bur. Stand., 49 (1952), pp. 409-436.

[7] L. N. TREFETHEN AND D. BAU, Numerical Linear Algebra, SIAM, 1997.
[8] J. Nocedal and S. J. Wright, Numerical optimization, Springer Series in Operations
Research. Springer-Verlag, New York, 1999.
[9] H. SATO AND T. IWAI, A Riemannian optimization approach to the matrix singularvalue decomposition, SIAM J. Optim., to appear.

[10] S. T. SMITH, optimization techniques on Riemannian manifolds, in Hamiltonianand gradient flows, algorithms and control, Fields Inst. Commun., 3. Amerl Math.
Soc., Providence, $RI$ , 1994, pp. 113-136.