3

I have an MxM matrix S whose entries are zero on the diagonal, and non-zero everywhere else. I need to make a larger, block matrix. The blocks will be size NxN, and there will be MxM of them.

The (i,j)th block will be S(i,j)I where I=eye(N) is the NxN identity. This matrix will certainly be sparse, S has M^2-M nonzero entries and my block matrix will have N(M^2-M) out of (NM)^2 or about 1/N% nonzero entries, but I'll be adding it to another NMxNM matrix that I do not expect to be sparse.

Since I will be adding my block matrix to a full matrix, would there be any speed gain by trying to write my code in a 'sparse' way? I keep going back and forth, but my thinking is settling on: even if my code to turn S into a sparse block matrix isn't very efficient, when I tell it to add a full and sparse matrix together, wouldn't MATLAB know that it only needs to iterate over the nonzero elements? I've been trained that for loops are slow in MATLAB and things like repmat and padding with zeros is faster, but my guess is that the fastest thing to do would be to not even build the block matrix at all, but write code that adds the entries of (the small matrix) S to my other (large, full) matrix in a sparse way. If I were to learn how to build the block matrix with sparse code (faster than building it in a full way and passing it to sparse), then that code should be able to do the addition for me in a sparse way without even needing to build the block matrix right?

Travis Bemrose
  • 379
  • 1
  • 12
  • How big are `N` and `M`? (Note : I think you'll have `N*(M^2-M)` nonzero entries instead of `M^2-M`. ) – BillBokeey Nov 02 '15 at 08:51
  • If only the diagonal is zero then that's not very sparse. There is an extra overhead associated with a sparse matrix that will cause it to require more memory in such cases. – IKavanagh Nov 02 '15 at 09:37
  • [This answer](http://stackoverflow.com/a/3293725/2433501) suggests that you need to be allocating under 25% of the sparse matrix to see an improvement in performance. But in reality, I think that the only way to be sure if it would make a difference in your situation is to implement a minimal version of whatever algorithm you are running and time it. Without knowing the exact operations that you are performing on your matrix I think that it is impossible for us to provide a conclusive answer. – zelanix Nov 02 '15 at 10:03
  • As they already said here: If your number of non-zeros is bigger than 25% of the entries, dont bother with sparse. – Ander Biguri Nov 02 '15 at 12:09
  • @BillBokeey Ohp... you're right. `S` has `M^2-M` nonzero entries, which will correspond 1 to N with the nonzero entries of my new matrix. At first I wanted to be able to handle up to `N~1,000` and `M~10,000`, but a quick calculation says that my full `NMxNM` matrix will have `1e14` entries, which would be `~100 TB`. – Travis Bemrose Nov 02 '15 at 20:16
  • @IKavanagh `S` is not sparse, but when I take those entries and spread them out over the larger matrix, it will be sparse. – Travis Bemrose Nov 02 '15 at 20:18
  • My block matrix will certainly have less than 25% nonzero entries. (Even with `N=4` and `M=6` it's 21%, and as I scale up, the percentage will get smaller and smaller). The question isn't how sparse is my block matrix. @zelanix the operation I'm doing is that I'll then add this sparse block matrix to a full `NMxNM` matrix. I can mess with some toy problems, but before I spent time studying how to write efficient sparse code, I was wondering if there's any advantage because it may not gain me anything. – Travis Bemrose Nov 02 '15 at 20:29
  • Well, I think regardless of performance, if you're seriously thinking about a matrix with 1e14 values, which, incidentally, will be much bigger than 100 TB ([this page](http://uk.mathworks.com/help/matlab/matlab_prog/memory-allocation.html) indicates 8 bytes per numeric element) you will have no choice but to create a sparse matrix! – zelanix Nov 02 '15 at 21:03
  • @zelanix Yeah, I knew it'd be bigger, which is why I used the ~ for 'on the order of'. The takeaway is that I'll have to be satisfied with much smaller cases, because I'm adding this sparse block matrix to a FULL matrix of that size. There's just no way to get around that (though if this does what I want, I'll certainly look for ways around it). – Travis Bemrose Nov 02 '15 at 22:25
  • I'm afraid that's just not going to be possible - you will need to re-engineer this. Matlab just can't create 1e7 x 1e7 matrices. See [here](http://uk.mathworks.com/matlabcentral/answers/91711-what-is-the-maximum-matrix-size-for-each-platform) for details - you're limited to 8TB if you have the available memory, including swap space. If you're creating 800TB of data structures (1e7 x 1e7 x 8) you will need 800TB of storage space too! – zelanix Nov 02 '15 at 22:45

2 Answers2

1

If you can keep a full NMxNM matrix in memory, dont bother with sparse operations. In fact in most cases A+B, with A full and B sparse, will take longer than A+B, where A and B are both full.

Felix Darvas
  • 507
  • 3
  • 5
1

From your description, using sparse is likely slower for your problem:

If you're adding a sparse matrix A to a full matrix B, the result is full and there's almost certainly no advantage to having A sparse.

For example: n = 12000; A = rand(n, n); B1 = rand(n, n); B2 = spalloc(n, n, n*n); B2 is as sparse as possible, that is, it's all zeros! On my machine, A+B1 takes about .23 seconds while A + B2 takes about .7 seconds.

Basically, operations on full matrices use BLAS/LAPACK library calls that are insanely optimized. Overhead associated with sparse is going to make things worse unless you're in the special cases where sparse is super useful.

When is sparse super useful?

Sparse is super useful when the size of matrices suggest that some algorithm should be very slow, but because of sparsity (+ perhaps special matrix structure), the actual number of computations required is orders of magnitude less.

EXAMPLE: Solving linear system A*x=b where A is block diagonal matrix: As = sparse(rand(5, 5)); for(i=1:999) As = blkdiag(As, sparse(rand(5,5))); end %generate a 5000x5000 sparse block diagonal matrix of 5x5 blocks

Af = full(As); b = rand(5000, 1);

On my machine, solving the linear system on the full matrix (i.e. Af \ b) takes about 2.3 seconds while As \ b takes .0012 seconds.

Sparse can be awesome, but it's only helpful for large problems where you can cleverly exploit structure.

Matthew Gunn
  • 4,451
  • 1
  • 12
  • 30