Я сделаю это, так как я немного разбираюсь в этом для своей собственной работы. Я разделю вопрос на части, чтобы вы могли смотреть только на те, которые вас интересуют:
1:
Независима ли регрессия для каждой цели (или выходной) в множественном выходе Ридж регрессии?
A1:
Я думаю, что типичная линейная регрессия с несколькими выходами с М выходами равна так же, как М независимая линейная регрессия с одним выходом. Я думаю, что это так, поскольку выражение для обычных наименьших квадратов для случая множественного выхода совпадает с выражением для (суммы) M независимых, единичных выходных случаев. Чтобы мотивировать это, давайте рассмотрим глупый, двумерный случай вывода без регуляризации.
Рассмотрим два входных вектора столбцов x 1 и x 2 с соответствующими весовыми векторами w 1 и ш 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 ш 1 T ) 2 + (у 2 - x 2 ш 2 T ) 2
Мы можем видеть, что это просто сумма квадратов ошибок двух независимых регрессий. Теперь, чтобы оптимизировать, мы дифференцируем по весам и устанавливаем в ноль. Поскольку y 1 не зависит w 2 , и наоборот для y 2 и w 1 , оптимизация может проводиться независимо для каждой цели.
Я рассмотрел один образец здесь для иллюстрации, но с большим количеством образцов ничего не меняется. Вы можете написать это для себя. Добавление штрафного члена в виде | w 1 | или | ш 2 | также не меняет это, так как все равно не будет зависимости w 2 от y 1 , или наоборот для y 2 и ш 1 .
Хорошо, вот какое доказательство, что вы получите 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
Использовать встроенную поддержку множественного вывода 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 такая же, как и наша ручная оценка. Тем не менее, встроенный работает на 1155 * намного быстрее, поскольку он использует матричную алгебру, чтобы решить все это один раз, в отличие от цикла, который я использовал здесь (для матричной формулировки регуляризации Риджа см. Вики о регуляризации Тихонова). Вы можете проверить это сами, раскомментировав магию %% timeit выше)
Q2:
Как мы можем использовать отдельные параметры альфа-регуляризации для каждой регрессии?
A2:
Склеарн Ридж принимает разные регуляризации для каждого выхода (цели). Например, продолжая приведенный выше код, чтобы использовать разные альфа для каждого вывода:
# 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 может быть намного быстрее.
В3: Как мне использовать GridSeachCV? Сдаю ли я матрицу возможных параметров регуляризации?
* * 1185 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', мы сохраняем необработанное значение для каждой цели.
Надеюсь, это поможет!