[HIGHS conif] Maps dual solution back to dual variables (cvxpy#2838)

* [HIGHS conif] Map dual values back to dual variables

* fix tests and dotsort

* fix CI

* Add cvar test

* fix dotsort canonicalization

* remove variable bounds handling

---------

Co-authored-by: Steven Diamond <steven@gridmatic.com>

Author: senhalil

cvxpy/atoms/dotsort.py

@@ -133,6 +133,11 @@ def get_data(self):
         """
         return None
 
+    def is_pwl(self) -> bool:
+        """Is the expression piecewise linear?
+        """
+        return True
+
     @staticmethod
     def _get_args_from_values(values: List[np.ndarray]) \
             -> Tuple[np.ndarray, np.ndarray]:

cvxpy/reductions/eliminate_pwl/canonicalizers/dotsort_canon.py

@@ -35,10 +35,10 @@ def dotsort_canon(expr, args):
     # subject to  x @ w_unique.T <= t + q.T
     #             0 <= t
 
-    t = Variable((x.size, 1), nonneg=True)
+    t = Variable((x.size, 1))
     q = Variable((1, w_unique.size))
 
     obj = sum(t) + q @ w_counts
     x_w_unique_outer_product = outer(vec(x, order='F'), vec(w_unique, order='F'))
-    constraints = [x_w_unique_outer_product <= t + q]
+    constraints = [x_w_unique_outer_product <= t + q, t >= 0]
     return obj, constraints

cvxpy/reductions/qp2quad_form/canonicalizers/__init__.py

@@ -30,6 +30,7 @@
 CANON_METHODS = {}
 CANON_METHODS[abs] = CONE_METHODS[abs]
 CANON_METHODS[cumsum] = CONE_METHODS[cumsum]
+CANON_METHODS[dotsort] = CONE_METHODS[dotsort]
 CANON_METHODS[maximum] = CONE_METHODS[maximum]
 CANON_METHODS[minimum] = CONE_METHODS[minimum]
 CANON_METHODS[sum_largest] = CONE_METHODS[sum_largest]

cvxpy/reductions/solvers/conic_solvers/highs_conif.py

@@ -20,6 +20,7 @@
 import cvxpy.settings as s
 from cvxpy.error import SolverError
 from cvxpy.reductions.solution import Solution, failure_solution
+from cvxpy.reductions.solvers import utilities
 from cvxpy.reductions.solvers.conic_solvers.conic_solver import (
     ConicSolver,
     dims_to_solver_dict,
@@ -104,14 +105,26 @@ def invert(self, results, inverse_data):
         if status in s.SOLUTION_PRESENT:
             opt_val = results["info"].objective_function_value + inverse_data[s.OFFSET]
             primal_vars = {
+                # inverse_data[HIGHS.VAR_ID]: ...
+                # I don't understand how the line below works, the other conif solvers have
+                # something similar to the commented line above for the "key".
                 HIGHS.VAR_ID: intf.DEFAULT_INTF.const_to_matrix(
                     np.array(results["solution"].col_value)
                 )
             }
             # add duals if not a MIP.
             dual_vars = None
-            if not inverse_data["is_mip"]:
-                dual_vars = {HIGHS.DUAL_VAR_ID: -np.array(results["solution"].row_dual)}
+            if not inverse_data['is_mip']:
+                # The dual values are retrieved in the order that the
+                # constraints were added in solve_via_data() below. We
+                # must be careful to map them to inverse_data[EQ_CONSTR]
+                # followed by inverse_data[NEQ_CONSTR] accordingly.
+                y = -np.array(results["solution"].row_dual)
+                dual_vars = utilities.get_dual_values(
+                    y,
+                    utilities.extract_dual_value,
+                    inverse_data[HIGHS.EQ_CONSTR] + inverse_data[HIGHS.NEQ_CONSTR])
+
             attr[s.NUM_ITERS] = (
                 results["info"].ipm_iteration_count
                 + results["info"].crossover_iteration_count

cvxpy/reductions/solvers/solving_chain.py

@@ -260,6 +260,7 @@ def construct_solving_chain(problem, candidates,
         # Canonicalize as a QP
         solver = candidates['qp_solvers'][0]
         solver_instance = slv_def.SOLVER_MAP_QP[solver]
+        # TODO should CvxAttr2Constr come after qp2symbolic_qp?
         reductions += [
             CvxAttr2Constr(reduce_bounds=not solver_instance.BOUNDED_VARIABLES),
             qp2symbolic_qp.Qp2SymbolicQp(),

cvxpy/tests/test_conic_solvers.py

@@ -2154,7 +2154,12 @@ class TestHIGHS:
         ],
     )
     def test_highs_solving(self, problem) -> None:
-        problem(solver=cp.HIGHS)
+        # HACK needed to use the HiGHS conic interface rather than 
+        # the QP interface for LPs.
+        from cvxpy.reductions.solvers.conic_solvers.highs_conif import HIGHS
+        solver = HIGHS()
+        solver.name = lambda: "HIGHS CONIC"
+        problem(solver=solver)
 
     @pytest.mark.parametrize(
         ["problem", "confirmation_string"],

cvxpy/tests/test_qp_solvers.py

@@ -507,6 +507,34 @@ def test_highs_warmstart(self) -> None:
             result2 = prob.solve(solver=cp.HIGHS, warm_start=False)
             self.assertAlmostEqual(result, result2)
 
+    def test_highs_cvar(self) -> None:
+        """Test problem with CVaR constraint from 
+        https://github.com/cvxpy/cvxpy/issues/2836
+        """
+        if cp.HIGHS in INSTALLED_SOLVERS:
+            # Generate data
+            num_stocks = 5
+            num_samples = 25
+            np.random.seed(1)
+            pnl_samples = np.random.uniform(low=0.0, high=1.0, size=(num_samples, num_stocks))
+            pnl_expected = pnl_samples.mean(axis=0)
+
+            # Prepare to solve
+            quantile = 0.05
+            w = cp.Variable(num_stocks, nonneg=True)
+            cvar = cp.cvar(pnl_samples @ w, 1 - quantile)
+            pnl = w @ pnl_expected
+
+            # Solve
+            objective = cp.Maximize(pnl)
+            constraints = [cvar <= 0.5]
+            problem = cp.Problem(objective, constraints)
+            problem.solve(
+                solver=cp.HIGHS,
+            )
+            assert problem.status == cp.OPTIMAL
+
+
     def test_piqp_warmstart(self) -> None:
         """Test warm start.
         """