Comment effectuer une régression de crête non négative?

La réponse plutôt anti-climatique à " Quelqu'un sait-il pourquoi? " Est que tout le monde ne se soucie pas suffisamment de mettre en œuvre une routine de régression de crête non négative. L'une des principales raisons est que les gens ont déjà commencé à mettre en œuvre des routines à filet élastique non négatif (par exemple ici et ici ). Le filet élastique inclut la régression de crête comme cas spécial (l'un a essentiellement défini la partie LASSO pour avoir une pondération nulle). Ces travaux sont relativement nouveaux et n'ont donc pas encore été intégrés à scikit-learn ou à un package similaire à usage général. Vous voudrez peut-être demander aux auteurs de ces articles le code.

ÉDITER:

Comme @amoeba et moi avons discuté des commentaires, la mise en œuvre réelle de ceci est relativement simple. Supposons que l'on ait le problème de régression suivant:

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

où et sont tous deux des normales standard telles que: . Remarquez que j'utilise des variables prédictives standardisées, je n'ai donc pas besoin de normaliser par la suite. Par souci de simplicité, je n'inclus pas d'interception non plus. Nous pouvons immédiatement résoudre ce problème de régression en utilisant une régression linéaire standard. Donc, dans R, cela devrait être quelque chose comme ceci: $x_1$ $x_2$ $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

Remarquez la dernière ligne. Presque toutes les routines de régression linéaire utilisent la décomposition QR pour estimer . Nous aimerions utiliser la même chose pour notre problème de régression de crête. À ce stade, lisez ce message de @whuber; nous mettrons en œuvre exactement cette procédure. En bref, nous allons augmenter notre matrice de conception originale avec une matrice diagonale et notre vecteur de réponse avec zéros. De cette façon, nous serons en mesure de ré-exprimer le problème de régression de crête d'origine as où le $\beta$ $X$ $\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{}$ symbolise la version augmentée. Vérifiez les diapositives 18-19 de ces notes aussi pour être complètes, je les ai trouvées assez simples. Donc, dans R, nous aimerions certains:

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

$\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$

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

qui, comme prévu, fonctionne à nouveau. Alors maintenant, nous voulons juste: où . Ce qui est tout simplement le même problème d'optimisation mais contraint pour que la solution ne soit pas négative. $\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

ce qui montre que la tâche de régression de crête non négative d'origine peut être résolue en reformulant comme un simple problème d'optimisation contraint. Quelques mises en garde:

J'ai utilisé des variables prédictives normalisées (pratiquement). Vous devrez vous-même rendre compte de la normalisation.
La même chose vaut pour la non normalisation de l'interception.
J'ai utilisé optiml' argument L-BFGS-B de . C'est le solveur R le plus vanillé qui accepte les limites. Je suis sûr que vous trouverez des dizaines de meilleurs solveurs.
En général, les problèmes de moindres carrés linéaires de contrainte sont posés comme des tâches d' optimisation quadratiques . C'est une exagération pour ce poste, mais gardez à l'esprit que vous pouvez obtenir une meilleure vitesse si nécessaire.
Comme mentionné dans les commentaires, vous pouvez ignorer la régression de crête en tant que partie de régression linéaire augmentée et encoder directement la fonction de coût de crête en tant que problème d'optimisation. Ce serait beaucoup plus simple et ce poste beaucoup plus petit. Par souci d'argument, j'ajoute également cette deuxième solution.
Je ne suis pas complètement la conversation en Python , mais essentiellement vous pouvez reproduire ce travail en utilisant des NumPy de linalg.solve et de SciPy optimiser les fonctions.
Pour choisir l'hyperparamètre etc., il vous suffit de faire l'étape CV habituelle que vous feriez dans tous les cas; rien ne change. $\lambda$

Code pour le point 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
la source

C'est quelque peu trompeur. La régression de crête non négative est triviale à mettre en œuvre: on peut réécrire la régression de crête comme une régression habituelle sur des données étendues (voir les commentaires sur stats.stackexchange.com/questions/203687 ), puis utiliser des routines de régression non négatives.

amoeba

Je suis d'accord qu'il est simple à mettre en œuvre (+1 à cela). (J'ai voté plus tôt le vôtre et le commentaire de Glen sur l'autre fil aussi). La question est de savoir pourquoi n'est pas mis en œuvre cependant, pas si c'est difficile. À ce sujet, je soupçonne fortement que la formulation directe de cette tâche NNRR un problème d'optimisation est encore plus simple que de le formuler d'abord comme une régression de données étendue, puis d'utiliser Quad. Programme. optimisation pour résoudre cette régression. Je n'ai pas dit cela dans ma réponse, car cela se risquerait dans la partie mise en œuvre.

usεr11852

Ou écrivez-le simplement en stan.

Sycorax dit Réintégrer Monica

Ah ok; J'ai compris le Q comme demandant principalement comment faire une crête non négative (et demandant seulement pourquoi il n'est pas implémenté en passant); J'ai même édité pour mettre cela dans le titre. En tout cas, comment faire cela me semble être une question plus intéressante. Si vous pouvez mettre à jour votre réponse avec des explications sur la façon de mettre en œuvre une crête non négative, je pense que cela sera très utile pour les futurs lecteurs (et je serai heureux de voter de manière positive :).

amoeba

Cool, je le ferai plus tard (je n'ai pas remarqué le nouveau titre, désolé pour ça). Je vais probablement donner l'implémentation en termes d'OLS / pseudo-observations donc nous répondons également à l'autre question.

usεr11852

Rappelons que nous essayons de résoudre:

{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$

est équivalent à:

{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$

avec un peu plus d'algèbre:

{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$

La solution en pseudo-python est simplement de faire:

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

voir: Comment peut-on faire des moindres carrés non négatifs clairsemés en utilisant régulariseurs de la forme ? $K$ $x^\top R_k x$

pour une réponse un peu plus générale.

Charlie Parker
la source

La ligne c = - A'y ne doit-elle pas lire c = A'y? Je pense que c'est correct, mais il faut noter que la solution est légèrement différente de scipy.optimize.nnls (newMatX, newVecY), où newMatX est une ligne X augmentée d'une matrice diagonale avec sqrt (lambda) le long de la diagonale et NewVecY est Y augmenté de zéros nvar. Je pense que la solution que vous mentionnez est la bonne ...

Tom Wenseleers

Comment effectuer une régression de crête non négative?

Réponses: