データ数が入力次元より少ない場合の重回帰分析

はじめに

この記事は古川研究室 Workout_calendar ４日目の記事です。
本記事は古川研究室の学生が学習の一環として書いたものです。

本記事では重回帰分析について述べます。
重回帰分析を扱った時にプログラムは回ったのに予測値がめちゃくちゃ！ってことはよくありますよね。その原因の一つに学習データの型(サイズ)があるかもしれません。重回帰では入力の次元数よりもデータ数が少ないと、出力値に大きな差が出ることがあります。
本記事では重回帰分析のアルゴリズムについて詳しく述べた後、入力が高次元かつデータ数が少ない場合に何が起こっているのか、またその解決法について述べます。

重回帰分析とは

問題設定

観測データとしてD次元のベクトルの入力$\mathbf x\in \mathbb R^D$と出力$y\in \mathbb R$のペア$\left\{(\mathbf{x}_n,y_n)\right\} _{n=1}^N $から$y=f(\bf{x})$となるような関数$f$を推定すること

重回帰分析を解く

\begin{eqnarray}
    f(\mathbf{x})=w^{(1)}x^{(1)}+w^{(2)}x^{(2)}+...+w^{(D)}x^{(D)} 
\end{eqnarray}

とモデルを定義します。

\begin{eqnarray}
    f(\mathbf{x})=\sum_{d=1}^{D}{w^{(d)}}{x^{(d)}} 
\end{eqnarray}

$\mathbf{w}=({w^{(1)}},{w^{(2)}},...,{w^{(D)}})^T,\mathbf{x}=(x^{(1)},x^{(2)},...x^{(D)})^T$とおくと

\begin{equation}
    f(\mathbf{x})=\mathbf{w}^T\mathbf{x} 
\end{equation}

観測データの出力とモデルの出力の二乗誤差の総和$J_\mathrm{LS}(\mathbf{w})$が最小となる$\mathbf{w}$を求めていきます。

\begin{align}
J_\mathrm{LS}(\mathbf{w})&=\sum_{n=1}^{N}(y_n-f(\mathbf{w}))^2\\
&=\sum_{n=1}^{N}(y_n-\mathbf{w}^T\mathbf{x}_n)^2\\
&=\sum_{n=1}^{N}(y_n-\mathbf{x}_n^T\mathbf{w})^2
\end{align}

\mathbf{y}=(y_1,y_2,...,y_N)^T,\mathbf{X}=\begin{pmatrix}x_1^{(1)} & x_1^{(2)} & \cdots & x_1^{(D)} \\
x_2^{(1)} & x_2^{(2)} & \cdots & x_2^{(D)} \\
\vdots & \vdots & \ddots & \vdots \\
x_N^{(1)} & x_N^{(2)} & \cdots & x_N^{(D)}&\end{pmatrix}

とおくと

\begin{eqnarray}
    J_\mathrm{LS}(\mathbf{w})&=&(\mathbf{y}-\mathbf{X}\mathbf{w})^T(\mathbf{y}-\mathbf{X}\mathbf{w})\\
    &=&(\mathbf{y}^T-\mathbf{w}^T\mathbf{X}^T)(\mathbf{y}-\mathbf{X}\mathbf{w})\\
    &=&\mathbf{y}^T\mathbf{y}-\mathbf{y}^T\mathbf{X}\mathbf{w}-\mathbf{w}^T\mathbf{X}^T\mathbf{y}+\mathbf{w}^T\mathbf{X}^T\mathbf{X}\mathbf{w}\\
    &=&\mathbf{y}^T\mathbf{y}-2\mathbf{w}^T\mathbf{X}^T\mathbf{y}+\mathbf{w}^T\mathbf{X}^T\mathbf{X}\mathbf{w}\\
\end{eqnarray}

上式を$\mathbf{w}$について偏微分し0とおくと

\begin{eqnarray}
    \frac{\partial}{\partial \mathbf{w}}J_\mathrm{LS}&=&\frac{\partial}{\partial \mathbf{w}}(\mathbf{y}^T\mathbf{y}-2\mathbf{w}^T\mathbf{X}^T\mathbf{y}+\mathbf{w}^T\mathbf{X}^T\mathbf{X}\mathbf{w})\\
    0&=&0-2\mathbf{X}^T\mathbf{y}+\{\mathbf{X}^T\mathbf{X}+(\mathbf{X}^T\mathbf{X})^T\}\mathbf{w}\\
    2\mathbf{X}^T\mathbf{X}\mathbf{w}&=&2\mathbf{X}^T\mathbf{y}\\
\mathbf{X}^T\mathbf{X}\mathbf{w}&=&\mathbf{X}^T\mathbf{y}\\
\mathbf{w}&=&(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}
\end{eqnarray}

となりパラメータ$\mathbf{w}$を求めることができました！

具体的な例を当てはめてみる

アイスクリーム屋の一日の売上をその日の気温と湿度より予測する重回帰問題を考えてみます。

# 必要なライブラリのimport
import numpy as np
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 入力データセット
N = 20
a1 = 1
a2 = 15
temp = np.random.uniform(22,35,N)
humi = np.random.uniform(0.0,0.9,N)
sales = a1*temp + a2*humi + np.random.rand(N)*2 - 1

# キャンバス作成
fig = plt.figure(figsize=(9, 9))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.set_title("example data: sales of ice")
ax.set_xlabel("humidity")
ax.set_ylabel("temperature")
ax.set_zlabel("sales of ice")
ax.scatter(temp, humi, sales,s=30, c='red')

# データセットを行列に
x_1 =temp [:,None]
x_2 =humi [:,None]
y = sales
X = np.concatenate([x_1,x_2],axis=1)

# 重回帰の解
w = np.linalg.inv(X.T@X)@X.T@y

new_samples_num = 30
new_temp,new_humi=np.meshgrid(np.linspace(22,35,new_samples_num),np.linspace(0.0,0.9))
new_temp=new_temp[:,:,None]
new_humi=new_humi[:,:,None]
new_X=np.concatenate([new_temp, new_humi], axis=2)
new_y=new_X[:,:,0]*w[0]+new_X[:,:,1]*w[1]

fig=plt.figure(1)
ax=Axes3D(fig)
ax.set_xlabel('temperature')
ax.set_ylabel('humidity')
ax.set_zlabel('sales')
ax.set_title('regression by linear model')
ax.scatter(X[:,0],X[:,1],y,c='red',label='trainig data $(x_{n1},x_{n2},y_n)$')
ax.plot_wireframe(new_X[:,:,0], new_X[:,:,1], new_y, color='b',label='estimated function $f_w(x)$')
plt.legend(fontsize=12)
plt.show()

気温と湿度からその日のアイスの売上を予測することができるようになりました！

入力の次元数がデータ数よりも多い場合

以下の重回帰モデルの解を見てみると、重回帰では$\mathbf{X}^T\mathbf{X}$が逆行列を持つ、すなわち正則であることが前提としてあることが分かります。

\begin{eqnarray}
\mathbf{w}&=&(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}
\end{eqnarray}

ここで入力の次元数$D$がデータ数$N$より多い場合を考えてみましょう。$(D>N)$
$\mathbf{X}$を特異値分解します。

\begin{eqnarray}
\mathbf{X}=U\Sigma V^T
\end{eqnarray}

ただし$\Sigma$は$N \times D$の行列で以下に示す性質を持っています。($U,V$は直交行列)

\begin{eqnarray}
\Sigma &=& [\Delta,\mathbf{0}] \\
\Delta &=& {\rm diag}(\sigma_1, \dots, \sigma_N)
\end{eqnarray}

このとき

\begin{eqnarray}
\mathbf{X}^T\mathbf{X}&=&(U\Sigma V^T)^TU\Sigma V^T\\
&=&V \Sigma^T U^T U \Sigma V^T\\
&=&V \Sigma^T \Sigma V^T\\
&=&V\begin{pmatrix} \Delta ^2 & \mathbf{0} \\
\mathbf{0} & \mathbf{0}
\end{pmatrix}V^T
\end{eqnarray}

となり$\mathbf{X}^T\mathbf{X}$は非正則となります。数学的には、「逆行列は求められない」で終わりなのですが、数値計算上では$\Sigma^T \Sigma$の要素が常識のように完全に$0$になるのではなく、$0$に近い非常に小さな値が入ります。それに対して無理やり逆行列を計算することになるので$\mathbf{w}$の要素にはとても大きな値が入ることになり、そのモデルから導かれた出力間には大きな差がでてしまいます。

解決法

$\mathbf{X}^T\mathbf{X}$が非正則である時にリッジ回帰という手法が使えます。(他にも過学習を防ぐ場合などにも使われます)
リッジ回帰では重みベクトル$\mathbf{w}$の大きさにペナルティを加えます。($\alpha>0$)

\begin{eqnarray}
    J_\mathrm{LS}(\mathbf{w})&=&\|\mathbf{y}-\mathbf{Xw}\|^2+\alpha\|\mathbf{w}\|^2
\end{eqnarray}

同様に$\mathbf{w}$で微分して0とおくと

\begin{eqnarray}
    \mathbf{w}=(\mathbf{X}^T\mathbf{X}+\alpha I)^{-1}\mathbf{X}^T\mathbf{y}
\end{eqnarray}

となりますね。
$\mathbf{X}$の特異値分解を考えると

\begin{eqnarray}
    \mathbf{X}&=&U\Sigma V^T\\
    \mathbf{X}^T\mathbf{X}+\alpha I &=& (U\Sigma V^T)^T U\Sigma V^T+\alpha I\\
    &=&V\Sigma^2V^T+\alpha VIV^T\\
    &=&V(\Sigma^2+\alpha I)V^T
\end{eqnarray}

となります。
$\mathbf{X}^T\mathbf{X}+\alpha I$の固有値は$\mathbf{X}$の特異値の2乗に正の値を足したものとなるため全て正となります。
よって固有値は0をとらず

\begin{eqnarray}
    \rm{det}(\mathbf{X}^T\mathbf{X}+\alpha I)>0
\end{eqnarray}

となり、$\mathbf{X}^T\mathbf{X}+\alpha I$は正則となって逆行列を持ちます！

まとめ

重回帰分析において入力の次元数に対してデータ数が足りない場合の解決法とそれで解決される理由を示しました。
皆さんも問題が発生した際は自分が扱っている手法が何を前提としているかを考えるのもいいかもしれません。

参考文献：パターン認識と機械学習(上・下)