cvxpy returns 'infeasible_inaccurate' for a feasible simple problem

[Emmanuel Sérié]

I try to solve the following problem.
The optimization variable is $x$, $\delta$ is an auxiliary variable, the parameters $A$ and $B$ do not have a definite sign and $\Omega$ is definite positive.




Here is the code that I tried for a minimal instance of the problem ($i\in \{0\}$ and $a \in \{0,1\}$).
The code show that the problem is feasible but cvxpy say that it is 'infeasible_innacurate".

import cvxpy as cp
import numpy as np
from scipy.linalg import block_diag


def test_problem():
    # Define the problem
    na = 2
    ni = 1

    A = np.array([[1757.4000007, 8786.99999824]])
    B = np.array([[-4292.5, 12675.49999661]])

    # Prepare the quadratic form
    # array N x N
    omega = np.array([[1.0]])
    # array AN x AN
    omega = block_diag(omega, omega)

    # variable of size AN
    xx = cp.Variable(na * ni)

    delta = cp.mul_elemwise(A.flatten(), xx) + B.flatten()
    cost = cp.quad_form(delta, omega)
    obj = cp.Minimize(cost)
    assert obj.is_dcp()

    # Constraints
    # cnstr_pos = [xx >= -1.44]
    cnstr_pos = [xx >= 0]
    cnstr_sum = [cp.sum_entries(cp.mul_elemwise(proj, xx)) == 1. for proj in [[1, 1]]]
    cnstr = cnstr_pos + cnstr_sum

    prob = cp.Problem(obj, cnstr)
    assert prob.is_dcp()

    # prob = cp.Problem(obj)

    # Verify that problem is actually feasible
    xx.value = np.array([[1.], [0]])
    optimal_value = obj.value
    assert optimal_value == 167095032.17051098
    assert [d.value for d in delta] == [-2535.0999993, 12675.49999661]
    assert [_cstr.value for _cstr in cnstr] == [True, True]

    # Solve the problem
    res = prob.solve(verbose=True)
    assert not np.isfinite(res)
    assert prob.status == 'infeasible_inaccurate'

 

Can somebody help me to understand why it fails. Is it a problem in my formulation or an instability of cvxpy?

Thanks in advance.

[Steven Diamond]

Your problem is badly scaled, which is difficult for the solvers to handle. I ran MOSEK on your problem, probably the most robust cone solver, and it failed. Try rescaling A and B. For example, I divided them both by 100 and MOSEK found the optimal value.

[Emmanuel Sérié]

Thank you for your response, effectively, the following modification in the code make the solution optimal (with ECOS backend):
 

   scale = 1.e-3
   cost = cp.quad_form(scale * delta, omega)

In fact, I had initially tested the following solution that does not work:

   scale = 1.e-6
   cost = scale * cp.quad_form(delta, omega)

Is this an expected behavior of cvxpy?

[Steven Diamond]

Yes this is expected behavior, though obviously unintuitive. They yield different cone program formulations.

[Emmanuel Sérié]

Thanks, this is good to know.
Good end of year!