1．は じ め に
近年の科学技術・計測技術の進展により，我々は大量 の高次元データを得ることができるようになった．この 状況が，データ中心科学の勃興の一因であることはいう までもない．得られた高次元データを有効に生かすには， 情報科学と，そのデータ解析の対象となる学問分野が， 図 1 のように緊密に融合した新しい枠組みと方法論が必 要不可欠である．本解説では，対象の学問分野として自 然科学を想定し，この枠組みをデータ駆動科学と呼ぶ． 我々はデータ駆動科学のキーテクノロジーが，情報科 学分野で大きな注目を集めているスパースモデリングで あると考える．スパースモデリングとは，データの出自 は問わず，与えられたデータが疎（スパース）に表現で きることを指導原理としてデータ解析を行うことや，与 えられたデータをスパースに表現できる基底を探ること を目的としたデータ解析の枠組みである．本稿では，ス パースモデリングによる潜在構造の抽出に着目し，自然 科学の広範囲の分野において有効であるデータ駆動科学 の枠組みを紹介する [岡田 14]． 2章では，データ駆動科学の指針を，図 2 に示す David Marrの三つのレベルと照らし合わせて議論する [乾 87, 川人 96, Marr 82, 岡田 14]．Marr は，この三つ のレベルが脳の情報処理だけでなく，情報処理課題を 実行するシステムを理解するために有用であることを主 張している．このことに基づくと，データ駆動科学を Marrの三つのレベルの視点で議論することは重要であ る．3 章では，具体的データ解析の課題として，スペク トル分解 [Nagata 12, 岡田 14] を取り扱う．最後にまと めとして，本稿の内容と，新学術領域研究「疎性モデリ ング」との関連を述べる [岡田 14, SpM-HD3 13]．2．Marr の三つのレベルによる記述
神 経 科 学 者 で あ る David Marr は， 自 身 の 著 書 “Vision”において，複雑な情報処理課題を実行する機 械を理解するには図 2 の三つのレベルを考えることが有 用であると述べ，脳・神経科学における研究指針を示し ている．一方，データ駆動科学においても，複雑な情報 処理課題を実行するので，データ駆動科学を Marr の三 つのレベルの視点で議論することは興味深い． 図 2 の第一のレベルの計算理論では，データ解析に対 応する情報処理の目標，その目標を達成するための方略， それらの目標と方略の適切さなどを議論する．図 1 に示 すように，この計算理論のレベルが，対象となる科学と データ駆動科学をつなぐ要となる．例えば，物質科学や 材料科学を取り扱うマテリアルズインフォマティクスで は，実験計測のデータや，物質の電子状態に関する大規 模計算のデータから，電気伝導度，誘電率，透磁率など を予測することがデータ解析の目標となる（図 2）．そこ では，これらのパラメータの物理的化学的な意味や機序 などの物質科学的な学問背景なしには，データ解析の目 標を設定することは不可能である．この具体例から，計 算理論を構築するには，データ解析に関する知識や技術 だけでなく，対象となる分野に対する深い知識が必要でスパースモデリングによるデータ駆動科学
Data-Driven Science by Sparse Modeling
永田 賢二
東京大学大学院新領域創成科学研究科Kenji Nagata Graduate School of Frontier Sciences, the University of Tokyo. [email protected]
岡田 真人
（同 上）Masato Okada [email protected]
Keywords:
sparse modeling, data-driven science, spectral deconvolution, Bayesian estimation, exchange Monte Carlo method.「データ中心科学」
あることがわかる． 次に計算理論の方略について説明する．一般には，デー タ解析の目標を完全に満たす条件がデータに含まれてい るとは限らない．ここで三次元テレビの計算理論を考え てみよう．ここでのデータは，右目と左目への入力に対 応する，2 台の離れたカメラで撮像された画像である． これら 2 枚の画像だけでは，画像中の視覚物体の奥行き を計算することはできない．その理由は，2 枚の画像の それぞれのピクセルがどのように対応しているかを，画 像情報からだけでは決められないからである．そこで Marrと Poggio は，この対応関係が画像上で滑らかに変 化すると仮定し，奥行き計算のための計算理論を提案し た [Marr 76]．この滑らかさが図 2 の方略である．また， その方略の適切さは，物理量が空間的に滑らかに変化す ることで担保されている．以上のように，これらの目標， 方略，適切さに関する知識は，対象となる科学の知見に 基づいていることがわかる． 第二のレベルは，計算理論の目標を実現する方略を， どのようにして実現するかを問う表現と，その表現をど のように実行するかを問うアルゴリズムの二つから構成 されている．機械学習は，この第二のレベルに対応する． 通常，一つの計算理論に対して，複数のアルゴリズムや 表現が存在する．これら複数のアルゴリズムの正否は， 予測誤差やアルゴリズムの実行時間により評価されるこ とが一般的である．しかし，いくら予測誤差が低くても， 与えられた計算理論の目標や方略に合致しないアルゴリ ズムであれば，それは何の意味ももたない．例えば，図 3のようなノイジーな観測データが与えられたとしよう． この観測データからノイズを除去するためには，多項式 フィッティングするアルゴリズムも考えられるし，スプ ライン補間も考えることができる．フーリエ基底で展開 するということも考えられる．3 章において，図 3 の観 測データの情報処理に関する計算理論を考察し，それに 基づくアルゴリズムを構築する． 第三のレベルは，第二のレベルであるアルゴリズム が，どのようにして物理的に実現されるかを問うハード ウェア実装のレベルである．計算に用いるハードウェア や，利用できる計算時間などの制約から，対象とするア ルゴリズムを実行することが困難な場合，その近似手法 を取り扱うケースは少なくないが，近似を行うことによ る計算結果に及ぼす影響を正確に知ることは重要であ る．近年，行列分解やマルコフ確率場モデルのハイパー パラメータ推定に関する研究から，ベイズ推論を行う際 の近似手法である変分ベイズ法が，計算理論の観点から 無視できない近似誤差を含む可能性があることが示唆さ れている [Nakajima 11, Nakanishi 14]．例えば，変分 ベイズ法では陽に変数選択の事前確率が導入されていな くても，変分ベイズ法の数理的な性質により変数選択が 行われやすい性質をもつことが知られている [Nakajima 11]．このような性質から類推すると，計算理論の実現 のために導入した変数選択の結果として変数が選択され たのか，近似アルゴリズムの副次的効果で変数選択が行 われたかの区別がつかなくなる可能性がある．この可能 性を避けるための最も確実な方法は，現実的な範囲で使 用可能な高速な計算機を用いて，モンテカルロ法などの サンプリング手法を用いて愚直にベイズ推論することで ある．そのためには，並列計算機や GPGPU のハードウェ ア上の特性を良く理解して，モンテカルロ法をいかに実 装するかのテクニックが必須である．この数値的な厳密 計算のみが，計算理論を数理的に定式化した枠組みのぜ ひを直接的に問える手段である．また，これは変分ベイ ズ法などの近似手法が，どの程度の精度が保証されるか に関する指針にもなる． 以上からデータ駆動科学においても，計算論的神経科 学の現状と同様に，これら三つのレベルを俯瞰しながら 研究を進めることが今後，重要になると予想している． 情報処理の目標として，与えられた時間内に結果を出す 必要がある場合は存在する．このとき，近似精度と計算 時間のトレードオフの観点で，最適な近似手法をあらか じめ設定しておくことが重要である．そのようなシステ ムの設計指針を得るためにも，データ駆動科学の三つの レベルを貫いて近似手法を評価することが必要である．
3．スペクトル分解
3･1 スペクトル分解の計算理論 図 3 のような多峰性スペクトルを，ガウス関数に代表 される単峰性の基底関数の線形和で表現するスペクトル 分解について議論する [Nagata 12]． （1） G(x;θ,K) = K k=1 akexp −(x−µk)2 2σ2 k ここで，K はガウス関数の個数であり，μkは k 番目のガ ウス関数の平均値であり，akは k 番目のガウス関数の大 きさである．σ2 kは k 番目のガウス関数の分散に対応する． ここでパラメータセットθをθ＝ {ak, μk, σk}Kk＝1と定義す る．この問題は，与えられたデータを少ない（スパース な）基底関数の線形和で表すという観点で，目的がスパー スである．そのため，スパースモデリングの一つの問題 として捉えることができる． 分光学においては，このような多峰性スペクトルのス 図 2 機械学習と David Marr の三つのレベル [乾 87, 川人 96, Marr 82]ペクトルを構成する個々の単峰性の基底関数は，物質を 構成する電子のエネルギー準位などに対応する．このよ うなスペクトルを観測し，そのデータを解析する目標は， 測定対象である物質の未知の電子状態を知ることであ る．スペクトルは離散的なエネルギー状態に対応してお り，複数のエネルギー状態が同時に観測されている状況 に対応する．図 3 の多峰性スペクトルを単峰性の基底関 数に分解し，それぞれの単峰性の基底関数の位置が個々 のエネルギー状態に対応すると考える．これが方略であ る．その分解の結果，その物質の電子状態が何個のエネ ルギー準位から構成され，それぞれのエネルギー準位に 含まれる電子の状態密度を知ることで，対象の電子状態 を知るという情報処理の目標が達成できる． 3･2 スペクトル分解の表現とアルゴリズム ここで述べた方略を数学的に表現したものが，式（1） となる．ここでガウス関数の個数 K は多峰性スペクトル のピーク数であり，それは測定対象のエネルギー準位の 数に対応する．ガウス関数の平均値 μkは k 番目のエネ ルギー準位の値であり，ガウス関数の大きさ akは k 番 目のエネルギー準位の電子の状態密度に対応する．ガウ ス関数の幅に相当するσkは，装置の観測誤差を表すと ともに，X 線などの入力で励起され空になった k 番目の エネルギー準位が，平衡状態へ緩和していくときの緩和 時間を表現している． 以上の考察から，スペクトル分解の計算理論の数理的 な表現は，図 3 のような多峰性スペクトルを構成する N 個のデータセット D＝ {xi, yi}Ni＝1から，多峰性スペクト ルのピーク数に対応する基底関数の個数 K と，それに伴う パラメータセットθをデータから適切に決めることである． 一般的にパラメータセットθの決定法は，最小二乗法 がよく用いられる．N 個のデータセット D ＝ {xi, yi}Ni＝1 に対して，以下の平均二乗誤差 E（θ, K）を最小にする パラメータθを求める手法である． E(θ,K) = 1 2N N i=1 (yi−G(xi;θ,K))2 （2） ここで，パラメータθだけでなく，ピークの個数 K に ついても，最小二乗法で決定することを考える．この手 法では，ピーク数 K を多くすることで，いくらでも誤差 関数 E（θ, K）を小さくすることが可能であるが，そう して得られたピークはノイズにまでフィットするため， 誤った結果を抽出してしまう． また，仮にピーク数が事前にわかっているとしても， 最小二乗法によりパラメータを最適に決定することは非 常に難しい．それは，図 4 に示すように，誤差関数 E（θ, K）が一般に局所解をもつため，初期値の選び方によっ ては，正しい解が得られなくなるためである．簡単な対 策として，局所解に陥ったら，別の初期値から計算し直 すことが考えられるが，そもそも得られた解が局所解か どうかを判別することが難しい．別の対策として，パラ メータの候補すべてについて誤差関数 E（θ, K）の値を計 算する方法が考えられるが，計算量が膨大になってしま う．仮に，一つのパラメータについて 10 個の候補を設 定したとしても，ピーク数が三つの問題ですら，パラメー タの数は 9 個となるため，109＝1 000 000 000 通りの誤 差関数の計算をしなくてはならない． 以上の考察から，最小二乗法ではピークの個数 K ま で決定できないことがわかり，スペクトル分解の表現と して，適切でないことがわかる．そこで，本稿ではベイ ズ推論に基づくスペクトル分解を提案する．ベイズ推論 では，観測データ y が生成された因果律を確率的に定式 化する．まず確率 p（K）に従って，その対象の物質のエ ネルギー準位の数 K が生成されると考える．次に，その Kに従い，パラメータセットθが条件付き確率 p（θ｜K） で与えられると考える．ここでは簡単のために p（K）と p（θ｜K）は一様であるとする．観測データ y は，式（1） に従う真の値に平均 0，分散 1 のガウスノイズεが重畳 されて観測されるとする． y =G(x;θ,K) + （3） これらの仮定より，観測データ y は，ピーク数 K とパラ メータセットθが与えられた元で，入力 x の条件付き確 率で生成される． 図 3 スペクトル分解 図 4 スペクトル分解における局所解構造
p(y|x,θ,K) = √1 2πexp − (y −G(x;θ,K))2 2 （4） それぞれのデータ（xi, yi）が独立に得られたとすると， パラメータθが与えられた元での，N 個のデータ D の 条件付き確率は， （5） （6） となり，p（D｜θ, K）は誤差関数 E（θ, K）をエネルギー とみなし，N を逆温度とみなした場合のボルツマン分布 に従う．ここで，これまで登場してきたすべての変数の 同時確率 p（D, θ, K）を，以下のように形式的に書き下 すことが可能となる． （7） （8） この同時確率が存在すれば，ベイズの定理と周辺化の手 続きを用いて，データ D と基底関数の個数 K が与えら れた場合の，パラメータセットの事後確率 p（θ｜D, K） を以下のように表現することができる． （9） （10） ここで，p（K）と p（θ｜K）は一様であることを用いた．式 （9）と式（10）より，事後確率 p（θ｜D, K）を最大化するθ をパラメータの推定値とする最大事後確率推定は，式 （2）の最小二乗法から得られるθと一致することがわかる． データ D が与えられた場合の，基底関数の個数 K の 事後確率 p（K｜D）は， p(D,K) = dθp(D,θ,K) ∝ dθexp(−NE(θ,K)) （11） p(K|D) = p(D,K)p(D) ∝ dθexp(−NE(θ,K)) （12） で与えられる．先ほどと同様に，D が与えられたもとで のピーク数 K の推定は，最大事後確率推定を用いて，式 （12）の事後確率 p（K｜D）を最大化する K を推定値とす る．また式（12）の右辺は統計力学の分配関数に相当す る．分配関数を用いて自由エネルギー F（K）を， （13） F(K) = −log dθexp(−NE(θ,K)) と定義する．対数の単調性から，事後確率 p（K｜D）を最 大化する K と自由エネルギー F（K）を最小化する K は 一致する． ここで，誤差関数 E（θ, K）がパラメータθの二次式と して展開できるとすると，式（13）は，ガウス積分を用 いることで （14）
F(K) = NE(ˆθ,K) +3K2 log N +const.
として表現される．ここで，ˆθは誤差関数 E（θ, K）を最 小にするパラメータθである．これを見ると，自由エネ ルギー F（K）は誤差関数の項 NE（ˆθ, K）と，第二項（3K/2） log Nのトレードオフにより成り立つことがわかる．式 （14）はベイズ情報量規準（BIC）として知られている． 式（13）と比較すると，パラメータθの積分をすること なく，θˆの導出のみで利用できる性質をもつ． しかしながら，スペクトル分解の場合，パラメータθ が階層的な構造をもつことや，モデルに含まれる非線形 性などに由来し，鞍点近似で仮定するパラメータθにお ける漸近正規性が成り立たず，自由エネルギー F（K）と BICは一致しないことが知られている．このようなクラ スの確率モデルは，特異モデルと呼ばれ，スペクトル分 解だけでなく，多層ニューラルネットワークや，混合正 規分布，隠れマルコフモデルなど，通常の機械学習で用 いられる確率モデルにも多く存在する．具体的には，特 異モデルにおける自由エネルギー F（K）はデータ数 N が十分多い極限において，以下の漸近型をもつ．
F(K) ≈ NE(ˆθ,K) +λlog N +O(loglog N)（15）
ここで，λは式（14）の係数 3K/2 と一致するとは限ら ず，モデルのピーク数 K だけでなく，データを生成して いる真のピーク数 K0などにも依存する係数である．特 異モデルにおいて係数λを解析的に導出するには，代 数幾何学の手法である特異点解消を用いる必要があり [Watanabe 01]，実行することが非常に困難である．また， 数値的に求めようとしても，複雑な多重積分の計算を要 するため，通常の数値積分では不十分である． 3･3 スペクトル分解のアルゴリズム： 交換モンテカルロ法の導入 3･2 節のスペクトル分解のベイズ推論を実行する際に は二つの困難が存在する．一つは，式（2）や式（10） の誤差関数 E（θ, K）が局所解をもつために，誤差関数 を最小にするパラメータ ˆθを求めることが難しいことで ある．もう一つは，式（12）や式（13）のθに関する数 値積分である． これらの問題を解決するため，交換モンテカルロ （EMC）法を用いる [Hukushima 96]．交換モンテカル ロ法はマルコフ連鎖モンテカルロ（MCMC）法の一つ であり，物性物理学のスピングラスと呼ばれる多峰性の エネルギーをもつ系を統計力学を用いて数値的に研究す る際に提案された手法である．ボルツマン分布の対応か ら，式（9）の事後確率 p（θ｜D, K）に逆温度β＝1/T を 導入し，以下の確率分布を考える． pβ(θ|D,K) ∝ exp(−NβE(θ,K)) （16） p(D|θ,K) = n i=1 p(yi|xi, θ,K) ∝ exp (−NE(θ,K)) p(D,θ,K) = p(D|θ,K)p(θ|K)p(K) ∝ exp (−NE(θ,K)) p(θ|D,K) = p(D,θ,K)p(D,K) ∝exp(−NE(θ,K)) −log p(θ|D,K) = NE(θ,K) +const.
温度 T が違う系を複数用意し，図 5 のように同時並列 にサンプリングを行う．具体的なアルゴリズムは，以下 の 2 種類の状態の更新を考え，これらを交互に実行する ことで定義される． （1）それぞれの分布におけるサンプリング メトロポリス法などの従来の MCMC 法により，そ れぞれの確率分布 pβ（l θl｜D, K）からのサンプリングを 並列に実行する． （2）二つの分布間における確率的な交換 上記の操作に加えて，適当なステップごとにサンプ ルθlとθl＋1を以下の確率 u で交換する． （17） （18） （19） これら複数の系をレプリカと呼ぶ．シミュレーテッドア ニーリングでは，温度Tを高温から低温に徐々に下げて（ア ニールして）いく．一方，EMC 法では，図 6 に示すように， 各温度で同時並列に行う MCMC 法の途中に，隣り合っ た温度間で確率的にレプリカを交換することで，アニー リングだけではなく，高温化の効果を取り入れ，局所解 から脱却し，効率的に最小解を探索することができる． 以下に示すように，複数温度でモンテカルロ法を行う EMC法の利点は，効率的に最小解が探索できた時点で， ピーク数の決定に必要な自由エネルギー F（K）も同時 に計算できる点である．ここで逆温度βでの自由エネル ギー f（K; β）を定義すると， （20） （21） （22） となる．式（22）は，自由エネルギーの逆温度βによる 微分が誤差関数 E（θ, K）の期待値になることを示して いる．EMC 法では，複数温度で MCMC 法を行ってい るので，複数温度での誤差関数 E（θ, K）の期待値はすで に求まっている．これらを用いて，図 7 に示すように， 式（21）の積分を数値的に行うことにより，自由エネル ギー F（K）を数値的に計算できる． 実際に交換モンテカルロ法を実行する際には，逆温度 βの設計が重要である．EMC 法のシミュレーションに おいて，逆温度βの設計は，どの程度交換が棄却されず に行われるかを支配しており，最適解への緩和効率を決 めることになる．我々は，この交換頻度に関する平均的 な挙動を解析的に導出することに成功し，その結果，驚 くべきことに，平均的な交換頻度は，式（15）の係数λ の関数として表現されることが明らかになった [Nagata 08]．これは，交換モンテカルロ法の挙動とベイズ推定 図 5 交換モンテカルロ法の概念図．各温度におけるエネルギー描像と各温度におけるフィッティングの例 図 6 交換モンテカルロ法のアルゴリズムの様子． 通常の MCMC 法の並列実行に加えて，確率的な状態の 交換を行う 図 7 交換モンテカルロ法による自由エネルギーの計算． 各棒グラフが，各温度ごとの誤差関数の期待値と して計算できる u = min(1,v) v = pβl+1(θl|D,K)pβl(θl+1|D,K) pβl(θl|D,K)pβl+1(θl+1|D,K)
=exp((βl+1−βl)(E(θl+1,K) −E(θl,K)))
f (K;β) = −log dθexp(−NβE(θ,K)) F(K) = f (K;1) = 1 0 dβ ∂f (K;β) ∂β ∂f (K;β) ∂β = dθNE(θ,K)exp(−NβE(θ,K)) dθexp(−NβE(θ,K))
におけるモデル選択の間に共通した性質を備えているこ とを表す重要な結果と解釈できる．また，この性質を利 用し，交換モンテカルロ法のシミュレーション結果によ り計算された交換頻度の平均値から係数λを導出するア ルゴリズムも開発している [Tokuda 13]． 3･4 スペクトル分解の適用例 図 8 は，これまで説明したスペクトル分解のベイズ推 論の枠組みを，図 3 のデータに関して適用した結果であ る．図 3 は人工データであり，K ＝ 3 を用いてデータは 生成されている．図 8（b），（c），（d）は，それぞれ，フィッ ティングに用いるピーク数を K ＝ 2，K ＝ 3，K ＝ 4 と変 化させたときの，θの推定値を用いて計算した結果であ る．図 8（a）は，EMC 法により数値的に求めた自由エ ネルギー F（K）のピーク数 K 依存性である．図に示す ようにピーク数が K ＝ 3 のときに自由エネルギーが最小 になり，K＝3 が正しく推定された． この様子を図 8（b），（c），（d）を用いて説明する． 自由エネルギーは，誤差関数 E（θ, K）の期待値に対応 するエネルギーと用いたモデルの複雑さを表すエントロ ピーの二つの項からなる．図 8（b）の K ＝ 2 では，デー タをフィットできずに誤差が大きく，エネルギーが高い． 一方，図 8（d）の K＝4 では，差は K＝3 と同等であるが， パラメータの次元数が高くなったことにより，誤差最小 を示すパラメータの密度が相対的に低くなり，結果とし てエントロピーが減少する．その結果，図 8（c）の K ＝ 3のときが，エネルギーとエントロピーのつり合いがと れたモデルとして選択されたことがうかがえる． 3･5 事後確率分布の様子 3･2 節で述べたように，スペクトル分解の問題におい て，BIC と自由エネルギー F（K）は一致しない．これは， データ数 N が大きいときにパラメータθ の事後確率 p（θ｜D, K）がガウス近似できないことに起因する．そ の様子を，実際のシミュレーション結果から事後確率を 可視化することで様子を見てみよう． 図 9 は，その様子を表した図である．図 9（a）が， 用いたデータであり，ピーク個数が二つのスペクトルに ノイズを重畳することで人工的に生成したデータであ る．このデータについて，フィッティングに用いるピー ク数を K＝2，K＝3 とした場合のベイズ推論を実行した． 図 9（b），（c）はそれぞれ，K ＝ 2，K ＝ 3 でのピーク位置 μkのサンプリング結果のラスタープロットを表す．スペ クトル分解のモデルでは，ピークに関する入替え対称性 があるため，K ＝ 2 では 2 か所，K ＝ 3 では 6 か所の部 分にサンプルが集まっていることがわかる．しかしなが ら，その形状は両者で大きく異なる．K ＝ 2 では，ある 1点回りでサンプリングされている．この場合，事後分 布をガウス近似でき，BIC と自由エネルギー F（K）は 漸近的に一致する．一方で，K＝3 では，十字型のサン プリングが 6 か所に点在し，全体として複雑な多様体を なしており，どの部分をピックアップしてもガウス近似 できないことが見て取れる．この状況では，自由エネル ギー F（K）と BIC が一致しない．今回は，ピーク位置 μkのパラメータだけに着目したが，実際には，ピーク強 度 akや幅 σkなどもあることを考慮すると，より一層複 雑な構造をしていることが予想される． この結果を見ると，図 9（b）のような状況のみ BIC を利用することができるように思われる．しかしながら， BICが利用できるかどうかを判別するにはフィッティン グに用いるピーク数 K だけでなく，真のピーク数 K0も 事前にわかっている必要がある．一般に K ＝ K0の場合， 自由エネルギー F（K）の係数λは，λ＝3K/2 となり， BICを利用することができる．一方，K ＞ K0の場合，λ＜ 3K/2となり BIC を利用すると誤った解釈を与える危険 性がある．このように，BIC を利用可能かを事前に判断 図 8 図 3 のベイズ推論の結果． （a）自由エネルギー F（K）の値であり，（b），（c），（d）はそれぞれ，ピーク数 Kが K ＝ 2， K ＝ 3， K ＝ 4 の場合のフィッティング結果である
することは困難であり，現状では，交換モンテカルロ法に より求める方法が精度の観点で優れている [Nagata 12]． 3･6 自然科学へのスペクトル分解の応用 以上のように，ベイズ推定によるスペクトル分解を用 いることで，ピークのパラメータθだけでなく，ピーク の個数 K も与えられたデータのみから推定することが可 能になる．現在，さまざまな自然科学分野においてスペ クトル分解の応用が進んでいる． 惑星科学の分野では，可視光の反射スペクトル解析 に，スペクトル分解が利用可能である．ここでの計算の 目標は，月や小惑星などにおける可視光の反射スペクト ルを計測することで，対象とする惑星の鉱物分布を知る ことである．カンラン石は，惑星の主要鉱物の一つであ り，鉄（Fe）成分をもつファイアライトとマグネシウム （Mg）成分をもつフォルステライトの連続固溶体である． [Sunshine 98]では，カンラン石の Mg/Fe の組成比に関 して，反射スペクトルが連続的に変化することに着目し， 反射スペクトルデータから Mg/Fe の組成比の抽出を試 みている．図 10（a），（b）内の実線は，それぞれ，フ ァイアライトとフォルステライトの反射スペクトルデー タである．それぞれ共通して，1.0 μm辺りに，下に凸 の複雑な吸収帯が存在することがわかる．この複雑な吸 収帯を下に凸のガウス関数の足し合せにより分解するこ とにより，複雑な吸収帯の定量的評価が可能になる．す なわち，方略として，Mg/Fe の組成比ごとにスペクトル 分解を適用し，組成比の変化に対して抽出されたガウス 関数の変化を調べることが考えられる． 具体的な定式は，以下のとおりである [Sunshine 98]． logR(x) = c0+c1/x + K k=1 akexp −(x −µk)2 2σ2 k （23） ここで，入力 x が波長（wavelength）を表し，それ に対する反射率 R（x）の自然対数 log R（x）が出力であ る．右辺第 1 項と第 2 項は，バックグラウンド成分であ り，第 3 項がガウス関数の線形和である．ここで，下に 凸のガウス関数により回帰するため，強度パラメータ ak について，ak 0 という制約がある． 図 10（a），（b）内の点線は，ベイズ推論により導 出された最適なピーク数におけるスペクトル分解の結 果である．それぞれの図の上部の点線が，各ガウス関 図 9 事後分布の形状． （a）はシミュレーションに用いたデータ．（b），（c）は，（a）のデータに対して，K ＝ 2，K ＝ 3 の設定で交換モンテカルロ法を実行した結果得られたピーク位置 μkのサンプリング結果．一般 に，真のピーク数 K0とピーク数 K について，K0 ＜ K のときに，ガウス近似ができなくなる 図 10 カンラン石における反射スペクトルのスペクトル分解． （a）（b）は，それぞれカンラン石の Mg/Fe のエンドメンバであるファイアライトとフォルステライトの反射スペクトルと ピーク分離結果．（c）分解されたピーク構造に基づいた三つのピーク位置と Mg/Fe 組成比の関係 [Sunshine 98]
数を表し，図中央の点線がバックグラウンド成分であ る．両者の図から，複雑な複合吸収帯を三つのガウス関 数で分離することが適切であることがわかる．これは， [Sunshine 98]において，惑星科学者の直感により恣意 的に決定したピーク数と同じである．また，図 10（c）は， Mg/Fe組成比と抽出された三つのガウス関数の中心位置 の関係を表した図である．両者には線形な関係があり， スペクトル分解の結果により Mg/Fe 組成比を抽出する ことの可能性を示唆している．
4．ま と め
本稿では，スパースモデリングによる潜在構造の抽出 に着目し，自然科学の広範囲の分野において有効である データ解析の方法論を紹介した．まず，David Marr の 三つのレベルを紹介し，データ駆動科学の方法論構築の 指針について議論し，その後，ベイズ推論によるスペク トル分解を紹介した． 通常，データ駆動科学の計算理論は，対象とするデー タの獲得者の目的や意図など，データの出自に依拠する． 一方，これまでのデータ駆動科学のアルゴリズムの応用 範囲の広さなどの汎用性の高さを考えると，多様に見え るデータの背後に，アルゴリズムを普遍的に適用できる 計算理論的な背景の存在が示唆される．近年のデータ科 学を背景とし，Science でもデータ科学が特集として取 り上げられている．その中に，天文学における高次元デー タ解析手法が，全く対象とスケールが異なる生命科学で も有効に働くことが報告されている [Reed 11]． この事例を受け流すことなく，分野を超えたアナロ ジー・普遍性への探究心とともに，多様な視点を貫く普 遍的な原理に基づく，新たな学問体系を提案する好機と 捉えることもできる．我々は，その普遍的な計算理論を 高次元データのスパース性であると考えた．データのス パース化を行った結果を，データ獲得者の学問的背景に 基づき判断することで，一連の枠組みにおける計算理論 の適切さを判断する．このように，自然科学を応用範囲 としたデータ駆動科学の計算理論の検証には，実験・計 測・情報科学の分野融合・連携が必須である．我々は， こうした幅広い自然科学の実験・計測研究者と情報科学 の連携により，計算理論とアルゴリズムのレベルをさま ざまな分野・対象においてループさせることで，データ 解析における普遍的な原理を明らかにし，高次元データ 駆動科学と称すべき新学術領域の創成を目指している [SpM-HD3 13]． 謝 辞 本研究は，科学研究費補助金・新学術領域研究 25120001， 25120009の援助を受けた．また，本稿の執筆にあたり， 東京大学の大道勇哉氏，中西（大野）義典氏，徳田 悟氏， 村田 伸氏には，内容に関して，深く議論させていただい た．特に，徳田 悟氏には，図 9 に関するシミュレーショ ンの実施など，多大なるご協力をいただいた．ここに， 感謝の意を表す．◇ 参 考 文 献 ◇
[Hukushima 96] Hukushima, K. and Nemoto, K.: Exchange Monte Carlo method and application to spin glass simulations,
J. phys. Soc. Jpn., Vol. 65, pp. 980-988（1996）
[乾 87] Marr, D. 著，乾 敏郎，安藤広志 訳：ビジョン，視覚の計 算理論と脳内表現，産業図書（1987）
[川人 96] 川人光男：脳の計算理論，産業図書（1996）
[Marr 76] Marr, D. and Poggio, T.: Cooperative computation of stereo disparity，Science, Vol. 194, 283-287（1976）
[Marr 82] Marr, D.: Vision, MIT Press（1982）
[Nagata 08] Nagata, K. and Watanabe, S.: Asymptotic behavior of exchange ratio in exchange Monte Carlo method, Neural
Networks, Vol. 21, pp. 980-988（2008）
[Nagata 12] Nagata, K. Sugita, S. and Okada, M.: Bayesian spectral deconvolution with the exchange Monte Carlo method, Neural Networks, Vol. 28, pp. 82-89（2012）
[Nakajima 11] Nakajima, S. and Sugiyama, M.: Theoretical analysis on Bayesian matrix factorization, J. Machine
Learning Research, Vol. 12, pp. 2579-2644（2011）
[Nakanishi 14] Nakanishi-Ohno, Y., Nagata, K., Shouno, H. and Okada, M.: Distribution estimation of hyperparameters in Markov random field models, J. Phys. A, Vol. 47, 045001（2014） [岡田 14] 岡田真人：スパースモデリングとデータ駆動科学，人工
知能学会全国大会，2014, 1F5-OS-06b-3（2014）
[Reed 11] Reed, S.: Is there an astronomer in the house?, Science, Vol. 331, pp.696-697（2011）
[SpM-HD3 13] 科学研究費補助金新学術領域研究：スパースモデリ ングの深化と高次元データ駆動科学の創成，http://sparse-modeling.jp/
[Sunshine 98] Sunshine, J. M. and Pieters, C. M.: Determining the composition of olivine from reflectance spectroscopy, J.
Geophysical Research, Vol. 103, pp. 13675-13688（1998） [Tokuda 13] Tokuda, S., Nagata, K. and Okada, M.: A numerical
analysis of learning coefficient in radial basis function network, IPSJ Trans. on Mathematical Modeling and Its
Applications, Vol. 7, pp. 15-21（2013）
[Watanabe 01] Watanabe, S.: Algebraic analysis for nonidentifiable learning machines, Neural Computation, Vol. 13, pp. 899-933 （2001） 2015年 1 月 29 日 受理
