R/spautor_matrixInv.R

In spmodel: Spatial Statistical Modeling and Prediction

#' Compute the covariance matrix inverse (and log-determinant) for autoregressive models
#'
#' @param spcov_params_val A \code{spcov_params} object
#' @param data_object The data object
#' @param dist_matrix_list The neighbor weight matrix \code{W} (despite the name,
#'   this is not actually a list for \code{spautor()}/\code{spgautor()} models;
#'   named for consistency with the analogous \code{splm()}/\code{spglm()} code)
#' @param randcov_params_val A \code{randcov_params} object
#' @param ldet Whether to also return the log-determinant
#'
#' @return A list with elements \code{SigInv} (the inverse covariance matrix)
#'   and \code{Sigldet} (its log-determinant, or \code{NULL} if \code{ldet} is
#'   \code{FALSE}), computed via the Sherman-Morrison-Woodbury identity when
#'   there is a nonzero independent error variance, since directly inverting
#'   the dense \code{de + ie} covariance matrix would discard the sparsity of
#'   the CAR/SAR precision matrix
#'
#' @noRd
spautor_cov_matrixInv <- function(spcov_params_val, data_object,
                                  dist_matrix_list, randcov_params_val, ldet = TRUE) {
  # floor de away from zero -- during optimization de can be proposed near 0,
  # which would make the CAR/SAR precision matrix (and its inverse) singular
  if (spcov_params_val[["de"]] < 0.001) {
    spcov_params_val[["de"]] <- 0.001
  }

  # MAKE IT CLEAR DIST_MATRIX_LIST IS NOT A LIST
  dist_matrix <- dist_matrix_list

  if (is.null(data_object$partition_matrix)) {
    # compute the dependent inverse portion
    SigInv_de <- spcov_matrixInv_de(spcov_params_val, dist_matrix, data_object$M)
    ldet_Sig_de <- -as.numeric(Matrix::determinant(SigInv_de)$modulus)
    # if there is no independent random error variance, stop with the inverse
    # and log determinant -- else find inverse and log determinant of
    # (de + ie)
    if (spcov_params_val[["ie"]] == 0) {
      # if the de is zero this is the inverse
      SigInv <- SigInv_de
      # if the de is zero this is the ldet
      ldet_Sig <- ldet_Sig_de
    } else {
      # finding the cholesky and inverse of middle smw (A^-1 + B^-1)^-1
      ## FOR RAND EFFECTS NEED TO SMW STEP THROUGH Z %*% I %*% Zt
      smw_mid <- SigInv_de
      diag(smw_mid) <- 1 / spcov_params_val[["ie"]] + diag(smw_mid)
      smw_mid_upchol <- chol(forceSymmetric(smw_mid))
      Inv_smw_mid <- chol2inv(smw_mid_upchol)
      ldet_smw_mid <- 2 * sum(log(diag(smw_mid_upchol)))
      # finding the inverse (A + B)^-1 = A^-1 - A^-1 (A^-1 + B^-1)^-1 A^-1
      SigInv <- SigInv_de - SigInv_de %*% Inv_smw_mid %*% SigInv_de
      # finding ldet (A + B)^-1 = A^-1 - A^-1 (A^-1 + B^-1)^-1 A^-1
      ldet_Sig <- ldet_Sig_de + data_object$n * log(spcov_params_val[["ie"]]) + ldet_smw_mid
    }

    # do random effects
    smwInv_rand_val <- smwInv_rand(SigInv, ldet_Sig, randcov_params_val, data_object$randcov_Zs)

    # HW steps
    hwInv_val <- hwInv(smwInv_rand_val$SigInv, smwInv_rand_val$Sigldet, data_object$observed_index)

    SigInv <- hwInv_val$SigInv
    Sigldet <- hwInv_val$Sigldet
  } else {
    # a partition factor multiplies the covariance element-wise by a 0/1
    # matrix (Hadamard product), which destroys the CAR/SAR sparsity pattern
    # the SMW shortcut above relies on -- fall back to building the dense
    # covariance matrix directly and inverting it via Cholesky
    # making a covariance matrix
    cov_matrix_val_full <- cov_matrix(
      spcov_params_val, dist_matrix, randcov_params_val,
      data_object$randcov_Zs, data_object$partition_matrix, data_object$M
    )

    # subsetting for only observed
    cov_matrix_val <- cov_matrix_val_full[data_object$observed_index, data_object$observed_index, drop = FALSE]

    ## find cholesky
    chol_cov_matrix_val <- chol(forceSymmetric(cov_matrix_val))

    ## store value
    SigInv <- chol2inv(chol_cov_matrix_val)
    Sigldet <- 2 * sum(log(diag(chol_cov_matrix_val)))
  }

  # store value
  if (ldet) {
    return(list(SigInv = SigInv, Sigldet = Sigldet))
  } else {
    return(list(SigInv = SigInv, Sigldet = NULL))
  }
}