Finding a loop used the Vietoris-Rips complex because it is the default. TDA4j offers several other ways to turn a point cloud into a filtered complex, and they trade off size, exactness and what they need to know about your data. This tutorial runs the same data, data/noisy-circle.csv, through five of them and compares what comes out.

Persistence takes the complex as an option for the three that need only the points, and any complex you build yourself as its input:

import scala.language.experimental.modularity
import org.appliedtopology.tda4j.*

val points = CSV.readPointCloud("_docs/tutorials/data/noisy-circle.csv")
val metricSpace = EuclideanMetricSpace(points)
val landmarks = LandmarkSelector.maxmin(metricSpace, 15).landmarks

val vr = Persistence(points)                                  // Vietoris-Rips, the default
val cech = Persistence(points, complex = Cech)
val alpha = Persistence(points, complex = AlphaShapes)
val sparse = Persistence(SparseRips(metricSpace, epsilon = 0.5, maxDimension = 1))
val witness = Persistence(Witness(metricSpace, landmarks, maxDimension = 1))

val diagrams = List(vr, cech, alpha, sparse, witness)
val loops = diagrams.map(_.dim(1).longest.get)   // the longest-lived loop of each: the table below
diagrams.map(_.bettiNumbers(0))                  // one connected component in each

What they cost is the size of the complex, the number of simplices the engine has to process:

VietorisRips(metricSpace, maxDimension = 1).iterator.size   // 24711

The five complexes

complex simplices the loop (birth, death)
Vietoris-Rips 24,711 (0.595, 1.707)
Čech 36,050 (0.297, 0.935)
alpha 325 (0.297, 0.935)
sparse Rips (ε = 0.5) 1,730 (0.708, 1.789)
witness (15 landmarks) 441 (0.278, 0.742)

All five find the loop, and all five find exactly one connected component. What differs is the cost and the units.

Vietoris-Rips: the default

The Vietoris-Rips complex joins any group of points whose pairwise distances are all at most r. It needs only distances, not coordinates, so it works on any distance matrix, and it has the fastest engines (the Ripser-style engine behind the command line and MATLAB defaults is specific to it). Its drawback is size: 60 points already give 24,711 simplices, and the count grows quickly with the number of points and the dimension you ask for.

Čech: the geometrically exact one

The Čech complex at scale r joins any group of points whose balls of radius r have a common point. It is the complex the "nerve theorem" speaks about directly: it has the same shape as the union of the balls. It also needs coordinates (it finds the smallest ball around each group), and it is the largest of the five (36,050 simplices here).

Notice the units. Čech scale is a radius, Vietoris-Rips scale is a diameter, so the two loops are born at scales that differ by a factor of two (0.595 against 0.297) and die at different places, √3 = 1.73 for Rips against 1 for Čech. Bars from different complexes must not be compared by their endpoints, only by their shape: here, one clear loop, far longer-lived than anything else.

Alpha: Čech without the bulk

The alpha complex is a subcomplex of the Delaunay triangulation of the points, chosen so that at each scale it has the same shape as the union of balls, which is to say the same shape as the Čech complex. So the two give identical bars, the same loop to six digits, from 325 simplices instead of 36,050. When you can use it, that is a hundredfold saving for nothing.

You can use it when your points have coordinates in a low-dimensional space (the Delaunay triangulation is what limits it: fine in two or three dimensions, rapidly harder beyond). The fast-alpha engine, described in the user guide, goes further by skipping the generic reduction.

Sparse Rips: Rips for many points

Sparse Rips (Cavanna, Jahanseir and Sheehy) keeps only a well-spread subset of the points at each scale, which makes the complex linear in the number of points. The price is approximation: bars are guaranteed correct only up to a factor 1 + ε. Here the loop is born at 0.708 instead of 0.595 and dies at 1.789 instead of 1.707, so it is clearly the same loop, with the endpoints shifted by the approximation. Smaller ε is closer to plain Rips and larger complexes. Use it when plain Vietoris-Rips is too big.

Witness: a few landmarks, all the points

The witness complex first picks a small set of landmark points (here 15, chosen by the furthest-point rule so that they cover the cloud evenly) and builds the complex on them only, using the other points as witnesses that decide which simplices exist. That makes it tiny (441 simplices) and very cheap, at the cost of thinking about landmarks. The loop is still there, at its own scale. See the user guide for choosing landmarks and the two witness variants.

Which one?

  • Start with Vietoris-Rips. It needs the least and is the fastest to compute with.
  • Low-dimensional coordinates, want the exact geometry: alpha. It is the same answer as Čech at a hundredth of the size.
  • Too many points for Vietoris-Rips: sparse Rips if you want to keep all the points and accept an approximation, witness if you are happy to work with landmarks.
  • Noisy data, or outliers: none of these is robust; see the next tutorial.

The whole script

import scala.language.experimental.modularity
import org.appliedtopology.tda4j.*

val points = CSV.readPointCloud("_docs/tutorials/data/noisy-circle.csv")
val metricSpace = EuclideanMetricSpace(points)
val landmarks = LandmarkSelector.maxmin(metricSpace, 15).landmarks

val vr = Persistence(points)
val cech = Persistence(points, complex = Cech)
val alpha = Persistence(points, complex = AlphaShapes)
val sparse = Persistence(SparseRips(metricSpace, epsilon = 0.5, maxDimension = 1))
val witness = Persistence(Witness(metricSpace, landmarks, maxDimension = 1))
val diagrams = List(vr, cech, alpha, sparse, witness)
val loops = diagrams.map(_.dim(1).longest.get)

val sizes = List(
  VietorisRips(metricSpace, maxDimension = 1),
  Cech(metricSpace, maxDimension = 1),
  AlphaShapes(points),
  SparseRips(metricSpace, epsilon = 0.5, maxDimension = 1),
  Witness(metricSpace, landmarks, maxDimension = 1)
).map(_.iterator.size)
javaaddpath('target/out/jvm/scala-3.9.0/tda4j/tda4j-0.5.0-SNAPSHOT-assembly.jar');   % from the repository root, after sbt assembly
import org.appliedtopology.tda4j.matlab.*;

points = readmatrix('_docs/tutorials/data/noisy-circle.csv');

% Each row: the value of the 'complex' option and any extra options it needs
settings = {
    'vr',          {}
    'cech',        {}
    'alpha',       {}
    'sparse-rips', {'sparseEpsilon', '0.5'}
    'witness',     {'numLandmarks', '15'}};

for s = 1:size(settings, 1)
    options = [{'maxDimension', '1', 'complex', settings{s,1}}, settings{s,2}];
    result = TDA4j.computeFromPoints(points, options);
    bars = result.toArray();
    loops = bars(bars(:,1) == 1, :);
    [~, longest] = max(loops(:,3) - loops(:,2));
    fprintf('%-12s %6d simplices, loop [%.3f, %.3f)\n', settings{s,1}, result.numCells(), loops(longest,2), loops(longest,3));
end

In MATLAB a complex is a value of the complex option rather than a function, and numCells() is the number of simplices the engine had to process. MatlabTabsSpec makes these same calls from Scala and checks every number quoted in the MATLAB tab; MATLAB itself is not run by the test suite.