カルマンフィルターは、観測できない確率変数(潜在変数)があるモデル(状態空間モデル)において、観測可能な確率変数に基づいて観測不可能な変数の状態を推定する手法です。カルマンフィルターは行列計算を用いて初等的に導出することも可能ですが、計算が結構複雑になる上に、逆行列の存在を仮定しなければならない部分もあったりと難点もあります。この記事では、ヒルベルト空間における直交射影の理論を用いてカルマンフィルターの漸化式を導出します。こちらのほうが理論がすっきりしており、数学の知識がある方には理解しやすいと思います。
問題設定
時刻k=0,1,2,⋯に対し、xkを観測できない状態量、ykを観測可能な観測値、ukを既知の入力ベクトルとし、それらが以下の条件を満たしているとします。 {xk+1=Akxk+Bkuk+ξkyk=Ckxk+Dkuk+ηk (k=1,2,⋯)
ここで、
- 状態変数xkおよびノイズξk はRn値確率変数、 観測ベクトルykおよび ノイズηk はRp値確率変数 、入力ベクトルuk∈Rm
- Ak∈Mn(R),Bk∈Mm,n(R),Ck∈Mn,p(R),Dk∈Mm,p(R)
- x0,ξ1,ξ2,⋯,η1,η2,⋯は独立な確率変数で、 E[ξk]=0, V[ξk]=Qk, E[ηk]=0, V[ηk]=Rk
- さらに、 Rkは正定値であるとする。
この設定の下で、既知のパラメータAk,Bk,Ck,Dk,Qk,Rkおよび観測値y1,y2,⋯,から、xkたちの最適な推定量を構成する方法を考えます。さらに、ここでは初期値x0も既知であるとします。
まずは、入力uk=0の場合を考えます。この場合、xk,ykの満たす方程式は、{xk+1=Akxk+ξkyk=Ckxk+ηk (k=1,2,⋯)となります。
直交射影による推定量の定義
ここでは、2つの確率変数ベクトルの「内積」もどきを定義し、直交射影を構成します。 以下では、行列Aに対しその転地行列をA′で表します。また、ほとんどいたるところ等しい確率変数は同一視することにします。
定義
Rd1値確率変数x、Rd2値確率変数yに対し、(x,y)および‖x‖を (x,y)=E[xy′], ‖x‖=(x,x)12 と定める。
ここで、(x,y)はスカラーではなく行列であることに注意します。 これに対しても、次に示すような内積に似た性質が成り立ちます。
命題1
Rd1値確率変数x,x1,x2、Rd2値確率変数y,y1,y2、A∈Mn,d1(R),B∈Mn,d2(R)に対し、次が成り立つ。
(1) (y,x)=(x,y)′
(2) (x1+x2,y)=(x1,y)+(x2,y), (x,y1+y2)=(x,y1)+(x,y2)
(3) (Ax,y)=A(x,y), (x,By)=(x,y)B′
(4) tr‖x‖=0 ⟺ x=0
(証明) (1) (y,x)=E[yx′]=E[xy′]′=(x,y)′
(2) (x1+x2,y)=E[(x1−x2)y′]=E[x1y′]+E[x2y′]=(x1,y)+(x2,y) であり、第2式も同様である。
(3) (Ax,y)=E[(Axy′]=AE[xy′]=A(x,y) であり、また(1)より、 (x,By)=(By,x)′={B(y,x)}′=(x,y)B′ である。
(4) tr‖x‖=trE[xx′]=trE[x′x]=E[|x|2]だから示すべきことを得る。 □
次に、Lk={A0x0+k∑j=1Ajyk;A0∈Mn,n(R),Aj∈Mn,p(R)}, L−1={0} と置きます。
上の命題より、Rn値確率変数全体はtr(⋅,⋅)を内積とするHilbert空間で、Lkはその閉部分空間となるので、直交射影の存在と一意性から次の命題が得られます。
命題2
Rn値確率変数xに対し、y∈Lkで tr‖y−x‖2=E[|y−x|2] を最小にするもの(xのLkへの直交射影)が一意的に存在する。
以下では、これをPkxと書く。特にP−1x=0である。
これを用いて、「観測値y1,⋯,yjに基づくxkの推定量」xk|jを次のように定義します。
定義
k,j=0,1,2,⋯に対し、ˆxk|jを ˆxk|j=Pjxk によって定める。
これは、x0およびy1,⋯,yjの線形結合で表されるxkの推定量のうち、平均2乗誤差を最小にするものということです。
入力がない場合のカルマンフィルタ
次に、上で構成した推定量xk|jを具体的に計算するための式を導きます。以下では、時刻が進むたびにy1,y2,⋯たちが順次観測されてゆき、ykが観測された時点で、その時の状態量xkを推定したい状況を考えます。つまり、xk|kを知りたいということです。
この時、時刻kにおいて新たに観測されたykと、 時刻k−1で求めたxk−1|k−1の値を使って、新たな推定量xk|kを計算するような漸化式が得られれば、計算のコストが省けてうれしいわけです。そのような漸化式を導くために、いくつかの補題を証明します。いずれも補題も、直交射影の性質から類推されるものです。
以下記述を簡単にするため、y0=x0と定めます。
補題1
ˆy∈LkおよびRn値確率変数xに対し、次は同値である。
(1) ˆy=Pkx
(2) j=0,1,⋯,kに対し、(ˆy−x,yj)=O
(証明) (1)⟹(2):背理法で示す。(1)を満たすˆyに対し、ある0≤j≤kが存在して(ˆy−x,yj)≠Oであるとする。この時、yj≠0であるから、tr‖yj‖2=E[|yj|2]>0である。そこで、 y=ˆy−(ˆy−x,yj)tr‖yj‖2yj と置くと、 tr‖y−x‖2=tr‖ˆy−x−(ˆy−x,yj)tr‖yj‖2yj‖2=tr‖ˆy−x‖2−tr((ˆy−x,yj)tr‖yj‖2yj,ˆy−x)−tr(ˆy−x,(ˆy−x,yj)tr‖yj‖2yj)+tr‖(ˆy−x,yj)tr‖yj‖2yj‖2=tr‖ˆy−x‖2−tr(ˆy−x,yj)(yj,ˆy−x)tr‖yj‖2−tr(ˆy−x,yj)(ˆy−x,yj)′tr‖yj‖2+tr(ˆy−x,yj)(ˆy−x,yj)′tr‖yj‖2=tr‖ˆy−x‖2−tr(ˆy−x,yj)(ˆy−x,yj)′tr‖yj‖2 であるが、(ˆy−x,yj)≠Oから上と同様にtr(ˆy−x,yj)(ˆy−x,yj)′>0となり、(1)の仮定に反する。よって、(2)が成り立つ。
(2)⟹(1):y∈Lkとすると、仮定より(ˆy−x,y−ˆy)=Oであるから、 tr‖y−x‖2=tr‖ˆy−x+y−ˆy‖2=tr‖ˆy−x‖2+tr‖y−ˆy‖2≥tr‖ˆy−x‖2 となり(1)が成り立つ。 □
上の補題により、次のような結果が得られます。
補題2
(1) 任意のA∈Mn,p(R)およびk=0,1,2,⋯に対し、 PkAx=APkx が成り立つ。
(2) 任意の0≤j≤kに対し、Pjξk=0が成り立つ。
(証明) (1) 任意のj=0,1,⋯,kに対し、 (APkx,yj)=A(Pkx,yj)=O だから上の補題により示すべきことを得る。
(2) (2) 任意のi=0,1,⋯,jに対し、yiはξi−1までの式で表されるから、ξkとは独立である。よって、(ξk,yi)=0となるから、上の補題により結論が従う。 □
次に、{yk}の「直交化」であるイノベーション過程を定義します。
定義
zk=yk−CkPk−1xk=yk−Ckˆxk|k−1 で定まる{zk}をイノベーション過程と呼ぶ。特に、P−1=Oにより、z0=y0=x0である。
次の補題により、これが{yk}の「直交化」になっていることがわかります。
補題3
イノベーション過程{zk}に対し、 (zi,zj)={Ci‖xi−Pi−1xi‖2C′i+Ri(i=j≥1)O(i≠j)Hk={k∑j=1Ajzk;A0∈Mn,n(R),Aj∈Mn,p(R)} である。
(証明) zk=Ckxk+ηk−CkPk−1xk=Ck(xk−Pk−1xk)+ηk であるから、i=jの時xi−Pi−1xiとηiの独立性により、 (zi,zi)=‖Ci(xi−Pi−1xi)‖2+‖ηi‖2=Ci‖xi−Pi−1xi‖2C′i+Ri となる。一方i<jとすると、ηjはηiおよびxi−Pi−1xiと独立であるから、 (zi,zj)=(Ci(xi−Pi−1xi)+ηi,Cj(xj−Pj−1xj))=(zi,Cj(xj−Pj−1xj)=(yi−CiPi−1xi,Cj(xj−Pj−1xj))=O となる。最後の等号は、補題1より(yi−CiPi−1xi,xj−Pj−1xj)=Oとなることを用いた。よって前半の主張が示された。後半の主張は帰納的に明らかである。 □
上の補題により、k≥1の時‖zk‖2は半正定値+正定値の形なので、正定値となることがわかります。そこで、ek=‖zk‖−1zkと定めると、{ek}k≥1は{yk}k≥1の「正規直交化」となり、次の補題が成り立ちます。
補題4
ek=‖zk‖−1zk (k≥1)とすると、Rn値確率変数xに対し、 Pjx=P0x+j∑i=1(x,ei)ei が成り立つ。
(証明) y=P0x+j∑i=1(xi,ei)eiと置くと、前の補題よりHjであr任意のei (1≤i≤j)に対し、 (x−y,ei)=(x,ej)−(y,ej)=(x,ei)−{(P0x,ej)+(x,ei)(ei,ei)}=0(x−y,x0)=(P0x−y,x0)=0 となる。これと前の補題により、任意のyi (0≤i≤j)に対し、(x−y,yi)=0となるから、補題1によりy=Pjxである。 □
以上の準備の下で、xk|kを計算するための漸化式が証明できます。
定理5
k≥1に対し、Gk=(xk,ek)‖zk‖−1と置くと、 ˆxk|k−1=Ak−1ˆxk−1|k−1ˆxk|k=ˆxk|k−1+Gk(yk−Ckˆxk|k−1) が成り立つ。
(証明) 補題2により、 ˆxk|k−1=Pk−1xk=Pk−1(Ak−1xk−1+ξk−1)=Ak−1Pk−1xk−1=Ak−1ˆxk−1|k−1 となるから最初の式が従う。
また、補題4により ˆxk|k=Pkxk=Pk−1xk+(xk,ek)ek=xk|k−1+(xk,ek)‖zk‖−1zk=xk|k−1+(xk,ek)‖zk‖−1(yk−Ckˆxk|k−1) となり2式目も従う。
上で登場したGkはカルマンゲインと呼ばれます。これを計算するための漸化式を求めます。
定理6
Pk,j=‖xk−ˆxk|j‖2とすると、 Pk,k−1=Ak−1Pk−1,k−1A′k−1+Qk−1Pk,k=(I−GkCk)Pk,k−1Gk=Pk,k−1C′k(CkPk,k−1C′k+Rk)−1 が成り立つ。
(証明) 前の定理を用いると、(xk−1−ˆxk−1|k−1)とξk−1の独立性により、 Pk,k−1=‖xk−ˆxk|k−1‖2=‖A(xk−1−ˆxk−1|k−1)+ξk−1‖2=‖A(xk−1−ˆxk−1|k−1)‖2+‖ξk−1‖2=Ak−1Pk−1,k−1A′k−1+Qk−1 となり最初の式が得られる。
また、ηkとxkの独立性および補題7から(xk,ηk)=0,(xk|k−1,xk−xk|k−1)=0となることに注意すると、 Gk=(xk,ek)‖zk‖−1=(xk,zk)‖zk‖−2=(xk,Ck(xk−xk|k−1)+ηk)‖zk‖−2=(xk,Ck(xk−xk|k−1))‖zk‖−2=(xk−xk|k−1,Ck(xk−xk|k−1))‖zk‖−2=‖xk−xk|k−1‖2C′k(zk,zk)−1=Pk,k−1C′k(CkPk,k−1C′k+Rk)−1 となり、3番目の式を得る。
さらにこれより、 Pk,k=‖xk−ˆxk|k‖2=‖xk−ˆxk|k−1−Gk(yk−Ckˆxk|k−1)‖2=‖xk−ˆxk|k−1−Gk(Ckxk−Ckˆxk|k−1+ηk)‖2=‖(I−GkCk)(xk−ˆxk|k−1)‖2+‖Gkηk‖2=(I−GkCk)Pk,k−1(I−GkCk)′+GkRkG′k=(I−GkCk)Pk,k−1−Pk,k−1C′kG′k+GkCkPk,k−1C′kG′k+GkRkG′k=(I−GkCk)Pk,k−1−Gk(CkPk,k−1C′k+Rk)G′k+GkCkPk,k−1C′kG′k+GkRkG′k=(I−GkCk)Pk,k−1 となって、2番目の式を得る。 □
以上により、xk|kを計算するための漸化式が得られました。定義により初期推定量x0|0はx0|0=P0x0|0=x0となることも合わせると、次のような結果になります。
カルマンフィルターによる状態推定(入力uk=0の場合)
Pk,k−1=Ak−1Pk−1,k−1A′k−1+Qk−1Pk,k=(I−GkCk)Pk,k−1Gk=Pk,k−1C′k(CkPk,k−1C′k+Rk)−1ˆx0|0=x0ˆxk|k=ˆxk|k−1+Gk(yk−Ckˆxk|k−1)ˆxk|k−1=Ak−1ˆxk−1|k−1
入力がある場合
次にukが一般の場合を考えます。この時は{x(1)0=x0,x(1)k+1=Akx(1)k+ξkx(1)1=0,x(2)k+1=Akx(2)k+Ckuk, {y(1)k=Bkx(1)k+ηky(2)k=Bkx(2)k+Dkukと定めると、xk=x(1)k+x(2)k, yk=y(1)k+y(2)kとなります。今x(2)kについては直接計算することができますので、x(1)kの値を推定できれば十分です。そこで、前節でxk→x(1)k,yk→y(1)kとして得られる推定量をˆx(1)k|jとし、ˆxk|j=ˆx(1)k|j+x(2)kと定めると、定理5のGkに対し、ˆxk|k=ˆx(1)k|k+x(2)k=ˆx(1)k|k−1+Gk(y(1)k−Ckˆx(1)k|k−1)+x(2)k=ˆxk|k−1+Gk{(yk−y(2)k)−Ck(ˆxk|k−1−ˆx(2)k|k−1)}=ˆxk|k−1+Gk(yk−Ckˆxk|k−1−Dkuk)ˆxk|k−1=ˆx(1)k−1|k+x(2)k=Ak−1ˆx(1)k−1|k−1+x(2)k=Ak−1(ˆxk−1|k−1−x(2)k−1)+Ak−1x(2)k−1+Ck−1uk−1=Ak−1ˆxk−1|k−1+Ck−1uk−1となります。以上をまとめると、次の結果になります。
カルマンフィルターによる状態推定(入力\bm u_kがある場合)
\begin{align*} &P_{k,k-1}=A_{k-1}P_{k-1,k-1}A_{k-1}’+Q_{k-1}\\ &P_{k,k}=(I-G_kC_k)P_{k,k-1}\\ &G_k=P_{k,k-1}C_k'(C_kP_{k,k-1}C_k’+R_k)^{-1}\\ &\hat{\bm x}_{0|0}=\bm x_0\\ &\hat{\bm{x}}_{k|k}=\hat{\bm{x}}_{k|k-1}+G_k(\bm{y}_k-C_k\hat{\bm{x}}_{k|k-1}-D_k\bm u_k)\\ &\hat{\bm{x}}_{k|k-1}=A_{k-1}\hat{\bm{x}}_{k-1|k-1}+C_{k-1}\bm u_{k-1}\end{align*}
P_{k,j}やG_kは入力があってもなくても同じ値になります。さらに、これらは観測データ\bm y_kがなくても事前に計算することができます。そのため、G_kをあらかじめ計算しておき、\bm y_kが観測された時点で\bm x_kを計算することにより、計算リソースを節約することができます。