今倉暁 | 数理人

スライドショーには JavaScript が必要です。

講演者：今倉暁（筑波大学）
題　目：周回積分型固有値解法とその並列ソフトウェア　
日　時：2014年11月24日(月,祝) 18:00–
場　所：早稲田大学西早稲田キャンパス102号室

オープニング

司会：お忙しい中お集まりいただき、ありがとうございます。数理人（すうりすと）セミナーというのを立ち上げて、若手で数理を盛り上げていこうということで、今年の3月ぐらいにみんなで一回集まって、それでその後ずっと頓挫してたんですが、また秋に始まったという感じです。
　この後、講演をしていただくわけなんですけれども、その講演は後でWeb上で公開する予定で、それは今ここにマイクがありますけれども、このマイクで拾って文字起こしをしようと思ってます。質問は適宜してもらって構わないというセミナーにしようと思っているんですが、その辺がまだ全然決まってませんので、その辺も含めてやってみて、悪いところはちょっとずつ改善していこうと思っています。

講演

司会：記念すべき第1回は、筑波大学の今倉先生に講演をお願いいたしました。だいたい75分ぐらい。では、すみませんが、お願いいたします。

今倉：録音するからといってやり方を変えるほど器用ではないので、どうしていいか分からないですけれども。1回目でこういう機会をいただきまして、ありがとうございます。じゃあ、始めます。
　今回のタイトルは『周回積分型固有値解法とその並列ソフトウェア』というタイトルで。今日はアカデミックな話ではあるんですが、僕が普段やっている数理的なところも話すんですが、それよりは僕が所属しているプロジェクトの話を少ししようかなと思ってます。そこでこの固有値解法のベースにしたソフトウェアの開発をやっているので、そこの部分が中心になるか分かりませんが、普段あまり話したことが無かったんで聞いてもらおうかなと思います。

　最初に今日何をやるかというのを今言いましたが、僕が所属しているCRESTプロジェクト、筑波大学の櫻井鉄也先生が代表をされている『ポストペタスケールに対応した階層モデルによる超並列固有値解析エンジンの開発』というタイトルのプロジェクトの紹介がまず一個。あと、その中で使っている固有値解法の概要とソフトウェアの紹介というのをやります。こっち（前者）はそういう漠然とした話で、こちら（後者）はちょっと内容に突っ込んだ話になってます。

　このプロジェクトのもともとの研究背景というのはポストペタスケールを目指そうという話で、ポストペタスケールに対応したソフトウェアを作りましょうという話なんですが。もともとは計算機の環境というのがすごい階層的に最近なってきていますよというのが背景にあります。ノードとかがすごい階層的に深くなってくる。ノードがそもそも独立にあったり、中でコアが複数あったりなんかするような感じになってます。加速器とかが下の方にくっついてたりという感じで、非常に階層的な構造を持っているというのが最近のスーパーコンピュータの特徴です。
　この階層的なアーキテクチャ上で性能を発揮するためには、アルゴリズムレベルでも階層性を持っていた方がいいだろうと。階層性を持ってないアルゴリズムでもうまく実装してやればいいんですが、それよりもアルゴリズムの段階からこの並列化を意識したようなものを扱った方が性能が出るだろうというのが研究背景にあります。

　我々の所属しているプロジェクトでは固有値解法に着目して、そのソフトウェアを作りましょうというのが大きな目標です。ご承知の方も多いかとは思いますが、固有値解法には大きく2つの分類があります。密行列の全固有値を求めるような相似変形を使っていって対角化や上三角化していくようなアルゴリズムと、射影法に基づくような疎行列向けのアルゴリズムの2つがあります。
　ですので、ソフトウェアの段階からそれら2つをどちらだけではなくて両方を開発しようということで、密行列の方がEigenExaというソフトウェアで、疎行列の方がz-Paresというようなソフトウェアというのを作ってます。あとは、これを作る上での並列化も含めた高性能な基盤技術の開発をしてたり、あと、アプリケーションを実際に使う上での利用技術というのを開発してたりというのをやってます。

　プロジェクトの意義は、我々が作っているのはアプリそのものではなくて、その基盤を支えるようなソフトウェア――プロジェクトの中ではエンジンというふうに言ってますが――、このエンジンを高性能なものを作っておく、提供することで、それぞれのアプリにこのエンジンを組み込むといろんなアプリケーションが加速されますよというのがプロジェクトの意義です。だから、これをやっておくといろんな範囲――特定の分野だけじゃなくていろんな分野――で貢献できるというのが基盤技術なので、そういうことになってます。

　ここから、プロジェクトがどういう構成でやられているかと現在までどんなことをやってきたかというのを簡単に。あと、2つのソフトウェアを作っているので、それぞれの概要というのを非常に簡単に。特にこちら（EigenExa）は僕はあまりそこまで詳しくないので、非常に簡単に紹介したいと思います。

　これはプロジェクトの成果報告っぽいスライドなんですが、実際そうなんですが。プロジェクトの進捗は平成24年度ぐらいに固有値解法ソフトウェア2つ、z-ParesとEigenExaというのが既にWeb上で公開されてます。フリーで公開されました。26年度現在、京コンピュータで性能評価というのをやって、高性能化というのをやってきてます。具体的に何をやっているかというと、大規模並列の時にボトルネックになるような部分の改善とか、アプリケーションに組み込んだときの高性能化なんかをやってたりします。
　平成27年度、来年度というのはより高性能化、というのは通信回避とかがやっぱり非常に重要ですよというのでこれをどんどんやっていくとか、高信頼性というのも少しやってます。精度保証までいくわけではないんですが、例えば耐故障性なんかというのも最近やったりしています。あとはアプリにいろいろ組み込んでやりましょうとかいうことです。今後はさらにいろんなアプリに組み込んで性能が出るようにして、インフラとして実現できたらいいなというふうに進めているという段階です。なので一応公開はできているんだけど、これの高性能化を今後もっとやっていかないといけないという段階になります。

　プロジェクトは、僕は櫻井グループに所属しているんですが、そこだけでやっているわけではなくて割と多くのグループが連携してやってます。疎行列向けの固有値解法の開発は、主に筑波大学の櫻井グループと名古屋大学の張グループが連携してやっているというかたちになってます。密行列向けのEigenExaの方は、理研AICSの今村先生の今村グループと電通大に移られた山本有作先生の山本グループが連携してやっているという状況です。あと、実際にアプリに組み込んで性能評価とかもしなきゃいけないので、アプリケーションのグループとして筑波大の藏増先生という方が基礎科学分野、QCDとかの分野のアプリをされている方なのでそこの分野と、あと、鳥取大の星先生が電子状態計算とかをされているグループと連携してやっているというかたちになります。基本的にはこの6グループが連携していてちょっと多めなんですが、やっているという状況になってます。

　主要メンバーはさっきのグループリーダーを含めて、これはたぶん学生を除いた分ぐらいはだいたい入れたと思うんですけど、筑波は多田野さんと僕と二村君。あと、産総研の池上さんと、東大の中務さんが去年ぐらいから入ってるんだと思います。あと、張先生の所は愛知県立大の曽我部さんと名古屋大の宮田さん。あと、今村グループの所に深谷君と廣田君が。山本先生の所は今移られたところで、たぶんいないんだと思います。もともと深谷君とか廣田君はこっち（山本グループ）に所属してました。

　あとはこちらだとあまり知られないかもしれませんが、井町さんはたぶん東大にいらっしゃったんですけど今は鳥取大に移られていますが、この辺のメンバーでやっているという状況です。

　細かいところは後でこちらの方はやりますが、基本的なところだけ、どんなのを作っているかというのをザッと紹介したいと思います。
　先ほども言いましたが、2つソフトウェアを作ってます。エンジンといってますが、ソフトウェアを作ってます。一つが疎行列固有値問題向けのz-Paresで、対象としている問題が次に紹介するEigenExaと違うんですが、一番違うのが疎行列向けというよりも部分固有対を求めるのに向いた射影法に基づく解法です。これは標準固有値問題、一般化固有値問題に対応している。で、実・複素どっちでもいけて、対称・非対称、エルミート・ノンエルミート、どれでもいけますよというような解法になってます。
　部分固有値というのなのでどこを求めるかというと、最大とか最小とかではなくて、ユーザーの指定する任意領域を求めていきますというような仕様になってます。
　基本的には、後でやっていきますが周回積分型固有値解法という解法に基づいてて、解法の名前でいうとBlock SS- RR法というのとHankel法の2種類がパラメータで選べるようになっているかと思います。一部のパラメータについては自動設定の機能がついてたりというような、ソフトウェアの使いやすさとしてはそういうのがついてたりします。この辺は後でもう一回ちゃんとやります。

　このアルゴリズム、ソフトウェアの実装の一番大きな特徴が、アルゴリズムレベルで階層的な構造を持っているということです。後で具体的にはやりますが、この解法のソフトウェアというかアルゴリズムの中で一番計算量の重いのが連立一次方程式、線形方程式を複数個解かなきゃいけない。各積分路に対して16個とかその程度の数の線形方程式を解かなきゃいけない。それがほとんどの計算時間を占めることになります。
　この部分に対してだいたい3階層ぐらいの階層性がある。並列の意味の階層性があるというふうに言ってます。一つは指定領域の固有値を求める解法なので、そもそも領域を独立に置いてあげれば、その2つの領域、複数の領域は全部並列に計算できる。何の通信とかデータのやり取りはなく、完全に独立に解くことができる。こういうのが一階層です。二つ目は、計算主要部の線形方程式の部分に対しては積分点ごとに解くので、ここは完全に独立に解くことができる。これが二つ目。三つ目は、これは完全に独立かというと怪しいんですが、線形方程式自体も並列の解法。厳密な意味で完全独立じゃないけど並列にできるような解法というのはいろいろありますので、それを並列に解くことができる。なので、この三つの階層というのがある。
　最初にも紹介しましたが、研究背景にもなっていたんですが、アーキテクチャの方にも階層性がある。だから、これらの階層に合わせてアーキテクチャの階層とマッチングさせるような感じで実装してあげれば、並列性がすごく高くなります。具体的には、ノード間の通信を比較的ローカルなところに限定してやることができるので、全体の通信が無くなるので性能が出ますねというのがアルゴリズムになってます。

　今、開発の状況がどうなっているか？　僕もこれは詳しくちゃんとフォローできてないんですが、だいたいこんな感じだと思います。2013年12月に初版かベータ版か分かりませんがリリースされてます。これが一番最初。Webで、おそらくこの段階だと京向けのバイナリとして公開されてた。今は完全にコードとして公開されているはずです。バージョンが若干上がってて、これがいつになったらver. 1になるのか僕は知らないんですが、まだやりたくないというか、ベータ版の段階でずっと更新されてます。
　主な成果は、このソフトウェアを使って、実際にはもうちょっと工夫も入っているんですが、企業の中で出てくるような振動解析の問題で非常に大規模な問題に適用してうまくいってますというような結果が出ています。このデータは僕の方で出せないので今日はお見せできないんですが、こういうような成果が出ているということです。
　今後は、実際にはこの解法をアプリに使うためにいろいろこのままでは……動くんですが、性能を出すにはもうちょっとアプリごとの工夫が要るので、その辺をやって、計算カーネル部分、具体的には線形方程式のところをアプリごとの問題でどう解くかというのをちょっと考えないといけないので、その部分のことをやっていって性能を上げていこうというのが今後の展開になると思います。この辺がz-Pares、疎行列向けの固有値解析エンジンの話です。

　もう一つのが密行列向けのEigenExaと呼ばれる、昔から理研の今村先生、昔、電通大にいた今村先生がずっと作られてたやつの最新版のやつでEigenExaというソフトです。これはさっきのと違って、全固有値を基本的には計算するものです。部分的にも計算可能なんですが、ほとんど演算量的には全固有値を求めるのと同じようなものになってしまうので、基本的には全固有値を対象にしています。問題設定は一番シンプルなやつです。実対称の標準固有値問題を求めるというような方法になってます。一番シンプルな問題です。
　今、実対称一般化固有値問題、 $A$ が対称で $B$ が正定値対称なやつには適用できるようにはしていってるんですが、非対称にいくのはそうそう難しいなと。エルミートもなんとかいけるかもしれないけど、非対称は当分いく気はないような感じで今はやってます。
　基本的なアルゴリズムは、基本的に対称行列の固有値問題の全固有値を求めるのに何をやるかというと、大体のやり方、オーソドックスなのはHouseholder法を使って三重対角化して、三重対角行列の固有値をQR法だったり分割統治法だったり何かの方法で解きましょうという方法なんですが、最近は帯行列をいったん作るという方法が実際やられてます。
　このEigenExaはもうちょっと特殊で、帯行列をした後に三重対角化をしないような方法をとります。帯のまま分割統治法を適用してやるというアルゴリズムがベースになってます。他のソフトウェアとはちょっとここが違うところです。

　もうちょっと言いますと、もともとの実対称な行列に対してオーソドックスなのは三重対角化してこれの固有対を求める。逆変換してあげて、固有ベクトルを求めてあげようというような方法をとるんです。ただ、ここの三重対角化の部分はデータ転送量が非常に大きくて、具体的にはBLASのLebel-3が使いにくい、使えない。Lebel-2クラスで実装されてしまうので、いったん帯行列を経由します。ここがLebel-3 BLASで組めて非常に速いんです。ということを使います。その後、三重対角化してあげて戻してあげようという方法をとるというのが、高性能計算の方では従来よくやられていた方法です。
　ただ、ここの逆変換の部分というのが実は高性能実装では非常に難しいらしいです。それを嫌がって、この辺がELPAとかDPLASMAというのがよく使われているソフトですが、ScaLAPACKはこいつで一番オーソドックスな方法をとるんですが、これは性能が出なくてもうちょっと改良したというのがいろんな海外で使われているソフトなんですが、EigenExaは帯にした後、三重対角化を経由せずに帯のまま分割統治法をやってあげる。でやると逆変換を回避できるので、高性能でいけますということをやっています。アルゴリズムレベルで若干違うということになってます。

　これまでの開発状況は、2013年8月にこちらはver. 1というふうにしてちゃんとリリースされてます。12月にver. 2が提供されていて、これはよく分からないんですが、「三重対角化版をやっぱり作ってみました」と言ってました。作って一応提供してます。これは比較のためで、実際は帯化のほうが使えるということになってます。ただ、比較用に一応作りました。
　細かく言いますと、この帯といっているのは実は五重対角に限定して実装されています。だからそんなに帯幅は広くないです。場合によってはもしかしたら七の方がいいケースもあるかもしれないけど、ハードによっては。現状は五で限定して実装されているところです。この五を七に変える実装は非常に難しくて、一から作るぐらい大変らしいので、そうそういじれないみたいです。
　こちらは全固有値を求めるタイプなので成果というと非常に分かりやすくて、どんだけでかい問題を解けたかということに尽きます。EigenExaは昨年に――今年かもしれないけど、いつかは忘れましたが――、100万次元の実対称密行列の対角化を実現しましたというのでプレスリリースされております。たぶん世界初でこれを実現しました。京のフルノードを使って約1時間弱、3,200秒かなんかそれぐらいだったと思いますが、100万次元の固有対が解けています。なので、京を全部使えれば、1時間あれば100万次元全部解けるみたいです。部分的に解くとかいらなくて、全部解く。ごり押しに全部求めればいいというようなことになってます。
　あとはアプリで理研の別の三好先生という方のチームと共同でデータ同化の分野で使われて、これもなんか従来法に比べて非常に速いということでプレスリリースされたそうです。ちょっと具体的には分からないんですが。
　今後は、これはずっと京でされてて京向けにチューニングされているので――今村先生が理研の方なので、そういうのもあるんですが――、もうちょっとBlue Geneとか別のハードの上で高性能化するとか、あと、途中の通信回避のところをもうちょっと工夫しましょうとかというのがされているそうです。こんな感じです。

　以上がプロジェクトの紹介です。これで20分ぐらいですね。残りはいつもの発表みたいな感じになっちゃうんですが。僕がやっているのは、櫻井チームに所属しているので密の方はあまり手を出してなくて、疎行列の周回積分型固有値解法のアルゴリズムの中身、数学的なところを中心にやってます。今日はそこの部分はちょっと触れますが、並列ソフトウェアとして実装する時にどんなことをやっているのかというのを若干触れます。それぐらいです。

　もともと解法の話なので、対象とする問題をもう一回整理しておきますと、対象としているのは一般化固有値問題です。アルゴリズム自体は非線形でもいけるんですが、今回は一般化固有値問題向けの話です。 $A$ と $B$ が複素の行列で問題なくて、複素で固有値。複素平面上の指定された領域Ω内部の固有値および固有ベクトルを全部求めたいというような問題設定になってます。こういうような問題を対象としているので、ユーザーが与えた指定された領域の固有値を求めるというのに対応できるという解法になっています。

　対象とする解法は周回積分型固有値解法という、2003年に筑波大の櫻井先生と、今、南山大にいらっしゃる杉浦先生が共同で作られた解法になってます。周回積分を使ってやることで対象の固有ベクトル成分を抜き出してあげましょうというような解法になっています。
　これは2003年の解法で、当然このまま使っているわけではなくて、その後いろんな改良法というのが提案されています。いろいろ挙げますと、これは非線形向けですが、いろんな解法というのが提案されてます。ソフトウェアで使われているのは、この中でBlock SS-Hankel法というやつとBlock SS-RR法――原著論文のBlock版拡張とRayleigh-Ritz拡張されたものです――の2種類というのが、現状のソフトウェアでは組み込まれてます。僕はこの辺のは最近作ったので、この辺も組み込んでほしいなとは思っているんですけど。実装している人が忙しそうなのであまり言えずにいたわけだけど、実はこっちの方が速いと思うんですが、まあいいです。

　今日の講演は、ソフトで使われている2つの解法の概略というのを紹介します。数学的なところよりもアルゴリズムがどうなっているかというのと、あと、それを使って組まれているソフトウェア、最終的にはz-Paresというのが公開されているんですが、それ以外にもMATLAB用のソフトウェアとかPETSc用のソフトウェアとかというのが公開されてます。なので、この辺もちょっとだけ。あと、これはこのソフトウェアと直接関係するわけじゃないんですけど、我々のプロジェクトで研究している内容の実用化のための関連技術として、さっき「精度保証じゃないけど」と言ってた耐故障性というのをやっているので、こんなのも。他にもいろいろやってるんですが、僕が出せそうなのがこれしかないのでこのデータを持ってきました。ここは僕やっているところなので出せるかなと思って持ってきました。

　流れは、まず解法の概略とソフトウェアを順番に説明していきたいと思います。

　基本的なアイデアは、これは何度かいろんな場で話してるかもしれませんが、元の一般化固有値問題の固有値を極に持つような有理関数、こんなものを考えます。こっち（ $u$ ）とこっち（ $v$ ）のベクトルは今別々に書いてますが、何でもいい、別々で書いてますが。最終的に使う場合はだいたい一緒に揃えて使います。この有理関数の極が元の問題の固有値に対応する。欲しいのは領域 $\Omega$ 内部の固有値なので、領域Ω内部の極だけ求めればいいということで、Cauchyの積分方式を使って周回積分をやってあげれば、ここの内部の極だけ出てきますねというのが原著のアイデアです。周回積分をやって、残りの成分を全部落としてあげましょうということをやります。

　アルゴリズム的に何をやるか。有理関数、解析関数の領域 $\Gamma$ 内部の極の計算というのはいろんなやり方があるんですが、安定性を考えると周回積分を使って出てきたこういう複素モーメント、関数の値にzのk乗を掛けて積分したやつを $\mu_k$ とおいて、その $k$ を0から順番にどんどん大きくしていって作るようなHankel行列を2つ用意します。これの一般化固有値問題を解いてやると、この極が出てきますよというような研究があります。他にもやり方はあるんだけど、このやり方が安定に求まるでしょうというふうにやっている。これはPeter Kravanjaの――櫻井先生もここに入ってますね――この論文を基に、周回積分型の固有値解法を作りましたという論文になってます。なので解法としては、元の一般化固有値問題をこのHankel行列の一般化固有値問題に帰着させてあげるというのがSS-Hankel法の大本です。

　「いろんな解法が提案されていますよ」と言いましたが、もう一回同じことになりますが、これは解析関数の極を求めるという立場で導出されたアルゴリズムです。これをもっと線形代数っぽく作り直して、Rayleigh-Ritzの技法を使ってもうちょっと安定的に高精度に求めるようにしたという解法がSS-RR法と現在呼んでいる解法になります。あと、Block化というふうなことで、あとで詳細は言いますが、重複固有値があった場合に安定的に求まるようにした方法がそれぞれに対して提案されていて、これがBlock SS-Hankel法とBlock SS-RR法――ここがBlockが抜けてますが――というのがあります。
　最近、我々がArnoldi法に基づいてRayleigh-RitzをArnoldiに拡張したような感じでアルゴリズムを組み直して提案した方法というのがこれです。あと、これの改良というのもちょっとやったりしてます。
　関連研究としては、似たようなアルゴリズムでFEASTというのがあります。これはややこしいんですが、アルゴリズムの名前であって、さらにソフトウェアの名前としてもFEASTという名前で公開されてます。z-Paresのライバルのソフトウェアというのが、このFEAST。またこれ（Beyn法）は非線形固有値問題向けに作られていて、基本的には非常に似たような解法になってます。z-Paresはこの（赤枠で囲んだ）2つの解法を中に持っているので、今回はこの2つの解法について紹介したいと思ってます。

　原著論文、SS-Hankel法について、もう一回アルゴリズムの話として紹介します。何をやっているかというと、こんなベクトル（ $v$ ）を考えます。インプットとしてn次元のベクトルを1本持ってくる。これは何でもいいというわけではなくて一応ルールはあって、欲しい固有ベクトルの成分を全部持ってます。どれかに直交しているとそいつは出てこないので、一応全部持っている。だいたい乱数ベクトルとか入れておけば大丈夫なので、乱数ベクトルをよく与えます。
　そういうベクトルを入れて、さっきの関数の部分を計算してやる。 $(zB-A)$ のインバースの $B$ 倍してこれに掛ける。 $z$ の $k$ 乗を掛けて積分したのをベクトル $s_k$ とします。それに、これは一緒じゃなくてもいいんですが、面倒くさいからだいたい一緒にしているんですが、ベクトルの転置をとって掛けて、スカラーで $\mu_k$ をとる。さっきの定義だとこの $v$ 以下はここ（積分の中）にあったんですが、ベクトルなので $z$ に依存しないので外に出している書き方をしていますが、こういうかたちで書けます。さっきと同じようにHankel行列を2つ作ってあげましょうと。一般化固有値問題を解けば、SS-Hankel法はできあがりますねというかたちになります。

　Block化といったのは、ここから何をやるかというと、入力のベクトルを複数本にします。ただそれだけです。1本だったやつを $L$ 本。パラメータが1個増えますが、 $L$ 本入れてあげる。こちら（左）から掛けるのも同じ $L$ 本やってあげると、これ（ $\mu$ ）のサイズが $L \times L$ の行列になるんです。こっち（左）が $L$ 本入ってきて、だからこいつ（ $V$ ）が $n \times L$ の行列で、 $n \times L$ の転置をかけるので、 $L \times L$ の行列。ちょっとややこしいですが、小っちゃい $L$ です。
　これ（ $\mu$ ）を並べた行列というのはHankel行列ではないんだけど、ブロックの単位で見るとHankel行列になるので、Block Hankel行列というふうに呼びます。表記を全く一緒にしちゃっているので申し訳ないんですが、何かここに記号を付ければよかったんですが、早くて間に合わなかったんですが。サイズが $M$ といってたのは、ここの $k$ をどんだけ増やしていくかですね。ここのサイズが $M$ です。 $M \times M$ のHankel行列。ここの $k$ の上げ方を同じまでいくと $L$ 倍されるので、 $LM \times LM$ のBlock Hankel行列が2個できて、これの固有値を解けばいいことになります。数学的にはほとんど変わらないことをやっているんですが、こっちでやると重複固有値もちゃんと拾ってこれるというような解法になってます。

　アルゴリズムの数学的な意味ではこれで完成で、あとはこれをどう実装するかということですが、まずやるのは当然この連続の積分ができないので数値積分をやりましょう。例えば $N$ 点台形則とかで近似しますとやると、こんな感じです。 $N$ 個の値で積分してあげましょう。その近似の意味で^（ハット）を付けてますが、こういうかたちで近似した行列の一般化固有値問題を解けばいいと。

　これでアルゴリズムは完成なんですが、実用上の工夫としてもうちょっと手を抜きます。何をやるかというと、Hankel行列のLow-Rank近似をします。低ランク近似をします。数値的に非常に不安定になってくるので数値的にはランク落ちをするので、いったんLow-Rank近似をしてやって、特異値の大きい方か、このしきい値をどう決めるんだというのはまた難しい問題なんですが、ゼロに十分小さい数値ランクのところまでで止めて、残りは捨ててあげるというように近似してあげて、せっかく特異値分解をしたんだから $\widehat{H}^m$ の方をひっくり返してこっち（左辺）に掛けてあげて、一般化固有値問題が最終的には標準固有値問題のかたちに落ちる。サイズがちょっと小っちゃくなった標準固有値問題になります。
　ですので、実際にやられているアルゴリズム、実装の上でやられているBlock SS-Hankel法というのは、対称の一般化固有値問題がBlock Hankel行列からやる、そのものではないんだけどLow-Rank近似でやって作られた標準固有値問題に帰着させる。これが実装の意味でのBlock SS-Hankel法になってます。この $d$ というのが数値ランクです。だから、固有値数よりも必ず大きくて、ただ $LM$ よりは小っちゃいような値でサイズを落としてるというかたちになります。

　アルゴリズムはこんな感じで、線形方程式を各積分点ごとに掛けて積分してあげる。それを $\mu$ を作るために $V$ 転置を全部掛けてあげる。一方、 $S$ の方でもそれぞれ並べたようなこんな行列を作ってあげるということをやります。Hankel行列を作ってそれのLow-Rank近似をしてあげて、この標準固有値問題の固有値を作る。最後に固有値の方はそのままなんですが、固有ベクトルの方はこの $\widehat{S}$ に $U_1$ が掛かったものが固有ベクトルになると。こんな感じで固有ベクトルを作り出すことができます。これがBlock SS-Hankel法のアルゴリズムです。

　もう一個の方がRayleigh-Ritzの方なんですが。やっていることは非常に似ているんですが大本のアイデアがちょっと違って、何をやっているかというと、さっき計算していた $S$ ―― $S$ というのはこいつ（スライド29の $\widehat{S}$ ）の連続でやったやつ。だから、ここ（ $\widehat{S}_k$ の数値積分）が連続の積分になってて並べたやつが $S$ です――がちゃんとインプットの行列ができていて固有値数よりも作ったサイズが大きければ、作った部分空間がちゃんと欲しかった固有ベクトルの部分空間に等しくなりますよ、連続の積分でやれば等しくなりますよということを言ってます。こんな定理があるので、これがあるんだからHankel行列とか作るのが面倒くさいので、そんなことをやらずにRayleigh-Ritzを使えばちゃんと固有ベクトルをこうやって出せますよねというのがアイデアです。

　アルゴリズムはすごい単純で、この $\widehat{S}$ の方、数値的にやるのでどうしても近似が入っちゃうんですが、 $S$ の^（ハット）が入ったやつの部分空間、レンジの空間に固有ベクトルが入るようにします。残差が同じ空間と直交するようにしましょう。で、Rayleigh-Ritzの完成です。
　これでいいんですが、こっちでもさっきと同じような工夫をしていて、基底が数値的にランク落ちするのでLow-Rank近似をしてあげよう。今度は基底の行列のLow-Rank近似をしてあげて、この $U_1$ でRayleigh-Ritzをやりましょうと。数値ランクの部分だけを取ってきた $U_1$ という行列です。この部分空間でRayleigh-Ritzをやりましょうというのが、実際にやられている方法です。Rayleigh-Ritzなので、何も考えずにこういうような元の $A$ と $B$ を使って一般化固有値問題がこういうふうにできますよというかたちになります。
　ですので、Hankelと違って、こっちは最終的には一般化固有値問題に落ちます。対象の一般化固有値問題はRayleigh-Ritzの技法によって得られる小規模な一般化固有値問題に落ちるのがSS-RR法の方法になります。

　アルゴリズムもさっきとほとんど同じでやっていることは一般的です。線形方程式を解いて積分をしてあげる。こっちは $\mu$ とかを作らずに $S$ をそのまま作って、これのLow-Rank近似をやってあげる。小規模な一般化固有値問題を解いて、固有対はRayleigh-Ritzのでこのように作ってあげればいいです。基底ので作った行列 $U$ に小規模問題の固有ベクトルを掛けてあげれば、元の固有ベクトルが出てくるのがアルゴリズムになっています。

　これでアルゴリズムの説明はひととおり終わりで、ここからちょっとソフトウェアの話に入ります。今までの紹介した両方の解法というのは、アルゴリズムの流れがそれぞれ若干違うんですがほとんどやっていることは一緒で、こういう3つのステップからアルゴリズムというのはできてます。
　最初は各積分点で連立一次方程式を解いてあげる。どんなのかというと、こいつです。 $(zB-A)$ のインバースかける $BV$ 、こいつです。 $L$ 本の右辺ベクトルがあるような係数行列が積分点ごとに行列 $B$ でシフトしてくるような感じなので、シフト線形方程式ではないんですが、こういうようなかたちで積分点ごとに係数行列が異なる、各積分点に $L$ 本の右辺ベクトルがあるような連立一次方程式を解く。で、積分計算、こいつを計算してあげて、何らかの行列の特異値分解をやって、標準・一般どっちだか分からないけど小規模な固有値問題を作って解いてあげるというような3つのステップになってます。
　当然といえば当然なんですが、計算コストというのはほとんどこの部分（ステップ１）の線形方程式の求解になってます。ここ（ステップ２）はもう足し算するだけだし、ここ（ステップ３）も小規模に落ちてるので、特異値数が小さいような領域を設定してやればここ（ステップ３）のサイズはすごく小さいので、ほとんどの計算時間はこの部分（ステップ１）になってきます。ですので、並列化とか計算コスト、並列化したときに最終的にどれぐらいの計算時間がかかるかというのは、ここの部分（ステップ１）がどれぐらいちゃんときれいに並列化できるかということになってきます。
　最初の方でも少し説明しましたが、この解法というのはここの部分で高い階層性があると。何かというと、そもそも積分領域が独立にできますよね。この辺りの複素平面上のある特定領域の固有値が欲しいとなったときに、その領域を分割して設定することもできます。そうすると、複数の領域に分割できるとその計算というのは独立にやることができる。
　二つ目が、数値積分というのが積分点ごとに完全独立な計算になっているので、これは独立に計算できますよねということです。三つ目は、最後の線形方程式もいろんな並列化の技術というのが既にありますので、こいつの階層はちょっと弱いんだけれども線形方程式も並列に解くことはできると。ということで、3つの階層の並列性を持っているので、最近の超並列計算機、京とかそういう非常に大規模の並列計算機では性能が出るでしょうということになってます。

　さっきも出しましたが同じ絵ですが、さっき言った積分領域と積分点、あと線形方程式自身の並列性というのをアーキテクチャの階層性に合わせてやることで、通信というのをローカルにできる。全体の通信というのを減らすことができて、アルゴリズム全体の性能が上がりますということになってます。
　これはソフトです。アルゴリズム開発者からの目線としてはこういうふうになってます。ソフトウェアが実際にそう組み込むかはまた別の話で、これに合わせてソフトというのは開発してますというのが現状です。

　もう一つあります。二つ目に、これが関連ソフトウェアとの比較で、どういうことをやっているかということです。一般化固有値問題を射影法で解くオーソドックスな方法というのは、Shift invert Arnoldiみたいな方法です。何をやるかというと、反復ごとに線形方程式を解いて空間を増やしていく。とやると、線形方程式の求解が逐次に現れる。だから、並列性が損なわれる。そのシフト点を複数置けるのでその意味で階層性はあるんですが、線形方程式の求解自体は逐次に出てきちゃうのであまり性能が出ませんというのが、従来から言われたArnoldi-basedの方法です。
　もう一個はFEASTです。FEASTは非常に似たような解法で、何が違うかというと……。（スライド３２に）戻りますと、FEASTはこいつ（block SS-RR法）とほとんど一緒なんですが、ここ（ $S$ ）が $S_0$ しか使わないやつです。それをすると何が起こるかというと、部分空間のサイズを作りたい分だけ線形方程式の右辺ベクトルをどんどん増やしていかなきゃいけない。このBlock SS-RR法はここでスカラー倍だけ動かしてベクトルの本数を増やせるので、右辺ベクトルの本数を、ここでいうパラメータを $M$ というのを入れた $M$ 分の1に抑えることができるというのが違います。 $M$ を1にするとFEASTとほとんど同じことができるというような解法になってます。そういう意味では、拡張した解法というふうに見ることができるかと思います。
　そういう意味で、線形方程式自体は独立に周回積分に基づく解法なので解けるんだけれども、右辺ベクトルが非常に多くなってしまうので、ここの部分（線形方程式の求解部分）が演算量が多くなってしまう。対して、我々が作っている解法は高次モーメント、さっきの $k$ のべき乗のところをうまく使うことで右辺ベクトルの本数を減らして、代わりに高次モーメントのスカラー倍の計算だけでベクトル本数を水増ししてやることができるわけで、線形方程式の計算量を削減することができるというのが売りです。
　そういう意味で、FEASTが今いろんな所で海外で使われているんだけど、こっちに切り替えてくださいねというような説明のスライドになってます。

　アルゴリズムは以上なんですが、それをどうソフトウェア化しているかということで、現在3つのソフトウェアというのが公開されてます。
　一つがMATLAB版で使う用の「試験的に使ってください」というための方法です。周回積分型固有値解法が射影法に基づく解法なので非常に性能が問題依存というのが大きいので、ユーザーが自分が持っている行列で使えるかというのはやっぱり試してみないと分からないんです。全固有値を求める場合は比較的そういうのが少ないので、問題サイズだけ分かっていればだいたいこれぐらいの性能というのは分かるんだけど、射影法の方はそうはいかないのでやっぱり使ってもらわなきゃいけない。そういう意味で、MATLABで使って簡単に使えるように公開されたソフトというのがSSEIGというのがあります。これは並列化はされてないです。シリアルで動く。MATLABが勝手に並列にしてくれるんだと思いますが、複数コアでたぶん動くんだと思いますが、MPIみたいなそういう通信はされてないです。
　二つ目がCISSというソフトウェアで、これはSLEPcという、なんていったらいいのか僕もよく分かってないんですが、分散並列環境でいろんなソフトというか比較的簡単に使えるようなパッケージがあって、その中に我々の解法が組み込まれています。そうすると、Arnoldi-basedの方法とかと比較が簡単にできる。これは分散環境でちゃんと動くので、そういう意味では非常に大規模は苦しいんですが、中規模ぐらいだったらこれで比較的簡単に動かすことができます。
　最後のが一番高性能なやつということで、FortranとMPIでしっかり組まれたやつというのが公開されてます。こいつが最終的なアプリのユーザーが使っているやつです。この（SSEIGとCISSの）2つは中規模問題とか、その使うための準備として使えるかなというのの検証用というかたちになります。公開ポリシーは完全にフリーで使ってもらってもいいし、改変しても構わないというかたちになってます。Webで現在公開されているというようになってます。

　一個一個のソフトウェアを簡単に紹介したいと思います。
　SSEIGはMATLABで使われているやつで、もう言っちゃいましたが逐次版です。MATLABの自動並列は有効になるはずです。MATLABで組まれてて検証用に使ってくださいと。対象の固有値問題はちゃんと全部に対応していて、実でも複素でも、標準でも一般でも、当然、密でも疎でも動きます。積分は本当はアルゴリズム的にはどんな積分則でもいいんですが、このMATLABは楕円領域に限定していて、あと $N$ 点台形則で組まれてます。領域も別に本当は何でもいいはずなんだけど、このソフトは楕円でしか動かないで、線形方程式の求解も本当は反復法を使って解くのがいいんだけど、MATLABはバックスラッシュで解いているというのが現状です。

　このスライドはたぶん本当はデモ用に作っているやつなのでもう結果が配置になっていますが、非常に簡単にサンプルコードとかも付いてて簡単に動かすことができます。MHD1280って「Matrix Market」にたぶん載ってる問題だと思います。1280次元の小っちゃい電磁学の問題で、全固有値はこんな感じで、複素平面上にバラバラとあるような問題です。この（赤枠の）部分のに領域を円として設定します。
　右辺ベクトルの本数、入力のベクトルの本数が16本とか、 $m$ 、高次のモーメントをどれぐらい作るかというパラメータと、積分点数というのを設定してやって、こういう固有値問題を……これ、結果が無い。ごめんなさい。そうか、デモ用だから結果が無いのか。ちゃんと求まりますよということになってます。

　二つ目がCISSというソフトです。これはSLEPcとかPETScが入っている環境でしか動きませんが、3階層あった分散並列全てにちゃんと対応して並列化されています。実・複素、標準・一般、疎・密、全部対応していて、こちらも楕円と $N$ 点台形則に限定されてます。線形方程式の求解は、PETScがいろんな線形方程式の求解のパッケージを含んでいるので、そこのライブラリーから選択して使ってくださいというかたちになってます。なので、疎行列向けのKrylov部分空間法とかを使うことができます。分散なので行列の格納形式というのがいろいろ考えられるんですが、これはPETScの行列格納形式に準拠しているというかたちになってます。

　これはWebで公開でされていて、PETScとかSLEPcがこんなライブラリーで中にアルゴリズムを含んでいますよという一覧なんですが、固有値問題だとKrylov-Schur法とJDとかCG-basedとかがあるんですが、その一つにCISSというのが組み込まれているというかたちになってます。なので、この辺とかも簡単に使うことができるので、比較とかするためには簡単に使えるような環境にはなってます。

　最後が、我々がイチ押しというか最終的に成果として作っているz-Paresというアルゴリズムです。当然、全ての階層の分散並列に対応してます。環境はFortranとMPIで組まれていて、今、Cのラッパーというのも作ろうとしているらしいんですが、需要が無いということで後回しになっているそうです。アプリの人がみんなFortranらしくて、実際にCのはあまり需要が……需要があれば作りますが、今のところはまだ無いみたいです。
　もちろん全部に対応していて、楕円と $N$ 点台形則でデフォルトは動くんですが、その他ユーザーが指定した領域とか積分則にも対応できます。というのは、積分点と重みをユーザーが全部指定するということができるので、自分で計算して放り込めばちゃんと動くようにはなってます。
　線形方程式の求解は何でもいいですというんですが、実は、後で言いますが、これはリバース・コミュニケーション・インタフェースというのを使ってて、「ここは自分で組んでくださいね」というかたちになってます。ただ、デフォルトとして疎行列向け直接法とかは使えるようにはなっていますが、自分で勝手に組み替えることが簡単にできるようになってます。
　行列の格納形式も完全にユーザーが定義することができます。z-Paresの中ではベクトルの分散は知らないまま動くようになってます。だから、ベクトルのどの部分をどのノードが持っているかは分からなくても動くような感じになってます。通信はAllReduceしかしないので、そこは正直いらないというかたちです。なので、その辺も全部ユーザーが定義することができます。

　アプリケーションごとにデータ分散形式というのは非常に異なるので全てに対応するのは不可能なので、そこはもうユーザーに丸投げして、こちらでは触らないようにしましょうというのがポイントです。そのためにリバース・コミュニケーション・インタフェースというのを使っています。これはいろんなソフトで既によく使われているやつらしいです。

　どんなのかというと、線形方程式の求解とか行列ベクトル積とかで行列のデータが欲しくなるんですが、そういうときは制御をいったんユーザーに返すということを毎回やります。ユーザーが書くサブルーチンみたいに使って、そこに毎回飛んでいきます。そこでユーザーが書いたコードで計算がされて、また中に戻っていくというかたちになります。だからユーザーはこの中に「ここからベクトルを1本、このベクトルに行列を掛けたベクトルをここに入れてください」という指示がされていて、自分の行列に合わせてそこは組んでくださいというかたちになってます。線形方程式も同じように自分でやってくださいと。ただし、簡単なやつだけはこっちで提供しているけど、性能を出したかったら自分でやってねというかたちになってます。while文で何度もこのz-Paresのメインルーチンを呼び出しながら、チョイチョイ返ってくるというかたちを繰り返すというかたちになってます。
　ARPACKとか他のFEASTソフトウェアとか、これも似たようなソフトなんですが、こういうのにも採用されている比較的よく知られた実装方法です。若干、実装は非常に複雑になるんですが、よく知られた方法ではあります。

　実装はこんな感じでややこしいんですが、これはwhile文で、タスクがfinishするまでずっと回ってくる。で、何度も呼び出される。こいつがメインサブルーチンで、何度も呼び出される。タスクが例えば「まとめてやってね」というのになったら、ここがグルグル回ってて、IF文で判定して、ここがユーザーでまとめてこう書くと、ここにユーザーが線形ソルバーを実装するとかいうかたちになってます。work1に初期ベクトルが来て、work2に計算されたベクトルを入れて返すとか、そういう指示がされます。「ここだけ自分で書いてね」というかたちになってます。

　簡単なテスト結果だけお見せします。実際のアプリでの結果というのはちょっと僕が出していいのか怪しかったので今日は持って来なくて、簡単な実験で。ELSES Matrix Libraryというのは一緒のグループの星先生が作られているライブラリーで、電子状態計算の行列をストックしているライブラリーです、一般化固有値問題でいっぱいの。
　一つ目が、実対称一般化固有値問題で9,000次元。パラメータはなんかこんな感じで入れて、線形方程式は反復法ではなくてMUMPSという疎行列向けの直接法のライブラリです。これは京で1プロセスから16プロセスまで、積分点が16点、これは32なんですが、対称の固有値問題の場合は上下の積分点を上だけで計算できるんです。下の方は複素対称で対称点を使って計算をしなくて済むので実質16点しかやってないんですが、ここの中位階層のみの並列化しかやってない、すごい単純なケースです。

　これが実験結果で、縦軸が計算時間、横軸がノード数で、ノード数を増やすと2倍に計算時間のほうが減る。これはなんでかというと、ほぼ当たり前なんですが、線形方程式は完全独立に解けてて、線形方程式は16個あるので、一人で解いたら16回解くのが、二人だと8回ずつ。だからここがリニアに完全に落ちる。そういうことで非常にスケーリングがいいですよという結果です。

　これは相対残差です。横軸は固有値の値で、この範囲の固有対を求める。縦軸が残差で、この辺がマイナス13乗、14乗で、この辺は10乗で、ところどころ精度が若干落ちているんですが、比較的いい精度で全ての固有対を探し出すことができているかたちになってます。

　これもほとんど同じで、もうちょっと大きい問題です。複素の問題で若干大きな問題で、同じような結果です。

　実験をやってみると、同じようにちゃんとリニアにスケールしますねというかたちになりました。

　精度も同じようによく出てる。全部がよく出てるわけじゃないんだけど、この辺のやつはいいけど、ところどころ悪いやつもいるんですが、ちゃんと精度は出てるという結果になってます。
　ここまでがソフトウェアの紹介です。もっと大きい問題で実アプリでちゃんと対も出てますというのもあるんですが、ちょっと今日はお見せできないので。

　最後が少しソフトの話からずれるんですが、最近、高性能計算とか大規模環境でやる時によく言われるテーマの一つで耐故障性というのがあります。何かというと、京でもノード数が8万個ぐらいあるんです。それぐらいあると、どこかがしょっちゅうコケるというのがあって。京はその辺が壊れにくい性能らしいんですが、「24時間連続稼働できます」とかいうのを売りにしているんです。逆に言うと24時間ぐらいしか動かないんです。1カ月とかずっと動かしっぱなしにすると、必ずどこかが故障する。それぐらいの頻度ではだいたい壊れるということになってます。
　ですので、最近言われているのは、どこかでノードが故障しても計算が継続できるようなアルゴリズムだったりハードウェアだったりを作らなきゃいけないねということが言われてます。その一環で、この周回積分型固有値解法の数学的なソフトの意味でアルゴリズムベースの耐故障性というのを研究してます。研究してますというか、もともとそういう性質があったので、耐保障性をちゃんとやろうという話になった。
　ベースになっているのは、もともとの解法の誤差解析です。誤差解析はこんな結果になってます。ちょっと式の定義を全然してないみたいですが、こっち（左辺）は残差みたいなものだと思ってください。指定された領域、欲しい固有対の残差みたいなやつ。こっち（右辺）がこの残差を上から抑えるパラメータ。こいつら（ $\alpha, \beta, \gamma$ ）はパラメータでちょっといったん忘れて、こういう $f$ の関数の自分の固有値の関数値と $LM+1$ ―― $LM$ というのは部分空間のサイズです――番目の関数の比で抑えられるというようなことをやっています。これが誤差解析の結果です。

　この $f$ という関数は何かというと、積分がどれぐらい厳密にできていますかというような指標です。連続積分だと $f$ の値は領域内部で1、外部で0というような関数の値になります。連続積分じゃなくて数値積分なので、それがなまってくるようなイメージで思ってください。この絵だと-1から1の所が領域内部、だいたい1ぐらいの部分。遠くの方だとだいたい0。間の辺はなんとなくなだらかに落ちていくような感じになってます。
　この（スライド５１の定理の）式がいっているのは、これ（右辺の分母）は自分自身の値です。これは領域内部にいるので、だいたい1の値を取ってます。こいつ（右辺の分子）は部分空間のサイズ足す1です。こいつが0に極めて近ければ、例えばこれ（分母）が1でこれ（分子）が10の-16乗だと、だいたいこいつの残差が10の16乗ぐらいまで落ちる。で、いい性能が出ますねということを言っている。ちょっと定数は無視しますが。
　この $\lambda$ は、この関数が大きい順にソートされているというかたちにしてください。なので、 $\lambda$ の $LM+1$ 番目が0になるぐらい部分空間を大きくとれば、こいつは0になるんだから精度がいいですねということ。だから、近接固有値があるとよく周回積分のやつは精度が悪いんじゃないかと言われるんですが、隣に固有値があっても $LM+1$ がこの辺にあるぐらいまで部分空間サイズをちゃんと大きくとれば精度は出ますよということを言っている。固有値が実際どこにあるかはあんまり関係なくて、部分空間をそれに合わせて大きくとりなさいということ。

　誤差が入った場合、故障がどこかで入った場合というのを考える。どうやって入るかというと、具体的には計算時間の一番重いところで誤差が入る可能性が高いので、線形方程式を解いてる時に特定のノード……特定の線形方程式の解が完全に壊れましたというケースを考えます。なので、解ベクトルが乱数ベクトルに置き換わりましたという状況です。
　そうすると何が変わるかというと、パラメータも若干変わるんですが一番大きく変わるのはこの部分で、 $LM+1$ だったのが $LM-L+1$ になります。どういうことかというと、部分空間のサイズが故障によって縮められるだけ。だから、別に大きくは壊れないんです。部分空間をちゃんと大きくとってれば、故障も想定してちゃんと大きくとっておけば、縮められたやつでも十分小っちゃくなるように部分空間を設定します。そうすると、解の精度にはほとんど影響しませんということ。あらかじめ情報を余分に持っておけば、若干失われてもちゃんと解は復元できますよということです。

　これが実験結果で、これは故障が無いやつです。これ、見にくくて申し訳ないんですが、部分空間を大きくとって、固有値が0.01から0.1刻みで大きくなっていくような感じです。領域が-1から1で、1.01という非常に領域外部に近接した固有値があるんですが、部分空間サイズをこれが30とか40ぐらい大きくとってあげると、ちゃんと-15乗ぐらいの精度で解が出てくる。外側は別に出てないですけど、これは別にいらないやつなのでどうでもよくて、欲しい固有対に対してはちゃんと精度が出てますねということを言っています。 $LM$ が例えば30だと、31番目のフィルター関数の値の大きな固有値を入れればいいので、だいたいこの辺りを入れます。31というのは3.01なので、だいたいこの辺。近いような値の残差の精度が出てますねということを言ってます。

　故障した場合はどうなるかというと、こんな感じです。赤が故障が無かったやつの30。緑は、部分空間サイズを増やさなかったらやっぱ精度が悪くなるんだけど、ちゃんと大きくとっておけば故障してもほとんど――若干、精度は落ちるんですが――精度が劣化せずに、故障が起こった時にも正しい固有対を計算することができるということです。
　これでこのアルゴリズムの耐故障性があるということをやっています。だから、なんかやったというわけではなくて、もともと解法がこういう性質を持っていたから良かったねということで、今ちょっとアピールしてるという段階です。

　今日の話はこれで終わりなんですが、今日の講演でこのまとめのスライドが若干正しくないような気もしますが、解法の数理的なところというよりは我々が所属しているプロジェクトの紹介になりました。ソフトとして、今3つのやつを公開してます。z-Paresは今ここで公開されていて今も改良中なんですが、一応たぶんソースでちゃんと公開されているので使ってみたい方は使ってみてください。以上です。

質疑

Ａ：どうもありがとうございました。（拍手）質問があればお願いします。
Ｂ：無限次元の周回積分を無限次元のままやろうとした人とかはいますか？　というのは、だいたい分点を取って周期積分に落としてそれから先につかまるんですけど、ここの部分の精度を落とさなければここがうまくいくとか……。
今倉：できるかは置いておいて、落とさなければ、このフィルターがさっき言った厳密にシャープなやつなんです。そうすると、部分空間サイズは中の固有値数と一致した段階で全固有値が必ず出るわけです。他の丸め誤差を無視すれば厳密に出るので、無駄な計算が無くなってうれしいです。
Ｂ： Hankel行列とかを使って適当に分点を取ってきて、それがあまりにも大き過ぎちゃうと、出てきた行列がランク落ちしちゃうみたいな感じになるのかな。
今倉：そうです。本当はフィルターが厳密にあれば必ずランク落ちするんです。そういうものなんです。だから最後の実用上の工夫でLow-Rank近似しているのは、そこを落とすためです。
　だから逆に言うと、数値的なフィルターでもランク落ちするぐらいちゃんと部分空間を広げなきゃいけない。ランク落ちしないということはまだ足りなくて、十分広げてランク落ちするぐらいになったらもういらないので、空間サイズとしてはピッタリですねという感じになります。
　ソフトのパラメータの自動設定があるというふうに言ったかと思いますが、それはそこの情報を使って特異値分解をした時にランク落ちしているか、最初、特異値が10のマイナス十何乗まできてるかを見て、きてなければ部分空間サイズが足りないのでもうちょっと増やしましょうということを内部でやっているんです。
Ｂ：ありがとうございます。
今倉：積分も厳密にやればいいというのはそのとおりで、前、相島君が教えてくれてたんですが、それを全固有値でやったやつがdqds 法に対応するそうです。僕はちょっとあまりフォローできてないんですが、そういう歴史的背景があるんです。
Ａ：これ、積分の分点数は結構……さっきの例はすごく楕円自体が小っちゃかったからかもしれないけど、すごい分点数が少なかったんですけど、例えばこれは単位円と書いてあるけれども、これは大きいけど分点ってどれぐらい？
今倉：分点は、これ、32点。
Ａ：あんまり大きくないね。
今倉：結局、分点を増やすとフィルターがシャープになるんです。シャープになってうれしいのは、ここのやつが早く落ちるんで、部分空間サイズを小っちゃくしても、lm＋1というのが固有値数よりちょっと増えたぐらいでも、そこの固有値の……うまく説明できないけど、フィルターがシャープだと、近接固有値の値でもちゃんとこいつが0に近くなってくれる。
Ａ：そうだけど……
今倉：それはうれしいんですが……
Ａ：近似解のほうが悪くなったりしないのかな。
今倉：近似解の精度はこいつで効くので、丸め誤差の話を無視してこういう積分の近似だけで言うと、この式でちゃんと表現できているので。部分空間を広げることと積分をシャープにやることが、だいたい等価になるんです。等価というか、意味の上では近い。あとはどっちが計算量的にうれしいかという感覚です。
　積分点を増やすのは線形方程式の本数というか行列の数を増やさなきゃいけなくて大変で、部分空間を増やすのは右辺を増やすだけでいいわけです。どっちが楽かというと、右辺ベクトルを増やす方が楽、計算量的には楽なのでだいたいそっちが使われるんですが、その加減がどっちがいいかというのは実は非常に難しくて、だいたいこれぐらいの32とかぐらいが使われているけど、場合によってはもうちょっと増やさなきゃいけないことももちろんあります。
Ａ：なるほど。あと、ごめんなさい。僕、途中で抜けたから、もしかしたら分かってないかもしれない。故障無しの場合と故障ありの場合と分けてるんだけど、故障したということは判定できるの？　つまりだから……
今倉：場合によってできるような方法もあるし、できないこともある。今どういう故障があり得るかというのは実は僕もそこまでちゃんとフォローできてないんだけど、今こういう研究も発展段階で難しいらしいんですが、想定しているのはビットフリップで中性子線点かなんか当たってビット上の0、1がひっくり返ったとか。だから「数字がちょっと変わっちゃったね」というケースで。
　どうやら非常に大きく値が変わった場合は判定ができるらしい。例えばInfとか桁が大きく上がった場合は、ハード的に判定が可能だということはあるらしいです。Infになったとか、Nanが。そういうのは判定できるらしい。ただ、丸め誤差ぐらいの感覚で微妙に変わっちゃったら判定は非常に難しい感じ。この耐故障は判定すらしなくていいんですよと言ってるんです。部分空間をもともとちゃんと大きくとっておけば、判定なんかしなくたって気にしなくていいので。
　判定できれば他にやり様が無いことはなくて、その壊れた積分点の計算を回避するとか。一番簡単なのはやり直すことですけど、やり直すのは計算でも絶対やりたくなくて。その積分点を足し算しないとかという方法もある。足し算しないと、結局これと同じ定理を使って解を解析できるというかたちになります。
　判定できるかは場合によるのかもしれないし、なかなか……。多数決を採って合わせて、チェックとかソフト的には簡単にできることを。それがいいのかはなかなか……。ちょっと難しいところです。僕もあまりフォローできないので。
Ａ：なんか比較的簡単な式で抑えられて、それでなんか可能性が分かったら、それはそれでちょっと面白いというか。他に何かありますか？
Ｃ：いいですか？　25ページの紙、いいですか？　そのところを僕、完璧に理解してなかったんですけど。それは要するに $\theta_i$ というのはもともとの求めたかった問題の固有値に一致してるんですか？
今倉：一致します。
Ｃ：一致してて、 $t_i$ はなんか行列を掛けてあげると、もともとの固有ベクトルになるんでしたっけ？
今倉：そうですね。 $M$ 本だからここの $s_0$ から $s_{M-1}$ まで並べたやつに $t_i$ を掛けると、固有ベクトルに戻ります。
Ｃ：さっきから言ってた部分空間というのは、これでいうと例えば $M \times M$ というのが部分空間の……
今倉：部分空間でいうと、これは分かりにくいんだけど、よくよく考えると部分空間というのはこの $s$ を並べたやつが部分空間なんです。
Ｃ： $s$ で張られる、全部の $s$ を並べて……
今倉： $s_0$ から $s_{M-1}$ で張られる空間が固有ベクトルを含むようになるので、そこで射影法をやっているという解釈はできます。なのでRayleigh-Ritzと結局似たようなことをやっているんだけど、導出は全然違うんだけど結局やっていることは非常に似てて。
Ｃ：部分空間ってずっとおっしゃられたその意味がちょっとまだよく分かってなくて。
今倉：部分空間射影法というのは、結局、欲しい固有ベクトルを何らかの部分空間に近似しておきますねという方法です。
Ｃ：なるほど。もともとは行列のサイズMaxだとありそうな感じで……
今倉：そうそうそう。 $M$ 次元空間でどこかに固有ベクトルがあるから、 $N$ 次元空間で考えなきゃいけないけどそれは大変なので、これでいうと $M$ 次元の部分空間に落として、その中で考えましょうと。
　実は、これを連続でやれば $s_0$ から $s_{m-1}$ で作った $M$ 個の $M$ 次元の部分空間の中に必ずちゃんと正しい解がある、固有ベクトルがあるというのが数学的に分かっているので、これで十分解ける。あとはこれを近似しちゃうので、実際、当然出ないんですが。
Ｃ：その $M$ というのはどこで決定されるの？
今倉：パラメータですね。そこが難しくて。結局、連続だと内部の固有値の数でいいんです。ただそれは事前には知らないので、大きくしていって、どこからって。これが連続できたとして探し方は、特異値分解をしてみて、こいつがランク落ちしたかしないか。ランク落ちしたということはサイズが十分足りましたということになる。
Ｃ：そうか。そう予想できるということですね。
今倉：そうです。
Ｃ：要するに $\Omega$ に依存して……
今倉：そうです。
Ｃ： $\Omega$ が大きければ大きいほどサイズが大きいことが期待されて……。
今倉：ただ、まあ、分からないので……。
Ｃ：まあ、分からないよね。
今倉：物理的にだいたい分かっている場合もあるらしいんですけど、基本的には分かってない。分からないときにどうするかという研究も我々のプロジェクトでもやってて、固有値の分布を推定するみたいな方法があって、どの辺に固有値がどれぐらいあるかというのを。簡単に言うと、粗っぽく固有値計算をやってみると。周回積分の方法、この方法の近いやり方で領域の円内の固有値の数を計算するという方法があるんです。数学的にはちゃんとできるんだけど、それをすごく粗っぽくやって数を事前に計算する方法もできることはできます。
Ａ：どのぐらい正しく出てくるの？
今倉：結構出るっぽい。けど、そこもだから前処理っぽく使うので、どれぐらい厳密にやっていいのかというのは難しくて。結局、周回積分を使ってやるので、その積分点をどうするか。結局4点ぐらいしか使わなかったり。線形方程式も解くのも大変なので、すごいラフに1桁、2桁しか解かないとかでやるので。どれぐらい合うかというと、まあまあ。逆に数は難しいけど、無いところを無いというのは比較的当たるらしいです。
Ａ：0になるというところだよね。
今倉：そうそう。こっち側は結構いくんだけど、「10個ですか？」「12個ですか？」というのは結構難しくて。
Ｃ：それは無理だと思うんですけど。
今倉：その辺も。ただ、現状のソフトにそこまではまだ組み込まれてないんだけど、一応、技術的にはそういうのもありです。
Ｃ：ランク落ちするまで広げていくイメージだったと思うんだけど、それって重複固有値があっても大丈夫なんですか？
今倉：こいつはダメで、ブロック化した場合に重複固有値も拾ってこれるというところですね。
Ｃ：なんで、それ大丈夫なの？
今倉：結局、ここのなんかフィルターのイメージのエラーということ。ベクトルを1本入れてフィルターをかけてアウトプットを出すんです。そのときに固有値が重複すると同じベクトル2本が自由に動く。それを2本入れておけば、ちゃんと両方再現できますねという、たぶんそういう技術です。
Ｃ：なるほど。（笑）
今倉：うさんくさいな。
Ｃ：分かりました。まあまあ、 $v$ を固有ベクトルに展開させるという式を見るというところで……。
今倉：そうね。ちゃんと標準形にバラしてやって、固有ベクトルに展開するということをすると分かる感じがする。
Ｄ：そのうち2本を一遍に中に入れていくイメージが……、 $L$ 本……
Ｅ：だから $L$ 本までだったら十分で、重複量が $L$ までだったら……
今倉：そうそう。だから本当に何も分からないと、 $L$ でどこまで大きくとらなきゃいけないかというと、固有値数分まで大きくとらないと全体像は見えなくて、そうするとこの $M$ って1にしなきゃいけないじゃないかっていう話もあって。そうすると、これってFEASTになりますよねという話もあって、なかなか難しいところなんだけど。実際ここまで重複するかと言われると……
Ａ：政治的な臭いがしました。
Ｃ：randとかで作ると、ほとんど重複しないですもんね。
今倉：係数行列のほうの意味なので。係数行列の領域 $\Omega$ 内部の固有値が2個重なっていると、 $L$ を2以上にしないと解けないということです。
Ｃ：数値データにおける重複固有値とは何だろうという問題はあるな。（笑）
今倉：そうね。何を重複とするかというのも分からない、難しいかもしれない。物理の連続の世界では重複してるけど、離散になった瞬間に重複なんかするわけないとかかもしれないし。
Ｃ：積分も近似でやっちゃってますしね。
今倉：もう訳分かんないんだけど。ここはだいたい大きくしちゃいます。あともう一個理由があって、 $M$ が無尽蔵に大きくはできないから。これが数学的な意味では積分点数より大きくはできないというのがあります。
Ｃ：そうなんですか。
今倉：そうです。だから、これは実はあんまり大きくはしたくなくて、だいたい8とか16ぐらいしか大きくできない。あと、数値的に $M$ をあまり大きくすると不安定になるというのも分かっているので、実はどんどんいくらでも大きくすりゃいいというもんでもなくて、そんなに大きくない。だいたいは $L$ の方で部分空間サイズは広げていって、 $M$ は8とか16辺りの決め打ちという感じが実用上の話です。
　演算量だけ考えると $M$ を増やした方が必ず少ないんですが、線形方程式の右辺ベクトル本数を極力下げたい。ただ、数字的に若干不安定になる。イメージとしてはKrylov列を量的にexplicitに作っているというイメージなので、 $M$ をあんまり大きくすると性能が落ちてしまう。数学的な意味で解の精度も落ちるのも分かっている。
Ｃ：先生、まだ時間ありますか？
Ａ：全然あります。
Ｃ：30ページに確か定義が書いてあったと思うんですけど。それが……
今倉：これは対角化可能な……
Ｃ：僕が理解力が足りてなくて、言葉の定義が分かってなくて。
今倉：ごめんなさい、そうですね。これ、部分空間のサイズです。こっち、定義してなかったですけど、 $\Omega$ の中の固有値の数。「これより大きくしてください」というのは当然で。もう一個は入力行列が、これ、固有ベクトルの成分全部をちゃんと持ってくださいね。もともと入れたやつが何かの固有ベクトルと直交していると当然取れるはずがなくて、乱数とか適当に入れて、欲しい固有ベクトルの成分を必ず持ってくださいということを言ってる。
Ｃ： $X_\Omega$ って何でしょう？
今倉：書いてないね。こいつだ。領域の $\Omega$ 内の固有値に対応する固有ベクトルで張る部分固有値、ごめんなさい、並べた行列。
Ｃ：なるほど、なるほど、分かりました、分かりました、分かりました。
今倉：ここまでやって、これを連続でやったやつが $S_k$ で、こいつが $S_k$ で、ここが $k$ 上で増えるところのやつがここの $k$ に対応してて、これを――いろいろnotationが足りないな――0から $M-1$ まで並べたやつが $X$ 。最終的にこの列ベクトルで張る部分空間でRayleigh-Ritzをやりましょうという解法です。
Ｃ：すみません、 $R(\cdot)$ って何だっけ？
今倉：レンジ空間。行列のレンジ。
Ｃ：レンジでいいんだ？　分かりました。
今倉：これ、「スパン」って書いてあったのを怒られて、「これはレンジにしなさい」って保國君に言われたので「レンジ」に変えました。（笑）列ベクトルで張る部分空間。だから、欲しい固有ベクトルで張る部分空間とこいつが等しくなるんだから、そこの部分空間で小規模問題でRayleigh-Ritzをやれば、ちゃんと厳密に出ますねということです。これはもちろん連続でやれば出る。これをどうせ近似でやるので、当然出ないけれども。
　ここが厳密にイコールになるのは対角化可能なときです。行列束が対角化可能。一般にはこっち向きかな。こっち向きの $\subset$ になります、たぶん。（笑）必ずなります。
Ｃ：全体でという意味ですよね。行列全体で対角化が可能なという意味？
今倉：そうです。対角化可能というのは、AとB、2つの行列が別々の正則な行列で挟んで両方対角にできた場合が、行列束が対角化可能。一般にはそうやって習いましたね。これだけたぶん、この定義は成り立つ。でも似たような話が別に一般の行列束でもなるので大丈夫なんですが、こっちの方が大きくなる。こっちより大きくなるので、 $S$ の方が大きくなるので、Rayleigh-Ritzをやれば出るということは変わらないんですが、イコールになるのは対角化可能なときだけと言っていいのかな。正確には $\Omega$ の中の固有値に対応する固有ベクトルのところで対角化が可能だったら大丈夫。
Ｃ：なるほどね。それが聞きたかったの。
今倉：ややこしいですね。
Ｃ：その領域の中で考えた行列という意味が正解。
今倉：ごめんなさい、そうです、そうです、そうです。
Ｃ：落ちた行列というか。
今倉：それが正しいです。行列全体はブロックが必要だけれども、Jordanブロックサイズが1だったらそのとおりです。そのとおりのご指摘で、そのとおりです。
Ｃ：よく分かりました。ありがとうございました。
Ａ：ちょうどいいぐらいの時間だね。では、どうもありがとうございました。（拍手）

クロージング

Ａ：では、これで今日はおしまいです。次回は後ろにいらっしゃる友枝先生に（笑）。2週間後にまたこの場所でやりますので、もしご予定がつけばいらしてください。どうもありがとうございました。（拍手）
Ｃ：この後、ご予定は？
Ａ：飲みに行ける人で飲みに行きましょう。

プロフィール

2011 年3 月名古屋大学大学院工学研究科博士課程修了。博士（工学）。筑波大学計算科学研究センター研究員を経て、現在、筑波大学システム情報系助教。大規模線形方程式、固有値問題およびそれに類する問題に対する数値解法の研究に従事。日本応用数理学会、情報処理学会各会員。

ホームページ：http://www.cs.tsukuba.ac.jp/~imakura/