第２５回助教の会

こんにちは，数理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

第２４回助教の会

2013年度2回目となる第24回助教の会では，数理第六研究室助教の冨岡亮太さんに「低ランク行列補完のためのマトロイド理論」というタイトルで発表して頂きました．

全ての場合の事象を観測することができず，部分的にしか情報が得られない場合は世の中に多数存在します．例えば，映画の評判に関するアンケートを例に取ると，各人が見た映画は異なっており，全てのユーザーに対して全ての映画の評価を得ることはできません． このような全ての情報を得ることができない場合に，手持ちの情報からどれだけの情報を原理的に復元できるかという問題は応用上も重要となります．

今回の講演における問題設定は，与えられたマスク（観測できる要素が特定される情報）の下で，観測できていない要素をどれだけ理論的に復元できるかというものでした．マスクは観測できる要素を１，観測できる要素を０とする行列表現以外にも，この行列の行と列に対応する頂点集合のうち，見えている要素に対応する枝を与えた２部グラフ表現も考えられます． この２部グラフの組み合わせ的な性質が解析において重要ということでした．


上記の問題は「部分的な観測から低ランク行列を復元できるか」という問題と数学的に表現でき，この必要条件と十分条件が今回の研究により求められました．必要条件は，行列のサイズを m * n，行列のランクを r とすると
(1) 観測された要素の数が r(m+n-r) 個以上であること，
(2) 各行，各列において r 個以上の観測があること，
(3) マスクの２部グラフから任意の r-1 本の枝を取り除いても連結であること，
でした．また，ランク r 以下の行列では任意の (r+1) * (r+1) 小行列式が零になるという性質を用いると，(r+1) * (r+1)の部分マスクで枝が一本のみ欠けている場合にその欠けている要素を復元でき，この復元された要素に対応する枝をマスクに加えます．今回の研究では上記の操作を繰り返してマスクが完全グラフになる場合を r-closable と呼んでおり，これが一意穴埋め可能の十分条件ということでした．

更に，今回の研究では matroid (マトロイド)という概念を用いて，部分的な復元という視点を導入しています．従来の穴埋めに関する研究では行列の穴埋めが成功するか失敗するかという観点にのみ着目されており，この点で先行研究から大きな進展が存在します．Matroid は何らかの意味で独立な部分集合の全体と定義されており，例えば，行列の列ベクトル集合のうち線形独立な列ベクトルの全体は vector matroid (ベクトルマトロイド)と呼ばれています．今回の講演では， determinantal matroid という行列の要素全体の集合のうちランク r で独立に選べる要素の集合を用いた理論が展開されました．この determinantal matroid の独立性に関する解析から，観測できない要素のうち復元できる要素とできない要素は，観測要素に従属か独立かという点で説明できることが今回の研究によりわかりました．この独立性 or 従属性はある非線形写像のヤコビアンから得ることができます．今回の結果には様々な応用例を考えることができ，最初の映画の例だと，追加の調査を行う際に新たに評価すべき映画を指定できるということが考えられます．

観測できる要素の箇所という情報だけから原理的に復元可能な要素とそうでない要素の区別を理論的に与えた今回の研究結果からは数理研究の持つ可能性の奥深さを強く感じました．


数理６研　森野佳生

第２２回助教の会

本年度の最終回，第２２回助教の会では数理情報学第２研究室助教の小林佑輔さんに発表していただきました．発表タイトルは「制約付き２マッチング問題に対するアルゴリズム」ということで，最適化と離散構造に関するお話でした．

組合せ最適化の分野における代表的・古典的問題のひとつに「最大マッチング問題」があります．与えられたグラフにおいて頂点の重複がないように枝の本数を最大化する問題です．この問題の解ではひとつの頂点につながる枝は高々１本ですが，２本まで許すことにして枝数の最大化を考えるのが「最大2-マッチング問題」です．この問題の解はサイクルとパスの集まりになりますが，最大マッチング問題に帰着されることが知られているので計算は多項式時間でできます．さらに，「長さ$k$以下のサイクルは使わない」という制約を入れたものを「制約付き最大2-マッチング問題」と呼びます．禁止サイクル長$k$の値が変わることによって問題の難しさが大きく変わってくることが知られています．

最大マッチング問題では，各枝に重みを与えることにより重みの和の最大化を考えることもできます．制約付き最大2-マッチング問題で各枝に重みを与えたものを「制約付き最大2-マッチング問題（重み付き）」と呼ぶことにします．重み無しの場合と同様に，制約付き最大2-マッチング問題（重み付き）においても禁止サイクル長$k$を変えることで問題の難しさが変わることが知られています（図１）．

図１　先行研究

 
$k=3$である場合については，計算時間評価の境目に当たる問題になっており，一般に多項式時間で解けるかどうかは分かっていません．グラフの次数を3以下に限った場合にはHartvigsen-Li (2007)により問題が多項式時間で解かれることが示されていますが，提案されているアルゴリズムは単純ではありません．今回の小林さんの研究では，グラフの次数を3以下に限った場合を考え，問題の構造をうまく抽出して利用することにより単純な局所探索アルゴリズムを得ています（図２）．得られたアルゴリズムの正当性は，離散凸解析とグラフに関する考察から導かれます．

図２　研究成果

もう少し詳しく見てみます（図３）．まず，枝集合の次数列（枝集合が与えられ場合の各頂点の次数を並べたもの）に注目します．制約付き最大2-マッチング問題（重み付き）を直接解くのではなく，次数列ごとに枝の重みの和を最大化することを考えます．ここで重要なのは，

・次数列から対応する重みの最大和を求める関数は次数列に関してM凹関数である
・グラフの次数列が指定されれば重みの最大和は多項式時間で計算可能である

という２点です．後者により，局所探索が実行可能であることが分かります．次数列が与えられた場合の重みの最大和の計算可能性はグラフに関する考察から得られますが，その際にグラフの次数が3以下であることが効いてきます．さらに前者と離散凸解析の結果から，次数列に関する局所探索により制約付き最大2-マッチング問題（重み付き）が多項式時間で解けることが保証されます．

図３　提案アルゴリズム

今回の小林さんの発表は，「個別の問題における計算可能性・計算時間」というミクロな視点と「効率的な計算を可能にする離散構造」というマクロな視点とを組み合わせた研究についてでした．次数列を経由することで，考えていた問題の離散構造が分かりやすい形で現れ，局所探索という単純なアルゴリズムで問題が解けることが示されました．ミクロとマクロをつなぐ中間的な手段を経由することの「うまさ」，妙味が感じられる内容で，とても興味深いものでした．


数理５研助教 廣瀬善大

Jokyo20130218 from y-kobayashi

第２１回助教の会

第21回目となる今回の数理助教の会は、国立情報学研究所　ERATO特任研究員の前原貴憲 (まえはらたかのり)さんに「行列の同時ブロック対角化問題」について発表していただきました。
タイトルにある「行列の同時ブロック対角化問題」とは n次の実正方行列がN個与えられた下で、共通の直交行列Pをうまくとってできるだけ細かいブロック対角行列に変形するという問題です。

この問題は、古くは、結晶物理に出てくるハミルトニアンの一般的な形を考える際にWignerが考えていました。その後も、物理では群論的手法として抽象的な議論がされてきました。しかし、このような問題は、工学における半正定値計画や、信号処理における独立成分分析、最近では、動的ネットワークの安定性解析など、いたるところに現れてきます。後者では、規模の大きい行列を数値的に扱うことを想定した　うまい同時ブロック対角化が議論されています。

理論的には one-by-one method として二つの可換な対称行列の同時対角化の証明が100年以上も前から知られていました。しかし、このような素朴な方法は、巨大な行列に対して数値的に対角化を行う場合には、誤差に弱いことが知られていました。 Bunse-Gerstner, Byers, Mehrmann (1990)は最適化の枠組みで、このような問題を取り扱っています。（Jacobi-like method）その後、JADEとよばれる方法によって、複数の行列の同時ブロック対角化が提案されましたが、厳密にうまくいくことが証明されたわけではありません。うまくいかない例もつくられています。


一方で、前原さんの研究では、このような問題に対して、数学的にしっかりとした理論を用いて、アルゴリズムを提案しています。前原さんのアイディアの数学的な基礎は1905年のSchurの研究にまでさかのぼります。（いわゆるSchur の補題として知られる結果。）　行列の組によって生成される＊代数 (star algebra) の交換子代数を考え、その元をランダムに一つとってきて、対角化する行列を見つけるというものです。行列の組がもつ対称性に注目し、行列自身を群の（有限次元）表現とみなすことで、既約分解の一般論を適用します。このような発想によって、one-by-one method の拡張を得たのが前原さんの一連の研究結果です。

スライドに出てくるArtin-Wedderburnの構造定理は、かなり抽象的な代数の話題ですが、これらの結果をうまく現実的な問題に適用できることが素晴らしいと思いました。ランダム行列理論を用いた誤差評価などは、数値的に統計的推定を扱う際にも役立ちそうで興味深かったです。 Jacobson根基が残るケース（半単純代数でないケース）の分解については、まだできていないとの話でしたが、工学的な応用のみならず、統計物理でよく出てくるHeisenberg群などの例もあるので、今後に期待したい所です。

また、この助教の会のブログをいつもご覧になっている皆さんの中には、坂下達哉さんの発表(第12回 2012年1月)の際の密度行列のテンソル積の既約分解の話 (同スライド p.23)が、前原さんの取り扱っている問題の特殊ケースになっていることに気付いた方もおられるかもしれません。発表の際にも指摘したのですが、坂下さんの研究では、並列計算を駆使して高精度の固有値計算を行っており、今回の前原さんの研究とも大いに関連があるように思いました。

simultaneous block diagonalization of matrices from tmaehara

数理４研　田中 冬彦

第２０回助教の会

今回のブログを担当する東大情報基盤センターの佐藤一誠です。

今回は、東京大学情報理工学系研究科　助教の田中冬彦さんに発表していただきました。

田中さんは、統計理論の研究だけでなく、理論物理・実験物理への統計理論の啓蒙や研究交流に日夜励んでいらっしゃいます。

今回の発表では、情報幾何のベイズ統計における応用例について発表して頂きました。

過去のデータから未来のデータについて予測を行う場合、ベイズ統計ではベイズ予測分布を用います。

最尤推定が１つのモデルパラメータを用いて予測を行うのに対し、

ベイズ予測分布では、モデルパラメータの事後分布を求め、事後分布で未来のデータの分布を平均化することで予測を行います。

さて、このベイズ予測分布はどのような意味で最適なのでしょうか？
ベイズ予測分布は、平均リスクを最小にするという意味で最適であることが、1975年にAitchisonによって示されています。

ただし、この最適性は事前分布の取り方に依存するという問題点があります。

つまり、どのような事前分布が望ましいかという疑問が新たに生じます。

特にパラメータに関する事前情報がない場合に使う事前分布（無情報事前分布）を
統計モデルのみから決めるという研究が海外を中心に盛んに行われています。
ただし、無情報事前分布というのは慣用的な用語で、デフォルトで使う事前分布といった意味だそうです。
これまで様々な文脈で無情報事前分布が提案されてきました。
今回注目する事前分布は優調和事前分布です。
よく知られている汎用的な無情報事前分布にJeffreys事前分布がありますが，
もし優調和事前分布が存在するならば，Jeffreys事前分布より予測の意味でも最適であることが知られています。

ここで重要なのは

・優調和事前分布が存在するか

・存在するならばどのようにして構成できるのか

です。

これらを理論的に示す道具として、統計モデルの情報幾何学が活躍します。

統計研究者の中でもあまり理解されてない情報幾何の利点として以下が挙げられます

１．記述の簡潔さ（“縮約“やEinstein規約）

２．パラメータの取り方に依存しない結果の抽出

３．望ましいパラメータの取り方を見つける

４．微分幾何学的な道具立ての活用

これらは、今回の発表の後半で、時系列モデルを扱う場合に確認することができます。

Seminar 130107 ver2_public_embed from tanafuyu

第１９回助教の会

第１９回助教の会は情報基盤センターで助教をされています佐藤一誠さんに話をしていただきました．タイトルは「“基礎”からのBayesian Nonparametrics-点過程と機械学習の数理-」ということで，ランダム測度からはじまり機械学習で広く使われている様々な確率モデルとの関係を概観していただきました．今回の話のキモは「フビニの定理」です．今日，自然言語処理や機械学習の分野では「ノンパラメトリック」なベイズモデルが広く使われています．ここでいうノンパラメトリックとは特定のパラメトリックモデルを仮定しない広いクラス（無限次元）のモデルです．ノンパラメトリックなモデルを考えると数学的に難しい部分が出てきますが，フビニの定理を通して眺めるとすっきりするという点は重要であったと思います．

まずはその雰囲気を概観してみましょう．普通のパラメトリックモデルでのベイズモデリングでは次の「ベイズの定理」から話が始まります：
$$p(\theta|x)=\frac{p(x|\theta)\pi(\theta)}{p(x)}.$$
ここで，$\pi(\theta)$はパラメータの事前的な確からしさを表わす事前分布で，$p(x|\theta)$はパラメータ$\theta$のもとでのデータ$x$に対する当てはまりの良さを表わす尤度です．こうして得られた$p(\theta|x)$を事後分布と呼びます．しかしノンパラメトリックモデルの場合，上のように密度関数で割ったりするという操作は必ずしも自明ではありません．そこで，以下のような別の表記を使ってみましょう：
$$\int \int h(x,\theta) p(\theta|x) d\theta p(x) dx = \int \int h(x,\theta) p(x|\theta) dx \pi(\theta) d\theta,$$
ただし$h$は任意の非負関数．ここで，積分の順序が交換されていることに注意してください（内側が$x$に関する積分か$\theta$に関する積分か）．これを「フビニの定理」と呼びます．フビニの定理は無限次元の世界でも（ある条件のもと）成り立ち，確率密度関数で割る操作を陽に行わないで，事後分布を導く一つの見方を与えます．

さて，話の本筋に入りたいと思います．まずはCompletely Random Measure (CRM) から話は始まります．CRMはランダムな非負測度であり，互いに素な集合上の測度は独立になるようなものです．例えば全世界でおきる交通事故の件数の分布を考えると分かりやすいでしょう．ある地域とそれ以外の地域の交通事故の発生件数は独立であるとすると，これはCRMになります．このCRMは自動的に無限分解可能分布の構造を持ち，そのためLevy Processとして表わすことができます．さらにCRMの非負性からガウス成分はなくポアソン成分のみが残ります（Levy-Ito Decomposition）．ちょっと話が難しくなりましたが，要約しますとCRMはポアソン分布とその飛躍の大きさに関する分布で表わせます．先ほどの交通事故の例ですと，交通事故の頻度がポアソン分布に従い，その時の損害金額が飛躍の大きさと考えることができます．ここで，各地域における頻度と飛躍の大きさを表現する関数をLevy measureと呼び，これをいろいろとモデリングすれば様々なCRMを導くことができます（下図）．

例えばLevy measureとしてガンマ分布を用いると (飛躍の大きさをガンマ分布とし，頻度はbase measureで与える)，対応するCRMはガンマ過程と呼ばれます．ベイズモデリングを考えると，ガンマ過程を事前分布とした場合の事後分布を求めたくなります．たとえば交通事故による損害金額の分布を推定する際に，事前分布にガンマ過程を用いて，実際観測されたデータから事後分布を構成することを考えます．その際，役に立つのが先ほど述べましたフビニの定理です．フビニの定理を用いて次々と積分の順序交換をしてゆくことにより事後分布が自然に得られるます．さらにその操作を通じて，Chinese Restaurant Processと呼ばれるサンプリング手法も自然に導かれるとのことです．

ガンマ過程は確率測度を与えませんが（積分して１にならない），正規化することにより確率測度が得られます．ガンマ過程を正規化したものは有名なDirichlet過程になります．Dirichlet過程に関してもフビニの定理をうまく使うことにより，事後分布が自然に求まります．最後に今回話していただいた，CRMからポアソン過程→ガンマ過程→Dirichlet過程という一連の繋がりを包括した概観図を載せておきます（下図）．

ベイズジアンノンパラメトリクスの分野は様々な応用があるだけでなく，昔の確率論から流れる堅い数学的なベースもあり，奥深い分野だと感じました．

数理第五研　鈴木大慈

第１８回助教の会

第１８回目の助教の会では，数理第５研究室の鈴木大慈さんに発表していただきました．発表のタイトルは”PAC-Bayesian Bound for Gaussian Process Regression and Multiple Kernel Additive Model”というもので，無理やり和訳すると「ガウス過程回帰を用いたマルチプルカーネル加法モデルに対するPAC-Bayesian 上界」とでも言うのでしょうか．本発表はConference on Learning Theory (COLT2012) という学習理論の国際会議での発表を，一般の数理の人向けに少しアレンジしたものでした．

発表タイトルの意味するところを読み解いていくと，研究成果としては，ノンパラメトリック回帰のためのスパース加法モデルを推定するために用いられるマルチプルカーネル学習に対し，PAC-Bayes 的な手法を用いて収束のレートを評価したというものです．その際，事前分布にガウス過程を導入し，これが理論評価のための１つのキーになっているそうです．

少し用語の説明も加味しながら研究成果をざっくりと解説すると，まずノンパラメトリック回帰とは，学部の授業で習う線形回帰Y=Xθのようなパラメータθを含む形をとらずに，（ある程度なめらかな関数を用いて）モデルを推定する方法です．やりたいのは，下図のように，サンプル点が与えられたときに適切な曲線をえがくことです．

スパース加法モデルは高次元ノンパラメトリック回帰において各説明変数ごとに（非線形）関数を割り当て，それらの和をデータに当てはめるモデルです．その際，各非線形関数を再生核ヒルベルト空間というある関数空間から適当にとってきた関数で表現する，いわゆるマルチプルカーネルが有用であるということでした．具体的にモデルを推定する方法としては，ある最適化問題を解くことになりますが，その際，正則化項を導入して少ないカーネル関数の和でモデルを表現することを目指します．これがマルチプルカーネル学習です（下図）．

この推定において，サンプル数nを増やすにつれて推定された関数は求めるべき関数に収束します．その収束性に関しては，Restricted Eigenvalue 条件を仮定することで高速なレートを達成できることが理論保証されていましたが，Restricted Eigenvalue 条件は強い仮定で，これを外すことが望まれていました．

鈴木さんはこの問題に対し，上に示した意味でのスパースな解による推定と，事前分布に（スケーリングした）ガウス過程のミクスチャーを導入した，ベイズ的マルチプルカーネル学習を考え，このタイプのマルチプルカーネル学習を行えば，既存の収束レート評価に必要だったRestricted Eigenvalue 条件を外した形で高速な収束レートを実現できることを示しました．発表ではいくつかの問題例に対する収束レートを示し，スケールミクスチャーの導入により高速なレートの評価ができることを強調してらっしゃいました．

鈴木さんの助教の会での発表は２回目でして，今回は数式や関数解析に登場する用語がかなり多くて専門的な内容になっているなあと思って１回目の鈴木さんの発表（第１回助教の会）のブログを見てみると，前回もかなり内容が盛りだくさんだったことを思い出しました．鈴木さんの手を抜かない性格と研究の迫力を感じさせられました．

数理3研助教 相島健助

PAC-Bayesian Bound for Gaussian Process Regression and Multiple Kernel Additive Model from Taiji Suzuki