Menerapkan regresi punggungan: Memilih kisi cerdas untuk

Saya menerapkan Ridge Regression dalam modul Python / C, dan saya telah menemukan masalah "kecil" ini. Idenya adalah bahwa saya ingin sampel derajat efektif kebebasan kurang lebih sama spasi (seperti plot pada halaman 65 pada "Elemen statistik Learning" ), yaitu, sampel:

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},$ mana

d_{i}^{2}

$d_i^2$ adalah nilai eigen dari matriks

X^{T} X

$X^TX$ , dari

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

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

d f (λ_{min}) = p

$\mathrm{df}(\lambda_{\min})=p$ . Cara mudah untuk mengatur batas pertama adalah untuk membiarkan

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

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

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

$\lambda_{\max} \gg d_i^2$ ), di mana

c

$c$ adalah konstanta kecil dan mewakili kira-kira tingkat kebebasan minimum yang ingin Anda sampel (misalnya

c = 0.1

$c=0.1$ ). Batas kedua tentu saja

λ_{min} = 0

$\lambda_{\min}=0$ .

Seperti judulnya, maka, saya harus sampel $\lambda$ dari $\lambda_{\min}$ ke $\lambda_{\max}$ dalam beberapa skala sehingga $\mathrm{df}(\lambda)$ adalah sampel (kurang-lebih), mengatakan, di $0.1$ interval dari $c$ ke $p$ ... apakah ada yang mudah cara untuk melakukan ini? Saya pikir menyelesaikan persamaan $\mathrm{df}(\lambda)$ untuk setiap $\lambda$ menggunakan metode Newton-Raphson, tetapi ini akan menambahkan iterasi yang terlalu banyak, khususnya ketika $p$ besar. Ada saran?

ridge-regression Néstor
sumber

Fungsi ini adalah penurunan fungsi rasional cembung

λ \geq 0

$\lambda \geq 0$ . Root, terutama jika dipilih di atas grid diad, harus sangat cepat ditemukan.

kardinal

@ kardinal, Anda mungkin benar. Namun, jika memungkinkan, saya ingin tahu jika ada beberapa kotak "default". Sebagai contoh, saya mencoba mendapatkan grid dengan melakukan

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

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

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

$s=(1,2,...,s_{max})$ , dan bekerja cukup baik untuk beberapa derajat kebebasan, tetapi sebagai

d f (λ) \to p

$df(\lambda)\to p$ , itu meledak. Ini membuat saya bertanya-tanya bahwa mungkin ada beberapa cara yang rapi untuk memilih kisi untuk

λ

$\lambda$ , yang saya tanyakan. Jika ini tidak ada, saya juga akan senang mengetahui (karena saya dapat meninggalkan metode Newton-Rapson dengan senang hati dalam kode saya mengetahui bahwa "tidak ada cara yang lebih baik").

Néstor

Untuk mendapatkan ide yang lebih baik tentang kesulitan potensial yang Anda temui, apa nilai tipikal dan terburuk dari

p

$p$ ? Apakah ada sesuatu yang Anda tahu apriori tentang distribusi eigen?

kardinal

@ cardinal, nilai-nilai khas

p

$p$ dalam aplikasi saya akan berkisar antara

15

$15$ hingga

40

$40$ , tapi saya ingin membuatnya sealami mungkin. Tentang distribusi nilai eigen, tidak terlalu banyak.

X

$X$ adalah matriks yang berisi prediktor di kolomnya, yang tidak selalu ortogonal.

Néstor

Newton-Raphson biasanya menemukan akar ke

10^{- 12}

$10^{-12}$ akurasi dalam

3

$3$ ke

4

$4$ langkah untuk

p = 40

$p=40$ dan nilai-nilai kecil

d f (λ)

$df(\lambda)$ ; hampir tidak pernah lebih dari

6

$6$ langkah. Untuk nilai yang lebih besar, kadang-kadang diperlukan hingga

30

$30$ langkah. Karena setiap langkah memerlukan perhitungan

O (p)

$O(p)$ , jumlah total perhitungan tidak penting. Memang, jumlah langkah tampaknya tidak tergantung pada

p

$p$ jika nilai awal yang baik dipilih (saya memilih salah satu yang akan Anda gunakan jika semua

d_{i}

$d_i$ sama dengan mean mereka).

whuber

Jawaban:

Ini jawaban yang panjang . Jadi, mari berikan versi cerita pendeknya di sini.

Tidak ada solusi aljabar yang bagus untuk masalah pencarian-akar ini, jadi kita perlu algoritma numerik.
Fungsi $\mathrm{df}(\lambda)$ memiliki banyak properti bagus. Kita dapat memanfaatkan ini untuk membuat versi khusus dari metode Newton untuk masalah ini dengan konvergensi monoton dijamin untuk setiap root.
Bahkan Rkode mati otak tidak ada upaya optimasi dapat menghitung kisi ukuran 100 dengan dalam beberapa detik. Kode yangditulis dengan hati-hatiakan mengurangi ini setidaknya 2-3 urutan besarnya. $p = 100\,000$ C

Ada dua skema yang diberikan di bawah ini untuk menjamin konvergensi monoton. Satu menggunakan batas yang ditunjukkan di bawah ini, yang tampaknya membantu menyelamatkan satu atau dua langkah Newton pada kesempatan.

Contoh : dan kotak seragam untuk derajat kebebasan ukuran 100. Nilai-nilai eigennya didistribusikan Pareto, karenanya sangat condong. Di bawah ini adalah tabel jumlah langkah Newton untuk menemukan setiap root. $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

Tidak akan ada solusi bentuk tertutup untuk ini , secara umum, tetapi ada adalah banyak struktur sekarang yang dapat digunakan untuk menghasilkan solusi yang sangat efektif dan aman menggunakan metode akar-menemukan standar.

Sebelum menggali terlalu dalam ke hal-hal, mari kita mengumpulkan beberapa sifat dan konsekuensi dari fungsi

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} \>.$

Properti 0 : adalah fungsi rasional . (Hal ini terlihat dari definisi.) Konsekuensi 0 : Tidak ada solusi aljabar umum akan ada untuk menemukan akar . Ini karena ada masalah penemuan akar polinomial setara dengan derajat dan jadi jika tidak terlalu kecil (yaitu kurang dari lima), tidak ada solusi umum yang akan ada. Jadi, kita membutuhkan metode numerik. $\df$ $\lambda$
$\df(\lambda) - y = 0$ $p$ $p$

Properti 1 : Fungsi yaitu cembung dan menurun pada . (Ambil turunannya.) Konsekuensi 1 (a) : Algoritma pencarian akar Newton akan berperilaku sangat baik dalam situasi ini. Biarkan menjadi derajat kebebasan yang diinginkan dan akar yang sesuai, yaitu, . Secara khusus, jika kita memulai dengan setiap nilai awal (jadi, $\df$ $\lambda \geq 0$
$y$ $\lambda_0$ $y = \df(\lambda_0)$ $\lambda_1 < \lambda_0$ ), maka urutan iterasi langkah-Newton akan konvergensecara monotonke solusi unik . Konsekuensi 1 (b): Selanjutnya, jika kita memulai dengan , makalangkahpertamaakan menghasilkan $\df(\lambda_1) > y$ $\lambda_1,\lambda_2,\ldots$ $\lambda_0$
$\lambda_1 > \lambda_0$ $\lambda_2 \leq \lambda_0$ , dari mana itu akan meningkat secara monoton ke solusi dengan konsekuensi sebelumnya (lihat peringatan di bawah). Secara intuitif, fakta terakhir ini mengikuti karena jika kita mulai di sebelah kanan akar, derivatif adalah "terlalu" dangkal karena busung dari dan langkah Newton pertama akan membawa kita ke suatu tempat ke kiri dari akar. NB Sejak adalah tidak di cembung umum untuk negatif , ini memberikan alasan yang kuat untuk memilih mulai ke kiri dari akar yang diinginkan. Jika tidak, kita perlu periksa bahwa langkah Newton tidak menghasilkan nilai negatif untuk perkiraan akar, yang dapat menempatkan kita di suatu tempat di bagian nonconvex dari . $\df$ $\df$ $\lambda$ $\df$
Konsekuensi 1 (c) : Setelah kami telah menemukan akar untuk beberapa dan kemudian mencari akar dari beberapa , menggunakan sehingga sebagai dugaan awal kami jaminan kami mulai di sebelah kiri dari root kedua. Jadi, konvergensi kami dijamin monoton dari sana. $y_1$ $y_2 < y_1$ $\lambda_1$ $\df(\lambda_1) = y_1$

Properti 2 : Ada batasan yang wajar untuk memberikan titik awal yang "aman". Menggunakan argumen konveksitas dan ketidaksetaraan Jensen, kita memiliki batasan berikut

\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} \>.$ Consequence 2: This tells us that the root

λ_{0}

$\lambda_0$ satisfying

d f (λ_{0}) = y

$\df(\lambda_0) = y$ obeys

\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$}$ So, up to a common constant, we've sandwiched the root in between the harmonic and arithmetic means of the

d_{i}^{2}

$d_i^2$ .

This assumes that $d_i > 0$ for all $i$ . If this is not the case, then the same bound holds by considering only the positive $d_i$ and replacing $p$ by the number of positive $d_i$ . NB: Since $\df(0) = p$ assuming all $d_i > 0$ , then $y \in (0,p]$ , whence the bounds are always nontrivial (e.g., the lower bound is always nonnegative).

Here is a plot of a "typical" example of $\df(\lambda)$ with $p = 400$ . We've superimposed a grid of size 10 for the degrees of freedom. These are the horizontal lines in the plot. The vertical green lines correspond to the lower bound in $(\star)$ .

Example dof plot with grid and bounds

An algorithm and some example R code

A very efficient algorithm given a grid of desired degrees of freedom $y_1, \ldots y_n$ in $(0,p]$ is to sort them in decreasing order and then sequentially find the root of each, using the previous root as the starting point for the following one. We can refine this further by checking if each root is greater than the lower bound for the next root, and, if not, we can start the next iteration at the lower bound instead.

Here is some example code in R, with no attempts made to optimize it. As seen below, it is still quite fast even though R is—to put it politely—horrifingly, awfully, terribly slow at loops.

# 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)))
}

Below is the final full algorithm which takes a grid of points, and a vector of the $d_i$ (not $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
}

Sample function call

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)

cardinal
sumber

Favoriting the question so I can refer back to this answer. Thanks for posting this detailed analysis, cardinal.

Macro

Amazing answer :-), thanks a lot cardinal for the suggestions AND the answer.

Néstor

In addition, a couple of methods exist that will calculate the complete regularization path efficiently:

The above are all R packages, as you are using Python, scikit-learn contains implementations for ridge, lasso and elastic net.

sebp
sumber

The ols function in the R rms package can use numerical optimization to find the optimum penalty using effective AIC. But you have to provide the maximum penalty which is not always easy.

Frank Harrell

A possible alternative according to the source below seems to be:

The closed form solution: $df(\lambda) = tr(X(X^{\top} X + \lambda I_{p})^{-1}X^{\top})$

Should you be using the normal equation as the solver or computing the variance-covariance estimate, you should already have computed $(X^{\top}X + \lambda I_{p})^{-1}$ . This approach works best if you are estimating the coefficients at the various $\lambda$ .

Source: https://onlinecourses.science.psu.edu/stat857/node/155

José Bayoán Santiago Calderón
sumber