# Linear Least Squares with Non Negativity Constraint

I am interested in the linear least squares problem: $$\min_x \|Ax-b\|^2$$

Without constraint, the problem can be directly solved. With an additional linear equality constraint, the problem can be directly solved too, thanks to a Lagrange multiplier.

However, I could not find a direct way to solve this problem with a linear inequality constraint. The problem belongs to quadratic programming, and the methods mentioned on Wikipedia involve an iterative procedure. I also know about Karush–Kuhn–Tucker conditions, but I cannot deal with them in this particular problem since the primal and dual feasability conditions, and the complementary slackness conditions, are numerous and not helpful in an abstract setting.

Let us assume the linear inequality constraint is indeed a simple enforcement of non-negativity: $$x\geq0$$

Is there a direct method which could be directly applied to this simpler case? The only hope I could find so far lies in the method explained in a technical report by Gene H. Golub and Michael A. Saunders, released in May 1969, called Linear Least Squares and Quadratic Programming, and which was linked in another question.

• Lagrange multipliers account for inequality constraints as well. – Jacob Oct 25 '12 at 21:54
• Then we talk about Karush–Kuhn–Tucker conditions. – Wok Oct 25 '12 at 21:57

You are asking about the nonnegative least squares problem. Matlab has a function to solve this type of problem. Here is a paper about an algorithm to solve nonnegative least squares problems:

Arthur Szlam, Zhaohui Guo and Stanley Osher, A Split Bregman Method for Non-Negative Sparsity Penalized Least Squares with Applications to Hyperspectral Demixing, February 2010

I also think the proximal gradient method can be used to solve it pretty efficiently.

Here are some details about how to solve this problem using the proximal gradient method. The proximal gradient method solves optimization problems of the form $$\text{minimize} \quad f(x) + g(x)$$ where $f:\mathbb R^n \to \mathbb R$ is convex and differentiable (with a Lipschitz continuous gradient) and $g:\mathbb R^n \to \mathbb R \cup \{\infty\}$ is convex. (We also require that $g$ is lower semi-continuous, which is a mild assumption that is usually satisfied in practice.) The optimization variable is $x \in \mathbb R^n$. The nonnegative least squares problem has this this form where $$f(x) = \frac12\| Ax - b \|_2^2$$ and $g$ is the convex indicator function for the nonnegative orthant: $$g(x) = \begin{cases} 0 & \quad \text{if } x \geq 0, \\ \infty & \quad \text{otherwise.} \end{cases}$$ The proximal gradient iteration is $$x^{k+1} = \text{prox}_{tg}(x^k-t \nabla f(x^k)).$$ Here $t > 0$ is a fixed step size which is chosen in advance, and which should satisfy $t \leq 1/L$ (where $L$ is a Lipschitz constant for $\nabla f$). The prox-operator of $g$ is defined by $$\text{prox}_{tg}(x) = \text{arg}\,\min_u \, g(u) + \frac{1}{2t} \| u - x \|_2^2.$$ For our specific choice of $g$, you can see that $\text{prox}_{tg}(x) = \max(x,0)$, the projection of $x$ onto the nonnegative orthant. The gradient of $f$ is given by $$\nabla f(x) = A^T (Ax - b).$$

So, it's as simple as that. There is also a line search version of the proximal gradient method which I recommend using. The line search version adjusts the step size $t$ adaptively at each iteration, and it often converges much faster than the fixed step size version.

• Thanks for the answer. I will have a look at the article. However, keep in mind I would be interested in a direct formula rather than an iterative procedure. If this is not possible, then be it, I will try the iterative procedure. – Wok Oct 25 '12 at 21:42
• I'm pretty sure the existence of this Matlab function (and the iterative algorithm given in the paper, although it solves a slightly harder problem) mean that no direct formula exists. Pretty sure you have to use an iterative procedure. – littleO Oct 25 '12 at 22:13
• @littleO, How would you solve it using PGM? Something like Projected Gradient Descent? – Royi Aug 21 '17 at 19:54
• @Royi I updated my answer to provide details about that. – littleO Aug 21 '17 at 20:13
• @littleO, As I wrote, this is exactly Projected Sub Gradient Descent as I wrote here - math.stackexchange.com/a/2383728/33. By the way, the best way to accelerate it is using Nesterov / FISTA approach. – Royi Aug 21 '17 at 20:29

Solving regular Least Squares problem using Gradient Descent is easy.
How will it be different with the constraint applied?
Well, you will have to project the solution of each iteration onto the set of the constraints.

Here is the code:

%% Solution by Projected Gradient Descent

vX = zeros([numCols, 1]);

for ii = 1:numIterations
vX = vX - ((stepSize / sqrt(ii)) * mA.' * (mA * vX - vB));
vX = max(vX, 0);
end


The full code and validation against CVX is in my Mathematics Q561696 GitHub Repository.