# Quadratic programming when the matrix is not positive definite

R package quadprog seems to be able to solve the quadratic programming problem only when the matrix $D$ is positive definite.

However, there is a case when the matrix $D$ is not positive definite. such as

\begin{eqnarray} \min(x^2 + y^2 - 6xy) \\ \text{subject to}\quad\quad x + y &\leq& 1,\\ 3x + y &\leq& 1.5,\\ x,y &\geq& 0. \end{eqnarray}

How can I solve this kind of problem?

• The problem might relate to the fact that if the quadratic is not positive definite it doesn't have a local minimum. In this case there should still be a global minimum, since the region is bounded. Aug 22 '13 at 22:48
• If $D$ is not PSD, then the problem is not convex. Any gradient descent algorithm will land you on a local minimum, more or less depending on its starting point. You might have to come up with an heuristic to decide when to stop the search. Aug 22 '13 at 22:54
• It's not all that hard to decide which segment of the boundary the minimum lies on. Then given that constraint, it's easy enough to cast it as a problem which does have a local minimum ... but @user603's suggestion of using a standard minimization algorithm like gradient descent can be quite useful as a general approach. Aug 22 '13 at 23:13

There are optimization routines specifically for local or global optimization of Quadratic Programming problems, whether or not the objective function is convex.

BARON is a general purpose global optimizer which can handle and take advantage of quadratic programming problems, convex or not.

CPLEX has a quadratic programming solver which can be invoked with solutiontarget = 2 to find a local optimum or = 3 to find a global optimum. In MATLAB, that can be invoked with cplexqp.

General purpose local optimizers which can handle linear constraints can also be used to find a local optimum. An example in R is https://cran.r-project.org/web/packages/trust/trust.pdf . Optimizers for R are listed at https://cran.r-project.org/web/views/Optimization.html .

In MATLAB, the function quadprog in the Optimization Toolbox can be used to find a local optimum.

In Julia, there are a variety of optimizers available.

"Any" gradient descent algorithm might not land you on anything, let alone dealing with constraints. Use a package developed by someone who knows what they are doing.

The example problem provided is easily solved to provable global optimality. Perhaps with the passage of more than 2 years it is no longer needed, or maybe being an example it never was, but in any event, the global optimum is at x = 0.321429, y = 0.535714

• +1. Lagrange multiplier methods to solve such problems are routinely taught in second-year Calculus classes. With them, one readily obtains $x=9/28$ and $y=15/28$ (which is attained along the boundary $3x+y=3/2$).
– whuber
Nov 15 '15 at 14:44

You can build a workaround by using nearPD from the Matrix package like so:
nearPD(D)\$mat.

nearPD computes the nearest positive definite matrix.

• +1 because it is a relatively straightforward approximate solution. (i don't recall seeing this question otherwise I would have given it myself in a comment.) Having said that, one essentially sets the negative eigenvalues to zero when using this technique and then reconstruct the original matrix; if the corresponding modes of variation are significant this approximation can be seriously flawed. Nov 15 '15 at 10:20
• Agree with last sentence in preceding comment. This is an excellent technique to use as long as you don't care in the slightest whether your answer is correct, or even in the right ballpark, city, or state. If your objective "Hessian" matrix is within "tolerance" away from being positive definite, this approach could actually be reasonable, otherwise, not. Nov 15 '15 at 12:49