【python】numpyで最小二乗法を実装（線形、多項式、正則化など）

はじめに

　最小二乗法をnumpyで実装してみた。

　理論背景についてはこちらを参照（外部リンク）。
mathtrain.jp

PRMLの線形回帰モデル（線形基底関数モデル） from Yasunori Ozaki

やるべきこと

　最小二乗法（正確には線形基底関数モデルによる回帰）は目的変数を説明変数の線形結合で表現しようというアイデア。面倒くさいことをすっ飛ばして言うと、次の式を解いて重みベクトル $\vec{w}$ を求めれば良い。この $\vec{w}$ は誤差関数（誤差のニ乗和）を最小化するようなパラメータになっている。

${ \displaystyle A = \left( \begin{array}{ccccc} 1 & x_{11} & x_{12} & \ldots & x_{1n} \\ 1 & x_{21} & x_{22} & \ldots & x_{2n} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{m1} & x_{m2} & \ldots & x_{mn} \end{array} \right) }$

${ \displaystyle \vec{w}=(A^TA)^{-1}A^T\vec{y} }$
　 $A$ は $n$ 次元の横ベクトルを $m$ 個縦に並べたもの・・・で、よくsklearnに入力しているnumpy配列と変わらない。1列目がぜんぶ1になっているのはバイアス項といって、要するに目的変数が（x-yグラフのy軸と考えて）上にずれたり下にずれたりするのを表現するもの。

　この重みベクトル $\vec{w}$ をどう使うか？　お察しの通り、未知の入力を（バイアス項を足して）この重みと掛け算してぜんぶ加算してあげると、推定された目的変数が出てきます・・・。

　でもこれだと線形関数しか回帰できないので、多項式を使って非線形関数にも対応させてあげる。別に難しいことはなく、入力ベクトルの次元を増やし、そこにそういった関数に入力を入れた値を突っ込んでやるだけ。

${ \displaystyle A = \left( \begin{array}{ccccc} 1 & x_{11} & x_{12} & \ldots & x_{1n} & f_1(x_{11} ) & \ldots & f_2(x_{11} ) &\ldots\\ 1 & x_{21} & x_{22} & \ldots & x_{2n} & f_1(x_{21} ) & \ldots & f_2(x_{21} ) &\ldots\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ldots & \vdots &\ldots\\ 1 & x_{m1} & x_{m2} & \ldots & x_{mn} & f_1(x_{m1}) & \ldots & f_2(x_{m1}) &\ldots \end{array} \right) }$
　簡単そうで素晴らしい。また、正則化を行いたい場合は
${ \displaystyle \vec{w}=(A^TA+\lambda I)^{-1}A^T\vec{y} }$
　これでいい（L2正則化)。別にL2が正則化のすべてではないけど、一番実装が楽そうなのはこれ。ここで $I$ は単位行列（行列の対角線上がぜんぶ1の奴）、 $\lambda$ は正則化パラメータで数字を大きくすればするほど正則化が強く働く。適正値はパラメタチューニングして探す必要がある。

　L2正則化は過学習を防止し、汎化性能を高める効果がある。どうしてL2正則化で過学習を防げるのか？　定性的な説明で済ませるが、要するにこれを足すことで $\vec{w}$ の長さ（L2ノルム）が大きくなりすぎないようにできる。最小二乗法で教師データへのフィッティングを極限まで高めようとすると、とにかく複雑な関数をたくさん使って正確にフィットさせよう・・・！　となるので、結果的に $\vec{w}$ の長さが大きくなる。こうなることを抑制すると、逆にできるだけ単純な関数の和で行こう・・・ってなる。そんな感じ。

実装

　sklearn風に最小二乗法のクラスを作り、fitとpredictが呼べるようにすると使いやすい。まずLSM(Least Squares Method)クラスを作ろう。
　パラメータは色々あると思うが、とりあえず「基底関数どうするか（なし、多項式）」と「L2正則化するか」と「L2正則化の $\lambda$ をどれくらいにするか」だけ決められるようにしておく。基底関数に関しては、先に関数のリストを作っておいて後から便利に使う方針を取ることにする。多項式は10次多項式とする。本当はどれくらい次数を上げると良いのかもチューニングするべきだと思う・・・。

import numpy as np
from itertools import product

class LSM:
    def __init__(self, func=None, l2=False, l2_lambda=1):
        self.func = func
        self.func_list = self._gen_func_list(func)
        self.l2 = l2
        self.l2_lambda = l2_lambda
    
    def _gen_func_list(self, func):
        if func is None:
            return [np.vectorize(lambda x:x)]
        elif func == "poly":
            return [self._gen_poly(i) for i in range(1, 11)]
        else:
            raise Exception

    def _gen_poly(self, i):
        return np.vectorize(lambda x:np.power(x, i))

　クロージャ使って綺麗に書こうとしたら微妙に禍々しくなった。でもまあ、こうやって下ごしらえしておけば後は楽である。あとは無心でfitとpredictを作るだけである。コツはndarrayで処理しようとすると色々面倒くさいので、素直にnp.mat型で扱うこと。

class LSM:
    # 中略
    def fit(self, X, y):
        A = np.mat(np.hstack([np.c_[np.ones(X.shape[0])]] + 
                             [f(X) for f in self.func_list]))
        A_t = A.T
        y_mat = np.mat(y).T
        if self.l2:
            lambda_I = self.l2_lambda*np.mat(np.identity(A.shape[1]))
            self.w = ((((A_t*A) + lambda_I ).I)*A_t*y_mat)
        else:
            self.w = (((A_t*A).I)*A_t*y_mat)

    def predict(self, X):
        X = np.mat(np.hstack([np.c_[np.ones(X.shape[0])]] + 
                             [f(X) for f in self.func_list]))
        lst = []
        for x in X:
            lst.append((x*self.w)[0,0])
        return np.array(lst)

　こんな感じ。あとは試しにsin関数（を適当な区間で切り取って更に誤差を足したもの）でも回帰させてみる。汚いけど、テストに使ったプログラムをぜんぶ載せてみる。

# coding: UTF-8
import numpy as np
from sklearn.metrics import mean_squared_error as mse

import matplotlib.pyplot as plt

def rmse(y1, y2):
    return np.sqrt(mse(y1,y2))

def scatter(x_train, y_train, x_test, y_test, y_preds, filename):
    fig = plt.figure()
    
    ax = fig.add_subplot(1,1,1)
    ax.scatter(x_train, y_train, color="b", label="train data")
    ax.plot(x_test, y_test, color="r", label="true line")
    ax.plot(x_test, y_preds, color="g", label="predicted line")
    ax.set_xlabel("rmse:{0:.6f}".format(rmse(y_test, y_preds)))
    ax.set_ylabel("")
    plt.legend()
    plt.savefig(filename)

def gen_sin_data(a=1, b=0, n=20, v=5, x_range=20):
    X = np.array(sorted((np.random.rand(n) - np.array([0.5]*n))*np.array([x_range]*n)))
    Y = np.sin(X)*a + np.array([b]*n) + np.random.normal(0, v, n)
    return X, Y

class LSM:
    def __init__(self, func=None, l2=False, l2_lambda=1):
        self.func = func
        self.func_list = self._gen_func_list(func)
        self.l2 = l2
        self.l2_lambda = l2_lambda
    
    def _gen_func_list(self, func):
        if func is None:
            return [np.vectorize(lambda x:x)]
        elif func == "poly":
            return [self._gen_poly(i) for i in range(1, 11)]
        else:
            raise Exception

    def _gen_poly(self, i):
        return np.vectorize(lambda x:np.power(x, i))

    def fit(self, X, y):
        A = np.mat(np.hstack([np.c_[np.ones(X.shape[0])]] + 
                             [f(X) for f in self.func_list]))
        A_t = A.T
        y_mat = np.mat(y).T
        if self.l2:
            lambda_I = self.l2_lambda*np.mat(np.identity(A.shape[1]))
            self.w = ((((A_t*A) + lambda_I ).I)*A_t*y_mat)
        else:
            self.w = (((A_t*A).I)*A_t*y_mat)

    def predict(self, X):
        X = np.mat(np.hstack([np.c_[np.ones(X.shape[0])]] + 
                             [f(X) for f in self.func_list]))
        lst = []
        for x in X:
            lst.append((x*self.w)[0,0])
        return np.array(lst)

def main():
    X_train, y_train = gen_sin_data(n=50, v=0.4)
    X_test, y_test = gen_sin_data(n=5000, v=0.0)

    X_train_v = np.array(np.mat(X_train).T)
    X_test_v = np.array(np.mat(X_test).T)

    for l in [0, 0.01, 0.1, 1]:
        print("linear reg  lambda:", l)
        lsm = LSM()
        lsm.fit(X_train_v, y_train)
        y_preds = lsm.predict(X_test_v)
        print("rmse:", rmse(y_test, y_preds))
        scatter(X_train, y_train, X_test,  y_test, y_preds, "lin_{0:.2f}.png".format(l))

        print("poly reg  lambda:", l)
        lsm = LSM(func="poly", l2=True, l2_lambda=l)
        lsm.fit(X_train_v, y_train)
        y_preds = lsm.predict(X_test_v)
        print("rmse:", rmse(y_test, y_preds))
        scatter(X_train, y_train, X_test,  y_test, y_preds, "poly_{0:.2f}.png".format(l))
        
if __name__ == '__main__':
    main()

　実行すると次のようにRMSEが表示される。

linear reg  lambda: 0
rmse: 0.723655771604
poly reg  lambda: 0
rmse: 0.311310262135
linear reg  lambda: 0.01
rmse: 0.723655771604
poly reg  lambda: 0.01
rmse: 0.311303335228
linear reg  lambda: 0.1
rmse: 0.723655771604
poly reg  lambda: 0.1
rmse: 0.311251369924
linear reg  lambda: 1
rmse: 0.723655771604
poly reg  lambda: 1
rmse: 0.311606815936

　線形、多項式でそれぞれ一番いいのを見せてみる。