xiangze's sparse blog

機械学習、ベイズ統計、コンピュータビジョンと関連する数学について

力学系の性質とその分類、予測への応用に関する論文紹介(2)

論文 物理 力学系 機械学習

機械学習の分野で力学系の知見を取り入れること、あるいは力学系研究において機械学習や統計的予測の手法を用いることは有用です。
前のエントリで紹介した２つの論文では時系列を生み出す力学系の方程式を求めるようなことはせず、その不変量から分類、予測を行っていました。
今回は式が既に知られている力学系の性質(分岐の様子)をパーティクルフィルタによる予測で再現することを試みた論文を紹介します。

Identifying dynamical systems with bifurcations from noisy partial observation(arxiv) (Phys. Rev. E 87, 042716)

概要

微分方程式の式の構造の予測をパーティクルフィルタを用いて行い、その分岐構造が再現されていることを見るもので、細胞分裂周期に関する２種類のモデルに対して双子実験の形でそれを行っています。

分子生物学において細胞分裂などの細胞内の特定の機能に関わる遺伝子、分子の種類は一般的には非常に多く、それらを取り入れた(確率)微分方程式のモデルは非常に複雑なものになります。
しかしながら外部から観測可能な現象は往々にして少数の自由度であり、本質的な現象も少数の変数のみで記述できる場合が多く、また現象理解の観点からも自由度の数は限られます。
そのような状況で数式に基づいたモデルが現象を表現するのに適切なものであるかを判定する基準としてこの論文では力学系の分岐に着目しています。非線形力学系の線形力学系で近似できない性質はパラメータの値変化に応じて現れる分岐現象にあると言えるからです。
この論文ではパーティクルフィルタを用いて力学系の挙動を学習することで既存の数理モデルの持つ分岐現象をとらえることができるか否かを検証しています。

手法

状態空間モデル
x_i^{t+1}=x_i^{t}+\Delta t f_i({x^t},s)
y_i^r=g_i(x_i^r)+\eta_i \phi_i^r
(xは隠れた変数、yは観測変数,t,rはそれぞれシミュレーション、観測の時間間隔、sは外部入力)
で推定したいモデルfが非線形力学系であり、その形状を
f_i({x_j^t},s)=\sum_n k_i^n f_i^n({x_j},s)
という線形の重ね合わせで表現しようとする場合を想定します。
論文では細胞分裂に関わる分子量の変動を表す２つのモデル Novak and Tyson modelとFerrell modelという２つの異なるモデルに対して数値実験を行っています。この２つのモデルはパラメータの変化に伴う分岐の構造がsaddle-node型分岐とsupercritical Hopf型分岐と異なっています(下図,論文のFig3 (a),(b))。


saddle-node型分岐(左)では固定点、リミットサイクル間の変化は不連続的で有限の大きさのリミットサイクルが突然現れるのに対し、supercritical Hopf型分岐(右)では連続的で固定点は極小の大きさのリミットサイクルに変化し、パラメータの変化に伴ってだんだん大きくなっていく。

Novak and Tyson modelとFerrell modelはそれぞれ9種類の分子を変数とした微分方程式であり、論文のSupplemental Materialに式が記載されています。
観測変数y_iとして論文では既存の実験結果に基づいてCdc2,Cyclinの２つの物質の量としています。
f_iの関数型は生化学でよく用いられるMichaelis-Menten型関数の組み合わせではなく、多項式
f_i({x_j^t},s)=k_i^1+k_i^2 x_1^t + k_i^3 x_2^t +\cdot\cdot\cdot+ k_i^{N_i} (x_D^t)^M
の重ね合わせで近似しており、モデルに対する事前知識をなるべく入れない状況を想定していると言えます。
パラメータの最適化としてEMアルゴリズムを用いており、そのExpectation stepでは時系列全体に対して初期値を変動させることを繰り返すようなパーティクルフィルタを用いて事後確率
p(X|Y,\theta)
を求めています。 \thetaは重ね合わせの係数 k_i^n,ノイズの分散 \sigma_i, \eta_iなどのパラメータを成分とするベクトルです。
Maximazation stepでは各イテレーションiに対し対数尤度の平均が極大となる場合
\theta_i=argmax_{\theta} [ Q(\theta) \ , Q(\theta,\theta_o)=_{p(X_o,Y_o,\theta_o)}]
すなわち
\frac{\par Q(\theta_i)}{\par \theta_i} =0
となるような条件から k_i^nなど各パラメータの最適値を求めています。
このアルゴリズム、式の詳細も論文のSupplemental Materialに詳細に記載されています。
実験で制御可能なのは外部入力sですが、固定点、リミットサイクルの分岐点をまたぐような3個のsの値に対してパーティクルフィルタでの学習を行い、分岐の様子の学習を行っています。
変数が多ければ多いほど時系列のフィッティングの精度は増しますが、未知のデータへの適合性が悪くなる過学習の現象が見られるので、ここではAIC, BICを多項式の次数Mごとにとってその値の極小値を採用しています。

結果

学習で得られた力学系が元の力学系の分岐構造を再現していることは入力sを変動させた場合のアトラクタの変化によって見ることができます。学習に用いる外部入力sの値を2種類だけにしてもTyson, Ferrell両モデルの分岐の特徴をとらえることができているのは特筆すべき点です。
(論文のFig. 3(g),(h) )

力学系の特徴的な性質が低次元であっても保存されていることはCdc2-Cyclinに相当する２つの変数の相図をプロットし、学習で得られた力学系の軌道を見るとわかります。図(論文のFig 4)ではヌルクライン(微分方程式で特定の変数の時間微分=0の条件を満たす曲線)も示されています。

対比としてTyson, Ferrell両モデルでCdc2、total Cyclin以外の変化速度の遅い変数を定数として固定した場合も重ね書きされています。
これは断熱近似と呼ばれ物理学や化学ではよく用いられる手法なのですが、ここではヌルクラインの形(図の破線)がもとのモデルの出力(図の散布点)と一致している学習して得られたモデルのもの(実線)とは全然違う形になってしまっています。特に紫色の線で示されているtotal Cyclinの乖離が大きく、これはtotal Cyclinの変化の速度が他の省略された変数に比べ十分に速くないことに起因すると考えられます。(このような時間スケールの近接は生物学において単純に物理学の知見を用いることが難しい理由の一つだと思います。)

Predicting Catastrophes in Nonlinear Dynamical Systems by Compressive Sensing

Phys. Rev. Lett. 106, 154101 arxiv

タイトルの通り圧縮センシングの手法を用いてHenon mapやLorenz,Rossler系のような少数自由度カオス力学系の方程式と動きの予測をしています。

本論文ではcatastropheすなわち時系列の性質の急激な変動とその予測に主眼がおかれています。これは力学系においてはパラメータの変化に伴う分岐に相当し、現実世界においては分岐点近傍においてノイズによる変動が不連続な変化によって増幅するような場合が想定されます。
このような力学系的見方からcatastropheの予測は局所的な線形予測では不十分で時系列全体を記述する数式の推定が必要という立場です。
そして時系列の挙動を微分方程式(力学系)として表現した場合その挙動に影響する項は少数であるという仮定をおけば疎なベクトルに対する圧縮センシングを用いることで高速で推定することが可能であるという研究です。
ここでは力学系の変数の時間微分であるベクトル場の係数をベクトル化したものに対して圧縮センシングを用いていて、変数値そのものの観測値からのずれを直接には用いていないことになります。
圧縮センシングとは
X = G \cdot a , X\in R^wを満たすようなv次元ベクトル a \in R^vが疎でありかつv>>w の場合に
argmin \sum_{i=1}^v ||a_i || , X = G \cdot a , X\in R^w
を線形計画法を用いて解くことができるというものです。求めるべきベクトルaはここでは微分方程式右辺の係数
\frac{dx}{dt}=[F(x)]_j =\sum_{l_1}^n \sum_{l_2}^n \cdots \sum_{l_m}^n (a_j)_{l_1, \cdots , l_m} x_1^{l_1} x_2^{l_2} \cdots x_m^{l_m}
を並べたものとしています。
解かれるべき問題X=GaはここではベクトルXを変数xの時間微分の観測値を並べたものとして
(\dot{x}(t_1), \cdots ,\dot{x}(t_w) )^T =(g(t_1) \cdots g(t_w))^T \cdot a_1
と書かれます。
行列Gを構成するベクトルg(t_i)は各t_iに対して変数 x_1(t_i), x_2(t_i), \cdots ,x_m(t_i) (x(t_i),y(t_i),z(t_i)) の値を列のラベル (l_1, \cdots l_m)に応じて代入することで得られます。
左辺Xはx_iごとに異なるものがあり、m次元力学系ではｍ個の式が得られることになります。論文での例はすべて３次元以下の場合です。

代数方程式で表されるHenon map,Lorenz系、Rossler系に対して適用し、もとの力学系に類似した分岐図が得られること(fig. 2,fig.3)、予測誤差のaの要素数に対する関係にしきい値があり、一定以上の係数を仮定すれば元の非線形力学系の動きが再現できると言えることを示しています(論文のfig. 4)。

最後に論文では多自由度系、確率的な系への適用について触れていますが、(aの要素数が組み合わせ的に増大するので)単純な適用はできず、ベイズ的手法を用いるのがよいのではないかとしています。

また弱点としてはレプリケーター方程式のような指数関数を含む力学系への適用は一見すると難しいように見えることなどが考えられます(レプリケーター方程式自体は変数変換により指数関数を含まない形にできます)。

感想

非線形力学系の挙動、とくに予測不可能性をもつカオス力学系の動きをどの程度パーティクルフィルタでとらえることができるかは個人的興味がありました。本論文の主張によれば少なくとも分岐の構造をとらえることに関しては十分な性能を持っていることが示されました。
生物学などのモデルをつくるという実際的な目的での有用性がしめされたことも見逃せません。
本論文では力学系の初期値を振り、対数尤度の最大値となる値を求めることでパラメータを決定していましたが、パラメータ空間内でパーティクルフィルタを実行する方法もあるのかなと思いました。生物学の文脈においては進化を模倣することに相当します。　

参考文献
    • Parameter Estimation of In Silico Biological Pathways with Particle Filtering Towards a Petascale Computing (pdf)

「データ同化入門」でも１章が割かれて説明されている研究です。分子生物学における複雑な情報伝達経路の予測、パラメータ探索をパーティクル１億個を用いた力技で行っています。