Quadratic programs with rank deficient positive semidefinite matrices - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2019-08-20T23:35:36Z https://scicomp.stackexchange.com/feeds/question/27066 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://scicomp.stackexchange.com/q/27066 3 Quadratic programs with rank deficient positive semidefinite matrices Bryson of Heraclea https://scicomp.stackexchange.com/users/24486 2017-06-06T09:26:38Z 2017-07-05T13:14:31Z <p>Let $A$ be a $n\times n$ square symmetric matrix. In addition, $A\succeq0$ and $\mathrm{rank}(A)&lt;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.</p> <p>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 <a href="http://sedumi.ie.lehigh.edu/" rel="nofollow noreferrer">SeDuMi</a> and with matrix $$A=\begin{bmatrix}10&amp;17&amp;25&amp;-5&amp;-9\\ 17 &amp; 29 &amp; 43 &amp; -9 &amp; -16\\ 25 &amp; 43 &amp; 65 &amp;-15 &amp;-26\\ -5 &amp; -9 &amp; -15 &amp; 5 &amp; 8\\ -9 &amp; -16 &amp; -26 &amp; 8 &amp; 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 <code>eig()</code> function gives following results:</p> <pre><code>&gt;&gt; eig(A) ans = -1.6661e-014 -4.4496e-016 4.2249e-015 5.0536 116.95 </code></pre> <p>SeDuMi (via <a href="https://yalmip.github.io/" rel="nofollow noreferrer">Yalmip</a>, excuted in Matlab 2007b) gives following error message:</p> <pre><code>Exiting: the solution is unbounded and at infinity; the constraints are not restrictive enough. </code></pre> <p>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?</p> https://scicomp.stackexchange.com/questions/27066/-/27076#27076 0 Answer by Rodrigo de Azevedo for Quadratic programs with rank deficient positive semidefinite matrices Rodrigo de Azevedo https://scicomp.stackexchange.com/users/20417 2017-06-07T13:56:29Z 2017-07-05T13:14:31Z <p>The following <a href="http://www.cvxpy.org" rel="nofollow noreferrer">CVXPY</a> script:</p> <pre><code>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 &gt;= 0 ] prob = Problem(objective, constraints) # solve optimization problem prob.solve() print "x =", x.value </code></pre> <p>produces the following vector</p> <pre><code>x = [[ 4.44444306e-01] [ 2.85606990e-11] [ 1.71645999e-10] [ 2.22221718e-01] [ 2.22222588e-01]] </code></pre> <p>which does satisfy the constraints. Perhaps the constraints are indeed "restrictive enough".</p> https://scicomp.stackexchange.com/questions/27066/-/27085#27085 4 Answer by Johan Löfberg for Quadratic programs with rank deficient positive semidefinite matrices Johan Löfberg https://scicomp.stackexchange.com/users/1662 2017-06-08T16:06:32Z 2017-06-08T16:38:13Z <p>To ensure this does not drown in the comments, I make it an answer.</p> <p>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.</p>