R/071_atoms_elementwise_rel_entr.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/elementwise/rel_entr.py
## RelEntr -- elementwise relative entropy: x*log(x/y)


RelEntr <- new_class("RelEntr", parent = Elementwise, package = "CVXR",
  constructor = function(x, y, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    y <- as_expr(y)
    shape <- sum_shapes(list(.shape(x), .shape(y)))
    obj <- .fast_new(RelEntr, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(x, y),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- sign: unknown (can be positive or negative) -----------------
method(sign_from_args, RelEntr) <- function(x) {
  list(is_nonneg = FALSE, is_nonpos = FALSE)
}

# -- curvature: convex --------------------------------------------
method(is_atom_convex, RelEntr) <- function(x) TRUE
method(is_atom_concave, RelEntr) <- function(x) FALSE

## CVXPY elementwise/rel_entr.py: rel_entr is smooth.
method(is_atom_smooth, RelEntr) <- function(x) TRUE

# -- monotonicity: not incr in either; decr in arg 2 -------------
method(is_incr, RelEntr) <- function(x, idx, ...) FALSE
method(is_decr, RelEntr) <- function(x, idx, ...) {
  idx == 2L  # decreasing in second argument (1-indexed: idx==2)
}

# -- domain: x >= 0, y >= 0 --------------------------------------
method(atom_domain, RelEntr) <- function(x) {
  list(.args(x)[[1L]] >= 0, .args(x)[[2L]] >= 0)
}

# -- numeric: x*log(x/y) -----------------------------------------
## Matches scipy.special.rel_entr(x, y)
method(numeric_value, RelEntr) <- function(x, values, ...) {
  xv <- as.matrix(values[[1L]])
  yv <- as.matrix(values[[2L]])
  ## Broadcast scalar to match shape (R doesn't broadcast ifelse/&)
  if (length(yv) == 1L && length(xv) > 1L) yv <- matrix(yv, nrow(xv), ncol(xv))
  if (length(xv) == 1L && length(yv) > 1L) xv <- matrix(xv, nrow(yv), ncol(yv))
  ## rel_entr(x, y) = x*log(x/y) when x > 0, y > 0
  ## rel_entr(0, y) = 0 when y >= 0
  ## rel_entr(x, 0) = Inf when x > 0
  result <- ifelse(xv > 0 & yv > 0,
                   xv * log(xv / yv),
                   ifelse(xv == 0 & yv >= 0, 0, Inf))
  result
}

# -- graph_implementation: stub -----------------------------------
method(graph_implementation, RelEntr) <- function(x, arg_objs, shape, data = NULL, ...) {
  cli_abort("graph_implementation for {.cls RelEntr} not yet implemented.")
}

# -- .grad: per-atom subgradient ----------------------------------
## CVXPY SOURCE: atoms/elementwise/rel_entr.py:69-93 (rel_entr._grad).
## rel_entr(x, y) = x * log(x/y); domain x > 0, y > 0.
## d/dx = log(x/y) + 1; d/dy = -x/y.
method(.grad, RelEntr) <- function(x, values, ...) {
  v0 <- values[[1L]]; v1 <- values[[2L]]
  if (min(v0) <= 0 || min(v1) <= 0) return(list(NULL, NULL))
  div <- v0 / v1
  cols <- as.integer(prod(.shape(x)))
  rows0 <- as.integer(prod(.arg_shape(x)))
  rows1 <- as.integer(prod(.arg_shape(x, 2L)))
  list(
    .elemwise_grad_to_diag(log(div) + 1, rows0, cols),
    .elemwise_grad_to_diag(-div,         rows1, cols)
  )
}

#' Relative Entropy: x*log(x/y)
#'
#' @param x An Expression
#' @param y An Expression
#' @returns A RelEntr atom
#' @export
rel_entr <- function(x, y) {
  RelEntr(x, y)
}