R/114_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 (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    A <- as_expr(A)
    ## Shape is always scalar
    shape <- c(1L, 1L)

    obj <- .fast_new(LogDet, 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 <- .args(x)[[1L]]
  if (length(.shape(A)) != 2L || .shape(A)[1L] != .shape(A)[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 v1.8.2 fix: log_det can be negative (log(det(A)) < 0 when 0 < det(A) < 1).
## Prior to v1.8.2, CVXPY incorrectly returned is_nonneg=TRUE.
method(sign_from_args, LogDet) <- function(x) {
  list(is_nonneg = FALSE, 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]]
  ## CVXPY v1.8.2 fix: guard against NULL (unset Variable)
  if (is.null(A)) return(NULL)
  ## 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(.args(x)[[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.")
}

# -- .grad: subgradient of log_det(X) -----------------------------
## CVXPY SOURCE: atoms/log_det.py:86-105 (log_det._grad).
## d/dX log det(X) = X^{-T} on the PSD interior. NULL if min eigenvalue
## <= 0 (outside the open PSD cone). Mirrors CVXPY's LA.eigvalsh
## domain check.
method(.grad, LogDet) <- 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))
  rows <- as.integer(prod(.arg_shape(x)))
  list(.dense_to_csc_vector(as.numeric(D), rows))
}

# ==================================================================
# 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)
}