R/is_sym_pos_def.R

Defines functions is_sym_pos_def

Documented in is_sym_pos_def

# ---------------------------------------------------------------------------- #

utils::globalVariables(c("fmt_mtx_rows"))

# ---------------------------------------------------------------------------- #
#' Check if a matrix is symmetric positive-definite
#'
#' Checks a matrix to determine if it is symmetric positive-definite.
#'
#' @param x A numeric matrix which should be checked to see if it is symmetric
#'   positive-definite.
#'
#' @return TRUE if the supplied matrix is symmetric positive-definite, FALSE
#'   otherwise.
#'
#' @examples
#' is_sym_pos_def(matrix(c(2.51, 2.01, 2.01, 1.74), 2, 2)) # => TRUE
#' is_sym_pos_def(matrix(c(1, 1, 1, 1), 2, 2)) # => FALSE
#'
#' @export
#'
is_sym_pos_def <- function(x) {

  # Note cli::cli_abort collapses adjacent spaces, making output hard to read -
  # is there a better way to display this in the event of an error?

  format_mtx_rows <- function(x) c("Matrix values:", capture.output(x))

  if (!test_matrix(x)) {
    cli_abort(c(
      "{.var x} must be a matrix",
      "x" = "You've supplied a {.cls {typeof(x)}}."
    ),
    class = "jute_error"
    )
  }

  if (!test_numeric(x)) {
    msg_lines <- format_mtx_rows(x)
    cli_abort(c(
      "{.var x} must be numeric",
      "x" = "You've supplied a {.cls {typeof(x)}} matrix.",
      setNames(msg_lines, rep("i", length(msg_lines)))
    ),
    class = "jute_error"
    )
  }

  # Use (base) S3 method base::isSymmetric to check for symmetry
  if (!isSymmetric(x)) {
    msg_lines <- format_mtx_rows(x)
    cli_abort(c(
      "{.var x} must be symmetric",
      "x" = "You've supplied a non-symmetric matrix.",
      setNames(msg_lines, rep("i", length(msg_lines)))
    ),
    class = "jute_error"
    )
  }

  # If a matrix is *symmetric* and positive-definite, then all of its
  # eigenvalues are positive
  if (length(which(eigen(x)$values <= 0.0)) > 0) {
    FALSE
  } else {
    TRUE
  }
}

# ---------------------------------------------------------------------------- #
toniprice/jute documentation built on Jan. 11, 2023, 8:23 a.m.