Positive semi-definite 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 48×48 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.