Back to NIMBLE Vignettes
Here we’ll show how to call the Matrix::tcrossprod function from the Matrix package inside a compiled nimbleFunction.
A function to construct a full (dense) matrix, from the component i (row indices), j (column indices), and x (non-zero values) sparse representation.
makeMatrix <- function(i, j, x, nr, nc, ijx) {
## ijx is an (optional, alternate) representation of {i,j,x} as a 3-column array
if(!missing(ijx)) { i <- ijx[,1]; j <- ijx[,2]; x <- ijx[,3] }
if(length(i) != length(j) | length(i) != length(x)) stop()
if(missing(nr)) nr <- max(i)
if(missing(nc)) nc <- max(j)
mat <- array(0, c(nr, nc))
for(ind in seq_along(i)) mat[i[ind], j[ind]] <- x[ind]
return(mat)
}
Try using makeMatrix, and also the standard tcrossprod function from the base package.
i <- c(1,3,7,7)
j <- c(1,1,3,4)
x <- c(1,2,3,4)
A <- makeMatrix(i, j, x)
A## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 0 0 0 0
## [3,] 2 0 0 0
## [4,] 0 0 0 0
## [5,] 0 0 0 0
## [6,] 0 0 0 0
## [7,] 0 0 3 4A %*% t(A) ## what tcrossprod(A) does## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 0 2 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 2 0 4 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 25base::tcrossprod(A) ## standard tcrossprod from base package## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 0 2 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 2 0 4 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 25
Now we’ll do the same tcrossprod calculation using the sparse representation and the Matrix::tcrossprod function.
library(Matrix)
Asparse <- Matrix::sparseMatrix(i, j, x = x)
Asparse## 7 x 4 sparse Matrix of class "dgCMatrix"
##
## [1,] 1 . . .
## [2,] . . . .
## [3,] 2 . . .
## [4,] . . . .
## [5,] . . . .
## [6,] . . . .
## [7,] . . 3 4Matrix::tcrossprod(Asparse)## 7 x 7 sparse Matrix of class "dsCMatrix"
##
## [1,] 1 . 2 . . . .
## [2,] . . . . . . .
## [3,] 2 . 4 . . . .
## [4,] . . . . . . .
## [5,] . . . . . . .
## [6,] . . . . . . .
## [7,] . . . . . . 25
Matrix::tcrossprod from NIMBLENow, we’re ready to build a wrapper function for the Matrix::tcrossprod function, then use nimbleRcall to register this as a nimbleFunction, that can be used inside other compiled nimbleFunctions.
This function below, R_sparse_crossprod is a regular R function wrapping a call to Matrix::tcrossprod. The key point here is that the input arguments are basic numeric vectors and the return value is a numeric array. These “basic” input and output forms are required for use of this function in other nimbleFunctions. But inside the body of this function, any regular R code is fine.
R_sparse_crossprod <- function(i, j, x) {
require(Matrix)
Asparse <- Matrix::sparseMatrix(i, j, x = x)
ans.dsCMatrix <- Matrix::tcrossprod(Asparse)
ans.dgTMatrix <- as(ans.dsCMatrix, 'dgTMatrix')
i <- ans.dgTMatrix@i + 1
j <- ans.dgTMatrix@j + 1
x <- ans.dgTMatrix@x
ijx <- cbind(i, j, x)
return(ijx)
}Check regular R execution of our R function:
ijx <- R_sparse_crossprod(i, j, x)
ijx## i j x
## [1,] 1 1 1
## [2,] 3 1 2
## [3,] 1 3 2
## [4,] 3 3 4
## [5,] 7 7 25makeMatrix(ijx = ijx)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 0 2 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 2 0 4 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 25
Now, we’re ready to call this R function R_sparse_crossprod from within a nimbleFunction.
library(nimble)Register use of the R_sparse_crossprod function, essentially just like any other nimbleFunction, using nimbleRcall. This is where we specify the forms of the (basic numeric) inputs and outputs to our R function.
See help(nimbleRcall).
nimble_sparse_crossprod <- nimbleRcall(
prototype = function(i = double(1), j = double(1), x = double(1)) {},
returnType = double(2),
Rfun = 'R_sparse_crossprod'
)For example, here we’ll define another nimbleFunction (e.g., this could represent a SGV density evaluation) which will internally use nimble_sparse_crossprod in order to call Matrix::tcrossprod.
This nimbleFunction will be R-executable, and also compilable and executable – just like any other nimbleFunction, but this one also happens to call out to arbitrary R code.
Rnf <- nimbleFunction(
run = function(i = double(1), j = double(1), x = double(1)) {
ijx <- nimble_sparse_crossprod(i, j, x)
returnType(double(2))
return(ijx)
}
)
Cnf <- compileNimble(Rnf)## compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compilation details.## compilation finished.Rnf(i, j, x)## i j x
## [1,] 1 1 1
## [2,] 3 1 2
## [3,] 1 3 2
## [4,] 3 3 4
## [5,] 7 7 25Cnf(i, j, x)## [,1] [,2] [,3]
## [1,] 1 1 1
## [2,] 3 1 2
## [3,] 1 3 2
## [4,] 3 3 4
## [5,] 7 7 25makeMatrix(ijx = Rnf(i, j, x))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 0 2 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 2 0 4 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 25makeMatrix(ijx = Cnf(i, j, x))## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 0 2 0 0 0 0
## [2,] 0 0 0 0 0 0 0
## [3,] 2 0 4 0 0 0 0
## [4,] 0 0 0 0 0 0 0
## [5,] 0 0 0 0 0 0 0
## [6,] 0 0 0 0 0 0 0
## [7,] 0 0 0 0 0 0 25
Some brief timing comparisons, between regular R execution of R_sparse_crosspord wrapper for Matrix::tcrossprod, compared against call-back from compiled execution though NIMBLE using the Cnf compiled nimbleFunction.
compare <- function(n) {
set.seed(0)
i <- ceiling(runif(n, 1, n))
j <- ceiling(runif(n, 1, n))
x <- 1:n
R.time <- system.time(R.ans <- R_sparse_crossprod(i, j, x))[3]
C.nimble.time <- system.time(C.nimble.ans <- Cnf(i, j, x))[3]
if(!identical(unname(R.ans), C.nimble.ans)) stop('error')
return(data.frame(n = n,
R.time = as.numeric(R.time),
C.nimble.time = as.numeric(C.nimble.time)))
}
p <- 6
nvals <- 10^(1:p)
as.data.frame(t(sapply(nvals, compare)))## n R.time C.nimble.time
## 1 10 0.001 0.001
## 2 100 0.001 0.001
## 3 1000 0.002 0.001
## 4 10000 0.004 0.005
## 5 100000 0.039 0.048
## 6 1000000 1.608 1.351