Speed and scalability improvements for graph multiresolution #3
Conversation
|
Hi, thanks for the contribution. Could you check and correct the errors that were raised in the TravisCI checks? They seem to be mostly errors in the code examples in the documentation, raising errors such as "NameError: name 'sparse' is not defined", "NameError: name 'extract_submatrix' is not defined", or "NameError: name 'block' is not defined". |
|
There's still a silly error in the doctest (it's one that I often forget is a problem): pygsp/pygsp/utils.py", line 261, in utils.py There should be a newline after (8,8) in the documentation |
|
My bad, it should be fine now. |
| if self.directed or not self.connected: | ||
| raise NotImplementedError('Focusing on connected non directed graphs first.') | ||
|
|
||
| start_nodes, end_nodes, weights = sparse.find(sparse.tril(self.W)) |
There was a problem hiding this comment.
Do you know if this is faster/slower/same as what's done in adj2vec(), lines 23-24 in pygsp/data_handling.py? Indeed, I rather prefer your design, than the use of adj2vec we're making now for computing graph gradients. I might make the necessary adaptions to replace it by this call of your create_incidence_matrix soon.
pygsp/operators/reduction.py
Outdated
| @@ -12,17 +12,26 @@ | |||
| logger = build_logger(__name__) | |||
|
|
|||
|
|
|||
| def graph_sparsify(M, epsilon, maxiter=10): | |||
| def graph_sparsify(M, epsilon, maxiter=10, fast=True): | |||
There was a problem hiding this comment.
I would set fast=False by default, so that the previous default behavior of this function remains unchanged. What do you think?
pygsp/operators/reduction.py
Outdated
| @@ -295,12 +300,22 @@ def kron_reduction(G, ind): | |||
| Graph structure or weight matrix | |||
| ind : list | |||
| indices of the nodes to keep | |||
| threshold: float | |||
| Threshold applied to the reduced Laplacian matrix to remove numerical | |||
| noise. (default: marchine precision) | |||
pygsp/operators/reduction.py
Outdated
| Lnew.eliminate_zeros() | ||
|
|
||
| # Enforces symmetric Laplacian | ||
| Lnew = (Lnew + Lnew.T) / 2. |
There was a problem hiding this comment.
Do we always want the Laplacian to be symmetric here? In the previous implementation, this line was called only under the conditional statement.
|
Thanks for the contribution. :) Could you look at @rodrigo-pena's comments so that we can merge it ? |
mdeff
force-pushed
the
master
branch
3 times, most recently
from
8321f45
to
4af195d
Compare
The motivation for these modification is to improve the scalability of the code to be able to handle much larger graphs.
With the baseline code, I couldn't compute graph pyramids for graphs with about 40,000 nodes, the peak memory consumption went above 500GB ! Turns out this was due to 2 main problems:
In the
kron_reductionfunction, the computation of the Schur complement was actually turning back the sparse input matrices into dense matrices. Plus, the computation of the Schur complement can be sped up by using a linear equation solver optimized to handle SDD matrices (symmetric diagonally dominant), which is the case of graph laplacians. There are even more efficient algorithms out there (starting with Cholesky decomposition) such as Koutis et al. (2011) arXiv:1102.4842 but I didn't find an out of the box python implementation and SuperLU seemed to do the job for me. I ended up rewriting bits ofkron_reductionto ensure that the sparse matrices were not turned into dense matrices when not necessary, and I addedpygsp.utils.splu_inv_dotas an optimized alternative to the use ofscipy.sparse.spsolve.In the
graph_sparsifyfunction, the computation of the resistance distances was done using blunt matrix inversion inpygsp.utils.resistance_distancewhich tries to invert it a potentially very large Laplacian matrix. After a kron reduction, this laplacian can have a lot of non zero elements, that's what was causing the worse of the memory consumption in my case. However, it turns out that the whole point of the Spielman-Srivastava spectral sparsification algorithm is to avoid having to compute this inverse, the algorithm only requires approximate distances and provides a procedure to compute them. I implementedpygsp.utils.approx_resistance_distanceto compute these distances based on the algorithm described in arxiv:0803.0929 . It still requires a way of solving a linear inverse system, where I used again my customizedpygsp.utils.splu_inv_dotbut it should be noted that there are much faster ways of doing this. Including the near linear time algorithm described in arXiv:1209.5821v3With these modifications, I can now compute graph multiresolutions in practice for graphs with about 50,000 nodes and 1e6 edges without running into memory issues on a standard desktop. The slowest point in the process is the sampling time for the edges in
graph_sparsify,scipy.stats.rv_discretetakes about 1s to sample 10,000 deviates, in my case the algorithm needed 16e6 which takes about 30min just to sample random numbers whereas the rest of the algorithm only takes minutes, so it's really stupid, but at least it eventually works.I tried to comment and document these modifications in the hope that they might be useful to others, everything seems to work for me but someone should take a close look to check that it doesn't break anything.