https://tsujimotter.hatenablog.com/entry/geometry-of-least-squares-method

「最小二乗法」とは、データ点をうまく表現するような直線を求めるための方法です。

統計学における一手法として計算的な側面が紹介されることが多いですが、その背景には実は幾何学的な面白い解釈があるという話を紹介したいと思います。

世の中には相関があるデータが山ほどあります。

ここで一つ、なにやら相関のありそうな変数 $x, y$ について、 $N$ 点データをとってみてグラフに表してみたときに、こんな関係があったとします。

このようなデータに対しては、データ点に近しいところを通る直線 $y = ax + b$ を「えいっ」と引きたくなりますよね。

ところで、この直線 $y = ax + b$ の係数 $a, b$ としては、どのようなものを選ぶのが妥当なのでしょうか？

ちなみに、エクセルさんはとても賢くて、散布図のデータ点を右クリックして、「近似曲線を追加」を選ぶと一発で最適な直線を引いてくれます。

今回の話は、このエクセルさんの背景にある数学のお話しです。

最小二乗法とは

妥当な $a, b$ の決定方法として、「誤差の二乗の総和を最小化」するというのが一つの有望な考え方です。

すなわち、直線とデータ点 $(x_1, y_1), \; (x_2, y_2), \; \ldots, \; (x_N, y_N)$ の差として

$\varepsilon_1 = ax_1 + b - y_1$
$\varepsilon_2 = ax_2 + b - y_2$

$\;\;\;\;\;\; \vdots$

$\varepsilon_N = ax_N + b - y_N$

という値（偏差）を考えます。図で表すと、緑色の波線の長さを正負の符号込みで考えていることになります：

この $\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_N$ たちの二乗の総和（の1/2）

$\displaystyle f(a, b) = \frac{1}{2}\sum_{i=1}^{N} \varepsilon_i^2 \tag{1}$

を最小化するような $a, b$ を妥当な係数としましょうという方針です。

二乗の総和を最小化するので、最小二乗法といいます。

なんで二乗をとるのかということについて簡単に触れておきます。単に $\varepsilon_i$ を足した $\sum_{i=1}^{N} \varepsilon_i$ を計算すると、偏差の $\pm$ が相殺されてほとんど消えてしまうからです。
では、絶対値でもいいじゃんと思うわけですが、二乗にしておくことで、後に微分するときに計算が簡単になるのです。

なぜ $1/2$ をするのかということについても、後で微分するときの事情です。

さて、式 $(1)$ の $f(a, b)$ を最小化するような $a, b$ を求めよという問題を考えましょう。

答えは簡単で、 $f$ を $a, b$ でそれぞれ偏微分して、 $= 0$ となるような $(a, b)$ の組を求めればよいです。

まず、 $f$ を $b$ で偏微分してみます。

$\displaystyle \begin{align*} \frac{\partial f}{\partial b} &= \frac{1}{2}\sum_{i=1}^{N} \frac{\partial }{\partial b} (ax_i + b - y_i)^2 \\ &= \frac{1}{2}\sum_{i=1}^{N} 2(ax_i + b - y_i) \\ &= b N + a \sum_{i=1}^N x_i - \sum_{i=1}^{N} y_i \end{align*}$

次に、 $f$ を $a$ で偏微分します。

$\displaystyle \begin{align*} \frac{\partial f}{\partial a} &= \frac{1}{2}\sum_{i=1}^{N} \frac{\partial }{\partial a} (ax_i + b - y_i)^2 \\ &= \frac{1}{2}\sum_{i=1}^{N} 2(ax_i + b - y_i) \cdot x_i \\ &= b \sum_{i=1}^N x_i + a \sum_{i=1}^N x_i^2 - \sum_{i=1}^{N} x_i y_i \end{align*}$

よって、 $\displaystyle \frac{\partial f}{\partial b} = 0, \; \frac{\partial f}{\partial a} = 0$ とすると

$\begin{cases} \displaystyle b N + a \sum_{i=1}^N x_i - \sum_{i=1}^{N} y_i = 0 \\ \displaystyle b \sum_{i=1}^N x_i + a \sum_{i=1}^N x_i^2 - \sum_{i=1}^{N} x_i y_i = 0 \end{cases} \tag{3}$

なる連立方程式が得られます。

一見複雑に見えますが、係数がややこしいだけで、ただの $(b, a)$ を変数とする連立一次方程式です。

行列で表すと、次のようになります：

$\displaystyle \begin{pmatrix} N & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i & \sum_{i=1}^N x_i^2 \end{pmatrix} \begin{pmatrix} b \\ a \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{N} y_i \\ \sum_{i=1}^{N} x_i y_i \end{pmatrix} \tag{4}$

この連立方程式を正規方程式といいます。何で「正規」なのかは後で分かります。

あとは連立方程式を解くだけなので、ここまでで最小二乗法は完了ということになりますが、係数行列の部分はいったい何の行列なのかということが疑問に湧くと思います。この辺をもう少し深掘りしてみましょう。

行列を使って解く

正規方程式の謎を解くには、最小二乗法をベクトルと行列を用いて表現し直すことが鍵です。

以下の4つの行列を用意します：

$\displaystyle \boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_N \end{pmatrix}, \;\;\;\; \mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix},$
$\displaystyle X = \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_N \end{pmatrix}, \;\;\;\; \displaystyle \boldsymbol{\theta} = \begin{pmatrix} b \\ a \end{pmatrix}$

すると、最初の偏差の定義 $\varepsilon_i = ax_i + b - y_i$ は、まとめて

$\boldsymbol{\varepsilon} = X\boldsymbol{\theta} - \mathbf{y} \tag{5}$

と表せます。 $X\boldsymbol{\theta}$ の部分が少しトリッキーですが、行列の積をよくよく考えると分かります。

$X\boldsymbol{\theta} = \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots \\ 1 & x_N \end{pmatrix} \begin{pmatrix} b \\ a \end{pmatrix} = \begin{pmatrix} ax_1 + b \\ ax_2 + b \\ \vdots \\ ax_N + b \end{pmatrix}$

これによって、最小二乗法はベクトルのノルムを用いて

$\displaystyle f(\boldsymbol{\theta}) = \frac{1}{2}\|\boldsymbol{\varepsilon}\|^2 = \frac{1}{2}\|X\boldsymbol{\theta} - \mathbf{y}\|^2 \tag{6}$

を最小化する $\boldsymbol{\theta} = (b, a)$ を求める問題だと言い換えられます。

これを $b, a$ で偏微分すればよいわけですが、まとめて

$\displaystyle \nabla_{\boldsymbol{\theta}} = \begin{pmatrix} \dfrac{\partial }{\partial b} \\ \dfrac{\partial}{\partial a} \end{pmatrix}$

という演算子を考えます（これをナブラ演算子といいます）。

ナブラ演算子を $f(\boldsymbol{\theta})$ に適用して

$\displaystyle \nabla_{\boldsymbol{\theta}} f(\boldsymbol{\theta}) = \begin{pmatrix} \dfrac{\partial f}{\partial b} \\ \dfrac{\partial f}{\partial a} \end{pmatrix} = \mathbf{0}$

を解けば最適な係数 $\boldsymbol{\theta} = (b, a)$ が得られるという寸法です。

ここでベクトルのノルムは内積によって次のように変形できることに注意します：

$\displaystyle f(\boldsymbol{\theta}) = \frac{1}{2}\|X\boldsymbol{\theta} - \mathbf{y}\|^2 = \frac{1}{2}(\boldsymbol{\theta}^T X^T - \mathbf{y}^T)(X\boldsymbol{\theta} - \mathbf{y})$

これを用いて以下のように式変形します。

$\begin{align*} f(\boldsymbol{\theta}) &= \frac{1}{2}(\boldsymbol{\theta}^T X^T - \mathbf{y}^T)(X\boldsymbol{\theta} - \mathbf{y}) \\ &= \frac{1}{2}\left\{ \boldsymbol{\theta}^T X^T X \boldsymbol{\theta} - \boldsymbol{\theta}^T X^T\mathbf{y} - \mathbf{y}^T X\boldsymbol{\theta} + \mathbf{y}^T\mathbf{y} \right\} \\ &= \frac{1}{2}\left\{ \boldsymbol{\theta}^T X^T X \boldsymbol{\theta} - 2\boldsymbol{\theta}^T X^T\mathbf{y} + \mathbf{y}^T\mathbf{y} \right\} \end{align*}$

最後の式変形は、 $\boldsymbol{\theta}^T X^T\mathbf{y}, \;\; \mathbf{y}^T X\boldsymbol{\theta}$ はスカラーであり、スカラーは転置操作に対して不変であることから分かります。

この最後の式に対して、ナブラによる微分を計算するには、次のような公式が活用できます：

$\nabla_{\boldsymbol{\theta}}(c) = \mathbf{0}$ （ $c$ は定数）
$\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta}^T \mathbf{a}) = \mathbf{a}$
$\nabla_{\boldsymbol{\theta}}(\boldsymbol{\theta}^T A \boldsymbol{\theta}) = (A + A^T) \boldsymbol{\theta}$

$X^T X$ は対称行列（ $(X^T X)^T = X^T X$ ）であることから

$\begin{align*} \nabla_{\boldsymbol{\theta}} f(\boldsymbol{\theta}) &= \frac{1}{2} \{ X^T X + (X^T X)^{T} \} \boldsymbol{\theta} - X^T \mathbf{y} \\ &= X^T X \boldsymbol{\theta} - X^T \mathbf{y} \end{align*} \tag{7}$

が得られます。

ここで $\nabla_{\boldsymbol{\theta}} f(\boldsymbol{\theta}) = \mathbf{0}$ から次の式が得られます：

$X^T X \boldsymbol{\theta} = X^T \mathbf{y} \tag{8}$

実は、この式 $(8)$ こそが正規方程式です。

実際、次のように計算できることから、式 $(4)$ が復元されます。

$\displaystyle X^T X = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_N \end{pmatrix} \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_N \end{pmatrix} = \begin{pmatrix} N & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i & \sum_{i=1}^N x_i^2 \end{pmatrix}$

$\displaystyle X^T \mathbf{y} = \begin{pmatrix} 1 & 1 & \cdots & 1 \\ x_1 & x_2 & \cdots & x_N \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_N \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{N} y_i \\ \sum_{i=1}^{N} x_i y_i \end{pmatrix}$

正規方程式はただの2元連立一次方程式なので、普通に解いたらよいのですが、 $X^T X$ が正則行列であれば（逆行列を持てば）

$\boldsymbol{\theta} = (X^T X)^{-1} X^T \mathbf{y}$

のように直接的に解を表示することができます。

幾何的な解釈

正規方程式

$X^T X \boldsymbol{\theta} = X^T \mathbf{y} \tag{8再掲}$

の意味は一体何なのでしょうか。ここからが今回の本題です。

そもそもこの方程式の解は、二乗誤差を最小化する係数ベクトル $\boldsymbol{\theta} = \begin{pmatrix} b \\ a \end{pmatrix}$ なので、これを $\hat{\boldsymbol{\theta}}$ と置くことにしましょう。

$X^T X \hat{\boldsymbol{\theta}} = X^T \mathbf{y}$

右辺を移項させると

$X^T X \hat{\boldsymbol{\theta}} - X^T \mathbf{y} = \mathbf{0}$

となります。左辺を $X^T$ で括ると

$X^T (X \hat{\boldsymbol{\theta}} - \mathbf{y}) = \mathbf{0} \tag{9}$

となります。実はこの等式は幾何的に解釈できます。

$X\boldsymbol{\theta} - \mathbf{y}$ は偏差を表すベクトルなのでした。

ここで $\mathbf{y}$ はデータ点の $y$ 座標、 $X\boldsymbol{\theta}$ は対応する直線の点の $y$ 座標がなす $N$ 次元ベクトルでした。

$X\boldsymbol{\theta}$ の図形的な解釈のために、行列を展開して考えます：

$\begin{align*} X\boldsymbol{\theta} &= \begin{pmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_N \end{pmatrix} \begin{pmatrix} b \\ a \end{pmatrix} \\ &= b \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix} + a \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix} \\ &= b \mathbf{x}_1 + a \mathbf{x}_2 \end{align*}$

ここで、 $\mathbf{x}_1 = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix}, \; \mathbf{x}_2 = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix}$ は、 $X$ の縦ベクトルです。

上の式から、 $X\boldsymbol{\theta}$ の全体は、（ $N$ 次元空間の中で）縦ベクトル $\mathbf{x}_1, \; \mathbf{x}_2$ が張る平面だということができます。ベクトル $\mathbf{x}_1, \; \mathbf{x}_2$ はこの平面の基底ベクトルです。

言い換えると、この平面は（正解データの $x$ 座標を固定して）パラメータ $\boldsymbol{\theta} = (b, a)$ を動かしたときに「直線モデルの $y$ 座標の組が実現できる領域」を表している平面だということになります。

モデルが実現できる平面なので、「モデル平面」と呼ぶことにしましょう。

このように解釈すると、式 $(9)$ の意味が見えてきます。 $X^T = \begin{pmatrix} \mathbf{x}_1^T \\ \mathbf{x}_2^T \end{pmatrix}$ とおくと、式 $(9)$ は

$\begin{pmatrix} \mathbf{x}_1^T \\ \mathbf{x}_2^T \end{pmatrix} (X \hat{\boldsymbol{\theta}} - \mathbf{y}) = \mathbf{0}$

と置き換えられます。これを成分ごとに分けると

$\begin{cases} \mathbf{x}_1^T (X \hat{\boldsymbol{\theta}} - \mathbf{y}) = 0 \\ \mathbf{x}_2^T (X \hat{\boldsymbol{\theta}} - \mathbf{y}) = 0 \end{cases}$

となります。

これは基底ベクトル $\mathbf{x}_1, \; \mathbf{x}_2$ と、偏差ベクトル $X \hat{\boldsymbol{\theta}} - \mathbf{y}$ の内積が、どちらも $0$ だということです。
つまり、これらのベクトルは直交しているということですね。

したがって、偏差ベクトル $X \hat{\boldsymbol{\theta}} - \mathbf{y}$ とモデル平面は直交することになります。

図で表すとこのようになります。

$\mathbf{y}$ は近似したい正解データを表すベクトルであり、モデル平面は（正解データの $x$ 座標を固定して）パラメータ $\boldsymbol{\theta} = \begin{pmatrix} b \\ a \end{pmatrix}$ を動かしたときに「直線の $y$ 座標の組が実現できる領域」を表している平面だということになります。

平面上の点で、正解データ $\mathbf{y}$ にもっとも「近い」点は、 $\mathbf{y}$ からモデル平面に下ろした垂線の足だと思うことができます。これこそが最小二乗法の解 $\hat{\boldsymbol{\theta}}$ による $X\hat{\boldsymbol{\theta}}$ だというわけです。