R/103_atoms_matrix_frac.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/matrix_frac.py
## MatrixFrac -- trace(X^T * P^{-1} * X)
##
## Also provides matrix_frac() convenience function with QuadForm shortcut.


MatrixFrac <- new_class("MatrixFrac", 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
    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 -----------------------------------------------------
## CVXPY: matrix_frac.py lines 87-99
method(validate_arguments, MatrixFrac) <- function(x) {
  X <- x@args[[1L]]
  P <- x@args[[2L]]
  if (length(P@shape) != 2L || P@shape[1L] != P@shape[2L]) {
    cli_abort("The second argument to {.fn matrix_frac} must be a square matrix.")
  }
  if (X@shape[1L] != P@shape[1L]) {
    cli_abort("The arguments to {.fn matrix_frac} have incompatible dimensions.")
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
## CVXPY: matrix_frac.py lines 101-104 -- returns tuple()
method(shape_from_args, MatrixFrac) <- function(x) c(1L, 1L)

# -- sign ---------------------------------------------------------
## CVXPY: matrix_frac.py lines 106-109 -- (True, False)
method(sign_from_args, MatrixFrac) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
## CVXPY: matrix_frac.py lines 111-119 -- convex, not concave
method(is_atom_convex, MatrixFrac) <- function(x) TRUE
method(is_atom_concave, MatrixFrac) <- function(x) FALSE

# -- monotonicity -------------------------------------------------
## CVXPY: matrix_frac.py lines 121-129 -- not monotone
method(is_incr, MatrixFrac) <- function(x, idx, ...) FALSE
method(is_decr, MatrixFrac) <- function(x, idx, ...) FALSE

# -- quadratic overrides ------------------------------------------
## CVXPY: matrix_frac.py lines 131-144
method(is_quadratic, MatrixFrac) <- function(x) {
  is_affine(x@args[[1L]]) && is_constant(x@args[[2L]])
}

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

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

# -- numeric ------------------------------------------------------
## CVXPY: matrix_frac.py lines 36-46
method(numeric_value, MatrixFrac) <- function(x, values, ...) {
  X <- values[[1L]]
  P <- values[[2L]]
  ## Hermitian form: X^H P^{-1} X for complex
  product <- if (is.complex(X) || is.complex(P)) {
    Conj(t(X)) %*% solve(P) %*% X
  } else {
    t(X) %*% solve(P) %*% X
  }
  if (length(dim(product)) == 2L) {
    matrix(sum(diag(product)), 1L, 1L)
  } else {
    matrix(product, 1L, 1L)
  }
}

# -- domain -------------------------------------------------------
## CVXPY: matrix_frac.py lines 48-51 -- P >> 0
method(domain, MatrixFrac) <- function(x) {
  list(PSD(x@args[[2L]]))
}

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

# -- graph_implementation -----------------------------------------
method(graph_implementation, MatrixFrac) <- function(x, arg_objs, shape, data = NULL, ...) {
  cli_abort("graph_implementation for {.cls MatrixFrac} not available; use Dcp2Cone canonicalization.")
}

# ==================================================================
# Convenience function
# ==================================================================

#' Matrix fractional function
#'
#' Computes \eqn{\mathrm{trace}(X^T P^{-1} X)}. If P is a constant matrix, uses
#' a QuadForm shortcut for efficiency.
#'
#' @param X A matrix expression (n by m)
#' @param P A square matrix expression (n by n), must be PSD
#' @returns An expression representing \eqn{\mathrm{trace}(X^T P^{-1} X)}
#' @export
matrix_frac <- function(X, P) {
  ## CVXPY: matrix_frac.py lines 147-153
  ## If P is a constant numeric matrix, shortcut via QuadForm
  if (is.matrix(P) && !inherits(P, "Expression") && !S7_inherits(P, Expression)) {
    invP <- solve(P)
    return(quad_form(X, (invP + t(invP)) / 2.0))
  }
  MatrixFrac(X, P)
}