Как выполнить неотрицательную ребристую регрессию?

10

Как выполнить неотрицательную ребристую регрессию? Доступно неотрицательное лассо scikit-learn, но для риджа я не могу навязать неотрицательность бета-версий, и действительно, я получаю отрицательные коэффициенты. Кто-нибудь знает, почему это?

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

Барон
источник
1
Здесь есть два совершенно ортогональных вопроса, и я хотел бы рассмотреть вопрос о том, можно ли выделить «я могу реализовать ребро в терминах наименьших квадратов» как отдельный вопрос.
Мэтью Друри

Ответы:

8

Весьма антиклиматический ответ на вопрос « Кто-нибудь знает, почему это так? » Заключается в том, что просто никто не заботится о том, чтобы реализовать процедуру отрицательной регрессии гребня. Одна из основных причин заключается в том, что люди уже начали внедрять подпрограммы неотрицательных эластичных сетей (например, здесь и здесь ). Эластичная сеть включает в себя регрессию гребня как особый случай (по существу, для детали LASSO задается нулевой вес). Эти работы являются относительно новыми, поэтому они еще не включены в scikit-learn или аналогичный пакет общего пользования. Возможно, вы захотите спросить авторов этих документов для кода.

РЕДАКТИРОВАТЬ:

Поскольку @amoeba и я обсуждали комментарии, фактическая реализация этого относительно проста. Скажем, у кого-то есть следующая проблема регрессии:

y=2x1x2+ϵ,ϵN(0,0.22)

где и оба являются стандартными нормалями, такими как: . Обратите внимание, что я использую стандартизированные переменные предиктора, поэтому мне не нужно нормализовать потом. Для простоты я также не включаю перехват. Мы можем немедленно решить эту проблему регрессии, используя стандартную линейную регрессию. Так что в R это должно быть примерно так:x 2 x pN ( 0 , 1 )x1x2xpN(0,1)

rm(list = ls()); 
library(MASS); 
set.seed(123);
N = 1e6;
x1 = rnorm(N)
x2 = rnorm(N)
y = 2 * x1 - 1 * x2 + rnorm(N,sd = 0.2)

simpleLR = lm(y ~ -1 + x1 + x2 )
matrixX = model.matrix(simpleLR); # This is close to standardised
vectorY = y
all.equal(coef(simpleLR), qr.solve(matrixX, vectorY), tolerance = 1e-7)  # TRUE

Обратите внимание на последнюю строку. Почти все процедуры линейной регрессии используют QR-разложение для оценки . Мы хотели бы использовать то же самое для нашей задачи регрессии гребня. На данный момент прочитайте этот пост @whuber; мы будем реализовывать именно эту процедуру. Короче говоря, мы будем расширять нашу оригинальную конструкцию матрицы с диагональная матрица и наш ответ вектор с нулей. Таким образом, мы сможем повторно выразить исходную проблему регрессии гребня как гдеX βXyp(XTX+λI) - 1 XTy( ˉ X T ˉ X ) - 1 ˉ X T ˉ y ¯λIpyp(XTX+λI)1XTy(X¯TX¯)1X¯Ty¯¯символизирует дополненную версию. Проверьте полноты слайдов 18-19 из этих заметок , я нашел их довольно простыми. Таким образом, в R мы хотели бы следующее:

myLambda = 100;  
simpleRR = lm.ridge(y ~ -1 + x1 + x2, lambda = myLambda)
newVecY = c(vectorY, rep(0, 2))
newMatX = rbind(matrixX, sqrt(myLambda) * diag(2))
all.equal(coef(simpleRR), qr.solve(newMatX, newVecY), tolerance = 1e-7)  # TRUE

и это работает. Итак, мы получили часть регрессии гребня. Мы могли бы решить по-другому, хотя, мы могли бы сформулировать это как задачу оптимизации, где остаточная сумма квадратов является функцией стоимости, а затем оптимизировать ее, т.е. . Конечно же, мы можем сделать это:minβ||y¯X¯β||22

myRSS <- function(X,y,b){ return( sum( (y - X%*%b)^2 ) ) }
bfgsOptim = optim(myRSS, par = c(1,1), X = newMatX, y= newVecY, 
                  method = 'L-BFGS-B')
all.equal(coef(simpleRR), bfgsOptim$par, check.attributes = FALSE, 
          tolerance = 1e-7) # TRUE

который, как и ожидалось, снова работает. Так что теперь мы просто хотим: где . Это просто та же самая проблема оптимизации, но ограниченная, так что решение неотрицательно.minβ||y¯X¯β||22β0

bfgsOptimConst = optim(myRSS, par = c(1,1), X=newMatX, y= newVecY, 
                       method = 'L-BFGS-B', lower = c(0,0))
all(bfgsOptimConst$par >=0)  # TRUE
(bfgsOptimConst$par) # 2.000504 0.000000

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

  1. Я использовал (практически) нормализованные предикторы. Вам нужно будет самостоятельно учесть нормализацию.
  2. То же самое касается и отсутствия нормализации перехвата.
  3. Я использовал optim«с L-BFGS-B аргумент. Это самый ванильный R решатель, который принимает границы. Я уверен, что вы найдете десятки лучших решателей.
  4. В общем случае линейные задачи наименьших квадратов представляют собой задачи квадратичной оптимизации . Это перебор для этого поста, но имейте в виду, что при необходимости вы можете получить лучшую скорость.
  5. Как упомянуто в комментариях, вы можете пропустить ребро-регрессию как часть расширенной линейной регрессии и напрямую закодировать функцию стоимости ребра в качестве задачи оптимизации. Это было бы намного проще, и этот пост значительно меньше. Ради аргумента я добавляю и это второе решение.
  6. Я не полностью разговорный Python , но по существу вы можете повторить эту работу, используя Numpy в linalg.solve и SciPy в Оптимизируйте функции.
  7. Чтобы выбрать гиперпараметр и т. Д., Вы просто делаете обычный CV-шаг, который вы сделали бы в любом случае; ничего не меняется.λ

Код для пункта 5:

myRidgeRSS <- function(X,y,b, lambda){ 
                return( sum( (y - X%*%b)^2 ) + lambda * sum(b^2) ) 
              }
bfgsOptimConst2 = optim(myRidgeRSS, par = c(1,1), X = matrixX, y = vectorY,
                        method = 'L-BFGS-B', lower = c(0,0), lambda = myLambda)
all(bfgsOptimConst2$par >0) # TRUE
(bfgsOptimConst2$par) # 2.000504 0.000000
usεr11852
источник
1
Это несколько вводит в заблуждение. Реализовать неотрицательную гребень тривиально: можно переписать гребень как обычную регрессию для расширенных данных (см. Комментарии к stats.stackexchange.com/questions/203687 ), а затем использовать процедуры неотрицательной регрессии.
амеба
Я согласен, что это просто реализовать (+1 к этому). (Я ранее проголосовал за ваш и комментарий Глена в другой ветке). Вопрос в том, почему это не реализовано, а не если это сложно. В связи с этим я сильно подозреваю, что непосредственно сформулировать эту задачу NNRR задача оптимизации еще проще, чем сначала сформулировать ее как расширяющуюся регрессию данных, а затем использовать Quad. Prog. оптимизация для решения этой регрессии. Я не сказал этого в своем ответе, потому что он рискнул бы в части реализации.
usεr11852
Или просто напишите это на стан.
Sycorax говорит восстановить Monica
Ах хорошо; Я понял, что Q в основном спрашивает, как сделать неотрицательный гребень (и только спрашивает, почему он не реализован попутно); Я даже отредактировал, чтобы поместить это в заголовок. В любом случае, как это сделать, мне кажется, более интересный вопрос. Если вы сможете обновить свой ответ объяснениями о том, как реализовать неотрицательный гребень, я думаю, что он будет очень полезен для будущих читателей (и я буду рад поднять голос :).
амеба
1
Круто, я сделаю это позже (я не заметил новый заголовок, извините за это). Я, вероятно, приведу реализацию в терминах OLS / псевдонаблюдений, поэтому мы ответим и на другой вопрос.
usεr11852
4

Пакет R glmnet, который реализует эластичную сеть и, следовательно, лассо и гребень, позволяет это. С помощью параметров lower.limitsи upper.limitsможно установить минимальное или максимальное значение для каждого веса в отдельности, поэтому, если вы установите нижние пределы на 0, он будет выполнять неотрицательную эластичную сеть (лассо / гребень).

Существует также оболочка Python https://pypi.python.org/pypi/glmnet/2.0.0

rep_ho
источник
2

Напомним, мы пытаемся решить:

minimizexAxy22+λx22s.t. x>0

эквивалентно:

minimizexAxy22+λxIxs.t. x>0

с еще немного алгебры:

minimizexxT(ATA+λI)x+(2ATy)Txs.t. x>0

Решение в псевдо-питоне просто сделать:

Q = A'A + lambda*I
c = - A'y
x,_ = scipy.optimize.nnls(Q,c)

см .: Как можно сделать разреженные неотрицательные наименьшие квадраты, используя регуляризаторов вида ?KxRkx

для чуть более общего ответа.

Чарли Паркер
источник
Должна ли строка c = - A'y не читать c = A'y? Я думаю, что это правильно, хотя следует отметить, что решение немного отличается от scipy.optimize.nnls (newMatX, newVecY), где newMatX - это строка X, дополненная диагональной матрицей с sqrt (лямбда) вдоль диагонали, а NewVecY - Y дополнен нвар нолями. Я думаю, что решение, которое вы упоминаете, является правильным ...
Том Wenseleers