はじめに

　スモールデータと機械学習、スパース回帰分析とパターン認識より、線形回帰（最小二乗法、LASSO）の理解を深めるために記事を作成します。n番煎じです。この記事では基本的にベクトルは縦ベクトルとして扱います。
　ついでに線形回帰の考え方を利用した実験計画法のサンプル選択方法の一つであるD最適基準の内容についても理解を出力します。

線形回帰

　データセット ($X, X \in \mathbb{R}^{N \times M}$, $N$: データ数, $M$: 特徴量の次元) 、および目的変数 ($y = {y_1 \ldots y_n}$) があります。

導出

$M$個の特徴量と$N$個のサンプルを持つデータセット$X (x_1 ... x_n) \in \mathbb{R}^{N \times M}$ および目的変数 ($y = {y_1 \ldots y_n}$) があります。

X
=\begin{pmatrix}
x_{1,1} \cdots x_{1,m} \\
\vdots \qquad \qquad \vdots \\
x_{n,1} \cdots x_{n,m}
\end{pmatrix}

　$X$ および$y$ に対してオートスケーリングの処理を入れることで、各特徴量の平均を0、分散を1に変換します。ここでは変換前の特徴量ごとの $X$ および $y$ の平均を $\mu_x, \mu_y$ 、分散を $\sigma^{2}_x, \sigma^{2}_y$ とします。

\begin{align}
(x_n - \mu_x)/ \sqrt{\sigma^{2}_x} \\
(y_n - \mu_y)/ \sqrt{\sigma^{2}_y}
\end{align}

この処理により、平均が0になるため、以下の線形回帰を考えるときに切片を考慮する必要がなくなります。(誤差: $\epsilon$, 係数: $\beta = \beta_1 \ldots \beta_m$)

y = X\beta  + \epsilon

また、誤差を考慮しない回帰直線からなる予測値 ($\hat{y}$) を以下のように記述します。

\hat{y} = X \beta

最小二乗法による回帰は、各次元（$X_1 \cdots X_M$）に対する係数($\beta_1 \cdots \beta_M$) を求める操作を行います。この時、$\beta$ は予測値 ($\hat{y} = X \beta$) と実測値 ($y$) の差 ($Q$) が小さくなるような値を選択します。また、$y$ と $\hat{y}$ の差は負の値をとる場合もあるため、2乗して計算することで、誤差の比較を行いやすくします。

\begin{align}
Q 
&= \sum_{i = 1}^N (y_i - \hat{y_i})^2 \\
&= (y - \hat{y})^\mathsf{T} (y - \hat{y}) \\
&= (y - X \beta)^\mathsf{T} (y - X \beta） \\
&= y^\mathsf{T} y - y^\mathsf{T} X\beta - \underbrace{(X\beta)^\mathsf{T} y}_{計算結果がスカラーになるため、 y^\mathsf{T} X\beta に変換} + (X\beta)^\mathsf{T} (X\beta) \\
&= y^\mathsf{T} y - 2y^\mathsf{T} X\beta  + (X\beta)^\mathsf{T} (X\beta)
\end{align}

$Q$ を最小にする$\beta$ を探すため、$\beta$ で偏微分します。

\begin{align}
\frac{\partial Q}{\partial \beta} 
&= 2X^\mathsf{T} X \beta - 2X^\mathsf{T} y　\\
&= 0 \\
X^\mathsf{T} X\beta 
&= X^\mathsf{T}y \\
\end{align}

$X^\mathsf{T}X$ が正則である場合、逆行列が存在するため、$(X^\mathsf{T}X)^{-1}$ を両辺に左からかけます。

\begin{align}
(X^\mathsf{T}X)^{-1}X^\mathsf{T} X\beta 
&= (X^\mathsf{T}X)^{-1}X^\mathsf{T}y \\
\beta 
&= (X^\mathsf{T}X)^{-1}X^\mathsf{T}y 
\end{align}

この操作により、$\beta$ をすでに得られている値から求めることができました。一方、$X$ がランク落ちしているなど、$X^\mathsf{T}X$ が正則にならないような条件の場合には、$\beta$ を求めることができません。また、$X$ がランク落ちまではいかなくても強い多重共線性がある場合、 $\beta$ が $X$ のわずかな数値の揺れに大きく影響され、不安定になります。特に実験などで得られる計測値はノイズを含むため、線形従属の関係にある特徴量を含む実験データで最小二乗法を組むと $\beta$ がノイズの影響を強く受けます。$X$ の特徴量間の多重共線性の程度を評価する指標として、VIF や $X$ の最大特異値 ($\sigma_{max}$)と最小特異値 ($\sigma_{min}$) の比である条件数 ($\kappa = \frac{\sigma_{max}}{\sigma_{min}}$) があります。このVIF および条件数が大きい（$e.g. VIF >10$）場合、行列のベクトル間の線形従属の関係が強くなります。

VIF

各特徴量($1\cdots m$)のVIF (Variance Inflation Factor, $VIF_1 \cdots VIF_m$) は以下の式で計算されます。 $$ VIF_m = \frac{1}{1-R_m^2} $$ $R_m^2$ は$x_m$ を$x_m$ 以外の特徴量で重回帰した際の決定係数($R^2$) です。 $$ \hat{x_m} = w_1 x_1 + \cdots + w_{m-1} x_{m-1} $$

$X$ に多重共線性の問題が生じている場合、VIF や PCA、後述する Lasso などを用いた変数削減により行列の悪条件を改善するテクニックがあります。

実装

sklearn.datasets.load_diabetesを使用して、numpyによる実装と sklearn による実装を比較します。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_diabetes

# load dataset
data = load_diabetes()

x = data.data
y = data.target
print(x.shape, y.shape)

def least_squares_regression(X, y):
    beta = np.linalg.inv(X.T @ X) @ X.T @ y
    return beta
    
# autoscaling
x_as = StandardScaler().fit_transform(x)
y_as = StandardScaler().fit_transform(y.reshape(-1, 1))

# calculate beta
beta_ls = least_squares_regression(x_as, y_as)

beta_ls_sklearn = LinearRegression().fit(x_as, y_as).coef_.T

plt.plot(beta_ls, 'o', label="numpy")
plt.plot(beta_ls_sklearn, 'x', label="sklearn")
plt.xlabel("features No.", fontsize=14)
plt.ylabel("beta", fontsize=14)
plt.legend()
plt.show()

$\beta$ の値が numpy によるスクラッチと sklearn で重なったため、実装は問題なさそうです。

# prediction y
y_pred = x_as @ beta_ls
plt.scatter(y_as, y_pred)
plt.xlabel("y_autoscaled", fontsize=14)
plt.ylabel("y_pred_autoscaled", fontsize=14)
plt.show()

Lasso (Least absolute shrinkage and selection operator) による変数選択

導出

　最小二乗法では、

 \underset{\beta \in \mathbb{R}^d}{\min} ||y - X \beta||^2_{2}

が最小となる $\beta$ を探索しました。Lasso では、L1 ノルムを加えて、$\beta_{lasso}$ がとりえる範囲を制限します。

 \underset{\beta_{lasso} \in \mathbb{R}^d}{\min} ||y - X \beta_{lasso}||^2_{2} + \lambda ||\beta_{lasso}||_1

わかりやすく表現するために、第1項に対し、最小二乗法の $\beta$ を導入します。

\begin{align}
y - X\beta_{lasso} 
&= y - X\beta + X\beta - X\beta_{lasso} \\
&= \underbrace{(y - X\beta)}_{最小二乗法の残差 e}　+ \underbrace{X(\beta - \beta_{lasso})}_{最小二乗法の解とlasso の解の差} \\
&= e + X(\beta - \beta_{lasso})
\end{align}

$||y - X\beta_{lasso}||^2_2$ を展開します。

\begin{align}
||y - X\beta_{lasso}||^2_2 
&= ||e + X(\beta - \beta_{lasso})||^2_2 \\
&= ||e||^2_2 + ||X(\beta - \beta_{lasso})||^2_2 + \underbrace{2eX(\beta - \beta_{lasso})}_{残差ベクトル e は X の空間に対して垂直になるため、eX = 0} \\
&= ||e||^2_2 +||X(\beta - \beta_{lasso})||^2_2
\end{align}

最小二乗法の残差 $e$ は求まっているため、最小化に関して、

||\beta - \beta_{lasso}||^2_2　+ ||\lambda \beta||_1

を考えます。計算の都合上、調整パラメーターを $2\lambda$ とします。また、第 1 項はL2 ノルムであるため、$\beta$ と $\beta_{lasso}$ を入れかえ、$\beta_{lasso}$ の各次元を最小化します。

 \underset{\beta_{lasso} \in \mathbb{R}^d}{\min} ||\beta_{lasso} - \beta||^2_{2} + \lambda ||\beta_{lasso}||_1　  = \sum_{j =1}^d \underset{\beta_{lasso_j} \in \mathbb{R}}{\min} ||\beta_{lasso_j} - \beta_j||^2_{2} + \lambda ||\beta_{lasso_j}||_1

以下の式を $\beta_{lasso_j}$ が正の場合 ($\beta_{lasso_j} \geq 0$) と負の場合 ($\beta_{lasso_j} \leq 0$)に分けて展開します。

正の場合 ($\beta_{lasso_j} \geq 0$)

\begin{align}
f(\beta_{lasso_j}) &= ||\beta_{lasso_j} - \beta_j||^2_{2} + \lambda ||\beta_{lasso_j}||_1 \\
&= (\beta_{lasso_j} - \beta_j)^2 + \lambda \beta_{lasso_j}　\\
f^{\prime}(\beta_{lasso_j}) 
&= 2\beta_{lasso_j} - 2\beta_j + \lambda \\
&= 0 \\
\beta_{lasso_j}  &= \beta_j - \frac{\lambda}{2} \geq 0 \\
\beta_j &\geq \frac{\lambda}{2}
\end{align}

そのため、$\beta_j \geq \frac{\lambda}{2}$ の時、$\beta_{lasso_j}$ は $\beta_j - \frac{\lambda}{2} $ で最小値をとります。

負の場合 ($\beta_{lasso_j} \leq 0$)

\begin{align}
f(\beta_{lasso_j}) &= ||\beta_{lasso_j} - \beta_j||^2_{2} + \lambda ||\beta_{lasso_j}||_1 \\
&= (\beta_{lasso_j} - \beta_j)^2 \underbrace{- \lambda \beta_{lasso_j}}_{\beta_{lasso_j} が負であるため}　\\
f^{\prime}(\beta_{lasso_j}) 
&= 2\beta_{lasso_j} - 2\beta_j - \lambda \\
&= 0 \\
\beta_{lasso_j}  &= \beta_j + \frac{\lambda}{2} \leq 0 \\
\beta_j &\leq - \frac{\lambda}{2}
\end{align}

そのため、$\beta_j \leq -\frac{\lambda}{2}$ の時、$\beta_{lasso_j}$ は $\beta_j + \frac{\lambda}{2} $ で最小値をとります。
まとめると、$\beta_{lasso_{j}}$ は以下の条件で最小になります。

\beta_{lasso_j} = S_{\lambda}(\beta_j) = \begin{cases}
\beta_j - \frac{\lambda}{2} & (\beta_j \geq \frac{\lambda}{2}) \\
    0 &  (\frac{\lambda}{2} \geq \beta_j \geq -\frac{\lambda}{2}) \\
\beta_j + \frac{\lambda}{2} & (\beta_j \leq \frac{\lambda}{2}) \\
\end{cases}

Lasso と Ridge の比較

Lasso は $\lambda$ をしきい値として、最小二乗法の $\beta$ を調節、または $0$ に変換することで、過学習を防ぐ機能があります。$\lambda$ を変化させた際の $\beta_{lasso}$ の変化をプロットしたものを解パスといいます。Lasso の解パスをRidge と比較すると、Lasso は基準以下の説明変数を $0$ に変換する一方、Ridge は説明変数の値が $0$ になりにくい特徴があります。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Ridge, Lasso
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_diabetes

# load dataset
data = load_diabetes()

x = data.data
y = data.target

x_as = StandardScaler().fit_transform(x)
y_as = StandardScaler().fit_transform(y.reshape(-1, 1))

n = 10
lambda_scale_ridge = np.logspace(0, 4 , n)
lambda_scale_lasso = np.linspace(0, 1 , n)
labels = data.feature_names
df_ridge_coef = np.zeros((len(labels), n))
df_lasso_coef = np.zeros((len(labels), n))
fig, ax = plt.subplots(1,2, figsize=(12, 8), tight_layout=True)
ax = ax.flatten()
for i in range(n):
    r = lambda_scale_ridge[i]
    model_ridge = Ridge(alpha=r)
    ridge = model_ridge.fit(x_as, y_as)
    ridge_coef_ = ridge.coef_.T
    df_ridge_coef[:, i] = ridge_coef_.flatten()

    l = lambda_scale_lasso[i]
    model_lasso = Lasso(alpha=l)
    lasso = model_lasso.fit(x_as, y_as)
    lasso_coef_ = lasso.coef_.T
    df_lasso_coef[:, i] = lasso_coef_.flatten()

ax[0].plot(lambda_scale_ridge, df_ridge_coef.T, label=labels)
ax[0].set_xlabel("lambda", fontsize=14)
ax[0].set_ylabel("beta_ridge", fontsize=14)

ax[1].plot(lambda_scale_lasso, df_lasso_coef.T, label=labels)
ax[1].set_xlabel("lambda", fontsize=14)
ax[1].set_ylabel("beta_lasso", fontsize=14)


plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

VIF による多重共線性の確認

　元のデータセットの VIF を確認し、多重共線性の程度を調べました。VIF を調べる目的は多重共線性の程度が弱くなるような変数選択であるため、先人の記事を参考に、VIFを算出し、最も高い VIF が目的の閾値以下になるまで再帰的に計算するコードを書きました。ついでに変数選択による汎化性能の変化をデータセットの内挿、外挿領域で見てみたかったため、データセットの $\pm 4 SD$ で一様分布を用いてデータを生成しました。生成したデータの目的変数を全ての説明変数または変数選択した説明変数による線形回帰で予測し、主成分分析で圧縮したプロットに予測値を色付けしました。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_diabetes
from sklearn.decomposition import PCA


def generate_diabetes_data():
   diabetes = load_diabetes()
   X_original = diabetes.data
   feature_names = diabetes.feature_names
   
   df_original = pd.DataFrame(X_original, columns=feature_names)
   
   # data generation
   n_samples = 10000
   generated_data = {}
   
   for col in df_original.columns:
       mean = df_original[col].mean()
       std = df_original[col].std()
       
       # 4 SD
       low = mean - 4 * std
       high = mean + 4 * std
       generated_data[col] = np.random.uniform(low, high, n_samples)
   
   df_generated = pd.DataFrame(generated_data)
   df_original.insert(loc = 0, column="target", value=diabetes.target)
   
   return df_original, df_generated


def visualize_pca_colormap(df_train ,target,  df_generated):
   y_train = target
   X_test = df_generated
   
   df_train_trimmed = df_train.drop(["s1", "s4"], axis=1)
   df_generated_trimmed = df_generated.drop(["s1", "s4"], axis=1)

   # Linear Regression
   lr = LinearRegression()
   lr.fit(df_train, y_train)
   pred_lr = lr.predict(X_test)
   
   # linear regression with df_train_trimmed, df_generated_trimmed
   lr_trimmed = LinearRegression()
   lr_trimmed.fit(df_train_trimmed, y_train)
   pred_lr_trimmed = lr_trimmed.predict(df_generated_trimmed)


   # pca with original df_train
   pca = PCA(n_components=2)
   pca.fit(df_train)
   coords_train = pca.transform(df_train)
   coords_test = pca.transform(X_test)
   exp_ratio = pca.explained_variance_ratio_

   # pca with trimmed, df_train_trimmed
   pca_trimmed = PCA(n_components=2)
   pca_trimmed.fit(df_train_trimmed)
   coords_train_trimmed = pca_trimmed.transform(df_train_trimmed)
   coords_test_trimmed = pca_trimmed.transform(df_generated_trimmed)
   exp_ratio_trimmed = pca_trimmed.explained_variance_ratio_

   # min, max value
   all_values = np.concatenate([y_train, pred_lr])
   vmin = all_values.min()
   vmax = all_values.max()
   
   # min, max value for trimmed    
   all_values_trimmed = np.concatenate([y_train, pred_lr_trimmed])
   vmin_trimmed = all_values_trimmed.min()
   vmax_trimmed = all_values_trimmed.max()

   # visualize
   fig, axes = plt.subplots(1, 2, figsize=(18, 8))
   axes = axes.flatten()
   
   if vmax_trimmed > vmax:
       vmax = vmax_trimmed
   if vmin_trimmed < vmin:
       vmin = vmin_trimmed
   
   def plot_data(ax, title, pred_values, exp_ratio):
       # plot generated data
       sc = ax.scatter(coords_test[:, 0], coords_test[:, 1], 
                       c=pred_values, cmap='jet', vmin=vmin, vmax=vmax,
                       s=15, alpha=0.4, edgecolors='none', marker='o',
                       label='Generated (Pred)')
       
       # plot original data
       ax.scatter(coords_train[:, 0], coords_train[:, 1], 
                  c=y_train, cmap='jet', vmin=vmin, vmax=vmax,
                  s=60, alpha=1.0, edgecolors='black', linewidth=1.2, marker='o',
                  label='Original (True)')
   
       ax.set_title(title, fontsize=14)
       ax.set_xlabel(f'PC1 ({exp_ratio[0]*100:.1f}%)')
       ax.set_ylabel(f'PC2 ({exp_ratio[1]*100:.1f}%)')
       ax.grid(True, linestyle='--', alpha=0.5)
       ax.legend(loc='upper right')
       
       cbar = plt.colorbar(sc, ax=ax)
       cbar.set_label('Target Value (True & Pred)', rotation=270, labelpad=15)

   # plot linear regression
   plot_data(axes[0], 'Linear Regression Prediction', pred_lr, exp_ratio=exp_ratio)
   
   # plot linear regression with trimmed data
   plot_data(axes[1], 'Linear Regression Prediction (drop: "s1", "s4")', pred_lr_trimmed, exp_ratio=exp_ratio_trimmed)

   plt.tight_layout()
   plt.show()
   print(f" min_linear_regression: {pred_lr.min()} \n min_lr_trimmed: {pred_lr_trimmed.min()}\n min_original: {y_train.min()}")

df_orig, df_gen_labeled = generate_diabetes_data()
visualize_pca_colormap(df_orig.iloc[:, 1:],df_orig.iloc[:, 0], df_gen_labeled)

def track_vif_history(df, threshold=10.0):
   # create vif_df 
   vif_history_df = pd.DataFrame(index=df.columns)
   
   df_current = df.copy()
   trial = 1
   
   while True:
       # add constant 
       df_with_const = add_constant(df_current)
       # calc. vif
       current_vifs = pd.Series(
           [variance_inflation_factor(df_with_const.values, i) 
            for i in range(df_with_const.shape[1])],
           index=df_with_const.columns
       ).drop('const')
       
       # add vif value to history dataframe
       col_name = f'vif_trial_{trial}'
       vif_history_df[col_name] = current_vifs
       
       # check max vif
       max_vif = current_vifs.max()
       max_feature = current_vifs.idxmax()
       
       if max_vif <= threshold:
           print(f"Trial {trial}: 完了 (Max VIF: {max_vif:.2f})")
           break
       
       print(f"Trial {trial}: {max_feature} を削除 (VIF: {max_vif:.2f})")
       
       # drop max vif feature
       df_current = df_current.drop(columns=[max_feature])
       trial += 1
       
   return vif_history_df, df_current



vif_history, df_current = track_vif_history(df_x, threshold=5)

s1、s4 が順次削除され、それに応じてs2やs3 のVIFが低下することがわかりました。

s1, s4 を削除することで、線形回帰の予測値が極端な値をとりにくくなりました。

実験計画法

　実験計画法（Design of Experiments, DOE）は、効率的に実験を行い、得られたデータから最大限の情報を引き出すための統計的な手法です。直行表などに代表される実験計画法と、順次実験条件を生成する適応的実験計画法があり、物質生産の最適化や寄与する変数の探索などに用いられています。実験計画法に関する解説記事は多々あるため、本記事では適応的実験計画法に用いられるD最適基準の意味についての自身の理解を出力します。

適応的実験計画法の概要

　適応的実験計画法では、以下の順でデータの取得が進められます

特徴量探索範囲の設定
設定探索範囲内の実験点生成
繰り返し
1. 実験による応答データ収集
2. 応答データを基にモデル作成
3. 目的の応答が得られる確率の高い実験点生成

3.2.でよいモデルを作るには、特徴量間に多重共線性がないようにする必要があります。D最適基準を用いることで、2.の操作を行う際に多重共線性の程度が小さくなるようなデータ点を選択することができます。

D最適基準の意義

　最小二乗法による線形回帰によって計算される $\beta$ の期待値（$\mathbb{E}[\beta]$）および共分散行列（$Cov[\beta]$）は以下のように計算できます。

\begin{align}
\mathbb{E}[\beta] 
&= \mathbb{E}[(X^\mathsf{T}X)^{-1}X^\mathsf{T} y] \\
&= (X^\mathsf{T}X)^{-1}X^\mathsf{T} \mathbb{E} [X \beta + \epsilon] \\
&= (X^\mathsf{T}X)^{-1}X^\mathsf{T} X \beta \\
&= \beta \\
\end{align}

$\beta = (X^\mathsf{T}X)^{-1}X^\mathsf{T} y = Cy$ と置くと

\begin{align}
Cov[\beta] 
&=  \mathbb{E}[(\beta - \mathbb{E}[\beta])(\beta - \mathbb{E}[\beta])^\mathsf{T}] \\
&= \mathbb{E}[(Cy - \beta)(Cy - \beta)^\mathsf{T}] \\
&= \mathbb{E}[(\underbrace{CX\beta}_{(X^\mathsf{T}X)^{-1}X^\mathsf{T}X = I} + C\epsilon - \beta)(CX\beta + C\epsilon - \beta)^\mathsf{T}] \\
&= \mathbb{E}[(C\epsilon)(C\epsilon)^\mathsf{T}] \\
&= \mathbb{E}[\epsilon \epsilon^\mathsf{T}]C C^\mathsf{T} \\
&= \sigma^2 C C^\mathsf{T} \\
C C^\mathsf{T} 
&= ((X^\mathsf{T}X)^{-1}X^\mathsf{T})((X^\mathsf{T}X)^{-1}X^\mathsf{T})^\mathsf{T} \\
&= ((X^\mathsf{T}X)^{-1}X^\mathsf{T})(X (X^\mathsf{T}X)^{-1}) \\
&= (X^\mathsf{T}X)^{-1}
\end{align}

そのため、$Cov[\beta] = \sigma^2 (X^\mathsf{T}X)^{-1}$ の体積（行列式: $det(Cov[\beta])$）を最小化することは、$\beta$ がとりえる空間を小さくすることと同じため、$\beta$ を精度よく推測することができます。$\sigma^2$ の推定には実験結果 ($y$) の値が必要なため、計算から除外すると、

\begin{align}
{Minimize} \ {det}((X^\mathsf{T}X)^{-1}) 
\end{align}

になるような $X$ を生成することで、データ取得後の解析が行いやすくなります。

おわりに

最小二乗法、LASSO、D最適基準を理解するために、自身の学びを出力しました。違っている点など見つかれば適宜修正します。ご指摘いただけると感謝です