Как многоплановая Ридж-регрессия работает в скиките?

Question

Я изо всех сил пытаюсь понять следующее:

Scikit-learn предлагает несколько выходных версий для регрессии гребня, просто передавая двумерный массив [n_samples, n_targets], но как это реализовать?

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

Правильно ли предположить, что каждая регрессия для каждой цели независима? В этих условиях, как я могу адаптировать это, чтобы использовать отдельные параметры альфа-регуляризации для каждой регрессии? Если я использую GridSeachCV, мне придется передать матрицу возможных параметров регуляризации, или как это будет работать?

Заранее спасибо - я искал часами, но ничего не смог найти по этой теме.

Answer 1

Я дам этому шанс, так как я немного разбираюсь в этом для своей собственной работы. Я разделю вопрос на части, чтобы вы могли смотреть только на те, которые вас интересуют:

Q1: Независима ли регрессия для каждой цели (или выхода) в множественной выходной регрессии Риджа?

A1: Я думаю, что типичная линейная регрессия с несколькими выходами с М выходами такая же, как М независимая линейная регрессия с одним выходом. Я думаю, что это так, поскольку выражение для обычных наименьших квадратов для случая множественного выхода такое же, как выражение для (суммы) M независимых, единичных выходных случаев. Чтобы мотивировать это, давайте рассмотрим глупый, двумерный случай вывода без регуляризации.

Рассмотрим два входных вектора столбцов x ₁ и x ₂ с соответствующими весовыми векторами w ₁ и w ₂.

Они дают нам два одномерных выхода: y ₁ = x ₁ w ₁ ^T + e ₁ и y ₂ = x ₂ w ₂ ^T + e ₂, где es - независимые ошибки.

Сумма квадратов ошибок выражается как:

e ₁ ² + e ₂ ² = (y ₁ - x ₁ w ₁ ^T) ² + (y ₂ - x ₂ w ₂ ^T) ²

Мы можем видеть, что это просто сумма квадратов ошибок двух независимых регрессий. Теперь, чтобы оптимизировать, мы дифференцируем по весам и устанавливаем в ноль. Поскольку y ₁ не зависит от w ₂, и наоборот для y ₂ и w ₁, оптимизация может проводиться независимо для каждой цели.

Я рассмотрел один образец здесь для иллюстрации, но ничего больше не меняется с большим количеством образцов. Вы можете написать это для себя. Добавление штрафного члена в форму | w ₁ | или | w ₂ | также не меняет это, так как все равно не будет никакой зависимости w ₂ от y ₁ или наоборот для y ₂ и w ₁.

Итак, это своего рода доказательство того, что вы получите C- (с щедрым профессором в этом). Поскольку это отражает sklearn, ручная реализация независимых регрессий и встроенная поддержка множественного вывода дадут тот же результат. Итак, давайте проверим это с помощью некоторого кода (я использую py2.7 с Jupyter):

Вещи, которые нам нужны

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

Настройте данные

## set up some test data
# N 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

Использовать встроенную поддержку множественного вывода

#%%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 такая же, как и наша ручная оценка. Тем не менее, встроенный намного быстрее, поскольку он использует матричную алгебру, чтобы решить все это один раз, в отличие от цикла, который я использовал здесь (для матричной формулировки регуляризации Риджа см. Вики о регуляризации Тихонова). Вы можете проверить это сами, раскомментировав магию %%timeit выше)

Q2: Как мы можем использовать отдельные параметры альфа-регуляризации для каждой регрессии?

A2: sklearn Ridge принимает различные регуляризации для каждого выхода (цели). Например, продолжая приведенный выше код, использовать разные альфы для каждого вывода:

# 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 независимых регрессий, каждая из которых имеет свою собственную альфа:

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() идет с факторизацией Холецкого (по крайней мере, для не разреженных данных), ища код для этого на github (_solve_cholesky в https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/linear_model/ridge.py), я вижу, что когда альфы выбираются отдельно для каждой цели, sklearn на самом деле подходит им отдельно. Я не знаю, является ли это присущим Cholesky или просто реализацией (у меня такое ощущение, что это последнее).

Однако для различных решателей, например на основе SVD (_solve_svd()), кажется, что код уже включает различные альфа-значения в матричную формулировку проблемы, поэтому весь процесс запускается только один раз. Это означает, что когда альфы выбираются отдельно для каждого выхода и когда имеется много выходов, решатель SVD может быть намного быстрее.

Q3: Как я могу использовать GridSeachCV? Сдаю ли я матрицу возможных параметров регуляризации?

A3: Я не использовал встроенный поиск по сетке, потому что он не совсем соответствовал моей проблеме. Однако, с учетом вышеизложенного, это легко реализовать; просто получите несколько CV-сгибов, используя sklearn.model_selection.KFold() или аналогичный, а затем тренируйтесь для каждого сгиба, используя разные альфа:

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)]

Я написал это довольно спешно, так что проверьте его внимательно. Обратите внимание, что нам нужно использовать метрический модуль, поскольку встроенный Ridge.score() усредняет по всем целям, чего мы здесь не хотим. Используя опцию multioutput='raw_values', мы сохраняем необработанное значение для каждой цели.

Надеюсь, это поможет!