3

我很难理解以下内容:

Scikit-learn 为 Ridge Regression 提供了多输出版本,只需传递一个 2D 数组 [n_samples, n_targets],但它是如何实现的呢?

http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html

假设每个目标的每个回归都是独立的是否正确?在这些情况下,我如何调整它以对每个回归使用单独的 alpha 正则化参数?如果我使用 GridSeachCV,我将不得不交出一个可能的正则化参数矩阵,或者这将如何工作?

在此先感谢-我一直在搜索几个小时,但找不到有关此主题的任何内容。

4

1 回答 1

5

我会试一试,因为我一直在为自己的工作进行一些研究。我会将问题分解为多个部分,因此您只能查看您感兴趣的部分:

Q1:多输出岭回归中每个目标(又名输出)的回归是否独立?

A1:我认为典型的 M 个输出的多输出线性回归与M 个独立的单输出线性回归是一样的。我认为是这种情况,因为多重输出案例的普通最小二乘表达式与 M 个独立的单输出案例(总和)的表达式相同。为了激发这一点,让我们考虑一个没有正则化的愚蠢的二元输出案例。

让我们考虑两个列向量输入 x 1x 2,以及相应的权重向量w 1w 2

这些给了我们两个单变量 输出,y 1 = x 1 w 1 T + e 1和 y 2 = x 2 w 2 T + e 2,其中 es 是独立误差。

误差平方和写成:

e 1 2 + e 2 2 = (y 1 - x 1 w 1 T ) 2 + (y 2 - x 2 w 2 T ) 2

我们可以看到这只是两个独立回归的平方误差之和。现在,为了优化,我们区分权重并将其设置为零。由于 y 1不依赖于w 2 ,对于 y 2w 1反之亦然,因此可以针对每个目标独立地进行优化。

我在这里考虑了一个示例以进行说明,但更多示例没有太大变化。你可以自己写出来。在表格中添加惩罚项 | w 1 | 或 | w 2 | 也不会改变这一点,因为w 2仍然不依赖于 y 1,反之亦然 y 2w 1

好的,这就是能让你获得 C- 的证明(有一位慷慨的教授)。就这反映 sklearn 而言,手动实现独立回归和内置的多输出支持将给出相同的结果。所以让我们用一些代码来检查一下(我正在使用带有 Jupyter 的 py2.7):

我们需要的东西

 import numpy as np
 import matplotlib.pyplot as plt
 from sklearn import linear_model, model_selection

设置数据

## set up some test data
# T samples, K features, M outputs (aka targets)
T = 1000
K = 100
M = 500

#get the samples from independent, multivariate normal
means_X = np.zeros(K)
cov_X = np.identity(K) 
X = np.random.multivariate_normal(means_X,cov_X,T)

#Make up some random weights. 
#Here I use an exponential form since that means some would be quite small, and thus regularization is likely to help
#However for the purposes of the example it doesn't really matter

#exponential weights
W = 2.0 ** np.random.randint(-10,0,M * K) 

#shape into a weight matrix correctly
W = W.reshape((K,M))

# get the ouput - apply linear transformation
Y = np.matmul(X, W)

# add a bit of noise to the output
noise_level = 0.1
noise_means = np.zeros(M)
noise_cov = np.identity(M) 

Y_nse = Y + noise_level * np.random.multivariate_normal(noise_means,noise_cov,T)

# Start with one alpha value for all targets 
alpha = 1

使用 sklearn 内置的多输出支持

#%%timeit -n 1 -r 1
# you can uncomment the above to get timming but note that this runs on a seperate session so 
# the results won't be available here 
## use built in MIMO support 

built_in_MIMO = linear_model.Ridge(alpha = alpha)
built_in_MIMO.fit(X, Y_nse)

为输出独立运行优化

# %%timeit -n 1 -r 1 -o
## manual mimo
manual_MIMO_coefs = np.zeros((K,M))

for output_index in range(M):
    
    manual_MIMO = linear_model.Ridge(alpha = alpha)
    manual_MIMO.fit(X, Y_nse[:,output_index]) 
    manual_MIMO_coefs[:,output_index] = manual_MIMO.coef_

比较估计值(不包括图)

manual_MIMO_coefs_T = manual_MIMO_coefs.T

## check the weights. Plot a couple
check_these_weights = [0, 10]
plt.plot(built_in_MIMO.coef_[check_these_weights[0],:],'r')
plt.plot(manual_MIMO_coefs_T[check_these_weights[0],:], 'k--')

# plt.plot(built_in_MIMO.coef_[check_these_weights[1],:],'b')
# plt.plot(manual_MIMO_coefs_T[check_these_weights[1],:], 'y--')

plt.gca().set(xlabel = 'weight index', ylabel = 'weight value' )
plt.show()

print('Average diff between manual and built in weights is %f ' % ((built_in_MIMO.coef_.flatten()-manual_MIMO_coefs_T.flatten()) ** 2).mean())


## FYI, our estimate are pretty close to the orignal too, 
plt.clf()
plt.plot(built_in_MIMO.coef_[check_these_weights[1],:],'b')
plt.plot(W[:,check_these_weights[1]], 'y--')
plt.gca().set(xlabel = 'weight index', ylabel = 'weight value' )
plt.legend(['Estimated', 'True'])

plt.show()

print('Average diff between manual and built in weights is %f ' % ((built_in_MIMO.coef_.T.flatten()-W.flatten()) ** 2).mean())

输出是(不包括这里的图)

Average diff between manual and built in weights is 0.000000 

Average diff between manual and built in weights is 0.000011 

所以我们看到内置的 sklearn 估计与我们的手动估计相同。然而,内置的更快,因为它使用矩阵代数来解决整个问题,这与我在这里使用的循环相反(关于 Ridge 正则化的矩阵公式,请参阅关于 Tikhonov 正则化的 wiki)。您可以通过取消注释上面的 %%timeit 魔法来自己检查)

Q2:我们如何为每个回归使用单独的 alpha 正则化参数?

A2:sklearn Ridge 对每个输出(目标)接受不同的正则化。例如继续上面的代码,为每个输出使用不同的 alpha:

# now try different alphas for each target.
# Simply randomly assign them between min and max range 
min_alpha = 0
max_alpha = 50
alphas = 2.0 ** np.random.randint(min_alpha,max_alpha,M)
built_in_MIMO = linear_model.Ridge(alpha = alphas)
built_in_MIMO.fit(X, Y_nse) 

如果将此与手动实现 M 个独立回归进行比较,每个回归都有自己的 alpha:

manual_MIMO_coefs = np.zeros((K,M))

for output_index in range(M):
    
    manual_MIMO = linear_model.Ridge(alpha = alphas[output_index])
    manual_MIMO.fit(X, Y_nse[:,output_index]) 
    manual_MIMO_coefs[:,output_index] = manual_MIMO.coef_

你得到相同的结果:

manual_MIMO_coefs_T = manual_MIMO_coefs.T

## check the weights. 
print('Average diff between manual and built in weights is %f ' % ((built_in_MIMO.coef_.flatten()-manual_MIMO_coefs_T.flatten()) ** 2).mean())

Average diff between manual and built in weights is 0.000000 

所以这些都是一样的。

但是,在这种情况下,性能很大程度上取决于求解器(就像@Vivek Kumar 直觉的那样)。

默认情况下,Ridge.fit() 使用 Cholesky 分解(至少对于非稀疏数据),查看 github 上的代码(https://github.com/scikit-learn/scikit-learn/blob中的_solve_cholesky /master/sklearn/linear_model/ridge.py),我看到当为每个目标单独选择 alpha 时,sklearn 实际上确实分别适合它们。我不知道这是 Cholesky 固有的还是只是实现的东西(我感觉是后者)。

但是,对于不同的求解器,例如基于 SVD 的 (_solve_svd()),代码似乎已经将不同的 alpha 合并到问题的矩阵公式中,所以整个事情只运行一次。这意味着当为每个输出单独选择 alpha 并且有许多输出时,svd 求解器可以快得多。

Q3:如何使用 GridSeachCV?我是否交出一个可能的正则化参数矩阵?

A3:我没有使用内置的网格搜索,因为它不太适合我的问题。但是,通过上述说明,实现这一点很简单;只需使用 sklearn.model_selection.KFold() 或类似方法获得一些 CV 折叠,然后使用不同的 alpha 对每个折叠进行训练:

from sklearn import metrics, model_selection
# just two folds for now
n_splits = 2
#logarithmic grid
alphas = 2.0 ** np.arange(0,10) 
kf = model_selection.KFold(n_splits=n_splits)

# generates some folds
kf.get_n_splits(X)

# we will keep track of the performance of each alpha here
scores = np.zeros((n_splits,alphas.shape[0],M))

#loop over alphas and folds
for j,(train_index, test_index) in enumerate(kf.split(X)):
    
    for i,alpha in enumerate(alphas):
        
        cv_MIMO = linear_model.Ridge(alpha = alpha)
        cv_MIMO.fit(X[train_index,:], Y_nse[train_index,:]) 
        cv_preds = cv_MIMO.predict(X[test_index,:])
        scores[j,i,:] = metrics.r2_score(Y_nse[test_index,:], cv_preds, multioutput='raw_values')

## mean CV score  
mean_CV_score = scores.mean(axis = 0)
# best alpha for each target
best_alpha_for_target = alphas[np.argmax(mean_CV_score,axis = 0)]

我写的很匆忙,所以请仔细检查。请注意,我们需要使用 metric 模块,因为内置的 Ridge.score() 对所有目标进行平均,我们不希望在这里。通过使用 multioutput='raw_values' 选项,我们保留每个目标的原始值。

希望有帮助!

于 2018-08-28T05:35:36.580 回答