R/111_atoms_matrix_frac.R

In CVXR: Disciplined Convex Optimization

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

## CVXPY SOURCE: atoms/matrix_frac.py
## MatrixFrac -- trace(X^T * P^{-1} * X)
##
## Also provides matrix_frac() convenience function with QuadForm shortcut.


MatrixFrac <- new_class("MatrixFrac", parent = Atom, package = "CVXR",
  constructor = function(X, P, id = NULL) {
    if (FALSE) new_object(S7_object())  ## S7 static-check guard
    if (is.null(id)) id <- next_expr_id()
    X <- as_expr(X)
    P <- as_expr(P)
    ## Shape is always scalar
    shape <- c(1L, 1L)

    obj <- .fast_new(MatrixFrac, S7_object(),
      id    = as.integer(id),
      .cache = new.env(parent = emptyenv()),
      args  = list(X, P),
      shape = shape
    )
    validate_arguments(obj)
    obj
  }
)

# -- validate -----------------------------------------------------
## CVXPY: matrix_frac.py lines 87-99
method(validate_arguments, MatrixFrac) <- function(x) {
  X <- .args(x)[[1L]]
  P <- .args(x)[[2L]]
  if (length(.shape(P)) != 2L || .shape(P)[1L] != .shape(P)[2L]) {
    cli_abort("The second argument to {.fn matrix_frac} must be a square matrix.")
  }
  if (.shape(X)[1L] != .shape(P)[1L]) {
    cli_abort("The arguments to {.fn matrix_frac} have incompatible dimensions.")
  }
  invisible(NULL)
}

# -- shape --------------------------------------------------------
## CVXPY: matrix_frac.py lines 101-104 -- returns tuple()
method(shape_from_args, MatrixFrac) <- function(x) c(1L, 1L)

# -- sign ---------------------------------------------------------
## CVXPY: matrix_frac.py lines 106-109 -- (True, False)
method(sign_from_args, MatrixFrac) <- function(x) {
  list(is_nonneg = TRUE, is_nonpos = FALSE)
}

# -- curvature ----------------------------------------------------
## CVXPY: matrix_frac.py lines 111-119 -- convex, not concave
method(is_atom_convex, MatrixFrac) <- function(x) TRUE
method(is_atom_concave, MatrixFrac) <- function(x) FALSE

# -- monotonicity -------------------------------------------------
## CVXPY: matrix_frac.py lines 121-129 -- not monotone
method(is_incr, MatrixFrac) <- function(x, idx, ...) FALSE
method(is_decr, MatrixFrac) <- function(x, idx, ...) FALSE

# -- quadratic overrides ------------------------------------------
## CVXPY: matrix_frac.py lines 131-144
method(is_quadratic, MatrixFrac) <- function(x) {
  is_affine(.args(x)[[1L]]) && is_constant(.args(x)[[2L]])
}

method(has_quadratic_term, MatrixFrac) <- function(x) {
  is_constant(.args(x)[[2L]])
}

method(is_qpwa, MatrixFrac) <- function(x) {
  is_pwl(.args(x)[[1L]]) && is_constant(.args(x)[[2L]])
}

# -- .mf_inv: the inverse CVXPY computes here ----------------------
## Every `LA.inv(P)` in matrix_frac.py -- numeric (:42,44), _grad (:68) and the
## `matrix_frac()` constructor shortcut (:149) -- is `numpy.linalg.inv`, which
## has NO conditioning test: it goes straight to LAPACK and raises only on
## EXACT singularity.
##
## Base R's `solve()` is different. It applies `tol = .Machine$double.eps`
## against the reciprocal condition number and refuses anything below it with
## "system is computationally singular". For a merely ILL-CONDITIONED P the two
## languages therefore disagree about whether a value exists at all, and CVXR
## raised where CVXPY returns a number.
##
## Found by revdepcheck, not by the suite: SLSEdesign's vignette computes a
## c-optimal design whose optimum is DEGENERATE -- the information matrix is
## rank-deficient in the directions the criterion does not constrain -- so the
## better the solver converges, the closer P gets to singular. At the optimum
## rcond(P) = 4.3e-18 and `solve(P)` refused, which aborted
## `value(problem@objective)` inside psolve() and took the whole solve with it.
##
## `tol = 0` restores numpy's behavior, and the numbers agree exactly:
##     R    solve(P, diag(n), tol = 0)  ->  matrix_frac = -5.23946e-14
##     numpy np.linalg.inv(P)           ->  matrix_frac = -5.2394597073595286e-14
## which is the right answer (mathematically 0, approached from below). CVXR
## 1.9.1 did not error on this problem only because it stopped further from the
## optimum, where it returned -4.08e-05 -- also negative, and four orders worse.
##
## DELIBERATE ADDITION over upstream: numpy inverts an ill-conditioned matrix
## silently, and the result can be numerically meaningless. R's refusal existed
## for a reason, so rather than simply deleting the check CVXR keeps the value
## (parity) and WARNS when the conditioning is below what base R would have
## accepted. The user gets the number CVXPY would give, plus the information R
## would have given, instead of one at the cost of the other.
##
## Returns NULL when P is exactly singular, which is `LA.LinAlgError`.
.mf_inv <- function(P, warn = TRUE) {
  n <- nrow(P)
  inv <- tryCatch(solve(P, diag(n), tol = 0), error = function(e) NULL)
  ## Warn only when a value is actually being RETURNED. An exactly singular P
  ## returns NULL and the caller reports that; warning "ill-conditioned" on the
  ## way to a "singular" error would be two messages for one problem.
  if (!is.null(inv) && isTRUE(warn)) {
    rc <- tryCatch(rcond(P), error = function(e) NA_real_)
    if (!is.na(rc) && rc < .Machine$double.eps) {
      cli_warn(c(
        "The second argument to {.fn matrix_frac} is ill-conditioned \\
         (reciprocal condition number {format(rc, digits = 3)}).",
        "i" = "Its inverse, and any value computed from it, may be numerically meaningless.",
        "i" = "Base {.fn solve} would refuse this matrix; the value is still returned, \\
               matching {.pkg CVXPY}."
      ), class = "matrixFracConditioning")
    }
  }
  inv
}

# -- numeric ------------------------------------------------------
## CVXPY: matrix_frac.py lines 36-46
method(numeric_value, MatrixFrac) <- function(x, values, ...) {
  X <- values[[1L]]
  P <- values[[2L]]
  P_inv <- .mf_inv(as.matrix(P))
  if (is.null(P_inv)) {
    ## Exactly singular. CVXPY's LA.inv raises LinAlgError here too; upstream
    ## has never decided what the value should be (matrix_frac.py:38 is a
    ## standing "# TODO raise error if not invertible?"), so keep the error.
    cli_abort(c(
      "Cannot evaluate {.fn matrix_frac}: the second argument is singular.",
      "i" = "{.fn matrix_frac} is {.code trace(t(X) %*% solve(P) %*% X)}, which requires an invertible {.arg P}."
    ))
  }
  ## Hermitian form: X^H P^{-1} X for complex
  product <- if (is.complex(X) || is.complex(P)) {
    Conj(t(X)) %*% P_inv %*% X
  } else {
    t(X) %*% P_inv %*% X
  }
  if (length(dim(product)) == 2L) {
    matrix(sum(diag(product)), 1L, 1L)
  } else {
    matrix(product, 1L, 1L)
  }
}

# -- domain -------------------------------------------------------
## CVXPY: matrix_frac.py lines 48-51 -- P >> 0
method(domain, MatrixFrac) <- function(x) {
  list(PSD(.args(x)[[2L]]))
}

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

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

# -- .grad: subgradient of matrix_frac(X, P) = X^T P^-1 X ---------
## CVXPY SOURCE: atoms/matrix_frac.py:53-91 (matrix_frac._grad).
## Two args: X (vector / matrix) and P (PSD matrix).
##   d/dX = (P^-1 + P^-T) X
##   d/dP = - (P^-1 X X^T P^-1)^T   (already symmetrised in CVXPY)
## NULL on both args when P is singular.
method(.grad, MatrixFrac) <- function(x, values, ...) {
  X <- as.matrix(values[[1L]])
  P <- as.matrix(values[[2L]])
  ## CVXPY SOURCE: matrix_frac.py:67-70 -- `try: P_inv = LA.inv(P) except
  ## LA.LinAlgError: return [None, None]`. `.mf_inv` is that `LA.inv`; see its
  ## header for why base `solve()` is not.
  ##
  ## `warn = FALSE`: the conditioning warning belongs to the VALUE, which is
  ## computed from the same P in the same solve and warns there. Warning again
  ## here would double every message, and upstream's _grad is silent.
  P_inv <- .mf_inv(P, warn = FALSE)
  if (is.null(P_inv)) return(list(NULL, NULL))

  ## d/dX = (P^-1 + P^-T) X
  DX <- (P_inv + t(P_inv)) %*% X
  rows_X <- as.integer(prod(.arg_shape(x)))

  ## d/dP = -(P^-1 X X^T P^-1)^T  (CVXPY's exact form, see lines 80-83)
  DP <- P_inv %*% X
  DP <- DP %*% t(X)
  DP <- DP %*% P_inv
  DP <- -t(DP)
  ## CVXPY then transposes again before raveling: t(DP).ravel(order='F')
  ## i.e. emits the *original* (untransposed) entries in column-major.
  rows_P <- as.integer(prod(.arg_shape(x, 2L)))

  list(
    .dense_to_csc_vector(as.numeric(DX), rows_X),
    .dense_to_csc_vector(as.numeric(t(DP)), rows_P)
  )
}

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

#' Matrix fractional function
#'
#' Computes \eqn{\mathrm{trace}(X^T P^{-1} X)}. If P is a constant matrix, uses
#' a QuadForm shortcut for efficiency.
#'
#' @param X A matrix expression (n by m)
#' @param P A square matrix expression (n by n), must be PSD
#' @returns An expression representing \eqn{\mathrm{trace}(X^T P^{-1} X)}
#' @export
matrix_frac <- function(X, P) {
  ## CVXPY: matrix_frac.py lines 147-153
  ## If P is a constant numeric matrix, shortcut via QuadForm
  if (is.matrix(P) && !inherits(P, "Expression") && !.s7_is(P, Expression)) {
    ## CVXPY SOURCE: matrix_frac.py:148-150 -- `invP = LA.inv(P)`, same
    ## conditioning semantics as the other two sites. See `.mf_inv`.
    invP <- .mf_inv(P)
    if (is.null(invP)) {
      cli_abort(c(
        "The second argument to {.fn matrix_frac} is singular.",
        "i" = "A constant {.arg P} is inverted eagerly to build a {.fn quad_form}."
      ))
    }
    return(quad_form(X, (invP + t(Conj(invP))) / 2.0))
  }
  MatrixFrac(X, P)
}