正規方程式(normal equation)は回帰式に対して最小二乗法(least square method)などを用いてパラメータの推定を行う際に出てくる数式です。$y = \beta_0 + \beta_1 x$のような単回帰の式の場合はシンプルな一方で、重回帰などを考えるにあたって行列表記を行うと難しそうに見えるので、当記事では行列に基づく正規方程式の導出の詳細を確認し、その直感的な理解に関して取りまとめを行いました。
内容の作成にあたっては「現代数理統計学（学術図書出版社）」の$11.1$節の「回帰モデル」と$11.2$節の「回帰モデルの推定」を主に参考にしました。
Contents
基本事項の確認
計画行列と回帰式
$$
\large
\begin{align}
y_i &= \beta_0 + \beta_1 x_i + \varepsilon_i \\
\varepsilon_i & \sim N(0,\sigma^2), \quad i.i.d.
\end{align}
$$
単回帰などを行うにあたっては、$x \in \mathbb{R}, y \in \mathbb{R}$のようなスカラー値を元に上記のように式を表す。同様に$k$次元の重回帰の式を考えると、下記のように表すことができる。
$$
\large
\begin{align}
y_i &= \beta_0 + \beta_1 x_{i1} + … \beta_k x_{ik} + \varepsilon_i \\
\varepsilon_i & \sim N(0,\sigma^2), \quad i.i.d.
\end{align}
$$
ここで上記の$2$式はどちらも$i$番目のサンプルに対応しており、さらに重回帰の式では複数変数を取り扱っており、要素を一つ一つ表すとなかなか表記が大変になる。
ここで表記をシンプルに行うにあたっては$n$対のサンプルに関する計画行列(design matrix)という行列を定義すると以下のように式の整理が行いやすい。
$$
\large
\begin{align}
y &= X \beta + \varepsilon \\
y &= \left(\begin{array}{c} y_{1} \\ \vdots \\ y_{n} \end{array} \right) \\
X &= \left(\begin{array}{cccc} 1 & x_{11} & \cdots & x_{1k} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & \cdots & x_{nk} \end{array} \right) \\
\beta &= \left(\begin{array}{c} \beta_{0} \\ \beta_{1} \\ \vdots \\ \beta_{k} \end{array} \right) \\
\varepsilon & \sim N_{n}(0,\sigma^2I_{n})
\end{align}
$$
ここで上記の$X$を計画行列(design matrix)、$\beta$を回帰係数ベクトルという。
残差ベクトル$e$と誤差の二乗和
「計画行列と回帰式」で表したような実際に得られた$y$に対して$X \beta$を用いて計算した予測値との差を考え、$e = y – X \beta$のように定義するとき、$e$は残差ベクトルという。
また、下記で表したように$e^{\mathrm{T}}e$は誤差の二乗和に一致する。
$$
\large
\begin{align}
e^{\mathrm{T}}e &= (y – X \beta)^{\mathrm{T}} (y – X \beta) \\
&= \sum_{i=1}^{n} (y_i – (\beta_{0} + \beta_{1}x_{i1} + … + \beta_{1}x_{ik}))^{2}
\end{align}
$$
正規分布の最尤推定＝最小二乗法
最小二乗法は「残差ベクトル$e$と誤差の二乗和」で確認した誤差二乗和の$e^{\mathrm{T}}e$が最小になるような$\beta$をパラメータの推定量と考える手法である。ここで「最小二乗法」によって導出される結果が「正規分布の最尤推定」の結果に一致することを以下に示す。
観測値を元に作成する計画行列$X \in \mathbb{R}^{n \times (k+1)}$と$X$に対応する$y \in \mathbb{R}^{n}$に対し、パラメータのベクトル$\beta$に関する尤度関数$L(\beta)$は$y$に関する同時確率密度関数を元に下記のように定義される。
$$
\large
\begin{align}
L(\beta) &= f(y|X,\beta,\sigma^2) \\
&= \frac{1}{(2 \pi \sigma^2)^{n/2}} \exp \left( -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) \right) \\
&= \exp \left( -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)} \right)
\end{align}
$$
上記に対して「$L(\beta)$が$\beta$に関して最大 $\iff$ $\log{L(\beta)}$が$\beta$に関して最大」である。このとき$\log{L(\beta)}$は下記のように表すことができる。
$$
\large
\begin{align}
\log{L(\beta)} &= \log \left[ \exp \left( -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)} \right) \right] \\
&= -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)}
\end{align}
$$
上記の$\beta$に関する最大化を考えることは$e^{\mathrm{T}}e = (y – X \beta)^{\mathrm{T}} (y – X \beta)$の$\beta$に関する最小化を考えることに一致することが式より確認できる。
よって、「最小二乗法」によって導出される結果が「正規分布の最尤推定」の結果に一致することが示される。
ベクトル・行列に関する微分の公式
・$a^{\mathrm{T}}b$の微分
$$
\large
\begin{align}
\left(\begin{array}{c} \frac{\partial}{\partial a_{1}} \\ \vdots \\ \frac{\partial}{\partial a_{q}} \end{array} \right) a^{\mathrm{T}}b &= \left(\begin{array}{c} \frac{\partial}{\partial a_{1}} \\ \vdots \\ \frac{\partial}{\partial a_{q}} \end{array} \right) \sum_{k=1}^{p} a_{k} b_{k} \\
&= \left(\begin{array}{c} b_{1} \\ \vdots \\ b_{q} \end{array} \right) = b
\end{align}
$$
・対称行列$C$に関して$a^{\mathrm{T}}Ca$の微分
$$
\large
\begin{align}
\left(\begin{array}{c} \frac{\partial}{\partial a_{1}} \\ \vdots \\ \frac{\partial}{\partial a_{q}} \end{array} \right) a^{\mathrm{T}}Ca &= 2 \left(\begin{array}{ccc} c_{11} & \cdots & c_{1q} \\ \vdots & \ddots & \vdots \\ c_{q1} & \cdots & c_{qq} \end{array} \right) \left(\begin{array}{c} a_{1} \\ \vdots \\ a_{q} \end{array} \right) \\
&= 2Ca
\end{align}
$$
詳しい導出は「現代数理統計学（学術図書出版社）」の章末課題$11.2$で取り扱った。
対称行列$X^{\mathrm{T}}X$の半正定値性
「特異値分解(SVD)」の「半正定値行列の固有値」で詳しく取り扱ったが、対称行列$X^{\mathrm{T}}X$は半正定値性を持つ。
正規方程式の導出
導出の詳細
「正規分布の最尤推定＝最小二乗法」の議論で、$\beta$に関する対数尤度$\log{L(\beta)}$は下記のように導出した。
$$
\large
\begin{align}
\log{L(\beta)} = -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)}
\end{align}
$$
上記は$\beta^{\mathrm{T}} X^{\mathrm{T}} y, y^{\mathrm{T}} X \beta$がどちらもスカラーであり$\beta^{\mathrm{T}} X^{\mathrm{T}} y = y^{\mathrm{T}} X \beta$が成立することに基づいて下記のように変形できる。
$$
\large
\begin{align}
\log{L(\beta)} &= -\frac{1}{2 \sigma^2}(y – X \beta)^{\mathrm{T}} (y – X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)} \\
&= -\frac{1}{2 \sigma^2}(y^{\mathrm{T}}y – y^{\mathrm{T}} X \beta – \beta^{\mathrm{T}}X^{\mathrm{T}}y + \beta^{\mathrm{T}} X^{\mathrm{T}} X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)} \\
&= -\frac{1}{2 \sigma^2}(y^{\mathrm{T}}y – 2 \beta^{\mathrm{T}}X^{\mathrm{T}}y + \beta^{\mathrm{T}} X^{\mathrm{T}} X \beta) – \frac{n}{2}\log{(2 \pi \sigma^2)}
\end{align}
$$
上記に対してベクトル$\displaystyle \beta = \left(\begin{array}{c} \beta_{0} \\ \vdots \\ \beta_{p} \end{array} \right)$での微分を考える。
$$
\large
\begin{align}
\left(\begin{array}{c} \frac{\partial}{\partial \beta_{0}} \\ \vdots \\ \frac{\partial}{\partial \beta_{p}} \end{array} \right) \log{L(\beta)} = 2 X^{\mathrm{T}} X \beta – 2 X^{\mathrm{T}}y
\end{align}
$$
上記が零ベクトルに一致する$\beta$を$\hat{\beta}$とおくと、$\hat{\beta}$は下記のように表すことができる。
$$
\large
\begin{align}
2 X^{\mathrm{T}} X \hat{\beta} – 2 X^{\mathrm{T}}y &= 0 \\
X^{\mathrm{T}} X \hat{\beta} &= X^{\mathrm{T}}y \\
\hat{\beta} &= (X^{\mathrm{T}} X)^{-1} X^{\mathrm{T}}y
\end{align}
$$
ここで$\hat{\beta}$の導出に用いた$X^{\mathrm{T}} X \beta = X^{\mathrm{T}}y$を正規方程式(normal equation)とよぶ。また、上記の$\hat{\beta}$は最尤推定量に一致する。
・参考
「現代数理統計学（学術図書出版社）」 章末課題$11.3$
正規方程式の直感的理解
正規方程式の直感的な理解を行うにあたって、$y_i = \beta_{0} + \beta_{1} x_i$の式に対し、そのまま偏微分を行うパターンと$y = X \beta$を用いるパターンの$2$通りで導出を行い、それぞれの式の確認を行う。
・$y_i = \beta_{0} + \beta_{1} x_i$の微分
$y_i = \beta_{0} + \beta_{1} x_i$に関しては誤差の関数$E(\beta_{0},\beta_{1})$は下記のように表される。
$$
\large
\begin{align}
E(\beta_{0},\beta_{1}) = \sum_{i=1}^{n} (y_{i} – (\beta_{0} + \beta_{1} x_i))^2
\end{align}
$$
ここで上記の$E(\beta_{0},\beta_{1})$に対し、$\beta_{0},\beta_{1}$で偏微分を行うとそれぞれ下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial E(\beta_{0},\beta_{1})}{\partial \beta_{0}} &= – 2 \sum_{i=1}^{n} (y_{i} – (\beta_{0} + \beta_{1} x_i)) \\
\frac{\partial E(\beta_{0},\beta_{1})}{\partial \beta_{1}} &= – 2 \sum_{i=1}^{n} (y_{i} – (\beta_{0} + \beta_{1} x_i)) x_i
\end{align}
$$
上記が$0$に一致すると考えるとき、それぞれの式は下記のように整理できる。
$$
\large
\begin{align}
\frac{\partial E(\beta_{0},\beta_{1})}{\partial \beta_{0}} &= 0 \\
– 2 \sum_{i=1}^{n} (y_{i} – (\beta_{0} + \beta_{1} x_i)) &= 0 \\
\sum_{i=1}^{n} (\beta_{0} + \beta_{1} x_i) &= \sum_{i=1}^{n} y_{i} \quad (1)
\end{align}
$$
$$
\large
\begin{align}
\frac{\partial E(\beta_{0},\beta_{1})}{\partial \beta_{1}} &= 0 \\
– 2 \sum_{i=1}^{n} (y_{i} – (\beta_{0} + \beta_{1} x_i)) x_i &= 0 \\
\sum_{i=1}^{n} (\beta_{0} x_i + \beta_{1} x_i^2) &= \sum_{i=1}^{n} x_{i}y_{i} \quad (2)
\end{align}
$$
・$y = X \beta$の微分
$y = X \beta$に関しては誤差の関数$E(\beta)$は下記のように表される。
$$
\large
\begin{align}
E(\beta) &= (y – X \beta)^{\mathrm{T}} (y – X \beta) \\
&= y^{\mathrm{T}}y – 2 \beta^{\mathrm{T}} X^{\mathrm{T}} y + \beta^{\mathrm{T}} X^{\mathrm{T}} X \beta
\end{align}
$$
ここで上記の$E(\beta)$に対し、ベクトル$\beta$で偏微分を行うと下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial}{\partial \beta} E(\beta) = 2 X^{\mathrm{T}} X \beta – 2 X^{\mathrm{T}} y
\end{align}
$$
上記が$0$に一致すると考えるとき、それぞれの式は下記のように整理できる。
$$
\large
\begin{align}
\frac{\partial}{\partial \beta} E(\beta) &= 0 \\
2 X^{\mathrm{T}} X \beta – 2 X^{\mathrm{T}} y &= 0 \\
X^{\mathrm{T}} X \beta &= X^{\mathrm{T}} y \quad (3)
\end{align}
$$
・式の対応
ここで$(1),(2)$式を$(3)$式の形式で表すことを考える。
$$
\large
\begin{align}
\left(\begin{array}{c} \displaystyle \sum_{i=1}^{n} (\beta_{0} + \beta_{1} x_i) \\ \displaystyle \sum_{i=1}^{n} (\beta_{0} x_i + \beta_{1} x_i^2) \end{array} \right) &= \left(\begin{array}{c} \displaystyle \sum_{i=1}^{n} y_{i} \\ \displaystyle \sum_{i=1}^{n} x_{i}y_{i} \end{array} \right) \\
\left(\begin{array}{c} \displaystyle n \beta_{0} + \beta_{1} \sum_{i=1}^{n} x_i \\ \displaystyle \beta_{0} \sum_{i=1}^{n} x_i + \beta_{1} \sum_{i=1}^{n} x_i^2 \end{array} \right) &= \left(\begin{array}{c} \displaystyle \sum_{i=1}^{n} y_{i} \\ \displaystyle \sum_{i=1}^{n} x_{i}y_{i} \end{array} \right) \\
\left(\begin{array}{c} n & \displaystyle \sum_{i=1}^{n} x_i \\ \displaystyle \sum_{i=1}^{n} x_i & \displaystyle \sum_{i=1}^{n} x_i^2 \end{array} \right) \left(\begin{array}{c} \beta_{0} \\ \beta_{1} \end{array} \right) &= \left(\begin{array}{ccc} 1 & \cdots & 1 \\ x_1 & \cdots & x_n \end{array} \right) \left(\begin{array}{c} y_{1} \\ \vdots \\ y_{n} \end{array} \right) \\
\left(\begin{array}{ccc} 1 & \cdots & 1 \\ x_1 & \cdots & x_n \end{array} \right) \left(\begin{array}{cc} 1 & x_1 \\ \vdots & \vdots \\ 1 & x_n \end{array} \right) \left(\begin{array}{c} \beta_{0} \\ \beta_{1} \end{array} \right) &= \left(\begin{array}{ccc} 1 & \cdots & 1 \\ x_1 & \cdots & x_n \end{array} \right) \left(\begin{array}{c} y_{1} \\ \vdots \\ y_{n} \end{array} \right) \\
X^{\mathrm{T}} X \beta &= X^{\mathrm{T}} y
\end{align}
$$
上記のように式の対応を考えることができる。よって、$X^{\mathrm{T}} X \beta = X^{\mathrm{T}} y$の形で表される正規方程式は直感的には左辺が予測に基づく式、右辺が実測値に基づく式であるとそれぞれ解釈しておくと良い。
・参考
最小二乗法・相関係数・決定係数
https://www.hello-statisticians.com/explain-terms-cat/regression1.html