Hierarchical Matrices

Gypsophila

REFERENCES

  1. A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices, by W. Hackbusch. Computing, 62:89–108, 1999.
  2. H²-matrix approximation of integral operators by interpolation, by W. Hackbusch and S. Börm. Applied Numerical Mathematics, 43:129–143, 2002.
  3. On generators of quasiseparable finite block matrices, by Y. Eidelman and I. Gohberg. Calcolo, 42:187–214, 2005.
  4. 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.
  5. 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.
  6. 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.

后续进展补充

  1. Adaptive multiplication of H²-matrices with block-relative error control, by S. Börm. Numerical Algorithms, 101:65–96, 2026.
  2. 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.

本文将上述文献按三条主线组织: 矩阵的分块与共享基底、拟可分/SSS 矩阵的顺序生成元、HSS 矩阵的树结构与快速直接求解。每篇分别介绍其贡献和可执行的计算步骤,再比较它们之间的联系。前六份研究覆盖 1999—2010 年,建立了今天许多层次矩阵算法的基本框架。文末另用两项已发表的后续工作说明近年来在误差控制和随机化构造方面的进展。

1. 层次矩阵核心思想

从一个低秩块开始

考虑 。直接保存全部条目需要 空间,一次矩阵向量乘法也需要 工作。层次矩阵试图利用的不是“许多条目等于零”,而是某些子块的行或列彼此近似相关。

设行指标集为 、列指标集为 。如果

那么这个块只需保存 个数,作用于向量时只计算

先把 个输入压缩成 个系数,再展开成 个输出,成本由 降到 。

这里的 表示共轭转置。后文按原文分别使用 、 或实矩阵的 ,含义均在相应章节说明。数值低秩还依赖容差:例如按谱范数取最小的 ,使

因此,“秩小”必须与精度、划分和问题规模一起讨论。整个矩阵可以满秩,只有适当选择的非对角块低秩。

为什么局部相互作用会低秩?这类低秩性的一个典型来源是积分算子。若目标簇 和源簇 相距较远,核 在对应区域内平滑,便可以用少量分离项近似:

离散化后就得到低秩矩阵块。邻近区域可能包含核的奇性,需要继续细分,最后用小稠密块直接保存 [1, §5;2, §§2–3]。另一个来源是消元。稀疏矩阵的逆、Schur 补通常是稠密的,但不同子区域之间的影响可能只经过少量接口变量。此时,条目虽然增加,跨区域相互作用的有效维数仍然很小。这也是 SSS、拟可分和 HSS 算法能够用于直接求解的原因之一[3, §§2–4;5, §1;6, §§1,5]。

低秩分解解决一个块,层次组织解决整个矩阵。 只找到几个低秩块还不够:需要决定哪些块一起存储、它们能否共享基底,以及运算之后怎样继续维持小秩。

分块树与顺序链

在树结构中,把指标集递归分成子簇。例如

相应的矩阵分成四块。符合低秩条件的块停止细分;其他块继续递归。这样,大范围但平滑的相互作用用大低秩块表示,局部复杂相互作用用小块表示。

在顺序结构中,把变量分为 个连续块,考察每一个切口:

这里的指标是块指标。若这些整个矩形子矩阵的秩都小,就可以用少量状态跨越切口传递信息,形成 SSS/拟可分表示。要求不是“每个小块 的秩小”;标量分块时,后一个条件对任何矩阵都成立,几乎不提供信息。[3, Lemma 3.1、Theorem 4.2;4, 原式(1)]

结构 被压缩的对象 保存的主要数据 关键计算组织
可接受分块中的低秩子块 各块各自的低秩因子、小稠密块 块递归与截断
可接受分块,并共享簇基底 叶基底、传递矩阵、耦合矩阵、近场块 向上压缩、块间作用、向下展开
拟可分/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]

主要工作与贡献

这篇论文提出 矩阵的基本框架:用层次块划分和块内低秩表示压缩稠密矩阵,并在这种格式中完成加法、乘法、求逆和 LU 分解。关键推进是,快速计算不再局限于 ,而是形成一套运算后仍返回压缩格式的近似矩阵算术。

原文具体研究两种格式。第一种非常简单,只递归细分对角块;第二种在对角线附近保留更细的分块,用于逼近积分算子。多维几何划分的一般理论留给后续工作,不能把后来完整的 矩阵理论全部归入这一篇。[1, §§1–2,5]

Step 1 用树和块划分定义存储格式

原文以 表示指标集,以 表示簇树,以 表示 的块划分,并定义

是对应块。原文的第一种模型把叶块细分到标量,因而“所有终止块的秩不超过 ”包括这些叶块;实际实现通常在小块处停止,直接采用稠密存储。[1, Definition 2.6、Footnote 5]

对于 ,第一种格式递归写为

其中 继续采用同样格式, 的秩不超过 。原文用 展开复杂度推导,一般固定 具有相同的 阶数。[1, 原式(5)–(8)、§2.5]

对于秩一块,

这里 为共轭转置。乘法可以直接在因子上完成:

这说明“大矩阵乘低秩矩阵”可以退化成少量矩阵向量乘法。

Step 2 用小矩阵分解完成截断

低秩格式通常不对加法封闭。若

则

秩可能由 增至 。原文以 表示截断回秩 的结果。[1, §2.4.2,原式(10)]

下面用 QR 加小 SVD 写出同一截断操作,便于直接实现。设待压缩块为 , 各有 列:

只保留前 个奇异值,得到

大维度只出现在瘦 QR 中,SVD 的尺寸至多为 。成本约为

而不是对 稠密块做完整 SVD。若正交分解和 SVD 精确,则

原文通过秩一因子的 Gram 矩阵说明这一思想;这里的 QR 写法是对该计算的等价整理。它也更清楚地说明如何按容差自适应选择保留秩。

Step 3 递归完成加法和乘法

加法按共同块划分逐块进行:稠密叶直接相加,低秩块先拼接因子,再截断。

乘法仍遵循普通分块公式:

区别在于执行方式。例如

第一项递归计算,第二项利用低秩因子计算;相加时再次截断。两个低秩块相乘满足

真正需要形成的中间矩阵只有 。非对角结果同样通过“作用于少量列、相加、再压缩”产生。[1, §3.4]

需要注意,截断改变了算术:先分别截断再相加,通常不同于相加后统一截断。因此,普通矩阵乘法的结合律不能不加分析地转移到固定秩近似算术中。

Step 4 用 Schur 补递归求逆

设 和 Schur 补均可逆,记

分块求逆公式为

原文算法依次执行:

  1. 递归计算 。
  2. 利用低秩因子计算 ,形成并截断近似 Schur 补 。
  3. 递归计算 。
  4. 把 代入上面的四个块,必要时再截断。

在第一种格式中, 本身低秩,因此 Schur 更新也是低秩更新。快速求逆的核心是让消元过程始终以层次块和低秩更新的形式进行。[1, §3.5,原式(17)]

这里要求递归主块可用;矩阵整体可逆并不自动保证每个指定主块可逆。原文公式也不构成对任意非正定矩阵均无需主元处理的保证。

例子:三对角矩阵的逆

把三对角矩阵按连续指标分成左右两半,跨越分界的 至多各有一个非零元,因此秩至多为一。消去左半后,

只修改右半矩阵边界处的一个对角元,Schur 补仍是三对角矩阵。递归下去,逆矩阵的对应非对角块也至多秩一。

因此,在递归消元可进行的条件下,第一种 格式能够精确表示并计算三对角矩阵的逆。这说明层次低秩并不总是近似性质;它也可能来自精确的代数结构。[1, Proposition 4.3]

第二种格式:在近场保留更细的块

相邻区间之间的核可能在接触点附近变化很快,把整个兄弟块立即压成低秩并不合适。原文为此引入邻接型 块:

只有靠近接触点的左下子块 继续采用 结构;上方非对角邻接块使用这种模式,下方邻接块使用转置的分块模式。这里转置指划分模式,不要求上下数值块相同。递归对角块与这两类邻接块组合成第二种 格式。[1, §5.1,原式(18)–(20)]

原文以对数核的分片常数配点离散为例。对实值问题可写成

这里把原文的对数核明确写成实值的 。在远场块中围绕 展开,令

于是

直接给出秩不超过 的因子。第二种分块保证每个远场块内有

原文由 Taylor 余项得到该离散模型的误差界

这里是矩阵诱导无穷范数;近场小块直接计算。该结论对应这一特定核和分块,不是所有 近似的通用误差界。[1, §5.3,Proposition 5.2]

复杂度与适用边界

对原文研究的平衡递归格式,固定块秩时有以下结论。[1, §§3.1–3.6,5.2]

运算 原文的渐近阶数
存储、矩阵向量乘法
截断加法
截断矩阵乘法、近似求逆、近似 LU

这是对给定格式的算术成本,而不是对任意矩阵都能以固定秩达到指定精度的承诺。求逆时若所需秩持续增长,或截断误差被病态问题放大,实际成本和求解精度都要重新评估。

这篇论文建立了完整的基本循环:分块低秩表示,执行普通代数公式,再把结果截断回压缩格式。 下一篇进一步减少了不同块、不同层之间重复保存的基底。

3. Hackbusch–Börm(2002):用插值构造嵌套的 基底

H²-matrix approximation of integral operators by interpolation. [2]

主要工作与贡献

这篇论文给出一种可直接实现的 构造:用多项式插值近似积分核,用插值多项式的父子关系产生传递矩阵。用户主要需要提供核函数的点值和叶层积分,而不必实现任意阶核导数。

其主要结果是:在合适的块稀疏性、叶大小和核光滑性条件下,构造可分析的积分算子近似,并把存储和矩阵向量乘法降到 。论文重点是构造与乘法,并没有在这里给出一般 求逆的完整算术。[2, §§1,3–5]

Step 1 从 Galerkin 矩阵和几何划分出发

考虑定义在 上的积分方程。这里 是嵌入空间维数,例如圆边界虽然是一维曲线,其包围盒和张量插值位于 ,故取 :

采用局部支撑基函数 ,离散方程为

其中

质量矩阵 通常稀疏,困难在于稠密的 。

对簇 ,取包含全部相关支撑的轴平行包围盒 。原文用

判断块 是否可接受。可接受则停止划分;否则继续细分,直到得到可直接计算的叶块。[2, 原式(4)、Table 1]

控制精度与块数的折中:更严格的分离条件通常让单块近似更容易,但产生更多块。它本身并不能替代核的光滑性假设。

Step 2 用核的点值构造低秩耦合

在 内取张量积插值点 ,对应 Lagrange 多项式记为 。在可接受块上,

定义

代入积分并分离变量,得到

只取决于簇 ,可以服务于所有与它相互作用的列簇; 则描述特定两簇之间的耦合。独立低秩因子由此变为共享的簇基底。[2, 原式(8)–(11)]

Step 3 利用插值的再现性得到精确嵌套

设 是 的子簇,各簇采用同一个多项式空间。父簇上的多项式可以被子簇插值精确再现:

定义传递矩阵

把上式与 积分,就得到

因此只需要显式保存叶簇的 ,内部节点只保存 。这里的嵌套关系来自多项式再现性,本身没有增加插值误差;核近似误差来自 Step 2。[2, 原式(12)–(14)]

Remark: 多项式次数、节点数和秩的记号

原文主体用 表示张量积基底数,附录又用次数为 、每方向 个 Chebyshev 点的记号。实现时应明确区分这两个量。若每方向取 个点,则每方向次数至多 ,实际基底数为 ;若次数取 ,则 。两种写法的渐近阶同为 ,但数组尺寸不能混用。

例如,附录的 Chebyshev 根节点为

先把这些点仿射映射到各坐标区间,再取张量积即可。

Step 4 三次遍历完成矩阵向量乘法

设要计算 。先令输出向量 ,并把所有目标系数置零。以下公式分别对应原文的 ForwardTransformation、Multiplication 和 BackwardTransformation。[2, §4.1,Tables 2–3]

第一次:从叶到根压缩源向量。

叶簇计算

内部节点利用已计算的子簇系数:

第二次:在可接受块之间传播作用。

所有目标系数初始为零,随后逐个可接受块累加

第三次:从根到叶展开。

将父簇贡献传给每个子簇:

到叶簇后计算

最后,把不可接受叶块的直接乘法

计入结果。近场块不能因为远场采用了插值而被省略。

三次遍历各有明确作用:把源数据变成系数,在小系数空间中处理相互作用,再把目标系数变回物理变量。父子间传递保证同一份信息不会在每层从头计算。

复杂度

令所有簇的基底维数为 。原文要求叶簇大小与 可比:

并要求每个簇在块划分中参与的块数由常数 控制。于是叶数为 ,二叉树节点数也是 。每个传递或耦合矩阵的计算量是 ,因此

叶基底和近场块也满足同一阶数。[2, 原式(15)–(17)、§4.1]

这个线性结论不是把 乘法当成 ,而是节点数量随叶大小相应减少。若叶大小固定为 ,直接计数应保留 项。

构造时,耦合矩阵只需核点值;叶基底和近场块还需数值积分。原文分别计入了插值、求积阶数和建树成本;典型几何建树本身可需要 ,不能把整个输入预处理都无条件称为线性。[2, §§2,4.2]

误差分析和变阶扩展

原文假设核渐近光滑,即离开 后有

配合几何可接受性和 Chebyshev 插值,得到

若导数常数具有原文 Remark 5.2 给出的阶乘型增长,并满足相应的 ,误差随阶数快速下降。积分方程解的误差还需要连续问题的椭圆性、离散空间稳定性和求积误差控制。[2, §5]

Remark: 插值误差常数包含什么?

原文 Lemma 5.1 中的常数为

其中阶乘分母来自插值余项,对数因子来自张量 Chebyshev 插值的稳定性,导数常数来自核的渐近光滑性。这里按原文在 的阶数范围理解。若

则阶乘增长被抵消,得到由 控制的几何衰减及较缓的对数因子。因此,只有几何分离而没有导数常数的控制,还不足以推出给定的收敛率。[2, 原式(19)–(20)、Remark 5.2]

论文 §6.3 还讨论变阶策略:小簇用较低阶,大簇用较高阶。如果子簇次数低于父簇次数,父多项式在子簇上的再插值便只是近似,嵌套关系会引入额外误差。适当地控制阶数随层级的增长,可使所有层的总成本保持 ;这与前面固定次数时的精确多项式嵌套应分开理解。

其计数可以这样理解:距叶层高度为 的节点约有 个,若每节点工作量只随 多项式增长,则

大多数细簇保持低阶,少数粗簇才使用高阶。原文单位圆单层势实验采用分片常数 Galerkin、 和 。从 到 个自由度,规模增加 倍,构造与乘法时间分别由 、 秒增至 、 秒。这些是原实现对近线性增长的实验支持,并非当前硬件上的性能预测。[2, §7,Table 4]

对于配点法、双层势等变体,主要变化落在叶基底的定义上,核点值耦合和多项式传递结构仍可复用。[2, §6]

这篇论文把“嵌套基底”从抽象结构变成了可计算的公式:叶上做积分,块间取核值,层间评价插值多项式。

4. Eidelman–Gohberg(2005):用最少的生成元描述跨切口相互作用

On generators of quasiseparable finite block matrices. [3]

主要工作与贡献

这篇论文主要是一篇带有新结果的理论综述。它系统回答三个问题:拟可分矩阵怎样由生成元表示,生成元的最小维数是多少,同一矩阵的不同生成元之间有什么关系;同时给出生成元约简、构造、乘法和求逆算法。原文特别指出,局部生成元的最小性判据 Theorem 4.5 是其中的新结果之一。[3, §§1,4]

它的重要性在于,快速算法的成本由生成元维数控制。仅知道矩阵存在低秩结构还不够,还要能识别并移除表示中的冗余状态。

本节沿用原文: 为矩阵, 为块数, 的大小为 ,整个矩阵的行列总数相等。这里的 与后文 SSS/HSS 论文中的总矩阵阶数不是同一含义。

Step 1 用输入、传递和输出生成元表示非对角块

下三角部分写为

其中

上三角部分类似:

其中

相邻块之间的空乘积为单位矩阵,故 。对角块独立保存为 。[3, §2,原式(3)–(4)]

例如,从第 个输入块到第 个输出块, 把输入变成低维状态, 逐站传递, 把状态变成输出。各切口允许具有不同状态维数 ,不必强制所有秩相同。

Step 2 由切口秩理解最小状态维数

原文定义累积生成元

于是每个下三角切口矩形都有

所以它的秩不会超过 。反过来,Theorem 4.2 证明:若

则一定可以找到 的生成元。也就是说,最小状态维数恰好等于跨越相应切口的相互作用秩,而不是某个人为选定的统一参数。[3, Lemma 3.1、Theorem 4.2]

从系统角度看, 说明过去的输入能产生哪些状态, 说明这些状态能影响哪些未来输出。若有状态不能被输入产生,或产生后不能影响任何输出,它就是冗余的。

原文 Theorem 4.5 把这种全局最小性化为小矩阵的局部条件。对一组阶数为 的候选生成元,最小性等价于以下满秩条件:

检验这些条件不需要形成所有大切口子矩阵。

Step 3 两次扫描,把已有生成元约简到最小

假设已有 ,但维数可能偏大。原文 Theorem 6.1 给出两阶段算法。

第一次,从左向右,使累积输入基底正交。

先做识别实际秩的紧致 LQ 或 SVD 分解

对 ,分解小矩阵

并按列拆成

同时更新

这一步删除输入端的冗余,并使新的 具有正交行。

第二次,从右向左,使累积输出基底正交。

先分解

对 ,计算

并拆成

最后令

若所有分解仅删除精确零奇异值,则所得生成元表示同一个矩阵,并达到最小阶数。若按容差删除小奇异值,则得到近似约简,误差需要结合传递过程分析,不能把局部阈值直接当成整个矩阵的相对误差。[3, §6]

局部正交关系还可写为

或

它们分别对应原文的右正规形和左正规形。正交化不仅为降秩服务,也能减少状态坐标选择引起的数值放大。

Remark: 如果输入是稠密矩阵,怎样构造生成元?

Theorem 6.4 不要求先有生成元。算法先压缩第一列下方的块:

然后在第 步拆分

把已经压缩的历史与新的一列拼接:

最后取 。每一步只携带旧切口的基底,不重复分解整个切口矩阵。在块大小和秩有界时,总成本为 ;从显式稠密输入开始,读取全部数据本身就可能需要二次工作。[3, Theorem 6.4]

Algorithm 6.5 将这里的分解替换为截断 SVD,允许指定最大生成元阶数或奇异值阈值。实现时必须把截断后的右因子按旧状态列与新输入列拆成 ,不能把整块右因子只当成 。

Step 4 两个方向的状态递归完成乘法

计算 时,下三角贡献通过正向状态得到:

上三角贡献通过反向状态得到:

最后

首尾不存在的项直接省略。若块边长不超过 、生成元阶数不超过 ,原文给出的成本界为

这里对块数 线性,但不能忽略 项。[3, Algorithm 7.1,原式(77)]

乘积仍可用生成元表示

设两个矩阵分别记为 。固定一个切口,把它们各分成 块,则

因此

原文 Algorithm 7.2 不仅给出秩界,还通过两次扫描直接计算乘积生成元。两个核心累积量为

为相容的空零矩阵,端点不存在的生成元按零处理。相应的对角块为

为写出其余生成元,定义

则

新的上下生成元阶数分别是两组旧阶数之和,因此运算后可以再用前面的约简算法压缩。这里 的两块按维数与乘积恒等式排列,修正了原文 Algorithm 7.2 中两块倒置的印误。固定块大小和秩时,整个乘法仍对块数线性。[3, Algorithm 7.2]

求解方程:把消元压缩为小状态更新

原文 §7.3 还指出:对按相同方块划分的可逆矩阵,逆矩阵具有相同的上下切口秩;其显式求逆算法进一步要求矩阵强正则,即所有领先主块可逆。

下面把该节的 Schur 更新整理成便于理解的 LU 求解形式;新增的横线记号只用于区分消元后的生成元。假定对角块均为方块,初始化

对 ,

最后

这对应 ,其中 为单位下三角, 的对角块为 ,非对角块满足

因而求解 不必形成 。前代 使用

回代 使用

首尾取空状态,省略不存在的项。实现中的 表示小线性系统求解,不需要显式形成小逆矩阵。

这些公式展示了普通消元与生成元算术的关系,但强正则性只保证代数上可进行,并不等价于数值稳定。接下来的 SSS 论文使用正交变换,正是为了把快速求解和稳定性结合起来。

不同的最小生成元只差状态坐标变换

若两组生成元都最小,则存在各切口上的可逆矩阵 ,使

代入 后,相邻的 和 逐项抵消,所以矩阵不变。这是原文 Theorem 5.4 的相似性结论。

最小性不可省略:非最小表示中可能存在彼此不同的冗余状态,未必能通过可逆变换互相得到。对于算法实现,比较两组生成元的具体数值通常没有意义;应该比较它们表示的矩阵、切口秩和数值条件。

5. Chandrasekaran 等人(2005):在 SSS 生成元上完成稳定求解

Some Fast Algorithms for Sequentially Semiseparable Representations. [4]

主要工作与贡献

这篇论文把顺序状态表示推进到三个计算问题:一般 SSS 矩阵的最小二乘与最小范数解;当 都以 SSS 形式给出时,直接构造 的结构化解 ;以及对运算后膨胀的生成元进行模型约简。

与上一节的无主元 Schur 消元相比,本文主要使用正交或酉变换,并允许矩形、奇异的对角块。重点不是重新证明低秩表示存在,而是在保持结构的同时,处理求解中的数值稳定性与秩亏问题。[4, §§1–3]

SSS 表示与拟可分表示的对应

原文把 矩阵分成 个块,,表示为

上标 为共轭转置。对矩形问题,可令 为 ,总行数与总列数不必相等。除尺寸相容外,对角块不要求可逆。[4, 原式(1)、§2]

拟可分论文 [3] SSS 论文 [4] 作用
对角块
向右传播的下三角作用
向左传播的上三角作用

因此两篇使用的基本代数结构是相同的。这里保留两套原文符号,避免把 [4] 的状态传递矩阵 与 [3] 的整个矩阵 混为一谈。

算法一:最小二乘需要左右两种变换

给定 ,目标是找到所有最小二乘解中范数最小的一个:

不能逐块见到方程就立即求解。例如

的解为 ,第一条方程不能在看见第二条之前独立决定 。原文因此先用右变换整理变量,再用左变换分离独立方程和不可消除的残差。[4, §2]

Step 1 用右侧酉变换消除下三角耦合。

下面沿用原文 §2.2 的最小生成元假设,使携带的下状态维数等于相应切口秩。若输入生成元有冗余,可以先做上一节的双向约简,或用秩揭示分解按实际秩调整携带状态。对角块允许秩亏,与生成元是否最小是两个不同的问题。

第一个块列的下方都以 结尾。做 QL 型分解

则

把 右乘到第一个块列,就能把部分列与所有后续行的耦合同时消为零。尚未处理完的列保留为小状态。

到第 步,把上一步保留的状态与当前变量合起来,对

进行同样的局部 QL 分解,并将相应变换施加到当前对角块和相关生成元。这里 是辅助因子,不是原来的 。[4, §§2.1–2.2]

累计这些右变换,得到

其中 为块上三角 SSS 矩阵。实际保存的是局部反射和新的生成元,不形成大稠密 。

Step 2 用左侧酉变换分离残差与独立行。

将上一步的生成元记为 。前一步尚未确定作用的方程用小矩阵 传递。第 步对

做带秩判断的正交行变换,使之成为

其中 和 都取满行秩,允许其中一组为空。第一步没有 ,只需处理 ;最后一步没有未来状态。三类行分别表示:

  • 系数全部为零,对应不可通过选取未知量消除的残差;
  • 当前可以确定的独立方程,保留在最终系统中;
  • 仍只通过未来变量起作用的方程,压成 继续传递。

同一个 必须作用于右端项。经过全部步骤,并把零行集中排列,可写成

其中 满行秩。[4, §2.3]

Step 3 在独立方程中取最小范数解,再恢复变量。

令 ,则

只需求 的最小范数解,最小残差为 。为看清“最小范数”的含义,概念上可继续作

其中 酉、 可逆下三角,于是

自由零空间分量取零便得到最小范数。实现中仍对上三角 SSS 系统使用结构化正交求解;原文在此调用作者先前发表于 HiPC 2002 的 SSS 求解方法,没有重新列出最后一步的全部局部递推。上式说明其目标,而不是要求形成稠密分解。右变换、左变换可按依赖顺序交织执行,原文称为一次顺序扫描的“延迟计算”。[4, §2]

在块大小和状态维数受控时,每步只分解小矩阵,工作量对块数线性。论文的矩形问题成本分析还保留了总行数、总列数和块数,不能把极高或极宽的矩形问题都简写成同一个 。[4, §2.4]

算法二:直接求出 的 SSS 表示

这里假设 为可逆方阵, 和 可以是矩形矩阵,并按行与 共形分块。 本身具有 SSS 表示,算法输出 的生成元,不逐个生成全部 。这个可逆性假设不要求原始对角块 各自可逆。[4, §3]

考虑当前第一个块行。若 的列数小于该块大小,先构造酉变换

再用右侧变换使

同时变换 的第一个块行,以及未知量的相应变量 。上方被隔离的方程只需解小三角系统,其生成元满足

这三个式子一次处理对角、下方状态和上方状态所代表的全部右端项,而不是按 的每一列分别求解。[4, §3.1]

局部剩余行的右端生成元随之更新。例如

这里 是同一个首对角块变换后的上下两部分; 不表示原来的第二个对角块。 也作同样形式的更新。已解变量对更远块行的影响都沿 的下三角生成元传播,原文把它们汇总为辅助状态 。若

则辅助状态更新为

初始 。保存这个小状态,便不必每消去一个块就更新整片尚未处理的 。[4, 原式(11)、§3.1]

若 的列数不小于当前块大小,块内没有可先释放的变量方向,就合并 的前两个行块和列块。与此同时, 与 只合并前两个行块,保留原有列划分,并把活动对角位置向右移一块。对 ,省略参数 后,新的局部对角块为

面向后续块的生成元可写成

对 ,原来的首对角块此时进入扩展的下状态,同时继续携带延迟修正 。随后继续压缩和消元,直到只剩一个小系统;递归返回后拆开 的行块,将活动对角位置左移,并逆向应用 恢复生成元。[4, §§3.2–3.3,3.6]

它与下一节 HSS 求解器有相同的局部动作:压缩、消去、合并、恢复。区别在于这里沿顺序链推进,HSS 则在树上合并兄弟节点。

算法三:运算后重新压缩 SSS 状态

结构化乘法和求解会增大生成元。本文 §3.7 用两次相反方向的 SVD 扫描恢复紧凑表示,与上一节的最小生成元约简相呼应。

这里只写上三角部分。定义切口基底

以及

上切口矩形为 。若这两个因子具有任意缩放,仅对其中一个截断可能严重放大误差。因此先正交化一侧,再从另一侧判断有效秩。

正向扫描。 用截断 SVD 计算

左因子列正交,并更新

首步的 为空。这样得到的累积列基底 列正交。

反向扫描。 对已完成正向处理的生成元重新使用无帽记号,由右向左计算

并更新

末端的 为空。下三角部分对 采用同类处理。[4, §3.7]

论文强调:若原表示的另一侧累积基底可能很大,第一遍应取零截断容差,只做不截断的正交化;待尺度得到控制后,第二遍再使用目标阈值。局部 SVD 的小残差本身并不保证最终整个矩阵具有同样大小的相对误差。

稳定性

最小二乘部分以正交变换构造算法,并通过实验比较其后向误差。对 的多右端结构化求解,原文 §3.5 明确指出,全面误差分析超出本文范围,主要提供归一化残差实验。因此,不能把整篇所有算法概括成“已证明对任意生成元都后向稳定”。

这篇论文的实际推进是:在顺序低秩结构上,建立能处理秩亏、多个结构化右端项和秩增长的计算流程;其中的正交化和状态约简,也是后面 HSS 直接求解的重要基础。

6. Chandrasekaran–Gu–Pals(2006):把顺序消元变成树上的隐式 ULV

A Fast ULV Decomposition Solver for Hierarchically Semiseparable Representations. [5]

主要工作与贡献

这篇论文针对 HSS 表示给出快速矩阵向量乘法和一般非对称线性系统求解器,并说明如何从普通矩阵构造数值 HSS 表示。对于非奇异矩阵 ,求解器隐式进行

其中 为酉矩阵, 为按消元次序排列的下三角矩阵。它把 SSS 的“压缩—消去—合并”放到近完全二叉树上,使多个互不依赖的块可以在同一层处理。[5, §§2–4]

这里用 表示全局因子,区别于 HSS 生成元 。所谓“隐式”,是说算法通过局部变换完成分解和求解,不需要先构造完整全局因子。

Step 1 明确 HSS 压缩的整个外部相互作用

设节点 对应指标集 ,补集为 。HSS 低秩条件针对

选择 张成第一个矩阵的列空间, 张成第二个矩阵共轭转置的列空间。于是对兄弟节点 ,

父节点 的基底满足

此处用节点字母简化原文的层号与节点号 等; 对应原文 。本文依据父子关系写出递推,修正原文 §2 个别下标与 排版不一致之处。[5, §2;6, 原式(1)]

实际保存的对象只有:

  • 叶节点的对角块 和基底 ;
  • 父子之间的传递矩阵 ;
  • 兄弟之间的耦合矩阵 。

内部节点的长基底不显式保存。一般非对称问题中,上下耦合不互为共轭转置, 的列数也可以不同。

例如有四个叶簇 ,左右父簇分别为 ,则一个跨父簇的块为

远处相互作用沿树路径通过小矩阵传递,不必为这个块另存一份完整基底。

Step 2 两个方向的遍历完成乘法

计算 时,向上计算输入系数

(叶节点)

(内部节点)

然后令 ,从根向下递推

叶节点恢复

根节点没有外部相互作用,无需保存外部基底或计算 ;从根到两个孩子时,直接取 、。

对于一个叶簇,它在每一层接收祖先的另一侧子树贡献;这些源区域恰好分割叶簇的补集,因此所有非对角作用既不遗漏,也不重复。[5, 原式(4)–(8)]

这个算法与 的三阶段乘法具有相同的信息复用思想。区别在于 HSS 的跨域作用由兄弟耦合组织,而几何 的块划分通常包含另外保留的近场关系。

Step 3 压缩整个块行,找出可以先求解的变量

考虑当前活跃节点 ,其块大小为 ,行基底 有 列。用 QL 型分解选择酉矩阵 ,使

上方有 行。由于 同时表示该节点对整个外部的作用,这些行对所有外部变量的系数都变为零。

接着选取酉矩阵 ,对当前对角块做右变换:

其中 为下三角。这一步可通过上方行块的 LQ 分解完成。整块列同时被变换,所以要更新

现在,前 条方程独立成为

这就是可以先求解并消去的局部部分。[5, §4.1]

Step 4 更新全部剩余右端项,再递归缩小系统

必须把已求出的 对所有未解方程的作用扣除。令

并把每个 补成

执行

仍具有 HSS 表示,只需更新当前节点的 ,高层传递和耦合关系不变。因此这次更新可以用前面的快速乘法完成,并利用已知零块节省工作。实现时还应复用跨层传递中已经计算的状态,才能保持所述成本。[5, §4.1;6, §4.5]

只在本地减去 不够:别的节点已经解出的变量也可能作用于当前剩余行。原文的 HSS 右端更新正是处理这种全局作用。

把已解坐标统一放在前面后,系统具有

所以剩余矩阵直接是 ,不需要另做一次普通 Schur 补矩阵更新;影响已经进入 。这一点与下一节保持对称性的算法不同。

Step 5 无法继续压缩时,合并兄弟节点

若当前块大小已经等于保留维数,便没有可消去的零耦合方向。对兄弟节点 形成

两个孩子合成一个新的活跃块,树少一层。新块的大小是两个孩子保留维数之和,而不再是原矩阵中这一大区域的全部变量数。[5, §4.2]

随后继续“压缩或合并”。到根部只剩一个小稠密系统,直接求解;恢复时把父节点解分给孩子,接回先前保存的 ,再计算

至此得到完整解。[5, §§4.1–4.3]

HSS 表示的构造

原文 §5 从完整非对角块行、块列的 SVD 出发,提取 。若基底列正交,则兄弟耦合可以写成

父节点的外部块行由两个孩子的外部块行删去内部兄弟作用后合并得到。利用孩子已经压缩的系数更新父级分解,再得到传递矩阵。对精确嵌套的基底,

近似构造时,应在孩子提供的压缩空间中生成父基底并控制截断,不能各层独立选取近似基底后就假定它们自动精确嵌套。

从一般显式稠密矩阵构造低秩 HSS,原文给出小秩情形的 工作量。所谓线性表示空间不包括原始稠密输入的 存储。[5, §5]

论文还给出两个特殊途径[5, §§5.1–5.2]:

  • 平滑核矩阵: 用分片、分尺度的 Chebyshev 展开形成基底和传递;近奇异对角线继续细分,生成元本身还可具有稀疏结构。不能把整个相邻半区都无条件看作同一个固定秩展开。
  • 稀疏矩阵: 用选择矩阵提取对节点外部仍有非零作用的行列,合并时删除只剩内部相互作用的坐标。合适的行列访问数据结构可使构造与非零数同阶,但这些接口基底未必已经最小,也未必具有很小的秩。

复杂度与稳定性

若各层秩均受 控制,叶大小与 同阶,则

存储与每次应用分解每个右端项求解

原文在统一秩 、叶大小 、 个右端项的模型下给出总成本上界

这里沿用原文该段的 表示秩;并不把它与前面插值次数混同。[5, §4.4]

同一节还分析秩随几何尺度增长的情形,二维和三维模型的分解成本可分别达到 、。因此,线性求解器依赖低秩结构,不能仅由问题来自 PDE 或积分方程就推出线性成本。

稳定性也依赖生成元的表示方式。原文给出的有利条件包括

或更一般地要求路径传递的范数增长受控;正交构造可以提供相应的 -范数控制。论文给出稳定性论断和数值实验,并明确将完整证明留待别处。[5, §4.6]

这篇论文最值得掌握的机制是:先把对外耦合压到少量坐标,再消去其余坐标;保留部分向父节点汇合,重复同一过程。

7. Xia 等人(2010):一般 HSS 树与可复用的显式分解

Fast algorithms for hierarchically semiseparable matrices. [6]

本节沿用该文的实矩阵记号, 表示转置。复数对应形式应使用共轭转置;下面的分解算法按原文详细讨论对称正定情形。

主要工作与贡献

这篇论文针对早期 HSS 算法的表示和实现作出进一步改进:

  1. 用后序单编号表示一般二叉 HSS 树,允许不等深的子树。
  2. 构造时复用已经提取的正交基底,减少后续压缩工作的尺寸。
  3. 给出普通 Cholesky 的 HSS 表示算法,以及对称正定情形下的显式 ULV/广义 HSS Cholesky 算法。
  4. 把后者的求解组织成前向和反向树遍历,使局部因子能够用于多个右端项。

这四点中,普通 Cholesky 的 HSS 形式与广义 HSS Cholesky 是不同算法。原文给前者二次成本,后者在小秩和合适叶尺度下达到线性成本。[6, §1.2、§§3–4]

Step 1 后序编号让数据依赖与程序执行一致

对四个叶节点,原文采用后序编号:叶 合成父节点 ,叶 合成父节点 ,最后合成根 。每个父节点总在两个孩子之后处理。

令节点 的孩子为 ,则

原文把兄弟之间的耦合简写为 ;例如 表示从节点 到节点 的耦合。

这里“一般完整二叉树”指每个内部节点恰有两个孩子,不要求所有叶在同一层。若有 个叶,节点总数为 。后序处理可以用栈实现:孩子完成后把小结果压栈,父节点弹出两个孩子的结果并合并。[6, §2,Theorem 2.1]

树不平衡本身不会破坏受控秩条件下的节点总量计数,但会影响遍历深度、并行性和误差传播。

Step 2 构造 HSS 时,不重复携带已经提取的正交基底

用 表示节点 的完整非对角块行。处理节点 时,计算

之后处理节点 ,其块行的一部分已写成 ,因此

右边最后一个矩阵具有正交行,所以不必带着它做下一次压缩;只需分解前面的较小系数矩阵。它们具有相同的非零奇异值,列空间的有效信息也相同。[6, §3.1]

当两个孩子完成后,父节点的外部块行可以在压缩坐标中拼成

父节点外部不再包含兄弟内部作用;被分离出来的内部系数形成 。块列对称地得到 。对全部节点按后序重复,就得到 Algorithm 1。[6, §3.2]

该算法仍需处理显式稠密输入,在固定小秩情形下为 。它改进的是重复压缩量、常数和数据访问方式,不是无条件绕过稠密输入的数据量。若 HSS 块在稀疏消元过程中逐步产生,则有机会避免先形成完整稠密矩阵。

普通 HSS Cholesky

以下主要考虑实对称正定矩阵 。对称构造可以选取满足以下关系的生成元:

其中 为兄弟节点。

原文 §4.1 先介绍普通 Cholesky。消去一个叶对角块,

对其兄弟,Schur 更新为

以及

更远节点也要沿树路径更新。例如在上述四叶后序编号下,记

则右侧子树中的更新包括

这与前面的 一起展示了 三类生成元的变化。结果可以继续用 HSS 形式表示,但连续消元会反复访问尚未处理的节点,原文给出的总成本为固定小秩下 。[6, §4.1]

这套方法对于取得部分消元后的 Schur 补很有用,但不是下面的线性成本算法。

广义 HSS Cholesky

核心变化是:在每个节点中,先隔离对外无作用的坐标,只消去这部分,把少量保留坐标传给父节点。全部直接数据交换只发生在父子之间。

Step 1 用同一个正交变换处理行和列。

设当前活跃块 为 , 为 。作

上方有 个坐标。由于对称性,它们的外部行耦合与列耦合同时变为零,合同变换还保持正定性。[6, 原式(5)–(6)]

Step 2 只分解无外部耦合的主块。

按消去和保留坐标划分

计算

于是

这是原文部分 Cholesky 的对称写法,与原式(7)–(8)等价。 是保留坐标上的小 Schur 补。正定性保证可消去主块和相应 Schur 补仍正定。

Step 3 合并两个孩子的保留坐标。

若父节点 的两个孩子已处理完,构造

父节点活跃维数为 。这里的 是消元后系统的生成元,不是原矩阵中父区域的完整大块。[6, 原式(9)]

重复上述步骤直到根节点 ,最后做

若某节点 ,就没有可消去坐标。此时取 、、,直接把它保留给父节点,跳过空的 运算。

显式因子

算法保存每个节点的 以及合并坐标所需的置换 。它们组成

是一系列局部正交变换、块三角矩阵和置换的递归乘积,不必是按原始变量顺序展开的普通下三角矩阵。显式是指这些可复用的局部因子被保存,而不是形成一个稠密 因子。[6, 原式(10)–(11)]

正交变换可以存成 Householder 向量,置换通常只需要索引和尺寸;大量右端项可以重复使用同一棵分解树。

因子求解:一次后序前代,一次逆后序回代

先解 。节点 收到当前右端向量 后,计算

保留 ,把 传给父节点;父节点拼接两个孩子的保留右端。到根部求解

再解 。根部先求

父节点把解按两个孩子的保留维数分开。设节点 收到的保留解为 ,则

内部节点继续把 分给两个孩子,叶节点得到原始变量的解。[6, §4.3,原式(14)–(16)]

此处分块尺寸必须是 ,不能误用已经消去的维数 。原文回代段的该处尺寸有排版问题;只有恰好取 时,两者数值才偶然相同。

复杂度

设 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 补、局部因子与父子传递

这并非一条“后者完全替代前者”的路线。几何 容许在近场保留很多局部关系;HSS 用更强的完整块行/块列低秩条件换取简洁的树上消元;SSS/拟可分则利用一个顺序上所有切口的小状态结构。合适的格式取决于实际秩分布与输入方式。

复杂度应分成构造、分解和应用。下面仅比较文献已讨论的任务的复杂度。 为矩阵总阶数;HSS 和固定阶 的叶大小取为与秩 同阶。

格式与算法 从一般显式稠密输入构造 已有表示的一次乘法 已有表示的分解/求解
[1] 的两种 格式 本文没有任意稠密输入的通用快速构造结论 固定秩下 固定秩下近似逆/LU 为
[2] 的核插值 借助核点值与局部积分构造;另计求积、建树 本篇不提供一般求逆算法
[3,4] 的拟可分/SSS 小秩和受控块大小下仍通常为二次工作 两次状态扫描,见下方尺寸公式 生成元上的局部消元、正交分解;固定块大小和秩时对块数线性
[5] 的 HSS 固定小秩下一般构造 分解 ,每右端项
[6] 的 HSS 改进的后序构造仍为固定小秩下 普通 HSS Cholesky 为固定小秩下 ;广义分解 ,每右端项

对 个大小约为 的 SSS 块,总阶数约为 。存储和乘法直接按生成元尺寸计数为

若 ,则是 ;若采用标量块而让 增长,通常仍有 项。这两种计数没有矛盾,区别在于块大小是否随状态维数一起变化。

对于树结构,若所有节点都保存 小矩阵,直接计数还应包含

或

真实秩分布比一个统一最大秩更能反映实际成本。固定秩下的线性结果,是在这个求和中进一步使用节点数与叶大小关系得到的。

算法的三种误差要分别看:

  • 离散化误差:连续积分方程或 PDE 被替换成有限维矩阵问题。
  • 压缩误差:用低秩矩阵 代替原离散矩阵。插值、截断 SVD、生成元约简都在这里产生误差。
  • 浮点求解误差:即使对 采用稳定算法,计算结果通常也对应一个邻近系统。若

    则对原矩阵的残差为

    所以,稳定分解 并不意味着压缩误差消失。令

    在一致的诱导范数下,若 ,标准扰动分析给出

    该式是对总扰动的直接估计:问题越病态,所需压缩精度通常越高。对称矩阵还要注意正定性;例如对称扰动满足 是保持正定的一个充分条件,单有“每个块压得很准”并不自动提供这个保证。

局部阈值不能直接当作全局误差。若互不重叠的块构成一个完整矩阵划分,那么

这种情况下,可以按块分配 Frobenius 误差预算。但嵌套基底上的一次截断可能同时影响多个块;随后经过传递、相加、Schur 更新,误差还会继续传播。它与只在互不重叠块上各做一次截断不同。

因此,一个可靠实现至少应区分:

  • 基底是否已正交化,阈值对应哪个尺度;
  • 阈值控制单个局部分解、一个矩阵块,还是整个算子;
  • 当前采用绝对误差还是相对误差;
  • 压缩后求解的残差,是相对于 还是原来的 。

对线性系统,常用的整体检查量是

它有助于判断原离散问题的求解效果,但不能单独替代对病态问题前向误差的分析。[4, §3.5;5, §4.6;6, §4.4]

与蝶形算法的区别

蝶形算法主要利用互补尺度上的低秩性:一侧簇变小,另一侧簇变大,二者的组合仍保持可压缩。这里的传统 通常利用几何分离后的平滑性,SSS/HSS 则要求规定的切口或完整外部块具有小秩。因此,对高频振荡问题,传统远场块或 HSS 块的秩可能随频率增长。此时应检查实际秩,判断是否需要方向性表示或互补低秩表示,而不是仅仅换一棵树就假定固定秩复杂度仍然成立。

9. 后续进展

前六篇说明了低秩结构存在时如何计算。后续研究进一步追问两个问题:能否只调用算子而不读取全部条目?能否为最终压缩结果提供更明确的整体或逐块误差保证?

下面两项工作用于说明这两条进展,不替代前面的六篇文献,也不作为截至 2026 年的完整研究综述。

  • 自适应 乘法:Börm 的后续工作 [7] 处理 的 近似。算法先用专门的树结构表示精确乘积,再通过秩揭示分解生成自适应嵌套基底,并压缩成新的 矩阵。其目标包括为结果中的可接受块控制相对误差,成本为。这项进展回应了早期截断算术的一个难点:若中间步骤不断截断,最终某个块的相对误差不容易直接追踪;先保留可分析的精确中间表示,再安排压缩,就能给误差预算更清楚的依据。[7]
  • 随机 HSS 构造:Amsel 等人的工作 [8] 只访问 和 的矩阵向量乘积,构造块秩至多为 的 HSS 近似 。其结论包括
    算子调用数额外运算
    以及相对于同一 HSS 结构中最佳秩 近似的期望误差保证:

    这里控制的是平方 Frobenius 误差的期望,并非每次随机运行都满足同一个确定性误差界;总时间还应加上实际算子调用的成本。这使早期“已有快速 HSS 表示之后如何求解”的问题,进一步延伸到“如何以少量算子访问得到有保证的表示”。[8]

Remark. 面对一个具体数值问题,可以按下面的顺序组织实现:

第一步,明确矩阵从哪里来。 如果有可快速求值的核和几何点集,插值 的基底可直接构造;如果只有显式稠密矩阵,要接受构造过程可能需要二次工作;如果只能调用 ,则需要与这种访问方式相配的随机化构造。

第二步,检查真正需要的低秩对象。 检查可接受块;拟可分/SSS 检查每个顺序切口;HSS 检查每个树节点的整个外部块行和块列。抽查几个孤立小块不足以确认格式适用。

第三步,先实现表示的应用。 让生成元或嵌套基底的乘法与小规模稠密矩阵一致,再进入求解。这个步骤能较早发现传递方向、共轭转置、兄弟耦合和端点状态的错误。

第四步,把压缩与消元分开检查。 先验证在不截断时的代数恒等式,再检查容差变化对秩、表示误差和原矩阵残差的影响。对广义 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.
Comments