On reading papers for what they don't say

A nice paper

I stumbled upon a review article by Jose Perea in the Notices of the AMS, Topological Time Series Analysis. It walks you through Takens’ embedding theorem, sliding-window point clouds, persistent homology, and then applies the whole apparatus to problems like detecting wheezes in breath recordings and scoring genes by their circadian rhythmicity. For someone who doesn’t know much about topology, the mathematics is clean and understandable, the applications are eclectic and honest. Overall a fun and interesting read.

What I couldn’t stop thinking about afterwards was Theorem 3.1, which is Takens’ embedding theorem dressed up for the sliding-window case. The theorem says, very roughly, that if you have a scalar measurement of a periodic dynamical system and you stack its time-shifted copies into vectors, the point cloud you get back is a topological copy of the underlying attractor. A circular orbit produces a loop. A Keplerian orbit at eccentricity 0.9 also produces a loop. They look different, but topologically they’re the same thing: one $\text{H}_1$ class, Betti number $\beta_1 = 1$. The shape of the loop depends on the waveform; its existence doesn’t.

Every period-finding method I use in my day job assumes something about the shape of the signal. Lomb-Scargle assumes a sinusoid. BLS assumes a boxcar transit. TLS adds a transit model with limb darkening and ingress/egress. While BLS and TLS aren’t really affected by waveform mismatches (transits tend to generally look the smae or very similar), GLS does suffer from a waveform mismatch from orbital eccentricity: an eccentric Keplerian RV is a genuine multi-harmonic Fourier series, with power spreading across P, P/2, P/3 and beyond as eccentricity rises, and the absolute amplitude of the fundamental collapsing. Cumming (2004) documents the sensitivity loss, which gets painful for e $\gtrsim$ 0.5. The field handles this with dedicated Keplerian periodograms. BLS and TLS work fine on eccentric transits; what those methods actually miss, per Hippke & Heller (2019), is small planets and shallow ingress/egress detail.

Here was a framework that on paper didn’t care about any of that. It would build a loop whether the waveform was a sine, a sawtooth, or an eccentric Keplerian spike. And when I went looking, I couldn’t find anyone who had tried it on astronomical data. The review cites applications to biology, to engineering, to video; none to astronomy. None to the specific problem of exoplanet period search.

So I built the thing.

Building it

What I built is a Python package, takens_periodogram, that does the natural thing. For each trial period P, construct the sliding-window delay embedding with lag $\tau = \frac{P}{m}$ in $m$ delay coordinates, compute the $\text{H}_1$ persistence of that cloud with Ripser, and return the max persistence as a function of trial period. Add the standard scaffolding: bootstrap false-alarm levels via phase-randomised surrogates, candidate detection, harmonic deduplication, NetCDF serialisation. Make it install with uv. Add injection-test notebooks for the canonical cases. The works.

It works on synthetic data in the sense that the code runs and the numbers are reproducible. Whether it finds periods turned out to be a different question.

The first wall

I ran a 16-case injection benchmark. Six circular-orbit periods (2, 3, 5, 15, 30, 80 days) at e = 0, plus P = 3 d at e = 0.5 and 0.7, each evaluated in both time-domain and phase-folded embedding. Clean white noise, a textbook Keplerian signal.

Detection rate: 2/16. 12.5 per cent.

And it wasn’t that the peaks were slightly off. They landed at completely wrong periods in 14 of 16 cases. Here’s one: a circular orbit at $P = 3$ d, the easiest regime you’d design a benchmark around. The algorithm places its top peak at about 7.86 d.

The landscape is noisy and multi-modal. The literature period sits in a local dip. This isn’t a signal-to-noise problem; the $\frac{K}{\sigma}$ ratio of the injection is plenty.

It gets worse with eccentricity. A 30-seed sweep at $P = 15$ d, $K = 10$ $\text{ms}^{-1}$, $N = 80$ observations shows $\pi(P)$ collapsing eight-fold from $e = 0$ to $e = 0.6$ in phase-folded mode, five-fold in time-domain. The first harmonic 2P overtakes the fundamental at $e \approx 0.3$ for time-domain, $e \approx 0.5$ for phase-folded:

The physical picture is clear. Eccentric orbits spend most of their time near apastron, so the periastron arc of the embedded loop becomes a sparsely sampled gap and the $\text{H}_1$ feature breaks. Folding at 2P doubles the periastron sampling, which is why 2P wins.

Over the next couple of months I tried eight variants. Time-domain and phase-folded embeddings, adaptive tolerance floors, adaptive embedding dimensions, a joint RV+photometry mode, peak refinement, FAL smoothing, a hybrid of time-domain and phase-folded. Each was motivated and each was honestly tested. None moved the detection rate above about 12 per cent.

Some of them helped with subsidiary problems, foe example, phase-folded embedding cleanly fixed a short-period NaN issue that time-domain suffered, but fixing that didn’t lift the detection rate. The real issue was that max $\text{H}_1$ persistence wasn’t a sharp enough period statistic. The contrast ratio between the true period and the largest spurious peak was $1.5-2\times$. That’s not going to cut it; you need more like $5\times$.

At some point I sat with the numbers for a while and accepted that the core periodogram was dead.

Rescue attempt: a diagnostic

I didn’t want to write the whole project off, so I looked at what had fallen out of the work that might survive the verdict on the periodogram itself. Two things had: a lifetime-ratio statistic (ℓ, comparing the top $\text{H}_1$ persistence to the second) and an epoch-split stability test I was calling S_norm. It was basically a discretised Wasserstein-1 distance on per-phase-bin means, comparing the first and second halves of the baseline. The intuition: a Keplerian signal is temporally stable (the phase-folded curve looks the same in both halves); stellar activity isn’t (spots grow, decay, migrate).

On 40 synthetic injection cases, S_norm looked really clean:

Signal class	S_norm median
Keplerian, $e \leq 0.5$	1.3
GP quasi-periodic	6.1
Evolving activity	4.3
White noise	0.9

Factor-2.5 separation between any planet-like and any activity-like signal, eccentricity-independent, 100 per cent flagging of activity with zero false positives on planets. If it held on real data it would have been a letter on its own. Planet/activity disambiguation is a genuinely live problem in precision RV surveys.

So I ran it on 19 HARPS systems where the literature verdict is known: seven confirmed planets (Pi Mensae, $\epsilon$ Indi A, GJ 887, $\nu^2$ Lupi, HD 114613, HD 192310, HD 38858), three retracted activity false positives (HD 41248, Kapteyn’s star, GJ 581), six controversial cases, three quiet calibrators.

Every one of the seven confirmed planets got flagged as “evolving” (activity) at the top GLS peak.

Pi Mensae: S = 229. $\epsilon$ Indi A: S = 25. HD 38858: S = 34. The median S for confirmed planets came in at 4.50. For retracted activity false positives it was 2.22. The diagnostic was inverted: confirmed planets looked more activity-like than known activity.

I spent an afternoon or two trying to rescue it. The formal RV uncertainties from the DRS pipeline are photon-limited, but real RV scatter on solar-type stars is floor-limited by granulation, oscillations, and instrumental systematics. So the S_norm denominator was too small by a factor of several, and the statistic blew up for any real signal regardless of its temporal stability. I swapped in an empirical-variance denominator (in-half sample variance instead of propagated formal error), per-instrument-per-half offset fitting to absorb the HARPS fibre-upgrade discontinuity, and Keplerian-subtracted residuals at the known planet period. The inversion went away:

class	median S_norm (rescued)
confirmed planets	1.49
activity FPs	1.16
controversial	1.88

But the factor-2.5 separation didn’t come back. 1.49 vs 1.16 is a 29 % gap, and the within-class spreads are much wider than that. Confirmed planets ranged from 1.40 to 11.17; activity FPs from 1.06 to 2.30. The synthetic-clean separation was an artefact of two things: formal errors being calibrated exactly to the injected noise, and the injected Keplerian being the only structure in the data. Real stars violate both.

That was the second burial.

A side quest

While I was working on the diagnostics, I got curious about a tangential question. Persistent homology is supposed to find topology in point clouds. Does it see the known structure in the exoplanet population itself. The Neptune Desert, the radius gap, the hot-Jupiter clump, that kind of thing.

I pulled the NASA Exoplanet Archive Composite catalog, subsampled down to 2000 planets via farthest-point sampling, and ran persistent homology up to $\text{H}_2$ on the 4D point cloud (log P, log M, log R, log insolation). Null: permute each coordinate independently across 20 surrogates, compare to the real distribution.

dim	data max persistence	null max persistence	p-value
H0	2.88	$2.48 \pm 0.42$	0.10
H1	0.33	$0.51 \pm 0.06$	1.00
H2	0.15	$0.28 \pm 0.03$	1.00

The population has comparable H0 persistence to the null and suppressed H1 and H2. I had to stop and think about why.

Persistent homology on an unweighted point cloud probes the topology of the support: where the distribution is connected, where it has holes, where it has enclosed voids. The Neptune Desert is not a hole in the support; it’s a density depression inside a continuous distribution. A Vietoris-Rips filtration on the unweighted cloud just doesn’t see that unless the depression is deep enough to disconnect the cloud, which it isn’t. So persistent homology doesn’t recover the known astrophysical structure, and not because the structure is absent. It’s because the structure isn’t topological in the sense the tool measures.

Clean result, but not one that rescued anything.

Going denser

The most sympathetic steelman of the periodogram failure was: “you tested on sparse RV. Of course it struggles. Try it on a signal where you have a lot of data per cycle.” Which is roughly what pulsating variable stars are. RR Lyrae, especially RRab pulsators, are non-sinusoidal (sawtooth in flux), well understood, and densely sampled by modern photometric surveys. If a waveform-agnostic periodogram was ever going to show its worth, it would be here.

I grabbed Kepler Q1 photometry of RR Lyrae itself: KIC 7198959, the prototype of the class. 1624 long-cadence points over 33 days, sampling 59 cycles of the literature period $P = 0.5668$ d. Ran the package across four configurations: phase-folded and time-domain, with embedding dimensions m = 3 and m = 5.

configuration	top peak (d)	offset	$\pi(P_\text{lit})$	$\pi(\text{runner-up})$	contrast
phase-folded, m = 3	0.998	76 %	0.43	1.45	0.29
phase-folded, m = 5	0.740	31 %	2.02	2.61	0.78
time-domain, m = 3	0.232	59 %	0.66	1.60	0.41
time-domain, m = 5	0.292	48 %	1.82	2.22	0.82

GLS on the same light curve hits the literature period at 0.5676 d, off by 0.13 %.

Every contrast ratio is below one. The literature period carries less $\text{H}_1$ persistence than the algorithm’s chosen peak in every configuration. Higher $m$ raises the persistence at the true period by roughly $4\times$, but raises spurious-peak persistence by a similar factor. The landscape gets more topologically interesting everywhere, not sharper where it should be.

So denser sampling didn’t rescue it either. If anything it did worse than the exoplanet tests: the contrast ratio of $0.29-0.82$ at RRab is below the $1.5-2\times$ of the exoplanet regime. Not the result I wanted.

Back to the math, with actual attention

At this point I was out of empirical rescue attempts and went back to Perea with a different question. Not “how do I tune this to work?” but “why does it not work?”

The paper’s Theorem 5.1 gives the semi-axes of the sliding-window ellipse for a pure sinusoid $\zeta(t) = \sin(Lt + \phi)$:

\[a^2,\ b^2 \;\propto\; (d+1) \pm \left|\frac{\sin(L(d+1)\tau)}{\sin(L\tau)}\right|\]

The ellipse becomes a perfect circle and $\text{mp}_1$ reaches its maximum exactly when the numerator $\sin(L(d+1)\tau)$ vanishes. That happens at $\tau = n\pi / (L(d+1)) = (n/2) \cdot P / (d+1)$ for every positive integer n. The paper’s “optimal $\tau$” ($\tau_d = 2\pi / (L(d+1))$) is just the $n = 2$ case.

My lag rule was $\tau = P/(d+1)$ (the natural choice: it makes the embedding window span one period). Plug that in and you get circular optima at

\[P \;=\; \tfrac{n}{2}\,P_\mathrm{true}\quad\text{for every } n = 1, 2, 3, \dots\]

P_true, P_true/2, 3P_true/2, 2P_true, 5P_true/2: every half-integer multiple. The optima are dense. Add harmonics (eccentric orbits, sawtooth RRab waveforms) and each harmonic contributes its own family. The comb gets denser still.

This is why the spurious peaks at RR Lyr landed at 0.232 $\approx$ P/2, 0.740 $\approx$ 4P/3, and 0.998 $\approx$ 7P/4. Those aren’t bugs. They’re optima.

The paper doesn’t say any of this, not because the paper is wrong, but because every one of its applications uses $\text{mp}_1$ as a feature at a single known $\tau$. Sweeping $\tau$ across a range of trial periods is just not a thing the paper does. The multi-modal optima structure is implicit in Theorem 5.1 (it’s three lines of trig) and nobody had written it down, including me until after the empirical failures forced me to look.

A quick check on the simplest possible case, a pure sinusoid at P = 1 with $\tau$ swept across the range, confirms the derivation directly:

$\text{mp}_1$($\tau$) sits at its maximum of about 4.0 across essentially the entire interior of the $\tau$ range, with sharp drops only at $\tau$ = 0 and $\tau$ = P/2 where sin(L$\tau$) itself vanishes and the ellipse collapses. The ten predicted circular optima at $\tau$ = n/22 for n = 1, 2, …, 10 all sit on a single plateau. A peak-finder picks off the boundary ones but the interior ones are indistinguishable.

The corresponding Perea-faithful diagnostic on RR Lyr, at d = 20 (the embedding dimension Perea recommends for waveforms with $\lesssim$ 10 harmonics), shows the same pattern on real data. $\pi$(P) = 4.59 $\pm$ 0.08 across 80 trial periods from 0.30 to 0.85 days, with the literature period sitting right in the flat region. The white-noise null shows how much loop structure is actually in the real signal. The framework clearly sees something (67 $\sigma$ above Gaussian noise). The phase-randomised null is nearly indistinguishable from the real signal. Phase-randomised RR Lyr is still periodic at the same frequency, just with shuffled phase relationships between harmonics, and its SW cloud has similar topology to the real one:

The framework picks out real signals from white noise at 67 $\sigma$. It doesn’t distinguish the coherent signal from its phase-shuffled counterpart. It doesn’t produce a swept landscape with a peak at the right period.

The punchline

What I had built was a sweep over $\tau$ whose output was $\text{mp}_1$. What the theory actually supports is an evaluation of $\text{mp}_1$ at a single $\tau$ picked to match an expected period. These have the same inputs and outputs. They’re not the same operation.

Max $\text{H}_1$ persistence of a sliding-window embedding is a periodicity statistic, not a periodogram. The property that makes the feature useful (robustness of $\text{mp}_1$ to the choice of $(d, \tau)$ near the optimum) is exactly what kills the periodogram interpretation. A periodogram needs a sharp peak at the right period; $\text{mp}_1$($\tau$) has a plateau of circular optima at every half-integer multiple of the true period.

Every application cited in Perea’s review uses $\text{mp}_1$ as a feature:

• SW1PerS (Perea, Deckard, Haase, Harer 2015): rank genes by circadian periodicity at the known 24-hour cycle.
• Wheeze detection (Emrani, Gentimis, Krim 2014): classify breath recordings.
• Chatter detection in machining (Khasawneh & Munch 2016): classify.
• Motion recognition from wearable sensors (Dirafzoon, Lokare, Lobaton 2016): classify.
• Video quasiperiodicity quantification (Tralie & Perea 2018): known-period scoring.

Not one of them searches for an unknown period. That pattern was the informative thing I didn’t read correctly going in.

What I take from this

Real data sooner, don’t focus on synthetic. The most expensive mistake was extensively validating every variant on clean injection tests, and spending weeks optimising hyperparameters on those tests, before putting real HARPS or Kepler data anywhere near the pipeline. When I finally did, the conclusions I’d drawn on synthetics got contradicted. Synthetic- first is attractive because it’s fast and controllable. That’s exactly what can make it a bad test bed: real-world noise structure, instrumental systematics, stellar jitter. These are exactly the things a clean injection test doesn’t contain.

Trust the early null. The core failure mode was visible in the first serious injection run. I remember thinking “that’s strange, but it’s probably the tolerance parameter.” Then “that’s strange, but it’s probably the embedding dimension.” Then “that’s strange, but the phase-folded variant will fix it.” At each step there was a local reason to keep going. In retrospect there wasn’t a global one.

Define the kill criterion before the first experiment. I kept moving the goalposts, which is what an invested and hopeful person does. Define a proper kill criterion before you invest too much time and effort into an idea, or at the very least heed the early warning signs and proceed with caution.

To finish off

There isn’t enough discourse online about failures in experimentation. There are enough self pat-in-the-backs and victory stories that create false expectations for the community. Moreover, I don’t want to lean on platitudes about resilience that turn a specific failure into a generic life lesson, or frame the failure as noble. That’s the positive-result framing in a mirror. Neither version is useful to the next person.

What I hope this post does instead: if you’re considering persistent homology on delay embeddings as a periodogram for astronomical time series, the repository contains the injection tests, the eccentricity sweep, the HARPS diagnostic benchmark, the RR Lyrae test, and the Perea-faithful d = 20 diagnostic. Decide for yourself whether to try anyway. If you have an idea that escapes the dense-optima plateau (sliding-window persistence landscapes instead of max persistence, supervised TDA features trained against a known catalog, something else), the scaffolding is a reasonable starting point. If you don’t, you’ve saved whatever fraction of my time you would otherwise have spent reaching the same conclusion.

Two things I won’t try to sell: that the failure was a learning experience in some deep sense, or that you should keep going regardless of the evidence. The first is self-congratulatory, the second is often wrong. Projects end for reasons. Writing them up is the useful part.

Repository: github.com/jvines/takens_periodogram.

---

References

• Takens, F. 1981, Detecting Strange Attractors in Turbulence, in Dynamical Systems and Turbulence, Warwick 1980, Lecture Notes in Mathematics 898, 366. The original delay-embedding theorem.
• Perea, J. A. 2019, Topological Time Series Analysis, Notices of the AMS 66(5), 686. doi:10.1090/noti1869. The review that launched this project.
• Perea, J. A., Deckard, A., Haase, S. B. & Harer, J. 2015, SW1PerS: Sliding windows and 1-persistence scoring; discovering periodicity in gene expression time series data, BMC Bioinformatics 16, 257. doi:10.1186/s12859-015-0645-6. The canonical $\text{mp}_1$-as-periodicity-score implementation.
• Bauer, U. 2021, Ripser: efficient computation of Vietoris–Rips persistence barcodes, Journal of Applied and Computational Topology 5, 391. doi:10.1007/s41468-021-00071-5. The homology engine the code is built on.
• Cumming, A. 2004, Detectability of extrasolar planets in radial velocity surveys, MNRAS 354(4), 1165. doi:10.1111/j.1365-2966.2004.08275.x. The canonical reference for LS sensitivity loss at high eccentricity.
• Hippke, M. & Heller, R. 2019, Optimized transit detection algorithm to search for periodic transits of small planets, A&A 623, A39. doi:10.1051/0004-6361/201834672. The TLS paper; documents what BLS actually misses and what adding a transit model with limb darkening and ingress/egress recovers.
• Emrani, S., Gentimis, T. & Krim, H. 2014, Persistent homology of delay embeddings and its application to wheeze detection, IEEE Signal Processing Letters 21(4), 459. doi:10.1109/LSP.2014.2305700. $\text{mp}_1$ as a feature for binary wheeze classification.
• Khasawneh, F. A. & Munch, E. 2016, Chatter detection in turning using persistent homology, Mechanical Systems and Signal Processing 70–71, 527. doi:10.1016/j.ymssp.2015.09.046. $\text{mp}_1$ as a feature for machining-chatter classification.
• Dirafzoon, A., Lokare, N. & Lobaton, E. 2016, Action classification from motion capture data using topological data analysis, IEEE Global Conference on Signal and Information Processing (GlobalSIP),
  1. doi:10.1109/GlobalSIP.2016.7906045. $\text{mp}_1$ as a feature for motion classification.
• Tralie, C. J. & Perea, J. A. 2018, (Quasi)Periodicity Quantification in Video Data, Using Topology, SIAM Journal on Imaging Sciences 11(2), 1049. doi:10.1137/17M1137642. $\text{mp}_1$ as a score at a known (or searched-by-grid) period for video recurrence.