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?

— Baron
źródło

Są tutaj dwa dość ortogonalne pytania, osobne pytanie rozważałbym wybranie „czy mogę zastosować grzbiet w kategoriach najmniejszych kwadratów”.

— Matthew Drury

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:

Użyłem (praktycznie) znormalizowanych zmiennych predykcyjnych. Będziesz musiał sam wyjaśnić normalizację.
To samo dotyczy nienormalizacji przechwytywania.
Ż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ń.
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ść.
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.
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 .
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

— usεr11852
źródło

Jest to nieco mylące. Nieujemna regresja kalenicy jest trywialna do wdrożenia: można przepisać regresję kalenicy jak zwykle regresję na rozszerzonych danych (patrz komentarze do stats.stackexchange.com/questions/203687 ), a następnie zastosować procedury regresji nieujemnej.

— ameba

Zgadzam się, że jest to łatwe do wdrożenia (+1 do tego). (Przegłosowałem wcześniej twój i komentarz Glena także w drugim wątku). Pytanie brzmi, dlaczego nie jest wdrażane, nie jeśli jest to trudne. W związku z tym mocno podejrzewam, że bezpośrednie sformułowanie tego zadania NNRR jest problemem prostszym niż sformułowanie go jako rozszerzenia regresji danych, a następnie użycie Quad. Żarcie. optymalizacja w celu rozwiązania tej regresji. Nie powiedziałem tego w mojej odpowiedzi, ponieważ przedsięwziąłoby to część implementacyjną.

— usεr11852

Lub po prostu napisz to w stanie.

— Sycorax mówi Przywróć Monikę

Ah, dobrze; Zrozumiałem Q jako głównie z pytaniem, jak wykonać nieujemny grzbiet (i tylko pytaniem, dlaczego nie zostało to zrealizowane przejściowo); Zredagowałem nawet, aby umieścić to w tytule. W każdym razie, jak to zrobić, wydaje mi się bardziej interesujące pytanie. Jeśli potrafisz zaktualizować swoją odpowiedź objaśnieniami, jak wdrożyć nieujemny grzbiet, myślę, że będzie bardzo przydatny dla przyszłych czytelników (i chętnie głosuję :).

— ameba

Fajnie, zrobię to później (nie zauważyłem nowego tytułu, przepraszam za to). Prawdopodobnie podam implementację pod kątem OLS / pseudo-obserwacji, więc odpowiadamy również na inne pytanie.

— usεr11852

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

— rep_ho
źródło

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

{minimize}_{x} {‖ A x - y ‖}_{2}^{2} + λ ‖ x ‖_{2}^{2} s.t. x > 0

$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda \| x \|^2_2 \,\,\,\,\text{s.t. }x>0$

jest równa:

{minimize}_{x} {‖ A x - y ‖}_{2}^{2} + λ x^{⊤} I x s.t. x > 0

$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda x^\top I x \,\,\,\,\text{s.t. }x>0$

z nieco większą algebrą:

{minimize}_{x} x^{T} (A^{T} A + λ I) x + {(- 2 A^{T} y)}^{T} x s.t. x > 0

$\text{minimize}_{x}\,\,\,\,x^{T}\left(A^{T}A+ \lambda I \right)x+\left(-2A^{T}y\right)^{T}x\,\,\,\,\text{s.t. }x>0$

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.

— Charlie Parker
źródło

Czy linia c = - Nie czyta c = A? Myślę, że jest to poprawne, choć należy zauważyć, że rozwiązanie nieco różni się od scipy.optimize.nnls (newMatX, newVecY), gdzie newMatX jest wierszem X powiększonym o macierz diagonalną z sqrt (lambda) wzdłuż przekątnej, a NewVecY to Y powiększony zerami nvar. Myślę, że wspomniane rozwiązanie jest prawidłowe ...

— Tom Wenseleers