Wdrażanie regresji grzbietu: Wybieranie inteligentnej siatki dla

17

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:

df(λ)=i=1pdi2di2+λ,
di2XTXdf(λmax)0df(λmin)=pλmax=ipdi2/cλmaxdi2cc=0.1λmin=0

Jak sugeruje tytuł, muszę próbkować λ z λmin do λmax w pewnej skali, tak że df(λ) 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 df(λ) dla każdego λ przy użyciu metody Newtona-Raphsona, ale to doda zbyt wiele iteracji, szczególnie gdy p jest duże. Jakieś sugestie?

Néstor
źródło
1
Ta funkcja jest malejącą wypukłą racjonalną funkcją λ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”). λ=log(s)λmax/log(smax)s=(1,2,...,smax)df(λ)pλ
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
kardynał
@ kardynalne, typowe wartości w mojej aplikacji byłyby w zakresie od do 40 , ale chcę, aby było tak ogólne, jak to możliwe. O rozkładzie wartości własnych, tak naprawdę niewiele. X jest macierzą, która zawiera predyktory w swoich kolumnach, które nie zawsze są ortogonalne. p1540X
Néstor
1
Newton-Raphson zwykle znajduje pierwiastki z dokładnością do 1012 w ciągu 3 do 4 kroków dla p=40 i małych wartości df(λ) ; prawie nigdy więcej niż 6 kroków. W przypadku większych wartości czasami wymagane jest do 30 kroków. Ponieważ każdy krok wymaga obliczeń O(p) , całkowita ilość obliczeń jest nieistotna. Rzeczywiście, liczba kroków nie wydaje się zależeć od p jeśli wybrano dobrą wartość początkową (wybieram tę, której użyłbyś, gdyby wszystkie d_i byłydi równe ich średniej).
whuber

Odpowiedzi:

19

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 df(λ) 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=100000C

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=100000

# 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

df(λ)=i=1pdi2di2+λ.

Właściwość 0 : jest racjonalną funkcjądfλ . (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(λ)y=0pp

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 pierwszydfλ0
yλ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ązaniey=df(λ0)λ1<λ0df(λ1)>yλ1,λ2,λ0
λ1>λ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λ2λ0dfdfλ, 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
y1y2<y1λ1df(λ1)=y1

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 .

p1+λpdi2df(λ)pidi2idi2+pλ.
λ0df(λ0)=y
()11pidi2(pyy)λ0(1pidi2)(pyy).
di2

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).di>0idipdidf(0)=pdi>0y(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(λ)p=400()

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.y1,yn(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 (di nie !).di2

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
1

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

  1. GPS
  2. glmnet
  3. gcdnet

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

wrz
źródło
1
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,
0

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

Rozwiązanie w postaci zamkniętej: df(λ)=tr(X(XX+λIp)1X)

Jeśli używasz równania normalnego jako solwera lub obliczasz oszacowanie wariancji-kowariancji, powinieneś już obliczyć . Podejście to działa najlepiej, jeśli szacujesz współczynniki dla różnych λ .(XX+λIp)1λ

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

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