R/114_atoms_lambda_sum_largest.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/lambda_sum_largest.py
## LambdaSumLargest -- sum of the k largest eigenvalues of a symmetric matrix
##
## Inherits from LambdaMax conceptually (in CVXPY).
## In S7, no parent delegation, so we replicate structure.
##
## Also provides lambda_sum_smallest() pure function: -lambda_sum_largest(-X, k)


LambdaSumLargest <- new_class("LambdaSumLargest", parent = Atom, package = "CVXR",
  properties = list(
    ## Numeric, not integer: CVXPY accepts a fractional k (lambda_sum_largest.py:40-54,
    ## 63-81; test_atoms.py::test_sum_largest "float k"). `as.integer(1.5)` silently
    ## gave k = 1 (claims-audit row C0214, repro notes/audit/repro/19_*).
    k = class_numeric
  ),
  constructor = function(A, k, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    A <- as_expr(A)
    k <- as.numeric(k)
    ## Shape is always scalar
    shape <- c(1L, 1L)

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

# -- validate -----------------------------------------------------
## CVXPY: lambda_sum_largest.py lines 31-38 (k > 0 is the only requirement on k)
method(validate_arguments, LambdaSumLargest) <- function(x) {
  A <- .args(x)[[1L]]
  if (length(.shape(A)) != 2L || .shape(A)[1L] != .shape(A)[2L]) {
    cli_abort("First argument to {.fn lambda_sum_largest} must be a square matrix.")
  }
  if (length(x@k) != 1L || is.na(x@k) || x@k <= 0) {
    cli_abort("Second argument to {.fn lambda_sum_largest} must be a positive number.")
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
method(shape_from_args, LambdaSumLargest) <- function(x) c(1L, 1L)

# -- sign ---------------------------------------------------------
## Same as lambda_max: (False, False)
method(sign_from_args, LambdaSumLargest) <- function(x) {
  list(is_nonneg = FALSE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
## Convex (like lambda_max)
method(is_atom_convex, LambdaSumLargest) <- function(x) TRUE
method(is_atom_concave, LambdaSumLargest) <- function(x) FALSE

# -- monotonicity -------------------------------------------------
method(is_incr, LambdaSumLargest) <- function(x, idx, ...) FALSE
method(is_decr, LambdaSumLargest) <- function(x, idx, ...) FALSE

# -- numeric ------------------------------------------------------
## CVXPY: lambda_sum_largest.py lines 40-54. Sum of the floor(k) largest
## eigenvalues plus (k - floor(k)) times the next one when k is fractional.
## .eigvalsh returns values in DECREASING order (CVXPY's eigvalsh is ascending
## and indexes from the end), so the k largest are the first k here.
method(numeric_value, LambdaSumLargest) <- function(x, values, ...) {
  A <- values[[1L]]
  evals <- .eigvalsh(A, only_values = TRUE)$values
  k_floor <- floor(x@k)
  k_frac <- x@k - k_floor
  result <- sum(evals[seq_len(min(k_floor, length(evals)))])
  if (k_frac > 0 && k_floor < length(evals)) {
    result <- result + k_frac * evals[k_floor + 1L]
  }
  matrix(result, 1L, 1L)
}

# -- get_data -----------------------------------------------------
## CVXPY: lambda_sum_largest.py lines 56-59 -- returns [self.k]
method(get_data, LambdaSumLargest) <- function(x) list(x@k)

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

# -- .grad: subgradient of the sum of the k largest eigenvalues --------
## CVXPY SOURCE: atoms/lambda_sum_largest.py:63-81 (_grad_matrices) composed with
## the inherited lambda_max._grad (atoms/lambda_max.py:62-84), real since CVXPY
## 1.9.0 (1.8.2 raised NotImplementedError; CVXR returned NULL until 1.9.2.9135,
## claims-audit row C0214). D = V_top V_top^T + k_frac * v_next v_next^T where
## V_top holds the eigenvectors of the floor(k) largest eigenvalues and v_next the
## next one. CVXPY uses LA.eigh (ascending) and takes the last columns; .eigvalsh
## is decreasing so the same vectors are the first columns. D is raveled
## column-major into the (n*n, 1) Jacobian column, as in .grad(LambdaMax).
method(.grad, LambdaSumLargest) <- function(x, values, ...) {
  A <- as.matrix(values[[1L]])
  n <- ncol(A)
  ev <- .eigvalsh(A, only_values = FALSE)
  k_floor <- floor(x@k)
  k_frac <- x@k - k_floor
  D <- matrix(0, n, n)
  if (k_floor > 0) {
    V_top <- ev$vectors[, seq_len(min(k_floor, n)), drop = FALSE]
    D <- D + V_top %*% t(V_top)
  }
  if (k_frac > 0 && k_floor < n) {
    v_next <- ev$vectors[, k_floor + 1L, drop = FALSE]
    D <- D + k_frac * (v_next %*% t(v_next))
  }
  rows <- as.integer(prod(.arg_shape(x)))
  list(.dense_to_csc_vector(as.numeric(D), rows))
}

# ==================================================================
# Convenience functions
# ==================================================================

#' Sum of largest k eigenvalues
#' @param A A square matrix expression
#' @param k Number of largest eigenvalues to sum (positive number; a fractional
#'   part weights the next eigenvalue, as in CVXPY)
#' @returns An expression representing the sum of the k largest eigenvalues
#' @export
lambda_sum_largest <- function(A, k) {
  LambdaSumLargest(A, k)
}

#' Sum of smallest k eigenvalues
#' @param A A square matrix expression
#' @param k Number of smallest eigenvalues to sum (positive number; a fractional
#'   part weights the next eigenvalue, as in CVXPY)
#' @returns An expression representing the sum of the k smallest eigenvalues
#' @export
lambda_sum_smallest <- function(A, k) {
  -lambda_sum_largest(-as_expr(A), k)
}