部分空間同定法の確率的解釈

はじめに

システム同定について学んでいたところ，部分空間同定法の確率的解釈という話題があった。今回読んだ論文は，上甲昌郎他著「局所線形アライメントによる非線形動的システムの学習法」(人工知能学会論文誌，2011)である。

cir.nii.ac.jp

本論文では部分空間同定法を確率的に解釈することによって，機械学習の手法が適用しやすくなり，非線形動的システムの学習にも活用できるという説明がなされていた。本記事は，部分空間同定法の確率的解釈について調べた内容である。

本記事を読むことで得られること

本記事を読むことで得られることは，主に以下の内容である。

正準相関分析に基づく部分空間同定法の導出
確率モデルを用いた，部分空間同定法の確率的解釈
状態空間モデル(前向きマルコフモデル・後ろ向きマルコフモデル)と正準相関分析の関係

はじめに
- 本記事を読むことで得られること
- 目次
対象とするシステムと目的
時系列データのデータセットの作成
- 特徴抽出に向けた時系列データの変形
- 過去の領域と未来の領域のデータ作成
正準相関分析に基づく部分空間同定法
部分空間同定法の確率的解釈
まとめと感想
- 参考サイト

対象とするシステムと目的

この論文で対象とするシステムは，以下のような線形状態空間モデルである(本論文の(5)式)。

$\begin{align} &\boldsymbol{x} (t+1) = \boldsymbol{A} \boldsymbol{x} (t) + \boldsymbol{w} (t) \\ \\ &\boldsymbol{y} (t) = \boldsymbol{C} \boldsymbol{x} (t) + \boldsymbol{v} (t) \\ \\ &\boldsymbol{x} (t) \in \mathbb{R}^d, \boldsymbol{y} (t) \in \mathbb{R}^p \\ \end{align}$

本論文では，非線形な動的システムの状態ベクトル $\boldsymbol{x}^{(t)}$ と，システム行列 $\boldsymbol{A}, \boldsymbol{C}$ を，データから推定することである。
なお本論文では，複数の状態空間モデルからなる混合状態空間モデルを想定しているが，本記事では1つの状態空間モデルを扱う。

時系列データのデータセットの作成

特徴抽出に向けた時系列データの変形

本論文では，離散時間のベクトル過程 $\boldsymbol{y}^{(\tau)} , \tau = 0, \pm 1, \cdots$ が得られているものとする。

状態ベクトルの推定に，部分空間同定法を用いるため，観測データから2種類の時系列データを作成する。時刻 $t$ を境に，過去の領域(領域"p"(past))と未来の領域(領域"f"(future))に分け，長さ $k$ の部分列を作成する。

$\begin{align} \check{\boldsymbol{y}}_k (t) \equiv \begin{bmatrix} \boldsymbol{y} (t-1) \\ \boldsymbol{y} (t-2) \\ \vdots \\ \boldsymbol{y} (t-k) \\ \end{bmatrix} \in \mathbb{R}^{pk} , \quad \boldsymbol{y}_k (t) &\equiv \begin{bmatrix} \boldsymbol{y} (t) \\ \boldsymbol{y} (t+1) \\ \vdots \\ \boldsymbol{y} (t+k-1) \\ \end{bmatrix} \in \mathbb{R}^{pk} \\ \\ \end{align}$

過去の領域のデータである $\check{\boldsymbol{y}}_k (t)$ は，要素の上から下に向かって，未来から過去に向かってデータが並んでいる。
一方で未来の領域のデータである $\boldsymbol{y}_k (t)$ は，要素の上から下に向かって，過去から未来に向かってデータが並んでいる。

過去の領域と未来の領域のデータ作成

有限長のデータ列から過去の領域と未来の領域のデータを作るには， $\tau = 1,...,N+2k-1$ の計 $N+2k-1$ のデータを準備したうえで，長さ $k$ の部分列のベクトルを列ベクトルとする行列 $\boldsymbol{Y}_p$ と $\boldsymbol{Y}_f$ を次のように構成する。

$\begin{align} \boldsymbol{Y}_p &\equiv \begin{bmatrix} \check{\boldsymbol{y}}(k+1) & \check{\boldsymbol{y}}(k+2) & \cdots & \check{\boldsymbol{y}}(N+k) \\ \end{bmatrix} \\ \\ &= \begin{bmatrix} \boldsymbol{y} (k) & \boldsymbol{y} (k+1) & \cdots & \boldsymbol{y} (N + k -1)\\ \boldsymbol{y} (k-1) & \boldsymbol{y} (k) & \cdots & \boldsymbol{y} (N + k - 2) \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{y} (1)& \boldsymbol{y} (2)& \cdots & \boldsymbol{y} (N) \\ \end{bmatrix} \in \mathbb{R}^{pk \times N} \\ \\ \boldsymbol{Y}_f &\equiv \begin{bmatrix} \boldsymbol{y}(k+1) & \boldsymbol{y}(k+2) & \cdots & \boldsymbol{y}(N+k) \\ \end{bmatrix} \\ \\ &= \begin{bmatrix} \boldsymbol{y} (k+1) & \boldsymbol{y} (k+2) & \cdots & \boldsymbol{y} (N + k)\\ \boldsymbol{y} (k+2) & \boldsymbol{y} (k+3) & \cdots & \boldsymbol{y} (N + k +1) \\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{y} (2k)& \boldsymbol{y} (2k+1)& \cdots & \boldsymbol{y} (N + 2k - 1) \\ \end{bmatrix} \in \mathbb{R}^{pk \times N} \\ \\ \end{align}$

このようにデータ行列を構成すると，過去の領域のデータはテプリッツ行列，未来の領域のデータはハンケル行列になっている。

なお本論文中では明記されていないが，未来の領域のデータ・過去の領域のデータはそれぞれ中心化，すなわち $E [ \boldsymbol{y}_k (t)] =0, E[ \check{\boldsymbol{y}}_k (t) ]=0$ になるように，データから平均値が引かれている。

正準相関分析に基づく部分空間同定法

ここからは，本論文における式展開を追ってみる。なお本論文では，

$\boldsymbol{\mu}_p$ : 過去の部分列のサンプル平均
$\boldsymbol{\mu}_f$ : 未来の部分列のサンプル平均
$\boldsymbol{\Sigma}_{pp}$ : 過去の部分列のサンプル分散共分散行列
$\boldsymbol{\Sigma}_{fp}$ : 過去の部分列と未来の部分列のサンプル共分散行列
$\boldsymbol{\Sigma}_{ff}$ : 未来の部分列のサンプル分散共分散行列

としているが，

$\boldsymbol{\mu}_p \equiv E[ \check{\boldsymbol{y}}(t) ]$

のように表しており，サンプル平均などの計算において $E[ \cdot]$ という演算子を使っていることに注意が必要である。

正準相関分析の適用

部分空間同定法では，過去の部分列と未来の部分列の相関係数が大きくなるように，線形変換の組 $\boldsymbol{a}_1, \boldsymbol{a}_2$ を求める。

「続々制御工学のこころ」を読む～第6章ビヨンド最小二乗法 ⑤補足：正準相関分析(理論編)～ - jiku logを参考にすると，相関係数の最大化問題は

$\begin{align} &\max_{\boldsymbol{a}_1, \boldsymbol{a}_2 } \boldsymbol{a}_1^T \boldsymbol{\Sigma}_{fp} \boldsymbol{a}_2 \\ \\ &\text{subject to } \boldsymbol{a}_1^T \boldsymbol{\Sigma}_{ff} \boldsymbol{a}_1 =1, \boldsymbol{a}_2^T \boldsymbol{\Sigma}_{pp} \boldsymbol{a}_2 =1 \\ \end{align}$

となる。
制約付き最大化問題であるため，ラグランジアンを導入し偏微分を実行すると，

$\begin{align} \boldsymbol{0} &= \boldsymbol{\Sigma}_{fp} \boldsymbol{a}_2 - \rho \boldsymbol{\Sigma}_{ff} \boldsymbol{a}_1 \\ \\ \boldsymbol{0} &= \boldsymbol{\Sigma}_{pf} \boldsymbol{a}_1 - \rho \boldsymbol{\Sigma}_{pp} \boldsymbol{a}_2 \\ \\ \Leftrightarrow & \begin{bmatrix} \boldsymbol{O} & \boldsymbol{\Sigma}_{fp} \\ \boldsymbol{\Sigma}_{pf} & \boldsymbol{O} \\ \end{bmatrix} \begin{bmatrix} \boldsymbol{a}_1 \\ \boldsymbol{a}_2 \\ \end{bmatrix} = \rho \begin{bmatrix} \boldsymbol{\Sigma}_{ff} & \boldsymbol{O} \\ \boldsymbol{O} & \boldsymbol{\Sigma}_{pp} \\ \end{bmatrix} \begin{bmatrix} \boldsymbol{a}_1 \\ \boldsymbol{a}_2 \\ \end{bmatrix} \\ \end{align}$

となり，本論文における(1)式のような，一般化固有値問題の式が得られる。この一般化固有値問題は，次式の特異値分解を実行することにより計算できる。

$\begin{align} \boldsymbol{W} \equiv \boldsymbol{\Sigma}_{ff}^{-1/2} \boldsymbol{\Sigma}_{fp} \boldsymbol{\Sigma}_{pp}^{-1/2} = \boldsymbol{U}_d \boldsymbol{S}_d \boldsymbol{V}^T \\ \\ \end{align}$

ここで，

$\begin{align} \boldsymbol{U}_d &= \begin{bmatrix} \tilde{\boldsymbol{a}}_1^1 & \tilde{\boldsymbol{a}}_1^2 & \cdots & \tilde{\boldsymbol{a}}_1^d \\ \end{bmatrix} \\ \\ \boldsymbol{V}_d &= \begin{bmatrix} \tilde{\boldsymbol{a}}_2^1 & \tilde{\boldsymbol{a}}_2^2 & \cdots & \tilde{\boldsymbol{a}}_2^d \\ \end{bmatrix} \\ \end{align}$

はそれぞれ左特異ベクトル・右特異ベクトルを並べた行列であり，
$\tilde{\boldsymbol{a}}_1^i = \boldsymbol{\Sigma}_{ff}^{1/2} \boldsymbol{a}_1^i$ ， $\tilde{\boldsymbol{a}}_2^i = \boldsymbol{\Sigma}_{pp}^{1/2} \boldsymbol{a}_2^i$ である。さらに，

$\begin{align} \boldsymbol{S}_d = \mathrm{diag} \begin{bmatrix} \sigma_1 & \sigma_2 & \cdots & \sigma_d \end{bmatrix} \\ \end{align}$

は，特異値を対角成分に持つ対角行列である。

$\boldsymbol{a}_1$ による未来の部分列の線形変換後の変数を $\boldsymbol{\alpha}$ ， $\boldsymbol{a}_2$ による過去の部分列の線形変換後の変数を $\boldsymbol{\beta}$ とすると， $\boldsymbol{\alpha}, \boldsymbol{\beta}$ はともに正準ベクトルと呼ばれ，

$\begin{align} \boldsymbol{\alpha} = \begin{bmatrix} \boldsymbol{a}_1^{1T} \\ \boldsymbol{a}_1^{2T} \\ \vdots \\ \boldsymbol{a}_1^{dT} \\ \end{bmatrix} \boldsymbol{y}_k (t) = \boldsymbol{U}_d^T \boldsymbol{\Sigma}_{ff}^{-1/2} \boldsymbol{y}_k (t) \\ \end{align}$

$\begin{align} \boldsymbol{\beta} = \begin{bmatrix} \boldsymbol{a}_2^{1T} \\ \boldsymbol{a}_2^{2T} \\ \vdots \\ \boldsymbol{a}_2^{dT} \\ \end{bmatrix} \check{\boldsymbol{y}}_k (t) = \boldsymbol{V}_d^T \boldsymbol{\Sigma}_{pp}^{-1/2} \check{\boldsymbol{y}}_k (t) \\ \end{align}$

となる。

状態ベクトルの定義

正準相関分析にもとづく部分空間同定法では，状態空間の基底は過去の部分列 $\check{\boldsymbol{y}}_k(t)$ と未来の部分列 $\boldsymbol{y}_k (t)$ の正準ベクトルで表される。このとき状態ベクトルは，以下のように定義される。

$\begin{align} \boldsymbol{x} (t) &\equiv \boldsymbol{S}_d^{1/2} \boldsymbol{\alpha}(t) \\ \\ & = \boldsymbol{S}_d^{1/2} \boldsymbol{U}_d^T \boldsymbol{\Sigma}_{ff}^{-1/2} \boldsymbol{y}_k (t) \in \mathbb{R}^d \\ \end{align}$

この式より，状態ベクトルは，「過去の部分列と未来の部分列の相関係数が最大になる」ように選ばれた正準ベクトルの各成分について，特異値によって重みづけしたベクトルであると解釈することができる。

なおこの定義より，状態ベクトル $\boldsymbol{x} (t)$ は以下の関係を満たす。

$\begin{align} E[\boldsymbol{x}(t) \boldsymbol{x}(t)^T ] &= \boldsymbol{S}_d^{1/2} E[\boldsymbol{\alpha}(t) \boldsymbol{\alpha}(t)^T ] \boldsymbol{S}_d^{1/2} \\ \\ &= \boldsymbol{S}_d^{1/2} \boldsymbol{U}_d^T \boldsymbol{\Sigma}_{ff}^{-1/2} E [ \boldsymbol{y}_k (t) \boldsymbol{y}_k (t)^T ] \boldsymbol{\Sigma}_{ff}^{-1/2} \boldsymbol{U}_d \boldsymbol{S}_d^{1/2} = \boldsymbol{S}_d \\ \end{align}$

未来の部分列と状態ベクトル・過去の部分列の関係

未来の部分列と状態ベクトルの関係

状態空間モデルの状態方程式と観測方程式を用いると，

$\begin{align} \boldsymbol{y}(1) \approx \boldsymbol{C} \boldsymbol{x}(1), \quad \boldsymbol{y}(2) \approx \boldsymbol{C} \boldsymbol{A} \boldsymbol{x}(1), \quad \boldsymbol{y}(3) \approx \boldsymbol{C} \boldsymbol{A}^2 \boldsymbol{x}(1), ... \\ \end{align}$

という関係が得られる。ただし $\approx$ は，「ノイズを除いて等しい」ということを表す。この関係式を用いると，

$\begin{align} \boldsymbol{y}_k(t) &= \begin{bmatrix} \boldsymbol{y}(t) \\ \boldsymbol{y}(t+1) \\ \vdots \\ \boldsymbol{y}(t+k-1) \\ \end{bmatrix} \approx \begin{bmatrix} \boldsymbol{C} \boldsymbol{x}(t) \\ \boldsymbol{C} \boldsymbol{A} \boldsymbol{x}(t) \\ \vdots \\ \boldsymbol{C} \boldsymbol{A}^{k-1} \boldsymbol{x}(t) \\ \end{bmatrix} = \begin{bmatrix} \boldsymbol{C} \\ \boldsymbol{C} \boldsymbol{A} \\ \vdots \\ \boldsymbol{C} \boldsymbol{A}^{k-1} \\ \end{bmatrix} \boldsymbol{x}(t) \equiv \boldsymbol{O}_k \boldsymbol{x}(t) \\ \end{align}$

となる。ただし $\boldsymbol{O}_k$ は拡大可観測行列である。ノイズも含めて表現すると，

$\begin{align} \boldsymbol{y}_k(t) = \boldsymbol{O}_k \boldsymbol{x}(t) + \boldsymbol{v}_k(t)\\ \end{align}$

となる。

なお，

$\begin{align} \boldsymbol{x} (t) = \boldsymbol{S}_d^{1/2} \boldsymbol{U}_d^T \boldsymbol{\Sigma}_{ff}^{-1/2} \boldsymbol{y}_k (t) \quad \Leftrightarrow \quad \boldsymbol{y}_k (t) = \boldsymbol{\Sigma}_{ff}^{1/2} \boldsymbol{U}_d \boldsymbol{S}_d^{-1/2} \boldsymbol{x} (t)\\ \end{align}$

であるので，

$\begin{align} \boldsymbol{O}_k = \boldsymbol{\Sigma}_{ff}^{1/2} \boldsymbol{U}_d \boldsymbol{S}_d^{-1/2} = \begin{bmatrix} \boldsymbol{C} \\ \boldsymbol{C} \boldsymbol{A} \\ \vdots \\ \boldsymbol{C} \boldsymbol{A}^{k-1} \\ \end{bmatrix} \end{align}$

となる。

未来の部分列と過去の部分列の関係

過去の部分列 $\check{\boldsymbol{y}}_k(t)$ と未来の部分列 $\boldsymbol{y}_k(t)$ の関係式を

$\begin{align} \boldsymbol{y}_k(t) = \boldsymbol{P} \check{\boldsymbol{y}}_k(t) + \boldsymbol{v}_k(t)\\ \end{align}$

と表す。過去の部分列 $\check{\boldsymbol{y}}_k(t)$ と未来の部分列 $\boldsymbol{y}_k(t)$ は，ともに観測されているので， $\boldsymbol{P}$ は最小二乗法で求めることができる。

$\begin{align} J(\boldsymbol{P}) &= E[ (\boldsymbol{y}_k(t) - \boldsymbol{P} \check{\boldsymbol{y}}_k(t) ) (\boldsymbol{y}_k(t) - \boldsymbol{P} \check{\boldsymbol{y}}_k(t) )^T ] \\ \\ &= \boldsymbol{\Sigma}_{ff} - \boldsymbol{\Sigma}_{fp} \boldsymbol{P}^T - \boldsymbol{P} \boldsymbol{\Sigma}_{pf} + \boldsymbol{P} \boldsymbol{\Sigma}_{pp} \boldsymbol{P}^T \\ \end{align}$

について， $\partial{J} / \partial{\boldsymbol{P}} = 0$ を解くと， $\hat{\boldsymbol{P}} = \boldsymbol{\Sigma}_{fp} \boldsymbol{\Sigma}_{pp}^{-1}$ なので，

$\begin{align} \boldsymbol{y}_k(t) = \boldsymbol{\Sigma}_{fp} \boldsymbol{\Sigma}_{pp}^{-1} \check{\boldsymbol{y}}_k(t) + \boldsymbol{v}_k(t)\\ \end{align}$

となる。

これらの関係から，状態ベクトル $\boldsymbol{x} (t)$ は，未来の部分列 $\boldsymbol{y}_k(t)$ を予測するために必要となる過去のデータを縮約したものである，と言える。

部分空間同定法の確率的解釈

ここからは，部分空間同定法に確率モデルを導入する。確率モデルを導入するには，

変数間のグラフィカルモデル
確率分布

を考える必要があるが，今回は確率分布として正規分布を考えているので，グラフィカルモデルを考える。

「正準相関分析に基づく部分空間同定法」におけるグラフィカルモデル

正準相関分析に基づく部分空間同定法では，正準相関分析の結果から状態ベクトル $\boldsymbol{x}(t)$ を導入した。この結果から，

状態ベクトル $\boldsymbol{x}(t)$ は，未来の部分列 $\boldsymbol{y}_k(t)$ を予測するために，過去の部分列 $\check{\boldsymbol{y}}_k(t)$ を縮約したベクトルである。

と解釈できる。この解釈に基づき，以下のようなグラフィカルモデルを考える。

確率モデル

このグラフィカルモデルから，以下の2つの条件付き分布を考える。

$\boldsymbol{y}_k(t)$ と $\boldsymbol{x}(t)$ の関係式 : $p(\boldsymbol{y}_k(t) \lvert \boldsymbol{x}(t) )$
$\check{\boldsymbol{y}}_k(t)$ と $\boldsymbol{x}(t)$ の関係式 : $p(\check{\boldsymbol{y}}_k(t) \lvert \boldsymbol{x}(t) )$

1つ目の $\boldsymbol{y}_k(t)$ と $\boldsymbol{x}(t)$ の関係式は，正準相関分析に基づく部分空間同定法によると

$\begin{align} \boldsymbol{y}_k(t) & = \boldsymbol{\Sigma}_{ff}^{1/2} \boldsymbol{U}_d \boldsymbol{S}_d^{-1/2} \boldsymbol{x} (t) + \boldsymbol{v}_k(t) \quad (= \boldsymbol{O}_k \boldsymbol{x}(t) + \boldsymbol{v}_k(t) ) \\ \\ \boldsymbol{x} (t) &\equiv \boldsymbol{S}_d^{1/2} \boldsymbol{\alpha}(t) \\ \end{align}$

となる。すなわち，未来の部分列(データ行列はハンケル行列)から得られる正準ベクトルを用いている。

2つ目の $\check{\boldsymbol{y}}_k(t)$ と $\boldsymbol{x}(t)$ の関係式は，1つ目の関係式とのアナロジーを考えると，過去の部分列(データ行列はテプリッツ行列)から得られる正準ベクトルを用いていればよい。

$\begin{align} \check{\boldsymbol{y}}_k(t) & = \boldsymbol{\Sigma}_{pp}^{1/2} \boldsymbol{V}_d \boldsymbol{S}_d^{-1/2} \boldsymbol{x} (t) + \check{\boldsymbol{v}}_k(t)\\ \\ \boldsymbol{x} (t) &\equiv \boldsymbol{S}_d^{1/2} \boldsymbol{\beta}(t) \\ \end{align}$

ただし $\check{\boldsymbol{v}}_k(t)$ は正規分布にしたがうノイズである。

これらの関係式から，以下の確率モデルが得られる。

$\begin{align} p(\boldsymbol{x}(t)) &= \mathcal{N}(\boldsymbol{0}, \boldsymbol{I}_d) \\ \\ p(\boldsymbol{y}_k(t) \lvert \boldsymbol{x}(t)) &= \mathcal{N}(\boldsymbol{F}_1 \boldsymbol{x}(t), \boldsymbol{\Psi}_1) = \mathcal{N} (\boldsymbol{\Sigma}_{ff}^{1/2} \boldsymbol{U}_d \boldsymbol{S}_d^{-1/2} \boldsymbol{x} (t), \boldsymbol{\Psi}_1) \\ \\ p(\check{\boldsymbol{y}}_k(t) \lvert \boldsymbol{x}(t)) &= \mathcal{N}(\boldsymbol{F}_2 \boldsymbol{x}(t), \boldsymbol{\Psi}_2) = \mathcal{N} (\boldsymbol{\Sigma}_{pp}^{1/2} \boldsymbol{V}_d \boldsymbol{S}_d^{-1/2} \boldsymbol{x} (t), \boldsymbol{\Psi}_2) \\ \end{align}$

ただし， $\boldsymbol{\Psi}_1, \boldsymbol{\Psi}_2$ はそれぞれ $\boldsymbol{v}_k(t), \check{\boldsymbol{v}}_k(t)$ の分散共分散行列である。

状態ベクトルの事後分布の推定

状態ベクトルから未来の部分列・過去の部分列を求める条件付き分布が得られたので，事後分布を用いればデータから状態ベクトルを推定することが可能になる。特に今回の場合，それぞれの変数が正規分布にしたがっているので，多変量正規分布の条件付き分布の公式が使える。この公式を用いると，過去の部分列が得られたときの状態ベクトルの事後分布の平均・分散共分散行列は

$\begin{align} E [ \boldsymbol{x}(t) \lvert \check{\boldsymbol{y}}_k(t) ] = \boldsymbol{S}_d^{1/2} \boldsymbol{V}_d^T \boldsymbol{\Sigma}_{pp}^{-1/2} \check{\boldsymbol{y}}_k(t), \quad V [ \boldsymbol{x}(t) \lvert \check{\boldsymbol{y}}_k(t) ] = \boldsymbol{I}_d - \boldsymbol{S}_d \\ \end{align}$

となる。

状態空間モデルとの関係

正準相関分析に基づく部分空間同定法は，未来の部分列と状態ベクトルの関係式を，拡大可観測行列を用いて表現した。過去の部分列と状態ベクトルの関係式についても，同様の表現を求めてみる。

正準相関分析を用いると，状態ベクトルとして2種類の正準ベクトルが得られる。これらの状態ベクトルを，以下のように定義する(なお"f"は「前向き」(forward)，"b"は「後ろ向き」(backward)を表す)。

$\begin{align} \boldsymbol{x}^f (t) \equiv \boldsymbol{S}_d^{1/2} \boldsymbol{\alpha}(t), \quad \boldsymbol{x}^b (t) \equiv \boldsymbol{S}_d^{1/2} \boldsymbol{\beta}(t) \\ \end{align}$

これらの状態ベクトルは，それぞれ以下の前向きマルコフモデル，後ろ向きマルコフモデルの状態ベクトルに相当する。

$\begin{align} \boldsymbol{x}^f (t+1) &= \boldsymbol{A} \boldsymbol{x}^f (t) + \boldsymbol{w}^f (t), \quad \boldsymbol{y}^f (t) = \boldsymbol{C} \boldsymbol{x}^f (t) + \boldsymbol{v}^f (t) \\ \\ \boldsymbol{x}^b (t-1) &= \boldsymbol{A}^T \boldsymbol{x}^b (t) + \boldsymbol{w}^b (t), \quad \boldsymbol{y}^b (t) = \bar{\boldsymbol{C}}^T \boldsymbol{x}^b (t) + \boldsymbol{v}^b (t) \\ \end{align}$

ただし $\bar{\boldsymbol{C}} = E[ \boldsymbol{x}^f (t+1) \boldsymbol{y}^f (t) ]$ である。

後ろ向きマルコフモデルは，Katayama "Subspace Methods for System Identification"の"4.8 Backward Markov Model" に紹介されている。後ろ向きマルコフモデルは，状態空間モデルの双対表現として導出される。詳細はリンク先を参照していただきたいが，導出のアウトラインは以下の通りである。

$\boldsymbol{\Pi} = E[ \boldsymbol{x}^f (t) \boldsymbol{x}^f (t) ^T ]$ とする。
システムノイズを $\boldsymbol{w}^b (t) \equiv \boldsymbol{\Pi}^{-1} \boldsymbol{x}^f (t) - \boldsymbol{A}^{T} \boldsymbol{\Pi}^{-1} \boldsymbol{x}^f (t+1)$ とする。
状態ベクトルを $\boldsymbol{x}^b (t) \equiv \boldsymbol{\Pi}^{-1} \boldsymbol{x}^f (t+1)$ とする。
状態ベクトルとシステムノイズ・観測ノイズの分散共分散行列を求めて，後ろ向きマルコフモデルも状態空間モデルのような性質を満たすことを確認する。

後ろ向きマルコフモデルを用いると，過去の部分列と状態ベクトルの関係式は以下のようになる。

$\begin{align} \check{\boldsymbol{y}}_k(t) = \boldsymbol{O}_k^b \boldsymbol{x}(t) + \check{\boldsymbol{v}}_k(t) , \quad \boldsymbol{O}_k^b = \begin{bmatrix} \bar{\boldsymbol{C}}^T \\ \bar{\boldsymbol{C}}^T \boldsymbol{A}^T \\ \vdots \\ \bar{\boldsymbol{C}}^T (\boldsymbol{A}^T)^{k-1} \\ \end{bmatrix} \\ \end{align}$

まとめと感想

今回は，上甲昌郎他著「局所線形アライメントによる非線形動的システムの学習法」(人工知能学会論文誌，2011)で紹介されていた部分空間同定法の確率的解釈についてまとめた。

単に正準相関分析を用いるだけでも，部分空間同定法の目的である状態ベクトルを導出することができる。しかしこれを確率モデルにすることで，ベイズ推定などの機械学習で良く用いられる手法が適用できるようになったり，本論文でも紹介されているような非線形化の拡張ができるようになったりする，といったメリットが理解できた。

また部分空間同定法は，拡大可観測行列など状態空間モデルともかかわりが深いが，後ろ向きマルコフモデルを導入することで，部分空間同定法の確率モデルとの関連性も理解ができた。

本論文はこの後，混合確率モデルを導入して非線形システムのシステム同定を行なっているのだが，今後はこの内容についても理解を深めていきたい。

本記事を最後まで読んでくださり，どうもありがとうございました。