I’m pleased to announce Ripserer.jl, a pure Julia library for computing persistent homology. It uses the same algorithm as Ripser, which is known to be one of the fastest pieces of software for this purpose.
The list of supported features currently includes the following:
Vietoris-Rips persistent homology.
Sublevel set persistent homology for multidimensional image and time series data.
Calculation of persistent homology with coefficients in any (possibly user defined) field
with the default of intgers mod p for a prime p.
Sparse distance matrix and thresholding support.
Computing and plotting representative cocycles of persistent cohomology.
For a quick demo, please checkout the Examples section in the docs. If you’re unfamiliar with persistent homology, I recommend starting with this blogpost.
The main benefit (in my opinion) of this project is its generic API, which allows you to easily add new features to Ripserer without having to maintain your own fork of the project. For example, cubical persistent homology is implemented in 190 lines of code, which include docstrings.
The packages are registered, but still in development, so keep that in mind if APIs change unexpectedly. All ideas, feedback, issues and pull requests are more than welcome!
Wah, very cool. Thank you for linking that blog post as well, I needed it. But yea, that piqued my interest. It would be worth adding more examples of where this can be used and implemented. Looking forward to trying it out on images and trajectories. Nice!
Nice. How does it compare with existing packages in Julia for persistence homology? I remember at least two packages implementing similar functionality.
Thanks! Sadly, I haven’t had the chance to use this stuff in the real world yet, but I’m planning on adding a few more realistic examples. Pull requests or issues with examples or ideas would be very welcome!
Yeah, there are quite a few projects out there. Here’s a quick overview of what I know about:
Ripser.jl: my wrapper of the original C++ program. Very bare-bones and a bit outdated, but does essentially the same thing. Being a C++ wrapper, things like cubical homology are out of the question. According to my benchmarks, it’s around 30% faster. I’m still working on getting that number lower and from what I’ve seen Ripserer will actually outperform it once Julia 1.5 comes out. I haven’t tested enough to be sure.
Sparips.jl: this is a preprocessor that comes with a different wrapper of Ripser. The preprocessor allows you to compute persistent homology of very large datasets. Stealing its code and adding it to Ripserer is on my TODO list.
ComputationalHomology.jl: haven’t used it yet, but from the looks of it, it does quite a few things besides persistent homology. From what I’ve heard, the persistent homology algorithm is orders of magnitude slower than Ripserer.
PersistentCohomology.jl: does essentially the same thing, but with a different algorithm, that is orders of magnitude slower.
Eirene.jl: Not really a package, but more like an interactive tool for computing persistent homology. I tried wrapping my head around the algorithm it uses and failed. Has heaps of features. From what I heard its actually pretty fast. I still haven’t taken the time to learn to use it.
Calculation of persistent homology with coefficients in any (possibly user defined) field with the default of intgers mod p for a prime p .
Do you also have plans to support computing over the integers? I know this is terrible computationally, but from a math point of this would be quite interesting.
I also noticed performance was not great with PersistentCohomology but could not quite figure out why. Is the difference at the algorithmic level? PersistentCohomology uses the persistent cocycle algorithm from the de Silva et al. paper, are there algorithms that are much more efficient?
Sadly, that would require a completely different algorithm and is out of the scope of this project. I think I read somewhere ComputationalHomology.jl can do that, but I haven’t tried it.
I really recommend reading Uli’s paper on Ripser or watching his talk. I think he explains it better than I could, but essentially, it comes down to three things:
The twist optimization. This allows you to skip the reduction of some columns by ordering simplices by dimension.
Emergent pairs. This skips computing some zero length persistence intervals, like the ones that appear on each triangle.
The simplicial complex and some parts of the reduction matrix are not stored in memory, but rather recomputed on the fly. It turns out that reducing memory accesses makes a huge difference in speed. The CPU cache has nightmares about persistent homology computations.
Feel free to hit me up if you have any other specific questions on the algorithm.
By the way, have you looked at computing circular coordinates from persistent cohomology? That would be a really cool thing to have, but I haven’t looked into it enough to be sure it’s feasible.
No, I’ve never looked into that. I was mainly playing with implementations to see if it is possible to get an algorithm that is in some sense “functorial”, in that if there is a discrete group action (say a mirror symmetry) that preserves the filtration, then it is possible to see the effect of the group action on the cornerpoints. I haven’t done that yet (and am not actively working on it at the moment), but I’ll be curious to know if you end up having that feature at some point.
I’m not sure I understand this correctly. Do you want to map a function (mapping a simplex to a different simplex) on the simplices of the filtration and observe what happens to the results?
I think something like that could be done by defining a new simplex type that applies the function in the constructor. In theory, things should just work after that.
Also, my package Grassmann.jl supports homology and cohomology computations, although it’s not specifically designed for persistent homology, and also defines an interface for Simplex, SubManifold and Chain types (along with the boundary and coboundary operations), which are the foundation of simplicial homology.
The idea is that if you have a finite group G of symmetries acting on a filtration of topological spaces (or on a filtration of simplicial complexes), then the homology groups are not in \mathbf{Vec}_\mathbb{K} (where \mathbb{K} is the base field), but really they are in \mathrm{Rep}_\mathbb{K}(G), the category of group representations of G, and cornerpoints in the persistence diagram are indecomposable representations of G, each drawn in a color that represents the isomorphism class of the indecomposable representation. Shameless self-plug: see the section “multicolored persistence” in this paper, figure 2 may help to get some intuition.
I think it’d be nice if eventually there would be an efficient way to compute these “multicolored persistence diagrams”. To come back to your question, in practice I guess it boils to something that needs as input a filtration of simplices with a group of permutation acting on them (so each simplex has an “orbit” of other simplices on which it is mapped by this group). If the group is for example \mathbb{Z}_2, then it is either mapped onto itself, or onto another simplex.
This is probably a very good intuition. I was using simple tuples for simplexes, but with a custom type (encoding for example the “orbit” that the simplex belongs to, or something like that) could actually be enough to recover the multicolored diagram.
Nice! I really should get into your package some time, looks like it can do some pretty cool things. Maybe it would be possible to somehow integrate with Ripserer… I wouldn’t however change the Simplex type I’m using in the main project because it’s very specialized for this specific algorithm.
@piever I skimmed through the paper and it that looks pretty cool! I’ll try to go through the paper some time soon. If you ever feel an urge to work on this kind of stuff again, let me know.
