対称固有値問題のフィルタ対角化法の実験と考察

全文

(1)Vol.2010-HPC-127 No.1 2010/10/13. 情報処理学会研究報告 IPSJ SIG Technical Report. が区間から離れた固有ベクトルはほとんど含まれないので，F は固有値が区間近傍にある固有ベクトルの組で張られた比較的小さい次元の不変部分空間への近似的な射影を与える．. 対称固有値問題のフィルタ対角化法の実験と考察. 十分多くのランダムなベクトルの組を，線形独立性を高めるために計量 B で正規直交化したベクトルの組をフィルタへの入力とする．フィルタからの出力ベクトルの組を計量 B. 村. 上. 弘†1. で特異値分解し，誤差の拡大を防ぐために特異値が相対的に小さい特異ベクトルは取り除「部分空間法」をく．除かれずに残った（B-正規直交である）特異ベクトルの組に対して，適用すれば，固有値が区間 I の近傍にある固有対の近似が得られる．. レゾルベントの線形結合をフィルタとするフィルタ対角化法により，対称固有値問題の固有値が指定区間内にある固有対だけを解く実験の例を紹介する．近似対をフィルタ通過後のベクトルの組から求める過程についても考察する．. 以上のフィルタ対角化法により得られた近似対の残差ベクトルを求めて，精度を改良する必要があれば，Rayleigh 逆反復もしくは Ritz 同時逆反復を数回程度適用する．. 2. レゾルベントの線型結合によるフィルタ作用素. Experiments and Considerations on the Filter Diagonalization Method for Symmetric Eigenvalue problems. 実対称定値一般固有値問題 Av = λBv に対応して，レゾルベント R(τ ) = (A−τ B)−1 B と恒等演算子 I の線形結合でフィルタ作用素 F = γ∞ I +. 2n. p=1. γp R(τp ) を作ると，固有値問題. の任意の固有対 (λ(ν) , v(ν) ) に対して，Fv(ν) = v(ν) f (λ(ν) ) が成立する．これは，元の固有値. Hiroshi Murakami†1. 問題の固有対の固有ベクトルがフィルタ作用素 F の固有ベクトルでもあり，その固有値は f (λ) になっていることを意味する．ここで f (λ) は λ の有理関数 f (λ) = γ∞ +. 2n. p=1. γp /(λ − τp ). でフィルタ作用素の伝達関数（transfer function）と呼ぶ．. By the filter diagonalization method whose filter is a linear combination of resolvents, some results of experiments are shown to solve only those pairs of a symmetric eigenproblem whose eigenvalues are in the specified interval. Considerations are also made on the process to calculate the approximated eigenpairs from the set of filtered vectors.. いまこの伝達関数の通過帯域を λ ∈ [a, b] として，それを標準区間 t ∈ [−1, 1] に対応させる線形変換を λ = L(t) =. b−a 2. t+. b+a 2. とする．さらに，規格化された座標 t による伝達関. 数 g(t) = f (λ) = f (L(t)) を定義する．g(t) は実軸上で有界な有理関数で，区間 [−1, 1] の特性関数（区間内で値が 1，それ以外では値が 0）をよく近似するように決める（特性関数の有理関数による近似法は，電子回路のアナログフィルタ理論で古くから用いられてきた手. 1. フィルタ対角化法の原理. 法である2),3),5) ）．. 行列 A，B が実対称で，B が正定値の場合の一般固有値問題 Av = λBv で，固有値が指. g(t) の逆数である A(t) = 1/g(t) を減衰関数（attenuation function）と呼ぶ．そうして，. 定区間 I = [a, b] にある固有対を求めるものとする．行列 A，B の次数 N は極めて大きい. 通過帯域（passband）|t| ≤ 1 では 1 ≤ A(t) ≤ Amax（この上限は tight とする）を要求し，. とする．. 少し離れた外側の阻止帯域（stopbands）|t| ≥ µ > 1 では Amin ≤ A(t) を要求する（図 1）．. フィルタ F は，固有値が区間内 [a, b] にある固有ベクトルは良く通過させるが，固有値が. 両者の中間の領域 1 < |t| < µ は遷移帯域（transitionbands）と呼ぶ．. 区間から離れた固有ベクトルは強く減衰させる性質を持つようにうまく構成された線形作. 典型的な 4 種類のフィルタとして Butterworth，Chebyshev，inverse-Chebyshev，ellip-. 用素であるとする．このフィルタに任意の入力ベクトルを通すと，出力ベクトルには固有値. tic がよく知られている3),5) ．これらはすべて t の偶関数で，種類に対応して有理関数 A(t) の関数形が決まる．形状条件の 3 つ組 (µ, Amax , Amin ) を指定すると，それを実現できる次数 n の最小値が決まる．その最小値以上の次数 n を指定すれば，条件を満たす関数 A(t) が. †1 首都大学東京数理情報科学専攻 Department of Mathematics and Information Sciences, Tokyo Metropolitan University. 完全に決まり，その逆数である g(t) の部分分数展開 g(t) = c∞ +. 1. 2n. p=1. cp /(t − tp ) の複素. c 2010 Information Processing Society of Japan .

(2) Vol.2010-HPC-127 No.1 2010/10/13. 情報処理学会研究報告 IPSJ SIG Technical Report 200. elliptic passband stopband. Attenuation[dB]. 150. 100. 50. 0. 0. 0.5. 1. 1.5. 2. t. 図 1 減衰率関数の形状パラメタ Fig. 1 Shape Parameters for Attenuation Function.. 図 2 楕円フィルタの減衰率関数 A(t) のグラフの例 Fig. 2 Sample Graphs of Attenuation Functions of Elliptic Filter. (µ=1.3, Amax =10[dB], Amin =100[dB], n=9.). 数の極と係数の数値は解析的な式に基づいて計算できる5),13)–15) ．フィルタ作用素 F のレゾルベントのシフト量 τp と線形結合係数 γp ，γ∞ は，λ=L(t) で，. β も極めて悪条件の行列となるので，縮小された固有値問題を数値的に良い精度で解くこと. g(t)=f (λ) であることから，τp =L(tp )，γp =L · cp ，γ∞ =c∞ と決まる（注記：文献13),14). が困難になる．そこで，Y の計量 B の特異値分解と特異値の閾値による切断による正則化. ではこの f (λ) と g(t) の極の係数 γp と cp の間の対応関係の記述に簡単な誤りがある）．任. を行って，より条件の良い部分空間の基底の組を作り，部分空間法を適用する必要がある．. 意の区間 [a, b] を通過帯域とするフィルタ作用素 F は，標準座標 t での g(t) の極とその係数. 3.1 部分空間の基底の構成法. をあらかじめ表にして用意しておけば，分点と係数の線形変換により簡単に構成できる．同. フィルタからの出力ベクトルの組を Y とする．Y の列で張られた部分空間内の固有解の. じ形状条件の 3 つ組に対して 4 種類の典型的フィルタのうちで，次数 n の最小値が最も小さ. 近似を求めるために部分空間法を用いると，縮小された行列 α = Y T AY と β = Y T BY を. いものは elliptic である．（楕円フィルタの減衰率関数 A(t) の例を µ=1.3，Amax =10[dB]，. 係数とする小規模の対称一般固有値問題 αz = λβz の固有対 (λ, z) に対応して，元の問題. Amin =100[dB] で，最小の次数 n=9 にとった場合のグラフを図 2 に示す．）. の近似対は (λ, Y z) となる．しかし，Y はその性質から一般には階数がかなり多く落ちた・落ちかけている条件の極め. 3. 部分空間法に対する考察. て悪いベクトルの組で，β = Y T BY は悪条件が 2 乗となり，β を（念を入れるならば対称. 元の一般固有値問題 Ax=λBx をフィルタ F の出力ベクトルの組 Y の張る部分空間内で. な対角ピボット交換付きの）Cholesky 法等で対称分解し一般固有値問題を標準固有値問題. 近似する部分空間法は，縮小行列を α=Y T AY ，β=Y T BY とおくと，元の問題の固有ベク. に帰着する通常の解法で解いても数値的困難が予想される. トル x が Y の張る部分空間に含まれていれば x=Y w で，縮小された固有値問題 αw=λβw. B 自身の条件が良い場合には，β の悪条件の原因は，部分空間の基底の組 Y の悪条件か. には固有対 (λ, w) がある．. ら来るので，基底のとり直しや切断正則化を行って線形独立性を改善し，新しい基底の組で縮小した係数行列の条件が悪くならないようにすることが考えられる．. しかしフィルタからの出力の組 Y は，その由来から通常は極めて悪条件であり，そのため. その一つの方法は，Y を（閾値に基づいた）有効階数が決定できる列ピボット交換付き. 2. c 2010 Information Processing Society of Japan .

(3) Vol.2010-HPC-127 No.1 2010/10/13. 情報処理学会研究報告 IPSJ SIG Technical Report. の Householder-QR 法で分解して，単位計量での正規直交基底の組 Q を作り，Q が張る部. 精度を最大で d 桁失なっているはずである．. 分空間に縮小された行列 α = QT AQ と β = QT BQ を係数とする小規模の縮小された一. 今回は，フィルタを 1 回だけ作用させる定式化にしているが，例えば 2 段階に分けて適. . 般固有値問題を解くことである．もしも縮小された係数行列 β があまり悪条件でなければ，. 用（文献例 11)）するのであれば，仮に 2 段階で実現されるトータルの減衰率特性が同等で. 通常の β の対称分解を用いた一般固有値問題の解法が使えることになる．. あっても，得られる不変部分空間の基底は近似品質が向上するはずである．第 2 回目のフィ. 今現在は，反復改良付きの古典的な B-計量の特異値分解による方法8) を用いて，Y の B-. ルタへの入力の組は，第 1 回目の入力の組であるランダムベクトルの組の B-正規直交化し. 特異ベクトルの組（B-正規直交である）Q を以下の方法で作成している．. たものではなくて，必要な固有ベクトルがほとんどの割合を占めるベクトルの組を B-正規. 計量 B の SVD を用いた正規直交化：. 直交化したものになるので，第 2 回目のフィルタ出力の B-SVD で得られる特異ベクトルの. (1). Y (0) = Y とする．. (2). Y (0) の自己 B-内積の行列 S (0) = Y (0) BY (0) を作る．対称行列 S (0) の固有値分解が S. (3). 特異値は，必要な固有ベクトルを多く含む特異ベクトルについては 1 に近いので，相殺によ T. (0). =U. (0). D. (0). U. (0) T. であるとして，Y. (1). =Y. (0). U. (0). る情報落ちの影響が少ない．. を作る．. 4. 丸め誤差の影響の拡大率の見積り. T. Y (1) の自己 B-内積の行列 S (1) =Y (1) BY (1) を作る（S (1) =D(0) となるはずであるが，誤差の影響でずれる）．対称行列 S (1) の固有値分解が. S (4). (1). =U. (1). D. (1). U. (1) T. であるとして，Y. (2). =Y. (1). U. (1). フィルタ F の出力ベクトルの組 Y から，B-正規直交である B-特異ベクトルの組 Q を作る．直交化の後に正規化した際の拡大係数の組を横ベクトル s とする．s は特異ベクトルと. を作る．. T. Y (2) の自己 B-内積の行列 S (2) = Y (2) BY (2) を作る（S (2) = D(1) となる. 対応させて特異値の逆数を並べたものである．. はずであるが，誤差の影響でずれる）．…. Q の列が張る部分空間に，Ax = λBx に対する部分空間法を適用すると，縮小された行. この反復により Y (j) の B-直交性が改善される．通常は，最初の直交化に 1 回の改良を加. 列 α = QT AQ を対角化する直交行列を U とすると，近似固有ベクトルの組は Z=QU で. えた段階の Y (2) で十分な精度が得られる．そこで，Y (2) の各列を計量 B での長さを 1 に. ある．. 正規化したものを Q とすると，Y の B-正規直交化が完了する．Q の列ベクトルが B-SVD. U の各要素をその絶対値に置き換えた行列を U として，横ベクトル ζ=sU を作ると，ζ. による特異ベクトルで，その特異値は対応する Y (2) の列ベクトル y の計量 B での長さ. . の各要素は，Y の条件の悪さに起因して元の固有値問題の近似固有ベクトルの組 Z の各列. yT By になる．B-SVD の後で，ある閾値よりも小さい特異値を持つ特異ベクトルは，部. に混入した丸め誤差の拡大倍率の目安になる．拡大倍率の値が非常に大きいものは，大きな. 分空間の基底ベクトルの組からは除外する．. 誤差を含む可能性のある解であるか，あるいは誤差により生じた偽の解である．もちろん，元の固有値問題の近似固有対 (λ, v) の品質は，残差ベクトル r = Av − λBv. この反復改良法は，Y の自己 B-内積である対称行列 S の条件数が，特異値の最大最小比の 2 乗になることから精度が低下するのを回避するのに有効である．但し，この反復で B-. のノルムで調べることができる．. 直交性の改良が有効にはたらくためには，対称行列の固有値分解として Householder 三重. 5. 特異値の切断の閾値. 対角化を経由する方法はだめで，Rutishauser の Jacobi 法1) を用いる必要がある（悪条件の問題に対して Jacobi 法の精度が良いこと等については，文献 4) が参考になる）．. 現在用いているフィルタ対角化法では，十分多くのランダムなベクトルの組を計量 B で. この計算で得られる特異値は，対応する B-特異ベクトルの Y の組に含まれていた合計量. 正規直交化したものを，フィルタへの入力としている．フィルタからの出力ベクトルの組 Y. を表す．一般には特異値の小さい特異ベクトルほど相対精度が失われている．最初に（それ. に対して，計量 B による特異値分解を行い，丸め誤差の影響の拡大を抑えるために，特異. ほど精密でなくても良い）B-正規直交化をしておいた入力ベクトル X が，フィルタにより. 値が非常に小さい特異ベクトルを棄却する処理「切断」を加える．. 多くの不用成分が除かれて縮んで Y となる過程で数値的相殺により有効精度が失われ，そ. 数値計算では，実際に実現可能な Amin の値は演算精度で制限される．いま演算に用いる. の後の Y の B-直交化の過程でも数値的相殺で有効精度が失われる．トータルとして数値的. 浮動小数点数のマシンイプシロンを M とすると，フィルタを丸め誤差を伴う演算で適用し. 相殺による情報落ちで特異値が 10−d の特異ベクトルは，不変部分空間の基底としての有効. 3. c 2010 Information Processing Society of Japan .

(4) Vol.2010-HPC-127 No.1 2010/10/13. 情報処理学会研究報告 IPSJ SIG Technical Report. た場合の現実の減衰率が −1 M より大きくなることは期待できない．このことから切断の特. 2. 異値の閾値は M よりも「ある程度」大きい値に設定しなければならない．. 1. 固有値が阻止帯域にある不要な固有ベクトルを除去するためには，特異値の閾値は. A−1 min. 0 -1. よりも「ある程度」大きい値に設定しなければならない．固有ベクトルが入力ベクトルの組 X に含まれている強度」の最小値に. A−1 max. LOG10(DELTA). 固有値が通過帯域にある必要な固有ベクトルを残すためには，特異値の閾値は「残したいを乗じた値よ. りも「ある程度」小さくなければならない．「残したい固有ベクトルが入力ベクトルの組 X に含まれている強度」は，入力 X に依存していて事前には分からない．フィルタ対角化の. -4 -5. -7. トルが全く含まれていないかあるいはごく僅かしか含まれていない場合である．そのような. -8. 稀な状況の生じる確率は乱数ベクトルの個数を増すと小さくなる．このような統計的分布を. -9. 持つデータに対する算法の振舞いの理論的な考察が必要で，今後は統計的な議論も展開する必要がある．少なくとも閾値は，最大の特異値 σ1 に. -3. -6. 算法にとって不都合な状況は，乱数で生成された入力ベクトル X の組に求めたい固有ベク. A−1 max. -2. を乗じた値よりも「ある程度」. 60. 65. 70. 75 80 85 EIGENVALUE. 90. 95. 100. 図 3 固有対の残差ノルム（偽の固有対の出現例）指定区間 [60, 100] Fig. 3 Norm of residual of eigenpairs (spurious pair appeared). 小さく設定する必要がある． c 以上を総合すると，r = Amax / min(Amin , −1 M ) とおいて，閾値を σ1 r とするとき，指. すると，最初の r 個の特異ベクトル q (1) , . . . , q (r) までは，いずれも固有値が通過帯域に. 数 c を 1/3∼2/3 の範囲程度にするのが良いように思われるが，適切な閾値をどう決めるかの指針については，今後もより深く検討すべきである．. ある固有ベクトル r 個の相対的に係数の大きな線形結合に，固有値が遷移帯域にある固有. 5.1 偽の固有対について. ベクトル 2 個がごく弱く混ざったものになる．そうして，q (r+1) は，固有値が遷移帯域にあ. 実対称定値一般固有値問題で，区間 [a, b] に固有値がある固有対をフィルタ対角化法によ. る固有ベクトル 2 個のある任意の割合の線形結合でほとんど表される．. り求めたとき，得られた近似対の組にまれに，固有方程式の残差が非常に大きい「偽の固有. そのような B-正規直交ベクトル (r+1) 個を用いて部分空間法の縮小行列 [α]i,j = T. q (i) Aq (j) を作ると，その最後の第 (r+1) 番目の行あるいは列は，対角要素を除いてごく弱い. 対」が混入することがある．多少の実験による経験からは，偽の対が発生する傾向が高まるのは，フィルタの特性パラ. 値が並ぶが，（固有値が遷移帯域にある 2 個の固有対をそれぞれ (λr+1 , u(r+1) )，(λr+2 , u(r+2) ). メタ µ が大きいとき，フィルタに通すランダムベクトルの本数を抑え気味にしたとき，特. とし，第 (r+1) 番目の特異ベクトルがほとんど q (r+1) = cr+1 u(r+1) + cr+1 u(r+1) と表さ. 異値の切断に用いる閾値を下げたとき，などである．. れるとすれば），第 (r+1) 番目の対角の値は αr+1,r+1 = c2r+1 λ(r+1) + c2r+1 λ(r+2) ，但し. c2r+1 + c2r+2 = 1 となる．第 (r+1) 番目の行および列の非対角要素は弱いので，部分空間法. 発生のモデルとして，以下の状況を考えてみた．通過帯域 [a, b] に固有値のある固有対は. はこの αr+1,r+1 に近い値を固有値として持つ．つまり λ(r+1) と λ(r+2) の任意の中間の値. r 個で，フィルタの遷移帯域には固有値のある固有値がたとえば 2 個あるとする．フィルタから出力されたベクトルの組を Y とする．ベクトルの組 Y に含まれるベクトル. になる．遷移帯域にある 2 つの固有値 λ(r+1) と λ(r+2) がそれぞれ区間 [a, b] の外部の両側. の個数は (r+1) 以上の数であるとする（たとえば (r+10) としてよい）．Y の各ベクトル. にあれば，中間の値は [a, b] に入る可能性が高く，区間 [a, b] に入る場合にはそれに近い固. は，固有値が通過帯域にある固有ベクトルの r 個の相対的に係数の大きな線形結合に，固有. 有値を持つ「偽の固有対」が出現することになる．. 値が遷移帯域にある固有ベクトルの 2 個が（フィルタにより減衰を受けて）ごく弱く混ざっ. 図 3 は，フィルタ対角化で偽の固有対が出現した例である．固有対の残差が異常に大きい. たものになる．Y を B-特異値分解して得られた特異対のうちで，ある閾値以上の特異値を. ことで容易に判別が可能である．. 持つものが (r+1) 個になったとする．. 4. c 2010 Information Processing Society of Japan .

(5) Vol.2010-HPC-127 No.1 2010/10/13. 情報処理学会研究報告 IPSJ SIG Technical Report. 5.2 部分空間の基底の構成法 2（考察中）. 参. 大規模 N 次の対称定値一般固有値問題 Ax=λBx を解くために，レゾルベントの線形結. 考. 文. 献. 1) Rutishauser, H.: The Jacobi method for real symmetric matrices, Numerische Mathematik, Vol.9, pp.1–10(1966). 2) 渡部和: 線形回路理論, 電子回路講座第 2 巻, 昭晃堂, 1971. 3) Daniels,R.W.: Approximation Methods for Electronic Filter Design, McGraw-Hill, 1974. 4) Demmel,J. and Vesselić, K.: Jacobi’s Method is More Accurate than QR, SIAM J. Matrix Anal. & Appl., Vol.13, No.4, pp.1204-1245, (1992). 5) Lutovac,M.D., To˘sić, D.Y. and Evans,B.L.: Filter Design for Signal Processing, §12.8, Prentice Hall, 2001. 6) Toledo, S. and Rabani, E.: Very Large Electronic Structure Calculations Using an Out-of-Core Filter-Diagonalization Method, J. Comput. Phys., Vol.180, No.1, pp.256–269 (2002). 7) Sakurai,T. and Sugiura,H: A projection method for generalized eigenvalue problems using numerical integration, J. Comp. Appl. Math., Vol.159, pp.119–128 (2003). 8) 村上弘: 非常に細長い大規模行列に対する記憶参照局所性が高い正規直交化法の実験, 日本応用数理学会論文誌, Vol.17, No.4,pp.399–454,(2007). 9) 村上弘: 帯対称定値一般固有値問題のフィルタ対角化法の実験, 情報処理学会研究報告, 2007-HPC-110(6), pp.31–36 (2007). 10) 村上弘: レゾルベントの線形結合によるフィルタ対角化法, 情報処理学会論文誌コンピューティングシステム, Vol.49, No.SIG2 (ACS21), pp.66–87 (2008). 11) Polizzi, E.: Density-matrix-based algorithm for solving eigenvalue problems, Phys.Rev.B, Vol.79, No.11, p.115112[6pages] (2009). 12) Ikegami,T., Tadano,H., Umeda,H. and Sakurai,T.: Hierarchical parallel algorithm to solve large generalized eigenproblems, HPCS2010 論文集, pp.107–114 (2010). 13) 村上弘: フィルタ対角化法の帯域通過フィルタの最適化, 情報処理学会研究報告, Vol.2010-HPC-124, No.3 (2010). 14) 村上弘: 楕円フィルタによる実対称定値一般固有値問題のフィルタ対角化法の実験, 情報処理学会研究報告, Vol.2010-HPC-125, No.1 (2010). 15) 村上弘: 固有値が指定された区間内にある固有対を解くための対称固有値問題用のフィルタの設計, 情報処理学会論文誌コンピューティングシステム (ACS31), Vol.3, No.3, pp.1–21 (2010 年 9 月), to appear. 16) Hiroshi Murakami: Filter Diagonalization Method by Resolvents for Symmetric Eigenproblems, Poster presentation at ICMS2010 on Sep.16th,2010 at Kobe University (2010). （ICMS2010 で会場配布される予定のソフトウェア DVD 集の一部として，フィルタ対角化法の参考ソースコードと説明のスライドを収納）．. 合によるフィルタを用いて部分空間の基底を構成する以下の方法を，現在検討中である（まだ実験しておらず，内容が誤っている可能性がある）．. X は B-正規直交化されたベクトルの組とし，それをフィルタ F を通して得られた出力ベクトルの組を Y = F X とする．いま，フィルタ作用素の固有値問題 Fy = φy をベクトル y を Y で張られた部分空間内に制約して解くことを考える．そこで，y = Y u と置くと，. F Y u = φY u となり，この式の両辺に，左から X T B を乗じれば，X T BF Y u = φX T BY u となる．（レゾルベントの定義 R(τ ) = (A − τ B)−1 B と A，B の対称性から BR(τ ) =. R(τ )T B で，F がレゾルベントの線形結合だから BF = F T B が成り立つので）左辺は X T BF Y u=X T F T BY u=(F X)T BY u=Y T BY u となる．縮小された行列を α=Y T BY ， β=X T BY とおくと，縮小されたフィルタ作用素の固有値問題は αu=φβu となる．ここで，係数行列の α は対称で，また β=X T BY =X T BF X=X T F T BX=Y T BX=(X T BY )T =β T であるから β も対称である．β の条件数は，Y の条件数の 2 乗ではなく 1 乗になる．（注：もしも β が負の固有値を持たないことを要求する場合は，フィルタの伝達関数が非負値でなければならない．すると c∞ が零でなければ，それを係数とする項を伝達関数あるいはフィルタ作用素の構成から省略することはできない．）. y が元の固有値問題 Ax = λBx の固有ベクトルならば，その固有値を λ とすると φ=f (λ) だから，縮小されたフィルタ作用素の固有値問題の固有値 φ は，y の固有値が通過帯域にあれば 1 付近で 1/Amax 以上の値を持ち，y の固有値が阻止帯域にあれば 0 付近で 1/Amin 以下の値を持つ（注：フィルタ作用素の c∞ の項を省略する場合には，限界値は c∞ だけ減少する）．そこで，縮小されたフィルタ作用素の固有値問題 αu=φβu の固有対 (φ, u) で固有値 φ が 1 付近にあるものを解いて，y=Y u とおくと，y は固有値が通過帯域にある固有ベクトルである．（u が β-正規直交に規格化された固有ベクトルの場合には，y を B-正規直 √ 交ベクトルにするには，y = (1/ φ) · Y u とすればよい）．そのような y を集めれば，固有値が通過帯域内にある固有ベクトルだけで張られた不変部分空間を近似する基底の組が得られるので，それを用いたサブスペース法により，元の固有値問題 Ax = λBx の近似対を得ることができる．そこで，縮小されたフィルタ作用素の固有値問題 αu=φβu は，β が悪条件であることを考慮して精度良く β=U DU T と固有値分解し，u=U z，γ=U T αU とおけば，方程式 γz=φDz となる．絶対値が閾値以下の D の対角要素を切断して 0 と置き，方程式をそれと整合するように修正して解く．必要な固有対はあくまでも固有値 φ が 1 付近のものだけである．. 5. c 2010 Information Processing Society of Japan .

(6)