１.　学会の概要

本会議は SIAM（米国応用数理学会）主催で3年に1回開催される線形代数の応用に関する国際会議である。通常は米国内で行われるが，今回は GAMM（ドイツ応用数学・力学会）との共催でドイツのデュッセルドルフで行われ，約 230 件の発表があった。この会議は線形計算の分野で第一線の研究者が多数参加し，最先端の研究内容が報告されるのが特徴である。今年も Demmel，Higham，van der Vorst，Saad，Trefethen，Hackbush など，著名な研究者が多数参加していた。

２.　目立ったテーマ

今回の会議では，(i) 行列の指数関数やp 乗根など行列の関数の効率的な計算法，(ii) 鞍点型の係数行列を持つ連立一次方程式の解法，特に反復法のための前処理法，(iii) 大規模な線形システムの簡約化手法，といったテーマが目立ち，それぞれ10件以以上の発表があった。(i)では，硬い連立常微分方程式を解くための exponential integrator で用いる行列の指数関数の計算法，複素積分を用いて各種の行列関数を統一的に計算するための手法などが報告されていた。(ii)では，ナビエ・ストークス方程式や制約付き最適化問題の KKT 条件など，鞍点型の係数行列が現れる典型的な問題に対し，問題の性質を利用して前処理行列を構築したり，その収束性を理論的に解析する発表が目立った。

この他に興味深かった発表としては，テンソルの分解に関する基調講演，Hierarchical matrix に関する基調講演，今夏公開予定の LAPACK 4.0 の紹介などがあった。また，今回は新しい試みとして線形計算の歴史に関するセッションが2個（計9講演）設けられ，Golub，Gutknecht，Strakos，Varga などが自分の研究分野を中心に歴史を振り返る講演を行った。このセッションは大盛況であり，DVD が発売予定となっている。また，2日目午後には，恒例の SIAG/LA Prize の授賞式があり，今回は対称三重対角行列の固有値・固有ベクトル計算のための MRRR アルゴリズムの開発に対して，Indergit S. Dhillon と Beresford N. Parlett が受賞した。

３.　興味深かった講演

特に興味深かった発表について，以下にまとめる。

(1) Tensor Decompositions: An Overview (Plenary talk)
　　Lieven de Lathauwer (CNRS)

3 階以上のテンソルについて，線形代数におけるランクの概念の拡張と，テンソルを低ランクのテンソルの和に分解する手法が紹介された。通常の行列 A は，ランク 1 の行列の和として表すことができる。ここで，（ランク 1 の行列を構成する）横ベクトルどうし，縦ベクトルどうしに直交性の条件を課すと，分解は一意に定まる。これが特異値分解である。特異値分解では，和を大きい方から k 番目の特異値までで打ち切ったとき，ランク k の行列による A の最良近似が得られるという性質がある。

これに対し，テンソルの場合には様子が異なる。与えられたテンソルをランク 1 のテンソルの和として表し，各方向のベクトルが直交系をなすようすることは可能である。これを HOSVD（higher Order Singular Value Decomposition）と呼ぶ。しかし HOSVD では，特異値分解と異なり，和を途中で打ち切った結果は一般に最良近似とはならない。一方，ランク 1 のテンソル k 個の和によって与えられたテンソルの最良近似を求めることができる。この場合，直交性という条件を課さなくても，分解は（適当な条件の下で）一意に定まる。これを PARAFAC と呼ぶ。ただし PARAFAC は，k を変えると最初から作り直す必要があり，特異値分解のように最良近似の性質を保ちつつ逐次的に項を追加していくことはできない。

直交性を課さなくても低ランクのテンソルによる分解が（最良近似の意味で）一意的に定まるという性質は，もともと直交性が本質的でない応用に対してはメリットとなる。この応用例として，ICA（独立成分分析）が挙げられる。ICA では，入力信号のベクトルを線形変換し，互いに独立な信号の和に分解する。このためには，変換後の信号の成分が互いに無相関という条件だけでなく，高次の統計量に対しても条件を課す必要がある。これはテンソルの近似的対角化に相当し，ヤコビ法に類似のアルゴリズムを用いて行うことができる。

感想： テンソルの分解というあまり馴染みのない分野についての紹介であったが，直交性という条件を課さなくても最良の分解が一意的に決まるのが興味深かった。また，「独立成分分析 = テンソルの分解」という見方により，独立成分分析が何をやっているのかがよく分かった気がした。

(2) Refinement Algorithms for Bandwidth Reduction
　　Jennifer Scott & John Reid (Rutherford Appleton Laboratory)

行列の帯幅を圧縮するアルゴリズムは，対称行列に対しては Reverse Cuthill Mckee（RCM）アルゴリズムをはじめとして多数提案されているが，非対称行列についてはあまり提案されていない。本報告では，RCM アルゴリズムを非対称行列に拡張し，さらに Lim らのアルゴリズム [2-1] により局所的な最適化を行う手法を開発した。本手法を，A+A^T に対して対称行列用帯幅圧縮アルゴリズムを適用する方法と比べたところ，多くの問題で顕著な改善が見られた。

参考文献
[2-1] A. Lim, B. Rodrigues and F. Xiao: “a Centroid-based Approach to Solve the Bandwidth Minimization Problem”, Proceedings of the 37th Hawaii International Conference on System Sciences, IEEE (2004).
[2-2] J. Reid and J. Scott: “Reducing the Total Bandwidth of a Sparse Unsymmetric Matrix”, Technical Report RAL-TR-2005-001, Rutherford Appleton Laboratory, 2005.

(3) Diagonalization Algorithms in Real Space Methods for Electronic Structure
　　Yousef Saad (Univ. of Minnesota)

著者らが10年以上に渡り開発してきた電子状態計算プログラム PARSEC の紹介。PARSEC は DFT-LDA（密度汎関数法の局所密度近似）と擬ポテンシャルに基づくプログラムであり，実空間での高精度差分法を用いて計算を行う。数学的には，非線形固有値問題である Kohn-Sham 方程式を，線形固有値問題の反復により解いており，線形固有値問題の求解が計算時間の大部分を占める。

この求解のアルゴリズムとして IRAM（Implicitly Restarted Arnoldi Method），AMLS（Automated Multi Level Substructuring）など様々な手法を試したが，IRAM は多数の電子軌道を求める場合には効率的でなく，AMLS は精度が十分でないという問題があった。そこで現在は，Lanczos 法を用いている。ただし，Lanczos 法では再直交化のコストが大きい。そこで，途中の線形固有値問題においては固有ベクトルの直交性は必ずしも重要でなく，占有軌道に対応する固有ベクトル群の張る空間のみが重要であることを利用し，部分的再直交化を行うに留めている。さらに，必要なエネルギー順位に属する固有ベクトルの成分を増幅するため，Chebyshev 多項式によるフィルタリングを行っている。

Si_525 H_296 の電子状態の計算（行列サイズ 290,000，占有軌道数 1,194）を例に取ると，これらの改良により，1997年時点でのプログラムでは CRAY T3E 48CPU を使って20時間かかった計算が，現在のプログラムでは 1CPU で3時間程度できるようになった。さらに，Si_9041 H_1860 など，巨大な系に対する計算も行っている。今後の課題は，途中の線形固有値問題において，固有ベクトルを求めるのを止め，固有空間のみで計算を行うようアルゴリズムを変えることである。

(4) From qd to LR and QR (historical session)
　　Martin Gutknecht (ETH) and Beresford N. Parlett (Univ. of California Berkeley)

Rutishauser により固有値計算のための qd 法が発見され，さらにLR 法，QR 法の発見へとつながっていく経緯を歴史的に振り返る講演。ただし時間の制約のため，実際には qd 法の発見までで話が終わってしまった。

A を行列，x，y を任意のベクトル，z を複素数のスカラーとするとき，y^T (zI-A)^{-1} x をローラン級数に展開すると，その係数に A の固有値のモーメントμ_νが出てくることは良く知られている。ここで，μ_{ν+1}/μ_ν は ν → ∞ のとき絶対値最大の固有値 λ_1 に収束することから，このようにして λ_1 の計算が可能となる。Stiefel は1953年ごろ，同様のアイディアでモーメントからすべての固有値を計算する方法を考案してみるよう Rutishauser に勧めた。Rutishauser はやってみたが，この方法では非常に悪条件の問題になってしまうことがわかった。

一方，モーメントμ_ν，μ_{ν+1}，… を順に反対角線に沿って並べてできる k×k の Hankel 行列の行列式を H_k^νとすると，|λ_{k+1}| < |λ_k|ならば，H_k^{ν+1}/H_k^ν はλ_1 λ_2 … λ_k に収束することが Hadamard（1892年）によって示されていた。さらに，H_k^νは2次同次式の形の漸化式を満たすため，H_k^0 がわかれば H_k^νは順に計算できる。さらに Rutishauser は，H_k^ν の比の差として定義される変数 q_k^ν，e_k^νを新たに定義し，漸化式をこれらに対する漸化式として書き直した。こうして得られたのが qd 法である。

感想： 行列のレゾルベントからモーメントが導かれ，そのモーメントを並べてできる Hankel 行列式の漸化式から qd 法が自然に導かれるところが非常に印象深かった。このあたりは新旧さまざまな固有値解法の基礎となっている部分であり，既に知られている結果を利用して新しい解法の安定性を解析することなどもできそうである。勉強しておく必要があると思った。

(5) Two Preconditioners for the Saddle Point Problems
　　Arie de Niet (Univ. of Groningen)

非圧縮流体の方程式を差分法あるいは有限要素法で離散化すると，鞍点型の係数行列（2×2のブロック構造を持ち，(1,1) ブロックが対称正定値，(2,1) ブロックが (1,2) ブロックの転置，(2,2) ブロックが 0 の行列）を持つ連立一次方程式が生じる。これに対する前処理として，Silvester & Wathen による SIMPLE(R) 前処理，Kay，Loghin & Wathen による前処理（KLW）などが提案されているが，著者らは新しい前処理方式として，grad-div 前処理と artificial compressibility 前処理の2つを提案した。

これらの前処理を施した係数行列はすべて，高い多重度の固有値1とそれ以外の固有値を持つが，固有値1に対する収束はクリロフ部分空間法では1反復で終了するため，収束性にほとんど影響を与えない。そこで，1 以外の固有値のみを使って実効的条件数を定義し，各前処理の比較を理論的に行った。その結果，提案する2つの前処理は，SIMPLE(R) や KLW に比べて少なくとも同程度に良いという結果を得た。Stokes 方程式，Oseen 方程式，海洋流の方程式の3つに対して数値実験を行ったところ，提案する前処理は反復回数でも実行時間でも SIMPLE(R) と KLW に比べて優ることが確認できた。また，反復回数は問題サイズによらず数回程度という結果であった。

感想： 鞍点型問題に対する解法の発表はこれ1つしか聞かなかったが，他にも多くの発表があり，解法の開発と解析が急速に進んでいるという印象を受けた。grad-div 前処理は逆行列をかける処理を含んでおり，そのためこの部分でまた前処理付き反復法が必要になるとのことで，面白いと思った。

(6) he Use of Total Least Squares Data Fitting in the Shape from Moments Problem
　　Mieke Schuermans (Katholieke Universiteit Leuven)

平面上に多角形があるとき，その頂点の位置をモーメントから推定する問題は，通常，モーメントを要素とする Hankel 行列 H，H^{<} を係数とする一般固有値問題 H^{<}x = λHx の求解に帰着される。しかし，一般にモーメントは誤差を含むため，その影響を低減することが必要となる。

そこで本研究では，係数行列の誤差まで考慮する total least squares problem として問題を定式化し，さらに，係数行列の誤差は Hankel 構造を保存する形で生じると仮定して問題を解いた。このような問題の解法としては，反復法が提案されている。モーメントに統計的誤差を加えた問題で数値実験を行った結果，本手法により解の精度が向上することが示された。

感想： Hankel 行列を係数とする一般固有値問題の扱いは，私も類似の問題に取り組んでいるところだったので，興味深く聞いた。モーメントの誤差が一般固有値問題の解にどの程度の影響を及ぼすかは，一般固有値問題の条件数に依存する。条件数が小さいならば本発表の手法が役立つと思うが，多角形の隣接する辺が（反）平行に近づくと，H と H^{<} が同じ方向の null space を持つようになって非正則な一般固有値問題に近づき，モーメントの誤差の影響が急速に拡大する。これを解決するには，QZ 法を適用する前に予め null space を同定し，デフレーションするなどの方法が必要ではないかと思う。

(7) Iterations for Matrix Roots: Convergence and Stability (Plenary talk)
　　Nicholas J. Higham (University of Manchester)

行列 A の p 乗根を求める問題は，独立成分分析の前処理や，マルコフ連鎖の確率推移行列を観測データから推定する問題など，多くの応用がある。p 乗根の計算手法としては，行列を Schur 形（対称行列ならば対角化）に変形し，上三角行列の p 乗根を漸化式により計算する方法 [7-1]，ニュートン法に基づく反復法，複素積分を用いる方法 [7-3] などが提案されてきた。

p>2 の場合，従来はニュートン法は収束性の問題によりあまり使われていなかったが，最近，新たな公式がいくつか開発されると共に，収束性の解析が進み，元の行列の固有値が複素平面上のある凸領域にすべて含まれるならば，収束が保証されることがわかってきた。さらに，行列を Schur 形に変換してからニュートン法を適用することで計算量を削減する Schur-Newton 法も提案されている。 p 乗根の応用として，企業の1ヶ月間における信用度の変化を表す確率推移行列 Pを，3ヶ月間における確率推移行列 P^3 から求める問題が挙げられる。この場合，求めた P は要素が非負で行の和がすべて 1 でなければならないという条件がある。ニュートン法では，P^3 の行の和がすべて 1 ならば，Pも自動的にその性質を受け継ぐことが示される。ただし，要素が非負という条件は自動的には満たされないため，得られた結果を修正することが必要である [7-5][7-6]。

参考文献
[7-1] Smith (2003)
[7-2] Bini, Higham and Meini (2005)
[7-3] Iannazzo (2006)
[7-4] Guo and Higham (2006)
[7-5] Israel et al. “Finding Generators for Markov Chains”, Mathematical Finance.
[7-6] Kreinin, Algo Research Quarterly.

感想： 対称正定値行列の平方根の計算が現在仕事上で必要になっているので，興味深く聞いた。標準的な解法は対角化を経由する方法であるが，通常の対角化では元の行列が持っている望ましい性質（単位行列に近いなど）を利用できず，アルゴリズムも複雑であり，また，固有ベクトルの直交性保証などの面倒な問題がある。これに対してニュートン法系統の反復法は，行列の乗算・加算だけで済むので極めて簡単であり，キャッシュの利用効率も高い。両者のアルゴリズムを実装し，様々な条件で性能を比較してみたい。

(8) Orthogonal Eigenvectors & Gram-Schmidt
　　Inderjit S. Dhillon (Univ. of Texas at Austin) and Beresford N. Parlett (Univ. of California Berkeley)

今年度の SIAG/LA Prize の受賞記念講演であり，対称三重対角行列の固有値・固有ベクトル計算のための MRRR アルゴリズムの基本的なアイディアと現状が紹介された。MRRR アルゴリズムの鍵となるのは，次の4つのアイディアである。

　(i) 三重対角行列Tを，修正コレスキー分解 LDL^T の形で保持し，計算を行うこと。
　(ii) シフトを行う際は，三重対角行列を経由せず，シフト後の行列に対する分解 L'D'L'^T を LDL^T から differential qd 変換により直接計算すること。
　(iii) 固有ベクトルを計算する際には，ツイスト分解 NDN^T を使うこと。ツイスト分解は，L'D'L'^Tから differential qd 変換により直接計算すること。
　(iv) representation tree を使い，固有値の階層的なクラスターから固有値を 1 個ずつ分離し，計算を行うこと。

(i)が望ましいのは，T の要素が微小に変動すると固有値は大きな相対的変動を被ってしまうが，L または D の要素が微小に変動しても，固有値の相対的変動は小さい場合が多いからである。このような性質を持つ行列の表現を，Relatively Robust Representation（RRR）と呼ぶ。修正コレスキー分解は，正定値の場合は RRR となることが保証され，不定値の場合でも，（固有ベクトルの計算で必要となる）絶対値最小の固有値に対しては，ほとんどの場合に RRR となる。これは，コレスキー分解の際に大きな element growth があっても成り立つ。この RRR の考え方をベースとして，あとは計算の過程で L，D の要素に小さい相対誤差しか入らないことを保証するよう，differential qd 変換を用いてアルゴリズムが構成される。

MRRR アルゴリズムの現状としては，分割統治法との性能比較と，極めて密集した固有値を持つ問題への対処法が示された。MRRR は分割統治法に比べ，演算量は小さいが，演算がすべて BLAS1 レベルという欠点がある。そのため，演算の大部分を BLAS3 で実行できる分割統治法と比べると，1万元程度の問題に対し，実行時間では接戦となるとのことであった。また，glued Wilkinson matrix など極めて密集した固有値を持つ問題に対しては，シフトを行っても固有値を数値的に分離できない場合があるが，これを解決するには，元の表現 LDL^T に微小な摂動を加えることが有効とのことであった。

感想： MRRR は対称三重対角行列の固有ベクトルの計算において直交化を（ほとんどの場合に）不要とし，演算量の大幅な削減と並列性の増大を達成したアルゴリズムであり，relatively robust representation，混合型誤差解析などの新しい概念も生み出していることから，受賞は極めて妥当だと思う。講演では，実際の数値例を見せながら上の4つのアイディアがなぜ必要かを順に説明していったが，論理展開は明快で，非常に納得できるものであった。特に，ツイスト分解を用いて逆反復法における初期ベクトルを最適に設定し，そのときの残差がε|λ| となることを示す部分は，理論による予言と数値結果とがぴたりと一致し，鮮やかであった。ただ，MRRR アルゴリズムには，極めて密集した固有値に対する処理以外にも，RRR の存在問題などいくつかの課題が残されているので，それらをまとめて未解決問題として提示してくれればもっとよかったと思う。

(9) Stability of Group-Theoretic Algorithms for Fast Matrix Multiplication (Plenary talk)
　　James Demmel (Univ. of California Berkeley)

n×n の正方行列どうしの乗算を O(n^3) より少ない演算回数で行うアルゴリズムは数多く発表されているが，小さいべきを持つアルゴリズムは定数が極めて大きく，通常使われるサイズの n に対しては実用的ではないとされている。

これに対し，最近，群論を用いたタイプの高速乗算アルゴリズムが提案されている。これは，入力行列 A，B をある群に埋め込み，群上での FFT を用いて乗算を行い，その結果を取り出すという 3 ステップからなるアルゴリズムである [9-1]。乗算の結果を取り出せるためには，群に特別な条件が必要となるが，最近，この条件を満たす群が発見され，実際にアルゴリズムの構築が可能となった [9-2]。

本研究では，このタイプの高速乗算アルゴリズムに対し，誤差解析を行った。その結果，計算結果の行列と真の結果との誤差のノルムは，c×n^d×ε×‖A‖‖B‖（c，d は定数，εはマシンイプシロン）で抑えられ，これらのアルゴリズムは後退誤差解析の意味で安定であることを示すことができた。

参考文献
[9-1] H. Cohn and C. Umans: “A Group-Theoretic Approach to Matrix Multiplication”, in Foundations of Computer Science. 44th Annual IEEE Symposium, pp. 438-449, 2003.
[9-2] H. Cohn, R. Kleinberg, B. Szegedy and C. Umans: “Group-Theoretic Algorithms for Matrix Multiplication”, in Foundations of Computer Science. 46th Annual IEEE Symposium, pp. 379-388, 2005.

(10) The Future of LAPACK and ScaLAPACK
　　James Demmel (Univ. of California Berkeley) and Sven Hammarling (NAG Ltd.)

線形計算ライブラリ (Sca)LAPACK の現状と今後の方向が，今夏と来夏の2回に分けて公開される LAPACK 4.0 の話題を中心に紹介された。今後のハードウェアは，並列性の向上（256コアのプロセッサなど）と，プロセッサ性能とメモリ性能とのギャップの指数的増大とが大きなトレンドとなる。そのため， LAPACK ではより良いアルゴリズムを開発と自動チューニングの採用によりこれらの動向に対応するとともに，より精度の高いアルゴリズムの開発，およびメニューの増加にも努めていくとのことであった。 LAPACK 4.0 で新規に提供されるルーチンとしては次が挙げられた（他にもあったがメモしたのは以下のみ）。

(i) 改良版の MRRR アルゴリズムによる対称三重対角行列の固有値・固有ベクトル計算。広く使われている分割統治法に比べ，演算量が大幅に削減される。ただし，分割統治法は BLAS3 を用いているため，1万元程度の問題では実行時間は同程度。また，MRRR アルゴリズムでは二分法の計算がネックとなるため，この部分を GPU を用いて高速化する研究も（LAPACK 4.0 とは別に）進行中である。Volkov の実装によると，現在，CPU のみを使う場合と比べて 100 倍の高速化が得られている。

(ii) 特別の高精度 BLAS を用いた連立一次方程式の解の反復改良アルゴリズム。条件数κが n^{1/2}ε以下であれば，正しい解が得られることを保証する。

(iii) Drmac と Vecelic の開発した新しいヤコビ法による特異値分解アルゴリズム。QR 法よりも高速で高精度。

(iv) Braman，Byers，Mathias の新しいマルチシフト QR 法による非対称ヘッセンベルグ行列の固有値計算アルゴリズム（前回の SIAG/LA Prize を受賞）。BLAS3 の利用と aggressive early deflation の導入により，従来法の 10 倍程度高速。また，ヘッセンベルグ行列への変換も，van de Geijn らによる BLAS3 の割合を増加したアルゴリズムにより，十数%の性能向上が図られる。

(v) 新しいマルチシフト QZ 法による非対称の一般固有値解法。

(vi) 対称非負定値行列のためのピボット選択付きコレスキー分解。

(vii) 一般化特異値問題の新しいルーチンと2次固有値問題のためのルーチン。

これらに加え，進化するアーキテクチャに適応するための自動チューニング機能が重要な目標として挙げられ，具体例として，ScaLAPACK におけるデータ分割の自動最適化などが挙げられていた。また，Cell や GPU など，単精度の計算であれば通常の CPU の数十倍の性能を達成できる専用プロセッサを効率的に利用するアルゴリズムの開発も今後の研究テーマとして挙げられていた。たとえば連立一次方程式の直接解法では，LU 分解を単精度で行い，反復改良のみを倍精度で行えば倍精度の解が得られるが， Cell での評価結果では，このアルゴリズムは全てを倍精度で計算する通常のアルゴリズムの 10 倍の性能を達成できたということである [10-1]。また，データ格納形式に関しては，対称行列を1次元配列に格納する従来法の代わりに，2×2 のブロックに分割して長方形の配列に詰め込む RFP（Rectangular Fully Packed）format が提案されていた。従来法に比べ，キャッシュを意識した格納形式になっているとのことである。

なお，LAPACK 4.0 では外部の研究者・ユーザからの様々な形での貢献を歓迎するので，詳しくは http://www.netlib.org/lapack-dev にある Contributor's Guide を参照して欲しいとのことであった。

参考文献
[10-1] J. Langou, J. Langou, P. Luszczek, J. Kurzak, A. Buttari and J. Dongarra: “Exploiting the Performance of 32 Bit Floating Point Arithmetic in Obtaining 64 Bit Accuracy (Revisiting Iterative Refinement for linear Systems)”, LAPACK Working Notes, No. 175, 2006.

(11) The Technique of Hierarchical Matrices (Plenary talk)
　　Wolfgang Hackbush (Max-Planck Institute)

n×n の行列に対する加算，乗算，逆転などの演算は，密行列であれば O(n^2) の領域と O(n^2) あるいは O(n^3) の計算量が必要であるが，近似を導入することで，これらをすべて O(n (log n)^q)（q は定数）で行おうというのが hierarchical matrix（以下 H 行列）の考え方である。近似は，行列をブロックに分割し，各ブロックに対して低ランクの行列による近似を行うか，あるいは再帰的にH行列による近似を行うことで達成される。

最も簡単な例としては，行列を 2×2 のブロックに分割して非対角ブロックを低ランクの行列で近似し，対角ブロックを再帰的に H 行列で近似する方法が考えられる。低ランクの行列による近似は特異値分解を用いて行う。この場合，2つの H 行列の和は O(n log n)，積は O(n (log n)^2) で計算できることが示される（ただし H 行列はこれらの演算について閉じていないため，結果を再びH行列で近似することが必要となる）。また，逆行列も，2×2 のブロック行列に対する逆行列の公式を使えば，元のブロックと対角ブロックの逆行列との和・差・積で表すことができる。ここで，対角ブロックの逆行列を再帰的に同じ方法で計算すれば，最終的に逆行列も O(n (log n)^2) で計算できることが示される。

一般には，行列の非ゼロ要素をすべて覆うように行列をいくつかのブロックの和に分割し，各ブロックを低ランクの行列あるいは H 行列で近似する。この場合にも，適当な条件の下で，上記の演算量が達成できる。近似による誤差は，楕円型方程式を境界要素法あるいは有限要素法で離散化してできる行列の場合，指数関数的に小さくできることが示される。

応用としては，反復法における前処理としての利用，シルベスター方程式や代数的リッカチ方程式など行列方程式の解法への利用，（複素積分を通じて）行列の指数関数の計算への応用，行列を陽的に生成することが不可能な超高次元空間での問題から生じる連立一次方程式への応用などが挙げられた。特に超高次元空間の問題は，行列のテンソル積構造を利用することで，1024^2048 元の問題を数分程度で解くことも可能とのことであった。

参考文献
[11-1] http://www.mis.mpg.de/scicomp
[11-2] http://www.hmatrix.org/index.html

感想： 3年前の同じ会議では，制御可能な誤差を許容することで，行列とベクトルの積や連立一次方程式の求解を O(n (log n)^q) の演算量で行うという発表が目立ったが，今回は同じ方向の研究が hierarchical matrix という形で統一されたように思う。H 行列に関する発表は本発表を含めて全部で 5 件あった。また，高速多重極子展開法も，行列をブロックに分け，各ブロックを低ランクの行列で近似して行列ベクトル積の演算量を削減する手法であるから，H 行列はこの発展形と見ることもできる。発表の中では，超高次元空間の問題への応用が特に興味深かった。たとえば物理における 3 次元空間での m 点相関関数の問題は 3m 次元空間での問題となり，一方向を N 点で離散化すると，行列サイズは N^{3m} となる。このような問題に H 行列の考え方が使えれば面白いと思う。

(12) Multishift Variants of the QZ Algorithm with Aggressive Early Deflation
　　Bo Kagstrom and Daniel Kressner (Umea University)

LAPACK 4.0 にも採用予定の Braman，Byers，Mathias の新しいマルチシフト QR 法の QZ 法への拡張。前者において高速化の主な要因となっている BLAS3 の利用と aggressive early deflation は，QZ 法にも自然に拡張できる。しかし，QZ 法特有の問題として，無限固有値の問題がある。無限固有値は，QZ 法を無限精度で行うならば，行列束 S-λT（S はヘッセンベルグ行列，T は上三角行列）における T の対角要素が 0 になることで検出できる。この場合，デフレーションを行うことで，無限固有値を分離できる。しかし，有限精度では，無限固有値が存在しても T の対角要素が十分小さくならない場合があり，このとき，デフレーションが行われないままに計算が進行し，計算誤差のために無限固有値が有限固有値に化けてしまうことがある。この問題の解決が大きな課題であり，本研究では，QZ 法の適用前に行列束 A-λB を GUPTRI（Generalized UPper TRIangular）標準形に変換し，行列束のクロネッカー構造を求めることで，事前に無限固有値のデフレーションを行っている。

本アルゴリズムでは，マルチシフトのシフト数 m，バルジ 1 個あたりが含むシフトの数 n_s，デフレーション・ウィンドウのサイズ n_w の3つののパラメータがあるが，これらを全数探索により最適化し，従来の QZ 法と性能を比較したところ，多くの問題で数倍の性能が得られた。

また，分散メモリ型計算機上での並列化も既に試みており，複数個（m/n_s 個）のバルジからなる列をさらに複数個導入し，パイプライン処理を行うことで，効率的に並列化ができる見込みを得た。

参考文献
[12-1] B. Kagstrom and D. Kressner: “Multishift Variants of the QZ Algorithm with Aggressive Early Deflation”, LAPACK Working Notes, No. 173, 2005.
[12-2] B. Adlerborn, B. Kagstrom and D. Kressner: “Parallel Variants of the Multishift QZ Algorithm with Advanced Deflation Techniques”, PARA06, 2006.

感想： QR 法の画期的な改良版とされている Braman らのアルゴリズムの拡張に関する報告であり，大変興味深かった。特に，無限固有値のデフレーション方法に関する考察は，モーメントに関する一般固有値問題における 0/0 の形の固有値の処理に通じる部分があり，参考になった。

本手法は（Bramanらによるマルチシフト QR 法と同様に）3 つの性能パラメータを含むため，今回の実験ではその最適値を全数探索により決めている。しかし，これはコストが大きいため，今後は実験計画法を用いた探索範囲の限定やモデルによる性能予測など，自動チューニング手法の適用が必要になると思われる。特に，デフレーション・ウィンドウのサイズ n_w については，問題により最適値のばらつきが大きく，統計的手法が必要になる。このように，本問題は自動チューニングの様々な手法を検証するための例題としても面白いと思う。

分散メモリ型計算機上での並列化に関しては，Braman らの QR 法についてはまだ報告されていないので，QZ 法が先行していると思われる。本研究のアプローチでは，単一プロセッサの場合のバルジ列をさらに複数導入することで並列性を抽出しているとのことだが，この場合，シフト数 m が非常に大きくなり（単一プロセッサの場合でも 100 以上が普通），シフトの計算に非常な時間がかかるように思われる。また，シフトの計算を行っている間は，パイプラインが停止してしまうように思われる。これらの点をどう解決しているか，調べてみたい。

(13) Accurate Eigenvalues of Sign-Regular Matrices
　　Plamen Koev (MIT)

Sign regular matrix とは，すべての k 次小行列式が，k が同じなら同じ符号を持つような行列である。ここでは特に，k 次小行列式の符号が (-1)^{k(k-1)/2} の場合を考える。このような行列は，totally positive matrix（すべての小行列式が正である行列）の列の順番を反転させたものになっている。

Koev は，前回のこの会議で，totally positive matrix に対して倍精度演算のみを使って任意の固有値λを誤差 O(ε)λで計算するアルゴリズムを発表した [13-1]。これは，そのような行列が要素が非負の二重対角行列の積として表現できること，二重対角行列の特異値が相対誤差の意味で高精度に計算できることに基づいている。対象とする行列には，ヒルベルト行列など極めて条件数の大きい行列も含まれているので，これは驚くべき結果である。今回の発表では，sign regular matrix が totally positive matrix の列の順番を反転させたものになっていることを利用し，前者に対しても同様に高精度な固有値計算ができることを示した。

なお，Koev は，本発表とは別に，楕円型偏微分方程式を四面体要素を用いた有限要素法により離散化してできる行列の固有値に対しても，倍精度演算のみを用いて小さい相対誤差で固有値を計算するアルゴリズムについて発表をしていた　[13-2]　が，残念ながらこちらは聞き漏らした。

参考文献
[13-1] P Koev: “Accurate Eigenvalues and SVDs of Totally Nonnegative Matrices”, SIAM J. Matrix Anal. Appl., Vol. 27, pp. 1-23 (2005).
[13-2] J. Demmel, P. Koev and S. Vavasis: “Accurate Eigenvalues for Elliptic Problems using Finite Elements”, in preparation.

(14) New Shifted Algorithm with Ensured Convergence to Singular Values of Bidiagonal Matrices
　　Yoshimasa Nakamura and Masashi Iwasaki (Kyoto Univ.)

離散 Lotka-Volterra 系に基づく新しい特異値計算アルゴリズム mdLVs と，特異ベクトル計算のための可積分系を用いた新しいコレスキー分解アルゴリズムの紹介。mdLVs では，qd 法と同様にシフトを用いて収束を加速するが，qd 法と異なり，適切なシフトを選べば，特異値への収束が理論的に保証されるという特長を持つ。また，最終的には 3 次収束をする。また，特異ベクトル計算のためのコレスキー分解についても，2 種類のシフトを導入することで，悪条件の行列に対するコレスキー分解が回避でき，高精度な固有ベクトルが計算できる。

数値実験では，多くの例題に対して，mdLVs が qd 法よりも高精度に特異値を計算できることが示された。また，特異ベクトルの計算においても，2 種類のシフトを適切に定めることにより，直交性を向上できることが示された。

感想： 日本発の独創的なアルゴリズムで，かつ理論的にも様々な特長を有するmdLVsに関する発表であり，非常に興味深かった。ただ，発表が最終日の最終時間帯ということで既に帰ってしまっていた参加者が多かったことと，発表内容が盛りだくさんで質問時間が取れなかったことが残念であった。

４.　写真

学会会場の University of Dusseldorf

フランク・ゲーリーのデザインしたビルなど，モダンな建築が多いライン川沿いのメディアセンター

近隣の町ケルンの大聖堂にて