pyhaya’s diary

プログラミング、特にPythonについての記事を書きます。Djangoや機械学習などホットな話題をわかりやすく説明していきたいと思います。

機械学習を原理から理解する 線形回帰

最近始めた機械学習の理論的な話の3つめです。相変わらず線形の話が続きます。過去記事を下に貼っておきます。

pyhaya.hatenablog.com
pyhaya.hatenablog.com

今回は、前回までの分類アルゴリズムとは異なり、回帰アルゴリズムについて書きます。

回帰とは

回帰というのは、データセットがあったときにパラメータと目的変数の間に成り立つ関係式を見つける操作です。実験とかで結果を理論曲線でフィッティングするイメージです。

...もう少し具体的な例を出すと、例えば、パラメータとして、過去1カ月分の降水量と気温、天候が与えられたときにある特定の日に収穫されたリンゴの糖度(目的変数)を予測するとかなると回帰を使います。降水量、気温、天候というデータを組み合わせることでリンゴの糖度を表現するわけです。

線形回帰

線形という表現からもわかるように、回帰の式に一次関数を使うのが線形回帰です。例えばパラメータが1つあって(xと書く)、そこからある変数の値(yと書く)を予測したいとします。いま、学習のためにいくつかデータが与えられていて、下のグラフのように分布しているとします。
f:id:pyhaya:20181129000333p:plain
これはおおむね線形に分布していることがわかります。これを一変数の一次関数で回帰したのがオレンジの直線になります。

回帰直線を求める

パラメータ

線形回帰をする際に求める必要があるのは次の2種類のパラメータです。一つは重みパラメータ \mathbf{w}そしてもう一つがバイアスbです。回帰に使うパラメータが1種類の時には重みパラメータは2次元平面上の直線の傾きを表します。


h(\mathbf{x}) = \mathbf{w}\cdot\mathbf{x}+b \tag{1}

ここで\mathbf{x}は与えられている学習データで、h(\mathbf{x})がそのデータから計算した目的変数の予測値です。ふつうは、バイアスも重みパラメータの中に取り込んだ下のような形を使います。


h(\mathbf{x}) = \mathbf{w}\cdot\mathbf{x} \tag{2}

ここで\mathbf{w}=(w_1, w_2,\cdots, b), \mathbf{x}=(x_1, x_2,\cdots, 1)

損失関数

回帰の際には、どれくらいうまく回帰ができているかを知るために損失関数を定義します。回帰によってこの損失関数を最小化することが目標になります。

線形回帰の際によく用いられる損失関数は二乗誤差です。二乗誤差は下のように定義されます。

\displaystyle
L_S(h)=\frac{1}{n}\sum_{i=1}^n(h(\mathbf{x_i})-y_i)^2 \tag{3}

nは学習データの数です。

損失関数を最小化する

損失関数をなるべく小さくするにはどのようにすればよいでしょうか?微分して0になる値を探します。

\displaystyle
\frac{d}{d\mathbf{x}}L_S(h)=\frac{2}{n}\sum_{i=1}^n(\mathbf{w}\cdot\mathbf{x_i}-y_i)\mathbf{x_i} =0 \tag{4}

ここで、2つ変数を導入します。

\displaystyle
A =\sum_{i=1}^n \mathbf{x_i}\mathbf{x_i}^T,\ \ \mathbf{b}=\sum_{i=1}^n y_i\mathbf{x_i} \tag{5}

Aは行列です。これを使うと、(4)式は下のように簡単に書くことができます。

 A\mathbf{w}=\mathbf{b} \tag{6}

つまり、パラメータはA逆行列が求まれば、


\mathbf{w}=A^{-1}\mathbf{b}\tag{7}

と求めることができます。

Pythonで書く

上の計算をPythonのコードに書き下すと下のようになります。

import matplotlib.pyplot as plt
import numpy as np

from sklearn.model_selection import train_test_split
from sklearn.datasets import make_regression


def LinearRegression(X, y):
    A = np.zeros((X.shape[1]+1, X.shape[1] + 1))
    b = np.zeros((X.shape[1]+1, 1))
    for i in range(len(X)):
        x = list(X[i])
        x.append(1)
        x = np.array(x).reshape((2, 1))
        A += x.dot(x.T)
        b += y[i] * x

    return np.linalg.inv(A).dot(b)


if __name__ == '__main__':
    X, y = make_regression(
        n_samples=50,
        n_features=1,
        noise=7,
    )

    w = LinearRegression(X, y)
    fit_x = np.linspace(-2, 2, 500)
    fit_y = fit_x * w[0, 0] + w[1, 0]
    plt.figure()
    plt.plot(X, y, 'o')
    plt.plot(fit_x, fit_y)
    plt.show()

9, 10行目でX.shape[1]+1のように+1しているのは、バイアスを重みパラメータの中に取り込むために、配列の要素数を1つ増やすためです。

このコードでは、簡単のために変数は1種類にしています。そのため計算結果は簡単に図示することができて、
f:id:pyhaya:20181129011247p:plain
うまくいっていることが見て取れます。

逆行列が存在しないとき

ここからは、少しだけ発展的な内容になります。

このアルゴリズムでは、計算に逆行列を使っています。しかし、逆行列は必ずしも存在するとは限りません。このようなときにはどのように計算すればよいか考えます。

逆行列が存在しないときというのは、線形代数の知識を使えば固有値に0を含むときであることがわかります。このようなときに行列Aを


A = VDV^\dagger \tag{8}

のように直交行列V固有値を対角成分にもつ対角行列Dを使って分解することができます。Dの対角成分に0が含まれなければ、Dのすべて対角成分をすべて逆数にした行列D^+がDの逆行列になり、実際、


VDV^\dagger VD^+V^\dagger=VDD^+V\dagger=VV^\dagger=1 \tag{9}

のように、真ん中から次々に単位行列に代わっていきます(直交行列Vは転置して複素共役をとったV^\daggerと掛けると単位行列になる性質を利用しています。 Wiki参照)。

対角成分に0が含まれていては、このような操作は行うことができません。この時にはD^+を次のように定義します。


D^+_{ii} = \begin{cases}1/D_{ii}\ \ \ \ (D_{ii}\neq 0)\\ 0 \ \ \ \ (\mathrm{else})\end{cases} \tag{10}

要するに、逆数とれるならとって、とれないなら0のままにしておくということです。このようなD^+によって計算されるA

A^+=VD^+V^\dagger \tag{11}

として逆行列の代わりに使えば、

\displaystyle
A\mathbf{w}=A(A^+\mathbf{b})=VDV^\dagger(VD^+V^\dagger)\mathbf{b} = \sum_{i:D_{ii}\neq0} \mathbf{v}_i\mathbf{v}_i^T\mathbf{b} \tag{12}

となります。A\mathbf{v}V中の列ベクトルです。最右辺を見てみると、A\mathbf{w}D_{ii}\neq 0であるようなiについての\mathbf{v}_iによって張られる空間への射影になっていることがわかります。このように空間を制限して考え直すことで、逆行列が存在しない場合においても回帰計算をすることが可能になります。