１. 学会の概要
本会議は SIAM（米国応用数理学会）主催で3年に1回開催される線形代数の応用に関する国際会議である。線形代数に関してはおそらく世界最大の会議であり，今回は12件の招待講演，75のミニシンポジウム（4講演／セッション），42の一般講演セッション（4講演／セッション）があった。講演の総数は450件程度である。内容は線形代数のほとんどあらゆる応用にわたり，線形計算はもちろん，制御理論，マルコフ連鎖，画像処理，グラフ理論などへの応用も含む。今回も Demmel，Higham，van Loanをはじめ，著名な研究者が多数参加していた。
２. 目立ったテーマ
今回の会議では，次のようなテーマが目立った。
３. 興味深かった講演
特に興味深かった発表と自分の発表について，以下にまとめる。なお，階層的テンソル分解に関するGrasedyckの講演，クリロフ部分空間法における基底の数値的直交性に関するPaigeの講演については，まだまとめられていない。両方とも大変興味深い講演だったので，今後勉強してまとめたい。
(1) Multi-Window Approach to Deflation in the QR Algorithm
Karen Braman (South Dacota School of Mines and Technology)
Hessenberg QRアルゴリズムにおける非常に効果的なデフレーション戦略として，Aggressive Early Deflation（AED）がある。AEDでは，反復途上の行列の右下のk×k部分行列（ウィンドウ）を取り，そのSchur分解を計算して，それに対応する直交変換（を拡大したもの）を行列全体に施す。すると，右下k×k行列は上三角行列になり，その左横に，スパイクと呼ばれる長さkの非ゼロ要素の列ができる。もしこのスパイクの要素のいくつかが十分小さければ，それらを0で置き換え，再び行列全体に直交変換の逆変換を施すことで，多くの副対角要素を一度に0にすることができる。
本発表では，このウィンドウを右下部だけでなく対角線に沿って複数取り，それぞれに対してSchur分解を計算するマルチウィンドウ型のAEDを提案した。この場合，右下以外のウィンドウに対しては，全体行列に直交変換を施すと，上三角行列の左と下の両方にスパイクができる。このスパイクのうち，k個以上の要素を0と見なせるならば，デフレーションを行えることが示された。
全体行列の中央部でデフレーションを行えれば，固有値問題が複数の小さい固有値問題に分割できるので，計算量の面から非常に有利である。数値実験の結果，このようなデフレーションはいつも可能とは限らないが，いくつかのテスト例題については大きな効果が得られることが示された。また，マルチウィンドウ型のAEDは，シフトの数，複数段のバルジ追跡をまとめてブロック化する際の段数，ウィンドウのサイズ，個数，位置など多くのパラメータを持つため，これらを最適化することの重要性が指摘された。
(2) Aggregation of the Compact WY Representations Generated by the TSQR Algorithm
Yusaku Yamamoto (Kobe University)
m×nの縦長行列Aに対し，そのQR分解を計算し，別の行列m×l行列Bに対してQTBを計算する問題を考える。このような問題は，ブロック化QR分解，対称行列の帯行列化など，様々な場面で現れる。この計算において，QR分解をマルチコアCPUで行い，QTBの計算をGPUなどのアクセラレータで行う場合，前者には並列粒度の大きいTSQRアルゴリズムを使うのが適切である。しかし，TSQRアルゴリズムで生成されるQ因子は，従来のハウスホルダーQR法で生成されるQ因子に比べてサイズが1/p（pはマルチコアプロセッサのコア数）と小さい複数の因子からなり，GPUでQTBの計算を高速に行うのが難しい。
そこで，TSQRアルゴリズムで生成されるQ因子を合成して，一つの大きなQ因子にすることが望まれる。このとき，計算の効率化のため，Q因子はコンパクトWY表現という形式で保持する必要がある。ところが，簡単な数学的議論により，TSQRアルゴリズムで生成されるQ因子は，従来法で生成されるQ因子と同じサイズのコンパクトWY表現では表せないことが示される。そこで本研究では，Aを上三角行列Rに変換するm×m直交行列は無限個存在し，応用の観点からは，多くの場合，そのうちのどれでも等価であることに着目した。そして，うまく直交行列Q'を選べば，AをRに変換するという性質と，従来法と同じサイズのコンパクトWY表現で表せるという性質が両立できることを示した。このようなQ'のコンパクトWY表現は，TSQR法のQ因子から，2mn2の演算量で計算できる。この計算を，コンパクトWY表現の合成と呼ぶ。
本手法を4コアXeon + Tesla C1060という環境で評価したところ，n=100，m=l=5120の場合に10%程度の高速化が達成できることが示された。また，本手法におけるQ'の直交性および残差A - Q'Rは，TSQR法の場合と同程度であり，ハウスホルダーQR法を用いた場合に比べて，より良い結果となった。今後の課題としては，具体的なアルゴリズムへの応用と理論的な誤差解析が挙げられる。
質疑：
Q1: 本手法で計算した行列Q' = I - YSYTは直交行列になるのか？（Julien Langou）
A1: （式を見せて直交行列になることを説明）
Q2: 列選択付きQR分解には適用できるか？
A2: TSQR法は列選択付きQR分解には対応していないので，現状ではできない。なお，（質問の趣旨とずれるが）アルゴリズム中で必要となる連立1次方程式の求解を列選択付きQR分解を用いて行えば，悪条件の問題の場合により良い精度が得られると考えている。
Q3: 対称行列の帯行列化のような両側変換への適用は？（van Loan）
A3: いったんコンパクトWY表現を構築してしまえば，従来法と同じにできる。
感想： 今回の発表は，今までこの会議で何回か発表した中では一番反応がよく，TSQRアルゴリズムの提唱者であるJulien Langou，UCBの学生，INRIAの研究者などからスライドが欲しいと言われた。Langouによると，TSQRで生成されるQ因子をまとめて大きくできないかという相談は，以前から何件か来ていたという。今回の手法では，実験的にはTSQRアルゴリズムと同等の精度が達成できているが，連立1次方程式を解くというプロセスを含むため，理論的な誤差解析が難しい。係数行列の条件数，および軸選択付きLU分解における成長因子（growth factor）について，良い上界を得ることがまだできていない。今後，さらに理論的解析を進める必要がある。
追記： 誤差解析がうまく行かずにしばらく放っている間に，Demmelたちのグループが誤差解析付きの類似のアルゴリズムを発表してしまった（参考文献[2-1]）。上記の発表を引用し，その中で使っているLU分解の部分を修正することで，後退安定性を証明しているようである。発表の反応が良かったので喜んでいたが，反応が良いということは，他のグループも同じ研究をやり始める可能性が高いということなので，もっと切迫感を持って急いで誤差解析までやり遂げて論文にすべきであった（2013. 12/16）。
参考文献
[2-1] Grey Ballard, James Demmel, Laura Grigori, Mathias Jacquelin, Hong Diep Nguyen and Edgar Solomonik: "Reconstructing Householder Vectors from Tall-Skinny QR", Technical Report No. UCB/EECS-2013-175, Electrical Engineering and Computer Sciences, University of California at Berkeley, October 26, 2013.
(3) A Generalized SVD for Collections of Matrices
Charles van Loan (Cornell University)
2つの行列X1∈Rm1×n，X2∈Rm2×n（ただしm1, m2≧n）に対し，X1 = U1Σ1VT，X2 = U2Σ2VTを満たす列直交行列U1∈Rm1×n，U2∈Rm2×n，n×n対角行列Σ1，Σ2，n×n直交行列Vが存在することが知られている。この分解を，一般化特異値分解（GVSD）と呼ぶ。本発表では，行列がX1，...，XNとN個存在し，それぞれが列フルランクである場合に対して，GSVDの拡張であるHigher order GSVD（HO GSVD）を提案した。
いま，もしすべてのXiが共通の行列Vを用いてXi = UiΣiVTと書けたとする。このとき，Ci=XiTXiとおくと，簡単な計算より，Vの各列は，任意のi，jについて，CiCj-1の固有ベクトルになっていなければならない。さらに，行列の形を対称化して，S=(1/N(N-1))Σi>j (CiCj-1+CjCi-1)とおくと，Vの各列はSの固有ベクトルになっていなければならないことがわかる。
そこで，行列Sの固有値問題について考えると，Sは対角化可能であり，固有値はすべて実数であり，1以上であることが示せる。さらに，ある固有値がλk=1を満たすとき，対応する固有ベクトルvkは，任意の2つの行列Xi，Xjのペアに対して，そのGSVDにおける右特異ベクトルになっていることが示せる。このような固有ベクトル群の張る空間を，common HO GSVD subspaceと呼ぶ。Common HO GSVD subspaceを計算することで，行の数が必ずしも等しくない複数の2次元データに対して，それらの間に共通な特徴を自動的に抽出することができる。λが1ではないが1に近い値であるとき，何が言えるかについては，今後の課題である。
感想： 行列の同時対角化に関してはこれまでに多くの研究があるが，いずれも同じ型の行列が多数ある場合を対象としている。本手法は，同時対角化ではなく，一部の共通な特徴のみを抜き出す手法であるものの，各行列の行の数が違うときにも適用できるという点で興味深かった。反復法でなく，固有値問題を1回解くだけで特徴抽出ができるのも興味深い。
(4) Divide and Conquer the CS Decomposition
Brian D. Sutton (Randolph-Macon College)
X=[X1T, X2T]Tを列直交行列とし，X1∈Rm1×n，X2∈Rm2×n（ただしm1, m2≧n）とする。このとき，X1，X2は，列直交行列U1，U2，C2 + S2=Iを満たすn×n対角行列C，S，およびn×n直交行列Vを用いて，X1 = U1CVT，X2 = U2SVTと表現できる。これを（2×1の）CS分解と呼ぶ。CS分解は，一般化特異値分解において入力行列が列直交行列である特別な場合であり，原理的には一般化特異値分解のルーチンを用いて計算できる。しかし，この方法では，Vの数値的直交性が保証できないという問題があった。
本研究では，CS分解に対して，2つのステップからなる解法を提案した。第1のステップでは，X1とX2の同時上2重対角化を行う。特異値分解の前処理の場合と同様に，まずX1，X2にそれぞれ左からハウスホルダー変換HL1(1)，HL2(1)をかけて第1列の第2要素以降を0にする。次に，右から共通のハウスホルダー変換HR(1)をかけて，X1の第1行の第3要素以降を0にする。このとき，Xの列直交性（左右からの直交変換によって保存される）により，X1の第1行の第2要素以降とX2の第1行の第2要素以降とは線形従属になっているので，X2の第1行の第3要素以降も自動的に0になる。これを繰り返すことで，X1，X2はそれぞれ上2重対角行列B1，B2に変換される。
第2のステップでは，B1，B2に対して分割統治法を適用し，対角化を行う。ただし，不用意に分割統治法を適用すると，数値誤差のために得られたCとSがC2 + S2=Iを満たさなくなる。そこで，B=[B1T, B2T]Tが列直交行列であることを利用し，B1，B2の要素をcosとsinで表現する。そして，行列要素どうしの加算，乗算などの演算を，三角関数の加法定理を逆に用いて，cos，sinの引数に対する演算として行うようにする。これにより，数値誤差があっても，C2 + S2=Iという性質は，丸め誤差レベルで成り立つ。本アルゴリズムは後退安定であることが示される。
感想： 行列要素をcos，sinで表現することにより，数値誤差があってもC2 + S2=Iを強制的に成立させるというアイディアは，非常に面白いと思った。なお，同時上2重対角化の説明において，X2の(1,1)要素が0の場合は，右からHR(1)をかけてもX2の第1行の第3要素以降は0にはならないが，この場合，実はX2の第1行の第3要素以降がもともと0になっている。したがって，X2の第1行の第3要素以降を0にするようにHR(1)を決めればよい。しかし，X2の(1,1)要素が0に非常に近い場合はどうするのだろうか。たとえばX1とX2の(1,1)要素どうしを比べて，絶対値が小さいほうの行を基にHR(1)を生成すれば，数値的安定性が保たれるのだろうか（これがvan Loanの質問の趣旨だったように思う）。まだ著者のウェブサイトにプレプリントが掲載されていないので，アルゴリズムの詳細はわからなかったが，掲載されたら読んでみたい。
(5) Hierarchical Tensor Decomposition and Approximation
Lars Grasedyck (RWTH Aachen)
(6) Backward Stability of Iterations for Computing the Polar Decomposition
Nicholas Higham and Yuji Nakatsukasa (University of Manchester)
任意のn×n行列Aは，ユニタリ行列Uと非負定値エルミート行列Hにより，A=UHと表現できる。これを極分解（polar decomposition）と呼ぶ。極分解は様々な応用を持つ。たとえば，Aの列ベクトルの張る空間の正規直交基底は，極分解を使うとQR分解より精度よく計算できる。また，最近では，極分解に基づくcommunication-avoiding型の固有値計算法も提案されている。
極分解に対しては，スケーリング付きNewton反復，動的重み付きHalley反復など，様々な反復法が提案されている。しかし，これらの反復法の安定性については，十分知られていなかった。
本発表では，ユニタリ因子に対する反復式Xk+1 = f(Xk) が，(i) 混合型の安定性を持ち，(ii) Xk の特異値を大きく減少させない，という2つの条件を満たす場合に，この反復が後退安定であることを示した。この解析によると，動的重み付きHalley反復は，列ピボット選択および行のスケーリング／ソーティング付きのハウスホルダーQR法を用いた場合に後退安定，スケーリング付きNewton反復は，逆行列を混合型の安定性を持つアルゴリズムで計算した場合に後退安定であることが示される。また，Newton-Schulz型の反復は，特定の条件の下でしか安定でないことも示される。
(7) Improving Performance and Robustness of Incomplete Factorization Preconditioners
Anshul Gupta (IBM T. J. Watson Research Center)
不完全コレスキー分解や不完全LU分解は，クリロフ部分空間法のための最も一般的な前処理法であるが，高性能計算の観点から見ると，十分な最適化を行わずに使われていることが多い。本講演では，IBMの疎行列計算パッケージWSMP（Watson Sparse Matrix Package）の開発者として知られるGupta氏が，疎行列直接法（sparse direct solver）の技法を利用した不完全分解の高性能実装法について報告した。
対称行列向けの疎行列直接法では，スーパーノードを利用して最適化を行う。スーパーノードとは，コレスキー因子において，ほぼ同じ非ゼロ構造を持つ連続した列の集合である。これらの列を一まとめにして消去演算を行うことで，間接参照を大きく削減し，かつ，密行列向けのBLASを利用できる。場合によっては，多少疎性を犠牲にしても，隣接する列集合をスーパーノード化したほうが性能上有利である場合も多い。これらのことは不完全分解の場合も同様であることが述べられた。ただし，不完全分解の場合は，ある閾値τに基づき0にする要素を決めるため，実際に数値的分解を行ってみないとスーパーノードを決められないのが，疎行列直接法と異なる点である。
不完全分解では，閾値τを変えることで，前処理行列としての精度と前処理にかかる演算量・メモリ量のトレードオフを調節できる。一方，反復法としてGMRES法を使う場合は，リスタート回数mを変えることで，収束性とメモリ量・反復当たりの演算量のトレードオフを調節できる。不完全分解による前処理を用いたGMRES法では，τとmとを同時に調節することで，所要メモリ量を一定に保ちつつ，計算時間を最適化できる。同じタイプの連立1次方程式を複数回解く場合は，自動チューニングも有効である。
(8) Hyperdeteminant and the Condition Number of a Multilinear System
Lek-Heng Lim (University of Chicago)
連立1次方程式の条件数は，その係数行列から最も近い非正則行列までの距離の逆数と解釈できる[8-1]。同じことが，テンソルを係数とする多重線形方程式の場合にも言える。この場合，非正則なテンソルとは，超行列式（hyperdeterminant）が0のテンソルである。超行列式が0という条件は，テンソルの特異値が0であることとも言い換えられる。また，対称なテンソルの場合，この条件は，固有値が0であることとも言い換えられる。テンソルをk次元の超行列と見なすと，テンソルの正則性は，そのすべての面に対して，同じガウス消去法あるいは同じハウスホルダー／ギブンス変換を適用しても変わらないことが示せる。
参考文献
[8-1] J. W. Demmel: "Applied Numerical Linear Algebra", SIAM, Philadelphia, 1997.
(9) Orthogonality and Stability in Large Sparse Matrix Iterative Algorithms
Chris Paige (McGill University) and Wolfgang Wiiling (W2 Actuarial & Math Services Ltd.)
(10) Prescribing the Behavior of the GMRES Method and the Arnoldi Method Simultaneously
Jurjen Duintjer Tebbens (Academy of Sciences of Czech Republic) and Gerard Meurant (Commissariat a l'Energie Atomique, Paris)
固有値計算のためのArnoldi法では，各ステップでHessenberg行列のサイズを1ずつ拡大する。第kステップでは(k+1)×kのHessenberg行列H(k)ができ，その最初のk行からなる正方行列の固有値（Ritz値）μ1(k)，…，μk(k)が，Aの固有値の近似値として採用される。これは，kが大きくなるにつれ，Ritz値が徐々にAの固有値の良い近似値になっていくという期待に基づいている。
ところが本発表では，これが必ずしも成り立たないことが示された。すなわち，{μ1(k)，…，μk(k)}（k=1, 2, …, n-1）という数列（の数列）は，Aの固有値（すなわち{μ1(n)，…，μn(n)}）と独立に，どんな値でも取りうることが示された。したがって，たとえば，最初のn-2ステップでは，Ritz値が徐々にAの固有値に近づいていき，第n-1ステップでこれらが全くでたらめな値を取る，というようなパターンも可能である。これは，任意のk=1, 2, …, nに対して，第k×k首座小行列の固有値が{μ1(k)，…，μk(k)}となるようなHessenberg行列を実際に構成することによって示される（第kステップまでに現れるRitz値の総数は，下側副対角が1のk×k Hessenberg行列の要素の自由度と等しい）。特に，このHessenberg行列の下側副対角要素は1に取ることができる。
さらに，この行列に対して対角行列による相似変換を行うと，首座小行列の固有値を変えずに，下側副対角要素を任意の値にすることができる。下側副対角要素は，各ステップにおけるGMRES法の残差に関係付けられるから，これは，GMRES法の収束履歴がRitz値の収束とは無関係でありうることを意味する。ただし，Ritz値が0になる場合にGMRES法が停滞することは知られており，その点についてだけは，両者に関係がある。
感想： Ritz値は元の行列の固有値に徐々に収束していくものと漠然と思っていたが，理論的にはどんな可能性もありうるのだと知って驚いた。導出も難しくなく，面白い結果だと思う。ただし，多くの問題では，Ritz値はやはり元の行列のよい近似になっている（そうでないとArnoldi法は役に立たない）。その点を説明できる理論があればと思う。
(11) SIGLA Prize: Communication-Avoiding Algorithm
James W. Demmel (UC Berkeley)
近年の計算機では，メモリアクセスやプロセッサ間通信などのデータ移動が，性能上の最大のボトルネックとなっている。そのため，プロセッサの演算性能を最大限に活用するには，データ移動量をできるだけ抑えたアルゴリズムが必要である。
本研究では，各プロセッサが主メモリに加えて容量M語の高速メモリを持つと仮定する。また，それらのプロセッサが複数個，ネットワークにより接続されていると仮定する。このとき，演算量がO(n3)となる密行列演算（行列乗算，LU分解，QR分解，特異値分解，テンソルの縮約など）に対して，プロセッサあたりのメモリアクセス量（主メモリから高速メモリへのデータ移動量）が，（プロセッサあたりの演算量）／M1/2という下界を持つことを示した。さらに，プロセッサあたりの通信量が，（プロセッサあたりの演算量）／M3/2という下界を持つことを示した。これらの下界は，高速メモリに格納する部分行列の形を長方形に限った場合には比較的容易に証明できる。これに対して本研究では，任意の3次元領域の体積と，その領域をxy，yz，zxの各平面に射影してできる図形の面積との間に成り立つ不等式を用いることで，どのように行列要素を抜き出しても，この下界を下回ることはできないことを示した。
この下界に照らしてLAPACKで使われているアルゴリズムを評価すると，多くの場合，現状のアルゴリズムは最適ではない。そこで，本研究では，様々な密行列演算に対して，下限を達成するアルゴリズムを新たに開発した。具体例として，2.5次元行列乗算とTSQR（Tall-and-Skinny QR）アルゴリズムが挙げられる。2.5次元行列乗算では，n=8,192の乗算を16,384ノードのBlueGeneで行う場合に，通信時間を95%削減できた。また，TSQRアルゴリズムでは，細長い行列のQR分解に対し，多くの場合に従来法の数倍の高速化が得られている。さらに，ピボット選択付きのLU分解に対しても，TSQRアルゴリズムと同様の原理に基づくcommunication-avoiding型のアルゴリズムが得られている。
感想： プロセッサの演算性能を引き出すにはメモリアクセスや通信量の削減が不可欠だということは，現在では常識である。しかし，どこまで削減できるかという下界を一般的に示し，その下界を達成するアルゴリズムを様々な行列計算に対して具体的に構築して見せた点は，非常に迫力があった。今回の学会でも，communication-avoiding型のアルゴリズムというのは，大きなトレンドになっている。こういうトレンドを作り出せるという意味で，Demmelはやはりすごいと思う。
(12) Solving Nonlinear Eigenvalue Problems in Electronic Structure Calculation
Chao Yang (Lawrence Berkeley National Laboratory)
密度汎関数理論に基づく電子状態計算では，Kohn-Sham方程式が基本方程式となる。これは，準粒子の1電子エネルギーと1電子波動関数を求める固有値問題であるが，演算子中のポテンシャルの項が，固有関数である1電子波動関数ψiに電子密度ρを通して依存するため，非線形固有値問題になっている。そのため，通常は，まずポテンシャルを固定して線形固有値問題を解いて1電子波動関数を求め，そこから電子密度を計算して再びポテンシャルを求める，という操作を，ポテンシャルが収束するまで繰り返す。しかし，この方法は，固有値問題を解く必要があることと，self-consistent計算の反復必ずしも収束するとは限らないという問題点がある。
一方，1電子波動関数をあらわに扱わない方法として，密度行列に基づく方法がある。この方法では，まず電子の分布関数における不連続性を緩和するために有限の温度（の逆数）βを導入し，密度行列 Σi |ψi>(2/(1+exp(β(Ei-μI))))<ψi| を基に，密度をはじめとするすべての物理量を計算する。密度行列は1電子ハミルトニアンHを用いて Σi 2/(1+exp(β(H-μI))) と書き直すことができ，特に密度ρはそのトレースとなる。したがって，密度行列のトレースを計算できれば，1電子波動関数を求めることなく（すなわち固有値問題を解くことなく），ポテンシャルから電子密度を求めることができ，そこからまたポテンシャルを求める，というself-consistent計算を行うことができる。
密度行列のトレースを計算するには，Fermi演算子 2/(1+exp(β(H-μI))) を複素平面上の周回積分として表現し，数値積分を行うことで，Hに関する部分分数展開の形に近似する。したがって，(H-zkI)-1 のトレースの計算が必要となる。そこで，この行列に対してLU分解を行い，LU因子から逆行列のトレースを効率的に計算するErisman & Tinneyのアルゴリズム[12-1]を用いる。一方，self-consistent計算の収束性を保証するには，ポテンシャルから次のポテンシャルを求める写像において，ヤコビアンのスペクトル半径が1より小さくなる必要がある。しかし，特に金属系では，系のサイズの増大と共にスペクトル半径が大きくなり，収束しにくくなることがわかっている。減速を行えば収束するものの，収束は遅くなる。これに対してはいくつかの改良方法があり，半導体や絶縁体では成功を収めているものの，金属系ではまだ良い方法が見つかっていない。
参考文献
[12-1] A. M. Erisman and W. F. Tinney: On Computing Certain Elements of the Inverse of a Sparse Matrix, Communications of the ACM, Vol. 18, No. 3, pp. 177--179 (1975).
感想： 大規模電子状態計算における最先端の研究が紹介され，興味深かった。逆行列のトレースという量は，dqdsアルゴリズムにおけるシフトの計算や，複素周回積分に基づく非線形固有値の解法でも出てくる量であり，Erisman & Tinneyのアルゴリズムも使ったことがある。同じ量が電子状態計算でも重要な役割を果たすのは興味深い。
(13) Hierarchical QR Factorization Algorithms for Multi Core Cluster Systems
Julien Langou (Univ. of Colorado, Denver), Jack Dongarra, Mathieu Faverge (Univ. of Tennessee), Thomas Herault (Univ. Paris-Sud) and Yves Robert (Ecole Normale Superieure de Lyon)
コレスキー分解，QR分解などの様々な行列分解に対し，行列をL×Lのタイル（ブロック）に分割してタイル単位で処理を行うタイルアルゴリズムが提案されている。タイルアルゴリズムはデータの局所性が高く，キャッシュメモリの有効利用や通信階数の削減などの点で優れており，次世代の密行列アルゴリズムとして注目されている。
M×N個のタイルからなる行列に対してタイルアルゴリズムによるQR分解を行う場合，第K段では，第K列の対角タイル以下にあるタイルの間で消去演算（2個のタイルをまとめて2L×Lのタイルとし，それに対してQR分解を行うことにより，下側のタイルを0にする）を行い，最終的に対角タイルを上三角行列，それより下のタイルをすべて0にする。消去演算は，この条件さえ満たせば，どのタイルでどのタイルを消去してもよい。したがって，第K列の消去演算は，葉の数がM-K+1の2分木により表現される。なお，各消去演算に伴い，第K+1列以降のタイルに対しても同じ直交変換による更新を行う必要がある。
本研究では，ハードウェアとしてマルチコアのクラスタを想定し，最適な2分木の構成法を示した。具体的には，2分木を3つの階層に分割し，最上位の2分木をノード間での消去に対応させる。一方，中間の2分木はノード内でのコア間の消去に対応させ，最下位の2分木はコア内での消去に対応させる。最初の2つの2分木は，並列性を重視し，なるべく高さが低い均等な2分木とする。一方，最後の2分木は，コア内でのキャッシュメモリの有効利用が目的のため，コアが担当する中で一番上のタイルで他のすべてのタイルを消去するタイプの2分木とする（これにより，三角のタイルで四角のタイルを消去する，効率的なTSカーネルが使える）。コアの割り当ては，ノード内・ノード間を合わせた全コアを長方形のテンプレートに並べ，これを敷き詰める形で2次元のサイクリック割り当てを行う。
上記の方針では，第K段のみに着目すると効率的な計算が行えるが，次の段との関係を考えると，パイプライン化ができないという問題がある。なぜなら，第K+1段における各コア内での消去では，再び各コアが持つ一番上のタイルで他のタイルを消去するが，このタイルは，すぐ左隣のタイルが第K段においてコア間の消去によって0にされ，それに対応する更新を受けて初めて消去に使える状態になるからである。したがって，第K段におけるコア内の消去がすべて完了するまで，第K+1段の消去が待たされてしまう。そこで，第K+1段の消去では，各コアが持つ上から2番目のタイルで消去を行うことにする。このタイルの左隣のタイルに対する消去は，第K段の最初に行われるので，対応する更新を受けてすぐ，このタイルは第K+1段目の消去に使えることになる。これにより，第K段目の消去と第K+1段目の消去を重ねてパイプライン的に実行できる。このような消去の順番をcoupling level treeと呼ぶ。
以上のアルゴリズムを，DAGuEスケジューリングツールを用いて実装し，4コアXeon×2のノードからなる64ノードのマシンで，300,000×4,480までの行列を用いて評価したところ，2,505GFLOPS（ピーク性能の57.5%）と極めて高い性能が得られた。これは，ScaLAPACKに比べて約10倍の性能である。
感想： 密行列のQR分解に関する最新の研究成果の報告であり，極めて興味深かった。実装に当たっては，消去の順番を記述してDAGuEに入力するだけで，適切なスケジューリングが決定され，ノード間通信を含むプログラムが自動的に生成されるというのもすごいと思う。これからは，OpenMPとMPIを使って自分で並列化を行い，手動でパイプラインを最適化するというような牧歌的なプログラミングは過去のものとなるのだろうか。
なお，coupling level treeという呼び名の適切さには疑問が残った。本発表における2分木は，本来，各列内での消去の順番を記述する手段として定義されている。ところが，coupling level treeと言っているのは，効率のよいパイプライン並列化を行うには，最下位の（コア内消去に対応する）2分木をKにつれてどのように変えていけばよいかという，その変形の方法なのである。したがって，実際に存在するのは，各Kにおける最下位の2分木だけであり，それ以外にcoupling level treeというものがあるわけではない。誤解を避けるため，別の名称にしたほうがよいのではないかと思う。
(14) Randomized Algorithms in Linear Algebra
Petros Drineas (Rensselaer Polytechnic Institute)
大規模な2次元データ A∈Rm×n （m≧n）から何らかの特徴抽出を行う場合，行列の特異値分解A=UΣVTがよく使われる。しかし，Aが巨大になると，特異値分解の計算は困難になる。また，最大特異値に対する特異ベクトル，たとえば左特異ベクトルは，元のデータ中の各列が持つ特徴をよく表してはいるものの，元のデータ中には存在しない合成された列ベクトルであり，しばしば解釈が困難なことがある。たとえば，整数を要素とする行列に対しても，左特異ベクトルの要素は一般に実数となるから，解釈をしにくいといったことが起こりうる。
そのため，特異ベクトルの代わりとして，Aの列ベクトルの中から，Aの特徴を最もよく表すベクトルを抜き出そうという問題が考えられる。さらに，その一般化として，CをAの列ベクトルをs本（s≦n）抜き出して並べたm×sの行列，Xをs×nの適当な行列とするとき，‖A-CX‖を最小化するようにCとXを定めるという問題が考えられる。これをCX分解という。この問題は，数値線形代数で研究されてきたColumn Subset Selection Problem（CSSP）と深い関係がある。行列ノルムとしてフロベニウスノルムを取った場合，X=C+A（C+はCの Moore-Penrose 一般逆行列）とするのが最適であることがわかるので，問題はCを決めることに帰着される。
この問題に対して，Aの各列をその2ノルムの2乗に比例する確率で抜き出し，Cの列と置き換える操作を繰り返すアルゴリズムを考えることができる。このアルゴリズムについて，誤差の期待値の上界を E[‖A-CC+A‖F2] ≦ ‖A-Ak‖F2 + (4k/s)1/2‖A‖F2 （1≦k≦rank(A)）と求めることができる。これは，本アルゴリズムの誤差が，特異値分解を使ったランクkの最良近似の誤差（右辺第1項）とランダム化による誤差（右辺第2項）の和で抑えられることを意味する。しかし，この第2項は ‖A‖F2 に比例する大きな量であり，かつ，s=k と取った場合でも0にならないという問題点がある。
そこで，別の形の誤差評価式を持つアルゴリズムを考える。いま，VTにおいて，大きいほうからk個の特異値に対応する行のみを抜き出してできるk×n行列をVkTとし，その第j列ベクトルを (VkT)(j) とする。このとき，Aの第j列を ‖(VkT)(j)‖22 に比例する確率でサンプリングすると，そのときの誤差は，確率0.9以上で ‖A-CC+A‖F2 ≦ (1+ε)‖A-Ak‖F2 となることが言える。ただし，ε=3k/s である。すなわち，このサンプリング方式では，（誤差の期待値の評価ではなく誤差の確率的評価になるが）相対誤差の意味での評価式を得ることが可能となる。‖(VkT)(j)‖22 を（jに関する和が1になるように）正規化したものを leverage score と呼ぶ。
ただし，leverage score を計算するには，一般に O(mn2) の演算量が必要である。そこで，前処理により，より少ない演算量で leverage score を計算する方法を考える。いま，Dを対角要素に+1，-1を等確率で並べたn×nのランダム対角行列，Hをn×nのウォルシュ・アダマール変換の行列とする。このとき，A'=ADH とすると，A' は（log n 倍の違いを無視すれば）一様な leverage score を持つことがわかる。そこで，A' の近似は，leverage score を計算しなくても，一様なサンプリングにより行える。サンプリングで得られた行列を C' とすると，確率0.9以上で ‖A'-C'C'+A'‖F2 ≦ (1+ε)‖A'-A'k‖F2 が成り立ち，さらに，DとHが直交行列であることに注意すると，‖A-C'C'+A‖F2 ≦ (1+ε)‖A-Ak‖F2 が言える。このアルゴリズムの演算量は，高速ウォルシュ・アダマール変換を利用すると，O(mn polylog(n) + mk2polylog(n)) にできる。
感想： CX分解を例題に，線形計算におけるランダマイズドアルゴリズムの設計法や誤差評価の方法が紹介され，わかりやすかった。特に，行列の摂動理論という古典的な結果と，確率論，高速JL変換などの道具が巧妙に組み合わされて，O(mn2) という演算量の壁を破るアルゴリズムが生み出されるところが興味深かった。伝統的な線形計算の考え方からすると，ランダム化を取り入れるのにはやや抵抗があるが，本研究が対象としている行列（ゲノム情報処理におけるデータマイニング）のように，行列自体がもともとサンプリングによって得られたものの場合には，行列に対する計算においてもサンプリングにより効率化を図るのは自然なアイディアだと思った。あと，たとえば特異値分解では行列の2重対角化という前処理を行うことで，演算量の削減を図るが，ランダマイズドアルゴリズムでも同様に，固有の前処理を行うことでサンプリングを容易にし，演算量を削減できるというのは面白いと思った。
(15) Accurate Computations for Rational Bernstein-Vandermonde and Said-Ball-Vandermonde Matrices
Jorge Delgado (Universidad de Zaragoza) and Juan M. Pena (Universidad de Zaragoza)
計算機援用の形状設計（Computer-Aided Geometric Design）で有用な基底関数として，Rational Bernstein基底およびSaid-Ball基底がある。これらの基底に対するcollocation行列A=(aij)=(bi(tj))（biは基底関数，tjは座標）に対する連立1次方程式の求解や固有値・特異値計算は形状設計において重要である。
これらのcollocation行列は，Totally Nonnegative（TN）であることが示されている。そこで，TN行列に対するKoevの高精度アルゴリズムを適用することが可能である。ただし，そのためには，行列を2重対角分解の形で表現し，各2重対角行列の行列要素を高精度で求めておく必要がある。
本研究では，上記の2つのcollocation行列に対し，2重対角分解の要素をt_jの関数として陽的に表す式を導出し，これらの要素を高精度に計算することに成功した。この結果に基づき，Koevのアルゴリズムを適用した結果，連立1次方程式の求解，固有値計算，特異値計算などを高い相対精度で行うことができた。
感想： TN行列に対する数値解法の研究は，高い相対精度で解を計算できるという理論的興味が先行し，実際の応用がなかなか追い付いてこないという弱点がある。本発表では，Rational Bernstein-Vandermonde行列とRational Said-Ball Vandermonde行列という2種類のTN行列が，計算機援用の形状設計という実際の応用に現れることを示した点に意義があった。ただ，これらの行列を係数とする連立1次方程式，あるいはこれらの行列に対する固有値・特異値計算が，形状設計という文脈においてどういう意義を持つのかについても，説明が欲しかった。本発表で出てきた行列がHessenberg行列，あるいはせめて帯行列であれば，我々の提案した高精度固有値計算法も使えるが，実際には密行列であったため，我々の解法との間にはまだギャップがあるのが残念であった。
(16) Accurate Jordan Structures of Irreducible Totally Nonnegative Matrices
Plamen Koev (San Jose State University)
行列のジョルダン標準形を数値計算により求めることは，一般には数値誤差のために不可能である。しかし，既約なTotally Nonnegative（TN）行列の場合，0以外の固有値はすべて単根であることが，Fallat & Gekhtmanによって示されている。したがって，固有値0のジョルダン構造のみを求めればよいことになる。本発表では，TN行列の2重対角分解を用いた高精度アルゴリズムにより，これが数値計算によって可能であることを示した。
提案手法では，2重対角行列の積として表現されたn×nのTN行列Aを入力とし，Neville eliminationを用いて，3重対角行列への相似変換を行う。ただし，入力行列が特異であるため，Neville eliminationの過程でピボットが0になることが生じうる。このとき，別の非ゼロ要素をピボットとして用いると，TN性が崩れてしまう（消去に用いる行列が2重対角でなくなる）。そこで，ピボットが0になったら，行列の行と列を削除するデフレーションを行い，固有値0を除去する。この操作は，TN性を保つことが示される。こうして3重対角行列に変換した後，dqds法などで（非ゼロの）固有値を計算する。同時に，デフレーションの回数を数えることで，Aのランクも計算できる。
一方，固有値0に対するジョルダン細胞の数は，n-rank(A)で計算できる。また，次数が2以上のジョルダン細胞の数は，rank(A)-rank(A2)で計算できる。ここで，A2の計算は，高精度性を保つため，やはりTN行列の2重対角表現を用いて行う。こうして，A3，A4，…を計算していくことで，すべての次数のジョルダン細胞の個数が計算でき，Aのジョルダン標準形が求められる。演算量は最大でO(n4)である。
(17) A QR Algorithm with Generator Compression for Structured Eigenvalue Computation
P. Boito (Univ. of Limoges), Y. Eidelman (Tel-Aviv Univ.,) L. Gemignani (Univ. of Pisa) and I. Gohberg (Tel-Aviv Univ.)
n次代数方程式のすべての解を求める方法として，方程式の係数からコンパニオン行列を生成し，その固有値を求める方法がある。コンパニオン行列は，単位行列の要素を1ずつ巡回的に右にずらした巡回行列と，第n列にだけ非ゼロ要素を持つ行列との和として書くことができ，巡回行列+ランク1行列という構造を持っている。このような行列は，quasi-separable行列（行列の狭義下三角部分あるいは狭義上三角部分から取った部分行列のランクがすべて1以下となる行列）の一種である。
一般に，quasi-separable行列は，generatorと呼ばれるO(n)個の変数を用いて表現できる。また，quasi-separable行列に対するQR法は，quasi-separable構造を保存することが知られている。そこで，本研究では，入力のコンパニオン行列からまずgeneratorを計算し（generator compression），generatorに対する反復式としてQR法の反復式を書き直した。これにより，QR法の1ステップをO(n)のメモリとO(n)の演算量で行うことを可能にした。また，陰的なシングルシフト，ダブルシフトを行うアルゴリズムも可能にした。
数値実験の結果，提案法は通常のQR法とほぼ同様に，平均3反復で1個の固有値を計算できた。したがって，全計算量はO(n2)となる。また，悪条件の問題であっても，後退誤差（数値解を真の解として持つ代数方程式の係数と，元の代数方程式の係数との差）は丸め誤差オーダーであることが示された。1000次の代数方程式のすべての解を求める場合の計算時間は，LAPACKのQR法を使った場合に168秒であるのに対し，本手法では6秒であった。
今後の課題としては，理論的な後退誤差解析，および一般固有値計算のためのQZ法への拡張が挙げられる。
感想： 構造を持つ行列に対する数値解法は近年盛んに研究されている。中でもquasi-separable行列は，コンパニオン行列など応用上重要な行列を含むクラスであり，注目度が高い。本研究では，コンパニオン行列の固有値問題に対し，1反復あたりの計算量がO(n)の陰的ダブルシフトQR法を提案しており，これはn次代数方程式がほぼO(n2)の計算量で解けることを意味する。数値実験の結果，収束性は固有値1個あたり3反復程度と安定しており，後退誤差も小さい。今後，解析が進んでアルゴリズムの後退安定性が証明されたら，古典的なデュラン-ケルナー法に代わって，本手法が高次代数方程式を解くための標準的な手法になるのではないだろうか。
通常のQR法の場合，後退安定性はハウスホルダー変換の後退安定性に基づいて証明される。しかし，generatorを変数に取った場合には，1反復の計算は（これらの変数を並べたベクトルに対する）直交変換にはならないので，後退安定性の証明は自明ではない。また，誤差解析を完成させるには，generatorの微小変化に対する固有値の摂動についても解析する必要がある。一方，一般固有値問題への適用については，quasi-separable行列の拡張という文脈において，右辺の行列をどんな行列に取るのが自然なのかに興味がある。
なお，コンパニオン行列のquasi-separable構造を利用して1反復当たりO(n)の計算量でQR反復を行うアルゴリズムは，BoitoのほかにもWatkinsやEidelmanが発表しており，ホットなテーマになっている。
(18) Fast and Memory Efficient Approximate Sparse Factorizations using Hierarchically Semi-Separable Structures
Xiaoe Li (Lawrence Berkeley National Laboratory)
スパースダイレクトソルバ（疎行列直接解法）は，反復法では解きにくい問題に対する解法として，広く使われているが，演算量が大きいことに加え，演算量に対するプロセッサ間の通信量やメモリアクセス量が大きいという問題がある。実際，3次元の問題で，行列サイズがNの場合，総演算量はO(N2)であり，演算1回あたり対するデータ移動量を表すByte/Flopも，通信に対してO(√p / N4/3)，メモリアクセスに対してO(√p / N2/3)と，密行列系の演算に比べて大きい。
そこで，データ量と演算量を圧縮するため，Hierarchical Semi-Separable（HSS）近似を用いた行列データの格納が研究されている。HSS近似とは，行列を2×2のブロックに分割して2つの非対角ブロックを低階数の行列で近似し，対角ブロックについては再帰的に分割を行うことで行列を近似する方法である。偏微分方程式の離散化から生じる行列に対しては，多くの場合，この近似が有効であることが示されている。低階数近似には，特異値分解を使う方法のほか，より演算量の少ない方法として，Rank-Revealing QR分解を使う方法，ランダムサンプリングを使う方法などが提案されている。m×nの行列をランクkの行列で近似する場合，それぞれの演算量はO(mn2)，O(mnk)，O(mnk)であり，ランダムサンプリングを使う方法は，行列ベクトル積が（高速多重極子展開法などにより）O(mn)より少ない演算量で計算できる場合は，より演算量を減らせることが示されている。
本研究では，マルチフロンタル法に基づく疎行列直接解法に対し，HSS近似による行列の格納方式を採用することで，演算量とデータ量を大きく削減できる可能性があることを示した。また，メモリアクセスに対するByte/FlopをO(p log p / N)に改善できることを示した。本手法を実装したソルバにより，変数が64億個の非等方ヘルムホルツ方程式を4096プロセッサで解いたところ，20分で解を求めることができた。
感想： 疎行列直接法は線形計算の中でも特に演算量が膨大になるアルゴリズムであり，それだけに，HSS近似やランダムサンプリングなど，最先端の手法をいち早く取り入れようとしているように感じた。本発表の次に行われた疎行列直接法ライブラリMUMPSに関する発表でも，違った手法により，やはり低階数近似を取り入れていた。
疎行列直接法は，本来は単なるピボット選択付きLU分解のはずだが，変数のオーダリングやフロンタル行列の抽出／依存関係解析に高度なグラフ理論的手法を用いたり，消去につれて変化する非ゼロ要素位置を少ないメモリで効率的に管理するために高度な動的データ構造を用いたり，消去の段階ごとに異なる並列化手法を採用したり，今回のように演算量削減のために階層的な低階数近似を取り入れたりと，多くの技法が組み合わされた大規模な総合技術になっている。そのため，専門家以外にはソルバの全貌を理解するのが難しくなりつつあり，残念ながら新規参入もしにくくなっている。
しかし一方で，疎行列直接法でなければ解けない問題はまだ多く，最近では，反復法の前処理としての利用など，活用の場所も広がっている。その際，ブラックボックスとしてでなく，やはり原理をある程度理解した上で使うことが望ましい。疎行列直接法について比較的最近の進展まで含めて書かれた本として，UMFPACKの開発者であるTimothy Davisの本[18-1]があるが，この本ではほとんど正定値対称の場合しか扱っていない。疎行列直接法の難しさは非対称の場合にあるので，そこを正面から扱った本が出て欲しいと思う。
参考文献
[18-1] Timothy A. Davis: "Direct Methods for Sparse Linear Systems", SIAM, Philadelphia, 2006.
４. 聴けなかった講演
今回は多数のセッションがパラレルで行われたため，聴きたかったのに聴けなかった講演が多かった。そのうちのいくつかについて，自分用のメモとして，アブストラクト集を基に内容をまとめておく。
(19) Efficient, Communication-Minimizing Algorithms for the Symmetric Eigenvalue Problem and the Singular Value Decomposition
Yuji Nakatsukasa and Nicholas J. Higham (University of Manchester)
極分解を用いて，漸近的に最適な通信量で対称行列の固有値分解および特異値分解を行う方法。演算量は標準的なアルゴリズムの3倍以内。
(20) LU Factorization with Panel Rank Revealing Pivoting and Its Communiation Avoiding Version
Amal Khabou (INRIA Sacley), James Demmel (UC Berkeley), Laura Grigori (INRIA Sacley) and Ming Gu (UC Berkeley)
ブロックLU分解において，パネルの分解をRank Revealing QR分解を用いて行うアルゴリズム。通常のピボット選択付きLU分解では解きにくい悪条件の問題に強い。
(21) 2.5D Algorithms for Parallel Dense Linear Algebra
Edger Solomonik and James Demmel (UC Berkeley)
余っているメモリを有効に利用することで，データ移動量を漸近的に削減する方法。行列乗算，LU分解，コレスキー分解，コレスキー-QR分解などに適用可能。
(22) Random Shadow Vectors in IDR(s): An Explanation of Its GMRES-like Convergence
Peter Sonneveld (The Delft University of technology)
IDR(s)法において，sを大きくしたときの収束特性がGMRES法に似てくることの理論的説明。
(23) Highly Accurate Numerical Linear Algebra via Rank Revealing Decomposition
Frolian M. Dopico (Univ. Carlos III, Spain)
Highly accurate algorithmとは，悪条件の問題に対しても丸め誤差レベルの前進誤差を持つアルゴリズムのこと。Rank revealing decompositionとは，A=XDY（ただしXとYは良条件の行列，Dは正則な対角行列）という分解であり，このような分解が可能であるとき，highly accurate algorithmを構築できる場合が多い。このような分解が可能な行列のクラスとしては，Vandermonde行列，Cauchy行列，graded行列などがあり，A自体ではなく分解の要素X，D，Yに対して計算を行うことで，特異値分解や対称固有値分解に対するhighly accurate algorithmが構築されてきた。ごく最近，連立1次方程式の解法や最小2乗問題に対してもアルゴリズムが提案されている。
感想： 近年，様々な構造を持つ行列に対して，高い相対精度で解を求める方法が提案されている。本講演は，これらをrank revealing decompositionという観点から統一的に扱おうとしているようで，非常に興味深い。テンソル分解のセッション（実は玉石混交だった）ではなくこちらに出るべきだった。Dopico教授は，上記の行列のほか，quasi-separable行列に対する解法も手掛けており，構造を持つ行列に対する数値解法について精力的に研究している。また，同じスペインのPena教授も，P行列（首座小行列がすべて正の行列）とその様々な部分クラスに対する高精度な数値解法を研究しており，スペインでは構造を持つ行列に対する解法の研究が盛んなようである。
(24) Stability of Numerical Algorithms with Quasiseparable Matrices
Pavel Zhlobich (Univ. of Edinburgh), Frolian M. Dopico (Univ. Carlos III, Spain) and Vadim Olshevsky (Univ. of Connecticut)
近年，quasi-separable行列に対する研究が盛んであり，QR分解，連立1次方程式の求解，QR反復などに対してO(n)のアルゴリズムが提案されている。しかし，その数値的安定性については十分な研究が行われていない。本発表では，quasi-separable行列に対するQR分解および連立1次方程式の解法について，初めて誤差解析の結果を与える。また，本発表では，dqdsアルゴリズムの一般化と，代数方程式の求解への応用についても触れる。
(25) Inclusion Theorems for Pseudospectra of Block Triangular Matrices
Michael Karow (TU Berlin)
行列Aのε-pseudospectrum σε(A)とは，行列Eが条件‖E‖2 < εを満たしつつあらゆる行列を動くときの，A+Eのスペクトルの和集合のことである。本発表では，Aが2×2のブロック行列であるとき，σε(A)がσf(ε)(A11) ∪σf(ε)(A22) とσg(ε)(A11) ∪σg(ε)(A22)によってバウンドされる（集合の間に包含関係が成り立つ）ことを示す。ここで，f(ε)とg(ε)はブロックA11，A12，A22のある関数である。
(26) The Polynomial Root Finding Problems and Quasiseparable Representations of Unitary Matrices
Yuli Eidelman (Tel-Aviv University)
コンパニオン行列はユニタリ行列 + ランク1行列という構造を持ち，かつ，ユニタリ行列部分はquasi-separableである。この構造を利用して，コンパニオン行列に対するQR法が1ステップ当たりO(n)の演算量で実行できることを示す。
(27) The Rotation of Eigenspaces of Perturbed Matrix Pairs
Ninoslav Truhar (University J. J. Strossmayer), Luka Grubisic (University of Zagreb), Suzaa Miodragovic (University J. J. Strossmayer) and Kresimir Veselic (Fernuniversitat in Hagen)
行列束A-λB（Aはエルミート行列，Bはエルミート正定値）の固有空間の相対的摂動に関する新しいsinΘ定理。
(28) Fast Computation of Zeros of a Polynomial
D. S. Watkins, J. L. Aurentz (Washington State Univ.) and Raf Vandebril (K. U. Leuven)
コンパニオン行列の構造を利用して，QR法の1反復をO(n)の演算量で行うアルゴリズムの提案。コンパニオン行列は，実質的に2×2の行列のn-1個の積として表現できる（Fiedler分解）。この構造を利用するQR法の変種を提案した。所要メモリ量は4nである。
感想： Watkinsの方法は，コンパニオン行列の分解に基づく方法なので，行列に対してどんな操作を行っているかがわかりやすい。それに対して，generatorを変数に取る方法は，漸化式を見ても，行列として見たときにどんな操作を行っているのかが理解しにくい。その意味で，通常のQR法のアルゴリズムに馴染んでいる人にとっては，Watkinsの方法のほうが理解しやすいのではないかと思う。
５. 雑感