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 (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)
print(paste("a = ", a+1, "b = ", b+1, "d = ", 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))
# The below line of code generates a plot using the 'lattice' package.
# For small values of p uncomment and run to see the plot
# levelplot(Q, col.regions= colorRampPalette(c("red","green","blue")))
# levelplot(Qf.OSI, col.regions= colorRampPalette(c("red","green","blue")))
# The below three lines of code generates a plot using the 'ggplot2' package.
# For small values of p uncomment and run to see the plot
# dd <- expand.grid(x = 1:ncol(Qf.OSI), y = 1:nrow(Qf.OSI))
# dd$col <- unlist(c(Qf.OSI))
# ggplot(dd, aes(x = x, y = y, fill = factor(col))) + geom_tile()
'''
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 H*Qf.OSI*Ht is finally a 48x48 matrix. What is an efficient way to generate the Q matrix?
```