I am repeatedly solving quadratic program,

$x^T Q x$ with time dependent linear constraints $Ax=b_t$.

Dimension of $x$ is around 10000 and there are around 50 constraints. I want to solve the system as fast as possible at each time $t$. I know that $b_t$ and $b_{t+1}$ are pretty close, say, in L2-norm. I also know, that the optimal solution $x^*|b_t$ is pretty close to $x^*|b_{t+1}$.

My question is how can I best employ this knowledge to speed up the computation.

So far, I am thinking about iterative minimization of the corresponding Lagrangian with providing the previous solution as an initial guess. I can also do a lot of precomputations and evaluate expensive inversions such as $Q^{-1}A^{T}$. The dependence of the solution on $b$ also leads me to some sensitivity analysis, but I’m clueless here.

My goal is obtaining the solution in 5 s, which as I am inexperienced in this area, don’t know if it’s even possible.

| cite | improve this question | | | | |
  • $\begingroup$ You've only got linear equality constraints- no inequalities? $\endgroup$ – Brian Borchers Dec 1 '18 at 18:52
  • $\begingroup$ Are your $Q$ and $A$ matrices sparse or dense? Is $Q$ symmetric and positive definite? $\endgroup$ – Brian Borchers Dec 1 '18 at 19:22
  • $\begingroup$ How much time are you willing to spend on one-time computations whose cost can be amortized over the many right hand sides? $\endgroup$ – Brian Borchers Dec 1 '18 at 19:42

Here's an approach for dense $Q$ and $A$ that does some relatively expensive preliminary computations that make the individual solutions very fast. I doubt that any iterative approach would be faster for the case where $Q$ and $A$ are dense. If $Q$ and $A$ are very sparse, you might be better off with an iterative approach.

We'll assume that $Q$ is symmetric (any quadratic form can easily be symmetrized) and positive definite (if not, then the problem could be unbounded.) In this case, the optimal solution is unique, and the Karush-Kuhn-Tucker conditions can be written as

$ \left[ \begin{array}{cc} Q & A^{T} \\ A & 0 \end{array} \right] \left[ \begin{array}{c} x \\ \lambda \end{array} \right]= \left[ \begin{array}{c} 0 \\ b \end{array} \right] $

This can be reduced to

$ AQ^{-1}A^{T}\lambda=-b $


$ x=-Q^{-1}A^{T}\lambda. $

Before solving any of these systems, we can precompute the Cholesky factorization of $Q$,

$ Q=R^{T}R $

We then use the Cholesky factorization to get

$ B=Q^{-1}A^{T} $

and let

$ C=AQ^{-1}A^{T}=AB $

Note that you want to avoid computing $Q^{-1}$, since this is much more expensive than computing $Q^{-1}A^{T}$.

We also compute the Cholesky factorization of $C$,

$ AQ^{-1}A^{T}=C=S^{T}S. $

Now, for each right hand side $b$, we use $S$ to solve


and then use $R$ to obtain


On my desktop machine, for a random dense problem with $n=10000$ and $m=50$, the setup computations take about 6 seconds, and then each individual solve runs in less than 0.1 seconds.

| cite | improve this answer | | | | |
  • $\begingroup$ Great. This works perfectly. In my case of Q I had to add slight regularization to the diagonal, even with all its Eigen values positive, probably due to numerical errors. Although it makes sense for my problem, since it is originally a regression. Thank you. $\endgroup$ – reicja Dec 3 '18 at 19:45

Your Answer

