A separable nonnegative quadratic program

Question

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).

Answer 1

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.

Answer 2

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.