Wdrażanie regresji grzbietu: Wybieranie inteligentnej siatki dla

Wdrażam regresję Ridge'a w module Python / C i natrafiłem na ten „mały” problem. Chodzi o to, że chcę próbować efektywne stopnie swobody mniej więcej równomiernie rozmieszczone (jak wykres na stronie 65 w „Elementach uczenia statystycznego” ), tj. Próbka:

d f (λ) = \sum_{i = 1}^{p} \frac{d_{i}^{2}}{d_{i}^{2} + λ},

$\mathrm{df}(\lambda)=\sum_{i=1}^{p}\frac{d_i^2}{d_i^2+\lambda},$

d_{i}^{2}

$d_i^2$

X^{T} X

$X^TX$

d f (λ_{max}) \approx 0

$\mathrm{df}(\lambda_{\max})\approx 0$

d f (λ_{min}) = p

$\mathrm{df}(\lambda_{\min})=p$

λ_{max} = \sum_{i}^{p} d_{i}^{2} / c

$\lambda_{\max}=\sum_i^p d_i^2/c$

λ_{max} ≫ d_{i}^{2}

$\lambda_{\max} \gg d_i^2$

c

$c$

c = 0.1

$c=0.1$

λ_{min} = 0

$\lambda_{\min}=0$

Jak sugeruje tytuł, muszę próbkować $\lambda$ z $\lambda_{\min}$ do $\lambda_{\max}$ w pewnej skali, tak że $\mathrm{df}(\lambda)$ jest próbkowany (w przybliżeniu), powiedzmy, w $0.1$ przedziałów od $c$ do $p$ ... czy jest to łatwy sposób? Myślałem, że rozwiązując równanie $\mathrm{df}(\lambda)$ dla każdego $\lambda$ przy użyciu metody Newtona-Raphsona, ale to doda zbyt wiele iteracji, szczególnie gdy $p$ jest duże. Jakieś sugestie?

ridge-regression

— Néstor
źródło

Ta funkcja jest malejącą wypukłą racjonalną funkcją

λ \geq 0

$\lambda \geq 0$ . Korzenie, szczególnie jeśli zostaną wybrane ponad siatką dyadic, powinny być bardzo szybkie do znalezienia.

— kardynał

@ kardynał, prawdopodobnie masz rację. Jednak, jeśli to możliwe, chciałbym wiedzieć, czy istnieje jakaś „domyślna” siatka. Na przykład próbowałem uzyskać siatkę, wykonując , gdzie , i działał całkiem dobrze dla niektórych stopni swobody, ale jak , to wybuchło. To mnie zastanowiło, że może istniał jakiś fajny sposób na wybranie siatki dla , o co pytam. Jeśli tak nie jest, chętnie bym to wiedział (ponieważ mógłbym szczęśliwie zostawić metodę Newton-Rapson w moim kodzie, wiedząc, że „nie ma lepszego sposobu”).

λ = l o g (s) λ_{m a x} / l o g (s_{m a x})

$\lambda=log(s)\lambda_{max}/log(s_{max})$

s = (1, 2, . . ., s_{m a x})

$s=(1,2,...,s_{max})$

d f (λ) \to p

$df(\lambda)\to p$

λ

$\lambda$

— Néstor

Aby lepiej zrozumieć potencjalne trudności, jakie napotykasz, jakie są typowe i najgorsze wartości ? Czy jest coś, co wiesz a priori o rozkładzie wartości własnych?

p

$p$

— kardynał

@ kardynalne, typowe wartości w mojej aplikacji byłyby w zakresie od do

, ale chcę, aby było tak ogólne, jak to możliwe. O rozkładzie wartości własnych, tak naprawdę niewiele.

jest macierzą, która zawiera predyktory w swoich kolumnach, które nie zawsze są ortogonalne.

p

$p$

15

$15$

40

$40$

X

$X$

— Néstor

Newton-Raphson zwykle znajduje pierwiastki z dokładnością do

10^{- 12}

$10^{-12}$ w ciągu

3

$3$ do

4

$4$ kroków dla

p = 40

$p=40$ i małych wartości

d f (λ)

$df(\lambda)$ ; prawie nigdy więcej niż

6

$6$ kroków. W przypadku większych wartości czasami wymagane jest do

30

$30$ kroków. Ponieważ każdy krok wymaga obliczeń

O (p)

$O(p)$ , całkowita ilość obliczeń jest nieistotna. Rzeczywiście, liczba kroków nie wydaje się zależeć od

p

$p$ jeśli wybrano dobrą wartość początkową (wybieram tę, której

gdyby wszystkie

d_{i}

$d_i$ równe ich średniej).

— whuber

To długa odpowiedź . Podajmy tutaj jego krótką wersję.

Nie ma dobrego algebraicznego rozwiązania tego problemu ze znalezieniem roota, więc potrzebujemy algorytmu numerycznego.
Funkcja $\mathrm{df}(\lambda)$ ma wiele fajnych właściwości. Możemy je wykorzystać, aby stworzyć specjalną wersję metody Newtona dla tego problemu z gwarantowaną konwergencją monotoniczną z każdym rdzeniem.
Nawet martwy mózg, Rnieobecny przy jakichkolwiek próbach optymalizacji, może obliczyć siatkę o wielkości 100 przy w ciągu kilku sekund. Starannie napisany kod zmniejszyłby to o co najmniej 2–3 rzędy wielkości. $p = 100\,000$ C

Istnieją dwa schematy podane poniżej, aby zagwarantować monotoniczną konwergencję. Jeden wykorzystuje granice pokazane poniżej, które wydają się pomagać czasami w zachowaniu kroku Newtona lub dwóch.

Przykład : i jednolita siatka dla stopni swobody wielkości 100. Wartości własne są rozkładem Pareto, a zatem bardzo wypaczonym. Poniżej znajdują się tabele liczby kroków Newtona, aby znaleźć każdy katalog główny. $p = 100\,000$

# Table of Newton iterations per root.
# Without using lower-bound check.
  1  3  4  5  6 
  1 28 65  5  1 
# Table with lower-bound check.
  1  2  3 
  1 14 85

Ogólnie rzecz biorąc, nie będzie dla tego rozwiązania w formie zamkniętej , ale istnieje wiele struktur, które można wykorzystać do tworzenia bardzo skutecznych i bezpiecznych rozwiązań przy użyciu standardowych metod wyszukiwania root.

Zanim zagłębimy się zbyt głęboko w rzeczy, zbierzmy pewne właściwości i konsekwencje funkcji

d f (λ) = \sum_{i = 1}^{p} \frac{d_{i}^{2}}{d_{i}^{2} + λ} .

$\newcommand{\df}{\mathrm{df}} \df(\lambda) = \sum_{i=1}^p \frac{d_i^2}{d_i^2 + \lambda} \>.$

Właściwość 0 : jest racjonalną funkcją $\df$ $\lambda$ . (Wynika to z definicji.)
Konsekwencja 0 : Nie będzie ogólnego rozwiązania algebraicznego dla znalezienia pierwiastka . Wynika to z faktu, że istnieje równoważny wielomianowy problem ze znalezieniem pierwiastka stopnia a więc jeśli nie jest bardzo małe (tj. Mniej niż pięć), nie będzie ogólnego rozwiązania. Potrzebujemy więc metody numerycznej. $\df(\lambda) - y = 0$ $p$ $p$

Właściwość 1 : Funkcja jest wypukła i maleje w . (Weź pochodne.) Konsekwencja 1 (a) : Algorytm znajdowania pierwiastka Newtona zachowa się bardzo ładnie w tej sytuacji. Niech będą pożądanymi stopniami swobody, a . Konsekwencja 1 (b) : Ponadto, jeśli zaczynamy od , wówczas pierwszy $\df$ $\lambda \geq 0$
$y$ $\lambda_0$ odpowiednim pierwiastkiem, tj. . W szczególności, jeśli zaczniemy od dowolnej wartości początkowej (więc, ), to sekwencja iteracji kroku Newtona zbiegnie się monotonicznie do unikalne rozwiązanie $y = \df(\lambda_0)$ $\lambda_1 < \lambda_0$ $\df(\lambda_1) > y$ $\lambda_1,\lambda_2,\ldots$ $\lambda_0$
$\lambda_1 > \lambda_0$ krok dałby , skąd monotonicznie wzrośnie do rozwiązania na podstawie poprzedniej konsekwencji (patrz zastrzeżenie poniżej). Intuicyjnie ten ostatni fakt wynika, ponieważ jeśli zaczniemy na prawo od nasady, pochodna jest „zbyt” płytka z powodu wypukłości więc pierwszy krok Newtona zabierze nas gdzieś na lewo od nasady. Uwaga: Ponieważ nie jest ogólnie wypukłe dla ujemnego $\lambda_2 \leq \lambda_0$ $\df$ $\df$ $\lambda$ , to stanowi silny powód, aby preferować rozpoczęcie od lewej strony żądanego katalogu głównego. W przeciwnym razie musimy dokładnie sprawdzić, czy krok Newtona nie spowodował ujemnej wartości szacowanego pierwiastka, co może umieścić nas gdzieś w niep wypukłej części . Konsekwencja 1 (c) : Po znalezieniu katalogu głównego dla niektórych a następnie szukamy katalogu głównego dla niektórych , używając taki sposób, że jako nasze początkowe przypuszczenie gwarantuje, że zaczynamy na lewo od drugiego katalogu głównego. Zatem nasza konwergencja jest odtąd monotoniczna. $\df$
$y_1$ $y_2 < y_1$ $\lambda_1$ $\df(\lambda_1) = y_1$

Właściwość 2 : Istnieją rozsądne granice zapewniające „bezpieczne” punkty początkowe. Używając argumentów wypukłości i nierówności Jensena, mamy następujące granice Konsekwencja 2 : To mówi nam, że root spełniający posłuszny Tak więc, do wspólnej stałej, rdzeń między środkami harmonicznymi i arytmetycznymi .

\frac{p}{1 + \frac{λ}{p} \sum d_{i}^{- 2}} \leq d f (λ) \leq \frac{p \sum_{i} d_{i}^{2}}{\sum_{i} d_{i}^{2} + p λ} .

$\frac{p}{1+ \frac{\lambda}{p}\sum d_i^{-2}} \leq \df(\lambda) \leq \frac{p \sum_i d_i^2}{\sum_i d_i^2 + p \lambda} \>.$

λ_{0}

$\lambda_0$

d f (λ_{0}) = y

$\df(\lambda_0) = y$

\begin{matrix} (⋆) & \frac{1}{\frac{1}{p} \sum_{i} d_{i}^{- 2}} (\frac{p - y}{y}) \leq λ_{0} \leq (\frac{1}{p} \sum_{i} d_{i}^{2}) (\frac{p - y}{y}) . \end{matrix}

$\frac{1}{\frac{1}{p}\sum_i d_i^{-2}}\left(\frac{p - y}{y}\right) \leq \lambda_0 \leq \left(\frac{1}{p}\sum_i d_i^2\right) \left(\frac{p - y}{y}\right) \>. \tag{$\star$}$

d_{i}^{2}

$d_i^2$

Zakłada się, że dla wszystkich . Jeśli tak nie jest, to ta sama granica obowiązuje, biorąc pod uwagę tylko dodatnie i zastępując liczbą dodatnich . NB : Ponieważ przy założeniu, że wszystkie , to , skąd granice są zawsze nietrywialne (np. Dolna granica jest zawsze nieujemna). $d_i > 0$ $i$ $d_i$ $p$ $d_i$ $\df(0) = p$ $d_i > 0$ $y \in (0,p]$

Oto wykres „typowego” przykładu z . Nałożyliśmy siatkę o rozmiarze 10 dla stopni swobody. Są to poziome linie na wykresie. Pionowe zielone linie odpowiadają dolnej granicy w . $\df(\lambda)$ $p = 400$ $(\star)$

Example dof plot with grid and bounds

Algorytm i przykładowy kod R.

Bardzo skutecznym algorytmem posiadającym siatkę pożądanych stopni swobody w jest sortowanie ich w malejącej kolejności, a następnie sekwencyjne znajdowanie pierwiastka każdego z nich, przy użyciu poprzedniego pierwiastka jako punktu początkowego dla następujących elementów 1. Możemy to udoskonalić, sprawdzając, czy każdy pierwiastek jest większy niż dolna granica dla następnego pierwiastka, a jeśli nie, możemy zamiast tego rozpocząć następną iterację od dolnej granicy. $y_1, \ldots y_n$ $(0,p]$

Oto przykładowy kod Rbez żadnych prób jego optymalizacji. Jak widać poniżej, wciąż jest dość szybki, choć R- mówiąc grzecznie - przerażająco, okropnie, strasznie wolno w pętlach.

# Newton's step for finding solutions to regularization dof.

dof <- function(lambda, d) { sum(1/(1+lambda / (d[d>0])^2)) }
dof.prime <- function(lambda, d) { -sum(1/(d[d>0]+lambda / d[d>0])^2) }

newton.step <- function(lambda, y, d)
{ lambda - (dof(lambda,d)-y)/dof.prime(lambda,d) }

# Full Newton step; Finds the root of y = dof(lambda, d).
newton <- function(y, d, lambda = NA, tol=1e-10, smart.start=T)
{
    if( is.na(lambda) || smart.start )
        lambda <- max(ifelse(is.na(lambda),0,lambda), (sum(d>0)/y-1)/mean(1/(d[d>0])^2))
    iter <- 0
    yn   <- Inf
    while( abs(y-yn) > tol )
    {
        lambda <- max(0, newton.step(lambda, y, d)) # max = pedantically safe
        yn <- dof(lambda,d)
        iter = iter + 1
    }
    return(list(lambda=lambda, dof=y, iter=iter, err=abs(y-yn)))
}

Poniżej znajduje się końcowy pełny algorytm, który pobiera siatkę punktów i wektor ( $d_i$ nie !). $d_i^2$

newton.grid <- function(ygrid, d, lambda=NA, tol=1e-10, smart.start=TRUE)
{
    p <- sum(d>0)
    if( any(d < 0) || all(d==0) || any(ygrid > p) 
        || any(ygrid <= 0) || (!is.na(lambda) && lambda < 0) )
        stop("Don't try to fool me. That's not nice. Give me valid inputs, please.")
    ygrid <- sort(ygrid, decreasing=TRUE)
    out    <- data.frame()
    lambda <- NA
    for(y in ygrid)
    {
        out <- rbind(out, newton(y,d,lambda, smart.start=smart.start))
        lambda <- out$lambda[nrow(out)]
    }
    out
}

Przykładowe wywołanie funkcji

set.seed(17)
p <- 100000
d <- sqrt(sort(exp(rexp(p, 10)),decr=T))
ygrid <- p*(1:100)/100
# Should take ten seconds or so.
out <- newton.grid(ygrid,d)

— kardynał
źródło

Podoba mi się pytanie, dzięki czemu mogę wrócić do tej odpowiedzi. Dziękuję za opublikowanie tej szczegółowej analizy, kardynał.

— Makro

Niesamowita odpowiedź :-), wielkie dzięki za sugestie ORAZ odpowiedź.

— Néstor

Ponadto istnieje kilka metod, które skutecznie obliczą pełną ścieżkę regularyzacji:

Powyżej są wszystkie pakiety R, ponieważ używasz Pythona, scikit-learn zawiera implementacje dla ridge, lasso i elastycznej siatki.

— wrz
źródło

olsFunkcji w R rmspakiecie można użyć optymalizacji numerycznej, aby znaleźć optymalne kary za pomocą skutecznej AIC. Ale musisz podać maksymalną karę, która nie zawsze jest łatwa.

— Frank Harrell,

Możliwą alternatywą według poniższego źródła wydaje się być:

Rozwiązanie w postaci zamkniętej: $df(\lambda) = tr(X(X^{\top} X + \lambda I_{p})^{-1}X^{\top})$

Jeśli używasz równania normalnego jako solwera lub oszacowanie wariancji-kowariancji, powinieneś już obliczyć . Podejście to działa najlepiej, jeśli szacujesz współczynniki dla różnych . $(X^{\top}X + \lambda I_{p})^{-1}$ $\lambda$

Źródło: https://onlinecourses.science.psu.edu/stat857/node/155

— José Bayoán Santiago Calderón
źródło