Ridge penalized GLMs using row augmentation?

Question

I've read that ridge regression could be achieved by simply adding rows of data to the original data matrix, where each row is constructed using 0 for the dependent variables and the square root of $k$ or zero for the independent variables. One extra row is then added for each independent variable.

I was wondering whether it is possible to derive a proof for all cases, including for logistic regression or other GLMs.

Answer 1

Ridge regression minimizes $\sum_{i=1}^n (y_i-x_i^T\beta)^2+\lambda\sum_{j=1}^p\beta_j^2$.

(Often a constant is required, but not shrunken. In that case it is included in the $\beta$ and predictors -- but if you don't want to shrink it, you don't have a corresponding row for the pseudo observation. Or if you do want to shrink it, you do have a row for it. I'll write it as if it's not counted in the $p$, and not shrunken, as it's the more complicated case. The other case is a trivial change from this.)

We can write the second term as $p$ pseudo-observations if we can write each "y" and each of the corresponding $(p+1)$-vectors "x" such that

$(y_{n+j}-x_{n+j}^T\beta)^2=\lambda\beta_j^2\,,\quad j=1,\ldots,p$

But by inspection, simply let $y_{n+j}=0$, let $x_{n+j,j}=\sqrt{\lambda}$ and let all other $x_{n+j,k}=0$ (including $x_{n+j,0}=0$ typically).

Then

$(y_{n+j}-[x_{n+j,0}\beta_0+x_{n+j,1}\beta_1+x_{n+j,2}\beta_2+...+x_{n+j,p}\beta_p])^2=\lambda\beta_j^2$.

This works for linear regression. It doesn't work for logistic regression, because ordinary logistic regression doesn't minimize a sum of squared residuals.

[Ridge regression isn't the only thing that can be done via such pseudo-observation tricks -- they come up in a number of other contexts]

Answer 2

Generalizing this recipe to GLMs indeed is not difficult as GLMs are usually fit using iteratively reweighted least squares. Hence, within each iteration one can subsitute the regular weighted least squares step with a ridge penalized weighted least squares step to get a ridge penalized GLM. In fact, in combination with adaptive ridge penalties this recipe is used to fit L0 penalized GLMs (aka best subset, ie GLMs where the total number of nonzero coefficients are penalized). This has been implemented for example in the l0ara package, see this paper and this one for details.

It's also worth noting that the fastest closed-form way of solving a regular ridge regression is using

lmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  solve(crossprod(X) + diag(lambdas), crossprod(X, y))[, 1]
}

for the case where n>=p, or using

lmridge_solve_largep = function (X, Y, lambda) (t(X) %*% solve(tcrossprod(X)+lambda*diag(nrow(X)), Y))[,1]

when p>n and for a model without intercept.

This is faster than using the row augmentation recipe, i.e. doing

lmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  qr.solve(rbind(X, diag(sqrt(lambdas))), c(y, rep(0, ncol(X))))
}

If you would happen to need nonnegativity constraints on your fitted coefficients then you can just do

library(nnls)

nnlmridge_solve = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x
}

which then gives a slightly more accurate result btw than

nnlmridge_rbind = function (X, y, lambda, intercept = TRUE) {
  if (intercept) {
    lambdas = c(0, rep(lambda, ncol(X)))
    X = cbind(1, X)
  } else { lambdas = rep(lambda, ncol(X)) }
  nnls(A=rbind(X,diag(sqrt(lambdas))), b=c(Y,rep(0,ncol(X))))$x 
}

(and strictly speaking only the solution nnls(A=crossprod(X)+diag(lambdas), b=crossprod(X,Y))$x is then the correct one).

I haven't yet figured out how the nonnegativity constrained case could be further optimized for the p > n case - let me know if anyone would happen to know how to do this... [lmridge_nnls_largep = function (X, Y, lambda) t(X) %*% nnls(A=tcrossprod(X)+lambda*diag(nrow(X)), b=Y)$x doesn't work]