2013年1月21日月曜日

第20回助教の会



今回のブログを担当する東大情報基盤センターの佐藤一誠です。
今回は、東京大学情報理工学系研究科 助教の田中冬彦さんに発表していただきました。

田中さんは、統計理論の研究だけでなく、理論物理・実験物理への統計理論の啓蒙や研究交流に日夜励んでいらっしゃいます。
今回の発表では、情報幾何のベイズ統計における応用例について発表して頂きました。

過去のデータから未来のデータについて予測を行う場合、ベイズ統計ではベイズ予測分布を用います。
最尤推定が1つのモデルパラメータを用いて予測を行うのに対し、
ベイズ予測分布では、モデルパラメータの事後分布を求め、事後分布で未来のデータの分布を平均化することで予測を行います。

さて、このベイズ予測分布はどのような意味で最適なのでしょうか?
ベイズ予測分布は、平均リスクを最小にするという意味で最適であることが、1975年にAitchisonによって示されています。
ただし、この最適性は事前分布の取り方に依存するという問題点があります。
つまり、どのような事前分布が望ましいかという疑問が新たに生じます。

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

ここで重要なのは

・優調和事前分布が存在するか
・存在するならばどのようにして構成できるのか

です。

これらを理論的に示す道具として、統計モデルの情報幾何学が活躍します。
統計研究者の中でもあまり理解されてない情報幾何の利点として以下が挙げられます

1.記述の簡潔さ(“縮約“やEinstein規約)
2.パラメータの取り方に依存しない結果の抽出
3.望ましいパラメータの取り方を見つける
4.微分幾何学的な道具立ての活用


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


2013年1月2日水曜日

第19回助教の会


第19回助教の会は情報基盤センターで助教をされています佐藤一誠さんに話をしていただきました.タイトルは「“基礎”からの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と呼ばれるサンプリング手法も自然に導かれるとのことです.

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

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

数理第五研 鈴木大慈

2012年11月21日水曜日

第18回助教の会


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

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

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

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

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

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

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

数理3研助教 相島健助

2012年8月24日金曜日

第17回助教の会

今年度5回目となる今回はお隣りのシステム情報学専攻から亀岡弘和先生をゲストとしてお招きして「スパース表現による音響信号処理」というタイトルでお話しいただきました。亀岡先生は2007年に博士の学位を取られてからNTTコミュニケーション科学基礎研究所に勤めていらっしゃいますが、2011年からは客員准教授として6号館に戻っていらっしゃいます。

今回は複素NMFと、複合自己回帰系という2つの手法に基づいた音源分離の話を紹介して頂きました。音源分離はその名前の通り、観測された音声・音響信号を少ない数の音源に分離する技術で、音楽の録音を楽器の音色ごとに分離したり、会議の録音を話者ごとに分離したりするなど様々な応用があります。

NMF (Nonnegative Matrix Factorization; 非負行列分解)は音源分離に頻繁に用いられる手法で、数学的には与えられた非負の要素からなる行列(非負行列)を2つの非負行列の積に分解するというふうに定式化されます。そもそもどうして非負の要素からなる行列を考えるかと言うと、音響信号のフーリエ変換の振幅は負の値を取らないからです。実際の音声信号の振幅スペクトログラムを図1右上に示します。振幅スペクトログラムは縦の長さがフーリエ変換の周波数の数、横の長さが音響信号の長さ(時間)の非負行列です。この図からわかるように、現実の音声信号は異なる周波数の音がばらばらに鳴っているのではなく、かなり規則的にそろって発生しています。このことから、「どの周波数の信号をどれだけ含んでいるか」という情報と「どのタイミングで音が鳴ったか」という情報のペアをひとつの音源とみなし、現実の音声信号を複数の音源の重ね合わせとして表現しようというのが、NMFを用いる動機です。多くの既存研究では振幅スペクトログラムに対してNMFを適用しています(図1右上)。NMFでは前者の「それぞれの音源はどの周波数の信号をどれだけ含んでいるか」という情報を行列Hで、後者の「それぞれの音源はどのタイミングでどの強さで鳴るか」という情報を行列Uで表現します(図1右下)。

図1 振幅スペクトグラムの分解表現

既存のNMFは振幅スペクトログラムを対象としているため、観測信号が複数の音源が線形に混ざったものであったとしても、振幅スペクトルの上では線形な混合にならないという問題点がありました。亀岡先生はこの点に注目し、混合の線形性が成り立つ複素スペクトルの上でNMFを定式化することに成功しました[1]。振幅スペクトルの代わりに複素スペクトルを対象にするためには信号の位相を考える必要があります。複素NMFは上述のペアがすべての時間点、周波数において位相の任意性を持っていると考え、これを含めて最適な信号の分解を計算するアルゴリズムになっています。複素NMFでは位相の情報も最適化の中で陽に扱っているので、従来の振幅スペクトログラムに基づくNMFの場合のように後処理で位相を求める必要がありません。

このアルゴリズムを使うと分解した1つ1つの音声信号の成分に対して、音量をコントロールしたりピッチを変更したりすることができます。これらは亀岡先生が作成された MusicFactorizer というツールのデモページで聞くことができます。

後半のお話は、複数の人間が同時に話しているような状況での(複数のマイクの)録音から個々の話者の音声を分離する、いわゆるカクテルパーティー問題が対象です。

ここでの音源とは人間の発する音声ですので、どのような周波数でも好きに組み合わせられるわけではありません。このような音声をモデル化するための有名な手法のひとつとして「自己回帰系」があります。P次の自己回帰系は、白色雑音を入力として、過去P次点前までの出力に依存して現在の出力が決まるような、時系列モデルです。従来の自己回帰系では、入力が白色雑音に限定されているため、P個の自己回帰フィルタ係数で音源の周波数特性を表現する必要がありました。一方、亀岡先生の提案している「複合自己回帰系」モデルは複数の白色とは限らないパワースペクトル密度と、複数の自己回帰フィルタを用意し、その組み合わせで音源を表現します。つまり、あらかじめ多くのパワースペクトル密度と多くの自己回帰フィルタを用意しておけばその組み合わせで各フレームの音源を簡単に表現することができるというわけです。その上、従来はフレームと呼ばれる信号の単位ごとに自己回帰フィルタを計算していたのですが、「複合自己回帰系」では用意しておくパワースペクトル密度や自己回帰フィルタを信号全体で最適化し、フレームをまたいで再利用することができるので無駄がない、というメリットもあります。

図2: 複合自己回帰系

このような信号源のモデルに残響のモデルを組み合わせることにより、現実的な環境で従来の方法よりうまく話者を分離できることが分かりました[2]。音声信号処理の伝統的な方法である自己回帰モデルと、少ない数の基底を組みあせて使うという行列分解モデルの考え方がうまく組み合わさって高い性能を出せるという点が非常に興味深いと感じました。

数理6研助教 冨岡亮太


[1] 亀岡弘和, 小野順貴, 柏野邦夫, 嵯峨山茂樹 (2008) 複素NMF: 新しいスパース信号分解表現と基底系学習アルゴリズム. 日本音響学会論文集.
[2] 濱村真理子, 亀岡弘和, 吉岡拓也, ルルージョナトン, 柏野邦夫 (2010) 複合自己回帰系による音声パワースペクトル密度モデルを用いたブラインド音源分離と残響除去. 日本音響学会講演論文集.

2012年7月12日木曜日

第16回助教の回

第16回目となる今回の数理助教の会は, 数理第4研究室の西遼佑さんに「個々の車の振る舞いによる合流渋滞の改善手法の研究」について発表していただきました。西さんは、日本の渋滞研究で有名な西成先生の研究室で学位を取得後、特任研究員として数理4研で研究活動を進めています。

日本で交通渋滞は、1年で11兆円もの経済損失につながっているとの試算があるそうです。そのため、渋滞緩和は重要な社会的課題の一つとなっています。 渋滞緩和では、道路などのインフラ整備(大規模)、信号やスピード制限などの規制を設ける方法、そして、個々の車の振る舞い(例:車間距離)に注目するという3つの切り口があります。インフラ整備や交通規制は大掛かりになりますので、個々の車の動きに注目した、低コストで渋滞緩和に効果を発揮する方法が望まれています。高速道路上でインターチェンジや合流部はゆるやかな上り坂(サグ部)、トンネルについで、渋滞の発生する場所ですが、日本で特に有名なのは小菅ジャンクション(JCT)だそうです。

交通流の研究は意外に歴史が古く、1930年代(日本は昭和初期)のアメリカで既に様々なデータ(速度 - 車間距離など)がとられていたようです。1990年代になると、物理学的な手法の1つとして、一つ一つの車を自己駆動粒子(Self Driven Particle; SDP)とみなし、交通流を多数の自己駆動粒子の集まりととらえて数理モデルで表現する研究が発展しました。また、統計物理の立場で、渋滞している状態と、そうでない状態を相転移ととらえるような研究もあるそうです。

西さんは、上のような渋滞現象の数理モデルでの記述から、さらに一歩進んで、実際に渋滞が起きるケースに注目し、渋滞緩和の具体的な方法を研究しています。私はこの点が非常に優れていると思いました。それでは、西さんたちが注目している小菅ジャンクション(JC)での渋滞緩和のアイディアをみていきましょう。

一般に車線変更が起きるような道路ではジッパー配置が理想的とされています。ジッパー配置とは、互い違いなど隣り合わずに並走している状態で、いったんジッパー配置になればスムーズに車線変更が可能となります。

ところが、下図のような、二つの道路がいったん合流して、すぐに分岐するような場所(例:小菅JCT)では、ジッパー配置になる前に車線変更してくる車は、後ろの車を減速させて、渋滞を発生させてしまいます。そこで、小菅JCTのような合流地点で、図の下側にあるようなオレンジライン(車線変更禁止ライン)を引くことを考えます。合流直後の車線変更を防ぎ、車線間相互作用を活用してジッパー配置を励起するのが狙いです。非常に安価であるというメリットがあります。

図の下側のように、オレンジラインに沿って上下車線の車がジッパー配置になることでスムーズに車線変更がおきます。 さて、このジッパー合流ですが、小菅JCTのように、合流部分が7~800メートルと短い区間では、オレンジラインの長さの調整が大変難しいそうです。また、実際に実験するのも大変コストがかかり危険を伴うため、まずは、数値シミュレーションで解析する必要があります。 西さんは、セルオートマータ(CA)を用いて、2車線の系で自動車がどのようにふるまうかを様々な条件のもと、シミュレーションしています。



 各車は、隣の車線に車があると減速したり、前が空いたら加速するといった複雑なふるまいをしますので、連続時間で表現するのは難しく、まずは、離散的な時間で1マスずつ車を動かしています。上のモデルでは、パラメータ a  (0<a<1) を動かした時にどれだけ短い距離でジッパー配置を達成できるかがわかりました。また、平均場近似による計算でも説明できます。他にもパラメータを調整して様々な結果がわかってきました。しかしながら、このようなシンプルな数理モデルでも、実際の状況に合うパラメータをどうやって推定すればいいかわからないという状況で、応用には様々な課題があるそうです。
私の感想としては、まず、実験的に検証しづらい現象をここまで数学的にモデル化できたというのが驚きです。 これは、他に先行研究のない全く新しい研究で、今後、日本が世界をリードできる分野の一つと感じました。また、1%でも渋滞緩和ができれば、経済効果としては単純計算で1年1000億円になります。道路での実証実験と実用化を目指し、数理モデルでのシミュレーションにある程度の予算と人材を投入すべきと思います。
発表では、既にある小菅JCTの課題をとらえて解決案を探っていましたが、今後の道路網整備や都市計画などへも、非常に有効な知見を与えられるのではないでしょうか。道路のような巨大なインフラでは、ちょっとした設計ミスが将来、多額の国家損失になりえますので、やはり、ある程度の研究予算を投入すべきと思いました。

数理4研 田中冬彦

2012年6月13日水曜日

第15回助教の会


第15回助教の会では公立はこだて未来大学の田中健一郎さんに発表していただきました.田中さんは2007年3月に東京大学大学院情報理工学系研究科数理情報学専攻にて博士号を取得後,東京海上日動火災保険株式会社に就職.2011年4月より公立はこだて未来大学システム情報科学部複雑系知能学科の助教として研究活動をされています.大学院時代の同期である数理4研助教の田中冬彦さんの仲介により今回の発表が実現しました.発表タイトルは「微分作用素に対するHillの方法の収束次数および明示的な誤差上界」ということで,微分方程式に対する数値解法であるスペクトル法,特にHillの方法に関する研究成果のお話でした.

物理的な現象や社会的な現象の中には微分方程式により記述されるものが多くあります.本発表で出てくる波の挙動はその代表的な例といえます.微分方程式の数値解法としてはスペクトル法の他に差分法や有限要素法がよく知られています.スペクトル法の特徴は,他の手法と比べて非常に精度よく解を近似できることです.田中さんの研究は,Hillの方法による近似誤差の収束次数と上界に関する理論的な結果を与えています.特に,誤差上界の評価については,スペクトルの真値を利用しないことで計算可能な上界が得られています.数値手法を実際に利用する上で理論的な評価は重要な指針となります.本発表の主結果(の一部)は次の通りです.

図1:主結果

具体的な問題設定は以下の通りです(図2).まず,対象は非線形波動方程式です.この非線形波動方程式を定常解の周りで線形近似することにより,摂動$\eta$を使って解$\zeta$の安定性を調べることが目的です.ここではさらに$\eta(x,t)=\Psi(x)\exp(\lambda t)$を仮定して,時刻に依存しない関数$\Psi$の微分方程式へと変形します.今,解の安定性は微分作用素$S_p$のスペクトル,特にその実数部分の正負で判定されます.つまり,スペクトルの実部が負であれば安定,正であれば安定ではないということです.以上より,問題は微分方程式$S_p\Psi = \lambda \Psi$を数値的に解くことに帰着し,あとはこの微分方程式を数値的に解けばよいことになります.そして田中さんは,Hillの方法による近似を考え,その収束次数と誤差上界を理論的に導出しています.

図2:非線形方程式の変形


Hillの方法はもともと1886年にHillにより提唱されましたが,有効な数値計算手法として認識されたのはごく最近になってからの話だそうです.Hillの方法による近似は
1. Floquet-Bloch分解
2. Fourier級数による近似
の2段階で行われます.まず,Floquet-Bloch分解の概要を図3に示します.Floquet-Bloch分解をなぜ行うかというと,それは作用素$S_p$が連続スペクトルをもち直接的な解析が難しいためです.Floquet-Bloch分解では,非周期性を指定するパラメータ$\mu$の値に応じて周期解と微分作用素$S_p^{\mu}$を考えます.パラメータ$\mu$の値を変化させることで作用素$S_p$のスペクトルがすべて現れ,$\mu$の値ごとに作用素$S_p^{\mu}$が点スペクトルだけをもつので解析がしやすくなります(図4).

図3:Floquet-Bloch分解(その1)
図4:Floquet-Bloch分解(その2)

Fourier級数による近似では,周期関数を使って周期解を無限級数展開し,$2N+1$項の有限打ち切りによる近似を行います(図5).Hillの方法により,微分方程式$S_p\Psi = \lambda \Psi$という問題がパラメータ$\mu$ごとに$2N+1$次正方行列の固有値問題に帰着されます.

図5:Fourier級数による近似


固有値問題の解ともとの問題$S_p\Psi = \lambda \Psi$のスペクトルとの誤差を評価するために,田中さんは誤差を「固有値問題の固有関数ともとの問題の固有関数との誤差」と「Fourier近似による誤差」の2つに分解しました(図6).後者についてはFourier解析の一般論から誤差の評価が得られます.前者の評価にはGershgorinの定理とDunford積分という2つの道具が使われています.Gershgorinの定理は固有値の存在する範囲を教える定理であり,Dunford積分は複素解析における留数の定理を作用素に対して考えたようなものとのことです.上界の評価においては,収束次数の評価で用いた手法に工夫を加え,スペクトルの真値を直接利用しないことで明示的な上界を与えています.以上より精度誤差の収束次数と上界に関する主結果を得ます.

図6:精度誤差評価のためのアイデア

例題について,近似誤差と田中さんの与えた上界との比較結果も紹介されました(図7).

図7:精度誤差と上界との比較

今回の田中さんの発表では「実際に計算できるもので評価する」という点が強調されていたように思います.誤差の理論的な評価は数値手法に関する重要な指標です.特に,実際に計算して使える指標の存在はやはり心強いものだと思います.大変興味深く面白いお話でした.


数理5研助教 廣瀬善大



2012年6月8日金曜日

第14回助教の会

今年度より助教の会に入会した林です.5月28日に行われた第14回助教の会では廣瀬善大さんに "Dimension Reduction Based on the Geometry of Dually Flat Spaces" (双対平坦空間の幾何学に基づいた次元削減)というタイトルで発表していただきました.

廣瀬さんは昨年度3月に東京大学の数理5研で博士号を取得され,今年度より同研究室で助教に着任されました.

今回のお話は廣瀬さんの博士課程での研究成果をまとめたもので,一言でいうと Least Angle Regression (LARS, 最小角回帰法) の幾何学的な拡張と具体的な確率モデルへの応用について発表されました.発表の前半ではLARSの解説とその拡張について,後半ではグラフィカルガウシアンモデルへの応用が紹介されています.

そもそものモチベーションとして,統計モデルを学習させるときに疎(スパース)な解を得たいというものがあります.例えば,線形回帰の問題では入力Xと出力yの組に対して線形な関係を仮定し,二乗誤差が最小となるような重みパラメータを求めます(最尤推定).この求められた重みの強さによってどの特徴量がyの予測に役立っているかを知ることができますが,重みが疎,すなわち多くの値が0で埋まっていればyと関係する特徴量を簡単に絞り込むことができます.一般的に線形回帰で求められた解は疎にならないため,L1ノルムによる正則化項を付け加える(LASSO)などなんらかの工夫が必要となります.

LARSは線形回帰の解を疎に求めるアルゴリズムとして2004年に提案されました.
  • Efron, Hastie, Johnstone and Tibshirani: Least Angle Regression (with discussion). Annals of Statistics, vol. 32 (2004), pp. 407-499.
LARSによるパラメータ推定の手順は一種の繰り返しアルゴリズムになっています.通常の最尤推定値を$\bar{y}$としたとき,LARSは予測に使う特徴量を一つづつモデルに加えていき,原点から$\bar{y}$にいたるまでのパスを計算します.具体的なアルゴリズムは以下になります.またアルゴリズムの具体例を図1に示します.
  1. すべての特徴量を使わない状況,すなわちすべての重みを0としてスタート.
  2. まだ使っていない特徴量のうち,相関が一番大きい特徴量を基底に加える.
  3. 基底に入っている特徴量と$\bar{y}-\hat{\mu}$との相関を考える.$\hat{\mu}$がその上を動くような直線として,これらの相関が互いに等しいままであるような直線をとる(角の二等分線に相当する直線).使っていない特徴量のうち,これらの相関と同じだけ大きい相関をもつ特徴量が現れたら,2.に移る.
  4. 推定値$\hat{\mu}$が$\bar{y}$に到達するまで2.と3.を繰り返す.
このアルゴリズムをk回目で止めるとk個の非ゼロな重みをもつパラメータを得ることができるため,ユーザは任意の疎な解を得ることができます.
図1:LARSによって得られるパス(赤線)

廣瀬さんはLARSのアイデアを情報幾何の世界で解釈しなおし,他の統計モデルにとって自然な拡張を行いました.LARSはパスを求める際,角の二等分線というユークリッド幾何的な性質を用いており,この性質は線形回帰と相性がいいものですが,一般の統計モデルにとってはそうとは限らないからです.LARSは角度という幾何的性質を用いたアルゴリズムですが,廣瀬さんはこの代わりに距離を用いてLARSの拡張を行いました.具体的にはKL擬距離という確率分布にとって自然な距離を導入し,一般化線形モデルやグラフィカルガウシアンモデル(GGM)に応用されました.(この記事では深くは立ち入りませんが,情報幾何とはユークリッド幾何を一般化したもので,ある確率分布がある一点として表現されるような空間を扱う学術体系です.)

発表の後半ではGGMへの応用が紹介されました.GGMは複数の変数が与えられたとき,ある2変数の間に(他の変数が与えられた下での)条件付き独立性があるかないかを調べるときに用いられるモデルです.具体的にはすべての変数は1つの多次元ガウス分布から生成されたと仮定し,推定された精度行列(共分散行列の逆行列)が0であれば対応する変数間には条件付き独立性がないと考えます.変数の数が増えるにしたがってその組み合わせの数も増加するため,ここでもやはり疎な解が歓迎されます.
図2:提案手法の実験結果
実データを用いた実験もいくつか紹介されました.その1つは学生88人の数学5科目の成績から科目間の条件付き独立性を調べたもので,図2は実際に提案手法によって得られた結果を示しています.2科目間の条件付き独立性の有無がグラフ上の辺によって表現されており,たとえば代数と解析の従属性が強いとされるなど興味深い結果となっています.シミュレーションでは提案手法が既存手法よりも優れた予測精度を達成したことが示されています.

今回の廣瀬さんのお話は,線形回帰モデルにのみ用いられてきたLARSアルゴリズムを一般の統計モデルに対しても使えるよう拡張されたという点で,非常に意義深いと感じました.自分も含め,最近は確率モデルや統計モデルを用いた研究はとても一般的になってきたため,ユーザの立場からもこのような一般化はありがたいのではないでしょうか.ただ,提案手法の計算量はそこそこかかってしまう(GGMの場合,変数の数の二乗のオーダ)と伺ったので,今後は計算の効率化にぜひ期待させていただきたいと思いました.

数理第6研究室 林浩平