Multiplying matrices with millions of entries

[ashokkrish]

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?

[woodward]

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.

http://www.johnmyleswhite.com/notebook/2011/10/31/using-sparse-matrices-in-r/