Persistent Topology for Peak Detection
This is a little story of how an abstract mathematical theory like persistent homology can be a practical tool for a usual task in industry.
Update: A related paper on this topic appeared at iDSC2020.
Motivation
At B+R Industrial Automation we have an open discussion culture and from time to time there are little take-away questions for people from computer science or mathematics. One such task was the following:
In order to identify the natural frequency of some mechanics attached to an ACOPOS servo drive one actuates the drive with a frequency sweep signal and records the response signal of the position. The (Fourier) spectrum of such a response signal could then look like in the picture below. The question is now how to detect the “main peaks” in this signal?
The original problem statement came from motion control, however, the underlying task carries more universality1: Given a real-valued signal resp. a time series, where are the “significant peaks”?
One of the first attempts would be to find maxima. The global maximum is not what we look for because we might be interested in this very significant local maximum at 10. However, all (or only some?) local maxima is probably also not what we look for: (i) the local maximum at 33 is quite insignificant compared to the local maximum at around 55 although the signal is higher and (ii) measured signals come with noise that cause many local maxima that are irrelevant for our analysis.
We could get rid of the noise by smoothing the signal. However, this often comes at costs: For instance the sharp peak at 10 would be quickly dampened and our impression on the “significance” of the peak would strongly altered. If possible we would always prefer to operate on the original data. But even if we would decide to smooth the signal, say, by a Gaussian smoothing kernel, what would be suitable standard deviation? Bringing bigger guns like scale space into the game appears to be inadequate to address the simple question originally raised.
Persistence
What we look for is a method to identify the “relative maxima” in a stable, noise-resilient fashion, and, moreover, quantify the significance (or dominance) of these maxima. If we look for the eight most significant peaks in the signal above then we would expect something similar to this:
What the plot above depicts are the eight most persistent 0-dimensional cycles in the sublevel set filtration of the function graph.
In simple words: Consider a water level starting above the global maximum that continuously lowers its level. Whenever the water level reaches a local maximum a new island is born. Whenever it reaches a local minimum two islands merge, and we consider the lower island to be merged into the higher island. This gives a precise notion of a lifespan of an island; it is the significance we mentioned above, or the so-called persistence in persistent homology. In this precise sense the above plot shows the eight most persistent islands/0-dimensional cycles on the original data.
The blue bars start at local maxima and grow downwards until the level where the corresponding island merges with another island. Only the island raised by the global maximum never dies; everything is eventually merged into it on the way to minus infinity. The blue bars are known as the barcode in persistent topology.
The persistence diagram is another representation of the features of persistent homology. It depicts every island as a point with birth- and death-level as the x- and y-coordinates. In our example the (0-dimensional) persistence diagram looks as follows:
An efficient algorithm
The above example is a simple case of the general problem of computing the persistent homology (over \(\mathbb{Z}_2\)) in a given dimension of a filtration of a simplicial complex. The classical approach to this is via matrix reduction (c.f. Gaussian elimination) by Edelsbrunner and Harer2 and takes \(O(n^3)\) time, where \(n\) denotes the length of the signal in the above example, or the number of simplices of the underlying simplcial complex in general. This is what my C++ library libstick or one of the several other available implementations do.
In the simple case of the sublevel set filtration of a real-valued signal, however, we can do much better, namely in the speed of sorting numbers, which is \(O(n \log n)\), and the algorithm is very simple, too:
class Peak:
def __init__(self, startidx):
self.born = self.left = self.right = startidx
self.died = None
def get_persistence(self, seq):
return float("inf") if self.died is None else seq[self.born] - seq[self.died]
def get_persistent_homology(seq):
peaks = []
# Maps indices to peaks
idxtopeak = [None for s in seq]
# Sequence indices sorted by values
indices = range(len(seq))
indices = sorted(indices, key = lambda i: seq[i], reverse=True)
# Process each sample in descending order
for idx in indices:
lftdone = (idx > 0 and idxtopeak[idx-1] is not None)
rgtdone = (idx < len(seq)-1 and idxtopeak[idx+1] is not None)
il = idxtopeak[idx-1] if lftdone else None
ir = idxtopeak[idx+1] if rgtdone else None
# New peak born
if not lftdone and not rgtdone:
peaks.append(Peak(idx))
idxtopeak[idx] = len(peaks)-1
# Directly merge to next peak left
if lftdone and not rgtdone:
peaks[il].right += 1
idxtopeak[idx] = il
# Directly merge to next peak right
if not lftdone and rgtdone:
peaks[ir].left -= 1
idxtopeak[idx] = ir
# Merge left and right peaks
if lftdone and rgtdone:
# Left was born earlier: merge right to left
if seq[peaks[il].born] > seq[peaks[ir].born]:
peaks[ir].died = idx
peaks[il].right = peaks[ir].right
idxtopeak[peaks[il].right] = idxtopeak[idx] = il
else:
peaks[il].died = idx
peaks[ir].left = peaks[il].left
idxtopeak[peaks[ir].left] = idxtopeak[idx] = ir
# This is optional convenience
return sorted(peaks, key=lambda p: p.get_persistence(seq), reverse=True)(I have received many mails about the license of the code. You may use the code under MIT license.)
The general theme for this algorithm is to track growing islands and merge them from time to time. A union-find data structure is what we look for in general. However, since our sets are disjoint intervals on the real line, we can replace the part after sorting by signals values to a loop with a simple, constant-time body instead of otherwise time-optimal \(O(\alpha(n))\) union-find operations.
Concluding remarks
Typically one considers the superlevel set filtration in which the water level rises and death levels are higher than birth levels, which then leads to the usual definition of persistence as death level minus birth level and dots in the persistence diagrams reside above the diagonal. However, the sublevel set filtration of one function is the superlevel set filtration of the negative function.
References and further reading
This article deals with algorithms and mathematical formalisms that belong to the field of computational topology. The introduction to this field is the book by Edelsbrunner and Harer2.
Examples for the 0-th and 1-st dimensional persistence diagram on two-dimensional input are given on the demonstration page of the libstick library. The getting started of libstick describes the general work flow from defining simplicial complexes to computing persistence diagrams.
Update. A demonstration and code of peak detection in two-dimensional data is given here.