はじめに
カルマンフィルタを解説する記事はたくさんあります。 詳しい理論や、細かい実装を知りたい場合は下記の記事などを参考にすると良いでしょう。
今回はTFPのdistributionsモジュールの中にある、比較的高レベルなAPIであるLinearGaussianStateSpaceModelというものを使い方の備忘録と、カルマンフィルタの意義の軽い説明です。特に状態観測器としての側面を理解することに重点を置きます。
カルマンフィルタは時系列モデルを扱う上で、学びやすく、そして学べることの多い素晴らしい題材です。
カルマンフィルタの意義
カルマンフィルタは結構多くの解説で「ノイズ除去に使える」という話がありますが、実はそれはカルマンフィルタの半分の魅力しか伝えられていません。
例えば下記の状態方程式
$$ \begin{align} x[t] &= Ax[t-1] \\ y[t] &= Cx[t] \end{align} $$
があったとしましょう。もしも $y[t]$ を観測できた場合には、$x[t]$ がどういう値であったのかを $x[t] = C^{-1}y[t]$ と計算することで把握できます。では下記の場合はどうでしょうか。
$$ \begin{align} x[t] &= Ax[t-1] \\ y[t] &= Cx[t] + v \end{align} $$
ここで $v$ とは観測時に生ずるノイズであり $v \sim \mathcal N({\bf 0}, \Sigma)$ という多変量正規分布から生起しているものとします。すると、先ほどのように $y[t]$ が観測できていたとしても、単純な式変形では
$$ x[t] = C ^ {-1} (y[t] - v) $$
という形式になってしまいます。ところがそもそも観測値 $y[t]$ に関して、ノイズ $v$ がどのようなものであったのかを知らないのですから、上記のように引き算を実行することはできずに終わってしまいます。
そこで我々は $x[t] = Ax[t-1]$ なる関係式を知っていることを思い出します。実は1個前の状態 $x[t-1]$ が分かっているならば、$x[t]$ を算出することはできてしまうのです。
しかしまてまて、1個前の$x[t-1]$ はどうやって知ればいいのだろうか？ そもそも1個前の観測 $y[t-1]$ が得られたところで同様の問題で $x[t-1]$ を算出できないのですから困ったものです。
……いや、困っていません。初期値 $x[0]$ さえ把握すれば $x[t] = A ^ {t}x[0]$ ではありませんか。解決解決。
そうだとすると、下記の場合はどうですか。
$$ \begin{align} x[t] &= Ax[t-1] + w\\ y[t] &= Cx[t] + v \end{align} $$
ここで $w$ も多変量正規分布に従うノイズです。さて、 $x[0]$ が分かれば $x[t]$ が分かるという理屈が通用しなくなりました。我々は観測できる $y[t]$ だけから内部の状態 $x[t]$ を知りたいというのです。そして行列 $C$ によっては、$y[t]$ の時系列と $x[t]$ の時系列が似ているという保証もありません。
例えば下記の図はどうでしょう。青と緑が観測時系列で、薄い赤と紫が対応する状態です。我々が実際に観測できるのは青と緑の時系列だけであることに注意してください。
図を見てみると、2つの観測データはいずれも似たような値を取っており、正の相関がかなり強そうです。ところが、その内部に隠れている状態の時系列は、負の相関を持っていることが見えます。
これは、状態遷移行列 $A$ が2つの状態変数が負の相関を持つような変換になっているのに対して、その符号を打ち消してしまうような観測行列 $C$ があるために観測された時には同じような値に見えているのです。
カルマンフィルタは状態方程式を "既知" とした場合に、観測 $y[t]$ から $x[t]$ を算出する方法（$C$ を使う方法）と、 $x[t-1]$ から $x[t]$ を算出する方法 （$A$を使う方法）を上手に織り交ぜた状態観測（推定）器だと言えます。制御の世界では確率的状態観測器として知られています（織り交ぜ具合を、ノイズを見積もりながら決めていくのである。そしてその織り交ぜ具合がカルマンゲインと呼ばれる係数で現れるのだ）。
TFPでのカルマンフィルタ
TFPでのカルマンフィルタは tfp.distributions.LinearGaussianStateSpaceModel クラスを利用することで実施できます。以下では実装のポイントを織り交ぜながら見ていきましょう。
モジュール
import numpy as np from matplotlib import pylab as plt import tensorflow.compat.v2 as tf import tensorflow_probability as tfp from tensorflow_probability import distributions as tfd from tensorflow import linalg as tfl plt.style.use("seaborn") tf.enable_v2_behavior()
まずは上記を使うのでインポートしておきます。普段あまり使わないものかもしれませんが、tf.linalg というモジュールが必要になります。ここでは tfl としておきましょう（TFP自体がもともと、TF本体にあったdistributions モジュールと linalg モジュールを取り出して使っているライブラリである）。
データの生成
さて、慣れている numpy でデータを生成してみましょう。
A = np.array([[1.2, 0.4], [-0.3, 0.5]]) C = np.array([[1.2, -0.2], [0.8, -0.6]]) def sample_data(): init_x = np.array([0., 1.0]) trans_cov = np.diag([0.05, 0.05]) obs_cov = np.diag([0.4, 0.4]) state_data = [] obs_data = [] x = init_x for _ in range(100): x = np.dot(x, A.T) + np.random.multivariate_normal([0, 0], trans_cov.T) y = np.dot(x, C.T) + np.random.multivariate_normal([0, 0], obs_cov.T) state_data.append(x) obs_data.append(y) return np.stack(state_data), np.stack(obs_data)
によってデータを作ることにします。ここで注意点は、（あとでも先でも良いが）基本的にデータの shape が (time_step, dims) になるようにしなければならないところです。数学では基本的に縦ベクトルとして記述されることが多いので、行列とベクトルの軸に注意してください。
ここでは 行列 $A$ と $C$ を数学の形式で定式化しておき、実際にプログラムで書くときにデータを横ベクトルとして扱うために転置することとします。このことは意識して置かなければなりません。なぜなら、TFPの高レベルAPIを使う時も、内部では横ベクトルとして計算しているが、モデルを組むときにはデータが縦ベクトルである前提で行列を渡さなければならないからです。
これによって生成されるデータが
state_data, obs_data = sample_data() plt.plot(obs_data) plt.plot(state_data, alpha=0.5) plt.legend(["obs1", "obs2", "state1", "state2"])
であったのでした。
TFPで線形状態空間モデルを作る
tfd.LinearGaussianStateSpaceModel を利用すれば簡単に作れます。
$$ \begin{align} x[t] &= Ax[t-1] + w\\ y[t] &= Cx[t] + v \end{align} $$
というモデルにおける、 $A$ と $C$ そして、ノイズ $w, v$ が従う正規分布の分散共分散行列を与えれば良いのです（加えて初期状態 $x[0]$ に関する事前分布も必要。定数を与えたい場合は平均に与えたい定数、分散を異様に小さくしとけばいいでしょう）
model_kf = tfd.LinearGaussianStateSpaceModel( num_timesteps=100, transition_matrix=tfl.LinearOperatorFullMatrix(A), transition_noise=tfd.MultivariateNormalFullCovariance( covariance_matrix=np.diag([0.05, 0.05]) ), observation_matrix=tfl.LinearOperatorFullMatrix(C), observation_noise=tfd.MultivariateNormalFullCovariance( covariance_matrix=np.diag([0.4, 0.4]) ), initial_state_prior=tfd.MultivariateNormalDiag( scale_diag=1 * tf.ones(2, dtype=tf.float64)) )
ここで注意点があります。行列を与えれば良いのだ！ということで、例えば transition_matrix=Aで渡したりしたくなってしまうのですが、これは後にエラーが出ます。numpy.array 及び、tf.eager_tensorが持っていないメソッドや属性を内部で使ってしまうため、tfd あるいは tfl のモジュールで表される行列を利用する必要があります。
また、既に述べた通り、ここでは $A, C$ また共分散行列も状態や観測が縦ベクトルである場合の形式で与えなければなりません。転置されたものを与えてしまうと違うモデルになってしまうので注意が必要です。
また、$A, C$ は自分たちが設計した機械製品であったりプラントであったり、何らかの同定可能な部分であるため既知とするのですが、遷移ノイズ transition_noise と観測ノイズobservation_noise は既知ではないケースも多いため、ここは値を試行錯誤的に決める場合も多いです。
例えば遷移ノイズを非常に小さく取るというのは状態の遷移が既知の $A$ によってのみほとんど支配されているはずだと考えていることに相当します。ある意味、設計したプラントが非常に正確であると考えていることになるでしょう。一方で観測ノイズを非常に小さく取るというのは、観測が物理的測定 $C$ によって正確に行われていると考えている場合に相当します。
カルマンフィルタの実行
作ったモデルで下記を唱えればフィルタリング後の時系列が得られます。
_, filt_state, _, _, _, _, _ = model_kf.forward_filter(
tf.convert_to_tensor(obs_data, dtype=tf.float64)
)
戻り値は7つあるので注意しましょう。フィルタリングされた時系列データは2番目に格納されています。3番目は共分散行列が入っているので、必要であればこれも取り出しましょう（例えば状態推定のばらつきを見たい場合）。ほかはカルマンフィルタのアルゴリズムの中で生ずる中間的なデータ達です。
plt.plot(obs_data, alpha=0.3) plt.plot(state_data) plt.plot(filt_state) plt.legend(["obs1", "obs2", "state1", "state2", "filt_state1", "filt_state2"])
見事に推定ができています。この推定は薄い青と緑の時系列だけからなされていることは驚くべきことでしょう。
…、いや驚くべきことでしょうか、思えば $A, C$ は既知なのですから、ノイズが小さいならば冒頭で述べた逆行列計算だけでなんとかなってしまうのでは…？と思えてきます（とはいえ、ノイズ除去の効果を、遷移ノイズと観測ノイズの設計で実際に実施できるのですから意義はある）。
追加実験
追加実験１：状態と観測の次元が異なるケース
さて、上記まででカルマンフィルタの基本的な用途を見ました。 次は、表題の通り、状態と観測の次元が異なるケースを見ます。特に、状態に対して観測が少ないケースです。 例えば以下のようなモデルです。
A = np.array([[1.2, 0.4, -0.1], [-0.3, 0.5, 0.1], [0.3, -0.2, 1.0]]) C = np.array([[1.2, -0.2, 0.2], [0.8, -0.6, 0.3]]) def sample_data(): init_x = np.array([0., 1.0, 1.0]) trans_cov = np.diag([0.05, 0.05, 0.05]) obs_cov = np.diag([0.4, 0.4]) state_data = [] obs_data = [] x = init_x for _ in range(100): x = np.dot(x, A.T) + np.random.multivariate_normal([0, 0, 0], trans_cov.T) y = np.dot(x, C.T) + np.random.multivariate_normal([0, 0], obs_cov.T) state_data.append(x) obs_data.append(y) return np.stack(state_data), np.stack(obs_data)
状態は３次元にもかかわらず、その３次元のデータを２次元でしか観測できない場合です（ $C$ が２行３列の行列である）。カルマンフィルタはこういった場合にでもちゃんと機能することを見ます。このとき観測されるデータと、状態をそれぞれ見ておきましょう（薄いほうが状態の時系列である）。カルマンフィルタはこの２本の時系列から、薄い３本の時系列を推定することができます。
ここでモデルを下記のように書きます。状態共分散と初期状態の事前分布を３次元にしておくことを忘れないようにしましょう。
model_kf = tfd.LinearGaussianStateSpaceModel( num_timesteps=100, transition_matrix=tfl.LinearOperatorFullMatrix(A), transition_noise=tfd.MultivariateNormalFullCovariance( covariance_matrix=np.diag([0.05, 0.05, 0.05]) ), observation_matrix=tfl.LinearOperatorFullMatrix(C), observation_noise=tfd.MultivariateNormalFullCovariance( covariance_matrix=np.diag([0.4, 0.4]) ), initial_state_prior=tfd.MultivariateNormalDiag( scale_diag=1 * tf.ones(3, dtype=tf.float64)) )
このモデルで同じように推定をしてみると
_, filt_state, _, _, _, _, _ = model_kf.forward_filter(
tf.convert_to_tensor(obs_data, dtype=tf.float64)
)
plt.figure(figsize=(12,8), dpi=100)
plt.plot(state_data)
plt.plot(filt_state)
plt.legend([
"state1", "state2", "state3",
"filt_state1", "filt_state2", "filt_state3"])
と得られるのです。とても先ほどの２本の観測時系列から得られたとは思えません。これが状態観測器としてのカルマンフィルタの重要な役割なのです。
よもや次元が違うために、ノイズが非常に小さければ、逆行列をたどるだけで良いのでは…？という発言は出ないでしょう。次元が欠落していても、３次元を復元できているのです（とは言いつつも、最小二乗法がそうであるように擬似逆行列の存在を知っていれば不思議ではない）。
ところが、カルマンフィルタを使っても、全く情報を復元できないケースもあるのです。
追加実験2: 不可観測系
なぜ先ほどのように、２次元から３次元を復元できたのでしょうか。それは $C$ で３次元から２次元に落ちたとしても、情報がごちゃまぜになって解像度が悪くなっただけであり、実際のところ情報が消えたわけではないからです。
例えば観測行列 $C$ が単位行列とした場合には全ての状態を単独で測定できるということに相当します。一方で 単位行列 $C$ からある1つの行を取り去ってしまった場合にはどうでしょうか。その成分は全く観測されないということになります。すなわちその成分に関する情報は、観測値からは一切得られないということです。
例えば下記のような観測行列 $C$ は状態 $x[t]$ に対して $y[t] = Cx[t]$ と作用し、2次元目の情報を完全に欠落してしまいます。
A = np.array([[1.2, 0.4, 0], [-0.3, 0.5, 0], [0.3, -0.2, 1.0]]) C = np.array([[1., 0., 0.], [0., 0., 1.]]) def sample_data(): init_x = np.array([0., 1.0, 1.0]) trans_cov = np.diag([0.05, 0.05, 0.05]) obs_cov = np.diag([0.4, 0.4]) state_data = [] obs_data = [] x = init_x for _ in range(100): x = np.dot(x, A.T) + np.random.multivariate_normal([0, 0, 0], trans_cov.T) y = np.dot(x, C.T) + np.random.multivariate_normal([0, 0], obs_cov.T) state_data.append(x) obs_data.append(y) return np.stack(state_data), np.stack(obs_data)
しかし、それでもカルマンフィルタは下記のように動作します。
カルマンフィルタが $x[t-1]$ から $A$ を使って $x[t]$ を推定するような動作も含んでいることを話しました。なので、 $A$ が状態遷移の中で色々な情報を混ぜあわせてくれるお陰で、ある瞬間は観測の立場で2番目の成分を完全に欠落しているとはいえど、常に立ち代わりで状態を見続けることができているのです。
ところが、この状態遷移の方でさえ、下記のように2番目の成分が他の成分から独立してしまった場合は打つ手がありません。
A = np.array([[0.9, 0., 0.4], [0., 1.0, 0.], [-0.3, 0., 1.0]]) C = np.array([[1., 0., 0.], [0., 0., 1.]]) def sample_data(): init_x = np.array([0., 1.0, 1.0]) trans_cov = np.diag([0.05, 0.05, 0.05]) obs_cov = np.diag([0.4, 0.4]) state_data = [] obs_data = [] x = init_x for _ in range(100): x = np.dot(x, A.T) + np.random.multivariate_normal([0, 0, 0], trans_cov.T) y = np.dot(x, C.T) + np.random.multivariate_normal([0, 0], obs_cov.T) state_data.append(x) obs_data.append(y) return np.stack(state_data), np.stack(obs_data)
これは状態遷移行列が第1成分と第3成分が値を混ぜあわしつつも、第2成分は自分自身の値しか寄与しません。したがって、その後観測行列では第２成分を知るすべが無いため、情報が一切なくなってしまいます。そのためカルマンフィルタは下記のように振る舞います。
このような状態方程式は、不可観測システムであると言います。 不可観測系に対してはどんなオブザーバを立てても推定はできません（当たり前だ。これはなんかシステムが時間発展してるらしいが、データを一度も見たことがありませんと言っているのだから）。
線形システムの場合は、システムが可観測であるか不可観測であるかを式で求めることができます。
最後に
さて、今回はTFPを使いながらカルマンフィルターの振る舞いを見ました。 もともとは単にTFPの備忘録でしたが、ふと制御理論が懐かしくなって、簡単にですが状態推定器としてのカルマンフィルタを概観してみたところです。
ここから発展していく方向性だけいくつか述べておきます。
まず線形カルマンフィルタの線形性を廃するのであれば、単に非線形を時間局所的に線形で近似できるとみなした拡張カルマンフィルタが使えます。この場合 $A, C$ の代わりに、それぞれヤコビ行列が用いられるだけで大筋は変わりません。
また、カルマンフィルタは正規分布の誤差を仮定しますが、この仮定を廃してしまったものがUncented Kalman Filterと呼ばれる手法です。また、Particle filterもこの種の仲間でしょう。
一方で、離散変数（カテゴリ変数）の時間発展を扱うような場合には隠れマルコフモデルというものが利用されます。ここでも正規分布の仮定は取り払われています。
さて、今回は状態方程式は既知の元、観測データだけから状態を推定するという試みをしました。一方で機械学習の立場で見れば、これは出来上がったモデルを使って潜在変数の事後分布を単に計算しているだけに過ぎません。時系列データに対して機械学習を用いる場合、むしろ観測データなどから、そのデータが従っている状態方程式を知りたい（すなわち $A$ や $C$ を推定したい）という問題に出くわします。
カルマンフィルタの場合はEMアルゴリズムと呼ばれる手法によって、潜在変数と状態方程式のパラメータを推定するということが一応定式化されます（あまり使われるところは見たことがないが）。
ちょうど隠れマルコフモデルではEMアルゴリズムの一種であるバウム・ウェルチのアルゴリズムがパラメータの推定に用いられます。また、学習結果を用いて潜在変数の時系列を観測データから推定する場合にはビタビアルゴリズム（動的計画法の一種）が用いられます。
またカルマンフィルタを状態観測器として見た場合には、その計算に繰り返し処理（例えば確率分布からのサンプリングなど）を必要としないことが大きく利きます。制御ではサンプリング周波数が100kHz以上あるような場合が多く、要0するに推定が10ms以内に終了しなければならないことが要求されます。
適切な制御設計をするためには、まずは線形システムの事自体を深く知る必要があるでしょう。（僕は全然しらんけど）
