R/013_utilities_psd_utils.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/utilities/psd_utils.R
#####

## CVXPY SOURCE: utilities/psd_utils.py
##
## Shared PSD/svec helpers.  New at CVXPY 1.9.0 (#3080/#3268): the triangle
## kind and sqrt(2) scaling a solver expects are solver *attributes*
## (`PSD_TRIANGLE_KIND` / `PSD_SQRT2_SCALING`, solver.py:85-90) threaded through
## the chain on `SolverInfo` (solving_chain.py:168-175), and there is ONE
## parameterized packer/unpacker instead of a bespoke pair per solver interface.
##
## CVXR ported this at 1.9.2 under ADR D_19.6, which reopened the D_19.1
## cluster-2 defer.  The decisive evidence was a bug class, not a diff size:
## CVXR had three hand-written expanders (clarabel/scs/mosek) that each
## re-derived the numpy-vs-R triangle-order question independently, and two got
## it wrong -- every PSD dual was transposed for n >= 3 in CRAN 1.9.1 (fixed in
## 893ca4b).  A single shared helper makes that error unrepresentable.

# -- TriangleKind --------------------------------------------------
## CVXPY SOURCE: psd_utils.py lines 23-26 (an `enum.Enum`).
## R has no enum type; a named list of the SAME string values is the
## idiomatic stand-in.  Compare with `identical()`, never with `==` on a
## possibly-NA value.

TriangleKind <- list(LOWER = "lower", UPPER = "upper")

# -- tri_to_full ---------------------------------------------------
## CVXPY SOURCE: psd_utils.py lines 28-79
##
## Expand a vectorized triangle to the full symmetric matrix.
##
## TRIANGLE ORDER (R differs from numpy -- do not "fix" this to match CVXPY's
## index set).  Solvers vectorize their triangle in COLUMN-MAJOR order
## (psd_utils.py:52-61).  numpy's `triu_indices`/`tril_indices` enumerate
## ROW-major, so CVXPY must use the OPPOSITE index set (psd_utils.py:67:
## `triu_indices` for LOWER, `tril_indices` for UPPER).  R's logical-mask
## assignment is ALREADY column-major, so R uses the mask matching the solver's
## own triangle.  Using the opposite mask transposes the layout and silently
## scrambles every PSD dual for n >= 3 (the two orders coincide only at n = 2).
##
## N/A FOR R: the batched `num > 1` path (psd_utils.py:63-79, returning a
## (num, n, n) array).  It exists for CVXPY's N-D PSD constraints; CVXR is 2-D
## (`ALLOW_ND_EXPR = FALSE`) and `num_cones(PSD)` is always 1L (psd.R:67), so a
## single triangle is the only case reachable here.
##
## Unlike the deleted per-solver expanders, this returns the (n, n) MATRIX that
## CVXPY returns, not `as.vector(full)`.  Callers that still need the flat
## column-major form wrap it in `as.vector()`.

tri_to_full <- function(tri_vec, n, triangle, sqrt2_scaling) {
  n <- as.integer(n)
  tri_dim <- (n * (n + 1L)) %/% 2L
  if (length(tri_vec) != tri_dim) {
    cli_abort(c(
      "{.arg tri_vec} has length {length(tri_vec)}, expected {tri_dim} for n = {n}.",
      "i" = "Batched (num > 1) triangles are not reachable in CVXR: expressions are 2-D."
    ))
  }
  full <- matrix(0, n, n)
  mask <- if (identical(triangle, TriangleKind$LOWER)) {
    lower.tri(full, diag = TRUE)
  } else {
    upper.tri(full, diag = TRUE)
  }
  full[mask] <- tri_vec
  full <- full + t(full)
  diag(full) <- diag(full) / 2
  if (isTRUE(sqrt2_scaling)) {
    off_diag <- !diag(TRUE, n)
    full[off_diag] <- full[off_diag] / sqrt(2)
  }
  full
}

# -- psd_format_mat ------------------------------------------------
## CVXPY SOURCE: psd_utils.py lines 82-147
##
## Sparse matrix M with M %*% vec(X, order = "F") = the scaled vectorized
## triangle of X: triangle extraction (scaled by sqrt(2) off the diagonal when
## `sqrt2_scaling`) composed with symmetrization (X + t(X)) / 2.
##
## N/A FOR R: the batched `num > 1` block-diagonal + de-interleaving permutation
## (psd_utils.py:131-147), for the same reason as `tri_to_full` above.

psd_format_mat <- function(constr, triangle, sqrt2_scaling) {
  n <- .arg_shape(constr)[1L]
  entries <- (n * (n + 1L)) %/% 2L

  mask <- if (identical(triangle, TriangleKind$LOWER)) {
    lower.tri(matrix(0, n, n), diag = TRUE)
  } else {
    upper.tri(matrix(0, n, n), diag = TRUE)
  }

  ## Column indices: the triangle's positions as column-major flat indices.
  ## `which()` on a logical matrix already enumerates column-major and returns
  ## them sorted, which is exactly CVXPY's
  ## `np.sort(np.ravel_multi_index(tri_idx, (n, n), order = "F"))`.
  col_arr <- which(mask)
  row_arr <- seq_len(entries)

  ## Value array: sqrt(2) off-diagonal (when scaled), 1.0 on the diagonal.
  val_mat <- matrix(0, n, n)
  val_mat[mask] <- if (isTRUE(sqrt2_scaling)) sqrt(2) else 1.0
  diag(val_mat) <- 1.0
  val_arr <- as.vector(val_mat)      # column-major
  val_arr <- val_arr[val_arr != 0]

  scaled_tri <- Matrix::sparseMatrix(
    i = row_arr, j = col_arr, 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_tri %*% symm_matrix
}