R/067_atoms_elementwise_huber.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/elementwise/huber.py
## Huber -- Huber loss function
##
## huber(x, M) = 2M|x| - M^2 if |x| >= M
##              |x|^2          if |x| <= M
## M is a parameter/constant stored as a property, NOT in args.


Huber <- new_class("Huber", parent = Elementwise, package = "CVXR",
  properties = list(
    M = class_any  # Constant or Parameter
  ),
  constructor = function(x, M = 1, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    x <- as_expr(x)
    ## M must be scalar, nonneg, and constant or parameter
    M_expr <- as_expr(M)
    if (!expr_is_scalar(M_expr)) {
      cli_abort("{.arg M} must be a scalar, got shape {.val {M_expr@shape}}.")
    }
    shape <- .shape(x)

    obj <- .fast_new(Huber, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(x),
      shape = shape,
      M     = M_expr
    )
    validate_arguments(obj)
    obj
  }
)

# -- parameters ---------------------------------------------------
## CVXPY: huber.py:52-54 — override needed because `M` is stored as a
## property, not in `args`. Without this, a parametric `M` would be
## invisible to parameters(problem) and miss out on chain-rule gradients.
method(parameters, Huber) <- function(x) {
  unique_list(c(
    unlist(lapply(.args(x), parameters), recursive = FALSE),
    parameters(x@M)
  ))
}

# -- validate: M must be nonneg, scalar, constant (CVXPY huber.py lines 94-98)
method(validate_arguments, Huber) <- function(x) {
  M <- x@M
  if (!(is_nonneg(M) && expr_is_scalar(M) && is_constant(M))) {
    cli_abort("{.arg M} must be a non-negative scalar constant or {.cls Parameter}.")
  }
  invisible(NULL)
}

# -- sign: always nonneg ----------------------------------------
method(sign_from_args, Huber) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

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

# -- monotonicity -------------------------------------------------
method(is_incr, Huber) <- function(x, idx, ...) is_nonneg(.args(x)[[idx]])
method(is_decr, Huber) <- function(x, idx, ...) is_nonpos(.args(x)[[idx]])

# -- quadratic ----------------------------------------------------
method(is_quadratic, Huber) <- function(x) {
  is_affine(.args(x)[[1L]])
}

method(has_quadratic_term, Huber) <- function(x) TRUE

# -- get_data -----------------------------------------------------
method(get_data, Huber) <- function(x) {
  list(x@M)
}

# -- numeric: 2 * SciPy Huber -------------------------------------
## CVXPY: 2 * scipy.special.huber(M, x)
## scipy.special.huber(delta, r) = { r^2  if |r| <= delta
##                                 { 2*delta*|r| - delta^2  if |r| > delta
## So CVXPY: 2*scipy.huber = our huber definition.
method(numeric_value, Huber) <- function(x, values, ...) {
  v <- values[[1L]]
  M_val <- as.numeric(value(x@M))
  ## huber(x, M) = x^2 if |x| <= M, else 2M|x| - M^2
  abs_v <- abs(v)
  ifelse(abs_v <= M_val, v^2, 2 * M_val * abs_v - M_val^2)
}

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

# -- .grad: per-atom subgradient ----------------------------------
## CVXPY SOURCE: atoms/elementwise/huber.py:100-118 (huber._grad).
## huber(x, M) is quadratic for |x| <= M and linear for |x| > M.
## Subgradient: 2 * sign(x) * min(|x|, M), elementwise diagonal.
method(.grad, Huber) <- function(x, values, ...) {
  rows <- as.integer(prod(.arg_shape(x)))
  cols <- as.integer(prod(.shape(x)))
  v <- values[[1L]]
  M_val <- as.numeric(value(x@M))
  min_val <- pmin(abs(v), M_val)
  grad_vals <- 2 * sign(v) * min_val
  list(.elemwise_grad_to_diag(grad_vals, rows, cols))
}

#' Create a Huber loss atom
#'
#' @param x An Expression
#' @param M Numeric threshold (default 1)
#' @returns A Huber atom
#' @export
huber <- function(x, M = 1) {
  Huber(x, M)
}