# How to perform non-negative ridge regression?

How to perform non-negative ridge regression? Non-negative lasso is available in scikit-learn, but for ridge, I cannot enforce non-negativity of betas, and indeed, I am getting negative coefficients. Does anyone know why this is?

Also, can I implement ridge in terms of regular least squares? Moved this to another question: Can I implement ridge regression in terms of OLS regression?

• There are two quite orthogonal questions here, I'd consider breaking out the "can I implement ridge in terms of least squares" as a separate question. – Matthew Drury Mar 25 '16 at 16:08

The rather anti-climatic answer to "Does anyone know why this is?" is that simply nobody cares enough to implement a non-negative ridge regression routine. One of the main reasons is that people have already started implementing non-negative elastic net routines (for example here and here). Elastic net includes ridge regression as a special case (one essentially set the LASSO part to have a zero weighting). These works are relatively new so they have not yet been incorporated in scikit-learn or a similar general use package. You might want to inquire the authors of these papers for code.

EDIT:

As @amoeba and I discussed on the comments the actual implementation of this is relative simple. Say one has the following regression problem to:

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

where $x_1$ and $x_2$ are both standard normals such as: $x_p \sim N(0,1)$. Notice I use standardised predictor variables so I do not have to normalise afterwards. For simplicity I do not include an intercept either. We can immediately solve this regression problem using standard linear regression. So in R it should be something like this:

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


Notice the last line. Almost all linear regression routine use the QR decomposition to estimate $\beta$. We would like to use the same for our ridge regression problem. At this point read this post by @whuber; we will be implementing exactly this procedure. In short, we will be augmenting our original design matrix $X$ with a $\sqrt{\lambda}I_p$ diagonal matrix and our response vector $y$ with $p$ zeros. In that way we will be able to re-express the original ridge regression problem $(X^TX + \lambda I)^{-1} X^Ty$ as $(\bar{X}^T\bar{X})^{-1} \bar{X}^T\bar{y}$ where the $\bar{}$ symbolises the augmented version. Check slides 18-19 from these notes too for completeness, I found them quite straightforward. So in R we would some like the following:

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


and it works. OK, so we got the ridge regression part. We could solve in another way though, we could formulate it as an optimisation problem where the residual sum of squares is the cost function and then optimise against it, ie. $\displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$. Sure enough we can do that:

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  which as expected again works. So now we just want :$ \displaystyle \min_{\beta} || \bar{y} - \bar{X}\beta||_2^2$where$\beta \geq 0$. Which is simply the same optimisation problem but constrained so that the solution are non-negative. 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  which shows that the original non-negative ridge regression task can be solved by reformulating as a simple constrained optimisation problem. Some caveats: 1. I used (practically) normalised predictor variables. You will need to account of the normalisation yourself. 2. Same thing goes for the non normalisation of the intercept. 3. I used optim's L-BFGS-B argument. It is the most vanilla R solver that accepts bounds. I am sure that you will find dozens of better solvers. 4. In general constraint linear least-squares problems are posed as quadratic optimisation tasks. This is an overkill for this post but keep in mind that you can get better speed if needed. 5. As mentioned in the comments you could skip the ridge-regression as augmented-linear-regression part and directly encode the ridge cost function as an optimisation problem. This would be a lot simpler and this post significantly smaller. For the sake of argument I append this second solution too. 6. I am not fully conversational in Python but essentially you can replicate this work by using NumPy's linalg.solve and SciPy's optimize functions. 7. To pick the hyperparameter$\lambda$etc. you just do the usual CV-step you would do in any case; nothing changes. Code for 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  • This is somewhat misleading. Non-negative ridge regression is trivial to implement: one can rewrite ridge regression as usual regression on extended data (see comments to stats.stackexchange.com/questions/203687) and then use non-negative regression routines. – amoeba Mar 25 '16 at 18:44 • I agree it is simple to implement (+1 to that). (I upvoted earlier yours and Glen's comment on the other thread too). The question is why is not implemented though, not if it is difficult. On that matter, I strongly suspect that directly formulating this NNRR task a optimisation problem is even simpler that first formulating it as an extending data regression and then using Quad. Prog. optimisation to solve this regression. I did not say this in my answer because it would venture in the implementation part. – usεr11852 Mar 25 '16 at 19:20 • Or just write it in stan. – Sycorax Mar 25 '16 at 19:24 • Ah, okay; I understood the Q as mainly asking how to do non-negative ridge (and only asking about why it's not implemented in passing); I even edited to put this into the title. In any case, how to do it seems to me to be a more interesting question. If you can update your answer with explanations on how to implement non-negative ridge, I think it will be very useful for future readers (and I will be happy to upvote :). – amoeba Mar 25 '16 at 20:01 • Cool, I will do it later (I didn't notice the new title, sorry about that). I will probably give the implementation in terms of OLS/pseudo-observations so we answer the other question too. – usεr11852 Mar 25 '16 at 20:17 R package glmnet that implements elastic net and therefore lasso and ridge allows this. With parameters lower.limits and upper.limits, you can set a minimum or a maximum value for each weight separately, so if you set lower limits to 0, it will perform nonnegative elastic-net (lasso/ridge). There is also a python wrapper https://pypi.python.org/pypi/glmnet/2.0.0 Recall we are trying to solve: $$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda \| x \|^2_2 \,\,\,\,\text{s.t. }x>0$$ is equivalent to: $$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+ \lambda x^\top I x \,\,\,\,\text{s.t. }x>0$$ with some more algebra: $$\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$$ The solution in pseudo-python is simply to do: Q = A'A + lambda*I c = - A'y x,_ = scipy.optimize.nnls(Q,c)  for a slightly more general answer. • Should the line c = - A'y not read c = A'y ? I think this is correct, though one should note that the solution is slightly different from scipy.optimize.nnls(newMatX, newVecY), where newMatX is X row augmented with a diagonal matrix with sqrt(lambda) along the diagonal and NewVecY is Y augmented with nvar zeros. I think the solution you mention is the correct one... – Tom Wenseleers Aug 29 '19 at 14:22 • this is incorrect, as I just discovered, and leads to a particularly large error for example in over-specified systems ($n \gt p$): we want to augment$A$with$\lambda I$and$y$with$0$, not simplify as above. The simplification is only correct in an unconstrained OLS settings. – Pierre D Feb 12 at 20:37 • @PierreD I don't understand what is wrong with the above. Why does it fail to you? If you provide your own answer that would help too or error messages. – Charlie Parker Feb 12 at 21:12 • It took me a while to figure out why solving via nnls() for$(X_{*}^\prime X_{*})\beta = X_{*}^\prime y_{*}$gives correct results, but incorrect for$(X^\prime X + \lambda I)\beta = X^\prime y$when the system is over-determined (more rows than columns). I think that's because the latter already does a least square projection (and no constraints). In summary: the two forms are equivalent if$n = p$or full OLS (no constraints on the coefficients). But if$n > p\$ and we want non-negative Ridge, then the two yield completely different results. In my answer I point to a full notebook. – Pierre D Feb 12 at 21:56
• @CharlieParker: I amended my own answer to show the difference. – Pierre D Feb 12 at 22:18

I know it's an oldie, but I put an example here of how to do this in Python either by using ElasticNet (approximation), or via nnls().

Basically, for ElasticNet, you can use:

from sklearn import linear_model as lm

eln = lm.ElasticNet(l1_ratio=0, fit_intercept=False)
act_alphas, coefs, dual_gaps = eln.path(X, y, alphas=alphas, positive=True)


but be forewarned that ElasticNet spews a lot of warnings for large values of alpha, as it doesn't converge.

Via nnls():

from scipy.optimize import nnls

def nnls_ridge(X, y, alpha):
"""return non-negative Ridge coefficients"""
p = X.shape[1]
Xext = np.vstack((X, alpha * np.eye(p)))
yext = np.hstack((y, np.zeros(p)))
coefs, _ = nnls(Xext, yext)
return coefs


Note

(Following the notations of this excellent post).

While it is true that in an unconstrained setting, the OLS formulation of Ridge:

$$(X_{*}^\prime X_{*})\beta = X_{*}^\prime y_{*}.$$

is equivalent to:

$$(X^\prime X + \lambda I)\beta = X^\prime y.$$

in a constrained setting (e.g. bounds = (0, np.inf), or using nnls()), the two are not equivalent, as $$(X^\prime X + \lambda I)$$ has already been reduced to a $$p \times p$$ matrix by (unconstrained) LS projection.

In particular, my initial attempt at nnls_ridge() is wrong:

def nnls_ridge_wrong(X, y, alpha):
p = X.shape[1]
Xmod = X.T @ X + alpha * np.eye(p)
ymod = X.T @ y
coefs, _ = nnls(Xmod, ymod)
return coefs


Example (using the functions described in the GitHub link at the beginning of this post):

X = np.arange(6).reshape((-1, 2))
y = np.array([8,1,1])
a = 0.001

>>> ols_ridge(X, y, a)
array([-8.56345838,  6.81837427])

>>> nn_ridge_path_via_elasticnet(X, y, [a])[1].T
array([[0.        , 0.45708041]])

>>> nnls_ridge(X, y, a)
array([0.        , 0.45714284])

>>> nnls_ridge_wrong(X, y, a)
array([0.        , 0.37663842])


For a $$X_{12,10}$$ and $$y_{12}$$ where all components are Normal $$N(0,1)$$, we see some drastic differences of $$R^2$$ for the nnls_ridge_wrong (bottom left pair of plots):