Hierarchical Matrices
REFERENCES
- A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices, by W. Hackbusch. Computing, 62:89–108, 1999.
- H²-matrix approximation of integral operators by interpolation, by W. Hackbusch and S. Börm. Applied Numerical Mathematics, 43:129–143, 2002.
- On generators of quasiseparable finite block matrices, by Y. Eidelman and I. Gohberg. Calcolo, 42:187–214, 2005.
- Some Fast Algorithms for Sequentially Semiseparable Representations, by S. Chandrasekaran, P. Dewilde, M. Gu, T. Pals, X. Sun, A.-J. van der Veen and D. White. SIAM Journal on Matrix Analysis and Applications, 27(2):341–364, 2005.
- A Fast ULV Decomposition Solver for Hierarchically Semiseparable Representations, by S. Chandrasekaran, M. Gu and T. Pals. SIAM Journal on Matrix Analysis and Applications, 28(3):603–622, 2006.
- Fast algorithms for hierarchically semiseparable matrices, by J. Xia, S. Chandrasekaran, M. Gu and X. S. Li. Numerical Linear Algebra with Applications, 17:953–976, 2010.
后续进展补充
- Adaptive multiplication of H²-matrices with block-relative error control, by S. Börm. Numerical Algorithms, 101:65–96, 2026.
- Quasi-optimal Hierarchically Semi-separable Matrix Approximation, by N. Amsel, T. Chen, F. Duman Keles, D. Halikias, Cameron Musco, Christopher Musco and D. Persson. SIAM Journal on Matrix Analysis and Applications, 2026.
本文将上述文献按三条主线组织:
1. 层次矩阵核心思想
从一个低秩块开始
考虑
设行指标集为
那么这个块只需保存
先把
这里的
因此,“秩小”必须与精度、划分和问题规模一起讨论。整个矩阵可以满秩,只有适当选择的非对角块低秩。
为什么局部相互作用会低秩?这类低秩性的一个典型来源是积分算子。若目标簇
离散化后就得到低秩矩阵块。邻近区域可能包含核的奇性,需要继续细分,最后用小稠密块直接保存 [1, §5;2, §§2–3]。另一个来源是消元。稀疏矩阵的逆、Schur 补通常是稠密的,但不同子区域之间的影响可能只经过少量接口变量。此时,条目虽然增加,跨区域相互作用的有效维数仍然很小。这也是 SSS、拟可分和 HSS 算法能够用于直接求解的原因之一[3, §§2–4;5, §1;6, §§1,5]。
低秩分解解决一个块,层次组织解决整个矩阵。 只找到几个低秩块还不够:需要决定哪些块一起存储、它们能否共享基底,以及运算之后怎样继续维持小秩。
分块树与顺序链
在树结构中,把指标集递归分成子簇。例如
相应的矩阵分成四块。符合低秩条件的块停止细分;其他块继续递归。这样,大范围但平滑的相互作用用大低秩块表示,局部复杂相互作用用小块表示。
在顺序结构中,把变量分为
这里的指标是块指标。若这些整个矩形子矩阵的秩都小,就可以用少量状态跨越切口传递信息,形成 SSS/拟可分表示。要求不是“每个小块
| 结构 | 被压缩的对象 | 保存的主要数据 | 关键计算组织 |
|---|---|---|---|
| 可接受分块中的低秩子块 | 各块各自的低秩因子、小稠密块 | 块递归与截断 | |
| 可接受分块,并共享簇基底 | 叶基底、传递矩阵、耦合矩阵、近场块 | 向上压缩、块间作用、向下展开 | |
| 拟可分/SSS | 每个顺序切口两侧的整个非对角矩形 | 输入、输出和状态传递生成元 | 沿顺序链扫描 |
| HSS | 每个树节点对应的整个非对角块行/块列 | 叶基底、父子传递、兄弟耦合、叶对角块 | 树上压缩、消元与合并 |
这些名称描述的结构相近,但不能不加条件地互换。用后来的常用术语说,仅要求每层两个兄弟非对角块低秩的表示称为 HODLR;HSS 还要求跨层共享嵌套基底,并压缩整个非对角块行/块列。HODLR 在这里仅用于帮助比较,并非两篇 HSS 原文采用的术语。[5, §2;6, Definitions 1.1–1.4]
若各层的低秩块都独立保存基底,在平衡二叉树上,每层通常需要
若父节点基底在子节点上的限制满足
就不必保存每层越来越长的
上式中的限制记号也可写成
这个计数同时解释了复杂度结论的条件:嵌套结构、受控的秩和适当的叶大小,三者缺一不可。 下一步是把这些想法变成具体算法。
2. Hackbusch(1999):为分块低秩矩阵建立近似算术
A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices. [1]
主要工作与贡献
这篇论文提出
原文具体研究两种格式。第一种非常简单,只递归细分对角块;第二种在对角线附近保留更细的分块,用于逼近积分算子。多维几何划分的一般理论留给后续工作,不能把后来完整的
Step 1 用树和块划分定义存储格式
原文以
对于
其中
对于秩一块,
这里
这说明“大矩阵乘低秩矩阵”可以退化成少量矩阵向量乘法。
Step 2 用小矩阵分解完成截断
低秩格式通常不对加法封闭。若
则
秩可能由
下面用 QR 加小 SVD 写出同一截断操作,便于直接实现。设待压缩块为
只保留前
大维度只出现在瘦 QR 中,SVD 的尺寸至多为
而不是对
原文通过秩一因子的 Gram 矩阵说明这一思想;这里的 QR 写法是对该计算的等价整理。它也更清楚地说明如何按容差自适应选择保留秩。
Step 3 递归完成加法和乘法
加法按共同块划分逐块进行:稠密叶直接相加,低秩块先拼接因子,再截断。
乘法仍遵循普通分块公式:
区别在于执行方式。例如
第一项递归计算,第二项利用低秩因子计算;相加时再次截断。两个低秩块相乘满足
真正需要形成的中间矩阵只有
需要注意,截断改变了算术:先分别截断再相加,通常不同于相加后统一截断。因此,普通矩阵乘法的结合律不能不加分析地转移到固定秩近似算术中。
Step 4 用 Schur 补递归求逆
设
分块求逆公式为
原文算法依次执行:
- 递归计算
。 - 利用低秩因子计算
,形成并截断近似 Schur 补 。 - 递归计算
。 - 把
代入上面的四个块,必要时再截断。
在第一种格式中,
这里要求递归主块可用;矩阵整体可逆并不自动保证每个指定主块可逆。原文公式也不构成对任意非正定矩阵均无需主元处理的保证。
例子:三对角矩阵的逆
把三对角矩阵按连续指标分成左右两半,跨越分界的
只修改右半矩阵边界处的一个对角元,Schur 补仍是三对角矩阵。递归下去,逆矩阵的对应非对角块也至多秩一。
因此,在递归消元可进行的条件下,第一种
第二种格式:在近场保留更细的块
相邻区间之间的核可能在接触点附近变化很快,把整个兄弟块立即压成低秩并不合适。原文为此引入邻接型
只有靠近接触点的左下子块
原文以对数核的分片常数配点离散为例。对实值问题可写成
这里把原文的对数核明确写成实值的
于是
直接给出秩不超过
原文由 Taylor 余项得到该离散模型的误差界
这里是矩阵诱导无穷范数;近场小块直接计算。该结论对应这一特定核和分块,不是所有
复杂度与适用边界
对原文研究的平衡递归格式,固定块秩时有以下结论。[1, §§3.1–3.6,5.2]
| 运算 | 原文的渐近阶数 |
|---|---|
| 存储、矩阵向量乘法 | |
| 截断加法 | |
| 截断矩阵乘法、近似求逆、近似 LU |
这是对给定格式的算术成本,而不是对任意矩阵都能以固定秩达到指定精度的承诺。求逆时若所需秩持续增长,或截断误差被病态问题放大,实际成本和求解精度都要重新评估。
这篇论文建立了完整的基本循环:分块低秩表示,执行普通代数公式,再把结果截断回压缩格式。 下一篇进一步减少了不同块、不同层之间重复保存的基底。
3. Hackbusch–Börm(2002):用插值构造嵌套的 基底
H²-matrix approximation of integral operators by interpolation. [2]
主要工作与贡献
这篇论文给出一种可直接实现的
其主要结果是:在合适的块稀疏性、叶大小和核光滑性条件下,构造可分析的积分算子近似,并把存储和矩阵向量乘法降到
Step 1 从 Galerkin 矩阵和几何划分出发
考虑定义在
采用局部支撑基函数
其中
质量矩阵
对簇
判断块
Step 2 用核的点值构造低秩耦合
在
定义
代入积分并分离变量,得到
Step 3 利用插值的再现性得到精确嵌套
设
定义传递矩阵
把上式与
因此只需要显式保存叶簇的
Remark: 多项式次数、节点数和秩的记号
原文主体用
例如,附录的 Chebyshev 根节点为
先把这些点仿射映射到各坐标区间,再取张量积即可。
Step 4 三次遍历完成矩阵向量乘法
设要计算
第一次:从叶到根压缩源向量。
叶簇计算
内部节点利用已计算的子簇系数:
第二次:在可接受块之间传播作用。
所有目标系数初始为零,随后逐个可接受块累加
第三次:从根到叶展开。
将父簇贡献传给每个子簇:
到叶簇后计算
最后,把不可接受叶块的直接乘法
计入结果。近场块不能因为远场采用了插值而被省略。
三次遍历各有明确作用:把源数据变成系数,在小系数空间中处理相互作用,再把目标系数变回物理变量。父子间传递保证同一份信息不会在每层从头计算。
复杂度
令所有簇的基底维数为
并要求每个簇在块划分中参与的块数由常数
叶基底和近场块也满足同一阶数。[2, 原式(15)–(17)、§4.1]
这个线性结论不是把
构造时,耦合矩阵只需核点值;叶基底和近场块还需数值积分。原文分别计入了插值、求积阶数和建树成本;典型几何建树本身可需要
误差分析和变阶扩展
原文假设核渐近光滑,即离开
配合几何可接受性和 Chebyshev 插值,得到
若导数常数具有原文 Remark 5.2 给出的阶乘型增长,并满足相应的
Remark: 插值误差常数包含什么?
原文 Lemma 5.1 中的常数为
其中阶乘分母来自插值余项,对数因子来自张量 Chebyshev 插值的稳定性,导数常数来自核的渐近光滑性。这里按原文在
则阶乘增长被抵消,得到由
论文 §6.3 还讨论变阶策略:小簇用较低阶,大簇用较高阶。如果子簇次数低于父簇次数,父多项式在子簇上的再插值便只是近似,嵌套关系会引入额外误差。适当地控制阶数随层级的增长,可使所有层的总成本保持
其计数可以这样理解:距叶层高度为
大多数细簇保持低阶,少数粗簇才使用高阶。原文单位圆单层势实验采用分片常数 Galerkin、
对于配点法、双层势等变体,主要变化落在叶基底的定义上,核点值耦合和多项式传递结构仍可复用。[2, §6]
这篇论文把“嵌套基底”从抽象结构变成了可计算的公式:叶上做积分,块间取核值,层间评价插值多项式。
4. Eidelman–Gohberg(2005):用最少的生成元描述跨切口相互作用
主要工作与贡献
这篇论文主要是一篇带有新结果的理论综述。它系统回答三个问题:拟可分矩阵怎样由生成元表示,生成元的最小维数是多少,同一矩阵的不同生成元之间有什么关系;同时给出生成元约简、构造、乘法和求逆算法。原文特别指出,局部生成元的最小性判据 Theorem 4.5 是其中的新结果之一。[3, §§1,4]
它的重要性在于,快速算法的成本由生成元维数控制。仅知道矩阵存在低秩结构还不够,还要能识别并移除表示中的冗余状态。
本节沿用原文:
Step 1 用输入、传递和输出生成元表示非对角块
下三角部分写为
其中
上三角部分类似:
其中
相邻块之间的空乘积为单位矩阵,故
例如,从第
Step 2 由切口秩理解最小状态维数
原文定义累积生成元
于是每个下三角切口矩形都有
所以它的秩不会超过
则一定可以找到
从系统角度看,
原文 Theorem 4.5 把这种全局最小性化为小矩阵的局部条件。对一组阶数为
检验这些条件不需要形成所有大切口子矩阵。
Step 3 两次扫描,把已有生成元约简到最小
假设已有
第一次,从左向右,使累积输入基底正交。
先做识别实际秩的紧致 LQ 或 SVD 分解
对
并按列拆成
同时更新
这一步删除输入端的冗余,并使新的
第二次,从右向左,使累积输出基底正交。
先分解
对
并拆成
最后令
若所有分解仅删除精确零奇异值,则所得生成元表示同一个矩阵,并达到最小阶数。若按容差删除小奇异值,则得到近似约简,误差需要结合传递过程分析,不能把局部阈值直接当成整个矩阵的相对误差。[3, §6]
局部正交关系还可写为
或
它们分别对应原文的右正规形和左正规形。正交化不仅为降秩服务,也能减少状态坐标选择引起的数值放大。
Remark: 如果输入是稠密矩阵,怎样构造生成元?
Theorem 6.4 不要求先有生成元。算法先压缩第一列下方的块:
然后在第
把已经压缩的历史与新的一列拼接:
最后取
Algorithm 6.5 将这里的分解替换为截断 SVD,允许指定最大生成元阶数或奇异值阈值。实现时必须把截断后的右因子按旧状态列与新输入列拆成
Step 4 两个方向的状态递归完成乘法
计算
上三角贡献通过反向状态得到:
最后
首尾不存在的项直接省略。若块边长不超过
这里对块数
乘积仍可用生成元表示
设两个矩阵分别记为
因此
原文 Algorithm 7.2 不仅给出秩界,还通过两次扫描直接计算乘积生成元。两个核心累积量为
为写出其余生成元,定义
则
新的上下生成元阶数分别是两组旧阶数之和,因此运算后可以再用前面的约简算法压缩。这里
求解方程:把消元压缩为小状态更新
原文 §7.3 还指出:对按相同方块划分的可逆矩阵,逆矩阵具有相同的上下切口秩;其显式求逆算法进一步要求矩阵强正则,即所有领先主块可逆。
下面把该节的 Schur 更新整理成便于理解的 LU 求解形式;新增的横线记号只用于区分消元后的生成元。假定对角块均为方块,初始化
对
最后
这对应
因而求解
回代
首尾取空状态,省略不存在的项。实现中的
这些公式展示了普通消元与生成元算术的关系,但强正则性只保证代数上可进行,并不等价于数值稳定。接下来的 SSS 论文使用正交变换,正是为了把快速求解和稳定性结合起来。
不同的最小生成元只差状态坐标变换
若两组生成元都最小,则存在各切口上的可逆矩阵
代入
最小性不可省略:非最小表示中可能存在彼此不同的冗余状态,未必能通过可逆变换互相得到。对于算法实现,比较两组生成元的具体数值通常没有意义;应该比较它们表示的矩阵、切口秩和数值条件。
5. Chandrasekaran 等人(2005):在 SSS 生成元上完成稳定求解
Some Fast Algorithms for Sequentially Semiseparable Representations. [4]
主要工作与贡献
这篇论文把顺序状态表示推进到三个计算问题:一般 SSS 矩阵的最小二乘与最小范数解;当
与上一节的无主元 Schur 消元相比,本文主要使用正交或酉变换,并允许矩形、奇异的对角块。重点不是重新证明低秩表示存在,而是在保持结构的同时,处理求解中的数值稳定性与秩亏问题。[4, §§1–3]
SSS 表示与拟可分表示的对应
原文把
上标
| 拟可分论文 [3] | SSS 论文 [4] | 作用 |
|---|---|---|
| 对角块 | ||
| 向右传播的下三角作用 | ||
| 向左传播的上三角作用 |
因此两篇使用的基本代数结构是相同的。这里保留两套原文符号,避免把 [4] 的状态传递矩阵
算法一:最小二乘需要左右两种变换
给定
不能逐块见到方程就立即求解。例如
的解为
Step 1 用右侧酉变换消除下三角耦合。
下面沿用原文 §2.2 的最小生成元假设,使携带的下状态维数等于相应切口秩。若输入生成元有冗余,可以先做上一节的双向约简,或用秩揭示分解按实际秩调整携带状态。对角块允许秩亏,与生成元是否最小是两个不同的问题。
第一个块列的下方都以
则
把
到第
进行同样的局部 QL 分解,并将相应变换施加到当前对角块和相关生成元。这里
累计这些右变换,得到
其中
Step 2 用左侧酉变换分离残差与独立行。
将上一步的生成元记为
做带秩判断的正交行变换,使之成为
其中
- 系数全部为零,对应不可通过选取未知量消除的残差;
- 当前可以确定的独立方程,保留在最终系统中;
- 仍只通过未来变量起作用的方程,压成
继续传递。
同一个
其中
Step 3 在独立方程中取最小范数解,再恢复变量。
令
只需求
其中
自由零空间分量取零便得到最小范数。实现中仍对上三角 SSS 系统使用结构化正交求解;原文在此调用作者先前发表于 HiPC 2002 的 SSS 求解方法,没有重新列出最后一步的全部局部递推。上式说明其目标,而不是要求形成稠密分解。右变换、左变换可按依赖顺序交织执行,原文称为一次顺序扫描的“延迟计算”。[4, §2]
在块大小和状态维数受控时,每步只分解小矩阵,工作量对块数线性。论文的矩形问题成本分析还保留了总行数、总列数和块数,不能把极高或极宽的矩形问题都简写成同一个
算法二:直接求出 的 SSS 表示
这里假设
考虑当前第一个块行。若
再用右侧变换使
同时变换
这三个式子一次处理对角、下方状态和上方状态所代表的全部右端项,而不是按
局部剩余行的右端生成元随之更新。例如
这里
则辅助状态更新为
初始
若
面向后续块的生成元可写成
对
它与下一节 HSS 求解器有相同的局部动作:压缩、消去、合并、恢复。区别在于这里沿顺序链推进,HSS 则在树上合并兄弟节点。
算法三:运算后重新压缩 SSS 状态
结构化乘法和求解会增大生成元。本文 §3.7 用两次相反方向的 SVD 扫描恢复紧凑表示,与上一节的最小生成元约简相呼应。
这里只写上三角部分。定义切口基底
以及
上切口矩形为
正向扫描。 用截断 SVD 计算
左因子列正交,并更新
首步的
反向扫描。 对已完成正向处理的生成元重新使用无帽记号,由右向左计算
并更新
末端的
论文强调:若原表示的另一侧累积基底可能很大,第一遍应取零截断容差,只做不截断的正交化;待尺度得到控制后,第二遍再使用目标阈值。局部 SVD 的小残差本身并不保证最终整个矩阵具有同样大小的相对误差。
稳定性
最小二乘部分以正交变换构造算法,并通过实验比较其后向误差。对
这篇论文的实际推进是:在顺序低秩结构上,建立能处理秩亏、多个结构化右端项和秩增长的计算流程;其中的正交化和状态约简,也是后面 HSS 直接求解的重要基础。
6. Chandrasekaran–Gu–Pals(2006):把顺序消元变成树上的隐式 ULV
A Fast ULV Decomposition Solver for Hierarchically Semiseparable Representations. [5]
主要工作与贡献
这篇论文针对 HSS 表示给出快速矩阵向量乘法和一般非对称线性系统求解器,并说明如何从普通矩阵构造数值 HSS 表示。对于非奇异矩阵
其中
这里用
Step 1 明确 HSS 压缩的整个外部相互作用
设节点
选择
父节点
此处用节点字母简化原文的层号与节点号
实际保存的对象只有:
- 叶节点的对角块
和基底 ; - 父子之间的传递矩阵
; - 兄弟之间的耦合矩阵
。
内部节点的长基底不显式保存。一般非对称问题中,上下耦合不互为共轭转置,
例如有四个叶簇
远处相互作用沿树路径通过小矩阵传递,不必为这个块另存一份完整基底。
Step 2 两个方向的遍历完成乘法
计算
然后令
叶节点恢复
根节点没有外部相互作用,无需保存外部基底或计算
对于一个叶簇,它在每一层接收祖先的另一侧子树贡献;这些源区域恰好分割叶簇的补集,因此所有非对角作用既不遗漏,也不重复。[5, 原式(4)–(8)]
这个算法与
Step 3 压缩整个块行,找出可以先求解的变量
考虑当前活跃节点
上方有
接着选取酉矩阵
其中
现在,前
这就是可以先求解并消去的局部部分。[5, §4.1]
Step 4 更新全部剩余右端项,再递归缩小系统
必须把已求出的
并把每个
执行
只在本地减去
把已解坐标统一放在前面后,系统具有
所以剩余矩阵直接是
Step 5 无法继续压缩时,合并兄弟节点
若当前块大小已经等于保留维数,便没有可消去的零耦合方向。对兄弟节点
两个孩子合成一个新的活跃块,树少一层。新块的大小是两个孩子保留维数之和,而不再是原矩阵中这一大区域的全部变量数。[5, §4.2]
随后继续“压缩或合并”。到根部只剩一个小稠密系统,直接求解;恢复时把父节点解分给孩子,接回先前保存的
至此得到完整解。[5, §§4.1–4.3]
HSS 表示的构造
原文 §5 从完整非对角块行、块列的 SVD 出发,提取
父节点的外部块行由两个孩子的外部块行删去内部兄弟作用后合并得到。利用孩子已经压缩的系数更新父级分解,再得到传递矩阵。对精确嵌套的基底,
近似构造时,应在孩子提供的压缩空间中生成父基底并控制截断,不能各层独立选取近似基底后就假定它们自动精确嵌套。
从一般显式稠密矩阵构造低秩 HSS,原文给出小秩情形的
论文还给出两个特殊途径[5, §§5.1–5.2]:
- 平滑核矩阵: 用分片、分尺度的 Chebyshev 展开形成基底和传递;近奇异对角线继续细分,生成元本身还可具有稀疏结构。不能把整个相邻半区都无条件看作同一个固定秩展开。
- 稀疏矩阵: 用选择矩阵提取对节点外部仍有非零作用的行列,合并时删除只剩内部相互作用的坐标。合适的行列访问数据结构可使构造与非零数同阶,但这些接口基底未必已经最小,也未必具有很小的秩。
复杂度与稳定性
若各层秩均受
原文在统一秩
这里沿用原文该段的
同一节还分析秩随几何尺度增长的情形,二维和三维模型的分解成本可分别达到
稳定性也依赖生成元的表示方式。原文给出的有利条件包括
或更一般地要求路径传递的范数增长受控;正交构造可以提供相应的
这篇论文最值得掌握的机制是:先把对外耦合压到少量坐标,再消去其余坐标;保留部分向父节点汇合,重复同一过程。
7. Xia 等人(2010):一般 HSS 树与可复用的显式分解
Fast algorithms for hierarchically semiseparable matrices. [6]
本节沿用该文的实矩阵记号,
主要工作与贡献
这篇论文针对早期 HSS 算法的表示和实现作出进一步改进:
- 用后序单编号表示一般二叉 HSS 树,允许不等深的子树。
- 构造时复用已经提取的正交基底,减少后续压缩工作的尺寸。
- 给出普通 Cholesky 的 HSS 表示算法,以及对称正定情形下的显式 ULV/广义 HSS Cholesky 算法。
- 把后者的求解组织成前向和反向树遍历,使局部因子能够用于多个右端项。
这四点中,普通 Cholesky 的 HSS 形式与广义 HSS Cholesky 是不同算法。原文给前者二次成本,后者在小秩和合适叶尺度下达到线性成本。[6, §1.2、§§3–4]
Step 1 后序编号让数据依赖与程序执行一致
对四个叶节点,原文采用后序编号:叶
令节点
原文把兄弟之间的耦合简写为
这里“一般完整二叉树”指每个内部节点恰有两个孩子,不要求所有叶在同一层。若有
树不平衡本身不会破坏受控秩条件下的节点总量计数,但会影响遍历深度、并行性和误差传播。
Step 2 构造 HSS 时,不重复携带已经提取的正交基底
用
之后处理节点
右边最后一个矩阵具有正交行,所以不必带着它做下一次压缩;只需分解前面的较小系数矩阵。它们具有相同的非零奇异值,列空间的有效信息也相同。[6, §3.1]
当两个孩子完成后,父节点的外部块行可以在压缩坐标中拼成
父节点外部不再包含兄弟内部作用;被分离出来的内部系数形成
该算法仍需处理显式稠密输入,在固定小秩情形下为
普通 HSS Cholesky
以下主要考虑实对称正定矩阵
其中
原文 §4.1 先介绍普通 Cholesky。消去一个叶对角块,
对其兄弟,Schur 更新为
以及
更远节点也要沿树路径更新。例如在上述四叶后序编号下,记
则右侧子树中的更新包括
这与前面的
这套方法对于取得部分消元后的 Schur 补很有用,但不是下面的线性成本算法。
广义 HSS Cholesky
核心变化是:在每个节点中,先隔离对外无作用的坐标,只消去这部分,把少量保留坐标传给父节点。全部直接数据交换只发生在父子之间。
Step 1 用同一个正交变换处理行和列。
设当前活跃块
上方有
Step 2 只分解无外部耦合的主块。
按消去和保留坐标划分
计算
于是
这是原文部分 Cholesky 的对称写法,与原式(7)–(8)等价。
Step 3 合并两个孩子的保留坐标。
若父节点
父节点活跃维数为
重复上述步骤直到根节点
若某节点
显式因子
算法保存每个节点的
正交变换可以存成 Householder 向量,置换通常只需要索引和尺寸;大量右端项可以重复使用同一棵分解树。
因子求解:一次后序前代,一次逆后序回代
先解
保留
再解
父节点把解按两个孩子的保留维数分开。设节点
内部节点继续把
此处分块尺寸必须是
复杂度
设 HSS 秩不超过
保存局部因子及一次前代、回代的成本分别为
这正是原文 Theorem 4.1 和 §4.3 的结论。若不选取合适叶尺度,应保留叶分解和内部节点的分别计数,而不是直接套用
| 方面 | 2006 年隐式 ULV [5] | 2010 年广义 HSS Cholesky [6] |
|---|---|---|
| 矩阵范围 | 一般可逆、可非对称 | 本文详细推导针对对称正定 |
| 局部变换 | 左右分别选择酉变换 | 同一正交变换用于合同变换 |
| 消元后矩阵 | 列变换后保留子系统 | 显式更新小 Schur 补 |
| 主要组织 | 以层为单位压缩、合并 | 后序处理一般二叉树 |
| 右端更新 | 利用 HSS 乘法传播已解变量的作用 | 父子间传递小向量 |
| 因子复用 | 论文重点是隐式分解与求解 | 显式保存局部因子和分解树 |
2010 年的改进还利用了对称正定性,不能把它的常数下降描述为在完全相同的一般非对称问题上无条件优于前者。[6, §4.5]
稀疏 PDE 求解
论文的一个应用方向是:先对稀疏 PDE 矩阵使用 nested dissection 和 multifrontal 消元,再用 HSS 压缩其中产生的稠密前沿矩阵或 Schur 补。也就是说,HSS 算法处理的是稀疏消元中的稠密中间对象。[6, §5]
若这些对象足够低秩,就可用 HSS 进行快速直接分解;若秩增长较快,也可以较粗地截断,再把近似因子用作预条件器。原文的三维高对比系数算例正属于后者,并不证明任意三维 PDE 都可精确地在线性时间内求解。
这一篇完成的推进,是把 HSS 结构变成更便于程序组织和重复求解的数据结构:后序构造,小块消元,父子合并,因子保存,双向回代。
8. 六篇文献对比
| 文献 | 主要推进 | 算法中最重要的中间量 |
|---|---|---|
| Hackbusch 1999 [1] | 建立分块低秩格式中的近似加、乘、逆与分解 | 低秩因子与截断后的 Schur 补 |
| Hackbusch–Börm 2002 [2] | 从核插值直接构造跨块、跨层共享的基底 | 叶基底 |
| Eidelman–Gohberg 2005 [3] | 解释生成元的最小维数、相似性及规范化 | 每个切口的状态与双向累积基底 |
| Chandrasekaran 等 2005 [4] | 在 SSS 表示上进行最小二乘、结构化求解和约简 | 局部酉变换、携带状态与延迟更新 |
| Chandrasekaran–Gu–Pals 2006 [5] | 在 HSS 树上进行一般非对称隐式 ULV 求解 | 已消去坐标与保留的对外作用坐标 |
| Xia 等 2010 [6] | 一般树、复用式构造及 SPD 显式分解和求解 | 小 Schur 补、局部因子与父子传递 |
这并非一条“后者完全替代前者”的路线。几何
复杂度应分成构造、分解和应用。下面仅比较文献已讨论的任务的复杂度。
| 格式与算法 | 从一般显式稠密输入构造 | 已有表示的一次乘法 | 已有表示的分解/求解 |
|---|---|---|---|
| [1] 的两种 |
本文没有任意稠密输入的通用快速构造结论 | 固定秩下 |
固定秩下近似逆/LU 为 |
| [2] 的核插值 |
借助核点值与局部积分构造;另计求积、建树 | 本篇不提供一般求逆算法 | |
| [3,4] 的拟可分/SSS | 小秩和受控块大小下仍通常为二次工作 | 两次状态扫描,见下方尺寸公式 | 生成元上的局部消元、正交分解;固定块大小和秩时对块数线性 |
| [5] 的 HSS | 固定小秩下一般构造 |
分解 |
|
| [6] 的 HSS | 改进的后序构造仍为固定小秩下 |
普通 HSS Cholesky 为固定小秩下 |
对
若
对于树结构,若所有节点都保存
真实秩分布比一个统一最大秩更能反映实际成本。固定秩下的线性结果,是在这个求和中进一步使用节点数与叶大小关系得到的。
算法的三种误差要分别看:
- 离散化误差:连续积分方程或 PDE 被替换成有限维矩阵问题。
- 压缩误差:用低秩矩阵
代替原离散矩阵。插值、截断 SVD、生成元约简都在这里产生误差。 - 浮点求解误差:即使对
采用稳定算法,计算结果通常也对应一个邻近系统。若
则对原矩阵的残差为
所以,稳定分解 并不意味着压缩误差消失。令
在一致的诱导范数下,若 ,标准扰动分析给出
该式是对总扰动的直接估计:问题越病态,所需压缩精度通常越高。对称矩阵还要注意正定性;例如对称扰动满足 是保持正定的一个充分条件,单有“每个块压得很准”并不自动提供这个保证。
局部阈值不能直接当作全局误差。若互不重叠的块构成一个完整矩阵划分,那么
这种情况下,可以按块分配 Frobenius 误差预算。但嵌套基底上的一次截断可能同时影响多个块;随后经过传递、相加、Schur 更新,误差还会继续传播。它与只在互不重叠块上各做一次截断不同。
因此,一个可靠实现至少应区分:
- 基底是否已正交化,阈值对应哪个尺度;
- 阈值控制单个局部分解、一个矩阵块,还是整个算子;
- 当前采用绝对误差还是相对误差;
- 压缩后求解的残差,是相对于
还是原来的 。
对线性系统,常用的整体检查量是
它有助于判断原离散问题的求解效果,但不能单独替代对病态问题前向误差的分析。[4, §3.5;5, §4.6;6, §4.4]
与蝶形算法的区别
蝶形算法主要利用互补尺度上的低秩性:一侧簇变小,另一侧簇变大,二者的组合仍保持可压缩。这里的传统
9. 后续进展
前六篇说明了低秩结构存在时如何计算。后续研究进一步追问两个问题:能否只调用算子而不读取全部条目?能否为最终压缩结果提供更明确的整体或逐块误差保证?
下面两项工作用于说明这两条进展,不替代前面的六篇文献,也不作为截至 2026 年的完整研究综述。
- 自适应
乘法:Börm 的后续工作 [7] 处理 的 近似。算法先用专门的树结构表示精确乘积,再通过秩揭示分解生成自适应嵌套基底,并压缩成新的 矩阵。其目标包括为结果中的可接受块控制相对误差,成本为 。这项进展回应了早期截断算术的一个难点:若中间步骤不断截断,最终某个块的相对误差不容易直接追踪;先保留可分析的精确中间表示,再安排压缩,就能给误差预算更清楚的依据。[7] - 随机 HSS 构造:Amsel 等人的工作 [8] 只访问
和 的矩阵向量乘积,构造块秩至多为 的 HSS 近似 。其结论包括算 子 调 用 数 额 外 运 算
以及相对于同一 HSS 结构中最佳秩 近似的期望误差保证:
这里控制的是平方 Frobenius 误差的期望,并非每次随机运行都满足同一个确定性误差界;总时间还应加上实际算子调用的成本。这使早期“已有快速 HSS 表示之后如何求解”的问题,进一步延伸到“如何以少量算子访问得到有保证的表示”。[8]
Remark. 面对一个具体数值问题,可以按下面的顺序组织实现:
第一步,明确矩阵从哪里来。 如果有可快速求值的核和几何点集,插值
第二步,检查真正需要的低秩对象。
第三步,先实现表示的应用。 让生成元或嵌套基底的乘法与小规模稠密矩阵一致,再进入求解。这个步骤能较早发现传递方向、共轭转置、兄弟耦合和端点状态的错误。
第四步,把压缩与消元分开检查。 先验证在不截断时的代数恒等式,再检查容差变化对秩、表示误差和原矩阵残差的影响。对广义 HSS Cholesky,要用不等叶大小、不等秩的例子检查保留坐标,避免
第五步,按实际用途选择精度。 多右端项问题重视因子复用;直接求解需要足够小的总扰动;预条件应用则可以容许更粗的近似,通过迭代修正误差。固定秩、固定误差和固定迭代次数是不同要求,应分别记录。
这些算法的共同计算原则可以明确地表述为:找到每个区域与外部交换信息所需的少量坐标,在这些坐标中传递和更新;只在局部形成稠密小矩阵,并在运算产生冗余时重新压缩。层次树和顺序链负责复用信息,低秩基底负责减少信息量,正交分解与误差预算负责让这种压缩可靠地服务于数值计算。
- Title: Hierarchical Matrices
- Author: Gypsophila
- Created at : 2026-11-04 15:04:19
- Updated at : 2026-11-04 20:30:27
- Link: https://chenx.space/2026/11/04/HierarchicalMat/
- License: This work is licensed under CC BY-NC-SA 4.0.