# 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 ?

 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 '14 at 21:43
• What is the dimension of C and d of Brian Borchers's answer? – user18168 Nov 9 '15 at 0:48
• If $A$ is of size $m$ by $n$ and $L=I$, then $C$ is of size $m+n$ by $n$ and $d$ is a vector of length $m+n$. Note that if $A$ is sparse then $C$ will also be sparse. – Brian Borchers Mar 8 at 0:37

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 '14 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 '14 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 '14 at 2:37

If the fortran code https://github.com/scipy/scipy/blob/master/scipy/optimize/nnls/nnls.f doesn't have a Tikhonov regularization option, then it would probably take some work to add this feature to scipy.

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).

• That was my best idea. However, my FORTRAN coding is too weak to do it. Maybe someone can code it and contribute with the Scipy community. – user3259573 Feb 1 '14 at 16:09
• I really only want to add any regularization to the NNLS. It was a point that a reviewer on my paper brought up. So I believe I would have to stick with python and NNLS. There are implementations in matlab, but now is too late to change (the paper is almost accepted). – user3259573 Feb 1 '14 at 16:24