Large discrepacy between optimal solution and DQCP solution

[literarycorner]

I started to experiment with quasiconvexity and DQCP by trying to solve a trivial problem.
$\min_{x} \quad \sqrt x$
$\textrm{subject to:} \quad 1 \leq x \leq 2$

The solution has to be x = 1 and p = 1.
However using cvxpy I get the following:
Optimal value of x: 1.4250764776108187
Optimal objective value: 1.193765671147742

Here is the code I used:

import cvxpy as cp

x = cp.Variable()

objective = cp.Minimize(cp.sqrt(x))

constraints = [x <= 2, x >= 1]

problem = cp.Problem(objective, constraints)

problem.solve(qcp=True)

print("Optimal value of x:", x.value)
print("Optimal objective value:", objective.value)

I am surprised by the big error. Is that to be expected when using cvxpy with DQCP or am I missing something? Can this problem be circumvented?

[rileyjmurray]

I just ran your example and got the same results, which I found the surprising.

Maybe the DQCP --> DCP transformations applied to this problem make it hard to quantify the error of an approximate solution. If that happens then CVXPY's bisection approach could erroneously terminate early, thinking that a solution has been obtained to fairly high accuracy. I wasn't involved in the development of CVXPY's DQCP support though, so I could be wrong about this. @akshayka may be able to provide insight.

[jcert]

I suspect there is something wrong with the DQCP --> DCP transformation.

Using problem.solve(qcp=True, verbose=True) we see that

===============================================================================
                                     CVXPY                                     
                                     v1.4.2                                    
===============================================================================
(CVXPY) Mar 27 10:57:12 PM: Your problem has 1 variables, 2 constraints, and 0 parameters.
(CVXPY) Mar 27 10:57:12 PM: It is compliant with the following grammars: DQCP
(CVXPY) Mar 27 10:57:12 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Mar 27 10:57:12 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
(CVXPY) Mar 27 10:57:12 PM: Your problem is compiled with the CPP canonicalization backend.
(CVXPY) Mar 27 10:57:12 PM: Reducing DQCP problem to a one-parameter family of DCP problems, for bisection.

********************************************************************************
Preparing to bisect problem

minimize 0.0
subject to var1 <= 2.0
           1.0 <= var1
           var1 <= var27
           SOC(reshape(var27 + 1.0, (1,), F), Vstack(reshape(var27 + -1.0, (1, 1), F), reshape(2.0 @ param17, (1, 1), F)))

Finding interval for bisection ...
initial lower bound: 0.000000
initial upper bound: 1.000000

(iteration 0) lower bound: 0.000000
(iteration 0) upper bound: 1.000000
(iteration 0) query point: 0.500000 
(iteration 0) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26992123), 27: array(1.67808185)}, dual_vars={5: 0.0, 9: 0.0, 25: -3.341629027804886e-18, 36: array([0., 0., 0.])}, attr={'solve_time': 1.6381e-05, 'setup_time': 1.8024e-05, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26992123, 1.67808185]), 'y': array([ 0.00000000e+00,  0.00000000e+00, -3.34162903e-18,  0.00000000e+00,
        0.00000000e+00,  0.00000000e+00]), 's': array([0.73007872, 0.2699214 , 0.40816045, 2.67808181, 0.67808188,
       0.99999997]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -0.0, 'res_pri': 1.7201348232776397e-07, 'res_dual': 3.341629027804886e-18, 'gap': 0.0, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 4.158514163920964e-18, 'solve_time': 0.016381, 'setup_time': 0.018024, 'lin_sys_time': 0.0045260000000000005, 'cone_time': 0.002337000000000001, 'accel_time': 0.0011030000000000002, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 5) lower bound: 0.000000
(iteration 5) upper bound: 0.031250
(iteration 5) query point: 0.015625 
(iteration 5) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301817), 27: array(1.70452178)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([0., 0., 0.])}, attr={'solve_time': 1.0330000000000001e-05, 'setup_time': 1.2543999999999999e-05, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301817, 1.70452178]), 'y': array([0., 0., 0., 0., 0., 0.]), 's': array([0.73698177, 0.26301837, 0.4415034 , 2.70452175, 0.70452181,
       0.03125   ]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -0.0, 'res_pri': 2.1227291938616778e-07, 'res_dual': 0.0, 'gap': 0.0, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 0.0, 'solve_time': 0.01033, 'setup_time': 0.012544, 'lin_sys_time': 0.0027009999999999985, 'cone_time': 0.0013240000000000005, 'accel_time': 0.0006540000000000003, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 10) lower bound: 0.000000
(iteration 10) upper bound: 0.000977
(iteration 10) query point: 0.000488 
(iteration 10) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301861), 27: array(1.70452235)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([0., 0., 0.])}, attr={'solve_time': 9.598e-06, 'setup_time': 9.688000000000001e-06, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301861, 1.70452235]), 'y': array([0., 0., 0., 0., 0., 0.]), 's': array([7.36981327e-01, 2.63018814e-01, 4.41503530e-01, 2.70452232e+00,
       7.04522385e-01, 9.76562471e-04]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -0.0, 'res_pri': 2.122716947698384e-07, 'res_dual': 0.0, 'gap': 0.0, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 0.0, 'solve_time': 0.009598, 'setup_time': 0.009688, 'lin_sys_time': 0.0027449999999999987, 'cone_time': 0.0013310000000000006, 'accel_time': 0.0006720000000000003, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 15) lower bound: 0.000000
(iteration 15) upper bound: 0.000031
(iteration 15) query point: 0.000015 
(iteration 15) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301861), 27: array(1.70452236)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([ 0.00000000e+00,  0.00000000e+00, -5.96941776e-22])}, attr={'solve_time': 9.879000000000001e-06, 'setup_time': 1.03e-05, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301861, 1.70452236]), 'y': array([ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,  0.00000000e+00,
        0.00000000e+00, -5.96941776e-22]), 's': array([7.36981327e-01, 2.63018814e-01, 4.41503530e-01, 2.70452233e+00,
       7.04522385e-01, 3.05175772e-05]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': 1.8217217291346015e-26, 'res_pri': 2.1227169430997754e-07, 'res_dual': 0.0, 'gap': 1.8217217291346015e-26, 'res_infeas': 0.0, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 5.501865187369163e-26, 'solve_time': 0.009879, 'setup_time': 0.0103, 'lin_sys_time': 0.0027449999999999987, 'cone_time': 0.0012950000000000006, 'accel_time': 0.0006800000000000003, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 20) lower bound: 0.000000
(iteration 20) upper bound: 0.000001
(iteration 20) query point: 0.000000 
(iteration 20) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301861), 27: array(1.70452236)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([0.00000000e+00, 0.00000000e+00, 1.86544305e-23])}, attr={'solve_time': 1.0069e-05, 'setup_time': 9.869e-06, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301861, 1.70452236]), 'y': array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
       0.00000000e+00, 1.86544305e-23]), 's': array([7.36981327e-01, 2.63018814e-01, 4.41503530e-01, 2.70452233e+00,
       7.04522385e-01, 9.53674288e-07]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -1.7790251261073571e-29, 'res_pri': 2.1227169393003956e-07, 'res_dual': 0.0, 'gap': 1.7790251261073571e-29, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 5.372915222042165e-29, 'solve_time': 0.010069, 'setup_time': 0.009869, 'lin_sys_time': 0.003066999999999999, 'cone_time': 0.0013420000000000005, 'accel_time': 0.0006910000000000003, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 25) lower bound: 0.000000
(iteration 25) upper bound: 0.000000
(iteration 25) query point: 0.000000 
(iteration 25) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301861), 27: array(1.70452236)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([0., 0., 0.])}, attr={'solve_time': 1.3476e-05, 'setup_time': 2.0329e-05, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301861, 1.70452236]), 'y': array([0., 0., 0., 0., 0., 0.]), 's': array([7.36981327e-01, 2.63018814e-01, 4.41503530e-01, 2.70452233e+00,
       7.04522385e-01, 2.98023215e-08]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -0.0, 'res_pri': 2.122716939300392e-07, 'res_dual': 0.0, 'gap': 0.0, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 0.0, 'solve_time': 0.013476, 'setup_time': 0.020329, 'lin_sys_time': 0.0034149999999999988, 'cone_time': 0.0018609999999999994, 'accel_time': 0.0009030000000000004, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

(iteration 30) lower bound: 0.000000
(iteration 30) upper bound: 0.000000
(iteration 30) query point: 0.000000 
(iteration 30) query was feasible. Solution(status=optimal, opt_val=0.0, primal_vars={1: array(1.26301861), 27: array(1.70452236)}, dual_vars={5: 0.0, 9: 0.0, 25: 0.0, 36: array([0., 0., 0.])}, attr={'solve_time': 1.4626999999999999e-05, 'setup_time': 1.7884e-05, 'num_iters': 25, 'solver_specific_stats': {'x': array([1.26301861, 1.70452236]), 'y': array([0., 0., 0., 0., 0., 0.]), 's': array([7.36981327e-01, 2.63018814e-01, 4.41503530e-01, 2.70452233e+00,
       7.04522385e-01, 9.31322547e-10]), 'info': {'status_val': 1, 'iter': 25, 'scale_updates': 0, 'scale': 0.1, 'pobj': 0.0, 'dobj': -0.0, 'res_pri': 2.1227169400601143e-07, 'res_dual': 0.0, 'gap': 0.0, 'res_infeas': nan, 'res_unbdd_a': nan, 'res_unbdd_p': nan, 'comp_slack': 0.0, 'solve_time': 0.014627, 'setup_time': 0.017884, 'lin_sys_time': 0.0036279999999999993, 'cone_time': 0.0020729999999999993, 'accel_time': 0.0009110000000000004, 'rejected_accel_steps': 0, 'accepted_accel_steps': 0, 'status': 'solved'}}}))

Bisection completed, with lower bound 0.000000 and upper bound 0.0000000
********************************************************************************

So the bisection is using the following problem to solve the original DQCP:

$$ \min_{x,z} 0 \quad \text{ s.t. } 1 \leq x \leq 2, \quad x \leq z, \quad \left|\left| \begin{bmatrix}z-1 \\ 2p\end{bmatrix} \right|\right|_2 \leq z+1, $$

where $p$ is the parameter being changed to find when the problem becomes infeasible. The last constraint makes is the same as $p^2-z\leq 0$, so for z large enough it will always be feasible.

By debugging the code I found that in one of the many executions of canonicalize_tree(self, expr) the expression var1 <= power(param17, 2.0) is to be added as a constraint but then it is further modified into

SOC(reshape(var27 + 1.0, (1,), F), Vstack(reshape(var27 + -1.0, (1, 1), F), reshape(2.0 @ param17, (1, 1), F)))

[PTNobel]

var1 <= power(param17, 2.0) is to be added as a constraint but then it is further modified into
SOC(reshape(var27 + 1.0, (1,), F), Vstack(reshape(var27 + -1.0, (1, 1), F), reshape(2.0 @ param17, (1, 1), F)))

This transformation is correct, the second-order constraint is $||(z-1, 2p)||_2 \leq z+1 \iff z^2 -2z +1 + 4p^2 \leq z^2 + 2z + 1 \iff p^2 \leq z$.

[jcert]

Ok, so I must be missing something. We have the problem

$$ \min_{x,z} 0 \quad \text{ s.t. } 1 \leq x \leq 2, \quad x \leq z, \quad p^2\leq z, $$

which will always be feasible, just pick $z$ large enough. Then how does the bisection work?

From debugging, the bisection search is reducing the interval for $p$ from [0,1] to [0,1/2] and ... [0,2**-n], since the problem is always feasible, and eventually stopping and returning some $x$ that satisfies 1 <= x <= 2 and $x \leq z$, but with the objective function being a constant there's nothing "pushing" $z$ to have $p^2\leq z$ binding at the optimum.

[jcert]

I did some more debugging and it seems that there really is something wrong with the DQCP --> DCP transformation, in particular how it deals with expressions of Parameters.

Running

import cvxpy as cp

v = cp.Parameter(value=1.0)

x = cp.Variable()
y = cp.Variable(pos=True)
objective_fn = -cp.sqrt(x) / y
problem = cp.Problem(cp.Minimize(objective_fn), [cp.exp(x) <= y, y<=v])
problem.solve(qcp=True, solver=cp.SCS, verbose=True)

assert problem.is_dqcp()
print("Optimal value: ", problem.value)
print("x: ", x.value)
print("y: ", y.value)

we obtain Optimal value: -0.008138005937416148, x: 6.622888237102037e-05, y: 1.0000131496100515, which looks fine.
Now, changing the constraint y<=v to y<=v**2 (which should yield the same solution) and running it again yields Optimal value: -0.42848064495897403, x: 0.5013713754755035, y: 1.6525269797085378 which is wrong for the problem, but matches the solution given in the DQCP example. Again we see that somewhere along the DQCP --> DCP transformation the constraint $y \leq v^2$ was pretty much removed by introducing an auxiliary variable $z$ and replacing $y \leq v^2$ with { $y\leq z$ , SOC($z+1$,($z-1$, $2v$) )}, that is the same as { $y\leq z$ , $v^2\leq z$} and is not equivalent to $y \leq v^2$. In fact, it will always be satisfied for a large $z$.

[rileyjmurray]

The idea that DQCP and DPP didn't play nice together sounded familar, so I did some digging. So far I've found this point in the release notes when DPP was added in CVXPY 1.1:
https://www.cvxpy.org/updates/index.html#known-issues.

The original example provided by @literarycorner didn't include Parameter objects. @jcert, based on the digging you've done in the DQCP canonicalization, do you think that the same problem with parameters is causing this bug and preventing use of DPP?

[jcert]

The original example provided by @literarycorner didn't include Parameter objects in its DQCP form, but the DQCP --> DCP transformation adds a Parameter to it. The modification of the DQCP example was only used to show more explicitly that expressions of Parameters are being incorrectly modified in a way that sometimes a constraint will get removed.

I believe that the handling of expressions of Parameters is the source of of the bug reported by @literarycorner.

By modifying canonicalize_tree into leaving expressions of parameters and constants as they are I got the correct result for @literarycorner 's problem and to the modification of the DQCP example

def canonicalize_tree(self, expr):
        """Recursively canonicalize an Expression."""
        # TODO don't copy affine expressions?
        if type(expr) == cvxtypes.partial_problem():
            canon_expr, constrs = self.canonicalize_tree(
              expr.args[0].objective.expr)
            for constr in expr.args[0].constraints:
                canon_constr, aux_constr = self.canonicalize_tree(constr)
                constrs += [canon_constr] + aux_constr
        else:
            canon_args = []
            constrs = []
            ###===> modified from here
            if hasattr(expr, 'curvature'):
                if expr.curvature == s.CONSTANT:
                    return expr, constrs
            ###<=== modified to here
            for arg in expr.args:
                #this is where things get weird 
                canon_arg, c = self.canonicalize_tree(arg)
                canon_args += [canon_arg]
                constrs += c
            canon_expr, c = self.canonicalize_expr(expr, canon_args)
            constrs += c
        return canon_expr, constrs

I don't know DPP, is there a reason for which expressions containing just Parameters and constants to be canonicalize? Or would this fix I found be reasonable?

[SteveDiamond]

This fix seems fine, at least at first glance. Can you add some test cases and try it out, including the original problematic case? Add the tests in test_dqcp.py

[chofchof]

The mysterious result of @literarycorner's DQCP code comes from the same reason as @ramonfmir's issue #2165.

After replacing the below two lines of "cvxpy/reductions/canonicalization.py" (line 109)

        if isinstance(expr, Expression) and (
                expr.is_constant() and not expr.parameters()):

into

        if isinstance(expr, Expression) and expr.is_constant():

I could get the following correct result.

Optimal value of x: 0.9999999408964901
Optimal objective value: 0.9999999704482446

[SteveDiamond]

The canonicalization of parameters is necessary for DPP. This particular solution will not work for that reason, but your explanation is very helpful. You could apply the solution to DQCP alone potentially though.

[SteveDiamond]

After further discussion it looks like there is a subtle bug due to the interactions between DPP and DQCP. IMO the best fix for now is to disable DPP when DQCP is active. But we will look into it further and prioritize a fix.

[SteveDiamond]

Fixed by #2424