回帰問題における最小二乗解とノルム最小解

はじめに

前回の記事（ノルム最小解の射影による導出） で、連立線形方程式の

• 最小二乗解
• ノルム最小解

についてまとめました。ここでは機械学習の分野でこれらが登場する場面として、線形回帰問題、およびその双対系であるカーネル回帰をとりあげます。
具体的には PRML の §3.1, §6.1 で登場する話題です。

線形回帰問題の定式化

基本的にPRMLと同じ表記を使います。

学習データとして、$N$ 個の入力データ $\{x_i\}$、出力データ $\{t_i\}$ の組がある。
入力データ $x$ をある関数セット $\phi_j\ (j=1,\dots,M)$ に通し、これを並べたもの
$$\phi(x) = \begin{pmatrix}
\phi_1(x) \\
\vdots \\
\phi_M(x)
\end{pmatrix} $$ を特徴量ベクトル¹と呼ぶことにします。関数セットのことは基底関数と呼びます。
線形回帰モデルでは、特徴量の線形和を出力とします。つまり、
$$\sum_{j=1}^M \phi_j(x_i) w_j \simeq t_i~~~(i=1, \cdots , N)$$ であるような係数 $w_j$ を学習します。係数をベクトルにしてまとめて書くと、
$$\Phi w \simeq t$$ と書けます。この行列 $\Phi$ は計画行列と呼ばれます。重要なので、きちんと書き下しておきます。
$$\Phi := \begin{pmatrix}
\phi(x_1)^T \\
\vdots \\
\phi(x_N)^T \\
\end{pmatrix}
=\begin{pmatrix}
\phi_1(x_1) &\cdots & \phi_M(x_1) \\
\vdots & & \vdots \\
\phi_1(x_N) &\cdots & \phi_M(x_N) \\
\end{pmatrix}$$ 計画行列は「データ数(N)×特徴量の数(M)」の行列であることをしっかり覚えておくことが重要です。

結局、線形回帰問題とは、$w$ についての「解けない」連立線形方程式
$$\Phi w = t$$ について（何らかの意味で）最良な解を求める問題であることが分かります。
以前の記事でまとめた最小二乗解、ノルム最小解（もしくは一般逆行列）がここで登場します。

最小二乗回帰

通常の回帰問題では、データ数 $N$ が多いため、特徴量ベクトル $\phi$ の次元 $M$ よりも大きいです。
このような場合は、最小二乗解を使えば二乗誤差を最小化する近似となります。
$$w = (\Phi^T\Phi)^{-1}\Phi^T t$$

最小二乗回帰では、右辺に出てくる $\Phi^T\Phi $ について解くことが計算の中心になります。
$$ \Phi^T\Phi = \sum_{i=1}^N \phi(x_i)\phi(x_i)^T $$ これは次数 $M$ の正方行列です。$M$ が小さい、すなわち、基底関数の数が少ないほど、少ない計算量で解けることが分かります。したがって、問題に応じて、そのような上手い基底関数を見つけることが重要になってきます。

多項式関数 $\phi_j(x) = x^j\ (j=0,1,\dots)$ を基底関数として最小二乗回帰を行う Python のコードは以下のようになります²。

# 多項式最小2乗回帰
class LeastSquare:
    # N: 多項式次数
    def __init__(self, N):
        self.N = N

    # データを多項式基底関数による表示に直す
    def __data_vec(self, x):
        return np.array([x**n for n in range(self.N + 1)])

    # 計画行列 X を返す
    def __design_mat(self, xv):
        return np.array([self.__data_vec(x) for x in xv])

    # 学習データを入力して係数wを決定する
    def fit(self, x, y):
        X = self.__design_mat(x)

        # solve (X^t X) w = X^t y        
        self.w = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, y))

    # 予測値を返す
    def predict(self, x):
        # 配列かどうかで分けている（ダサい）
        try:
            X = self.__design_mat(x)
        except TypeError:
            X = self.__data_vec(x)
        return np.dot(X, self.w)

入力データとして正弦波にガウス乱数を足したものを与えて、5次多項式で回帰した例を示します。青線がデータ生成に用いた正弦波、水色の点がそれに乱数を足した学習データ、赤線が回帰曲線です。

カーネル回帰

カーネル回帰では、特徴量ベクトル $\phi(x)$ の次元 $M$ がデータ数 $N$ よりも大きいことを仮定します。

$M>N$ のときは、方程式 $\Phi w = t$ は不定で、ノルム最小解
$$ w = \Phi^T(\Phi\Phi^T)^{-1} t $$ が使えます。ただし、ノルム最小解は二乗誤差を0とするため、回帰問題にそのまま適用すると過学習になってしまいます。そこで、(L2)正則化項を入れて、
$$ w = \Phi^T(\Phi\Phi^T + \lambda I)^{-1} t $$ を使います。こうすることで、誤差を一定以下にしつつ、$||w||^2$ を最小化する解が得られます³。

カーネル回帰では右辺に現れる行列 $\Phi\Phi^T$ について解くことが計算の中心になります。この行列はグラム行列と呼ばれ、その行列要素は特徴量ベクトルの内積です。
$$ K_{ij} := (\Phi\Phi^T)_{ij} = \phi(x_i)^T\phi(x_j) $$ これは次数 $N$ の正方行列です。$N$ が小さい、すなわち、データ数が少ないほど、少ない計算量で解けます、したがって、データ数が少ない場合にカーネル回帰は適用しやすくなります。

カーネルトリック

この記事の主題からは逸れるのですが、カーネル回帰では基底関数 $\phi(x)$ を明示的に作らずにすべての計算を済ませることができます。その概要は以下の通りです。

カーネル回帰では基底関数 $\phi(x)$ を与える代わりに、その内積 $\phi(x_i)^T\phi(x_j)$ を関数として与えます。
$$ k(x, x') := \phi(x)^T\phi(x') $$

以下に示すように、（驚くべきことに）基底関数そのものではなく、その内積さえ与えれば、カーネル回帰の解は決まってしまいます。

まず、グラム行列は、カーネル関数を使って
$$ K_{ij} = k(x_i, x_j) $$ と求められます。次に、$ w = \Phi^T a$ とおいて、$w$ の代わりに $a$ を求めるようにすると、
$$ a = (\Phi\Phi^T + \lambda I)^{-1} t = (K + \lambda I)^{-1}t $$ となります。$a$ を使うと、$x$ における回帰予測値は、
$$ \phi(x)^T w = \phi(x)^T \Phi^T a = \sum_{j=1}^N \phi(x)^T \phi(x_j)a_j = \sum_{j=1}^N k(x, x_j) a_j $$ となるので、$\phi(x)$ を使わずにカーネル関数だけで得られます。ただし、学習に使った入力データセット $\{x_i\}$ を保持しておく必要があります。

以上に述べたカーネル回帰を実装すると、以下のようになります。
カーネル関数として、ガウスカーネル $\exp(-\beta(x-x')^2)$ を使っています。

# ガウスカーネル回帰
class KernelRegression:
    # beta: ガウスカーネルパラメータ
    # lam : 正則化係数
    def __init__(self, beta, lam):
        self.beta = beta
        self.lam = lam

    # gauss kernel
    def __kernel(self, x1, x2):
        return np.exp(-self.beta * (x1 - x2)**2)

    # グラム行列を返す
    def __gram_mat(self, x):
        N = x.size
        K = np.fromfunction(lambda i, j: self.__kernel(x[i], x[j]), 
                            (N, N), dtype=int)
        return K

    # 学習データを入力して係数aを決定する
    def fit(self, x, y):
        # 学習データを保持する
        self.x = x

        # solve (K + lambda I) a = y
        K = self.__gram_mat(x)
        self.a = np.linalg.solve(K + self.lam * np.eye(x.size), y)

    # 予測値を返す
    def predict(self, x):
        # \sum_j k(x, x_j) a_j
        K = np.fromfunction(lambda i, j: self.__kernel(x[i], self.x[j]),
                            (x.size, self.x.size), dtype=int)
        return K.dot(self.a)

先ほどと同じ問題で、カーネル回帰を実行した例を示します。ガウスカーネルのパラメータを $\beta=1.0$, 正則化係数を $\lambda=0.1$ としました。

まとめ

数式やコードが増えてしまいましたが、この記事で強調したかったのは、具体的な計算方法ではなく、線形回帰問題を線形方程式の一般的な問題として捉えたときの、最小二乗解とカーネル法の解の間の関係を対比的に見ることです。

• 最小二乗回帰 : $N>M$ が前提、最小二乗解 $(\Phi^T\Phi)^{-1}\Phi^Tt$、$M$ 次正方行列 $\Phi^T\Phi$ が主役
• カーネル回帰 : $N<M$ が前提、ノルム最小解 $\Phi^T(\Phi\Phi^T)^{-1}t$、$N$ 次正方行列 $\Phi\Phi^T$ が主役

実は主成分分析とカーネル主成分分析にも似たような関係があるのですが、それについてはまた別の記事で書ければと思います。

---

1. ベクトルとスカラーは面倒くさいので特に区別せずに書いています。添え字のない$\phi$, $w$, $t$, $a$ はベクトルと思ってください。 ↩

2. コード上の変数名は数式と一致させていません。 ↩

3. 最小二乗解にも正則化項を入れることができますが、触れませんでした。ここで載せた正弦波の例のように、データ数に比べて特徴量ベクトルの次元が十分小さい場合、もともと表現力が小さいので、過学習になりづらいからです。一方、カーネル回帰の場合は、もともとが表現力が高く、学習データに完全にフィットする解を作れてしまうため、正則化が必須になります。 ↩