Bagaimana cara melakukan regresi ridge non-negatif?

Jawaban yang agak anti-iklim untuk " Apakah ada yang tahu mengapa ini? " Adalah bahwa tidak ada yang cukup peduli untuk menerapkan rutinitas regresi ridge non-negatif. Salah satu alasan utama adalah bahwa orang-orang sudah mulai menerapkan rutin jaring elastis non-negatif (misalnya di sini dan di sini ). Jaring elastis mencakup regresi punggungan sebagai kasus khusus (satu dasarnya mengatur bagian LASSO memiliki bobot nol). Karya-karya ini relatif baru sehingga belum dimasukkan dalam scikit-learn atau paket penggunaan umum serupa. Anda mungkin ingin menanyakan kode ini kepada pembuat makalah ini.

EDIT:

Seperti @amoeba dan saya bahas di komentar implementasi sebenarnya ini relatif sederhana. Katakanlah seseorang memiliki masalah regresi berikut untuk:

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

di mana dan keduanya normals standar seperti: $x_1$ $x_2$ . Perhatikan saya menggunakan variabel prediktor standar sehingga saya tidak harus menormalisasi setelahnya. Untuk kesederhanaan saya juga tidak menyertakan intersep. Kita dapat segera menyelesaikan masalah regresi ini dengan menggunakan regresi linier standar. Jadi di R itu harus seperti ini: $x_p \sim N(0,1)$

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

Perhatikan baris terakhir. Hampir semua rutinitas regresi linier menggunakan dekomposisi QR untuk memperkirakan . Kami ingin menggunakan hal yang sama untuk masalah regresi ridge kami. Pada titik ini baca posting ini oleh @whuber; kami akan menerapkan prosedur ini dengan tepat . Singkatnya, kita akan menambah matriks desain asli $\beta$ $X$ dengan matriks diagonal dan vektor respons kamidengannol. Dengan cara itu kita akan dapat kembali mengungkapkan masalah regresi ridge aslisebagai mana $\sqrt{\lambda}I_p$ $y$ $p$ $(X^TX + \lambda I)^{-1} X^Ty$ $(\bar{X}^T\bar{X})^{-1} \bar{X}^T\bar{y}$ $\bar{}$ melambangkan versi augmented. Lihat slide 18-19 dari catatan ini juga untuk kelengkapannya, saya menemukannya cukup mudah. Jadi dalam R kami ingin beberapa hal berikut:

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

dan itu berhasil. OK, jadi kami mendapat bagian regresi ridge. Kita dapat menyelesaikannya dengan cara lain, kita dapat memformulasikannya sebagai masalah optimisasi di mana jumlah sisa kuadrat adalah fungsi biaya dan kemudian mengoptimalkannya, yaitu. $\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$ . Cukup yakin kita bisa melakukan itu:

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

yang seperti yang diharapkan lagi berfungsi. Jadi sekarang kita hanya ingin: mana $\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

yang menunjukkan bahwa tugas regresi ridge non-negatif yang asli dapat diselesaikan dengan merumuskan kembali sebagai masalah optimisasi terbatas yang sederhana. Beberapa peringatan:

Saya menggunakan (secara praktis) variabel prediktor yang dinormalisasi. Anda harus memperhitungkan normalisasi sendiri.
Hal yang sama berlaku untuk non normalisasi intersep.
Saya menggunakan optim 's L-BFGS-B argumen. Ini adalah pemecah vanilla R paling yang menerima batas. Saya yakin Anda akan menemukan lusinan pemecah yang lebih baik.
Secara umum kendala linear kuadrat-terkecil masalah diajukan sebagai tugas optimasi kuadrat . Ini adalah kerja keras untuk posting ini tetapi perlu diingat bahwa Anda bisa mendapatkan kecepatan yang lebih baik jika diperlukan.
Seperti disebutkan dalam komentar, Anda dapat melewati regresi ridge sebagai bagian augmented-linear-regression dan secara langsung menyandikan fungsi biaya ridge sebagai masalah optimisasi. Ini akan menjadi jauh lebih sederhana dan posting ini jauh lebih kecil. Demi argumen saya menambahkan solusi kedua ini juga.
Saya tidak sepenuhnya berbicara dalam Python tetapi pada dasarnya Anda dapat meniru pekerjaan ini dengan menggunakan linalg.solve NumPy dan fungsi optimalisasi SciPy .
$\lambda$

Kode untuk poin 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
sumber

Ini agak menyesatkan. Regresi ridge non-negatif adalah sepele untuk diterapkan: seseorang dapat menulis ulang regresi ridge sebagai regresi biasa pada data yang diperluas (lihat komentar untuk stats.stackexchange.com/questions/203687 ) dan kemudian gunakan rutinitas regresi non-negatif.

amoeba

Saya setuju itu mudah diterapkan (+1 untuk itu). (Saya membenarkan komentar Anda sebelumnya dan komentar Glen di utas lainnya juga). Pertanyaannya adalah mengapa tidak diterapkan, tidak jika itu sulit. Dalam hal ini, saya sangat curiga bahwa secara langsung merumuskan tugas NNRR ini masalah optimasi bahkan lebih sederhana yang pertama merumuskannya sebagai perpanjangan data regresi dan kemudian menggunakan Quad. Prog. optimisasi untuk menyelesaikan regresi ini. Saya tidak mengatakan ini dalam jawaban saya karena itu akan berani di bagian implementasi.

usεr11852

Atau tulis saja dalam stan.

Sycorax berkata Reinstate Monica

Ah, oke; Saya mengerti Q terutama bertanya bagaimana melakukan punggungan non-negatif (dan hanya bertanya mengapa itu tidak diterapkan secara sepintas); Saya bahkan mengedit untuk memasukkan ini ke dalam judul. Bagaimanapun, bagaimana melakukannya menurut saya menjadi pertanyaan yang lebih menarik. Jika Anda dapat memperbarui jawaban Anda dengan penjelasan tentang bagaimana menerapkan punggungan non-negatif, saya pikir itu akan sangat berguna bagi pembaca masa depan (dan saya akan senang untuk mengungguli :).

amoeba

Keren, saya akan melakukannya nanti (saya tidak melihat judul baru, maaf tentang itu). Saya mungkin akan memberikan implementasi dalam hal OLS / pengamatan semu jadi kami juga menjawab pertanyaan lainnya.

usεr11852

Ingatlah bahwa kami mencoba untuk memecahkan:

{memperkecil}_{x} {‖ SEBUAH x - y ‖}_{2}^{2} + λ ‖ x ‖_{2}^{2} st x > 0

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

setara dengan:

{memperkecil}_{x} {‖ SEBUAH x - y ‖}_{2}^{2} + λ x^{⊤} saya x st x > 0

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

dengan beberapa aljabar:

{memperkecil}_{x} x^{T} ({SEBUAH}^{T} SEBUAH + λ saya) x + {(- 2 {SEBUAH}^{T} y)}^{T} x st 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$

Solusi dalam pseudo-python adalah melakukan:

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

$K$ $x^\top R_k x$

untuk jawaban yang sedikit lebih umum.

Charlie Parker
sumber

Haruskah baris c = - A'y tidak membaca c = A'y? Saya pikir ini benar, meskipun orang harus mencatat bahwa solusinya sedikit berbeda dari scipy.optimize.nnls (newMatX, newVecY), di mana newMatX adalah baris X yang ditambah dengan matriks diagonal dengan sqrt (lambda) sepanjang diagonal dan NewVecY adalah Y ditambah dengan nol nvar. Saya pikir solusi yang Anda sebutkan adalah yang benar ...

Tom Wenseleers

Bagaimana cara melakukan regresi ridge non-negatif?

Jawaban: