Generate a large block diagonal covariance matrix with exponential decay

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.

• I'm voting to close this question as off-topic because it belongs on cs.stackexchange.com – dfhwze Jul 20 at 8:15
• I think this question should go to SO. @dfhwze there is no [r] tag in cs.stackexchange.com I dont think it will be very useful. – Grada Gukovic Jul 21 at 3:31
• @GradaGukovic Initially I posted this on SO they sent me over to CodeReview. There someone voted to my question here so I have the same question posted on 3 locations. Not sure where is the right place for this question. Any help with code is highly appreciated. – Ashok Jul 21 at 5:44
• Rather than toggling between sites, perhaps it should stay then. Hopefully some expert in this subject can help you out. – dfhwze Jul 21 at 8:02

    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.

• Thank you Peter Taylor for your code and insight! The code works great. I just have a quick question. The Q matrix is symmetric so essentially Q[a, b] = Q[b, a]. I wonder if we can save time in not generating the entire matrix rather generate only the lower triangular matrix and copy it to the upper triangular matrix. Would appreciate your comments or code for this? Much appreciated. – Ashok Jul 23 at 20:05
• @Ashok, the only way to know whether that would be faster is to try it. I suspect that it would be slower, because the inner loop only has about 6 floating point operations, but accessing both triangles of the matrix at the same time wouldn't be good for cache coherence. Also, in case it wasn't clear, the code was to make a point about readability. If you break the calculation into chunks of ncols by ncols then I think each chunk is symmetric and can be used multiple times, but should fit into L2 cache. – Peter Taylor Jul 24 at 6:20

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

• Thank you @minem the code works great! – Ashok Jul 29 at 21:29