二村 保徳

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

講演者：二村保徳（筑波大学）
題　目：モンテカルロ法に基づく超大規模疎行列の固有値分布計算
日　時：2015年4月28日(火) 18:00-
場　所：早稲田大学西早稲田キャンパス63号館数学応用数理会議室（102号室）

○司会　じゃあ、きょうは数理人の第6回目ということで、最近ちょっと来れてなかったので申しわけないんですが、櫻井先生のグループの最近助教になられた二村さんにご講演をお願いしました。タイトルは「モンテカルロ法に基づく超大規模疎行列の固有値分布計算」。櫻井研で、CRESTのチームでソフトウエアとか開発しているんですけど、ほぼ二村君が実装しているという感じで、内部の一番頑張っている方です。
　じゃあ、よろしくお願いします。
○二村　ご紹介ありがとうございます。筑波大の二村です。

　では、始めます。アウトラインはこのようになっていまして、最初、固有値問題のおさらいみたいな、どんな方が来られるかよくわからなかったので、そんな話と、その後にメインの話になる、固有値数のモンテカルロ法による確率的推定法、それから、途中で小規模問題での数値実験を挟みまして、それらの方法の実用化と、さらに並列化の話、そして、超大規模な問題での数値実験という流れでいきたいと思います。このように時間もある程度長いので、各章の切れ目で質問とかあれば、どんどんいただければと思います。

　今回の講演では、標準固有値問題を扱います。ご存じのとおり、この場合、 $n$ 個の固有値が存在しまして、今回、簡単な、 $A$ が実対称行列の場合を扱いますので、この場合は固有値が全て実数で、相異なる固有値の固有ベクトルが互いに直交します。今回は行列 $A$ が超大規模な疎行列である場合を踏まえて話を進めていきます。

　固有値問題といいましても、実際の応用ではさまざまな状況がありまして、行列の大きさがどこから大規模かというのもあるんですけれども、小規模のものから大規模のものまで。行列 $A$ の要素の、ゼロが多い疎行列だったり、ゼロがほとんどない密行列だったりというバリエーションがあります。
　それから、固有値問題として、固有値と固有ベクトルを求めるという話なんですけど、固有値だけが欲しいという場合もあったりします。
　また、実際全ての固有値を求める場合だけではなくて、絶対値の最大最小とか、単純に最大最小から幾つ欲しいといった要求があったり、ある特定の区間内のものを全部欲しいですとか、ある値から近いものを幾つか欲しいといった需要があります。

大ざっぱなんですけど、固有値問題の解法を大きく分けますと、直交変換に基づく方法と射影法の２つに大ざっぱに分けることができます。
　こちらの方法は、主に直交変換を使って、もとの行列 $A$ を三重対角──三重対角だったり五重対角だったりすることもあるんですけども、などの簡単な形に変換しまして、その変換後の行列の性質を利用して、QR法とか分割統治法とか二分法、逆反復法の組み合わせなどの反復計算を使って固有値を求めます。
　一方で、射影法と呼ばれる方法は、これはちょっと抽象的なんですけど、２つの部分空間 $\mathcal{K}$ 、 $\mathcal{L}$ というものをつくって、 $u$ が $\mathcal{K}$ に属して、 $A u - \lambda u$ がもう１つの部分空間 $\mathcal{L}$ に直交するような、そんな条件を満たす $\lambda, u$ を近似固有対とする方法です。
　部分空間の生成法や部分空間の選び方によって、さまざまな解法がありまして、Lanczos法だとか、Jacobi-Davidson法だとか、私がソフトウエアを実装していますが、櫻井・杉浦法といった方法などが属します。

　それぞれの方法群にはこのような特徴がありまして、直交変換に基づく方法は行列の要素を直接変換することを行います。射影法はなるべく行列要素を変換しない前提で進めていきます。ですので、こちらの方法は、密行列に対して主に使われます。一方、射影法は疎行列に対して使われます。
　また、こちらの方法は、計算量としては求める固有値に依存しない項が基本的に支配的なので、密行列の中でも特に全固有値を求める場合に有効です。一方、射影法のほうは、計算量として求める固有値の数に依存する項が支配的になりますので、特に一部の固有値を求める場合に有効になります。

　本講演では、指定したある区間内の固有値とそれに対応する固有ベクトルを求める場合について考えます。
　指定した値付近の固有値と対応する固有ベクトルを求める手法としまして、Shift-and-invert Lanczos法やJacobi-Davidson法などが挙げられます。これも指定した区間内にある値を指定しまして、その付近を求めるというふうにして、指定した区間内の固有ベクトルを求めます。
　一方で、直接指定区間内、厳密には複素平面上の領域を考えるんですけども、その固有値、固有ベクトルを求める方法として、先ほど述べました櫻井・杉浦法とか、FEAST、こちら、Referenceなんですが、などの方法があります。

　ある区間内の固有値を求めるんですけども、その区間内の固有値の数というのがかなり多くなってきますと、解法の中でグラム・シュミット直交化などの計算量が固有値数 $m$ の二乗に比例する部分の計算コストがきいてきて、この部分がつらくなってきます。素朴にその場合どうしたらいいかと考えると、求めたい区間を複数の小区間に分けて、それぞれの区間の固有値計算を個別に行うことが考えられます。さらに、異なる区間の計算は独立に、当たり前ですけど、行うことができますので、並列計算は可能になります。

　先ほど申しましたように、各小区間にもとのある区間内の固有値を求めるとした場合に、こんなふうにバッテンで固有値をあらわしていますけど、こんなふうに固有値が分布しているとしまして、それぞれの小区間に、各区間の固有値数が均等になるようにすると、計算量が均等になりやすいので、大体この状況だと、例えば何となくこんなふうに小区間に分けたいと。ですが、固有値の大体の分布がわかるとうれしいんですけども、当たり前なんですけども、固有値がわからなくて固有値問題を解こうとしているので、それがわかるならば苦労しないと。とはいえ、粗い近似精度でもいいので、何かしら少ない計算コストで固有値分布を計算することはできないものかという考えに至ります。

　実際、固有値の分布を求める方法として挙げられるのがこちらの３つの手段で、よく知られているGershgorin定理の利用で、固有値をGershgorin円と言われる円の計算をして、それで、固有値がどの辺にあるかという当たりをつけるんですけども、ちょっと問題点として、いつもそうとは限らないですけど、Gershgorin円が大き過ぎる場合があるという問題が考えられます。
　また、行列要素の絶対値をとって計算したりするので、行列が陽に与えられていない場合には使用が困難になると考えられます。
　もう１つがSylvester慣性則の利用も考えられるんですけども、こちらは行列 $A$ の修正コレスキー分解が必要になりますので、こちらも行列が容易に与えられない場合というのは難しくなります。
　ここでさっきから行列が陽に与えられない場合という話をしているんですけど、大規模問題だと、行列掛けるベクトルという作用だけしか得られないというか、行列を陽につくってしまうと、メモリが膨大に必要になったりする場合があって、そういう場合にはちょっと使用が難しい形になります。
　最後に、ここに出しています行列のTraceのモンテカルロ計算に基づく手法というのが今回のメインテーマでして、こちらはほかの方法に対して、行列が容易に与えられなくても使えるんですけども、弱点としては、モンテカルロ法なので結果は不正確、著しく不正確な場合も考えられます。今回はこちらの話をしていきます。
　ここまで何か。大丈夫ですね。

　じゃあ、ここから行列Traceのモンテカルロ計算に基づく固有値分布の計算法の概要です。

　これはちょっとどんなふうに話が進んでいくかという話なんですけども、固有値分布の近似法というのは、このような行列値関数、これはスカラーのときの表現なんですけども、Traceをとることによって固有値数を与える、こういった矩形の行列値関数を考えて、その行列値関数を行列多項式や有理関数で近似します。その近似した行列値関数のTraceを求めるのにモンテカルロ計算を使います。

　最初の話というのは、この $h$ というものを使って話が進んでいくんですけど、まず定義としまして、 $\Lambda$ は固有値を対角に持つ対角行列です。ここで $U$ というのが固有ベクトル $u_1$ から $u_n$ を並べた直交な正方行列と定義します。そうしますと、ご存じのとおり、行列Ａはこのように固有値分解されまして、 $U$ は直交行列になります。
　ですので、この $h(t)$ というのは、スカラーが入った場合に、スカラーが区間Ａの中に入った場合は1、それ以外の場合はゼロといった関数を考えまして、形式的に行列値関数として扱って、 $A$ を入れたもののTraceを考えますと、固有値分解がありますので、こういうふうに表現できまして、それで、この対角行列の中にある固有値の部分だけが残りますので、 $u_i u_i^T$ の中で対応する固有値が区間内にある部分だけのTraceの総和をとることになりまして、それは結局、固有値数 $m$ になるということになります。

　矩形関数 $h$ を近似して、実際使えるようにするんですけど、１つは多項式近似でやる方法が提案されていまして、*1の論文に載っています。
　もう１つが周回数値積分による近似を行う方法でして、こちらは私が過去に研究した内容になります。今回はこちらの方面でお話ししていきます。

　まず、ここで複素平面を考えまして、こちらを $a$ 、こちらを点 $b$ としまして、区間 $[a,b]$ がこちらにありまして、この $a$ と $b$ を通るような単純閉曲線 $\Gamma$ を考えます。そうしまして、このような周回積分を考えますと、 $z$ を変数としてこの $Gamma$ の積分をしますと、このように固有値分解を使って式変形しまして、それぞれ各固有ベクトルと各固有値の和の形でこれを書き直しますと、このようにあらわせます。そうすると、コーシーの積分公式で、この関数の極が $\lambda$ になりまして、留数が1になりますので、最終的にこの式というのがこの区間 $[a,b]$ 内の固有値数 $m$ となります。なので、この $h(A)$ はこのようにあらわせます。ただし、周回路が $a,b$ を通りますので、固有値が $a$ や $b$ になっている場合は、ここの部分のイコールは成り立たない話になります。

　この周回積分を実際に計算するために、数値積分で近似します。例えば台形則とかガウス・ルジャンドル則などです。一般にこのようにして周回積分を近似します。これで $w_j$ が重みで、 $z_j$ が積分路上の積分点です。
　こちらは、先ほどの行列値関数 $h(A)$ の有理関数近似としても見てとることができます。

　さて、先ほど、多項式で近似したり、有理関数で近似したりという話になってくるわけなんですけども、こちら、行列が大規模な場合というのは、ここに逆行列が出てきますので、行列多項式や逆行列を陽に計算することは非現実的ですね。ですが、これのTraceを計算するという話なので、そこで、行列Traceを求めるモンテカルロ系の方法を用いて計算量を削減することを考えます。
　こちらは、行列のTraceを求めるモンテカルロ系の方法というのは結構、このようないろいろな文献で扱われています。

　今の状況を整理するんですけども、多項式近似の場合は、実対称行列の実係数多項式は実対称行列になりますので、実対称行列のTraceを近似するという話になります。
　一方、数値周回積分の場合、固有値はこの実軸上にしかないので、計算上いろいろ都合が良いので、積分点を複素共役の対でとることが多いです。
　実際、そのような場合を考えると、ある $z$ というのが積分点だとして、それに対応する重みが $w$ だとして、先ほどのところに出ていましたように、この和を計算するので、それぞれ共役対になった部分和というのがこういうふうにあらわせるわけなんですけど、これが実対称行列になるので、結局、それらの和を計算することになっていますので、数値周回積分の場合も、共役対で積分点をとっていれば、実対称行列の議論に落ちます。なので、ここから話の問題設定としまして、実対称行列のTraceを近似する。もちろん、ただ行列 $A$ が陽に与えられていたら、他に対角を足し合わせるだけなので、もちろんそんなのは当たり前です。なので、陽には得られない行列のTraceを近似するという話を進めていきます。

　行列Traceのモンテカルロ計算なんですけども、フレームワークとしてはこういう形になっています。実対称行列のTraceをあるベクトル、 $v_1$ から $v_s$ のs個のベクトルを使って、このような形で近似するのが行列Traceのモンテカルロ計算の方法になります。ここで $v_i$ というのは、ある種の乱数でつくられたベクトルで、基本的には $v_i^T A v_i$ の期待値が $tr(A)$ になるようにつくられます。
　実際こういうふうな形になっていますので、行列 $A$ にベクトル $v$ を掛けるという式だけでできていますので、行列が陽に与えられない場合でも、行列 $A$ とベクトルの積が計算できれば使えるという話になります。

　この行列Traceを推定するモンテカルロ系の方法は、このように4種類ぐらい挙げられまして、それぞれHutchinson estimator、Gaussian estimator、Normalized Rayleigh商 estimator、Unit vector estimatorが挙げられます。こちらは、こちらの論文で詳しく述べられています。
　それぞれについて順番に説明していきたいと思います。

　Hutchinson estimatorというのを最初に紹介します。これが一番古くから使われている方法で、1990年、Hutchinsonによって提案された方法でして、各要素が等確率で1か-1のベクトルを使います。こちらが最初の提案に使われて、広く用いられていると言われています。こちらというのが、この $v^T A v$ というのを１つのサンプルでの分散を評価しますと、このようにあらわせまして、 $A$ のフロベニウスノルムの二乗引く対角、これ、 $a_{ii}$ というのが $A$ の対角要素です。 $a_{ii}$ の二乗を引いているので、これは非対角要素の二乗和ということになります。なので、非対角要素が小さいほど分散が相対的に小さくなります。この1、-1のベクトルを考えるとすぐわかるんですけど、 $A$ が対角行列の場合は、ベクトル1本で厳密にTraceが得られます。

　次はGaussian estimatorについてです。こちらは各要素が互いに独立に標準正規分布に従うベクトルです。標準正規分布というのは、平均がゼロで分散が1の正規分布のことです。これはちょっとおもしろい性質がありまして、こういったベクトルというのは、直交行列を掛けても各要素が互いに独立な標準正規分布に従うベクトルになります。なので、行列 $A$ の固有値分解を考えると、結局、こういうベクトルを並べた行列に $v$ を掛けたやつも、統計的には $v$ と同じ性質を持ちますので、こんなふうに書きまして、分散は固有値のみで決まります。どんな固有ベクトルを持っていても、固有値が同じなら、同じような統計的な振る舞いをすることが知られています。

　
　次はNormalized Rayleigh商 estimatorなんですけども、これは今までとこの後述べるやつに対して１つ抽象度が高い存在なんですけども、これは単純に各要素が互いに独立で、平均ゼロのある確率分布に従う乱数ベクトルを使う方法です。ただし、ここで $v$ はこのように正規化されているものとします。先ほど最初に言いましたHutchinson estimatorは、これの一例になります。ここで確率分布を特に規定しないので、分散が何になるか、確率分布に依存するとしか言いようがなくて、この中では１つ抽象的な存在です。

　
　最後に紹介しますのがUnit vector estimatorで、これはある種Deterministicというか、これはまず、1から $n$ の値の中で１つ、一様に乱択しまして、それの単位ベクトル $e_i$ をこの $v$ として使う方法になります。なので、結局、ある単位ベクトル $e_i$ を使って $e_i^T A e_i$ とやると、明らかに $i$ 番目の対角要素が出てきますね。なので、それを適当に、ランダムに各対角要素をとって足して、サンプルの数で割っている値の倍数ということになります。そういうふうにやっても期待値は $A$ のTraceになります。こんな方法なので、非対角要素は全く見ていないので、この分散というのは対角成分がどのぐらいばらついているかで決まります。対角要素が全て等しい場合というのは、当たり前なんですけど、ベクトル1本で真値が出ます。
　さらに、こちら、対角成分がばらついている場合はなかなかうまく出ないんですけども、 $v$ の要素をフーリエ変換などでシャッフルしてあげると、見かけ上の行列の対角成分を均等にするような処理をすることで、精度を上げるといった試みも提案されています。分散はこのようにあらわせます。

　
　これは本当に当たり前の話なんですが、一応述べるんですけど、ここまでは区間内の固有値の数の推定の話だったんですけど、分布を見たいということなので、素朴にある意味ヒストグラムみたいなものを計算すればいいので、もとの区間の $[a,b]$ というのを、特に事前情報がないので、基本的には等間隔に分割しまして、それぞれの小区間で固有値の数を求めて固有値分布とします。
　ここまでで1つ終わったんですけど、何か不明点とかありますでしょうか。

　
　じゃあ、進みます。先ほどの行列Traceの固有値数の確率的推定法の実際やったらどうなるかというのを見るために数値実験をします。

　
　まず、指定した区間内の固有値の数を推定します。ここでTrace estimatorは、全てHutchinson estimatorを使っています。Matrix Marketの行列を使っています。
　ここでまず先に、数値積分をやっている方法ですので、単純に行列のTraceをexactに求めて、積分点数がどういうふうに影響するかというのを見ます。その後に、こちらのTrace estimatorを使って、固有値数がどのように得られるかということを見ます。最後に、指定した区間の固有値分布、ヒストグラムみたいなものの計算を実アプリの行列で行います。

　
　まず、ちょっと行列の数が少ないんですけど、Matrix MarketのPLAT1919というのとBCSST13というものを使います。こちら、区間の $a$ と $b$ はこのようになっています。中の固有値、正確な固有値の数はそれぞれこのようになっています。

　
　まず、Traceをexactに計算して、数値積分の変化による挙動を見ます。

　
こちらは真値40なんですけども、最初は55とかの値が出ているんですけど、積分点をふやすとだんだん近づいてきますね。これはふやしていくと、ちょっとUnder estimateするような結果になっています。

　次、もう１つの行列なんですけども、これは下から近づいてきまして、真値11に大体30ぐらいで、これはかなりスケールが狭いんですけども、かなり合ってきています。

　今の結果の積分点数 $N=16$ のときの結果というのを今回真値としまして、サンプルベクトル数 $s$ の増加でどのように固有値数の推定が変わっていくのかというところを見ます。

　これがPLAT1919の結果なんですけども、この真値、40.6ぐらいのところでして、かなりもたついているんです。また、1回近づいたと思ったらまた離れたと、こういった挙動をします。全体的にこのサンプルの数 $s$ が20ぐらいからこのように、大体プラマイ1ぐらいの範囲には入っています。
○質問者　ごめんなさい。その前にお調べになった $N$ が16というのは何でしたっけ。
○二村　積分点の数ですね。
○質問者　今のグラフの横軸は。
○二村　横軸はサンプルベクトルの数ですね。こちらはTraceの推定をした……
○質問者　ああ、そうか、そうか。なるほど。
○質問者　このときの赤い線は何ですか。
○二村　これは向かっている真値ですね。なので、今、この問題でやっていて、もともと周回積分をやっているので、連続の周回積分だと、厳密にこの40という値になるわけなんですけども、それを数値積分で近似していて、その数値積分の積分点数を上げていって、どんなふうに、まず、周回積分にいかに近づくかというのを見ていて、大体こんなものでいいかということで16ということにして、16というのが今の40.5とか、そのぐらいなので、じゃあ、この積分、16点の数値積分に対してTraceの推定、１つ近似して、もう１つ近似しているので、これは１つ目の近似の結果ですね。さらにモンテカルロ法を使って近似するというのがこの話です。
○質問者　縦軸の意味も1回目と2回目で変わっているということですか。カウントが非整数というのは……。
○二村　いや、一緒ですね。非整数というのは、整数に丸めてもいいんですけど、ありのままを見せているだけという形ですね。もちろん丸める前でどのぐらい近いのかというのをやっぱり見たいじゃないですか、40.6なのか40.9なのか。なので、非整数で出しています。
○質問者　すみません。それって、300と書いてあるのは、サンプルベクトルを300個発生させて、それを１回やったときの値。
○二村　そうです、1回で。
○質問者　300個を何回もやっているわけではない。
○二村　やっていないです。1回だけです。もちろん1回だけなので、複数やった結果も見るべきとは思います。

　もう１つの行列は、今度はこのような振る舞いをしていまして、ここは真値がもとの11にかなり近い値に対して、最初14ぐらいからスタートして、ちょっとOver estimateした状態で、ずっと続いているといった結果になっています。実際のMatrix Marketではこんなような振る舞いをするというお話です。

　最後は密度汎関数計算という、ナノシミュレーションのアプリケーションで標準固有値問題があらわれますので、そちらであらわれる行列でテストをしてみました。行列のサイズは17万次元で、この問題は最初から1,020個の固有値数が必要な問題ですね。そのような固有値が存在する区間の固有値分布を推定してみました。
　この問題というのが、最初からずっと言っている行列が陽に与えられない問題でして、行列ベクトル積の関数としてのみ与えられる問題で、こういった問題に適用したくて実際やってきた感じです。

　パラメーターですが、数値積分の積分点数は8点で、サンプルベクトルの数は20と非常に少ない数を使って、固有値数を複数の小区間でやって、固有値分布を計算したのがこの結果になります。

　これ、横軸が、すみません、小区間の左から見た番号です。等間隔に小区間に分割しています。縦軸は固有値の数です。これはヒストグラムで、黒い線が推定値の値で、赤い線というのが厳密な固有値の数です。
○質問者　黒いのは整数で、赤いのは実数になっているんですか。
○二村　そうですね。赤いのが……
○質問者　黒いのは実数の値で、整数になりますよね。
○二村　赤いのが整数ですね。
○質問者　赤いのが整数。
○二村　赤が答えみたいな、赤ペンみたいな、そういうイメージ。
　こうやって見ていただくとわかるように、結構いい感じで出ているんじゃないかと思います。ここは多分一番差が大きくて、３ぐらいですかね、ずれているんですけど、大体、なかなか精度よく出ているので。これを使って固有値問題の解法のパラメーターを設定するには、悪くないんじゃないかなといった結果になります。

　これは先ほどのこの図を左から足し上げていった結果なので、こっちのほうが少しわかりやすい感じになるんですけども、左から足し上げていくと、このようなグラフになりまして、大体、これはスケールは結構大きいんですけども、100はいかないけど、大体50ぐらいですね、Over estimateしているような形になります。
　実際、この1,020番目の固有値がどこにあるかみたいなことも知りたかったりするので、こういった使い方をして、その番号の固有値がどこにあるか、最初がどこにあるか知らないとちょっと厳しいんですけども、そんなふうにして調べていくことができます。

　こちら、ちょっと手前みそなんですけども、もともとこの方法は非対称の固有値問題をもうちょっと一般の話でやっていまして、その論文を過去に私、書かせていただきました。さらに、非線形の固有値問題を扱った論文もございます。
　ここまでで、インターバルなんですけど、大丈夫ですかね。
○質問者　さっきの4種類あった中で、結局どれを使っていたんだっけ。
○二村　これは1個目の1、-1の。
○質問者　それは何か理由があるんですか。
○二村　はい。それしか使っていなかったのは、この実験をやっていた当時にそれしか知らなかったからで、結構ほかのやつは論文としては新しくて、2011年とか。
○質問者　今回対象としている行列がフィルターがかかった行列で、固有値がほぼゼロ、1になるとかという情報があるんだけど、それはガウシアンのやつと相性がいいんですか、悪いんですか。固有値のみに依存するという、ゼロ、1というのは分布が大きいという判断ですか。
○二村　ちょっとそれはこの話の範囲を超えるんですけど、ただ、一応そういう話も考えていて。確率不等式でどのぐらい押さえられるかというのがあって、それを使うと、行列のTraceの真値がもしわかっていれば、変な話なんですけど、かなりタイトなものが得られるという、Gaussの場合だと。ということで、結局、Traceがわかっていると、中の数がわかっていないといけないので、それの例えば下限とかをわかっていれば、そういう性質が使えたりするんですけど。ちょっと、今示している分散の値の話だけだと議論はそこまで行かないんですけど、確率不等式の話とかもし出すと、いろいろおもしろい結果が得られる。Hutchinson以外のestimatorでもいろいろやってみたいと思っています。

　では、次に進みます。

　先ほど、この周回数値積分による方法というのが、数値積分とTraceのモンテカルロ計算をやることで、このようにして固有値数を求めるわけなんですけども、ここを見ていただくとわかるように、行列のInverseが出ていまして。なので、 $z$ が、 $j$ が1から $N$ までで、右辺としてこの右辺 $v_i$ の $i$ が1から $s$ まであるので、結局、 $s$ 掛ける $N$ 個の線形方程式を解かないといけないという難しさがあります。
　例えばこの線形方程式をLU分解を使った直接法で解くことを考えると、そもそもこれ、LU分解しているぐらいなんだから、修正コレスキー分解をやって最初に言っていたSylvester慣性則を使えばいいじゃんという、という話になります。なので、基本的にこの線形方程式を解くときは、Krylov部分空間法の反復法を使うことを考えます。

　今回、この行列 $A$ が実対称としていますので、これ、 $z I - A$ というのは複素対称になります。なので、複素対称行列向きの解法を使うことを考えます。複素対称行列向きのKrylov部分空間法は、例えばCOCG法とかCOCR法、QMR-SYM法などが挙げられます。今回、COCG法を使って話を進めていくことにします。ここでCOCG法を選んだ理由というのが、ベクトル和の計算が一番少ないからです。

　こちらが前処理なしのCOCG法のアルゴリズムです。これの問題として、複素対称行列Ｃに関する線形方程式、 $C x = b$ を解くアルゴリズムとして表現しています。小文字の太字はベクトルで、小文字の太字じゃないギリシャ文字がスカラーで、これは $C$ だけですけど、大文字が行列になります。この後のアルゴリズムの説明でも全部そういったノーテーションでやっていきます。

　これを見ていただくとわかるように、毎回、反復ごとに行列 $C$ とベクトルの行列ベクトル積があらわれます。こちらのCOCG法のKrylov部分空間反復法では、前処理を行わない場合ですと、行列 $C$ に関する疎行列ベクトル積が計算の主要部になります。こういった例外もあります。

　今回解きたい線形方程式というのは、ある右辺ｖで固定しますと、 $z_j I - A$ で、 $z_j$ が、 $j$ ｊが1から $N$ まで変わるような複数の方程式を解くという話になっています。この $z_j I - A$ というのは、この行列 $A$ を $z_j$ だけずらした、符号も変わっていますけど、固有値をずらした行列になっていますので、シフト線形方程式と呼ばれています。
　こういった問題を複数解きたいので、ここでKrylov部分空間のシフト不変性に着目します。Krylov部分空間はこのようにしていまして、このように行列 $A$ に対してスカラー掛ける単位行列 $I$ が足し込まれていても変わらないといった性質になります。

　
　このシフト不変性を利用して、効率的に複数のシフト方程式を解くという解法が幾つも提案されていまして、例えばShifted CGとかShifted COCG、Shifted COCRなどが挙げられます。

　
　先ほどCOCGを使うという話でしたので、Shifted COCGについて見ていきます。Shifted COCGというのは、後でアルゴリズムを示すんですけども、シフト線形方程式、 $\sigma_j$ というのが複素のスカラーなんですけども、このようなシフト線形方程式群をそれぞれに対して、それぞれCOCG法を適用したものと数学的に等価のアルゴリズムになります。
　ですが、ここで特筆すべきなのは、 $N$ のシフト線形方程式を解くのに、1回反復当たり、 $C$ に関する行列ベクトル積がたった1回で済むという性質があります。

　
　こちらがShifted COCGのアルゴリズムを示していまして、こちらが、省略してあるのが、先ほど示したCOCGのループの中身の部分になります。この中で１回行列ベクトル積を行っています。
　一方で、各 $j = 1$ から $N$ に対して、このように各シフト方程式の解ベクトルを漸化式で更新している部分があります。こちらは全部スカラーの計算で、ベクトルの計算はこちらに出ているのがわかると思います。なので、この $j$ に依存する部分というのは、一切行列ベクトル積は出てこないというアルゴリズムになります。ですが、この近似解というのが、COCG法を全てのシフト線形方程式に対して適用したのと数学的に等価になります。

　これを図示すると、COCG法をそれぞれのシフト線形方程式に別個に使った場合というのは、例えばこの青の部分が行列ベクトル積とか内積の計算で、赤の部分がベクトルスカラー積、ベクトル和の計算だとすると、こんなふうにあらわせるんですけども、これはShifted COCG法を使うことによって、1回分の行列ベクトル積と1個分の行列ベクトル積と内積だけでよくて、それ以外の部分だけこうあるという感じになります。

　このShifted Krylov部分空間法というのが万能というわけではなくて、基本的にKrylov部分空間法では前処理というのをやって収束性を改善していきます。ただし、一般に前処理を適用しますと、さっきから使っていましたKrylov部分空間法のシフト不変性は保たれないので、Shifted COCG法のシフト不変性を用いた解法では前処理の適用が一般には困難になります。

　さらに話を拡張しまして、今回、さっき、ベクトル1本だけで話していましたけど、もともとは複数のサンプルベクトルに対しての線形方程式を解いてあげる必要があるので、この複数あるというのをうまく使うことを考える話があります。というのが、この複数の辺ベクトルをまとめて扱うことで、全体、トータルの反復回数を減らすことが可能な、いつもとは言わないんですけど、いわゆるBlock Krylov部分空間法を使います。
　こちらのBlock Krylov部分空間法にもシフト不変性がありますので、シフト方程式向けの解法が幾つか提案されています。例えばShifted Block CGなどがあります。
　話はここまでですね。何かわからない点とかありますでしょうか。

　では、次は並列化の話をしていきます。超大規模問題を解くという話をしていますので、並列計算というのは避けて通れないわけです。

　並列計算を大ざっぱに分けますと、共有メモリ型並列と分散メモリ型並列があります。分散メモリ型並列の下で共有メモリ型並列が動きますので、今回は共有メモリ型並列というのは扱わないで、分散メモリ型並列の話を本当に少しだけします。
　分散メモリ型並列で今回のような計算をすることを考えますと、行列やベクトルのデータを複数の計算モードに分散しておきます。それは計算量の問題だったり、メモリの制約の問題だったりして、そういうことを行う必要があります。実際、計算にはインタラクションがありますので、逐一必要になったときに、ほかのノードにあるデータを通信ライブラリーを使って通信を行って、データを得る形になります。

　Shifted COCG法の計算カーネルは、この疎行列ベクトル積、内積、ベクトル和、ベクトルスカラー積です。それぞれ並列実装について、当たり前の話かもしれないですけど、一応述べておきます。

　並列実装する場合はこんな感じになります。疎行列ベクトル積は、基本的に近接通信になります。というか、超大規模並列計算をする場合、近接通信になるように、もとの微分方程式を離散化しようとするので、多くの場合近接通信になります。
　一方で、内積というのがそれぞれのプロセスごとの部分和を計算しまして、通信ライブラリーでそれらを足し合わせるということをします。ここの計算というのは、疎行列ベクトル積の近接通信と違って、全てのプロセスがかかわりますので、これは結構オーバーヘッドが大きくて、なるべく避けたいという事情があります。今回のテーマにはあんまり関係ないんですけども、なるべく内積を減らすようなKrylov部分空間法のアルゴリズムを考えるみたいなテーマもあります。
　最後に、ベクトル和、スカラー積、これは一切通信なしで計算を進めることができます。

　最後に、ちょっと並列計算からそれちゃうんですけども、先ほど線形方程式を求めるということを話していましたけど、本当に計算したいのは解ベクトル、 $x_{\sigma_j}$ そのものではなくて、 $v^T A v$ みたいなのを計算しているので、左から、線形方程式で言うところの、右辺ベクトル $b$ との内積を計算したいという、もともとそういう話です。なので、この解ベクトルを求めるということよりも、このスカラー量、内積の値を直接計算することを考えます。

　これは先ほどのShifted COCG法のアルゴリズムなんですけども、ここで、シフト線形方程式の解ベクトルの更新のところを見てみますと、先ほど示しましたように、ベクトルの漸化式があらわれています。簡単な話なんですけども、このベクトルの漸化式というのを抜き出してきまして、ここで示しているんですけども、この漸化式の両辺に $b^T$ を左から掛けるので、それぞれ $b$ との内積をとりますと、このように、これはそのままなんですけども、こちらは $b$ と残差 $r$ が直交するという性質を利用して、ここの項を消しているんですけども、こんなふうな漸化式。本当はもっと簡単になるんですけど、大体こんな感じにして漸化式を組めます。こちらはベクトルの漸化式だったんですけども、スカラーの漸化式にして使うことができます。

　なので、一切 $x_k^{\sigma_j}$ とか $p_k^{\sigma_j}$ という、（ $r_k$ は必要なんですけど）それぞれのシフト方程式に対して用意していたベクトルというのを、計算も要らないのはそうなんですけども、メモリも必要ないというようなことになります。なので、計算量、メモリ要求量の大幅な削減ができます。これはこのシフト方程式の数が本当にすごくたくさんある場合、数千とか数万とか数十万とかあるような場合、ものすごくこれがきいてきます。

　さらに言いますと、さっき言ったような超大量のシフト方程式をまとめて解くことを考える場合というのが、特に、先ほど話しましたBlock Krylov部分空間法への一般化をするとそうなるんですけど、漸化式の係数の計算というのが実は計算の大部分を占めるような状況になります。というのが、漸化式の係数って、 $\alpha_k$ 、 $\beta_k$ とかのことです。ここの部分ですね。この辺の計算ですね。これがShifted COCG、スカラーのうちは全然問題ないんですけど、Block Krylovという、右辺をまとめて解くので、まとめて解く右辺の数というのが、いつも $L$ と言っているんですけど、このブロックサイズ $L$ の $L \times L$ の正方行列になるんですね、これが一般化されると。なので、 $L \times L$ の正方行列の掛け算を結構たくさんやる。これは $n$ がものすごくたくさんあると結構ばかにならなくて、しかも並列計算できないので、冗長にやらないといけないところなので、かなり重くなるんですね。ですが、今回スカラーの漸化式にしたことで、Block Krylovはスカラーではないんですけど、漸化式にしたときに、この計算というのが $j$ に対して並列に行うことができるんですね。なので、ここを並列にやって、最後に結果を集めるという。残差が十分小さくなったときに、最後に結果を集めるということをすることができるので、それを使うととてもよいという話になります。

　最後、大規模問題の数値実験の話です。ここまで何かございますでしょうか。

（注：まだPublishされていない結果のため、スライド、および文字起こしは非公開になります。）

　こちらで発表は最後なんですけども、モンテカルロ法に基づく固有値分布計算手法とその実用化と並列化について述べました。電子状態計算にあらわれる超大規模の問題で、固有値解法を走らせる前のプロファイリングとして十分高速ではないかということを確認しました。
　今後は、先ほどちょっと質問の答えの中で出たんですけど、どのぐらい真値に近いかというのが興味のあるところだと思うんですけど、これはある程度確率不等式で押さえられるような話がありますので、この辺を研究していきたいなと考えています。
　以上です。ご質問等あればお願いいたします。
○質問者　数値実験の２つ目の図、結果をもう一回見せてもらっていいですか。
○二村　これですか。これは元データがあります。拡大しますか。
○質問者　これは、細か過ぎてあんまりよくわかっていないんですけど、幅になっているわけじゃなくて、線が単純に結ばれているというだけですか。これはヒストグラムみたいになっているはずなんですよね。
○二村　これはヒストグラムなんです。
○質問者　何で下のほうはないんですか。ちょっとちんぷんかんぷん。
○二村　ヒストグラムというのは、書き方として、大分前のスライドなんですけども……。
○質問者　あれ、だから、要は区間になっているわけじゃないですよね。
○二村　こうやって出しているんですよね。これがもっと少ない場合で、これを出しているんです。これが今、めちゃ細かくなって、今みたいになっているので。ここは下、続いていないですね。下までいっていないものになります。
○質問者　ああ、なるほど。そうか、そうか。わかりました。
○質問者　これ、細かく区間を区切って、合計を足し合わせたのと、全体を一区間と思って計算したときで、それはずれますか。
○二村　ずれますね。
○質問者　どっちが正しいとかというのは。
○二村　細かいほうがいいはずですね。
○質問者　積分点数を同じぐらい大量にとっても。例えば積分を無視したとして、モンテカルロの意味だけで……
○二村　同じ数の、前、いっぱいあったほうがいいという着想に至ったんですけど、その根拠を、すみません、ちょっと忘れてしまったので。もちろん積分点が──どういう積分をとるかもありますけどね、全体をとったときに。こういうふうに積分点を置いたら、置いたほうが精度が悪そうですよね。端のところがやっぱり──等間隔に置けばですよ。端のところが結構疎になってしまう。
○質問者　ちょっと計算量とかとりあえず無視したとして、区間があって、全体で１つの区間を囲って積分した値と、それを単純に２で割って、２つ円を描いて、それぞれ求めて足し合わせた値というのはどれぐらいずれる。積分の精度と言ってしまえばそれまで。
○二村　そうですね。
○質問者　フィルターがなまって、外側の部分を取り込んで二重にカウントしていたりしたらとか思ったんだけど、そういうわけではないのか。
○二村　フィルターの裾を考えて、微妙に区間をふやしながらやってもらうやり方はあるんですけど。
○質問者　積み上げていって、最後にここのところで50個ぐらいずれたというのが、全体で囲ったときもそれぐらいずれるのかという。やってみないとわからない。
○二村　ちょっとまた積分の誤差とモンテカルロの誤差が両方あるので、どう考えていいんですかね。その辺はおもしろいですけどね。
○質問者　どういう問題が苦手なんですか。
○二村　問題ですか。それは例えばモンテカルロの意味で。
○質問者　そうですね。
○二村　モンテカルロの意味で言いますと、苦手かどうかは、どのestimatorを使うかというところに来るんじゃないですかね。なので、Hutchinson、今回使ったやつだと、対角要素に比べて非対角要素がでかいと、かなり嫌なんですよね。非対角要素の中に埋もれてしまってというんですかね。ということが起きるので、嫌だというのはありますね。Gaussianだとその限りではないんですけど。Unit vectorを使った場合だと、例えば非対角要素は関係ないので、対角要素がどれだけ平均が、平均というのは、分散が小さいかだけで決まるという話になりますね。非対角要素といっても、行列多項式の意味でとか、行列の数値積分した後の要素に有理関数の行列がどうなるかという話なので、もとの行列 $A$ ではないので、また実際にやるときにわかるかというのは、そんな簡単じゃないから。
○質問者　そうすると、じゃあ、今目指しているのはやっぱり、高速にということは一番主眼としてあるということなんですかね。
○二村　はい。この流れだとそうですね。ただ、最後にも言いましたように、どれぐらい真値に近いかというのを確率不等式で押さえるというのも興味あるところになりますね。それは大規模の問題でも、やった後に、このぐらいは合っているかもよというのは、確率的ですけど、一応出せるので。
○質問者　なるほど。それはおもしろそうですね。
○質問者　どれぐらい合っていてほしいものですか。どう使うかは、例えば今想定しているのは多分、櫻井・杉浦法の前処理として使ったとして、どれぐらい合っていれば満足できるでしょう。
○二村　例えば100だったのが、200は超えてほしくないかなという。
○質問者　でも、それぐらいの、1割ぐらいずれていても大したことはない。
○二村　1割は全然余裕だと思うんですね。倍とかになっていても、まだ許せるかなというぐらい。倍でも、もとの数にもよりますよね。5万だったら1万とか。そのサイズは見ることはないですけど、もうちょっと小さい値を、固有値数が見たいですけど。
○質問者　それは大き過ぎてestimateが外れると、SS法のどこで問題となるんでしょう。
○二村　Under estimateは困りますね。Over estimateはまだいいんですよね。
○質問者　そうですよね。Over estimateだったら別に、さほど大きな問題にはならない。
○二村　ちょっと計算時間がかかるだけで、答えは出るあれなんですけど、Under estimateは怖いですね。
○質問者　じゃあ、今、例えば2倍になっちゃったのが怖いとおっしゃっているのは、時間的な意味合いでということですね。
○二村　そうです。
○質問者　じゃあ、別に計算できないと言っているわけじゃない。
○二村　そうですね。Under estimateだったらもっと怖いですね。今回この方法は、Over estimateかUnder estimateか多分わからない方法になるので、恐らくかなり困難だと思うので、そこは難しいところではあるんですけど。
○質問者　なるほど。
○二村　何もわからないよりははるかにましというところからスタートはしているので。
○質問者　SS法でUnder estimateだった場合の対処ってあるんですか。
○二村　一応ありますね。Under estimateにしていて、それでももとのパラメーターを、部分空間サイズをどこまでとるかというようなパラメーターなんですけど、それの最大値をUnder estimateに従ってあまりにも小さくすると、無理ですね。
○質問者　じゃ、そっちを解消しちゃうと。何回もやってあげれば大丈夫。
○二村　はい。どのぐらいまでふやすかみたいな、それはソフトウエアの話になっちゃうんですけど。無限にふやせばいいじゃんということは別に、数学的にはそういう話がありますけど、具体的には何か上限を決めなきゃいけないので、そこがあまりに小さいと無理ですね。
○質問者　なるほど。どのEstimatorを選ぶかというところもまた、何かパラメーターは、人が選ぶのではなくて、自動化的な意味ではあるんですか。
○二村　行列の要素が得られるならできると思うんですけど。それは今回のテーマでは無理という話から来ているので、そういう場合はもうアプリケーションで、これでやっていて、今まで小さい問題でうまくいっているから、こういう傾向があるからこれを使うかなというぐらいに、どうしても本当に実際やるときはそういう形になると思いますね。
○質問者　なるほど。
○質問者　Unit vector estimatorは、この使い方は厳しいんじゃないかなと思って。だから、固有値の値みたいな──違う。固有値じゃないのか。対角要素を返すのか。厳密には固有値じゃないのか。
○二村　そうですね。
○質問者　違うね。関係ないか。
○二村　何かの変換をかけて、直交変換をかけて、対角要素をなるべく均等化するようなアプローチもある。これだけでなく、これを出発点にして、いろんな方法があり得るという。
○質問者　今、ソフトウエアを作成されていて、実際に、理想というのはかなり難しいですけど、できばえとしてご自身ではどのくらいだとお考えなんですか。
○二村　この推定の結果のできばえ。
○質問者　そうです、そうです。もう既に実装されているわけですよね。
○二村　はい。一応実装はしています。たくさん、分布を計算するんじゃなくて、結局、櫻井・杉浦法って周回積分で、これも１つの周回積分の方法じゃないですか。なので、途中まで同じ計算でできるけど、櫻井・杉浦法のサンプルベクトル、ここに出ていたやつなので、部分空間をとってくるもとのベクトルをこのサンプルベクトルにすれば、櫻井・杉浦法の副産物としてこの推定値が出るので。その副産物を使ってパラメーターを参考に決めるというのは実装していますね。そのベクトルにHutchinson vectorを使っているんですけど、これでいいのかというのはまた……
○質問者　また１つ議論。
○二村　議論はあるんです。
○質問者　それを使ったときに苦手となるようなものが出てきてしまうというのは、行列の性質から何かあったりするんですか。
○二村　そうですね。Traceを計算する行列が非対角優位だと、大きいとなるというのと、あと、何か嫌な理由はあるかな。
○質問者　それって対角要素との相対的なものということ。
○二村　はい。1、-1のベクトルで挟むと、同じ足というか、要素のところ同士で対角をとってきますよね。なので、対角に当たったところは、1、-1のベクトルなんですけど、同じインデックスのところなので、1でも-1でも、1掛ける1と-1掛ける-1になるので、必ず1になって、対角は絶対残るんですよね。それに対して非対角は、1、-1でサンプルをふやしていくと、非対角要素が例えば $a_{ij}$ だとすると、 $a_{ij}$ 足す $a_{ij}$ 、足す $-a_{ij}$ 足す $-a_{ij}$ と、大体平均としてプラマイゼロになる。でも、非対角要素が相対的に対角要素に対してでかいと、それだけぶれたときの影響が大きいですよね。なので、非対角要素が対角要素に比べて相対的に小さいほうがいいという話になります。それは分散の計算をすると、こういうふうになって、分散の式にもあらわれてくる。
○質問者　幾つか数値実験、この間も見せていただきましたけど、そのときに使われているテスト行列というか、使っている行列というのは、全部の固有値が単純であるとか、もしくは重複固有値を含む含まないとか、あらかじめわかっているものを使っているんですか。
○二村　重複度に関してはほとんどわかっていないんですけども、数値的にこの問題とかは、数値計算すると、数値的にある程度縮重している問題──これじゃないですね。最後の大規模な数値実験の１番目とかは、数値的に結構縮重しているというのがあるのと、あと、こういうシリコンがあるんですけど、その結晶構造が、同じものが周期的に続いているみたいな構造をしていまして、そういうのって数学的に縮重するような問題になるみたいで。実際これで見ると、数値的に縮重しているという……
○質問者　要するに、重複している固有値がある可能性があるというか、多分あるという。
○二村　きっとあるという。
○質問者　その場合は、Exactな固有値の数というのは、グラフに出ていたらどういう意味で書いていらっしゃるんですか。
○二村　Exactは、もともとこのアプリケーションの中のコードに実装されている固有値解法で、十分残差が小さくなるまで求めた結果になっていまして、縮重している固有値も求まるような方法になっているので、それに関しては大丈夫と思います。
○質問者　大丈夫というのは、例えば3重複だった場合に、それを3個とカウントするのか、1個とカウントするのか。
○二村　３個とカウントするようにしています。全部この発表の中での固有値数というのは、重複しているものをちゃんと数えた結果になっていますし、方法としても、縮重は、重複していても大丈夫な方法になっています。
○質問者　需要として固有値の数を知りたい場合に、それは縮重している場というのは、3重複を1と考えるのか、それともばらして数えるのが需要としてあるということでいいんですか。それとも、一応同じ固有値だった場合は1個とカウントするほうが望ましいといえば望ましいのか。
○二村　重複していても、それぞれ別なものとして見たいですね。特に、求めた固有ベクトルの構成が大事になるアプリケーションなので、それはやっぱりちゃんと分離して、重複している分の固有ベクトルも直交化したいですよね。という意味では、別々に数えてあげたほうがいい話になってくるかと思いますね、この問題ですと。
○質問者　わかりました。ありがとうございます。
○質問者　アプリで、下から10番目と言われたら、10というのは、ばらばらなやつの10番目──ばらばらじゃないや。ごめん。重複値も含む形の10番目を探してくださいという意味。
○二村　はい、そうです。
○質問者　10位タイ的な。
○質問者　10位タイになるかもしれないけども。あと、これ、非対称にいったときは、どこに困難さが出てくるんですか。
○二村　非対称といいますと。
○質問者　非対称の固有値問題にいったときは、このアルゴリズムのどこに難しくなってくるのか、特に何も考えなくていけちゃうのか。
○二村　非対称の場合、今回対称にやったのは、先の話が対称ばっかりだったということなんですけど、非対称で、このTrace estimatorが素朴に使えるかどうか、ちゃんと確認というか、自分で証明していなくて。というのは、Hutchinsonは非対称でできるのは証明したんですけど、結構対称の議論しかしていなくて、Hutchinsonはいけるのはわかっていて、Gaussianもできると書いてあって、この辺とか、Unit vectorはいいんですけど。Gaussianが非対称でできるよと書いてあるだけで、ちょっと証明とかは書いていなかったので。この辺は気持ち悪いところがあるので、それを対称にしたというのは、非対称になったときに、この辺が──この辺ですね、特に。どの辺までVarianceが幾つになるかというのが、一緒になるかというのが、そこまでちゃんと絞り込めていなくて。なので、非対称の難しさというか、まだこいつが全部使えるかどうかわからないという。この辺は研究の余地はあると思います。
　あとは何か、非対称とか。重複でも別に問題なくできることはわかっていますし。
○質問者　対角化は関係ない。
○二村　対角化が可能かどうかは関係ないです。
○司会　ほかの人、よろしいですか。
　じゃあ、きょうはありがとうございました。（拍手）

プロフィール

2014年3月筑波大学大学大学院システム情報工学研究科コンピュータサイエンス専攻博士後期課程修了。博士（工学）。筑波大学システム情報系研究員を経て、2014年10月より現在、筑波大学システム情報系助教。大規模連立一次方程式および固有値問題の並列解法、並列数値計算ソフトウェアの研究に従事。日本応用数理学会、情報処理学会各会員。
ホームページ：http://www.cs.tsukuba.ac.jp/~futamura/