6
$\begingroup$

I am trying to create a Block Toeplitz matrix likes this: $\boldsymbol U = \begin{pmatrix} U(0) & U(1) & \cdots & U(n)\\ U(1)^T & U(0) & \cdots & U(n-1)\\ \vdots & \vdots & \vdots & \vdots\\ U(n)^T & U(n-1)^T & \cdots & U(0) \end{pmatrix}$ where $U(i)$, $i = 1, \dots, n$ is a $k \times k$ matrix.

I found the toeplitz() function in R only deals with the case when $U(i)$s are scalars. Currently I only know to use loops to create the matrix $\boldsymbol U$, which is not efficient in R. If someone have any suggestions to vectorize the code, could you please help me out.

For example, I am interested in creating a matrix like this:

Let $U(0) = \begin{pmatrix} 0.8 & 1 \\ 0 & 0.8 \end{pmatrix}$, $\Phi = \begin{pmatrix} 0.9 & 1 \\ 0 & 0.8 \end{pmatrix}$, and let $U(i) = \Phi^i U(0)$ for $i = 1, \dots, n$.

The R code I wrote is,

U0 <- matrix(c(0.8, 0, 1, 0.8), 2, 2) # create U0
phi <- matrix(c(0.9, 0, 1, 0.8), 2, 2) # create phi
n <- 100 # n is the number of U(i)
k <- dim(U0)[1]
U <- matrix(NA, k*n, k*n) # create an empty matrix to receive value
for (i in 1:n) {
   for (j in 1:n) {
     if (i >= j) {
       U[((i-1)*k+1):(i*k), ((j-1)*k+1):(j*k)] <- phi^(i-j) %*% U0
     } else { 
       U[((i-1)*k+1):(i*k), ((j-1)*k+1):(j*k)] <- t(phi)^(j-i) %*% U0
     }
   }
 }

Is there anyway to avoid the loops?

Thanks.

$\endgroup$
5
$\begingroup$

This question might be of interest to statisticians and data scientists because of the multidimensional array handling issues it presents.

Loops aren't the problem. Your goal likely is to perform the operation efficiently, regardless of how it might be done, and perhaps to do it in an obviously parallelizable way.

Exploit R's ability to manage multidimensional arrays. (It inherits this from FORTRAN, which even in its original design handled up to seven dimensions.) Move groups of data en masse into an array and then reshape it as needed. These native operations will be fast, even if some looping needs to be done.

Although it uses short high-level loops (each executing $n+1$ times), the following solution will create a "block Toeplitz" array with $n=999$ and $2\times 2$ blocks--an array with four million elements--within a half second and small use of overhead storage. It is based on the observation that the Toeplitz pattern can be constructed from the string

$$U_n, U_{n-1}, \ldots, U_1, U_0, U_1^\prime, U_2^\prime, \ldots, U_n^\prime$$

by taking the last $n+1$ elements to form the first column of output, shifting those elements to the left to form the second column of output, and so on, until the last column of output is formed from the first $n+1$ elements of this string.

The code uses n to represent $n+1$ (because R indexing starts at 1 instead of 0). The string of matrices itself, as a $k\times (2n-1)$ array, is stored in a temporary array strip.

At the end, the $k\times k\times n\times n$ array that is created by "blasting" parts of the strip down the columns of the output X is permuted to get its values stored in the right order (I was too lazy to work out the optimal storage configuration at the outset) and then recast as a $kn\times kn$ array (which involves no movement of data). These column-wise operations are each highly vectorized and capable of massive parallelism.

The initial value of U contains obvious patterns to check that the output is correct.

#
# Create n square matrices of dimension k by k.
#
U <- list(matrix(1:4, 2), matrix(5:8, 2), matrix(9:12, 2))
#U <- lapply(1:1000, function(i) matrix(-3:0 + 4*i, 2))

system.time({
  k <- min(unlist(lapply(U, dim)))
  n <- length(U)
  #
  # Create the "strip".
  #
  strip <- array(NA, dim=c(k,k,2*n-1))
  for (i in 1:n) strip[,,i] <- U[[n+1-i]]
  if (n > 1) for (i in 2:n) strip[,,n+i-1] <- t(U[[i]])
  #
  # Assemble into "block-Toeplitz" form.
  #
  X <- array(NA, dim=c(k,k,n,n))
  #
  # Blast the strip across X.
  #
  for (i in 1:n) X[,,,i] <- strip[,,(n+1-i):(2*n-i)]
  X <- matrix(aperm(X, c(1,3,2,4)), n*k)
})
$\endgroup$
1
$\begingroup$

The following function that takes as argument a list of blocks. However, this may be not the best solution: there still one loop present. And it needs some more work, since it doesn't do the transpose of the blocks lower diagonal (in my case I have symmetric matrices).

toeplitz.block <- function(blocks) {
    l <- length(blocks)
    m.str <- toeplitz(1:l)

    res <- lapply(1:l,function(k) {
        res <- matrix(0,ncol=ncol(m.str),nrow=nrow(m.str))
        res[m.str == k] <- 1
        res %x% blocks[[k]]
    })

    Reduce("+",res)
}  
$\endgroup$
  • 2
    $\begingroup$ Welcome. There is no penalty for answering old questions here, so well done! There is also no need for apologies (which I have edited out). The philosophy of the site is to cut to the chase and provide answers to questions, you can have a look at the tour to familiarise yourself with the philosophy of the site. $\endgroup$ – Antoine Vernet Jul 26 '16 at 14:42
  • 4
    $\begingroup$ Note that loops are not necessarily inefficient in R as the apply family are just loops internally. $\endgroup$ – mdewey Jul 26 '16 at 14:56
  • 1
    $\begingroup$ @mdewey:+1. Especially apply functions can routinely be as slow as or even slower than for-loops. C wrapper functions like colMeans etc. can make a world difference but standard apply's are not as efficient as many people think. $\endgroup$ – usεr11852 says Reinstate Monic Jul 26 '16 at 15:13
  • 1
    $\begingroup$ @usεr11852, that's right. I use *apply mostly because they're concise & easier to read & follow later. $\endgroup$ – gung - Reinstate Monica Jul 26 '16 at 16:17
  • $\begingroup$ @gung: I completely agree; this why people should primarily use them. My comment was just a slight expansion on @mdewey's. I routinely use *apply functions too. $\endgroup$ – usεr11852 says Reinstate Monic Jul 26 '16 at 17:39

Not the answer you're looking for? Browse other questions tagged or ask your own question.