R/106_atoms_log_det.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/log_det.py
## LogDet -- log-determinant of a PSD matrix: log(det(A))


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

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

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

# -- shape --------------------------------------------------------
## CVXPY: log_det.py lines 56-59 -- returns tuple()
method(shape_from_args, LogDet) <- function(x) c(1L, 1L)

# -- sign ---------------------------------------------------------
## CVXPY: log_det.py lines 61-64 -- (True, False)
## Note: This is CVXPY's convention. Mathematically log_det can be negative
## for 0 < det < 1, but CVXPY reports is_nonneg=TRUE. We replicate.
method(sign_from_args, LogDet) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
## CVXPY: log_det.py lines 66-73 -- NOT convex, IS concave
method(is_atom_convex, LogDet) <- function(x) FALSE
method(is_atom_concave, LogDet) <- function(x) TRUE

# -- monotonicity -------------------------------------------------
## CVXPY: log_det.py lines 75-83 -- not monotone
method(is_incr, LogDet) <- function(x, idx, ...) FALSE
method(is_decr, LogDet) <- function(x, idx, ...) FALSE

# -- numeric ------------------------------------------------------
## CVXPY: log_det.py lines 36-48
## CVXPY uses np.linalg.slogdet which handles complex Hermitian matrices.
## R's determinant() does NOT support complex matrices, so we use
## eigenvalue decomposition: log(det(A)) = sum(log(eigenvalues(A)))
## for Hermitian/symmetric PSD matrices.
method(numeric_value, LogDet) <- function(x, values, ...) {
  A <- values[[1L]]
  ## Take Hermitian part for numerical stability (conj transpose, not just transpose)
  symm <- (A + Conj(t(A))) / 2
  if (is.complex(symm)) {
    ## For complex Hermitian matrices, eigenvalues are real
    ev <- Re(eigen(symm, symmetric = FALSE, only.values = TRUE)$values)
    if (all(ev > 0)) {
      matrix(sum(log(ev)), 1L, 1L)
    } else {
      matrix(-Inf, 1L, 1L)
    }
  } else {
    ## Real symmetric case: use determinant() for numerical stability
    det_result <- determinant(symm, logarithm = TRUE)
    if (isTRUE(det_result$sign == 1)) {
      matrix(as.numeric(det_result$modulus), 1L, 1L)
    } else {
      matrix(-Inf, 1L, 1L)
    }
  }
}

# -- domain -------------------------------------------------------
## CVXPY: log_det.py lines 107-110 -- A >> 0
method(domain, LogDet) <- function(x) {
  list(PSD(x@args[[1L]]))
}

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

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

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

#' Log-determinant
#'
#' Computes log(det(A)) for PSD matrix A.
#'
#' @param A A square PSD matrix expression
#' @returns An expression representing log(det(A))
#' @export
log_det <- function(A) {
  LogDet(A)
}