学会会場の CAAS とその付近
１. 学会の概要
本会議は，1996年から隔年で行われている，固有値計算に関する国際会議である。参加者は30〜50名と比較的小規模であるが，固有値計算に関する第一線の研究者が多数参加し，レベルの高い講演が多い。また，パラレルセッションがなく，1講演の時間が25分〜50分と長いこと，コーヒーブレイク，昼食，エクスカージョンなどの時間がたっぷり取られていることなどから，参加者どうしの議論が十分にできるのが特色である。
今回は，クロアチアのドブロブニクで開催された。ドブロブニクは，「アドリア海の真珠」として世界的に知られ，中世の街並みが残る美しい町である。また，クロアチアは，固有値計算のためのヤコビ法の理論的解析に関しては，世界をリードする実績がある。我々は最近，ヤコビ法に基づく超並列向けのソルバを開発しており，それに関連して収束性についての解析も少し行っていたことから，今回，参加することにした。
２. 目立ったテーマ
多項式固有値問題と有理式固有値問題に関する講演が全部で6件あり，関心の高さが伺えた。これらの固有値問題は，線形化と呼ばれる手法によって一般固有値問題に帰着できるが，新しい線形化の手法や，線形化の手法によって条件数がどう変わるかなどのテーマが議論されていた。
ほとんどの発表は固有値問題をテーマとしていたが，Riccati方程式やSylvester方程式の数値解法，行列関数の計算法，線形計算におけるランダマイズドアルゴリズムなど，関連するテーマの講演もあった。
３. 興味深かった講演
特に興味深かった発表と自分の発表について，以下にまとめる。
(1) A Restarted Induced Dimension Reduction method to approximate eigenpairs of large unsymmetric matrices
R. Astudillo and M. B. van Gijzen (Delft University of Technology)
IDR(s)法に基づく固有値解法の提案。大規模疎行列の固有値計算でよく使われるArnoldi法では，クリロフ部分空間の正規直交基底を作っていくことで，行列AのArnoldi分解 AVm = VmHm + vm+1hm+1,memT を計算し，Hmの固有値を求めることで A の固有値の近似値を求める。ここで，数値的安定性を別にすれば，基底は別に直交していなくても，同様の分解を作ることができ，Hmから A の近似固有値を求めることができる。
IDR(s)法の場合，基底 wiは，\({\bf w}_{i+1}=(A-\mu_{j+1}I)\left({\bf w}_i-\sum_{l=1}^s c_l {\bf w}_{i-l}\right)\)という式で計算される。ただし，μj+1は任意のスカラーであり，clは，右側の括弧の中のベクトルが，あるs次元線形部分空間Pに直交するように決められる。この式より，wi+1 は\({\bf w}_{i-s}, \ldots, {\bf w}_i, A{\bf w}_{i-s}, \ldots, A{\bf w}_i\) の線形結合で書けることがわかるが，帰納法より，\(A{\bf w}_{i-s}, \ldots, A{\bf w}_{i-1}\)は {wl} の線形結合で書けることがわかるので，結局，wi+1 は Awi と，wl（l≦i）の線形結合で書けることになる。これより，IDR(s)法の基底ベクトルについても，Arnoldi分解と同じ形の式 AWm = WmHm + wm+1hm+1,memT が成り立つことがわかる。これをIDR分解と呼ぶ。このヘッセンベルグ行列の固有値を求めることで，A の近似固有値が求められる。
Arnoldi分解との違いとして，IDR分解におけるヘッセンベルグ行列では，上三角部分が帯行列になる。これは，漸化式にwi-s 以降のベクトルしか現れないことに対応する。また，数値的安定性のため，基底ベクトルを直前のs本のベクトルと正規直交化することが推奨される。これは，ヘッセンベルグ行列の帯構造を壊さない。また，直交化の演算量は，m に比例する程度で済む。これがIDR分解の利点である。
欲しい固有値だけを求めるには，2つの方法がある。1つは，SorensenらによってArnoldi法のために提案された陰的リスタートを使う方法である。具体的には，Hm に対して，不要な固有値をシフトとするQR反復を適用することで，不要な基底ベクトルを分離する。もう1つは，任意パラメータ μj を調節することで，IDR法の残差多項式（のうち，パラメータに依存する部分）の絶対値が，不要な固有値を含む領域で小さくなるように調節することである。これには，チェビシェフ多項式が利用できる。
感想： 本発表はポスター発表であるが，IDR(s)法に基づく固有値計算の原理がわかりやすく説明され，興味深かった。IDR(s)法に基づく固有値計算は，Gutknecht & Zemke も行っているが，そちらでも，基底ベクトルの直交化は直前のs本に限っているようである。なぜそれでうまく行くのか，知りたいと思う。
(2) The Fiedler matrix and the LR algorithm
Carla Ferreira (Universidade do Minho, Portugal) and Beresford N. Parlett
代数方程式の解法として，コンパニオン行列の固有値問題に基づく解法がある。これに対して，Fiedler [2-1] は，n次のコンパニオン行列 Cnが Cn=An-1・・・A1A0 とn個の行列の積に分解できることを示した。ここで \(A_k = I_{n-k-1}\oplus \begin{pmatrix}-a_k & 1 \\ 1 & 0\end{pmatrix}\oplus I_{k-1}\)であり，akは，代数方程式のk次の項の係数（ただし最高次の係数は1）である。Fiedlerはさらに，A0, A1, ..., An-1の並び順を任意に変えてできるn!個の行列はすべて相似であることを示した。したがって，代数方程式の解を求めるには，これらのn!個の固有値問題のどれを解いてもよい。
特に，偶数番目のAiだけを掛け合わせた行列をAE，奇数番目のAiだけを掛け合わせた行列をAOとおくと，これらはそれぞれ2×2の対角ブロックを持つブロック対角行列であり，それらの積AEAOは5重対角行列となることが示せる。これをFiedlerコンパニオン行列と呼ぶ。
そこで，この行列にLR法を適用することで，代数方程式の根を求める方法を提案した。LR法は帯幅を保存することから，この方法は1反復あたりの計算量がO(n)となる。安定性の理論的証明はないものの，FrancisシフトとWilkinsonシフトを使って計算を行ったところ，多くの代数方程式に対し，安定に計算を行うことができた。不安定になるのは，LR分解において要素の絶対値が大きくなる場合であるが，その場合はシフトを取り直して計算することで対処する。
さらに，Fiedlerコンパニオン行列は，逆行列もまた5重対角になるという著しい性質を持つ [2-2]。このことを利用すると，反復の途中で，逆行列に対する反復に切り替えるといったことが自由にできる。これは，安定性の向上や収束の加速に役立つのではないかと期待される。
感想： コンパニオン行列に対するO(n2)の固有値解法は，最近，ホットな研究テーマであり，[2-3] など，いくつもの解法が提案されている。その中で，本発表の解法は，(i) コンパニオン行列の分解，(ii) 因子の並べ替え，(iii) その結果得られる5重対角行列に対するLR法の適用，という3つの要素からなっており，非常にわかりやすい点に特徴がある。ただし，最近の研究では，代数方程式の係数の絶対値が大きい場合，Fiedlerコンパニオン行列の後退誤差は，普通のコンパニオン行列の後退誤差に比べて大きくなるという結果もある [2-4]。今後，種々の解法間の比較が進むことが期待される。
コーヒーブレイクの際に発表者の Carla Ferreira と話をしたところ，精度と安定性の向上のために，5重対角行列用のLR法をdqds型に書き換えようとしたが，うまく行っていないとのことであった。そこで，5重対角は難しいが，下2重対角が1個と上2重対角が複数個の積の行列に対してならば，dqds法を拡張でき，各2重対角行列の要素が正ならば，小さい相対誤差で全固有値を計算できると話したところ，興味を持ってもらえた。論文を送ることにした。
参考文献
[2-1] M. Fiedler: "A note on companion matrices". Linear Algebra and Its Applications, Vol. 372, pp. 325-331 (2003).
[2-2] G. Strang: "Fast transforms: Banded matrices with banded inverses", Proc. of National Academy of Sciences, Vol. 107, No. 28, pp. 12413-12416 (2010).
[2-3] J. L. Aurenz, R. Vandebril and D. S. Watkins: "Fast computation of the zeros of a polynomial via factorization of the companion matrix", SIAM Journal on Scientific Computing, Vol. 35, No. 1, pp. A255-A269 (2013).
[2-4] F. D. Teran, F. M. Dopico and J. Perez: "Backward stability of polynomial root-finding using Fiedler companion matrices", submitted.
(3) On the convergence of the cyclic and quasi-cyclic block Jacobi methods
Erna Begović and Vjeran Hari (University of Zagreb, Croatia)
固有値計算のためのヤコビ法は，3重対角化に基づく方法に比べて計算量は多いものの，並列性が高いこと，ある条件下では小さい相対誤差で固有値を求められるなどの優れた特徴を持つ。ヤコビ法の拡張として，三角行列の特異値分解を計算するための Kogbetliantz 法，一般固有値問題を解くための Falk-Lamgemeyer 法などがある。一方，高性能計算機上でキャッシュメモリを有効に利用して演算を高速化するための技法としてブロック化があり，ヤコビ法系統の解法に対しても，ブロック化アルゴリズムが重要となっている。
本研究では，巡回型のブロックヤコビ法およびその拡張について，大域的収束性を議論するための一般的な理論を構築した。そのための主な道具は，Henrici と Zimmermann によって導入された巡回ヤコビ演算子の理論の，ブロックへの一般化である。巡回型ヤコビ演算子とは，あるステップでの行列 A の非対角要素を並べてできるベクトルを a，1巡回後の行列 A' に対するベクトルを a' とするとき，a から a' への写像のことである。この写像は a に依存するので非線形であるが，行列 F(a) を用いて a'=F(a)a と書ける。したがって，行列 F(a) の2ノルムが，a によらない1未満の上界を持つならば，ヤコビ法の大域的収束性が示せる。
本発表では，（一般化された）巡回型ブロックヤコビ法 A(k+1)=Fk-1A(k)Fk が次の3つの条件：
を満たせば，大域的収束性を持つことを示した。なお，通常のヤコビ法では，Ukはもともとユニタリ行列なので，(ii)は自明に成り立つ。また，(iii)は UBC（Uniformly Bounded Cosine）条件と呼ばれ，[3-2] で導入された条件である。これは，（ブロック版でない）巡回型ヤコビ法の大域的収束性の解析において，Henrici らによって導入された，Givens 回転のコサインが0より大きい下界を持つという条件に対応している。
また，非対角ブロックの個数が（上三角部分だけを数えて）Nの場合，ピボット戦略はN!通り存在するが，本研究ではそれらの間に，行・列ブロックの番号の置換により同値なもの，消去順序の反転により同値なもの，弱同値なものなど，いくつかの同値関係を定義し，ある戦略に対して大域的収束性が証明されれば，それと同値な戦略についても大域的収束性が成り立つことを示した。特に，これまで知られていた行（列）方向の巡回（左の列から右の列に，列内では上から下に）だけでなく，「左の列から右の列に進むが，各列内での消去順序は任意」というピボット戦略でも，大域的収束性が成り立つことが示された。
感想： ブロックヤコビ法の大域的収束性については，長らく理論的証明がなかったが，最近，クロアチアの研究者によって急速に研究が進んでいる。最初の理論的成果は Drmać [3-2] によるもので，行（列）方向巡回型のブロックヤコビ法において，変換に用いるユニタリ行列が UCB 条件を満たせば大域的収束性が保証されること，また，UCB 条件を満たすユニタリ行列は，2×2ブロックの固有値問題を解き，その固有ベクトル行列に適当な列の置換を行えば実際に構成できることを示した。本研究では，この結果が，拡張されたブロックヤコビ法，及び，より広いクラスのピボット戦略について成り立つことが示された。論文 [3-1] はかなり長く，込み入っているが，ヤコビ法の研究をするには読む必要があると思われる。
参考文献
[3-1] V. Hari: "Convergence to zero of off-diagonal part in block Jacobi-type methods", submitted.
[3-2] Z. Drmać: "A Global Convergence Proof for Cyclic Jacobi Methods with Block Rotations", SIAM Journal on Matrix Analysis and Applications, Vol. 31, No. 3, pp. 1329-1350 (2009).
(4) Accurate eigenvalue decomposition of arrowhead matrices, rank-one modification of diagonal matrices and applications
Nevena J. Stor, Ivan Slapničar (University of Split, Croatia) and Jesse Barlow
A=D+uuT という，対角+ランク1の形をした行列の固有値問題は，分割統治法などで重要である。本研究では，この問題に対して，前進安定な解法，すなわち，計算された固有値・固有ベクトルが真の固有値，固有ベクトルに対して小さい相対誤差を持つ解法を提案した。また，対角+ランク1行列とほぼ等価な arrowhead 行列に対する解法も提案した。
解法の基本的なアイディアは，A に対してシフト付き逆反復法を適用することである。Sherman-Morrison 公式より，対角+ランク1行列（をシフトした行列）の逆行列は，再び対角+ランク1行列になる。しかも，実際に式を書いてみればわかる通り，逆行列の要素は，1つの要素を除き，スカラーの掛け算と割り算のみで計算できるので，高い相対精度での計算が可能である。残りの1つの要素については，内積を計算する必要があるが，この部分だけは4倍精度で計算して精度を保つ。なお，シフトした結果，D のある対角要素が 0 になることがあるが，その場合にも，類似の公式が存在する。固有値 λi が求まったら，固有ベクトル x を，標準的な分割統治法と同様に，公式 x = (D-λiI)-1u により計算する。
本アルゴリズムは，現在使われている Gu & Eisenstat による対角+ランク1の固有値問題の解法と同様，固有値ごとの並列性を持つ。しかも，Gu & Eisenstat の方法が後退安定性しか持たないのに対し，本アルゴリズムは前進安定なので，これまでの解法を置き換えることが期待できる。
感想： 本発表はポスター発表である。比較的単純なアイディアで前進安定性を実現できる新しい解法ということで，大変興味深かった。縮重に非常に近い場合まで含めて数値実験を行い，理論通りの結果が得られているということである。密行列の固有値問題に適用する場合は，3重対角化が前進安定性を崩してしまうので，本解法の特色がそのまま生かせないのは残念である。しかし，最近エクサスケール向けの固有値解法 Eigen Exa [4-2] などで使われている5重（一般には2m+1重）対角行列の固有値解法では，対角+ランク1行列の固有値問題を，ランク1の摂動を2回（m回）に渡って加えることにより解く。その際，精度を保つためには，対角+ランク1行列の固有値問題を高精度に解く必要があるので，本解法を使うメリットがあるかもしれない。
後で気が付いたが，前進安定性が意味を持つためには，対角+ランク1という行列の表現法が，D および u に微小な相対的摂動が加わったとき，固有値・固有ベクトルも微小な相対的摂動しか受けないという性質を持つ（すなわち，Relatively Robust Representation になっている）必要があると思う。このことは証明されているのだろうか。
参考文献
[4-1] N. J. Stor, I. Slapničar and J. Barlow: "Accurate eigenvalue decomposition of arrowhead matrices and applications", to appear in Linear Algebra and Its Applications.
[4-2] EigenExa, http://www.aics.riken.jp/labs/lpnctrt/EigenExa.html.
(5) Solving eigenvalue problems with a googol of unknowns
Daniel Kressner (École polytechnique fédérale de Lausanne, Switzerland)
量子力学における多電子問題では，高次元空間における偏微分方程式の固有値問題を解く必要がある。また，確率的ネットワークの解析でも，高次元の固有値問題が生じる。これらの固有値問題はテンソル構造を持っており，効率的な求解のためには，その構造を活用した数値解法の開発が必要である。本発表では，低ランク近似に基づくテンソルの表現法と，それを利用した固有値問題の解法が紹介された。
A を各方向のサイズが n の d 次元テンソルとし，その表現法，及び低ランク近似の方法を考える。よく知られた手法として，CP（Candecomp/Parafac）分解とTucker分解（Higher Order SVD）がある。前者はランク1のテンソルの和としてテンソルを表現する方法であり，後者は，ランク d のコアテンソルを考え，その各添字に対し，2次元の直交行列との縮約を行うことで，元の d 次元テンソルを表現する方法である。前者は，ランクを r で打ち切ったとき，表現に用いるパラメータが dnr 個と少なくて済むが，頑健な分解アルゴリズムが知られていない。一方，後者は特異値分解を用いて安定に計算できるが，コアテンソルの各方向のサイズを r としたとき，必要なパラメータ数が dnr + rd 個と多い。そこで，両者の特長を併せ持った表現，すなわち，パラメータ数が d の指数オーダーにならず，かつ，特異値分解を用いて安定に計算できる表現法が求められる。
これに応えるのが，2009年頃より研究が盛んになってきた Tensor Train 形式である。これは，物性物理で用いられる Matrix Product States と同じものである。Tensor Train による A の表現は，\(A(i_1, \ldots, i_d) = \sum_{\alpha_1, \ldots, \alpha_{d-1}} G_1(i_1,\alpha_1)G_2(\alpha_1,i_2,\alpha_2)\cdots G_d(\alpha_{d-1},i_d)\) と書かれる。すなわち，2個の2階のテンソル \(G_1, G_d\)，d-2個の3階のテンソル\(G_2, \ldots, G_{d-1}\)で A を表す。必要なパラメータ数は dnr2 個と少ない。なお，右辺の和を陽的に計算してしまうと，r に関して指数オーダーの計算量と格納領域が必要となるが，この和は計算せずに，\(G_1, G_2, \ldots, G_d\) のみをデータとして持つ。この形式では，テンソルに対する様々な基本的演算が，d に関して線形，r に関して多項式オーダーの計算量で行える。また，Tensor Train を拡張した形式として，Tensor Network があり，これは，CP分解，Tucker分解，Tensor Train，Hierarchical Tucker など，これまでテンソルの表現に用いられてきたほとんどの形式をその特別な場合として含む。
テンソルの低ランク表現を用いると，高次元の固有値問題を数値的に解く際に，近似解の探索空間を低ランクのテンソルの空間に制限して解くことができる。これにより，計算量と記憶容量が大きく削減できる。低ランクの空間で探索を行うための手法としては，Riemannian optimization がある。また，Tensor Train 形式において，各 Gi に対して順番に最適化を行い，これを反復する Alternating Least Squares という手法も有効である。
参考文献
[5-1] L. Grasedyck, D. Kressner and C. Tobbler: "A literature survey of low-rank tensor approximation techniques" , arXiv:1302.7121 (2013).
[5-2] W. Hackbusch, "Tensor Spaces and Numerical Tensor Calculus", Springer, 2012.
(6) On error resilience of a complex moment-based eigensolver
Tetsuya Sakurai, Yasunori Futamura and Akira Imakura (University of Tsukuba)
Rayleigh-Ritz 型の Sakurai-Sugiura 法（CIRR法）の持つ耐故障性についての解析。CIRR法では，n次一般固有値問題 Ax=λBx の，複素平面上の円Γ内にある固有値を求めるために，n×L（Lは適当な自然数）の乱数行列Vを用意し，\(S_k \equiv \frac{1}{2\pi i}\oint_{\Gamma}z^k (zB-A)^{-1}BV\)，\(S\equiv[S_0, S_1, \ldots, S_{M-1}]\)（Mは適当な自然数）として，n×ML 行列 S に対して Rayleigh-Ritz 法を適用する。このとき，ML がΓ内の固有値の個数以上ならば，Γ内の固有値をすべて求めることができる。さらに，円Γの中心に近い順に，全固有値に番号を付けるとき，i 番目の固有対の残差は，\(|f(λ_{LM+1})/f(λ_i)|\)に比例することが示される [6-1]。ただし，f はこの複素周回積分（を数値積分で計算したもの）に対応するフィルタ関数である。
さて，ある積分点 zj での計算に，浮動小数点数のビットフリップなどによる大きな誤差が混入し，そこでの計算結果が \((z_j B-A)^{-1}BV + E_j\)となったとする。ここで，Ej は誤差行列で，そのランクは K≦L であるとする。また，Im(Ej) の基底を \({\bf e}_1, {\bf e}_2, \ldots, {\bf e}_K\) とする。このとき，Ej の影響で S にも誤差 E が混入するが，容易にわかるように，E の各列は，{ei} の線形結合になる。したがって，S の第 i 列を si とすると，S+E の第 i 列 \(\tilde{\bf s}_i\) は\(\tilde{\bf s}_i={\bf s}_i + \sum_{k=1}^K c_{ik}{\bf e}_k\) と書ける。これより，\(\{\tilde{\bf s}_i\}\) から {ei} を消去し，Im(S+E) の基底ベクトルで {ei} を含まないものを ML-K 本作ることができる。これらが張る ML-K 次元空間で Rayleigh-Ritz 法を行うことができたとすると，誤差 E の影響を受けない固有対が ML-K 個求まる。実際には，（誤差 ek の係数 cik はわからないから）これに残りの K 本の基底（{ei} を含むもの）を加えた空間 Im(S+E) で Rayleigh-Ritz 法を行うことになるが，これは，拡大した空間での Rayleigh-Ritz 法であり，固有対の精度が良くなることはあっても，悪くなることはない。以上より，大きな誤差 Ej の混入があっても，ML-K 個の固有対は正しく求まることがわかる。誤差解析によれば，このとき，i 番目の固有対の残差は，\(|f(λ_{LM-K+1})/f(λ_i)|\)に比例することが示される [6-1]。
感想： エクサフロップスマシンでは，耐故障性が重要な課題になると言われるが，本研究はアルゴリズムレベルでの耐故障性を実現しようとする試みであり，非常に興味深かった。同じくアルゴリズムレベルでの耐故障性を図った疎行列固有値解法として，Multiple Explicitly Restarted Arnoldi Method（MERAM）がある。これは，初期ベクトルが異なる複数の Arnoldi 法を同時に走らせてできる複数の Krylov 部分空間の和空間に対して Rayreigh-Ritz 法を適用する方法である。誤差のために1つの Arnoldi 法が全く変な結果を返しても，それは部分空間が無駄に拡張されるだけで，残りの Arnoldi 法が生成した空間で，十分精度の良い解が得られるという原理である。障害による結果を，無駄な空間拡張として解釈し，空間の拡張により固有対の精度が良くなることはあっても悪くなることはないということを利用して，耐故障性を実現している点は，本手法と共通する。ただし，MERAM は空間の拡張方法に無駄が多いので，効率は CIRR 法のほうが優れていると思われる。
参考文献
[6-1] A. Imakura, L. Du and T. Sakurai: "Accuracy analysis on the Rayleigh-Ritz type of the contour integral based eigensolver for solving generalized eigenvalue problems", Technical Report of Department of Computer Science, University of Tsukuba (CS-TR), CS-TR-14-23 (2014).
[6-2] N. Emad, S. Petiton and G. Edjlali: "Multiple explicitly restarted Arnoldi method for solving large eigenproblems", SIAM Journal on Scientific Computing, Vol. 27, No. 1, pp. 253-277 (2005).
(7) On quadratic convergence of the block Jacobi method in the presence of multiple eigenvalues
Yusaku Yamamoto (The University of Electro-Communications)
ブロックヤコビ法の大域的収束性については，Drmać [3-2]，Hari [3-1] らにより，行（列）方向の巡回型のピボット戦略の場合を中心に，収束証明が確立されている。一方，局所的収束性については，Drmać [3-2] により，固有がすべて単純である場合について，2次収束性が証明されている。また，Drmaćは，重複固有値の場合にも，ブロックの区分けを反復途中で変更することで，2次収束を実現できることを示唆している。しかし，分散メモリで並列化を行う場合には，各ブロックをノードに割り当てるので，途中で区分けを変更することは，通信や負荷分散などの点で望ましくない。そのため，重複固有値の場合でも，区分けの変更なしに2次収束を達成できる手法が望まれる。本研究では，2×2のブロックに対する直交変換の行列を修正することで，これを実現する手法を提案した。
提案手法を説明するため，まず，固有値がすべて単純で，その最小間隔が d の場合の2次収束性の証明を振り返る。ブロックヤコビ法における2×2の直交変換は，ピボットブロックの番号を (X,Y) とするとき，\(\begin{pmatrix}A_{XX}^{(k+1)} & A_{XY}^{(k+1)} \\ A_{YX}^{(k+1)} & A_{YY}^{(k+1)}\end{pmatrix} = \begin{pmatrix}P_{XX}^{(k)} & P_{XY}^{(k)} \\ P_{YX}^{(k)} & P_{YY}^{(k)}\end{pmatrix}^{\top} \begin{pmatrix}A_{XX}^{(k)} & A_{XY}^{(k)} \\ A_{YX}^{(k)} & A_{YY}^{(k)}\end{pmatrix} \begin{pmatrix}P_{XX}^{(k)} & P_{XY}^{(k)} \\ P_{YX}^{(k)} & P_{YY}^{(k)}\end{pmatrix}\)，あるいは \(\tilde{A}^{(k+1)} = \left(\tilde{P}^{(k)}\right)^{\top} \tilde{A}^{(k)} \tilde{P}^{(k)}\)と書ける。ここで，直交行列\(\tilde{P}^{(k)}\)は，\(\tilde{A}^{(k)}\)を対角化するように取る。いま，A(k)の非対角要素の絶対値がすべて δ（≪ d）以下ならば，sinΘ 定理により，PXY(k)，PYX(k) の2ノルムが O(δ/d) になるように（\(\tilde{P}^{(k)}\)の列を並べ替えて）できる。すると，ブロック間の回転角（直交変換のCS分解により定義される角度）は O(δ/d) となるから，ある非対角ブロックの更新における，他の非対角ブロック（既にノルムが O(δ) となっている）からの寄与は O(δ2/d) となる。したがって，巡回により一度消去されたブロックは，再び非ゼロになることがあっても，そのノルムは高々 O(δ2/d) に留まる。これが2次収束の理由である。
一方，重複した固有値λがあり，かつ，それらに対応する対角要素がブロック間にまたがっている場合は，sinΘ 定理による回転角の抑え込みができない。そのため，上記の証明が崩れてしまう。
そこで本研究では，重複固有値に対応する\(\tilde{P}^{(k)}\)の対角ブロック\(\tilde{P}_{\lambda}^{(k)}\)に注目した。これは，上記の2×2分割で言うと，2つのブロック行・列にまたがる対角ブロックである。このブロックを取り出し，極分解を行うことにより，2次収束を妨げる不要な回転成分を抽出し，それを\(\tilde{P}^{(k)}\)から取り除くことにした。その結果，ブロック間の回転角は，O(δ/d) に抑えられる。この場合，\(\tilde{A}^{(k+1)}\)は厳密な対角行列ではなくなるが，0 となるべき要素の大きさは O(δ2/d) になることが示せるので問題ない。こうして，重複固有値がある場合にも2次収束を達成するアルゴリズムが構成できる。本手法は，ブロックの区切りを変更する必要がないため，分散メモリにも問題なく適用できる。簡単な数値実験により，提案手法がうまく働くことを示した。
感想： 本発表はポスター発表である。扱っているのがかなり細かい問題であり，かつ，アイディアも非常に単純なので，IWASEP で発表するのはどうかと思ったが，Parlett，Drmać，Hari など，ヤコビ法に詳しい方々からは，非常に好意的な評価をいただいた。たぶん，今時ヤコビ法の収束性を研究している人など（クロアチアを除いては）いないので，励まして下さったのだろう。特に Hari は，重複固有値でもブロックの区分けを変更せずに済む点を，非常に気に入って下さった。私のアブストラクトを見て，自分で詳しい計算をされたようで，転置記号の付け忘れなども指摘して頂いた。
Hari からは，浮動小数点での演算に関しても，注意を頂いた。極分解が，ヤコビ法の持つ高い相対精度を壊す恐れがあるので，最初から，\(\tilde{P}_{\lambda}^{(k)}\)がなるべく対角行列に近くなるよう，注意して Givens 回転による消去を行うべきだという指摘である。我々は，ヤコビ法の高精度性にではなく，むしろ超並列計算機での効率性に着目してソルバを開発していたため，その点には注意が及ばなかった。実際，特異値問題ではなく固有値問題を解いている点で，既に精度面でのペナルティがあり，さらに，\(\tilde{A}^{(k)}\)の対角化にヤコビ法ではなくLAPACKの分割統治法を用いている点で，相対誤差の意味での高精度性は崩れてしまっている。そのため，極分解を気にせずに使っていたが，世の中ではヤコビ法と言えば高精度性を期待する場合も多いので，相対誤差の意味での高精度性を崩さずに重複固有値に対処する方法も検討してみたい。
Ipsen は，熱心にポスターを見て下さったが，いきなり，「なぜ極分解なのか。QR分解では駄目なのか」と質問されたのでびっくりした。確かに，\(\tilde{P}_{\lambda}^{(k)}\)はほとんど直交行列になっているので，QR分解でも近似的に回転成分を抽出でき，オーダーの評価に関しては（議論は少し難しくなるが）同じ結論が導けるのではないかと思う。他の講演での質問を聞いていても感じたが，短時間で内容を理解して本質的な質問をする能力はすごいと思った。
４. わからなかったが興味のある講演
面白そうだが難しくてわからなかった講演について，以下にまとめる。
(8) Variational characterization of eigenvalues of a non-symmetric eigenvalue problem governing elastoacoustic vibrations
Markus Stammberger and Heinrich Voss (Hamburg University of Technology)
穴の開いた弾性体と，その中に満たされた液体とからなる系について，小さい振幅の振動モードを記述する線形固有値問題の解析。弾性体の変位と液体の圧力とを変数に取ると，得られる偏微分方程式は，両者の間の境界条件のために非自己随伴型となるが，自己随伴型と似た性質を持つ。たとえば，固有値は実で，非負で，可算個であり，かつ，Rayleigh汎関数（Rayleigh商の一般化）を定義すると，固有値についてミニマックス型の特徴付けができる。また，固有関数は直交系をなすように取れる。これらの性質は離散化後も成り立ち，数値解法で利用できる。
(9) On structured pencils arising in Sonneveld methods
Jens-Peter M. Zemke (Technical Universitat Hamburg-Harburg)
IDR(s)法およびIDRStab法（GBi-CGSTAB(s,L)法）に基づく固有値解法では，特殊な形をした係数行列を持つ一般固有値問題 Ax = λBx が生じる。ここで，A はヘッセンベルグ行列において上三角部分がブロック3重対角である行列，B は上三角行列である。この問題を通常のQZ法で解くと，この非ゼロ構造が壊れてしまうため，これに適した効率的なソルバが必要である。(1, 2s-1) 型 quasi-separable 行列に対する Zhlobich の dqds 法 [9-3] が有望だと思われる。また，IDR法の特殊性のために，この一般固有値問題は無限固有値を持ち，それを予めデフレーションすることができる。また，残りの固有値を高い信頼性で計算するための手法を提案する。
感想： 一般固有値問題の導出とデフレーションの話は，説明が速過ぎて付いて行けなかったが，係数行列の形は面白いと思った。密行列系の解法では，特殊な非ゼロ構造を壊さずに，連立1次方程式や固有値問題を解く技法が多数研究されているが，そのような技法が大規模疎行列の解法でも必要となるのは興味深い。疎行列と密行列の研究者の連携により，より効率的な解法が作れるのではないかと思われる。なお，本発表では，谷尾・杉原のGBi-CGSTAB(s,L)法についても言及されていた。
参考文献
[9-1] M. H. Gutknecht and J-P. M. Zemke: "Eigenvalue computations based on IDR", SIAM Journal on Matrix Analysis and Applications, Vol. 34, No. 2, pp. 283-311 (2013).
[9-2] O. Rendel, A. Rizvanolli and J-P. M. Zemke: "IDR: a nw generation of Krylov subspace methods?", Linear Algebra and Its Applications, Vol. 439, No. 4, pp. 1049-1061 (2013).
[9-3] P. Zhlobich: "Differential qd algorithm with shifts for rank-structured matrices", SIAM Journal on Matrix Analysis and Applications, Vol. 33, No. 4, pp. 1153-1171 (2012).
(10) Low-rank updates to matrix functions
Ana Susnjara (EPFL SB MATHICSE ANCHP) and Daniel Kressner
exp(A) や sign(A) などの行列関数において，A に低ランクの摂動を加えた場合に関数を効率的に再計算する手法についての研究。テンソル型の多項式クリロフ部分空間と有理クリロフ部分空間を使う。本発表はポスター発表である。
５. 雑感
６. 写真
学会会場の CAAS とその付近
ドブロブニクの旧市街。
夜8時過ぎまで明るく，10時近くになっても大勢の観光客がいた。
旧市街を囲む城壁と塔，塔から見た旧市街の全景。
エクスカーションで行った Trsteno の植物園と途中の風景
学会出張報告のページに戻る
Topに戻る