スパース主成分分析とは

主成分分析は、[1]によれば、『しかし、各主成分の寄与率が小さいとき、（例えば 10 個の特徴量が持つ情報量を 5 個の主成分で 6 割しか説明できないよう場合であれば、）データ全体を解釈するために、より多くの主成分を選ぶ必要がある。より多くの主成分を選ぶと、その主成分の解釈が難しくなる。この主成分に対してもスパース推定というアプローチをすることができる。主成分を推定する際に、主成分がなるべく 0 になるように推定するアプローチである。』

つまり、
次元圧縮の際、主成分分析は二乗誤差最小化を使うが、スパース主成分分析は、L1正則化(Lasso回帰)を行う。これにより、各主成分から解釈に不要な変数を取り除くことができる。[2]

主成分分析[2]


\hat{V} = \underset{V_k}{\operatorname{argmin}} \sum_{i=1}^n　||x_i - \hat{x_i} ||^2 \quad s.t. \quad \hat{x_i} = V_k V_k^T x_i.

主成分分析についてよくまとまった記事がありましたので、詳しく知りたい方はどうぞ。

スパース主成分分析[3]


V(\hat{A}, \hat{B}) = \underset{A,B}{\operatorname{argmin}} \{ \sum_{i=1}^n　||x_i - AB^Tx^i ||^2 + \lambda_2 \sum_{j=1}^k　|| \beta_j ||^2 + \sum_{j=1}^k \lambda_{1,j} || \beta_j ||_1 \}.

レポジトリーはこちらです。

スパース推定(Lasso回帰)について


E(x) = ||y - Ax||^2 + \lambda \sum_{i=1}^n|w|_1

ここで、過学習を説明する、過学習とは、訓練データを用いた場合は、最適解ではあるが、
未知のデータに適用した場合は、最適解にはならないというものである。

そこで、通常の最小二乗法に対して、正則化項を加えたものを(Lasso回帰、Ridge回帰)という。
一般にL1正則化(Lasso回帰)にスパース回帰に向いているが、L2正則化(Ridge回帰)は向いていない。

上の式の意味は
正則化項で答えの候補を絞り込み
観測データから候補を出す

この二つのバランスさせる。

Lasso回帰、Ridge回帰についての詳しい解説はこちら

正則化項の決定

スパース主成分分析は正則化項を自由に設定できる。
正則化項を変更することで、過学習を回避できる。
正則化項の決定には、ベイズ情報量規準(以下、bicと記す)、赤池情報量基準(以下、aicと記す)を使う。

[5]によれば、『予測ではなく妥当なモデルの構造を知りたい』ときに使う。

以下[6]より、


AIC = -2 \sum_{i=1}^n logp(X_i| W^*) + 2d


BIC =  -2 \sum_{i=1}^n logp(X_i| W^*) + dlogn

実装

import numpy as np
import numpy.random as random
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib as mpl
import seaborn as sns
import scipy as sp
from sklearn.linear_model import LassoLarsIC
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_friedman1
from sklearn.decomposition import SparsePCA

X, y = make_friedman1(n_samples=200, n_features=30, random_state=0)

#標準化
scaler = StandardScaler()
scaler.fit(X)
X = scaler.transform(X)

spca_std = SparsePCA(n_components=5, random_state=0)
spca_std.fit(X)
X_spca_std = spca_std.transform(X)

主な出力

#データから抽出されたスパース成分。
print(spca_std.components_)
#各イテレーションにおけるエラーのベクトル
print(spca_std.error_)
#想定される構成要素の数(既存バグでエラー、予測変換でも出てこない)
#print(spca_std.n_components_)
print(spca_std.n_components)
#イテレーションの実行回数
print(spca_std.n_iter_)
#トレーニングセットから推定された、各特徴の経験的平均値。X.mean(axis=0)と同じです。
print(spca_std.mean_)

正則化項の決定

# defauoltはalpha=1で出す。スパース性が低いとalpha=0が算出されるので注意が必要
model_bic = LassoLarsIC(criterion='bic')
model_bic.fit(X, y)
alpha_bic_ = model_bic.alpha_
print(alpha_bic_)
model_aic = LassoLarsIC(criterion='aic')
model_aic.fit(X, y)
alpha_aic_ = model_aic.alpha_
print(alpha_aic_)

代入

spca_bic = SparsePCA(n_components=5, random_state=0, alpha = alpha_bic_)
spca_bic.fit(X)
X_spca_bic = spca_bic.transform(X)
print(X_spca_bic.shape)

spca_aic = SparsePCA(n_components=5, random_state=0, alpha = alpha_aic_)
spca_aic.fit(X)
X_spca_aic = spca_aic.transform(X)
print(X_spca_aic.shape)

結果

print('defoult a=1                 :',np.mean(spca_std.components_ == 0))
#defoult a=1                 : 0.6933333333333334
print('bic a=0.02871277337904591   :',np.mean(spca_bic.components_ == 0))
#bic a=0.02871277337904591   : 0.12666666666666668
print('aic a=0.022802510590416925  :',np.mean(spca_aic.components_ == 0))
#aic a=0.022802510590416925  : 0.07333333333333333

情報量規準を用いて正則化項を変更したことで、主成分が0に近づいたことが分かる。

参考文献

[1]

[2]

[3]

[4]

[5]

[6]