Butterfly Algorithm | ALL YOU NEED IS LOVE

REFERENCES

A Fast Butterfly Algorithm for the Computation of Fourier Integral Operators, by Emmanuel Candès, Laurent Demanet, and Lexing Ying. Multiscale Modeling & Simulation, 7(4), 1727–1750, 2009.

Butterfly Factorization, by Yingzhou Li, Haizhao Yang, Eileen R. Martin, Kenneth L. Ho, and Lexing Ying. Multiscale Modeling & Simulation, 13(2), 714–732, 2015.

Interpolative Decomposition Butterfly Factorization, by Qiyuan Pang, Kenneth L. Ho, and Haizhao Yang. SIAM Journal on Scientific Computing, 42(2), A1097–A1115, 2020.

许多积分变换和特殊函数变换在离散后，都要求计算一个稠密矩阵与向量的乘积

如果输入和输出各有个数，直接计算需要次运算。Fourier 变换可以借助 FFT 降到；但对相位更复杂的振荡核，或者由特殊函数组成的矩阵，通常没有相同的精确分解。另一方面，这类矩阵往往整体满秩，也不能用一个小规模低秩分解代替。

蝶形算法利用的是另一种结构：目标区间越小，允许与它配对的源区间就越大；在这些尺度互补的区间对上，矩阵块可以保持低数值秩。 算法逐层改变两个区间的尺度，用少量局部系数传递整个求和。它所压缩的是不同尺度上的相互作用，而不是整个矩阵的全部列空间。[1, §2；2, §1.1]

本文围绕上述三篇论文介绍这条思路及其发展。Candès–Demanet–Ying 从 Fourier 积分算子的相位出发，证明局部低秩并给出显式插值算法；Li 等人把互补低秩矩阵组织成一串稀疏矩阵，并提供两种黑箱构造方法；Pang–Ho–Yang 则用分层插值分解，让矩阵条目可以直接访问时的构造成本也接近线性。

这三篇工作的差别，主要在于局部表示怎样获得、层间信息怎样复用，以及整个表示怎样低成本地建立。下面先解释它们共有的计算结构，再分别展开具体算法。文中的“原式”、节号和定理号均指对应参考文献；各篇的规模符号不强行统一。这里讨论的进展以所选文献中的 2020 年 IDBF 为终点。

1. 蝶形算法基本思想

整体低秩与互补低秩

普通低秩分解把矩阵写成

这样可以先计算，再乘，成本为。但以 Fourier 矩阵为例，经过归一化后它是酉矩阵，所有奇异值都相同；若要求很小的相对算子误差，就不能通过保留少量奇异向量来压缩整个矩阵。

蝶形算法先分别为行索引和列索引建立二叉树。本节采用根层为、叶层为的编号，并暂取、每个叶子含一个索引。如果位于行树的第层，位于列树的第层，那么对应矩阵块的大小为

互补低秩性质要求：这些块在给定精度下都有低秩近似，而且所需秩随规模增长得很慢，理想情况下与无关。这里“互补”指两个树层号之和固定，而不是行、列索引集合互为补集。[2, §1.1]

以为例：

行树层号	列树层号	每个矩阵块的大小	这一层的块数

每层都覆盖整个矩阵，但块的形状从高窄变成方形，再变成矮宽。不能只检查中间层的方块低秩；快速递归需要相邻层上相应的矩阵块也具有低秩结构。

层次分块示意图，其中的每一块都是低秩的。[3, figure 1]

这与通常依靠几何远离来压缩非对角块的分层矩阵方法有所不同。振荡核即使在远离对角线的位置，也可能需要很高的秩；蝶形算法通过一个区间缩小、另一个区间扩大的尺度配对，控制局部振荡自由度。

一个 Fourier 核的例子

考虑

并令、。注意到

于是有精确恒等式

前三个因子各自至多依赖一个变量。真正把与耦合在一起的是最后的残差相位。如果两个区间的宽度和（即与可取的最大值）满足

那么的变化范围有界，不随全局频率范围增大，于是残差相位可以用少量项逼近

每一项都是“只依赖的函数”乘“只依赖的函数”。这就给出了局部低秩近似。全局核可以剧烈振荡，但去掉两端可以分别处理的相位后，局部剩余的耦合并不复杂。这一计算是理解 [1] 相位分解的最简单例子；实际算法使用更稳定、便于递归的 Chebyshev 插值，而不直接保存高次单项式系数。[1, §§1.1,3]

低秩块之间的信息传递

根据上面的分析，当满足低秩互补性质时，在对应的块上可以被低秩逼近，现在需要将这些低秩块上的信息快速地组合以得到最终的结果。设某个互补区间对上有

由内源点产生的部分和可以写成

因此，对当前目标区间而言，整组源点只需用个系数表示。若每次都重新对内所有源点求和，这仍然很贵。蝶形算法真正省下的工作来自下面的递归。

令是的父节点，是的两个子节点。上一层已经知道、的系数，而

下一层只需把这两组系数转换成适合的一组系数：

是小规模传递矩阵；这个统一记号只用于说明，不替代三篇论文各自的公式。空间目标越来越局部，源点覆盖范围越来越大，直至每个目标点收集全部源点的贡献。在二维四叉树中，每次相应地合并四个源子节点。[1, §2]

这一过程解释了两个关键事实：

每层的状态总量保持同一数量级。 行树有个节点，列树有个节点，所以互补节点对共有个。
层数只有对数级。 若每个节点对用个系数、每次用稠密小矩阵传递，则一层需要工作，共需。

“蝶形”指这种分裂、合并、再分裂的交叉连接组织。它与 FFT 都具有分层计算图，但一般蝶形算法中的小矩阵来自局部近似，不必是 FFT 中的精确二点运算，也不要求存在一个可以直接调用的 FFT。

最后还剩一个问题：这些小矩阵本身怎样得到？ 如果为了寻找每个低秩块的基底而遍历全部条目，预计算仍可能需要二次工作。下面三篇论文，正是从不同角度回答这个问题。

符号约定

文献	规模符号的含义	树层编号	局部表示
Candès–Demanet–Ying [1]	是二维网格每一方向的点数，输入总量为	根层，叶层	个 Chebyshev 点，
Li 等人 [2]	为矩阵，输入总量为	根层，叶层，中层	数值秩的 SVD 基底
Pang–Ho–Yang [3]	，为初始叶宽	根层，叶层	秩参数、实际骨架秩的 ID

特别是，[1] 的与 [2,3] 的都是相对于输入总量的近线性复杂度，不能仅凭的幂次判断它们谁更快。

2. Candès–Demanet–Ying（2009）：从残差相位得到显式蝶形算法

A Fast Butterfly Algorithm for the Computation of Fourier Integral Operators. [1]

主要问题与贡献

本文考虑二维离散 Fourier 积分算子 (Fourier Integral Operators, FIO)

其中

在时光滑，并且关于频率具有一次齐次性：

直接计算需要次运算。作者给出的算法在固定精度下需要次运算。[1, 原式(1.1)，§1]

这篇论文的主要贡献是为一般齐次 FIO 核找到了可分析、可直接计算的局部表示：先剥离主振荡，只对残差核做插值；由此得到显式的局部基函数和传递公式，避免逐块进行昂贵的代数低秩分解。共享的张量积插值还降低了运算常数与存储需求。

Step 1: 用极坐标消除频率原点带来的相位困难

对于，原文引入

这里最后的指数表示平面单位方向向量。映射后的点位于。对零频率采用的连续延拓，并指定；所有映射后的点共同组成集合。再定义

一次齐次性将径向依赖变成显式的线性因子，角向函数只在单位圆上求值。因此可以在新的参数域上光滑表示，原问题变为

这里仅对离散源点重新编号，没有把数据插值到均匀极坐标网格，也没有引入连续积分换元的 Jacobian。也可将已包含在此式中的零频贡献单独计入。两棵四叉树分别建立在与上；中的点通常并不均匀。[1, §1.1，原式(1.2)，§3]

Step 2: 分离主相位，只逼近剩余耦合

取目标方盒、源方盒，中心分别为，边长分别为。简写两个中心为，定义残差

于是有

残差中，单独随或变化的部分被消掉，只留下混合变化。由混合导数有界可得

因此，只要采用互补尺度

就有。全局相位的振荡次数随增大，但这个局部残差的相位范围保持有界。[1, 原式(1.3)–(1.5)]

严格的精度依赖还需要更强的正则性。原文 Theorem 3.1 假设关于联合变量实解析，并具有统一的解析导数界；对充分大的和充分小的，残差核具有

项的分离近似。这个秩界与无关，但依赖相位的解析性常数。不能仅由就断言任意光滑核都有相同的对数精度依赖。[1, §3.1]

Step 3: 在较小的一侧进行 Chebyshev 插值

在每个坐标方向上取个 Chebyshev–Lobatto 点。原文使用的参考节点为

将二维张量积节点平移、缩放到方盒，得到；对应的二维 Lagrange 基函数记为。目标盒内类似地定义。下标实际代表一对一维节点下标，总节点数为

因为，至少有一侧的边长不超过。原文 Theorem 3.3 说明：选择

时，在较小一侧对残差核做张量积 Chebyshev 插值，可以获得所需局部精度。随后把剥离的相位乘回去，得到两种用途不同的表示。

源盒较小时：用等效源表示部分和。 当，

因而

其中

可以理解为放在上的等效源强度。它不是原输入在 Chebyshev 点上的函数值，而是原源点数据经过带相位插值压缩得到的系数。[1, 原式(3.8)、(3.9)]

目标盒较小时：用局部势的采样值表示部分和。 当，

此时令

就得到

这时系数近似的是目标 Chebyshev 点上的部分和。前半程传递等效源，后半程传递局部势值；两种系数虽然使用同一个符号，含义并不相同。[1, 原式(3.10)、(3.11)]

Step 4. 沿双树执行，并在中间层切换表示

下面给出实际计算所需的完整递归。为方便书写，设为偶数。表示目标盒的父盒，表示源盒的四个子盒。

初始化。 取为目标树根盒，对每个源叶盒直接计算

这一阶段只在源叶子内访问原输入。[1, 原式(4.1)]

前半程：合并等效源。 对，遍历行树第层与列树第层的盒对，令

具体而言，将子盒上的等效源乘上新目标中心对应的相位，用插值矩阵转到父源盒的 Chebyshev 节点，再消去这些节点上的相位。这里使用的是新目标盒中心，并不是旧中心。[1, 原式(4.2)]

中间层：等效源转为势值。 在，两个盒的边长均为，两种表示都适用。计算

这一步是一个的小核矩阵乘法，将等效源在目标节点上求值。它不是单纯重新命名或排列系数。[1, 原式(4.3)]

后半程：插值局部势，再合并源贡献。 对，计算

这次是在目标侧插值：先对每个源子盒的势值消去主相位，从父目标盒插值到子目标盒，再乘回相位并求和。相位中心是各个，不能提前统一换成。[1, 原式(4.4)]

结束。 当为源树根盒时，它已经包含全部源点。对每个目标叶盒中的实际点，输出

至此不再需要访问原输入；所有源贡献都已包含在当前系数中。[1, 原式(4.5)]

这套公式之所以只需要较低的计算和存储成本，关键不只是秩与无关，还在于传递矩阵中的插值部分具有张量积结构。将一组个系数整理成一个数组，二维插值可以分两次一维乘法完成：

首尾的相位因子是逐项乘法；中间的插值从一般稠密矩阵所需的降为。父子盒之间的相对几何只有固定几种，因此与盒的位置和层号无关，可以预先计算并反复使用。不过，相位乘子仍依赖具体盒对；共享的是标准插值矩阵，不能把整套 FIO 传递矩阵都视为相同。[1, §4, pp.1743–1744]

二维每层至多有

个盒对。由前后递归及中间层切换得到成本

第二项来自中间层的小核矩阵作用，不能省略其因子。代入定理给出的保守插值阶数，得到

在固定精度下，以上就是。执行时仅保留相邻两层的系数，存储为；不必保存所有层、所有盒对的稠密传递矩阵。[1, §4]

原文所说的“无需预计算”，应理解为无需逐盒进行昂贵的核相关低秩分解。构造树、生成 Chebyshev 节点、保存少量公共插值矩阵仍然是必要的。实践中，当一个盒内的真实点数少于等效源数时，直接保留真实点可能更便宜；原文因此从较粗的源层开始，并在较粗的目标层结束递归。[1, §5]

振幅、多右端与数值证据

若核还带有振幅，且它具有低秩表示

则将输入改为，分别执行相位核的蝶形作用，再乘并相加即可。保持固定精度下近线性复杂度，还要求与无关。总体成本仍含这些额外通道的工作；作者建议同时处理它们，以共享相位计算。这里要求振幅本身具有低成本的分离表示，不能对任意振幅直接宣称常数不变。[1, §4, pp.1744–1745]

原文 Table 5.1 的二维 FIO 实验，在时给出：

每方向节点数	每盒系数数	抽样相对误差

这些结果说明，增加插值阶数可以明显改善精度，但张量积节点数与传递成本也随之上升。误差是对随机输入、在个目标点上与直接计算比较得到的，并不是整个矩阵的谱范数误差，也不是统一的浮点后向稳定性证书。[1, 原式(5.1)，Table 5.1]

这篇论文建立的是一条完整的分析型路线：利用相位证明互补低秩，再用显式插值把它实现为无需昂贵低秩预计算的快速求和。 它为一般 FIO 提供了清晰的数学基础；代价是张量积插值阶数可能高于矩阵真正需要的数值秩。这正是下一篇论文希望进一步压缩的部分。[2, §1.2]

3. Li 等人（2015）：把蝶形计算写成稀疏矩阵乘积

Butterfly Factorization. [2]

主要问题与贡献

这篇论文从一个具有互补低秩性质的矩阵出发，不要求知道其解析相位。作者希望构造

这就是 butterfly factorization，简称 BF。在局部秩固定的条件下，它由个稀疏因子组成，每个因子只有个非零元素，因此存储和一次作用都为。[2, 原式(1)、(7)]

此前的蝶形算法已经包含多层传递。本文的贡献是把它明确组织成可保存、可重用的稀疏矩阵分解，并给出两种不依赖显式相位的构造接口：

已有快速计算和的程序；
可以在时间内访问任意矩阵条目。

第一种接口还允许压缩已有快速算法或若干算子的复合；第二种接口适用于有快速条目评价公式的变换。局部基底通过 SVD 获得，有机会比固定阶数的解析插值使用更小的秩。[2, §§1.2,3]

Step 1: 从中间层构造初始分解

沿用原文记号，令为行树第层的第个节点，为列树第层的第个节点，并记

在中间层，两棵树各有

个节点，每个的大小约为。首先对全部中层块建立秩近似。

原文不是直接把正交奇异向量放在两侧，而是采用带奇异值权重的形式。若某块的截断 SVD 为

则定义

于是

注意右侧因子的下标为。奇异值保留在两侧，可以让后续递归继续辨别不同方向对原矩阵的贡献大小；这些因子本身并非全部正交。这里只对实际保留的非零奇异值取倒数，不是求整个矩阵的逆。[2, 原式(3)，§3.1]

固定目标块，把所有横向拼接：

类似地构造。所有则放入一个加权置换矩阵，把按源块排列的局部系数转换为按目标块排列。由此得到

这一分解的作用可以直接按块理解：

负责中间的重排与对角缩放，不是一般稠密的全局耦合矩阵。若每块秩均为，则有列，为阶矩阵；但只有个非零元素。[2, 原式(8)、(9)]

中间层的低秩计算:

方式一：用快速矩阵作用产生结构化随机采样。
对一个普通低秩块，随机 SVD 可以先生成两个瘦高 Gaussian 矩阵，计算

对做带主元 QR，得到列空间和行空间的近似正交基。随后只需对小核心矩阵

进行 SVD，再把小矩阵的奇异向量乘回两侧基底。原文 Algorithm 2.1 使用个随机向量，为附加过采样量。[2, §2.1]
对整个 BF，不能逐块读出再执行上述步骤。作者把随机输入组织成块对角矩阵，每个对角块有列。于是

一次批量作用同时得到许多块的采样，而总共只需个全矩阵输入向量。中间结果随后按分块处理；没有必要形成任何完整的矩阵。[2, §3.1]

Remark: 怎样用已有采样恢复小核心？

上面的描述了投影后的核心，但实现批量黑箱构造时，不应为每个块额外调用许多次完整的。一种与现有采样兼容的整理方式是：利用

解一个小规模最小二乘问题，例如

表示伪逆；实际可用 QR 或 SVD 求解。此式是对采样数据流的实现说明，并非原文单独编号的公式。它要求采样能辨认保留的子空间，且相应小问题不过度病态。若增加拟合用的随机探针，也应按相同的块对角方式组织，并保持常数批数；最终近似可用少量独立的全局测试向量校验。核心恢复始终在小矩阵上进行，从而保持原文次全矩阵作用的规模。

方式二：随机抽取少量行列。
如果任意条目可直接求值，则对每个中层块抽取少量行，做带主元 QR 找到重要列；再抽取少量列，做对应的行选择。交替更新几次后，用选出的完整列和完整行构造基底。
随后选择小索引集，在交叉子矩阵上拟合核心：

即

再对做 SVD，就获得该块的近似 SVD。固定时，每块只需读取数量级的条目，而不是全部个条目。[2, Algorithm 2.2]
这种抽样方法依赖行列信息能够被少量样本捕获。原文特别讨论了与坐标向量不相干的情形，并明确表示没有为 Algorithm 2.2 给出一般的误差与成功概率定理。因此“具有互补低秩”本身不能保证任意抽样策略都可靠。[2, §2.2]

Step 2: 向两端递归压缩中层因子

中间层的还含有数量级的数据，不能直接作为最终存储。接下来利用下一层的互补低秩性质继续分解。

以左侧为例。把按目标树的两个子节点分成上半部与下半部，记为和。再把来自两个相邻源节点的因子并在一起：

这两个矩阵各有列，分别对应“目标盒缩小一半、源盒扩大一倍”后的相互作用。由于新的矩阵块、仍低秩，这些拼接基底在相应近似精度下也可以重新压缩到约列：

这里的大小为。例如，对上半部拼接矩阵作截断 SVD ，令、；下半部同理。这样继续把奇异值留在待压缩的一侧。[2, 原式(5)、(13)–(16)]

把所有这些小分解按树组织起来，就得到

右侧对做相同操作：

将它们代入中层分解，便得到本节开头的完整 BF。每一层只保存连接相邻尺度的小矩阵，最后的只与常数规模的物理叶节点相连。[2, §§3.2.1–3.2.2]

Step 3 应用保存的稀疏因子

对新输入，按照从右到左的顺序计算：

用把输入映射为最细层局部系数；
依次作用，到达中层；
用重排并缩放；
依次作用，向目标细层传播；
用输出目标值。

这些步骤不再需要评价核函数或做 SVD。多个右端也可以共同通过相同的局部矩阵。但最终速度仍受局部秩、因子布局和小矩阵乘法效率影响，不能只由稀疏因子的数量判断。

构造成本、存储与数值结果

令表示一次调用已有或程序的成本。在固定秩下，原文的复杂度为：

构造接口	中层构造	递归压缩	总构造成本
快速作用			若，则为
任意条目访问

第二阶段为什么仍需？在第层，许多待压缩基底仍有行，该层工作为；从中层求和，

因此，在上述访问模型及低秩采样条件下，2015 年 BF 的构造为，或在快速矩阵作用接口下达到典型的；构造尚未达到。[2, §3.3，原式(20)]

若同时保存全部中层因子，峰值内存为。条目访问版本还可以逐个生成并压缩中层块行，把构造内存降到；不能把默认中层存储的峰值与最终分解存储混为一谈。最终 BF 的存储和作用均为，若显式保留秩依赖，内部传递通常带有因子。[2, §3.3，Table 1]

原文 Table 2 的一维 FIO 实验，在固定秩时给出：

矩阵阶数	构造时间，分钟	一次作用时间，秒	抽样相对误差

这里是原论文的历史 MATLAB 实验，用于说明构造和作用是两种不同的成本；不能与另一篇论文在不同硬件、不同核函数上的绝对时间直接比较，也不能将其直接求值基线等同于已经预存矩阵的最优 BLAS 乘法。误差在随机选取的个输出点上计算，并非矩阵谱范数误差。[2, §4，原式(21)，Table 2]

本文还展示了只通过快速作用程序来压缩复合算子的能力。对具有适当 FIO 结构的复合，新的矩阵仍可能保持互补低秩，但所需秩会增加；这并不意味着任意两个互补低秩矩阵的乘积都自动具有同样的低秩界。[2, §4, Example 3]

2015 年 BF 的核心推进是：用接近实际数值秩的代数基底，把蝶形作用变成通用的稀疏矩阵链，并将黑箱构造成本显著降低。 对反复使用同一算子的任务，这种预计算可以换取较便宜的每次作用；但若只作用少数几次，构造仍可能成为主要负担。

4. Pang–Ho–Yang（2020）：用骨架索引实现近线性构造

Interpolative Decomposition Butterfly Factorization. [3]

主要问题与贡献

2015 年 BF 的瓶颈在于：先在中间层生成许多较长的基向量，再把它们向两侧压缩。虽然最终表示很小，中途仍处理了量级的数据。

2020 年的 interpolative decomposition butterfly factorization，简称 IDBF 改变了构造顺序。它从“行叶子对列根”和“行根对列叶子”两端开始，保留代表性行列及其插值矩阵，再向中间推进。下一层所需的信息仍可通过原矩阵条目获得，不必先形成中层的大基底。[3, §§1–2]

其完整分解仍是一串“左侧插值—中间骨架—右侧插值”因子。沿用原文的，但用方括号上标明确表示递归操作次数，可以写为

这种表示与原文式(1)、(25)具有相同结构，操作次数由树的停止准则决定。方括号编号用于避开原文末层下标的一处偏移，下面将给出具体计数。现在是行、列插值及重排因子；右侧的已经是作用因子，不应再照搬 2015 年公式给它额外加一个伴随。

在互补秩有界、每个条目可用工作准确评价、且有足够有效的预选采样点时，原文给出固定参数下的构造、作用和存储成本。这是一项关于构造方式的推进，不是把稀疏因子重新命名为 ID。[3, §1，原式(1)，§2.4]

Step 1: 用少量采样行找到列骨架

先考虑一个数值秩约为的矩阵。插值分解的基本形式是

是一小组实际列索引，称为列骨架，表示其余列如何由骨架列线性组合。与 SVD 的区别在于，左侧保留的是原矩阵的列，而不是另外生成的一组正交基向量。

对精确秩为的矩阵，做带列置换的秩揭示 QR：

令对应置换后的前列，对应其余列，则

实际通过三角求解计算，无需显式形成。在数值低秩情形，上述等式换成近似。传入的是试验秩上限，实际保留的骨架数记为，不超过及矩阵的可用维数。[3, Lemma 1，原式(2)–(6)]

如果，直接对整个高瘦矩阵做 QR 仍然很贵。作者只选取一个含约个行索引的集合；为与总递归次数区别，这里将原文的采样集合记为。对

做相同计算，得到列索引和插值矩阵。是过采样参数。固定时，计算只需运算和存储，成本不再与全部行数成正比。[3, §2.1，原式(7)]

这里有一个重要的实现区别：计算不要求立刻把完整的读出来。 可以先只保存列索引，以后需要哪一个骨架条目，再向原矩阵请求该条目。这是 IDBF 能避免巨大中间矩阵的基础。

对矮宽矩阵，则对其伴随执行相应的列 ID，得到行 ID：

其中为实际行索引，为行插值矩阵，成本为。[3, 原式(11)]

Remark: 采样点为什么采用 mock-Chebyshev 分布？

少量行能够决定整个列空间，需要这些行充分反映矩阵的变化。如果只在不具代表性的位置采样，就可能完全漏掉某些方向。

原文比较了均匀随机采样与 mock-Chebyshev 采样。后者是在已有的等距离散点中，寻找距离理想 Chebyshev 节点最近的一组点。例如，对有个有序等距行节点的区间，以个 Chebyshev–Lobatto 节点为示例，可以先把

映射到对应的行索引范围，再取最近的可用索引；实现时需处理重复索引与采样数超过可用点数的情形。这个式子只是最近邻选点的具体说明，节点仍来自原离散集合。原文没有限定必须采用这一组包含端点的 Chebyshev 规则。

这种分布在区间端点附近较密，对许多来自光滑物理核的局部插值更有效。这些点是用来发现骨架的采样点，并不意味着最终的骨架列就是 Chebyshev 点。 最终哪些列保留，仍由带主元 QR 决定。[3, §2.1]

原式(9)用 QR 对角项相对首项的大小判断截断。对复数问题应比较其模，例如把

作为出现可忽略方向的信号。但采样矩阵上的截断不能无条件保证整个的误差；其有效性依赖采样是否覆盖了真实子空间。原文引言明确要求矩阵具有合适的预选代理点，数值部分也强调局部容差与全局误差并不直接相等。[3, §§1,2.1,3]

Step 2: 从行叶子与列叶子提取两侧骨架

原文称这一步为 leaf-root complementary skeletonization，LRCS，即叶—根互补骨架化。这里用表示当前处理的互补低秩矩阵。

将行索引与列索引分别分成叶集合

开始时。

先压缩行。 对每个矮宽块做行 ID：

令

则

再压缩列。 对每个做列 ID：

令、，得到

这就是原文式(14)的内容。两端的都容易保存，中间的却可能仍很大。算法的关键是：此时只记录的行列索引，不显式计算全部。 后续若请求，把骨架索引追溯回原矩阵并评价一个条目即可。[3, §2.2]

这一点是 ID 比任意低秩基底更适合本算法的原因。经过普通基底变换，中间矩阵的一个条目可能需要一长串求和；经过骨架选择，中间矩阵仍是原矩阵的子矩阵，保留了廉价条目访问能力。

Step 3: 分裂矩阵，并合并相邻骨架叶子

仅仅对同一个反复做相同的叶压缩并不够。一个叶子里已经只剩代表性行列，继续压缩通常没有空间。要继续降低自由度，必须同时改变它与另一侧区间的配对关系。

作者引入 matrix splitting with complementary skeletonization，MSCS。先按行、列树根节点的两个孩子，把当前矩阵分成四块：

分别对每个做 LRCS：

不同象限可选择不同的骨架索引，因为它们代表的相互作用不同。不能把某个叶子针对一个象限选出的骨架，未经检查就用于另一个象限。

为了说明这些局部分解怎样组装，给所有空块显式补零。原文式(16)等价于

于是。的交叉位置用于匹配两侧不同的排列，并非可以随意省略的装饰。对输入，这一结构具体计算的是

每个中间块仍只由原矩阵的骨架索引描述。[3, §2.3，原式(15)、(16)]

下一次递归之前，在每个子问题内合并相邻的骨架叶子：

新叶子至多包含两组各约个的代表索引，因而为下一轮提供了大小的候选集合。下一轮先缩小配对区域，再通过 LRCS 重新选择骨架。[3, Figures 4–6，§§2.2–2.4]

因此，一次递归循环由根部的象限分裂、各象限的 LRCS，以及为下一轮准备的骨架合叶组成。合叶本身不额外增加一组插值因子；重新压缩发生在下一轮的 LRCS 中。这种尺度重组使后续处理始终局限在小候选集合中，而不需要重新访问整块矩阵。

Step 4: 递归到中间尺度，只在最后形成小核心

记第一次 MSCS 为

对的每个非零核心块继续做 MSCS，再按已有排列组装，得到

所以

重复操作，得到

这里的每个都表示带有固定排列结构的一组核心块，不是把所有核心填成一个大稠密矩阵。[3, §2.4]

停止位置可以直接从树来判断。初始每侧有个叶子；完成次“分象限—骨架化—合叶”后：

$核心块数每个核心的可压缩树叶数$

当为偶数时，经过次操作，每个核心只剩一个叶子。此时核心的实际行列数为，就可以直接评价其全部条目并保存；无需再分裂。原文的稀疏因子链就是这一递归的组装形式。[3, §2.4, p.A1106]

对新向量的作用顺序为：先做，再做最终核心，最后做。与 BF 一样，构造完成后的作用阶段只使用保存的因子。

从实现角度，可以把整个构造概括为以下流程：

保存当前子问题的行列索引及二叉树关系，并提供从当前索引到原矩阵条目的访问函数。
将当前核心按根节点分成四个象限，分别用采样行 ID、列 ID 提取两侧骨架。
保存小插值矩阵及必要的排列；核心只保存新的骨架索引。
合并相邻骨架叶子，将形成的子问题交给下一轮。
当树不再具有可继续配对的层级时，直接生成最后的小核心块。

这五步说明了 IDBF 与“把所有互补块分别做 ID”的区别。后者仍可能反复读取大量条目；IDBF 通过骨架索引的继承，使每轮只处理上一轮留下的少量候选。

复杂度

完成次操作后，所有核心的叶数相加为

因此，虽然子问题越来越多，每个子问题的树也相应变小，每轮需要处理的局部骨架问题总数仍为。合叶后，一个待压缩方向只有个候选，代理采样也只有个点；局部 QR 和三角求解需要工作，保存的插值矩阵有个条目。

固定过采样参数，按原文的叶宽与秩计数，有

原文隐藏在中的主要参数依赖由此可见。减少每个局部问题的候选规模，并保持核心条目的直接访问，是构造加速的原因。[3, §§1,2.4]

若把作为独立变化的参数，还应单独计入初始物理叶层。一个更完整的计数写法是

这两式是把原文流程中的端点工作显式补入后的计数；在通常的固定叶宽、固定秩设置下，都恢复近线性结论。实际局部秩还要截断在可用行列数以内，不能机械地为一个只有列的块保留个独立骨架。

这里的“kernel-independent”也有明确含义：构造无需知道核的解析相位公式，但仍需索引的几何或顺序信息、可靠的少量采样，以及快速条目访问。若每个特殊函数条目都要通过长递推或昂贵求解才能获得，其生成成本就必须加入上述复杂度；仅仅拥有一个返回条目的函数，并不等于拥有的条目访问。[3, §§1,2.1]

奇异核 — 在外层分块中使用 IDBF

实际高频积分方程的矩阵可能并非整体互补低秩，尤其是在对角线附近含有奇异或近奇异相互作用时。原文后两个例子采用 HODLR 型外层结构：对矩阵分块，将适合的非对角块交给 IDBF，继续递归处理对角块，直到最小对角块直接保存。[3, §3, Examples 4–5]

示意地写为

其中使用蝶形表示，再作同样分块。这里“非对角块可压缩”指它们具有适当的互补低秩结构，并非整个非对角块必须是普通低秩矩阵。

外层层级与每个块内部的蝶形层级会叠加。设外层深度为，止于固定大小的对角叶。若第个外层有个规模为的蝶形块，在固定秩条件下，总工作约为

这解释了为什么原文 EFIE、CFIE 实验给出的是量级，而不是整张矩阵直接套用一次 IDBF 的。这些实验展示的是矩阵压缩与快速作用，为迭代求解或进一步构造求解器提供工具；它们本身不等于完成了一个近线性的直接求逆算法。

数值实验结果

论文使用单线程 MATLAB 实现，设置初始叶宽、过采样参数，并比较不同容差与秩参数。前三类实验包括的设置；后两个积分方程例子使用。它的五类实验如下：

实验	被作用的矩阵或变换	所验证的重点
Example 1	相位为的一维 FIO	mock-Chebyshev 与随机采样的精度差别，秩参数和局部容差的影响
Example 2	构成的 Schlömilch 展开求值矩阵	特殊函数条目可快速生成时的构造与作用
Example 3	非均匀 Fourier 变换矩阵	非均匀点集上的适用性
Example 4	二维 EFIE 阻抗矩阵	非对角块采用 IDBF 的分层作用
Example 5	二维 CFIE 阻抗矩阵	另一类高频积分方程矩阵上的分层作用

其中 FIO 相位使用

它与 2015 年实验中的具体相位参数不同。作者为简化展示，直接对整个核使用 IDBF，而未先作极坐标或按原点附近区域分开处理；原文将附近的非光滑性列为局部秩随规模略有增加的原因。[3, 原式(28)，Example 1]

实验主要支持以下几项结论。

第一，采样点选择会直接影响精度。 对文中的 FIO，mock-Chebyshev 采样明显优于均匀随机采样。后续例子因此采用前者。这是针对所测试核与采样设置的经验结论，并不是确定性采样对任意矩阵都优于随机采样。[3, Figures 9–11]

第二，局部容差与最终误差不同。 论文比较了和等设置，但明确指出最终误差不一定达到：局部截断误差会通过层级传播，测试秩不足也会限制精度。因此，不能将设置直接写成整个变换达到了双精度舍入误差。[3, §3, p.A1108]

第三，理论计数与有限规模测量需要分开。 前三个例子报告的构造和作用曲线主要处于型增长附近，并讨论了较大规模时非零元数量趋近的现象；后两个外层分块例子报告型增长。固定秩模型下的不意味着所有实验都已经测得同样的增长规律。[3, Figures 9–16]

与前两篇类似，误差采用抽取个输出点后的相对向量误差。它适合检查一次作用的实际精度，但不能作为任意输入下的完整算子误差证明。[3, 原式(26)]

IDBF 最重要的贡献由此可以概括为：用双侧骨架化保留原条目访问能力，再通过根部切分和叶部合并，使构造全程只处理近线性数量的小问题。 它解决了 2015 年方法先形成较大中层表示的构造负担，同时把采样质量、局部秩和误差传播作为仍需关注的条件保留下来。

5. 三篇论文之间的关系

为了横向比较，本节单独用表示输入向量的长度。因此，2009 年二维问题中，另外两篇中。下表的近线性复杂度均以固定局部秩、固定精度及相应访问条件为前提。

方法	构造需要的信息	核心构造方式	主要存储
2009 分析型蝶形 [1]	可评价的相位，以及残差的解析或插值结构	显式调制、共享插值；没有逐块代数预压缩	流式执行的系数为
2015 BF，条目接口 [2]	任意条目访问，以及有效的低秩采样	中层随机 SVD，再向两端压缩；构造	最终因子；逐块构造可达到同阶峰值
2015 BF，作用接口 [2]	快速程序	结构化随机探针；典型构造	最终因子；构造峰值
2020 IDBF [3]	廉价条目访问、索引结构、有效代理采样点	从两端向中层递归骨架化；构造	因子

第一行的是按需计算相位、只保存相邻层系数时的存储，不能与“把全部相位和传递矩阵预存下来”的另一种实现混用。其余行的存储主要指可反复使用的分解。上述比较对应 [1, §4]、[2, §3.3]、[3, §§1–2.4]。

这也说明，2020 年 IDBF 并不在所有条件下替代 2009 年方法：如果相位结构明确，显式插值能够避免寻找局部骨架，并共享标准插值矩阵；如果只有快速作用程序而不能廉价访问条目，2015 年的作用接口仍有独立价值。应当根据可获得的信息选择构造路线。

三篇文献分别给出三项实质推进：

2009 年：让局部低秩成为可直接计算的解析结构。 残差相位给出互补尺度的依据，两种插值表示给出完整传递公式，张量积结构则把理论上的低秩转化为较便宜的实际计算。它回答的是：面对复杂振荡相位，为什么可以压缩，以及怎样避免昂贵地寻找基底。
2015 年：让蝶形成为通用、可保存的代数分解。 中层构造与双侧递归把矩阵变成稀疏因子链；随机化接口使已有快速程序和条目查询都能用于构造。它回答的是：没有可直接利用的相位公式时，怎样建立一个紧凑、适合重复作用的表示。
2020 年：让表示的构造本身也保持近线性工作量。 ID 保留原矩阵索引，MSCS 维持跨尺度骨架的可访问性；整个构造不再先处理大量长基向量。它回答的是：怎样在找到这些局部表示的过程中，也避免付出超近线性的代价。

因此，这三篇论文的递进关系不能简单概括为“插值换成 SVD，再把 SVD 换成 QR”。低秩分解只是局部工具；决定整体复杂度的是分解顺序、候选集合的大小，以及下一层信息能否廉价取得。

6. 精度、伴随与实际效率

局部低秩误差不等于完整变换误差

对 2009 年的局部核近似，如果

那么直接用精确源数据求一次局部部分和，可得

这是一个点态绝对误差界。若输出中存在强烈相消，它并不能直接转化为相对输出误差；递归中的源系数还带有此前各层的近似误差。[1, §2]

BF 与 IDBF 同样需要区分局部截断与完整乘积。用两个因子说明，

因而误差不只取决于的大小，也取决于传递因子的范数。层数增加、插值矩阵较大、骨架选择不充分时，局部误差都可能被放大。这一简单恒等式解释了为什么“小块都达到容差”还不是全局证明；[3, §3] 对分层 ID 的最终精度也明确保留了这一限制。

阅读数值结果时，应分别看局部近似条件、完整变换的误差测量，以及实际浮点程序的稳定性。三篇论文的抽样误差支持其测试实例上的有效性，但不能自动推广为所有输入下的谱范数保证或后向稳定性定理。

保存同一套因子，可以直接得到它的伴随

设构造出的近似矩阵为

则同一表示的伴随为

因此，只需反向遍历计算图，并将每个局部作用改成共轭转置，就能作用。这是矩阵乘积的直接性质，不需要为伴随独立寻找另一套骨架。对 2009 年按需生成的插值传递，也可以按相同原则反向执行其固定参数下的线性计算图。

这里得到的是同一个近似算子的伴随。它与对真实单独构造另一份近似不同；它也不等于逆矩阵。若真实是酉矩阵，压缩后的仍可能有正交性缺陷。忽略实际乘法中的舍入，若

则由可推出

这说明恢复接近恒等的往返作用，首先需要控制完整算子误差；仅仅使用相同因子的伴随，还不能保证往返误差达到机器精度。这个估计是一般矩阵代数推论，不是三篇论文另行证明的快速变换稳定性结论。

常数、预计算与交叉点

同样写成，三条路线的实际代价仍有明显区别。

2009 年方法在二维保存个局部系数，传递成本主要为；BF 和 IDBF 的内部小矩阵通常带有秩的平方成本，IDBF 构造还包含秩的立方成本。更少的系数未必意味着更少的总时间：解析插值可能有更大的秩，却拥有共享的张量积运算；代数方法可能秩较小，却要进行更多不规则的数据重排与小矩阵操作。

预计算还需结合使用次数评价。若构造成本为，每次作用为，同一算子使用次，则总成本为

这与只比较单次作用时间是不同的问题。多个右端可以共享相位、索引和局部矩阵，但批量计算的速度收益仍要结合实现测量；三篇论文没有给出统一的多右端交叉点。

对于一个新的变换，最有用的检查顺序是：先确认互补尺度上的数值秩，再确认核条目或快速作用接口的成本，随后检查局部表示能否跨层复用，最后比较完整构造、作用、存储和误差。只有把这些环节放在同一计算流程中，才能判断蝶形结构是否真正带来优势。

共同思想与后续问题

蝶形算法将一次全局稠密求和拆成许多尺度互补的局部相互作用。在每一层，每个相互作用只保留少量系数；当目标区域缩小、源区域扩大时，通过小规模传递重新组织这些系数。它获得快速作用的原因，是局部表示可以跨尺度接续，而不是矩阵在某一个分块尺度上看起来低秩。

三篇论文分别把这一思想推进到解析构造、通用代数分解和近线性骨架构造。它们也共同指出了实际算法仍需面对的几个问题：局部秩是否随规模或精度增长，少量采样是否可靠，误差如何穿过整张计算图，以及所节省的算术是否足以抵消插值、重排和预计算的常数。

从这组文献中最值得保留的方法，是把“证明结构”“建立表示”和“作用表示”分开研究，再检查三者是否真正衔接。明确的相位结构可以降低局部表示成本；保留真实骨架可以降低后续信息获取成本；保存小传递矩阵可以降低重复作用成本。蝶形算法的进展，正发生在这些环节被更紧密地组织起来的时候。

Remark: 原文定位、记号整理与核查说明

本文以三篇所附期刊 PDF 为主要依据，保留各文的主要符号，用原文节号、公式号和图表号定位结论。少量基础推导用于解释计算流程，不作为原论文的新增定理。复杂度始终说明局部秩与访问条件；实验数据均标明其误差测量对象。

为避免继承原文的排版或编号问题，作了以下整理。

[1, p.1743] 将“其余步骤”写成，但中层有个盒对，且紧接着的总复杂度也包含。本文按中层计数写为。
[3, 原式(24)] 左侧组装矩阵的最后一个对角块应为，右侧的最后一个对角块应为。本文采用原式(16)的正确装配关系，显式标出零块，避免复制这两个笔误。
[3, §2.4] 原式(25)中从到的含端点因子编号，与正文“共次递归”及叶数计数存在一层偏移。本文用方括号上标记录实际 MSCS 次数，并以“每核心只剩一个可压缩树叶”为停止准则。
[3, p.A1100] 采样 QR 中的只具有采样行数，不能直接当作完整矩阵的列骨架。本文始终将骨架定义为置换选出的原列，并将 QR 对角阈值解释为秩选择信号。

文稿还对残差相位消去、两侧调制插值、张量积收缩、中层奇异值缩放、采样核心恢复、行列 ID、MSCS 块排列与伴随次序做了小规模代数核验。这些检查用于排除符号、维数和转置错误，不是对三套完整快速算法的重新实现，也不替代原论文的误差分析。

本文没有将“局部低秩”“抽样误差较小”或“使用同一分解的伴随”提升为全局后向稳定性结论，也没有用不同年代、不同硬件的计时拼接跨论文速度排名。