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. Feb 6, 2014 at 21:43
• What is the dimension of C and d of Brian Borchers's answer?
– user18168
Nov 9, 2015 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. Mar 8, 2019 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 \\ \lambda Lx_{0} \end{array} \right]$$.

• Nice +1 I would upvote but I'm too lazy to register.
– k20
Feb 1, 2014 at 16:45
• Thanks Brian. I implemented that, and it runs. I have one more question: How do I choose the value of lambda? Feb 4, 2014 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. Feb 4, 2014 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$. Jun 16, 2020 at 15:12
• @tBuLi fixed. Thanks for spotting that. Jun 16, 2020 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. Feb 1, 2014 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). Feb 1, 2014 at 16:24