Finding black hole spins efficiently during a numerical binary evolution

Abstract: The dynamics of a binary black hole system depend on its masses and spins. For a binary at finite separation, it is not possible to define these quantities in an unambiguous way; however, there are several reasonable definitions that reduce to the expected values in the limit of infinite separation. Approximate Killing vector (AKV) spin is one of the spin definitions used in the numerical relativity code SpEC. AKV spin requires finding approximate Killing vectors on an apparent horizon, which reduces to a generalized eigenvalue problem of size $\mathcal{O}(L<sup>2)$, and a direct solve has time complexity $\mathcal{O}(L<sup>6)$, where $L$ is the highest spherical harmonic mode used to represent the apparent horizon. This scaling means that the cost of computing AKV spins increases rapidly as we simulate systems at higher resolutions, especially those with high spin or mass ratios. We describe a new algorithm for computing AKV spins that is much faster than the current algorithm.

• The paper replaces a full dense generalized eigensolve with ARPACK shift-invert iteration and cached LU factorization, reducing the eigensolve from about 13,100 seconds to 18 seconds at spherical-harmonic resolution L = 100.
• The method uses the shear norm’s non-positive generalized eigenvalue spectrum to guarantee that any positive shift selects the desired three spin eigenpairs, while remaining accurate at condition numbers as high as 10¹².
• The algorithm delivers roughly a two-order-of-magnitude eigensolve improvement and about a one-order-of-magnitude end-to-end speedup, shifting the main performance bottleneck to parallelizable matrix assembly.

Motivation and context

Numerical relativity (NR) simulations of binary black hole (BBH) mergers provide the ground truth for waveform models used by LIGO-Virgo-KAGRA, and next-generation detectors will require simulations at substantially higher resolution. In codes such as SpEC, black hole masses and spins are not inputs but are measured from the evolving solution; in particular, the approximate Killing vector (AKV) spin is computed on each apparent horizon at every time step. The AKV construction reduces to a generalized eigenvalue problem $Hz = \lambda D^2 z$ discretized in spherical harmonics up to mode $L$, producing dense matrices of size $N = (L-1)^2 - 1$. Chaudhary et al. observed that this calculation had become a severe bottleneck: LAPACK's dense solver dggev computes all $\mathcal{O}(L^2)$ eigenpairs with theoretical cost $\mathcal{O}(L^6)$ — measured closer to $\mathcal{O}(L^7)$, plausibly due to memory pressure from the $\mathcal{O}(L^4)$ matrix elements — and runs single-threaded while the rest of the simulation waits. For $L = 100$, a single spin evaluation took over three hours.

The AKV spin definition

The AKV method seeks a divergence-free vector field $\phi^A = \epsilon^{AB}D_B z$ tangent to the horizon that minimizes the shear norm $\|\sigma\|^2 = \oint z H z\, dA$, where $L$0 is a fourth-order operator involving the horizon Laplacian and Ricci scalar. Minimizing subject to a norm constraint yields the generalized eigenvalue problem $L$1; after removing the shared constant kernel ($L$2), the three smallest-magnitude eigenvalues generalize the rotational Killing vectors of a round sphere to a deformed horizon. Spin magnitudes follow from normalizing the eigenfunctions via the Kerr-based condition $L$3 — an assumption that the dynamical horizons are approximately Kerr, which the authors note is the most robust normalization they know but is exact only for Kerr.

The new algorithm

The key observation is that only three eigenpairs are needed, so an iterative eigensolver suffices. The authors use ARPACK with a shift-invert transform $L$4, which converts the desired small-magnitude eigenvalues into large, well-separated ones. Because the chosen shear norm makes $L$5 positive semidefinite and $L$6 negative semidefinite, all generalized eigenvalues are non-positive, so any positive shift $L$7 is guaranteed to select the correct eigenpairs — a structural property the authors point out is absent for other norms such as Cook–Whiting's. They adopt $L$8 as a default balancing conditioning against convergence speed; notably, ARPACK remained accurate even when the condition number reached $L$9, and for large $N = (L-1)^2 - 1$0 it was more accurate than dggev near zero eigenvalues.

Each ARPACK iteration requires solving $N = (L-1)^2 - 1$1. Rather than inverting explicitly, the algorithm performs one LU factorization of $N = (L-1)^2 - 1$2, caches it, and applies back-substitution per iteration (~30 iterations typically). Although LU factorization shares the $N = (L-1)^2 - 1$3 theoretical scaling of dggev, its heavy use of BLAS Level 3 routines makes it two orders of magnitude faster in practice.

Performance results

On Caltech HPC hardware, dggev runtime grows from 0.07 s at $N = (L-1)^2 - 1$4 to about 13,100 s at $N = (L-1)^2 - 1$5, whereas ARPACK plus LU takes only 4.7 s plus 18 s respectively at $N = (L-1)^2 - 1$6 — roughly two orders of magnitude faster for the eigensolve itself. Including matrix assembly, which both algorithms require and which scales as $N = (L-1)^2 - 1$7–$N = (L-1)^2 - 1$8, the end-to-end speedup is closer to one order of magnitude, with matrix generation now the dominant bottleneck. Extrapolating measured scalings, LU factorization would overtake matrix generation only at $N = (L-1)^2 - 1$9 on HPC-class systems (or $\mathcal{O}(L^2)$0 on consumer hardware with less optimized BLAS), well beyond anticipated needs. A useful diagnostic appendix shows that without shift-invert, finding the three smallest eigenvalues of a test diagonal matrix required 31,302 matrix-vector products versus 430 for the largest; with shift-invert, only 38 inverse applications were needed, and the advantage grows with matrix size.

Limitations and open questions

The authors are candid that their original goal — a fully matrix-free algorithm using SpEC's operator-application functions $\mathcal{O}(L^2)$1 and $\mathcal{O}(L^2)$2 directly with GMRES or BiCGSTAB inner solves — failed in practice. Simple preconditioners worked most of the time but converged slowly for highly deformed horizons near merger, and since robustness was preferred over peak speed for an embedded diagnostic, they retained explicit matrix construction. A promising middle ground they identify but do not implement is reusing a fixed matrix as a preconditioner across many time steps, though rapidly changing spins near merger may force frequent rebuilds. Matrix generation is also embarrassingly parallel and could be accelerated in next-generation codes like SpECTRE. An additional unstated dependency is that the guarantee on the shift $\mathcal{O}(L^2)$3 relies on the specific norm choice; extending the approach to other quasilocal spin definitions would require revisiting this argument.

Conclusion

This paper replaces a dense generalized eigensolve with an ARPACK shift-invert scheme backed by a cached LU factorization, reducing AKV spin computation from hours to seconds at high resolution and shifting the bottleneck to embarrassingly parallelizable matrix assembly. The change removes a scaling obstacle that would otherwise have grown prohibitive as higher-resolution simulations of high-spin and extreme-mass-ratio binaries become necessary for next-generation gravitational-wave data analysis.