R/235_reductions_solvers_qp_solvers_qp_solver.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/reductions/solvers/qp_solvers/qp_solver.R
#####

## CVXPY SOURCE: reductions/solvers/qp_solvers/qp_solver.py
## QpSolver -- base class for QP-path solvers
##
## Converts ConeMatrixStuffing data (A, b, c, P, dims) into standard QP form:
##   minimize    0.5 x' P x + q' x
##   subject to  A_eq x  = b_eq
##               F_ineq x <= g_ineq
##
## The KEY difference from ConicSolver: QpSolver NEGATES the inequality rows
## of AF (the combined constraint matrix), producing F = -AF[ineq]. This sign
## flip is what makes all QP solver dual conventions match CVXPY -- eliminating
## the ad-hoc per-solver sign hacks that were needed when QP solvers inherited
## from ConicSolver.


# -- QpSolver class ----------------------------------------------
## CVXPY SOURCE: qp_solver.py lines 35-149

QpSolver <- new_class("QpSolver", parent = Solver, package = "CVXR",
  properties = list(
    SUPPORTED_CONSTRAINTS = new_property(class_list,
      default = list(Zero, NonNeg)),
    REQUIRES_CONSTR = new_property(class_logical, default = FALSE)
  ),
  constructor = function(MIP_CAPABLE = FALSE, BOUNDED_VARIABLES = FALSE) {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = MIP_CAPABLE,
      BOUNDED_VARIABLES = BOUNDED_VARIABLES,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg),
      REQUIRES_CONSTR = FALSE
    )
  }
)

# -- reduction_accepts ------------------------------------------
## QP solvers accept only Zero + NonNeg constraints.

method(reduction_accepts, QpSolver) <- function(x, problem, ...) {
  if (!is.list(problem) || is.null(problem[["constraints"]])) return(FALSE)
  constrs <- problem[["constraints"]]
  if (length(constrs) == 0L && x@REQUIRES_CONSTR) return(FALSE)
  all(vapply(constrs, function(c) {
    S7_inherits(c, Zero) || S7_inherits(c, NonNeg)
  }, logical(1L)))
}

# -- reduction_apply --------------------------------------------
## CVXPY SOURCE: qp_solver.py lines 79-149
##
## Converts ConeMatrixStuffing output into QP standard form:
##   Equality:   A_eq x = b_eq    (first len_eq rows, negate b)
##   Inequality: F_ineq x <= g_ineq (next len_ineq rows, NEGATE F)
##
## The sign flip F = -AF[ineq] is the KEY that makes all QP solver
## dual conventions match CVXPY -- no per-solver sign hacks needed.

method(reduction_apply, QpSolver) <- function(x, problem, ...) {
  data <- problem  ## data list from ConeMatrixStuffing
  inv_data <- list()

  inv_data[[SOLVER_VAR_ID]] <- data[["x_id"]]
  cone_dims <- data[[SD_DIMS]]
  inv_data[[SD_DIMS]] <- cone_dims

  ## Validate: only Zero + NonNeg (no SOC, ExpCone, PSD, etc.)
  if (length(cone_dims@soc) > 0L || cone_dims@exp > 0L ||
      length(cone_dims@psd) > 0L || length(cone_dims@p3d) > 0L) {
    cli_abort(c(
      "QP solver {.fn solver_name} cannot handle conic constraints.",
      "i" = "Problem has non-QP cones. Use a conic solver instead."
    ))
  }

  ## Split into eq/ineq
  constr_map <- group_constraints(data[["constraints"]])
  inv_data[[SOLVER_EQ_CONSTR]] <- constr_map[["Zero"]]
  inv_data[[SOLVER_NEQ_CONSTR]] <- constr_map[["NonNeg"]]

  n <- length(data[[SD_C]])
  len_eq <- cone_dims@zero
  len_ineq <- cone_dims@nonneg
  AF <- data[[SD_A]]
  bg <- data[[SD_B]]

  ## KEY SIGN FLIP -- matches CVXPY qp_solver.py lines 122-132
  ## ConeMatrixStuffing produces: A*x + b = 0 (Zero), A*x + b >= 0 (NonNeg)
  ## QP form: A_eq*x = b_eq, F_ineq*x <= g_ineq
  ## Equality: A_eq = AF[1:len_eq, ],  b_eq = -bg[1:len_eq]
  ## Inequality: F = -AF[ineq, ],  g = bg[ineq]  (NEGATE F!)
  A_eq <- if (len_eq > 0L) {
    AF[seq_len(len_eq), , drop = FALSE]
  } else {
    Matrix::sparseMatrix(i = integer(0), j = integer(0),
                         dims = c(0L, n))
  }
  b_eq <- if (len_eq > 0L) -bg[seq_len(len_eq)] else numeric(0)

  F_ineq <- if (len_ineq > 0L) {
    -AF[(len_eq + 1L):(len_eq + len_ineq), , drop = FALSE]
  } else {
    Matrix::sparseMatrix(i = integer(0), j = integer(0),
                         dims = c(0L, n))
  }
  g_ineq <- if (len_ineq > 0L) {
    bg[(len_eq + 1L):(len_eq + len_ineq)]
  } else {
    numeric(0)
  }

  ## Build solver data: P, q, A_eq, b_eq, F_ineq, g_ineq
  solver_data <- list()
  solver_data[[SD_P]] <- data[[SD_P]]
  solver_data[["q"]] <- data[[SD_C]]
  solver_data[["A_eq"]] <- A_eq
  solver_data[["b_eq"]] <- b_eq
  solver_data[["F_ineq"]] <- F_ineq
  solver_data[["g_ineq"]] <- g_ineq
  solver_data[[SD_DIMS]] <- cone_dims

  ## Pass through MIP indices
  solver_data[["bool_idx"]] <- data[["bool_idx"]]
  solver_data[["int_idx"]] <- data[["int_idx"]]

  inv_data[[SD_OFFSET]] <- data[[SD_OFFSET]]
  inv_data[["is_mip"]] <- length(data[["bool_idx"]] %||% integer(0)) > 0L ||
    length(data[["int_idx"]] %||% integer(0)) > 0L

  list(solver_data, inv_data)
}

# -- print ------------------------------------------------------

method(print, QpSolver) <- function(x, ...) {
  cat(sprintf("QpSolver(%s)\n", solver_name(x)))
  invisible(x)
}