# A separable nonnegative quadratic program

I have spent quite some time trying to solve the following quadratic program: $$\min \sum_{i=1}^n (\frac{1}{2}x_i^TQx_i+c_i^Tx_i), \quad \mathrm{s.t. } \quad x_i\ge 0 \quad \forall i,$$ where $n$ is very large and $Q$ is a sparse positive semi-definite matrix (given at the end of this post).

Although the sub-problem for each $i$ can be easily solved (numerically), since $n$ is very large, it is infeasible to loop and solve the sub-problems separately.

I would like to know if there is any special method to deal with this type of problem.

--

The matrix $Q$ in my particular problem is given by $Q=A^TA+B^TB$, where $A$ and $B$ are defined as follows.

Denote $e=(1,1,\ldots,1)^T\in\mathbb{R}^d$, then $$A=\begin{bmatrix} e & & & \\ & e & & \\ & & \ddots & \\ & & & e \end{bmatrix}$$ and $$B=\begin{bmatrix}\mathrm{diag}(e) & \mathrm{diag}(e) & \cdots & \mathrm{diag}(e)\end{bmatrix}$$ where $e$ appears $d$ times in $A$ and $d$ times in $B$ (i.e. $A,B\in\mathbb{R}^{d\times d^2}$).

For example, if $d=3$ we have:

$$A=\begin{bmatrix} 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \end{bmatrix}$$

$$B=\begin{bmatrix} 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 & 0\\ 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1 \end{bmatrix}$$ and thus $$Q=\begin{bmatrix} M & I & I\\ I & M & I\\ I & I & M \end{bmatrix}$$ where $$M=\begin{bmatrix} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{bmatrix},\quad I = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix}.$$

The Matlab code for creating $Q$:

d = 16;
one = ones(d,1);
A = kron(eye(d),one');
B = repmat(diag(one),1,d);
Q = A'*A + B'*B;


Thank you in advance for any suggestions.

Update: according to the comments, it was not clear why I do not want to solve the subproblems separately using a loop. So here are some clarifications.

If we do a loop and solve each subproblem at each iteration, then it will take too much time if the number of subproblems is huge. For example, if solving a subproblem takes 1 second, then solving 220512 subproblems will take more than 60 hours.

We use some iterative method to solve the subproblems, then instead of updating each $x_i$ the ones after the others (using the loop), can we update them simultaneously? We can reformulate the problem as: $$\min\quad 1/2x^T\mathrm{diag}(Q,Q,...,Q)x + c^Tx, \quad \mathrm{s.t. } \quad x\ge 0,$$ where $x=(x_1,...,x_n)$ and $c=(c_1,...,c_n)$ and then use a sparse QP solver. I tried using Matlab's solver (an example is given here), however, the sparse matrix $\mathrm{diag}(Q,Q,...,Q)$ takes ~9GB of memory and thus slows down everything. Moreover, in using (naively) a sparse QP solver, the special structure of $\mathrm{diag}(Q,Q,...,Q)$ is not exploited (although it is considered to be sparse by the solver, I think we can exploit even more than that).

• Is there anything coupling the subproblems? From what you've written, I don't see why you can't solve them separately (i.e., you have $n$ small problems, not one large problem). For the small problems, look into Newton-Krylov methods for solving the subproblems, for which you only need to compute matrix-vector products (which you can do without constructing the matrices $A$ or $B$, let alone $Q$). Commented Jun 27, 2014 at 21:39
• @ChristianClason: please see the update.
– f10w
Commented Jun 27, 2014 at 23:14
• The point, which Brian also made, is that there is no mathematical reason why solving the small problems all at once as one big problem should take less time than solving them in sequence; in fact, quite the opposite. Commented Jun 28, 2014 at 8:41
• Dear @ChristianClason. Thank you very much for your interest in my question, and please forgive for this late reply. I have been involved in another problem and have to put all my effort into it. When I finish I will come back to the above question. Thank you (and sorry) again. Best regards.
– f10w
Commented Jul 17, 2014 at 22:12

With $d=16$, the $Q$ matrix is just 256 by 256. Thus the individual subproblems are quite small.

Your $Q$ matrices are singular, so the optimal solution to each of the subproblems is likely to be nonunique. Do you care which solutions you end up with? For example, are minimum norm solutions required?

You write that "since n is very large, it is infeasible to loop and solve the sub-problems separately." Solving one very large problem shouldn't be any more efficient than solving the subproblems separately. If it is, then this indicates that overhead in the solver is excessive and you might want to look for a better solver.

However, if you want to better understand the performance of the QP solver that you're using, then you might try some computational experiments. You could for example truncate the sum to $n=10000$ terms, then solve 1 subproblem at time, and compare this with solving 10 subproblems at a time and then 100 subproblems at a time to see whether you get any speedup from combining the subproblems. If this does help things, then you can always adjust the number of subproblems that you solve at a time to match the available memory.

If you've got access to a cluster of computers, then you could easily make efficient use of as many processors as you can get by breaking the sum up into parts and parceling these out to the individual processors.

• Hi Brian. Thanks for the answer. a) No I do not care which solution I will get. b) Yes I see your point but what I meant by "it is infeasible to loop and solve the sub-problems separately" is that if we do a loop and solve each subproblem at each iteration, then it will take too much time if the number of subproblems is huge. For example, if solving a subproblem takes 1 second, then solving 220512 subproblems will take more than 60 hours. c) Thus, what I meant by "not to solve the subproblems separately" is very similar to what you suggested... (continued below)
– f10w
Commented Jun 27, 2014 at 20:04
• ...to what you suggested about gathering a number of subproblems and then solve them at once. I though that, if we use some iterative method to solve the subproblems, then instead of updating each $x_i$ the ones after the others (using the loop), can we update them simultaneously? We can reformulate the problem as $$1/2x^T\mathrm{diag}(Q,Q,...,Q)x + c^Tx$$ where $x=(x_1,...,x_n)$ and $c=(c_1,...,c_n)$ and then use a QP solver. However,... (continued below)
– f10w
Commented Jun 27, 2014 at 20:12
• However, as I mentioned in a comment (to the other answer), the sparse matrix $\mathrm{diag}(Q,Q,...,Q)$ (I said $Q$ in that comment, it was a typo) takes ~9GB of memory and thus slows down everything. Moreover, in using (naively) a sparse QP solver, the special structure of $\mathrm{diag}(Q,Q,...,Q)$ is not exploited (although it is considered to be sparse by the solver, I think we can exploit even more than that).
– f10w
Commented Jun 27, 2014 at 20:17
• The main structure-exploiting tricks you could use would be 1) using an LU decomposition (with pivoting) of $Q$ to determine an LU decomposition of your block diagonal matrix, 2) using matrix-free methods plus a good preconditioner. 1) is going to be roughly equivalent to solving each subproblem separately, which is why BrianBorchers is telling you that solving the entire problem all at once won't be faster than solving subproblems separately. 2) is what ChristianClason is suggesting, which might be helpful, if you can find a good implementation out there; it is still an active research area. Commented Jun 27, 2014 at 22:01
• Thanks @GeoffOxberry. What I meant by "separately" should read "separately using a loop", which I would like to avoid. I think this is a question of implementation: the subproblems can be solved separately, but how to do that in a parallel manner?
– f10w
Commented Jun 27, 2014 at 23:37

Before venturing for specialize algorithms, have you tried a good QP solver as MOSEK, CPLEX or GUROBI?

Given the extreme sparsity of your matrices, you could be able to solve this problems.

• In my probelm: n = 220512 and d = 16. I tried to use Matlab QP solver, like this example: mathworks.fr/fr/help/optim/ug/…. However, the sparse (not full) matrix $Q$ takes ~9GB of memory and thus slows down everything (although it takes only 1.1Gb when saved to file), so I guess there is a better a way to do.
– f10w
Commented Jun 27, 2014 at 16:07