we know there is a fused version of LASSO. Fused LASSO adds a further regularizer demanding the smoothness of \beta. More details could be found here

I am wondering why I cannot find something similar for Ridge regression. We can also add a 'smoothness' term for Ridge regression:

$\sum_{k}||y-X\beta||^2 + \lambda_{1} * ||\beta||^2 + \lambda_{2} * ||D*\beta||^2$

where $D$ is the matrix that represents the space-relationship between individual $\beta$'s

Why I cannot find any research work/R-libraries doing this kind of regression? Is there any flaw of this formula?

  • 3
    $\begingroup$ No flaw. Indeed, I wrote about this possibility at the end of an answer at stats.stackexchange.com/a/164546/919. That answer shows that you don't need any new software to do these calculations: you can carry them out with your ordinary least squares routine. Note, too, that for given $\lambda_i$ your formula can be expressed as a standard Ridge regression problem simply by a linear change of variables in the model matrix $X$. $\endgroup$ – whuber Feb 7 '18 at 22:06
  • 2
    $\begingroup$ I recall reading a paper that used this kind of approach. They had a linear ordering of the coefficients, and encouraged them to be smooth by penalizing $\sum_i (\beta_i - \beta_{i-1})^2$ (which corresponds to a particular choice of your $D$ matrix). They solved for the coefficients by augmenting the data matrix as @whuber suggests. Unfortunately, this was a while ago so I can't remember the paper, but it's out there somewhere... $\endgroup$ – user20160 Feb 7 '18 at 22:37

Yes, there definitely is such a thing in the $\ell_2$ (ridge) case, typically used for some form of smoothness. First, note that like ordinary ridge, we have a closed form solution:

You can combine the two penalty terms, by noting that the first ($\lambda_1 \|\beta\|_2^2$) can be written as $\lambda_1 \beta^T I \beta$ where $I$ is the $p\times p$-identity matrix, while the second ($\lambda_2 \|D\beta\|_2^2$) can be written as $\lambda_2\beta^T D^TD\beta$. If we put these together, we have a combined penalty of the form $\beta^T(\lambda_1 I + \lambda_2 D^TD)\beta$.

This suggests that we should consider penalizations of the form $\beta^T\Omega \beta$ for some strictly positive-definite matrix $\beta$ as our general case. (We can sometimes get away with non-negative definite, but let's not worry about that).

In matrix form, the problem is then

$$ \text{arg min}_{\beta} \|y - X\beta\|_2^2 + \beta^T\Omega \beta$$

Taking the gradient with respect to $\beta$ gives the stationarity conditions:

$$ -X^T(y - X\beta) + \Omega \beta = 0 $$

(I dropped the factor of two that would appear in both terms). We can solve this for $\beta$ directly to get

$$ \hat{\beta}_{\text{GenRidge}} = (X^TX + \Omega)^{-1}X^Ty$$

If we set $\Omega = \lambda I$, we recover "standard" ridge regression.

As @whuber notes in another question (https://stats.stackexchange.com/a/164546/61353), if we don't want to solve this using the analytic solution, we can use a normal OLS routine if we augment our data set with some "pseudo-observations."

(There are some interesting connections to posterior means in exponential families here if you like that sort of stuff; very roughly, we can think of the penalty as being equivalent to some prior and since everything is normal, the MAP and posterior mean coincide).

In particular, instead of adding zero-response observations at $\begin{pmatrix} \sqrt{\lambda} I \end{pmatrix}$, we can add zero-response observations at $L$ where $L$ is the Cholesky factor of $\Omega = L^TL$.

Then, $X^* = \begin{pmatrix} X \\ L \end{pmatrix}$ and $y^* = \begin{pmatrix} y \\ 0 \end{pmatrix}$ so OLS on $(X^*, y^*)$ gives

$$ \hat{\beta} = ((X^*)^TX^*)^{-1}(X^*)^Ty^* = \left(\begin{pmatrix} X^T & L^T\end{pmatrix}\begin{pmatrix} X\\ L \end{pmatrix}\right)^{-1}\begin{pmatrix} X^T & L^T\end{pmatrix}\begin{pmatrix} y \\ 0 \end{pmatrix} = \left(X^TX + L^TL\right)^{-1}X^Ty = (X^TX + \Omega)^{-1}X^Ty$$

as desired.

So now we can solve the problem for general $\Omega$ -- when would we actually use such a thing?

Like the related $\ell_1$-construction (the "generalized lasso"), the most commonly used case is when $\Omega = D^TD$ where $D$ is the pairwise difference operator and $X = I$. This gives a denoising procedure very similar to total variation denoising.

A similar construction is also used for fitting spline models where we want to force the spline term to be smooth in the sense of having small second derivative everywhere. If we take $\Omega$ as a 3-banded matrix with elements $(-1, 2, -1)$, we are penalizing (a finite difference approximation to) the second derivative, yielding a "smooth" fit.

This idea can be extended further to arbitrary senses of "smoothness," e.g., spatial or temporal. If we keep the $X = I$ assumption, we're essentially denoising on an arbitrary grid (the temporal case recovers a grid which is just a chain). If we allow general $X$, the applications aren't quite as obvious, but we could use a model of this form to do, e.g., a time-varying coefficient model.

Consider a hypothetical study where we have a set of pregnant women who were exposed (one time) to a potentially beneficial, but rare, environmental factor during their pregnancy. For each subject, we have a measure of the degree of the exposure, various standard measures of maternal health, and some measure of the baby's health at birth.

The simplest thing to do would be to set up a standard GLM with maternal health and the level of exposure as terms linear predictor, but this ignores the fact that the mothers were exposed at different points in their gestation. If we believe that time-of-exposure is important, we can use a smoothly varying coefficient model where the effect of exposure at different times is assumed to vary smoothly. In this case, our model looks something like:

$$\text{argmin}_{\beta} \left\|y - \begin{pmatrix} X_{\text{health}} & X_{\text{exposure}}\end{pmatrix} \begin{pmatrix} \beta_{\text{health}} \\ \beta_{\text{exposure}}\end{pmatrix}\right\|_2^2 + \lambda \beta_{\text{exposure}} \Omega \beta_{\text{exposure}}$$

where $\Omega$ is the smoothing matrix described above, and $X_{\text{exposure}}$ is a $n \times T$-matrix which gives subject $n$'s exposure in column $t$ if she was exposed during time unit $t$ and is 0 otherwise.

The penalty term here will force the elements of $\beta_{\text{exposure}}$ to be smooth with respect to time.

(This is not a good study design, or a particularly good model for this type of data, but hopefully you get the idea.)

For a spatial example, consider building an image classifier where the pixels of the images are the columns of $X$. We would expect pixels which are nearby in the raw image to have similar information and so we could use a smoothing penalty to enforce spatial smoothness. (Here we would use a 2D version of the finite difference approximation as $\Omega$)

I can't think of examples where $\Omega$ does not encode some sort of smoothness, but I'm sure there are some in the literature.

Finally, note that you can do both sparsity and smoothness in the same model: see, e.g., the neuroscience example in Section 6 of https://arxiv.org/pdf/1309.2895.pdf.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.