4
$\begingroup$

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?

Improve this question
$\endgroup$
11
  • $\begingroup$ Does adding a small multiple of the identity helps? $\endgroup$
    – nicoguaro ♦
    Commented Jun 6, 2017 at 13:01
  • $\begingroup$ 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? $\endgroup$
    – Bryson of Heraclea
    Commented Jun 6, 2017 at 14:17
  • $\begingroup$ 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. $\endgroup$
    – nicoguaro ♦
    Commented Jun 6, 2017 at 14:23
  • 2
    $\begingroup$ 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 $\endgroup$
    – Johan Löfberg
    Commented Jun 8, 2017 at 15:49
  • 2
    $\begingroup$ 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) $\endgroup$
    – Johan Löfberg
    Commented Jun 8, 2017 at 15:55

2 Answers 2

Reset to default
5
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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... $\endgroup$
    – Bryson of Heraclea
    Commented Jun 9, 2017 at 7:41
  • 1
    $\begingroup$ 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 $\endgroup$
    – Johan Löfberg
    Commented Jun 9, 2017 at 9:03
0
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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$. $\endgroup$
    – Bryson of Heraclea
    Commented Jun 8, 2017 at 8:15
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bryson of Heraclea
    Commented Jun 8, 2017 at 8:20
  • $\begingroup$ @BrysonofHeraclea You are correct. I have edited my answer. $\endgroup$
    – Rodrigo de Azevedo
    Commented Jul 5, 2017 at 13:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.