add sparse cholesky functionality

nCompiler/R/Rexecution.R

@@ -196,6 +196,8 @@ makeReturnVector <- function(fillValue, length, recycle) {
 #' @export
 #' 
 nEigen <- function(x, symmetric, valuesOnly = FALSE) {
+  if(inherits(r, 'dCHMsimpl'))
+    stop("`nEigen` not implemented for sparse Cholesky factor input")      
   res <- eigen(x = x, symmetric = symmetric, only.values = valuesOnly)
   ans <- EigenDecomp$new()
   ans$values <- res$values
@@ -257,6 +259,8 @@ nEigen <- function(x, symmetric, valuesOnly = FALSE) {
 #'      return(singularValues)
 #'  })
 nSvd <- function(x, vectors = 'full') {
+  if(inherits(r, 'dCHMsimpl'))
+    stop("`nSvd` not implemented for sparse Cholesky factor input")      
   n <- nrow(x)
   p <- ncol(x)
   if(vectors == 'full') { # vectors = 2, when converted to int
@@ -284,6 +288,8 @@ nSvd <- function(x, vectors = 'full') {
 #' @export
 #' 
 nDiag <- function(x, ...) {
+  if(inherits(x, 'dCHMsimpl'))
+    stop("`nDiag` not implemented for sparse Cholesky factor input")      
   diag(x, ...)
 }
 
@@ -299,21 +305,42 @@ nDiag <- function(x, ...) {
 #' @export
 #' 
 nChol <- function(x) {
+  if(inherits(r, 'dCHMsimpl'))
+    stop("`nChol` not meaningful for sparse Cholesky factor input")      
   chol(x)
 }
 
+#' Compute the sparse Cholesky decomposition of a matrix for further operations
+#' 
+#' @details This function finds a sparse Cholesky factor using a fill-reducing
+#' `P` (AMD reordering by default), so the factorization is actually:
+# \eqn{P A P^\top = L L^\top}.
+#' 
+#' @param x a positive definite matrix
+#'
+#' @author nCompiler development team
+#'
+#' @export
+#'
+#' @details
+#'
+#' @examples TODO
+#' 
+sparseCholFactor <- function(x) {
+  Matrix::Cholesky(x, LDL = FALSE)  
+}
+
 #' Compute the log-determinant of a matrix
 #' 
 #' In a \code{nFunction}, \code{nLogdet} is identical to \code{logdet}
 #'
-#' @details This function is similar to R's \code{\link{diag}} function, but 
-#'   can be used in a nFunction and compiled using \code{nCompile}.  
-#' 
-#' @param x a square matrix
+#' @param x a square matrix or a sparse Cholesky factor
 #'
 #' @export
 #' 
 nLogdet <- function(x) {
+    if(inherits(x, 'dCHMsimpl'))
+       return(sum(log(diag(expand1(chol, "L")))))
     ldet <- determinant(x, logarithm = TRUE)
     ifelse(ldet$sign >= 0, ldet$modulus, NaN)
   }
@@ -410,27 +437,40 @@ asDense <- function(x) {
 #' 
 #' @export
 nMul <- function(x, y) {
+  if(inherits(x, "dCHMsimpl"))
+    # P^{top} L* x (see ?Matrix:::solve))
+    return(solve(ch, expand1(ch, "L") %*% x, system = "Pt"))  
   x %*% y
 }
 
+## This would only work if Matrix is loaded.
+## setMethod("%*%", c(x = "dCHMsimpl", y = "ANY"),
+##          function(x,y)
+##            solve(ch, expand1(ch, "L") %*% x, system = "Pt"))  # P^{top} L* x (see ?Matrix:::solve))
+
 #' Wrapper for solve
 #'
 #' @export
 nSolve <- function(a, b) {
-  solve(a,b)
+  solve(a,b)  # This works fine if `a` is of type `dCHMsimpl`.
 }
 
 #' Wrapper for forwardsolve
 #'
 #' @export
 nForwardsolve <- function(l, x) {
+  ## TODO: add error trap for non-sparse Cholesky factor input
+  if(inherits(r, 'dCHMsimpl'))
+    stop("`nForwardsolve` not implemented for sparse Cholesky factor input")      
   forwardsolve(l,x)
 }
 
 #' Wrapper for backsolve
 #'
 #' @export
 nBacksolve <- function(r, x) {
+  if(inherits(r, 'dCHMsimpl'))
+    return(solve(r, solve(ch, r, system = "Lt"), system = "Pt"))  # P^{top} U*^{-1} x
   backsolve(r,x)
 }
 

nCompiler/R/compile_aaa_operatorLists.R

@@ -1106,6 +1106,15 @@ assignOperatorDef(
   )
 )
 
+assignOperatorDef(
+  'sparseChol',
+  list(
+    labelAbstractTypes = list(
+      handler = 'sparseChol'
+    )
+  )
+)
+
 assignOperatorDef(
   c('nChol'),
   list(
@@ -1120,7 +1129,7 @@ assignOperatorDef(
   c('nLogdet'),
   list(
     labelAbstractTypes = list(
-      handler = 'UnaryReduction'
+      handler = 'nLogdet'
     )
   )
 )

nCompiler/R/compile_labelAbstractTypes.R

@@ -522,19 +522,38 @@ inLabelAbstractTypesEnv(
 )
 
 inLabelAbstractTypesEnv(
-  nChol <- function(code, symTab, auxEnv, handlingInfo) {
+  sparseChol <- function(code, symTab, auxEnv, handlingInfo) {
     inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
     argType <- code$args[[1]]$type
     if(inherits(argType, 'symbolSparse')) {
+      # Cholesky factor of a sparse matrix is a collection of information
+        code$type <- symbolNC$new(
+          name = code$name, type = 'sparseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE
+      )
+    } else {   # Fill in here when handle dense too.
+      # Cholesky factor of a dense matrix is a dense matrix (i.e., same type)
+      type <- setReturnType(handlingInfo, argType$type)
+      nDim <- setReturn_nDim(handlingInfo, argType$nDim)
+      code$type <- symbolBasic$new(type = type, nDim = nDim)
+    }
+    invisible(inserts)
+  }
+)
+
+inLabelAbstractTypesEnv(
+  nChol <- function(code, symTab, auxEnv, handlingInfo) {
+    inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+    argType <- code$args[[1]]$type
+    #if(inherits(argType, 'symbolSparse')) {
       # Cholesky factor of a sparse matrix is a collection of matrices
       # TODO: do we need to specify arguments for the initializer?
-      code$type <- symbolSparseCholesky$new(name = code$name)
-    } else {
+    #  code$type <- symbolSparseCholesky$new(name = code$name)
+    #} else {
       # Cholesky factor of a dense matrix is a dense matrix (i.e., same type)
       type <- setReturnType(handlingInfo, argType$type)
       nDim <- setReturn_nDim(handlingInfo, argType$nDim)
       code$type <- symbolBasic$new(type = type, nDim = nDim)
-    }
+    #}
     invisible(inserts)
   }
 )
@@ -1477,7 +1496,8 @@ inLabelAbstractTypesEnv(
     # product between two input vectors, which returns a matrix
     resDim <- 2
     # return a dense object if any arguments are dense, o/w return sparse
-    if(all(sapply(code$args, function(a) inherits(a$type, 'symbolSparse')))) {
+    if(all(sapply(code$args, function(a) inherits(a$type, 'symbolSparse'))) ||
+       inherits(code$args[[1]]$type, 'symbolSimplicialLLT') && inherits(code$args[[2]]$type, 'symbolSparse')) {
       code$type <- symbolSparse$new(nDim = resDim, type = returnType)
     } else {
       code$type <- symbolBasic$new(nDim = resDim, type = returnType)
@@ -1601,6 +1621,41 @@ inLabelAbstractTypesEnv(
   }
 )
 
+inLabelAbstractTypesEnv(
+    # nChol
+  chol <- function(code, symTab, auxEnv, handlingInfo) {
+    insertions <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+    argType <- code$args[[1]]$type
+    if(inherits(argType, 'symbolSparse')) {
+      code$type <- symbolSimplicialLLT$new(name = code$name)
+    } else
+        ## Put dense chol handling here
+      stop(exprClassProcessingErrorMsg(
+        code,
+        'sparse Cholesky factor only supported for sparse matrices.'
+      ), call. = FALSE)
+    invisible(NULL)
+  }
+)
+
+inLabelAbstractTypesEnv(
+  nLogdet <-
+    function(code, symTab, auxEnv, handlingInfo) {
+      if(length(code$args) != 1)
+        stop(exprClassProcessingErrorMsg(
+          code,
+          'nLogdet called with argument length != 1.'
+        ),
+        call. = FALSE)
+
+      inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+
+      code$type <- symbolBasic$new(nDim = 0,
+                                   type = setReturnType(handlingInfo, "double"))
+      inserts
+    }
+)
+
 inLabelAbstractTypesEnv(
   nEigen <- function(code, symTab, auxEnv, handlingInfo) {
     # determine object's natural type

nCompiler/R/symbolTable.R

@@ -629,3 +629,25 @@ symbolCppVar <- R6::R6Class(
     }
   )
 )
+
+symbolSimplicialLLT <- R6::R6Class(
+  classname = "symbolSimplicialLLT",
+  inherit = symbolBase,
+  portable = TRUE,
+  public = list(
+    interface = FALSE,  
+    initialize = function(...) {
+      super$initialize(..., interface = self$interface)
+      self$type <- "simplicialLLT"
+      self  
+    },
+    print = function() {
+      writeLines(paste0(self$name, ': symbolSimplicialLLT'))
+    },
+    genCppVar = function() {
+      cppVarFullClass$new(name = self$name,
+                          baseType = "Eigen::SimplicialLLT",  #  "std::shared_ptr",
+                          templateArgs = list("Eigen::SparseMatrix<double>"))
+    }
+  )
+)

nCompiler/R/typeDeclarations.R

@@ -268,6 +268,9 @@ typeDeclarationList <- list(
                            type = "double") {
     nSparseType(scalarType = type, nDim = 1)
   },
+  simplicialLLT = function() {
+    symbolSimplicialLLT$new()
+  },
   nList = function(type) {
     elementSym <- argType2symbol(type)
     symbolNlist$new(elementSym = elementSym)

nCompiler/R/zzz_NC_Predefined.R

@@ -125,3 +125,14 @@ OptimResultList <- nClass(
     message = 'RcppCharacterVector'
   )
 )
+
+#' @export
+sparseCholFactor <- nClass(
+    classname = 'sparseCholFactor',
+    predefined = quote(system.file(file.path("include","nCompiler", "predef"), package="nCompiler") |>
+                       file.path("sparseCholFactor_nC")),
+    Cpublic = list(
+        llt = 'simplicialLLT'
+    ) 
+)
+

nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h

@@ -211,19 +211,37 @@ Eigen::Tensor<Scalar, TensorExpr::NumDimensions> binaryOp(
     return binaryOp<OP_, Scalar>(xEval, y);
 }
 
-// forward declaration of nCompiler struct to store Sparse Chol. decompositions
-class SparseCholesky;
+// forward declaration of class to store sparse Cholesky decompositions
+// class SparseCholesky;
+class sparseCholFactor;
+
+// TODO: Not sure if IsSparseCholFactor stuff should be here or in tensorOperations_chol.h.
+// I think we need to use `IsSparseCholFactor` in `IsEvaluatedType`.
 
 /**
- * Template meta programming check to see if Class is an Eigen::SparseCholesky
+ * Template meta programming check to see if Class is an Eigen::SparseCholFactor
  *
  * @tparam Class type to inspect
  */
 template<typename Class>
-struct IsSparseCholesky : std::is_base_of<
-    SparseCholesky,
+struct IsSparseCholFactor : std::is_base_of<
+    std::shared_ptr<sparseCholFactor>,
     Class
 > { };
+
+/**
+  * Returns true if template Class has N dimensions
+  *
+  * Intended to be used as a template metaprogramming aid.
+  *
+  * @tparam Class Type to inspect, restricted to std::shared_ptr<sparseCholFactor> by SFINAE
+  * @tparam N Number of dimensions to test for
+  */
+template<typename Class, int N>
+constexpr typename std::enable_if<IsSparseCholFactor<Class>::value, bool>::type
+  HasNumDimensionsN() {
+    return N == 2;  // matrices are inherently 2-dimensional
+  }
 
 /**
  * Template meta programming check to see if Class is an Eigen::SparseMatrix
@@ -394,7 +412,7 @@ std::conditional<
 template<typename Class>
 struct IsEvaluatedType : std::conditional<
     IsSparseType<Class>::value || IsTensor<Class>::value ||
-    IsSparseCholesky<Class>::value || std::is_arithmetic<Class>::value ||
+    IsSparseCholFactor<Class>::value || std::is_arithmetic<Class>::value ||
     IsMap<Class>::value || IsTranspose<Class>::value,
     std::true_type,
     std::false_type
@@ -510,6 +528,7 @@ TENSOR_SPMAT_OP(!=, nCompiler::logical_neq)
      return N == 1;  // vectors are inherently 1-dimensional
    }
 
+
 /**
  * Implicitly convert a Tensor expression input to an Eigen::Tensor object
  *
@@ -747,6 +766,7 @@ Eigen::Tensor<Scalar, 2> asDense(SpMatExpr &x) {
     return res;
 }
 
+/*
 template<typename SparseCholType = SparseCholesky, typename Scalar>
 SparseCholType nChol(const Eigen::SparseMatrix<Scalar> &x) {
     Eigen::SimplicialLLT<Eigen::SparseMatrix<Scalar>> llt(x);
@@ -758,6 +778,7 @@ SparseCholType nChol(const Eigen::SparseMatrix<Scalar> &x) {
     res.P = P;
     return res;
 }
+*/
 
 /**
  * Compute the Cholesky decomposition for a symmetric matrix stored as an
@@ -965,7 +986,7 @@ Eigen::SparseMatrix<Scalar> nDiagonal(Xpr x, Index nrow, Index ncol) {
   * DiagIO class that uses conversion operators and overloaded
   * assignment operators to provide additional functionality (i.e., assignment
   * and extraction of diagonal entries) in a way that is compatible with
-  * nCompiler/R-like syntax.  Without the added flexibility, nCompiler would
+  * nCompiler/R-like syntax.  Without the added flexibility, nCompiler wonuld
   * need to perform additional, possibly complex, modification of an nFunction's
   * AST in order to change R-like syntax into Eigen-like syntax.
   *
@@ -1330,6 +1351,7 @@ Eigen::Tensor<typename RHS::Scalar, RHS::NumDimensions> nBacksolve(
     return triangularsolve<Eigen::Upper>(U, b);
 }
 
+
 /**
  * Matrix multiplication x %*% y when both inputs are matrix-like objects, i.e.,
  * rank 2 Eigen::Tensor objects, or Tensor expressions
@@ -1593,7 +1615,6 @@ Scalar nLogdet(const Xpr & x) {
     return qrdecomp.logDeterminant();
 }
 
-// TODO: implement nLogdet for sparse matrices
 
 // This is drafted but not yet used.
 template<typename Scalar >