tICAにおける一般化固有値問題の解き方

2023年5月2日 18:51

tICA（Time-structure Independent Components Analysis or Time-lagged Independent Components Analysis）は様々な時間スケールの情報が混在している入力値から時間スケール毎の独立成分を抽出するための手法です。生体分子シミュレーションで得られるtrajectory(軌跡)はまさに様々な時間スケールの情報が混在しており，またその中から特定の時間スケールの情報を得たいという場面があり，tICAが有効な手法となることが期待されます。tICA自体については，参考文献１（tICA提案者による解説記事）で詳しく説明されていますので，興味ある方はご参照ください。
MDシミュレーション結果に対してtICAを実行するためには，式(1)の一般化固有値問題を解くアルゴリズムを実装する必要があります。

$$
\begin{align*}
\mathbf{C}(\tau)\mathbf{U}&=\mathbf{C}(0)\mathbf{U}\boldsymbol\Lambda
\end{align*}\tag{1}
$$

本記事では，参考文献２に基づいて式(1)を解くためのアルゴリズムについて解説します。

式(1)内における太文字表記された行列の定義はそれぞれ以下の通りです。
$${\mathbf{C}(\tau)}$$は時間遅延の分散共分散行列(time-lagged covariance matrix)であり，

$$
\begin{align*}
\mathbf{C}(\tau)&:=\langle \mathbf{r}(t)\mathbf{r}(t+\tau)^{\rm T}\rangle\\
\mathbf{r}(t)&= \tilde{\mathbf{r}}(t)-\langle\tilde{\mathbf{r}}(t)\rangle
\end{align*}\tag{2}
$$

と定義されます。$${\langle\cdots\rangle}$$はアンサンブル平均を表し，MD計算の場合は時間平均を代用して値を求めます。$${\tilde{\mathbf{r}}(t)}$$は，適当に併進・回転操作された解析対象部分の配置になります。
$${\mathbf{U}, \boldsymbol\Lambda}$$はそれぞれ一般化固有値方程式の固有ベクトル$${\{\mathbf{u}_i\}}$$，固有値$${\{\lambda_i\}}$$を行列としてまとめたものになります。

$$
\begin{align*}
\mathbf{U}&:=\begin{bmatrix}\mathbf{u}_1&\cdots&\mathbf{u}_N\end{bmatrix}\\
\boldsymbol\Lambda&:={\rm diag}(\lambda_1,\cdots,\lambda_N)
\end{align*}\tag{3}
$$

固有ベクトル行列の分解

固有ベクトルを並べた行列である$${\mathbf{U}}$$を複数の行列に分解することにより，式(1)を標準固有値方程式に帰着させることがポイントになります。天下り的ですが，

$$
\begin{align*}
\mathbf{U}&=\mathbf{W}\boldsymbol\Sigma^{-1}\mathbf{V}
\end{align*}\tag{4}
$$

と分解することで標準固有値方程式に帰着できます。
ここで，$${\mathbf{W}, \boldsymbol\Sigma}$$はそれぞれ$${\mathbf{C}(0)}$$の固有ベクトル$${\{\mathbf{w}_i\}}$$，固有値$${\{\sigma_i\}}$$を行列としてまとめたものになり，主成分分析に対応します。一方，$${\mathbf{V}}$$は標準固有値問題における固有ベクトル$${\{\mathbf{v}_i\}}$$を行列としてまとめたものになります。
式(4)を式(1)に代入すると，

$$
\begin{align*}
\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1}\mathbf{V}&=\mathbf{C}(0)\mathbf{W}\boldsymbol\Sigma^{-1}\mathbf{V}\boldsymbol\Lambda\\
\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1}\mathbf{V}&=\mathbf{W}^{\rm T}\mathbf{C}(0)\mathbf{W}\boldsymbol\Sigma^{-1}\mathbf{V}\boldsymbol\Lambda\\
&=\mathbf{V}\boldsymbol\Lambda\ \ \ \ (\because \mathbf{W}^{\rm T}\mathbf{C}(0)\mathbf{W}=\boldsymbol\Sigma)\\
\end{align*}\tag{5}
$$

式(5)は行列$${\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1}}$$に対する標準固有値方程式になっています。

行列$${\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1}}$$は$${\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1/2}}$$が正方行列，$${\boldsymbol\Sigma^{-1/2}}$$が対角化行列であることを利用することにより，

$$
\begin{align*}
\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1}&=\left(\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1/2}\right)\boldsymbol\Sigma^{-1/2}\\
&=\boldsymbol\Sigma^{-1/2}\left(\mathbf{W}^{\rm T}\mathbf{C}(\tau)\mathbf{W}\boldsymbol\Sigma^{-1/2}\right)\\
&=\boldsymbol\Sigma^{-1/2}\mathbf{W}^{\rm T}\langle \mathbf{r}(t)\mathbf{r}(t+\tau)^{\rm T}\rangle\mathbf{W}\boldsymbol\Sigma^{-1/2}\\
&=\left\langle\left\{\boldsymbol\Sigma^{-1/2}\mathbf{W}^{\rm T}\mathbf{r}(t)\right\}\left\{\boldsymbol\Sigma^{-1/2}\mathbf{W}^{\rm T}\mathbf{r}(t+\tau)\right\}^{\rm T}\right\rangle\\
&=:\left\langle\mathbf{y}(t)\mathbf{y}(t+\tau)^{\rm T}\right\rangle
\end{align*}\tag{6}
$$

となり，$${\mathbf{y}(t):=\boldsymbol\Sigma^{-1/2}\mathbf{W}^{\rm T}\mathbf{r}(t)}$$に対する時間遅延の分散共分散行列として表現することができます。

AMUSEアルゴリズム

参考文献２によりますと，上記の解き方はAMUSEアルゴリズムと称されるらしいです。
アルゴリズムの処理内容は以下の通りです。

$${\{\tilde{\mathbf{r}}(t_1),\cdots,\tilde{\mathbf{r}}(t_M)\}}$$を用意
$${\langle\tilde{\mathbf{r}}(t)\rangle=\frac{1}{M}\sum_{i=1}^M\tilde{\mathbf{r}}(t_i)}$$を計算
$${\mathbf{r}(t)=\tilde{\mathbf{r}}(t)-\langle\tilde{\mathbf{r}}(t)\rangle}$$を計算
$${\mathbf{C}(0)=\langle \mathbf{r}(t)\mathbf{r}(t)^{\rm T}\rangle＝\frac{1}{M-1}\sum_{i=1}^M\mathbf{r}(t_i)\mathbf{r}(t_i)^{\rm T}}$$を計算
$${\mathbf{C}(0)}$$を対角化することで$${\mathbf{W}, \boldsymbol\Sigma}$$を算出
固有値が小さい主成分モードを除去$${\mathbf{W}\rightarrow\tilde{\mathbf{W}}, \boldsymbol\Sigma\rightarrow\tilde{\boldsymbol\Sigma}}$$
$${\mathbf{y}(t_i)=\tilde{\boldsymbol\Sigma}^{-1/2}\tilde{\mathbf{W}}^{\rm T}\mathbf{r}(t_i), i=1,\cdots,M}$$を計算
$${\mathbf{C}^{\mathbf{y}}(m*\Delta t):=\langle \mathbf{y}(t)\mathbf{y}(t+m*\Delta t)^{\rm T}\rangle＝\frac{1}{M-m-1}\sum_{i=1}^{M-m}\mathbf{y}(t_i)\mathbf{y}(t_{i+m})^{\rm T}}$$を計算
$${\mathbf{C}_{\rm sym}^{\mathbf{y}}(m*\Delta t):=\frac{1}{2}\left\{\mathbf{C}^{\mathbf{y}}(m*\Delta t)+\mathbf{C}^{\mathbf{y}}(m*\Delta t)^{\rm T}\right\}}$$を計算（対称行列化）
$${\mathbf{C}_{\rm sym}^{\mathbf{y}}(m*\Delta t)}$$を対角化することで$${\mathbf{V}, \boldsymbol\Lambda}$$を算出
$${\mathbf{U}=\tilde{\mathbf{W}}\tilde{\boldsymbol\Sigma}^{-1}\mathbf{V}}$$を計算
tIモード行列$${\mathbf{G}=\mathbf{C}(0)\mathbf{U}}$$を計算
tIコンポーネント（$${\mathbf{r}(t)}$$をtIモードで分解した際の係数）$${\mathbf{z}(t)=\mathbf{U}^{\rm T}\mathbf{r}(t)}$$を計算

参考文献

この記事が気に入ったらサポートをしてみませんか？