Spherical harmonic transform

REFERENCES

1. Rapid application of the spherical harmonic transform via interpolative decomposition butterfly factorization by James Bremer, Ze Chen, Haizhao Yang, SISC 2021.
2. Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series by Richard Mikael Slevinsky, ACHA 2019.
3. Fast algorithms for spherical harmonic expansions by Vladimir Rokhlin, Mark Tygert, SISC 2006.
4. Fast algorithms for spherical harmonic expansions, II by Mark Tygert, JCP 2008.
5. Fast algorithms for spherical harmonic expansions, III by Mark Tygert, JCP 2010.
6. A fast spherical harmonics transform algorithm by Reiji Suda, Masayasu Takami, Math.Comp 2001.

球谐函数（spherical harmonics）是单位球面 上的一组完备正交函数，可以看作 Fourier 级数在球面几何上的自然推广。取球坐标 ，其中 为从北极起算的极角（余纬）， 为经度，相应的球面面积元为 。球谐函数 是球面 Laplace–Beltrami 算子 的特征函数，满足

这里 称为次数（degree）， 称为阶数（order）。将球面上的特征值问题分离变量，可以得到球谐函数的基本形式

其中 是关联 Legendre 函数， 经过了 上的 归一化；不同文献可能采用不同的整体相位约定，尤其是负阶 的定义，但这不影响变换的主要计算结构。由球谐函数的正交完备性，任意 都可以在 意义下表示为

这种正交展开使球谐函数成为球面函数逼近与谱方法的重要基础：光滑球面函数往往具有快速衰减的球谐系数；而对于球面 Poisson、Helmholtz 等偏微分方程，上述特征值关系还可用于简化谱空间中的算子表示。与此同时，球谐函数沿经度和纬度呈现的不同结构，也决定了快速球谐变换（spherical harmonics transform, SHT）的主要计算困难：球谐函数在经度方向上已经是 Fourier 模式，可以直接使用 FFT；但在纬度方向上，不同阶数对应不同的关联 Legendre 函数族，无法把整个计算简单地交给一个固定的 FFT。真正需要加速的是这些随阶数变化的稠密一维变换。

本文是关于上面六篇论文的简单综述，它们给出了目前最有效的几种进行快速球谐变换的算法。这六篇文献的差别，主要在于它们选择了什么中间表示，以及如何利用这个表示中的结构：Suda–Takami 把递推生成的函数拆成低相对次数的两部分，再快速插值；Rokhlin–Tygert 将跨阶基底连接转化为结构化特征向量变换；Tygert II 将特征向量的含义改为特殊节点上的函数值；Tygert III 用蝶形压缩代替 II 中的特征向量算法；Slevinsky 从解析的正交旋转出发，构造到 Fourier 系数的稳定连接；Bremer–Chen–Yang 则利用函数的振荡区域与非振荡区域，把矩阵的构造和作用一并加速。

阅读时需要始终区分三个问题：一次变换需要多少运算、为反复使用它需要做多少预计算，以及误差结论针对哪一个计算步骤。 这三个问题的答案并不总是相同。以下公式除明确标注的统一说明符号外，均按对应论文的记号书写；文中的“原式”与节号用于定位附件 PDF，而不作为本文公式的连续编号。

1. 共同问题与两条路线

先用一个仅在本节及最后比较中使用的带限参数 表示最大球谐次数。记 为在 上关于 归一化的关联 Legendre 函数，满足

这里 只是统一说明的占位符，不取代各篇论文中的 、 或 。相位因子随原文约定吸收到系数中。一个带限球面函数可以写成

球谐系数一共有 个。先固定 ，括号内的计算就是 associated Legendre transform，以下简称 ALT。如果有 个纬度节点，每个 的直接 ALT 需要 次运算，对所有 求和便得到 。经度方向则是在 条纬线上分别做长度为 的 FFT，合计 。因此，快速球谐变换的主要任务是把纬向部分从三次复杂度降到接近数据量的复杂度。[1, §1；3, §1；6, §1]

下文统一用合成(synthesis)表示“球谐系数到函数值”，用分析(analysis)表示“函数值到球谐系数”。各文对 forward/inverse 的命名并不一致；例如 Bremer–Chen–Yang 把合成称作 forward。因此不能仅凭这两个英文词判断矩阵乘法的方向。[1, §3.1；3, §4–5；5, §6]

这六篇文献中的使用的算法大致可以分成值域方法与系数域方法两条技术路线：

路线	固定阶数时希望先得到	本组文献中的代表	主要结构
节点值路线	在某组纬度节点上的关联 Legendre 展开值	Suda–Takami [6]；Tygert II [4]；Tygert III [5]；Bremer–Chen–Yang [1]	分裂递推与插值、采样矩阵的正交特征向量、互补低秩、振荡区域分块
系数连接路线	公共低阶基底的系数，再转为正弦、余弦或 Fourier 系数	Rokhlin–Tygert [3]；Slevinsky [2]	Sturm–Liouville 算子、半可分矩阵、解析 Givens 旋转、Fourier 积分算子、Gamma 比值连接

这里的两条路线并不对应两种互斥的压缩算法。蝶形分解是一种工具，既可以压缩“系数到节点值”的矩阵，也可以压缩“一个基底的系数到另一个基底的系数”的矩阵。 Tygert III 和 Bremer–Chen–Yang 属于前一种用法，Slevinsky 属于后一种用法。同样，Rokhlin–Tygert 与 Tygert II 都使用快速特征向量变换，但它们对角化的算子不同，所得特征向量的含义也不同。

这一划分有助于澄清文献之间的实际继承关系：Rokhlin–Tygert 已经走过“先转成三角系数”的路线，Slevinsky 进一步改善的是连接的构造、正交性利用和压缩方式；Tygert III 继承的是 II 的节点框架，替换了其中的核心谱变换；Bremer–Chen–Yang 直接面向公共 Gauss–Legendre 节点上的矩阵，特别处理长期存在的预计算瓶颈。

正交性、求积权重与逆变换

设 、 是足够高阶的 Gauss–Legendre 节点和权重，在一个固定阶数 下令

由正交性可知精确求积将给出

因而合成是 ，分析是 。如果改用加权节点值 ，相应矩阵

满足 。当 为方阵时，它是正交矩阵；为长方阵时，它是保持系数二范数的等距嵌入。后一情形下， 是左逆， 是投影到可表示节点数据的投影算子，不能把它写成整个节点空间上的恒等算子。[4, §2.3、§3.1；5, §3.3、§5；1, 原式(25)]

Remark: 为什么奇阶关联 Legendre 函数也能使用 Gauss–Legendre 精确求积？

关联 Legendre 函数可写为

其中 是次数为 的多项式，归一化常数包含在 中。当 为奇数时，单个 通常不是代数多项式；但同阶函数的乘积

是次数不超过 的多项式。如果 ，使用 个 Gauss–Legendre 节点就足以精确计算这些内积。这才是上式离散正交性的依据。这个说明也避免继承 [1, §3.1.1] 对单个关联 Legendre 函数多项式次数的过于简略的表述。

符号约定

文献	原文主要带限或规模符号	关联 Legendre 记号	阅读时特别注意
Suda–Takami [6]	为最高次数； 为经纬方向采样数	理论主要用 ；实验用	实验归一化的平方积分为 ，与其余几篇常用的单位范数不同
Rokhlin–Tygert [3]	全局最高次数为 ；局部奇偶连接用		奇阶降到 ，正偶阶降到 ； 单独处理
Tygert II [4]	全局最高次数为 ；单个奇偶分支长度为		随阶数变化， 为公共节点； 是权重
Tygert III [5]	核心测试用阶数 与分支长度		是加权值， 是展开系数；不能丢掉平方根权重
Slevinsky [2]	全局最高次数为 ；局部公式也使用 作次数下标		相邻阶连接 ；最终降到 或
Bremer–Chen–Yang [1]	全局最高次数为		是半带限尺度；矩阵行列规模不能直接与 [3] 的 等同

下面不把这些符号强行统一。每一篇开始时先重新说明局部记号；只在最后的横向比较中使用 表示共同的最高次数尺度。

2. Suda–Takami（2001/2002）：分裂 Legendre 函数与快速插值

A fast spherical harmonics transform algorithm, by Reiji Suda and Masayasu Takami. [6]

从高次数展开转向局部函数值表示

Suda–Takami 考虑的主要问题是固定阶 的关联 Legendre 展开求值。沿用原文符号， 为截断次数， 为多项式次数， 为纬向采样点，需要计算

这就是原文式(3)。经度方向已经可以使用 FFT，因此能否加速球谐变换，主要取决于如何降低全部 对应的这些求值成本。采用与求积一致的规范化后，分析步骤可通过求值矩阵的加权转置实现，因而一套求值算法也提供构造分析算法的基础。[6, §1.1, pp.704–705]

本文的出发点是：一个部分和所需的自由度由其次数区间的宽度决定，应当用相应数量的函数值表示它，而不必始终保留全局规模的数据。 困难在于，高次数区间的部分和并不是低次多项式。若直接剥去公共因子、再用多项式系数组织递归，中间表示可能很不稳定。作者因此引入两条分裂 Legendre 分量，使高次数区间也能写成低次多项式与适当基准函数的乘积。随后，各层部分和都通过精心选择的采样值表示，层间传递则由快速插值完成。[6, §§1.2,2.1]

这里的快速工具是 FMM 加速的 Cauchy 核求和。正文并没有对关联 Legendre 矩阵构造 SVD 或 ID 蝶形分解；其核心结构首先来自三项递推给出的精确恒等式。[6, §§2.1–2.4]

为何一次快速插值还不够？假如只在 个选点上直接计算原展开，再插值到全部目标点，虽然减少了求值点数，最初那批直接求值对所有阶累加仍是三次成本。在作者采用的 个纬向点设置下，这只能把相应主项常数降到原来的 。真正的突破必须发生在选点值的计算本身：它们也要由更小的部分和递归生成，才能逐层削减总工作量。[6, §§1.1,2.3]

Step 1 用两条分裂分量表示高次数部分和

关联 Legendre 函数满足

由于这是关于次数的二阶递推，给定相邻的两个基准函数 和 ，其他次数的函数可以表示为这两个函数的多项式线性组合。原文 Theorem 1 将其写成

这分别是式(6)、(7)，组合使用这两个公式可以给出

这说明我们可以将一个阶数为的关联Legendre函数分裂成两个阶数为与的关联Legendre函数的组合，其中的 称为分裂点。除恒等于零的特殊情形外，

当 时，对应的 恒为零。由此可见，若待求和的次数位于 ，取 后，所需多项式次数随区间宽度增长，而不随区间所在的绝对次数增长。这是分治能够减少采样点数的代数依据。[6, §2.1, p.706]

不同子区间可以选择不同分裂点，但合并它们之前，必须把表示迁移到共同的分裂点。原文定义

并给出式(8)：

第一个下标 是原分裂点，第二个下标 是新分裂点；上标 表示从原通道 向新通道 传递。比值在分母非零处定义，相关关系来自函数恒等式。关键是 不依赖求和次数 ，所以同一个 点值变换可以迁移整个部分和，而无需逐项重算。[6, §2.1]

分裂点的选择还关系到稳定性。作者指出， 的分裂及 的迁移与较可靠的递推方向一致；反向移动则可能引入不稳定的次数递推。虽然把分裂点放在区间中部可以进一步降低多项式次数，实际实现仍选择区间下端点 ，因为误差对分裂稳定性十分敏感。[6, §§2.1,2.6]

例如，原文将高次数部分和写为

取 后，每条通道只含一个固定基准函数，变化都集中在低次多项式中。进一步二分高次数部分和时，两个子区间可以各自选择新的基准函数；待返回父节点，再用上述 将它们变换回共同表示。因而，“分裂”并不是近似地丢掉高频，而是在每一层精确地更换描述同一个部分和的两个基准函数。[6, §2.4]

Step 2 用带缩放的 Cauchy 插值传递节点值

设低次多项式 在采样点 上的值已知，目标是在另一组点 求值。定义

Lagrange 插值可以写成

在节点固定时，这就是“输入对角缩放—Cauchy 矩阵作用—输出对角缩放”：

首尾两步具有线性成本，中间一步可由 FMM 加速。若源点数和目标点数分别为 、——这里沿用原文插值小节的局部记号——其成本为 ，其中 是 FMM 精度参数。与采样点重合的目标点直接保留已有值，不参与奇异分母的计算。[6, §2.2, pp.707–708]

分裂分量具有形式 ，因此只需在两端缩放中再纳入 ：

因此，算法实际传递的是带有基准函数的值，避免了显式形成可能不适宜数值计算的多项式系数。不过，具有快速插值公式并不自动意味着插值稳定；采样点的选择是这条路线的另一个关键环节。

Step 3 稳定选点，再沿次数区间递归合并

对以 为基准的插值，原文引入矩阵

若矩阵元素很大，插值便可能通过大项相消产生较小结果。作者用 Cramer 法则说明：在候选节点中，若采样矩阵的行列式绝对值达到全局最大，则 。这说明稳定的节点集合具有存在性依据，但寻找全局最优集合是昂贵的组合优化。[6, §2.3, p.708]

实际算法在候选点集合 上采用加权贪心选择：从 开始，每次选取使 最大的点 ，再令 。对于两条分裂分量，作者不分别选择两套节点，而是从

开始，共同选择使

最大的节点，再同时更新两条权函数。这一 max-min 规则试图兼顾两个基准函数，避免分别取点后将采样规模约增一倍。必须区分：最大行列式集合给出的是存在性论证，而上述贪心规则的稳定性在本文中由实验支持，没有获得同样的全局保证。[6, §2.3, pp.708–709]

完整计算由两个递归过程组成：一个返回普通部分和，另一个返回两条分裂部分和。对每个次数区间，先选择适合该区间的分裂点和较小采样集合，再将次数区间二分；两个子结果迁移到同一分裂点后相加，最后通过快速插值传递到父节点需要的点上。两条分量共享递归过程，并不是将每一层展开成四个相互独立的子问题。当继续插值或二分不再划算时，转为直接求值。[6, §2.4, pp.709–711]

在固定 FMM 精度下，每层所有子区间所处理的点数总和为 ，递归深度为 ，所以单个阶的作用成本为 ；对全部 求和，再加经向 FFT，得到 的球谐变换作用成本。这个结论假设分割点、分裂点、采样集合、函数值、插值缩放系数、迁移系数及 FMM 树均已生成。原文明确给出相关预计算的总成本为 ，并把降低预计算及存储需求列为后续问题。[6, §2.5, p.711; §3, p.714]

预计算昂贵也有具体原因：贪心选点需要反复检查候选集合，插值及迁移还需要生成许多与节点相关的系数。作用时可把这些操作变为现成数据上的乘加，但首次使用仍要承担建设成本。本文没有给出足以据此宣称完整计划具有某个近二次存储界的分析；其贡献应落在重复作用的加速，而完整计划的成本和存储仍须单独评估。[6, §§2.5,3]

误差方面，FMM 的 是内部近似参数，不能直接视为最终球谐变换的统一误差。作者明确承认尚无整个算法的解析近似误差界。因此，应将它评价为具有清楚作用复杂度分析、并通过选点和递推方向改善稳定性的快速方法；完整的统一精度保证仍未完成。[6, §§2.5,3]

这一缺口与具体计算过程有关：Cauchy 求和的误差还要经过输入输出缩放、分裂点迁移以及多层合并。即使单次插值矩阵的元素较小，也仍需估计这些操作共同作用时误差如何传播。作者通过适应每个部分和的采样集合，避免所有区间被迫使用同一种不合适的表示，但尚未把这一设计转化为全过程的误差定理。[6, §3, p.714]

实验所支持的结论

实验使用 Gauss–Legendre 节点，并将函数规范化为

这一规范化与后面 Tygert 使用的单位 规范化有所不同。快速法与直接法均利用对称性；比较中的 lazy direct 还跳过低于精度相关阈值的极小函数值，快速算法内部的直接计算也采用这一处理。[6, §2.6, pp.712–713]

表1在 FMM 展开阶数为 时， 的快速时间为 秒，lazy direct 为 秒，尚未体现加速；到 ，两者分别为 秒和 秒，加速约 倍，但报告的相对无穷范数误差为 。因此，摘要中的交叉规模需要结合所要求的精度和所选直接基准理解。[6, 表1, p.712]

这些时间由抽样的 个阶估计全部模态工作量，不包含预计算，也不是含经向 FFT 的完整端到端计时；误差来自随机向量下与直接结果的比较。作者补充，全模态测试的误差可能比表中抽样结果大 – 倍。这些实验说明分裂函数和稳定选点可以使 FMM 分治取得实际加速，但不足以支持所有规模、所有输入均达到机器精度的判断。[6, §2.6, p.713]

3. Rokhlin–Tygert（2006）：低阶连接与半可分谱变换

Fast algorithms for spherical harmonic expansions, by Vladimir Rokhlin and Mark Tygert. [3]

从球谐变换到跨阶连接问题

Rokhlin 与 Tygert 的出发点仍然是在经向 Fourier 变换已经可以交给 FFT 的前提下，尝试同时加速各个角阶的 associated Legendre 变换。他们没有直接压缩函数在纬向节点上的采样矩阵，而是先研究不同角阶的正交基之间能否快速换基。原文以 表示通常的 associated Legendre 函数，以 表示关于 归一化后的函数，满足

球谐展开截断到次数 ，共 个系数。分离经向变量之后，需要计算的基本对象是

本文默认使用等角中点纬度，其对应的 节点为 。这些节点适合执行 DCT 和 DST，因而只要能快速把 associated Legendre 系数变成相应的三角级数系数，就能完成纬向变换。[3, §1、§2.1，式 (1.4)–(1.8)、(2.5)–(2.6)，pp.1904–1908]

作者选择的中间基是阶数 和 的 associated Legendre 函数。这个选择来自端点行为：正偶数阶 的函数包含因子 ，可以在相同次数奇偶的 中有限展开；奇数阶 的函数包含 ，可以在相应的 中有限展开。因此，这一步是有限维函数空间内的精确换基，并没有用低阶函数近似替代高阶函数。 则是普通 Legendre 变换，可以单独使用已有的普通 Legendre 快速转换工具，不应放进要求端点消失的阶数 基中。[3, Lemmas 2.1–2.4、§2.3，pp.1908–1913]

Step 1 把微分算子写成对角加半可分矩阵

关键数学工具是 associated Legendre 方程的自伴形式。原文定义

不同 的微分算子只有一个乘法项不同。例如，。以 均为偶数且 为例，原文用 表示 （注意这里的上标不是平方而是表示），得到 矩阵

由于基函数已经对角化 ，该矩阵具有分解

更强的结论是， 的非对角部分并不是一般稠密矩阵，而是精确的秩一半可分结构：

其中

因而可以由对角矩阵与一个秩一半可分矩阵相加得到（Diagonal-Plus-Rank-One, DPR1），于是任何严格位于对角线一侧的连续矩形块都可以写成一个列向量与一个行向量的外积。这里的秩一是解析恒等式，不是规定容限下观察到的数值低秩。[3, Lemma 3.2，式 (3.9)–(3.14)，p.1916]

这个结论把原本需要逐项生成的稠密矩阵，压缩成了一个对角向量和两条生成元。直接应用 时，半可分部分可以通过前缀和、后缀和在线性成本内计算。不过，球谐变换所需的是 的特征向量矩阵，而不是乘以 本身，因此仅仅发现半可分结构还没有完成快速变换。真正需要衔接的下一步，是让特征分解也保留这种递归可压缩性。这一区分对于理解本文的复杂度来源很重要：线性成本的矩阵乘法，通常不能单独推出线性成本的特征基变换。

这一结构与三项递推的关系也很清楚。在固定次数奇偶的低阶基中，乘法 由三对角矩阵表示；与之互逆的乘法 对应三对角逆所具有的半可分结构。原文先在无限维乘法算子层面论证，再取有限主子块，因此不能把 直接等同于未经边界修正的有限三对角截断矩阵的逆。其余三个奇偶组合只改变基函数的起始次数；奇数阶 则改用 和系数 。 的对应修正为零，可以直接处理。[3, §3.1，式 (3.15)–(3.35)，pp.1916–1918]

Step 2 用特征向量实现高阶到低阶的连接

设目标函数在低阶基中的展开为

将它代入 ，再与低阶基取内积，就得到

所以， 的相应特征向量就是高阶基函数的连接系数。由于两套基均关于 正交归一，这些连接列也正交归一。原文将它们补成实正交矩阵 ，并把所需物理次数对应的列排在前面。[3, Lemma 3.1，式 (3.5)–(3.8)，pp.1915–1916]

具体地，若

则原文给出的关系是

这里 已补零到长度 。真实高阶空间通常比低阶空间小，因此实际换基只使用 的部分列；余下特征向量构成正交补。这个细节解释了为什么可以调用方阵的正交谱变换，同时仍只保留有限个高阶系数。[3, Lemma 3.3，式 (3.36)–(3.40)，p.1919]

从数值角度看，这种中间表示的吸引力在于能量规范始终清楚：连接前后的函数相同，两套基都使用同一个 内积，因而物理系数向量的二范数在换基中不变。与先求导得到高次幂系数、再重新组合的思路相比，它不需要显式引入尺寸悬殊的单项式系数。与此同时，正交补也提醒我们：任意低阶系数未必属于目标高阶空间，分析方向是取出相应物理分量；只有当输入确实在目标空间中时，其余分量才严格为零。这是有限维带限变换的空间约束，不应被忽略。

Step 3 快速作用谱连接，再转为三角系数

到这里，问题转化为快速应用 的特征向量矩阵。原文调用 Chandrasekaran–Gu 的半可分分治谱算法及相关的 Gu–Eisenstat 工具，给出 、 的作用成本 。其基本原因是：分治后的谱合并可归结为对角矩阵的低秩修正，秩一修正的特征向量分量具有缩放 Cauchy 核形式，因而可以借助快速求和完成各层合并。本文自身的主要贡献在于建立 associated Legendre 换基与这个结构化谱问题之间的联系；它没有在正文重新展开底层 Cauchy/FMM 合并程序的全部实现细节。[3, Observation 2.11、Remark 2.12，p.1914]

具体说，半可分矩阵按连续索引分为两块时，连接左右子块的非对角块仍为秩一；在子问题的特征基中，这种小秩耦合仍然是小秩耦合。分治算法由此可以先处理两个子问题，再通过结构化谱更新合并。合并矩阵虽然仍然稠密，但其条目不再是任意数据，而由两套实谱点之间的 Cauchy 核和少量缩放量控制。快速方法保存的是这些可快速应用的合并表示，再逐层将输入传递上去。最后把所有合并步骤逆序转置，就能得到相应的转置作用，而不必重新压缩另一张不相关的矩阵。

因此，本文同时使用了精确代数结构和规定精度的解析近似：半可分表示以及连接系数的特征向量解释是精确的，而底层快速核求和允许误差容限 。将整套算法称为“完全不使用近似的代数快速变换”并不准确；但将它理解为一般稠密采样矩阵的黑箱低秩拟合，也会掩盖其最重要的结构来源。

高阶到低阶的连接还不是完整变换。阶数 的展开接到 ，阶数 的展开接到 ，利用

此处 是第二类 Chebyshev 多项式，与前面的特征向量矩阵 不同。原文将 Chebyshev 系数到低阶 associated Legendre 系数的连接记为 ，反方向记为 ，，符号 区分奇偶块。[3, §2.2–2.3，pp.1908–1913]

这些连接矩阵也有显式解析结构。例如，若

则 、。记

矩阵 的非零项为

这种显式的差指标、和指标及对角缩放结构，使作者可以调用 Alpert–Rokhlin 的快速 Legendre–Chebyshev 转换工具。八个低阶连接矩阵的作用成本均为 。阶数 的连接矩阵为长方形，因为对应的多项式子空间满足端点约束，不能将其当成任意 Chebyshev 系数上的方阵逆。[3, 式 (2.21)、(2.26)–(2.27)、Lemma 2.6、Observation 2.10，pp.1909–1913]

合成连接中的三角形零结构也有直接的多项式解释：一个固定次数的基函数，在同类 Chebyshev 展开中不可能产生比自身更高的次数。非零区域则由上述 Gamma 比值给出，可以按其解析变化规律快速组织计算。本文没有逐阶累积大量低阶到高阶的升阶矩阵，而是一次通过 把高阶系数送到公共低阶基，再调用少数几种低阶连接。这样，不同角阶可以共享低阶前端；真正需要随 改变的是结构化谱连接部分。

完整流程、贡献与实际限制

分析方向的步骤因而是：先用镜像节点拆分奇偶部分，再通过 DCT/DST 得到相应三角系数；应用 得到阶数 或 的系数；最后应用 ，取出目标阶所需的系数。原文在奇偶拆分时使用 的规范。合成则逆序执行：高阶系数补零、应用 、应用 ，最后做三角变换。经向 FFT 在外层处理全部 。[3, §4–5，pp.1921–1923]

这里还需要区分“逆变换”和“直接取整个合成矩阵的转置”。高阶到低阶的正交连接确实可以用转置逆转，但 Chebyshev 系数和 associated Legendre 系数所对应的自然内积不同，所以原文分别构造 与 ；网格值到三角系数也有自己的归一化因子。只有把这些尺度全部纳入，才能得到正确的分析流程。相同原因也决定了不同阶段应分别检查误差，不能仅凭 就判定完整系数—值变换的全部误差性质。

固定阶的总作用成本是 ，全部球谐变换为 。但这并不是论文代码的首次使用成本。作者为简化实现，采用了每个阶 的预计算，并明确说明理论上可依据引用技术降为固定精度下的 ；表格没有实现后一结果。若为所有阶分别建立该版计划，计划总量仍可达到 。[3, Remark 4.1、§6，pp.1922、1926]

此外，中间低阶空间的维数主要由最大次数决定。当 很大时，实际非零高阶系数已经很少，原文仍通过补零后的较大正交矩阵完成连接。这不改变上述复杂度上界，却会影响有限规模的常数，也解释了后续直接采用本阶节点与本阶维数的分解为何值得研究。算法比较因此既要看渐近阶数，也要看每个角阶实际经过了多大的中间空间。

数值实验使用双精度 Fortran，但目标只有约六位正确数字，测试的是单阶、单个次数奇偶块。作者的直接参考时间是 稠密矩阵向量积时间的三分之一，用以估计当时 semi-naive 方法的成本，并非直接运行同一完整 SHT 软件所得。表中交叉约在 至 之间。例如 Table 1 在 时，快速作用为 秒，直接参考为 秒，预计算为 秒；到 ，三者分别为 、 和 秒。表中的相应合成相对 RMS 误差约为 和 。[3, Tables 1–8、Observation 6.2，pp.1924–1926]

这些结果说明，在线算法确实开始呈现低于二次的逐模增长，并在大规模上取得作用优势，但计划的摊销需要很多次重复使用。它们也不能给出所有角阶在同一完整球面变换中的统一最坏误差，因为表格只选取了几个代表性的 与次数奇偶组合。由于不同文献使用的截断参数、目标精度、直接程序和计时边界并不相同，这里的交叉规模更适合说明该实现的历史表现，不宜直接当成所有快速球谐方法的普遍门槛。

这篇工作的价值，是用精确的微分算子结构把球谐变换接入快速正交谱作用，避免不稳定的单项式中间坐标，并为 associated Laguerre、spheroidal 等函数的推广提供共同框架。其局限也应按原文保留：正交连接矩阵本身条件良好，并不等于低阶连接、节点缩放、谱计划构造合在一起已有统一的机器精度后向稳定性保证；约六位精度的实验也不能用来证明这一更强结论。就实际使用而言，较大的作用常数与二次逐模预计算，使这一版本更适合用来理解结构和验证大规模渐近优势。[3, §6–7，pp.1926–1927]

4. Tygert（2008，II）：三对角采样变换与公共网格插值

Fast algorithms for spherical harmonic expansions, II, by Mark Tygert. [4]

从跨阶连接转向本阶节点上的变换

Tygert 的第二篇工作保留了快速结构化特征向量作用这个工具，但改变了需要对角化的矩阵。2006 年的矩阵表示微分算子，特征向量给出不同角阶之间的连接系数；2008 年的矩阵来自乘法 的三项递推，特征向量直接给出同一角阶的基函数在特殊节点上的归一化值。因此，本文不必先经过阶数 或 ，而是先在最适合该阶的节点上完成变换，再将值插值到所有阶公用的网格。[4, §1、§2.4、§3，pp.4262–4263、4267–4269、4273–4275]

原文把球谐展开截断到次数 ，系数记为 ，公用纬向节点为 的 个 Gauss–Legendre 根。正半轴的公共节点记为 。固定 后，先按关于 的奇偶性拆成

由于 ，第一支是偶函数，第二支是奇函数；这里说的是函数奇偶，不一定是次数 的奇偶。实际带限变换中两支长度可以相差一。原文在叙述单支算法时用 表示其维数，下面也沿用这个局部记号。[4, 式 (4)–(13)、(83)–(87)，pp.4261–4262、4272–4273]

Step 1 从乘法 的递推构造三对角矩阵

在同一个奇偶子序列内，associated Legendre 函数满足

其中起始次数 没有第一项。沿用原文此处表示递推次数的 ，系数为

据此构造两个 实对称三对角矩阵：

使用 的作用，是把两个奇偶子空间分开，使每个矩阵只对应正半轴节点，矩阵规模也随之减半。[4, Lemma 2.12、式 (39)–(44)，p.4267]

令 为 的正根， 为 的正根。定义

把递推式代入这些根后，最后一项因下一阶函数为零而消失，于是

根平方互异，所以对应于不同特征值的归一化特征向量正交。这样， 都是实正交矩阵，且采样矩阵与特征向量矩阵之间只差对角缩放。[4, Lemma 2.13，式 (45)–(53)，pp.4267–4268]

这是正交多项式、Gauss 求积与 Jacobi 矩阵之间经典联系在本文中的具体应用。一般读者可以这样理解：三项递推描述“乘以坐标”如何改变展开系数，Gauss 节点则恰好使截断后的最后一项消失，所以在这些特殊坐标处，整列函数值自动满足一个有限维特征向量方程。归一化以后，不需要再额外证明整张采样矩阵的复杂正交恒等式，实对称矩阵的谱定理已经提供了这一性质。本文加速的正是这个有限维特征基变换，而不是在每个节点上独立递推计算所有多项式值。[4, §2.4，pp.4267–4269]

Step 2 通过正交谱变换取得本阶节点值

原文定义对角矩阵

若 、，则合成和分析分别为

这也是全文最核心的矩阵分解。 的行对应基函数，列对应节点，故合成用 ； 才是与系数具有相同欧氏范数的加权值，原始采样映射 本身不正交。后来的 Tygert III 改用显式乘平方根求积权的形式，正是为了直接呈现这里的正交结构。[4, Lemma 2.14、式 (49)–(59)、(89)–(90)，pp.4268–4269、4273；5, Remark 5.1]

从这一分解还可以看出计划阶段需要固定哪些信息：除了三对角矩阵的递归谱数据，还包括节点顺序、特征向量的符号以及恢复原始函数值的对角尺度。一般特征值程序返回的归一化特征向量具有任意列符号，而这里的符号已经由 的定义确定；若不作一致化，得到的仍是正交矩阵，却不再是所需的函数采样变换。这是把抽象快速谱程序接入特殊函数算法时必须处理的一步。

为了快速应用 、 及其转置，本文调用 Gu–Eisenstat 三对角分治谱算法。其合并机制可以用一个标准秩一修正说明：对子块完成谱分解后，合并问题成为 ；新特征值满足

相应的归一化特征向量项具有 的形式。这里帽号仅为说明所加。这个矩阵是缩放 Cauchy 矩阵，可由 FMM 快速作用；每一层处理的总维数与 成正比，递归深度与 成正比，因此得到固定计算精度下的 作用。稳定求解 secular equation、归一化和近重合情形处理则由所调用的谱算法承担。[4, §2.2，尤其 Remark 2.5，pp.4263–4265]

Step 3 用 Christoffel–Darboux 插值连接公共网格

仅在 、 上得到值还不够，因为这些节点依赖 ，而经向 FFT 需要所有角阶共享同一组纬圈。本文利用 Christoffel–Darboux 恒等式完成这一步。偶支的恒等式为

这里直接使用时要求 。它由三项递推相乘相减、求和后的望远镜消去得到。再将其中一个变量取为 的根，就只剩一个乘积。[4, Lemma 2.15，式 (60)–(62)，pp.4269–4270]

原文的特殊节点系数与公共节点系数分别为

其中 已合并正负节点的贡献，为通常单个正 Gauss–Legendre 节点权重的两倍。由精确求积和再生核性质，本阶根到公共根的插值为

反向则为

反向式来自公共 Gauss 求积，并非对前一个长方形插值矩阵求普通逆。两式都由预计算的左右对角因子和核 组成，因此再次可以使用 FMM，固定精度成本为 。节点重合时应直接取值或使用连续极限，不能计算形式上的 。[4, 式 (32)、(37)、(65)、(71)、(91)–(92)，pp.4266、4270–4271、4274]

插值之所以可以这样进行，不是因为任意节点数据都满足这两条关系，而是因为输入来自指定的有限维 associated Legendre 空间。在这个空间内，正交展开先把函数写成再生核与自身的积分；适当的 Gauss 求积对该乘积精确，于是积分变成有限和；最后使用 Christoffel–Darboux 公式把含有 项的再生核化成一个 Cauchy 分母。这个顺序同时解释了权重的来源和多项式次数条件，也说明为什么无需显式形成两个稠密评价矩阵再相乘。

对球面分析而言，原始节点数据可以先理解为相应带限离散问题的数据；超出带限空间的分量需要按所选离散投影或截断规范处理。上述精确等式描述的是带限部分的转换，不意味着有限节点能够无损恢复任意连续函数。同时，公式中看似很大的 Cauchy 项与很小的边界阶函数值可能相互补偿，因此左右缩放和近场处理是数值实现的一部分，不能在证明代数恒等式之后随意重排为缺少尺度保护的浮点计算。

奇支完全类似：将最高阶 改成 ，次高阶 改成 ，将 改成 ，并使用

对应的双向公式见 (66)、(72)。离散正交关系还给出 、，将这里的求积权与前面的对角归一化联系起来。所以本文的 FMM 有两种用途：加速谱分治的 Cauchy 合并，以及加速不同节点集之间的 Cauchy 插值。[4, §2.1、§2.2、§3.2；5, Remark 5.1]

Remark: 求积规范、后续勘误与两支维数

2008 年式 (35)、(36) 的右端漏掉了端点权因子，这一点已由 Tygert III 的 Remark 3.8 明确纠正。对于相应次数范围内的偶多项式 ，正确关系是

以及

其中 。前一式允许 ，后一式允许 。由于 associated Legendre 函数的乘积已经含有 ，前述以函数值写出的插值式仍使用 ，不能再重复乘一次端点权。[4, Lemma 2.10、式 (32)–(36)，p.4266；5, Remark 3.8]

原文为同时叙述两支的反向插值，采用共同充分条件 。若按最高次数 精确分配维数，则偶函数支为 ，奇函数支为 。由同一 Gauss 求积的次数检查，偶支的条件可放宽为 ，奇支仍为 。这是对原文公式应用范围的补充说明，不能将两支总是设为同一个长度。[4, Lemma 2.17，p.4271]

Prüfer 变换与预计算

若在每个节点都从低次数递推到最高次数，根、权和插值因子的生成仍会带来二次工作。本文为此使用 Sturm–Liouville 方程的 Prüfer 变换。对 ，在 的区间内定义连续相位 和幅度 ：

零点对应 ，而相位和幅度满足

当相位单调时，可以以 为自变量积分，从一个零点连续推进到下一个零点。这里的 是局部 Prüfer 相位，不是球面余纬。对于 ，取 、；在根处还有 ，因此幅度也直接给出所需求积系数。[4, Lemma 2.18、式 (75)–(82)、(93)，pp.4272、4274]

作者使用高精度显式预测—校正 ODE 积分器，生成零点、权和插值所需的边界阶函数值，并报告这些辅助预计算随 线性增长。高阶情况下 在转折点外会变号，因此不能把上述实 Prüfer 公式无条件用于整个区间；其零点定位是在相应振荡区进行，并结合原 ODE 处理函数值。[4, §3.3，pp.4274–4275]

这里的相位并不承担球谐数据的快速作用，也没有把待变换函数改写成一个可直接调用 FFT 的振荡积分。它服务于计划阶段：相邻零点的相位差固定，零点搜索由此避免在整个区间反复试探；幅度方程则把权重计算与根定位一并完成。在线变换仍由前面的正交谱作用和 Cauchy 插值构成。将这两个阶段分开，有助于判断算法的实际收益来自哪里，也避免把辅助节点生成的线性成本误认成全部计划都已经线性。

完整成本与稳定性应如何评价

合成的完整流程是“系数经过 和 到本阶节点值，经过 Cauchy 插值到公共节点值，恢复正负节点与各个 Fourier 模，再做经向 FFT”。分析逆序进行。按原文固定计算精度的口径，每个奇偶块的谱作用为 、插值为 ，全部球谐变换为 。更细的容限依赖中，谱计划为 ，已有计划后的谱作用为 ，插值为 。[4, §2.1–2.2、Remark 2.3、§3.1–3.2]

但是，实际代码仍使用逐模 的谱预计算，且所有计划均用软件四倍精度；在线作用才是双精度。作者明确说明，采用引用文献中的更完整方案原则上可以使谱计划近线性，但本篇未实现。表中 的本阶奇次数合成耗时 秒，稠密乘法为 秒，计划却需 秒；到 ，快速作用为 秒，计划为 秒。后者对应的稠密参考 秒由较小规模外推。网格插值另行计时，不能把该谱作用速度直接当作完整 SHT 的速度。[4, Tables 1–6、Remarks 4.1–4.2，pp.4275–4278]

若对所有角阶分别建立这种二次谱计划，将各支的 相加，仍然得到关于全局带限参数 的三次总计划成本；辅助节点和插值因子的线性预计算并不能消去这一项。相反，若谱计划真正换成原文所引用的近线性构造，总计划才有可能与近二次的完整作用相匹配。该区别对多次重复变换尤其重要：已有计划后的在线优势，与一次或少数几次使用时的总耗时优势，是两个需要分别检验的问题。

稳定性方面，数学矩阵 的正交性提供了良好基础，但谱计划的误差分析还含最小特征值间隔 ：原文给出的分量绝对误差尺度为 ，并说明省略了小的 次幂。高精度计划正是为改善最终准确度而采用的措施。此外， 把加权值还原为原始值，也会改变误差大小。数值表使用“输出误差的 RMS 除以输入 RMS”，例如上述合成的误差从 增至 ；它不是全规模统一的机器精度后向误差保证。[4, §2.2、Remarks 2.2–2.3、Table 1、§4，pp.4264–4265、4275–4278]

相较于 2006 年的方法，本文的直接三对角采样分解减少了中间换基，使较小规模上的在线变换更有竞争力，并把“快速变换”和“公共网格适配”清楚分开。代价是计划仍重、实现依赖高精度预计算，而且这种依赖多项式三项递推的构造不能直接覆盖 2006 年可处理的所有 spheroidal 情形。其贡献应归结为一个更高效的分解和在线算法，而不是已经解决了所有参数下的低成本计划与严格全流程稳定性问题。[4, Remark 4.2，p.4278]

5. Tygert（2010，III）：ID 蝶形与深度优先预计算

Fast algorithms for spherical harmonic expansions, III, by Mark Tygert. [5]

一种保留原有球谐流程的替代方案

Tygert III 延续第二篇文章的球谐变换框架：经度方向由 FFT 处理，纬向先在依赖阶数的专用节点上完成正交变换，再通过快速插值连接到公共球面网格。本文替换的是其中的核心模态变换，即用 ID 蝶形代替第二篇采用的三对角谱分治。作者希望保持快速的正向及逆向作用，同时使压缩预计算更容易实现，并避免该压缩步骤为了接近机器精度而依赖扩展精度算术。[5, §1, p.6181; §5, p.6188]

另一个实际问题是矩阵元素的生成方式。关联 Legendre 函数可以在一个固定节点上按次数递推，廉价生成整条采样向量，但任意一个高次数元素未必可以被廉价地独立查询。本文因此不仅讨论如何压缩矩阵，还设计了深度优先的预计算次序，避免为保留已经计算出的元素而存储完整稠密矩阵。[5, §§1,4]

Step 1 先建立加权正交坐标

原文用

表示规范化关联 Legendre 函数，其中 是通常的 Legendre 多项式；这一定义不额外包含 。固定 时， 在 上关于通常积分内积正交归一。 的奇偶性决定函数关于 的奇偶性，因此可以分别处理次数 和 两组。[5, §3.3, 式(19), p.6185]

对长度为 的偶支，取 的正根

对奇支，取 的正根

对应权重为

权重的上述显式定义见原文式(26)、(27)，使用这些权重的 Gauss–Jacobi 求积规则见式(29)、(30)，已经包含将正负节点配对后的因子。它们随 变化，不能直接当作公共 Gauss–Legendre 网格上的权重。[5, §3.3, p.6186]

本文加速的两个变换分别为

以上编号沿用原文。 是系数， 是经过权重平方根缩放的节点值。这个缩放具有实质作用：由 Gauss–Jacobi 求积的精确性，

奇支同样成立；奇函数在中心节点 的值为零，所以无需在该方阵变换中增加中心项。因此，两组加权矩阵均为实正交矩阵，逆变换就是转置。III 还专门指出第二篇部分求积公式漏写了权函数因子，并说明未加权节点值与这里的正交坐标之间需要一个对角变换。[5, Lemma 3.7、Remark 3.8, p.6186; §5、Remark 5.1, p.6188]

Step 2 在互补尺度的子块中构造 ID

蝶形分解并不要求整个振荡矩阵具有低秩。其出发点是：当次数范围较宽时限制节点范围，当节点范围较宽时限制次数范围，相应矩形子块可能仍只有较少的有效自由度。 原文对一般 矩阵 的假设是：任意至多包含 个元素的连续矩形子块，都能在给定精度下以秩 近似。这就是文中的 rank property，也就是互补尺度矩形块的低秩性质。[5, §3.2, p.6184]

每个局部近似采用插值分解

其中 由原矩阵的 个实际列组成， 中含有一个 单位子矩阵。这里 暂时表示待压缩的局部矩阵。保留实际列使后续层可以继续从已有骨架列中选取代表列，而无需引入一套全新的全局基。[5, §3.1, 式(6)]

原文给出的 ID 存在性结果控制了插值系数及残差；实际使用的高效版本允许插值系数绝对值至多为 ，并对一个 局部矩阵给出

其中 是 的第 个奇异值，局部维数 与球谐阶的 是不同语境的记号。这说明在局部矩阵确实容易被低秩逼近时，可以选择受控的插值表示；它本身并没有证明关联 Legendre 子块的奇异值一定快速衰减。实现中，各块的秩按照所要求的精度自适应选择。[5, Remark 3.2, 式(8); Remark 3.4]

Step 3 合并、切分并重新压缩，形成蝶形图

设 是等高、左右相邻的两个矩形块，已有

把两块左右拼接，再沿行索引切成上下两半：

这样得到的 宽度加倍、高度减半，元素总数与原来的单个块相同，因而仍处于低秩假设适用的尺度。对已经选出的骨架列做同样操作，

然后只需把这两组 个候选列再次压缩为约 个代表列：

例如，上半块最终具有表示

这就是原文式(11)–(18)的合并、切分、重压缩步骤。重复进行后，最初的高窄块逐层变成矮宽块，局部矩阵的面积大体不变，各层以小插值矩阵相连。[5, §3.2, pp.6184–6185]

算法从约 的块开始，共有 层。每层约有 个 传递矩阵，因此在上述秩假设下，作用和分解存储为

若固定精度下有效秩保持较小，且 与 同阶，便得到所希望的近线性作用。转置只需反向通过同一计算图，并转置其中的小矩阵；不需要为逆方向重新构造一套蝶形分解。[5, §3.2、Remark 3.3]

Step 4 以深度优先次序建立计划，再接入原有球谐流程

自然的广度优先建设方式会先生成并处理全部最高块，再进入下一层。如果单个矩阵元素不能被廉价重算，为避免重复生成数据，这种做法可能需要保存 个元素。本文改用深度优先次序：已生成的相邻块一旦能够合并，就立即合并、切分并重压缩；只有充分处理了已有左侧块，才生成下一个新的最高块。[5, §4, p.6187]

已经处理完的块只保留小插值矩阵及骨架列索引，不再保留长的骨架列值；等待合并的块才需要暂存这些值。由于每次处理都将列的高度减半，同时活跃的长列数据受到控制，预计算峰值内存可降至与蝶形表示相同量级的

这一改进解决了预计算的内存问题，但没有把其时间也降到近线性。本文压缩预计算仍需 运算，实验中其常数受有效秩影响。[5, §4; Remark 6.2; §7, p.6192]

生成所需采样数据时，作者使用保持奇偶性的三项递推

其中

最初两项省略不存在的低次项。在一个固定节点上，递推生成所有相关次数的函数值，并通过追踪指数处理极小值，避免下溢。按式(35)、(36)的行列方向，这生成的是求值矩阵的一整行；可以对其转置建设蝶形分解，从而满足深度优先算法所需的廉价整列生成条件，而不影响双向作用。[5, §3.3, 式(31)–(34); §5, p.6188]

整个 SHT 仍需要第二篇的附加步骤：专用节点及权重由 Prüfer 坐标下的微分方程积分生成，专用节点与公共 Gauss–Legendre 节点之间则通过对角缩放的 Cauchy 矩阵及 FMM 插值连接。III 没有取消这些步骤，只替换其中的正交模态变换。[5, §5; 4, §3.2, 式(91)–(92)]

原文全局展开的最高次数为 。在固定精度、有效秩适当受控的前提下，对全部阶及奇偶块累加，再计入插值和 FFT，可得到预期的 全变换作用成本。相应地，若逐模态使用本文的二次压缩预计算，总建图成本仍为三次量级；若保留全部蝶形表示，存储为二次乘对数的量级。这些全局成本是对单模态分析求和的结果，应与本文实际测量的单模态时间区分。[5, §§3.2,5,7]

理论边界与实验结果

这篇工作的证据有三个不同层次：加权节点矩阵的正交性是精确数学事实；在给定秩性质下，ID 蝶形的快速作用与存储结构有明确分析；关联 Legendre 矩阵是否对所有规模和阶数都具有足够强的统一秩性质，本文尚未证明。作者在结论中明确承认这一缺口，并将数值实验视为进一步建立理论的依据。[5, Remark 6.3, p.6191; §7, pp.6191–6192]

这里的秩性质在原文中以局部块的 Frobenius 范数近似表述，ID 则给出相应的谱范数残差控制。局部压缩可靠与全局成本足够低是两个相关但不同的问题：预计算能够根据误差要求增加秩，以保持局部近似质量；若某些块需要越来越高的秩，作用常数及增长率就会受到影响。本文的自适应建设使这类困难可在实际秩和存储量中观察到，却没有预先排除它在未测试规模出现的可能。[5, §§3.1–3.2; §7]

实验对各块自适应选秩，首层块宽取 。例如偶支 从 增至 时，平均秩仅从 增至 ，但最大秩由 增至 。前者支持压缩在平均意义上有效，后者也提醒读者不能把实验解释成所有子块已被同一小常数严格界住。[5, 表1, p.6189; Remark 6.1]

在单核 GHz Intel Core 2 Duo 上， 的偶支前向与逆向作用分别约为 秒、 秒，而节点权重预计算约为 秒，压缩预计算约为 秒。这里的前向是系数到加权节点值，不能与全文开头“值到系数”的 forward SHT 命名混为一谈，也不能将这些局部作用时间理解为完整球谐变换时间。[5, 表1–2, p.6189; §6, pp.6190–6191]

误差测试使用单位 范数的随机输入。表中的 是蝶形输出与稠密乘法输出的最大分量差， 则是连续进行蝶形正向和逆向后的往返误差。例如 时，两者分别为 与 ，说明前向对稠密吻合并不能替代完整往返检验。[5, 表4, p.6190; §6]

最后，论文所称“不需要扩展精度”应准确落在蝶形压缩及其作用上。报告实验中，生成节点和权重的 Prüfer 积分仍使用了扩展精度，以多保留约两位有效数字；递推和 ID 压缩则使用双精度。本文确实降低了核心压缩过程对高精度算术的依赖，但没有以全流程均未使用高精度的实验来证明这一点。[5, §5, p.6188; §6, p.6191]

6. Slevinsky（2019）：正交连接与双 Fourier 系数变换

Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series, by Richard Mikael Slevinsky. [2]

Slevinsky 的出发点是，球谐展开与双 Fourier 展开各自适合不同的计算任务。球谐函数对角化球面 Laplace–Beltrami 算子，便于在系数空间处理微分算子；双 Fourier 展开则容易借助 FFT 在等角网格上求值，便于计算局部非线性，也能与已有的球面函数表示相衔接。因此，文献[2]把问题表述为两种展开之间的快速基转换，首先研究连接矩阵本身的条件数与代数结构，再设计它的快速作用方法。这一表述也避免了将特定 Gauss 求积规则作为变换的组成部分。[2, §1]

两级转换与共同的正交空间

沿用原文， 表示球谐展开的最大次数， 是余纬， 是经度。原文式(4)把球谐函数写成

这里的相位因子属于本文的符号约定：对于非负整数 ，它等于 。固定 后， 在 中正交归一化。下文只需讨论 的连接，负阶由对称性处理。算法利用关联 Legendre 函数与 Gegenbauer（ultraspherical）多项式的关系

其中 为 Pochhammer 符号。当 为偶数时，这些函数是代数多项式，能够有限地展开到 ；当 为奇数时，它们是 乘代数多项式，能够有限地展开到 。因此，原文图1中的第一步是将各个高阶层转换到阶数0或1，第二步才将这两类低阶展开转换为余弦或正弦展开。[2, §2, 式(4)–(5), §3]

第一步具有一个直接影响数值稳定性的性质。若 是两个正交归一化函数族之间的有限连接矩阵，其列给出源基函数在目标基中的展开，那么 Parseval 等式给出

源空间通常是目标空间的一个真子空间，所以 一般是具有正交列的矩形矩阵。对来自源空间的数据，转置能够恢复原系数；对任意目标系数， 是投影，而非整个目标空间上的恒等映射。这既说明了连接问题的条件数为1，也说明了“反向转换”的适用范围。尤其不能据此认为任意双 Fourier 系数与有限球谐系数之间存在双射，因为完整转换的第二步还涉及不同权重下的基，且一般的双 Fourier 数据未必满足相应的球面正则性约束。[2, Theorem 2.2及其后的讨论]

Step 1 用解析 Givens 旋转构造相邻阶连接

为了构造高阶到低阶的连接矩阵，作者首先解决相邻同奇偶阶数 的问题。取 个源函数，并以 个目标函数表示它们，即

这些连接系数并不需要通过数值积分生成。Gegenbauer 多项式的参数提升关系对应一个上三角带状矩阵，乘以 的关系则对应另一个带状矩阵。组合这两个关系，得到“上三角带状矩阵的逆乘带状矩阵”的表示。带状矩阵逆的严格上三角部分具有半可分结构，而基函数的奇偶性使连接矩阵出现棋盘状的零元素。在正交归一化之后，可以用一串只相隔两个坐标的 Givens 旋转消去其下方非零元素。[2, Theorem 2.4]

具体地，令 只改变第 与第 个坐标，其非平凡的二阶子矩阵为

其中正弦和余弦都有显式公式：

连接矩阵的有限维分解于是为

这里用 明确表示列数，从而使旋转个数与矩阵尺寸一致；其正弦、余弦公式就是原文式(11)。应用时，先在输入末尾补上两个零，再从右向左应用旋转；转置作用则反向遍历并转置各个旋转，最后保留相应坐标。每次相邻阶数转换只需 次运算。[2, Definition 2.3, Theorem 2.4]

Remark: 相邻阶连接所用的 Gegenbauer 恒等式与连接系数

原文的推导从下面两个恒等式出发：

转回归一化关联 Legendre 基后，原文式(10)的连接系数为

从左侧逐列作 Givens 消元时，列正交性保证已经处理过的行不会留下额外非零项，最终得到矩形恒等阵。这就是显式旋转分解的来源。上述有限截断统一以 列计数；原文定理中的 尺寸与只写到 的乘积上限存在一项偏差，按列数计数时应包含全部 个旋转。

这一精确分解也提供了可验证的稳定实现。作者分别从两个正量比值计算 与 ，避免通过 求出很小的分量。Algorithm 2.5进一步规定先用整数运算形成分子、分母，再转为浮点数，进行除法与平方根运算；在文中给出的整数范围内，得到每个旋转 、 的界。由正交变换的标准误差分析，未压缩旋转连接及其转置最小二乘具有2范数后向稳定性，无需通过扩展精度来稳定这些步骤。[2, Algorithm 2.5及其后讨论]

不过，仅有相邻阶数的线性成本还不够。若每个阶数都逐级走到0或1，全部球谐层的作用仍需 运算。解析旋转把应保存的参数降到 ，并减少平方根计算，却没有自动把完整变换降到近二次。因此，论文继续压缩这些旋转链所表示的高阶连接矩阵。[2, §2.1, §3.1]

Step 2 用蝶形压缩高阶到低阶的连接矩阵

这里压缩的是连接系数，而非关联 Legendre 函数的节点取值。原文采用的秩性质是：如果连续矩形子块所含条目数不超过 ，那么在指定精度下，该子块可由较小的秩 近似。蝶形算法同时利用不同长宽比例下的这种结构。它从狭长列条带开始，用插值分解保留少量实际列；随后合并相邻条带，再将合并后的矩阵沿行分成上下两半，并对两块中的旧骨架列重新压缩。不断重复后，行范围逐渐缩小、列范围逐渐增大，而每层只传递少量骨架自由度。[2, §3.1]

如果层数为 、每层有 个 插值矩阵，则单个连接矩阵及其转置的作用成本为

通常令 与 同阶，可写成 。实现还利用三类对称性：正负阶数共享结构，偶数阶与奇数阶分别转换，同一连接矩阵再按次数奇偶重排以移除棋盘零。插值矩阵中的恒等子阵也作为置换直接作用。这些安排不会改变主要复杂度，却明显影响实际运算量和存储常数。[2, §2.3, §3.1]

连接矩阵的正交性还要求对蝶形构造作一个调整。完整列彼此正交，直接对这些列进行第一次插值分解不会产生压缩；因此，作者跳过这一层，从更合适的子块开始。后续若某块几乎满秩，其插值矩阵也接近置换，此时应避免把它作为一般稠密矩阵乘法处理。这反映了一个实际原则：决定是否值得压缩的是当前子块的秩与尺寸，而不是整个变换已经被冠以快速算法的名称。[2, §3.1, §3.3]

作者对这一压缩的理论解释来自 Fourier–Bessel 表示。沿用原文该小节的局部下标，偶阶连接系数定义为

这个小节沿用原文式(14)的下标方向：第一个下标 标记高阶源函数，第二个下标 标记零阶目标函数；它与前面式(10)先写目标行、后写源列的顺序不同。因为归一化 Legendre 多项式的 Fourier 变换满足

其中 是球 Bessel 函数，所以原文式(16)将连接系数写为

这个公式将连接算子分成关联 Legendre、Fourier 和球 Bessel 三部分。Fourier 核在有限相空间矩形中的近似秩受矩形面积控制，因此关键是证明只需考虑有限的有效波数范围。小波数处利用 的幂次界；大波数处利用 在端点的消失，以及 Gegenbauer 连接公式带来的 Fourier 衰减界。Theorem 3.5据此给出显式的 与 ，使区间 以外的贡献为 。其中 ；当 与 同尺度时，。小阶数情形另有半可分结构可用。[2, Lemmas 3.3–3.4, Theorem 3.5, Remark 3.6]

Remark: Fourier–Bessel 表示的有效频带与两端尾积分

对于上述偶阶连接，原文 Theorem 3.5给出的边界为

第一个界利用

控制小波数区。第二个界先将高参数 Gegenbauer 多项式连接到适合权函数 的族，再利用球 Bessel 函数的大波数界及端点消失阶数控制尾积分。原文证明分别把 与 的贡献界到 ，所以将两尾合并时严格得到至多 ；若希望总截断误差不超过 ，在两个边界的定义中分别使用 即可。这个常数调整不改变相应的复杂度讨论。[2, Theorem 3.5及其证明, pp.595–597]

这些界针对连接系数的积分表示。它们说明哪些波数对指定精度有贡献，并为随后讨论Fourier核的可压缩性提供依据；它们不是直接作用算法中的数值积分规则，也没有给出全部蝶形因子的浮点误差界。

上述积分是解释压缩结构的分析工具，实际算法不对每个输入计算这一双重积分。文献[2]通过这一 Fourier 积分算子分析给出整个球谐系数转换 的执行复杂度结论。应同时明确，Theorem 3.5直接陈述的是有效波数截断；它不提供多层插值分解在浮点运算中的统一全局后向误差界。论文后面对平均数值秩的统计，以及压缩算法的误差实验，承担的是不同层面的验证作用。[2, 摘要, §3.1.1, §4–5]

Step 3 通过共享的非对角低秩转换得到三角系数

第一阶段之后，所有阶数只需调用两类固定的三角连接矩阵。定义原文式(19)中的函数

低阶关联 Legendre 到三角函数的系数可以由这个 Gamma 比值直接表示：

这分别是原文式(22)、(24)。反向转换也有显式系数公式，见原文式(23)、(25)；余弦反向公式的零阶项及对角项需要用极限或单独的对角公式处理。计算 时，作者使用移位后的渐近展开，并通过

把较小参数移到适合渐近展开的区间，避免分别形成两个可能很大的 Gamma 值再相除。[2, Algorithm 3.8, 式(19)–(25)]

去掉可分离的行、列缩放后，余弦连接的主要核为

它远离主对角线时足够平滑，适合通常的分层非对角低秩表示：近对角块直接保存，远对角块用少量 Chebyshev 插值点近似。原文将 且 的正方形称为分离块；这里用于连接矩阵的次数索引区域，故 。令 ，以 表示在 上的次数 的 Chebyshev 插值多项式。对任意固定的 ，当 时，Theorem 3.10引用 Keiner 的结果，给出

这里 是连接核，普通字母 是误差界中的常数，两者沿用原文的区分。几何衰减意味着，达到固定精度所需的插值阶数不必随矩阵规模增长；更一般地，它只随 增长。[2, Definition 3.9, Theorem 3.10, 式(26)–(27)]

实际实现采用 Chebyshev 重心插值。对一个有 行、 列的远块，在行变量中取 个插值节点，储存 的核样本矩阵 ，再以 的重心插值矩阵 恢复各行，即 。作用时先计算 ，再乘 ，成本为 。分层组织这些块就得到文中的 HODLR/FMM 式快速转换。它利用的是非振荡核的远场平滑性，与第一阶段的互补低秩结构有所不同。两类低阶转换由所有 共享，因而无需为每个阶数重复建立第二阶段计划。[2, §3.2]

预计算、thin plan与数值表现

完整流程可以概括为：预计算相邻阶旋转，用其生成并压缩高阶到0或1阶的连接矩阵，建立共享的低阶到三角系数变换；作用时依次执行这两级转换，需要节点值时再接 FFT 或 DCT/DST。反向计算使用第二级的反向连接以及第一级连接矩阵的转置。这里第二级的逆一般不能直接用转置替代。[2, §3]

预计算仍是本方法的重要开销。原文摘要和结论给出的总预计算量级为 。为改善实际成本，§3.3提出只对间隔约 的阶数建立蝶形，其中 是平均数值秩；其余层先通过若干次相邻阶 Givens 转换到已建立计划的阶数，再调用对应蝶形。这样形成的 thin plan 减少昂贵蝶形的数量，同时利用相邻阶转换的线性成本。文中关于相对完整计划降低约 构造成本的计数，结合了实验支持的 ；它没有把整套预计算降为近二次。[2, §3.3, §5]

这个取舍可以用一组相邻阶数来理解。单层蝶形约需 运算，而把阶数 、 的一组展开都转换到阶数 ，总旋转成本为 。适当选择已计划阶数之间的间隔，就能用较便宜的短旋转链替代一部分昂贵计划，而不使执行成本失去快速算法的优势。间隔过小，节省的预计算有限；间隔过大，许多输入需要经过较长的旋转链。原文实验采用固定的64作为间隔参数，并未声称这一数值对全部规模最优。[2, §3.3–4]

从软件组织看，高阶连接矩阵可以通过正交旋转逐层更新，再交给蝶形算法压缩，不必同时保存全部未经压缩的矩阵。第二阶段只保存所有阶数共享的两个三角连接计划。这样可以降低构造时的临时存储，然而所有已压缩层及其插值因子仍要占据大量内存。论文指出，沿用Tygert的蝶形存储组织，单层可采用 存储，整体为 量级；这类渐近计数仍需结合秩参数和实际字节数判断是否能留在主存。[2, §4]

论文使用 Julia 的 FastTransforms.jl 与 LowRankApprox.jl，在2.8GHz Intel Core i7、16GB内存的机器上测试。输入球谐系数先逐列归一化，经过转换到双 Fourier 系数再转回，统计最大列2范数误差，并对三次试验平均。未压缩 Givens 方案的误差增长经验上约为 ，thin 蝶形方案约为 。这说明压缩带来额外误差，但不能把这些随机输入下的趋势当成最坏情形误差定理。[2, §4, 图4–5]

还应把误差来源分开理解。正交连接使原始变阶问题良条件，解析旋转控制了生成和直接应用这些连接的舍入误差；蝶形插值则另外引入可调的近似误差，并伴随多层矩阵乘法的舍入。原文对插值分解的稳定性援引了经验依据，没有将其全部推导为与Givens相同的全局后向稳定结论。因此，该文标题中的稳定性贡献有明确的数学核心，快速版本的整体精度则需要结合压缩容差与数值实验评价。[2, §2.1–2.2, §4–5]

实际存储也限制了渐近优势的显现： 时 thin plan 已需约26GB，超过物理内存后，执行时间混入了磁盘换页开销。作者明确指出，测试范围中的执行趋势更接近 ，尚未呈现 的最终渐近行为。因而，这篇工作的核心贡献是把良条件的正交连接、解析可得的稳定旋转和快速压缩组织到同一个基转换框架中，并显著改善预计算的实际可用性；低预计算、低存储以及完整快速算法的统一舍入误差控制，仍应分别评价。[2, §4–5, 图5–6]

7. Bremer–Chen–Yang（2021）：转折区域分块与 IDBF

Rapid application of the spherical harmonic transform via interpolative decomposition butterfly factorization, by James Bremer, Ze Chen, and Haizhao Yang. [1]

Bremer、Chen与Yang关注的是另一项直接影响实际使用的问题：已有蝶形方法即使能够很快地作用变换矩阵，往往仍需先生成或扫描接近完整的稠密矩阵。对于每个角向阶数，这意味着二次预计算；全部球谐层合起来，预计算仍为立方量级。文献[1]的目标因此同时包括矩阵构造和单次作用。其方法依据关联 Legendre 函数的振荡、非振荡及转折结构划分矩阵，并结合按需条目评价，避免把“快速作用”建立在过于昂贵的稠密预处理上。[1, §1]

变换定义、求积结构与奇偶分解

原文的规模参数是 ，最大球谐次数为 ，这一点与直接以最大次数为参数的记法不同。其展开为

其中 表示固定 下的 归一化关联 Legendre 函数。纬向采用 点 Gauss–Legendre 求积，经向采用 个均匀节点

对每个 ，先计算

再沿经度用 FFT 形成 。经向工作的总成本是 ，真正需要加速的是固定 的关联 Legendre 变换，原文简称 ALT。本文把“系数到节点值”称为 forward，把“节点值到系数”称为 inverse。[1, 式(1)–(6), §3.1]

设 是固定 的 评价矩阵， 是 Gauss–Legendre 权重构成的对角矩阵，则其离散正交关系为

其中 表示反向 ALT。这里即使 为奇数，求积仍然精确：两个同阶关联 Legendre 函数的乘积包含因子 ，整体是次数不超过 的多项式，而 点 Gauss–Legendre 求积能够精确积分到 次。因此，一旦能够快速作用 ，只需把同一计算图转置，再结合权重，就能够计算反向变换，无须另建一个无关的快速逆算子。[1, §3.1.1, 式(25)]

另外，，而 Gauss 节点关于零点对称。作者据此把 为奇数、偶数的系数分别处理，只在半数节点上计算，再通过符号对称恢复另一半。所需的两个矩阵分别为

这种分解在进入压缩之前就减少了重复计算，并使转折结构能够在一个半球上清楚地表达。[1, §3.1.2, 式(26)–(29)]

Step 1 由转折曲线划分振荡区与非振荡区

对于 ，Legendre 评价矩阵具有适合标准蝶形算法的振荡结构；高阶情形则不同。为了理解其区别，可以把关联 Legendre 方程写成标准形。令

则由关联 Legendre 微分方程得到

当 时，方括号中的系数在

处变号。这就是原文式(30)使用的转折点。靠近极点的一侧对应非振荡、快速衰减的行为，靠近赤道的一侧对应振荡行为。由于 随次数 改变，它在“节点—次数”矩阵中形成一条曲线，将矩阵分成性质不同的区域。[1, §1, §3.2]

区别并不只是某些条目较小。非振荡区域中的一个大块可以具有通常的低数值秩；振荡区域中的大块一般不是低秩，但在适当的长宽配对下，其子块可能具有互补低秩；穿过转折曲线的块则同时包含两种行为，不适合直接套用纯振荡区域的压缩模型。因此，本文先依据曲线划分计算对象，再对不同区域采用相应的压缩工具。上述 应理解为标准形方程的转折点，而不是未加权 自身的精确拐点。[1, §3.2–3.3]

具体的层次分块从奇、偶 ALT 矩阵的近方形初始块开始。完全处于转折曲线一侧的块停止细分；与曲线相交的块继续作 二分，直到其行数或列数小于给定阈值 。这种过程只在转折曲线附近持续细分，而不会把整个矩阵统一分到最细层。按原文的计数，第 层有 个相关块，其边长为 ，层数为 ；划分本身只需线性量级的工作。[1, §3.2, 图2]

Step 2 按需评价条目并构造插值分解

本文使用的基本工具是插值分解，简称 ID。对一个近似秩为 的矩阵 ，列 ID 选择一组实际列的索引 ，使

这里 是骨架列， 是将所有列表示为骨架列组合的插值矩阵，并在相应位置含有恒等子阵。其意义是：后续计算只需保留少量有代表性的原始列及它们的插值关系，而不必保存全部矩阵元素。[1, §2.1, 式(9)]

通常对完整矩阵做秩揭示分解仍然昂贵。作者采用抽样版本：先在行方向取少量代表点，只评价这些行与全部列构成的子矩阵，再用带列主元的 QR 选出骨架列，解小型三角系统构造插值关系。代表行既可随机选取，也可取离散网格上接近 Chebyshev 点的 mock-Chebyshev 点。若精度不足，则增加用于采样和压缩的秩参数。对于 的块，当孤立条目可以按需快速评价时，列 ID 的构造工作约为 ，存储为 ；行 ID 则相应交换两维。[1, §2.1]

这种方法的成本优势来自“只评价需要的条目”。若先把整个 矩阵全部生成，再进行同样的 QR，其存储与生成成本并不会消失。同样，少量样本决定的是一个压缩模型；局部 QR 的截断容差不应直接等同于整套 ALT 对所有输入的最终误差上界。文中将采样选择、秩参数和最终作用误差分别进行数值检验，正是因为这几个因素共同决定压缩质量。[1, §2.1, §5]

Step 3 对三类矩阵块采用相应的快速方法

振荡块采用 IDBF。 对一个 的振荡块 ，分别对行、列索引建立深度为 的二叉树。互补低秩性质指：行树第 层的一个节点，与列树第 层的一个节点组成的子矩阵，能够以低秩近似。随着行区间缩小，允许的列区间相应扩大；这类配对覆盖了蝶形计算所需的不同长宽比例。[1, §2.2]

IDBF对这些子块实施行、列插值分解，将 近似表示为一串稀疏因子：

中间的 连接两侧的骨架自由度，外围各层负责在不同分辨率之间转换。单次作用就是依次乘这些稀疏矩阵；转置作用则按相反次序应用转置因子。在固定块秩较小时，每层保存和传递的数据远少于原始稠密块。原文对这一分解的构造、作用和存储均给出 的界，其中 是相关数值秩的上界。[1, 式(17), §3.3, 式(34)]

非振荡块采用抽样低秩分解。 在这一侧，函数远离转折点后快速衰减。作者先缩小块的有效范围，排除机器精度以下的可忽略部分；完全可忽略的块不再保留。对剩余部分，Algorithm 1通过少量行、列抽样构造一个近似 SVD。[1, §2.3, §3.3]

它首先在抽到的行、列上作带主元 QR，选出代表列和代表行，再对相应的完整列、完整行正交化，得到 与 。随后只读取一个小交叉子矩阵 ，拟合核心矩阵

对 做小规模截断 SVD，再乘回两侧正交基，即得到 。因而，本方法确实包含 SVD，但主要的 SVD 运算发生在低维核心上。采用固定过采样参数时，一个边长 、秩上界为 的非振荡块，其构造成本为 ，单次作用和存储为 。[1, Algorithm 1, 式(35)]

采样数量通常略大于最终保留秩，使选中的行、列更可能覆盖主要数值子空间。原文用过采样参数 控制这一数量，在其试验中 已有较好的表现。这是抽样质量的经验选择，并不是任意低秩矩阵都能由同样少的随机行列可靠恢复的无条件保证。相应地，振荡区ID中的预设秩参数、非振荡区SVD的截断秩，以及理论中的 含义不同，不能把它们作为同一个固定常数使用。[1, §2.1, §2.3, §5]

转折附近的小块直接作用。 与转折曲线相交的块经过反复细分后已经足够小，作者对它们保存所需条目并直接做矩阵向量乘法，同时也裁掉可忽略的部分。它们只占转折曲线附近的窄区域，并不意味着算法在高阶情形退回整个稠密 ALT。按原文式(33)，这些块的总工作与存储有 的上界，低于主导的分层振荡块成本。[1, §3.2–3.3]

特殊函数条目评价为何是预计算复杂度的前提

上述压缩算法不断要求评价一些不连续的、由抽样和骨架选择决定的 。普通的次数三项递推善于按顺序生成一整列函数值，但若为了取得某个孤立条目，每次都从低次数递推到高次数，单条目成本便会随 增长。这样即使最终只保存少量条目，压缩构造的总工作仍可能过大。

为此，文献[1]明确依赖 Bremer 此前提出的关联 Legendre 数值评价算法，提供运行时间不随 degree、order 增长的条目评价接口。这个接口与 IDBF 的抽样构造相配合，才使得“需要多少条目”能够转化为“花费多少构造工作”的近线性计数。本文并未自行重新推导这一特殊函数算法；综述其成本时，应将这种按需 条目访问列为计算模型中的必要条件。[1, §1]

这一点也说明，与已有矩阵压缩方法的区别不仅在最终保存了多少个因子。构造阶段能否快速取得选中的条目，会直接决定预计算是否仍然需要二次甚至更高的工作。本文把层次分块、IDBF和特殊函数求值接口一起设计，才避免了先完整生成每个 ALT 的步骤。[1, §1, §3]

建立计划时，每个固定阶数所需的工作依次是确定奇偶列、生成转折分区、裁剪可忽略区域，并为各类块保存相应因子或直接条目。输入系数到来之后，这些几何判断和矩阵分解都不再重复：只需提取对应子向量，应用块因子并累加结果，再恢复另半球和完成经向FFT。反向作用沿同一块划分作转置累加，并接入求积权重。计划与输入无关，因此重复变换可以摊销构造成本；本文的复杂度分析进一步要求，即使把首次构造算进去，也保持接近线性的固定阶成本。[1, §3–4]

成本的已知部分与秩增长猜想

原文用 同时界住非振荡块的普通数值秩，以及振荡块在互补分块下的数值秩。在第 层，振荡块数量为 ，每块规模约为 ，所以对所有层求和得到

非振荡块的分解更便宜，转折叶块也不改变这个主导项。因此，Theorem 4.1给出固定精度下，一个奇或偶 ALT 的时间、存储复杂度均为

其中时间包括分解构造和单次作用，而不是只包括已有计划之后的乘法。[1, §4, Theorem 4.1]

尚未闭合的是 的一致增长界。Conjecture 4.2根据数值实验猜想 ；在这一猜想下，单个 ALT 的成本为 ，全部角向阶数加上 FFT 后，整个球谐变换为 。论文把这一结果称为接近渐近最优，并明确说明需要进一步证明相应的块秩界。因此，对该工作的准确评价应同时保留其含预计算的近线性目标、已经给出的条件复杂度定理，以及仍由实验支持的秩增长部分。[1, Conjecture 4.2, §6]

数值证据与稳定性的实际含义

论文的 MATLAB 代码在单个3.2GHz CPU核上运行，测试规模为 ，阶数取 。ID 容差为 ，转折块阈值 ，非振荡低秩分解的秩参数为30。对ID采样的比较显示，mock-Chebyshev 点比纯随机点更准确；后续试验将其自适应秩参数设为150。这些参数说明论文主要检验的是给定容差下的低常数快速作用，而不是采用机器精度容差后对所有参数给出统一保证。[1, §5, 图8]

其误差定义也需要与上述设置一起理解。前向测试在非零区域随机取256个输出坐标，比较快速结果 与直接矩阵乘法结果 ：

反向测试类似地抽取256个系数坐标。图9中非零阶数的误差大致处于 至 ，随规模略有增长； 的特殊情形接近机器精度。这样的测试支持指定样本和容差下的数值稳定性，但不是全矩阵谱范数误差，也不直接等价于对所有输入的统一误差界。[1, 式(36)–(37), 图9]

图10将构造和作用时间与 、 参考曲线比较，图11进一步与显式矩阵生成和直接矩阵向量乘法对照，表明大规模下有明显节省。不过，这里的主要实验对象是固定 的 ALT，不能据此直接给出完整 SHT 与其他高度优化软件之间的统一性能排序。类似地，离散正交性保证精确算子的加权转置关系，局部压缩与浮点舍入仍会引入误差；本文没有建立覆盖全部抽样、截断和多层因子作用的全局后向稳定定理。[1, §5–6]

这篇工作的贡献在于把转折几何纳入压缩结构，并使矩阵条目生成、计划构造与单次作用的复杂度相互一致。它展示了直接节点变换如何同时降低 setup 和 apply 的可能性，也清楚地留下了进一步研究的重点：证明随阶数和规模变化的块秩界，改进代表点选择及低秩分解，并把局部近似容差与完整变换的误差联系起来。[1, §6]

两条路线的对比

六篇工作可以用两个计算流程来概括。节点值路线完成

系数连接路线完成

节点值路线直接贴近“给定球谐系数，在求积网格上求值；给定这些值，求回球谐系数”的需求。在 Gauss–Legendre 网格上，分析公式由求积正交性直接给出，数据和球面积分也自然匹配。代价是某些中间节点会随 改变：Tygert II/III 必须额外处理专用节点到公共节点的插值。Bremer–Chen–Yang 直接压缩公共节点上的矩阵，省掉了这层转换，却必须面对转折曲线穿过采样矩阵时的非均匀结构。[4, §3.2–4；5, §5；1, §3.2]

系数连接路线把“如何表示函数”和“在哪些点求值”分成两个阶段。当后续工作本来就需要三角系数，或希望在等角网格上反复使用 FFT 时，这种中间表示很有价值；转换得到的 Fourier 系数自身也是可复用的数据。它同时需要承担高阶到低阶连接、低阶到三角系数连接两层工作。球面函数的极点正则性和奇偶约束也会限制所得 Fourier 系数，并非任意一个双变量 Fourier 系数数组都对应同一带限下的光滑球谐展开。[3, §3–5；2, §1、§3]

因此，“需要节点值”与“需要 Fourier 系数”是选择路线时比“是否使用蝶形”更靠前的问题。两条路线都能最终求值；这里比较的是哪一种中间表示最符合后续计算，而不是宣称某一种表示只能用于某一类网格。

其次，相同工具在不同路线中承担不同任务：

• 特征向量方法的区别。 Rokhlin–Tygert 对角化的是写在低阶基底中的 Sturm–Liouville 算子。其特征向量回答“高阶基函数如何在低阶基底中展开”。Tygert II 对角化的是同奇偶基底中的 乘法算子。其特征向量回答“基函数在相应求积节点上取什么值”。前者输出连接系数，后者输出加权采样值；这也是 [3] 与 [4] 最重要的结构区别。
• FMM 的区别。 Suda–Takami 的 FMM 加速分裂分量的 Lagrange 插值；Tygert II 的 FMM 既服务于快速特征向量算法中的 Cauchy 型乘法，又服务于 Christoffel–Darboux 插值；Rokhlin–Tygert 的低阶到三角连接也利用可快速作用的结构化核。FMM 在这里主要解决的是“如何快速求和”，它本身不自动证明外层递推、节点选择或整个变换的稳定性。[6, §2.2；4, §2.1–2.2、§3.2；3, §2.3–2.4]
• 蝶形分解的区别。 Tygert III 压缩加权节点采样矩阵；Slevinsky 压缩跨阶正交连接矩阵；Bremer–Chen–Yang 只对自适应分块之后的振荡区域使用 IDBF，另外两类区域采用其他处理。三者的低秩结构发生在不同的矩阵中，不能把一篇论文对其矩阵的理论解释直接移植成另一篇的统一秩定理。[5, §3–5；2, §3.1；1, §3.2–4]

这也解释了一个容易造成困惑的现象：Slevinsky 给出了 Fourier 积分算子与截断频带的理论分析，Bremer–Chen–Yang 仍明确保留秩增长猜想。二者不构成直接矛盾。除了矩阵不同，还需要分别检查所证明的频带结论、互补分块方式、对所有阶数的统一性，以及由局部近似误差得到整个压缩算子误差的步骤。[2, Theorem 3.5、Remark 3.6；1, §4]

另外，复杂度需要把预计算与作用分开。下表用 表示统一的最高球谐次数尺度，隐藏常数精度要求下的部分因子；各文的原始 并不相等。“作用”指计划已建立之后的一次变换。表中的全球复杂度由所有阶数的纬向变换加上经度 FFT 构成，不能与某一张单阶 ALT 实验表直接混用。

文献	预计算的关键结论	计划建立后的作用复杂度	应当保留的条件
Suda–Takami [6]	论文给出全部阶数 预计算	固定插值精度下	复杂度分析不等于完整递归误差界；实际稳定选点是贪心近似
Rokhlin–Tygert [3]	实现为每阶二次、合计 ；另指出更快理论预计算可用		利用结构化特征向量与低阶连接算法；实验目标约六位精度
Tygert II [4]	理论快速谱预计算与实际实现分开；实际每阶 ，合计三次量级	固定精度参数下	特征向量误差与谱间隔有关；实际计划采用四倍精度
Tygert III [5]	每个奇偶模态压缩仍为 ；深度优先主要降低峰值内存	条件性近于 ，须保留秩因子	单个 奇偶模态矩阵为 ；有效秩的良好行为主要由实验支持
Slevinsky [2]	总预计算仍为 量级；thin plan 减少需要建立和保存的连接	作者给出	理论解释来自 FIO 与分层压缩；不应额外附加原文没有证明的统一全局误差结论
Bremer–Chen–Yang [1]	预计算、作用、存储一起获得秩相关上界	全阶为 ；秩猜想下为	必须有固定精度的 条目评价，并将有效秩控制到所有相关阶数及子块

对于最后一行， 表示比较时采用的相关有效秩上界；固定阶数的原定理先给出 。把它推广为全阶表达式时，需要对参与的阶数一并控制秩。论文给出的 是猜想，不能省略这一步就称整个算法已经具有无条件的 保证。[1, Theorem 4.1、Conjecture 4.2]

这个表格说明，早期几篇工作的主要成果已经是快速作用，而 Bremer–Chen–Yang 特别推进了快速建立变换。Slevinsky 的重点则是可解释、可准确生成的系数连接和到 Fourier 表示的完整链条。只按作用复杂度中对数的次数排序，会忽略实际应用中很重要的预计算时间、内存以及目标精度。

对于数值稳定性而言，精确算子的正交性只是起点。如果高阶与低阶基底在同一个 内积下都正交归一，系数连接本身就保持二范数；如果节点采样先乘上正确的平方根求积权重，节点路线也有同样的等距结构。因此不能笼统地说“系数路线稳定，节点路线不稳定”。两条路线都能从良态的精确算子出发，真正的区别在于实现是否充分利用了这个结构，以及近似产生在哪里。例如 Tygert III 的精确加权采样矩阵是正交矩阵，但把它近似成许多 ID 因子的乘积后，局部因子一般不再正交。Slevinsky 的解析 Givens 因子则直接保持旋转结构，未压缩连接有明确的后向稳定分析；其快速算法还包含 ID 蝶形与低秩连接，因而需要区分这些后续步骤与解析旋转核心的保证。[5, Lemma 3.7、§5；2, §2.1–2.2、Algorithm 2.5]

Remark: 为什么每一块误差小，不直接等于整个变换误差小？

设一个变换表示为 ，计算时使用 。一个精确的望远镜展开为

由此可见，误差不仅取决于 ，还取决于其前后因子乘积的范数。若各因子都接近正交，这个累积较容易控制；若只是完整乘积在精确算术下正交，中间因子却含有较大的插值系数或缩放，就需要额外分析。这个基本事实解释了为什么 ID 条目界、单块逼近容差、完整矩阵正交性和整个程序后向稳定性是不同层次的命题。

同样，节点值与加权节点值之间还隔着 或 。系数的二范数对应球面 范数，加权节点二范数也有明确的求积意义；未经加权的最大分量误差则是另一种量。报告“若干位准确数字”时，必须同时说明比较对象和误差范数。

各篇的证据层次可以更具体地归纳如下。

文献	可直接依赖的稳定性结构	不宜扩大解释的地方
[6]	顺稳定方向使用递推；最大体积选点提供插值系数控制的存在性依据	实际 max-min 贪心不保证全局最大体积；原文没有完整算法的解析误差界
[3]	同一内积下的正交连接；快速结构化谱算法	六位左右精度的实验不能当作双精度机器精度验证
[4]	对称三对角谱分解；加权采样的正交性；求积与 CD 恒等式	谱间隔和容差进入误差分析；实验的四倍精度预处理不能省略
[5]	精确加权矩阵正交；一个压缩图可转置用于反向计算	ID 压缩与完整流程的误差主要由数值实验支持，前向误差不能代替往返误差
[2]	相邻阶解析 Givens 旋转及其生成算法具有明确的后向稳定分析	该结论不能原封不动地覆盖全部蝶形近似、三角连接和最后 FFT
[1]	精确离散正交性 ；反向使用带权转置，压缩精度由实验检验	抽样输出误差、局部容差和秩猜想不构成完整 SHT 的统一后向稳定定理

这里对结论范围的限制并不否定这些算法的数值价值，而是为了准确区分六篇论文中的精确数学性质、快速近似步骤和实验验证结果。

需要注意的是，六篇论文的实验跨越二十年，硬件、编程语言、精度要求、对照实现和测试对象都不同。Suda–Takami 把少量阶数的测试按比例汇总来估计整体纬向工作；Rokhlin–Tygert 使用稠密矩阵乘法的代理基准；Tygert II/III 主要测试单个奇偶模态；Slevinsky 测试球谐系数与 Fourier 系数的往返；Bremer–Chen–Yang 的主要图也是固定阶数 ALT，并在随机抽取的坐标上计算相对误差。它们并不是在共同平台上对同一份全球数据做同一精度的测试。[6, §2.6；3, §6；4, §4；5, §6；2, §4；1, §5]

从这些实验可以提取“某个结构能否在大规模下起作用”“预计算是否成为瓶颈”“秩随规模是否增长缓慢”等证据，但不能据此得出某篇算法在今天的所有设置中最快。尤其需要避免把相对均方根误差、归一化输入的最大分量误差、系数往返二范数误差和抽样相对误差放在同一列中直接比较数字大小。

如果要重新比较，至少应固定同一带限、同一采样或输出表示、同一误差门槛和同一精度环境，分别测量计划时间、单次合成、单次分析、峰值内存及多次执行后的总时间。这样的比较是由现有论文揭示的必要评估口径，不是这些论文已经完成的共同实验。

预计算耗时是整个算法不能忽视的一部分，应用次数决定昂贵预计算能否被摊薄。设建立计划的时间为 ，一次快速作用为 ，一个具有相同输入输出和精度要求的直接方法为 。做 次变换时，快速方法的总时间是

如果 ，它开始节省总时间的条件为

上式假设直接法没有另计的准备成本；若两种方法都需要准备，应把分子换成两者准备时间之差，并按同一计时范围比较。这个简单关系说明了为何“长期固定带限反复使用”和“每次改变带限或阶数”会给出不同的选择。前一种情况容许花较长时间建立一个紧凑而快速的计划；后一种情况对计划成本敏感，Bremer–Chen–Yang 同时加速构造与作用的设计就更值得关注。还需要考虑内存：一个理论上能摊薄时间的计划，如果大到触发换页，实际作用时间也会改变。Slevinsky 的大规模实验明确展示了这种影响。[2, §4]

总的来说，这两条路线各有优劣：

• 节点值路线的长处是与求积和物理网格衔接直接。 在公共 Gauss 网格上，合成、带权转置和球面积分之间的关系清楚；它也允许直接利用关联 Legendre 函数在物理空间中的递推、衰减区和振荡区。Suda–Takami 的分裂插值和 Bremer–Chen–Yang 的转折分块，都是沿着这种函数值结构开展的。其困难是不同阶数的行为不均匀，端点附近有很大的尺度变化，转折区域和专用节点插值都需要专门处理。Tygert II/III 还说明，得到一个漂亮的正交采样矩阵，并不自动解决节点生成、权重生成和低成本压缩计划的问题。
• 系数连接路线的长处是把正交基之间的代数关系显式利用起来。 同一个函数在两个正交基中的系数变化具有清楚的范数解释，Slevinsky 的相邻阶旋转更使这种关系成为可准确生成的局部操作。低阶到三角基底的转换又有共享的 Gamma 比值核，可以跨多个 复用。其困难是高阶到公共低阶的连接仍然是大规模工作，直接串接旋转会回到三次复杂度；即使用蝶形加速，生成和保存许多连接也可能代价很大。此外，低阶关联 Legendre 基与正余弦基对应的自然积分权不同，最后这一步不能仅靠“正交连接取转置”一语带过。
• 从这六篇论文中得不到一条在所有指标上占优的路线。 若终点就是 Gauss 网格上的分析与合成，节点值路线与该数据接口衔接直接，值得优先评估，实际速度仍取决于具体实现；若 Fourier 系数本身有后续用途，系数连接路线的中间产物便有独立价值。若带限经常改变，计划复杂度更重要；若同一计划执行很多次，则作用时间和存储更重要；若需要严格追踪舍入误差，解析正交因子提供了很有价值的起点，但仍应把后续压缩误差单独纳入分析。这些判断来自计算流程本身，不依赖跨年代性能数字的直接拼接。

共同技术特征与主要贡献

这组六篇论文并不是沿着一条直线不断替换旧算法。它们反复表明，一个快速变换的构造通常需要同时解决三个层次的问题。

第一层是选择保留结构的表示。 把高次数函数拆成两支低相对次数分量、把微分算子写入低阶基底、把乘法算子写成三对角矩阵、或者把相邻阶连接写成平面旋转，都在减少“无结构的稠密矩阵”这一表面描述带来的困难。真正的进展往往发生在压缩之前。

第二层是找到合适的快速作用形式。 对 Cauchy 核使用 FMM，对准确的半可分结构使用快速谱算法，对互补尺度上的近似低秩使用蝶形，对远离对角线的光滑核使用分层低秩。不同结构需要不同的压缩组织方式；矩阵整体低秩、普通非对角低秩和互补低秩并不是同一个概念。

第三层是让表示和压缩本身可以低成本、稳定地建立。 这包括稳定递推方向、节点和权重生成、Givens 参数的准确计算、深度优先的内存调度、只保存部分连接的 thin plan，以及可以随机访问的特殊函数条目评价。忽略这一层，会得到一个每次作用很快、却难以实际使用的算法。

因此，这些工作的核心贡献可以分别概括为：Suda–Takami 建立了稳定分裂递推与快速插值的分治方法；Rokhlin–Tygert 把跨阶连接识别为结构化谱问题；Tygert II 把本阶节点采样与三对角特征向量联系起来；Tygert III 用蝶形与内存调度改善这一采样框架；Slevinsky 从解析正交连接出发构造到 Fourier 表示的完整快速链条；Bremer–Chen–Yang 利用区域结构和随机条目访问进一步处理预计算瓶颈。

综合来看，六篇论文分别利用了递推分裂、算子谱结构、正交连接、快速插值和蝶形压缩等工具，但其预计算成本、作用复杂度与稳定性依据并不相同。这些结论仅对应各篇所研究的球谐或关联 Legendre 变换，不能在缺少额外分析的情况下视为其他函数族或其他变换的定理。

Remark: 阅读与公式核对说明

本文以所附六篇期刊 PDF 为依据。参考文献中的原文节号、公式号和定理号用于区分“作者明确证明”“在结构或秩假设下成立”“实验观察”和“本文根据流程作出的比较”。没有把某篇论文的单阶实验时间扩展为实际测得的完整 SHT 时间，也没有将局部稳定性推广为整个近似程序的统一后向稳定性。

数学整理中特别保留了以下更正：Tygert II 的求积公式按 Tygert III Remark 3.8 补回权函数；Slevinsky 的有限 Givens 乘积用列数 消除截断下标不一致；关联 Legendre 离散正交性由同阶乘积的多项式性质解释；Bremer–Chen–Yang 的转折点按标准形微分方程理解。需要变更或增补原记号的地方均在对应正文中说明，不以改换符号掩盖差异。

文稿还对相邻阶旋转连接、半可分生成元、三对角特征向量与采样值关系、正反 CD 插值，以及低阶到正余弦的显式公式做了独立的小规模数值核对。这些检查用于排除符号、归一化、转置和系数抄写错误，不构成对六个完整快速算法的重新实现，也不替代原论文的理论证明。