逐次最小二乗法

最小二乗法の逐次解法 †

最小二乗法で代表されるシステムのパラメータ推定手法について、紹介する。特に、時系列データーのように、入出力データーが逐次的に得られる場合のパラメーター推定の方法を示す。

線形離散システムのパラメーター推定 †

線形システムの次数n を既知と仮定すると，その入出力関係は次式のような差分方程式で表される．

xt = Σaixt-i +Σbiut-i i=1,n　　　　　　eq1
yt = xt + vt                          eq2
但し、utは入力、ytは出力
      Vtは平均値ゼロの離散値確率雑音。入力と無相関　E[utvt] = 0

システムの入出力の観測値ut とyt からシステムのパラメータ(ai,bi)i=1～n　を推定することを離散システムの同定問題という。

↑

ベクトル形式で表現 †

2式を1式に代入して整理する。

yt = Σaiyt-i +Σbiut-i + rt
但し　rt =Σaivt (a0 = 1)　　　　　　　　eq3

上式は次のようなベクトル形式で表現できる。

yt = z't θ + rt
z't = [-yt-1,....,-yt-n, ut-1,...., ut-n]
θ' = [a1,...., an, b1,...., bn]

↑

誤差最小化問題と推定値 †

パラメータを推定するため，次のような二乗誤差評価を最小にすることを考える。データ数はNとしよう。

J(θ) =Σρ^(N-t)rt^2 =Σρ^(N-t)(yt - z't θ)^2
ただし，0 < ρ <= 1 は忘却係数である．

忘却の効果を小さくすれば，ゆっくり変動する時変パラメータの推定にも適用できる。

J を最小にする必要条件はJ(θ)の一階偏微分ベクトル = 0 である．すなわち，

∂J/∂θ= -2Σρ^(N-t)ztyt + 2Σρ^(N-t)ztz't θ = 0

これより、θの最小二乗推定値は

θ*=[Σztz't]^(-1)・Σρ^(N-t)ztyt

上式が成り立つためには，上記の逆行列が存在する必要があるが、一般に入力信号が十分多くの周波数成分を含んでいれば成立する。

↑

逐次推定の考え方と方法 †

データ数がNの場合のθの最小二乗推定値をθ*(N)と置く。パラメータ推定のための評価規範を、下記のように置いてみよう。

J(θ;N) =(1/N)Σρ^(N-t)rt^2 =(1/N)Σρ^(N-t)(yt - z't θ)^2

これは、次のように書きかえられる。

J(θ;N) =C(N)-2θ'b(N)+θ'R(N)θ
ただし　b(N)=(1/N)Σρ^(N-t)ytzt
       R(N)=(1/N)Σρ^(N-t)ztz't
       C(N)=(1/N)Σρ^(N-t)(yt)^2

ここでJ(θ;N) をθ に関して微分して0とおく。すると、正規方程式(normal equation) と呼ばれるθに関する連立一次方程式が得られる。この解が推定値である。

R(n)θ*(N)=b(N)　：正規方程式

行列R(N) が正定値であるかどうかが問題になるが、正定値である場合には、逆行列を用いてパラメータ推定値を求めることができ、この同定法は一括処理最小二乗法(batchleastsquares method)、あるいはオフライン最小二乗法(off-line least-squares method) と呼ばれる。

前節で説明した最小二乗法をもとに、再帰的最小二乗法(RLS アルゴリズム)を導出してみよう。

正規方程式は次のように、分解表現できる。

θ*(N)=R(n)^(-1)・b(N)
     =[Σρ^(N-t)ztz't]^(-1)・(Σρ^(N-t)ytzt)　　　　　　５式

まず、行列P(N) を

P(N)=[Σρ^(N-t)ztz't]^(-1)　　　　　　　　　　　　　　　　６式

とおき、これを共分散行列と呼ぶ。すると、

P(N)^(-1) =Σρ^(N-t)ztz't
          =P(N-1)^(-1)+zNz'N　　　　　　　　　　　　　　７式

が得られる。同様にして

b(N)=Σρ^(N-t)ytz't + yNz'N　　　　　　　　　　　　　　　　8式

６から８式を５式に代入して

θ*(N)=P(N)・(Σρ^(N-t)ytzt)  
     =P(N)・[P(N-1)^(-1)θ*(N-1)+yNzN]
     =P(N)・{[P(N)^(-1)-znz'n]θ*(N-1)+yNzN]
     =θ*(N-1)-P(N)zn[z'Nθ*(N-1)-yN]
     =θ*(N-1)+P(N)zn[yN-z'Nθ*(N-1)]                    ９式

９式と７式が、再帰的最小二乗法(RLS アルゴリズム)である。しかし、P(N)の逆行列をオンラインで計算することは困難である。そこで、逆行列補題(matrix inversion lemma) を用いて、(1.17)式をオンライン計算が可能な形式に変形しよう。

逆行列補題(matrix inversion lemma)

逆行列補題:ある正則行列A に対して次式が成立する。
(A+ BC)^(−1) = A^(−1) −A^(−1)B[I +CA^(−1)B]^(−1)CA^(−1)

7式に逆行列補題を適用すると、次式が得られる

P(N)=[P(N-1)^(-1)+zNz'N]^(-1)
    =P(N-1)-P(N-1)(zNz'N)P(N-1)/{1+z'NP(N-1)zN}     10式

9式も右辺第二項に含まれるP(N)zNを上式を代入することで、さらに簡単に修正できる。

 P(N)zN={P(N-1)-P(N-1)(zNz'N)P(N-1)/[1+z'NP(N-1)zN]}zN
       =  P(N-1)zN/[1+z'NP(N-1)zN]　　　　　　　　　　１１式

上のP(N)zNを９式に代入して整理する。

θ*(N)=θ*(N-1)+{P(N-1)zN/[1+z'NP(N-1)zN]}・[yN-z'Nθ*(N-1)]　１２式

この１２式と１０式が、最終的な逐次最小二乗法の式である。

↑

推定値の漸近的性質 †

データ数が多くなるにつれて，確率極限

N-->無限大の時、plim θ*=　θ

が成立するとき，θ* は一致推定量であるという。

証明： ρ = 1の場合、すなわち忘却しない場合で考えよう。

plim θ*=　｛plim (1/N)Σztz't｝^(-1) plim (1/N)Σzt[z't θ + rt]
       = θ +｛plim (1/N)Σztz't｝^(-1)plim (1/N)Σzt rt

が得られるので、一致推定量であるためには，ベクトルzt と式誤差rt は無相関でなければならない。

plim (1/N)Σzt rt = 0

が成り立つ必要がある。

一般に上式が成立することが期待できないので，観測雑音のレベルが高い場合，最小二乗推定は一致推定量でなく，偏り（バイアス）を持つ

↑

補助変数法によるパラメータの一致推定方法 †

ベクトルzt と同サイズの補助変数ベクトルmt を導入し，次のように補助変数法による推定値 θ** を得る。

θ**=[Σmtz't]^(-1)・Σρ^(N-t)mtyt

前節と同様に式変換することで、補助変数ベクトルmt が次のような条件を満足すれば，パラメータ推定値は一致推定量であることが言える。

１．下記の逆行列が存在する。
```
[Σmtz't]^(-1)
```
2.補助変数ベクトルmt と式誤差rt は無相関である．すなわち，
```
plim (1/N)Σmt rt = 0
```

この２つの条件を満たす補助変数mtの選び方は、いろいろある。

遅延された入力を利用する方法
```
mt = [ut-d,...,ut-d-n, ut-1,...., ut-n]
```
遅延された出力を利用する方法
```
mt = [yt-1-d,....., yt-n-d, ut-1,...... ut-n]
```
- ただし，d の選び方は，明確ではない.

推定された出力を利用する方法
```
mt = [x*t-1,....,x*t-n, ut-1,....., ut-n]
ただし、
x*t = Σa*ix*t-i +Σb*iut-i i=1,n
```
- a*i 及びb*i は何らかの方法で得たパラメータの推定値である（一致推定量でなくてもいい）

最小二乗法および補助変数法は，データを蓄えておき，データ行列の逆行列を計算して推定値を得るバッチ処理である．このような方法は繰り返し計算を必要としないことからオフライン法とも呼ばれている.

実時間での予測や，異常の検出・診断には不向きである．

↑

逐次推定アルゴリズム †

N 個のデータから得られる推定値をθN とする．そこで，θN を過去の推定値及び時刻N でのデータより求めることを考える．新しくデータが得られる度に直前の推定値を修正していく方式（逐次計算式）が実現されれば，オンライン推定や実時間推定への適用が可能である．

最小二乗法と補助変数法に関する逐次推定アルゴリズムは以下のようになる．

ただし，PN は2n　x　2n の行列で，LN は2n x 1 のベクトルである．上式において，逐次最小二乗法は次のように与えられる．

また，逐次補助変数法は次のように与えられる．

推定値は，1 時刻前の推定値に，式誤差の推定値に比例した修正項を加えることによって生成される．したがって，バッチ処理に比べて記憶容量が少なくすむ．
逐次計算では逆行列の計算を必要としない．
一般に，忘却係数が小さいほど，忘却機能が強く働くので，推定速度が速く，パラメータの変動にも対応できるが，システムに関する情報を速く忘却するので，観測雑音に弱く，推定値が振動的になる．通常，時不変システムのパラメータ推定において，N = 1 としても，差し支えないが，忘却係数を時間の経過とともに1 に近づくように設定した方が，パラメータの推定値が速く真値へ収束する．
逐次アルゴリズムを実行するためには初期値の設定が必要である．通常，次のように選べばよい．