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.

[1] http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.nnls.html

[2] http://icses2012.pwr.wroc.pl/article/34.pdf

  • $\begingroup$ 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. $\endgroup$
    – Alexander
    Feb 6 '14 at 21:43
  • $\begingroup$ What is the dimension of C and d of Brian Borchers's answer? $\endgroup$
    – user18168
    Nov 9 '15 at 0:48
  • $\begingroup$ 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. $\endgroup$ Mar 8 '19 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]$


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

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

  • $\begingroup$ 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. $\endgroup$ Feb 1 '14 at 16:09
  • $\begingroup$ 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). $\endgroup$ Feb 1 '14 at 16:24

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.