Issues mixing sparse and dense matrices

[James-Thorson-NOAA]

I'm getting a variety of confusing error messages when exploring sparse matrices in RTMB.

For example:

• Diagonal(x=advector(1:3)) fails while diag(x=advector(1:3)) works
• I had some issues with elementwise product of a dense and sparse matrix (not shown)
• I had some issues with Matrix( i=, j=, x=) or Matrix( i=, p=, x=) giving an error message about an unsupported class dTMatrix even when passing argument repr = "C"
• I can't seem to do M = sparseMatrix(...); M@x <- vec assignment in package functions;
• The script below has output "Error: evaluation nested too deeply: infinite recursion / options(expressions=)?"

I'm sorry that I don't have better intuition about how to diagnose these class and coersion errors in RTMB. Any thoughts?

library(Matrix)
library(RTMB)
set.seed(1)
n <- 4
Msparse <- rsparsematrix(n, n, nnz = floor(n^(3/2)) )
x = rnorm(n)
Mdense = matrix( rnorm(n^2), n, n)

f <- function(x) {
  Mprod = Mdense * Msparse
  sum((Mprod %*% x)^2)
}
F <- MakeTape(f, x )
obj = MakeADFun(f, x)

[kaskr]

@James-Thorson-NOAA Thanks for pointing out these issues! It would be a huge work to make all Matrix functions work in RTMB. My original intention was to add a general converter from Matrix to RTMB - it's now added (AD) with demo (example(AD)). So you can do

• AD(Diagonal(3)) and then assign the x slot
• Elementwise products are now fixed
• Matrix and sparseMatrix can be used as well if wrapped in AD and then assigned to later.
• The infinite recursion in your example is actually caused by a method in the Matrix package. We could probably come up with a fix to never go into that method, but this is not top priority as you could simply add Mprod = as.matrix(Mprod):

selectMethod("%*%", c("dgeMatrix", "advector"))
Method Definition:

function (x, y) 
x %*% (if (length(dim(y)) == 2L) as.matrix else as.vector)(y)
<bytecode: 0x55f8f1de2168>
<environment: namespace:Matrix>

Signatures:
        x           y         
target  "dgeMatrix" "advector"
defined "Matrix"    "ANY"

[James-Thorson-NOAA]

Thanks for looking! This seems very promising, and I can confirm that the demo(AD) is working as stated.

However, I am still encountering some errors that I can replicate by modifying my original script as follows:

f <- function(x) {
  # Error in Diagonal(n, x = x) : 'x' has unsupported class "advector"
  #D = AD(Diagonal(n, x=x))

  # WORKS
  D = AD(Diagonal(n=n))
  D@x = x

  # WORKS
  #M = D

  # ERROR for either option: Error in as.vector(data) : no method for coercing this S4 class to a vector
  M = D %*% Msparse
  #M = crossprod(D, Msparse)

  sum(M %*% x)
}
F <- MakeTape(f, x )
obj = MakeADFun(f, x)

FWIW, I'm trying to explore converting EcoState to use a sparse diet matrix, to hopefully speed up the 100s to 1000s of dBdt function calls when taping the ODE solver.

[James-Thorson-NOAA]

OK, after more playing around, I can see that I needed another AD coersion:

f <- function(x) {
  # Error in Diagonal(n, x = x) : 'x' has unsupported class "advector"
  #D = AD(Diagonal(n, x=x))

  # WORKS
  D = AD(Diagonal(n=n))
  D@x = x

  # WORKS
  #M = D

  # ERROR for either option: Error in as.vector(data) : no method for coercing this S4 class to a vector
  M = D %*% AD(Msparse)

  sum(M %*% x)
}
F <- MakeTape(f, x )
obj = MakeADFun(f, x)

I'm gonna keep playing around, and closing the issue for now. Sorry for the extra messages!