R/287_reductions_solvers_qp_solvers_highs_qpif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/qp_solvers/highs_qpif.py
## HiGHS QP solver interface for LP/QP problems
##
## HiGHS solves: minimize 0.5 x'Qx + L'x  s.t. lhs <= Ax <= rhs,
##               lower <= x <= upper
## Accepts ONLY Zero (equality) and NonNeg (inequality) constraints.
## Inherits from QpSolver -- uses QpSolver.apply() for sign-correct data.
##
## Key differences from OSQP:
##   - Q matrix: full symmetric dgCMatrix (NOT upper-triangle)
##   - Dual variables: NEGATE ALL row_duals (matching CVXPY highs_qpif.py line 97)
##   - Status codes: integer codes from R highs (not string enum names)
##   - NOT MIP capable in QP path (MIQP not supported -- CVXPY highs_qpif.py line 37)


# -- HiGHS status map ---------------------------------------------------------
## Reuse HIGHS_STATUS_MAP from highs_conic_solver.R (defined there to avoid
## duplicate, since both files are loaded)

# -- HiGHS_QP_Solver class ----------------------------------------------------
## CVXPY SOURCE: highs_qpif.py lines 30-130

HiGHS_QP_Solver <- new_class("HiGHS_QP_Solver", parent = QpSolver,
  package = "CVXR",
  constructor = function() {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    .fast_new(HiGHS_QP_Solver, S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = TRUE,
      PSD_TRIANGLE_KIND = NA_character_,
      PSD_SQRT2_SCALING = NA,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg),
      REQUIRES_CONSTR = FALSE
    )
  }
)

method(solver_name, HiGHS_QP_Solver) <- function(x) HIGHS_SOLVER

# -- solve_via_data ------------------------------------------------------------
## CVXPY SOURCE: highs_qpif.py lines 122-170
## Receives QP data from QpSolver.apply(): P, q, A_eq, b_eq, F_ineq, g_ineq

method(solve_via_data, HiGHS_QP_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                                       solver_opts = list(), ...) {
  .require_solver_package(HIGHS_SOLVER)

  ## CVXPY SOURCE: highs_qpif.py:123 (warm_start, solver_cache plumbed via ...).
  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]

  L_vec <- data[["q"]]
  nvars <- length(L_vec)

  ## Q matrix (quadratic objective) -- full symmetric dgCMatrix for HiGHS.
  ## DIFFERS from OSQP (upper-triangle) and Clarabel/SCS (dsCMatrix).
  ##
  ## We symmetrize Q below (`(Q + t(Q))/2`). Rationale + a known deficiency:
  ##  - A quadratic form only sees the symmetric part of Q, so (Q + t(Q))/2 is
  ##    the EXACT canonical Hessian (not an approximation): x'Qx is identical.
  ##  - HiGHS >= 1.14 added a strict Hessian-symmetry check and REJECTS a Q
  ##    whose two triangles differ by fp-epsilon (e.g. a dense
  ##    quad_form(M' x, Sigma)), aborting in hi_new_solver() with "Square
  ##    Hessian contains N non-symmetries" (reproduced; CVXPY #3301). HiGHS
  ##    1.12 (current CRAN/supported target) has no such check. Symmetrizing
  ##    satisfies both; it is a no-op on 1.12.
  ##  - DEFICIENCY / FUTURE WORK (in OUR highs interface, not CVXR): CVXPY
  ##    never materializes a square Hessian -- it sends only the upper triangle
  ##    with HessianFormat.kTriangular (highs_qpif.py:199-204). CVXR cannot
  ##    mirror that today: the R `highs` package hardcodes the "square" Hessian
  ##    format in EVERY public path (model_set_hessian, R/highs.R); the
  ##    "triangular" format exists in C++ (highs_interface.cpp) but is reachable
  ##    only via the non-exported highs:::model_set_hessian_. The proper fix is
  ##    to expose a triangular Hessian in the R highs API, then drop this
  ##    symmetrize and pass triu(Q) to match CVXPY exactly. Until then,
  ##    symmetrize. See #3301 and test-infeasible-duals-v19.R.
  if (!is.null(data[[SD_P]])) {
    Q <- methods::as(methods::as(data[[SD_P]], "generalMatrix"), "CsparseMatrix")
    ## Force exact symmetry (see note above): canonical Hessian, no-op on
    ## highs 1.12, required by highs >= 1.14's symmetry check.
    Q <- (Q + Matrix::t(Q)) / 2
  } else {
    Q <- NULL
  }

  ## Stack A_eq and F_ineq into combined HiGHS constraint matrix
  ## Bounds: lhs = [b_eq, -inf*ones], rhs = [b_eq, g_ineq]
  A_eq <- data[["A_eq"]]
  b_eq <- data[["b_eq"]]
  F_ineq <- data[["F_ineq"]]
  g_ineq <- data[["g_ineq"]]

  len_eq <- nrow(A_eq)
  len_ineq <- nrow(F_ineq)

  if (len_eq > 0L && len_ineq > 0L) {
    A <- rbind(A_eq, F_ineq)
    lhs <- c(b_eq, rep(-Inf, len_ineq))
    rhs <- c(b_eq, g_ineq)
  } else if (len_eq > 0L) {
    A <- A_eq
    lhs <- b_eq
    rhs <- b_eq
  } else if (len_ineq > 0L) {
    A <- F_ineq
    lhs <- rep(-Inf, len_ineq)
    rhs <- g_ineq
  } else {
    A <- Matrix::sparseMatrix(i = integer(0), j = integer(0),
                               dims = c(0L, nvars))
    lhs <- numeric(0)
    rhs <- numeric(0)
  }

  ## Ensure A is dgCMatrix for highs
  if (!is.null(A) && !inherits(A, "dgCMatrix")) {
    A <- methods::as(A, "dgCMatrix")
  }

  ## Variable bounds: default to (-Inf, Inf) -- no MIP handling in QP path
  lower <- data[[LOWER_BOUNDS]] %||% rep(-Inf, nvars)
  upper <- data[[UPPER_BOUNDS]] %||% rep(Inf, nvars)

  ## All continuous (QP path is NOT MIP-capable)
  types <- rep(1L, nvars)

  ## Build HiGHS control
  ctrl <- highs::highs_control()
  ctrl$log_to_console <- verbose
  ## CVXPY SOURCE: highs_qpif.py:208 sets only log_to_console. R-SPECIFIC: the
  ## R `highs` 1.14 persistent-solver API installs a log callback whose console
  ## output is gated by `output_flag` (master switch, default TRUE), so
  ## log_to_console alone no longer silences HiGHS. Mirror CVXPY's intent
  ## (quiet unless verbose) by also setting output_flag.
  ctrl$output_flag <- verbose

  ## Apply user-specified solver options
  for (opt_name in names(solver_opts)) {
    ctrl[[opt_name]] <- solver_opts[[opt_name]]
  }

  ## CVXPY SOURCE: highs_qpif.py:175-248.
  ## Persistent-solver flow (highs_model -> hi_new_solver -> optional
  ## hi_solver_set_solution -> hi_solver_run -> getter calls) replaces
  ## the one-shot highs::highs_solve() call so that warm-start can feed
  ## the prior solution via hi_solver_set_solution(), mirroring CVXPY's
  ## `solver.setSolution()` at highs_qpif.py:232.  Requires highs >= 1.14.
  model <- highs::highs_model(
    Q       = Q,
    L       = L_vec,
    lower   = lower,
    upper   = upper,
    A       = A,
    lhs     = lhs,
    rhs     = rhs,
    types   = types,
    maximum = FALSE,
    offset  = 0
  )
  solver <- highs::hi_new_solver(model)
  highs::hi_solver_set_options(solver, ctrl)

  ## CVXPY SOURCE: highs_qpif.py:228-232 (warm-start primal/dual feed-in).
  cache_key <- HIGHS_SOLVER
  if (warm_start && !is.null(solver_cache) &&
      exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    old_status <- HIGHS_STATUS_MAP[[as.character(cached$result$status)]]
    nrow_A <- if (is.null(A)) 0L else nrow(A)
    if (!is.null(old_status) && old_status %in% SOLUTION_PRESENT &&
        length(cached$result$solver_msg$col_value) == nvars &&
        length(cached$result$solver_msg$row_value) == nrow_A) {
      tryCatch({
        prior <- cached$result$solver_msg
        highs::hi_solver_set_solution(
          solver,
          col_value   = prior$col_value,
          row_value   = prior$row_value,
          col_dual    = prior$col_dual,
          row_dual    = prior$row_dual,
          value_valid = isTRUE(prior$value_valid),
          dual_valid  = isTRUE(prior$dual_valid)
        )
      }, error = function(e) NULL)
    }
  }

  ## CVXPY SOURCE: highs_qpif.py:235-243 (run + collect result fields).
  highs::hi_solver_run(solver)
  solution <- highs::hi_solver_get_solution(solver)
  info <- highs::hi_solver_info(solver)
  result <- list(
    primal_solution = solution[["col_value"]],
    objective_value = info[["objective_function_value"]],
    status          = highs::hi_solver_status(solver),
    status_message  = highs::hi_solver_status_message(solver),
    solver_msg      = solution,
    info            = info
  )

  ## CVXR DEVIATION (deliberate; CVXPY's QP path has NO dual ray -- its
  ## highs_qpif.py:117 calls failure_solution(status, attr) with no duals, and
  ## only the conic highs_conif.py:348-349 captures getDualRay()).
  ##
  ## Why CVXR needs it here too: CVXPY 1.9's named-solver preference
  ## (solving_chain.py:313-321) puts the CONIC instance first unless the
  ## objective is quadratic, so an infeasible LP with solver="HIGHS" reaches
  ## highs_conif.py.  CVXR's .solve_as_qp() (solving_chain.R:170-178) instead
  ## only prefers conic when setdiff(conic_solvers, qp_solvers) is non-empty --
  ## naming HIGHS makes both sets {HIGHS}, so the same LP lands HERE.  Wiring
  ## the ray only into highs_conif.R would leave the certificate unreachable
  ## on CVXR's default path for a named-HiGHS LP.  The routing divergence is
  ## tracked separately; propagating the certificate on both paths is correct
  ## regardless of which one a given problem takes.
  if (identical(HIGHS_STATUS_MAP[[as.character(result$status)]], INFEASIBLE)) {
    result$dual_ray <- tryCatch(highs::hi_solver_get_dual_ray(solver),
                                error = function(e) NULL)
  }

  ## CVXPY SOURCE: highs_qpif.py:247-248 (cache for next warm-start).
  if (!is.null(solver_cache)) {
    assign(cache_key,
           list(solver = solver, result = result),
           envir = solver_cache)
  }

  ## Store len_eq for dual splitting (used by reduction_invert).
  result$.len_eq <- len_eq
  result
}

# -- reduction_invert ----------------------------------------------------------
## CVXPY SOURCE: highs_qpif.py lines 76-120
## Dual sign: negate ALL row_duals -- matching CVXPY highs_qpif.py line 97.

method(reduction_invert, HiGHS_QP_Solver) <- function(x, solution, inverse_data, ...) {
  attr_list <- list()

  ## Map status via integer status code.
  ## Same deliberate divergence from upstream's `s.UNKNOWN` fallback as on the
  ## conic side -- see the comment at highs_conif.R's reduction_invert and
  ## notes/upstream_report_cvxpy_highs_unknown.md.
  status_code <- solution$status
  status <- .status_from_map(HIGHS_STATUS_MAP, status_code)

  ## Timing and iteration info
  info <- solution$info
  if (!is.null(info)) {
    num_iters <- (info$simplex_iteration_count %||% 0L) +
                 (info$ipm_iteration_count %||% 0L) +
                 (info$qp_iteration_count %||% 0L) +
                 (info$crossover_iteration_count %||% 0L)
    if (num_iters > 0L) attr_list[[RK_NUM_ITERS]] <- num_iters
  }

  if (status %in% SOLUTION_PRESENT) {
    ## Objective value
    opt_val <- solution$objective_value + inverse_data[[SD_OFFSET]]

    ## Primal variables
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution$primal_solution

    ## Dual variables: negate ALL row_duals
    ## CVXPY SOURCE: highs_qpif.py line 97: y = -np.array(results["solution"].row_dual)
    if (!is.null(solution$solver_msg) &&
        isTRUE(solution$solver_msg$dual_valid)) {
      raw_dual <- solution$solver_msg$row_dual
      y <- -raw_dual  # negate ALL
      len_eq <- solution$.len_eq

      eq_dual <- if (len_eq > 0L) {
        get_dual_values(
          y[seq_len(len_eq)],
          extract_dual_value,
          inverse_data[[SOLVER_EQ_CONSTR]]
        )
      } else {
        list()
      }

      ineq_dual <- if (len_eq < length(y)) {
        get_dual_values(
          y[(len_eq + 1L):length(y)],
          extract_dual_value,
          inverse_data[[SOLVER_NEQ_CONSTR]]
        )
      } else {
        list()
      }

      dual_vars <- c(eq_dual, ineq_dual)
    } else {
      dual_vars <- list()
    }

    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    ## Infeasibility certificate (the dual ray), mapped exactly as the
    ## optimal-case duals above: negate, then split A_eq rows from F_ineq rows
    ## at .len_eq.  See the CVXR DEVIATION note in solve_via_data() for why the
    ## QP path carries this at all when CVXPY's does not, and the R-SPECIFIC
    ## guard note in highs_conif.R for why has_dual_ray must be checked.
    dual_vars <- list()
    ray <- solution$dual_ray
    if (status == INFEASIBLE && !is.null(ray) && isTRUE(ray$has_dual_ray)) {
      y <- -as.numeric(ray$dual_ray)
      len_eq <- solution$.len_eq %||% 0L
      if (length(y) >= len_eq) {
        eq_dual <- if (len_eq > 0L) {
          get_dual_values(y[seq_len(len_eq)], extract_dual_value,
                          inverse_data[[SOLVER_EQ_CONSTR]])
        } else {
          list()
        }
        ineq_dual <- if (len_eq < length(y)) {
          get_dual_values(y[(len_eq + 1L):length(y)], extract_dual_value,
                          inverse_data[[SOLVER_NEQ_CONSTR]])
        } else {
          list()
        }
        dual_vars <- c(eq_dual, ineq_dual)
      }
    }
    failure_solution(status, attr_list, dual_vars)
  }
}

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

method(print, HiGHS_QP_Solver) <- function(x, ...) {
  cat("HiGHS_QP_Solver()\n")
  invisible(x)
}