Rank Degenerate Non Negative Least Squares

Question

I'm following an algorithm in the book "Solving Least Squares Problems" by Lawson and Hanson (#15 in Siam's Classics in Applied Mathematics) for solving non-negative least squares.

$$\begin{array}{ll} \text{minimize} & \|Ex - f\|\\ \text{subject to} & x \ge 0\end{array}$$

It's basically using a pivoting scheme to try and find which variables can be set free and which have to be clamped to $0$. It forms a temporary matrix $E_p$, where for free variables, it takes the corresponding column in $E$, and for clamped variables, it takes a column of 0s, and computes the least squares solution for $E_p x \cong f$.

Only this isn't really the least squares solution, in the sense of having minimum residual and norm. It's a basic solution, where the fewest number of variables are set to be free as possible. In cases where the matrix is rank degenerate it won't "spread out" the answer as a proper least squares solution would.

For instance:

$$E = \left( \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 2 & 3 \\ 1 & 2 & 3 \\ \end{array} \right), \qquad f = \left( \begin{array}{ccc} 1 \\ 1 \\ 1 \end{array} \right)$$

produces the vector $x = \left( \begin{array}{ccc} 0 & 0 & \frac{1}{3} \end{array} \right)$ when I use my algorithm, but its unconstrained least squares solution is $x = \left( \begin{array}{ccc} \frac{1}{14} & \frac{2}{14} & \frac{3}{14} \end{array} \right)$. For comparison, Octave's lsqnonneg also returns $x = \left( \begin{array}{ccc} 0 & 0 & \frac{1}{3} \end{array} \right)$.

Both vectors, when multiplied by $E$, produce $f$ again, so they're both correct in the sense that they have minimum residual. But for my purposes the smaller norm of the second would be better.

Is there a way to "relax" this basic solution to get back a minimum norm solution? Note that the matrix isn't necessarily square. Or a different algorithm to use? I'd like to stick to pivoting algorithms over iterative gradient descent, but otherwise I'm pretty flexible.

Answer 1

You can solve the following damped system, where $\epsilon$ is a small number with respect to $\|E\|$, and $I$ is the identity matrix.

\begin{align} \text{minimize} & \quad \left\|\left( \begin{array}{c} E\\ \epsilon\, I \end{array} \right) x- \left( \begin{array}{c} f\\ 0 \end{array} \right)\right\|\nonumber\\ {\rm subject\; to} & \quad x\ge0 \end{align}

As an example, for your proposed problem, using Lawson and Hanson NNLS program (or MATLAB lsqnonneg) to solve the problem with $\epsilon=10^{-4}$ gives the desired minimum-norm solution to 11 digit accuracy.

Answer 2

You can try some post-processing after calculating of $x$: compute the null space of $E$, then solve $$ \min \|x + v \| $$ subject to $Ev=0$ and $x+v\ge 0$.