Memory efficient way to zero out the diagonal of a sparse matrix in R

user_4476155 0
Viewed 256

I'd like to zero out the diagonal of a sparse matrix in R. My brute force way is explicitly setting it to zero, but this seems inefficient. Is there a more efficient way?

require(Matrix)
A <- as(rsparsematrix(nrow = 1e7, ncol = 1e7, nnz = 1e4), "sparseMatrix")
diag(A) <- 0
A <- drop0(A)  # cleaning up

Clarification and resolution: my initial worry was that Matrix inflates the sparse matrix with actual zeros on the diagonal. This turns out not to be the case (in the end, although in the interim it is, see comment below). To see this, consider what would happen if we were to set the diagonal to one:

A <- as(rsparsematrix(nrow = 1e7, ncol = 1e7, nnz = 1e4), "sparseMatrix")
format(object.size(A), units = "Mb")

[1] "38.3 Mb"

diag(A) <- 1
format(object.size(A), units = "Mb")

[1] "152.7 Mb"

The many non-zero elements we have added use up O(n) memory, where n is the dim of the matrix. However, with diag(A) <- 0 we get:

diag(A) <- 1
format(object.size(A), units = "Mb")

[1] "38.3 Mb"

Namely, Matrix already handles this situation efficiently.

sparse-matrix r
1 Answers

You can find the non-zero entries really quickly:

ij <- which(A != 0, arr.ind = TRUE)

# Subset to those on the diagonal:

ij <- ij[ij[,1] == ij[,2],,drop = FALSE]

# And set those entries to zero:

A[ij] <- 0

Edited to add:

As the revision to the original question says, this doesn't save much memory in the end, but it is much faster. The diag(A) <- 0 statement takes about 3.2 seconds on my computer, whereas these 3 steps take about 0.2 seconds. Here's how to do the timing:

library(microbenchmark)
microbenchmark(A <- as(rsparsematrix(nrow = 1e7, ncol = 1e7, nnz = 1e4), "sparseMatrix"),
{A <- as(rsparsematrix(nrow = 1e7, ncol = 1e7, nnz = 1e4), "sparseMatrix"); diag(A) <- 0},
{A <- as(rsparsematrix(nrow = 1e7, ncol = 1e7, nnz = 1e4), "sparseMatrix");ij <- which(A != 0, arr.ind = TRUE);ij <- ij[ij[,1] == ij[,2],,drop = FALSE];A[ij] <- 0}, times = 10)

When I run it, I see median timings of 137 ms for the matrix creation and nothing else, 3351 ms for creation plus the diag(A) call, and 319 ms for creation followed by my code.

It also saves a lot of memory in intermediate steps, which can be seen using memory profiling: Rprof(memory=TRUE); run code ; Rprof(NULL); summaryRprof().

user_25543300
Related