diff --git a/cvxpy/cvxcore/python/canonInterface.py b/cvxpy/cvxcore/python/canonInterface.py index ebf2945e6b..ae4bbfae0a 100644 --- a/cvxpy/cvxcore/python/canonInterface.py +++ b/cvxpy/cvxcore/python/canonInterface.py @@ -292,6 +292,13 @@ def get_problem_matrix(linOps, default_canon_backend = get_default_canon_backend() canon_backend = default_canon_backend if not canon_backend else canon_backend + # DIFFENGINE is a reduction-layer replacement for ConeMatrixStuffing and + # has no lin_ops backend of its own. When code paths reach this matrix + # builder directly (e.g. tests that bypass the stuffing reduction), fall + # through to the CPP backend so the lin_ops pipeline still works. + if canon_backend == s.DIFFENGINE_BACKEND: + canon_backend = s.CPP_CANON_BACKEND + if canon_backend == s.CPP_CANON_BACKEND: from cvxpy.cvxcore.python.cppbackend import build_matrix return build_matrix(id_to_col, param_to_size, param_to_col, var_length, constr_length, linOps) diff --git a/cvxpy/reductions/dcp2cone/canonicalizers/perspective_canon.py b/cvxpy/reductions/dcp2cone/canonicalizers/perspective_canon.py index 973b84f24f..8be83140de 100644 --- a/cvxpy/reductions/dcp2cone/canonicalizers/perspective_canon.py +++ b/cvxpy/reductions/dcp2cone/canonicalizers/perspective_canon.py @@ -39,10 +39,22 @@ def perspective_canon(expr, args, solver_context: SolverInfo | None = None): prob_canon = chain.apply(aux_prob)[0] # grab problem instance # get cone representation of c, A, and b for some problem. - q = prob_canon.q.toarray().flatten()[:-1] - d = prob_canon.q.toarray().flatten()[-1] - Ab = prob_canon.A.toarray().reshape((-1, len(q) + 1), order="F") - A, b = Ab[:, :-1], Ab[:, -1] + # Extract cone representation q, d, A, b. + # ParamConeProg stores q as sparse (n+1, n_params+1) tensor with d + # embedded in the last row, and A as a flattened tensor with b embedded. + # DiffengineConeProgram stores q, d, A, b as separate concrete arrays. + if hasattr(prob_canon, 'd') and not hasattr(prob_canon.q, 'toarray'): + # DiffengineConeProgram: concrete matrices, d stored separately. + q = prob_canon.q + d = prob_canon.d + A = prob_canon.A.toarray() if hasattr(prob_canon.A, 'toarray') else prob_canon.A + b = prob_canon.b + else: + # ParamConeProg: sparse tensor with embedded offsets. + q = prob_canon.q.toarray().flatten()[:-1] + d = prob_canon.q.toarray().flatten()[-1] + Ab = prob_canon.A.toarray().reshape((-1, len(q) + 1), order="F") + A, b = Ab[:, :-1], Ab[:, -1] # given f in epigraph form, aka epi f = \{(x,t) | f(x) \leq t\} # = \{(x,t) | Fx +tg + e \in K} for K a cone, the epigraph of the diff --git a/cvxpy/reductions/dcp2cone/cone_matrix_stuffing.py b/cvxpy/reductions/dcp2cone/cone_matrix_stuffing.py index c6ade2ec8a..3411b97750 100644 --- a/cvxpy/reductions/dcp2cone/cone_matrix_stuffing.py +++ b/cvxpy/reductions/dcp2cone/cone_matrix_stuffing.py @@ -324,6 +324,53 @@ def split_solution(self, sltn, active_vars=None): value, var.shape, order='F') return sltn_dict + def apply_restruct_mat(self, restruct_mat, restruct_mat_op): + """Apply restructuring matrix to parametric tensor A. + + Parameters + ---------- + restruct_mat : list + List of sparse matrices or linear operators (unused for this path). + restruct_mat_op : LinearOperator + Block diagonal linear operator for restructuring. + + Returns + ------- + ParamConeProg + New program with restructured A tensor. + """ + if restruct_mat_op is not None: + unspecified, _ = np.divmod(self.A.shape[0] * self.A.shape[1], + restruct_mat_op.shape[1], dtype=np.int64) + reshaped_A = self.A.reshape(restruct_mat_op.shape[1], + unspecified, order='F').tocsr() + restructured_A = restruct_mat_op(reshaped_A).tocoo() + # Because of a bug in scipy versions < 1.20, `reshape` + # can overflow if indices are int32s. + restructured_A.row = restructured_A.row.astype(np.int64) + restructured_A.col = restructured_A.col.astype(np.int64) + restructured_A = restructured_A.reshape( + np.int64(restruct_mat_op.shape[0]) * (np.int64(self.x.size) + 1), + self.A.shape[1], order='F') + else: + restructured_A = self.A + return ParamConeProg( + self.q, + self.x, + restructured_A, + self.variables, + self.var_id_to_col, + self.constraints, + self.parameters, + self.param_id_to_col, + P=self.P, + formatted=True, + lower_bounds=self.lower_bounds, + upper_bounds=self.upper_bounds, + lb_tensor=self.lb_tensor, + ub_tensor=self.ub_tensor, + ) + def split_adjoint(self, del_vars=None): """Adjoint of split_solution. """ @@ -343,6 +390,65 @@ def split_adjoint(self, del_vars=None): return var_vec +def lower_and_order_constraints(constraints): + """Lower equality/inequality constraints and reorder by cone type. + + Converts Equality -> Zero, Inequality -> NonNeg, and normalizes + SOC/PowCone/ExpCone axes. Returns constraints ordered as: + Zero, NonNeg, SOC, PSD, ExpCone, PowCone3D, PowConeND. + + Parameters + ---------- + constraints : list + The problem constraints to lower and reorder. + + Returns + ------- + ordered_cons : list + The lowered and reordered constraints. + cons_id_map : dict + Mapping from constraint id to constraint id (identity map). + """ + cons = [] + for con in constraints: + if isinstance(con, Equality): + con = lower_equality(con) + elif isinstance(con, Inequality): + con = lower_ineq_to_nonneg(con) + elif isinstance(con, SOC) and con.axis == 1: + con = SOC(con.args[0], con.args[1].T, axis=0, + constr_id=con.constr_id) + elif isinstance(con, PowCone3D) and con.args[0].ndim > 1: + x, y, z = con.args + alpha = con.alpha + con = PowCone3D(x.flatten(order='F'), + y.flatten(order='F'), + z.flatten(order='F'), + alpha.flatten(order='F'), + constr_id=con.constr_id) + elif isinstance(con, PowConeND) and con.axis == 1: + alpha = con.alpha.T + W = con.W.T + con = PowConeND(W, con.z.flatten(order='F'), + alpha, + axis=0, + constr_id=con.constr_id) + elif isinstance(con, ExpCone) and con.args[0].ndim > 1: + x, y, z = con.args + con = ExpCone(x.flatten(order='F'), y.flatten(order='F'), + z.flatten(order='F'), + constr_id=con.constr_id) + cons.append(con) + + constr_map = group_constraints(cons) + ordered_cons = constr_map[Zero] + constr_map[NonNeg] + \ + constr_map[SOC] + constr_map.get(PSD, []) + \ + constr_map.get(SvecPSD, []) + constr_map[ExpCone] + \ + constr_map[PowCone3D] + constr_map[PowConeND] + cons_id_map = {con.id: con.id for con in ordered_cons} + return ordered_cons, cons_id_map + + class ConeMatrixStuffing(MatrixStuffing): """Construct matrices for linear cone problems. @@ -384,52 +490,19 @@ def stuffed_objective(self, problem, extractor): def apply(self, problem): inverse_data = InverseData(problem) - # Lower equality and inequality to Zero and NonNeg. - cons = [] - for con in problem.constraints: - if isinstance(con, Equality): - con = lower_equality(con) - elif isinstance(con, Inequality): - con = lower_ineq_to_nonneg(con) - elif isinstance(con, SOC) and con.axis == 1: - con = SOC(con.args[0], con.args[1].T, axis=0, - constr_id=con.constr_id) - elif isinstance(con, PowCone3D) and con.args[0].ndim > 1: - x, y, z = con.args - alpha = con.alpha - con = PowCone3D(x.flatten(order='F'), - y.flatten(order='F'), - z.flatten(order='F'), - alpha.flatten(order='F'), - constr_id=con.constr_id) - elif isinstance(con, PowConeND) and con.axis == 1: - alpha = con.alpha.T - W = con.W.T - con = PowConeND(W, con.z.flatten(order='F'), - alpha, - axis=0, - constr_id=con.constr_id) - elif isinstance(con, ExpCone) and con.args[0].ndim > 1: - x, y, z = con.args - con = ExpCone(x.flatten(order='F'), y.flatten(order='F'), z.flatten(order='F'), - constr_id=con.constr_id) - cons.append(con) + # Lower equality and inequality constraints and reorder by cone type. + ordered_cons, cons_id_map = lower_and_order_constraints(problem.constraints) + inverse_data.cons_id_map = cons_id_map + inverse_data.constraints = ordered_cons + # Need to check that intended canonicalization backend still works. - lowered_con_problem = problem.copy([problem.objective, cons]) + lowered_con_problem = problem.copy( + [problem.objective, ordered_cons]) canon_backend = get_canon_backend(lowered_con_problem, self.canon_backend) # Form the constraints extractor = CoeffExtractor(inverse_data, canon_backend) params_to_P, params_to_c, flattened_variable = self.stuffed_objective( problem, extractor) - # Reorder constraints to Zero, NonNeg, SOC, PSD, EXP, PowCone3D, PowConeND - constr_map = group_constraints(cons) - ordered_cons = constr_map[Zero] + constr_map[NonNeg] + \ - constr_map[SOC] + constr_map.get(PSD, []) + \ - constr_map.get(SvecPSD, []) + constr_map[ExpCone] + \ - constr_map[PowCone3D] + constr_map[PowConeND] - inverse_data.cons_id_map = {con.id: con.id for con in ordered_cons} - - inverse_data.constraints = ordered_cons # Batch expressions together, then split apart. expr_list = [arg for c in ordered_cons for arg in c.args] params_to_problem_data = extractor.affine(expr_list) diff --git a/cvxpy/reductions/dcp2cone/diffengine_cone_program.py b/cvxpy/reductions/dcp2cone/diffengine_cone_program.py new file mode 100644 index 0000000000..46b1f7b282 --- /dev/null +++ b/cvxpy/reductions/dcp2cone/diffengine_cone_program.py @@ -0,0 +1,373 @@ +""" +Copyright, the CVXPY authors + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" +from __future__ import annotations + +import numpy as np +import scipy.sparse as sp +from sparsediffpy import _sparsediffengine as _diffengine + +from cvxpy.expressions.variable import Variable +from cvxpy.reductions.dcp2cone.cone_matrix_stuffing import ConeDims, ParamConeProg +from cvxpy.reductions.matrix_stuffing import ( + extract_lower_bounds, + extract_mip_idx, + extract_upper_bounds, +) +from cvxpy.reductions.solvers.nlp_solvers.diff_engine.converters import convert_expr +from cvxpy.reductions.solvers.nlp_solvers.diff_engine.helpers import ( + build_var_dict, + normalize_shape, + to_dense_float, +) +from cvxpy.reductions.utilities import group_constraints + + +def build_capsule(objective_expr, constraint_exprs, inverse_data, params=None, verbose=False): + """Build a C diff engine problem capsule from CVXPY expressions. + + Parameters + ---------- + objective_expr : Expression + The objective expression to convert. + constraint_exprs : list[Expression] + Flat list of constraint expressions to convert. + inverse_data : InverseData + Variable/parameter offset information. + params : list[Parameter] or None + CVXPY Parameter objects, for registration and initial value setting. + verbose : bool + Whether the C engine should print output. + + Returns + ------- + capsule : PyCapsule + The C diff engine problem capsule. + n_vars : int + Total number of scalar decision variables. + param_dict : dict + Mapping {param_id: C parameter capsule}. + """ + var_dict, n_vars = build_var_dict(inverse_data) + + # Build parameter dict inline (separate from NLP path's build_param_dict). + param_dict = {} + if params: + for param_id, offset in inverse_data.param_id_map.items(): + if param_id not in inverse_data.param_shapes: + continue + d1, d2 = normalize_shape(inverse_data.param_shapes[param_id]) + param = next(p for p in params if p.id == param_id) + p = to_dense_float(param.value) + param_dict[param_id] = _diffengine.make_parameter( + d1, d2, offset, n_vars, p.flatten(order='F')) + + c_obj = convert_expr(objective_expr, var_dict, n_vars, param_dict) + c_constraints = [convert_expr(e, var_dict, n_vars, param_dict) for e in constraint_exprs] + + capsule = _diffengine.make_problem(c_obj, c_constraints, verbose) + + if param_dict and params: + _diffengine.problem_register_params(capsule, list(param_dict.values())) + theta = np.concatenate([ + np.asarray(p.value, dtype=np.float64).flatten(order='F') + for p in params + ]) + _diffengine.problem_update_params(capsule, theta) + + return capsule, n_vars, param_dict + + +class DiffengineConeProgram(ParamConeProg): + """A cone program with matrices extracted via the diffengine. + + Duck-type compatible with ParamConeProg. On first solve, stores the sparsity + pattern of A and P. When parameters are present, re-evaluates the + converted C expression trees on subsequent solves via apply_parameters(). + + minimize q'x + d + [(1/2)x'Px] + subject to cone_constr(A*x + b) in cones + """ + + def __init__( + self, + x: Variable, + A: sp.spmatrix, + b: np.ndarray, + q: np.ndarray, + d: float, + P, + constraints: list, + inverse_data, + formatted: bool = False, + lower_bounds=None, + upper_bounds=None, + capsule=None, + parameters=None, + quad_obj: bool = False, + ) -> None: + self.A = A + self.b = b + self.q = q + self.d = d + self.x = x + self.P = P + + self.constraints = constraints + self.constr_map = group_constraints(constraints) + self.cone_dims = ConeDims(self.constr_map) + + self.inverse_data = inverse_data + self.formatted = formatted + self.lower_bounds = lower_bounds + self.upper_bounds = upper_bounds + + self._capsule = capsule + self.quad_obj = quad_obj + self._restruct_mat = None + self.parameters = list(parameters) if parameters else [] + self.param_id_to_size = {p.id: p.size for p in self.parameters} + + # TODO what should we do about parametric bounds? + # Duck-type compatibility (always None for diffengine path). + self.lb_tensor = None + self.ub_tensor = None + + @property + def variables(self): + return [self.inverse_data.id2var[vid] + for vid in self.inverse_data.var_offsets] + + @property + def var_id_to_col(self): + return self.inverse_data.var_offsets + + @property + def id_to_var(self): + return self.inverse_data.id2var + + @property + def param_id_to_col(self): + return self.inverse_data.param_id_map + + def is_mixed_integer(self) -> bool: + return self.x.attributes['boolean'] or self.x.attributes['integer'] + + def apply_parameters(self, id_to_param_value=None, zero_offset: bool = False, + keep_zeros: bool = False, quad_obj: bool = False): + """Return problem matrices, re-evaluating if parameters are present. + + When no parameters exist, returns the stored matrices directly. + When parameters exist, updates the C DAG with current parameter values + and re-evaluates objective/constraints at x=0. + """ + if not self.parameters or self._capsule is None: + # Non-parametric: return stored matrices directly. + A = self.A + if quad_obj and self.P is not None: + return self.P, self.q, self.d, A, self.b + return self.q, self.d, A, self.b + + # Build theta vector from current parameter values. + if id_to_param_value is not None: + parts = [np.asarray(id_to_param_value[p.id], + dtype=np.float64).flatten(order='F') + for p in self.parameters] + else: + parts = [np.asarray(p.value, + dtype=np.float64).flatten(order='F') + for p in self.parameters] + theta = np.concatenate(parts) + + _diffengine.problem_update_params(self._capsule, theta) + + # Re-evaluate at x0 = 0. + n_vars = self.inverse_data.x_length + x0 = np.zeros(n_vars, dtype=np.float64) + + d = float(_diffengine.problem_objective_forward(self._capsule, x0)) + q = _diffengine.problem_gradient(self._capsule).copy() + + if self.constraints: + b_vec = _diffengine.problem_constraint_forward(self._capsule, x0) + jac_data, jac_indices, jac_indptr, jac_shape = \ + _diffengine.problem_jacobian(self._capsule) + # TODO: have the diffengine return CSC directly to avoid this conversion. + A = sp.csr_matrix( + (jac_data, jac_indices, jac_indptr), + shape=(jac_shape[0], n_vars)).tocsc() + else: + b_vec = np.array([], dtype=np.float64) + A = sp.csc_matrix((0, n_vars)) + + b = np.atleast_1d(b_vec) + + # Apply cached restructuring matrix if present. + if self._restruct_mat is not None and self._restruct_mat is not False: + A = self._restruct_mat @ A + b = np.asarray(self._restruct_mat @ b).flatten() + + # Update stored matrices. + self.A, self.b, self.q, self.d = A, b, q, d + + if quad_obj: + duals = np.zeros(b.shape[0], dtype=np.float64) + h_data, h_indices, h_indptr, h_shape = \ + _diffengine.problem_hessian(self._capsule, 1.0, duals) + P_csr = sp.csr_matrix((h_data, h_indices, h_indptr), shape=h_shape) + self.P = P_csr.tocsc() + return self.P, q, d, A, b + return q, d, A, b + + def apply_restruct_mat(self, restruct_mat, restruct_mat_op=None): + """Apply restructuring matrix to concrete A, b matrices. + + Materializes the block-diagonal restructuring matrix R and caches it + so that parametric re-evaluations can re-apply it. + + Parameters + ---------- + restruct_mat : list + List of sparse matrices or linear operators forming a block diagonal. + restruct_mat_op : LinearOperator or None + Unused for DiffengineConeProgram (uses restruct_mat directly). + + Returns + ------- + DiffengineConeProgram + New program with restructured A and b. + """ + R = None + if restruct_mat: + sparse_mats = [] + for mat in restruct_mat: + if sp.issparse(mat): + sparse_mats.append(sp.csc_matrix(mat)) + elif callable(mat): + eye = sp.eye_array(mat.shape[1], format='csc') + sparse_mats.append(sp.csc_matrix(mat(eye))) + else: + eye = sp.eye_array(mat.shape[1], format='csc') + sparse_mats.append(sp.csc_matrix(mat @ eye)) + + R = sp.block_diag(sparse_mats, format='csc') + new_A = R @ self.A + new_b = np.asarray(R @ self.b).flatten() + else: + new_A, new_b = self.A, self.b + + new_prog = DiffengineConeProgram( + self.x, new_A, new_b, self.q, self.d, self.P, + self.constraints, self.inverse_data, + formatted=True, + lower_bounds=self.lower_bounds, + upper_bounds=self.upper_bounds, + capsule=self._capsule, + parameters=self.parameters, + quad_obj=self.quad_obj, + ) + if R is not None: + new_prog._restruct_mat = R + return new_prog + + def split_solution(self, sltn, active_vars=None): + """ + Splits the solution into individual variables. + """ + from cvxpy.reductions import cvx_attr2constr + if active_vars is None: + active_vars = [v.id for v in self.variables] + sltn_dict = {} + for var_id, col in self.var_id_to_col.items(): + if var_id in active_vars: + var = self.id_to_var[var_id] + value = sltn[col:var.size + col] + if var.attributes_were_lowered(): + orig_var = var.leaf_of_provenance() + value = cvx_attr2constr.recover_value_for_leaf( + orig_var, value, project=False) + sltn_dict[orig_var.id] = np.reshape( + value, orig_var.shape, order='F') + else: + sltn_dict[var_id] = np.reshape( + value, var.shape, order='F') + return sltn_dict + + @classmethod + def from_problem(cls, problem, ordered_cons, inverse_data, quad_obj): + """Build a DiffengineConeProgram by evaluating expressions at x=0. + + Uses the sparsediffpy diffengine to compute A, b, q, d, P matrices. + + Parameters + ---------- + problem : Problem + The CVXPY problem (post Dcp2Cone, with affine constraints and + linear/quadratic objective). + ordered_cons : list + Ordered constraints (Zero, NonNeg, SOC, PSD, ExpCone, ...). + inverse_data : InverseData + Inverse data for the problem. + quad_obj : bool + Whether the objective is quadratic. + """ + # Build the diff engine problem capsule. + expr_list = [arg for c in ordered_cons for arg in c.args] + params = problem.parameters() + capsule, n_vars, _ = build_capsule( + problem.objective.expr, expr_list, inverse_data, params=params, verbose=True) + + # Create flattened variable with MIP info. + boolean, integer = extract_mip_idx(problem.variables()) + x = Variable(n_vars, boolean=boolean, integer=integer) + + # Evaluate at x0 = 0. + x0 = np.zeros(n_vars, dtype=np.float64) + + if quad_obj: + _diffengine.problem_init_derivatives(capsule) + else: + _diffengine.problem_init_jacobian(capsule) + + d = float(_diffengine.problem_objective_forward(capsule, x0)) + q = _diffengine.problem_gradient(capsule).copy() + + if expr_list: + b_vec = _diffengine.problem_constraint_forward(capsule, x0) + jac_data, jac_indices, jac_indptr, jac_shape = \ + _diffengine.problem_jacobian(capsule) + # TODO: have the diffengine return CSC directly to avoid this conversion. + A = sp.csr_matrix( + (jac_data, jac_indices, jac_indptr), + shape=(jac_shape[0], n_vars)).tocsc() + else: + b_vec = np.array([], dtype=np.float64) + A = sp.csc_matrix((0, n_vars)) + + P = None + if quad_obj: + duals = np.zeros(b_vec.shape[0], dtype=np.float64) + h_data, h_indices, h_indptr, h_shape = \ + _diffengine.problem_hessian(capsule, 1.0, duals) + P_csr = sp.csr_matrix((h_data, h_indices, h_indptr), shape=h_shape) + P = P_csr.tocsc() + + lower_bounds = extract_lower_bounds(problem.variables(), n_vars) + upper_bounds = extract_upper_bounds(problem.variables(), n_vars) + + return cls(x, A, np.atleast_1d(b_vec), q, d, P, + ordered_cons, inverse_data, + lower_bounds=lower_bounds, upper_bounds=upper_bounds, + capsule=capsule, parameters=params, quad_obj=quad_obj) diff --git a/cvxpy/reductions/dcp2cone/diffengine_matrix_stuffing.py b/cvxpy/reductions/dcp2cone/diffengine_matrix_stuffing.py new file mode 100644 index 0000000000..2fe4ae16cd --- /dev/null +++ b/cvxpy/reductions/dcp2cone/diffengine_matrix_stuffing.py @@ -0,0 +1,124 @@ +""" +Copyright, the CVXPY authors + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" +from __future__ import annotations + +import numpy as np + +import cvxpy.settings as s +from cvxpy.constraints import SOC, ExpCone +from cvxpy.problems.objective import Minimize +from cvxpy.reductions import InverseData, Solution +from cvxpy.reductions.dcp2cone.cone_matrix_stuffing import lower_and_order_constraints +from cvxpy.reductions.matrix_stuffing import MatrixStuffing + + +class DiffengineMatrixStuffing(MatrixStuffing): + """Construct matrices for cone problems using the diffengine (C autodiff) backend. + + Sibling of ConeMatrixStuffing that uses sparsediffpy's C autodiff engine + to extract A, b, q, d, P matrices directly, bypassing the parametric + tensor pipeline entirely. + """ + CONSTRAINTS = 'ordered_constraints' + + def __init__(self, quad_obj: bool = False): + self.quad_obj = quad_obj + + def accepts(self, problem): + from cvxpy.reductions import cvx_attr2constr + from cvxpy.reductions.utilities import are_args_affine + + valid_obj_curv = (self.quad_obj and problem.objective.expr.is_quadratic()) or \ + problem.objective.expr.is_affine() + return (type(problem.objective) == Minimize + and valid_obj_curv + and not cvx_attr2constr.convex_attributes(problem.variables()) + and are_args_affine(problem.constraints) + and problem.is_dpp()) + + def apply(self, problem): + from cvxpy.reductions.dcp2cone.diffengine_cone_program import ( + DiffengineConeProgram, + ) + + inverse_data = InverseData(problem) + + ordered_cons, cons_id_map = lower_and_order_constraints(problem.constraints) + + inverse_data.cons_id_map = cons_id_map + inverse_data.constraints = ordered_cons + inverse_data.minimize = type(problem.objective) == Minimize + + new_prob = DiffengineConeProgram.from_problem( + problem, ordered_cons, inverse_data, self.quad_obj) + + return new_prob, inverse_data + + def invert(self, solution, inverse_data): + """Retrieves a solution to the original problem.""" + var_map = inverse_data.var_offsets + con_map = inverse_data.cons_id_map + # Flip sign of opt val if maximize. + opt_val = solution.opt_val + if solution.status not in s.ERROR and not inverse_data.minimize: + opt_val = -solution.opt_val + + # Remap dual variables before SOLUTION_PRESENT check so duals + # are available for infeasible/unbounded certificates. + dual_vars = {} + if solution.dual_vars is not None: + for old_con, new_con in con_map.items(): + con_obj = inverse_data.id2cons[old_con] + shape = con_obj.shape + dual_value = solution.dual_vars.get(new_con) + # TODO rationalize Exponential. + if dual_value is not None: + if shape == () or isinstance(con_obj, (ExpCone, SOC)): + dual_vars[old_con] = dual_value + else: + dual_vars[old_con] = np.reshape(dual_value, shape, + order='F') + + primal_vars = {} + if solution.status not in s.SOLUTION_PRESENT: + return Solution(solution.status, opt_val, primal_vars, dual_vars, + solution.attr) + + # Split vectorized variable into components. + x_opt = list(solution.primal_vars.values())[0] + for var_id, offset in var_map.items(): + shape = inverse_data.var_shapes[var_id] + size = np.prod(shape, dtype=int) + primal_vars[var_id] = np.reshape(x_opt[offset:offset+size], shape, + order='F') + + # Remap duals again for optimal solutions (overwrite with full mapping). + if solution.dual_vars is not None: + for old_con, new_con in con_map.items(): + con_obj = inverse_data.id2cons[old_con] + shape = con_obj.shape + # TODO rationalize Exponential. + if shape == () or isinstance(con_obj, (ExpCone, SOC)): + dual_vars[old_con] = solution.dual_vars[new_con] + else: + dual_vars[old_con] = np.reshape( + solution.dual_vars[new_con], + shape, + order='F' + ) + + return Solution(solution.status, opt_val, primal_vars, dual_vars, + solution.attr) diff --git a/cvxpy/reductions/solvers/conic_solvers/conic_solver.py b/cvxpy/reductions/solvers/conic_solvers/conic_solver.py index b9842357b4..9b2073561c 100644 --- a/cvxpy/reductions/solvers/conic_solvers/conic_solver.py +++ b/cvxpy/reductions/solvers/conic_solvers/conic_solver.py @@ -271,45 +271,12 @@ def format_constraints(cls, problem, exp_cone_order): else: raise ValueError("Unsupported constraint type.") - # Form new ParamConeProg + # Apply restructuring via polymorphic method. if restruct_mat: - # TODO(akshayka): profile to see whether using linear operators - # or bmat is faster - restruct_mat = as_block_diag_linear_operator(restruct_mat) - # this is equivalent to but _much_ faster than: - # restruct_mat_rep = sp.block_diag([restruct_mat]*(problem.x.size + 1)) - # restruct_A = restruct_mat_rep * problem.A - unspecified, _ = np.divmod(problem.A.shape[0] * problem.A.shape[1], - restruct_mat.shape[1], dtype=np.int64) - reshaped_A = problem.A.reshape(restruct_mat.shape[1], - unspecified, order='F').tocsr() - restructured_A = restruct_mat(reshaped_A).tocoo() - # Because of a bug in scipy versions < 1.20, `reshape` - # can overflow if indices are int32s. - restructured_A.row = restructured_A.row.astype(np.int64) - restructured_A.col = restructured_A.col.astype(np.int64) - restructured_A = restructured_A.reshape( - np.int64(restruct_mat.shape[0]) * (np.int64(problem.x.size) + 1), - problem.A.shape[1], order='F') + restruct_mat_op = as_block_diag_linear_operator(restruct_mat) else: - restructured_A = problem.A - new_param_cone_prog = ParamConeProg( - problem.q, - problem.x, - restructured_A, - problem.variables, - problem.var_id_to_col, - problem.constraints, - problem.parameters, - problem.param_id_to_col, - P=problem.P, - formatted=True, - lower_bounds=problem.lower_bounds, - upper_bounds=problem.upper_bounds, - lb_tensor=problem.lb_tensor, - ub_tensor=problem.ub_tensor, - ) - return new_param_cone_prog + restruct_mat_op = None + return problem.apply_restruct_mat(restruct_mat, restruct_mat_op) def invert(self, solution, inverse_data): """Returns the solution to the original problem given the inverse_data. diff --git a/cvxpy/reductions/solvers/nlp_solvers/diff_engine/converters.py b/cvxpy/reductions/solvers/nlp_solvers/diff_engine/converters.py index 1fe6c9688b..b1d5f96179 100644 --- a/cvxpy/reductions/solvers/nlp_solvers/diff_engine/converters.py +++ b/cvxpy/reductions/solvers/nlp_solvers/diff_engine/converters.py @@ -105,6 +105,15 @@ def convert_expr(expr, var_dict, n_vars, param_dict=None): d1, d2 = normalize_shape(expr.shape) return _diffengine.make_parameter(d1, d2, -1, n_vars, c.flatten(order='F')) + # Constant atom: no variables/parameters, so it must have a concrete value. + # Handles atoms like floor(NegExpression(Constant(5.0))) after EvalParams. + # Check variables()/parameters() before .value to avoid triggering + # numeric() on atoms that don't support it (e.g. cached SymbolicQuadForm). + if not expr.variables() and not expr.parameters() and expr.value is not None: + c = to_dense_float(expr.value) + d1, d2 = normalize_shape(expr.shape) + return _diffengine.make_parameter(d1, d2, -1, n_vars, c.flatten(order='F')) + # Recursive case: atoms atom_name = type(expr).__name__ children = [convert_expr(arg, var_dict, n_vars, param_dict) for arg in expr.args] @@ -115,6 +124,11 @@ def convert_expr(expr, var_dict, n_vars, param_dict=None): C_expr = convert_matmul(expr, children, var_dict, n_vars, param_dict) elif atom_name == "multiply": C_expr = convert_multiply(expr, children, var_dict, n_vars, param_dict) + elif atom_name in ("QuadForm", "SymbolicQuadForm"): + from cvxpy.reductions.solvers.nlp_solvers.diff_engine.registry import ( + convert_quad_form, + ) + C_expr = convert_quad_form(expr, children, n_vars) elif atom_name in ATOM_CONVERTERS: C_expr = ATOM_CONVERTERS[atom_name](expr, children) else: diff --git a/cvxpy/reductions/solvers/nlp_solvers/diff_engine/registry.py b/cvxpy/reductions/solvers/nlp_solvers/diff_engine/registry.py index 098eb304bb..06a64e1be4 100644 --- a/cvxpy/reductions/solvers/nlp_solvers/diff_engine/registry.py +++ b/cvxpy/reductions/solvers/nlp_solvers/diff_engine/registry.py @@ -45,6 +45,36 @@ def convert_conv(expr, children): return _diffengine.make_convolve(children[0], children[1]) +def convert_concatenate(expr, children): + """Convert Concatenate along an existing axis via hstack/vstack. + + Concatenate semantics: + ndim=1, axis=0 → flat concat, same as hstack on 1D + ndim=2, axis=0 → stack rows, same as vstack + ndim=2, axis=1 → stack cols, same as hstack + """ + axis = expr.axis + arg_shapes = [arg.shape for arg in expr.args] + ndims = {len(s) for s in arg_shapes} + if len(ndims) != 1: + raise NotImplementedError( + "Concatenate across inputs of mixed dimensionality is not " + "supported by the diffengine backend." + ) + ndim = ndims.pop() + + if ndim == 1 and axis == 0: + return _diffengine.make_hstack(children) + if ndim == 2 and axis == 0: + return _diffengine.make_vstack(children) + if ndim == 2 and axis == 1: + return _diffengine.make_hstack(children) + raise NotImplementedError( + f"Concatenate with ndim={ndim}, axis={axis} is not supported " + "by the diffengine backend." + ) + + def convert_div(expr, children): """Convert x / c by multiplying x by the elementwise reciprocal of c. @@ -104,8 +134,39 @@ def convert_rel_entr(expr, children): return _diffengine.make_rel_entr(children[0], children[1]) -def convert_quad_form(expr, children): - """Convert scalar quadratic form x.T @ P @ x.""" +def _diagonal_weights(P_val): + """If P_val is (numerically) diagonal, return its diagonal as a 1D + float64 array; otherwise return None.""" + if sparse.issparse(P_val): + P_val = P_val.toarray() + P = np.asarray(P_val, dtype=np.float64) + if P.ndim != 2 or P.shape[0] != P.shape[1]: + return None + w = np.diag(P) + if np.allclose(P, np.diag(w)): + return w + return None + + +def convert_quad_form(expr, children, n_vars): + """Convert quadratic form x.T @ P @ x. + + Scalar output (shape ()) uses the native `make_quad_form` C node. + + Vector-shaped SymbolicQuadForm is lowered using existing atoms when P + is diagonal (every CVXPY canonicalizer that produces a vector-shaped + SymbolicQuadForm today builds it with P = eye or P = α·eye): + + - `block_indices is None`, P = diag(w): + y = w ⊙ x**2 + - `block_indices` given with uniform block size m, P = α·I: + y_j = α · Σ_{i ∈ I_j} x_i**2 + (gather → square → axis-0 sum → scalar-mult by α) + + A genuinely non-diagonal vector quad_form still raises + NotImplementedError — that case needs a native block quad_form node + in SparseDiffPy. + """ P = expr.args[1] if not P.is_constant(): @@ -118,27 +179,98 @@ def convert_quad_form(expr, children): "is not supported by the diff engine." ) - if not isinstance(P_val, sparse.csr_matrix): - P_val = sparse.csr_matrix(P_val) - + if expr.size > 1: + w = _diagonal_weights(P_val) + if w is None: + raise NotImplementedError( + f"Non-diagonal vector-shaped quad form (shape {expr.shape}) " + "is not supported by the diffengine backend. Needs native " + "block quadform in SparseDiffPy." + ) + + block_indices = getattr(expr, "block_indices", None) + child = children[0] + uniform_alpha = float(w[0]) if np.allclose(w, w[0]) else None + + if block_indices is None: + # Case A: y_j = w_j * x_j**2 (elementwise, K == N). + x_sq = _diffengine.make_power(child, 2.0) + if uniform_alpha is not None and np.isclose(uniform_alpha, 1.0): + out = x_sq + elif uniform_alpha is not None: + alpha_node = _diffengine.make_parameter( + 1, 1, -1, n_vars, + np.array([uniform_alpha], dtype=np.float64)) + out = _diffengine.make_param_scalar_mult(alpha_node, x_sq) + else: + d1, d2 = normalize_shape(expr.shape) + w_node = _diffengine.make_parameter( + d1, d2, -1, n_vars, + np.asarray(w, dtype=np.float64).flatten(order="F")) + out = _diffengine.make_param_vector_mult(w_node, x_sq) + else: + # Case B: y_j = α · Σ_{i ∈ I_j} x_i**2 (requires uniform blocks). + if uniform_alpha is None: + raise NotImplementedError( + "quad_form with block_indices and non-uniform P diagonal " + "is not supported by the diffengine backend." + ) + block_lens = {len(b) for b in block_indices} + if len(block_lens) != 1: + raise NotImplementedError( + "quad_form with non-uniform block_indices sizes is not " + "supported by the diffengine backend." + ) + m = block_lens.pop() + K = len(block_indices) + flat_idx = np.concatenate( + [np.asarray(b, dtype=np.int32) for b in block_indices] + ).astype(np.int32) + gathered = _diffengine.make_index(child, m, K, flat_idx) + sq = _diffengine.make_power(gathered, 2.0) + summed = _diffengine.make_sum(sq, 0) + if np.isclose(uniform_alpha, 1.0): + out = summed + else: + alpha_node = _diffengine.make_parameter( + 1, 1, -1, n_vars, + np.array([uniform_alpha], dtype=np.float64)) + out = _diffengine.make_param_scalar_mult(alpha_node, summed) + + # Match the declared output shape (convert_expr's final dim check + # compares C dims to normalize_shape(expr.shape)). + d1_out, d2_out = normalize_shape(expr.shape) + d1_c, d2_c = _diffengine.get_expr_dimensions(out) + if (d1_c, d2_c) != (d1_out, d2_out): + out = _diffengine.make_reshape(out, d1_out, d2_out) + return out + + # Scalar output: native quad_form node. + if sparse.issparse(P_val): + P_val = P_val.toarray() + P_arr = np.asarray(P_val, dtype=np.float64) + P_csr = sparse.csr_matrix(P_arr) return _diffengine.make_quad_form( children[0], - P_val.data.astype(np.float64), - P_val.indices.astype(np.int32), - P_val.indptr.astype(np.int32), - P_val.shape[0], - P_val.shape[1], + P_csr.data.astype(np.float64), + P_csr.indices.astype(np.int32), + P_csr.indptr.astype(np.int32), + P_csr.shape[0], + P_csr.shape[1], ) def convert_reshape(expr, children): - """Convert reshape. C-order via transpose(F-reshape(transpose(x))). + """Convert reshape. - Identity: reshape(x, (m, n), C) == transpose(reshape(transpose(x), (n, m), F)). + F-order (column-major) uses make_reshape directly. + C-order (row-major) is decomposed as: transpose(reshape(transpose(x), (n, m), F)) + since reshape(x, (m, n), C) == transpose(reshape(transpose(x), (n, m), F)). """ d1, d2 = normalize_shape(expr.shape) if expr.order == "F": return _diffengine.make_reshape(children[0], d1, d2) + # C-order: transpose input, F-reshape to swapped dims, transpose output transposed = _diffengine.make_transpose(children[0]) reshaped = _diffengine.make_reshape(transposed, d2, d1) return _diffengine.make_transpose(reshaped) @@ -215,7 +347,6 @@ def convert_diag_mat(expr, children): node = _diffengine.make_diag_mat(children[0]) # C produces (n, 1) but CVXPY shape is (n,) which normalizes to (1, n) # TODO add support for producing (1, n) directly in C and remove this reshape - # TODO also raise error that the k should be zero, since that's the only supported case n = expr.args[0].shape[0] return _diffengine.make_reshape(node, 1, n) @@ -248,7 +379,8 @@ def convert_upper_tri(_expr, children): # Reductions "Sum": convert_sum, # Bivariate - "QuadForm": convert_quad_form, + # QuadForm and SymbolicQuadForm are dispatched directly in convert_expr + # (they need n_vars for creating constant nodes in the vector case). "quad_over_lin": convert_quad_over_lin, "rel_entr": convert_rel_entr, # Elementwise univariate with parameter @@ -279,6 +411,7 @@ def convert_upper_tri(_expr, children): # Horizontal/vertical stack "Hstack": convert_hstack, "Vstack": convert_vstack, + "Concatenate": convert_concatenate, # 1D full convolution "conv": convert_conv, "convolve": convert_conv, diff --git a/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py b/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py index 5bdfee649c..667c4ea8bd 100644 --- a/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py +++ b/cvxpy/reductions/solvers/nlp_solvers/nlp_solver.py @@ -245,6 +245,7 @@ def update_params(self, problem) -> None: raise ValueError("update_params called but problem has no parameters. " "This is a bug and should be reported.") + theta = np.concatenate([ np.asarray(p.value, dtype=np.float64).flatten(order='F') for p in problem.parameters() diff --git a/cvxpy/reductions/solvers/solving_chain.py b/cvxpy/reductions/solvers/solving_chain.py index dbebca54b9..21af7d3fb6 100644 --- a/cvxpy/reductions/solvers/solving_chain.py +++ b/cvxpy/reductions/solvers/solving_chain.py @@ -35,7 +35,11 @@ from cvxpy.reductions.solvers.constant_solver import ConstantSolver from cvxpy.reductions.solvers.qp_solvers.qp_solver import QpSolver from cvxpy.reductions.solvers.solver import Solver, expand_cones -from cvxpy.settings import COO_CANON_BACKEND, DPP_PARAM_THRESHOLD +from cvxpy.settings import ( + COO_CANON_BACKEND, + DPP_PARAM_THRESHOLD, + SCIPY_CANON_BACKEND, +) from cvxpy.utilities.solver_context import SolverInfo from cvxpy.utilities.warn import warn @@ -208,6 +212,10 @@ def _build_solving_chain( if total_param_size >= DPP_PARAM_THRESHOLD: canon_backend = COO_CANON_BACKEND + # Resolve None to the configured default so the DIFFENGINE check below fires. + if canon_backend is None: + canon_backend = "DIFFENGINE" + # --- Canonicalization reductions (problem_form + solver_context) --- use_quad = True if solver_opts is None else solver_opts.get('use_quad_obj', True) @@ -234,10 +242,17 @@ def _build_solving_chain( if solver_instance.SOC_DIM3_ONLY and SOC in cones: reductions.append(SOCDim3()) - reductions += [ - ConeMatrixStuffing(quad_obj=quad_obj, canon_backend=canon_backend), - solver_instance, - ] + if canon_backend == "DIFFENGINE" and problem._max_ndim() > 2: + canon_backend = SCIPY_CANON_BACKEND + + if canon_backend == "DIFFENGINE": + from cvxpy.reductions.dcp2cone.diffengine_matrix_stuffing import ( + DiffengineMatrixStuffing, + ) + stuffing = DiffengineMatrixStuffing(quad_obj=quad_obj) + else: + stuffing = ConeMatrixStuffing(quad_obj=quad_obj, canon_backend=canon_backend) + reductions += [stuffing, solver_instance] return SolvingChain(reductions=reductions, solver_context=solver_context) diff --git a/cvxpy/settings.py b/cvxpy/settings.py index 9965319f9b..8a9e8e4218 100644 --- a/cvxpy/settings.py +++ b/cvxpy/settings.py @@ -194,6 +194,7 @@ RUST_CANON_BACKEND = "RUST" CPP_CANON_BACKEND = "CPP" COO_CANON_BACKEND = "COO" # 3D COO sparse tensor backend: O(nnz) operations for large parameters +DIFFENGINE_BACKEND = "DIFFENGINE" # Default canonicalization backend, pyodide uses SciPy DEFAULT_CANON_BACKEND = CPP_CANON_BACKEND if sys.platform != "emscripten" else SCIPY_CANON_BACKEND