# How does one do sparse non-negative least squares using $K$ regularizers of the form $x^\top R_k x$?

I want to solve:

$$J_{R_K,L1}(x) = ||Ax - y ||^2 + \sum^K_{k} \lambda_k x^\top R_kx + \alpha \| x \|_1, x>0$$

of course, in the case where $R_1 = I$ we get non-negative Elastic Net regularization:

$$J_{L2,L1}(x) = ||Ax - y ||^2 + \beta ||x||^2_2 + \alpha \| x \|_1, x>0$$

without the L1 norm its trivial to solve the problem:

$$J_{R_K,L1}(x) = ||Ax - y ||^2 + \sum^K_{k} \lambda_k x^\top R_k x, x>0$$

as long as one one can write it in a quadratic form $\frac{1}{2} x^\top Q x + c^\top x$ and use for non-negative least squares. This can be done by using the fact $||Ax - y ||^2 = x^\top A^T A x -2 y^\top A x + y^\top y$ and factoring $x^\top$ and $x$ from the left and right respectively:

$$\frac{1}{2} x^\top Q x + c^\top x = x^{\top} \left( A^\top A + \sum^K_{k} \lambda_k R_k \right)x + (-A^\top y)^\top x$$

and then plug in to any standard non-negative least squares software. This is really simple because one just needs to change the $A^\top A$ design matrix.

However, when one includes the L1, since (I assume) L1 cannot be written in a quadratic form $x^\top L_1 x$ then one requires different methods (like sub-gradient methods or proximity operators). Furthermore, the problem gets more complicated because we need to include the non-negativity constraint which I guess one can just include $max(x,0)$ in the objective function:

$$\frac{1}{2} x^\top Q x + c^\top x = x^{\top} \left( A^\top A + \sum^K_{k} \lambda_k R_k \right)x + (-A^\top y)^\top x + \sum^D_{d} \max(0,x_d)$$

and use sub-gradients or proximity operators.

So my question is:

1. How do we solve this problem mathematically? Do we just have to sort of re-do the mathematics that this paper offers? Or just do sub-gradients or proximity operators? Or is there something simpler? If we do solve it any of these ways do have to implement everything from scratch or is there optimized code we may re-use?
2. More importantly, my intuition tells me that there must be a way to change the design matrix and use non-negative elastic nets solvers like the numpy/scipy one that already exist. I want to do that mainly because I want to re-use optimized code so that the methods run fast.

I am sure that with enough patience I could do step 1 (or do sub-gradient methods or proximity operators). However, is it possible for me to avoid these complicated maths AND implementation and re-use optimized code for non-negative Elastic Net that already exists?

Recall we are trying to solve:

$$\text{minimize}_{x}\,\,\,\,\left\Vert Ax-y\right\Vert _{2}^{2}+\sum_{k}\lambda_{k}x^{T}R_{k}x+\alpha\left\Vert x\right\Vert _{1}\,\,\,\,\text{s.t. }x>0$$

This problem is a perfect example “not seeing the forest for the trees.” The solution is quite simple if one noticed that we have:

$x>0$, thus, $\left\Vert x\right\Vert _{1}=\sum_{i}x_{i}=1^{T}x$ This little observation turns the original problem into a classic non-negative quadratic program:

$$\text{minimize}_{x}\,\,\,\,x^{T}\left(A^{T}A+\sum_{k}\lambda_{k}R_{k}\right)x+\left(\alpha1-2A^{T}y\right)^{T}x\,\,\,\,\text{s.t. }x>0$$

which can be solved via a host of algorithms including projected gradient, active-set methods, etc.

so the solution in pseudo-python is simply to do:

Q = A'A + sum( lambda[k], R[k] )
c = alpha 1 - 2 A'y
x,_ = python_maths.QR_non_negative(Q,c)


Answer credit to: professor Reza Borhani on Quora