参考リンク
- Mukherjee A, Zhu J. Reduced Rank Ridge Regression and Its Kernel Extensions. Stat Anal Data Min. 2011;4(6):612–622. doi:10.1002/sam.10138
- RRRR - Python implementation of Reduced Rank Ridge Regression
簡単に！
「縮小ランク回帰」と「リッジ回帰」を組み合わせて、説明変数の共線性と目的変数の低次元構造に同時に適合するようにした、線形モデル回帰の手法。
Reduced Rank Regression（縮小ランク回帰）とは？
一般化線形モデル
問題設定として、説明変数 $X_{N \times P}$により、目的変数 $Y_{N \times Q}$を予測する（回帰する）ことを考えます。
ここでよく用いられるのが、次の式で表される一般化線形モデルです。$$Y=XB+E$$ここで、$B_{P \times Q}$は係数行列であり、$E_{N \times Q}$はランダムな誤差行列を表しています。
この最小二乗法（OLS）による解は以下のようになります。$$\hat{B}_{OLS}=(X^TX)^{-1}X^TY$$
今回は、一般化線形モデルを最小二乗法で解くときの２つの問題点に着目することにします。
- 目的変数 $Y_{N \times Q}$の真の次元が $Q$ より小さいときに、性能が落ちる。（$Y$のランク落ち）
- 説明変数 $X_{N \times P}$の行同士の相関が高いときに、性能が落ちる。（$X$の共線性）
次元削減
このように、高次元のデータの中に潜在的な低次元の構造がある場合には、次元削減が有効です。有名なのでいえば、主成分分析（PCA）が挙げられるでしょう。主成分分析によって次元削減した後に回帰をするモデルを、主成分回帰（PCR）と呼んだりします。
主成分回帰のように、説明変数の中から低次元の因子（Factor）を選んで目的変数を回帰するようなモデルは、線形因子モデル（Linear Factor Model）と呼ばれています。例えば、独立成分回帰、偏最小二乗回帰、正準相関分析などがあります。
これとは別の次元削減の手法として、縮小ランク回帰（Reduced Rank Regression）があります。縮小ランク回帰では、係数行列 $B$ のランクを制限しながら誤差関数を最小化することで、潜在的な低次元の構造を仮定しつつ回帰することができます。
参考：縮小ランク回帰
正則化
２つめの問題点である説明変数 $X$ の共線性に対しては、正則化（regularization）がしばしば行われます。有名な手法に、リッジ回帰（Ridge Regression）やLASSO回帰（LASSO Regression）があります。LASSO回帰は、主にスパース性の推定に用いられ、特徴量選択の手法としても用いられています。リッジ回帰は、説明変数の共線性による不良設定問題に対応するために広く用いられています。
Reduced Rank Ridge Regression （縮小ランクリッジ回帰）への拡張
2つの制約を誤差関数に導入する
以上を踏まえて、２つの問題に対処するために、次の２つの制約を二乗誤差に追加することにします。
- リッジ回帰の正則化項
- $\hat{B}$ のランク制約
したがって、次のように誤差の最小化を行うことで、$\hat{B}(\lambda,r)$ を得ることができます。
$$\underset{ \lbrace B:rank(B) \leq r \rbrace}{argmin} \Vert Y-XB \Vert_F^2 + \lambda \Vert B \Vert_F^2$$
ただし、$r \leq \min \lbrace P, Q \rbrace$ とし、$\Vert \cdot \Vert_F^2$ はフロベニウスノルム（行列ノルム）とします。
これを縮小ランク回帰の枠組みで考えるために、
X_{(N+P)\times P}^* =
\left( \begin{array}{c} X \\ \sqrt{\lambda}I \end{array}\right), \
Y_{(N+P)\times Q}^* =
\left( \begin{array}{c} Y \\ 0 \end{array}\right)
とすることで、誤差関数を次のように表現します。
$$\underset{ \lbrace B:rank(B) \leq r \rbrace}{argmin} \Vert Y^*-X^*B \Vert_F^2$$
また、リッジ回帰の推定値 $ \hat{Y_R}^* = X^* \hat{B_R}^*$ を用いて、正規直交性から、
\Vert Y^*-X^*B \Vert_F^2 = \Vert Y^*-\hat{Y}_R^* \Vert_F^2 + \Vert \hat{Y}_R^*-X^*B \Vert_F^2
と変形できるので、第一項は $B$ に依存しないため、誤差関数は次のように表現できます。
\underset{ \lbrace B:rank(B) \leq r \rbrace}{argmin} \Vert \hat{Y}_R^*-X^*B \Vert_F^2
線形問題による解を導出する
ここで、次のように特異値分解が与えられると仮定します。
\hat{Y}_R^* = \sum_{i=1}^{\tau} \sigma_i u_i v_i^T
ここで、$\sigma_i$ は特異値、$u_i, v_i$ は左右の特異値ベクトルであり、$\tau$ は $\hat{Y}_R^*$ のランク（正則化により $Q$ となっている）となります。
これは、フロベニウスノルムにおけるエッカート・ヤングの定理（参考）として考えることができ、最適なランク $r$ による近似は次のようになります。
\hat{Y}_r^* = \sum_{i=1}^{r} \sigma_i u_i v_i^T
これを用いて、最適な係数行列 $\hat{B}(\lambda,r)$ を表現しましょう。
\hat{Y}_r^* = \sum_{i=1}^{r} \sigma_i u_i v_i^T = \left(\sum_{i=1}^{\tau} \sigma_i u_i v_i^T\right) \left(\sum_{j=1}^{r} u_j v_j^T\right) = \hat{Y}_r^* P_r = X^* \hat{B}_R^* P_r = X^* \hat{B}(\lambda,r)
とすることで、解として $\hat{B}(\lambda,r) = \hat{B}_R^* P_r $ が得られます。
さらに、リッジ回帰の最小二乗解を用いて以下のように表現することができます。
\hat{B}(\lambda,r) = \hat{B}_R^* P_r = (X^T X + \lambda I)^{-1}X^T Y P_r \\
\hat{Y}(\lambda,r) = X \hat{B}(\lambda,r) = X (X^T X + \lambda I)^{-1}X^T Y P_r = \hat{Y}_{\lambda} P_r
ここで、$\hat{Y}_{\lambda}$ は $\lambda$ におけるリッジ回帰の推定解となり、「リッジ回帰の推定解を$r$次元の空間に射影することで最適解を得ている」と解釈することができます。そのため、$r=Q$では、単純なリッジ回帰の解が最適解となることになります。
パラメータのチューニング
パラメータ $\lambda,r$ のチューニングが、縮小ランクリッジ回帰の性能を左右することになります。参照論文では、K-fold 交差検証 を行って、最適なパラメータを決定していました。
パッケージを試してみる
データ
用意されていたデモ用のデータを使用しましょう。
import numpy as np
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import PredefinedSplit
import reduced_rank_regressor as RRR
import matplotlib.pyplot as plt
%matplotlib inline
N_PARAMETERS_GRID_SEARCH = 20
Data_path = "data/"
# Load Data
trainX = np.loadtxt(Data_path+"trainX.txt")
testX = np.loadtxt(Data_path+"testX.txt")
validX = np.loadtxt(Data_path+"validX.txt")
trainY = np.loadtxt(Data_path+"trainY.txt")
testY = np.loadtxt(Data_path+"testY.txt")
validY = np.loadtxt(Data_path+"validY.txt")
# Inspection of Data
f,ax = plt.subplots(1,2,figsize=(10,10))
ax[0].imshow(np.corrcoef(trainX), cmap='jet')
ax[0].set_title('trainX')
ax[1].imshow(np.corrcoef(trainY), cmap='jet')
ax[1].set_title('trainY')
plt.show()
相関行列から、目的変数 $Y$ には低次元の構造があることがわかります。
交差検証
制約ランクと正則化の強さについてハイパーパラメータの交差検証を行い、最適化を行います。
# Cross-validation setup. Define search spaces
rank_grid = np.linspace(1,min(min(trainX.shape),min(trainY.shape)), num=N_PARAMETERS_GRID_SEARCH)
rank_grid = rank_grid.astype(int)
reg_grid = np.power(10,np.linspace(-20,20, num=N_PARAMETERS_GRID_SEARCH+1))
parameters_grid_search = {'reg':reg_grid, 'rank':rank_grid}
valid_test_fold = np.concatenate((np.zeros((trainX.shape[0],))-1,np.zeros((validX.shape[0],))))
ps_for_valid = PredefinedSplit(test_fold=valid_test_fold)
# Model initialisation
rrr = RRR.ReducedRankRegressor()#rank, reg)
grid_search = GridSearchCV(rrr, parameters_grid_search, cv=ps_for_valid,
scoring='neg_mean_squared_error')
print("fitting...")
grid_search.fit(np.concatenate((trainX,validX)), np.concatenate((trainY,validY)))
# Display the best combination of values found
print(grid_search.best_params_)
means = grid_search.cv_results_['mean_test_score']
means = np.array(means).reshape(N_PARAMETERS_GRID_SEARCH, N_PARAMETERS_GRID_SEARCH+1)
print(grid_search.best_score_)
[output]
fitting...
{'rank': 316, 'reg': 1e-20}
-75126.47521541138
交差検証の結果について視覚化しましょう。最適ランクの部分(316付近)でスコアの順位が高くなっていることがわかります。
# Show CV results
scores = [x for x in grid_search.cv_results_['rank_test_score']]
scores = np.array(scores).reshape(N_PARAMETERS_GRID_SEARCH, N_PARAMETERS_GRID_SEARCH+1)
f,ax = plt.subplots(1,1,figsize=(10,10))
cbar = ax.imshow(scores, cmap='jet')
ax.set_title('Test score Ranking')
ax.set_xlabel('Regression')
ax.set_ylabel('Rank')
ax.set_xticks(list(range(len(reg_grid))))
ax.set_yticks(list(range(len(rank_grid))))
ax.set_xticklabels([str("%0.*e"%(0,x)) for x in reg_grid],rotation=30)
ax.set_yticklabels([str(x) for x in rank_grid])
f.colorbar(cbar)
plt.show()
推論
交差検証で求めた最適なハイパーパラメータを元に推論をします。
# Train a model with the best set of hyper-parameters found
rrr.rank = int(grid_search.best_params_['rank'])
rrr.reg = grid_search.best_params_['reg']
rrr.fit(trainX, trainY)
# Testing
Yhat = rrr.predict(testX).real
MSE = (np.power((testY - Yhat),2)/np.prod(testY.shape)).mean()
print("MSE =",MSE)
f,ax = plt.subplots(1,2,figsize=(10,10))
ax[0].imshow(np.corrcoef(testY), cmap='jet')
ax[0].set_title('testY')
ax[1].imshow(np.corrcoef(Yhat), cmap='jet')
ax[1].set_title('Yhat')
plt.show()
[output]
MSE = 0.1508996355795968
MSE（平均二乗誤差）は小さくなり、$\hat{Y}$ はよく予測できていると思われます。
まとめ
縮小ランクリッジ回帰は低次元である目的変数に適合する能力が高く、現実のデータ構造に広く適用できる可能性がありそうです。また、パッケージは、scikit-learn に準じて実装されているので、使い勝手がよさそうです。