From 2cc7511a3b62383859ae939dbe490fa5ecec9565 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sat, 2 May 2026 15:16:32 -0700 Subject: [PATCH 01/13] Make partial draft of sparse Cholesky, starting with factorization. --- nCompiler/R/Rexecution.R | 23 +++++++++++++++++++ nCompiler/R/compile_aaa_operatorLists.R | 9 ++++++++ nCompiler/R/compile_labelAbstractTypes.R | 23 +++++++++++++++---- nCompiler/R/symbolTable.R | 19 +++++++++++++++ .../ET_ext/post_Rcpp/tensorOperations.h | 12 ++++++++++ 5 files changed, 82 insertions(+), 4 deletions(-) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index ff7f3234..87b9c5fc 100644 --- a/nCompiler/R/Rexecution.R +++ b/nCompiler/R/Rexecution.R @@ -302,6 +302,29 @@ nChol <- function(x) { 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) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky +} + +## diag(Cholesky(x)) +## result <- solve(ch, b, system = "A"), or "L" + #' Compute the log-determinant of a matrix #' #' In a \code{nFunction}, \code{nLogdet} is identical to \code{logdet} diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index 0cac8c90..0c231939 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1106,6 +1106,15 @@ assignOperatorDef( ) ) +assignOperatorDef( + 'sparseCholFactor', + list( + labelAbstractTypes = list( + handler = 'sparseCholFactor' + ) + ) +) + assignOperatorDef( c('nChol'), list( diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 15d69c3b..d0cf0a73 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -525,16 +525,16 @@ inLabelAbstractTypesEnv( nChol <- function(code, symTab, auxEnv, handlingInfo) { inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) argType <- code$args[[1]]$type - if(inherits(argType, 'symbolSparse')) { + #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) } ) @@ -1601,6 +1601,21 @@ inLabelAbstractTypesEnv( } ) +inLabelAbstractTypesEnv( + sparseCholFactor <- 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 + stop(exprClassProcessingErrorMsg( + code, + 'sparse Cholesky factor only supported for sparse matrices.' + ), call. = FALSE) + invisible(NULL) + } +) + inLabelAbstractTypesEnv( nEigen <- function(code, symTab, auxEnv, handlingInfo) { # determine object's natural type diff --git a/nCompiler/R/symbolTable.R b/nCompiler/R/symbolTable.R index 5d22de7f..0d8b44c7 100644 --- a/nCompiler/R/symbolTable.R +++ b/nCompiler/R/symbolTable.R @@ -629,3 +629,22 @@ symbolCppVar <- R6::R6Class( } ) ) + +symbolSimplicialLLT <- R6::R6Class( + classname = "symbolSimplicialLLT", + inherit = symbolBase, + portable = TRUE, + public = list( + initialize = function(...) { + super$initialize(...) + self$type = "simplicialLLT" # "simplicialLLT_ptr" ? + self + }, + print = function() print(self$type), + genCppVar = function() { + cppVarFullClass$new(name = self$name, + baseType = "std::shared_ptr", + templateArgs = "Eigen::SimplicialLLT >") + } + ) +) diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h index e4607f03..9f027151 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h @@ -211,6 +211,8 @@ Eigen::Tensor binaryOp( return binaryOp(xEval, y); } +// TODO: presumably delete this SparseCholesky stuff + // forward declaration of nCompiler struct to store Sparse Chol. decompositions class SparseCholesky; @@ -747,6 +749,7 @@ Eigen::Tensor asDense(SpMatExpr &x) { return res; } +/* template SparseCholType nChol(const Eigen::SparseMatrix &x) { Eigen::SimplicialLLT> llt(x); @@ -758,6 +761,15 @@ SparseCholType nChol(const Eigen::SparseMatrix &x) { res.P = P; return res; } +*/ + +// Should this use SparseMatrix? +std::shared_ptr>> sparseCholFactor(const Eigen::SparseMatrix &x) { + //Eigen::SimplicialLLT> llt(x); + std::shared_ptr>> llt(new Eigen::SimplicialLLT>); + llt->compute(x); + return llt; +} /** * Compute the Cholesky decomposition for a symmetric matrix stored as an From 6c99fbc9d78e6182ef2a1cd0e483d12a531fe262 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Mon, 4 May 2026 08:33:41 -0700 Subject: [PATCH 02/13] Add draft of sparseCholLegdet. --- nCompiler/R/Rexecution.R | 21 ++++ nCompiler/R/compile_aaa_operatorLists.R | 9 ++ nCompiler/R/compile_labelAbstractTypes.R | 9 ++ .../ET_ext/post_Rcpp/tensorOperations.h | 12 ++ .../test-tensorOperations_sparse.R | 105 ++++++++++++++++++ 5 files changed, 156 insertions(+) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index 87b9c5fc..c641285b 100644 --- a/nCompiler/R/Rexecution.R +++ b/nCompiler/R/Rexecution.R @@ -322,6 +322,27 @@ sparseCholFactor <- function(x) { Matrix::Cholesky(x) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky } +sparseCholLogdet <- function(chol) { + return(sum(log(diag(expand1(chol, "L"))))) +} + +sparseCholSolve <- function(ch, x) { + ## E.g. for quadratic form with sparse covariance + return(solve(ch, x)) + # solve(ch, solve(ch, solve(ch, solve(ch, x, system="P"), system = "L"), system = "Lt"), system = "Pt")) +} + +sparseCholBacksolve <- function(ch, x) { + ## E.g., for generating MVN draws with sparse precision. + return(solve(ch, solve(ch, x, system = "Lt"), system = "Pt")) # P^{top} U*^{-1} x +} + +sparseCholMult <- function(ch, x) { + ## E.g., for gnerating MVN draws with sparse covariance. + solve(ch, expand1(ch, "L") %*% x, system = "Pt") # P^{top} L* x (see ?Matrix:::solve) +} + + ## diag(Cholesky(x)) ## result <- solve(ch, b, system = "A"), or "L" diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index 0c231939..faeb1c9b 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1115,6 +1115,15 @@ assignOperatorDef( ) ) +assignOperatorDef( + 'sparseCholLogdet', + list( + labelAbstractTypes = list( + handler = 'sparseCholLogdet' + ) + ) +) + assignOperatorDef( c('nChol'), list( diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index d0cf0a73..0657615a 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -1616,6 +1616,15 @@ inLabelAbstractTypesEnv( } ) +inLabelAbstractTypesEnv( + sparseCholLogdet <- function(code, symTab, auxEnv, handlingInfo) { + if(exists('paciorek')) browser() + insertions <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) + code$type <- symbolBasic$new(nDim = 0, type = 'double') + } +) + + inLabelAbstractTypesEnv( nEigen <- function(code, symTab, auxEnv, handlingInfo) { # determine object's natural type diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h index 9f027151..46ed94e3 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h @@ -766,11 +766,23 @@ SparseCholType nChol(const Eigen::SparseMatrix &x) { // Should this use SparseMatrix? std::shared_ptr>> sparseCholFactor(const Eigen::SparseMatrix &x) { //Eigen::SimplicialLLT> llt(x); + // auto llt = std::make_shared>>() std::shared_ptr>> llt(new Eigen::SimplicialLLT>); llt->compute(x); return llt; } +// Use Scalar as return type? +double sparseCholLogdet(const std::shared_ptr>>& llt) { + double logdet; // = llt->matrixL().diagonal().array().log().sum(); // + Eigen::SparseMatrix L; + L = llt->matrixL(); // Can't chain diagonal() on a triangular view, so instantiate the L matrix. + logdet = L.diagonal().array().log().sum(); + return logdet; +} + + + /** * Compute the Cholesky decomposition for a symmetric matrix stored as an * Eigen::Tensor object, or derived from a Tensor expression object diff --git a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R index e9e017e6..2c567fa2 100644 --- a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R +++ b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R @@ -270,3 +270,108 @@ expect_equal(cAdd2_force_dense_mixed(x = M_sparse, y = M2), M + M2) expect_equal(cAsSparse(x = M), M_sparse) expect_equal(cAsDense(x = M_sparse), M) expect_equal(cAdd2_force_dense_unnecessary(x = M, y = M2), M + M2) + + + +sparseCholFactor <- function(x) { + Matrix::Cholesky(x, LDL = FALSE) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky +} + +sparseCholLogdet <- function(chol) { + return(sum(log(diag(expand1(chol, "L"))))) +} + +sparseCholSolve <- function(ch, x) { + return(solve(ch, x)) + # solve(ch, solve(ch, solve(ch, solve(ch, x, system="P"), system = "L"), system = "Lt"), system = "Pt")) +} + + +sparseCholBacksolve <- function(ch, x) { + return(solve(ch, solve(ch, x, system = "Lt"), system = "Pt")) # P^{top} U*^{-1} x +} + +## Seems ok in dist. +sparseCholMult <- function(ch, x) { + solve(ch, expand1(ch, "L") %*% x, system = "Pt") # P^{top} L* x (see ?Matrix:::solve) +} + + +set.seed(1) +m <- 10 +n <- 10 +# proportion of non-zero entries +p <- .1 + +# set up a matrix with only a few random non-zero entries that will cause a permutation +M <- matrix(data = 0, nrow = m, ncol = n) +M[sample(x = length(M), size = p * length(M))] <- runif(n = p * length(M), 0, .2) +diag(M) <- 1 +M[lower.tri(M)] <- t(M)[lower.tri(M)] +M_sparse <- as(Matrix::Matrix(M, sparse = TRUE), 'generalMatrix') + +x <- rnorm(10) +rch <- chol(M) + +sum(log(diag(rch))) +backsolve(rch, forwardsolve(t(rch), x)) +forwardsolve(t(rch), x) +t(rch) %*% x + +ch <- Cholesky(M_sparse, LDL=FALSE) + + +ch <- sparseCholFactor(M_sparse) +sparseCholLogdet(ch) +sparseCholSolve(ch, x) # Correct +sparseCholForwardsolve(ch, x) + + +y1 <- solve(ch, rhs, system = "P") # y1 = P rhs (apply permutation) +y2 <- solve(ch, y1, system = "L") # y2 = L^{-1} y1 (forward solve) +y2a <- solve(ch, y1, system = "D") # y2 = L^{-1} y1 (forward solve) +y3 <- solve(ch, y2a, system = "Lt") # y3 = L'^{-1} y2 (back solve) +b <- solve(ch, y3, system = "Pt") # b = P' y3 (undo permutation) + + +m <- 10000 +out2 <- out1 <- matrix(0,m, 10) +set.seed(1) +for(i in 1:m) { + x <- rnorm(10) + out1[i,] <- t(rch) %*% x + out2[i,] <- sparseCholMult(ch, x)[,1] +} +plot(M,cov(out1), type = 'p') +points(M,cov(out2),col='red') + +m <- 10000 +out2 <- out1 <- matrix(0,m, 10) +set.seed(1) +for(i in 1:m) { + x <- rnorm(10) + out1[i,] <- backsolve(rch, x) + out2[i,] <- sparseCholBacksolve(ch, x) +} +plot(solve(M),cov(out1), type = 'p') +points(solve(M),cov(out2),col='red') + + +myfun <- nFunction( + fun = function(Q = 'nSparseMatrix') { + ch <- sparseCholFactor(Q) + return(0) + }, returnType = 'numericScalar') +cmyfun <- nCompile(myfun) + +myfun(m) +cmyfun(sm) + + +myfun <- nFunction( + fun = function(Q = 'nSparseMatrix') { + ch <- sparseCholFactor(Q) + out <- sparseCholLogdet(ch) + return(out) + }, returnType = 'numericScalar') +cmyfun <- nCompile(myfun) From 28e214d08591e663c4441694bce8c0c1710b6286 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Mon, 11 May 2026 11:49:09 -0700 Subject: [PATCH 03/13] Fully handling of sparse Cholesky, with sparseCholLogdet as first example operator. --- nCompiler/R/Rexecution.R | 2 +- nCompiler/R/compile_aaa_operatorLists.R | 4 +-- nCompiler/R/compile_labelAbstractTypes.R | 27 ++++++++++++++++--- nCompiler/R/symbolTable.R | 13 +++++---- nCompiler/R/typeDeclarations.R | 3 +++ nCompiler/R/zzz_NC_Predefined.R | 11 ++++++++ .../ET_ext/post_Rcpp/tensorOperations.h | 20 -------------- 7 files changed, 49 insertions(+), 31 deletions(-) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index c641285b..2b95a2f9 100644 --- a/nCompiler/R/Rexecution.R +++ b/nCompiler/R/Rexecution.R @@ -318,7 +318,7 @@ nChol <- function(x) { #' #' @examples TODO #' -sparseCholFactor <- function(x) { +sparseChol <- function(x) { Matrix::Cholesky(x) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky } diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index faeb1c9b..333d9661 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1107,10 +1107,10 @@ assignOperatorDef( ) assignOperatorDef( - 'sparseCholFactor', + 'sparseChol', list( labelAbstractTypes = list( - handler = 'sparseCholFactor' + handler = 'sparseChol' ) ) ) diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 0657615a..7c6f56e8 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -521,6 +521,26 @@ inLabelAbstractTypesEnv( } ) +inLabelAbstractTypesEnv( + sparseChol <- function(code, symTab, auxEnv, handlingInfo) { + inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) + argType <- code$args[[1]]$type + if(exists('paciorek')) browser() + 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) @@ -1602,12 +1622,14 @@ inLabelAbstractTypesEnv( ) inLabelAbstractTypesEnv( - sparseCholFactor <- function(code, symTab, auxEnv, handlingInfo) { + # 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 + } else + ## Put dense chol handling here stop(exprClassProcessingErrorMsg( code, 'sparse Cholesky factor only supported for sparse matrices.' @@ -1618,7 +1640,6 @@ inLabelAbstractTypesEnv( inLabelAbstractTypesEnv( sparseCholLogdet <- function(code, symTab, auxEnv, handlingInfo) { - if(exists('paciorek')) browser() insertions <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) code$type <- symbolBasic$new(nDim = 0, type = 'double') } diff --git a/nCompiler/R/symbolTable.R b/nCompiler/R/symbolTable.R index 0d8b44c7..329c9ac3 100644 --- a/nCompiler/R/symbolTable.R +++ b/nCompiler/R/symbolTable.R @@ -635,16 +635,19 @@ symbolSimplicialLLT <- R6::R6Class( inherit = symbolBase, portable = TRUE, public = list( + interface = FALSE, initialize = function(...) { - super$initialize(...) - self$type = "simplicialLLT" # "simplicialLLT_ptr" ? + super$initialize(..., interface = self$interface) + self$type <- "simplicialLLT" self }, - print = function() print(self$type), + print = function() { + writeLines(paste0(self$name, ': symbolSimplicialLLT')) + }, genCppVar = function() { cppVarFullClass$new(name = self$name, - baseType = "std::shared_ptr", - templateArgs = "Eigen::SimplicialLLT >") + baseType = "Eigen::SimplicialLLT", # "std::shared_ptr", + templateArgs = list("Eigen::SparseMatrix")) } ) ) diff --git a/nCompiler/R/typeDeclarations.R b/nCompiler/R/typeDeclarations.R index 603ac6d4..37ab1657 100644 --- a/nCompiler/R/typeDeclarations.R +++ b/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) diff --git a/nCompiler/R/zzz_NC_Predefined.R b/nCompiler/R/zzz_NC_Predefined.R index 4998a11e..eff5389a 100644 --- a/nCompiler/R/zzz_NC_Predefined.R +++ b/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' + ) +) + diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h index 46ed94e3..e0dfba2e 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h @@ -763,26 +763,6 @@ SparseCholType nChol(const Eigen::SparseMatrix &x) { } */ -// Should this use SparseMatrix? -std::shared_ptr>> sparseCholFactor(const Eigen::SparseMatrix &x) { - //Eigen::SimplicialLLT> llt(x); - // auto llt = std::make_shared>>() - std::shared_ptr>> llt(new Eigen::SimplicialLLT>); - llt->compute(x); - return llt; -} - -// Use Scalar as return type? -double sparseCholLogdet(const std::shared_ptr>>& llt) { - double logdet; // = llt->matrixL().diagonal().array().log().sum(); // - Eigen::SparseMatrix L; - L = llt->matrixL(); // Can't chain diagonal() on a triangular view, so instantiate the L matrix. - logdet = L.diagonal().array().log().sum(); - return logdet; -} - - - /** * Compute the Cholesky decomposition for a symmetric matrix stored as an * Eigen::Tensor object, or derived from a Tensor expression object From 5253e4ad012b900fb6a7e43def52fb88174dc98b Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Tue, 12 May 2026 10:41:38 -0700 Subject: [PATCH 04/13] Overload nLogdet for sparse chol factor and remove `sparseCholLogdet`. --- nCompiler/R/Rexecution.R | 10 +++----- nCompiler/R/compile_aaa_operatorLists.R | 11 +-------- nCompiler/R/compile_labelAbstractTypes.R | 23 +++++++++++++------ .../test-tensorOperations_sparse.R | 5 ---- 4 files changed, 20 insertions(+), 29 deletions(-) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index 2b95a2f9..e132abb9 100644 --- a/nCompiler/R/Rexecution.R +++ b/nCompiler/R/Rexecution.R @@ -322,9 +322,6 @@ sparseChol <- function(x) { Matrix::Cholesky(x) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky } -sparseCholLogdet <- function(chol) { - return(sum(log(diag(expand1(chol, "L"))))) -} sparseCholSolve <- function(ch, x) { ## E.g. for quadratic form with sparse covariance @@ -350,14 +347,13 @@ sparseCholMult <- function(ch, x) { #' #' 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) } diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index 333d9661..3e0bb5ba 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1115,15 +1115,6 @@ assignOperatorDef( ) ) -assignOperatorDef( - 'sparseCholLogdet', - list( - labelAbstractTypes = list( - handler = 'sparseCholLogdet' - ) - ) -) - assignOperatorDef( c('nChol'), list( @@ -1138,7 +1129,7 @@ assignOperatorDef( c('nLogdet'), list( labelAbstractTypes = list( - handler = 'UnaryReduction' + handler = 'nLogdet' ) ) ) diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 7c6f56e8..60b92867 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -525,7 +525,6 @@ inLabelAbstractTypesEnv( sparseChol <- function(code, symTab, auxEnv, handlingInfo) { inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) argType <- code$args[[1]]$type - if(exists('paciorek')) browser() if(inherits(argType, 'symbolSparse')) { # Cholesky factor of a sparse matrix is a collection of information code$type <- symbolNC$new( @@ -1639,12 +1638,22 @@ inLabelAbstractTypesEnv( ) inLabelAbstractTypesEnv( - sparseCholLogdet <- function(code, symTab, auxEnv, handlingInfo) { - insertions <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) - code$type <- symbolBasic$new(nDim = 0, type = 'double') - } -) - + 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) { diff --git a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R index 2c567fa2..38216081 100644 --- a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R +++ b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R @@ -272,15 +272,10 @@ expect_equal(cAsDense(x = M_sparse), M) expect_equal(cAdd2_force_dense_unnecessary(x = M, y = M2), M + M2) - sparseCholFactor <- function(x) { Matrix::Cholesky(x, LDL = FALSE) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky } -sparseCholLogdet <- function(chol) { - return(sum(log(diag(expand1(chol, "L"))))) -} - sparseCholSolve <- function(ch, x) { return(solve(ch, x)) # solve(ch, solve(ch, solve(ch, solve(ch, x, system="P"), system = "L"), system = "Lt"), system = "Pt")) From 60986bf3d586fe507b53c24ed9d7df2f9f64166b Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 13 May 2026 09:51:46 -0700 Subject: [PATCH 05/13] Overload nMul for sparse chol factor. --- nCompiler/R/Rexecution.R | 21 +++++------ nCompiler/R/compile_labelAbstractTypes.R | 3 +- .../ET_ext/post_Rcpp/tensorOperations.h | 35 ++++++++++++++----- .../test-tensorOperations_sparse.R | 4 +-- 4 files changed, 41 insertions(+), 22 deletions(-) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index e132abb9..cb1bc361 100644 --- a/nCompiler/R/Rexecution.R +++ b/nCompiler/R/Rexecution.R @@ -319,14 +319,7 @@ nChol <- function(x) { #' @examples TODO #' sparseChol <- function(x) { - Matrix::Cholesky(x) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky -} - - -sparseCholSolve <- function(ch, x) { - ## E.g. for quadratic form with sparse covariance - return(solve(ch, x)) - # solve(ch, solve(ch, solve(ch, solve(ch, x, system="P"), system = "L"), system = "Lt"), system = "Pt")) + Matrix::Cholesky(x, LDL = FALSE) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky } sparseCholBacksolve <- function(ch, x) { @@ -352,7 +345,7 @@ sparseCholMult <- function(ch, x) { #' @export #' nLogdet <- function(x) { - if(inherits(x, 'dCHMsimpl') + if(inherits(x, 'dCHMsimpl')) return(sum(log(diag(expand1(chol, "L"))))) ldet <- determinant(x, logarithm = TRUE) ifelse(ldet$sign >= 0, ldet$modulus, NaN) @@ -450,14 +443,22 @@ 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 diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 60b92867..e15491bb 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -1496,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) diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h index e0dfba2e..6b56963a 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h @@ -211,21 +211,37 @@ Eigen::Tensor binaryOp( return binaryOp(xEval, y); } -// TODO: presumably delete this SparseCholesky stuff +// forward declaration of class to store sparse Cholesky decompositions +// class SparseCholesky; +class sparseCholFactor; -// forward declaration of nCompiler struct to store Sparse Chol. decompositions -class SparseCholesky; +// 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 -struct IsSparseCholesky : std::is_base_of< - SparseCholesky, +struct IsSparseCholFactor : std::is_base_of< + std::shared_ptr, 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 by SFINAE + * @tparam N Number of dimensions to test for + */ +template +constexpr typename std::enable_if::value, bool>::type + HasNumDimensionsN() { + return N == 2; // matrices are inherently 2-dimensional + } /** * Template meta programming check to see if Class is an Eigen::SparseMatrix @@ -396,7 +412,7 @@ std::conditional< template struct IsEvaluatedType : std::conditional< IsSparseType::value || IsTensor::value || - IsSparseCholesky::value || std::is_arithmetic::value || + IsSparseCholFactor::value || std::is_arithmetic::value || IsMap::value || IsTranspose::value, std::true_type, std::false_type @@ -512,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 * @@ -969,7 +986,7 @@ Eigen::SparseMatrix 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. * @@ -1334,6 +1351,7 @@ Eigen::Tensor nBacksolve( return triangularsolve(U, b); } + /** * Matrix multiplication x %*% y when both inputs are matrix-like objects, i.e., * rank 2 Eigen::Tensor objects, or Tensor expressions @@ -1597,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 diff --git a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R index 38216081..eb77effe 100644 --- a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R +++ b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R @@ -365,8 +365,8 @@ cmyfun(sm) myfun <- nFunction( fun = function(Q = 'nSparseMatrix') { - ch <- sparseCholFactor(Q) - out <- sparseCholLogdet(ch) + ch <- sparseChol(Q) + out <- nLogdet(ch) return(out) }, returnType = 'numericScalar') cmyfun <- nCompile(myfun) From ac509f44070a6e25e677d677b25f5768e1ae0ef4 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 14 May 2026 08:30:02 -0700 Subject: [PATCH 06/13] Actually add tensorOperations_chol.h with C++ implementations for sparse chol operations. --- .../ET_ext/post_Rcpp/tensorOperations_chol.h | 113 ++++++++++++++++++ 1 file changed, 113 insertions(+) create mode 100644 nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h new file mode 100644 index 00000000..7f9c03c1 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h @@ -0,0 +1,113 @@ +#ifndef TENSOROPERATIONS_CHOL_H +#define TENSOROPERATIONS_CHOL_H + + +std::shared_ptr sparseChol(const Eigen::SparseMatrix& x) { + std::shared_ptr ch = nClass_builder()(); + ch->llt.compute(x); + if (ch->llt.info() != Eigen::Success) + throw std::runtime_error("Cholesky factorization failed"); + return ch; +} + +// This overloads nLogdet, on top of use for standard matrices with calculation +// done via the SVD. +double nLogdet(std::shared_ptr ch) { + Eigen::SparseMatrix L = ch->llt.matrixL(); + return L.diagonal().array().log().sum(); +} + +// Should we implement nLogdet for operation directly on sparse matrix? + +template +Eigen::Tensor nSolve(std::shared_ptr ch, const RHS & b + ) { + // explicit Eigen::Tensor types for inputs + typedef Eigen::Tensor bTensor; + // evaluate arguments, if necessary + const auto & b_eval = eval(b); + // initialize storage for solution, given problem dimensions + auto bdim = b.dimensions(); + bTensor res = Eigen::Tensor(bdim[0], bdim.size() > 1 ? bdim[1]: 1); + // map tensor objects to Eigen::Matrix types + auto bmap = matmap(b_eval); + auto resMap = matmap(res); + // solve linear system + resMap = ch->llt.solve(bmap); + return res; +} + +// In progress +template +Eigen::Tensor nBacksolve(std::shared_ptr ch, const RHS & b + ) { + // explicit Eigen::Tensor types for inputs + typedef Eigen::Tensor bTensor; + // evaluate arguments, if necessary + const auto & b_eval = eval(b); + // initialize storage for solution, given problem dimensions + auto bdim = b.dimensions(); + bTensor res = Eigen::Tensor(bdim[0], bdim.size() > 1 ? bdim[1]: 1); + // map tensor objects to Eigen::Matrix types + auto bmap = matmap(b_eval); + auto resMap = matmap(res); + // solve linear system + resMap = ch->llt.permutationPinv() * ch->llt.matrixU().triangularView().solve(bmap); + return res; +} + +// Multiply L by matrix. +template< + typename Ypr, + typename std::enable_if< + HasNumDimensionsN(), + Ypr + >::type* = nullptr, + typename ResultType = typename std::conditional< + IsSparseType::value, + Eigen::SparseMatrix, + Eigen::Tensor + >::type +> +ResultType nMul(std::shared_ptr ch, const Ypr & y) { + // evaluate arguments, if necessary + const auto & yeval = eval(y); + // map inputs and initialize output + auto ymap = matmap(yeval); + ResultType res(ymap.rows(), ymap.cols()); + // map and multiply! + matmap(res) = ch->llt.permutationPinv() * (ch->llt.matrixL() * ymap); + return res; +} + +// Multiply L by vector. +// TODO: is there any way for us to return a vector instead of a matrix? +template< + typename Ypr, + typename std::enable_if< + HasNumDimensionsN(), + Ypr + >::type* = nullptr, + typename ResultType = typename std::conditional< + IsSparseType::value, + Eigen::SparseMatrix, + Eigen::Tensor + >::type +> +ResultType nMul(std::shared_ptr ch, const Ypr & y) { + // evaluate arguments, if necessary + const auto & yeval = eval(y); + // map inputs and initialize output + auto ymap = matmap(yeval); + // initialize output + ResultType res(ymap.size(), 1); + // map and multiply! + matmap(res) = ch->llt.permutationPinv() * (ch->llt.matrixL() * ymap); + return res; +} + + + +// Add function for dense chol factor too. + +#endif From fbbd93a10afdef2d88757ebf864dabe6eb857e9c Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 14 May 2026 08:33:00 -0700 Subject: [PATCH 07/13] Omit stray test/check code in test-tensorOperations_sparse.R. --- .../test-tensorOperations_sparse.R | 99 ------------------- 1 file changed, 99 deletions(-) diff --git a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R index eb77effe..43233ea7 100644 --- a/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R +++ b/nCompiler/tests/testthat/tensorOps_tests/test-tensorOperations_sparse.R @@ -271,102 +271,3 @@ expect_equal(cAsSparse(x = M), M_sparse) expect_equal(cAsDense(x = M_sparse), M) expect_equal(cAdd2_force_dense_unnecessary(x = M, y = M2), M + M2) - -sparseCholFactor <- function(x) { - Matrix::Cholesky(x, LDL = FALSE) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky -} - -sparseCholSolve <- function(ch, x) { - return(solve(ch, x)) - # solve(ch, solve(ch, solve(ch, solve(ch, x, system="P"), system = "L"), system = "Lt"), system = "Pt")) -} - - -sparseCholBacksolve <- function(ch, x) { - return(solve(ch, solve(ch, x, system = "Lt"), system = "Pt")) # P^{top} U*^{-1} x -} - -## Seems ok in dist. -sparseCholMult <- function(ch, x) { - solve(ch, expand1(ch, "L") %*% x, system = "Pt") # P^{top} L* x (see ?Matrix:::solve) -} - - -set.seed(1) -m <- 10 -n <- 10 -# proportion of non-zero entries -p <- .1 - -# set up a matrix with only a few random non-zero entries that will cause a permutation -M <- matrix(data = 0, nrow = m, ncol = n) -M[sample(x = length(M), size = p * length(M))] <- runif(n = p * length(M), 0, .2) -diag(M) <- 1 -M[lower.tri(M)] <- t(M)[lower.tri(M)] -M_sparse <- as(Matrix::Matrix(M, sparse = TRUE), 'generalMatrix') - -x <- rnorm(10) -rch <- chol(M) - -sum(log(diag(rch))) -backsolve(rch, forwardsolve(t(rch), x)) -forwardsolve(t(rch), x) -t(rch) %*% x - -ch <- Cholesky(M_sparse, LDL=FALSE) - - -ch <- sparseCholFactor(M_sparse) -sparseCholLogdet(ch) -sparseCholSolve(ch, x) # Correct -sparseCholForwardsolve(ch, x) - - -y1 <- solve(ch, rhs, system = "P") # y1 = P rhs (apply permutation) -y2 <- solve(ch, y1, system = "L") # y2 = L^{-1} y1 (forward solve) -y2a <- solve(ch, y1, system = "D") # y2 = L^{-1} y1 (forward solve) -y3 <- solve(ch, y2a, system = "Lt") # y3 = L'^{-1} y2 (back solve) -b <- solve(ch, y3, system = "Pt") # b = P' y3 (undo permutation) - - -m <- 10000 -out2 <- out1 <- matrix(0,m, 10) -set.seed(1) -for(i in 1:m) { - x <- rnorm(10) - out1[i,] <- t(rch) %*% x - out2[i,] <- sparseCholMult(ch, x)[,1] -} -plot(M,cov(out1), type = 'p') -points(M,cov(out2),col='red') - -m <- 10000 -out2 <- out1 <- matrix(0,m, 10) -set.seed(1) -for(i in 1:m) { - x <- rnorm(10) - out1[i,] <- backsolve(rch, x) - out2[i,] <- sparseCholBacksolve(ch, x) -} -plot(solve(M),cov(out1), type = 'p') -points(solve(M),cov(out2),col='red') - - -myfun <- nFunction( - fun = function(Q = 'nSparseMatrix') { - ch <- sparseCholFactor(Q) - return(0) - }, returnType = 'numericScalar') -cmyfun <- nCompile(myfun) - -myfun(m) -cmyfun(sm) - - -myfun <- nFunction( - fun = function(Q = 'nSparseMatrix') { - ch <- sparseChol(Q) - out <- nLogdet(ch) - return(out) - }, returnType = 'numericScalar') -cmyfun <- nCompile(myfun) From efca3f0f98dccdc12b86e3d39035478eafcf63fa Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 14 May 2026 08:41:22 -0700 Subject: [PATCH 08/13] Clean up R execution for sparse chol factor input to various lin alg operators. --- nCompiler/R/Rexecution.R | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/nCompiler/R/Rexecution.R b/nCompiler/R/Rexecution.R index cb1bc361..e13dc874 100644 --- a/nCompiler/R/Rexecution.R +++ b/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,6 +305,8 @@ nDiag <- function(x, ...) { #' @export #' nChol <- function(x) { + if(inherits(r, 'dCHMsimpl')) + stop("`nChol` not meaningful for sparse Cholesky factor input") chol(x) } @@ -318,24 +326,10 @@ nChol <- function(x) { #' #' @examples TODO #' -sparseChol <- function(x) { - Matrix::Cholesky(x, LDL = FALSE) # unlike Matrix::chol, this does permutation and returns representation of the Cholesky +sparseCholFactor <- function(x) { + Matrix::Cholesky(x, LDL = FALSE) } -sparseCholBacksolve <- function(ch, x) { - ## E.g., for generating MVN draws with sparse precision. - return(solve(ch, solve(ch, x, system = "Lt"), system = "Pt")) # P^{top} U*^{-1} x -} - -sparseCholMult <- function(ch, x) { - ## E.g., for gnerating MVN draws with sparse covariance. - solve(ch, expand1(ch, "L") %*% x, system = "Pt") # P^{top} L* x (see ?Matrix:::solve) -} - - -## diag(Cholesky(x)) -## result <- solve(ch, b, system = "A"), or "L" - #' Compute the log-determinant of a matrix #' #' In a \code{nFunction}, \code{nLogdet} is identical to \code{logdet} @@ -465,6 +459,9 @@ nSolve <- function(a, b) { #' #' @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) } @@ -472,6 +469,8 @@ nForwardsolve <- function(l, x) { #' #' @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) } From 47c77473a522f2063f8f733a8bba04bcd74a1964 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 14 May 2026 09:12:40 -0700 Subject: [PATCH 09/13] Add predefined C++ for sparseCholFactor class. --- .../sparseCholFactor_copyFiles.txt | 0 .../sparseCholFactor_cppContent.cpp | 43 +++++++++++++++++++ .../sparseCholFactor_filebase.txt | 1 + .../sparseCholFactor_hContent.h | 24 +++++++++++ .../sparseCholFactor_manifest.txt | 7 +++ .../sparseCholFactor_post_cpp_compiler.txt | 0 .../sparseCholFactor_preamble.cpp | 6 +++ 7 files changed, 81 insertions(+) create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_copyFiles.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_cppContent.cpp create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_filebase.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_manifest.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_post_cpp_compiler.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_preamble.cpp diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_copyFiles.txt b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_copyFiles.txt new file mode 100644 index 00000000..e69de29b diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_cppContent.cpp b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_cppContent.cpp new file mode 100644 index 00000000..1c18aa87 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_cppContent.cpp @@ -0,0 +1,43 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __sparseCholFactor_CPP +#define __sparseCholFactor_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "sparseCholFactor_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] + + sparseCholFactor::sparseCholFactor ( ) { +RESET_EIGEN_ERRORS +} + +// [[Rcpp::export(name = "sparseCholFactor_new")]] + SEXP new_sparseCholFactor ( ) { +RESET_EIGEN_ERRORS +return CREATE_NEW_NCOMP_OBJECT(sparseCholFactor);; +} + +// [[Rcpp::export(name = "set_CnClass_env_sparseCholFactor_new")]] + void set_CnClass_env_sparseCholFactor ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(sparseCholFactor, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_sparseCholFactor_new")]] + Rcpp::Environment get_CnClass_env_sparseCholFactor ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(sparseCholFactor);; +} + +NCOMPILER_INTERFACE( +sparseCholFactor, +NCOMPILER_FIELDS(), +NCOMPILER_METHODS() +) +#endif diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_filebase.txt b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_filebase.txt new file mode 100644 index 00000000..0f6dfd79 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_filebase.txt @@ -0,0 +1 @@ +sparseCholFactor_c_ diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h new file mode 100644 index 00000000..7c455795 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h @@ -0,0 +1,24 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __sparseCholFactor_H +#define __sparseCholFactor_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include + +class sparseCholFactor : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + sparseCholFactor ( ) ; + Eigen::SimplicialLLT> llt; +}; + + SEXP new_sparseCholFactor ( ) ; + + void set_CnClass_env_sparseCholFactor ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_sparseCholFactor ( ) ; + +#include + +#endif diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_manifest.txt b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_manifest.txt new file mode 100644 index 00000000..ffb07a3d --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1778520082.54572, class = c("POSIXct", +"POSIXt")), packet_name = "sparseCholFactor", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles" +), files = list(preamble = "sparseCholFactor_preamble.cpp", cppContent = "sparseCholFactor_cppContent.cpp", + hContent = "sparseCholFactor_hContent.h", filebase = "sparseCholFactor_filebase.txt", + post_cpp_compiler = "sparseCholFactor_post_cpp_compiler.txt", + copyFiles = "sparseCholFactor_copyFiles.txt")) diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_post_cpp_compiler.txt b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_post_cpp_compiler.txt new file mode 100644 index 00000000..e69de29b diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_preamble.cpp b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_preamble.cpp new file mode 100644 index 00000000..61410333 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_preamble.cpp @@ -0,0 +1,6 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +// #define NCOMPILER_USES_TBB +#define NCOMPILER_USES_NLIST +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE From 29d2a8c79731a018e9055b912fdf27578a4e892b Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sat, 30 May 2026 13:35:42 -0700 Subject: [PATCH 10/13] Add dense chol factorization and related methods, akin to sparseChol. --- nCompiler/R/compile_aaa_operatorLists.R | 9 ++++ nCompiler/R/compile_labelAbstractTypes.R | 48 ++++++++++++------- nCompiler/R/symbolTable.R | 24 +++++++++- nCompiler/R/typeDeclarations.R | 3 ++ nCompiler/R/zzz_NC_Predefined.R | 11 +++++ .../ET_ext/post_Rcpp/tensorOperations.h | 28 +++++++++++ .../ET_ext/post_Rcpp/tensorOperations_chol.h | 10 ++-- 7 files changed, 111 insertions(+), 22 deletions(-) diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index 3e0bb5ba..31978769 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1115,6 +1115,15 @@ assignOperatorDef( ) ) +assignOperatorDef( + 'denseChol', + list( + labelAbstractTypes = list( + handler = 'denseChol' + ) + ) +) + assignOperatorDef( c('nChol'), list( diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index e15491bb..263facda 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -521,39 +521,57 @@ inLabelAbstractTypesEnv( } ) + +## TODO: merge these next two to become `Cholesky`. inLabelAbstractTypesEnv( 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) + } else { + code$type <- symbolNC$new( + name = code$name, type = 'denseCholFactor', NCgenerator = sparseCholFactor, 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 = denseCholFactor, 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 { + if(inherits(argType, 'symbolSparse')) { + + + } 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) } ) @@ -1629,11 +1647,7 @@ inLabelAbstractTypesEnv( 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) + code$type <- symbolLLT$new(name = code$name) invisible(NULL) } ) diff --git a/nCompiler/R/symbolTable.R b/nCompiler/R/symbolTable.R index 329c9ac3..9c11e635 100644 --- a/nCompiler/R/symbolTable.R +++ b/nCompiler/R/symbolTable.R @@ -646,8 +646,30 @@ symbolSimplicialLLT <- R6::R6Class( }, genCppVar = function() { cppVarFullClass$new(name = self$name, - baseType = "Eigen::SimplicialLLT", # "std::shared_ptr", + baseType = "Eigen::SimplicialLLT", templateArgs = list("Eigen::SparseMatrix")) } ) ) + +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? + } + ) +) diff --git a/nCompiler/R/typeDeclarations.R b/nCompiler/R/typeDeclarations.R index 37ab1657..d8e421a6 100644 --- a/nCompiler/R/typeDeclarations.R +++ b/nCompiler/R/typeDeclarations.R @@ -271,6 +271,9 @@ typeDeclarationList <- list( simplicialLLT = function() { symbolSimplicialLLT$new() }, + LLT = function() { + symbolLLT$new() + }, nList = function(type) { elementSym <- argType2symbol(type) symbolNlist$new(elementSym = elementSym) diff --git a/nCompiler/R/zzz_NC_Predefined.R b/nCompiler/R/zzz_NC_Predefined.R index eff5389a..8a7e3761 100644 --- a/nCompiler/R/zzz_NC_Predefined.R +++ b/nCompiler/R/zzz_NC_Predefined.R @@ -136,3 +136,14 @@ sparseCholFactor <- nClass( ) ) +#' @export +denseCholFactor <- nClass( + classname = 'denseCholFactor', + predefined = quote(system.file(file.path("include","nCompiler", "predef"), package="nCompiler") |> + file.path("denseCholFactor_nC")), + Cpublic = list( + llt = 'LLT' + ) +) + + diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h index 6b56963a..68f1e78a 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations.h @@ -215,6 +215,8 @@ Eigen::Tensor binaryOp( // class SparseCholesky; class sparseCholFactor; +class denseCholFactor; + // TODO: Not sure if IsSparseCholFactor stuff should be here or in tensorOperations_chol.h. // I think we need to use `IsSparseCholFactor` in `IsEvaluatedType`. @@ -291,6 +293,32 @@ struct IsSparseType : std::conditional< std::false_type >::type { }; +/** + * Template meta programming check to see if Class is an Eigen::DenseCholFactor + * + * @tparam Class type to inspect + */ +template +struct IsDenseCholFactor : std::is_base_of< + std::shared_ptr, + 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 by SFINAE + * @tparam N Number of dimensions to test for + */ +template +constexpr typename std::enable_if::value, bool>::type + HasNumDimensionsN() { + return N == 2; // matrices are inherently 2-dimensional + } + + /** * Template meta programming check to see if Class is an Eigen::Tensor type. * diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h index 7f9c03c1..c0d9aadd 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h @@ -6,10 +6,11 @@ std::shared_ptr sparseChol(const Eigen::SparseMatrix& std::shared_ptr ch = nClass_builder()(); ch->llt.compute(x); if (ch->llt.info() != Eigen::Success) - throw std::runtime_error("Cholesky factorization failed"); + throw std::runtime_error("sparse Cholesky factorization failed"); return ch; } + // This overloads nLogdet, on top of use for standard matrices with calculation // done via the SVD. double nLogdet(std::shared_ptr ch) { @@ -17,6 +18,7 @@ double nLogdet(std::shared_ptr ch) { return L.diagonal().array().log().sum(); } + // Should we implement nLogdet for operation directly on sparse matrix? template @@ -37,7 +39,7 @@ Eigen::Tensor nSolve(std::shared_ptr Eigen::Tensor nBacksolve(std::shared_ptr ch, const RHS & b ) { @@ -56,6 +58,7 @@ Eigen::Tensor nBacksolve(std::shared_p return res; } + // Multiply L by matrix. template< typename Ypr, @@ -80,6 +83,7 @@ ResultType nMul(std::shared_ptr ch, const Ypr & y) { return res; } + // Multiply L by vector. // TODO: is there any way for us to return a vector instead of a matrix? template< @@ -107,7 +111,5 @@ ResultType nMul(std::shared_ptr ch, const Ypr & y) { } - -// Add function for dense chol factor too. #endif From 6c4cbfea51bcb0f0cf2653201fca709f5ba7c53c Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sat, 30 May 2026 13:36:30 -0700 Subject: [PATCH 11/13] Add cpp files for denseChol. --- .../post_Rcpp/tensorOperations_denseChol.h | 113 ++++++++++++++++++ .../denseCholFactor_copyFiles.txt | 0 .../denseCholFactor_cppContent.cpp | 43 +++++++ .../denseCholFactor_filebase.txt | 1 + .../denseCholFactor_hContent.h | 24 ++++ .../denseCholFactor_manifest.txt | 7 ++ .../denseCholFactor_post_cpp_compiler.txt | 0 .../denseCholFactor_preamble.cpp | 6 + 8 files changed, 194 insertions(+) create mode 100644 nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_copyFiles.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_cppContent.cpp create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_filebase.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_hContent.h create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_manifest.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_post_cpp_compiler.txt create mode 100644 nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_preamble.cpp diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h new file mode 100644 index 00000000..ef315c62 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h @@ -0,0 +1,113 @@ +#ifndef TENSOROPERATIONS_DENSECHOL_H +#define TENSOROPERATIONS_DENSECHOL_H + + +std::shared_ptr denseChol(const Eigen::Tensor &x) { + auto x_map = matmap(x); + std::shared_ptr ch = nClass_builder()(); + ch->llt.compute(x_map); // ch->llt(x) ? + if (ch->llt.info() != Eigen::Success) + throw std::runtime_error("dense Cholesky factorization failed"); + return ch; +} + +// This overloads nLogdet, on top of use for standard matrices with calculation +// done via the SVD. + +double nLogdet(std::shared_ptr ch) { + return ch->llt.matrixL().nestedExpression().diagonal().array().log().sum(); +} + + + +template +Eigen::Tensor nSolve(std::shared_ptr ch, const RHS & b + ) { + // explicit Eigen::Tensor types for inputs + typedef Eigen::Tensor bTensor; + // evaluate arguments, if necessary + const auto & b_eval = eval(b); + // initialize storage for solution, given problem dimensions + auto bdim = b.dimensions(); + bTensor res = Eigen::Tensor(bdim[0], bdim.size() > 1 ? bdim[1]: 1); + // map tensor objects to Eigen::Matrix types + auto bmap = matmap(b_eval); + auto resMap = matmap(res); + // solve linear system + resMap = ch->llt.solve(bmap); + return res; +} + + +template +Eigen::Tensor nBacksolve(std::shared_ptr ch, const RHS & b + ) { + // explicit Eigen::Tensor types for inputs + typedef Eigen::Tensor bTensor; + // evaluate arguments, if necessary + const auto & b_eval = eval(b); + // initialize storage for solution, given problem dimensions + auto bdim = b.dimensions(); + bTensor res = Eigen::Tensor(bdim[0], bdim.size() > 1 ? bdim[1]: 1); + // map tensor objects to Eigen::Matrix types + auto bmap = matmap(b_eval); + auto resMap = matmap(res); + // solve linear system + resMap = ch->llt.matrixU().solve(bmap); + return res; +} + +// Multiply L by matrix. + +template< + typename Ypr, + typename std::enable_if< + HasNumDimensionsN(), + Ypr + >::type* = nullptr, + typename ResultType = typename std::conditional< + IsSparseType::value, + Eigen::SparseMatrix, + Eigen::Tensor + >::type +> +ResultType nMul(std::shared_ptr ch, const Ypr & y) { + // evaluate arguments, if necessary + const auto & yeval = eval(y); + // map inputs and initialize output + auto ymap = matmap(yeval); + ResultType res(ymap.rows(), ymap.cols()); + // map and multiply! + matmap(res) = ch->llt.matrixL() * ymap; + return res; +} + +// Multiply L by vector. +// TODO: is there any way for us to return a vector instead of a matrix? + + +template< + typename Ypr, + typename std::enable_if< + HasNumDimensionsN(), + Ypr + >::type* = nullptr, + typename ResultType = typename std::conditional< + IsSparseType::value, + Eigen::SparseMatrix, + Eigen::Tensor + >::type +> +ResultType nMul(std::shared_ptr ch, const Ypr & y) { + // evaluate arguments, if necessary + const auto & yeval = eval(y); + // map inputs and initialize output + auto ymap = matmap(yeval); + // initialize output + ResultType res(ymap.size(), 1); + // map and multiply! + matmap(res) = ch->llt.matrixL() * ymap; + return res; +} + +#endif diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_copyFiles.txt b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_copyFiles.txt new file mode 100644 index 00000000..e69de29b diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_cppContent.cpp b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_cppContent.cpp new file mode 100644 index 00000000..ce134ed0 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_cppContent.cpp @@ -0,0 +1,43 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __denseCholFactor_CPP +#define __denseCholFactor_CPP +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include +#include "denseCholFactor_c_.h" +using namespace Rcpp; +// [[Rcpp::plugins(nCompiler_Eigen_plugin)]] +// [[Rcpp::depends(RcppParallel)]] +// [[Rcpp::depends(nCompiler)]] +// [[Rcpp::depends(Rcereal)]] + + denseCholFactor::denseCholFactor ( ) { +RESET_EIGEN_ERRORS +} + +// [[Rcpp::export(name = "denseCholFactor_new")]] + SEXP new_denseCholFactor ( ) { +RESET_EIGEN_ERRORS +return CREATE_NEW_NCOMP_OBJECT(denseCholFactor);; +} + +// [[Rcpp::export(name = "set_CnClass_env_denseCholFactor_new")]] + void set_CnClass_env_denseCholFactor ( SEXP env ) { +RESET_EIGEN_ERRORS +SET_CNCLASS_ENV(denseCholFactor, env);; +} + +// [[Rcpp::export(name = "get_CnClass_env_denseCholFactor_new")]] + Rcpp::Environment get_CnClass_env_denseCholFactor ( ) { +RESET_EIGEN_ERRORS +return GET_CNCLASS_ENV(denseCholFactor);; +} + +NCOMPILER_INTERFACE( +denseCholFactor, +NCOMPILER_FIELDS(), +NCOMPILER_METHODS() +) +#endif diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_filebase.txt b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_filebase.txt new file mode 100644 index 00000000..d7a9ce54 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_filebase.txt @@ -0,0 +1 @@ +denseCholFactor_c_ diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_hContent.h b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_hContent.h new file mode 100644 index 00000000..4b42f518 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_hContent.h @@ -0,0 +1,24 @@ +/* OPENER (Do not edit this comment) */ +#ifndef __denseCholFactor_H +#define __denseCholFactor_H +/* BODY (Do not edit this comment) */ +#ifndef R_NO_REMAP +#define R_NO_REMAP +#endif +#include + +class denseCholFactor : public interface_resolver< genericInterfaceC >, public loadedObjectHookC { +public: + denseCholFactor ( ) ; + Eigen::LLT llt; +}; + + SEXP new_denseCholFactor ( ) ; + + void set_CnClass_env_denseCholFactor ( SEXP env ) ; + + Rcpp::Environment get_CnClass_env_denseCholFactor ( ) ; + +#include + +#endif diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_manifest.txt b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_manifest.txt new file mode 100644 index 00000000..0a117f3e --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_manifest.txt @@ -0,0 +1,7 @@ +list(saved_at = structure(1780167848.49078, class = c("POSIXct", +"POSIXt")), packet_name = "denseCholFactor", elements = c("preamble", +"cppContent", "hContent", "filebase", "post_cpp_compiler", "copyFiles" +), files = list(preamble = "denseCholFactor_preamble.cpp", cppContent = "denseCholFactor_cppContent.cpp", + hContent = "denseCholFactor_hContent.h", filebase = "denseCholFactor_filebase.txt", + post_cpp_compiler = "denseCholFactor_post_cpp_compiler.txt", + copyFiles = "denseCholFactor_copyFiles.txt")) diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_post_cpp_compiler.txt b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_post_cpp_compiler.txt new file mode 100644 index 00000000..e69de29b diff --git a/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_preamble.cpp b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_preamble.cpp new file mode 100644 index 00000000..61410333 --- /dev/null +++ b/nCompiler/inst/include/nCompiler/predef/denseCholFactor_nC/denseCholFactor_preamble.cpp @@ -0,0 +1,6 @@ +#define NCOMPILER_HANDLE_EIGEN_ERRORS +#define NCOMPILER_USES_EIGEN +// #define NCOMPILER_USES_TBB +#define NCOMPILER_USES_NLIST +#define USES_NCOMPILER +#define NCOMPILER_USES_NCLASS_INTERFACE From bcdf8aaefb80564e4aee3431308337824fe69e7b Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sat, 30 May 2026 13:38:48 -0700 Subject: [PATCH 12/13] Add comment about include insertion in predef code for chol factors. --- nCompiler/R/zzz_NC_Predefined.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/nCompiler/R/zzz_NC_Predefined.R b/nCompiler/R/zzz_NC_Predefined.R index 8a7e3761..bbb1c700 100644 --- a/nCompiler/R/zzz_NC_Predefined.R +++ b/nCompiler/R/zzz_NC_Predefined.R @@ -128,6 +128,9 @@ OptimResultList <- nClass( #' @export sparseCholFactor <- nClass( + # manually insert: + #include + # 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")), @@ -138,6 +141,9 @@ sparseCholFactor <- nClass( #' @export denseCholFactor <- nClass( + # manually insert: + #include + # 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")), From b83c4a3ba1d43a16f458a110d5bedf726a70bab0 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sat, 30 May 2026 14:18:59 -0700 Subject: [PATCH 13/13] Add Cholesky with overloading to replace {sparse,dense}Chol. --- nCompiler/R/compile_aaa_operatorLists.R | 10 +++++++++ nCompiler/R/compile_labelAbstractTypes.R | 22 ++++++++++++++++--- nCompiler/R/zzz_NC_Predefined.R | 2 +- .../post_Rcpp/tensorOperations_denseChol.h | 9 ++++++++ ...s_chol.h => tensorOperations_sparseChol.h} | 12 ++++++++-- .../sparseCholFactor_hContent.h | 2 +- 6 files changed, 50 insertions(+), 7 deletions(-) rename nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/{tensorOperations_chol.h => tensorOperations_sparseChol.h} (90%) diff --git a/nCompiler/R/compile_aaa_operatorLists.R b/nCompiler/R/compile_aaa_operatorLists.R index 31978769..bc8f57a5 100644 --- a/nCompiler/R/compile_aaa_operatorLists.R +++ b/nCompiler/R/compile_aaa_operatorLists.R @@ -1106,6 +1106,16 @@ assignOperatorDef( ) ) +assignOperatorDef( + 'Cholesky', + list( + labelAbstractTypes = list( + handler = 'Cholesky' + ) + ) +) + +## TODO: remove these next two as redundant with Cholesky. assignOperatorDef( 'sparseChol', list( diff --git a/nCompiler/R/compile_labelAbstractTypes.R b/nCompiler/R/compile_labelAbstractTypes.R index 263facda..851073de 100644 --- a/nCompiler/R/compile_labelAbstractTypes.R +++ b/nCompiler/R/compile_labelAbstractTypes.R @@ -521,8 +521,24 @@ 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: merge these next two to become `Cholesky`. +## TODO: remove these next two as redundant with Cholesky inLabelAbstractTypesEnv( sparseChol <- function(code, symTab, auxEnv, handlingInfo) { inserts <- recurse_labelAbstractTypes(code, symTab, auxEnv, handlingInfo) @@ -533,7 +549,7 @@ inLabelAbstractTypesEnv( ) } else { code$type <- symbolNC$new( - name = code$name, type = 'denseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE + name = code$name, type = 'denseCholFactor', NCgenerator = denseCholFactor, isArg = FALSE ) } invisible(inserts) @@ -547,7 +563,7 @@ inLabelAbstractTypesEnv( argType <- code$args[[1]]$type if(inherits(argType, 'symbolSparse')) { code$type <- symbolNC$new( - name = code$name, type = 'sparseCholFactor', NCgenerator = denseCholFactor, isArg = FALSE + name = code$name, type = 'sparseCholFactor', NCgenerator = sparseCholFactor, isArg = FALSE ) } else { code$type <- symbolNC$new( diff --git a/nCompiler/R/zzz_NC_Predefined.R b/nCompiler/R/zzz_NC_Predefined.R index bbb1c700..a1395a8b 100644 --- a/nCompiler/R/zzz_NC_Predefined.R +++ b/nCompiler/R/zzz_NC_Predefined.R @@ -129,7 +129,7 @@ OptimResultList <- nClass( #' @export sparseCholFactor <- nClass( # manually insert: - #include + #include # in the hContent file AFTER the sparseCholFactor class declaration classname = 'sparseCholFactor', predefined = quote(system.file(file.path("include","nCompiler", "predef"), package="nCompiler") |> diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h index ef315c62..a498f31c 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_denseChol.h @@ -1,7 +1,16 @@ #ifndef TENSOROPERATIONS_DENSECHOL_H #define TENSOROPERATIONS_DENSECHOL_H +std::shared_ptr Cholesky(const Eigen::Tensor &x) { + auto x_map = matmap(x); + std::shared_ptr ch = nClass_builder()(); + ch->llt.compute(x_map); // ch->llt(x) ? + if (ch->llt.info() != Eigen::Success) + throw std::runtime_error("dense Cholesky factorization failed"); + return ch; +} +// TODO: remove as redundant. std::shared_ptr denseChol(const Eigen::Tensor &x) { auto x_map = matmap(x); std::shared_ptr ch = nClass_builder()(); diff --git a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_sparseChol.h similarity index 90% rename from nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h rename to nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_sparseChol.h index c0d9aadd..57d58a57 100644 --- a/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_chol.h +++ b/nCompiler/inst/include/nCompiler/ET_ext/post_Rcpp/tensorOperations_sparseChol.h @@ -1,7 +1,15 @@ -#ifndef TENSOROPERATIONS_CHOL_H -#define TENSOROPERATIONS_CHOL_H +#ifndef TENSOROPERATIONS_SPARSECHOL_H +#define TENSOROPERATIONS_SPARSECHOL_H +std::shared_ptr Cholesky(const Eigen::SparseMatrix& x) { + std::shared_ptr ch = nClass_builder()(); + ch->llt.compute(x); + if (ch->llt.info() != Eigen::Success) + throw std::runtime_error("sparse Cholesky factorization failed"); + return ch; +} +// TODO: remove as redundant std::shared_ptr sparseChol(const Eigen::SparseMatrix& x) { std::shared_ptr ch = nClass_builder()(); ch->llt.compute(x); diff --git a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h index 7c455795..05d65f3e 100644 --- a/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h +++ b/nCompiler/inst/include/nCompiler/predef/sparseCholFactor_nC/sparseCholFactor_hContent.h @@ -19,6 +19,6 @@ class sparseCholFactor : public interface_resolver< genericInterfaceC +#include #endif