I have a PEPS tensor network where each tensor has 5^4 =625 entries, of which only 15 are nonzero. With 100 tensors in the network, there’d be 5^{100} entries in the contracted network, which is much more memory than I have at hand. However, I expect the vast majority of entries in the final contracted tensor to be 0, with less than 10^5 nonzero values.
While doing the contraction, I currently run out of memory much faster than I run out of computation time. It seems like I should be able to use the sparsity of my tensors to at least avoid running out of memory, even if I’m unable to utilize the sparsity to speed up the calculation.
How is ITensor’s sparse array functionality?
I know from this post that block sparse arrays are already supported with generation via quantum numbers, and that in the future general block sparse arrays will be available, but my tensor doesn’t have an obvious block sparse structure (although with so many 0s, it’s possible that there’s some opportunistic block sparsity that could be leveraged).
I know there is a SparseArraysBase with a DOK-based sparse array under development, but it seems like it isn’t integrated with the rest of the ecosystem yet. In addition, using this would save on memory but wouldn’t save on computation time (if I understand DOKs properly).
Julia’s SparseArrays library doesn’t support multidimensional arrays as far as I can tell, and it looks like ITensor needs access to the internal details of the arrays underlying tensors anyways, so it seems like the best option right now is to figure out how to use SparseArrayDOKs. Using some sort of SVD truncation scheme could be possible given that I have a finite system with nonperiodic boundary conditions - I’ve heard of treating one of the boundary rows as an MPS and then marching that state across the 2D lattice, truncating each time I contract in a new row. However, with so few entries it seems like I might actually sacrifice a lot of accuracy with that approach.
Are there other considerations I should keep in mind? What would the recommended approach be? Thanks in advance.
Good question. @JoeyT1994 was recently trying something like that (contracting a PEPS with very sparse tensors) using the SparseArrayDOK type defined in SparseArraysBase.jl, I think it worked but in the end for his use case it didn’t improve over doing standard boundary MPS contraction using dense tensors. You can use it as a backend for ITensorBase.jl, which is our “next-gen” replacement for ITensors.jl which is basically built to be able to easily wrap our new array libraries including SparseArraysBase.jl (or any Julia AbstractArray). However, use those new packages at your own risk, they are still under active development and are primarily for internal use right now. Note that tensor contractions of SparseArrayDOK could use some optimizing (for example we matricize the tensors and perform matrix multiplications, but don’t take advantage of optimized sparse matrix multiplication backends for that, and additionally it could be helpful to call out to sparse tensor libraries like Finch.jl but we haven’t investigated that very much), and also we don’t have tensor/matrix factorizations defined for SparseArrayDOK yet. So if you go that route, you may have to add those features yourself or wait until we do it if you aren’t able to do it yourself.
In terms of algorithms for dealing with networks like that, I’m not aware of much work on that. You can of course try standard methods like boundary MPS to approximately contract the network, and maybe use something like a rank revealing sparse QR to perform the truncations (I think SVD wouldn’t work since I’m not aware of a sparse SVD, though maybe one exists). I think that is something that would be considered an area of research where the best approaches aren’t yet known.
Thanks for the prompt and helpful response! When you say that you “matricize the tensors and perform matrix multiplications, but don’t take advantage of optimized sparse matrix multiplication backends for that”, does this mean:
that while the ‘latent’ storage of arrays in memory might be sparse, they get converted into a dense form before and after contractions? If that is the case, it seems like I will hit a bottleneck at the final contraction no matter what, and am better off trying boundary MPS.
that you keep the DOKs storage type for the entire computation, but during the computation itself you stubbornly iterate through each term in the matrix product rather than intelligently skipping rows or columns where possible (as you could with a CSR-type sparse array)
I imagine that you mean the second one, but just making sure.
Yes, I mean the second one. It would be very easy to convert to a sparse matrix CSC/CSR format when we matricize, which would accelerate the matrix multiplication, we just haven’t implemented it yet.
Yeah this is an interesting question. I did have a similar tensor network (TN) at some point (full of sparse tensors) that I wanted to contract. In my case the final result was just a scalar (no open indices on the TN).
How to do this is very much an open research question tbh.
I tried two tactics:
Working in a dense format. I Gauged the TN so that the tensors were no longer sparse but the overall TN was unchanged then I ran boundary MPS from left to right with truncation to get an approximate contraction.
I implemented my own “hacky” book-keeping for contracting sparse ITensors and started contracting the tensor network with this book-keeping.
Both schemes worked but the book-keeping still got out of hand for the sparse approach because the intermediate tensors were sparse but still had like 10^10 entries. If you think this wouldn’t be the case (intermediate contractions have < 10^5 entries) this could work well.
The dense format worked but the answer became approximate, which means if I wanted to resolve the contraction to high precision I had difficulties