第２３回助教の会

2013年度初回となる第23回助教の会では，数理第六研究室特任助教の森野佳生（もりのかい）さんに，「振動が失われたネットワークにおける効率的な大域的振動の回復」というタイトルで発表して頂きました．森野さんは合原研でこの3月に博士を取得され，4月に特任助教として数理六研に着任されました．非線形動力学や疾患の数理モデルを専門として研究されています．

講演の内容はまだ論文投稿中とのことですので，後日公開させて頂きます．

 ========

論文が出版されたということですので，2014年4月11日に内容を追加しました．

今回の講演では，振動が止まってしまったネットワークにおいて，外部から振動可能な素子を加えることによって再びネットワークの振動を引き起こす，という現象の数理モデルが扱われました．このモデルは再生医療を動機としたものになっています．

再生医療とは，機能を失った臓器に対して，何らかの細胞を移植することによって元々の組織の治癒能力を高め機能を回復する手段のことです．近年では，移植の際の副作用が小さい手法として，シート状に並べられた細胞を元の臓器に貼り付けることによって臓器の機能を回復する手法が研究されています．例えば，心筋細胞から作られているシートでは発火が起こりますが，シートの各部分の発火のタイミングは，時間の経過に伴って同期するということが知られています．今回の講演では，各細胞の振動が Stuart-Landau (SL) oscillator と呼ばれるモデルで表せると仮定し，振動がどのように時間発展していくのかを数理的に解析しています．すなわち，再生医療の視点から見ると，弱った機能が回復する過程をモデリングして解析しています．

各細胞の振動が SL oscillator で表されるとすると，細胞が active か inactive か（正常かどうか）を一つのパラメータの正負によって区別することができます．今回は，振動子は active なものと inactive なものの2種類のみであると仮定し，各振動子の間で相互に影響を及ぼしあう場合を考えています．これらの仮定の下では，inactive な振動子の割合がある閾値を超えると，全体の振動が止まってしまうことが先行研究により確認されています．

さて，振動が止まってしまった状態から，図のように active な振動子を加えて振動を回復させることはできるでしょうか？

本研究の1つ目の成果は，振動を回復するために外部に active な振動子をどれくらい加えれば良いのか，を理論的に解析したことです．振動子の数は大量にありますので，単純な計算では解析することができません．そこで本研究では，「同じ性質を持つグループに属する振動子は同期している」という妥当な仮定を加えたうえで，全体の振動が0かどうかの判定を原点の安定性の判定に置き換える，という工夫をして解析をしています．2つ目の成果は，効率的な振動子の加え方を解析したことです．理論的な解析によって，新たに加える振動子は active なものに優先して付加していけばよいという結果を得ています．最後に3つ目の成果として，2つ目の成果として述べた最適な方法で振動子を加えていったときに，振動の回復が起こる過程を初期状態によって4つに分類しています．

今回の講演では，再生医療という具体的な動機から，数理的なモデル化とその理論的な解析に至るまで詳細に話して頂きました．「現実的な問題を表現すること」と「理論的に解析しやすい形にすること」とを両立することがモデル化の難しさであり，面白さでもあると感じました．

第３０回助教の会

2014年度初回の助教の会では，数理２研の小林が発表しました．開始以来だんだんと長文になってきていたこのブログですが，担当者の負担が大きいこともあり，本年度から本人が書く形式にすることになりました．今後ともよろしくお願いします．

さて，今回は「円板形領域損傷モデルにおける 最大流最小カット定理と高速アルゴリズム」というタイトルで発表いたしました．ネットワークの信頼性の指標として用いられる「グラフの連結度」は，グラフ理論や最適化の分野で盛んに研究されています．連結度は，「何本の辺（点）を取り除くとグラフが連結でなくなるか？」を表す量であり，辺や点を取り除くことは，ネットワークにおいてリンクやノードが故障することを意味しています．

今回は，道路網のように平面上に埋め込まれているグラフを考えます．そして，個々のリンクやノードの故障を考える代わりに，自然災害や事故のように，損傷が地理的な広がりを持った領域で起こるモデルを考えています．具体的には，損傷がおこる領域は円板形で表され，損傷がおこるとその領域内のリンクやノードがすべて使えなくなるという設定です．

このモデルにおいて，いくつの損傷が起こると連結性が壊れるかを求める問題（最小カット問題）と，その双対にあたる問題（最大流問題）を扱ったのが今回の発表の内容です．本研究では，これらの問題に対して，最大最小定理を与え，高速なアルゴリズムを構築しています．

同じ内容で講演したものがこのページにありますので，興味のある方はご覧下さい．

数理第二研究室 小林佑輔

第２９回助教の会

久々の開催となりました第２９回助教の会では，本年度秋に数理４研に着任された松井千尋さんに「確率過程の可解性と代数構造」というタイトルで発表していただきました．松井さんはもともと物理学の出身で，今回の発表はこれまで行ってきた研究を物理に馴染みのない人でも分かるように構成したものとのことです．

統計物理では粒子同士の相互作用といったミクロな規則から統計的手法を用いて例えば超流動といったマクロな現象を導き出すということがテーマとなっています．例えば磁性体のモデルであるイジングモデルの場合，それまで揃っていたスピンの向きが系の温度を上げるにつれて熱揺らぎの効果が大きくなっていき，ある温度を境にスピンの向きが全く無秩序になるという臨界現象が知られています．このような臨界現象では，モデルの詳細よりむしろ「現象がモデルのスケール（例えば格子の数）に対してどのように依存性をもつか」といったマクロな特徴量が本質的となり，そのような量を導出することが物理現象を分類・解析するうえで重要な問題となります．

このような解析では近似や漸近的な議論が主に用いられますが，代数的によい構造をもつ「可積分系」とよばれる一部のモデルでは厳密解を求めることができ，そのような可積分系に対して松井さんはこれまで研究を行ってきました．今回話して頂いたのは，非対称単純排他過程(ASEP)をUq(sl2)代数とよばれる構造を用いて解くというものです．ASEPというのは１次元格子上に配置された粒子のモデルで，１個の格子（サイト）には１個（あるいは一般に$l$個）の粒子しか収められないという制約のもとでの粒子の動きを表現するものです．これはもともとRNAの動きのモデルとして提案されたものですが，渋滞のモデルとして聞いたことのある人も多いのではないでしょうか．

さて，ASEPにおける状態の時間変化は遷移行列を用いて表すことができます．具体的には，２サイト$(i,i+1)$間での粒子あり/なしの組み合わせそれぞれの確率を表すベクトル
$$
|P \rangle= |(\mbox{なし・なし}), (\mbox{なし・あり}), (\mbox{あり・なし}), (\mbox{あり・あり})\rangle
$$
に対して，その時間変化はある行列$M_{i,j}$をかけたものと表すことができ，それらを全てのサイト$i$について重ね合わせたものが全体としての状態変化となります．

この系の定常状態を求めるには
$$
\frac{\mathrm{d}}{\mathrm{d}t}|P\rangle= M |P\rangle
$$
の固有状態を求めればいいということになりますが，これを１サイトに収容可能な粒子数を$l\ge 2$個とした場合に拡張するのは簡単ではありません．

そこで用いるのが$M_{i,j}$のもつ代数的な性質です．$M_{i,j}$は
$$
e_i^2=(q+q^{-1})e_i\\
e_i e_{i+1}e_i=e_i\\
e_i e_j=e_j,e_i,\,\mbox{if } |i-j|>1
$$
という３つの代数関係を満たす行列$e_i$(Temperley-Lieb (TL)代数生成子)の直交変換で表すことができることが知られていて，このTL生成子を３状態以上に拡張したものをうまく組み合わせることで，一般の$l$における遷移行列$M$を表現することができます．

$M$をTL生成子で表すことのメリットは，TL生成子$e_i$に対して演算子$X$が可換となる(つまり，$e_i X=X e_i$となる)ような$X$が知られているということで，初期の定常状態$|P\rangle$に対して$X|P\rangle$もまた定常状態となることから次々と定常状態を求めていくことができます．このような性質をもつ$X$が冒頭に出てきたUq(sl2)代数生成子で，松井さんはこのことを用いて「粒子が右に進む確率が左に進む確率より大きい，粒子は計$n$個で境界外への出入りはない」という設定のもとで，「粒子のほとんどは右側$n/l$サイトに集中する，それより$r$サイト左に粒子が存在する確率は指数的に減衰し，その指数はサイトあたりの最大粒子数$l$に比例する」という非常にシンプルな結果を導出しました．

また発表の最後には，別の初期条件・境界条件に関する議論もありましたが，こちらは公表前の話題ということでここでは割愛させて頂きます．

私自身も研究で確率過程の挙動を調べることがしばしばありますが，基本的に近似と漸近論で済ます場合がほとんどなので，このように一見複雑な系でも代数構造を利用することで厳密な挙動が求まるというのは大変興味深い結果でした．

Multi-state extension of asymmetric simple exclusion process from Chihiro Matsui

第２８回助教の会

今回は数理情報学第１研究室の本多淳也さんに発表していただきました。タイトルは「多腕バンディット問題における漸近最適戦略について」です。本多さんは数理4研で修士をされ、その後数理1研で博士号を今年3月に所得し春から1研で助教をされています。本多さんは情報理論と統計的機械学習の両方に興味をお持ちで、今回の話は主に後者の問題を扱っています。

多腕バンディット問題というのは、ギャンブラーがスロットマシンで金儲けを狙う場合に、
$k$台ある性能の異なるマシンの中から$n$回($n>> k$)プレイする中で各マシンの報酬の期待値を推定し、それに応じて各台のプレイ回数を報酬ができるだけ大きくなるように決める問題です。仮定としては各マシンの報酬がある確率分布の集合に属すること、またその集合（例えばベルヌーイ分布、または正規分布であること）は既知であることです。応用例はギャンブルの他にいくつもあるそうです。

直感的には、期待値の一番高いマシンをできるだけ多くプレイし、その他のマシンはプレイしないことが理想となります。

 

 
ここで報酬を大きくするために、regretと呼ばれる’損失量’を小さくしたいのですが、最適戦略では各々のマシンを少なくとも$O(\log n)$回プレイすることが必要であることが知られています。この理論限界を達成する戦略をとれば、一番良いマシンを$n-O(k\log n)$回プレイできることになります。これを漸近的に達成するアルゴリズムを開発、解析することが多腕バンディット問題の目標となります。$\log n$の係数部はKL divergence を期待値制約のもとで最小化した量を用いてあらわされますが、これは分布$F$が期待値$\mu$以上の分布とどれくらい紛らわしくないかを表す指標で、これが小さい場合マシンは相対的にプレイされやすくなります。

発表ではまず先行研究の紹介と比較をしていただきました。既に提案されているUCB(Upper Confidence Bound)戦略と本多さんと竹村彰通先生が2010年に提案されたDMED(Deterministic Minimum Empirical Divergence)戦略、また近年再発見されたThompson samplingの3手法について計算量、性能（相対的に$O(\log n)$の項がどの程度小さくなるか）と解析の容易さの基準から比較しています。確率分布の変数（期待値、分散など）の数が$m$個の場合を$m$パラメータ問題と呼びます。確率分布の集合やコンパクト性、またパラメータの個数に応じても問題の難しさが変わるようです。

DMED戦略では、大まかには理論限界よりも少ない回数をプレイした台を探し、それがあればそれを次に選び、なければ現在の期待値の予想値が最大のものを選びます。UCBよりも計算が速い場合が多く、解析が容易（なため難しい確率分布モデルへも理論が拡張できる）ことが特長です。

確率分布を変えたときに漸近最適性の解析法や証明が異なり、特に個人的に興味深かったのは、一番自然に頭に浮かぶ正規分布の2パラメータ問題（期待値と分散が未知）については、複数パラメータで非コンパクトなこともあり解析が非常に難しいとされてきたことです。

regretを大きくしてしまう意味で問題となるのが、実際は期待値の大きいマシンが偶然期待値が小さいものと判断されてしまい、結果的になかなかそのマシンがプレイされないために期待値最大のマシンの推定を間違えている状態が続く状況です。発表後半ではこのような状況が起こる確率とそれから抜けるための待ち時間などを大偏差原理という理論を用いて評価することで解析を進め、この結果DMED戦略がコンパクト分布に対し漸近最適であること、更に非コンパクトな半有界サポートモデルでも同じ戦略で漸近最適であることを示されました。

実際の発表では、幸か不幸か導入部から山のように質問があり、スライド12枚目で2時間が経過する状況になりました。スライド後半では、なぜ正規分布の解析が難しいのかなども説明していただく予定だったようですが、結果的に数理的な内容は詳しい話をしていただく時間がありませんでした。質問の内容は多様でしたが特に、確率分布が時間定常でない場合はどうなるか、またDMEDとUCBの直感的違いなどについて議論になりました。

Asymptotically optimal policies in multiarmed bandit problems from j_honda

第２７回助教の会

こんにちは、第２７回数理助教の会ブログ担当の松島です。

先日の数理助教の会では、中務佑治さんに発表していただきました。

中務さんは今年の６月から数理第七研究室で助教に就任されている方で、主に数値線形代数を研究されています。今回は、最新の研究の成果である「Spectral divide-and-conquer algorithms for matrix eigenvalue problems and the SVD」について発表していただきました。

今回の発表で取り扱っている問題は対称正方行列における固有値分解(Symmetric eigenvalue decomposition)、そして一般の行列における特異値分解(SVD/Singular Value Decomposition) の問題です。

固有値分解では対称正方行列$A$が与えられた場合に、

$$ A = V^{\top} \Lambda V $$

となる直交行列$V$と対角行列$\Lambda$を求めます。

一方、特異値分解では$m$行$n$列の行列$A$が与えられた場合に、

$$ A = U^{\top} \Sigma V $$

となる$m$行$m$列の直交行列$U$、（各成分が負でない）$m$行$n$列の対角行列$\Sigma$、そして$n$行$n$列の直交行列$V$を求めるという問題です。

自分自身機械学習の研究をしていて、機械学習分野における固有値分解や特異値分解を用いたアルゴリズムというのもとてもなじみがありますが、固有値分解・特異値分解は一般に数理的なアプローチを行う非常に広範にわたる分野で使われている基本的な問題だと思います。なので固有値分解・特異値分解の新しいアルゴリズムというのは、それらがより具体的な問題解決のための基盤技術であるという意味でも非常に重要です。

特に最近では、多くの応用で、大きなサイズの行列を固有値分解・特異値分解するという場面が増えています。現代の計算機にはメモリの階層構造があり、高速にアクセス可能なメモリほど容量が小さくなっていくのですが、一度に高階のメモリに乗り切らないほど大きなサイズの行列を扱う場合、高階のメモリに行列の要素を移すコミュニケーションコストの方が、数値演算のコストよりも高くなってしまう事になります。そのため、大きい行列を扱う場合は、コミュニケーションコストを最小化しつつ、数値演算のコストも低く押さえることが望まれます。中務さんが発表されたのはそのような要望に応える新しいアルゴリズムの提案です。

今まで大きなサイズの行列に対して、コミュニケーションコストを最小におさえられる操作としては、行列積、コレスキー分解、QR分解が知られています。中務さんのアイディアは、これらの操作をうまく利用してやる事でコミュニケーションコストを最小に押さえられる固有値分解・特異値分解のアルゴリズムを与えようという事です。アルゴリズムの流れを固有値分解の場合で説明しますと、まず$A$の正の固有値と負の固有値がだいたい同じくらいになるように$A-\sigma I $とした後にQR分解を用いた極分解で$A-\sigma I =UH $における$U$を求めます(1.)。次に$(I+U)/2$の列空間とその直交補空間の正規直交基底に対応する行列$V_+$、$V_-$を求めます(2.)。これらを用いて$A$をブロック対角化することができるので(3.)、以上の操作をブロックに再帰的に施す事によって対角化を得る(4.)、という流れになっています。行列$A$の極分解とは

$$A=UH$$

となる直交行列$U$と半正定値行列$H$を求める事で、特に$A$が対称正方行列の場合は、$A$の全ての固有値を、符号に応じて$1$と$-1$に変えてできる行列$\mathrm{sign} (A)$に対応します。

アルゴリズムの内容もさることながら、このような数学的に定義された分解を正確に求めるためのアルゴリズムに関してはその評価も非常に重要になってきます。今、計算機上で表現されている行列$A$に対して固有値分解のアルゴリズムが直交行列$\hat{V}$と対角行列$\hat{\Lambda}$を出力したとします。このとき、$A$と$\hat{V}{}^{\top}\hat{\Lambda}\hat{V}$は非常に値が近い行列ではありますが、完全に一致するということは一般には期待できません。そのため、対象としているアルゴリズムが出力する$\hat{V}$および$\hat{\Lambda}$はどのような意味で$A$の真の固有値分解に近いのか、という問いに答えなければなりません。中務さんは、提案されているアルゴリズムが入力$A$に対して出力した$\hat{V}$および$\hat{\Lambda}$に関して、

$$A-\hat{V}{}^{\top}\hat{\Lambda}\hat{V} = e \|A\|_2$$

がいつでも成り立つ事を数学的に証明しています。ここで$e$はオーダーが$2^{-53}$程度の行列で、この値は標準的な浮動小数点計算で生じる誤差の大きさです。これはつまり、アルゴリズムが出力した$V$および$\hat{\Lambda}$が真の分解となる行列$\hat{A}=\hat{V}{}^{\top}\hat{\Lambda}\hat{V}$と、与えられた行列$A$が近いという事を示しており、特に後方安定性（Backward Stability）といわれています。このような後方安定性の議論は数値線形代数の分野においては多く見られる重要な議論のようですが、少し分野が違う自分にとっては目新しいもので、とても興味深かったです。

他にも、極分解を効率よく求めるためのポイントや特異値分解の場合への適用法、さらに極分解を効率よく求める方法の話など、興味深い内容がありましたが、ここでは割愛させていただきます。詳しくはスライドのほうをご覧ください。

Spectral divide-and-conquer algorithms for eigenvalue problems and the SVD from yuji_nakatsukasa

個人的には、先ほどの後方安定性の話の他にも数値線形代数の分野では良く知られているが自分には目新しい話というのがいくつも出てきて、このような話を共有させていただきとても勉強になりました。一方で、メモリ容量の問題で大規模な計算を行う場合計算アルゴリズム自体を再考しなければいけない、という現象はまさに自分が機械学習の分野で取り扱っていることと同じことでした。大規模データを扱うことのむずかしさと面白さというのは、「ビッグデータ」という言葉とともに多くの実社会に関する面で強調されている昨今ですが、数理工学の基礎的分野においても広く議論の的になっているとても普遍的な問題なのだということを感じました。

数理六研　松島

[参考文献]

Yuji Nakatsukasa and Nicholas J. Higham, Stable and efficient spectral divide and conquer algorithms for the symmetric eigenvalue decomposition and the SVD, SIAM Journal on Scientific Computing, Vol. 35(3), pp. A1325-A1349, 2013. pdf file. Codes available at MATLAB Central File Exchange

第２６回助教の会

第２６回助教の会では数理情報学第６研究室の松島慎さんに発表していただきました．松島さんは昨年度博士号を取得し，現在ポスドクとして数理６研に所属しています．今回の助教の会での発表のタイトルは” Scaling up Machine Learning algorithms for classification”でして， ざっくり言うと，分類問題に対する機械学習アルゴリズムを高速化したことについて話していただけました．これは博士論文にまとめられた成果（の一部）ということです．

機械学習と言っても数理的にはある種の最適化問題を解くことで，本研究では二値分類問題に対しサポートベクターマシンで最適化問題を定式化します．要するにやりたいのは与えられたサンプル点に対し下図（１つ目）のような直線を引くことで，その直線は下図（２つ目）にある最適化問題を解くことで得られます．

この分野では，この最適化問題を近似でよいので速く解くことが求められるのですが，そのための有力なアプローチが，双対問題に対し座標降下法（Coordinate descent method）を用いるというものです．ただし，原型となるアルゴリズムでは，データが巨大化するにつれ，同時にアクセスすべきデータも増えるため，高速化のためにはデータを分割してキャッシュメモリに移し，データを局所的に参照して計算できるようにすることが不可欠になります．今回の松島さんの発表では，その目的に適うように既存のアルゴリズムを改良して実装し，実際に高速化できていることを示してくれました．

研究の核心である実装部分の話についてもう少し詳しく述べると，巨大データを扱う際には2010年にBlock Minimization という技術が提案されており，これによりハードディスク上の巨大データを小分けにしてメモリに読み込ませることで，メモリ上で高速に処理できるようになったのですが，この技術ではメモリ上で処理している間ハードディスクからの読み込みが停止するため，その時間が無駄になっていました．これに対し松島さんの方法では，できるだけ読み込みと計算処理の両方が常に動き続けるように改良されています．視覚的に分かりやすく説明しようとすると，下図のReaderとTrainerが両方同時に動くようにしたということです．

この説明については，下のYoutubeの動画が非常によくできていますので是非ご覧ください（最下部に添付の発表スライドにも埋め込まれています）．

　　　　　

話は単純ですが，これを効率的に行うには，もとの座標降下法の長所である（１）反復で更新する変数は一つでよい，（２）一次収束性，（３）反復中で不必要なサンプル点を除去できる，といった性質を損なわないようにする必要があります．松島さんはそのメリットを維持した上でスケーラブルな実装を与え，実際に既存手法より数倍高速な実験例を示し，それだけでなく，既存手法では収束していなかったものでも本手法で現実的な時間で収束する例も紹介してくれました．また，松島さんの枠組みは，サポートベクターマシンに限らず広く応用できるので，最近は他の問題にも焦点を当てて研究を行っているそうです．

今回の話は数学的に難しい話は少なく，サポートベクターマシンや座標降下法などは学部生にも理解できる話なのでもっと周知する機会があってもいいような気もしました．発表を思い出すと、松島さんのやわらかい物腰のおかげか，いつも以上に質問やコメントが多かったように思います．近年，自分の研究する数値計算の分野でも，キャッシュメモリをうまく活用するアルゴリズム設計が重要視されているので興味深く感じました．そして研究の歴史に着目してみると，双対問題に座標降下法を用いることが提唱されたのが2008年で，Block Minimization が2010年です．日進月歩でこの手の研究が精力的に行われていることが伺えますね．この分野の中での松島さん自身の今後のご活躍を楽しみに思いました．

数理3研助教 相島健助

　　　　　

Scaling up Machine Learning Algorithms for Classification from smatsus

第２５回助教の会

こんにちは，数理6研助教の冨岡です．今回は数理3研の助教の相島健助さんに「固有値問題に対する射影法について」というタイトルで行列の固有値を求めるための Lanczos 法と，Jacobi-Davidson 法という2つの方法のレビューをして頂きました．相島さん自身の研究の話もあったのですが，未発表ということで，ここでは割愛させて頂きます．

量子力学における物理量の計算，Google のページランクや，システムの安定性解析など，多くの実世界の複雑な問題では，非常に大きな行列の固有値および固有ベクトルを求めることが必要になることがあります．ここで，固有値を計算したい$n\times n$行列を$A$（ここでは簡単のために対称と仮定します）として

$$ A x = \lambda x $$

を満たすような $\lambda$ と $x$ をそれぞれ，$A$の固有値，固有ベクトルと呼びます．固有値や固有ベクトルを求めるには，行列$A$が小さい場合にはQR法などの直接法が用いられますが，行列が大きい場合には反復法が有効であることが知られています．

反復法としてもっとも代表的な方法はべき乗法です．べき乗法では適当な $x_0$ を初期ベクトルとして，

$$ x_{n+1} = A x_{n} $$

のように反復計算をします．このような反復を十分長く続けると，ベクトル$x_n$が$A$の絶対値最大の固有値に対応する固有ベクトルに収束することが知られています．この方法は非常に単純であり行列Aとの掛け算さえ計算できればよいので広く活用されていますが，収束に要する反復回数が大きくなってしまうことが問題です．

べき乗法は1本のベクトルを更新して行くのに対して，射影法は$m (>1)$次元の部分空間を更新して行く方法です．具体的には

1. 何らかの方法で部分空間の正規直交基底 $V_m$ を生成
2. 小さな行列 $V_m^TAV_m$ の固有値$\theta$，固有ベクトル$y$を求める．
3. $V_my$ がもとの行列$A$の固有ベクトルになっていれば終了，そうでなければ何らかの方法で部分空間$V_m$を更新

というような手続きを行います．ここでのポイントは，$m$が小さければ小さいほど，2で解くべき固有値問題は簡単になる代わりに，より多くの反復が必要になるというトレードオフがあるという点です．逆に言えば，$m$を増やすことで，べき乗法（$m=1$）に比べれば1回あたりの計算は大変になるものの，少ない回数の反復で固有値，固有ベクトルが得られるということです．具体的にどのように部分空間を作り，どのように更新するのかはアルゴリズムによって異なります．

射影法のひとつであるLanczos法は部分空間$V_m$として，適当な$x$を初期ベクトルとする Krylov部分空間 ${\rm span}(x, Ax, A^2x,\ldots,A^{m-1}x)$を用います．さらに，上記3で収束しなかった場合，$x:=V_my$として再度Krylov部分空間を構築します．

射影法のもう1つの手法 Jacobi-Davidson法は適当な正規化された初期ベクトルを$V_1$として，各反復ごとに1本の基底を追加しつつ$V_m$を更新します．各反復で追加する基底は修正方程式の解として得られます．修正方程式は一見複雑な形をしていますが，実はその解はRayleigh商反復という，また別の反復法の更新式と等価であることを示すことができます．また，$m$が大きくなると2で解くべき固有値問題が重くなるため，あらかじめ定めた$m$まで基底を増やすごとに，初期ベクトル$V_1=V_my$として$m=1$にリセットすることが，実用上有効であるとのことでした．

ふだん固有値や特異値を計算することは多いのですが，実際に内部でどのようなアルゴリズムが動いているのかはあまり考えずに使っていることが多いので，ユーザの立場としても様々なトレードオフを理解した上で使うことが重要であると感じました．個人的には最近注目されているランダム射影に基づく方法[1]と上記の手法の関係が気になります．

[1] N. Halko, P. G. Martinsson, and J. A. Tropp (2011) Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Reviews, 53(2): 217-288.

なお，今回の発表の内容は来週，道後温泉で開かれる数値解析シンポジウム(NAS2013) でより詳しく発表するということです．

Jokyonokai130531 from nwpmq516