R/237_reductions_solvers_conic_solvers_scs_conif.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: reductions/solvers/conic_solvers/scs_conif.py
## SCS solver interface
##
## SCS uses lower-triangular svec format for PSD constraints.
## Convention: A*x + s = b, s in K


# -- SCS status map ------------------------------------------------
## CVXPY SOURCE: scs_conif.py lines 127-135

SCS_STATUS_MAP <- list(
  "1"  = OPTIMAL,
  "2"  = OPTIMAL_INACCURATE,
  "-1" = UNBOUNDED,
  "-6" = UNBOUNDED_INACCURATE,
  "-2" = INFEASIBLE,
  "-7" = INFEASIBLE_INACCURATE,
  "-4" = SOLVER_ERROR,
  "-3" = SOLVER_ERROR,
  "-5" = SOLVER_ERROR
)

# -- tri_to_full (SCS): expand lower-tri svec to full matrix ------
## CVXPY SOURCE: scs_conif.py lines 45-76
## Scales off-diagonal by 1/sqrt(2) as per SCS spec.
## Returns a numeric vector of length n*n (column-major).

scs_tri_to_full <- function(lower_tri, n) {
  full <- matrix(0, n, n)
  ## SCS tracks "lower triangular" in a way that corresponds to numpy upper tri
  ## (confusing but correct, per SCS docs)
  full[upper.tri(full, diag = TRUE)] <- lower_tri
  full <- full + t(full)
  diag(full) <- diag(full) / 2
  ## Scale off-diagonal by 1/sqrt(2)
  off_lower <- lower.tri(full)
  off_upper <- upper.tri(full)
  full[off_lower] <- full[off_lower] / sqrt(2)
  full[off_upper] <- full[off_upper] / sqrt(2)
  as.vector(full)  # column-major
}

# -- SCS psd_format_mat -------------------------------------------
## CVXPY SOURCE: scs_conif.py lines 171-206
## Lower-triangular extraction with sqrt(2) off-diagonal scaling
## + symmetrization. Returns a sparse matrix (entries x n^2).

scs_psd_format_mat_fn <- function(constr) {
  n <- constr@args[[1L]]@shape[1L]  # constr.expr.shape[0]
  entries <- (n * (n + 1L)) %/% 2L

  ## Row indices: 0 to entries-1
  row_arr <- seq(0L, length.out = entries)

  ## Column indices: lower triangle positions as column-major flat indices
  ## np.tril_indices(n) -> R: which(lower.tri(m, diag=TRUE), arr.ind=TRUE)
  idx <- which(lower.tri(matrix(0, n, n), diag = TRUE), arr.ind = TRUE)
  ## Convert to 0-based column-major flat indices (already sorted)
  col_arr <- (idx[, 2L] - 1L) * n + (idx[, 1L] - 1L)

  ## Value array: sqrt(2) for off-diagonal, 1.0 for diagonal
  val_mat <- matrix(0, n, n)
  val_mat[lower.tri(val_mat, diag = TRUE)] <- sqrt(2)
  diag(val_mat) <- 1.0
  val_arr <- as.vector(val_mat)  # column-major
  val_arr <- val_arr[val_arr != 0]  # nonzero entries

  scaled_lower_tri <- Matrix::sparseMatrix(
    i = row_arr + 1L, j = col_arr + 1L, x = val_arr,
    dims = c(entries, as.integer(n * n))
  )

  ## Symmetrization matrix: (M + M^T) / 2
  nn <- as.integer(n * n)
  K <- matrix(seq_len(nn) - 1L, nrow = n, ncol = n, byrow = TRUE)  # C-order
  row_symm <- c(seq_len(nn) - 1L, as.vector(K)) + 1L      # 1-based
  col_symm <- c(seq_len(nn) - 1L, as.vector(t(K))) + 1L   # 1-based
  val_symm <- rep(0.5, 2L * nn)

  symm_matrix <- Matrix::sparseMatrix(
    i = row_symm, j = col_symm, x = val_symm,
    dims = c(nn, nn)
  )

  scaled_lower_tri %*% symm_matrix
}

# -- dims_to_solver_dict_scs --------------------------------------
## CVXPY SOURCE: scs_conif.py lines 36-42
## SCS 3.0+ uses 'z' instead of 'f' for zero cone.

dims_to_solver_dict_scs <- function(cone_dims) {
  cones <- dims_to_solver_dict(cone_dims)
  ## SCS 3.x renamed 'f' to 'z'
  if (requireNamespace("scs", quietly = TRUE)) {
    scs_ver <- utils::packageVersion("scs")
    if (scs_ver >= "3.0.0") {
      cones[["z"]] <- cones[["f"]]
      cones[["f"]] <- NULL
    }
  }
  cones
}

# -- SCS_Solver class ---------------------------------------------
## CVXPY SOURCE: scs_conif.py lines 116-155

SCS_Solver <- new_class("SCS_Solver", parent = ConicSolver, package = "CVXR",
  constructor = function() {
    new_object(S7_object(),
      .cache = new.env(parent = emptyenv()),
      MIP_CAPABLE = FALSE,
      BOUNDED_VARIABLES = FALSE,
      SUPPORTED_CONSTRAINTS = list(Zero, NonNeg, SOC, ExpCone,
                                   PSD, PowCone3D),
      EXP_CONE_ORDER = c(0L, 1L, 2L),
      REQUIRES_CONSTR = TRUE
    )
  }
)

method(solver_name, SCS_Solver) <- function(x) SCS_SOLVER

## Override PSD format for SCS (lower-triangular + sqrt(2))
method(solver_psd_format_mat, SCS_Solver) <- function(solver, constr) {
  scs_psd_format_mat_fn(constr)
}

# -- SCS extract_dual_value ----------------------------------------
## CVXPY SOURCE: scs_conif.py lines 218-233
## PSD constraints use lower-tri svec format; expand to full.

scs_extract_dual_value <- function(result_vec, offset, constraint) {
  if (S7_inherits(constraint, PSD)) {
    dim <- constraint@args[[1L]]@shape[1L]
    lower_tri_dim <- (dim * (dim + 1L)) %/% 2L
    new_offset <- offset + lower_tri_dim
    lower_tri <- result_vec[(offset + 1L):(offset + lower_tri_dim)]
    full <- scs_tri_to_full(lower_tri, dim)
    list(value = full, new_offset = new_offset)
  } else {
    extract_dual_value(result_vec, offset, constraint)
  }
}

# -- SCS invert ----------------------------------------------------
## CVXPY SOURCE: scs_conif.py lines 235-278
## Parses SCS-specific result format.

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

  ## SCS result format: solution is the raw list from scs::scs()
  ## SCS 3.x: status_val in info, SCS 2.x: statusVal in info
  info <- solution[["info"]]
  status_val <- info[["status_val"]]
  if (is.null(status_val)) status_val <- info[["statusVal"]]
  status <- SCS_STATUS_MAP[[as.character(status_val)]]
  if (is.null(status)) status <- SOLVER_ERROR

  ## Timing attributes
  solve_time <- info[["solve_time"]]
  if (is.null(solve_time)) solve_time <- info[["solveTime"]]
  setup_time <- info[["setup_time"]]
  if (is.null(setup_time)) setup_time <- info[["setupTime"]]
  if (!is.null(solve_time)) attr_list[[RK_SOLVE_TIME]] <- solve_time / 1000
  if (!is.null(setup_time)) attr_list[[RK_SETUP_TIME]] <- setup_time / 1000
  attr_list[[RK_NUM_ITERS]] <- info[["iter"]]

  if (status %in% SOLUTION_PRESENT) {
    primal_val <- info[["pobj"]]
    opt_val <- primal_val + inverse_data[[SD_OFFSET]]
    primal_vars <- list()
    primal_vars[[as.character(inverse_data[[SOLVER_VAR_ID]])]] <- solution[["x"]]

    ## Dual variables: split at zero cone boundary
    y <- solution[["y"]]
    zero_dim <- inverse_data[[SD_DIMS]]@zero
    if (zero_dim > 0L) {
      eq_dual <- get_dual_values(
        y[seq_len(zero_dim)],
        scs_extract_dual_value,
        inverse_data[[SOLVER_EQ_CONSTR]]
      )
    } else {
      eq_dual <- list()
    }
    if (zero_dim < length(y)) {
      ineq_dual <- get_dual_values(
        y[(zero_dim + 1L):length(y)],
        scs_extract_dual_value,
        inverse_data[[SOLVER_NEQ_CONSTR]]
      )
    } else {
      ineq_dual <- list()
    }
    dual_vars <- c(eq_dual, ineq_dual)
    Solution(status, opt_val, primal_vars, dual_vars, attr_list)
  } else {
    failure_solution(status, attr_list)
  }
}

# -- SCS solve_via_data --------------------------------------------
## CVXPY SOURCE: scs_conif.py lines 304-354

method(solve_via_data, SCS_Solver) <- function(x, data, warm_start = FALSE, verbose = FALSE,
                                               solver_opts = list(), ...) {
  if (!requireNamespace("scs", quietly = TRUE)) {
    cli_abort("Package {.pkg scs} is required but not installed.")
  }

  dots <- list(...)
  solver_cache <- dots[["solver_cache"]]

  args <- list(A = data[[SD_A]], b = data[[SD_B]], c = data[[SD_C]])
  cone <- dims_to_solver_dict_scs(data[[SD_DIMS]])

  ## Parse solver options
  opts <- solver_opts
  scs_ver <- utils::packageVersion("scs")
  if (scs_ver >= "3.0.0") {
    if (is.null(opts[["eps_abs"]])) opts[["eps_abs"]] <- 1e-5
    if (is.null(opts[["eps_rel"]])) opts[["eps_rel"]] <- 1e-5
  } else {
    if (is.null(opts[["eps"]])) opts[["eps"]] <- 1e-4
  }

  ## Pass P for QP path -- SCS expects symmetric sparse (dsCMatrix)
  if (!is.null(data[[SD_P]])) {
    args[["P"]] <- Matrix::forceSymmetric(Matrix::triu(data[[SD_P]]), uplo = "U")
  }

  ## Warm-start: pass previous x/y/s as initial point
  ## CVXPY SOURCE: scs_conif.py lines 327-331
  cache_key <- SCS_SOLVER
  initial <- NULL
  if (warm_start && !is.null(solver_cache) && exists(cache_key, envir = solver_cache)) {
    cached <- get(cache_key, envir = solver_cache)
    initial <- list(x = cached$x, y = cached$y, s = cached$s)
  }

  ## Call SCS
  result <- scs::scs(
    A = args$A, b = args$b, obj = args$c, P = args[["P"]], cone = cone,
    initial = initial,
    control = c(list(verbose = verbose), opts)
  )

  ## Cache result for future warm-starts (only on optimal)
  ## CVXPY SOURCE: scs_conif.py lines 352-353
  if (!is.null(solver_cache)) {
    status <- SCS_STATUS_MAP[[as.character(result$info$status_val)]]
    if (!is.null(status) && status == OPTIMAL) {
      assign(cache_key, result, envir = solver_cache)
    }
  }

  result
}

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