R/095_atoms_quad_form.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/quad_form.py
## QuadForm -- quadratic form x^T P x


QuadForm <- new_class("QuadForm", parent = Atom, package = "CVXR",
  constructor = function(x, P, id = NULL) {
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    P <- as_expr(P)
    ## Shape is always scalar (1, 1)
    shape <- c(1L, 1L)

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

# -- validate -----------------------------------------------------
method(validate_arguments, QuadForm) <- function(x) {
  P <- x@args[[2L]]
  xarg <- x@args[[1L]]
  ## P must be square
  if (P@shape[1L] != P@shape[2L]) {
    cli_abort("{.arg P} must be square, got shape ({P@shape[1L]}, {P@shape[2L]}).")
  }
  ## x must be a vector compatible with P
  n <- P@shape[1L]
  if (xarg@shape[1L] != n) {
    cli_abort("{.arg x} must have {n} rows to match P, got {xarg@shape[1L]}.")
  }
  ## P must be symmetric/hermitian (CVXPY quad_form.py line 57)
  if (!is_hermitian(P)) {
    cli_abort("Quadratic form matrices must be symmetric/Hermitian.")
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
method(shape_from_args, QuadForm) <- function(x) c(1L, 1L)

# -- sign: depends on P definiteness ------------------------------
method(sign_from_args, QuadForm) <- function(x) {
  list(is_nonneg = is_atom_convex(x), is_nonpos = is_atom_concave(x))
}

# -- curvature: depends on P --------------------------------------
## CVXPY: convex iff P is constant + PSD; concave iff P constant + NSD
method(is_atom_convex, QuadForm) <- function(x) {
  P <- x@args[[2L]]
  is_constant(P) && is_psd(P)
}

method(is_atom_concave, QuadForm) <- function(x) {
  P <- x@args[[2L]]
  is_constant(P) && is_nsd(P)
}

# -- log-log curvature (CVXPY quad_form.py lines 76-84) ----------
method(is_atom_log_log_convex, QuadForm) <- function(x) TRUE
method(is_atom_log_log_concave, QuadForm) <- function(x) FALSE

# -- monotonicity -------------------------------------------------
method(is_incr, QuadForm) <- function(x, idx, ...) FALSE
method(is_decr, QuadForm) <- function(x, idx, ...) FALSE

# -- quadratic analysis -------------------------------------------
method(is_quadratic, QuadForm) <- function(x) TRUE
method(has_quadratic_term, QuadForm) <- function(x) TRUE
method(is_pwl, QuadForm) <- function(x) FALSE

# -- get_data -----------------------------------------------------
method(get_data, QuadForm) <- function(x) list()

# -- numeric ------------------------------------------------------
method(numeric_value, QuadForm) <- function(x, values, ...) {
  xv <- values[[1L]]
  Pv <- values[[2L]]
  ## Use Hermitian form x^H P x for complex, crossprod for real
  ## R's crossprod(xv, y) = t(xv) %*% y, NOT Conj(t(xv)) %*% y
  if (is.complex(xv) || is.complex(Pv)) {
    val <- Conj(t(xv)) %*% Pv %*% xv
    matrix(Re(as.vector(val)), 1L, 1L)
  } else {
    matrix(as.numeric(crossprod(xv, Pv %*% xv)), 1L, 1L)
  }
}

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

# -- decomp_quad: eigendecomposition for conic canonicalization ----
## CVXPY SOURCE: atoms/quad_form.py lines 186-251
## Returns list(scale, M1, M2) where P = scale * (M1 %*% t(M1) - M2 %*% t(M2))
decomp_quad <- function(P, cond = NULL) {
  P <- as.matrix(P)
  eig <- .eigvalsh(P, only_values = FALSE)
  w <- eig$values
  V <- eig$vectors

  if (is.null(cond)) {
    cond <- 1e6 * .Machine$double.eps
  }

  scale <- max(abs(w))
  if (scale == 0) {
    w_scaled <- w
  } else {
    w_scaled <- w / scale
  }

  maskp <- w_scaled > cond
  maskn <- w_scaled < -cond

  if (any(maskp) && any(maskn)) {
    cli_warn("Forming a nonconvex expression {.fn quad_form}(x, indefinite).")
  }

  ## Scale each column of V by sqrt of corresponding eigenvalue
  M1 <- if (any(maskp)) {
    t(t(V[, maskp, drop = FALSE]) * sqrt(w_scaled[maskp]))
  } else {
    matrix(numeric(0), nrow = nrow(P), ncol = 0L)
  }

  M2 <- if (any(maskn)) {
    t(t(V[, maskn, drop = FALSE]) * sqrt(-w_scaled[maskn]))
  } else {
    matrix(numeric(0), nrow = nrow(P), ncol = 0L)
  }

  list(scale = scale, M1 = M1, M2 = M2)
}

#' Quadratic form x^T P x
#'
#' When \code{x} is constant, returns \code{t(Conj(x)) \%*\% P \%*\% x}
#' (affine in \code{P}).
#' When \code{P} is constant, returns a \code{QuadForm} atom (quadratic in \code{x}).
#' At least one of \code{x} or \code{P} must be constant.
#'
#' @param x An Expression (vector)
#' @param P An Expression (square matrix, symmetric/Hermitian)
#' @param assume_PSD If TRUE, assume P is PSD without checking (only when P is constant).
#' @returns A QuadForm atom or an affine Expression
#' @export
quad_form <- function(x, P, assume_PSD = FALSE) {
  ## CVXPY SOURCE: quad_form.py lines 254-275
  x <- as_expr(x)
  P <- as_expr(P)
  ## Dimension checks
  if (length(P@shape) != 2L || P@shape[1L] != P@shape[2L]) {
    cli_abort("Invalid dimensions for arguments to {.fn quad_form}: P must be square.")
  }
  n <- P@shape[1L]
  if (x@shape[1L] != n) {
    cli_abort("Invalid dimensions for arguments to {.fn quad_form}: x has {x@shape[1L]} rows, P has {n}.")
  }
  if (is_constant(x)) {
    ## x constant: x^H P x is affine in P
    Conj(t(x)) %*% P %*% x
  } else if (is_constant(P)) {
    if (assume_PSD) {
      P <- psd_wrap(P)
    }
    QuadForm(x, P)
  } else {
    cli_abort("At least one argument to {.fn quad_form} must be non-variable.")
  }
}