NNLS convergence proof

Question

I was trying to understand why the nonnegative least squares (NNLS) algorithm of Lawson & Hanson converges to a solution of $$ \min f_{0}(x) = \| Cx - d \|^{2}, \\ \text{s.t.} \ x \geq 0. $$ Pseudocode can be found here: https://en.wikipedia.org/wiki/Non-negative_least_squares or in the references below.

Useful references to me were:

• Lawson & Hanson, "Solving Least Squares Problems"
• Palomo & Guerrero-Garcia, "Solving a sequence of sparse least squares problems"
• Bro & De Jong, "A fast non-negativity-constrained least squares algorithm" (http://xrm.phys.northwestern.edu/research/pdf_papers/1997/bro_chemometrics_1997.pdf)

All the references above just mention that the algorithm converges because after every iteration the residual norm strictly decreases and therefore we have a different partition of the indices of $x$ into sets $P$ and $R$ everytime, giving only a finite number of possibilities. But I'm having difficulty in proving that the residual norm decreases in each step.

All I see so far are the following:

1. At the end of each iteration $x$ is nonnegative
2. $x$ restricted to the indices in the set $P$ is an ordinary least squares solution for the problem $\min \| C^{(P)}z^{(P)} - d \|^{2}$
3. A lemma proofed in the book of Lawson & Hanson states that the variable index $j$ that enters $P$ at the beginning of each outer iteration will yield $s_{j} > 0$, where $s$ is the ordinary least squares solution to $\min \| [C^{(P)}, C_{*,j}] z - d \|^{2}$

I think the case where the inner loop is not executed is manageable to proof by using the uniqueness of least squares estimator for full column rank matrices:

Assume $x^{(i)}$ is the solution vector after iteration $i$ with sets $P = P^{(i)}$, $R = R^{(i)}$ and index $j \in R$ will be added to $P$ at beginning of iteration $i+1$. Let $s$ the OLS solution of $\min \| [C^{(P)}, C_{*,j}] z - d \|^{2}$. If we had $f_{0}(x^{(i)}) < f_{0}(s)$, then $\tilde{s} = \left( x^{(i)}, 0 \right)^{T} $ would yield $f_{0}(\tilde{s}) = f_{0}(x^{(i)}) < f_{0}(s)$ which contradicts minimality of $s$. Remains the case where $f_{0}(x^{(i)}) = f_{0}(s)$. Because of uniqueness of OLS minimizer we would get $s = x^{(i)}$ which is not possible since $s_{j} > 0 = x^{(i)}_{j}$ (j was in set $R^{(i)}$, meaning that $x$ at this position is 0) which leaves only $f_{0}(s) < f_{0}(x^{(i)})$ and $x^{(i+1)} = s$ which proves the claim.

But for the case where the inner loop is executed (possibly multiple times) I don't see how the proof could go.

Answer 1

Replying to this since I also got confused from the description in the textbook and it took me a bit to figure out the argument:

Suppose $x$ is the current solution (at some iteration of the algorithm) with active set $P$ (i.e., $x_i >0$ for $i \in P$ and $0$ otherwise). Here is what happens afterwards.

1. Let $w = C^\top (d-Cx)$ the current Lagrange dual multipliers. If all of the $w_i \leq 0$, then we terminate. Instead let $j$ be such that $w_j > 0$ is the largest. Then we add this $j$ to the active set $P$, i.e., we have $P' = P\cup\{j\}$.
2. Next, we solve the least squares problem $\lVert{C\tilde{x}-d}\rVert^2_2$ subject to the constraint that $\tilde{x}_i = 0$ for all $i \notin P'$. Let us call the solution $z$ (keeping the notation $x$ for the solution at the start of the iteration). This $z$ has the property (cf. point 3 in the question) that $z_j > 0$, where $j$ is the index that was selected in step 1. There are two things that could happen afterwards.
3. • Suppose all $z_i > 0$, then we go back to step 1. (where we update $x \leftarrow z$ and $P \leftarrow P'$).
   • Suppose some $z_i \leq 0$. Then we let $\alpha = \max_{i \in P': z_i \leq 0 } x_i/(x_i - z_i)$, update $x \leftarrow \alpha z + (1-\alpha)x$ and go back to step 1 (with the corresponding update to the active set $P$.

Here is what we need to prove: No matter what path we followed in step 3, we will get a feasible solution with a strictly improved residual norm. Let us note first that:

$$ \lVert{Cz-d}\rVert^2_2 < \lVert{Cx-d}\rVert^2_2$$ This follows from the argument in the question and demonstrate that if the 1st part of 3. occurs, then we improve the residual norm. It remains to handle the second case.

We note the following: $\alpha \in (0,1]$ since for all $i \in P': z_i \leq 0$ it holds that $x_i > 0$ and $z_i \leq 0$. Furthermore, it is clear by the definitions that $\alpha z + (1-\alpha)x = x + \alpha(z-c)$ has all elements $\geq 0$. For the residual norm decrease, we use the convexity of the norm (triangle inequality):

$$ \begin{aligned} \lVert{C(\alpha z + (1-\alpha)x)-d}\rVert_2 &= \lVert{\alpha(Cz-d) + (1-\alpha)(Cz-d)}\rVert_2 \\ & \leq \alpha \lVert{Cz-d}\rVert_2 + (1-\alpha) \lVert{Cx-d}\rVert_2 \\ & < \lVert{Cx-d}\rVert_2. \end{aligned} $$ In the above inequality we crucially relied on the fact that $\alpha >0$. It follows that also $\alpha z + (1-\alpha)x$ has smaller residual norm than $x$.