Generate a large block diagonal covariance matrix with exponential decay

Question

I am implementing Kalman filtering in R. Part of the problem involves generating a really huge error covariance block-diagonal matrix (dim: 18000 rows x 18000 columns = 324,000,000 entries). We denote this matrix Q. This Q matrix is multiplied by another huge rectangular matrix called the linear operator, denoted by H.

I am able to construct these matrices but it takes a lot of memory and hangs my computer. I am looking at ways to make my code efficient or do the matrix multiplications without actually creating the matrices exclusively.

library(lattice)
library(Matrix)
library(ggplot2)

nrows <- 125

ncols <- 172

p <- ncols*nrows

#--------------------------------------------------------------#
# Compute Qf.OSI, the "constant" model error covariance matrix #
#--------------------------------------------------------------#

  Qvariance <- 1
  Qrho <- 0.8

  Q <- matrix(0, p, p) 

  for (alpha in 1:p)
  {
    JJ <- (alpha - 1) %% nrows + 1
    II <- ((alpha - JJ)/ncols) + 1
    #print(paste(II, JJ))

    for (beta in alpha:p)
    {
      LL <- (beta - 1) %% nrows + 1
      KK <- ((beta - LL)/ncols) + 1

      d <- sqrt((LL - JJ)^2 + (KK - II)^2)
      #print(paste(II, JJ, KK, LL, "d = ", d))

      Q[alpha, beta] <-  Q[beta, alpha] <-  Qvariance*(Qrho^d)
    }
  } 

  # dn <- (det(Q))^(1/p)
  # print(dn)

  # Determinant of Q is 0
  # Sum of the eigen values of Q is equal to p

  #-------------------------------------------#
  # Create a block-diagonal covariance matrix #
  #-------------------------------------------#

  Qf.OSI <- as.matrix(bdiag(Q,Q))

  print(paste("Dimension of the forecast error covariance matrix, Qf.OSI:")); print(dim(Qf.OSI))

It takes a long time to create the matrix Qf.OSI at the first place. Then I am looking at pre- and post-multiplying Qf.OSI with a linear operator matrix, H, which is of dimension 48 x 18000. The resulting HQf.OSIHt is finally a 48x48 matrix. What is an efficient way to generate the Q matrix? The above form for Q matrix is one of many in the literature. In the below image you will see yet another form for Q (called the Balgovind form) which I haven't implemented but I assume is equally time consuming to generate the matrix in R.

Answer 1

    JJ <- (alpha - 1) %% nrows + 1
    II <- ((alpha - JJ)/ncols) + 1

That looks likely to be buggy. I would guess that a is supposed to be an encoding for a pair (row, col), but in that case the same base should be used for the %% and the /.

I would also suggest that if you can't use 0-indexed matrices then you do the offset to 1-based when you access the matrices, and keep the values you manipulate 0-based. See how much simpler this is:

for (rowa in 0:(nrows-1))
{
  for (cola in 0:(ncols-1))
  {
    a = rowa * ncols + cola
    for (rowb in 0:(nrows-1))
    {
      for (colb in 0:(ncols-1))
      {
        b = rowb * ncols + colb
        d = sqrt((rowa - rowb)^2 + (cola - colb)^2)
        Q[a+1, b+1] <- Qvariance * (Qrho^d)
      }
    }
  }
}

Incidentally, since Qvariance is multiplied into every single element you could pull that out and post-multiply the final \$48 \times 48\$ matrix instead.

---

Now, elimination of the matrix. We have \$(AB)_{i,j} = \sum_k A_{i,k} B_{k,j}\$, so $$(HQH^T)_{i,j} = \sum_k H_{i,k}(QH^T)_{k,j} = \sum_k H_{i,k} \sum_l Q_{k,l} H^T_{l,j} = \sum_k \sum_l H_{i,k} H_{j,l} Q_{k,l}$$ which allows you to restructure the code so as to avoid creating \$Q\$ in memory. However, it is at the cost of using the naïve algorithm for matrix multiplication, and your matrices are large enough that R is probably using a sub-cubic algorithm. So what you might want to do is to instead break it down into chunks: e.g. of size nrows \$\times\$ nrows. I don't know enough R to be certain, but I expect that its index range notation allows you to do this quite cleanly.

Following up on some comments, we can expand \$k = r_1 C + c_1\$, \$l = r_2 C + c_2\$ where \$C\$ is ncols, and get $$(HQH^T)_{i,j} = \sum_{r_1} \sum_{c_1} \sum_{r_2} \sum_{c_2} H_{i,r_1 C + c_1} H_{j,r_2 C + c_2} Q_{r_1 C + c_1,r_2 C + c_2} \\ = \sigma \sum_{r_1=1}^R \sum_{r_2=1}^R \sum_{c_1=1}^C \sum_{c_2=1}^C H_{i,r_1 C + c_1} H_{j,r_2 C + c_2} \rho^{\sqrt{(r_1-r_2)^2 + (c_1-c_2)^2}} $$

Let \$Q^{(\delta)}\$ be a symmetric \$C \times C\$ matrix with \$Q^{(\delta)}_{i,j} = \rho^{\sqrt{\delta^2 + (i-j)^2}}\$. Then $$(HQH^T)_{i,j} = \sigma \sum_{r_1=1}^R \sum_{r_2=1}^R \sum_{c_1=1}^C \sum_{c_2=1}^C H_{i,r_1 C + c_1} Q^{(|r_1-r_2|)}_{c_1,c_2} H_{j,r_2 C + c_2} \\ HQH^T = \sigma \sum_{r_1=1}^R \sum_{r_2=1}^R H_{1..48,r_1 C .. (r_1+1)C} Q^{(|r_1-r_2|)} H_{1..48,r_2 C..(r_2+1)C}^T $$

and the sum can be regrouped by \$|r_1 - r_2|\$ to calculate each \$Q^{(\delta)}\$ only once. When calculating \$Q^{(\delta)}\$ you can exploit the symmetry without worrying too much about cache coherence, because the whole of \$Q^{(\delta)}\$ should fit in L2 cache.

Answer 2

Maybe if you have enough RAM this will be faster: (only the Q creation)

alpha <- matrix(rep(1:p, p), p, p)
JJ <- (alpha - 1L) %% nrows + 1L
II <- ((alpha - JJ)/ncols) + 1L
LL <- t(JJ)
KK <- t(II)
d <- sqrt((LL - JJ)^2 + (KK - II)^2)
Q2 <- Qvariance*(Qrho^d)
all.equal(Q, Q2)
# TRUE