R/025_interface_matrix_utilities.R

In CVXR: Disciplined Convex Optimization

#####
## DO NOT EDIT THIS FILE!! EDIT THE SOURCE INSTEAD: rsrc_tree/interface/matrix_utilities.R
#####

## CVXPY SOURCE: interface/matrix_utilities.py
## R-side matrix interface utilities

# -- Convert values to canonical R matrix form -------------------------
## CVXPY equivalent: intf.convert() / DEFAULT_INTF.const_to_matrix()

#' Convert a value to a numeric matrix or sparse matrix
#'
#' Normalizes R values so the rest of CVXR can assume a consistent type, as
#' CVXPY's `Constant()` does through `const_to_matrix` (constant.py:51-54;
#' numpy_interface/sparse_matrix_interface.py:40-51, ndarray_interface.py):
#' every sparse input becomes a double CSC matrix, bool and int are cast to
#' double, the shape is kept and explicit zeros are kept.
#'
#' Here any `sparseMatrix` or `sparseVector` becomes a general double
#' `dgCMatrix` (symmetric, triangular, diagonal, triplet, row-compressed,
#' pattern and logical classes included); a dense `Matrix` becomes a base
#' matrix; a base logical, integer, double or complex value is stored as
#' double (complex kept), a scalar or vector as a column.
#'
#' @param val A logical, integer, double or complex scalar, vector or matrix,
#'   or a Matrix object
#' @returns A base matrix or a dgCMatrix
#' @keywords internal
intf_convert <- function(val) {
  if (inherits(val, "sparseVector")) {
    ## Matrix's own coercion to the n x 1 sparse column (a pattern vector gives
    ## an ngCMatrix); normalized by the sparseMatrix branch below.
    val <- methods::as(val, "CsparseMatrix")
  }
  if (inherits(val, "sparseMatrix")) {
    ## General double CSC. The three-step chain is Matrix's own replacement
    ## for the deprecated as(., "dgCMatrix") on logical and pattern classes;
    ## it expands one-triangle (dsC, dtC) and diagonal storage, materializes
    ## a unit diagonal's ones, and sums duplicate triplets. No drop0: CVXPY
    ## keeps explicit zeros.
    if (!inherits(val, "dgCMatrix")) {
      val <- as(as(as(val, "dMatrix"), "generalMatrix"), "CsparseMatrix")
    }
    return(val)
  }
  if (inherits(val, "Matrix")) {
    ## Dense Matrix class (dgeMatrix, lgeMatrix, ...): a base matrix, then
    ## the storage cast below.
    val <- as.matrix(val)
  }
  if (!(is.numeric(val) || is.logical(val) || is.complex(val)) || is.null(val)) {
    cli_abort("Cannot convert object of class {.cls {class(val)}} to a CVXR matrix.")
  }
  ## Logical and integer become double; `storage.mode<-` keeps the dim.
  if (!is.complex(val)) storage.mode(val) <- "double"
  d <- dim(val)
  if (is.null(d) || length(d) == 1L) {
    ## Scalar -> 1 x 1, vector or 1-d array -> column (Fortran order, as CVXPY).
    dim(val) <- c(length(val), 1L)
  } else if (length(d) != 2L) {
    cli_abort("Cannot convert an array with {length(d)} dimensions to a CVXR matrix; constants are at most two-dimensional.")
  }
  val
}

# -- Shape of a value --------------------------------------------------
## CVXPY equivalent: intf.shape()

#' Get the shape of a value as an integer vector c(nrow, ncol)
#' @param val A numeric scalar, vector, matrix, or Matrix object
#' @returns Integer vector of length 2
#' @keywords internal
intf_shape <- function(val) {
  if (is.matrix(val) || inherits(val, "Matrix")) {
    return(c(nrow(val), ncol(val)))
  }
  if (length(val) == 1L) {
    return(c(1L, 1L))
  }
  ## Vector
  c(length(val), 1L)
}

# -- Sign of a value --------------------------------------------------
## CVXPY equivalent: intf.sign()

#' Determine the sign of a numeric value
#' @param val A numeric matrix or sparse matrix
#' @returns List with \code{is_nonneg} (logical) and \code{is_nonpos} (logical)
#' @keywords internal
intf_sign <- function(val) {
  ## Handle empty values (0-element)
  if (length(val) == 0L) {
    return(list(is_nonneg = TRUE, is_nonpos = TRUE))
  }
  ## Complex values have no sign ordering
  ## CVXPY SOURCE: constant.py _compute_attr() sets nonneg=nonpos=False for complex
  if (is.complex(val)) {
    return(list(is_nonneg = FALSE, is_nonpos = FALSE))
  }
  if (inherits(val, "sparseMatrix")) {
    ## Pattern matrices (ngCMatrix etc.) have no @x slot -- treat as all zeros
    if (!.hasSlot(val, "x") || length(val@x) == 0L) {
      return(list(is_nonneg = TRUE, is_nonpos = TRUE))
    }
    vals <- val@x
  } else {
    vals <- as.numeric(val)
  }
  ## NaN -> treat as unknown sign (neither nonneg nor nonpos)
  if (anyNA(vals)) {
    return(list(is_nonneg = FALSE, is_nonpos = FALSE))
  }
  mx <- max(vals)
  mn <- min(vals)
  list(is_nonneg = (mn >= 0), is_nonpos = (mx <= 0))
}

# -- Is sparse? --------------------------------------------------------
## CVXPY equivalent: intf.is_sparse()

#' Check if a value is a sparse matrix
#' @param val A value
#' @returns Logical
#' @keywords internal
intf_is_sparse <- function(val) {
  ## sp.issparse: a sparseVector is sparse too (intf_convert stores it as a
  ## sparse column), so Constant()'s .sparse flag matches its stored value.
  inherits(val, "sparseMatrix") || inherits(val, "sparseVector")
}

# -- Hermitian / symmetric check ---------------------------------------
## CVXPY SOURCE: interface/matrix_utilities.py:276-299 is_hermitian (CVXPY 1.9.3)
##
## Tolerance is numpy's `np.allclose` (rtol 1e-5, atol 1e-8, FALSE on NaN),
## through `.np_allclose` (zzz_R_specific/utility.R).  Through CVXR 1.9.2.9154
## these tests used `all.equal` / `isSymmetric` at 1e-10, a MEAN relative
## difference that rejected what CVXPY accepts (an off-diagonal pair differing
## by 1e-7 made quad_form abort), accepted what CVXPY rejects (a 1e12 pair off
## by 1e-3 was skew-symmetric), and called a NaN matrix symmetric so that
## is_psd crashed in eigen() (claims audit C0487, decision 22, 1.9.2.9155).
##
## Dense matrices are compared entry by entry.  Before transposing, rows 1, 2,
## n-1 and n are compared with the matching columns, as base R's
## isSymmetric.matrix does: those comparisons are a subset of the full test,
## so the answer is unchanged, and an asymmetric 1000 x 1000 matrix is
## rejected in microseconds instead of after a 16 MB transpose.

#' Check if a matrix is symmetric (and Hermitian for real case)
#' @param val A numeric matrix or sparse matrix
#' @returns List with \code{is_symmetric} (logical) and \code{is_hermitian} (logical)
#' @keywords internal
intf_is_hermitian <- function(val) {
  if (is.complex(val)) {
    ## "symmetric" always means "real symmetric" in CVXPY (the comment at
    ## matrix_utilities.py:281-287), so is_symmetric is FALSE for complex.
    is_symm <- FALSE
    if (inherits(val, "sparseMatrix")) {
      is_herm <- intf_is_sparse_symmetric(val, complex = TRUE)
    } else {
      is_herm <- .dense_allclose_transpose(val, Conj)
    }
    return(list(is_symmetric = is_symm, is_hermitian = is_herm))
  }
  if (inherits(val, "sparseMatrix")) {
    is_symm <- intf_is_sparse_symmetric(val, complex = FALSE)
  } else {
    is_symm <- .dense_allclose_transpose(val, identity)
  }
  list(is_symmetric = is_symm, is_hermitian = is_symm)
}

## `np.allclose(val, cj(val.T))` for a dense matrix, with the row/column
## pre-check described above.  `cj` is Conj for the Hermitian test and
## identity for the symmetric one.
.dense_allclose_transpose <- function(val, cj) {
  if (inherits(val, "Matrix")) val <- as.matrix(val)
  n <- nrow(val)
  if (n > 1L) {
    for (i in unique(c(1L, 2L, n - 1L, n))) {
      if (!.np_allclose(val[i, ], cj(val[, i]))) return(FALSE)
    }
  }
  .np_allclose(val, cj(t(val)))
}

## CVXPY SOURCE: interface/matrix_utilities.py:315-366 is_sparse_symmetric
##
## Structure first, values second: after summing duplicates and dropping
## explicit zeros, the strict lower and strict upper triangles must hold the
## same positions (mirrored), and then the paired values must be allclose
## (conjugated for the Hermitian test).  A stored 1e-9 facing an absent entry
## is therefore NOT symmetric here even though the dense test would accept it,
## exactly as upstream.  The complex branch is unreachable in R today (Matrix
## has no complex sparse class) and is kept for isomorphism.
##
## Same test, no sorting.  Upstream lexsorts the two triangles and compares
## them.  A canonical compressed-column matrix (row indices increasing within
## each column, no duplicates, no explicit zeros) is a sorted position list
## already, and so is its transpose, so "the lower positions mirror the upper
## positions" is `identical(pattern(m), pattern(t(m)))`, and entry k of
## `t(m)@x` is the mirror of entry k of `m@x`.  The value comparison is
## restricted to the strict lower triangle so that the `rtol * |b|` side is
## the upper mirror, as in upstream's `np.allclose(vl, vu)`.  Measured on a
## 1000 x 1000 matrix with 2e4 nonzeros: the lexsort form took 411 us / 1.8 MB
## against 144 us for the old Matrix::isSymmetric; this form is in the table
## under "Decision 22" in notes/audit/claims_audit_item1.md.

#' Check if a sparse matrix is symmetric (or Hermitian)
#' @param m A square sparse matrix
#' @param complex If \code{TRUE}, test \code{m == Conj(t(m))}
#' @returns Logical scalar
#' @keywords internal
intf_is_sparse_symmetric <- function(m, complex = FALSE) {
  if (nrow(m) != ncol(m)) cli_abort("m must be a square matrix")
  m <- .csc_canonical(m)
  mt <- Matrix::t(m)
  ## Entry (i, j) in the lower triangle must mirror entry (j, i) in the upper.
  if (!(identical(m@p, mt@p) && identical(m@i, mt@i))) return(FALSE)
  lower <- m@i > rep.int(seq_len(ncol(m)) - 1L, diff(m@p))   # strict lower triangle
  vl <- m@x[lower]
  vu <- mt@x[lower]
  if (complex) .np_allclose(vl, Conj(vu)) else .np_allclose(vl, vu)
}

## `sp.coo_array(m); m.sum_duplicates(); m.eliminate_zeros()` for a Matrix
## object: a general dgCMatrix with duplicates summed and explicit zeros
## dropped.  Symmetric and diagonal storage classes keep one triangle only,
## so the matrix is expanded to a general one first; a dgCMatrix is already
## canonical up to explicit zeros.
.csc_canonical <- function(m) {
  if (!inherits(m, "dgCMatrix")) {
    m <- as(as(as(m, "dMatrix"), "generalMatrix"), "CsparseMatrix")   # sums duplicates
  }
  if (any(m@x == 0)) m <- Matrix::drop0(m)
  m
}

# -- Skew-symmetric check ---------------------------------------------
## CVXPY SOURCE: interface/matrix_utilities.py:302-312 is_skew_symmetric
##
## `np.allclose(constant + constant.T, 0.0)`: the reference is 0, so only
## atol = 1e-8 applies.  Same row/column pre-check as the symmetric test.

#' Check if a matrix is skew-symmetric (A + A^T == 0)
#' @param val A numeric matrix or sparse matrix
#' @returns Logical scalar
#' @keywords internal
intf_is_skew_symmetric <- function(val) {
  ## Non-square: error (matches CVXPY)
  if (nrow(val) != ncol(val)) {
    cli_abort("Skew-symmetry requires a square matrix.")
  }
  ## Complex matrices: skew-symmetric not defined in CVXPY (returns False)
  if (is.complex(val)) return(FALSE)
  if (inherits(val, "sparseMatrix")) return(intf_is_sparse_skew_symmetric(val))
  if (inherits(val, "Matrix")) val <- as.matrix(val)
  n <- nrow(val)
  if (n > 1L) {
    for (i in unique(c(1L, 2L, n - 1L, n))) {
      if (!.np_allclose(val[i, ] + val[, i], 0)) return(FALSE)
    }
  }
  .np_allclose(val + t(val), 0)
}

## CVXPY SOURCE: interface/matrix_utilities.py:369-413 is_sparse_skew_symmetric
##
## As is_sparse_symmetric, with the diagonal included on both sides (a
## diagonal entry mirrors itself, so it must be allclose to 0 after doubling)
## and the paired values summed against 0.  Same pattern-of-the-transpose
## form as above.

#' Check if a real sparse matrix satisfies A + t(A) == 0
#' @param A A square sparse matrix
#' @returns Logical scalar
#' @keywords internal
intf_is_sparse_skew_symmetric <- function(A) {
  if (nrow(A) != ncol(A)) cli_abort("m must be a square matrix")
  A <- .csc_canonical(A)
  At <- Matrix::t(A)
  ## Entry (i, j) with i >= j must mirror entry (j, i) with j <= i.
  if (!(identical(A@p, At@p) && identical(A@i, At@i))) return(FALSE)
  lower <- A@i >= rep.int(seq_len(ncol(A)) - 1L, diff(A@p))   # lower triangle with diagonal
  .np_allclose(A@x[lower] + At@x[lower], 0)
}

# -- PSD check via eigenvalues -----------------------------------------
## CVXPY equivalent: utilities/linalg.py::is_psd_within_tol

#' Check if a symmetric matrix is PSD within tolerance
#' @param val A symmetric numeric matrix
#' @param tol Eigenvalue tolerance
#' @returns Logical
#' @keywords internal
intf_is_psd <- function(val, tol = EIGVAL_TOL) {
  ev <- .eigvalsh(val, only_values = TRUE)$values
  all(ev >= -tol)
}