In R Is it possible to do matrix operations on matrices that contain functions? -

i working mathematical model functions contained within matrices simulate dynamics of biological populations. requires calling function update matrix, carrying out matrix multiplication project dynamics of system next time step. i'm wondering if possible embed function directly within matrix , skip explicitly updating matrix elements.

the model looks this: state of system in matrix x1, eg

x <- c(0,10) 

the system changes according matrix a

a <- matrix(data = c(0, na,0.75,0.75),nrow =2,byrow = t)      [,1] [,2] [1,] 0.00   na [2,] 0.75 0.75 

where element na function of element 2 in x1, this

f.a.12 <- function(x.2){1 + (1 - x.2/1000)} 

so system simulated this: update marix based on status of x

a[1,2] <- f.a.12(x[2]) 

iterate model updated matrix:

a%*%x       [,1] [1,] 19.9 [2,]  7.5 

this gets iterated 1000 times, updating before each multiplication

the real model uses larger matrices containing multiple functions. there r package allows me embed functions directly within matrix, such as

a.with.fx <- matrix(data = c(0,    f.a.12,                              0.75, 0.75),   nrow =2,byrow = t) 

and carry out regular matrix operations, eg

a.with.fx%*%x 

without having explicitly assign values in function of x in each iteration? guess require function modified %*% operation necessary lookup on own.

i think might trying break problem down in way ends not clean or standard. instead of thinking matrices contain scalar functions, think function returns matrix.

you can trivially do:

update <- function(x) matrix(c(0,.75, f.a.12(x[2]), .75), 2) 

and abuse of notation:

`%**%` <- function(f, x) f(x) %*% x x <- update %**% x 

note %*% s4 generic, it's bit of pain overload it.

alternatively, little algebra can split a constant , function of x, , handle complexity way:

a = [ [  0, 2 - x[2]/1000],       [.75, .75]            ]   = [[0, 2],[.75, .75]] + [[0, - x[2]/1000],[0, 0 ] ]   = c + g(x) 

so update can rewritten:

x = c * x + g(x) * x              

where g matrix-valued function, yielding slower cleaner:

g <- function(x) matrix( c(0,0, -x[2]/1000), 0), 2,2) c <- matrix(c(0,.75,2,.75),2) x <- c(0,10)  > (x <- c %*% x + g(x) %*% x)        [,1]   [1,] 19.9   [2,]  7.5   

we can little better composing g(x) %*% x 1 function :

f <- function(x)  c( -x[2]^2 / 1000, 0   )  

giving:

> x <- c(0,10) > (x <- c %*% x + f(x) )      [,1] [1,] 19.9 [2,]  7.5 

here's couple benchmarks:

microbenchmark(two_step={a <- update(x); %*% x},                   split=(c %*% x + g(x) %*% x),                 composed= c %*% x + f(x)  ) unit: microseconds      expr   min     lq    mean median     uq    max neval  two_step 3.608 4.0900 4.85369 4.3730 4.9290 13.089   100     split 5.587 6.0400 7.72511 6.4900 7.3785 53.047   100  composed 2.266 2.4835 2.78195 2.6815 2.9745  6.697   100