Skip to main content

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

3.1 はじめに​

システムの解析、予測、制御、シミュレーション及び異常診断のため、システムの動特性を知る必要がある。

システム同定とは

システムの入出力信号の観測値に基づいて、ある評価規範のもとでそのシステムをもっともよく記述する数学モデルを決定することをシステム同定（System Identification）という。

このように、統計学に基づいた同定法は大量かつ複雑なデータ処理を必要とする。最小二乗法で代表されるシステムのパラメータ推定手段は実システムに使える。

システム同定の概要

3.2 最小二乗法の導出​

グラフに示すような線形システムの入出力関係を表す数式モデルの記述について考える。線形システムの次数nnを既知と仮定すると、入出力関係は次式のような差分方程式で表される。

xt=−∑i=1naixt−i+∑i=1nbiut−ix_t=-\sum_{i=1}^{n}a_ix_{t-i}+\sum_{i=1}^{n}b_iu_{t-i}

ただし、utu_tとxtx_tはそれぞれシステムの入出力信号の離散時刻ttにおけるサンプル値であり、aia_iとbib_iはシステムモデルのパラメータである。

ここで、出力信号xtx_tには、平均値ゼロの離散値確率雑音vtv_tが加わると仮定する。

yt=xt+vty_t=x_t+v_t
無相関の仮定

一般に、観測雑音vtv_tはシステムの入力信号utu_tとは無相関である。すなわち、E[utvt]=0E[u_tv_t]=0。ただし、E[⋅]E[\cdot]は期待値を意味する。

我々の目的は、システムの入力utu_tと出力の観測値yty_tからシステムのパラメータを推定することである。

ここで、yty_tの式を差分方程式に代入すると：

yt=−∑i=1naiyt−i+∑i=1nbiut−i+rty_t=-\sum_{i=1}^{n}a_iy_{t-i}+\sum_{i=1}^{n}b_iu_{t-i}+r_t

が得られる。ただし、観測雑音による誤差は：

rt=∑i=0naivt−i(a0=1)r_t=\sum_{i=0}^{n}a_iv_{t-i}\quad(a_0=1)

ベクトル形式での表現​

yt=ztTθ+rty_t=z_t^T\boldsymbol{\theta}+r_t ztT=[−yt−1,…,−yt−n,ut−1,…,ut−n]z_t^T=\begin{bmatrix} -y_{t-1},\dots,-y_{t-n},u_{t-1},\dots,u_{t-n} \end{bmatrix} θ=[a1⋮anb1⋮bn]\boldsymbol{\theta}=\begin{bmatrix} a_1\\ \vdots \\a_n \\ b_1 \\ \vdots \\ b_n \end{bmatrix}

二乗誤差評価関数​

システムのパラメータを推定するため、次のような二乗誤差評価を導入する。

J=∑t=1NρN−trt2=∑t=1NρN−t(yt−ztTθ)2J=\sum_{t=1}^{N}\rho^{N-t}r_t^2=\sum_{t=1}^{N} \rho^{N-t}(y_t-z_t^T\boldsymbol{\theta})^2

ただし、0<ρ≤10<\rho\leq1は忘却係数である。

忘却係数の役割

NNを現在時刻とみれば、N−tN-tにすると、現在の時間より前の過去のtt個時間の意味である。時間ttが大きいほど、遠くなる過去値が現在に対する影響を意味する。よって、忘却係数が小さくすると、推定値への寄与（貢献）を小さくする役割を果たしており、一種の忘却機能を持っている。この忘却の効果のため、ゆっくり変動する時変パラメータの推定にも適用できる。

最小二乗推定値の導出​

JJを最小にする必要条件は∂J/∂θ=0\partial J/\partial\theta=0である。すなわち：

∂J∂θ=∂∂θ(∑t=1NρN−t(yt−ztTθ)2)=−∑t=1NρN−t2zt(yt−ztTθ)=−2∑t=1NρN−tztyt+2∑t=1NρN−tztztTθ=0 \begin{aligned} \frac{\partial J}{\partial \theta}&= \frac{\partial}{\partial \theta}\left( \sum_{t=1}^{N} \rho^{N-t}(y_t-z_t^T\boldsymbol{\theta})^2 \right)\\ &=-\sum_{t=1}^{N}\rho^{N-t}2z_t(y_t-z_t^T\boldsymbol{\theta})\\ &=-2\sum_{t=1}^{N}\rho^{N-t}z_ty_t+2\sum_{t=1}^{N}\rho^{N-t}z_tz_t^T\theta = 0 \end{aligned}
ベクトル微分

ベクトルAAとベクトルxxの積をベクトルxxで微分する： ∂∂x(Ax)=AT\frac{\partial}{\partial x}(Ax)=A^T よって： ∂∂θ(ztTθ)=zt\frac{\partial}{\partial \theta}(z_t^T\theta)=z_t

より、θ\thetaの最小二乗推定値が得られる（JJが最小にするθ^\hat\theta）：

θ^=(∑t=1NρN−tztztT)−1∑t=1NρN−tztyt\hat\theta=\left( \sum_{t=1}^{N}\rho^{N-t}z_tz_t^T \right)^{-1}\sum_{t=1}^{N}\rho^{N-t}z_ty_t

上式が成り立つためには、(∑t=1NρN−tztztT)−1\left( \sum_{t=1}^{N}\rho^{N-t}z_tz_t^T \right)^{-1}が存在しなければならない。

逆行列の存在条件

一般に、入力信号u(t)u(t)が十分多くの周波数成分を含んでいれば、上記の逆行列が存在する。

3.2.1 推定値の漸近的性質​

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

plimN→∞θ^=θ\text{plim}_{N\to\infty}\hat\theta=\theta

が成立するとき、θ^\hat\thetaは一致推定値であるという。

ここで、忘却係数ρ=1\rho=1、すなわち過去のデータを忘却しない場合について考察する。θ^\hat\thetaの式にyt=ztTθ+rty_t=z_t^T\theta+r_tを代入すると：

plimN→∞θ^=(plimN→∞1N∑t=1NztztT)−1(plimN→∞1N∑t=1Nzt{ztTθ+rt})=θ+(plimN→∞1N∑t=1NztztT)−1(plimN→∞1N∑t=1Nztrt) \begin{aligned} \text{plim}_{N\to\infty}\hat\theta &=\left( \text{plim}_{N\to\infty} \frac{1}{N}\sum_{t=1}^{N}z_tz_t^T \right)^{-1} \left( \text{plim}_{N\to\infty} \frac{1}{N} \sum_{t=1}^{N}z_t\{z_t^T\theta+r_t\} \right)\\ &=\theta+\left( \text{plim}_{N\to\infty} \frac{1}{N} \sum_{t=1}^{N}z_tz_t^T \right)^{-1}\left(\text{plim}_{N\to\infty} \frac{1}{N} \sum_{t=1}^{N}z_tr_t\right) \end{aligned}

が得られる。ベクトルztz_tと誤差rtr_tは無関係でなければならない、すなわち：

plimN→∞1N∑t=1Nztrt=0\text{plim}_{N\to\infty}\frac{1}{N}\sum_{t=1}^{N}z_tr_t=0

という条件が成り立たなければならない。

バイアスの問題

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

3.3 逐次最小二乗法（Recursive Least Squares）​

前に導出した最小二乗法は、データを蓄えておき、データ行列の逆行列を計算して推定値を得るバッチ処理である。このような方法は繰り返し計算を必要としないことからone-shot法またはオフライン法とも呼ばれている。

バッチ処理の限界

システムの特性を実時間での監視や、異常の検出と診断には不向きである。

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

逐次最小二乗法のアルゴリズム​

θ^N=θ^N−1+LNεN\hat\theta_N=\hat\theta_{N-1}+L_N\varepsilon_N
積分の形

過去のデータ + 比例係数 × 偏差

εN=yN−zNTθ^N−1\varepsilon_N=y_N-z_N^T\hat\theta_{N-1}
式誤差（偏差）

観測値 - 目標値

LN=PN−1zNρN+zNTPN−1zNL_N=\frac{P_{N-1}z_N}{\rho_N+z_N^TP_{N-1}z_N}
比例係数

ゲイン係数

PN=1ρN[PN−1−PN−1zNzNTPN−1ρN+zNTPN−1zN]P_N=\frac{1}{\rho_N}\left[P_{N-1}-\frac{P_{N-1}z_Nz_N^TP_{N-1}}{\rho_N+z_N^TP_{N-1}z_N}\right]

ただし、PNP_Nは2n×2n2n\times2nの行列で、LNL_Nは2n×12n\times1のベクトルである。

逐次最小二乗法の特徴​

  1. 推定値θ^N\hat\theta_Nは、1時刻前の推定値θ^N−1\hat\theta_{N-1}に、式誤差yN−zNTθ^N−1y_N-z_N^T\hat\theta_{N-1}に比例した修正項を加えることによって生成される。したがって、バッチ処理に比べて記憶容量が少なくすむ。また、逐次計算では逆行列の計算を必要としない。

  2. 上述の逐次アルゴリズムを実行するためには初期値の設定が必要である。通常、次のように選べばよい：

θ^0=0,P0=αI\hat\theta_0=0, \quad P_0=\alpha I

ただし、α=103∼105\alpha=10^3 \sim 10^5

忘却係数の設定​

忘却係数ρN\rho_Nも推定アルゴリズムの性能を大きく左右する。

忘却係数の特性

忘却係数が小さいほど、忘却機能が強く働くので、推定速度が速く、パラメータの変動にも対応できるが、システムに関する情報を速く忘却するので、観測雑音に弱く、推定値が振動的になる。

時不変システムの場合​

通常、時不変システムのパラメータ推定において、ρN=1\rho_N=1としても差し支えないが、忘却係数を時間の経過とともに1に近づくように設定したほうが、パラメータの推定値が速く真値へ収束する。

ρN=(1−0.01)ρN−1+0.01,ρ0=0.95\rho_N=(1-0.01)\rho_{N-1}+0.01, \quad \rho_0=0.95

すなわち、最初のうちに忘却係数を小さくすることによって推定速度を上げるが、推定値がある程度その真値に近づいたら、雑音の影響を軽減するために、忘却係数を1とすることによって推定速度を落とす。

時変システムの場合​

一方、時変システムのパラメータ推定において、ρN=0.97∼0.999\rho_N=0.97\sim0.999と、1よりやや小さい定数にすればよい。忘却係数が小さいほど、時変パラメータへの追従速度が速くなるが、雑音の影響を受けやすくなる。

時変システムの制限

いずれにしても、雑音のレベルが高い場合、時変システムのパラメータ推定は困難である。

3.4 補助変数法によるパラメータ一致推定​

ベクトルztz_tと同じサイズの補助変数ベクトルmtm_tを導入し、次のように補助変数法による推定値を得る。

θ^IV=(∑t=1NρN−tmtztT)−1∑t=1NρN−tmtyt\hat\theta_{IV}=\left( \sum_{t=1}^{N}\rho^{N-t}m_tz_t^T \right)^{-1} \sum_{t=1}^{N}\rho^{N-t}m_ty_t

そこで、忘却係数ρ=1\rho=1の場合（過去のデータを忘却しない）における補助変数推定値の漸近的性質について考察する。上式の推定値にyt=ztTθ+rty_t=z_t^T\theta+r_tを代入して、確率極限をとると：

plimN→∞θ^IV=(plimN→∞1N∑t=1NmtztT)−1(plimN→∞1N∑t=1Nmt{ztTθ+rt})=θ+(plimN→∞1N∑t=1NmtztT)−1(plimN→∞1N∑t=1Nmtrt) \begin{aligned} \text{plim}_{N\to\infty}\hat\theta_{IV} &= \left(\text{plim}_{N\to\infty}\frac{1}{N}\sum_{t=1}^{N} m_tz_t^T \right)^{-1}\left( \text{plim}_{N\to\infty} \frac{1}{N}\sum_{t=1}^{N}m_t\{ z_t^T\theta+r_t \} \right)\\ &=\theta+\left(\text{plim}_{N\to\infty}\frac{1}{N}\sum_{t=1}^{N} m_tz_t^T \right)^{-1}\left( \text{plim}_{N\to\infty} \frac{1}{N}\sum_{t=1}^{N}m_tr_t\right) \end{aligned}

一致推定の条件​

補助変数ベクトルmtm_tが次のような条件を満足すれば、パラメータ推定値は一致推定値である：

  1. 逆行列(plimN→∞1N∑t=1NmtztT)−1\left(\text{plim}_{N\to\infty}\frac{1}{N}\sum_{t=1}^{N} m_tz_t^T \right)^{-1}が存在する。

  2. 補助変数ベクトルmtm_tと式誤差rtr_tは無相関である。すなわち： plimN→∞1N∑t=1Nmtrt=0\text{plim}_{N\to\infty}\frac{1}{N}\sum_{t=1}^{N}m_tr_t=0

補助変数ベクトルの構成方法​

(1) 遅延された入力を利用​

mt=[ut−1−d,…,ut−n−d,ut−1,…,ut−n]Tm_t=\left[ u_{t-1-d},\dots,u_{t-n-d},u_{t-1},\dots,u_{t-n} \right]^T

ただし、d≥nd\geq n

(2) 遅延された出力を利用​

mt=[yt−1−d,…,yt−n−d,ut−1,…,ut−n]Tm_t=[y_{t-1-d},\dots,y_{t-n-d},u_{t-1},\dots,u_{t-n}]^T

ただし、d≥nd\geq n

(3) 推定された出力を利用（Bootstrap法）​

mt=[−x^t−1,…,−x^t−n,ut−1,…,ut−n]Tm_t=\left[ -\hat x_{t-1},\dots,-\hat x_{t-n},u_{t-1},\dots,u_{t-n} \right]^T

ただし、推定された出力： x^t=−∑i=1na^ix^t−i+∑i=1nb^iut−i\hat x_t=-\sum_{i=1}^{n}\hat a_i\hat x_{t-i}+\sum_{i=1}^{n}\hat b_iu_{t-i}

ここで、a^i,b^i\hat a_i, \hat b_iは何らかの方法で得たパラメータの推定値である（一致推定値でなくてもよい）。このように、システムの出力とパラメータを交互に推定することを**Bootstrap法（靴のひもかけ）**という。

3.4.1 逐次推定アルゴリズム​

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

θ^N=θ^N−1+LNεN\hat\theta_N=\hat\theta_{N-1}+L_N\varepsilon_N εN=yN−ϕNTθ^N−1\varepsilon_N=y_N-\phi_N^T\hat\theta_{N-1} LN=PN−1ψNρN+ϕNTPN−1ψNL_N=\frac{P_{N-1}\psi_N}{\rho_N+\phi_N^TP_{N-1}\psi_N} PN=1ρN[PN−1−PN−1ψNϕNTPN−1ρN+ϕNTPN−1ψN]P_N=\frac{1}{\rho_N}\left[P_{N-1}-\frac{P_{N-1}\psi_N\phi_N^TP_{N-1}}{\rho_N+\phi_N^TP_{N-1}\psi_N}\right]

ただし、PNP_Nは2n×2n2n\times2nの行列で、LNL_Nは2n×12n\times1のベクトルである。

上式において：

逐次最小二乗法の場合： ϕN=zN,ψN=zN\phi_N=z_N, \quad \psi_N=z_N

逐次補助変数法の場合： ϕN=zN,ψN=mN\phi_N=z_N, \quad \psi_N=m_N

Bootstrap補助変数法の注意点​

Bootstrap補助変数法に関して、下記の点に留意しておくべきである。

推定された出力を使って補助変数を構成する場合、推定されたシステムの安定性をチェックしなければならない。

補助変数ベクトルは次のように構成される：

mN=[−x^N−1,…,−x^N−n,uN−1,…,uN−n]Tm_N=\left[ -\hat x_{N-1},\dots,-\hat x_{N-n},u_{N-1},\dots,u_{N-n} \right]^T

ただし： x^N=−∑i=1naˉi,N−1x^N−i+∑i=1nb^i,N−1uN−i\hat x_N=-\sum_{i=1}^{n}\bar a_{i,N-1}\hat x_{N-i}+\sum_{i=1}^{n}\hat b_{i,N-1}u_{N-i}

ここで：

  • b^i,N−1\hat b_{i,N-1}は逐次推定アルゴリズムから得られた時刻N−1N-1における推定値
  • aˉi,N−1\bar a_{i,N-1}は安定性を考慮して修正された推定値

安定性チェック​

時刻NNにおいて、推定出力x^N\hat x_Nを計算するとき、上記システムの極が安定でなければ、値が発散してしまう。したがって、システムの極：

zn+aˉ1,N−1zn−1+⋯+aˉn,N−1=0z^n+\bar a_{1,N-1}z^{n-1}+\dots+\bar a_{n,N-1}=0

の解が単位円の内部にあれば（複素数の半径が1より小さい）、システムが安定である。

安定性修正

a^i,N−1(i=1,2,…,n)\hat a_{i,N-1}(i=1,2,\dots,n)が安定な極を構成していれば、aˉi,N−1=a^i,N−1(i=1,2,…,n)\bar a_{i,N-1}=\hat a_{i,N-1}(i=1,2,\dots,n)とする。そうでなければ、a^i,N−1(i=1,2,…,n)\hat a_{i,N-1}(i=1,2,\dots,n)を適当な安定パラメータに置き換える。

射影による修正手法​

複雑な手法として、以下のような「射影」による修正手法もある：

アルゴリズム：

  • Step a:
    • If a^i,N−1(i=1,2,…,n)\hat a_{i,N-1}(i=1,2,\dots,n) stable then
      • aˉi,N−1=a^i,N−1(i=1,2,…,n)\bar a_{i,N-1}=\hat a_{i,N-1}(i=1,2,\dots,n)
      • Stop
    • Else
      • aˉi,N−1=aˉi,N−2+0.5∗(a^i,N−1−aˉi,N−2)\bar a_{i,N-1}=\bar a_{i,N-2}+0.5*(\hat a_{i,N-1}-\bar a_{i,N-2})
    • Endif
  • Step b:
    • If aˉi,N−1(i=1,2,…,n)\bar a_{i,N-1}(i=1,2,\dots,n) stable then
      • Stop
    • Else
      • aˉi,N−1=aˉi,N−2+0.5∗(aˉi,N−1−aˉi,N−2)\bar a_{i,N-1}=\bar a_{i,N-2}+0.5*(\bar a_{i,N-1}-\bar a_{i,N-2})
    • Endif
  • Step c: Go to Step b
修正の範囲

上述の修正は、推定出力の計算だけのために行われることに注意されたい。逐次推定アルゴリズムに対しては、推定されたシステムの安定性の如何に関わらず、このような修正をしない。

最初のうち、パラメータ推定値が揺れるので、最初の50〜100ステップの推定を最小二乗法で行い、そのあとBootstrap補助変数法に切り替えたほうがよい。

3.6 逐次最小二乗法の導出​

最小二乗法による二乗誤差評価∂J/∂θ=0\partial J/\partial\theta=0のとき、最小二乗法推定値θ^\hat{\theta}により：

θ^=(∑t=1NρN−tztztT)−1∑t=1NρN−tztyt\hat{\theta}=\left(\sum_{t=1}^N\rho^{N-t}z_t z_t^T\right)^{-1}\sum_{t=1}^N\rho^{N-t}z_t y_t

が得られる。このとき、逆行列を仮定すると：

PN=(∑t=1NρN−tztztT)−1P_N=\left(\sum_{t=1}^N\rho^{N-t}z_t z_t^T\right)^{-1}

よって：

忘却係数の影響

忘却係数ρ\rhoが時不変であれば、ρt=ρN\rho_t=\rho_Nである。これによって、忘却係数ρ\rhoが小さくなると、昔のことを忘れてPNP_Nが大きくなる。

PN−1=(∑t=1NρNN−tztztT)=(ρNN−1z1z1T+ρNN−2z2z2T+⋯+ρNzN−1zN−1T)+ρN0zNzNT=ρN∑t=1N−1ρNN−1−tztztT+zNzNT=ρNPN−1−1+zNzNT \begin{aligned} P_N^{-1}&=\left(\sum_{t=1}^N\rho_N^{N-t}z_t z_t^T\right)\\ &=\left(\rho_N^{N-1}z_1z_1^T+\rho_N^{N-2}z_2z_2^T+\dots+\rho_{N}z_{N-1}z_{N-1}^T\right)+\rho_{N}^0z_{N}z_{N}^T\\ &=\rho_N\sum_{t=1}^{N-1}\rho_N^{N-1-t}z_t z_t^T+z_{N}z_{N}^T\\ &=\rho_NP_{N-1}^{-1}+z_Nz_N^T \end{aligned}

Woodbury逆行列補題の適用​

逆行列の補題（Woodbury恒等式）

(A−1+CTB−1D)−1=A−ACT(DACT+B)−1DA(A^{-1}+C^TB^{-1}D)^{-1}=A-AC^T(DAC^T+B)^{-1}DA

A=PN−1ρN,B=1,CT=zN,D=zNTA=\frac{P_{N-1}}{\rho_N}, \quad B=1, \quad C^T=z_N, \quad D=z_N^T

したがって：

PN=1ρNPN−1−1ρNPN−1zN[1+zNTPN−1ρNzN]−1zNTPN−11ρN=1ρN[PN−1−PN−1zNzNTPN−1ρN+zNTPN−1zN] \begin{aligned} P_N&=\frac{1}{\rho_N}P_{N-1}-\frac{1}{\rho_N}P_{N-1}z_N\left[1+z_N^T\frac{P_{N-1}}{\rho_N}z_N\right]^{-1}z_N^TP_{N-1}\frac{1}{\rho_N}\\ &=\frac{1}{\rho_N}\left[ P_{N-1} - \frac{P_{N-1}z_Nz_N^TP_{N-1}}{\rho_N+z_N^TP_{N-1}z_N} \right] \end{aligned}

PNP_Nを最小二乗推定値θ^\hat\thetaに代入すると：

θ^N=(∑t=1NρNN−tztztT)−1∑t=1NρNN−tztyt=PN[ρN∑t=1N−1ρNN−1−tztyt+zNyN] \begin{aligned} \hat\theta_N&=\left(\sum_{t=1}^N\rho_N^{N-t}z_t z_t^T\right)^{-1}\sum_{t=1}^N\rho_N^{N-t}z_t y_t\\ &=P_N\left[\rho_N\sum_{t=1}^{N-1}\rho_N^{N-1-t}z_t y_t+z_{N}y_{N}\right] \end{aligned}

整理すると：

θ^N=PN[ρNPN−1−1θ^N−1+zNyN]=PN[[PN−1−zNzNT]θ^N−1+zNyN]=θ^N−1+PNzN[yN−zNTθ^N−1] \begin{aligned} \hat\theta_N&=P_N\left[ \rho_N P_{N-1}^{-1}\hat\theta_{N-1}+z_Ny_N \right]\\ &=P_N\left[ [P_N^{-1}-z_Nz_N^T]\hat\theta_{N-1}+z_Ny_N \right]\\ &=\hat\theta_{N-1}+P_Nz_N\left[y_N-z_N^T\hat\theta_{N-1}\right] \end{aligned}
ゲインの意味

LN=PNzNL_N=P_Nz_NでPNP_Nが大きくなると、LNL_Nが大きくなる。偏差係数が大きくなるので、偏差に対応する変化が大きいため、前のことを忘れた。

PNzN=LN=1ρN[PN−1−PN−1zNzNTPN−1ρN+zNTPN−1zN]zN=PN−1zNρN+zNTPN−1zN \begin{aligned} P_Nz_N&=L_N\\ &=\frac{1}{\rho_N}\left[ P_{N-1} - \frac{P_{N-1}z_Nz_N^TP_{N-1}}{\rho_N+z_N^TP_{N-1}z_N} \right]z_N\\ &=\frac{P_{N-1}z_N}{\rho_N+z_N^TP_{N-1}z_N} \end{aligned}

3.7 射影と最小二乗法の幾何学的意味​

3.7.1 ベクトル空間と部分空間​

実数集合R\mathbb{R}または複素数集合C\mathbb{C}などは、加減乗除の演算、いわゆる四則演算が定義されている。このような集合を**体（field）**という。体FF上のベクトル空間XXは以下のように定義される。

定義 3.7.1：線形ベクトル空間​

x,y∈Xx, y\in Xおよびα∈F\alpha\in Fについて、x+y∈Xx+y\in Xとαx∈X\alpha x\in Xが定義され、以下の公理（axiom）が満足されるとき、XXを体FF（R\mathbb{R}あるいはC\mathbb{C}）上の線形ベクトル空間という。

線形空間の公理

FF上のベクトル空間XXの部分空間（subspace）SSはXXの空ではない部分集合であり、SS自身も同じFF上のベクトル空間である。

3.7.2 行列の階数と次元​

定義 3.7.2：一次独立性​

Cm\mathbb{C}^m上のベクトルの組{x1,…,xk}\{ x_1,\dots,x_k \}は、その**一次結合（線形結合）**が0であるとき： α1x1+⋯+αkxk=0\alpha_1x_1+\dots+\alpha_kx_k=0 α1=α2=⋯=αk=0\alpha_1=\alpha_2=\dots=\alpha_k=0であれば、**一次独立（linearly independent）**であるという。

逆に、すべてが0ではないあるスカラー組α1,…,αk\alpha_1,\dots,\alpha_kによって、{x1,…,xk}\{ x_1,\dots,x_k \}の一次結合がゼロであれば、{x1,…,xk}\{ x_1,\dots,x_k \}は**一次従属（linearly dependent）**であるという。

行列表現による理解

xi=(ai,bi,…,ni)Tx_i=(a_i,b_i,\dots,n_i)^Tがベクトル、αi\alpha_iがスカラーで： α1x1+⋯+αkxk=0⇒[x1  …  xk][α1⋮αk]=Ax=0\alpha_1x_1+\dots+\alpha_kx_k=0 \Rightarrow [x_1\;\dots\;x_k]\begin{bmatrix}\alpha_1\\\vdots\\\alpha_k\end{bmatrix}=Ax=0

ベクトルxix_iを展開したら、一次結合（線形結合）が行列とベクトルの掛け算がゼロという形になる。

  • この線形変換が成り立つベクトルが、ベクトルα\alphaが零ベクトルしかいないと、**一次独立（線形独立）**である。つまり、ベクトル組xxでの行列の階数が次元数と同じであり、線形変換した空間の次元が圧縮されない。
  • 逆にこの線形変換が成り立つベクトルが、ベクトルα\alpha以外に存在すると、**一次従属（線形従属）**である。つまり、ベクトル組xxでの行列の階数が次元数より小さく、線形変換した空間の次元が圧縮された。

定義 3.7.3：部分空間​

ベクトル組{x1,…,xk}∈Cn\{ x_1,\dots,x_k \}\in\mathbb{C}^nの一次結合の全体は、部分空間となり、{x1,…,xk}\{ x_1,\dots,x_k \}によって張られる**部分空間（subspace spanned by vectors）**と呼ばれる。以下のように表される：

span{x1,…,xk}={∑i=1nβixi,  βi∈C}\text{span}\{ x_1,\dots,x_k \}=\left\{ \sum_{i=1}^{n} \beta_i x_i, \;\beta_i\in\mathbb{C}\right\}

{x1,…,xn}\{ x_1,\dots,x_n \}が一次独立（線形独立）で、かつb∈span{x1,…,xn}b\in \text{span}\{ x_1,\dots,x_n \}であれば、bbはx1,…,xnx_1,\dots,x_nによる唯一の一次結合（線形結合）である。

基底ベクトルの特性

基底ベクトルはお互いに独立であれば、基底ベクトルで張られる空間に存在するベクトルは、当然唯一の座標がある。

X∈Rn×mX\in\mathbb{R}^{n\times m}の列ベクトルによって張られる部分空間R(X)\mathcal{R}(X)をその列空間という。

列空間と行空間​

R(X)={Xの列ベクトル{x1,…,xM}で張られる部分空間}\mathcal{R}(X) =\{ X \text{の列ベクトル}\{x_1,\dots,x_M\} \text{で張られる部分空間} \}

R(XT)={Xの行ベクトルで張られる部分空間}\mathcal{R}(X^T)=\{ X \text{の行ベクトルで張られる部分空間}\}

行列XT∈RN×MX^T\in\mathbb{R}^{N\times M}について、階数rrは、独立な列の最大個数あるいは独立な行の最大個数である。

dim⁡R(X)=dim⁡R(XT)=r\dim\mathcal{R}(X)=\dim\mathcal{R}(X^T)=r

左零空間​

αT(x1,…,xk)=xTA=0\alpha^T(x_1,\dots,x_k)=x^TA=0

左零空間： N(XT)={n∈RN∣XTn=0}\mathcal{N}(X^T)=\{n\in\mathbb{R}^N| X^Tn=0\}

dim⁡N(XT)=N−r\dim\mathcal{N}(X^T)=N-r

列空間と左零空間が直交している： R(X)  ⊥  N(XT),R(X)⊥=N(XT),N(XT)⊥=R(X)\mathcal{R}(X)\;\bot\;\mathcal{N}(X^T), \quad \mathcal{R}(X)^{\bot}=\mathcal{N}(X^T), \quad \mathcal{N}(X^T)^{\bot}=\mathcal{R}(X)

列空間と左零空間の和空間は、空間全体である： R(X)⊕N(XT)=RN\mathcal{R}(X)\oplus \mathcal{N}(X^T)=\mathbb{R}^N

零空間（核空間）​

N(X)={h∈RM∣Xh=0}\mathcal{N}(X)=\{h\in\mathbb{R}^M| Xh=0\}

dim⁡N(X)=M−r\dim\mathcal{N}(X)=M-r

R(XT)  ⊥  N(X),R(XT)⊥=N(X),N(X)⊥=R(XT)\mathcal{R}(X^T)\;\bot\;\mathcal{N}(X), \quad \mathcal{R}(X^T)^{\bot}=\mathcal{N}(X), \quad \mathcal{N}(X)^{\bot}=\mathcal{R}(X^T)

R(XT)⊕N(X)=RM\mathcal{R}(X^T)\oplus \mathcal{N}(X)=\mathbb{R}^M

空間の関係

3.7.3 射影と最小二乗法​

直線への射影＝ベクトルの内積​

yTx=xTy=∣∣x∣∣⋅∣∣y∣∣cos⁡θy^Tx=x^Ty=||x||\cdot||y||\cos\theta

射影​

yyがxxへの射影yr=axy_r=axとおくと（aaがスカラー、単なる比例を表す）：

(y−yr)Tx=(y−ax)Tx=0,a=yTxxTx=xTyxTx(y-y_r)^Tx=(y-ax)^Tx=0, \quad a=\frac{y^Tx}{x^Tx}=\frac{x^Ty}{x^Tx}
射影の意味

つまり、ベクトルがベクトルの線形結合である。元のyyベクトルをxxベクトルへの方向の成分を除くと、垂直成分しかないので、内積（射影）はゼロである。

yr=hx=xTyxTxx=[(xTx)−1xxT]y=[x(xTx)−1xT]y\boldsymbol{y}_r=h\boldsymbol{x}=\frac{\boldsymbol{x}^T\boldsymbol{y}}{\boldsymbol{x}^T\boldsymbol{x}}\boldsymbol{x}=\left[(\boldsymbol{x}^T\boldsymbol{x})^{-1}\boldsymbol{x}\boldsymbol{x}^T\right]\boldsymbol{y}=\left[\boldsymbol{x}(\boldsymbol{x}^T\boldsymbol{x})^{-1}\boldsymbol{x}^T\right]\boldsymbol{y}

射影と最小二乗法​

モデル：

y=yr+ynyr=XhX=[x1,⋯ ,xM]∈RN×Mh=[h1,⋯ ,hM]T\begin{aligned} &\boldsymbol{y}=\boldsymbol{y}_r+\boldsymbol{y}_n \\ &\boldsymbol{y}_{r}= Xh \\ &X= [\boldsymbol{x}_1,\cdots,\boldsymbol{x}_M]\in \mathbb{R}^{N\times M} \\ &h= [h_1,\cdots,h_M]^T \end{aligned}

XXの列空間R(X)\mathcal{R}(X)上の任意のベクトルXgXgについて考えると、yny_nがXgXgに垂直すれば：

(Xg)Tyn=(Xg)T(y−yr)=gTXT(y−Xh)=gT(XTy−XTXh)=0\begin{aligned} (\boldsymbol{X}\boldsymbol{g})^{T}\boldsymbol{y}_{n}& =(\boldsymbol{X}\boldsymbol{g})^T(\boldsymbol{y}-\boldsymbol{y}_r) \\ &=\boldsymbol{g}^T\boldsymbol{X}^T(\boldsymbol{y}-\boldsymbol{X}\boldsymbol{h}) \\ &=\boldsymbol{g}^T(\boldsymbol{X}^T\boldsymbol{y}-\boldsymbol{X}^T\boldsymbol{X}\boldsymbol{h})=0 \end{aligned}

よって、すべてのggについて成り立つために：

XTy=XTXh=XTyrX^Ty=X^TXh=X^Ty_r

部分空間への射影

つまり、部分空間XXに存在する任意のベクトルggに対して、あるベクトルyyのベクトル成分yny_nと垂直すれば、yyが部分空間XXへの射影はyry_rである。

すなわち： h=(XTX)−1XTyh=(X^TX)^{-1}X^Ty

まとめ​

本章では線形システムのパラメータ推定について以下の重要な手法を学習した：

  1. 最小二乗法: バッチ処理によるパラメータ推定の基礎
  2. 逐次最小二乗法: リアルタイム推定のための逐次アルゴリズム
  3. 補助変数法: 一致推定を実現するための手法
  4. 射影理論: 最小二乗法の幾何学的解釈

これらの手法は、システム同定における実用的なパラメータ推定の基盤となる重要な理論である。特に逐次最小二乗法は、オンライン推定や適応制御において広く応用されている。