はじめに

今回は、私がツイログから最近掘り出してきて、勉強している黒木玄さんの最小二乗法についてのノートブックを紹介したいと思います。紹介といっても、私がこのノートブックを理解するのに使ったものを紹介するというほうが正しいですが、、、
線形代数をきちんと勉強している方は、黒木玄さんのノートブックだけで理解できるのだと思います（羨ましい）。数学的な解釈や間違い等あればお手数ですがどなたか指摘していただければと思います。

引用だけでは、私が記事を書く意味がないので最尤推定について少し触れておきます。最尤推定の特別な場合が最小二乗法にあたるということを簡単にですが示したいと思います。

最小二乗法は直交射影の言い換え

黒木玄さんの書いた丁寧な説明とJuliaの実装が以下から見れます。まずはこちらをご覧ください。今回触れるのは直交射影ですが、是非続きも読んでみてください。

こちらの画像は直交射影の説明の時に使われていたものです。上記ノートブックと併せてご活用ください。

ということで、ノートブックと画像見ても何のことかわからないよという、私のような方のためにいろいろ用意してきました。まずは画像を理解するところから始めましょう。画像の伝えたいことを式に表しているのがノートブックです。

イメージ（像、Im）とカーネル(核、Ker)

図をみてどこがイメージでどこがカーネルなのかを確認しましょう。

射影について雰囲気をつかむ

ある二つのベクトルの内積が0のとき、なにがいえるでしょうか？

イメージ（像、Im）について雰囲気をつかむ

例えばベクトル $v$ に行列 $X$ を掛けたときどうなるでしょうか？そうです、ベクトルが伸び縮みしたり、反転したりしますね。
今ベクトルだけを考えましたが、空間全体も同様に伸び縮みし、反転すると考えられますね。この変換後の空間のことを写像 $X$ のイメージといい、 $ImX$ と書きます。以下のリンクの動画がわかりやすいです。

カーネル（核、Ker）について雰囲気をつかむ

例えばこんなことを考えてみましょう。ある空間上に平行四辺形があります。これを写像で写すと伸び縮みしたり反転したりするといったことは先ほど学びました。でも、写像によっては伸びたり反転したりではなく、「つぶれる」といったこともあります。平行四辺形の場合でいうと写像で写した後のイメージは直線になるといった感じです。下の記事をご覧ください。平行四辺形ではないですが、四角形が直線上につぶれている様子が見れます。

このとき写像 $X$ で写した点が全て原点に集まるような集合があります。そのような集合をカーネルといいます。上記リンクの例で考えてみましょう。このとき写像は $A$ になります。

A = 
\begin{pmatrix}
2 & 3 \\
1 & 1.5
\end{pmatrix} \quad
\boldsymbol{x} = 
\begin{pmatrix}
x_{1} \\
x_{2}
\end{pmatrix}

カーネルとは写した点が全て原点に集まるような集合のことでしたね。つまり $Ax = 0$ となるような $x$ はどのようなものがあるかを考えればいいわけです。例えば、$(x_{1}, x_{2}) = (0, 0)$ ですね。あとは、$(x_{1}, x_{2}) = (3, -2)$ も条件を満たします。以上のことは実際に計算を行ってみればわかると思います。

ここまでの説明を動画にしてあるのが以下のリンクです。

$$$$

まとめ

ここで直交射影の画像に戻ります。ここまでくれば割と画像内のことがわかってくると思います。

まず、$(X^{T}X)^{-1}$ が存在する（しない）場合というのは、イメージがつぶれる（つぶれない）くらいに考えてもらっていいと思います。

ちゃんと考えるなら、$(X^{T}X)^{-1}$ が存在するということは、$X^{T}X$ が正則であるということ。正則であるということは、$X^{T}X$ の列ベクトルは一次独立であるということ。$X^{T}X$ の列ベクトルが一次独立であるならば、$X$ の列ベクトルもまた一次独立である。$X$ が一次独立であるならば、$X$ のランクは次元数と等しくなり、また $ImX$ のランクとも等しくなる。（次元定理/階数・退化次数の定理とか言うやつ）つまりカーネルが、0ベクトルしか持たないということになる。つまりイメージがつぶれないということになる。（自信ない）

また、$R^{r}$ は列ベクトルの次元数で、ここでは $\beta \in R^{2}$ ということになります。写像 $X$ は $X \in R^{3 \times 2}$ になります。$R^{2}$ の列ベクトルが横に3つ並んでいると考えてください。そうなると $X$ と $\beta$ の積は $X\beta \in R^{3}$ になりますね。 $R$ についてはそんな感じに考えてください。

最後に、$\beta$ が赤い直線上のどこにいてもいいということについて説明して、画像についての解説を終わります。直線上のどこにいてもいいということは、解がいくつもあるということです。これは、$(X^{T}X)^{-1}$ が存在しない、つまり$X$ が一次従属であるということからわかります。
じゃあ、これは解けないのかという話なんですがそうではなく、$L_{2}$ノルムを最小化するような$\beta$ など条件を増やすみたいです。
以下のリンクに逆行列が存在しない場合の話が書いています。詳しくはこちらをご覧ください。

さて、長くなりましたがこれで、なんで平面から直線になっているの？とか $ImX$ への直交射影ってどこ？とかがわかるようになったと思います。雰囲気つかんだらぜひノートブックのほうに戻ってみてください。さっきとはまた違ったように見えると思います。

最小二乗法は正規分布を用いた最尤推定

ではここからは線形代数はお休みです。ここからは一変数で話を進めていきたいと思います。
次は最小二乗法が最尤推定の特別な場合であるということを示します。最尤推定を勉強するときに必ず出てくる話題ですね。

みんなが知っている残差を使った求め方

観測値 $\boldsymbol{y} = [y_{1},,,y_{n}]$ が与えられているとする。

この観測値をよく表すように直線を引く。そしてその直線と各観測値との誤差、二乗誤差を最小化するというのが最小二乗法ですね。
今回帰直線を以下のようにあらわすとき、二乗誤差は次のように表されます。二乗誤差の係数は計算を行いやすくするためのものですね。

\begin{align}
回帰直線：y_{n} &= ax_{n} + b　\\
二乗誤差：J(a, b) &= -\frac{1}{2} \Sigma_{n=1}^{N} (y_{n} -ax_{n} - b)^ 2 
\end{align}

さて、ここまでの話が分からない方は下の解説が詳しいです。

最尤推定

次のような回帰式を考えます。ここで$\epsilon$ は正規分布で、各 $\epsilon$ は独立に取り出されます。

y_{n} = ax_{n} + b + \epsilon_{n} \quad \epsilon ～ iid N(0, \sigma^{2}) \tag{1}

ここで $y_{n}$ が得られる確率は次のとおりです。

\begin{align}
y_{n} &= ax_{n} + b + \epsilon_{n} \\
\epsilon_{n} &= y_{n} -ax_{n} - b \\
\end{align}

\begin{align}
P(y_{n}|\theta) \Rightarrow 
\dfrac{1}{\sqrt{2\pi\sigma^{2}}}\exp(-\dfrac{(y_{n} -ax_{n} - b)^ 2}{2\sigma^ 2}) \tag{2}
\end{align}

$\theta = [a, x_{n}, b, \sigma^{2}]$ です。最後の矢印は以下の式を利用しています。

\begin{align}
P(\epsilon_{n}|\sigma^2) \Rightarrow 
\dfrac{1}{\sqrt{2\pi\sigma^{2}}}\exp(-\dfrac{\epsilon_{n}^{2}}{2\sigma^ 2}) \\
\end{align}

今、観測値 $\boldsymbol{y} = [y_{1},,,y_{n}]$ が与えられているとする。説明変数 $\boldsymbol{x} = [x_{1},,,x_{n} ]$ を使って観測値 $y$ が得られる確率は次のようになります。

\begin{align}
P(\boldsymbol{y}|\theta) &= P(y_{1}|\theta) \times P(y_{2}|\theta) \times ... \times P(y_{n}|\theta) \\
&= \Pi_{n=1}^{N} P(y_{n}|\theta) \\
&= \Pi_{n=1}^{N} \dfrac{1}{\sqrt{2\pi\sigma^{2}}}\exp(-\dfrac{(y_{n} -ax_{n} - b)^ 2}{2\sigma^ 2}) \\
\end{align}

得られた確率 $P(\boldsymbol{y}|\theta)$ が最大になるようなパラメータを求めます。通常、対数を取って計算しやすくします。

\begin{align}
L(\theta|\boldsymbol{y}) &= \Sigma_{n=1}^{N} \log( P(y_{n}|\theta)) \\
&= -\frac{N}{2} \log(2\pi) -\frac{N}{2} \log(\sigma^{2}) - \frac{1}{2\sigma^{2}} \Sigma_{n=1}^{N} (y_{n} -ax_{n} - b)^ 2  \\
\end{align}

ここで、指数部分であった第三項に注目してください。この項は $\sigma^{2}$ を定数とみると、二乗誤差の項と同じですね。なぜ定数とみれるかというと、パラメータの$a, b$ を求めるとき偏微分を行いますが、微分対象以外の変数はすべて定数になるからです。

\begin{align}
最小二乗法&：-\frac{1}{2} \Sigma_{n=1}^{N} (y_{n} -ax_{n} - b)^ 2 \\
最尤推定&：- \frac{1}{2\sigma^{2}} \Sigma_{n=1}^{N} (y_{n} -ax_{n} - b)^ 2  
\end{align}

まとめ

よって最小二乗法は最尤推定でノイズを正規分布としたモデルを設計したものと同じということになります。なので、最小二乗法を使っているときは暗に誤差を正規分布と仮定していることになります。

ちなみに線形代数から最小二乗法を求めることを最初にしました。あの場合も暗に正規分布が仮定されています。通常 $Ax=b$ の形は解をもつとは言えないですが、誤差を許すことによってその誤差を最小化するような $x$ を求めようとしています。$||Ax-b||^{2}$ を最小化するということは暗に誤差の存在を認めていることになります。この話についてもこちらのリンクがわかりやすく説明してくれています。

Juliaのバックスラッシュ

以下の式は一度試してもらいたいです。Aが正則じゃないことに注目

# Ay = b
f(x) = 2 + 3x + randn()
x = 0:0.1:10
A = hcat( f.(x).^0, f.(x) )
b = rand(length(x))
y = A \ b

> julia
julia> \.jlをコピペ

参考

• https://www.math.nagoya-u.ac.jp/~larsh/teaching/F2011_LA/lecture7.pdf