# Generalized Projection of a Matrix onto the Non Negative Orthant

I'd like to solve the following problem

$$\min_{Y \geq 0} \quad \operatorname{Tr} \left( (Y-Z)H(Y-Z)^T \right)$$

where $Z$ and $H$ are given matrices. $Y$ and $Z$ are $n \times r$ matrices, where $n$ may be large but $r$ is very small, and $H$ is $r \times r$ positive definite, but not sparse or structured. Assume I can do all the works on $H$ (Cholesky, eigendecomposition, etc)

This problem could be interpreted as a generalized projection, where the distance is somewhat warped by $H$ (but here $H$ is not the Hessian of some function, since it is tiny).

Is there any way of directly solving this problem (no iterative algorithms)? If not, is there a super fast iterative method I could use? I suppose the eigenvectors of $H$ should be helpful somehow.

An iterative method (also proposed below) would be to do a projected gradient algorithm, since the gradient is not that expensive and the projection is super cheap. But, because this is one subproblem in another iterative algorithm, I'd like to avoid anything with, say, more than $r$ steps. (Conj. gradient for example takes at most $r$ steps to solve an $r\times r$ linear system.)

Thanks for any tips!

• Why not solve the unconstrained quadratic program without the nonnegativity constraints and then check if the solution is admissible? – Rodrigo de Azevedo Mar 16 '18 at 22:56
• I agree that's a good first step, but in general the solution will not be admissible, for general $H$ and $Z$. (I need this in an iterative algorithm, where $H$ keeps changing.) – Y. S. Mar 16 '18 at 23:26
• I think this can be vectorized into operations on vectors, isn't it? – Royi Mar 16 '18 at 23:31
• I guess your $H$ is positive definite, isn't it? – Surb Mar 16 '18 at 23:44
• Yes $H$ is PD. Actually it's $Y^TY+I$ and I know what $Y$ is. – Y. S. Mar 17 '18 at 0:46

We have the following quadratic program in $\mathrm X \in \mathbb R^{n \times r}$

$$\min_{\mathrm X \geq \mathrm O} \quad \mbox{tr} \left( (\mathrm X - \mathrm A) \, \mathrm Q \, (\mathrm X - \mathrm A)^\top \right)$$

where $\mathrm Q \in \mathbb R^{r \times r}$ is positive definite and, thus, has a Cholesky decomposition of the form $\rm Q = R^\top R$.

$$\mbox{tr} \left( (\mathrm X - \mathrm A) \, \mathrm Q \, (\mathrm X - \mathrm A)^\top \right) = \mbox{tr} \left( (\mathrm X - \mathrm A) \,\mathrm R^\top \mathrm R \, (\mathrm X - \mathrm A)^\top \right) = \| (\mathrm X - \mathrm A) \,\mathrm R^\top \|_{\text{F}}^2$$

Vectorizing, we obtain

$$\| (\mathrm X - \mathrm A) \,\mathrm R^\top \|_{\text{F}}^2 = \cdots = \| \left(\mathrm R \otimes \mathrm I_n\right) \mbox{vec} (\mathrm X - \mathrm A) \|_2^2$$

which is a very standard convex quadratic function. Can you take it from here?

• Sorry, I don't see an obvious way to solve this over the constraint $X\geq 0$. I had tried something like replacing $U = RX$ where $H = R^TR$, but the issue is that $X\geq 0$, not $RX$. Basically, if I wanted to do $\min_{X\geq 0} \|X-A\|_F^2$ I would just take $X_{ij} = \max\{0,A_{ij}\}$. But with the "warped" distance function it feels less simple. – Y. S. Mar 17 '18 at 0:48
• What you have is a non-negative least squares problem in $\mbox{vec}(\rm X) \geq 0_{nr}$. Take a look at this MATLAB documentation. – Rodrigo de Azevedo Mar 17 '18 at 8:05
• Oh cool! It looks like they are still using an iterative algorotihm, but their active set / lagrange multipliers thing is kind of interesting. Thanks for the pointer! – Y. S. Mar 17 '18 at 17:10

As written in comments and done by Rodrigo, using vectorization one could reform the problem into regular Constrained Least Squares Problem with inequality constraints.
Since the equation includes a Matrix a closed form solution isn't feasible directly.
One method to do so is utilizing the Projected Gradient Method.

Yet, since the the Gradient of the original problem is simple it would be better to use it directly:

$$\frac{d}{ d Y } \operatorname{Tr} \left( \left( Y - Z \right) H \left( Y - Z \right)^{T} \right) = \left( Y - Z \right) \left( {H}^{T} + {H} \right)$$

The projection is simply clipping Negative Values.

A MATLAB Code for that is given by:

mY = zeros([numRows, numCols]);
vObjValPgd(1) = hObjFun(mY);

for ii = 2:numIterations
% vG = ((mY - mZ) * mH.') + ((mY - mZ) * mH);
vG = (mY - mZ) * (mH.' + mH);
mY = mY - (stepSize * vG);

mY = max(mY, 0); %<! Projection onto Non Negative Orthant

vObjValPgd(ii) = hObjFun(mY);
end


The result is verified against CVX:

Convergence is very fast (And can be accelerated using Accelerated Gradient Descent Methods) and operations are efficient.

The code is available (Including validation by CVX) at my StackExchange Mathematics Q2694373 GitHub Repository.

• Ok so here you are using an iterative method (projected gradient method)? Cool, that's a good option. But I was wondering if there's something faster. (e.g. if $H = I$ then it's just one step and lightning fast.) I was hoping knowing the eigenvectors and Cholesky factorization of $H$ can do some kind of preconditioning to make something super fast. – Y. S. Mar 17 '18 at 1:00
• No, I don't think so. As I wrote, since you have LS with a Matrix (Not the Identity) I don't see a way to directly solve this. I'd say that using Gradient Descent with acceleration the result will be very very fast as like a direct closed form solution. – Royi Mar 17 '18 at 1:06
• Ok, that is my suspicion as well. I want to keep the question open in case there is something surprising, but I do feel a bit pessimistic about it. Thanks for looking into it! – Y. S. Mar 17 '18 at 1:24