レイリー・リッツ法

レイリー・リッツ法(レイリー・リッツほう、: Rayleigh–Ritz method)は、固有値問題に対する数値近似解法の一つ。レイリー卿とヴァルター・リッツに名をちなむ。物理学上の境界値問題の解法として考案された。

固有値と固有ベクトルの近似をともなう全ての問題に応用でき、分野によってしばしば別名で呼ばれる。量子力学では、系を構成する粒子はハミルトニアンを用いて記述されるが、リッツ法(英語版)では試行波動関数を用いて最低エネルギー固有値に対応する固有関数を近似する。有限要素法の文脈では、数学的に等価なアルゴリズムが一般にリッツ・ガラーキン法(英語版)と呼ばれる。機械工学および構造工学では固有振動モードおよび共鳴周波数を近似する手法としてレイリー・リッツ法およびリッツ法という用語が用いられることが多い。

名称

本手法は1908年から1909年にかけてヴァルター・リッツが発表したもので、リッツ法と呼ぶべきであるという主張もある[1][2]。A. W. Leissa[1]によれば、レイリー卿は1911年にリッツの業績を顕彰する論文を書いたが、彼自身が書籍他の刊行物において本手法をすでに何度も用いていたと述べている。後に異論も出たもののこの主張にくわえ、射影に単一ベクトルを用いる自明な場合、本手法はレイリー商の計算に帰着するという事実もあり、異論もあるもののレイリー・リッツ法という名称が現在まで用いられている。S.Ilanko[2]リヒャルト・クーラントを引いて、レイリー卿とヴァルター・リッツがそれぞれ独立に、偏微分方程式境界値問題変分問題の等価性を活用し、有限のパラメータを決定すればよい極値問題で変分法を置き換えるというアイデアを独自に考案したとする。詳細については、リッツ法(英語版)の項を参照されたい。皮肉なことに、後にこの手法はより単純でより一般的な正射影を用いるよう改良され、ボリス・ガラーキン(英語版)に名を因んでガラーキン法(英語版)もしくはリッツ・ガラーキン法と呼ばれる。

行列の固有値問題への適用

数値線形代数において、レイリー・リッツ法は一般的[3]にサイズ N {\displaystyle N} 正方行列 A C N × N {\displaystyle A\in \mathbb {C} ^{N\times N}} についての固有値問題

A x = λ x {\displaystyle A\mathbf {x} =\lambda \mathbf {x} }
の近似解を得るために用いられる。まず、行列 A {\displaystyle A} はより小さなサイズの行列へと射影される。射影は、その列が正規直交系をなす射影行列 V C N × m {\displaystyle V\in \mathbb {C} ^{N\times m}} により行われる。行列版のレイリー・リッツ法は最も単純で、以下のように書き下せる。

  1. m × m {\displaystyle m\times m} 行列 V A V {\displaystyle V^{*}AV} を計算する。ここで、 V {\displaystyle V^{*}} V {\displaystyle V} 複素共役転置行列とする。
  2. 固有値問題 V A V y i = μ i y i {\displaystyle V^{*}AV\mathbf {y} _{i}=\mu _{i}\mathbf {y} _{i}} を解く。
  3. リッツベクトル x ~ i = V y i {\displaystyle {\tilde {\mathbf {x} }}_{i}=V\mathbf {y} _{i}} およびリッツ値 λ ~ i = μ i {\displaystyle {\tilde {\lambda }}_{i}=\mu _{i}} を計算する
  4. リッツ対と呼ばれる ( λ ~ i , x ~ i ) {\displaystyle ({\tilde {\lambda }}_{i},{\tilde {\mathbf {x} }}_{i})} をを元の行列 A {\displaystyle A} の固有値問題の近似解として出力する。

もし、行列 V C N × m {\displaystyle V\in \mathbb {C} ^{N\times m}} の列を正規直交基底として張られる線型部分空間が行列 A {\displaystyle A} の固有ベクトルに近い k m {\displaystyle k\leq m} 個のベクトルを含んでいれば、上記のレイリー・リッツ法はそれら固有ベクトルをよく近似する k {\displaystyle k} 個のリッツベクトルを与える。各リッツ対の精度は、容易に計算できる量 A x ~ i λ ~ i x ~ i {\displaystyle \|A{\tilde {\mathbf {x} }}_{i}-{\tilde {\lambda }}_{i}{\tilde {\mathbf {x} }}_{i}\|} により評価できる。

最も簡単な m = 1 {\displaystyle m=1} の場合、 N × m {\displaystyle N\times m} 行列 V {\displaystyle V} 単位列ベクトル v {\displaystyle v} 、行列 V A V {\displaystyle V^{*}AV} レイリー商 ρ ( v ) = v A v / v v {\displaystyle \rho (v)=v^{*}Av/v^{*}v} と一致するスカラーとなり、固有値問題の唯一の解は y 1 = 1 , μ 1 = ρ ( v ) {\displaystyle y_{1}=1,\mu _{1}=\rho (v)} 、唯一のリッツベクトルは v {\displaystyle v} それ自体となる。したがって、 m = 1 {\displaystyle m=1} の場合、レイリー・リッツ法はレイリー商の計算に帰着する。

別の有用なレイリー商とのつながりとして、各レイリー対 ( λ ~ i , x ~ i ) {\displaystyle ({\tilde {\lambda }}_{i},{\tilde {\mathbf {x} }}_{i})} に対して μ i = ρ ( v i ) {\displaystyle \mu _{i}=\rho (v_{i})} が成り立ち、したがってリッツ値は対応するレイリー商の理論から導かれるいくつかの性質をもつことが上げられる。たとえば、 A {\displaystyle A} エルミート行列のとき、そのレイリー商は(したがってリッツ値も)実数値をとり、 A {\displaystyle A} の最小固有値と最大固有値の間の閉区間におさまる。

行列

A = [ 2 0 0 0 2 1 0 1 2 ] {\displaystyle A={\begin{bmatrix}2&0&0\\0&2&1\\0&1&2\end{bmatrix}}}
の固有値は 1 , 2 , 3 {\displaystyle 1,2,3} であり、それぞれに対応する固有ベクトルは以下のとおりである。
x λ = 1 = [ 0 1 1 ] , x λ = 2 = [ 1 0 0 ] , x λ = 3 = [ 0 1 1 ] . {\displaystyle \mathbf {x} _{\lambda =1}={\begin{bmatrix}0\\1\\-1\end{bmatrix}},\quad \mathbf {x} _{\lambda =2}={\begin{bmatrix}1\\0\\0\end{bmatrix}},\quad \mathbf {x} _{\lambda =3}={\begin{bmatrix}0\\1\\1\end{bmatrix}}.}
ここで、
V = [ 0 0 1 0 0 1 ] {\displaystyle V={\begin{bmatrix}0&0\\1&0\\0&1\end{bmatrix}}}
とすると、
V A V = [ 2 1 1 2 ] {\displaystyle V^{*}AV={\begin{bmatrix}2&1\\1&2\end{bmatrix}}}
の固有値は 1 , 3 {\displaystyle 1,3} でありそれぞれに対応する固有ベクトルは以下のとおりとなる。
y μ = 1 = [ 1 1 ] , y μ = 3 = [ 1 1 ] {\displaystyle \mathbf {y} _{\mu =1}={\begin{bmatrix}1\\-1\end{bmatrix}},\quad \mathbf {y} _{\mu =3}={\begin{bmatrix}1\\1\end{bmatrix}}}
したがってリッツ値は 1 , 3 {\displaystyle 1,3} 、リッツベクトルは以下のように求まる。
x ~ λ ~ = 1 = [ 0 1 1 ] , x ~ λ ~ = 3 = [ 0 1 1 ] {\displaystyle \mathbf {\tilde {x}} _{{\tilde {\lambda }}=1}={\begin{bmatrix}0\\1\\-1\end{bmatrix}},\quad \mathbf {\tilde {x}} _{{\tilde {\lambda }}=3}={\begin{bmatrix}0\\1\\1\end{bmatrix}}}
この例で与えられた V {\displaystyle V} を用いると、リッツベクトルは A {\displaystyle A} の固有ベクトルのうち2つと完全に一致しており、リッツ値も3つの固有値のうち2つと完全に一致している。この例において近似解が厳密解と一致したのは、行列 V {\displaystyle V} 列空間が2つの固有ベクトル x λ = 1 {\displaystyle \mathbf {x} _{\lambda =1}} および x λ = 3 {\displaystyle \mathbf {x} _{\lambda =3}} によって張られる線形部分空間と一致していたからと説明される。

行列の特異値問題への適用

数値線形代数において、打ち切り特異値分解問題に対してもレイリー・リッツ法を適用することができる。これにより、サイズ M × N {\displaystyle M\times N} の行列 M C M × N {\displaystyle M\in \mathbb {C} ^{M\times N}} の与えられた線形部分空間内の左特異ベクトルおよび右特異ベクトルを近似的に求める問題を固有値問題に帰着させることができる。

正規行列を用いる場合

特異値 σ {\displaystyle \sigma } と、それに対応する左特異ベクトル u {\displaystyle u} と右特異ベクトル v {\displaystyle v} M v = σ u {\displaystyle Mv=\sigma u} および M u = σ v {\displaystyle M^{*}u=\sigma v} により定義される。左右どちらかの特異ベクトルおよび対応する特異値の集合を近似的に求めるには、エルミート正規行列 M M C N × N {\displaystyle M^{*}M\in \mathbb {C} ^{N\times N}} または M M C M × M {\displaystyle MM^{*}\in \mathbb {C} ^{M\times M}} のどちらか小さい方に対してナイーブにレイリー・リッツ法を適用すればよい。もう片方の特異ベクトルは単純に u = M v / σ {\displaystyle u=Mv/\sigma } または v = M u / σ {\displaystyle v=M^{*}u/\sigma } のように行列をかけて特異値で割れば得られる。しかし、割り算は特異値がゼロもしくはそれに近いとき不安定となる。

別のアプローチとして、サイズ N × N {\displaystyle N\times N} の正規行列 A = M M C N × N {\displaystyle A=M^{*}M\in \mathbb {C} ^{N\times N}} に対する、 N × m {\displaystyle N\times m} 正規直交行列 W C N × m {\displaystyle W\in \mathbb {C} ^{N\times m}} を用いたレイリー・リッツ法を適用すると、 m × m {\displaystyle m\times m} 行列

W A W = W M M W = ( M W ) M W {\displaystyle W^{*}AW=W^{*}M^{*}MW=(MW)^{*}MW}
に対する固有値問題を解くことになるが、これは N × m {\displaystyle N\times m} 行列 M W {\displaystyle MW} についての特異値問題とみなせることを活用する方法がある。この見方により左右両方の特異ベクトルを同時に得ることが以下のようにして可能となる。

  1. N × m {\displaystyle N\times m} 行列 M W {\displaystyle MW} を計算する。
  2. Thin-SVD[訳語疑問点] M W = U Σ V h {\displaystyle MW=\mathbf {U} \Sigma \mathbf {V} _{h}} を解き、 N × m {\displaystyle N\times m} 行列 U {\displaystyle \mathbf {U} } m × m {\displaystyle m\times m} 対角行列 Σ {\displaystyle \Sigma } m × m {\displaystyle m\times m} 行列 V h {\displaystyle \mathbf {V} _{h}} を得る。
  3. リッツ左特異ベクトル U = U {\displaystyle U=\mathbf {U} } とリッツ右特異ベクトル V h = V h W {\displaystyle V_{h}=\mathbf {V} _{h}W^{*}} を計算する。
  4. リッツ特異トリプレットと呼ばれる三つ組 U , Σ , V h {\displaystyle U,\Sigma ,V_{h}} が、元の行列 M {\displaystyle M} の打ち切り特異値分解問題の行列 W {\displaystyle W} の列空間上における近似解である。

このアルゴリズムは、固有値問題ソルバ(例:LOBPCG(英語版))から出力された行列 W {\displaystyle W} に対する後処理として実行することができる。

行列

M = [ 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 0 0 0 0 ] {\displaystyle M={\begin{bmatrix}1&0&0&0\\0&2&0&0\\0&0&3&0\\0&0&0&4\\0&0&0&0\end{bmatrix}}}
の正規行列は
A = M M = [ 1 0 0 0 0 4 0 0 0 0 9 0 0 0 0 16 ] {\displaystyle A=M^{*}M={\begin{bmatrix}1&0&0&0\\0&4&0&0\\0&0&9&0\\0&0&0&16\\\end{bmatrix}}}
また特異値は 1 , 2 , 3 , 4 {\displaystyle 1,2,3,4} 、対応するThin-SVDは以下のように求まる。
A = [ 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 ] [ 4 0 0 0 0 3 0 0 0 0 2 0 0 0 0 1 ] [ 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 ] , {\displaystyle A={\begin{bmatrix}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\\0&0&0&0\end{bmatrix}}{\begin{bmatrix}4&0&0&0\\0&3&0&0\\0&0&2&0\\0&0&0&1\end{bmatrix}}{\begin{bmatrix}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\end{bmatrix}},}
ここで左側の行列の列は行列 A {\displaystyle A} の左特異ベクトルの完全集合をなし、真ん中の対角行列の対角要素は特異値、右側の行列の列は右特異ベクトルの転置となっている(ただし、転置しても元と変わらない)
[ 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 ] = [ 0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 ] {\displaystyle {\begin{bmatrix}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\end{bmatrix}}^{*}\quad =\quad {\begin{bmatrix}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\end{bmatrix}}}
ここで、特異値1、2に対応する2つの右特異ベクトル厳密解

[ 0 1 1 0 0 0 0 0 ] {\displaystyle {\begin{bmatrix}0&1\\1&0\\0&0\\0&0\end{bmatrix}}}

により張られる空間を列空間とする行列、

W = [ 1 / 2 1 / 2 1 / 2 1 / 2 0 0 0 0 ] {\displaystyle W={\begin{bmatrix}1/{\sqrt {2}}&1/{\sqrt {2}}\\1/{\sqrt {2}}&-1/{\sqrt {2}}\\0&0\\0&0\end{bmatrix}}}

を導入する。

上記アルゴリズムのステップ1に従い、次の行列を得る。

M W = [ 1 / 2 1 / 2 2 2 0 0 0 0 ] {\displaystyle MW={\begin{bmatrix}1/{\sqrt {2}}&1/{\sqrt {2}}\\{\sqrt {2}}&-{\sqrt {2}}\\0&0\\0&0\end{bmatrix}}}
さらに、ステップ2に従い、thin-SVD M W = U Σ V h {\displaystyle MW=\mathbf {U} {\Sigma }\mathbf {V} _{h}} を解くと以下を得る。
U = [ 0 1 1 0 0 0 0 0 0 0 ] , Σ = [ 2 0 0 1 ] , V h = [ 1 / 2 1 / 2 1 / 2 1 / 2 ] {\displaystyle \mathbf {U} ={\begin{bmatrix}0&1\\1&0\\0&0\\0&0\\0&0\end{bmatrix}},\quad \Sigma ={\begin{bmatrix}2&0\\0&1\end{bmatrix}},\quad \mathbf {V} _{h}={\begin{bmatrix}1/{\sqrt {2}}&-1/{\sqrt {2}}\\1/{\sqrt {2}}&1/{\sqrt {2}}\end{bmatrix}}}
したがって、 Σ {\displaystyle \Sigma } から特異値として2と1が得られ、 U {\displaystyle \mathbf {U} } から対応する左特異ベクトル u {\displaystyle u} として [ 0 , 1 , 0 , 0 , 0 ] {\displaystyle [0,1,0,0,0]^{*}} [ 1 , 0 , 0 , 0 , 0 ] {\displaystyle [1,0,0,0,0]^{*}} が得られた。これら二つのベクトルは W {\displaystyle W} の列空間を張っており、与えらえれた W {\displaystyle W} による近似解が厳密解と一致する理由を説明する。

最後に、ステップ3に従い V h = V h W {\displaystyle V_{h}=\mathbf {V} _{h}W^{*}} を計算すると以下を得る。

V h = [ 1 / 2 1 / 2 1 / 2 1 / 2 ] [ 1 / 2 1 / 2 0 0 1 / 2 1 / 2 0 0 ] = [ 0 1 0 0 1 0 0 0 ] {\displaystyle \mathbf {V} _{h}={\begin{bmatrix}1/{\sqrt {2}}&-1/{\sqrt {2}}\\1/{\sqrt {2}}&1/{\sqrt {2}}\end{bmatrix}}\,{\begin{bmatrix}1/{\sqrt {2}}&1/{\sqrt {2}}&0&0\\1/{\sqrt {2}}&-1/{\sqrt {2}}&0&0\end{bmatrix}}={\begin{bmatrix}0&1&0&0\\1&0&0&0\end{bmatrix}}}
したがって右特異ベクトル v {\displaystyle v} [ 0 , 1 , 0 , 0 ] {\displaystyle [0,1,0,0]^{*}} および [ 1 , 0 , 0 , 0 ] {\displaystyle [1,0,0,0]^{*}} である。このうち前者が M v = σ u {\displaystyle Mv=\sigma u} を満たすことは以下のように確かめられる。
[ 1 0 0 0 0 2 0 0 0 0 3 0 0 0 0 4 0 0 0 0 ] [ 0 1 0 0 ] = 2 [ 0 1 0 0 0 ] {\displaystyle {\begin{bmatrix}1&0&0&0\\0&2&0&0\\0&0&3&0\\0&0&0&4\\0&0&0&0\end{bmatrix}}\,{\begin{bmatrix}0\\1\\0\\0\end{bmatrix}}=\,2\,{\begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}}}
また、 M u = σ v {\displaystyle M^{*}u=\sigma v} を満たすことも以下のように確かめられる。
[ 1 0 0 0 0 0 2 0 0 0 0 0 3 0 0 0 0 0 4 0 ] [ 0 1 0 0 0 ] = 2 [ 0 1 0 0 ] {\displaystyle {\begin{bmatrix}1&0&0&0&0\\0&2&0&0&0\\0&0&3&0&0\\0&0&0&4&0\end{bmatrix}}\,{\begin{bmatrix}0\\1\\0\\0\\0\end{bmatrix}}=\,2\,{\begin{bmatrix}0\\1\\0\\0\end{bmatrix}}}
このように、行列 W {\displaystyle W} の列空間が右特異ベクトルの厳密解により張られる空間と一致するとき、それら右特異ベクトルと、対応する左特異ベクトルおよび特異値の厳密解が得られる。任意の行列 W {\displaystyle W} に対しては、レイリー・リッツ法の意味で最適な特異値分解が近似解として求まる。

関連項目

  • リッツ法(英語版)
  • レイリー商
  • アーノルディの反復法(英語版)

出典

  1. ^ a b Leissa, A.W. (2005). “The historical bases of the Rayleigh and Ritz methods”. Journal of Sound and Vibration 287 (4–5): 961–978. Bibcode: 2005JSV...287..961L. doi:10.1016/j.jsv.2004.12.021. https://www.sciencedirect.com/science/article/abs/pii/S0022460X05000362. 
  2. ^ a b Ilanko, Sinniah (2009). “Comments on the historical bases of the Rayleigh and Ritz methods”. Journal of Sound and Vibration 319 (1–2): 731–733. Bibcode: 2009JSV...319..731I. doi:10.1016/j.jsv.2008.06.001. 
  3. ^ Trefethen, Lloyd N.; Bau, III, David (1997). Numerical Linear Algebra. SIAM. p. 254. ISBN 978-0-89871-957-4. https://books.google.com/books?id=JaPtxOytY7kC 

外部リンク

  • Course on Calculus of Variations, has a section on Rayleigh–Ritz method.