R/decomp_cov.R

In jfrench/bayesutils: Utilities for Bayesian Analysis

#' Decompose covariance matrix
#'
#' \code{decomp_cov} decomposes a covariance matrix
#' \code{v}.  If \code{A = decomp_cov(v)}, then
#' \code{tcrossprod(A, A) == v}.
#'
#' The \code{"chol"} method is the fastest but least
#' stable method.  The \code{"eigen"} method is slower, but more
#' stable.  The \code{"svd"} method is the slowest method,
#' but should be the most stable.
#'
#' @param v An \eqn{N \times N} covariance matrix.
#' @param method The method used to decompose \code{v}.
#'   Valid options are \code{"chol"}, \code{"eigen"}, or
#'   \code{"svd"}.
#'
#' @return Returns an \eqn{N \times N} matrix.
#'
#' @author Joshua French
#' @export
#' @examples
#' # generate data
#' n = 100
#' coords = matrix(runif(n*2), nrow = n, ncol = 2)
#' d = as.matrix(dist(coords))
#' # create covariance matrix
#' v = 3 * exp(-d/2) + 0.1 * diag(n)
#'
#' # decompose v using the three methods
#' d1 = decomp_cov(v, "chol")
#' d2 = decomp_cov(v, "eigen")
#' d3 = decomp_cov(v, "svd")
#'
#' # verify accuracy of decompositions
#' all.equal(v, tcrossprod(d1))
#' all.equal(v, tcrossprod(d2), check.attributes = FALSE)
#' all.equal(v, tcrossprod(d3), check.attributes = FALSE)
decomp_cov <- function(v, method = "eigen") {
  #sanity check
  if (!(is.matrix(v)  || inherits(v, "Matrix"))) {
    stop("v should be a matrix or Matrix")
  }
  if (nrow(v) != ncol(v)) {
    stop("v must be a square numeric matrix")
  }
  if (!is.element(method, c("chol", "eigen", "svd"))) {
    stop("method must be 'chol', 'eigen', or 'svd'")
  }
  if (method == "eigen") {
    eigenv <- eigen(v)
    return(eigenv$vectors %*% diag(sqrt(pmax(eigenv$values,0))))
  } else if (method == "chol") {
    return(t(chol(v)))
  }else if (method == "svd") {
    svdv <- svd(v)
    return(tcrossprod(svdv$u %*% diag(sqrt(svdv$d)), svdv$v))
  }
}

#' @export
#' @rdname decomp_cov
decomp.cov = decomp_cov