# Quadratic programs with rank deficient positive semidefinite matrices

Let $A$ be a $n\times n$ square symmetric matrix. In addition, $A\succeq0$ and $\mathrm{rank}(A)<n$. This means that all eigenvalues are non-negative, but also that there are some zero eigenvalues. I want to perform quadratic optimization of the functional $$\frac{1}{2}x^HAx+b^Hx$$ under convex conditions $x_i\geq0$ and $y^Hx=0$, where $b,y$ some fixed vectors.

The problem arises when, due to numerical rounding errors, the matrix has some very small negative spurious eigenvalues. Then, then functional becomes unbounded. I have tried this with SeDuMi and with matrix $$A=\begin{bmatrix}10&17&25&-5&-9\\ 17 & 29 & 43 & -9 & -16\\ 25 & 43 & 65 &-15 &-26\\ -5 & -9 & -15 & 5 & 8\\ -9 & -16 & -26 & 8 & 13 \end{bmatrix},\quad b=-\begin{bmatrix}1\\1\\1\\1\\1 \end{bmatrix},\quad y=\begin{bmatrix} 1\\1\\1\\-1\\-1 \end{bmatrix}$$ If you calculate the rank with MATLAB you will get $\mathrm{rank}(A)=2$. However, the numerical calculation with the help of the eig() function gives following results:

>> eig(A)

ans =

-1.6661e-014
-4.4496e-016
4.2249e-015
5.0536
116.95


SeDuMi (via Yalmip, excuted in Matlab 2007b) gives following error message:

Exiting: the solution is unbounded and at infinity;
the constraints are not restrictive enough.


Do you have any idea, how to effectively address this problem numerically? I have tried diagonalizing $A$ and replacing the small negative eigenvalues with arbitrary small positive numbers, but this seems, well ... arbitrary, and I fear that it might produce numerical errors in the optimization. What do you think?

• Does adding a small multiple of the identity helps? – nicoguaro Jun 6 '17 at 13:01
• I have tried replacing the negative eigenvalues with a small fixed number (e.g. a small multiple of the largest positive eigenvalue), which worked. Adding a multiple of the identity would be another option. However, the question is how to choose the size of the perturbation, so that the end result is not perturbed too much. Or is there a way to compensate for the perturbation a posteriori, to recover the exact solution? – Bryson of Heraclea Jun 6 '17 at 14:17
• If you add a multiple of the identity all the eigenvalues are shifted by the scaling constant. Then, I suppose that you can recover the "original" values that way. – nicoguaro Jun 6 '17 at 14:23
• To begin with, using an SDP solver (SeDuMi) to solve a simple convex QP is a bit overkill. Anyway, are you using a recent version of SeDuMi and YALMIP. I've tried it on my machine, and YALMIP+SeDuMi easily finds a solution without any issues. YALMIP has various safeguards to detect the low-rank but psd structure here (despite eig etc saying it is indefinite). YALMIP will exploit the rank-2 structure when setting up the resulting SOCP and thus the low-rank should be a non-issue – Johan Löfberg Jun 8 '17 at 15:49
• BTW, you're not using SeDuMi here. That error messsage is from quadprog. To use SeDuMI you explicitly have to select SeDuMi as a solver, as YALMIP picks the most simple solver possible (in your case, you have a quadratic program, and you have the QP solver quadprog installed, hence quadprog is used) – Johan Löfberg Jun 8 '17 at 15:55

To ensure this does not drown in the comments, I make it an answer.

The solver used is not SeDuMi, as claimed in the question. The solver used is quadprog, and that solver (or more specifically, a severely outdated version of it), apparently had numerical issues on this particular instance. A recent version of quadprog, or any reasonably robust solver, solves this problem without issues.

• Dear Johan, thank you very much for your reply, I didn't expect the creator of Yalmip to see this post! You are right, I forgot to set the solver settings to SeDuMi. When I used it, the problem was solved alright. By the way, since SeDuMi is "overkill", which solver would you recommend? I have a 2014 version of Yalmip... – Bryson of Heraclea Jun 9 '17 at 7:41
• First, don't use a 2014 version...For quadratic programming, you have excellent commercial solvers free for academia such as mosek, gurobi, mosek, but there are alternatives as listed here yalmip.github.io/tags/#quadratic-programming-solver – Johan Löfberg Jun 9 '17 at 9:03

The following CVXPY script:

from cvxpy import *
import numpy as np

# optimization variables
x = Variable(5)

# matrix A and vectors b and c
A = np.array([[10, 17, 25, -5, -9],
[17, 29, 43, -9,-16],
[25, 43, 65,-15,-26],
[-5, -9,-15,  5,  8],
[-9,-16,-26,  8, 13]])
b = np.ones(5)
c = np.array([ 1, 1, 1,-1,-1])

# build optimization problem
objective = Minimize( 0.5 * quad_form(x, A) - b.T * x )
constraints = [ c.T * x == 0, x >= 0 ]
prob = Problem(objective, constraints)

# solve optimization problem
prob.solve()
print "x =", x.value


produces the following vector

x = [[  4.44444306e-01]
[  2.85606990e-11]
[  1.71645999e-10]
[  2.22221718e-01]
[  2.22222588e-01]]


which does satisfy the constraints. Perhaps the constraints are indeed "restrictive enough".

• You are right that the unconstrained functional is unbounded from below. If we diagonalize, transform variables and complete the squares we can show that it is the sum of squares and a linear term (sth. like $d_1(\tilde{x}_1-c_1)^2+d_2(\tilde{x}_2-c_2)^2-x_3-x_4-x_5+c_0$. However, this functional is bounded below under the constraint $x_i\geq 0$. I insist that the problem arises due to negative eigenvalues. If e.g. there exists $d_3<0$ (very small) we would have $d_1(\tilde{x}_1-c_1)^2+d_2(\tilde{x}_2-c_2)^2-|d_3|(\tilde{x}_3-c_3)^2-x_4-x_5+c_0$ which becomes unbounded even if $x_i\geq0$. – Bryson of Heraclea Jun 8 '17 at 8:15
• A further indication that negative eigenvalues are to blame is that when I slightly perturb the matrix (I replaced the negative eigenvalues with 2.6e-14), so that numerically it has only positive eigenvalues, SeDuMi is able to solve the problem and gives the same solution as your CVX Python code. Nevertheless, I want to thank you very much for your reply, since it cleared up a misunderstanding and made me realize that other numerical routines (like CVX in Python) might be more stable. – Bryson of Heraclea Jun 8 '17 at 8:20
• @BrysonofHeraclea You are correct. I have edited my answer. – Rodrigo de Azevedo Jul 5 '17 at 13:16