R/125_constraints_second_order.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/constraints/second_order.R
#####

## CVXPY SOURCE: constraints/second_order.py
## SOC -- Second-Order Cone constraint

#' Create a Second-Order Cone Constraint
#'
#' Constrains \eqn{\|X_i\|_2 \le t_i} for each column or row \eqn{i},
#' where \code{t} is a vector and \code{X} is a matrix.
#'
#' @param t A CVXR expression (scalar or vector) representing the upper bound.
#' @param X A CVXR expression (vector or matrix) whose columns/rows are bounded.
#' @param axis Integer, 2 (default, column-wise) or 1 (row-wise). Determines whether columns (\code{2})
#'   or rows (\code{1}) of \code{X} define the individual cones.
#' @param constr_id Optional integer constraint ID.
#' @returns An \code{SOC} constraint object.
#' @export
SOC <- new_class("SOC", parent = Cone, package = "CVXR",
  properties = list(
    axis = new_property(class = class_integer, default = 2L)
  ),
  constructor = function(t, X, axis = 2L, constr_id = NULL) {
    t <- as_expr(t)
    axis <- as.integer(axis)

    ## CVXPY SOURCE: second_order.py lines 40-41
    ## t must not be a matrix (len(t.shape) >= 2 means genuinely 2D)
    if (expr_is_matrix(t) || !is_real(t)) {
      cli_abort("Invalid first argument: {.arg t} must be a vector and real.")
    }

    ## CVXPY SOURCE: second_order.py lines 43-49
    ## Check t has one entry per cone
    t_size <- expr_size(t)
    X_shape <- X@shape
    X_ndim <- length(X_shape)
    if ((X_ndim <= 1L && t_size > 1L) ||
        (X_ndim == 2L && t_size != X_shape[axis]) ||
        (X_ndim == 1L && axis == 1L)) {
      cli_abort(
        "Argument dimensions {.val {t@shape}} and {.val {X_shape}}, with axis={axis}, are incompatible."
      )
    }

    ## CVXPY SOURCE: second_order.py line 52 -- flatten 0-d t
    if (expr_is_scalar(t)) {
      t <- reshape_expr(t, c(1L, 1L))
    }

    if (is.null(constr_id)) constr_id <- next_expr_id()
    shape <- t@shape

    ## Both t and X go into args; dual variables auto-created for both
    args <- list(t, X)
    dvars <- lapply(args, function(a) Variable(a@shape))

    new_object(S7_object(),
      id    = as.integer(constr_id),
      .cache = new.env(parent = emptyenv()),
      args  = args,
      dual_variables = dvars,
      shape = shape,
      .label = "",
      axis  = axis
    )
  }
)

# -- expr_name ----------------------------------------------------
## CVXPY SOURCE: second_order.py line 55-56

method(expr_name, SOC) <- function(x) {
  sprintf("SOC(%s, %s)", expr_name(x@args[[1L]]), expr_name(x@args[[2L]]))
}

# -- is_dcp -------------------------------------------------------
## CVXPY SOURCE: second_order.py lines 158-164

method(is_dcp, SOC) <- function(x) {
  .all_args(x, is_affine)
}

# -- is_dgp -------------------------------------------------------
## CVXPY SOURCE: second_order.py lines 166-167

method(is_dgp, SOC) <- function(x) FALSE

# -- get_data -----------------------------------------------------
## CVXPY SOURCE: second_order.py lines 117-124

method(get_data, SOC) <- function(x) {
  list(x@axis, x@id)
}

# -- num_cones ----------------------------------------------------
## CVXPY SOURCE: second_order.py lines 126-129

method(num_cones, SOC) <- function(x) {
  expr_size(x@args[[1L]])
}

# -- .cone_size (internal helper) ---------------------------------
## CVXPY SOURCE: second_order.py lines 131-140

.soc_cone_size <- function(x) {
  X <- x@args[[2L]]
  X_shape <- X@shape
  ## Our shapes are always length 2 in R; CVXPY checks len(X.shape) == 0
  ## For scalar X: X_dim = 1
  if (expr_is_scalar(X)) {
    X_dim <- 1L
  } else {
    ## X.shape[self.axis] -> cone dimension = shape[3 - axis]
    X_dim <- X_shape[3L - x@axis]
  }
  1L + X_dim
}

# -- cone_sizes ---------------------------------------------------
## CVXPY SOURCE: second_order.py lines 148-156

method(cone_sizes, SOC) <- function(x) {
  rep(.soc_cone_size(x), num_cones(x))
}

# -- constr_size --------------------------------------------------
## CVXPY SOURCE: second_order.py lines 142-146

method(constr_size, SOC) <- function(x) {
  .soc_cone_size(x) * num_cones(x)
}

# -- residual -----------------------------------------------------
## CVXPY SOURCE: second_order.py lines 58-115
## SOC projection residual: ||(t,X) - proj(t,X)||

method(residual, SOC) <- function(x) {
  t_val <- value(x@args[[1L]])
  X_val <- value(x@args[[2L]])
  if (is.null(t_val) || is.null(X_val)) return(NULL)

  ## Ensure matrix form
  if (!is.matrix(t_val)) t_val <- matrix(as.numeric(t_val), ncol = 1L)
  if (!is.matrix(X_val)) X_val <- matrix(as.numeric(X_val), ncol = 1L)

  ## Reduce axis=2 to axis=1 by transposing X
  if (x@axis == 2L) X_val <- t(X_val)

  ## Track if X was promoted from 1D
  promoted <- (nrow(X_val) == 1L && ncol(X_val) > 1L && length(x@args[[2L]]@shape) <= 1L)
  if (nrow(X_val) == 1L && !promoted) {
    ## Already 2D-ish, keep as is
  }

  ## t as vector
  tv <- as.numeric(t_val)

  ## Compute row-wise norms of X
  norms <- sqrt(rowSums(X_val^2))
  n_cones <- length(tv)

  ## Initialize projections to zero (default: t <= -||x|| case)
  t_proj <- numeric(n_cones)
  X_proj <- matrix(0, nrow = nrow(X_val), ncol = ncol(X_val))

  ## Case 1: t >= ||x|| -> proj = (t, X)
  mask1 <- tv >= norms
  t_proj[mask1] <- tv[mask1]
  X_proj[mask1, ] <- X_val[mask1, , drop = FALSE]

  ## Case 2: -||x|| < t < ||x|| -> proj = 0.5*(t/||x|| + 1)*(||x||, x)
  mask2 <- abs(tv) < norms
  avg_coeff <- 0.5 * (1 + tv / norms)
  if (any(mask2)) {
    X_proj[mask2, ] <- avg_coeff[mask2] * X_val[mask2, , drop = FALSE]
    t_proj[mask2] <- avg_coeff[mask2] * norms[mask2]
  }

  ## Compute residuals: ||(t,X) - proj(t,X)||_2 per cone
  Xt <- cbind(X_val, tv)
  Xt_proj <- cbind(X_proj, t_proj)
  resid <- sqrt(rowSums((Xt - Xt_proj)^2))

  if (promoted || n_cones == 1L) resid[1L] else resid
}

# -- save_dual_value ----------------------------------------------
## CVXPY SOURCE: second_order.py lines 172-183
## Uses C-order reshape (CR-1)

method(save_dual_value, SOC) <- function(x, val) {
  cone_size <- .soc_cone_size(x)
  ## CVXPY: np.reshape(value, (-1, cone_size)) in C-order
  val_mat <- .reshape_c_order(val, length(val) %/% cone_size, cone_size)
  t_dual <- val_mat[, 1L]
  X_dual <- val_mat[, -1L, drop = FALSE]

  ## CVXPY: if len(self.args[1].shape) == 0 -> scalar X
  if (expr_is_scalar(x@args[[2L]])) {
    X_dual <- X_dual[1L, 1L]
  } else if (x@axis == 2L) {
    X_dual <- t(X_dual)
  }

  value(x@dual_variables[[1L]]) <- t_dual
  value(x@dual_variables[[2L]]) <- X_dual
  invisible(x)
}

# -- dual_cone ----------------------------------------------------
## CVXPY SOURCE: second_order.py lines 185-196
## SOC is self-dual

method(dual_cone, SOC) <- function(x, ...) {
  args <- list(...)
  if (length(args) == 0L) {
    SOC(x@dual_variables[[1L]], x@dual_variables[[2L]], axis = x@axis)
  } else {
    SOC(args[[1L]], args[[2L]], axis = x@axis)
  }
}