From 5a444a0c24a3d4fe04a655cc5b2c553baecf14b7 Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Thu, 18 Apr 2024 23:35:08 +0200 Subject: [PATCH 1/7] move subproblem to single file --- src/ncopt/sqpgs/cvxpy_subproblem.py | 168 +++++++++++++++++++++++++++ src/ncopt/sqpgs/main.py | 171 +--------------------------- tests/test_subproblem.py | 14 +-- 3 files changed, 179 insertions(+), 174 deletions(-) create mode 100644 src/ncopt/sqpgs/cvxpy_subproblem.py diff --git a/src/ncopt/sqpgs/cvxpy_subproblem.py b/src/ncopt/sqpgs/cvxpy_subproblem.py new file mode 100644 index 0000000..28ca74a --- /dev/null +++ b/src/ncopt/sqpgs/cvxpy_subproblem.py @@ -0,0 +1,168 @@ +from typing import List + +import cvxpy as cp +import numpy as np + +from ncopt.sqpgs.defaults import DEFAULT_OPTION + +CVXPY_SOLVER_DICT = {"osqp": cp.OSQP, "clarabel": cp.CLARABEL} + + +class CVXPYSubproblemSQPGS: + def __init__( + self, + dim: int, + p0: np.ndarray, + pI: np.ndarray, + pE: np.ndarray, + assert_tol: float, + solver: str = DEFAULT_OPTION.qp_solver, + ) -> None: + """ + dim : solution space dimension + p0 : number of sample points for f (excluding x_k itself) + pI : array, number of sample points for inequality constraint (excluding x_k itself) + pE : array, number of sample points for equality constraint (excluding x_k itself) + """ + + self.dim = dim + self.p0 = p0 + self.pI = pI + self.pE = pE + + self.assert_tol = assert_tol + + self.d = cp.Variable(self.dim) + self._problem = None + self._qp_solver = CVXPY_SOLVER_DICT.get(solver, cp.OSQP) + + @property + def nI(self) -> int: + return len(self.pI) + + @property + def nE(self) -> int: + return len(self.pE) + + @property + def has_ineq_constraints(self) -> bool: + return self.nI > 0 + + @property + def has_eq_constraints(self) -> bool: + return self.nE > 0 + + @property + def problem(self) -> cp.Problem: + assert self._problem is not None, "Problem not yet initialized." + return self._problem + + @property + def status(self) -> str: + return self.problem.status + + @property + def objective_val(self) -> float: + return self.problem.value + + @property + def setup_time(self) -> float: + return self.problem.solver_stats.setup_time + + @property + def solve_time(self) -> float: + return self.problem.solver_stats.solve_time + + def solve( + self, + L: np.ndarray, + rho: float, + D_f: np.ndarray, + D_gI: List[np.ndarray], + D_gE: List[np.ndarray], + f_k: float, + gI_k: np.ndarray, + gE_k: np.ndarray, + ) -> None: + """ + This solves the quadratic program + + Parameters + + L : array + Cholesky factor of Hessian approximation + rho : float + parameter + D_f : array + gradient of f at the sampled points + D_gI : list + j-th element is the gradient array of c^j at the sampled points. + D_gE : list + j-th element is the gradient array of h^j at the sampled points. + f_k : float + evaluation of f at x_k. + gI_k : array + evaluation of inequality constraints at x_k. + gE_k : array + evaluation of equality constraints at x_k. + + Updates + self.d: Variable + search direction + + self.lambda_f: array + KKT multipier for objective. + + self.lambda_gE: list + KKT multipier for equality constraints. + + self.lambda_gI: list + KKT multipier for inequality constraints. + + """ + + d = self.d + z = cp.Variable() + if self.has_ineq_constraints: + r_I = cp.Variable(gI_k.size, nonneg=True) + if self.has_eq_constraints: + r_E = cp.Variable(gE_k.size, nonneg=True) + + objective = rho * z + (1 / 2) * cp.sum_squares(L.T @ d) + + obj_constraint = f_k + D_f @ d <= z + constraints = [obj_constraint] + + if self.has_ineq_constraints: + ineq_constraints = [gI_k[j] + D_gI[j] @ d <= r_I[j] for j in range(self.nI)] + constraints += ineq_constraints + objective = objective + cp.sum(r_I) + + if self.has_eq_constraints: + eq_constraints_plus = [gE_k[j] + D_gE[j] @ d <= r_E[j] for j in range(self.nE)] + eq_constraints_neg = [gE_k[j] + D_gE[j] @ d >= r_E[j] for j in range(self.nE)] + constraints += eq_constraints_plus + eq_constraints_neg + objective = objective + cp.sum(r_E) + + problem = cp.Problem(cp.Minimize(objective), constraints) + problem.solve(solver=self._qp_solver, verbose=False) + + assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE} + self._problem = problem + + # Extract dual variables + duals = problem.solution.dual_vars + self.lambda_f = duals[obj_constraint.id] + + if self.has_ineq_constraints: + self.lambda_gI = [duals[c.id] for c in ineq_constraints] + else: + self.lambda_gI = [] + + if self.has_eq_constraints: + self.lambda_gE = [ + duals[c_plus.id] - duals[c_neg.id] + for c_plus, c_neg in zip(eq_constraints_plus, eq_constraints_neg) + ] + else: + self.lambda_gE = [] diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py index dcabfc7..0130717 100644 --- a/src/ncopt/sqpgs/main.py +++ b/src/ncopt/sqpgs/main.py @@ -14,12 +14,12 @@ import time from typing import List, Optional, Tuple, Union -import cvxpy as cp import numpy as np import torch from ncopt.functions import ObjectiveOrConstraint from ncopt.plot_utils import plot_metrics, plot_timings +from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS from ncopt.sqpgs.defaults import DEFAULT_ARG, DEFAULT_OPTION from ncopt.utils import compute_batch_jacobian_vmap, get_logger @@ -141,7 +141,9 @@ def solve(self): pE = np.repeat(pE_, self.dimE) ############################################################### - self.SP = SubproblemSQPGS(self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"]) + self.SP = CVXPYSubproblemSQPGS( + self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"] + ) E_k = np.inf # for stopping criterion x_kmin1 = None # last iterate @@ -503,168 +505,3 @@ def stop_criterion(gI, gE, g_k, SP, gI_k, gE_k, V_gI, V_gE, nI_, nE_): val5 = np.maximum(val5, np.max(SP.lambda_gE[j] * gE_vals[j])) return np.max(np.array([val1, val2, val3, val4, val5])) - - -# %% - -CVXPY_SOLVER_DICT = {"osqp": cp.OSQP, "clarabel": cp.CLARABEL} - - -class SubproblemSQPGS: - def __init__( - self, - dim: int, - p0: np.ndarray, - pI: np.ndarray, - pE: np.ndarray, - assert_tol: float, - solver: str = DEFAULT_OPTION.qp_solver, - ) -> None: - """ - dim : solution space dimension - p0 : number of sample points for f (excluding x_k itself) - pI : array, number of sample points for inequality constraint (excluding x_k itself) - pE : array, number of sample points for equality constraint (excluding x_k itself) - """ - - self.dim = dim - self.p0 = p0 - self.pI = pI - self.pE = pE - - self.assert_tol = assert_tol - - self.d = cp.Variable(self.dim) - self._problem = None - self._qp_solver = CVXPY_SOLVER_DICT.get(solver, cp.OSQP) - - @property - def nI(self) -> int: - return len(self.pI) - - @property - def nE(self) -> int: - return len(self.pE) - - @property - def has_ineq_constraints(self) -> bool: - return self.nI > 0 - - @property - def has_eq_constraints(self) -> bool: - return self.nE > 0 - - @property - def problem(self) -> cp.Problem: - assert self._problem is not None, "Problem not yet initialized." - return self._problem - - @property - def status(self) -> str: - return self.problem.status - - @property - def objective_val(self) -> float: - return self.problem.value - - @property - def setup_time(self) -> float: - return self.problem.solver_stats.setup_time - - @property - def solve_time(self) -> float: - return self.problem.solver_stats.solve_time - - def solve( - self, - L: np.ndarray, - rho: float, - D_f: np.ndarray, - D_gI: List[np.ndarray], - D_gE: List[np.ndarray], - f_k: float, - gI_k: np.ndarray, - gE_k: np.ndarray, - ) -> None: - """ - This solves the quadratic program - - Parameters - - L : array - Cholesky factor of Hessian approximation - rho : float - parameter - D_f : array - gradient of f at the sampled points - D_gI : list - j-th element is the gradient array of c^j at the sampled points. - D_gE : list - j-th element is the gradient array of h^j at the sampled points. - f_k : float - evaluation of f at x_k. - gI_k : array - evaluation of inequality constraints at x_k. - gE_k : array - evaluation of equality constraints at x_k. - - Updates - self.d: Variable - search direction - - self.lambda_f: array - KKT multipier for objective. - - self.lambda_gE: list - KKT multipier for equality constraints. - - self.lambda_gI: list - KKT multipier for inequality constraints. - - """ - - d = self.d - z = cp.Variable() - if self.has_ineq_constraints: - r_I = cp.Variable(gI_k.size, nonneg=True) - if self.has_eq_constraints: - r_E = cp.Variable(gE_k.size, nonneg=True) - - objective = rho * z + (1 / 2) * cp.sum_squares(L.T @ d) - - obj_constraint = f_k + D_f @ d <= z - constraints = [obj_constraint] - - if self.has_ineq_constraints: - ineq_constraints = [gI_k[j] + D_gI[j] @ d <= r_I[j] for j in range(self.nI)] - constraints += ineq_constraints - objective = objective + cp.sum(r_I) - - if self.has_eq_constraints: - eq_constraints_plus = [gE_k[j] + D_gE[j] @ d <= r_E[j] for j in range(self.nE)] - eq_constraints_neg = [gE_k[j] + D_gE[j] @ d >= r_E[j] for j in range(self.nE)] - constraints += eq_constraints_plus + eq_constraints_neg - objective = objective + cp.sum(r_E) - - problem = cp.Problem(cp.Minimize(objective), constraints) - problem.solve(solver=self._qp_solver, verbose=False) - - assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE} - self._problem = problem - - # Extract dual variables - duals = problem.solution.dual_vars - self.lambda_f = duals[obj_constraint.id] - - if self.has_ineq_constraints: - self.lambda_gI = [duals[c.id] for c in ineq_constraints] - else: - self.lambda_gI = [] - - if self.has_eq_constraints: - self.lambda_gE = [ - duals[c_plus.id] - duals[c_neg.id] - for c_plus, c_neg in zip(eq_constraints_plus, eq_constraints_neg) - ] - else: - self.lambda_gE = [] diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py index 8280a00..bf1e075 100644 --- a/tests/test_subproblem.py +++ b/tests/test_subproblem.py @@ -2,21 +2,21 @@ import numpy as np import pytest -from ncopt.sqpgs.main import SubproblemSQPGS +from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS @pytest.fixture -def subproblem_ineq() -> SubproblemSQPGS: +def subproblem_ineq() -> CVXPYSubproblemSQPGS: dim = 2 p0 = 2 pI = np.array([3]) pE = np.array([], dtype=int) assert_tol = 1e-5 - subproblem = SubproblemSQPGS(dim, p0, pI, pE, assert_tol) + subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol) return subproblem -def test_subproblem_ineq(subproblem_ineq: SubproblemSQPGS): +def test_subproblem_ineq(subproblem_ineq: CVXPYSubproblemSQPGS): D_f = np.array([[-2.0, 1.0], [-2.04236205, -1.0], [-1.92172864, -1.0]]) D_gI = [np.array([[0.0, 2.0], [0.0, 2.0], [1.41421356, 0.0], [1.41421356, 0.0]])] subproblem_ineq.solve( @@ -34,17 +34,17 @@ def test_subproblem_ineq(subproblem_ineq: SubproblemSQPGS): @pytest.fixture -def subproblem_eq() -> SubproblemSQPGS: +def subproblem_eq() -> CVXPYSubproblemSQPGS: dim = 2 p0 = 2 pI = np.array([], dtype=int) pE = np.array([4, 4]) assert_tol = 1e-5 - subproblem = SubproblemSQPGS(dim, p0, pI, pE, assert_tol) + subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol) return subproblem -def test_subproblem_eq(subproblem_eq: SubproblemSQPGS): +def test_subproblem_eq(subproblem_eq: CVXPYSubproblemSQPGS): D_gE = [ np.array([[1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0], [1.0, 0.0]]), np.array([[0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]), From e0e10f3bbdb2e3950bf8f380fe347dd2c6e5e9c3 Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Fri, 19 Apr 2024 14:29:56 +0200 Subject: [PATCH 2/7] direct call initial --- src/ncopt/sqpgs/osqp_subproblem.py | 304 +++++++++++++++++++++++++++++ 1 file changed, 304 insertions(+) create mode 100644 src/ncopt/sqpgs/osqp_subproblem.py diff --git a/src/ncopt/sqpgs/osqp_subproblem.py b/src/ncopt/sqpgs/osqp_subproblem.py new file mode 100644 index 0000000..3baa9f8 --- /dev/null +++ b/src/ncopt/sqpgs/osqp_subproblem.py @@ -0,0 +1,304 @@ +from typing import List + +import numpy as np +import osqp +from scipy import sparse + + +class OSQPSubproblemSQPGS: + def __init__( + self, + dim: int, + p0: np.ndarray, + pI: np.ndarray, + pE: np.ndarray, + assert_tol: float, + ) -> None: + """ + dim : solution space dimension + p0 : number of sample points for f (excluding x_k itself) + pI : array, number of sample points for inequality constraint (excluding x_k itself) + pE : array, number of sample points for equality constraint (excluding x_k itself) + """ + + self.dim = dim + self.p0 = p0 + self.pI = pI + self.pE = pE + + self.assert_tol = assert_tol + + self._problem = None + self.P, self.q, self.inG, self.inh, self.nonnegG, self.nonnegh = self._initialize() + + @property + def nI(self) -> int: + return len(self.pI) + + @property + def nE(self) -> int: + return len(self.pE) + + @property + def has_ineq_constraints(self) -> bool: + return self.nI > 0 + + @property + def has_eq_constraints(self) -> bool: + return self.nE > 0 + + @property + def problem(self) -> osqp.interface.OSQP: + assert self._problem is not None, "Problem not yet initialized." + return self._problem + + @property + def status(self) -> str: + return self.problem.status + + @property + def objective_val(self) -> float: + return self.problem.value + + @property + def setup_time(self) -> float: + return self.problem.solver_stats.setup_time + + @property + def solve_time(self) -> float: + return self.problem.solver_stats.solve_time + + def _initialize(self): + """ + The quadratic subrpoblem we solve in every iteration is of the form: + + min_y 1/2 * y.T @ P @ y + q.T @ y subject to G @ y <= h + + variable structure: y = (d,z,rI,rE) with + d = search direction + z = helper variable for objective + rI = helper variable for inequality constraints + rI = helper variable for equality constraints + + This function initializes the variables P,q,G,h. + The entries which change in every iteration are then updated in self.update() + + G and h consist of two parts: + 1) inG, inh: the inequalities from the paper + 2) nonnegG, nonnegh: nonnegativity bounds rI >= 0, rE >= 0 + """ + + dimQP = self.dim + 1 + self.nI + self.nE + + P = np.zeros((dimQP, dimQP)) + q = np.zeros(dimQP) + + inG = np.zeros((1 + self.p0 + np.sum(1 + self.pI) + 2 * np.sum(1 + self.pE), dimQP)) + inh = np.zeros(1 + self.p0 + np.sum(1 + self.pI) + 2 * np.sum(1 + self.pE)) + + # structure of inG (p0+1, sum(1+pI), sum(1+pE), sum(1+pE)) + inG[: self.p0 + 1, self.dim] = -1 + + for j in range(self.nI): + inG[ + self.p0 + 1 + (1 + self.pI)[:j].sum() : self.p0 + + 1 + + (1 + self.pI)[:j].sum() + + self.pI[j] + + 1, + self.dim + 1 + j, + ] = -1 + + for j in range(self.nE): + inG[ + self.p0 + 1 + (1 + self.pI).sum() + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1, + self.dim + 1 + self.nI + j, + ] = -1 + inG[ + self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1, + self.dim + 1 + self.nI + j, + ] = +1 + + # we have nI+nE r-variables + nonnegG = np.hstack( + (np.zeros((self.nI + self.nE, self.dim + 1)), -np.eye(self.nI + self.nE)) + ) + nonnegh = np.zeros(self.nI + self.nE) + + return P, q, inG, inh, nonnegG, nonnegh + + def _update(self, H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k): + self.P[: self.dim, : self.dim] = H + self.q = np.hstack((np.zeros(self.dim), rho, np.ones(self.nI), np.ones(self.nE))) + + self.inG[: self.p0 + 1, : self.dim] = D_f + self.inh[: self.p0 + 1] = -f_k + + for j in range(self.nI): + self.inG[ + self.p0 + 1 + (1 + self.pI)[:j].sum() : self.p0 + + 1 + + (1 + self.pI)[:j].sum() + + self.pI[j] + + 1, + : self.dim, + ] = D_gI[j] + self.inh[ + self.p0 + 1 + (1 + self.pI)[:j].sum() : self.p0 + + 1 + + (1 + self.pI)[:j].sum() + + self.pI[j] + + 1 + ] = -gI_k[j] + + for j in range(self.nE): + self.inG[ + self.p0 + 1 + (1 + self.pI).sum() + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1, + : self.dim, + ] = D_gE[j] + self.inG[ + self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1, + : self.dim, + ] = -D_gE[j] + + self.inh[ + self.p0 + 1 + (1 + self.pI).sum() + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1 + ] = -gE_k[j] + self.inh[ + self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() : self.p0 + + 1 + + (1 + self.pI).sum() + + (1 + self.pE).sum() + + (1 + self.pE)[:j].sum() + + self.pE[j] + + 1 + ] = gE_k[j] + + return + + def solve( + self, + H: np.ndarray, + rho: float, + D_f: np.ndarray, + D_gI: List[np.ndarray], + D_gE: List[np.ndarray], + f_k: float, + gI_k: np.ndarray, + gE_k: np.ndarray, + ) -> None: + """ + This solves the quadratic program. + + Updates + self.d: array + search direction + + self.lambda_f: array + KKT multipier for objective. + + self.lambda_gE: list + KKT multipier for equality constraints. + + self.lambda_gI: list + KKT multipier for inequality constraints. + + """ + + # update params + self._update(H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + + A = np.vstack((self.inG, self.nonnegG)) + u = np.hstack((self.inh, self.nonnegh)) + l = np.ones_like(u) * (-np.inf) + + problem = osqp.OSQP() + + problem.setup(P=sparse.csc_matrix(self.P), q=self.q, A=sparse.csc_matrix(A), l=l, u=u) + + res = problem.solve() + # TODO: assert status, where is it stored in OSQP? + self._problem = problem + + primal_solution = res.x + + self.d = primal_solution[: self.dim] + self.z = primal_solution[self.dim] + + self.rI = primal_solution[self.dim + 1 : self.dim + 1 + self.nI] + self.rE = primal_solution[self.dim + 1 + self.nI :] + + # TODO: pipe through assert_tol + assert len(self.rE) == self.nE + assert len(self.rI) == self.nI + # assert np.all(self.rI >= -1e-5) , f"{self.rI}" + # assert np.all(self.rE >= -1e-5), f"{self.rE}" + + # extract dual variables = KKT multipliers + dual_solution = res.y + lambda_f = dual_solution[: self.p0 + 1] + + lambda_gI = list() + for j in np.arange(self.nI): + start_ix = self.p0 + 1 + (1 + self.pI)[:j].sum() + lambda_gI.append(dual_solution[start_ix : start_ix + 1 + self.pI[j]]) + + lambda_gE = list() + for j in np.arange(self.nE): + start_ix = self.p0 + 1 + (1 + self.pI).sum() + (1 + self.pE)[:j].sum() + + # from ineq with + + vec1 = dual_solution[start_ix : start_ix + 1 + self.pE[j]] + + # from ineq with - + vec2 = dual_solution[ + start_ix + (1 + self.pE).sum() : start_ix + (1 + self.pE).sum() + 1 + self.pE[j] + ] + + # see Direction.m line 620 + lambda_gE.append(vec1 - vec2) + + self.lambda_f = lambda_f.copy() + self.lambda_gI = lambda_gI.copy() + self.lambda_gE = lambda_gE.copy() + + return From bc46c5d444f3ec0f15fc6d18526c1c8a27571a6d Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Fri, 19 Apr 2024 15:29:48 +0200 Subject: [PATCH 3/7] one test fails due to dual infeasible --- src/ncopt/sqpgs/cvxpy_subproblem.py | 5 +++- src/ncopt/sqpgs/main.py | 18 +++++++++---- src/ncopt/sqpgs/osqp_subproblem.py | 34 ++++++++++-------------- tests/test_subproblem.py | 40 +++++++++++++++++++++++++++++ 4 files changed, 71 insertions(+), 26 deletions(-) diff --git a/src/ncopt/sqpgs/cvxpy_subproblem.py b/src/ncopt/sqpgs/cvxpy_subproblem.py index 28ca74a..3d0f6a1 100644 --- a/src/ncopt/sqpgs/cvxpy_subproblem.py +++ b/src/ncopt/sqpgs/cvxpy_subproblem.py @@ -121,7 +121,7 @@ def solve( """ - d = self.d + d = cp.Variable(self.dim) z = cp.Variable() if self.has_ineq_constraints: r_I = cp.Variable(gI_k.size, nonneg=True) @@ -150,6 +150,9 @@ def solve( assert problem.status in {cp.OPTIMAL, cp.OPTIMAL_INACCURATE} self._problem = problem + # Extract primal solution + self.d = d.value.copy() + # Extract dual variables duals = problem.solution.dual_vars self.lambda_f = duals[obj_constraint.id] diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py index 0130717..a5a3e37 100644 --- a/src/ncopt/sqpgs/main.py +++ b/src/ncopt/sqpgs/main.py @@ -21,6 +21,7 @@ from ncopt.plot_utils import plot_metrics, plot_timings from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS from ncopt.sqpgs.defaults import DEFAULT_ARG, DEFAULT_OPTION +from ncopt.sqpgs.osqp_subproblem import OSQPSubproblemSQPGS from ncopt.utils import compute_batch_jacobian_vmap, get_logger @@ -141,9 +142,12 @@ def solve(self): pE = np.repeat(pE_, self.dimE) ############################################################### - self.SP = CVXPYSubproblemSQPGS( - self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"] - ) + if self.options["qp_solver"] == "osqp": + self.SP = OSQPSubproblemSQPGS(self.dim, p0, pI, pE, self.assert_tol) + else: + self.SP = CVXPYSubproblemSQPGS( + self.dim, p0, pI, pE, self.assert_tol, self.options["qp_solver"] + ) E_k = np.inf # for stopping criterion x_kmin1 = None # last iterate @@ -235,10 +239,14 @@ def solve(self): # Subproblem solve ############################################## t1 = time.perf_counter() - self.SP.solve(np.linalg.cholesky(H), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + if isinstance(self.SP, OSQPSubproblemSQPGS): + self.SP.solve(H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + else: + self.SP.solve(np.linalg.cholesky(H), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + + d_k = self.SP.d.copy() t2 = time.perf_counter() - d_k = self.SP.d.value.copy() # compute g_k from paper g_k = ( self.SP.lambda_f @ D_f diff --git a/src/ncopt/sqpgs/osqp_subproblem.py b/src/ncopt/sqpgs/osqp_subproblem.py index 3baa9f8..f614082 100644 --- a/src/ncopt/sqpgs/osqp_subproblem.py +++ b/src/ncopt/sqpgs/osqp_subproblem.py @@ -56,18 +56,6 @@ def problem(self) -> osqp.interface.OSQP: def status(self) -> str: return self.problem.status - @property - def objective_val(self) -> float: - return self.problem.value - - @property - def setup_time(self) -> float: - return self.problem.solver_stats.setup_time - - @property - def solve_time(self) -> float: - return self.problem.solver_stats.solve_time - def _initialize(self): """ The quadratic subrpoblem we solve in every iteration is of the form: @@ -253,10 +241,22 @@ def solve( problem = osqp.OSQP() - problem.setup(P=sparse.csc_matrix(self.P), q=self.q, A=sparse.csc_matrix(A), l=l, u=u) + problem.setup( + P=sparse.csc_matrix(self.P), + q=self.q, + A=sparse.csc_matrix(A), + l=l, + u=u, + eps_abs=1e-05, + eps_rel=1e-05, + verbose=False, + polish=1, + ) res = problem.solve() - # TODO: assert status, where is it stored in OSQP? + + assert res.info.status not in ["unsolved"] + print(res.info.status) self._problem = problem primal_solution = res.x @@ -267,12 +267,6 @@ def solve( self.rI = primal_solution[self.dim + 1 : self.dim + 1 + self.nI] self.rE = primal_solution[self.dim + 1 + self.nI :] - # TODO: pipe through assert_tol - assert len(self.rE) == self.nE - assert len(self.rI) == self.nI - # assert np.all(self.rI >= -1e-5) , f"{self.rI}" - # assert np.all(self.rE >= -1e-5), f"{self.rE}" - # extract dual variables = KKT multipliers dual_solution = res.y lambda_f = dual_solution[: self.p0 + 1] diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py index bf1e075..5799f7d 100644 --- a/tests/test_subproblem.py +++ b/tests/test_subproblem.py @@ -3,6 +3,7 @@ import pytest from ncopt.sqpgs.cvxpy_subproblem import CVXPYSubproblemSQPGS +from ncopt.sqpgs.osqp_subproblem import OSQPSubproblemSQPGS @pytest.fixture @@ -61,3 +62,42 @@ def test_subproblem_eq(subproblem_eq: CVXPYSubproblemSQPGS): ) assert subproblem_eq.status == cp.OPTIMAL assert np.isclose(subproblem_eq.objective_val, 1.0) + + +def test_subproblem_consistency(): + dim = 2 + p0 = 2 + pI = np.array([3]) + pE = np.array([], dtype=int) + assert_tol = 1e-5 + subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol) + subproblem2 = OSQPSubproblemSQPGS(dim, p0, pI, pE, assert_tol) + D_f = np.array([[-2.0, 1.0], [-2.04236205, -1.0], [-1.92172864, -1.0]]) + D_gI = [np.array([[0.0, 2.0], [0.0, 2.0], [1.41421356, 0.0], [1.41421356, 0.0]])] + + subproblem.solve( + L=np.eye(2, dtype=float), + rho=0.1, + D_f=D_f, + D_gI=D_gI, + D_gE=[], + f_k=1.0, + gI_k=np.array([-1.0]), + gE_k=np.array([], dtype=float), + ) + + subproblem2.solve( + H=np.eye(2, dtype=float), + rho=0.1, + D_f=D_f, + D_gI=D_gI, + D_gE=[], + f_k=1.0, + gI_k=np.array([-1.0]), + gE_k=np.array([], dtype=float), + ) + + assert np.allclose(subproblem.d, subproblem2.d) + assert np.allclose(subproblem.lambda_f, subproblem2.lambda_f) + assert np.allclose(subproblem.lambda_gI, subproblem2.lambda_gI) + assert np.allclose(subproblem.lambda_gE, subproblem2.lambda_gE) From ccbd52b42390704852ce64c8376eb432ef4fd3d4 Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Mon, 22 Apr 2024 13:58:28 +0200 Subject: [PATCH 4/7] add tikohonov, and subproblem test --- src/ncopt/sqpgs/cvxpy_subproblem.py | 2 +- src/ncopt/sqpgs/main.py | 14 +++++- src/ncopt/sqpgs/osqp_subproblem.py | 11 ++++- tests/test_subproblem.py | 68 +++++++++++++++++++++++++++-- 4 files changed, 87 insertions(+), 8 deletions(-) diff --git a/src/ncopt/sqpgs/cvxpy_subproblem.py b/src/ncopt/sqpgs/cvxpy_subproblem.py index 3d0f6a1..80f82c8 100644 --- a/src/ncopt/sqpgs/cvxpy_subproblem.py +++ b/src/ncopt/sqpgs/cvxpy_subproblem.py @@ -5,7 +5,7 @@ from ncopt.sqpgs.defaults import DEFAULT_OPTION -CVXPY_SOLVER_DICT = {"osqp": cp.OSQP, "clarabel": cp.CLARABEL} +CVXPY_SOLVER_DICT = {"osqp-cvxpy": cp.OSQP, "clarabel": cp.CLARABEL} class CVXPYSubproblemSQPGS: diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py index a5a3e37..f67d663 100644 --- a/src/ncopt/sqpgs/main.py +++ b/src/ncopt/sqpgs/main.py @@ -238,11 +238,21 @@ def solve(self): ############################################## # Subproblem solve ############################################## + _reg_H = 1e-04 t1 = time.perf_counter() if isinstance(self.SP, OSQPSubproblemSQPGS): - self.SP.solve(H, rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + self.SP.solve(H + _reg_H * np.eye(self.dim), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) else: - self.SP.solve(np.linalg.cholesky(H), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) + self.SP.solve( + np.linalg.cholesky(H + _reg_H * np.eye(self.dim)), + rho, + D_f, + D_gI, + D_gE, + f_k, + gI_k, + gE_k, + ) d_k = self.SP.d.copy() t2 = time.perf_counter() diff --git a/src/ncopt/sqpgs/osqp_subproblem.py b/src/ncopt/sqpgs/osqp_subproblem.py index f614082..7f6254d 100644 --- a/src/ncopt/sqpgs/osqp_subproblem.py +++ b/src/ncopt/sqpgs/osqp_subproblem.py @@ -4,6 +4,14 @@ import osqp from scipy import sparse +# see: https://osqp.org/docs/interfaces/status_values.html +OSQP_ALLOWED_STATUS = [ + "solved", + "solved inaccurate", + "maximum iterations reached", + "run time limit reached", +] + class OSQPSubproblemSQPGS: def __init__( @@ -255,8 +263,7 @@ def solve( res = problem.solve() - assert res.info.status not in ["unsolved"] - print(res.info.status) + assert res.info.status in OSQP_ALLOWED_STATUS, f"OSQP results in status {res.info.status}." self._problem = problem primal_solution = res.x diff --git a/tests/test_subproblem.py b/tests/test_subproblem.py index 5799f7d..de62890 100644 --- a/tests/test_subproblem.py +++ b/tests/test_subproblem.py @@ -64,7 +64,7 @@ def test_subproblem_eq(subproblem_eq: CVXPYSubproblemSQPGS): assert np.isclose(subproblem_eq.objective_val, 1.0) -def test_subproblem_consistency(): +def test_subproblem_consistency_ineq(): dim = 2 p0 = 2 pI = np.array([3]) @@ -76,7 +76,7 @@ def test_subproblem_consistency(): D_gI = [np.array([[0.0, 2.0], [0.0, 2.0], [1.41421356, 0.0], [1.41421356, 0.0]])] subproblem.solve( - L=np.eye(2, dtype=float), + L=np.eye(dim, dtype=float), rho=0.1, D_f=D_f, D_gI=D_gI, @@ -87,7 +87,7 @@ def test_subproblem_consistency(): ) subproblem2.solve( - H=np.eye(2, dtype=float), + H=np.eye(dim, dtype=float), rho=0.1, D_f=D_f, D_gI=D_gI, @@ -101,3 +101,65 @@ def test_subproblem_consistency(): assert np.allclose(subproblem.lambda_f, subproblem2.lambda_f) assert np.allclose(subproblem.lambda_gI, subproblem2.lambda_gI) assert np.allclose(subproblem.lambda_gE, subproblem2.lambda_gE) + + +def test_subproblem_consistency_ineq_eq(): + dim = 4 + p0 = 2 + pI = np.array([1]) + pE = np.array([1]) + assert_tol = 1e-5 + subproblem = CVXPYSubproblemSQPGS(dim, p0, pI, pE, assert_tol) + subproblem2 = OSQPSubproblemSQPGS(dim, p0, pI, pE, assert_tol) + D_f = np.array( + [ + [1.83529234, 2.51893663, -0.87507966, 0.53305111], + [0.7042857, -0.19426588, 1.26820232, 0.59255224], + [-0.87356341, -0.24994689, -0.82073493, -1.18734854], + ] + ) + + D_gI = [ + np.array( + [ + [-1.13249659, 0.67854141, -0.0316317, -1.37963152], + [-1.64858759, 0.65296873, -0.72343526, 0.60976315], + ] + ) + ] + + D_gE = [ + np.array( + [ + [1.05532136, -0.12589961, 0.49469938, 0.22879848], + [2.10668334, -0.00816628, -0.43333072, 0.22656999], + ] + ) + ] + + subproblem.solve( + L=np.eye(dim, dtype=float), + rho=0.1, + D_f=D_f, + D_gI=D_gI, + D_gE=D_gE, + f_k=1.0, + gI_k=np.array([-1.0]), + gE_k=np.array([-1.0]), + ) + + subproblem2.solve( + H=np.eye(dim, dtype=float), + rho=0.1, + D_f=D_f, + D_gI=D_gI, + D_gE=D_gE, + f_k=1.0, + gI_k=np.array([-1.0]), + gE_k=np.array([-1.0]), + ) + + assert np.allclose(subproblem.d, subproblem2.d) + assert np.allclose(subproblem.lambda_f, subproblem2.lambda_f) + assert np.allclose(subproblem.lambda_gI, subproblem2.lambda_gI) + assert np.allclose(subproblem.lambda_gE, subproblem2.lambda_gE) From 39cf7dae6f47ddbab67e15308bfd99763e67d35f Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Mon, 22 Apr 2024 16:24:49 +0200 Subject: [PATCH 5/7] handle device properly --- src/ncopt/functions/main.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/ncopt/functions/main.py b/src/ncopt/functions/main.py index 6bb32f3..3872be5 100644 --- a/src/ncopt/functions/main.py +++ b/src/ncopt/functions/main.py @@ -28,12 +28,17 @@ def __init__( self.prepare_inputs = prepare_inputs self.is_differentiable = is_differentiable - # If no device is provided, set it to the same as the first model parameter - # this might fail for distributed models + # If no device is provided, set it to the same as the model parameters # if model has no parameters, we set device to cpu if not self.device: if sum(p.numel() for p in model.parameters() if p.requires_grad) > 0: - self.device = next(model.parameters()).device + devices = set([p.device for p in model.parameters()]) + if len(devices) == 1: + self.device = devices.pop() + else: + raise KeyError( + "Model parameters lie on more than one device. Currently not supported." + ) else: self.device = torch.device("cpu") From d9f37faead7b1d9876390d34748fe85d0dec0be8 Mon Sep 17 00:00:00 2001 From: fabian-sp Date: Tue, 23 Apr 2024 15:54:23 +0200 Subject: [PATCH 6/7] tikhonov regularization as option --- src/ncopt/sqpgs/defaults.py | 1 + src/ncopt/sqpgs/main.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/ncopt/sqpgs/defaults.py b/src/ncopt/sqpgs/defaults.py index 5cbedfd..377c77c 100644 --- a/src/ncopt/sqpgs/defaults.py +++ b/src/ncopt/sqpgs/defaults.py @@ -33,6 +33,7 @@ class Dotdict(dict): "num_points_gI": 3, "num_points_gE": 4, "qp_solver": "osqp", + "reg_H": 1e-03, } DEFAULT_ARG = Dotdict(_defaults) diff --git a/src/ncopt/sqpgs/main.py b/src/ncopt/sqpgs/main.py index f67d663..c29f6ee 100644 --- a/src/ncopt/sqpgs/main.py +++ b/src/ncopt/sqpgs/main.py @@ -238,7 +238,7 @@ def solve(self): ############################################## # Subproblem solve ############################################## - _reg_H = 1e-04 + _reg_H = self.options["reg_H"] t1 = time.perf_counter() if isinstance(self.SP, OSQPSubproblemSQPGS): self.SP.solve(H + _reg_H * np.eye(self.dim), rho, D_f, D_gI, D_gE, f_k, gI_k, gE_k) From 253ee9dda54ab6bb81e19d6be022bbac5b43aef3 Mon Sep 17 00:00:00 2001 From: fabian-sp <51907567+fabian-sp@users.noreply.github.com> Date: Wed, 24 Apr 2024 21:28:05 +0200 Subject: [PATCH 7/7] Update src/ncopt/functions/main.py Co-authored-by: Philipp Schiele <44360364+phschiele@users.noreply.github.com> --- src/ncopt/functions/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ncopt/functions/main.py b/src/ncopt/functions/main.py index 3872be5..9916d00 100644 --- a/src/ncopt/functions/main.py +++ b/src/ncopt/functions/main.py @@ -32,7 +32,7 @@ def __init__( # if model has no parameters, we set device to cpu if not self.device: if sum(p.numel() for p in model.parameters() if p.requires_grad) > 0: - devices = set([p.device for p in model.parameters()]) + devices = {p.device for p in model.parameters()} if len(devices) == 1: self.device = devices.pop() else: