Computing the Skyscraper Invariant

Topological Data Analysis seeks to compute topological invariants from geometric data. In a typical setting, one constructs a simplicial complex $X$ filtered by a map $f:X\to {ℝ}^d$. The sublevel sets $X_{\le \alpha }≔\{x\in X|f(x)\le \alpha \}$ assemble to a functor $X_{\le \bullet }:{ℝ}^d\to \text{SCpx}$ via inclusion. Composition with homology ${𝐇}_{\ast }(-;𝐤)$ over a fixed field $𝐤$ yields a functor $V:{ℝ}^d\to {\text{Vect}}_{𝐤}$ - a Multiparameter Persistence Module.

The most important use case in practice are density-scale bifiltrations:

The analog of a bar in the multiparameter case is the Interval Module.

The Skyscraper Invariant arises from the Harder-Narasimhan filtration [31] of persistence modules for a specific set of stability conditions. We will explain how one could arrive at considering HN filtrations even without prior knowledge of geometric invariant theory [54]. This theory was originally developed for persistence modules in [27] and [26].

Let $V$ be a multiparameter persistence module. A sensible idea to extend the barcode is the following: search for an element with a "maximal lifetime" (an idea developed further in [52]). By "lifetime", we mean for any element $v\in V_{\alpha }$, with $\alpha \in {ℝ}^d$, the set

of parameter values under which $v$ persists and "maximal" refers here to its volume. To assure that every lifetime has a finite volume, we will assume that every module is bounded, which can always be realised by restricting the module to a bounded subset of ${ℝ}^d$.

After finding the vector $v$ of maximal lifetime, the submodule $⟨v⟩$ of $V$ generated by $v$ is an interval module and $L_v$ is the support $\mathrm{supp}⟨v⟩$ of $⟨v⟩$. Hence, if we compute the quotient $V/⟨v⟩$, and repeat this process, we obtain a filtration of $V$ in terms of interval modules.

The element of maximal lifetime is not compatible with direct sums of modules.

If we instead consider the minimal lifetime of an element, the construction becomes additive, meaning the desired vector lies in a direct summand:

This enables the use of decomposition algorithms (e.g. [20]) to speed up the search. Now, a priori, this element is not well defined, since in every non-zero persistence module we can find a sequence of elements whose lifetime approaches a set of empty volume. To fix this, we choose a concrete parameter value $\alpha \in {ℝ}^d$ and demand that the element be in $V_{\alpha }$.

Since we are only considering elements at $\alpha$, we can now restrict our computations to $⟨V_{\alpha }⟩$ – and this module is simpler than $V$:

In particular, $⟨V_{\alpha }⟩$ is more decomposable than $V$. Still, the choice of an element is not unique: there could be multiple vectors whose lifetime has the same (minimal) volume. This is the first hint that one should not only consider single elements in $V_{\alpha }$, but instead consider submodules of arbitrary dimension generated at a parameter value $\alpha$.

Let $U\subset V_{\alpha }$ be any sub vector space. The integral of the Hilbert function $\underset{¯}{dim}⟨U⟩$ over ${ℝ}^d$ would be the natural generalisation of the volume of the lifetime. However, using this definition would never make a subvector space of dimension higher than $1$ be the candidate for minimal lifetime. We need to average over the dimension of $U$.

The maximal submodule of highest slope $⟨U_1⟩\subset ⟨V_{\alpha }⟩$ is semistable. We define $⟨U_2⟩$ as the maximal submodule of highest slope in $⟨V_{\alpha }⟩/⟨U_1⟩$ and continue iteratively. This produces a sequence of semistable modules ${(⟨U_i⟩)}_{i\in [\mathrm{\ell }]}$. By defining inductively $F^i≔U_i+F^{i-1}$, we arrive at the Harder-Narasimhan-filtration [31]

of $⟨V_{\alpha }⟩$ by [26, Theorem (A)]. The semistable modules $⟨U_i⟩$ are by definition isomorphic to the factors $⟨F^i⟩/⟨F^{i-1}⟩$. The sequence $\mu (⟨U^i⟩)$ of slopes is decreasing in $i$, but this is not obvious. For a formal treatment, see [26] and [28, §2.4].

Because of the existence of modules like Section 1, the factors of the Harder-Narasimhan filtration at $\alpha$ should generally not all be intervals. Surprisingly, they did turn out to be intervals for all of our examples of density-scale bifiltrations of point sets in ${ℝ}^2$.

Then, in particular, all factors of the Harder-Narasimhan filtration are intervals by additivity.

Let $0\hookrightarrow F^1\hookrightarrow \mathrm{\dots }\hookrightarrow F^{\mathrm{\ell }}=V_{\alpha }$ be the HN filtration at $\alpha$.

This value only depends on the Hilbert function of the factors $\underset{¯}{dim}⟨U_j⟩$. Since $s_V^0(\alpha ,\beta )={dim}_{𝐤}{⟨V_{\alpha }⟩}_{\beta }=\mathrm{rank}(V_{\alpha \to \beta })$, the Skyscraper Invariant filters the rank invariant.

For each $\theta \in ℝ$, $s_V^{\theta }(\bullet ,\bullet )$ is a functor from ${({ℝ}^d)}^{\text{op}}\times {ℝ}^d$ to ${(\mathbb{Z}\cup \left\{+\mathrm{\infty }\right\})}^{\text{op}}$. This allows for an erosion distance $d_{\text{E}}$ [56, 55] under which the Skyscraper Invariant is stable [26, Theorem 3.7] with respect to the interleaving distance [14].

Given $s_V^{\theta }(\alpha ,\beta )$ for all $\alpha ,\beta$ in a grid $𝒢$, we can define the functor ${\widehat{s}}_V^{\theta }(\bullet ,\bullet )$ over ${({ℝ}^d)}^{\text{op}}\times {ℝ}^d$ via a projection to the next lower element in $𝒢$.

After computing the HN filtrations, we need to store the Hilbert functions of the filtration factors and their associated slopes. The Hilbert functions of each factor $⟨U_i⟩$ can be stored as a list of ${dim}_{𝐤}{U}_i$ intervals, ${\{W_{i,j}\}}_{j\in J}$, each being assigned the slope of $⟨U_i⟩$. For each $\alpha \in 𝒢$, this results in a list of at most $𝔱(V)$ intervals. To compute $s_V^{\theta }(\alpha ,\beta )$ from this list, we count how many $(W_{i,j})$ with an associated slope greater than $\theta$ contain $\beta$. The latter can be done in $𝒪(\log (m))$ time where $m$ is the number of relations of the interval. Since $m\le n$, This puts the query time in $𝒪(𝔱(V)\log (n))$.

The data structure in which persistence modules come is that of a presentation. This is a set of generators and a set of relations among the generators. They are collected in a matrix where each row corresponds to a generator and each column to a relation. The rows and columns are endowed with the parameter value – or degree – of the corresponding generator or relation, making this a graded matrix or monomial matrix [49]. If we write $A^G$ for the free module generated by a tuple of degrees $G$, then a graded matrix is nothing but a linear map $M:A^R\to A^G$ and it presents $V$ iff $V\stackrel{\sim }{\to }{A}^G/\mathrm{Im}(M)$.

Each finitely encoded persistence modules [51] has a finite minimal presentation which can be extended to a minimal resolution. The degrees of the generators, relations, higher relations etc. in a minimal resolution are independent of the chosen resolution and form multi-subsets $b_0(V),b_1(V),\mathrm{\dots },b_d(V)\subset {ℝ}^d$ – the graded Betti numbers of $V$. We denote the category of bounded finitely presented persistence modules by ${\mathrm{Pers}}_{\text{fp}}^{\text{b}}({ℝ}^d)$.

Computing a presentation for any induced submodule can be done by computing the kernel of a graded matrix in the category of persistence modules. For $d=2$, this can be done with the LWKR-algorithm (Lesnick, Wright, Kerber, and Rolle [43, 40]) and for $d=3$ there is an unpublished algorithm [5] by Bender, Gäfvert, and Lesnick. When $d$ is larger, we rely on classical machinery, like Schreyer’s algorithm [23], which is too slow in practice. We will therefore focus on the $2$-parameter case for this paper.

See [28, Proposition 3.1] for a chase in the following diagram:

$◀$

We condense the ideas from Section 1 into an algorithm.

Denote by $n=|G|+|R|$ the input size, $𝔱(V)$ the thickness of $V$, $k$ the maximal thickness of a uniquely generated indecomposable submodule, and by $T_{\text{hn-f}}$, $T_{\text{ker}}$, and $T_{\text{dec}}$ the times to compute the HN filtration, kernels, and decomposition for $d$-parameter modules.

Correctness follows from additivity (Equation 1) and the fact that we only need the ordered list of Hilbert functions of factors. At each step there are at most $𝔱(V)/k$ indecomposable submodules of thickness $k$ generated at $\alpha$. $◀$

Consider a minimal free resolution $0\to F_d\to \mathrm{\dots }\to F_0\to V\to 0$ which will be of maximal length $d$ by Hilbert’s Syzygy theorem. In our bounded case, it will even be exactly of length $d$ by duality of free and injective resolutions (first proved by Miller in [50] and rediscovered in [3] for the computation of Multiparameter Persistent Cohomology). Each $F_i$ has a basis with degrees at $b_i(V)$. By exactness, for each $\alpha \in {ℝ}^d$ we have

Let $B\in {ℝ}^d$ be large enough that $\mathrm{supp}(V)\subset {\prod }_{j\in [d]}(-\mathrm{\infty },B_j]$, then

As a function of $B$, Equation 5 is a polynomial and constant, so we can set $B=0$. $◀$

Instead of considering $⟨V_{\alpha }⟩$, we assume without loss of generality, that $V$ is uniquely generated at $0$. Recall that $𝔱(V)=\max \underset{¯}{dim}V$, and choose any basis ${𝐤}^{𝔱(V)}\stackrel{\sim }{\to }{V}_0$. We want to find the sub vector space $U\subset {𝐤}^{𝔱(V)}$ for which $⟨U⟩\subset V$ has the highest slope. An exhaustive search would involve computing $\mu (⟨U⟩)$ with Algorithm 1 for all sub vector spaces $U\subset V_0$. This is clearly possible if we work over a finite field, but also in the infinite case, the finite presentation of $V$ allows only a finite number of non-isomorphic submodules. Since all software uses finite fields, we fix a prime power $q$ and $𝐤={𝔽}_q$.

Denote by $𝒫V_0$ the set of 1-dimensional subspaces of $V_0$ and define $W_{\text{max}}\in \underset{W\in 𝒫V_0}{\arg \max }\mu (⟨W⟩)$

Consider any $U\subset V_0$ of dimension $k$. By assumption, there is a $1$-dimensional subspace $T\in 𝒫U\subset 𝒫V_0$ which is not in $V_k$. Therefore, .

$◀$

Every $k$-dimensional vector space over ${𝔽}_q$ contains exactly $(q^k-1)/(q-1)$ subspaces of dimension $1$, so the premise of Section 3 is trivially satisfied if . $V_0$ has $\sim q^{𝔱{(V)}^2/4}$ subspaces, but only $\sim q^{𝔱(V)}$ of these are $1$-dimensional.

Correctness follows from Section 3 and Section 3. Algorithm 1 computes a kernel for each subspaces of ${𝔽}_q^{𝔱(V)}$ , of which there are ${\sum }_i{\left(\genfrac{}{}{0pt}{}{𝔱(V)}{i}\right)}_q\sim q^{{(𝔱(V)/2)}^2+𝒪(𝔱(V))}$. $◀$

By calling Algorithm 3 recursively on the quotient by the submodule of highest slope we can compute the whole HN filtration at $\alpha$ in $𝒪(𝔱(V)\cdot T_{\text{ker}}\cdot q^{𝔱{(V)}^2/4+𝒪(𝔱(V))})$ time.

Fix two integers $N,N^{\prime }>0$ and a field $𝐤$. The space of $N\times N^{\prime }$ matrices with coefficients in $𝐤$ is denoted ${𝐤}^{N\times N^{\prime }}$ and we write $\otimes$ for the Kronecker product.

For the purpose of this paper, we present a simple random algorithm introduced in [29, Theorem 6.7] for $𝐤=\mathbb{Q}$. We extend this algorithm to any field using [34, Lemma 5.3]

Fix $W\in {\mathrm{Pers}}_{\text{fp}}^{\text{b}}({ℝ}^d)$ and a number $\epsilon >0$. Assume that the support of $W$ is contained in a bounded cube $B\subset {ℝ}^d$ and let $𝒢≔(\epsilon {\mathbb{Z}}^d)\cap B$. We denote by $V:𝒢\to \text{Vect}$ the restriction of $W$ to $𝒢$. Using the Persistence-Algebra library one can compute the structure maps of $V$ in time $n^3{(1/\epsilon )}^d+n^2{(1/\epsilon )}^{2d}$.

Let $\alpha \in 𝒢$, and let $\{\alpha ,{\beta }_1,\mathrm{\dots },{\beta }_k\}≔\{\beta \in 𝒢\mid \beta \ge \alpha \}$. Define $p_0≔dim{V}_{\alpha }$ and $q_0≔{\sum }_{i=1}^{k}dim{V}_{{\beta }_i}$. Consider the following subspace of ${𝐤}^{q_0\times p_0}$:

Cheng also observed in the above Theorem that when $(p,q)≔(p_0,q_0)$, either $\mathrm{\ell }=1$ or [15, Corollary 3.4]. They thus devised Algorithm 4.

By [26] (as detailed in [28, Corollary 4.5]), applying for every $\alpha \in 𝒢$ Algorithm 4 to the restriction of $V$ to $𝒢$ yields an $\epsilon$-approximation of $s_V$. By Section 4 with $\mathrm{\ell }≔k{p}_0{q}_0$ and $N≔p_0{q}_0$, line 3 in Algorithm 4 takes $𝒪(|𝒢|{(p_0{q}_0)}^{9+\eta /2})=𝒪(𝔱{(V)}^{18+\eta }{\epsilon }^{-10d-\eta })$ for any $\eta >0$. The recursive step adds a factor $𝔱(V)$, whence the desired complexity. $◀$

By harnessing the block structure of (6), we are able to improve the complexity in Section 4 for the subspaces appearing in Algorithm 4.

Further empirical improvements can be obtained by computing $i(p,q)$ in Section 4 for small integers $p,q>0$ such that $\frac{p}{q}\approx \frac{p_0}{q_0}$. Indeed, the complexity of computing $F^{i(p,q)}$ is heavily dependent on the size of $pq$. For simplicity, assume that $\frac{p}{q}>\frac{p_0}{q_0}$. If $F^{i(p,q)}\ne 0$, then we can replace $U$ by $F^{i(p,q)}$ in Algorithm 4. Otherwise, we know that there are no submodule $U\subset V$ such that $\mu (U)>\frac{p}{p+q}$. We keep choosing better approximations $\frac{p}{q}$ of $\frac{p_0}{q_0}$ (so with higher values of $p,q$) until either we find a proper submodule in the HN filtration of $V$ at $\alpha$, or we can guarantee that $V$ is semistable at $\alpha$.

Finally, given $p,q>0$ and a matrix $A\in {𝐤}^{p\times q}\otimes {𝒜}_{\alpha }$, our implementation uses the block structure of $A$ to partially reduce it: Since $A={({\lambda }_{ij}{V}_{\alpha \to {\beta }_{i\text{ modulo }k}})}_{ij}$ for some $\mathrm{\Lambda }\in {𝐤}^{pk\times q}$, we column-reduce the first $q$ block rows of $A$ by reducing $\mathrm{\Lambda }$ and each ${(V_{\alpha \to {\beta }_i})}_{1\le i\le k}$.

These two empirical improvements led to significant speed-ups (see [28, Table 3]).

For many computations on persistence modules, it is often enough to perform them only at those parameter values that share each of their coordinates with a generator or relation (cf. the barcode template [43, 44]). We call the set of these parameters, $\text{Grid}(V)\subset {ℝ}^d$, the induced grid of $V$. It partitions ${ℝ}^d$ into rectangular cells.

Consider any cell $C$ in the grid on the left in Figure 8, which is slightly coarser than the induced grid. For all $\alpha \in C$, the submodules $⟨V_{\alpha }⟩$ will look qualitatively similar: Each has a presentation with the same underlying matrix. We will use this fact to avoid recomputing decompositions and HN filtrations for every $\alpha \in C$ on an $\epsilon$-grid.

For $J\subset [d]$, denote by $F_J\subset {ℝ}^d$ the hyperplane spanned by the coordinate axes indicated by $J$, with the origin at $\alpha$. The graded Betti numbers, which lie in some $F_J$ for , will be shifted to align with $\alpha \in C$ when $\alpha$ changes (Figure 9).

Full proof in [28, Lemma 6.5]

$◀$

We describe this computation for $d=2$, where after ignoring the degree-$2$ part, the slope-polynomials are linear and so the minimisation diagram is convex-polygonal.

The thickness-$6$ module from Table 1 produces the following subdivision.

We ran Algorithm 5 for all modules used in Table 1. In each case, as in Figure 11, the first grid cell lies entirely in one face of the subdivision.

Correctness follows from Section 5 and Section 5. The runtime is the same as for Algorithm 3, except that kernel computation is in $𝔱{(V)}^3$ for 2-parameter modules and we compute a lower envelope. That can be done in expected time $n\log n$, putting the runtime in

$◀$

For each $i\in I$, there is a highest-slope submodule $W^{\{i\}}\subset ⟨V_{\alpha }⟩$ associated with the face $S_i$. By computing the quotient $V^{\{i\}}≔⟨V_{\alpha }⟩/W^{\{i\}}$ and recursively calling Algorithm 5 on the pair $(V^{\{i\}},S_i)$, we obtain a nested sequence of subdivisions of the grid cell $C$ and an associated tree of filtration factors. Glueing these for all grid cells, we obtain a data structure from which we can compute the exact value of $s_V^{\theta }(\alpha ,\beta )$ by locating $\alpha$ in the nested subdivisions [28, Section 6.4].

For all $\beta$ in a region of this nested sequence of subdivisions, the module $V$ has equivalent HN filtrations. This is one implication of Theorem 1.

In practice, where one wants to compute filtered vectorisations of the Skyscraper Invariant, it is not a problem to restrict the module to an $\epsilon$-grid $𝒢$ and use Algorithm 2. The real power of Theorem 1 and Algorithm 5 lies in the fact that it is a recipe to avoid recomputations: For any direct summand $V_i$ of the input module, the induced grid $\text{Grid}(V_i)$ is typically coarser than $𝒢$ in some places and finer than $𝒢$ in others. We therefore overlay the two grids to form the restriction or parallel grid and use Algorithm 5 only for those cells in each $\text{Grid}(V_i)$ that intersect $𝒢$. This technique enables fast computation of the Skyscraper Invariant at high resolution.

To generate indecomposable uniquely generated submodules of non-zero thickness we used synthetically generated noisy samples of a torus embedded in ${ℝ}^3$ (cf Section 1) and random presentation matrices. We test the exhaustive search Algorithm 3 and its filter Section 3. For thickness $9$, Algorithm 3 always timed out.

We conclude that the filter (Section 3) for Algorithm 3 is barely helpful.

After restricting the torus generated modules to equidistant grids, we have used Cheng’s algorithm (Algorithm 4) over ${𝔽}_2$ and $\mathbb{Q}$. Since we can only generate modules from bifiltrations over ${𝔽}_2$, we interpreted them as $\mathbb{Q}$-modules when possible.

We observe a great variability of runtimes for Cheng’s algorithm depending on whether the heuristic introduced at the end of Section 4 successfully computes and certifies the correctness of all the HN filtrations. The random algorithm from [30] performs computations on square matrices with $𝔱{(V)}^2{\epsilon }^{-d}r$ rows for a large enough $r>0$. For $𝐤=\mathbb{Q}$ (but not for $𝐤={𝔽}_2$), we have observed in our experiments that $r=1$ is large enough.

Table 2 shows the strong dependence of Algorithm 4 on the size of the grid to which the uniquely graded module has to be restricted. Algorithm 3 does not require this restriction but, as visible in Table 1, it is highly dependent on the thickness of the input.

The asymptotic runtimes of the algorithms suggest that for practical grids $(>50\times 50)$, Cheng’s algorithm would outperform the exhaustive search not before $𝔱(V)\ge 20$. At that point, the number of arithmetic operations is too large for any computer to handle.

We used benchmark data from [39, 61] to compare the parallel grid scan algorithm with Algorithm 2.

The parallel grid scan should not show the same quadratic dependency on the grid size of Algorithm 2. Table 4 confirms this and shows that our algorithm computes the Skyscraper Invariant in a practical timeframe for relevant data sizes.

The Multiparameter Landscape [60] is directly defined by the ranks of the diagonal maps $V_{\alpha }\to V_{\alpha +(t,t)}$. Therefore, we can filter the Landscape by replacing the rank with the Skyscraper Invariant. As an example, we will use the IHC stained sample LargeHypoxicRegion1 of a tumour, for which the authors of [61] analyse the location of immune cells using Multiparameter Landscapes.

Using multipers [46], we computed its density-alpha bifiltration with function-delaunay [63], a minimal presentation of the first persistent homology group thereof with mpfree [40], and restricted it to a $300\times 300$ grid. Computing the Skyscraper Invariant on this $35360\times 25553$ matrix took around 5 minutes.

Since the filtered Landscapes are stable [26, Corollary 3.12], their differences are too, revealing elements of long lifetime.

With our algorithms and the induced filtered Landscapes practitioners have a new invariant for topological data analysis that is robust to noise, finer than the rank invariant, efficiently computable, and is interpretable directly from its visual representation.

To incorporate invariants into Machine Learning pipelines, the authors of [58] define a framework which needs to be extended to include the Skyscraper Invariant. Instead of Landscapes, one can also filter other vectorizations of the rank invariant. By changing the stability condition [26], we obtain different HN filtrations and invariants. Can our methods be adapted to this general setting?

For each $\theta \in ℝ$, the Skyscraper invariant can be Möbius inverted [57], just like the rank invariant [48]. What are the properties of this path?