Jak wykonać nieujemną regresję kalenicy?

Jak wykonać nieujemną regresję kalenicy? Lasso nieujemne jest dostępne w scikit-learn, ale dla grzbietu nie mogę wymusić nieujemności bety i rzeczywiście otrzymuję współczynniki ujemne. Czy ktoś wie, dlaczego tak jest?

Czy mogę również zastosować grzbiet w kategoriach zwykłych najmniejszych kwadratów? Przeniesiono to do innego pytania: Czy mogę wdrożyć regresję kalenicy w odniesieniu do regresji OLS?

Raczej anty-klimatyczna odpowiedź na „ Czy ktoś wie, dlaczego tak jest? ” Polega na tym, że po prostu nikomu nie zależy na wdrożeniu nieujemnej regresji grzbietu. Jednym z głównych powodów jest to, że ludzie już zaczęli wdrażać nieujemne elastyczne procedury sieciowe (na przykład tutaj i tutaj ). Elastyczna siatka zawiera regresję kalenicy jako specjalny przypadek (jeden zasadniczo ustawia część LASSO na zerową wagę). Prace te są stosunkowo nowe, więc nie zostały jeszcze włączone do scikit-learn lub podobnego pakietu ogólnego zastosowania. Możesz poprosić autorów tych artykułów o kod.

EDYTOWAĆ:

Ponieważ @amoeba i ja omawialiśmy te komentarze, faktyczna implementacja tego jest względnie prosta. Powiedzmy, że mamy następujący problem z regresją:

$y = 2 x_1 - x_2 + \epsilon, \qquad \epsilon \sim N(0,0.2^2)$

gdzie $x_1$ i $x_2$ są standardowymi normami, takimi jak: $x_p \sim N(0,1)$. Zauważ, że używam standardowych zmiennych predykcyjnych, więc nie muszę później normalizować. Dla uproszczenia nie dołączam również przechwytywania. Możemy natychmiast rozwiązać ten problem regresji za pomocą standardowej regresji liniowej. Więc w R powinno być coś takiego:

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

Zwróć uwagę na ostatni wiersz. Prawie wszystkie procedury regresji liniowej wykorzystują rozkład QR do oszacowania$\beta$. Chcielibyśmy użyć tego samego w naszym problemie z regresją kalenicy. W tym momencie przeczytaj ten post przez @whuber; będziemy wdrażać dokładnie tę procedurę. Krótko mówiąc, będziemy powiększać naszą oryginalną matrycę projektową$X$ z $\sqrt{\lambda}I_p$ macierz diagonalna i nasz wektor odpowiedzi $y$ z $p$zera W ten sposób będziemy mogli ponownie wyrazić pierwotny problem regresji kalenicy$(X^TX + \lambda I)^{-1} X^Ty$ tak jak $(\bar{X}^T\bar{X})^{-1} \bar{X}^T\bar{y}$ gdzie $\bar{}$symbolizuje wersję rozszerzoną. Sprawdź też kompletność slajdów 18–19 z tych notatek , uważam je za dość proste. Więc w R chcielibyśmy niektóre z następujących:

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

i to działa. OK, więc mamy część regresji grzbietu. Moglibyśmy rozwiązać w inny sposób, moglibyśmy sformułować go jako problem optymalizacji, w którym rezydualna suma kwadratów jest funkcją kosztu, a następnie zoptymalizować względem niej, tj.$\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$. Rzeczywiście możemy to zrobić:

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

który zgodnie z oczekiwaniami znów działa. Teraz chcemy tylko: gdzie . Jest to po prostu ten sam problem optymalizacji, ale ograniczony, aby rozwiązanie nie było ujemne.$\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$$\beta \geq 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

co pokazuje, że pierwotne nieujemne zadanie regresji grzbietu można rozwiązać, zmieniając formułę jako prosty ograniczony problem optymalizacji. Niektóre zastrzeżenia:

1. Użyłem (praktycznie) znormalizowanych zmiennych predykcyjnych. Będziesz musiał sam wyjaśnić normalizację.
2. To samo dotyczy nienormalizacji przechwytywania.
3. Że stosuje się optim„a L-BFGS-B argumentu. Jest to najbardziej waniliowy solver R, który akceptuje granice. Jestem pewien, że znajdziesz dziesiątki lepszych rozwiązań.
4. Zasadniczo problemy liniowe najmniejszych kwadratów są przedstawiane jako kwadratowe zadania optymalizacji . To przesada dla tego postu, ale pamiętaj, że w razie potrzeby możesz uzyskać większą prędkość.
5. Jak wspomniano w komentarzach, można pominąć regresję kalenicy jako część regresji rozszerzonej i liniowej i bezpośrednio zakodować funkcję kosztu kalenicy jako problem optymalizacji. Byłoby to o wiele prostsze, a ten post znacznie mniejszy. Dla argumentu dołączam również to drugie rozwiązanie.
6. Nie jestem w pełni konwersacyjny w Pythonie, ale w zasadzie możesz replikować tę pracę, używając linalg.solve NumPy i funkcji optymalizacji SciPy .
7. Aby wybrać hiperparametr itp., Po prostu wykonaj zwykły krok CV, który zrobiłbyś w każdym przypadku; nic się nie zmienia.$\lambda$

Kod dla punktu 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

Pakiet R glmnet, który wprowadza elastyczną siatkę, a zatem pozwala na to lasso i grzbiet. Dzięki parametrom lower.limitsi upper.limitsmożesz ustawić minimalną lub maksymalną wartość dla każdej masy osobno, więc jeśli ustawisz dolne limity na 0, wykona ona nieujemną siatkę elastyczną (lasso / grzbiet).

Istnieje również opakowanie Pythona https://pypi.python.org/pypi/glmnet/2.0.0

Przypomnijmy, że próbujemy rozwiązać:

jest równa:

z nieco większą algebrą:

Rozwiązaniem w pseudo-python jest po prostu:

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

patrz: W jaki sposób wykonuje się rzadkie nieujemne najmniejsze kwadraty za pomocą w postaci ?$K$$x^\top R_k x$

dla nieco bardziej ogólnej odpowiedzi.