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,34 @@ assignOperatorDef(
   )
 )
 
+assignOperatorDef(
+  'Cholesky',
+  list(
+    labelAbstractTypes = list(
+      handler = 'Cholesky'
+    )
+  )
+)
+
+## TODO: remove these next two as redundant with Cholesky.
+assignOperatorDef(
+  'sparseChol',
+  list(
+    labelAbstractTypes = list(
+      handler = 'sparseChol'
+    )
+  )
+)
+
+assignOperatorDef(
+  'denseChol',
+  list(
+    labelAbstractTypes = list(
+      handler = 'denseChol'
+    )
+  )
+)
+
 assignOperatorDef(
   c('nChol'),
   list(
@@ -1120,7 +1148,7 @@ assignOperatorDef(
   c('nLogdet'),
   list(
     labelAbstractTypes = list(
-      handler = 'UnaryReduction'
+      handler = 'nLogdet'
     )
   )
 )

nCompiler/R/compile_labelAbstractTypes.R

@@ -521,14 +521,67 @@ inLabelAbstractTypesEnv(
   }
 )
 
+inLabelAbstractTypesEnv(
+  Cholesky <- function(code, symTab, auxEnv, handlingInfo) {
+    inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+    argType <- code$args[[1]]$type
+    if(inherits(argType, 'symbolSparse')) {
+        code$type <- symbolNC$new(
+          name = code$name, type = 'sparseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE
+      )
+    } else { 
+        code$type <- symbolNC$new(
+          name = code$name, type = 'denseCholFactor', NCgenerator = denseCholFactor, isArg = FALSE
+      )
+    }
+    invisible(inserts)
+  }
+)
+
+## TODO: remove these next two as redundant with Cholesky
+inLabelAbstractTypesEnv(
+  sparseChol <- function(code, symTab, auxEnv, handlingInfo) {
+    inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+    argType <- code$args[[1]]$type
+    if(inherits(argType, 'symbolSparse')) {
+        code$type <- symbolNC$new(
+          name = code$name, type = 'sparseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE
+      )
+    } else { 
+        code$type <- symbolNC$new(
+          name = code$name, type = 'denseCholFactor', NCgenerator = denseCholFactor, isArg = FALSE
+      )
+    }
+    invisible(inserts)
+  }
+)
+
+## Akin to `Matrix::Cholesky` in returning collection of information for further computation.
+inLabelAbstractTypesEnv(
+  denseChol <- function(code, symTab, auxEnv, handlingInfo) {
+    inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo)
+    argType <- code$args[[1]]$type
+    if(inherits(argType, 'symbolSparse')) {
+        code$type <- symbolNC$new(
+          name = code$name, type = 'sparseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE
+      )
+    } else { 
+        code$type <- symbolNC$new(
+          name = code$name, type = 'denseCholFactor', NCgenerator = denseCholFactor, isArg = FALSE
+      )
+    }
+    invisible(inserts)
+  }
+)
+
+## Akin to R `chol` in returning U matrix.
 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 {
       # Cholesky factor of a dense matrix is a dense matrix (i.e., same type)
       type <- setReturnType(handlingInfo, argType$type)
@@ -1477,7 +1530,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 +1655,37 @@ 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
+      code$type <- symbolLLT$new(name = code$name)
+    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,47 @@ 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", 
+                          templateArgs = list("Eigen::SparseMatrix<double>"))
+    }
+  )
+)
+
+symbolLLT <- R6::R6Class(
+  classname = "symbolLLT",
+  inherit = symbolBase,
+  portable = TRUE,
+  public = list(
+    interface = FALSE,  
+    initialize = function(...) {
+      super$initialize(..., interface = self$interface)
+      self$type <- "LLT"
+      self  
+    },
+    print = function() {
+      writeLines(paste0(self$name, ': symbolLLT'))
+    },
+    genCppVar = function() {
+      cppVarFullClass$new(name = self$name,
+                          baseType = "Eigen::LLT", 
+                          templateArgs = list("Eigen::MatrixXd"))  # Eigen::Matrix<double>?
+    }
+  )
+)

nCompiler/R/typeDeclarations.R

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

nCompiler/R/zzz_NC_Predefined.R

@@ -125,3 +125,31 @@ OptimResultList <- nClass(
     message = 'RcppCharacterVector'
   )
 )
+
+#' @export
+sparseCholFactor <- nClass(
+  # manually insert:
+  #include <nCompiler/ET_ext/post_Rcpp/tensorOperations_sparseChol.h>
+  # in the hContent file AFTER the sparseCholFactor class declaration
+    classname = 'sparseCholFactor',
+    predefined = quote(system.file(file.path("include","nCompiler", "predef"), package="nCompiler") |>
+                       file.path("sparseCholFactor_nC")),
+    Cpublic = list(
+        llt = 'simplicialLLT'
+    ) 
+)
+
+#' @export
+denseCholFactor <- nClass(
+  # manually insert:
+  #include <nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h>
+  # in the hContent file AFTER the denseCholFactor class declaration
+    classname = 'denseCholFactor',
+    predefined = quote(system.file(file.path("include","nCompiler", "predef"), package="nCompiler") |>
+                       file.path("denseCholFactor_nC")),
+    Cpublic = list(
+        llt = 'LLT'
+    ) 
+)
+
+