R/103_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 (FALSE) new_object(S7_object())  ## S7 static-check guard
    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 <- .fast_new(QuadForm, 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 <- .args(x)[[2L]]
  xarg <- .args(x)[[1L]]
  ## P must be square
  if (.shape(P)[1L] != .shape(P)[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 <- .shape(P)[1L]
  if (.shape(xarg)[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.
## CVXPY v1.9.0 #3142: inside quad_form_dpp_scope (QP-solver path), allow a
## param-affine P provided x is param-free -- enables DPP for quad_form(x, P)
## with a parametric P. See quad_form.py:62-99.
##
## .check_dpp_args mirrors quad_form.py _check_dpp_args(): x param-free
## (avoid quadratic-in-params) and P param-affine (DPP requirement).
.quad_form_check_dpp_args <- function(x) {
  is_param_free(.args(x)[[1L]]) && is_param_affine(.args(x)[[2L]])
}

method(is_atom_convex, QuadForm) <- function(x) {
  P <- .args(x)[[2L]]
  if (quad_form_dpp_scope_active()) {
    return(.quad_form_check_dpp_args(x) && is_psd(P))
  }
  is_constant(P) && is_psd(P)
}

method(is_atom_concave, QuadForm) <- function(x) {
  P <- .args(x)[[2L]]
  if (quad_form_dpp_scope_active()) {
    return(.quad_form_check_dpp_args(x) && is_nsd(P))
  }
  is_constant(P) && is_nsd(P)
}

## CVXPY quad_form.py: quad_form is smooth (a polynomial).
method(is_atom_smooth, QuadForm) <- function(x) TRUE

# -- 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 -------------------------------------------------
## CVXPY v1.8.2 fix: per-argument monotonicity (previously returned FALSE
## unconditionally, losing DCP composition information).
## idx=1 -> x (1-based; CVXPY idx=0), idx=2 -> P (CVXPY idx=1)
method(is_incr, QuadForm) <- function(x, idx, ...) {
  if (idx == 1L) {
    ## nabla_x f = 2Px: nonneg when (x>=0, P>=0) or (x<=0, P<=0)
    (is_nonneg(.args(x)[[1L]]) && is_nonneg(.args(x)[[2L]])) ||
      (is_nonpos(.args(x)[[1L]]) && is_nonpos(.args(x)[[2L]]))
  } else if (idx == 2L) {
    ## d f / d P_ij = x_i * x_j: nonneg when x has definite sign
    is_nonneg(.args(x)[[1L]]) || is_nonpos(.args(x)[[1L]])
  } else {
    FALSE
  }
}
method(is_decr, QuadForm) <- function(x, idx, ...) {
  if (idx == 1L) {
    ## nabla_x f = 2Px: nonpos when (x>=0, P<=0) or (x<=0, P>=0)
    (is_nonneg(.args(x)[[1L]]) && is_nonpos(.args(x)[[2L]])) ||
      (is_nonpos(.args(x)[[1L]]) && is_nonneg(.args(x)[[2L]]))
  } else {
    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.")
}

# -- .grad: per-atom subgradient ----------------------------------
## CVXPY SOURCE: atoms/quad_form.py:133-144 (QuadForm._grad).
## d/dx (x^T P x) = (P + P^T) x  (the symmetrised quadratic gradient).
## CVXPY returns a 1-element list [D]; the second arg (P) is always a
## Constant/Parameter in practice, whose grad is empty so the walker's
## per-arg loop never accesses grad_self[[2]]. We mirror the 1-element
## return exactly to keep parity.
method(.grad, QuadForm) <- function(x, values, ...) {
  xv <- as.numeric(values[[1L]])
  P  <- as.matrix(values[[2L]])
  D  <- (P + Conj(t(P))) %*% xv
  rows <- as.integer(prod(.arg_shape(x)))
  list(.dense_to_csc_vector(as.numeric(D), rows))
}

## NOTE: CVXPY also has SymbolicQuadForm in quad_form.py:147-172, whose
## _grad raises NotImplementedError. CVXR does not (yet) have a
## SymbolicQuadForm counterpart, so no port is needed here.

# -- 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(.shape(P)) != 2L || .shape(P)[1L] != .shape(P)[2L]) {
    cli_abort("Invalid dimensions for arguments to {.fn quad_form}: P must be square.")
  }
  n <- .shape(P)[1L]
  if (.shape(x)[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.")
  }
}