R/105_atoms_quad_over_lin.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/quad_over_lin.py
## QuadOverLin -- sum(x^2) / y, convex, domain y > 0
## CVXPY: inherits AxisAtom (supports axis/keepdims reduction)


QuadOverLin <- new_class("QuadOverLin", parent = AxisAtom, package = "CVXR",
  constructor = function(x, y, axis = NULL, keepdims = FALSE, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    y <- as_expr(y)
    if (!is.null(axis)) axis <- as.integer(axis)
    keepdims <- as.logical(keepdims)

    ## Shape: reduce x along axis, then divide by scalar y -> same reduced shape
    if (is.null(axis)) {
      shape <- c(1L, 1L)
    } else {
      shape <- .axis_shape(.shape(x), axis, keepdims)
    }

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

# -- validate: y must be scalar -----------------------------------
## CVXPY: quad_over_lin.py lines 144-152
method(validate_arguments, QuadOverLin) <- function(x) {
  y <- .args(x)[[2L]]
  if (!all(.shape(y) == c(1L, 1L))) {
    cli_abort("{.arg y} must be scalar in {.cls QuadOverLin}.")
  }
  ## CVXPY: quad_over_lin.py line 149-150 -- second arg cannot be complex
  if (is_complex(y)) {
    cli_abort("The second argument to {.cls QuadOverLin} cannot be complex.")
  }
  ## Axis validation (from AxisAtom)
  if (!is.null(x@axis)) {
    .validate_axis(x@axis, 2L)
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
## CVXPY: Inherits AxisAtom.shape_from_args -> reduce first arg along axis
method(shape_from_args, QuadOverLin) <- function(x) {
  .axis_shape(.arg_shape(x), x@axis, x@keepdims)
}

# -- sign: always nonneg ------------------------------------------
method(sign_from_args, QuadOverLin) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature: convex --------------------------------------------
method(is_atom_convex, QuadOverLin) <- function(x) TRUE
method(is_atom_concave, QuadOverLin) <- function(x) FALSE

## CVXPY quad_over_lin.py: quad_over_lin is smooth.
method(is_atom_smooth, QuadOverLin) <- function(x) TRUE

# -- monotonicity -------------------------------------------------
## CVXPY: quad_over_lin.py lines 134-143
## Increasing in x when x >= 0, decreasing when x <= 0
## Decreasing in y (always)
method(is_incr, QuadOverLin) <- function(x, idx, ...) {
  if (idx == 1L) is_nonneg(.args(x)[[1L]]) else FALSE
}
method(is_decr, QuadOverLin) <- function(x, idx, ...) {
  if (idx == 1L) is_nonpos(.args(x)[[1L]]) else TRUE
}

# -- log-log: convex (CVXPY quad_over_lin.py lines 124-132) ------
method(is_atom_log_log_convex, QuadOverLin) <- function(x) TRUE
method(is_atom_log_log_concave, QuadOverLin) <- function(x) FALSE

# -- domain: y >= 0 -----------------------------------------------
method(atom_domain, QuadOverLin) <- function(x) {
  list(.args(x)[[2L]] >= 0)
}

# -- quadratic analysis -------------------------------------------
## CVXPY: quad_over_lin.py lines 154-167
method(is_quadratic, QuadOverLin) <- function(x) {
  is_affine(.args(x)[[1L]]) && is_constant(.args(x)[[2L]])
}

method(has_quadratic_term, QuadOverLin) <- function(x) {
  is_constant(.args(x)[[2L]])
}

method(is_qpwa, QuadOverLin) <- function(x) {
  is_pwl(.args(x)[[1L]]) && is_constant(.args(x)[[2L]])
}

# -- get_data -----------------------------------------------------
method(get_data, QuadOverLin) <- function(x) {
  list(x@axis, x@keepdims)
}

# -- numeric ------------------------------------------------------
## CVXPY: quad_over_lin.py lines 51-61
method(numeric_value, QuadOverLin) <- function(x, values, ...) {
  xv <- values[[1L]]
  yv <- as.numeric(values[[2L]])
  ## For complex x, |x|^2 = Mod(x)^2, not x^2
  squared <- if (is.complex(xv)) Mod(xv)^2 else xv^2
  if (is.null(x@axis)) {
    matrix(sum(squared) / yv, 1L, 1L)
  } else if (x@axis == 2L) {
    res <- apply(squared, 2L, sum) / yv
    if (x@keepdims) matrix(res, nrow = 1L) else matrix(res, nrow = 1L)
  } else {
    res <- apply(squared, 1L, sum) / yv
    if (x@keepdims) matrix(res, ncol = 1L) else matrix(res, ncol = 1L)
  }
}

# -- graph_implementation: stub -----------------------------------
method(graph_implementation, QuadOverLin) <- function(x, arg_objs, shape, data = NULL, ...) {
  cli_abort("graph_implementation for {.cls QuadOverLin} not yet implemented.")
}

# -- .grad: per-atom subgradient ----------------------------------
## CVXPY SOURCE: atoms/quad_over_lin.py:73-103 (quad_over_lin._grad).
## Two args: X (vector / matrix) and y (scalar > 0).
## d/dX = 2X / y; d/dy = -||X||^2_2 / y^2. NULL on y <= 0 or with axis.
## Override .grad directly (rather than using .axis_grad) because the
## gradient is genuinely two-arg and the axis case is unsupported per
## CVXPY (returns [None, None]).
method(.grad, QuadOverLin) <- function(x, values, ...) {
  if (!is.null(x@axis)) return(list(NULL, NULL))
  X <- values[[1L]]
  y <- as.numeric(values[[2L]])
  if (y <= 0) return(list(NULL, NULL))

  ## d/dX = 2X / y, raveled column-major to match CVXPY's order='F'.
  DX <- 2 * as.numeric(X) / y
  ## d/dy = -||X||^2 / y^2, scalar.
  Dy <- -sum(as.numeric(X)^2) / (y^2)

  rows_X <- as.integer(prod(.arg_shape(x)))
  list(
    .dense_to_csc_vector(DX, rows_X),
    Matrix::sparseMatrix(
      i = 1L, j = 1L, x = Dy,
      dims = c(1L, 1L), repr = "C"
    )
  )
}

#' Sum of squares divided by a scalar
#'
#' @param x An Expression
#' @param y An Expression (scalar, positive)
#' @param axis NULL (all), 1 (row-wise), or 2 (column-wise)
#' @param keepdims Logical: keep reduced dimensions?
#' @returns A QuadOverLin atom
#' @export
quad_over_lin <- function(x, y, axis = NULL, keepdims = FALSE) {
  QuadOverLin(x, y, axis, keepdims)
}

#' Sum of squares (= quad_over_lin(x, 1))
#'
#' @param x An Expression
#' @param axis NULL (all), 1 (row-wise), or 2 (column-wise)
#' @param keepdims Logical: keep reduced dimensions?
#' @returns A QuadOverLin atom
#' @export
sum_squares <- function(x, axis = NULL, keepdims = FALSE) {
  QuadOverLin(x, Constant(1), axis = axis, keepdims = keepdims)
}