数理解析研究所講究録
第 1981 巻 2016 年 127-142

127

シュティーフエル多様本上のニュートン法とその収束性解析
東京理科大学工学部経営工学科

佐藤寛之

東京理科大学理学部数理情報科学科

相原研輔

Hiroyuki Satol
Department of Management Science, Tokyo Univerisity of Science

Kensuke Aihara
Department of Mathematical Information Science, Tokyo University of Science
概要

ユークリッド空間における最適化アルゴリズムは多数あるが，その中でも制約条件なしの問題に対しては，最急降下法やニュートン法，あるいは共役勾配法など，性質の良いものが知られている．一方，制約条件付きの最適化問題に対しては，これらのアルゴリズムにより生成される点列は制約条
件を満たさないため，制約条件付き問題に特化した他の最適化手法が用いられることが多い．しかし，制約条件付きの最適化問題であっても，その制約条件を満たす点全体がリーマン多様体をなす場合には，ニュートン法などをその多様体上に一般化したアルゴリズムを用いることができる．
本稿では，特異値分解と密接な関連を持つシュティーフェル多様体上の最適化問題を扱い，この問題に対するニュートン法について詳しく議論する．特に，ニュートン法の各反復において現れるニュートン方程式が線形方程式であることに着目し，この方程式に対する線形反復解法を提案する．また，ニュートン方程式を求解する際の精度とニュートン法自体の収束率の関係についても議論する．

序論

1

本稿では次の最適化問題を考える．

maximizesubject to$tr(U^{T}AVN)$ ,$U^{T}U=V^{T}V=I_{p},$(1.1)$U\in \mathbb{R}^{m\cross p},$

ここで，正の整数 $m,$ $n,p$ は $m\geq n\geq p$ を満たし，$V\in \mathbb{R}^{n\cross p}.$$A\in \mathbb{R}^{m\cross n},$$N=diag(\mu_{1}, \mu_{2}, \ldots, \mu_{p})$,$\mu_{1}>\mu_{2}>\cdots>\mu_{p}>0$ である．

数値線形代数の観点からは，問題 (1.1) を解くことは $A$ の打ち切り特異値分解 [4] を求めることと等価である．すなわち，問題 (1.1) の最適解を
の列ベ
とすると，
$(U_{*}, V_{*})$$1_{e}$-mail:

hsato[AT]ms.kagu.tus.ac.jp$U_{*},$$V_{*}$

128
クトルがそれぞれ $A$ の大きい方から$p$個の特異値に属する左特異ベクトル，右特異ベクトルであることが示されている [7].

最適化の観点からは，この問題はユークリッド空間

における $U^{T}U=V^{T}V=$$\mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$という制約条件付きの最適化問題であり，このままではニュートン法や共役勾配法などの無制約最適化手法 [5] を適用することができない．しかし，シュテイーフエル多様体と$I_{p}$

呼ばれる多様体$St(p, n)=\{Y\in \mathbb{R}^{n\cross p}|Y^{T}Y=I_{p}\}$を考えると，問題 (1.1) の制約条件は $U\in St(p, m)$ , $V\in St(p, n)$ と書けるから，問題 (1.1)
は次のように積多様体 $St(p, m)\cross St(p, n)$ 上の無制約最適化問題と見なすことができる．

minimize $F(U, V)=-tr(U^{T}AVN)$ ,subject to$(U, V)\in St(p, m)\cross St(p, n)$

(1.2)

.

本稿では，[7] で提案された，この問題に対するニュートン法を紹介するとともに，一般の に対するニュートン方程式の難しさについて述べ，それを解決する工夫について議$p$

論する．具体的には，複雑な行列方程式の形で表されるニュートン方程式を等価な連立一次方程式に変換し，クリロフ部分空間法の一種である共役残差法 [3] を適用する．さらに，
ニュートン法の各反復におけるニュートン方程式の解の精度とニュートン法の収束率の関係についても議論する．最後に，提案手法の収束性を数値実験により検証する．

2

シュティーフエル多様体上のニュートン法
ユークリッド空間$\mathbb{R}^{N}$における制約条件なしの問題

minimizesubject to$f(x)$ ,$x\in \mathbb{R}^{N}$に対するニュートン法では，最適解の十分近くの初期点$\eta\in \mathbb{R}^{N}$$x_{0}$

から始め，

$k$

回目の反復ではに対するニュートン方程式(2.1)$\nabla^{2}f(x_{k})\eta=-\nabla f(x_{k})$を解き，その解を$\eta_{k}$として次の点 $x_{k+1}$ を(2.2)$x_{k+1}=x_{k}+\eta_{k}$の勾配，ヘッセ行列であり，は探索方向と呼ばれる．この手法を $St(p, m)\cross St(p, n)$ 上に拡張し，問題 (1.2) に適用でによって計算する．ここで，$\nabla f$および$\nabla$2f はそれぞれ$f$$\eta_{k}$

きるようにしたい．ここでは，[7] で提案されたアルゴリズムを紹介する．
まず，ニュートン方程式 (2.1) に現れる勾配やヘッセ行列は目的関数のみの情報しか有していないため，これらを多様体 St $(p, m)\cross St(p, n)$ 上に拡張する必要がある．そのために，
リーマン計量を多様体 $St(p, m)\cross St(p, n)$ に導入する必要があるが， $St(p, m)\cross St(p, n)$129の標準内積から誘導される計量の部分多様体であるため，$St(p,m)\crossSt(p,n)$
上の各点 $(U, V)$ における接空間を考えることができる．すなわち，は$\mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$$\mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$$T_{(u,v)}(St(p, m)\cross St(p, n))$において，, $(\xi’, \eta’)\rangle_{(U,V)}=tr(\xi^{T}\xi’)+tr(\eta^{T}\eta’)$ ,, $(\xi’, \eta’)\in T_{(U,V)}$ $(St(p, m)\cross St(p, n))$(2.3)
によって定まる内積を定義し，これをすべての点 $(U, V)$ について集めた族を$St(p, m)\cross St(p, n)$ のリーマン計量とすることができる．このリーマン計量を与えることで，St $(p, m)\cross St(p, n)$ はリーマン多様体となり，その幾何学的構造がもたらす情報を含むもの
として目的関数 $F$ の勾配 grad F やヘシアン Hess F を定義することができる．実際，リーマン計量により点 $(U, V)$ における St $(p, m)\cross St(p, n)$ の法空間を定義することができ，接空間$T_{(u,v)}(St(p, m)\cross St(p, n))$ への直交射影 $P_{(U,V)}$ :$\langle(\xi, \eta)$$(\xi, \eta)$$\rangle_{(U,V)}$$\mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}arrow T_{(U,V)}(St(p, m)\cross St(p, n))$が$P_{(U,V)}(X, Y)=(X-Usym(U^{T}X), Y-Vsym(V^{T}Y))$により定まる．ここで， $sym(.$ は行列の対称部分を返す写像，すなわち $sym(Z)=(Z+$m p$Z^{T})/2$ である．目的関数 $F$ の $(U, V)\in St(p, m)\cross St(p, n)$ における勾配 grad F はにおけるユークリッド勾配を直交射影することで計算でき，$)$$\cross$$\mathbb{R}$$\cross$$\mathbb{R}^{n\cross p}$

grad F $(U, V)=P_{(U,V)}(\nabla\overline{F}(U, V))$

(2.4)

上で $\overline{F}(U, V):=-tr(U^{T}AVN)$ と定義される関数となる．ここで， は$\overline{F}|_{St(p,m)\crossSt(p,n)}=F$である．すなわち，である．さらに，勾配 grad F の微分をさらに射影することで， $(U, V)\in St(p, m)\cross St(p, n)$ におけるヘシアン Hess F $(U, V)$ の $(\xi, \eta)\in$$\tau_{(u,v)(St}(p, m)\cross$ St $(p, n))$ の作用を$\overline{F}$$\mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$

Hess F $(U, V)[(\xi, \eta)]=P_{(U,V)}(D($ grad F $)(U, V)[(\xi, \eta$(2.5)と表すことができる．

(2.4) および (2.5) から勾配とヘシアンを具体的に計算すると， $(U, V)\in$ St $(p, m)\cross$
St $(p, n)$ , $(\xi, \eta)\in T_{(u,v)}(St(p, m)\cross St(p, n))$ に対して，

grad F $(U, V)=(US_{1}-AVN, VS_{2}-A^{T}UN)$ ,Hess F $(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$(2.6)となる．ここで，$S_{1}=sym(U^{T}AVN) , S_{2}=sym(V^{T}A^{T}UN)$である．これらを用いると，ユークリッド空間上のニュートン法で現れるニュートン方程式 (2.1) を St $(p, m)\cross St(p, n)$ 上に拡張した方程式

Hess F $(U_{k}, V_{k})[(\xi_{k}, \eta_{k})]=-$ grad F $(U_{k}, V)$130を具体的に書き下すことができる．簡単のため反復回数$k$を省略すると，問題 (1.2) に対するニュートン方程式は(2.7)$\{\begin{array}{l}\xi S_{1}-A\eta N-U sym (U^{T}(\xi S_{1}-A\eta N))=AVN-US_{1},\eta S_{2}-A^{T}\xi N-V sym (V^{T}(\eta S_{2}-A^{T}\xi N))=A^{T}UN-VS_{2}\end{array}$と書ける．
一般の に対してニュートン方程式 (2.7) を解くのは一見すると困難で，それを解決$p$するための工夫についての議論が本稿の主題の一つであり，次節以降で詳しく述べる．
におけるニュートン方程式が解けたとする．その解をここでは，第 反復で点$(\xi_{k}, \eta_{k})\in T_{(U_{k},V_{k})}(St(p, m)\cross St(p, n))$ とすると，ユークリッド空間における点列の更新式(2.2) を $St(p, m)\cross St(p, n)$ 上に拡張しておけば，ニュートン法の 1 反復が完成する．一般が定義する半直線上から探索方向にユークリッド空間における最適化手法では，点に次の点を定めることになるが，多様体 $St(p, m)\cross St(p, n)$ 上では半直線の代わりに，現に接する曲線を考えれば良い．そこで，行列をにおいて探索方向在の点$k$$(U_{k}, V_{k})$$x_{k}$$(U_{k}, V_{k})$$\eta_{k}$$(\xi_{k}, \eta_{k})$引数としてその QR 分解の $Q$ 成分を返す写像 $qf()$ 用いて$(U_{k+1}, V_{k+1})=(qf(U_{k}+\xi_{k}), qf(V_{k}+\eta_{k}))$とすれば適切な更新式になることが示されている [1, 7]. これは， $(U_{k}, V_{k})+(\xi_{k}, \eta_{k})=$を $St(p, m)\cross St(p, n)$ に引き戻すような写像，すなわち
レトラクションと呼ばれる写像を考えていることになる．なお，一般のリーマン多様体上の最適化については [1] が詳しい．
ここまでの議論をアルゴリズムの形でまとめておく．
$(U_{k}+\xi_{k}, V_{k}+\eta_{k})\in \mathbb{R}^{m\cross p}\cross \mathbb{R}^{n\cross p}$

アルゴリズム 1 問題 (1.2) に対するニュートン法
1: 初期点 $(U_{0}, V_{0})\in St(p, m)\cross St(p, n)$ を選ぶ．
2: for $k=0$ , 1, 2, . . . do
3$\cdot$:ニュートン方程式$\{\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},\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))$について解く．ここで，$S_{1_{\}}k}$sym $(U_{k}^{T}AV_{k}N)$ , $S_{2,k}=sym(V_{k}^{T}A^{T}U_{k}N)$ である．
4:

次の点を $(U_{k+1}, V_{k+1})=(qf(U_{k}+\xi_{k}), qf(V_{k}+\eta_{k}))$ により計算する．

5: end for3

シュティーフエル多様体上のニュートン方程式の解法
本節では問題 (1.2) に対するニュートン方程式 (2.7) の解法について議論する．

$=$

131
ニュートン方程式の標準的な線形方程式への変換

3.1

方程式 (2.7) はそれ自体が複雑な形をしているが，さらに $(\xi, \eta)\in T_{(U,V)}(St(p, m)\cross$

St $(p, n))$ という条件があるため，これも考慮する必要がある．シュティーフェル多様体
St $(p, n)$ 上の点 における接空間 $T_{Y}St(p, n)$ は$Y$$T_{Y}St(p, n)=\{\xi|\xi^{T}Y+Y^{T}\xi=0\}$と表せるため [2],$(\xi, \eta)\in T_{(U,V)}(St(p, m)\cross St(p, n))$

は

$\xi$

TU UT $=\eta^{T}V+V^{T}\eta=0$ を
$+$

$\xi$

満たさなければならない．この条件を扱いやすくするため，本稿では [2] で示されている接空間の別の表現(3.1)$T_{Y}St(p, n)=\{YB+Y_{\perp}C|B\in Skew(p), C\in \mathbb{R}^{(n-p)\cross p}\}$を用いる．ここで，列である．また，Skew (p) は$Y\perp\in \mathbb{R}^{n\cross(n-p)}$

は $Y^{T}Y_{1}=0$ かつ $Y_{\perp}^{T}Y_{\perp}=I_{n-p}$ を満たす任意の行

次反対称行列全体の集合を表す．なお，$\xi\in T_{Y}St(p, n)$ に対して $\xi=YB+Y_{\perp}C$ となる $B\in Skew(p)$ ,$p$$Y_{1}$$C\in \mathbb{R}^{(n-p)\cross p}$を固定すると，が一意に存在する．

以後， $U^{T}U\perp=0,$ $U_{\perp}^{T}U\perp=I_{m-p},$ $V^{T}V\perp=0,$
よび$V\perp\in \mathbb{R}^{n\cross(n-p)}$$V_{\perp}^{T}V\perp=I_{n-p}$を満たす$U\perp\in \mathbb{R}^{m\cross(m-p)}$を固定しておく．式 (3.1) から， $(\xi, \eta)\in T_{(u,v)}(St(p, m)\cross St(p, n))$ は(3.2)$(\xi, \eta)=(UB_{x}+U_{\perp}C_{x}, VD_{x}+V_{\perp}E_{x})$$E$$D_{x}\in Skew(p)$ ,と書ける．ここで，
配 grad F $(U, V)$ もまた $(U, V)$ における接ベクトルであるので，$\in \mathbb{R}^{(n-p)\cross p}$$C_{x}\in \mathbb{R}^{(m-p)\cross p},$$B_{x},$$\mathbb{R}^{(m-p)\cross p},$お。

$E_{g}\in \mathbb{R}^{(n-p)\cross p}$$B_{g},$である．さらに，勾$D_{9}\in$

Skew(p),$C_{g}\in$を用いてgrad F $(U, V)=(UB_{g}+U_{1}C_{g}, VD_{g}+V_{\perp}E_{g})$と書ける．したがって，ニュートン方程式 (2.7) を解くことは，ヘシアンの作用によって$(-B_{g}, -C_{g}, -D_{g}, -E_{9})$に移されるような $(B_{x}, C_{x}, D_{x}, E_{x})$ を求めることと等価である．そのためにはヘシアンの作用を詳しく解析する必要があるが，ヘシアンは線形性を持つため適当な基底の下での表現行列がわかれば良いので， $(B_{x}, C_{x}, D_{x}, E_{x})$ の成分を 1 列に並べは反対称行列であるため，狭義下てベクトル表示することを考える．ただし，B。と三角部分の成分のみが独立であると考えることにする．そこで，行列$D_{x}$$X=(x_{i,j})\in \mathbb{R}^{m\cross n}$

に対して作用素 vec $()$ をvec(X) $=[x_{1,1}x_{2,1}\cdots x_{m,1}|x_{1,2}x_{2,2}\cdots x_{m,2}|\cdots|x_{1,n}x_{2,n}\cdots x_{m,n}]^{T}\in \mathbb{R}^{mn}$により定義する．同様に，反対称行列 $Z=(z_{i,j})\in Skew(p)$ に対して作用素 veck $()$ を

veck $(Z)=[z_{2,1}z_{3,1}\cdots z_{p,1}|z_{3,2}z_{4,2}\cdots z_{p_{)}2}|\cdots|z_{p,(p-1)}]^{T}\in \mathbb{R}^{p(p-1)/2}$

132と定める．すなわち，veck (Z) は $Z$ の狭義下三角部分の成分のみを左の列から順に縦に並べて得られる $p(p-1)/2$ 次元ベクトルである．これらの作用素を用いて$x:=\{\begin{array}{l}veck(B_{x})vec(C_{x})veck(D_{x})vec(E_{x})\end{array}\}, g:=-\{\begin{array}{l}veck(B_{g})vec(C_{g})veck(D_{g})vec(E_{g})\end{array}\}$とすることで，ニュートン方程式 (2.7) を$Hx=g$(3.3)と書くことができる．ここで， $H$ は Hess F $(U, V)$ の表現行列である．なお，線形方程式(3.3) の次元 $K:=p(m+n-p-1)$ は接空間 $T_{(U,V)}(St(p, m)\cross St(p, n))$ の次元と一致してといることに注意する．すなわち，(3.3) は吸 $u,v$ ) $(St(p, m)\cross St(p, n)$ をベクトル空間$\mathbb{R}^{K}$同一視したときの標準的な線形方程式である．

3.2

目的関数のヘシアンの表現行列とその性質

方程式 (3.3) の解法を探るために，ヘシアンの表現行列 $H$ を詳しく調べる必要がある．
によりそのためには，まず $B_{H},$ $D_{H}\in Skew(p)$ ,$C_{H}\in \mathbb{R}^{(m-p)\cross p},$ $E_{H}\in \mathbb{R}^{(n-p)\cross p}$

Hess F $(U, V)[(\xi, \eta)]=(UB_{H}+U_{\perp}C_{H}, VD_{H}+V_{\perp}E_{H})$(3.4)とし， $(B_{x}, C_{x}, D_{x}, E_{x})$ が $(B_{H}, C_{H}, D_{H}, E_{H})$ にどのように変換されるかを調べれば良い．

式 (2.6) から，$\{\begin{array}{l}UB_{H}+U_{\perp}C_{H}=\xi S_{1}-A\eta N-U sym (U^{T}(\xi S_{1}-A\eta N)) ,VD_{H}+V_{\perp}E_{H}=\eta S_{2}-A^{T}\xi N-V sym (V^{T}(\eta S_{2}-A^{T}\xi N))\end{array}$(3.5)

を得る．さて，$W:=U^{T}AV, W_{1}:=U^{T}AV_{1}, W_{\perp}.

とおく．式 (3.2) より $U^{T}\xi=B_{x},$に左からを乗じると，:=U_{\perp}^{T}AV, W_{\perp\perp}:=U_{\perp}^{T}AV_{\perp}$$U_{\perp}^{T}\xi=C_{x}$であることに注意して，(3.5) の第 1 式の両辺$U^{T}$$B_{H}=skew(B_{x}S_{1}-U^{T}A(VD_{x}+V_{\perp}E_{x})N)=skew(B_{x}S_{1}-WD_{x}N-W_{\perp}E_{x}N)$(3.6)となる．ここで，skew $()$ は行列の反対称部分を返す写像，すなわち skew $(Z)=(Z-Z^{T})/2$を乗じると，
である．同様に (3.5) の第 1 式の両辺に左から$U_{\perp}^{T}$$C_{H}=C_{x}S_{1}-U_{\perp}^{T}A(VD_{x}+V_{\perp}E_{x})N=C_{x}S_{1}-W_{\perp}.D_{x}N-W\perp\perp E_{x}N$

133
を得る．(3.5) の第 2 式についても同様に両辺に左から$V^{T}$

および$V_{\perp}^{T}$を乗じることで，$D_{H}=skew(D_{x}S_{2}-V^{T}A^{T}(UB_{x}+U_{\perp}C_{x})N)\cdot=skew(D_{x}S_{2}-W^{T}B_{x}N-W_{\perp}^{T}.C_{x}N)$

,

(3.7)

$E_{H}=E_{x}S_{2}-V_{\perp}^{T}A^{T}(UB_{x}+U_{\perp}C_{x})N=E_{x}S_{2}-W_{\perp}^{T}B_{x}N-W_{\perp\perp}^{T}C_{x}N$となる．後は，これらの式の両辺に vec や veck を作用させれば表現行列 $H$ の作用の様子がわかる．

ここで，クロネッカー積
に対して$\otimes$を導入する．すなわち，行列$X=(x_{ij})\in \mathbb{R}^{m\cross n},$ $Y\in \mathbb{R}^{p\cross q}$$X\otimes Y=\{\begin{array}{lll}x_{11}Y \cdots x_{1n}Y\vdots \vdots x_{m1}Y \cdots x_{mn}Y\end{array}\}$である．このとき，$X\in \mathbb{R}^{m\cross p},$$Y\in \mathbb{R}^{p\cross q},$$Z\in \mathbb{R}^{q\cross n}$に対して，$vec(XYZ)=(Z^{T}\otimes X)vec(Y)$であることが知られている [8]. さらに， $(i,j)$ 成分のみ 1 で，その他の成分はすべてあるような $m\cross n$ 行列を
と表し，

$0$

で

$E_{i,j}^{(m\cross n)}$$D_{p}= \sum_{p\geq i\geq j\geq 1}(E_{p(j-1)+i,j(p-(j+1)/2)-p+i}^{(p^{2}\cross p(p-1)/2)}-E_{p(i-1)+j,j(p-(j+1)/2)-p+i}^{(p^{2}\cross p(p-1)/2)})$とおくと，反対称行列 $S\in Skew(p)$ に対して$vec(S)=D_{p}veck(S)$ ,veck $(S)= \frac{1}{2}D_{p}^{T}vec(S)$が成り立つことが容易に証明される．
これらを用いると，式 $(3.6)-(3.7)$ により

veck$(B_{H})= \frac{1}{2}D_{p}^{T}$($(S_{1}\otimes I_{p})D_{p}$veck $(B_{x})-(N\otimes W)D_{p}$ veck $(D_{x})-(N\otimes W_{\perp})$ vec

vec $(C_{H})=(S_{1}\otimes I_{m-p})$ vec $(C_{x})-(N\otimes W_{\perp}.)D_{p}$ veck
veck $(D_{H})= \frac{1}{2}D_{p}^{T}(-(N\otimes W^{T})D_{p}$ veck$(D_{x})-(N\otimes W_{\perp\perp})$$(B_{x})-(N\otimes W_{\perp}^{T}.)$$-(S_{2}\otimes I_{p})D_{p}veck(D_{x}))$

vec

vec

$(E_{x})$$(E_{x})$,$(C_{x})$,$vec(E_{H})=-(N\otimes W_{\perp}^{T})D_{p}veck(B_{x})-(N\otimes W_{\perp\perp}^{T})vec(C_{x})+(S_{2}\otimes I_{n-p})vec(E_{x})$と計算できる．したがって$\mathbb{D}:=\{D_{p} I_{p(m-p)} D_{p} I_{p(n-p)}\}$

),

134
とすると，表現行列 $H$ は$H=\mathbb{D}^{T}\{\begin{array}{llll}\frac{1}{2}S_{1}\otimes I_{p} 0_{p^{2}\cross p(m-p)} -\frac{1}{2}N\otimes W -\frac{1}{2}N\otimes W_{1}0_{p(m-p)\cross p^{2}} S_{1}\otimes I_{m-p} -N\otimes W_{\perp} -N\otimes W_{\perp\perp}-\frac{1}{2}N\otimes W^{T} -\frac{1}{2}N\otimes W_{\perp}^{T} \frac{1}{2}S_{2}\otimes I_{p} 0_{p^{2}\cross p(n-p)}-N\otimes W_{\bullet\perp}^{T} -N\otimes W_{\perp\perp}^{T} 0_{p(n-p)\cross p^{2}} S_{2}\otimes I_{n-p}\end{array}\}\mathbb{D}$

(3.8)

と書ける．

一方，方程式 (3.3) の右辺も簡単に書き下せて，$g:=\{\begin{array}{l}veck(B_{g})vec(C_{g})veck(D_{9})vec(E_{g})\end{array}\}=\{\begin{array}{l}veck(skew(WN))vec(W_{\perp}.N)veck(skew(W^{T}N))vec(W_{\perp}^{T}N)\end{array}\}$となる．これで，(3.3) の係数行列 (ヘシアンの表現行列) も右辺ベクトル (勾配の列ベ
クトル表示) も明らかになったので， $H$ が正則であれば少なくとも理論的には$x=H^{-1}g$と解くことができる．これで veck
, veck $(D_{x})$ , vec $(E_{x})$ が求まるので，そこか, vec
らを構成でき，元のニュートン方程式 (2.7) の解を$(B_{x})$$(C_{x})$$B_{x},$ $C_{x},$ $D_{x},$ $E_{x}$$(\xi, \eta)=(UB_{x}+U_{1}C_{x}, VD_{x}+V_{\perp}E_{x})$と計算できる．ただし，実際に $H$ を陽に構成せずとも，任意のベクトル に対して $Hx$
さえ計算できれば線形方程式 (3.3) を解くことができ，次節ではその考えに基づき共役残差法による求解について述べる．
$x$

ところで，共役残差法を線形方程式に適用する際の前提として，その方程式の係数行列
が対称行列であるという条件が必要である．実は，(3.3) の係数行列 $H$ , すなわち (3.8) は対称ではない．しかし，容易にわかるように， $H$ の左から

$D’:=\{2I_{p(p-1)/2} I_{p(m-p)} 2I_{p(p-1)/2} I_{p(n-p)}\}$を乗じると， $D’H$ は対称行列である．このことは，接空間 $T_{(U,V)}(St(p, m)\cross St(p, n))$ におけるヘシアン Hess F $(U, V)$ の線形作用素としての対称性からも次のように導かれる．まを点 $(U, V)$ における接ベクトルとし，ず，,$(\xi_{1}, \eta_{1})$$(\xi_{2}, \eta_{2})$$(\xi_{1}, \eta_{1})=(UB_{1}+U_{\perp}C_{1}, VD_{1}+V_{\perp}E_{1})$,$(\xi_{2}, \eta_{2})=(UB_{2}+U_{\perp}C_{2}, VD_{2}+V_{\perp}E_{2})$135$D_{2}\in Skew(p)$ ,と書かれているとする．ここで，である．リーマン計量 (2.3) に基づいてこの 2 つの接ベクトルの内積を計算すると，$C_{1},$ $C_{2}\in \mathbb{R}^{(m-p)\cross p},$ $E_{1},$$B_{1},$ $B_{2},$ $D_{1},$$E_{2}\in$$\mathbb{R}^{(n-p)\cross p}$$\langle(\xi_{1}, \eta_{1}) , (\xi_{2}, \eta_{2})\rangle_{(U,V)}=tr(\xi_{1}^{T}\xi_{2})+tr(\eta_{1}^{T}\eta_{2})$$=tr(B_{1}^{T}B_{2}+C_{1}^{T}C_{2})+tr(D_{1}^{T}D_{2}+E_{1}^{T}E_{2})$$=2$. veck $(B_{1})^{T}veck(B_{2})+vec(C_{1})^{T}vec(C_{2})$$+2$ .

veck $(D_{1})^{T}veck(D_{2})+vec(E_{1})^{T}vec(E_{2})$$=\{\begin{array}{l}2\cdot veck(B_{1})vec(C_{1})2\cdot veck(D_{1})vec(E_{1})\end{array}\}\{\begin{array}{l}veck(B_{2})vec(C_{2})veck(D_{2})vec(E_{2})\end{array}\}$$=\{\begin{array}{l}veck(B_{1})vec(C_{1})veck(D_{1})vec(E_{1})\end{array}\}D’\{\begin{array}{l}veck(B_{2})vec(C_{2})veck(D_{2})vec(E_{2})\end{array}\}$となる．一方，ヘシアン Hess F $(U, V)$ はリーマン計量から定まる内積につぃて対称な線形作用素として定義されるので，$\langle(\xi_{1}, \eta_{1})$, Hess F$(U, V)[(\xi_{2}, \eta_{2})]\rangle_{(U,V)}=\langle$

Hess F$(U, V)[(\xi_{1}, \eta_{1} (\xi_{2}, \eta_{2})\rangle_{(U,V)}$が成り立つ．したがって，$\{\begin{array}{l}veck(B_{1})vec(C_{1})veck(D_{1})vec(E_{1})\end{array}\}D’H\{\begin{array}{l}veck(B_{2})vec(C_{2})veck(D_{2})vec(E_{2})\end{array}\}=[H\{\begin{array}{l}veck(B_{1})vec(C_{1})veck(D_{1})vec(E_{1})\end{array}\}]^{T}D’\{\begin{array}{l}veck(B_{2})vec(C_{2})veck(D_{2})vec(E_{2})\end{array}\}$となる．これは任意の$(\xi_{1\eta_{1}})$, $(\xi_{2}, \eta_{2})\in T_{(u,v)}(St(p, m)\cross St(p, n))$ について成り立つから，$D’H=H^{T}D’$(3.9)

を得る．式 (3.9) は行列 $D’H$ が対称行列であることを示している．したがって，線形方程式 (3.3) の両辺に左から $D’$ を乗じることで，係数行列が対称行列であるような方程式$D’Hx=D’g$(3.10)に書き直すことができる．

3.3

共役残差法による線形方程式の求解$H’:=D’H,$ $g’:=D’g$ とおいて，線形方程式 (3.3) および (3.10) と等価かつ，係数行列が対称行列である$H’$$H’x=g’$(3.11)

136が得られた．この線形方程式に対して共役残差法を適用したい．まず，共役残差法のアルゴリズムを記述しておく．
アルゴリズム 2 線形方程式 (3.11) に対する共役残差法$q_{0}=H’p_{0}$
とする．
を選び
–77J
, $r_{0}=g’-H’x_{0},$1: 初期値$\backslash$$\{\llcorner g$$p_{0}=r_{0},$$x_{0}$

2: for $l=0$ , 1, 2, . . . do

3:
4:
5:

$\alpha_{l}=(r_{l}, H’r_{l})/(q_{l}, q_{l})$.

$x_{l+1}=x_{l}+\alpha_{l}p_{l}.$$r_{l+1}=r_{l}-\alpha_{l}q_{l}.$$\beta_{l}=(r_{l+1}, H’r_{l+1})/(r_{l}, H’r_{l})$ .
6:$p_{l+1}=r_{l+1}+\beta_{l}p_{l}.$7:$q_{l+1}=H’r_{l+1}+\beta_{l}q_{l}.$8:9: end for

共役残差法においては H’ は行列-ベクトル積の形でしか現れないことに注意する．したがって，行列 H’ そのものを (3.8) に基づいて計算する必要はなく，任意のベクトル に対して $H’x$ が計算できさえすれば良い．
, Hess F $(U, V)[(\xi, \eta)]\in T_{(u,v)}(St(p, m)\cross St(p, n))$をそれぞれ式 (3.2), (3.4) のように表したとすると，式 $(3.6)-(3.7)$ により$x$$(\xi, \eta)$$B_{H},$ $C_{H},$ $D_{H},$ $E_{H}$は既に$B_{x},$ $C_{x},$ $D_{x},$ $E_{x}$で表せているため， $x$ から $Hx$ への変換，すなわち$\{\begin{array}{l}veck(B_{x})vec(C_{x})veck(D_{x})vec(E_{x})\end{array}\}\mapsto(B_{x}, C_{x}, D_{x}, E_{x})\mapsto(B_{H}, C_{H}, D_{H}, E_{H})\mapsto\{\begin{array}{l}veck(B_{H})vec(C_{H})veck(D_{H})vec(E_{H})\end{array}\}$は計算できる．ところで，変換 $y\mapsto D’y$ は， の一部の成分を 2 倍するだけであるから容易に計算できるため，結局 $H’:=D’H$ による行列-ベクトル積の計算は次のアルゴリズ
ムのように書ける．
$y$

アルゴリズム 3 $H’$ による行列-ベクトル積の計算
1: Given:$x\in \mathbb{R}^{K}.$2: $x=[veck(B)^{T}$$\mathbb{R}^{(m-p)\cross p},$3:

$B_{H},$$C_{H},$$vec(C)^{T}$$E\in \mathbb{R}^{(n-p)\cross p}$$D_{H},$$E_{H}$$veck(D)^{T}$$vec(E)^{T}]^{T}$となるように $B,$ $D\in Skew(p)$ , $C\in$を構成する．

を$\{\begin{array}{l}B_{H}=skew(B_{x}S_{1}-U^{T}A(VD_{x}+V_{\perp}E_{x})N) ,C_{H}=C_{x}S_{1}-U_{\perp}^{T}A(VD_{x}+V_{\perp}E_{x})N,D_{H}=skew(D_{x}S_{2}-V^{T}A^{T}(UB_{x}+U_{\perp}C_{x})N) ,E_{H}=E_{x}S_{2}-V_{\perp}^{T}A^{T}(UB_{x}+U_{\perp}C_{x})N\end{array}$により計算する．
4: $H’x=[2\cdot veck(B_{H})^{T}$$vec(C_{H})^{T}$$2\cdot veck(D_{H})^{T}$$vec(E_{H})^{T}]^{T}$137
これで，共役残差法を適用する準備が整った．なお，共役残差法では通常，相対残差ノルムがあらかじめ定めた定数 $\epsilon(>0)$ 未満になったところで線形方程式の十分な近似解が得られたと判断し，反復を終了する．すなわち，今回の場合であれば，(3.12)$\frac{\Vert g’-H’x_{l}\Vert}{\Vert g’\Vert}<\epsilon$

を $H’x=g’$ に対する共役残差法の停止条件とする．ただし，初期値は $x_{0}=0$ とする．以上の事柄の数値実験による検証は 5 節で行う．また，共役残差法の停止条件として，(3.12)
より良い条件の提案へと繋がる議論を次節で行う．

ニュートン方程式の解の精度とニュートン法の収束率4

一般のリーマン多様体 $M$ 上のニュートン法の局所的収束性について，次の定理が成り立つ [1, 7].

定理 4.1. リーマン多様体 $M$ 上の滑らかな目的関数 $f$ を最小化する最適化問題において，が存在するような $x_{*}\in M$ があると仮定する．このとgrad f $(x_{*})=0$ かつ Hess f
き，x、の $M$ における近傍 が存在して，任意のを初期点としてこの問題に対す$(x_{*})^{-1}$$\mathcal{U}$$x_{0}\in \mathcal{U}$るニュートン法により生成される点列$\{x_{k}\}$は$x_{*}$に 2 次収束する．

前節で述べたニュートン法の解法は，ニュートン方程式 (2.7) の代わりに標準形の線形方程式 (3.11) を解くというものであるが，数値誤差がなく厳密に解けるとすれば，理論的には (2.7) を解くことと等価である．したがって，定理 4.1 より，提案手法は局所的に 2 次収束する点列を生成する．

一方，ニュートン法の各反復におけるニュートン方程式を非厳密に解いた場合の収束性も議論されている．特に，リーマン多様体 $M$ 上では次の定理が成り立つ [1].
定理 4.2. 定理 4.1 で述べたリーマン多様体 $M$ 上のニュートン法の各反復において，ニュートン方程式

Hess f $(x_{k})[\eta_{k}]=-$ grad f$(x_{k})$(4.1)$(x_{k})+\epsilon_{k}$

(4.2)

の代わりに

Hess f $(x_{k})[\eta_{k}]=-$ grad fを満たすように探索方向であを計算したとする．ここで， $x_{k}\in M,$が存在するような $x_{*}\in M$ があると仮定する．
る．また，grad f $(x_{*})=0$ かつ Hess f
をなる $M$ の局所座標系とし，この局所座標系を用いた表現を により表すとする．さらに，定数 $\kappa>0$ が存在して$\eta_{k}$$\eta_{k},$$\epsilon_{k}\in T_{x_{k}}M$$(x_{*})^{-1}$$\wedge$$x_{*}\in \mathcal{U}’$$(\mathcal{U}’, \phi)$$\Vert\hat{\epsilon}_{k}\Vert\leq\min\{\Vert\overline{gradf}(\hat{x}_{k})\Vert^{\theta}, \kappa\}\Vert\overline{gradf}(\hat{x}_{k})\Vert$が任意ので成り立つとする．ここで， $\theta>0$ である．このとき， の $M$ における近傍が存在して，任意のを初期点としてこの非厳密なニュートン法により生成される点列は x、に超 1 次収束 $( \min\{1+\theta, 2\} 次収束 )$ する．特に， $\theta=1$ のときは 2 次$k$$\mathcal{U}$$x_{*}$$x_{0}\in \mathcal{U}$$\{x_{k}\}$

収束する．

138
$\epsilon_{k}$を残差と見なすと，方程式 (4.2) を満たす$\eta_{k}$はニュートン方程式 (4.1) の近似解であると考えられる．

定理 4.2 は本稿で考えている問題 (1.2) に適用できる．しかし，本稿では St $(p, m)\cross St(p, n)$の局所座標系を考えることなく，式 (3.1) に基づいて接空間における基底を適当に取ることで，その基底についての接ベクトルの成分表示やヘシアンの表現行列を用いてニュートン方程式を (3.11) に変換し，これを共役残差法で近似的に解くという方法を提案した．したがって，その差異についての議論は必要ではあるが，基本的には定理 4.2 の証明と同様にして，(3.11) を近似的に解いた場合のニュートン法の収束率を議論できる．すなわち，
各反復

$k$

における $H_{k}’x=g_{k}’$ の近似解

$x_{k}$

が(4.3)$H_{k}’x_{k}=g_{k}’+\epsilon_{k}$を満たすとし，定数$\kappa,$$\theta>0$

に対して

(4.4)

$\Vert\epsilon_{k}\Vert\leq\min\{\Vert g_{k}’\Vert^{\theta}, \kappa\}\Vert g_{k}’\Vert$が成り立つとき，アルゴリズム 1 においてニュートン方程式の解 ( k, 科の代わりに，(4.3)$\xi$$\eta$から構成した $(\xi_{k,\eta_{k}})\in T_{(u_{k},v_{k})}(St(p, m)\cross St(p, n))$ を用いると，局所的に$\min\{1+\theta, 2\}$ 次収束する点列が生成される．

を満たす$x_{k}$

特に， $\theta=1$ のとき，式 (4.4) は$\frac{||\epsilon_{k}||}{||g_{k}’||}\leq\min\{\Vert g_{k}’\Vert, \kappa\}$のときは (3.12) と同様に相対残差ノルムが定数 未満になったら共役が小さくなりが成り立残差法を終了するが，ニュートン法の反復が進んで未満になるまで共役残差法の反復をすることに
つようになると，相対残差ノルムがとなり，$\Vert g_{k}’\Vert>\kappa$$\kappa$$\Vert g_{k}’\Vert<\kappa$$\Vert g_{k}’\Vert$$\Vert g_{k}’\Vert$

が大きいうちは線形方程式 (3.11) を粗く解いても良いことになる．
なる．すなわち，
このように (3.11) を厳密に解かなくてもニュートン法の局所的な 2 次収束性は保たれると$\Vert g_{k}’\Vert$

いうわけである．

5

数値実験

本節では，3 節で提案した手法を用いた数値実験例を紹介する．実験は，Intel Core i74790 CPU 3.60 GHz, 16.0 GB RAM を搭載した PC により，MATLAB $R2014b$ を用いて行った．

まず，列直交行列, . . . , 1) としてを乱数を用いて生成し， $\Sigma_{*}=diag(n,$ $n-$とした．また，目的関数 $F$ を定義する際に必要な定数の対のそれぞれ左から 列
角行列 $N$ は， $N=diag(p,p-1, \ldots, 1)$ とした．以上から，
行列に対して， $(U_{opt}, V_{opt})\in$
目までを取り出してできる $m\cross p$ 行列 Uopt および$St(p, m)\cross St(p, n)$ は問題 (1.2) の最適解である．さらに，標準正規分布に従う乱数を各成および $n\cross p$ 行列」監 p を用いて，ニュートン法の初期点を分に持つ $m\cross p$ 行列$1$$U_{*}\in \mathbb{R}^{m\cross n},$ $V_{*}\in \mathbb{R}^{n\cross n}$$A=U_{*}\Sigma_{*}V_{*}^{T}$$U_{*},$$n\cross p$$V_{*}$$V_{opt}$$R_{m\cross p}$$\cross$$(U_{0}, V_{0})=(qf(U_{opt}+0.001R_{m\cross p}), qf(V_{opt}+0.001R_{n\cross p}))$

$p$

139
により構成した．すなわち，最適解に摂動を加えた点を初期点として実験を行った．これは，ニュートン法が大域的収束性を持たず，任意の初期点に対しては収束性が保証できな
いためである．なお，最適解の情報が利用できない場合に最適解の近くの点を得るためには，たとえば共役勾配法などの大域的収束性を持つ最適化手法を適用すれば良く，実用的には最適化の途中で共役勾配法からニュートン法に切り替えるというハイブリッドな方法が考えられる [7]. ニュートン法の反復回数は 3 回とし，各反復におけるニュートン方程式とした (3.12) を用い，最大反復回数は 500を解くための共役残差法の停止条件は回とした．問題のサイズは $m=3000,p=5$ と固定し， の値は 100, 500, 1000, 1500, 2000
としてそれぞれの場合で実験を行った．各問題に対するニュートン法についての勾配の相対誤差ノルム grad Fの収束履歴を表 1 に示す．表
grad F$\epsilon=10^{-12}$$n$$(U_{k}, V_{k})\Vert_{(U_{k},V_{k})}/\Vert$$\Vert$

1 より，

$n$

$(U_{0}, V_{0})\Vert_{(U_{0},V_{0})}$の値によらず勾配が滑らかに$0$に収束している様子がわかる．

表 1: $n=100$ , 500, 1000, 1500, 2000 での勾配の相対誤差ノルムの収束履歴

また， $n=500$ としたときのニュートン法の各反復におけるニュートン方程式に対する共役残差法の収束履歴を図 5.1 に，の相対誤差ノルム, 目的関数値の相対誤差 $|F(U_{*}, V_{*})-F(U_{k}, V_{k})|/|F(U_{*},$
勾配の相対誤差
ノルム grad Fの収束履歴を図 5.2-5.4 に示す．図grad F$k=1$ , 2, 3 は，ニュートン法の第 反復で現れるニュートン方程式に
5.1 において，対する共役残差法の，第 反復での残差ベクトルを表す．
$U_{k},$$\Vert U_{*}-U_{k}\Vert_{F}/\Vert U_{*}\Vert_{F},$$V_{k}$$V_{k}\Vert_{F}/\Vert V_{*}\Vert_{F}$$\Vert$$\Vert V_{*}-$$V_{*}$$(U_{k}, V_{k})\Vert_{(U_{k},V_{k})}/\Vert$$(U_{0}, V_{0})\Vert_{(U_{0},V_{0})}$$k$$r_{k,l},$$l$

これらの結果から，提案手法は最適解に 2 次収束する点列を生成できていることがわかる．特に，ニュートン方程式を解く際の共役残差法による相対残差ノルムも滑らかに収束している様子が確認できる．

6

結論と今後の課題
本稿では問題 (1.2) に対するニュートン法の各反復で現れるニュートン方程式を連立一次方程式に変換し，その共役残差法による効率的な解法を提案した．また，数値実験により提案手法の有効性を検証した．さらに，共役残差法による線形方程式の求解の精度とニュートン方程式の収束率の関係についても議論した．この議論を活かし，より効率的なアルゴリズムを提案することが今後の課題である．また，問題 (1.2) は特異値分解と関係があるが，[6] で示されているように，目的関数を変えると近似的同時特異値分解に対応
する St $(p, m)\cross St(p, n)$ 上の最適化問題が得られ，この問題に対しても本稿の議論は同様に成り立つと期待される．他にもシュティーフェル多様体上の最適化問題として定式化される応用問題は多くあり，本稿の結果を利用して，それらに対する解法アルゴリズムも提
案したい．

140

Number of iterations

図 5.1: 共役残差法による$\Vert r_{k,l}\Vert_{2}/\Vert g_{k}’\Vert_{2}$の収束履歴 $(k=1,2,3)$

Number of iterations

図 5.2:$\Vert U_{*}-U_{k}\Vert_{F}/\Vert U_{*}\Vert_{F}$および$\Vert V_{*}-V_{k}\Vert_{F}/\Vert V_{*}\Vert_{F}$

の収束履歴

141

Number of itel.ations

図 5.3: $|F(U_{*}, V_{*})-F(U_{k}, V_{k})|/|F(U_{*}, V_{*})|$ の収束履歴

Number of itcrations

図 5.4:$\Vert$grad F$(U_{k}, V_{k})\Vert_{(u_{k},v_{k})}/\Vert$grad F$(U_{0}, V_{0})\Vert_{(U_{0},V_{0})}$

の収束履歴

142
参考文献
[1] P.-A. Absil, R. Mahony, and R. Sepulchre, optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, 2008.
[2] A. Edelman, T.A. Arias, and S.T. Smith, The geometry of algorithms with orthogonality constraints, SIAM J. Matrix Anal. Appl., 20, 303-353, 1998.

[3] S.C. Eisenstat, H.C. Elman, and M.H. Schultz, Variational iterative methods for nonsymmetric systems of linear equations, SIAM J. Numer. Anal., 20, 345-357, 1983.
[4] G.H. Golub and C.F. Van Loan, Matrix Computations, 4th ed., The Johns Hopkins
University Press, Baltimore, 2013.

[5] J. Nocedal and S.J. Wright, Numerical optimization, 2nd ed., Springer, New York,2006.
[6] H. Sato, Joint singular value decomposition algorithm based on the Riemannian trustregion method, JSIAM Letters, 7, 13-16, 2015.
[7] H. Sato and T. Iwai, A Riemannian optimization approach to the matrix singularvalue decompositon, SIAM J. Optim., 23, 188-212, 2013.
[8] J.R. Schott, Matrix Analysis for Statistics, 2nd ed., Wiley, New York, 2005.