# 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 '19 at 0:37

## 2 Answers

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 \\ \lambda 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
• Excellent answer. I do think there is one small mistake: in case you include a prior you also have to multiply that by $\lambda$, so the last entry in $d$ should be $\lambda L x_0$. – tBuLi Jun 16 '20 at 15:12
• @tBuLi fixed. Thanks for spotting that. – Brian Borchers Jun 16 '20 at 15:20

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