Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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

share|improve this question
    
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

2 Answers 2

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]$.

share|improve this answer
    
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

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

share|improve this answer
    
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 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 at 16:24

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.