# Tikhonov regularization in the non-negative least square - NNLS (python:scipy)

I am working on a project that I need to add a regularization into the NNLS algorithm. Is there a way to add the Tikhonov regularization into the NNLS implementation of scipy [1]?

[2] talks about it, but does not show any implementation. Sklearn has an implementation, but it is not applied to nnls.

-
Could you please comment on the size of your problem? I need to apply NNLS to a problem with 10^5 unknowns, but I do not know if active set methods are feasible for such a large problem. –  Alexander Feb 6 at 21:43

If what you want is to solve

$\min \| Ax - b \|_{2}^{2} + \lambda^{2} \| x \|_{2}^{2}$

subject to

$x \geq 0$,

then this is easily implemented. Construct a matrix

$C=\left[ \begin{array}{c} A \\ \lambda I \end{array} \right]$

and a vector

$d=\left[ \begin{array}{c} b \\ 0 \end{array} \right]$.

Then use your nonnegative least squares solver on

$\min \| Cx - d \|_{2}^{2}$

subject to

$x \geq 0$.

You can easily extend this to problems of the form

$\min \| Ax-b \|_{2}^{2} + \lambda^{2} \| L (x-x_{0}) \|_{2}^{2}$

by letting

$C=\left[ \begin{array}{c} A \\ \lambda L \end{array} \right]$

and

$d=\left[ \begin{array}{c} b \\ Lx_{0} \end{array} \right]$.

-
Nice +1 I would upvote but I'm too lazy to register. –  k20 Feb 1 at 16:45
Thanks Brian. I implemented that, and it runs. I have one more question: How do I choose the value of lambda? –  user3259573 Feb 4 at 1:13
There are lots of methods for selecting the regularization parameter. If you know the noise level in $b$, then you can use it as a basis for selecting $\lambda$ (pick the largest lambda that still results in statistically adequate fit to the data.) A simple heuristic that is commonly used in practice is the L-curve criterion- plot $\| Ax - -b \|$ vs. $\| x \|$, and look for a value of $\lambda$ that gives a "corner" solution that is pareto optimal. In practice, the choice of $\lambda$ is often simply subjective- what makes the solution look good. –  Brian Borchers Feb 4 at 2:37
Edit: I looked at the .pdf. I guess you want to find the smallest $\lambda$ so that $\text{argmin}_x ||Ax-b||^2 + \lambda^2 ||x||^2$ is positive. I think the $\text{argmin}_x$ for any given $\lambda$ has a more or less closed form solution, and the .pdf mentions bisection search on $\lambda$ until you find the smallest one so that the best $x$ is positive. This seems pretty straightforward, in the sense that it shouldn't need any weird decompositions or messy fortran code.
Edit 2: I noticed that this question is tagged convex-optimization. If you want to take that approach literally, then you could try cvx (http://cvxr.com/cvx/) or cvxpy (https://github.com/cvxgrp/cvxpy).