• 検索結果がありません。

レゾルベントの線形結合をフィルタに用いたエルミート定値一般固有値問題のフィルタ対角化法

N/A
N/A
Protected

Academic year: 2021

シェア "レゾルベントの線形結合をフィルタに用いたエルミート定値一般固有値問題のフィルタ対角化法"

Copied!
16
0
0

読み込み中.... (全文を見る)

全文

(1)情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). レゾルベントの線形結合をフィルタに用いた エルミート定値一般固有値問題のフィルタ対角化法 村上 弘1,a) 受付日 2013年7月29日, 採録日 2013年11月15日. 概要:レゾルベントの線形結合をフィルタとして用いる「フィルタ対角化の方法」により,係数行列 A と B が実対称で B が正定値である一般固有値問題 Av = λBv の指定された実区間に固有値がある固有対す べての良い近似が得られることをこれまでの研究で示してきた.今回は,係数行列 A と B が(複素)エル ミートで B が正定値である一般固有値問題に対しても,これまでの方法を自然に拡張した同様の「フィル タ対角化の方法」により,指定された実区間に固有値がある固有対のすべての良い近似が得られることを 示す. キーワード:一般固有値問題,フィルタ対角化,エルミート. Filter Diagonalization Method for a Hermitian Definite Generalized Eigenproblem by Using a Linear Combination of Resolvents as the Filter Murakami Hiroshi1,a) Received: July 29, 2013, Accepted: November 15, 2013. Abstract: In previous researches, we have shown that for the generalized eigenproblem Av = λBv whose coefficient matrices both of A and B are real symmetric and B is also positive definite, the filter diagonalization method whose filter is a linear combination of resolvents gives good approximations of all those eigenpairs whose eigenvalues are in the specified real interval. In this research, we will show that for the generalized eigenproblem whose coefficient matrices both of A and B are complex hermitian and B is also positive definite, the filter diagonalization method in a similar manner by a natural extension to the previous one can also give good approximations of all those eigenpairs whose eigenvalues are in the specified real interval. Keywords: generalized eigenproblem, filter diagonalization, hermitian. 1. はじめに いま行列 A と B はエルミートで,B は正定値でもある とする.このとき(複素)エルミート定値一般固有値問題. 張により固有ベクトルは一般には実ベクトルにはならない が,固有値 λ は実対称の場合と同様に必ず実数になるので, 必要な固有対の固有値の存在範囲を実数の区間を用いて指 定することができる.. Av = λBv は,実対称定値一般固有値問題(A と B が実. この一般固有値問題(GEVP)に付随するベクトルの間. 対称で,B が正定値の場合)の係数行列を実対称の場合か. の B-エルミート内積を u, v = uH Bv とする.固有ベク. らエルミートの場合へ自然に拡張したものである.この拡. トルの全体は必ず空間全体を張り,この内積に関して正規. 1. 直交基底(B-正規エルミート直交基底)となるようにとれ. a). 首都大学東京大学院数理情報科学専攻 Department of Mathematics and Information Sciences, Tokyo Metropolitan University, Hachioji, Tokyo 192–0397, Japan mrkmhrsh@tmu.ac.jp. c 2014 Information Processing Society of Japan . る.ここで,右肩の H はエルミート転置(共役転置)の 操作を表している.そのように選んだ固有ベクトル全部を. 57.

(2) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). 並べた行列を V ,固有ベクトルに対応する順に固有値を並. 伝達関数 f (λ) が実関数であれば,BF はエルミートに. べた対角行列を Λ とすれば,AV =BV Λ,V H BV =I を満. なる.実際,c∞ は実数であり,任意の極の添字 p にはそ. たす.. れに対応する添字 q が唯一存在して,τp =τq ,γp =γq とな. 以上のことから,実対称の問題に対する手順からエルミー. ることを用いると,. トの問題に対する手順が B-内積や転置を B-エルミート内 積やエルミート転置にそれぞれ置き換えるなどにより,自. (BF)H = B c∞ I +. 然な拡張として得られることが期待できる.そこで本研究. = B c∞ I +. ではこの方針に基づいて,実対称定値一般固有値問題に対. = B c∞ I +. してすでに得ているフィルタ対角化法の手順 [4], [6], [7], [8]. = B c∞ I +. をできるだけ自然な形で模倣した拡張により,エルミート. . (A − τp B)−1 B }H. p. γp B H (A − τp B)−H B H. p. γp B(A − τp B)−1 B. q. γq B(A − τq B)−1 B.  . . q. γq (A − τq B)−1 B }. = BF. そうしてその拡張の結果として得られた,エルミート定値 が,実対称定値の場合と同様にうまく機能することを,数. p { γp B. = B{ c∞ I +. 定値一般固有値問題に対する手順を得ることを試みた [9]. の場合の一般固有値問題に対するフィルタ対角化の手順. . により示せる. フィルタ F あるいは伝達関数 f (λ) を定義するパラメー. 値実験でいくつかの例について確認した.. タである次数 2n,実係数 c∞ ,複素シフト τp ,複素結合係. 2. フィルタとその伝達関数. 数 γp は,伝達関数の形状に対して課す制約条件を満たす. 大規模な N 次の(複素)エルミート定値一般固有値問題. Av=λBv に対して,実区間 [a, b] を通過域とするフィルタ を用いたフィルタ対角化法により「区間 [a, b] に固有値が ある固有対」だけを選択的に求めることにする. い ま 複 素 数 の シ フ ト τ の レ ゾ ル ベ ン ト を R(τ ) =. (A−τ B)−1 B と定義する.フィルタ作用素 F として,シフト  が相異なるレゾルベントの線形結合 F = c∞ I + p γp R(τp ) を採用することにする.そのとき,元の固有値問題の固有 値 λ の固有ベクトル v に対しては Fv = f (λ)v が成立す  る.ここで λ の有理関数 f (λ) = c∞ + p=1 γp /(λ − τp ) は伝達関数と呼び,固有値が λ である固有ベクトルのフィ ルタ作用素による伝達率(拡大率)を与える.エルミート 定値一般固有値問題の固有値は必ず実数なので,フィルタ の設計を行う際には伝達関数 f (λ) の振舞いを λ が実数の 場合についてだけ考慮すればよい. 以上のことから,一般固有値問題のフィルタを用いた解 法では,実対称定値の場合に用いたのと同じ伝達関数をエ ルミート定値の場合にも使えることが分かる.固有値が実 数に限られることから,レゾルベントのシフトがすべて虚 数であれば,固有値の分布状況によらずに作用素 F の有界 性が保証され,伝達関数 f (λ) も実軸上では有界な関数と なる.後で見るように,互いに複素共役の虚数をシフトと して持つレゾルベントどうしはそれらの作用を与える連立. 1 次方程式の係数が互いのエルミート転置になるので,そ の対称性を利用すると LU 分解が共用できて分解の計算量 を減らせる.このことから伝達関数 f (λ) を実関数(実数 の引数に対して実数値を与える関数)となるようにして, 虚数のシフトが必ず複素共役対で現れるようにするのが有 利である.そこで本論文では,フィルタの伝達関数 f (λ) を実軸上に極を持たない(すべての極が虚数である)実関 数に限ることにする*1 .. c 2014 Information Processing Society of Japan . ように決める. これまでの「実対称定値一般固有値問題」に対して用い てきたフィルタの設計と構成の手法(フィルタ特性に対す る要求を満たすための伝達関数の極の個数と配置,極の係 数の決定法)は,今回の「複素エルミート定値一般固有値 問題」に対してもそのまま適用できる.. 2.1 フィルタの設計 まず,λ の正規化座標 t を λ ∈ [a, b] から t ∈ [−1, 1] へ の 1 次変換により定義する.そのとき t ∈ [−1, 1] は通過 域(passband)に,μ ≤ |t| は阻止域(stopbands)に,中 間の 1 < |t| < μ は遷移域(transitionbands)に,それぞ れ対応する.ただし,μ は 1 よりも(少し)大きい値を持 つパラメータである.フィルタの「種類」とは伝達関数. f (λ) = g(t) の「関数族」のことである. 伝達関数 g(t) の値の分布に制約条件: i) 通過域におけ る g(t) の tight な最小値は gpass で,逆に g(t) が gpass 以上 になるのは通過域に限る; ii) 阻止域における |g(t)| の上限 値は gstop である;を課する(図 1 参照) .条件 i) は g(t) ≥. gpass であることと t が通過域 [−1, 1] にあることとが同値 になることを意味する. 典型的な 4 種類のフィルタであるバターワース型,チェ ビシェフ型,逆チェビシェフ型,楕円型に対しては,制約 条件の 3 つ組(μ,gpass ,gstop )を指定すると,その条件 を満たせる整数である次数 n の最小値が決まり,n をその *1.  重複する極を持たない伝達関数 f (λ) = c∞ + p γp /(λ − τp ) が 実関数であるためには,係数 c∞ は実数で,実数の極の係数は実 数で,虚数の極にはその複素共役である極が必ず含まれていて, 互いに複素共役な極どうしの係数もまた互いに複素共役であるこ とが必要十分である.つまり任意の添字 p に対応する唯一の添字 q があり,τp =τq ,γp =γq となる(伝達関数 f (λ) の有界性を仮 定すると,極は虚数のみで個数が偶数 2n で,前半の添字が 1 か ら n までの極は虚部が正,後半の添字が n + 1 から 2n までの極 は前半の複素共役とできる).. 58.

(3) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). くと,疎行列 C を係数とする連立 1 次方程式の組 CZ=BX をベクトルの組 Z について解くことに帰着する.連立 1 次 方程式の組に対して係数行列 C の疎性を利用できる標準 的な解法としては,たとえば複数個のベクトルの組に対す る(ブロック)クリロフ部分空間法の系統の反復解法,あ るいはそれに不完全 LU 分解などの前処理を加えたものな どが考えられる.また行列 A と B が偏微分方程式の離散 化近似に由来するものであれば,問題の構造に密着したた とえばマルチグリッド法などの高速な解法が適用できる可 能性がある. 図 1. フィルタの伝達関数の絶対値 |g(t)| の概形図. Fig. 1 Conceptual graph of the filter’s transfer function |g(t)|.. 3.1 シフトが複素共役のレゾルベントの LU 分解の共用 3.1.1 これまでの「実対称定値」の場合. 値以上に設定すると g(t) が完全に決まる(同様に制約条件. これまでの実対称定値一般固有値問題の場合には,フィ. の 3 つ組を(n,μ,gpass )で与える場合には,条件を満た. ルタ F の伝達関数 f (λ) が実軸上に極を持たない実関数と. せる gstop の最小値が決まり,gstop をその最小値以上に設. すると,レゾルベントのシフトは必ず複素共役な虚数の対. 定すれば g(t) が完全に決まる.あるいはまた制約条件の 3. として現れ,互いに複素共役なシフトのレゾルベントに対. つ組を(n,gpass ,gstop )で与える場合には,条件を満た. するフィルタの線形結合係数もまた互いに複素共役になる.. せる μ の最小値が決まり,μ をその最小値以上に設定すれ ば g(t) が完全に決まる) .関数 g(t) が決まれば,有理関数. そうして実ベクトルの組 X にフィルタ F を作用させる計  算 FX = k γk R(τk )X の総和の中で,互いに複素共役で. である伝達関数 f (λ)=g(t) の持つすべての極の位置と各極. ある任意のシフトの対を τ ,τ とすれば,それに対応して. の係数が決まり,それにより通過域が λ ∈ [a, b] であるフィ. レゾルベントの結合係数は γ ,γ とかける.シフトが虚数. ルタ F のパラメータが完全に決まる.フィルタの設計法. のときにはレゾルベントは正則で存在して,実対称問題の. と構成法の詳細は文献 [5] に記述した(ただし,用いてい. 場合には A と B が実であることから複素共役対称性によ. る記号や定義は,本論文のものとは多少異なっている) .. り,R(τ ) = R(τ ) であることが容易に分かる.よって,実. 3. レゾルベントの作用の計算法について. ベクトル X へのフィルタ F の作用の計算の中で,対応す る 2 つのレゾルベントを含む項の和によるフィルタの作用. フィルタの作用を実現するには,N 次(列)ベクトルを. への寄与は,{γR(τ ) + γ R(τ )}X = 2Re {γR(τ )X} とな. m 個並べた組 X に対して,複素数のシフト τp のレゾルベ. り,互いに複素共役なシフトを持つ 2 つのレゾルベントの. ント R(τp )=(A − τp B). −1. B を各 p について作用させる計. 項を作用させた結果の和は,実際には片方の項を作用させ. 算が必要である.それは C=A − τp B とおけば,連立 1 次. た結果だけを求めてその実部として求められるので演算量. 方程式の組 CZ=BX をベクトルの組 Z について解くこと. と記憶参照量を容易に半分にできた.. に帰着する(この右辺の BX はシフトの値にはよらない量. 3.1.2 今回の「エルミート定値」の場合. である).行列 A と B がエルミートであっても τp が虚数 ならば行列 C はエルミートではないことに注意する.. 今回のエルミート定値一般固有値問題の場合にも,フィ ルタの伝達関数を実軸上に極を持たない実関数に限定する. いま A と B がともに帯行列で半帯幅が h であると,C. と,やはりレゾルベントのシフトは必ず複素共役な虚数の. もまた半帯幅が h となり,連立 1 次方程式は帯行列専用の. 対の形で現れ,互いに複素共役であるシフトのレゾルベン. 行ピボット交換付きの LU 分解法を用いることで安定に解. トに対するフィルタの線形結合係数もまた互いに複素共役. ける.その場合に分解で得られる L の下半帯幅は h で,U. になる.しかし以下に示すように,複素エルミートの場合. の上半帯幅は 2h 以下になる.. には実対称の場合と比べて,互いに複素共役の虚数をシフ. 元の一般固有値問題の固有値は実数なので,シフト τ が. トとして持つレゾルベントの性質の利用は多少複雑になる.. 虚数のときには行列 C は正則である.そこで行ピボット交. いまシフト τ に対応して行列 C(τ ) ≡ A − τ B を定義する. 換を省いて計算を行えば(その場合には LU 分解の L の下. と,A,B がエルミートであることから,関係 C ( τ ) = C(τ )H. 帯幅と U の上帯幅がともに h になる)演算量と記憶量を. が成り立つことが容易に分かる.つまり複素共役のシフト. 節約できるが,数値的な精度が悪化する可能性がある.. に対応する行列はエルミート転置になる.. 係数 A と B が帯行列ではなくて一般的な疎行列の場合で. そ こ で ,レ ゾ ル ベ ン ト を 複 素 ベ ク ト ル の 組 X に 作. あっても,複素数 τ のシフトを持つレゾルベント R(τ ) をベ. 用させる計算を見ていく.シフト τ のレゾルベントは. クトルの組 X に作用させる計算は,同様に C=A − τ B とお. R(τ ) = C(τ )−1 B であるから,組 X にシフト τ のレゾル. c 2014 Information Processing Society of Japan . 59.

(4) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). ベントを作用させる計算 V ≡ R(τ )X は,あらかじめシフ. 有対の個数よりも少し大きい値を個数 m として選ぶ.そ. トにはよらない右辺の組 X  ≡ BX を作っておき,連立 1. うして(通常は)乱数で生成した m 個以上のベクトルの組. 次方程式の組 C(τ )V = X  を解いてベクトルの組 V を求. から B-正規エルミート直交化により,m 個の(列)ベクト. めることに帰着できる.同様に,組 X にシフト τ のレゾ. ルの組 X を作成する(組 X は B-正規エルミート直交系,. ルベントを作用させる計算 W ≡ R(τ )X は,連立 1 次方程. つまり X H BX=Im である) .この X をフィルタを作用さ. . 式の組 C(τ )W = X を解いてベクトルの組 W を求めるこ. せる前のベクトルの組とする. 次に組 X の各列ベクトル x(k) ,k=1, 2, . . ., m に対して. とに帰着できる. ま ず 最 初 に ,シ フ ト τ に 対 応 す る 連 立 1 次 方 程 式. フィルタ F を適用して得られる列ベクトル y(k) =Fx(k) ,. C(τ )V = X  を係数行列の行選択付き LU 分解を用い. k=1, 2, . . ., m を対応順に並べて,m 個の列ベクトルの組 Y. て直接法により解いて V を求めるものとする.いま P を. を作成する.この Y がフィルタ F を X に作用させた後の. 行ピボット選択に対応する置換の行列として,係数行列. ベクトルの組 Y = FX である.. C(τ ) を LU 分解した結果が P C(τ ) = LU であるとすれば,. 阻止域におけるフィルタの伝達率の大きさが無視できる. V = U −1 L−1 P X  となり,これは置換操作,前進消去と後. かあるいは数値丸め誤差と比べて同等程度に小さいと仮定. 退代入により計算できる.次に,シフト τ に対応する連立. すると,フィルタを作用させた後の各ベクトル y(k) は固有. 1 次方程式の組 C(τ )W = X  も直接法で解いて W を求め. 値がフィルタの通過域または遷移域にある固有ベクトルの. る場合には,この係数行列 C(τ ) は LU 分解を新たに計算. 線形結合であると見なせるので,Y の張る空間はそのよう. しなくても,すでに得ている C(τ ) の LU 分解の結果のエ. な固有ベクトルで張られた部分空間の良い近似になる.. ルミート転置をとれば C(τ )P. T. H. H. =U L. T. となることが分 −H. かるので,そのことを用いれば W = P L. U. −H. . 計算効率の観点からは,F を X に作用させる処理の内. X とな. 側で行われる行列をベクトルに乗じる計算や連立 1 次方程. る(あるいは W の複素共役 W = P T L−T U −T X  を計算. 式を解く前進後退代入などの計算の際には,組 X に含ま. して,それから W を求めることもできる) .これも前進消. れる列ベクトル 1 つずつに対して F を適用するのではな. 去,後退代入と置換操作により計算できる.. くて,列ベクトルをある程度の個数分まとめて処理するブ. 以上のことから,レゾルベントの作用を与える連立 1 次. ロック化を行う方がトータルの記憶の参照や移動の回数が. 方程式を LU 分解を用いた直接法で解く場合には,互いに. 減らせるという面では効率的である.ただし,途中の作業. 複素共役であるシフトに対応する連立 1 次方程式の係数行. 用記憶量が増加することや潜在的な並列度が減ることなど. 列は互いのエルミート転置であるという性質を用いて,片. の面も考慮してバランスをとる必要がある.. 方で得られた LU 分解の結果を両方で共有して再利用すれ ば,LU 分解に掛かる演算量や記憶参照量を半減できる.. 5. レイリー・リッツ法に与える基底の構成. しかし,実対称の場合とは異なり,LU 分解後の連立 1. フィルタを作用させた後のベクトルの組 Y の張る空間. 次方程式の解の組 V と W を前進消去と後退代入を用い. は, 「通過域 [a, b] あるいは遷移域に固有値のある固有ベ. て求める計算の過程と結果は共用も削減もできない.複. クトル」の張る空間を近似している.そこで,レイリー・. 素ベクトルの組 X へのフィルタ F の作用の計算の中で,. リッツ法を用いることで Y の張る空間内で元の一般固有. シフトが τ と τ である 2 つのレゾルベントの線形結合係. 値問題(GEVP)の近似固有対を求めようとする.しかし. 数をそれぞれ γ と γ とするとき,対応する 2 つのレゾル. フィルタの持つ伝達特性の性質から,固有値が阻止域ある. ベントを含む項の和によるフィルタの作用への寄与は,. いは通過域から離れた遷移域にある固有ベクトルは強く減. {γR(τ ) + γ R(τ )}X = γV + γ W となる(注意:さらに,. 衰を受けるので,Y の数値的な階数は通常は多く落ちてい. 実対称の場合とは異なり,たとえ X を実ベクトルの組と. る(しかもベクトルの個数 m を増せばそれだけより多くの. して選べた場合であっても,X にフィルタを作用させた結. 階数が落ちる傾向がある).そのような数値的な条件の悪. 果 FX は実ベクトルの組にはならない) .. い基底 Y に対してレイリー・リッツ法を直接適用すると,. 付記:レゾルベントの作用を与える連立 1 次方程式を直. 得られる結果もまた数値的に不安定になる.. 接法ではなくて反復法を用いて解く場合にも,もしも係数. そこであらかじめ一種の正則化として,Y の張る空間の. 行列の不完全 LU 分解を前処理に用いるのであれば,複素. うちから数値的な線形独立性の高い基底 Z を選び出す処理. 共役であるシフトを持つ連立 1 次方程式の係数行列が互い. を行う.いま通過域 [a, b](のある近傍)に固有値がある不. のエルミート転置であるという性質を用いることで同様. 変部分空間を I とする.そのとき,Y の張る空間の中に制. に,不完全 LU 分解を共用できるであろう.. 限された B-正規エルミート直交基底 Z を適切に構成して,. 4. フィルタの作用前と作用後のベクトルの組 いま,フィルタの通過域または遷移域に固有値がある固. c 2014 Information Processing Society of Japan . Z の張る空間が I の良い近似となるようにできる. そのようにして構成した基底 Z にレイリー・リッツ法を 適用して得られるリッツ対を元の GEVP の近似固有対とす. 60.

(5) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). る.具体的には,Z は B-正規エルミート直交(Z H BZ = I ) であるので,次数の小さいエルミート行列  A = Z H AZ を. 来は縮退しているはずの固有値が縮退がとけて分散してし. 作り,その両側ユニタリ変換による対角化(固有値分解) を A = U DU H とすると,リッツ対は ZU の列ベクトル. たかどうかを知ることは難しいので,可能ならば摂動の限. まっている可能性がある.もちろん数値的には縮退であっ 界を見積もり,値が近接している固有値は縮退している可. (リッツベクトル)とそれに対応する D の対角要素(リッ. 能性があるものとして 1 つの塊のように扱う必要がある.. ツ値)の対として得られる.またこの方法により得られる 近似固有値(リッツ値)はエルミート行列  A の固有値な. 結局実用的には,数値で求めた固有値 φ の値を gpass 未満 のものから(あるいは安全上のマージンをとるために gpass. ので必ず実数になる(近似計算で固有値を求める場合でも. よりも少しだけ小さい値,たとえば 0.95gpass などから開始. 実数が得られるようにできる) .. して)減少順に探索してゆき,φ の分布の中に比較的大き な空隙(すきま)がある場所を見つけることができたら,. 5.1 不変部分空間 I の近似基底 Z の構成 フィルタ作用素 F がレゾルベントの線形結合の場合には,. それよりも上側にある φ の値はすべて近傍として含める, などとするのが安全であろう.区間の上限である 1 につい. 元の GEVP:Av=λBv の固有対を (λ, v) とすると,ベクト. ても誤差の影響を考慮して,1 を超えている φ の値も近傍. ル v は同時にまた F の固有ベクトルにもなり,Fv=f (λ)v. 内にあるとして含めるようにする.近傍に取り入れた最小. も満たす.ただし,f (λ) は伝達関数である.その逆の命題. の φ の値が,取り入れなかった最大の φ の値から十分に離. 「F の任意の固有ベクトルは元の GEVP の固有ベクトルで. れていて,しかも 1 に比べて微小でなければ,条件の良い. ある」は成り立たないが, 「F の(任意の)不変部分空間は. 近傍の設定に成功したことになる.もしも近傍として微小. 元の GEVP の不変部分空間である」は成り立つ.. な φ の値を含めてしまうと,あたかも最初から gpass を微. この性質をうまく利用すると, 「[a, b] のある近傍」に固. 小に設定していたのと同じことになり,その後の計算での. 有値がある元の GEVP の不変部分空間を I とするとき,. 精度低下の原因になるから,近傍に取り入れる φ の値の大. 以下の性質を満たすベクトルの組 Z が構成できる:. きさには下限を設けておくのが良い.たとえば (1/2)gpass. ( 1 ) 組 Z のベクトルは Y の張る空間内にある.. よりも小さい φ は近傍には含めない,などとする(ここで. ( 2 ) Z の張る空間は,不変部分空間 I の良い近似である.. は gpass は,用いたフィルタの通過域における伝達関数の. ( 3 ) 組 Z は B-正規エルミート直交系である.. 「実際の」下限値である.なお,たとえば典型的なフィルタ. そのように構成された組 Z に対してレイリー・リッツ. の場合に非負係数 c∞ が零でないときにそれが微小である. 法を適用すると,元の GEVP の固有対であって固有値が. として零と見なして省略する近似を採用する場合には,そ. 「[a, b] のある近傍」にあるものすべての近似をいっせいに. れにともなって伝達関数の値は一様に c∞ だけ減少するか. 求めることができる. 組 Z の構成は次のようにする.まず次節の方法(5.2 節) を用いることで,ベクトル v を Y の張る空間内に制限して. ら,伝達関数の最大値である 1 や通過域における伝達率の 最小値 gpass などは,微小な正数 c∞ だけ値を減じたもの として扱うことになる) .. 得られる作用素 F の近似固有対 (φ, v) の中から,固有値 φ. 注:固有対の抜け落ちをなるべく避けるためには,区間. の値が区間 [gpass , 1](通過域における伝達関数の値域)の. [gpass , 1] をマージンをとって少し広げて F の近似固有対の. 「適切な近傍」にあるものだけをすべて求める.そのよう. 固有値の取捨選択の閾値を下げて F の近似固有対を多め. なベクトル v は B-正規エルミート直交系にとれて,それ. に取り入れると,その結果として最終的に得られる(元の. らを集めて並べたものを組 Z にすればよい.. 固有値問題の)フィルタ対角化法による近似固有対には,. 補足的考察. 通過域 [a, b] から少し離れた遷移域に固有値を持つものが. ここで区間 [gpass , 1] の「適切な近傍」とは,固有値 φ の. 含まれる可能性がある.そのような本来は不必要な固有対. 分布を(誤差の影響を考慮したうえで)はっきりと分離す. はその固有値を調べて選別して除くことができる.固有値. るような近傍のことである.数値計算で求めた固有値 φ の. の値がきわどい位置にある近似対は精度を高める改良を加. 値は近似の誤差や数値演算の丸め誤差,連立 1 次方程式の. えたあとで選別する方が安全であろう.. 解計算の不完全性などに起因する誤差によって摂動を受 ける.たとえば gpass と比べて小さい φ の値でも,それが. gpass に非常に近い値であれば安全のためには含めておく必 要がある.また,固有ベクトルの組が不変部分空間となる. 5.2 F の近似固有対の解法 ベクトル v を Y の張る空間内に制限して得られる F の 近似固有対 (φ, v) は,以下の方法で求めている.. ためには,固有値が縮退している固有ベクトルがあるとき. 固有値方程式 Fv = φ v に対して,v が Y の張る空間内. にはそれらをすべてを含むかあるいはまったく含まないか. にあるという条件の式 v=Y u を課してやり,さらに両辺. のどちらかであり,けっして部分的に含まれることがない. に左側から X H B を乗じると,次数 m の小規模な GEVP:. ようにしなければならないが,数値計算では摂動により本. α u = φ β u が得られる(左側から任意のベクトルの組で. c 2014 Information Processing Society of Japan . 61.

(6) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). はなくて X H B を乗じることによる利点は,それにより. 列の数値的な条件が良い.そこでエルミート行列用のヤコ. 得られる小規模な GEVP の係数の行列 α,β が以下に示. ビ対角化法を用いてまずエルミート行列 β をユニタリ変換. すような良い性質を持つことにある.さらに,後の 7 章. で固有値分解する.これにより得られる固有値はすべて実. で行っている分析の結果もこの手法の妥当性を補強する. 数であり,固有ベクトルを固有値の値の減少順に対応して. H. H. 材料である).ここで α=X (BF)Y ,β=X BY である.. (負の固有値のものは末尾に置いて)並べた β の固有値分. いま,F はレゾルベントの線形結合で,伝達関数 f (λ) は. 解を β=QH DQ とする(このヤコビ対角化法は,2 次のユ. 実関数としたので,BF はエルミートである.このことか. ニタリ変換を繰り返し適用して対角行列に収束させること. H. H. H. H. H. ら,まず α = X (BF)Y = X (BF) Y = X F BY =. でエルミート行列の固有値分解を求める方法であり,通常. (FX)H BY = Y H BY により行列 α は半正定値エルミー. の実対称行列に対して 2 次の直交変換(ヤコビ回転)を繰. H. H. H. H. トで,次に β = X BY = X BFX = X (BF) X = H. H. H. H. X F BX = (FX) BX = Y BX = β. H. により行列 β. り返して対角行列に収束させる方法の自然な拡張になって いる(文献 [2] および付録 A.3)).. もエルミートであることが分かる.さらに,f (λ) が正値. ユニタリ行列 Q による変換 u ≡ QH w により方程式は. 関数ならば β は正定値であること,f (λ) が非負値関数な. Gw = φ Dw となる(ただし,G ≡ Q α QH で,行列 G は. らば β は半正定値であること,f (λ) の負の値の下限値が. エルミートである) .. 微小ならば β が負の固有値を持つ場合でもその絶対値は微. ここで数値計算上の安定化のために,きわめて小さい正. 小であること,を容易に示すことができる(付録 A.1 に記. の値 thres を閾値として導入し,thres 未満の β の固有値. 述), (ただし,β が理論上は正定値や半正定値になる場合. (D の固有値)を切断する方法で「方程式の正則化」を行. でも,数値計算で求めた固有対には微小な負の固有値を持. う(たとえば thres = 100 εmac などとする.ここで εmac. つものが混入しうる.これは絶対値が小さい固有値を持つ. は使用する精度での浮動小数点数のマシンイプシロンを表. 固有対は数値誤差の影響により精度が出ないためである.. す*2 ).閾値 thres 未満である D の対角要素に対応する行. また β が不定値になる場合でも負の固有値がすべて微小. および列をそれぞれ G,D,w から省くことで,切断された  w=φ w  はエルミー   が得られる(行列 G 固有値方程式 G D. であれば,数値計算としては実質的に半正定値と同様に見 値の固有ベクトルについての関係を切断処理で棄てる.そ.  は正の対角行列である).さらに,対角スケー トであり,D  −1/2 G  −1/2 z により H ≡ D D  −1/2 とおいて  ≡D ル変換 w. れにより最終的な計算結果が大きく影響されない定式化に. 得られるエルミート行列 H の標準固有値問題:Hz = φ z. なっていればそれでよい) . たとえば典型的な 4 種類のフィルタであるバターワース. を再びエルミート行列用のヤコビ対角化法を用いて解く.  −1/2 z  := D その固有対 (φ, z) のベクトル z から逆変換 w. 型,チェビシェフ型,逆チェビシェフ型,楕円型のうちで.  を得て,切断された行の自由度に零を により対応する w. 最初の 2 種類バターワース型,チェビシェフ型は f (λ) が. 補って w を得て,さらに u := QH w とする.. なせる.そこで後で見るように,β の大きさが微小な固有. 正値の実関数なので,β は正定値エルミートになる.フィ. このようにして構成された u から v := Y u で作った v. ルタの種類が逆チェビシェフ型あるいは楕円型の場合には. をすべて集めたものが,元の一般固有値問題の不変部分空. f (λ) は負の値を阻止域においてだけとり,負の値の大きさ. 間を近似する空間の基底となる.ただし,この v はそれ. の最大値は gstop である(もしも c∞ が非零のときにフィ. だけではまだ B-正規エルミート直交系にはなっていない.. ルタ作用素 F から項 c∞ I を省略する近似を採用した場合. いま上述のように,エルミート行列 H の標準固有値問題. にはその伝達関数の負の値の大きさの最大値は 2gstop とな. をヤコビ対角化法を用いて解いたとすると,算法の性質か. る) .よって,gstop の値が微小ならば β に負の固有値があ. ら得られた固有ベクトル z は自然に正規エルミート直交系. る場合でもその大きさは微小であることがいえる.. になる.すると(β の微小な固有値を閾値で切断した影響. フィルタの性質から自然に導かれることであるが,通常. を無視すれば)u は β-正規エルミート直交系になることが. は Y の特異値には(フィルタで濾過するベクトルの個数を. 分かる.そのときには上の v(i) := Y  u(i) の代わりにスカ   ラ倍の因子を付け加えたもの v(i) := 1/ φ(i) ·Y u(i) を. 十分大きくすればそれだけ)微小な値が多くなり,係数 α, 程式 α u = φ β u は以下(5.3 節)で記述する方法により丁. 作ると,この新しい v(i) が B-正規エルミート直交系にな    ることが示せる.ただし,この規格化因子 1/ φ(i) が. 寧に解いている.. 示唆しているように,もしも通過域にある固有ベクトルの. β はともに悪条件になる.そのことを考慮して,固有値方. 伝達率に相当する φ(i) に小さい値があると数値計算では誤. 5.3 「縮小された F の固有値方程式」α u=φ βu の解法 行列 Y は一般には非常に悪条件となる可能性があるが, そのような行列 Y がそれぞれ α の表式中には 2 回,β には. 1 回,積の形で入っていることから,β の方が α よりも行. c 2014 Information Processing Society of Japan . *2. マシンイプシロン mac は 1 よりも大きい最小の浮動小数点数値 と 1 との差であり,たとえば IEEE 754 規格の 2 進 64 bit 倍精 度であれば約 2.22×10−16 である.Fortran90 言語では組み込 み関数 EPSILON を用いることで簡単に表せる.. 62.

(7) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). がある.. 「Span(Y ) 内. Av=λBv,v=Y u;. (λ, v). の GEVP」. ↓. ↓. 乱数ベクトルから計量 B の正規エルミート直交化で作っ. 「Span(Y ) 内. F v=φv,v=Y u,φ=f (λ);. (φ, v). たベクトルの組 X は,ほとんどの場合には一般的(generic). の F の SEVP」. . . な状況にあるので,N  ≤ m であれば計量 B での Y の数. α u=φ β u;. (φ, u). 値的階数はちょうど N  になる.実際には m の値をちょう. 「縮小された F. α ≡ Y H BY ,β ≡ X H BY. の GEVP」. 図 2 対応関係の図. Fig. 2 Illustration of correspondences.. ど N  にはしないで余裕を持たせて少し大きめの値にする のが良い.なぜならば乱数を用いて一般的(generic)な状 況を作り出しているが,その数値的な generic の程度(一 種の条件数)があまり良くない状況に遭遇する確率を下げ. 差が拡大して B-正規エルミート直交性が(規格化だけで. られるからである.. はなくて直交性も)崩れる可能性がある.フィルタの伝達. 「固有値がフィルタの遷移域にある固有対」が多く存在. 関数の gpass の値が 1 の付近にある場合(たとえば典型的. している場合で N  ≤ m とはならずに N ≤ m < N  であ. な 4 種類のフィルタで gpass = 0.5 と設定する)には,そ. る場合には(たとえ伝達関数が急峻に変化して遷移域が狭. のような困難は生じない.しかし,もしも伝達関数の特性. くても,その遷移域に固有値が多く分布するとそのような. 値 gpass が微小であるような特性の良くないフィルタを用. 状況が生じうる),区間 [a , b ] にある固有値 λ を重複度も. いた場合には,数値丸め誤差の拡大によって B-正規エル. 込めて数えて,それらの固有値に対応する伝達率 f (λ) を. ミート直交性が崩れる可能性があり,そのような場合には. その絶対値の減少順に並べたとき,第 m 番目よりも後の. 崩れを補正するための B-正規エルミート直交化を追加す. 伝達率が(伝達関数の最大値 1 に比べて) 「数値的に微小で. る必要が出てくる.. 零と見なせる」のであるか否かによって,近似対の精度が. 本章の 5.1 節,5.2 節,5.3 節に出てきた方程式の対応関 係を簡潔にまとめたものを図 2 に記す.. 6. 濾過するベクトルの個数 m の決め方 フィルタで濾過するベクトルの個数 m の決め方につい て,多少考察をしてみる. フィルタの通過域 [a, b] に両側の遷移域を加えて広げた . . 左右されるであろう. 以上の考察により, 「阻止域に固有値がある固有ベクト ル」がフィルタによってほぼ完全に除去されると見なせる 場合には,フィルタで濾過するベクトルの個数 m として 「阻止域を除いた区間 [a , b ] に固有値がある固有対」の個 数 N  よりも少し大きい値を採用すればよいことになる. 濾過するベクトルの個数 m の実際の決め方としては,た. 区間を [a , b ] とする.これは実軸から阻止域を除いた区. とえば以下のような方法が考えられる.. 間である.以下では,固有値が [a, b] にある固有対の個数. ( 1 ) なんらかの根拠か経験に基づいて,N  の値かその上. を N とし,固有値が [a , b ] にある固有対の個数を N  と. 限値を知っている場合には,それに基づいて個数 m. する.. を N  < m を満たせる値で少し大きめに設定する.そ. m 個の B-正規エルミート直交ベクトルの組を与えてそ. のように決めた m の値を用いて以降の計算を行う場. れをフィルタで濾過される前のベクトルの組 X とし,フィ. 合に,必要となる記憶量や演算量が過大ではないと判. ルタによって濾過された後のベクトルの組を Y とする.い. 断できたら計算を進める(過大の場合には,元の区間. まフィルタの阻止域における伝達率の大きさの上限 gstop. [a, b] を分割して扱うことを検討する).. が数値的に微小であり,零と見なして無視する近似*3 が十. ( 2 ) N  の値についての予備知識がないか,推定値に誤差. 分によく成立するならば,計量 B での Y の数値的な階数. がともなう場合には,濾過するベクトルの個数 m の値. . は N を超えない. . を(m の値にともない必要となる記憶量や演算量が過. もしも m < N であるならば,Y の m 個のベクトルの. 大にはならない範囲で)少しずつ様子を見ながら段階. 線形結合では全部で N  個ある「固有値が [a , b ] にある固. 的に増やして計算を行うことができる(増やした場合. 有ベクトル」を表すことができない.すると必要としてい. には計算を最初から完全にやり直すのではなくて,そ. る「固有値が [a, b] にある固有ベクトル」もまた,うまく近. れ以前に決めた X の列は保ちながら,新たな列を X. 似ができないかあるいは近似の精度が十分に出ない可能性. に追加して,追加した後の X 全体が B-正規エルミー. *3. いまフィルタの伝達率の最大値が 1 付近に規格化されていると すれば,阻止域における伝達率の大きさの上限値 gstop が浮動小 数点数のマシンイプシロン(IEEE 754 規格の 2 進浮動小数点数 64 bit 倍精度では約 2.22×10−16 )の程度の 10−15 や 10−14 な どであれば丸め誤差とほぼ同程度で零と見なせる.また gstop が マシンイプシロンよりもある程度大きい値(たとえば 10−12 や 10−8 など)でも,最初から計算をあたかもその程度の演算精度 で行っていると思えば,やはり零と見なせることになる.. c 2014 Information Processing Society of Japan . ト直交性を保つようにできる.そうして追加を行った. X の列だけにフィルタ F を作用させて対応する Y の 列を追加する.このように処理を追加的に行うことが 可能である) .一般的な状況にある X にフィルタを作 用させて作った Y の数値的な階数が十分に落ちてい ない場合や,得られた近似固有対の残差が十分に小さ. 63.

(8) 情報処理学会論文誌. Vol.7 No.1 57–72 (Mar. 2014). コンピューティングシステム. くない場合には,m を増やす(m を増やし続けて,つ. フィルタで濾過される m 個のベクトルの組 X は N ×m. いにある限界を超えたら,区間 [a, b] の再分割を検討. 行列であり,X は B-正規エルミート直交系になっているの. する).. で X H BX=Im である.いま考えている固有値問題の固有. ( 3 ) エルミート定値一般固有値問題 Av=λBv では,指定. ベクトル全体は全空間を張る基底なので,X を固有ベクト. した区間に固有値がある固有対の個数をシルベスター. ルで展開した係数を並べた N ×m 行列 C が一意に存在し. の慣性律を用いて求めることができる.これは 2 分法. て X=V C と書ける(本論文の他の箇所では記号 C をレゾ. などの区間縮小法で固有値を求めるのに用いられてき. ルベントに対応する連立 1 次方程式の係数行列を表すのに. た方法である.いま実数のシフト σ を与えてエルミー. すでに用いているが,適切な意味を連想させる英大文字の. ト行列 C := A − σB を作り,エルミート行列用の修正. 記号が足らないので,本章では同じ記号 C をこの展開係数. H. コレスキー法による C の分解が C =: LDL. となっ. の行列を表すためだけに用いることにする) .X が B-正規. たとき,σ よりも小さい固有値の個数は実対角行列 D. エルミート直交系であることから C H C = Im なので C は. の負の対角要素の個数に等しい.このことを用いて区. フルランクである.そのとき,Y = FX = FV C = V ΦC. 間 [a , b ) に固有値がある固有対の個数 N  は σ を b. と書ける.. . および a にした場合のそれぞれの分解で得られた対. すると, 「縮小された F の固有値問題」の係数である m. 角行列 D の負の対角要素の個数の差として求められ. 次のエルミート行列はそれぞれ α = Y H BY = C H Φ2 C ,. る.ただし,区間の両端のシフトにおける修正コレス. β = X H BY = C H ΦC と書けることが分かる.いま,フィ. キー法の分解過程において破綻が生じなかったと仮定. ルタで濾過するベクトルの個数 m は十分に大きく選ばれて. する(破綻が生じる可能性もあるが,今のように区間. いて,フィルタの伝達特性の性質から第 m 番目までよりも. 内の固有対の個数の「上限」を求めさえすればよい用. 後の伝達率の大きさは無視できるほど十分に小さいと「仮. 途であれば,破綻を生じたシフトの値を区間を広げる. 定」する(たとえば阻止域における伝達関数の大きさの上. 方向に少しだけ移動して計算をやり直して,それで破. 限値 gstop が無視できるほど小さいと見なせるならば,固. 綻が生じなければそれを真の区間の端のシフトの代替. 有値が通過域または遷移域にある固有対すべての個数より. として少し広げた区間内の固有値の個数を求めればそ. も大きい値として m を選べばよい) .この「仮定」の下で,  までで「切断する近似」 行列 Φ をその m 次の首位行列 Φ. . れは N の上限になる.ただし,数値的な安定性につ いては検討を要する) .この方法はたとえば A と B が 幅の狭い帯行列であれば高速に実行しうる(修正コレ. を導入する.それに対応して C をその最初の m 行目まで H Φ  を用いると,α ≈ C 2 C , で切断した m 次正方行列 C. スキー法の代わりに Bunch-Kaufman の算法 [1], [3] な. H Φ C  と近似できる.さらに行列 C  の条件があまり β≈C. どを用いると分解計算の数値安定性を保証できるが,. 悪くないという仮定もおくことにする(X が乱数を元にし. 対称性を保つピボット交換にともなって分解の帯幅が 次第に広がるため,計算量や記憶量の面から実施が困. て(B-正規エルミート直交化により)作ったベクトルの組  はそれほど悪条件にはならないはず であれば,普通は C. 難である可能性がある) .. である.もしもそうでなかったときには乱数列を変更して. 7. 「縮小された F の固有値問題」α u=φ βu の固有ベクトル展開を用いた分析 この章では,縮小された固有値方程式 αu=φβu の解法 について理解を深めるために,固有ベクトル展開を用いた 分析を行う. いまフィルタの伝達関数 f (λ) は実関数であると仮定す. 最初からベクトルを作り直すかあるいはベクトルの個数 m を増やせば,条件が改善できると期待できる) .すると「縮 小された F の固有値問題」α u = φ β u は,この Φ を「切  u = w とおくと,対角行列を係数 断する近似」の下では C.  w となるので,その第 2 w = φ Φ とする一般固有値問題 Φ  の第 j 対角要素の値であることが分か j 番目の固有値は Φ  は存在するが既知の量ではないから, る.ただし,C や C. る.固有ベクトルは B-正規エルミート直交系となるよう. 固有値方程式 α u = φ β u はなんらかの数値的手段によっ. に規格化されているとする.そうして元の固有値問題の固. て直接解いて近似解を求めることになる.その際,この固. 有対に番号が付けられていて,固有ベクトルの伝達率の. 有値方程式の固有値 φ の値が微小である近似対は丸め誤差. 大きさが(広義)減少順になっているとする.いま固有対. による摂動を受けるため精度良く求めることはできないで. の番号順に固有ベクトルを並べた N 次の複素行列を V と. あろう.しかし今の場合には,固有対 (φ, u) は固有値 φ が. し,固有対の番号順に固有値を並べた N 次の実対角行列. [gpass , 1](の近傍)のものだけを解けばよいのである.す. を Λ とし,固有対の番号順に固有ベクトルの伝達率を並べ. るとフィルタの設計の際に通過域での最低保証伝達率 gpass. た N 次の実対角行列を Φ とする.すると,V H BV = IN ,. として(たとえば 0.5 程度の)微小ではない値が設定され. AV = BV Λ,FV = V Φ である.ここで IN は N 次の単. ているのであれば,数値計算上の困難さは低い(もしもそ. 位行列を表す.. うではなくて,フィルタの特性が悪くて gpass の大きさが. c 2014 Information Processing Society of Japan . 64.

(9) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). 1 に比べて微小な値であれば,丸め誤差の影響で困難が生. CentOS 6.4 for x86 64.. じる可能性がある) .. 近似対の精度の評価方法. 付記:上記では,元の固有値問題 Av = λBv に対して. F の固有値方程式 Fv = φ v を経由することでフィルタの. 近似対 (λ, v) のベクトル v が B-正規化されていれば,残 √ 差ベクトル r ≡ (A − λB)v の B −1 -ノルム Δ ≡ rH B −1 r. 伝達特性を考慮に入れた近似解の構成法を分析した.比較. は近似固有値の誤差上界を与え,λ から距離 Δ 以内に真の. のために,同様の分析を F の固有値方程式を経由しないで. 固有値がある(拡張された「Wilkinson の上界」).ただし. より直接的な方法で近似解を構成する場合についても行っ. 通常は(ベクトルの近似が悪ければ特に)固有値の誤差評. てみることにする.. 価としては過大になる(なお「真の残差」が非常に小さくな. いま元の固有値問題 Av = λBv に対して,固有ベク. る場合には,残差を求める計算は相殺により数値の有効桁. トル v を Y の張る空間内に制限する近似によって得ら. を多く失わせるので,得られた残差は丸め誤差に強く影響. れる方程式 AY u = λBY u を考えて,その残差と X の. されたものになる.そのため残差を求める計算は通常より. 標準のエルミート内積をとって方程式を作ってみること. も高い精度で演算しなければ,計算で求めたノルムの値は. にする.するとそれは小さい次数 m のエルミート行列. 本来の値よりも過大になる.しかし,今回の実験ではすべ. H. H. γ = X AY ,β = X BY を係数とする一般固有値問題 γ u = λ β u になる(行列 β は前と同様にエルミートであ. ての計算は IEEE 754 の 64 bit 倍精度演算だけを用いた) . そのほかにも,逆反復による近似固有値の変動の大きさ. り,また γ = C H Λ Φ C が示せるから,行列 γ もエルミー  u = w とおくと, トである) .Φ の「切断近似」の下で,C. も固有値精度の良い推定材料になる.. Φ  w = λΦ  w にな 対角行列を係数とする一般固有値問題 Λ. 8.1 例題 1(楕円フィルタ,伝達関数の次数 2n = 24). り,その第 j 番目の固有値は一応「数学的には」元の固有値  問題の固有値 λ(j) に等しいことが分かる.ただし,C や C. 係数行列の設定. は既知の量ではないので,係数行列がともにエルミートで. h = 50 の帯行列である.行列要素の添字 p と q は 1 か. 係 数 行 列 A,B は 次 数 が N = 1000,000 で 半 帯 幅 が. ある一般固有値問題 γ u = λ β u は実際には何らかの数値. ら N の範囲であるとする.添字が |p − q| > h である帯. 的手段を用いて直接解くことになるが,その場合に φ(j) の. 外の行列要素は値が零である.対角要素の値はそれぞれ. 値が微小であるとそれだけ丸め誤差の影響によって対応す 際には,γ u = λ β u だけを単独で解いた場合には,以前の. ap,p = p − 1,bp,p = p + 1 と設定し,対角線の下側にあ √ 1 1 る帯内の要素の値をそれぞれ ap,q = p+q+1 + p−q+1 −1, √ 1 1 bp,q = (p−q+1)2 + (p+q+1)2 −1 と設定した.行列のエル. F の固有値問題を経由した場合とは異なり,得られたそれ. ミート性から対角の上側にある帯内の要素の値はそれぞれ. ぞれの近似対 (λ, u) がはたしてどの番号 j に対応するもの. 転置の位置にある要素の値の複素共役である.. るλ. (j). の近似値の精度が失われる可能性がある(しかも実. (j). なのか,また対応する φ. の値の大きさがどうなっている. フィルタの設定. のか,も分からないのである.もしも対応する φ の値が微. 用いたフィルタは楕円型で,伝達関数の次数は 2n = 24. 小ならば,それは通過域から離れたところに固有値がある. で,gpass = 0.5,μ = 1.1 と設定したところ,阻止域での伝. フィルタで強く減衰させられたベクトルであるから,その. 達率の大きさの最大値は gstop = 3.93×10−11 となった.グ. 数値的な精度はすでに失われていて誤差を多く含むことを. ラフ(図 3)は,基準座標 t による伝達関数の大きさ |g(t)|. 意味する) .すると状況次第では,そのような近似固有値 λ. の対数値を左右の対称性から右半分だけ描いたものである.. は元の固有値問題の真の固有値から任意に離れた値となっ. さらに係数 c∞ = 3.93×10−11 (n が偶数の場合には c∞ は. てしまい,そのことにより本来の問題とは無縁な固有値を. 正で今の gstop と同じ値になる)は微小であるとして,c∞. 持つ残差の大きい「偽の固有対」が生じる可能性がある.. を省略する近似を用いた.使用したフィルタの構成パラ. 8. 数値実験例 使用した計算機システム. メータの詳細については表 A·1 に記述した. 固有値の区間の指定と濾過するベクトルの個数の決定 求める固有値の区間 [a, b] は a = 1.00000 00300,b =. 実験に用いた計算機システムの仕様は以下のとおり. 1.00000 00400 と設定した.またフィルタの阻止域の端は. である.CPU:intel Core i7-3930 K(3.2 GHz,6 コア,. a = 1.00000 00295,b = 1.00000 00405 である(μ = 1.1. L3=12 MB,AVX 命令,ソケット LGA2011,設定 Turbo-. なので,区間 [a , b ] の幅は区間 [a, b] の幅の 1.1 倍になる) .. ;主記憶:計容量 64 GB Boost=OFF,HyperThread=OFF). シルベスタの慣性律に基づいた修正コレスキー分解を用. (8 Gbyte モジュールを 8 本)クアドチャネル DDR3-1600. いた事前の調査で,フィルタの通過域 [a, b] に固有値があ. (PC3-12800,CL10) ;コンパイラ:intel Fortran v14.0.0. る固有対は 80 個であること,また区間 [a , b ] に固有値を. for x86 64 (コンパイルオプション-fast -openmp);数. 持つ固有対は 89 個であることが判明した.そこで, (必要. 値と演算:IEEE 754 二進浮動小数点 64 bit 倍精度;OS:. 最小個数である)m = 89 個のベクトルを乱数により生成. c 2014 Information Processing Society of Japan . 65.

(10) 情報処理学会論文誌. 図 3. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). 例題 1:正規座標 t でのフィルタの伝達関数の大きさ |g(t)|. 図 5 例題 1:φ の値の分布(m = 89,閾値 0.45). Fig. 5 Exam1: Distribution of φ (m = 89,threshold 0.45).. (右半分). Fig. 3 Exam1: Magnitude of transfer function |g(t)| in normalized coordinate t (right-half part).. 図 6 例題 1:近似対の残差のノルム (区間 [1.00000 00300, 1.00000 00400]) 図 4. 例題 1:行列 β の固有値分布. Fig. 4 Exam1: Distribution of eigenvalues of β.. Fig. 6 Exam1: Norms of residuals of approximated pairs (whose values are in [1.00000 00300, 1.00000 00400]).. し,それを用いてフィルタ対角化を実施した.. る近似対 80 個が得られた.. 実行結果. 近似対の改良. フィルタ作用素 F の近似固有値 φ は,縮小された固有. フィルタ対角化法で得られた近似対をレイリー商逆反復. 値問題 αu = φ βu を,β の固有値の切断による正則化の. で 2 回の改良を行った.グラフ(図 6)の中で,折れ線 IT0. ための閾値を thres = 10−12 として解いた.計算で求めた. はフィルタ対角化法の各近似対の固有値に対して残差のノ. β の固有値を値の減少順にプロットした対数グラフを図 4. ルムを対数でプロットしたものであり,折れ線 IT1,IT2. に示す(最後の 2 個の負の固有値は青色の点でプロットし. はそれぞれ逆反復法で 1 回,2 回改良した後の近似対の残. た).切断による正則化で自由度は m = 89 から 86 に減. 差のノルムの対数プロットである.残差のノルムは固有値. 少した.得られた 86 個の固有値 φ を値の減少順に(番号. の誤差の(通常は過大な)限界を与える.. を 1 番から 86 番まで付けて)並べると φ80 = 0.50005 53,. グラフから,フィルタ対角化法の近似対の固有値の誤差. φ81 = 0.38662 97,φ82 = 0.00350 89,などとなった.それ. は,1.5×10−13 以下で,固有値の相対精度は約 13 桁程度. ら 86 個の固有値の分布に対して,まず gpass = 0.5 未満で. であることが分かる.逆反復法 1 回目の改良では約 1 桁程. かつ分布中の広い間隙内にある値 0.45 を上下を分離する閾. 度改善されたが,倍精度計算でもあり,2 回目の改良では. 値に設定した.するとこの閾値の上側には必要な固有値の. あまり改良されなかった.. 個数にちょうど等しい 80 個の φ の近似値があった(図 5) .. 経過時間. これから(既述の方法により)不変部分空間を近似する 80. 計算は OpenMP で並列化し,CPU 内の 6 コアを用いて. 個の B-正規直交基底を構成し,それにレイリー・リッツ法. 6 スレッドで実行した(フィルタ対角化の中でフィルタの. を適用すると,元の一般固有値問題の固有値が [a, b] にあ. 適用に出てくるレゾルベントの作用は連立 1 次方程式を直. c 2014 Information Processing Society of Japan . 66.

(11) 情報処理学会論文誌. 図 7. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). 例題 2:正規座標 t でのフィルタの伝達関数の大きさ |g(t)|. 図 8. 例題 2:行列 β の固有値分布. Fig. 8 Exam2: Distribution of eigenvalues of β.. (右半分). Fig. 7 Exam2: Magnitude of transfer function |g(t)| in normalized coordinate t (right-half part).. 接法により行ピボット選択を行う LU 分解を用いて解いた が,その際にはシフトが複素共役のレゾルベントの組に対 しては LU 分解の結果を共用する方法を用いて計算してい る) .フィルタ対角化法で近似対 80 個を求めるまでの経過 時間は全部で 117.83 秒であり,そのうちで 12 次のフィル タを m = 89 個のベクトルに作用させる部分の経過時間は. 93.86 秒であった.また,フィルタ対角化法で得た近似対 80 個に対して,逆反復法による改良を 2 回繰り返すのに要 した経過時間は全部で 220.40 秒であった.. 8.2 例題 2(楕円フィルタ,伝達関数の次数 2n = 12). 図 9 例題 2:φ の値の分布(m = 301,閾値 0.45). Fig. 9 Exam2: Distribution of φ (m = 301, threshold 0.45).. 係数行列の設定 例題 2 の行列 A,B は次数が N = 1000, 000,半帯幅が. 個であることが分かった.そこで, (必要最小個数である). h = 50 で,それぞれ例題 1 とまったく同一のものである.. m = 301 個のベクトルを乱数により生成し,それを用いて. フィルタの設定. フィルタ対角化を実施してみた.. フィルタの種類は楕円型で,今度は伝達関数の次数を. 実行結果. 2n = 12 と例題 2 よりも小さくして,gpass = 0.5 で μ = 3.0. フィルタ作用素 F の固有値問題 αu = φ βu は β の正. と設定した.そのとき阻止域での伝達関数の大きさの最大値. 則化に用いる固有値の切断閾値を thres = 10−12 と設定し. は gstop = 2.55×10. −12. となった.グラフ(図 7)は,基準座. て解いた.計算で求めた β の固有値を減少順にプロット. 標 t による伝達関数の大きさ |g(t)| の対数値を左右の対称性. した対数グラフを図 8 に示す(257 番以降の負の値は青. −12. 色の点でプロットした).切断による正則化後に得られた. から右半分だけ描いたものである.係数 c∞ = 2.55×10. は微小なので省略する近似を用いた.用いたフィルタの構. 228 個の固有値 φ を減少順に並べると,φ100 = 0.50000 47,. 成パラメータの詳細については表 A·2 に記述した.. φ101 = 0.36785 11 などから,gpass = 0.5 未満の φ の分布. 固有値の区間の指定と濾過するベクトルの個数の決定. 中にある広い間隙中の値である 0.45 を φ の値の分布の上. 求めるべき固有値の範囲の区間 [a, b] は例題 1 の場合. 下を分離する閾値に設定した.するとこの閾値の上側には. と異なり,a = 1.00000 01000,b = 1.00000 01100 とし. 必要な固有値の個数にちょうど等しい 100 個の φ の近似値. . た.フィルタの阻止域の端点は a = 1.00000 00850 と . . . b = 1.00000 01250 になる(μ = 3.0 なので,区間 [a , b ] の 幅は区間 [a, b] の幅の 3.0 倍になる) . シルベスタの慣性律に基づいた修正コレスキー分解を用 いた事前の調査により,区間 [a, b] に固有値を持つ固有対 . . は 100 個であり,区間 [a , b ] に固有値を持つ固有対は 301. c 2014 Information Processing Society of Japan . があった(図 9). 不変部分空間を近似する B-正規直交ベクトルとして 100 個の基底を構成し,それにレイリー・リッツ法を適用する ことで,フィルタ対角化法による近似対が 100 個得られた. 例題 1 の場合と同様に,グラフ(図 10)には,フィルタ 対角化法で得られた近似対の残差のノルム(IT0) ,および. 67.

(12) 情報処理学会論文誌. コンピューティングシステム. Vol.7 No.1 57–72 (Mar. 2014). スレッド数を 6 から 3 にしたことにより演算に利用される. CPU コアの数が半分に減少したので,フィルタ対角化法で 求めた 100 個の近似対を逆反復で 2 回改良するのにかかっ た経過時間は,スレッド数が 6 の場合に比べて約 2 倍に増 えて 498.93 秒になった.. 8.3 切断正則化のための閾値に対する考察 5.3 節ですでに述べたように,不変部分空間を近似する 空間の B-正規直交基底を構成する際には,小規模の一般 固有値問題 αu = φ βu を係数行列が悪条件なのでまず正 則化を施してから解いている.正則化は閾値を設定して β 図 10 例題 2: 近似対の残差のノルム (区間 [1.00000 00300, 1.00000 00400]). Fig. 10 Exam2: Norms of residuals of approximated pairs (whose values are in [1.00000 00300, 1.00000 00400]).. の固有値を切断することで行う.最終的に得られる不変部 分空間の基底となるベクトルの個数は,この正則化で切断 されずに残った β の固有値の個数以下になる. いま k を元の GEVP の必要な固有値の固有対の個数と するときに,切断後の自由度が k 以上となるためには,閾. 逆反復法を 1 回,2 回適用して得られた近似対の残差のノ. 値 thres は β の最大側から第 k 番目の固有値よりも小さく. ルムをそれぞれ IT1,IT2 として対数でプロットした.グ. なければならず,また相対的に十分に小さい値にするべき. ラフから,フィルタ対角化法の近似対の精度はすでに 13∼. である.ただし閾値の絶対的な値を,マシン・イプシロン. 14 桁程度あるので,倍精度計算であるためか,逆反復法に. 程度以下のような小さい値にすることは適切ではないであ. よる精度の改良があまりないことが分かる.. ろう(閾値を下げれば正則化後の自由度が増えて近似に参. 経過時間. 加する空間の自由度が増やせるが,過度に下げると正則化. この例題 2 では,μ = 3.0 とフィルタの遷移域が広くて,. を行う意味が失われて,β の固有値の分布にもよるが,丸め. そこに固有値が入る固有対が多くあり,そのためフィルタ. 誤差の影響が拡大されて数値不安定となる可能性がある) .. で濾過するベクトルの個数 m が多くなった.行列の次数. ベクトルの組 X からフィルタ操作によりベクトルの組 Y. 6. が N = 10 ということもあって,現状のプログラムのま. を作る過程で「不要な固有値の固有ベクトル」を数値の相. までは 6 組の共役なシフトに対するレゾルベントの作用. 殺で減衰させる際に生じる丸め誤差を Y は含む.いま仮に. を OpenMP による 6 スレッドで同時並行に計算すると,. 相殺による減衰が起きなかったとした場合の β = X H BY. 実メモリの容量が不足してしまった(その原因は,6 対の. の潜在的な固有値の上限は 1 の程度なので,相殺の丸め誤. 相異なるシフトの共役対についてそれぞれ同時並行的に,. 差による摂動の影響の大きさはマシン・イプシロン(の定. m = 301 個の異なる右辺ベクトルに対応する連立 1 次方程. 数倍)程度であろうかと思われる.であれば,たとえ実際. 式を LU 分解後の前進消去後退代入などの際に m 組をま. の β の固有値がどれも(X が「必要な固有値の固有ベクト. とめていっせいに解いたので一時的な作業用記憶が多くな. ル」を含む割合が少ないために)1 よりもかなり小さい値. りすぎたためである).そのために 301 個のベクトルを濾. になっていたとしても,単純にそれに比例させて切断の閾. 過するのにかかった経過時間が 12,043 秒となり,フィルタ. 値を下げることは適切ではないと考えられる.. 対角化全体の経過時間は 13,077 秒となった.フィルタ対. 以上のことから,閾値 thres の設定の方針として,. 角化法で得られた 100 個の近似対を逆反復を用いて 2 回改. ( 1 ) β の最大側から第 k 番目の固有値よりも十分に小さい,. 良するのにかかった経過時間は 266.52 秒であった(この場. ( 2 ) マシン・イプシロンよりはある程度大きい,. 合にも,シフトが複素共役のレゾルベントの組に対しては. の両方の条件を満たすような値にするのが良いであろうと. LU 分解の結果を共用する方法を用いて計算をしている).. いうことになる(あるいは,k がいくつであるのかを事前. 使用スレッド数を減らせば,プログラムをまったく書き. に知らない場合には,簡単に条件 ( 2 ) を満たす値にしてお. 直さずに使用する実メモリの容量を減らすことができる.. いて,その値により条件 ( 1 ) も満たされることを期待する. そこでプログラム外部の環境変数 OMP NUM THREADS を設. というやり方をとる) .. 定することで許容する最大スレッド数を 6 から 3 に減ら. 表 1 と表 2 はそれぞれ例題 1 と例題 2 について,表の. して,プログラムを変更せずに実行してみると,フィルタ. 各列に行列 β の固有値の切断の閾値 thres を 10−6 ,10−8 ,. 対角化法全体の経過時間が大幅に短縮されて 537.46 秒と. 10−10 ,10−12 ,10−14 と変えて設定して計算を実行した結果. なった.これは実メモリの容量不足が解消されて仮想記憶. を掲げたものである(計算には IEEE 754 の 64 bit 倍精度. 機構によるオーバヘッドが減少したからである.一方で,. 浮動小数点数を用いている) .例題 1 の場合は,必要な固有. c 2014 Information Processing Society of Japan . 68.

図 1 フィルタの伝達関数の絶対値 | g(t) | の概形図 Fig. 1 Conceptual graph of the filter’s transfer function | g(t) | .
Fig. 2 Illustration of correspondences.
図 6 例題 1 :近似対の残差のノルム
図 9 例題 2 : φ の値の分布( m = 301 ,閾値 0.45 ) Fig. 9 Exam2: Distribution of φ (m = 301, threshold 0.45).
+4

参照

関連したドキュメント

Abstract. Recently, the Riemann problem in the interior domain of a smooth Jordan curve was solved by transforming its boundary condition to a Fredholm integral equation of the

When the matrices A and B are of small to moderate sizes, the Tikhonov minimization problem (1.4) is typically simplified by first computing the Generalized Singular Value

We study the classical invariant theory of the B´ ezoutiant R(A, B) of a pair of binary forms A, B.. We also describe a ‘generic reduc- tion formula’ which recovers B from R(A, B)

Due to Kondratiev [12], one of the appropriate functional spaces for the boundary value problems of the type (1.4) are the weighted Sobolev space V β l,2.. Such spaces can be defined

But in fact we can very quickly bound the axial elbows by the simple center-line method and so, in the vanilla algorithm, we will work only with upper bounds on the axial elbows..

We also discover continued fractions whose approximants generate every term in diagonals and branches of the Stern–Brocot tree.. Introduction

Spectral properties of sign symmetric matrices are studied.A criterion for sign symmetry of shifted basic circulant permutation matrices is proven, and is then used to answer

Note that any type II component of the singular locus is associated to exactly one reduced critical 3412 embedding.. However, if both A and B are nonempty, then we do not have a type