Sequentially multiply matrices from a list by a vector at time t-1 (recursively)

Question

I am trying to multiply a list of matrices, in the order they appear within the list starting with matrix 1, by an initial vector and then recursively; so matrix 2 from the list multiplied by that resulting vector. I have tried various iteration of lapply and map but am unable to project forward and perform this recursively. More explicitly: A[[1]] % * % allYears[,1], then A[[2]] % * % allYears[,2],.....,A[[4]] % * % allYears[,4] which produces the final 5th column in "allYears". Below is the sample code with the known error in the for loop at A[[i]] indexing as i is not explicitly referenced.

A <- lapply(1:4, function(x)  # construct list of matrices
  matrix(c(0, 0, 10, 
           rbeta(1, 5, 4), 0, 0, 
           0, rbeta(1, 10, 2), 0), nrow=3, ncol=3, byrow=TRUE, ))

n <- c(1000, 100, 10)  # initial vector of abundances

nYears <- 4  # define the number of years to project over

allYears <- matrix(0, nrow=3, ncol=nYears+1)  # build a storage array for all abundances

allYears[, 1] <- n  # set the year 0 abundance   

for (t in 1:(nYears + 1)) {   # loop through all years
  allYears[, t] <- A[[i]] %*% allYears[, t - 1]
}

I have corrected to nYears +1. Thanks. – user11057680 Jan 20 '22 at 17:26 — user11057680, Jan 20 '22 at 17:26

score 3 · Accepted Answer · answered Jan 20 '22 at 17:24

3

Based on the description, perhaps we need to loop over the sequence - i.e. length of A is 4, whereas the number of columns of 'allYears' is 5. Create an index from 2 to ncol of 'allYears', then loop over the sequence of that index, extract the corresponding element of 'A' based on the sequence whereas we get the allYears previous column

i1 <- 2:(nYears + 1)
 for(t in seq_along(i1)) {
    allYears[,i1[t]] <-  A[[t]] %*% allYears[,i1[t]-1]
 
 }

-output

> allYears
     [,1]      [,2]      [,3]       [,4]      [,5]
[1,] 1000 100.00000 817.24277 2081.08322  333.6702
[2,]  100 261.46150  55.44237  423.22095 1244.6680
[3,]   10  81.72428 208.10832   33.36702  355.5175

answered Jan 20 '22 at 17:24

akrun

874,273
37
540
662

1

that would be it, thank you kindly. – user11057680 Jan 20 '22 at 17:29
I'm not sure this is right? `A[[1]] %*% n` is (100, 569, 78) ... – Ben Bolker Jan 20 '22 at 18:53
@BenBolker There was no `set.seed` in OP's example. So, the values you got may be different – akrun Jan 20 '22 at 18:54
1

aha. Thanks .... – Ben Bolker Jan 20 '22 at 18:55

score 2 · Answer 2 · edited Jan 20 '22 at 20:29

2

Alternative, possibly too clever, solution:

Construct a list of (A[[1]], A[[1]] %*% A[[2]], ... )

Alist <- Reduce("%*%", A, accumulate=TRUE)

Multiply each of these by the initial value

vlist <- lapply(Alist, "%*%", n)

Combine:

do.call(cbind, vlist)
          [,1]      [,2]       [,3]     [,4]
[1,] 100.00000 856.66558 4864.20044  486.420
[2,] 569.23739  56.92374  543.02451 3307.690
[3,]  78.55101 553.62548   55.36255  445.619

@MikaelJagan points out that this can be done in fewer steps:

do.call(cbind, 
        rev(Reduce("%*%", rev(A), init = n, right = TRUE, accumulate = TRUE)))

or (in recent versions of R)

(A
  |> rev()
  |> Reduce(f = "%*%", init = n, right = TRUE, accumulate = TRUE)
  |> rev()
  |> do.call(what = cbind)
)

(The last step could be replaced by |> unlist() |> matrix(nrow = length(n)).)

edited Jan 20 '22 at 20:29

Mikael Jagan

9,012
2
17
48

answered Jan 20 '22 at 18:57

Ben Bolker

211,554
25
370
453

Maybe you want `matrix(unlist(Reduce("%*%", A, init = n, accumulate = TRUE)), nrow = length(n))`. – Mikael Jagan Jan 20 '22 at 19:40
You're too fast for me! – Ben Bolker Jan 20 '22 at 19:43
I have a feeling I've made a mistake, though... – Mikael Jagan Jan 20 '22 at 19:48
I think I've got it now... I've written this comment three times... By default, `Reduce` does `init %*% A[[i]]`. In order to get what we want, we need to pass `right = TRUE` _and_ reverse the order of the input and output. Hence `matrix(unlist(rev(Reduce("%*%", rev(A), init = n, right = TRUE, accumulate = TRUE))), nrow = length(n))`. – Mikael Jagan Jan 20 '22 at 20:19

Sequentially multiply matrices from a list by a vector at time t-1 (recursively)

2 Answers2

Linked