R/112_atoms_tr_inv.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/tr_inv.py
## TrInv -- trace of the inverse of a PSD matrix: tr(X^{-1})


TrInv <- new_class("TrInv", parent = Atom, package = "CVXR",
  constructor = function(X, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    X <- as_expr(X)
    ## Shape is always scalar
    shape <- c(1L, 1L)

    obj <- .fast_new(TrInv, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(X),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- validate -----------------------------------------------------
## CVXPY: tr_inv.py lines 52-57
method(validate_arguments, TrInv) <- function(x) {
  A <- .args(x)[[1L]]
  if (length(.shape(A)) != 2L || .shape(A)[1L] != .shape(A)[2L]) {
    cli_abort("The argument to {.fn tr_inv} must be a square matrix, got shape ({A@shape[1L]}, {A@shape[2L]}).")
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
## CVXPY: tr_inv.py lines 59-62 -- returns tuple()
method(shape_from_args, TrInv) <- function(x) c(1L, 1L)

# -- sign ---------------------------------------------------------
## CVXPY: tr_inv.py lines 64-67 -- (True, False)
method(sign_from_args, TrInv) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
## CVXPY: tr_inv.py lines 69-77 -- convex, not concave
method(is_atom_convex, TrInv) <- function(x) TRUE
method(is_atom_concave, TrInv) <- function(x) FALSE

# -- monotonicity -------------------------------------------------
## CVXPY: tr_inv.py lines 79-87 -- not monotone
method(is_incr, TrInv) <- function(x, idx, ...) FALSE
method(is_decr, TrInv) <- function(x, idx, ...) FALSE

# -- numeric ------------------------------------------------------
## CVXPY: tr_inv.py lines 37-50
method(numeric_value, TrInv) <- function(x, values, ...) {
  A <- values[[1L]]
  ## Check symmetry
  if (norm(A - t(A), "F") >= 1e-8) {
    return(matrix(Inf, 1L, 1L))
  }
  ## Take symmetric part for numerical stability
  symm <- (A + t(A)) / 2
  eig_vals <- .eigvalsh(symm, only_values = TRUE)$values
  if (min(eig_vals) <= 0) {
    return(matrix(Inf, 1L, 1L))
  }
  matrix(sum(1 / eig_vals), 1L, 1L)
}

# -- domain -------------------------------------------------------
## CVXPY: tr_inv.py lines 111-114 -- X >> 0
method(domain, TrInv) <- function(x) {
  list(PSD(.args(x)[[1L]]))
}

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

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

# -- .grad: subgradient of tr(X^-1) -------------------------------
## CVXPY SOURCE: atoms/tr_inv.py:89-108 (tr_inv._grad).
## d/dX tr(X^-1) = -X^{-T} X^{-T} = -(X^{-T})^2.
## NULL if min eigenvalue <= 0 (outside the open PSD cone).
method(.grad, TrInv) <- function(x, values, ...) {
  X <- as.matrix(values[[1L]])
  ev <- .eigvalsh(X, only_values = TRUE)$values
  if (min(ev) <= 0) return(list(NULL))
  D <- t(solve(X))
  D <- -D %*% D
  rows <- as.integer(prod(.arg_shape(x)))
  list(.dense_to_csc_vector(as.numeric(D), rows))
}

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

#' Trace of matrix inverse
#'
#' Computes \eqn{\mathrm{tr}(X^{-1})} for PSD matrix X.
#'
#' @param X A square PSD matrix expression
#' @returns An expression representing \eqn{\mathrm{tr}(X^{-1})}
#' @export
tr_inv <- function(X) {
  TrInv(X)
}