# 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.

``````    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?``````

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/

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