More efficient way to create Block Toeplitz matrix in R? 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.
 A: 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)
})

A: 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)
}  

