Hermite対称な定値一般固有値問題のフィルタ対角化法について

全文

(1)情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. Hermite 対称な定値一般固有値問題のフィルタ対角化法について村上弘1,a). 概要：係数行列 A と B が Hermite 対称で B が正定値である一般固有値問題 Av=λBv では固有値 λ はすべて実数である．この種類の問題に対して，指定した実区間に固有値がある固有対の近似だけを求めるフィルタ対角化法を示す．. キーワード：フィルタ対角化法，固有値問題, Hermite 対称定値. Filter Diagonalization for Complex Hermitian Definite Generalized Eigenproblems Hiroshi Murakami1,a). Abstract: Eigenvalues λ are all reals for a generalized eigenproblem Av=λBv whose coefficient matrices A and B are complex Hermitian and B is positive definite. For this kind of problems, a filter diagonalization method is presented which solves only those approximations of eigenpairs whose eigenvalues are in the specified real interval. Keywords: filter diagonalization, eigenproblem, complex Hermitian definite. 1. はじめに. し，固有ベクトルは実対称の問題の場合には B を計量とする正規直交条件 V T BV =I を，Hermite 対称の問題の場合. 係数行列 A と B が Hermite 対称で B が正定値の定値一. には V H BV =I を満たすようにできる．ここで右肩に H. 般固有値問題（GEVP）Av=λBv は，係数行列が実対称の. をつけたものは通常の Hermite 転置を表す．そうして問. 場合の問題を Hermite 対称の場合に拡張したものになって. 題の線形空間に付随するベクトル同士の自然な内積 hu, vi. いる．実対称の場合と同じく，この問題の固有値 λ はすべ. は，実対称の場合は uT Bv であり，Hermite 対称の場合は. て実数である．そのことから，求めたい固有対の固有値の. uH Bv である．. 存在範囲を実数の区間で指定できる．但し，固有ベクトル. これらのことから，実対称の場合の解法に現れる転置操. は実対称の場合はすべて実ベクトルにとれたが，Hermite. 作の箇所を Hermite 転置に置き換えることで，Hermite 対. 対称の場合には一般には実ベクトルにはできない．. 称の場合の問題が同様に扱えると期待できる．そこで本研. どちらの場合にも固有ベクトルの全体は固有値に重複が. 究では，既に構築済みの実対称の定値一般固有値問題に対. あっても必ず全空間を張る．そこで固有ベクトルの全体を. するフィルタ対角化の方法の構造をなるべく保ちながら拡. 順番に並べた行列を V ，それに対応して固有値を並べた対. 張し，Hermite 対称の場合を扱える定式化を得ることを試. 角行列を Λ とすれば，一般固有値問題は AV =BV Λ に帰. みた．. 1. a). 首都大学東京・数理情報科学専攻 Department of Mathematics and Information Sciences, Tokyo Metropolitan University [email protected]. ⓒ2012 Information Processing Society of Japan. 2. フィルタとその伝達関数大規模な N 次の Hermite 対称な定値一般固有値問題. 1.

(2) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. （GEVP）：Av=λBv の「区間 [a, b] に固有値がある固有対」だけを，通過域 [a, b] のフィルタを利用してフィルタ対角化法により求めるとする．シフトが複素数 τ のレゾルベントを R(τ )=(A − τ B). γp は，伝達関数の形状に対して課す制約条件を満たすように決める．実対称の問題の場合における通過域 [a, b] のフィルタの. −1. B. と定義する．フィルタ作用素 F として，シフトが相異なる ∑ レゾルベントの線形結合 F=c∞ I + p γp R(τp ) を採用す. 設計と構成の方法（フィルタの特性に対して課された条件を満たすための伝達関数の極の個数と極の配置，極の留数の値を決める方法）は，Hermite 対称の問題の場合にもそ. る．そのとき，元の固有値問題の固有値 λ の固有ベクトル. のまま利用できる．. v に対しては Fv=f (λ)v が成立する．ここで λ の有理関 ∑ 数 f (λ)=c∞ + p=1 γp /(λ − τp ) は，フィルタによる固有. 2.1 フィルタの設計. 値が λ の固有ベクトルの伝達率を与えるので伝達関数と呼. まず λ の正規化座標 t を λ ∈ [a, b] から t ∈ [−1, 1] へ. ばれる．Hermite 対称な定値一般固有値問題の固有値は実. の線形変換で定義する．そのとき t ∈ [−1, 1] は通過域. 数であるので，フィルタの設計の際には伝達関数 f (λ) の. （passband）に，µ ≤ |t| は阻止域（stopbands）に，中間の. 振る舞いは λ が実軸の上だけで調べれば良い．このことか. 1 < |t| < µ は遷移域（transitionbands）に，それぞれ対応. ら，フィルタの伝達関数 f (λ) が実対称の場合と同一のも. する．但し，µ は 1 より（少し）大きい値を持つパラメタ. のを Hermite 対称の問題においても使用できる．固有値は. である．フィルタの種類とは伝達関数 f (λ) = g(t) の関数. 実数だけなので，すべてのシフトを虚数にとれば F の有界. 族の違いである．. 性は保証され，伝達関数 f (λ) も実軸上の有界な関数になる．さらに後でみるように，伝達関数 f (λ) を実関数にすると，シフトが互いに複素共役なレゾルベント間の対称性. に対して実数値を与える関数）に限定する． ∑ 伝達関数 f (λ)=c∞ + p γp /(λ − τp ) （但し τp. の値は全て相異なるとする）が実関数となるには，. gstop. 係数 c∞ は実数で，実数の極はその留数は実数で，. 0. 虚数の極は複素共役の対として含まれ，複素共役な極に対する留数もまた複素共役であることが必. 要十分である．つまり実数の極の場合も含めて，. STOP. PASS −µ. −1. t. STOP 1. STOP. そこで本論文では，伝達関数 f (λ) は実関数（実数の引数. gpass |g(t)|. を利用することで計算量が削減できて有利である（後述）．. PASS. 1. µ. 図 1 フィルタの伝達関数の絶対値 |g(t)| の概形図．. 任意の添字 p に対してある添字 q が唯一決まり，. τp =τq ，γp =γq となる．（伝達関数 f (λ) の有界性. 伝達関数 g(t) の値の分布に制約条件： i) 通過域における. を仮定する場合には，極は虚数のみであり，添字. g(t) の tight な最小値は gpass で，逆に g(t) が gpass 以上に. の個数は偶数 2n で，添字 1 から n までの極は例え. なるのは通過域に限る； ii) 阻止域における |g(t)| の上限値. ば虚部が正のものとし，それらの複素共役が添字. n + 1 から 2n までの極となるように並べられる．）伝達関数 f (λ) が実関数の場合は，BF は Hermite 対称になる．実際，c∞ は実数で，任意の極の添字 p にはそれに. は gstop である；を課す（図 1 参照）．条件 i) は g(t) ≥ gpass と通過域 [−1, 1] に t があることが同値であることを意味している．典型的な 4 種類のフィルタ Butterworth，Chebyshev，. 対応する添字 q が唯一存在して，τp =τq ，γp =γq となるこ. inverse-Chebyshev，elliptic に対しては，制約条件の 3 つ. とを用いると， ∑ (BF)H = B c∞ I + p { γp B (A − τp B)−1 B }H ∑ = B c∞ I + p γp B H (A − τp B)−H B H ∑ = B c∞ I + p γp B(A − τp B)−1 B ∑ = B c∞ I + q γq B(A − τq B)−1 B ∑ = B{ c∞ I + q γq (A − τq B)−1 B }. 組（µ，gpass ，gstop ）を指定すると，その条件を満たせる次. = BF として示せる．. 数 n の最小値が決まり，n をその値以上に設定すると g(t) が完全に決まる．（同様に制約条件の 3 つ組を（µ，gpass ，. n）で与える場合は，条件を満たせる gstop の最小値が決まり，gstop をその値以上に設定すれば g(t) が完全に決まる．あるいはまた制約条件の 3 つ組を（gpass ，gstop ，n）で与える場合には，条件を満たせる µ の最小値が決まり，µ をその値以上に設定すれば g(t) が完全に決まる．）関数 g(t) が決まると，有理関数である伝達関数 f (λ)=g(t). フィルタ F あるいは伝達関数 f (λ) を定義するパラメタ. のすべての極の位置と各極の持つ留数が決まり，それによ. である次数 2n，実係数 c∞ ，複素シフト τp ，複素結合係数. り通過域が λ ∈ [a, b] であるフィルタ F のパラメタが完全. ⓒ2012 Information Processing Society of Japan. 2.

(3) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. に決まる．フィルタの設計法と構成法の詳細は文献 [6] に記述した（但し，そこで用いた記号や定義は本論文のものと多少異なる）．. 3. レゾルベントの作用の計算法について. 3.1 シフトが複素共役なレゾルベントの間の LU 分解の共用. 実対称の問題の場合には，任意の実ベクトルの組 X に（伝達関数 f (λ) が実関数の）フィルタ F を作用させる計算の中で，複素共役なシフトの添字を p，q とするとき，複素共役性から関係 { γp R(τp ) + γq R(τq ) }X = 2Re {γp R(τp )X}. フィルタの作用を実現するには，N 次（列）ベクトルの. が成り立つので複素共役対のシフトを持つ 2 個のレゾルベ. m 個の組である X に対して，複素数のシフト τp のレゾル. ントの項の作用の和は，片方の項の作用だけの計算で済ま. ベント R(τp )=(A − τp B)−1 B を各 p について作用させる. せることで，計算量と記憶参照量が半分にできていた．. 計算が必要である．それには C=A − τp B とおいて，連立. Hermite 対称の問題では実対称の場合に比べて複素共役. 1 次方程式の組 CZ=BX をベクトルの組 Z について解く. 性の利用は多少複雑になる．いま A，B が Hermite 対称で. ことに帰着する（右辺 BX はシフト τp によらないので，. あるとして，任意の複素ベクトルの組 X に（伝達関数 f (λ). X 0 := BX を計算して作っておくのが便利である）．行列. が実関数の）フィルタ F を作用させた結果を求めることに. A，B が Hermite 対称のとき，τp の値が実数でなければ行. する．互いに複素共役なシフトの添字を p，q とするとき，. 列 C は Hermite 対称にはならないことに注意する．. τp =τq で γp =γq である．いま C(τp ) ≡ A − τp B と定義す. ここでは A，B が帯行列でその半帯幅を h とすると，C. れば，添字 p のレゾルベントは R(τp )=C(τp )−1 B である. も半帯幅が h となり，連立 1 次方程式は帯行列専用の行ピ. から，レゾルベントの作用 V := R(τp )X の計算は，まず. ボット交換付きの LU 分解法を用いて安定に解ける．その. X 0 := BX とおいて，連立 1 次方程式の組 C(τp )V =X 0 を. 場合の LU 分解の L の下半帯幅は h で，U の上半帯幅は. 解けばよい．これを LU 分解で解く場合について以下で考. 2h 以下となる．. 察する．. 固有値は実数しかないので，シフト τ が実数でなければ. いま C(τp ) の LU 分解が C(τp ) =: LU と構成できたと. 行列 C は必ず正則であることから，行ピボット交換を省い. する．すると分解された連立 1 次方程式の組 LU V =X 0 を. て計算を行う（その場合の LU 分解の L の下帯幅と U の. 前進消去と後退代入で解いて V が得られる．同様にシフト. 上帯幅は共に h となる）と計算量と記憶量を節約できる可. が複素共役である τq のレゾルベントの作用 W := R(τq )X. 能性があるが，数値的な精度が悪化する可能性がある．. の計算には連立 1 次方程式 C(τq )W =X 0 を解けばよいが，. 固有ベクトルを逆反復法で改良する計算を行う場合に. C(τq ) = A − τp B = C(τp )H なのでその係数行列の LU. は，近似固有値である λ は実数であるので，それに対応し. 分解は先ほど構成した C(τp ) の LU 分解の因子を用いて. て Hermite 対称行列 C=A − λB を係数とする連立 1 次方. C(τq )H =: U H LH と表せる（実際 U H は下三角行列で，LH. 程式を解くことになるが，これもむしろ C の Hermite 対称. は上三角行列である）．すると U H LH W =X 0 を前進消去と. 性を用いて修正コレスキ分解 C=LDLH を用いるよりも，. 後退代入で解いて W が得られる（あるいは U H LH W =X 0. 行ピボット交換付きの帯行列用の LU 分解を用いて解いた. のかわりに，それの複素共役 U T LT W =X 0 を前進消去と. 方が計算量と記憶量は多くなるが数値精度を保つ為には望. 交代代入で解いて W を得た後で，その複素共役として W. ましいであろう．. を求める方が，多少効率的であろう）．結局，フィルタ F. 係数の Hermite 対称行列 A，B が帯行列ではなくて一般. の複素ベクトル X への作用の計算の中で，添字 p，q の互. 的な疎行列の場合にも，複素数 τ のシフトを持つレゾルベ. いに複素共役なシフトを持つ 2 つのレゾルベントの項の和. ント R(τ ) をベクトルの組 X に作用させる計算は，同様に. の寄与は，{γp R(τp ) + γq R(τq )}X = γp V + γq W となる. C=A − τ B とおいて，疎行列 C を係数とする連立 1 次方. （注：実対称の場合とは異なり，Hermite 対称の場合には，. 程式の組 CZ=BX をベクトルの組 Z について解くことに. 一般的にはたとえ X が実ベクトルの組でも，フィルタを. 帰着できる．行列 C の疎性を利用できる連立 1 次方程式の. 作用した結果は実ベクトルの組にはならない）．. 組に対する標準的解法としては，たとえば複数個のベクト. 以上のことから，互いに複素共役であるシフトを持つ 2. ルの組に対する（ブロック）クリロフ部分空間法の系統の. 個のレゾルベントの作用の計算には，片方の LU 分解の結. 反復解法，あるいはそれに不完全 LU 分解などの前処理を. 果を共有してもう片方で再利用することで，LU 分解にか. 加えたものなどが考えられる．また行列 A，B が偏微分方. かる計算量と記憶参照を半分に減らせることがわかる．. 程式の離散化近似に由来するものであれば，問題の構造に. 付記：シフト τ が実数の場合には C(τ )=A − τ B も Her-. 密着したマルチグリッド法などの高速解法が適用できる可. mite 対称である．B は正定値だから，「もしも」シフト τ. 能性がある．. が最小固有値よりも小さい実数であれば，Hermite 対称行列用の修正コレスキ分解 C(τ ) =: LDLH の計算はピボット交換なしでも安定に行える．. ⓒ2012 Information Processing Society of Japan. 3.

(4) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. る．すなわち Z は B-Hermite 正規直交（Z H BZ = I ）であるから，小さい次数の Hermite 対称行列 f A = Z H AZ. 4. フィルタへの入力とその出力いま整数 m を，フィルタの通過域または遷移域に固有値がある固有対全部の個数よりもいくらか大きい数にとる．そうしてまず m 個以上の乱数ベクトルの組に対して. B-Hermite 正規直交化を行い，B-Hermite 正規直交系である m 個のベクトルの組 X を作成する（つまり X BX=Im H. を作り，そのユニタリ変換による対角化（固有値分解）をf A = U DU H とするとき，Ritz 対は ZU の列ベクトル（Ritz ベクトル）とそれに対応する D の対角要素（Ritz 値) の対として得られる．この方法による近似固有値（Ritz 値）は Hermite 対称行列 f A の固有値であるから実数である．. である）．この X をフィルタへの入力ベクトルの組として用いる．このようにして得られた組 X の各ベクトル x(k) ，. 5.1 不変部分空間 I の近似基底 Z の構成フィルタ作用素 F がレゾルベントの線形結合の場合には，. k=1, 2, . . ., m のそれぞれに対してフィルタ F を適用し，得. 元の GEVP：Av=λBv の固有対を (λ, v) とすると，ベクト. られた結果 y. ル v は同時にまた F の固有ベクトルにもなり，Fv=f (λ)v. (k). (k). =Fx. ，k=1, 2, . . ., m を対応順に並べて，. m 個のベクトルの組 Y （= FX ）を作る．この Y がフィルタからの出力ベクトルの組になる．数値丸め誤差と比べて阻止域でのフィルタの伝達率の大きさが無視できる程度に小さければ，各出力ベクトルは. も満たす．ただし，f (λ) は伝達関数である．逆の命題「F の任意の固有ベクトルは元の GEVP の固有ベクトルである」は成り立たないが，「F の（任意の）不変部分空間は元の GEVP の不変部分空間である」は成り立つ．. y(k) は固有値がフィルタの通過域と遷移域にある固有ベク. この性質をうまく利用すると，「[a, b] のある近傍」に固. トルの線形結合とみなせるから，Span(Y ) はそのような固. 有値がある元の GEVP の不変部分空間を I とするとき，. 有ベクトルで張られた部分空間である．計算効率の観点からは，X に F を作用させる操作の内側で行われる行列ベクトル積や連立 1 次方程式を解く際に，. 以下の性質を満たすベクトルの組 Z を構成できる：. ( 1 ) Z は Y の張る空間内にある． ( 2 ) Z の張る空間は不変部分空間 I をよく近似する．. 組 X に含まれる各列ベクトルを 1 つずつ独立に F を適用. ( 3 ) Z は B-Hermite 正規直交系である．. するよりも，ある程度の個数をまとめてブロック化して処. そのように構成された組 Z に対して Rayleigh-Ritz 法を適. 理するとトータルの記憶参照回数が減らせるなどの利点が. 用すると，元の GEVP の指定区間 [a, b] のある近傍に固有. ある．但し，作業用の記憶の量をそれだけ多く確保する必. 値があるすべての固有対の近似を一斉に求めることがで. 要がある．. きる．. 5. Rayleigh-Ritz 法に与える基底の組の構成. Z の構成法は，まず次節の方法（節 5.2）を用いて組 Y の張る空間内で作用素 F の近似固有対 (φ, v) で（伝達関数. フィルタの出力ベクトルの組 Y が張る空間は，「通過域. f の [a, b] の値域である）[gpass , 1] の「適切な近傍」に固有. [a, b] の近傍に固有値がある固有ベクトル」で張られてい. 値 φ の値があるものだけをすべて求める．そのようなベク. る．そこで，Y の張る空間の内で元の GEVP の近似固有対. トル v は B-Hermite 正規直交系にとれ，それらを集めて. を Rayleigh-Ritz 法を用いて求めようとする．しかしフィ. 並べたものを組 Z にすればよい．. ルタの阻止域での性質により，固有値が阻止域の近傍にあ. ここで [gpass , 1] の「適切な近傍」とは，固有値. る固有ベクトルは非常に強く減衰をうけるので，Y は通常. φ の分布を（誤差の影響を考慮した上で）はっき. は数値的階数が多数落ちている行列である（しかも入力ベ. りと分離するような近傍のことである．数値計算. クトルの個数 m を増せばそれだけ，より多くの階数が落. で求めた固有値 φ の値は近似の誤差や数値演算の. ちる）．そのような数値的条件の悪い基底 Y に対して直接. 丸め誤差，連立 1 次方程式の解計算の不完全性な. に Rayleigh-Ritz 法を適用すると，得られる結果もまた数. どにより誤差を含んでずれる．たとえば gpass と. 値的に非常に不安定になる．. 比べて小さい φ の値でも，それが gpass に非常に. そこで一種の正則化として，Y の張る空間の内から数値. 近い値であれば安全のためには含めておく必要が. 的線形独立性の高い基底 Z をあらかじめ選び出す処理を. ある．また，固有ベクトルの組が不変部分空間に. 行う．いま固有値が [a, b]（のある近傍）にある不変部分. なるためには，縮退している固有値を持つ固有ベ. 空間を I とする．そのとき，Span(Y ) の中に制限された. クトルがあればそれらの全てを含むか全く含まな. B-Hermite 正規直交基底 Z を，Span(Z) が I を良く近似. いかのどちらかであって部分的に含まれることが. するように構成できる．. ないようにしなければならないが，数値計算では. そのような基底 Z を用いて Rayleigh-Ritz 法により元の. 摂動によって本来は縮退していたはずの固有値の. GEVP の Ritz 対をすべて求めて，それを近似固有対とす. 縮退がほどけて分散してしまっている可能性があ. ⓒ2012 Information Processing Society of Japan. 4.

(5) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. る．もちろん数値的には縮退であったかどうかを. え理論上は正定値や半正定値であっても絶対値の小さい負. 知ることは難しいので，摂動の限界を可能なら見. の固有値を持つ固有対が混入しうる．そこで後で見るよう. 積り，値が近接している固有値は縮退している可. に固有値の絶対値が非常に小さい固有ベクトルについての. 能性があるとして扱う必要がある．結局実用的に. 関係は切断操作により捨てる．それで最終的な計算結果が. は数値で求めた固有値 φ の値を gpass 未満のもの. 大きく影響を受けない定式化になっていればそれで良い．). から（あるいはマージンをとるために gpass より. たとえば典型的な 4 種類のフィルタ Butterworth，Cheby-. も少しだけ小さい値，例えば 0.95gpass などから開. shev，inverse-Chebyshev，elliptic では f (λ) が非負の実関. 始して）減少順に探索してゆき，φ の分布の中に. 数なので，β は半正定値 Hermite 対称行列になる．ただし，. 比較的大きな空隙（すきま）のある場所を見つけ. フィルタの係数 c∞ はフィルタの種類が inverse-Chebyshev. たならば，それよりも上側にある φ はすべて近傍. あるいは elliptic の場合でさらに n が偶数のときだけは零. として含めるとするのが安全であろう．上側の上. ではない正の値を持ち，普通はその大きさが微小なので零. 限 1 についても誤差の影響を考慮して 1 を越えた. として省く近似が有利で，そのような近似を入れた場合に. φ の値も近傍に含めるようにする．近傍に取り入. は β は大きさが微小な負の固有値を持つ可能性がある．. れた最小の φ が，取り入れなかった最大の φ との. フィルタの性質から自然に導かれることであるが，通常. 間に十分距離があり，かつ 1 に比べてあまり小さ. は Y の特異値には（入力ベクトルの個数を十分大きくとれ. い値でなければ，条件の良い近傍の設定に成功し. ばそれだけ）微小な値が多く，係数 α，β は共に悪条件に. たことになる．もしも近傍に微小な φ の値が含ま. なる．そのことを考慮して，固有値方程式 α u = φ β u は. れてしまうと，あたかも gpass が最初から微小で. 以下（節 5.3）の方法を用いて特別丁寧に解く．. あったのと同じことになり，その後の計算で精度低下の原因になるから，近傍に取り入れる φ の値. 5.3. 「縮小された F の固有値方程式」α u=φ βu の解法. には下限も設けておくのが良いであろう．例えば. Y が積の形で α の表式中には 2 回，β には 1 回，それぞ. (1/2)gpass よりも小さい φ は近傍には含めないな. れ入っているので，α よりも β の方が行列の数値的な条件. ど．（ここでは gpass は，用いたフィルタの通過域. が良い．そこで Hermite 対称な行列 β を Hermite 対称行. に於ける伝達関数の「実際の」下限値である．た. 列用の Jacobi 法を用いてユニタリ変換で固有値分解する．. とえば典型的なフィルタの場合で c∞ を零に置き. これにより得られる固有値は実数であり，固有値の減少順. 換えた近似を採用している場合には，その分を減. に（負のものは後に）固有ベクトルを並べた β の固有値分. じた値を用いるものとする．）. 解を β=QH DQ とする（ここで固有値分解用いる Jacobi 法は Hermite 対称行列に 2 次のユニタリ変換を作用させて. 5.2 F の近似固有対の解法. 対角行列に収束させることで固有分解を求める方法で，実. 組 Y が張る部分空間内における F の近似固有対 (φ, v). 対称行列の場合の 2 次の直交変換（Jacobi 回転）を繰り返して対角行列に収束させる方法の拡張になっている [1]）．. は以下の方法で求めている．固有値方程式 Fv = φ v に対して，v が組 Y の張る部. ユニタリ行列 Q による変換 u ≡ Qw により G ≡ Q α QH. 分空間内にあると言う条件の式 v=Y u を課し，さらに. とおくと，方程式は Gw = φ Dw となる（行列 G は Hermite. 両辺に左側から X H B を乗じると，次数 m の小規模な. 対称である）．. GEVP：α u = φ β u が得られる．ここで α=X (BF)Y ， H. β=X BY である．いま，F がレゾルベントの線形結合で， H. 伝達関数 f (λ) は実関数としたので，BF は Hermite 対称であ. きわめて小さい閾値（たとえば 100 εMAC ）をと置き，値が未満の D の対角要素に対応する行と列を G，D，w b w=φ bw b b から省くことにより，切断された固有値方程式 G D. X F BY = (FX) BY = Y BY により行列 α は半正. b は Hermite 対称である）．さらに，対角を得る（行列 G bD b −1/2 と b 1/2 z により H ≡ D b −1/2 G b ≡D スケール変換 w. 定値 Hermite 対称で，つぎに β = X H BY = X H BFX =. おいて得られる Hermite 対称行列 H の標準固有値問題：. X H (BF)H X = X H F H BX = Y H BX = β H により行列 β が Hermite 対称であることがわかる．さらに，f (λ) が正. Hz = φz を Hermite 行列用の Jacobi 法を用いて解く． b 1/2 z b := D その固有対 (φ, z) のベクトル z から逆変換 w. 値関数ならば β は正定値であること，f (λ) が非負値関数な. b を得て，切断された行の自由度に零をにより対応する w. らば β は半正定値であること，f (λ) の負の値の下限値の絶. 補って w を得て，さらに u := Qw を得て，v := Y u を得. 対値が微小ならば β に負の固有値が存在すればその絶対値. る．この v を集めたものが不変部分空間の近似空間の基底. は微小であること，を示すこともできる．（ただし実際に. となる．但しこのままでは，まだ B-Hermite 正規直交系に. は，β の固有対を数値計算で求めると，誤差の影響により. はなっていない．. る．このことから，まず α = X H (BF)Y = X H (BF)H Y = H. H. H. H. 固有値の絶対値が小さい固有対は精度が出ないので，たと ⓒ2012 Information Processing Society of Japan. 注：H に Hermite 対称行列用の Jacobi 法を適. 5.

(6) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. 用して求めた場合には固有ベクトル z は自然に. トの作用の計算を完全に省略できて，計算量を半分に減ら. Hermite 正規直交系となるから，u が β-Hermite 正. せたが，Hermite 対称の問題の場合に於いてフィルタ F を. 規直交系になる．そのときには上の v. 複素ベクトル x に作用させる計算では，複素共役なシフト τ. (i). (i). := Y u. のかわりに単にスカラ倍による規格化を付け加え √ た v(i) := (1/ φ(i) )·Y u(i) を作るだけでこの新し. と τ を持つ 2 つのレゾルベントの作用を与える連立 1 次方程式を解く際に LU 分解を用いるのであれば，連立 1 次方程. い v(i) は B-Hermite 正規直交系になることが示せ √ る．但し，規格化因子 (1/ φ(i) ) が現れることか. 式の係数行列の間の関係 (A − τ B) = (A − τ B)H を利用す. ら容易に想像されるように，もしも通過域にある. れるならば，共役のシフト τ の係数は (A − τ B) =: U H LH. 固有ベクトルの伝達率に相当する φ. という LU 分解を持つので分解の計算を片方だけ行って共. (i). に小さい値. があると数値計算では誤差の拡大から B-Hermite 正規直交性が（規格化だけではなくて直交性も）. ることで，シフト τ で係数が (A − τ B) =: LU と LU 分解さ. 有することで計算と記憶場所が節約できる．. Hermite 対称の問題に場合にも実対称の問題の場合と. 崩れうる．フィルタの伝達関数の gpass の値が 1. まったく同様に，フィルタの入力ベクトルの組 X と出力. の付近にある（例えば典型的な 4 種類のフィルタ. ベクトルの組 Y をフィルタの伝達関数の特性パラメタに照. で gpass = 0.5 にとるなど）のであれば，そのよう. らして分析することで，B-Hermite 正規直交基底 Z を，Y. な困難は生じないが，もしもそうではなくて gpass. が張る空間の中で構成して，Z の張る空間が不変部分空間. が微小な値であるような特性の悪いフィルタを用. I の良い近似となるようにできる．ただし不変部分空間 I. いたならば，B-Hermite 正規直交性が数値丸め誤. は元の GEVP の指定した区間 [a, b]（の適切な近傍）にあ. 差の拡大により崩れるので，せめてそれを補正す. る固有値を持つ固有ベクトルを含むものである．そうして. るための B-Hermite 正規直交化を行う必要が生. 構成された基底 Z を Rayleigh-Ritz 法に適用することで，. じうる．. 指定された固有値を持つ固有対の近似が一斉に得られる．. 「Span(Y ) 内の. Av=λBv，v=Y u；. (λ, v). 一般固有値問題」. ↓. ↓. 「Span(Y ) 内の. F v=φv，v=Y u，φ=f (λ)；. (φ, v). F の固有値問題」. l. l. 「縮小された F の. α u=φ β u；. (φ, u). 固有値問題」. α ≡ Y H BY ，β ≡ X H BY. 参考文献 [1]. [2]. 図 2 対応関係の図. [3]. 6. まとめ行列 A，B が実対称で B が正定値の場合の一般固有値. [4]. 問題に対するフィルタ対角化法の計算手順をほぼ忠実に模倣して拡張することにより，行列 A，B が Hermite 対称で. [5]. B が正定値の場合の一般固有値問題に対するフィルタ対角化法の計算手順を導いた．実対称と Hermite 対称の両方の問題で共通する，固有値. [6]. がすべて実数であることを用いると，フィルタの伝達関数. f (λ) を実関数（実数の引数に対して関数値が実数である. [7]. 関数）となるように選ぶならば，実対称の問題の場合に用いたフィルタの設計と構成の方法は Hermite 対称の問題にもそのまま利用できる．そのように構成されたフィルタ F. [8]. は実対称の問題では BF が実対称で，Hermite 対称の問題では BF が Hermite 対称となる．伝達関数として実関数を用いる場合には，実対称の問題の場合にはフィルタ F を実ベクトル x に作用させる計算で. [9]. Tibor Fiala: Unitary Jacobi-Method for the Eigenvalue Problem of an Arbitrary Normal Matrix, Ann. Univ. Sci. Budapest. Sect. Comput. Vol.3 (1982), pp.119-125 (1983), MathSciNet. Toledo, S. and Rabani, E.: Very Large Electronic Structure Calculations Using an Out-of-Core FilterDiagonalization Method, J. Comput. Phys., Vol.180, No.1 (2002), pp.256–269. Sakurai, T. and Sugiura, H: A Projection Method for Generalized Eigenvalue Problems Using Numerical Integration, J. Comp. Appl. Math., Vol.159 (2003), pp.119– 128. 村上弘: レゾルベントの線形結合によるフィルタ対角化法, 情報処理学会論文誌：コンピューティングシステム, Vol.49, No.SIG2 (ACS21) (2008), pp.66–87. Polizzi, E.: Density-Matrix-Based Algorithm for Solving Eigenvalue Problems, Phys. Rev. B, Vol.79, No.11 (2009), p.115112[6pages]. 村上弘: 固有値が指定された区間内にある固有対を解くための対称固有値問題用のフィルタの設計, 情報処理学会論文誌：コンピューティングシステム (ACS31), Vol.3, No.3 (2010), pp.1–21. 村上弘: フィルタで濾過されたベクトルの組から不変部分空間の直交基底の組を近似構成するフィルタ対角化法, 情報処理学会研究報告, Vol.2011-HPC-129, No.1 (2011), pp.1–8. 村上弘: 対称一般固有値問題のフィルタ作用素を用いた不変部分空間の近似構成, 先進的計算基盤システムシンポジウム SACSIS2011 論文集 (電子版), (2011), pp.332–339．村上弘: 対称一般固有値問題のフィルタ作用素を用いた不変部分空間の近似構成, 情報処理学会論文誌：コンピューティングシステム (ACS35), Vol.4, No.4 (2011), pp.1–14.. は複素共役なシフト τ と τ を持つ 2 つのレゾルベントの項. γR(τ ) と γ R(τ ) の間の複素対称性から片方のレゾルベン ⓒ2012 Information Processing Society of Japan. 6.

(7) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. 付. 録. の減少順に並べたときに，第 m + 1 番目の伝達率が（伝達. A.1 入力ベクトルの個数 m の決め方の考察入力ベクトルの個数 m の決め方について，多少詳しく考察をしてみる．フィルタの通過域 [a, b] に両側の遷移域を加えて広げた区間を [a0 , b0 ] とする．これは阻止域以外の区間である．以下では，固有値が [a, b] にある固有対の個数を N とし，固有値が [a , b ] にある固有対の個数を N とする． 0. 0. 0. m 個の B-Hermite 正規直交ベクトルを入力ベクトルの組 X とし，X にフィルタを作用させて得られる出力ベクトルの組を Y とする．いま阻止域に於けるフィルタの伝達率の大きさの上限 gstop が数値的に微小で零とみなせる近似が十分良く成り立つとすれば，Y の計量 B での数値的な階数は N を超えない． 0. フィルタの伝達率の最大値は 1 となるように規格化しているので（注：恒等演算子の項を省く近似を入れた場合には 1 より少し小さい），阻止域における伝達率の大きさの上限値 gstop が浮動小数点数のマシンイプシロン（IEEE 754 の倍精度 64bit 浮動小数点数では 2.22×10−16 ）の程度の 10−15 や. 10. −14. などであれば丸め誤差と同程度でほぼ零と. 見なせる．またマシンイプシロンよりもある程度大きな値 10−12 や 10−8 などであっても，最初からその程度の演算精度を用いて計算をしていると思うならば，ほぼ零と見なせる値であると思うことができよう．もしも m < N 0 であれば，N 0 個ある「固有値が [a0 , b0 ] にある固有ベクトル」のすべては，Y の m 個のベクトルの線形結合で表すことはできない．するといま必要としている「固有値が [a, b] にある固有ベクトル」も，うまく近似できないかまたは近似の精度が十分にできない可能性がある．乱数ベクトルから計量 B の正規直交化で作った X は，ふつうは一般的（generic）な状況にあるので，N 0 ≤ m であれば Y の計量 B での数値的階数はちょうど N になる． 0. 実際には m の値をちょうど N にせずに，いくらか大き 0. い値とするのがよい．その理由は，乱数を用いることで一般的（generic）な状況を作り出してはいるが，その数値的な generic の程度（一種の条件数）があまり良くない状況が生じる確率をできるだけ低くするためである．「固有値がフィルタの遷移域にある固有対」が多く存在している場合で（伝達関数の振る舞いが急峻であっても遷移域に固有値が密集して分布していればそのような状況が生じうる），N ≤ m とはならずに N ≤ m < N となっ 0. 0. た場合には，区間 [a0 , b0 ] にある固有値 λ を重複を込めて考え，それらの固有値に対応する伝達率 f (λ) をその絶対値. ⓒ2012 Information Processing Society of Japan. 関数の最大値 1 に比べて）「数値的に微小で零と見なせる」かどうかによって，近似対の精度が左右されるであろう．以上の考察に従えば，フィルタにより「阻止域に固有値がある固有ベクトル」がほぼ完全に除去されているとみなせるならば，m として「阻止域ではない区間 [a0 , b0 ] に固有値がある固有対」の個数 N 0 よりも少し大きい値を採用すればよいことになる．実際の入力ベクトルの個数 m の決め方としては，たとえば以下のような方法が考えられる．. ( 1 ) 何らかの根拠か経験に基づいて，N 0 の値かその上限値が既知の場合には，それに基づいて N 0 < m を満たす少し大きめの m を設定する．そうしてこの m の値を用いても以降の計算で必要な記憶量や演算量が過大にはならないと判断したら計算を進める（過大である場合には，元の区間 [a, b] を分割して扱うことを検討する）．. ( 2 ) N 0 の値について予備知識が無いか，推定値に誤差が伴う場合には，入力ベクトルの個数 m の値を（それに伴い必要となる記憶量や演算量が過大にはならない範囲で）様子を見ながら少しずつ段階的に増やして計算を行うことができる（増やす場合には計算を最初から完全にやり直すのではなくて，m を増やす前の部分の. X の列は保ちながら，X に新たな列を追加して，追加後の X 全体が B-Hermite 正規直交性を保つようにできる．そうして X の追加した列にあらたにフィルタ F を作用させて Y の対応して追加する部分の列を作る．このように追加的に処理を行うことが可能である）．一般的な状況にある X にフィルタを作用させて作った Y の数値的な階数が十分に落ちていない場合や，得られた近似固有対の残差が十分に小さくない場合には，m を増やす（m を増やし続けてついにある限界を超えたら，区間 [a, b] の再分割を検討する）．. ( 3 ) Hermite 対称な定値一般固有値問題 Av=λBv では，指定した区間に固有値がある固有対の個数は Sylvester の慣性律を用いて求められるこれは固有値を 2 分法などの区間縮小法で求めるのに用いられてきた方法である．いま実数のシフト σ を与えて Hermite 対称行列 C := A − σB を作り，Hermite 対称行列用の修正. Cholesky 法で C =: LDLH と分解できたとすると，σ よりも小さい固有値の個数は実対角行列 D の負の対角要素の個数に等しい．このことを用いると区間 [a0 , b0 ) に固有値がある固有対の個数 N 0 は σ が b0 と a0 の場合の分解で得られた対角行列 D の負の符号数の差として求められる．ただし，区間の両端のシフトにおいて分解が破綻せずにできたと仮定する（破綻が生じる可能性もあるが，区間内の固有対の個数の「上限」が求まれば良いだけの用途であれば，破綻が起きたシフ. 7.

(8) 情報処理学会研究報告. Vol.2012-HPC-134 No.1 2012/6/1. IPSJ SIG Technical Report. トの値を区間を広げる方向に少し移動して計算をやり直してみて破綻が起きなければそれを代替に用いるとそれで求めた符号数の差は上限を与える．但し数値安定性については検討を要する）．この方法はたとえば. A，B が幅の狭い帯行列であれば高速に実行しうる. A.2 「縮小された F の固有値問題」α u=φ βu の固有ベクトル展開を用いた分析以下では，縮小された固有値方程式 αu=φβu の解法についてより良く理解するため，固有ベクトル展開を用いた分析を行う．いま f (λ) は実数値をとると仮定している．あらかじめ伝達率 φ=f (λ) が（広義）減少順となるように，元の固有値問題の固有対には番号が付けられているとする．そうして第 j 列に第 j 番目の固有ベクトルを並べた N 次の正方行列を V とする．固有ベクトルは B-Hermite 正規直交条件で規格化されている，つまり V H BV =IN であるとする，ここで IN は N 次単位行列である．そうして固有ベクトルに対応するように元の固有値問題の固有値を並べた N 次実対角行列を Λ とすると AV =BV Λ である．さらに第 j 番目の固有ベクトルの伝達率 φ. (j). を第 j 番目の対角要素と. する N 次実対角行列を Φ とおくと FV =V Φ である．フィルタへの m 個の入力ベクトルの組 X は N ×m 行列で，X は B-Hermite 正規直交化されているので X BX=Im で H. ある．いま X の固有ベクトルによる展開係数を表す N ×m 行列 C が存在して X=V C と書けて，X の B-Hermite 正規直交性から C H C=Im であるから C はフルランクである．そうして Y = FX = FV C = V ΦC とかける．すると，「縮小された F の固有値問題」の係数である m 次の Hermite 対称行列はそれぞれ α = Y H BY = C H Φ2 C ，. β = X H BY = C H ΦC と書けることがわかる．いま，m が十分に大きくて，フィルタの伝達特性の性質から第 m + 1 番目以降の伝達率 φ. (j). の大きさは無視できるほど十分に小. さいと仮定する（たとえば阻止域に於ける伝達関数の大きさの上限値 gstop が無視しうるほど小さいと見なせるなら. b の第 j 対角要素の値 φ(j) になることが分目の固有値は Φ b は既知の量ではないから，固有値方かる．ただし，C や C 程式 α u=φ β u をなんらかの数値的手段で直接解いて近似解を求めることになる．その際には固有値 φ が微小な値を持つ近似解については丸め誤差の影響を受けて良い精度では求められないであろう．しかし今の場合には固有値 φ が. [gpass , 1]（の近傍）である固有対 (φ, u) だけを解くのであるから，フィルタの設計の際に gpass として（たとえば 0.5 程度の）微小ではない値が設定できていれば，数値計算上の困難の程度は少ない．（もしもそうではなくて，gpass の大きさが 1 に比べて微小な値であるような特性の良くないフィルタを用いると，丸め誤差の影響で困難が生じる可能性がある．）付記：いま比較のために，上記と同様にして，元の固有値問題 Av=λBv を v ∈ Span(Y ) に制限する近似で得られる方程式 AY u=λBY u の残差と X とで標準の Hermite 内積で方程式を作ってみる．するとそれは小さい次数 m の行列 γ=X H AY ，β=X H BY を係数とする Hermite 対称な一般固有値問題 γ u=λ β u になる（行列 β は前と同様に Hermite 対称であり，また γ=C H Λ Φ C が示せるので，行列 γ も Hermite 対称である）．Φ の切断近似の下で， b u=w とおくと，対角行列を係数とする一般固有値問題 C. bΦ b w=λ Φ b w になり，その第 j 番目の固有値は一応「数学 Λ 的には」元の固有値問題の固有値 λ(j) に等しいことが分か b は既知の量ではないので，Hermite 対る．ただし，C や C 称な一般固有値問題 γ u=λ β u は実際には何らかの数値的手段により直接解くことになるが，その場合に φ(j) の値が微小であればそれだけ丸め誤差の影響により対応する λ(j) の近似値の精度が失われる可能性がある（しかも実際には，. γ u=λ β u だけを単独で解いても，得られた各近似対 (λ, u) がどの番号 j と対応するか，また対応する φ(j) の値が何であるかはこの場合にはわからないのである）．極端な状況に於いては，この近似固有値 λ は元の固有値問題の固有値から大きく離れた（無縁な）値になる可能性もある．. ば，m の値を固有値が通過域あるいは遷移域にある固有対の個数よりも大きくすればよい）．そのとき，行列 Φ をそ b までで切断する近似を行う．それにの m 次の首位行列 Φ 対応して C をその最初の m 行目までで切断した m 次正方 bC b と近似で b ，β ≈ C bH Φ b2 C b を用いると，α ≈ C bH Φ 行列 C. b があまり悪条件ではないと仮定するきる．さらに行列 C （X が乱数を元にして（B-Hermite 正規直交化を施して） b はあまり悪条件に作ったベクトルの組であれば，通常は C はならないはずである．もしもそうでなければ乱数をとり直すかもしくは個数 m を増やせば改善することが期待できる）．すると「縮小された F の固有値問題」α u=φ β u は， b u=w とおくと，対角行列を係この Φ の切断近似の下で C. b 2 w=φ Φ b w となり，その第 j 番数とする一般固有値問題 Φ ⓒ2012 Information Processing Society of Japan. 8.

(9)