Use the SCS box cone API

CHANGELOG.md

@@ -7,8 +7,13 @@ All notable changes to this project will be documented in this file.
 ### Added
 
 - Document how to add a new QP solver to the library
+- Example with (box) lower and upper bounds
 - Test case where `lb` XOR `ub` is set
 
+### Changed
+
+- SCS: use the box cone API when lower/upper bounds are set
+
 ## [2.0.0] - 2022/07/05
 
 ### Added

examples/box_inequalities.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2016-2022 Stéphane Caron and the qpsolvers contributors.
+#
+# This file is part of qpsolvers.
+#
+# qpsolvers is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Lesser General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# qpsolvers is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with qpsolvers. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Test one of the available QP solvers on a small problem with box inequalities.
+"""
+
+import random
+
+from time import perf_counter
+
+import numpy as np
+
+from qpsolvers import available_solvers, print_matrix_vector, solve_qp
+
+M = np.array([[1.0, 2.0, 0.0], [-8.0, 3.0, 2.0], [0.0, 1.0, 1.0]])
+P = np.dot(M.T, M)  # this is a positive definite matrix
+q = np.dot(np.array([3.0, 2.0, 3.0]), M)
+A = np.array([1.0, 1.0, 1.0])
+b = np.array([1.0])
+lb = -0.5 * np.ones(3)
+ub = 1.0 * np.ones(3)
+
+if __name__ == "__main__":
+    start_time = perf_counter()
+    solver = random.choice(available_solvers)
+    x = solve_qp(P, q, A=A, b=b, lb=lb, ub=ub, solver=solver)
+    end_time = perf_counter()
+
+    print("")
+    print("    min. 1/2 x^T P x + q^T x")
+    print("    s.t.   A * x == b")
+    print("         lb <= x <= ub")
+    print("")
+    print_matrix_vector(P, "P", q, "q")
+    print("")
+    print_matrix_vector(A, "A", b, "b")
+    print("")
+    print_matrix_vector(lb.reshape((3, 1)), "lb", ub, "ub")
+    print("")
+    print(f"Solution:\n\n    x = {x}\n")
+    print(f"Found in {1e6 * (end_time - start_time):.0f} [us] with {solver}")

qpsolvers/solve_qp.py

@@ -143,11 +143,14 @@ def solve_qp(
     if isinstance(G, ndarray) and G.ndim == 1:
         G = G.reshape((1, G.shape[0]))
     check_problem_constraints(G, h, A, b)
-    G, h = concatenate_bounds(G, h, lb, ub)
     kwargs["initvals"] = initvals
     kwargs["verbose"] = verbose
     try:
-        return solve_function[solver](P, q, G, h, A, b, **kwargs)
+        if (lb is not None or ub is not None) and solver == "scs":
+            return solve_function["scs"](P, q, G, h, A, b, lb, ub, **kwargs)
+        else:  # all other solvers, bounds or not
+            G, h = concatenate_bounds(G, h, lb, ub)
+            return solve_function[solver](P, q, G, h, A, b, **kwargs)
     except KeyError as e:
         raise SolverNotFound(
             f"solver '{solver}' is not in the list "

qpsolvers/solvers/__init__.py

@@ -24,8 +24,7 @@
 
 from numpy import ndarray
 
-from .typing import CvxoptReadyMatrix
-from .typing import DenseOrCSCMatrix
+from .typing import CvxoptReadyMatrix, DenseOrCSCMatrix
 
 available_solvers = []
 dense_solvers = []
@@ -310,6 +309,8 @@
             Optional[DenseOrCSCMatrix],
             Optional[ndarray],
             Optional[ndarray],
+            Optional[ndarray],
+            Optional[ndarray],
             float,
             float,
             bool,

qpsolvers/solvers/scs_.py

@@ -20,17 +20,15 @@
 
 """Solver interface for SCS"""
 
-from typing import Optional
+from typing import Any, Dict, Optional
 from warnings import warn
 
-from numpy import hstack, ndarray
+import numpy as np
 from numpy.linalg import norm
 from scipy import sparse
 from scs import solve
 
-from .typing import DenseOrCSCMatrix
-from .typing import warn_about_sparse_conversion
-
+from .typing import DenseOrCSCMatrix, warn_about_sparse_conversion
 
 # See https://www.cvxgrp.org/scs/api/exit_flags.html#exit-flags
 __status_val_meaning__ = {
@@ -49,17 +47,19 @@
 
 def scs_solve_qp(
     P: DenseOrCSCMatrix,
-    q: ndarray,
+    q: np.ndarray,
     G: Optional[DenseOrCSCMatrix] = None,
-    h: Optional[ndarray] = None,
+    h: Optional[np.ndarray] = None,
     A: Optional[DenseOrCSCMatrix] = None,
-    b: Optional[ndarray] = None,
-    initvals: Optional[ndarray] = None,
+    b: Optional[np.ndarray] = None,
+    lb: Optional[np.ndarray] = None,
+    ub: Optional[np.ndarray] = None,
+    initvals: Optional[np.ndarray] = None,
     eps_abs: float = 1e-7,
     eps_rel: float = 1e-7,
     verbose: bool = False,
     **kwargs,
-) -> Optional[ndarray]:
+) -> Optional[np.ndarray]:
     """
     Solve a Quadratic Program defined as:
 
@@ -70,7 +70,8 @@ def scs_solve_qp(
             \\frac{1}{2} x^T P x + q^T x \\\\
         \\mbox{subject to}
             & G x \\leq h                \\\\
-            & A x = b
+            & A x = b                    \\\\
+            & lb \\leq x \\leq ub
         \\end{array}\\end{split}
 
     using `SCS <https://github.com/cvxgrp/scs>`_.
@@ -89,6 +90,10 @@ def scs_solve_qp(
         Linear equality constraint matrix.
     b :
         Linear equality constraint vector.
+    lb:
+        Lower bound constraint vector.
+    ub:
+        Upper bound constraint vector.
     initvals :
         Warm-start guess vector (not used).
     eps_abs : float
@@ -117,13 +122,13 @@ def scs_solve_qp(
     that SCS behaves closer to the other solvers. If you don't need that much
     precision, increase them for better performance.
     """
-    if isinstance(P, ndarray):
+    if isinstance(P, np.ndarray):
         warn_about_sparse_conversion("P")
         P = sparse.csc_matrix(P)
-    if isinstance(G, ndarray):
+    if isinstance(G, np.ndarray):
         warn_about_sparse_conversion("G")
         G = sparse.csc_matrix(G)
-    if isinstance(A, ndarray):
+    if isinstance(A, np.ndarray):
         warn_about_sparse_conversion("A")
         A = sparse.csc_matrix(A)
     kwargs.update(
@@ -133,24 +138,24 @@ def scs_solve_qp(
             "verbose": verbose,
         }
     )
-    data = {"P": P, "c": q}
-    cone = {}
+    data: Dict[str, Any] = {"P": P, "c": q}
+    cone: Dict[str, Any] = {}
     if initvals is not None:
         data["x"] = initvals
     if A is not None and b is not None:
         if G is not None and h is not None:
             data["A"] = sparse.vstack([A, G], format="csc")
-            data["b"] = hstack([b, h])
+            data["b"] = np.hstack([b, h])
             cone["z"] = b.shape[0]  # zero cone
-            cone["l"] = h.shape[0]  # positive orthant
-        else:  # A is not None and b is not None
+            cone["l"] = h.shape[0]  # positive cone
+        else:  # G is None and h is None
             data["A"] = A
             data["b"] = b
             cone["z"] = b.shape[0]  # zero cone
     elif G is not None and h is not None:
         data["A"] = G
         data["b"] = h
-        cone["l"] = h.shape[0]  # positive orthant
+        cone["l"] = h.shape[0]  # positive cone
     else:  # no constraint
         x = sparse.linalg.lsqr(P, -q)[0]
         if norm(P @ x + q) > 1e-9:
@@ -159,6 +164,16 @@ def scs_solve_qp(
                 "q has component in the nullspace of P"
             )
         return x
+    if lb is not None or ub is not None:
+        n = P.shape[1]
+        cone["bl"] = lb if lb is not None else np.full((n,), -np.inf)
+        cone["bu"] = ub if ub is not None else np.full((n,), +np.inf)
+        zero_row = sparse.csc_matrix((1, n))
+        data["A"] = sparse.vstack(
+            (data["A"], zero_row, -sparse.eye(n)),
+            format="csc",
+        )
+        data["b"] = np.hstack((data["b"], 1.0, np.zeros(n)))
     solution = solve(data, cone, **kwargs)
     status_val = solution["info"]["status_val"]
     if status_val != 1: