Multiplying matrices with millions of entries

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.


    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?

This is a big matrix, 2.5Gb. I think you will need to use sparse matrices, which requires that most of the cells are zero. Probably there is a package around somewhere for that. Otherwise you will need to decompose the problem somehow.

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.