　はじめまして、りょーつといいます。高専出身の大学院2年生です。研究の専門は力学や機構学で、Qiitaでは主に制御工学や数学に関する記事を書いています。統計学の知識を少し使ってモータのパラメータ同定をしてみたので、それらについてまとめようと思います。本記事はこちらの記事の内容をより一般化したものであり、実験データの自己相関を考慮したものとなります。統計学に強い友人から教わったものを自分なりに理解して復習のためにもまとめようと思いました。
　数学的にも洗練された内容で面白いなあと思ったのでぜひ読んでみてください。

目次

1.はじめに
 2.数学的準備
 3.モータ差分方程式と共分散行列
 4.一般化最小二乗法によるパラメータ同定
 5.おわりに

1. はじめに

　今回の記事では、一般化最小二乗法によるモータのパラメータ同定についてまとめます。理論と実験の2本構成を予定しており、今週は理論的背景についてしっかりまとめる予定です。数学的な部分がたくさん出てきますが、丁寧に書くので追ってみていただけると嬉しいです。また、本記事の内容は今月中旬に開催されるキャチロボ2025で使用予定です。

2. 数学的準備

2.1 最小二乗法　

　まずは通常の線形なシステムに対する最小二乗法について復習します。通常の最小二乗法は$n$個のデータから$m$個のパラメータを推定する手法で、実験データの出力と$\boldsymbol{y} \in \mathbb{R}^n$、入力$X \in \mathbb{R}^n \times \mathbb{R}^m$からパラメータ$\boldsymbol{\beta} \in \mathbb{R}^m$を推定します。各実験データの誤差をまとめたベクトルを$\boldsymbol{\varepsilon} \in \mathbb{R}^n$とすると、システムは以下のように表されます。

\boldsymbol{y}
=
X \boldsymbol{\beta} + \boldsymbol{\varepsilon}
\tag{1}

最小二乗法は(1)式における誤差$\boldsymbol{\varepsilon}$のノルムが極小になるようなパラメータ$\boldsymbol{\beta}$を推定する手法であり、その解は以下のように与えられます。導出過程は以前の記事に詳しく書いたので、そちらを参照ください。

\boldsymbol{\beta}
=
{(X^T X)}^{-1} X^T \boldsymbol{y}
\tag{2}

ここで注意しないといけないのが、通常の最小二乗法は誤差$\boldsymbol{\varepsilon}$が自己相関を持たないという条件下で使用する必要があります。
　では$\boldsymbol{\varepsilon}$の中身が

\boldsymbol{\varepsilon}
=
{\begin{bmatrix}
\varepsilon_1 \\
\varepsilon_2 \\
\vdots \\
\varepsilon_n \\
\end{bmatrix}}
\tag{3}

のように構成されていたとします。誤差には(きっと)バラツキがあるので、同じ入力条件で実験したとしても同じ出力が得られるとは限りません。そこで分散を考えてみましょう。
　自己相関が無いというのは、これらの実験データの誤差はそれぞれ独立して発生するということを意味しており、ことなるデータの誤差の共分散が0になるということを意味しています。共分散は記号"Cov"を使って表すことが多いので、これを式で表すと以下のようになります。

\rm{Cov}(\varepsilon_i, \varepsilon_j) = 0 \ \ \ \ \ (i \neq j)
\tag{4}

ちなみに分散は記号"Var"を使うことが多いので、以下のように定義しておきます。

\rm{Cov}(\varepsilon_i, \varepsilon_i) = \rm{Var}(\varepsilon_i)
\tag{5}

データの誤差どうしの共分散を全て書き下すのは大変なので、(4)式の結果を行列$\Omega$にまとめてみましょう。

\Omega
=
{\begin{bmatrix}
\rm{Var}(\varepsilon_1) & \rm{Cov}(\varepsilon_1, \varepsilon_2) & \cdots &  \rm{Cov}(\varepsilon_1, \varepsilon_n) \\
\rm{Cov}(\varepsilon_2, \varepsilon_1) & \rm{Var}(\varepsilon_2) & \cdots & \rm{Cov}(\varepsilon_2, \varepsilon_n)\\
\vdots & \vdots & \ddots & \vdots \\
\rm{Cov}(\varepsilon_n, \varepsilon_1) & \rm{Cov}(\varepsilon_n, \varepsilon_2) & \cdots & \rm{Var}(\varepsilon_n)\\
\end{bmatrix}}
\tag{6}

(6)式のように定義される行列$\Omega$を共分散行列といいます。最小二乗法を使用するためには、(4)式を満たす必要があるので、この時の共分散行列$\Omega$は対角行列になります。

　なお、誤差の期待値(平均値)が零であると仮定した場合、

E[\boldsymbol{\varepsilon}] = \boldsymbol{0}
\tag{7}

分散と共分散はそれぞれ誤差の積の期待値として以下のように表すことができます。

\rm{Var}(\varepsilon_i) 
=
E[(\varepsilon_i - \bar{\varepsilon_i})^2]
=
E[{\varepsilon_i}^2]
\tag{8}

\rm{Cov}(\varepsilon_i, \varepsilon_j) 
=
E[(\varepsilon_i - \bar{\varepsilon_i})(\varepsilon_j - \bar{\varepsilon_j})]
=
E[{\varepsilon_i}{\varepsilon_j}]
\tag{9}

(8)式と(9)式の関係を(6)式に代入して整理すると、共分散行列$\Omega$も期待値を使ってシンプルに記述することができます！

\Omega
=
\begin{bmatrix}
E[{\varepsilon_1}^2] & E[{\varepsilon_1}{\varepsilon_2}] & \cdots & E[{\varepsilon_1}{\varepsilon_n}] \\
E[{\varepsilon_2}{\varepsilon_2}] & E[{\varepsilon_2}^2] & \cdots & E[{\varepsilon_2}{\varepsilon_n}] \\
\vdots & \vdots & \ddots & \vdots \\
E[{\varepsilon_n}{\varepsilon_1}] & E[{\varepsilon_n}{\varepsilon_2}] & \cdots & E[{\varepsilon_n}^2]
\end{bmatrix}
=
E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T]
\tag{10}

本記事では、共分散行列$\Omega$が非対角の場合の最小二乗法について考えていきます。

2.2 共分散行列の性質

　共分散行列$\Omega$に関するいくつかの性質について導出しておきます。まずは共分散行列$\Omega$が正定値対称行列であることを示します。共分散行列$\Omega$が対称行列であることは(10)式より明らかですね。

\Omega^T
=
E\big[ \big( \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T \big)^T \big]
=
E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T]
=
\Omega
\tag{11}

では正定値性について証明しましょう。共分散行列を$\Omega \in \mathbb{R}^n \times \mathbb{R}^n$としたとき、ある非零の実ベクトル$\boldsymbol{x} \in \mathbb{R}^n$に対して以下の関係がなりたつことを正定値性と言います。

\boldsymbol{x}^T \Omega \ \boldsymbol{x}
>
0
\tag{12}

正定値性は(12)式に(10)式を代入し、内積の性質と行列の結合則を利用することで示すことができます。

\boldsymbol{x}^T E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T] \ \boldsymbol{x}
=
E[\boldsymbol{x}^T \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T \boldsymbol{x}]
=
E[(\boldsymbol{x}^T \boldsymbol{\varepsilon} )(\boldsymbol{\varepsilon}^T \boldsymbol{x})]
=
E[(\boldsymbol{x}^T \boldsymbol{\varepsilon} )(\boldsymbol{x}^T \boldsymbol{\varepsilon} )]
=
E[(\boldsymbol{x}^T \boldsymbol{\varepsilon} )^2]
>
0
\tag{13}

言葉で書くと、「あるスカラー量$\boldsymbol{x}^T \boldsymbol{\varepsilon} $を二乗した物の期待値は0以上になる」ということから(12)式を証明しています。

　正定値対称行列には、固有値$\lambda$が全て実数かつ正になるという特徴があります。仮に共分散行列$\Omega$の固有値$\lambda$に対応する固有ベクトルを$\boldsymbol{v}$としたとき、(12)式の正定値性より以下が成り立ちます。

\boldsymbol{v}^T \Omega \ \boldsymbol{v}
=
\boldsymbol{v}^T \lambda \ \boldsymbol{v}
=
\lambda \boldsymbol{v}^T \boldsymbol{v}
=
\lambda {\|\boldsymbol{v}\| }^2
>
0
\tag{13}

固有ベクトル$\boldsymbol{v}$は非零ベクトルなので$\lambda > 0$が示せます。

2.3 共分散行列の逆行列の性質　

　共分散行列$\Omega$の逆行列$\Omega^{-1}$に関する特性についても考えてみましょう。結論から言うと、共分散行列の逆行列$\Omega^{-1}$も、正定値対称行列となります。まずは対称性を示しましょう。逆行列の定義より、以下がなりたちます。

\Omega \Omega^{-1} = I
\tag{14}

なお、単位行列を$I$と置いています。(14)式全体の転置を取り、(11)式を代入してみると、逆行列の対称性が示せます。

(\Omega \Omega^{-1})^T 
=
(\Omega^{-1})^T \Omega^T
=
(\Omega^{-1})^T \Omega
=
I

\therefore
(\Omega^{-1})^T 
=
\Omega^{-1}
\tag{15}

さらに共分散行列の逆行列$\Omega^{-1}$の正定値性について考えてみましょう。共分散行列$\Omega$を固有値分解したとき、対角行列が以下のように表されるとします。固有値分解は対角化の逆の操作によって元の行列を表す操作です。

\Omega
=
P \Lambda P^T
\tag{16}

\Lambda
=
\begin{bmatrix}
\lambda_1 & 0 & \cdots & 0\\
0 & \lambda_2 & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \lambda_n
\end{bmatrix}
=
\rm{diag}(\lambda_i)
\tag{17}

なお共分散行列は対称行列なので$P^T = P^{-1}$となることが知られています。このとき、共分散行列の逆行列$\Omega^{-1}$は以下のように表されることが知られています。

\Omega^{-1}
=
P \Lambda^{-1} P^T
\tag{18}

\Lambda^{-1}
=
\begin{bmatrix}
\frac{1}{\lambda_1} & 0 & \cdots & 0\\
0 & \frac{1}{\lambda_2} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \frac{1}{\lambda_n}
\end{bmatrix}
=
\rm{diag} \biggl(\frac{1}{\lambda_i} \biggr)
\tag{19}

実際に検算してみると結合則により

\Omega \Omega^{-1}
=
P \Lambda P^T P \Lambda^{-1} P^T
=
P \Lambda \Lambda^{-1} P^T
=
P \Lambda \Lambda^{-1} P^T
=
P P^T
=
I
\tag{19}

たしかになりたちますね。ということは(18)、(19)式より、共分散行列の逆行列$\Omega^{-1}$の固有値も全て正であることが読み取れます。最後に「固有値が全て正の対称行列であるならば正定値対称行列である」ということを示せば十分ですね。(12)式と同様にある非零の実ベクトル$\boldsymbol{x} \in \mathbb{R}^n$に対して以下の関係を示したいです。

\boldsymbol{x}^T \Omega^{-1} \ \boldsymbol{x}
=
\boldsymbol{x}^T P \Lambda^{-1} P^T \boldsymbol{x}
>
0
\tag{20}

ここで新たなベクトル$\boldsymbol{y} \in \mathbb{R}^n$を

\boldsymbol{y}
=
P^T \boldsymbol{x}
=
\begin{bmatrix}
y_1\\
y_2\\
\vdots\\
y_n
\end{bmatrix}
\tag{21}

と定義して(20)式の左辺を展開すると、

\boldsymbol{x}^T \Omega^{-1} \ \boldsymbol{x}
=
\boldsymbol{y}^T \Lambda^{-1} \boldsymbol{y}
=
\begin{bmatrix}
y_1&
y_2&
\cdots&
y_n
\end{bmatrix}
\begin{bmatrix}
\frac{1}{\lambda_1} & 0 & \cdots & 0\\
0 & \frac{1}{\lambda_2} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \frac{1}{\lambda_n}
\end{bmatrix}
\begin{bmatrix}
y_1\\
y_2\\
\vdots\\
y_n
\end{bmatrix}
=
\frac{{y_1}^2}{\lambda_1} 
+
\frac{{y_2}^2}{\lambda_2} 
+
\cdots
+
\frac{{y_n}^2}{\lambda_n} 
=
\sum_{i=1}^n \frac{{y_i}^2}{\lambda_i} 
\tag{22}

となります。前述のとおり、固有値$\lambda_i$は全て正なので以下のように(20)式が証明できます。

\boldsymbol{x}^T \Omega^{-1} \ \boldsymbol{x}
=
\sum_{i=1}^n \frac{{y_i}^2}{\lambda_i}
>
0
\tag{23}

2.4 共分散行列の平方根　

　ここでは共分散行列$\Omega$の平方根$\Omega^{\frac{1}{2}}$についてまとめます。平方根$\Omega^{\frac{1}{2}}$の定義は以下のとおりです。

\Omega^{\frac{1}{2}} \big( \Omega^{\frac{1}{2}} \big)^T
=
\Omega
\tag{24}

少し天下り的ですが、$\Omega$が正定値対称行列である場合に限り、(24)式を満たす平方根は以下のように表されます。

\Omega^{\frac{1}{2}}
=
P \Lambda^{\frac{1}{2}} P^T
\tag{25}

\Lambda^{\frac{1}{2}}
=
\begin{bmatrix}
\sqrt{\lambda_1} & 0 & \cdots & 0\\
0 & \sqrt{\lambda_2} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \sqrt{\lambda_n}
\end{bmatrix}
=
\rm{diag}(\sqrt{\lambda_i})
\tag{26}

本当になりたつのかどうか確かめてみましょう。(25)(26)式を(24)式に代入し、左辺を展開してみます。結合則より

\Omega^{\frac{1}{2}} \big( \Omega^{\frac{1}{2}} \big)^T
=
P \Lambda^{\frac{1}{2}} P^T \big( P \Lambda^{\frac{1}{2}} P^T \big)^T
=
P \Lambda^{\frac{1}{2}} P^T P \Lambda^{\frac{1}{2}} P^T
=
P \Lambda^{\frac{1}{2}} \Lambda^{\frac{1}{2}} P^T
\tag{27}

ですね。ここで

\Lambda^{\frac{1}{2}} \Lambda^{\frac{1}{2}} 
=
\begin{bmatrix}
\sqrt{\lambda_1} & 0 & \cdots & 0\\
0 & \sqrt{\lambda_2} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \sqrt{\lambda_n}
\end{bmatrix}
\begin{bmatrix}
\sqrt{\lambda_1} & 0 & \cdots & 0\\
0 & \sqrt{\lambda_2} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \sqrt{\lambda_n}
\end{bmatrix}
=
\begin{bmatrix}
\lambda_1 & 0 & \cdots & 0\\
0 & \lambda_2 & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \lambda_n
\end{bmatrix}
=
\Lambda
\tag{28}

であることから、

(25)式の妥当性が示されます。

\Omega^{\frac{1}{2}} \big( \Omega^{\frac{1}{2}} \big)^T
=
P \Lambda^{\frac{1}{2}} \Lambda^{\frac{1}{2}} P^T
=
P \Lambda P^T
=
\Omega
\tag{29}

同様の手順で、共分散行列の逆行列の平方根$\Omega^{-\frac{1}{2}}$を定義することも可能です。検算は省略します。

\Omega^{-\frac{1}{2}} \big( \Omega^{-\frac{1}{2}} \big)^T
=
\Omega^{-1}
\tag{30}

\Omega^{-\frac{1}{2}}
=
P \Lambda^{-\frac{1}{2}} P^T
\tag{31}

\Lambda^{-\frac{1}{2}}
=
\begin{bmatrix}
\frac{1}{\sqrt{\lambda_1}} & 0 & \cdots & 0\\
0 & \frac{1}{\sqrt{\lambda_2}} & \cdots & 0 \\
\vdots & \dots & \ddots & \vdots\\
0 & 0 & \cdots & \frac{1}{\sqrt{\lambda_n}}
\end{bmatrix}
=
\rm{diag} \bigg( \frac{1}{\sqrt{\lambda_i}} \bigg)
\tag{32}

2.5 共分散の性質　

　最後に共分散に関する法則を導出しておきます。具体的には以下の式を導出します。

\rm{Cov}(\rho \varepsilon_i + \mu_j, \varepsilon_k)
=
\rho \ \rm{Cov}(\varepsilon_i, \varepsilon_k) 
+
\rm{Cov}(\mu_j, \varepsilon_k) 
\tag{33}

なお、$\rho$は一定のスカラー値であり、$\varepsilon_i$、$\mu_j$、$\varepsilon_k$はバラつきのあるデータとします。
　同じ実験条件で$L$個のデータを取った際の(33)式左辺の共分散は以下のように定義されます。なお、$\bar{x}$は$x$の平均値を意味するものとします。

\rm{Cov}(\rho \varepsilon_i + \mu_j, \varepsilon_k)
=
\frac{1}{L}
\sum_{l=1}^L \big( (\rho \varepsilon_{il} + \mu_{jl}) - \overline{(\rho \varepsilon_{il} + \mu_{jl})} \big)(\varepsilon_{kl} - \bar{\varepsilon_{kl}})
\tag{34}

(34)式の右辺を丁寧に展開していくと以下のようにして(33)式が導出できます。

\rm{Cov}(\rho \varepsilon_i + \mu_j, \varepsilon_k)
=
\frac{1}{L}
\sum_{l=1}^L \big( (\rho \varepsilon_{il} - \overline{\rho \varepsilon_{il}}) + (\mu_{jl} - \overline{\mu_{jl}}) \big)(\varepsilon_{kl} - \bar{\varepsilon_{kl}})

=
\frac{1}{L}
\sum_{l=1}^L  (\rho \varepsilon_{il} - \overline{\rho \varepsilon_{il}})(\varepsilon_{kl} - \bar{\varepsilon_{kl}})
+
\frac{1}{L}
\sum_{l=1}^L (\mu_{jl} - \overline{\mu_{jl}}) (\varepsilon_{kl} - \bar{\varepsilon_{kl}})

=
\rho \frac{1}{L}
\sum_{l=1}^L  ( \varepsilon_{il} - \overline{ \varepsilon_{il}})(\varepsilon_{kl} - \bar{\varepsilon_{kl}})
+
\frac{1}{L}
\sum_{l=1}^L (\mu_{jl} - \overline{\mu_{jl}}) (\varepsilon_{kl} - \bar{\varepsilon_{kl}})

=
\rho \ \rm{Cov}(\varepsilon_i, \varepsilon_k) 
+
\rm{Cov}(\mu_j, \varepsilon_k) 
\tag{35}

この式はパラメータ同定を行う際に共分散行列$\Omega$を考える上でとても重要になります。

3. モータ差分方程式と共分散行列

　モータの運動を表す差分方程式は、固有の定数$a$と$b$を用いて以下のように表されることが知られています。$\omega$を角速度、$u$入力電圧とします。同じ形式の差分方程式は以前に書いた記事でも登場しているため、詳しい導出などはそちらを参照ください。

\omega_{t+1}
=
a \omega_{t}
+
b u_t
\tag{36}

(36)式に対して測定値の誤差$\varepsilon$を考慮し、$t = 1$から順に$N$個並べると以下のようになります。

\omega_{2}
=
a \omega_{1}
+
b u_1
+
\varepsilon_1
\ \ \ \ \ (t = 1)
\tag{37}

\omega_{3}
=
a \omega_{2}
+
b u_2
+
\varepsilon_2
\ \ \ \ \ (t = 2)
\tag{38}

\vdots

\omega_{N+1}
=
a \omega_{N}
+
b u_{N}
+
\varepsilon_{N}
\ \ \ \ \ (t = N)
\tag{39}

分かりやすくするために(37)～(39)式を行列にまとめてみましょう。以下のように行列とベクトルを定義します。

\boldsymbol{y}
=
\begin{bmatrix}
\omega_2\\
\omega_3\\
\vdots\\
\omega_{N+1} \\
\end{bmatrix}
\in \mathbb{R}^N
\tag{40}

X
=
\begin{bmatrix}
\omega_1 & u_2\\
\omega_2 & u_3\\
\vdots & \vdots\\
\omega_{N} & u_{N} \\
\end{bmatrix}
\in \mathbb{R}^N \times \mathbb{R}^2
\tag{41}

\boldsymbol{\beta}
=
\begin{bmatrix}
a\\
b\\
\end{bmatrix}
\in \mathbb{R}^2
\tag{42}

\boldsymbol{\varepsilon}
=
\begin{bmatrix}
\varepsilon_1\\
\varepsilon_2\\
\vdots\\
\varepsilon_{N} \\
\end{bmatrix}
\in \mathbb{R}^N
\tag{43}

このとき、差分方程式は(2)式と同じ形式で書くことができます。念のため再掲しておきます。

\boldsymbol{y}
=
X \boldsymbol{\beta} + \boldsymbol{\varepsilon}
\tag{44}

本記事では(44)式の誤差$\boldsymbol{\varepsilon}$に自己相関があると考えて最小二乗法を解きたいと思います。例えば(37)(38)式に注目してください。これらには共通するデータとして$\omega_2$が含まれています。このように$t = T$のときのデータと$t = T+1$のときのデータには似た傾向が出やすいと予想できます。そこで、(43)式の誤差のデータには以下の関係が成り立つと仮定してみましょう。

\varepsilon_{t+1}
=
\rho \varepsilon_{t} + \mu_{t+1}
\tag{45}

$\rho$は比例係数、$\mu_t$は他のデータに起因しない独立誤差とします。言い換えると、$t=2$のときの誤差$\varepsilon_2$は$t=1$のときの誤差を引きずった項$\rho \varepsilon_{1}$と新しく生成された独立誤差$\mu_2$の重ね合わせで表現できる考えます。$\mu$は独立誤差なので、$\varepsilon$との共分散は0となります。

\rm{Cov}(\varepsilon_i, \mu_j)
=
0
\tag{46}

　また、もう1つ大胆な仮定として、各時刻における誤差の分散は一定値$\gamma$になると仮定します。

\rm{Var}(\varepsilon_1)
=
\rm{Var}(\varepsilon_2)
=
\cdots
=
\rm{Var}(\varepsilon_N)
=
\gamma
\tag{47}

(45)式、(46)式のような仮定を置くことで、誤差$\boldsymbol{\varepsilon}$の共分散も求めることができます。たとえば$t = 1$と$t = 2$のときの共分散は(33)式を使って以下のように計算できます。

\rm{Cov}(\varepsilon_1, \varepsilon_2)
=
\rm{Cov}(\varepsilon_1, \rho \varepsilon_1 + \mu_2)
=
\rho \rm{Cov}(\varepsilon_1,  \varepsilon_1) + \rm{Cov}(\varepsilon_1,  \mu_2)
=
\rho \rm{Var}(\varepsilon_1)
=
\rho \gamma
\tag{48}

(48)式を一般化すると次のような関係が得られます。

\rm{Cov}(\varepsilon_t, \varepsilon_{t - k})
=
\rm{Cov}(\rho \varepsilon_{t-1} + \mu_t, \varepsilon_{t - k})

=
\rho \rm{Cov}(\varepsilon_{t-1}, \varepsilon_{t - k})

=
\rho^2 \rm{Cov}(\varepsilon_{t-2}, \varepsilon_{t - k})

\vdots

=
\rho^k \rm{Cov}(\varepsilon_{t-k}, \varepsilon_{t - k})

=
\rho^k \rm{Var}(\varepsilon_{t-k})
=
\rho^k \gamma

\therefore
\rm{Cov}(\varepsilon_t, \varepsilon_{t - k})
=
\rho^k \gamma
\tag{49}

以上のことから、モータの差分方程式における共分散行列$\Omega$は(6)式に(47)式と(49)式を代入することで以下のように求まります。

\Omega
=
\begin{bmatrix}
\gamma & \rho \gamma & \cdots & \rho^{N - 1} \gamma\\
\rho \gamma & \gamma & \cdots & \rho^{N - 2} \gamma\\
\vdots & \vdots & \ddots & \vdots\\
\rho^{N - 1} \gamma & \rho^{N - 2} \gamma & \cdots & \gamma\\
\end{bmatrix}

\therefore
\Omega
=
\gamma
\begin{bmatrix}
1 & \rho  & \cdots & \rho^{N - 1} \\
\rho  & 1 & \cdots & \rho^{N - 2} \\
\vdots & \vdots & \ddots & \vdots\\
\rho^{N - 1}  & \rho^{N - 2}  & \cdots & 1\\
\end{bmatrix}
\in \mathbb{R}^N \times \mathbb{R}^N
\tag{50}

なお(50)式の共分散行列のように、対角項が同じ要素で構成される行列を"Toeplitz 行列"と呼びます。

4. 一般化最小二乗法によるパラメータ同定

4.1 一般化最小二乗法の計算式

　本記事で考えているモータのモデルは自己相関があるために通常の最小二乗法が使用できませんでした。本章では、適切な変数変換により、自己相関のないモデルを構築し、通常の最小二乗法と同様の手順でパラメータ同定を行う手法についてまとめます。

　モータのモデルに登場した(50)式の共分散行列$\Omega$にも平方根行列が存在するため、平方根行列$\Psi \in \mathbb{R}^N \times \mathbb{R}^N$を以下のように定義します。

\Psi \Psi^T
=
\Omega^{-1}
\tag{51}

ここで(44)式の両辺に左から$\Psi^T$をかけて以下のような変数変換をしてみましょう。

\boldsymbol{\widetilde{y}}
=
\Psi^T\boldsymbol{y}
\tag{52}

\widetilde{X}
=
\Psi^T X
\tag{53}

\boldsymbol{\widetilde{\varepsilon}}
=
\Psi^T\boldsymbol{\varepsilon}
\tag{54}

\boldsymbol{y}
=
X \boldsymbol{\beta} + \boldsymbol{\varepsilon}
\to
\boldsymbol{\widetilde{y}}
=
\widetilde{X} \boldsymbol{\beta} + \boldsymbol{\widetilde{\varepsilon}}
\tag{55}

このような変数変換をした場合、(54)式で与えられる誤差$\boldsymbol{\widetilde{\varepsilon}}$は独立誤差になることが知られています。

　いきなりそんなことを言われても信じられないのでしっかり確認していきましょう。(55)式の共分散行列$\widetilde{\Omega}$を(10)式に基づいて定義し、その行列が対角行列になれば$\boldsymbol{\widetilde{\varepsilon}}$が独立誤差(各誤差間の共分散が0)であると考えられます。定義に基づいて丁寧に式変形をしていくと

\widetilde{\Omega}
=
E[\boldsymbol{\widetilde{\varepsilon}} \ \boldsymbol{\widetilde{\varepsilon}}^T]

=
E[\Psi^T \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T \Psi]

=
\Psi^T E[ \boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^T ] \Psi

=
\Psi^T \Omega \ \Psi

=
\Psi^T \big( \Omega^{-1}\big) ^{-1} \ \Psi

=
\Psi^T \big( (\Psi \Psi^T)^{-1}\big) ^{-1} \ \Psi

=
\Psi^T (\Psi^T)^{-1} \Psi^{-1} \Psi

\therefore
\widetilde{\Omega}
=
I
\tag{56}

といった具合に共分散行列$\widetilde{\Omega}$は対角行列どころか単位行列$I$になってしまうことが分かりました。つまり、(55)式に対しては通常の最小二乗法が適用できそうです。
　通常の最小二乗法の解は(2)式のように与えられるため、(2)式を(55)式に適用することでパラメータを計算してみましょう。

\boldsymbol{\beta}
=
{(\widetilde{X}^T \widetilde{X})}^{-1} \widetilde{X}^T \boldsymbol{\widetilde{y}}

=
{\big( (\Psi^T X)^T (\Psi^T X)\big)}^{-1} (\Psi^T X)^T (\Psi^T X) \boldsymbol{y}

=
{\big( X^T \Psi \Psi^T X \big)}^{-1}  X^T \Psi \Psi^T X \boldsymbol{y}

\therefore
\boldsymbol{\beta}
=
{\big( X^T \Omega^{-1} X \big)}^{-1}  X^T \Omega^{-1} X \boldsymbol{y}
\tag{57}

ここで(57)式の変数を整理します。(50)式の共分散行列$\Omega$は分散$\gamma$と定数$\rho$で構成されています。ここで変数を分離するために$\hat{\Omega}(\rho)$を以下のように定義します。

\Omega (\gamma, \rho) = \gamma \ \hat{\Omega}(\rho)
\tag{58}

\hat{\Omega}(\rho)
= 
\begin{bmatrix}
1 & \rho  & \cdots & \rho^{N - 1} \\
\rho  & 1 & \cdots & \rho^{N - 2} \\
\vdots & \vdots & \ddots & \vdots\\
\rho^{N - 1}  & \rho^{N - 2}  & \cdots & 1\\
\end{bmatrix}
\in \mathbb{R}^N \times \mathbb{R}^N
\tag{59}

(59)式を(57)式に代入することで、$\boldsymbol{\beta}$を計算する式に含まれる未知量が$\rho$のみであることが分かります。

\boldsymbol{\beta}
=
{\bigg( X^T \frac{1}{\gamma} \hat{\Omega}^{-1} X \bigg)}^{-1}  X^T \frac{1}{\gamma} \hat{\Omega}^{-1} X \boldsymbol{y}

=
{\big( X^T \hat{\Omega}^{-1} X \big)}^{-1}  X^T \hat{\Omega}^{-1} X \boldsymbol{y}
\tag{60}

　最後に未知量$\rho$を計算することでパラメータ同定が完了します。ただし、$\boldsymbol{\beta}$の計算に$\rho$が、$\rho$の計算に$\boldsymbol{\beta}$が必要になるので、数値計算で$\rho$が収束するまで追い込む必要があります。

　未知量$\rho$は(45)式のように定義されました。そこで(45)式を添え字$t$について並べてみると以下のようになります。

\varepsilon_2
=
\rho \varepsilon_1 + \mu_2
\tag{61}

\varepsilon_3
=
\rho \varepsilon_2 + \mu_3
\tag{62}

\vdots

\varepsilon_N
=
\rho \varepsilon_{N-1} + \mu_N
\tag{63}

(61)～(63)式の内容をベクトルにまとめてみましょう。以下のようにベクトル$\boldsymbol{\varepsilon}_{t+1}$、$\boldsymbol{\varepsilon}_t$、$\boldsymbol{\mu}$を定義します。

\boldsymbol{\varepsilon}_{t+1}
=
\begin{bmatrix}
\varepsilon_2 \\
\varepsilon_3 \\
\vdots \\
\varepsilon_N \\
\end{bmatrix}
\in \mathbb{R}^{N-1}
\tag{64}

\boldsymbol{\varepsilon}_t
=
\begin{bmatrix}
\varepsilon_1 \\
\varepsilon_2 \\
\vdots \\
\varepsilon_{N-1} \\
\end{bmatrix}
\in \mathbb{R}^{N-1}
\tag{65}

\boldsymbol{\mu}
=
\begin{bmatrix}
\mu_2 \\
\mu_3 \\
\vdots \\
\mu_N \\
\end{bmatrix}
\in \mathbb{R}^{N-1}
\tag{66}

すると(61)～(63)式は(1)式と同じ形式で以下のようにまとめることができます。

\boldsymbol{\varepsilon}_{t+1}
=
\boldsymbol{\varepsilon}_t \rho
+
\boldsymbol{\mu}
\tag{67}

ベクトル$\boldsymbol{\varepsilon}_{t+1}$、$\boldsymbol{\varepsilon}_t$はデータ、$\boldsymbol{\mu}$は独立誤差(と定義した)であるため、未知量$\rho$は最小二乗法の解(2)式より以下のように計算できます。

\rho
=
\big( {\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_t \big)^{-1} {\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_{t+1}

\therefore
\rho
=
\frac{{\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_{t+1}}{{\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_t}
\tag{68}

(68)式の$\boldsymbol{\varepsilon_t}$と$\boldsymbol{\varepsilon}_{t+1}$は$\boldsymbol{\varepsilon}$の要素で構成されています。さらに$\boldsymbol{\varepsilon}$は(1)式から得られる$\boldsymbol{\beta}$の関数であるため、(68)は$\rho$を$\boldsymbol{\beta}$から求める式だと考えられます。

4.2 数値計算のながれ　

　前節で紹介したとおり、一般化最小二乗法は数値計算によって解がもとまります。具体的には未知量$\rho$が収束するように追い込みます。本節では計算のながれについてざっと説明します。なお、計算の収束条件などは(イマイチ理解していないので)省略します。きっとMatlabとかを使えばうまく計算してくれるでしょう。

　数値計算は以下の手順で行います。
① $\rho=0$を仮定し、通常の最小二乗法(2)式を使って$\boldsymbol{\beta}$を求めます。

\boldsymbol{\beta}
=
{(X^T X)}^{-1} X^T \boldsymbol{y}
\tag{a}

② 求めた$\boldsymbol{\beta}$を使って(1)式から$\boldsymbol{\varepsilon}$を計算し、(64)～(68)式をもとに$\rho$を推定します。

\boldsymbol{\varepsilon}
=
\boldsymbol{y} - X \boldsymbol{\beta}
\tag{b}

\rho
=
\frac{{\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_{t+1}}{{\boldsymbol{\varepsilon}_t}^T \boldsymbol{\varepsilon}_t}
\tag{c}

③ 推定した$\rho$を使って(59)式のように行列$\hat{\Omega}$を構成します。

\hat{\Omega}(\rho)
= 
\begin{bmatrix}
1 & \rho  & \cdots & \rho^{N - 1} \\
\rho  & 1 & \cdots & \rho^{N - 2} \\
\vdots & \vdots & \ddots & \vdots\\
\rho^{N - 1}  & \rho^{N - 2}  & \cdots & 1\\
\end{bmatrix}
\tag{d}

④ (60)式を解くことで、再度パラメータ$\boldsymbol{\beta}$を推定します。

\boldsymbol{\beta}
=
{\big( X^T \hat{\Omega}^{-1} X \big)}^{-1}  X^T \hat{\Omega}^{-1} X \boldsymbol{y}
\tag{e}

⑤ 求めた$\boldsymbol{\beta}$を使って再度②の手順で$\rho$を推定します。あとは$\rho$の値が収束するまで②～④の手順を繰り返し、収束した際の$\boldsymbol{\beta}$が一般化最小二乗法の解となります。

5. おわりに

　今回の記事では、一般化最小二乗法によるモータのパラメータ同定の理論についてまとめました。線形代数の知識をフル活用したような内容で、これまでに学んだ数学が生かされてるなあと感じながら記事を書いていました。
　今回はしれっとスルーしましたが、平方根$\Psi$の具体的な成分表示なども付録という形でまとめられたらいいなと考えています。
　来週の記事では通常の最小二乗法と一般化最小二乗法のパラメータ推定結果の比較とかができればいいなと考えています。キャチロボ本番が近づいてきたのでそろそろロボットの制御とかやらないといけませんが、記事の毎週投稿は今後も続けていく予定です。
　今週も最後まで読んでいただきありがとうございました！