last day (14 days later) » 

3:02 AM
2
A: How to calculate a large matrix of zeros and ones in R quickly and efficiently

李哲源 Zheyuan LiA great question!! I will take a step-by-step approach to show you how we can optimize this. Step 1: Vectorization If I have understood you correctly, we can use: outer(t, t, "<=") + 0L ## operator `"<="` returns TRUE / FALSE logical, ## add an integer 0 to make a 0 / 1 binary integer matri...

 
I dont think using sparse matrix will save storage space if it the matrix is not actually sparse - (may be less effcient??)
 
@user20650 If we do both row and column pivoting, we end up with a triangular matrix. This means at least we save half of the storage costs.
 
Sorry, im not sure i understand, what do you mean by "If we do both row and column pivoting, we end up with a triangular matrix". (but if half a sparse matrix is non-zero, then i think this may be more inefficient to store than a full matrix)
 
@user20650 Have a look at this matrix: Z <- outer(tt, tt, "<=") + 0L; it is a lower triangular matrix, with all lower triangular elements being 1. Then Z[piv, piv] gives the desired matrix. I am not saying that we will work this way by full pivoting, but this shows clearly that zeros and non-zeros in the matrix are half-half. So using a sparse matrix saves half of the storage.
@user20650 You are right that a triangular matrix is not really seen as a sparse matrix. And I have not benchmarked my two methods. Maybe I can do it now.
@user20650 shit... it seems that outer is faster
@user20650 Would you like to have a test too, on your machine? Mine is an old one, having been criticised a lot for performance. Try x <- runif(5000); system.time(f(x, TRUE)); system.time(f(x, FALSE)).
 
Re the sparse matrix: if half are zero then it does not mean that it size of the object is half the size. Try object.size(f(t, FALSE)) ; object.size(f(t, TRUE)). Also, a problem i find with outer is it can run into memory problems as the vector increases in length more than a loop will.
 
3:03 AM
It is getting interesting now
 
outer: 0.65 secs, sparse: 3.1 secs. I doubt your lappie is as old as mine!!
 
hahaha, I have 1.61 for outer, 3.5 for sparse
If you use a loop, how would you do it?
You said we can compute half a triangular in a comment, but how?
 
well i hadnt thought that far ahead.....
 
"outer" is known to have memory problem. do you know why? Does it use more storage than just the final matrix?
Well anyway, looks like I have to delete large chunk of my answer
As it turns out doing de-optimization
 
haha.. id keep the answer though, its interesting
 
3:13 AM
@user20650 Well, already updated.
People can still find it in history
I reckon I took one hour on writing that extra bit.
Damn it
So what is going on with matrix size?
I found that triangular matrix only saves 1/3 of storage rather than near 1/2
How did you measure that "outer" takes more memory than a "for" loop?
 
re outer: i seem to remember trying to use it on some large(ish) vectors ./or matrices and hitting memory constraints. I think i then used a loop ... this was all some time ago...
how are you saving the lower.tri matrix? As a vector?
 
No, I was just comparing object.size(f(t, TRUE)) with object.size(f(t, FALSE))
Maybe the "dgCMatrix" class requires some extra costs?
 
ahh..okay, yes i have definitely seen a several things on web about the costs of sparse matrices when they are not sparse
 
Great. I just learnt a lesson
 
3:30 AM
my thought on the calcs.. we only need to calculate the upper or lower, so just use
X <- matrix(0, length(x), length(x))
for(i in 1:5) {
X[i, -(1:i)] <- x[i] <= x[-(1:i)]
}
as we also dont need to calculate the diagonals
 
I am a bit of confused
The matrix is not a symmetric one, why can you arrange computation like this?
 
confused is my default
 
hahaha, well let me actually try it
 
You can arrange like this as if there is a one in the upper then the lower must be zero
to store both seems wasteful
 
With set.seed(0); x <- runif(5)
I don't quite understand the resulting matrix
it is definitely not outer(x, x, "<=") + 0L
 
3:39 AM
its not the same... it only calculates the upper triangle. Which is all that is needed to to fill in the rest of the matrix (if it is required)
Y <- matrix(0, length(x), length(x))
for(i in 1:length(x)){
Y[i, ] <- (x[i] <= x)
}
Y
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 1 0
[2,] 1 1 1 1 1
[3,] 0 0 1 0 0
[4,] 0 0 1 1 0
[5,] 1 0 1 1 1

X <- matrix(0, length(x), length(x))
for(i in 1:5) {
X[i, -(1:i)] <- x[i] <= x[-(1:i)]
}
X
[,1] [,2] [,3] [,4] [,5]
[1,] 0 0 1 1 0
[2,] 0 0 1 1 1
[3,] 0 0 0 0 0
[4,] 0 0 0 0 0
[5,] 0 0 0 0 0
So the upper of X is the same as the upper of Y. We knopw the diagonals are one, so we dont need to calculate these. And we know the lower triangle will be zero when the upper is one, and vice versa so we dont really need to calculate these either
 
Ah I think I got it
 
Just to cut down on the number of calculations that are being done : not storage
 
outer(x, x, "<=") + outer(x, x, ">") gives a matrix will all elements of 1
 
but could actually just save the upper triangle as a vector...
 
Yes, true
 
3:46 AM
okay, cheers for sharing, gtg
 
Why not post as answer?
what is "gtg" short for?
Generally I like using C code for performance. If you implement your idea in C, it will be very efficient
oh.... "gtg" = "got to go"?
Sorry to get you back here.
 

last day (14 days later) »