R/057_atoms_affine_upper_tri.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/atoms/affine/upper_tri.R
#####

## CVXPY SOURCE: atoms/affine/upper_tri.py
## UpperTri -- strict upper triangle of a square matrix as a vector


UpperTri <- new_class("UpperTri", parent = AffAtom, package = "CVXR",
  constructor = function(x, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    n <- .shape(x)[1L]
    ## strict upper triangle has n*(n-1)/2 entries
    entries <- (n * (n - 1L)) %/% 2L
    shape <- c(entries, 1L)

    obj <- .fast_new(UpperTri, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(x),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

method(validate_arguments, UpperTri) <- function(x) {
  arg <- .args(x)[[1L]]
  if (.shape(arg)[1L] != .shape(arg)[2L]) {
    cli_abort("{.cls UpperTri} requires a square matrix, got shape ({arg@shape[1L]}, {arg@shape[2L]}).")
  }
  invisible(NULL)
}

method(shape_from_args, UpperTri) <- function(x) {
  n <- .arg_shape(x)[1L]
  c((n * (n - 1L)) %/% 2L, 1L)
}

# -- log-log: affine (CVXPY upper_tri.py) ------------------------
method(is_atom_log_log_convex, UpperTri) <- function(x) TRUE
method(is_atom_log_log_concave, UpperTri) <- function(x) TRUE

method(numeric_value, UpperTri) <- function(x, values, ...) {
  m <- values[[1L]]
  ## CVXPY uses np.triu_indices(n, k=1, m) which returns row-major order:
  ##   (0,1), (0,2), (0,3), (1,2), (1,3), (2,3) for a 4x4 matrix.
  ## R's upper.tri() + m[mask] returns column-major order:
  ##   (1,2), (1,3), (2,3), (1,4), (2,4), (3,4) for a 4x4 matrix.
  ## We need to extract indices sorted by row then column.
  idx <- which(upper.tri(m), arr.ind = TRUE)
  idx <- idx[order(idx[, 1L], idx[, 2L]), , drop = FALSE]
  matrix(m[idx], ncol = 1L)
}

method(graph_implementation, UpperTri) <- function(x, arg_objs, shape, data = NULL, ...) {
  list(upper_tri_linop(arg_objs[[1L]]), list())
}

#' Extract strict upper triangle of a square matrix
#'
#' @param x An Expression (square matrix)
#' @returns An UpperTri atom (column vector)
#' @export
upper_tri <- function(x) {
  UpperTri(x)
}

#' Reshape a vector into an upper triangular matrix
#'
#' Inverts \code{\link{upper_tri}}. Takes a flat vector and returns an
#' upper triangular matrix (row-major order, matching CVXPY convention).
#'
#' @param expr An Expression (vector).
#' @param strict Logical. If TRUE, returns a strictly upper triangular matrix
#'   (diagonal is zero). If FALSE, includes the diagonal. Default is FALSE.
#' @returns An Expression representing the upper triangular matrix.
#' @export
vec_to_upper_tri <- function(expr, strict = FALSE) {
  expr <- as_expr(expr)

  ## Flatten to column vector if needed
  if (!expr_is_vector(expr)) {
    cli_abort("{.fn vec_to_upper_tri} requires a vector input.")
  }
  ## If expr is (n, 1) shape, get the length
  ell <- expr_size(expr)

  ## Compute n from the triangular number
  if (strict) {
    ## n * (n-1) / 2 == ell
    n <- as.integer((sqrt(8 * ell + 1) + 1) / 2)
    if ((n * (n - 1L)) %/% 2L != ell) {
      cli_abort("Vector length {ell} is not a strict triangular number.")
    }
  } else {
    ## n * (n+1) / 2 == ell
    n <- as.integer((sqrt(8 * ell + 1) - 1) / 2)
    if ((n * (n + 1L)) %/% 2L != ell) {
      cli_abort("Vector length {ell} is not a triangular number.")
    }
  }

  ## Flatten expr to (ell, 1) if not already
  if (.shape(expr)[1L] != ell || .shape(expr)[2L] != 1L) {
    expr <- vec(expr)
  }

  ## Build sparse coefficient matrix P
  ## P maps: (P @ expr).reshape((n, n)) is upper triangular
  ## CVXPY: row_idx = n * row + col for (row, col) in triu_indices(n, k)
  ## This places each vector entry at the correct position in the flattened nxn matrix
  k <- if (strict) 1L else 0L
  ## Get upper triangular indices (row-major order, matching CVXPY)
  idx <- which(upper.tri(matrix(0, n, n), diag = !strict), arr.ind = TRUE)
  idx <- idx[order(idx[, 1L], idx[, 2L]), , drop = FALSE]
  row_0 <- idx[, 1L] - 1L  # 0-based row
  col_0 <- idx[, 2L] - 1L  # 0-based col

  ## P_rows: flattened position in row-major nxn matrix (0-based)
  P_rows <- n * row_0 + col_0
  P_cols <- seq_len(ell) - 1L  # 0-based column indices
  P_vals <- rep(1, ell)

  ## Build sparse matrix (convert to 1-based for R)
  P <- Matrix::sparseMatrix(
    i = P_rows + 1L,
    j = P_cols + 1L,
    x = P_vals,
    dims = c(as.integer(n * n), as.integer(ell))
  )

  ## P %*% expr gives (n*n, 1) column, reshape to (n, n) in F-order then transpose
  ## CVXPY: reshape(P @ expr, (n, n), order='F').T
  ## In R: reshape_expr(..., c(n, n), order = "F") then transpose
  result <- Constant(P) %*% expr
  result <- reshape_expr(result, c(n, n), order = "F")
  Transpose(result)
}


# -- upper_tri_to_full ----------------------------------------------
## CVXPY SOURCE: atoms/affine/upper_tri.py lines 151-176
## IMPLEMENTED AT: src/RcppConv.cpp (wrapper generated into lin_ops/RcppExports.R)
##
## Returns an (n^2 x n*(n+1)/2) sparse matrix A such that `(A %*% v)` reshaped
## as (n, n) is symmetric. This is CVXPY's isomorphic home for the function; the
## implementation is native because the R generated wrapper cannot live anywhere
## but lin_ops/RcppExports.R, which is Rcpp's output path.
##
## HISTORY, because it cost a day to find. An R implementation of this function
## lived in `reductions/cvx_attr2constr.R` -- the wrong file -- from Phase 5b
## (ec9bb97) until 2026-08-14, and silently SHADOWED the C++ routine, which
## predates it (9d0a95a, Phase 0) and exists precisely to make this fast.
## `copy_r_source.R` prefixes R/ files by load order, so `260_..._cvx_attr2constr`
## loaded after `018_lin_ops_RcppExports` and simply overwrote the binding. No
## error, no warning: the native routine sat compiled, registered and
## unreachable for the whole of that time.
##
## The two were verified equivalent before the R copy was removed: identical
## sparse matrices for n = 1..20, and both satisfy the real contract that
## `(A %*% v)` reshaped is symmetric. The C++ one is 67x faster and allocates
## 10x less (2.38us / 2.14KB against 159.33us / 21.92KB at n = 12).
##
## `scripts/validate_rcpp_bridge.R` fails on exactly this shape now, so a
## duplicate definition cannot bury a native routine again.
##
## NOT PORTED from the same CVXPY file: `batched_upper_tri_to_full`
## (upper_tri.py:183) has no CVXR caller.