[论文解读] Rotation method for accelerating multiple-spherical Bessel function integrals against a numerical source function

Zachary Slepian, Yin Li+2 more|arXiv (Cornell University)|Nov 29, 2019

Radio Astronomy Observations and Technology参考文献 4被引用 8

一句话总结

本文提出了一种名为'旋转方法'的新算法，通过利用三角函数的积化和差恒等式，将多重球贝塞尔函数（sBF）积分重新表述，从而实现 $N_{\rm g}\log N_{\rm g}$ 复杂度的1D FFT加速。该方法相较原始的 $N_{\rm g}^{p+1}$ 方法，实现了高达1000–3000倍的加速，显著提升了宇宙学相关函数与协方差矩阵计算的效率。

ABSTRACT

A common problem in cosmology is to integrate the product of two or more spherical Bessel functions (sBFs) with different configuration-space arguments against the power spectrum or its square, weighted by powers of wavenumber. Naively computing them scales as $N_{ m g}^{p+1}$ with $p$ the number of configuration space arguments and $N_{ m g}$ the grid size, and they cannot be done with Fast Fourier Transforms (FFTs). Here we show that by rewriting the sBFs as sums of products of sine and cosine and then using the product to sum identities, these integrals can then be performed using 1-D FFTs with $N_{ m g} \log N_{ m g}$ scaling. This "rotation" method has the potential to accelerate significantly a number of calculations in cosmology, such as perturbation theory predictions of loop integrals, higher order correlation functions, and analytic templates for correlation function covariance matrices. We implement this approach numerically both in a free-standing, publicly-available extsc{Python} code and within the larger, publicly-available package exttt{mcfit}. The rotation method evaluated with direct integrations already offers a factor of 6-10$ imes$ speed-up over the naive approach in our test cases. Using FFTs, which the rotation method enables, then further improves this to a speed-up of $\sim$$1000-3000 imes$ over the naive approach. The rotation method should be useful in light of upcoming large datasets such as DESI or LSST. In analysing these datasets recomputation of these integrals a substantial number of times, for instance to update perturbation theory predictions or covariance matrices as the input linear power spectrum is changed, will be one piece in a Monte Carlo Markov Chain cosmological parameter search: thus the overall savings from our method should be significant.

研究动机与目标

• 解决在宇宙学中对多重球贝塞尔函数积分与数值源函数进行计算时存在的计算瓶颈问题。
• 克服针对 $p$ 个配置空间变量的原始积分方法所导致的 $N_{\rm g}^{p+1}$ 复杂度问题。
• 实现快速傅里叶变换（FFTs）在sBF积分中的高效应用，这些积分原本与标准FFT加速不兼容。
• 提供一种数值稳定、高精度的方法，用于大规模结构巡天中相关函数与协方差矩阵的计算。
• 支持计算密集型任务，如需要反复重新计算积分的马尔可夫链蒙特卡洛（MCMC）参数估计。

提出的方法

• 利用三角恒等式将球贝塞尔函数重写为正弦与余弦的和。
• 应用积化和差的三角恒等式，将多配置空间变量的积分转化为单个差值与和的积分。
• 将问题重新表述为在配置空间平面上进行45°旋转，从 $(a,b)$ 变换为 $(a-b, a+b)$，从而实现1D FFT的应用。
• 即时执行生成的1D FFT，避免使用大型查表，支持内存内计算。
• 利用旋转变换将高维sBF积分转化为一系列 $N_{\rm g}\log N_{\rm g}$ 复杂度的1D FFT。
• 在公开的Python包中实现该方法，并集成至mcfit框架中，以提升广泛可用性。

实验结果

研究问题

• RQ1 能否在数值源函数存在的情况下，将多重球贝塞尔函数积分的加速效果超越原始的 $N_{\rm g}^{p+1}$ 复杂度？
• RQ2 尽管sBF积分具有固有的多维结构，是否仍可将其重述为适合1D FFT处理的形式？
• RQ3 与直接数值积分及现有基于查表的方法相比，旋转方法可实现多大的加速比？
• RQ4 在高精度宇宙学应用（如红移巡天的协方差矩阵估计）中，该方法表现如何？
• RQ5 该方法能否扩展至更高阶的sBF阶数 $\ell$，同时保持数值精度？

主要发现

• 旋转方法通过启用基于FFT的计算，将多重sBF积分的计算复杂度从 $N_{\rm g}^{p+1}$ 降低至 $N_{\rm g}\log N_{\rm g}$。
• 在测试案例中，该方法相较直接积分实现6–10倍加速，结合FFT后最高可达1000–3000倍加速。
• 与原始方法相比，数值结果一致性良好，足以确保在DESI分析等应用中系统误差可忽略。
• 该方法避免使用大型查表，支持高效的内存内计算，减少I/O瓶颈。
• 该方法在 $\ell = 0$ 和 $\ell = 2$ 多极矩中已表现有效，覆盖了2PCF协方差矩阵的关键组成部分。
• 在更高 $\ell$ 时，由于被积函数动态范围过大，存在挑战，但该方法对当前规模巡天应用仍具可行性。