[科普中国]-稀疏矩阵算法

基本概念

在科学与工程领域中进行数值计算时经常需要处理大型的稀疏矩阵，这些矩阵的特征一是阶数很大，二是含有大量的零元，使得无法使用传统的稠密矩阵存储方案和算法。由于矩阵阶数太大，若存储整个矩阵，将会占用大量的存储空间，甚至会超出计算系统的存储容量，因此需要设计专门的稀疏存储方案，尽可能的减少零元素的存储。

稀疏矩阵算法是以稀疏矩阵作为核心数据结构的算法。与稠密矩阵算法相比，稀疏矩阵算法的最大特点是通过只存储和处理非零元素从而大幅度降低存储空间需求以及计算复杂度，代价则是必须使用专门的稀疏矩阵压缩存储数据结构，因而在计算过程中引入了大量的离散间接寻址操作。稀疏矩阵算法是典型的不规则算法，计算访存比很低，并且计算过程中的访存轨迹与稀疏矩阵的稀疏结构相关，很难发掘计算过程中的时空局部性，因此在传统的基于Cache的处理器上稀疏矩阵算法的计算效率很低。为了提高稀疏矩阵算法的计算效率，需要从稀疏存储数据结构和稀疏矩阵算法两方面对现有算法进行改进。
按照应用领域的不同，稀疏矩阵算法分为两类，一类是非数值计算算法，典型代表是图搜索算法，包括宽度优先搜索等核心算法；另一类则是数值计算算法，典型代表是稀疏线性方程组求解，包括稀疏矩阵向量乘、稀疏三角方程组求解以及稀疏矩阵分解等核心算法。

算法图遍历算法图遍历算法是最基本的图算法之一，由指定节点开始，按照一定规则遍历图结构中所有的连通节点，包括宽度优先搜索(Breadth First Search,BFS)和深度优先搜索(Depth First Search, DFS)等核心算法。
作为最基本的图遍历算法，宽度优先搜索算法代表了图遍历算法的计算特性，具有非常重要的研究意义。一方面，BFS算法是最短路径、邻接查询、可达性查询等算法的实现基础，广泛应用于图分割、信念传播统计以及网络社区结构发现等领域;另一方面，BFS算法作为典型的数据密集型算法，体现了数据密集型应用对高性能计算系统的需求，广泛应用于大规模并行计算系统的数据处理能力评测，已经成为Parboil, Rodinia和Graph500等基准测试程序的核心算法。
在实际应用中，图的规模在不断增大，相应的，对图的存储和处理开销不断增加，有效地实现大规模并行BFS算法具有十分重要的意义。

稀疏线性方程组求解法稀疏线性方程组的求解是对自然科学和社会科学中许多实际问题进行数值模拟时的关键技术之一。在高层建筑、桥梁、水坝、防洪堤的结构设计中，需对变形与应力情况进行模拟；在油气资源探测与分析、数值天气预报、飞行器的动力学分析中，需利用流体力学方程组进行模拟;在进行恒星大气分析与核爆实验时，常需利用辐射流体力学与粒子统计平衡等规律进行模拟。在对这些问题进行分析模拟时，通常利用偏微分方程建立数学模型。在对偏微分方程的离散求解过程中，稀疏线性方程组求解算法扮演着十分重要的角色。在许多不以偏微分方程建模的问题中，稀疏线性方程组求解同样发挥了重要的作用。在空中交通控制、电力线路中的最优电流问题中，需利用数学规划求解;在对采纳某项政策时在某给定条件下对国内、国际多个区域的相应经济指标进行预测时，需利用CGE模型进行分析；在可靠性分析、排队网络分析与计算机系统性能评估中，常利用具有大量状态的离散Markov链进行模拟。在这些问题的求解中，稀疏线性方程组的求解都占有重要位置，并且往往是整个计算过程中的性能瓶颈，稀疏线性方程组的高效求解是计算数学和工程应用中十分重要的课题之一。

解稀疏线性方程组的方法包括直接法(direct method)与迭代(iterative method)两类。直接法指在不考虑计算舍入误差的情况下，通过包括矩阵分解和三角方程组求解等有限步的操作求得方程组的精确解，因此又称精确法;迭代法指给定一个初始解向量，通过一定的计算构造一个向量列(一般通过逐次迭代得到一系列逼近精确值的近似解)，向量列的极限为方程组理论上的精确解。迭代法对存储空间的需求低，在求解高阶非病态稀疏线性方程组方面具有一定优势。然而，迭代法不适合求解病态问题，性能因问题而异，并且面临精度控制、收敛速度慢或不收敛等问题。与迭代法相比，直接法的通用性好，求解结果精度高，性能稳定。当矩阵分解结果能够被多次后续计算重用以及多右端项时，直接法的优势尤其明显。在有限元分析、模拟电路瞬态仿真等应用领域的商用软件均采用直接法求解器作为标准的稀疏线性方程组求解器。但直接法的缺点在于对存储资源要求较高，无法处理高阶稀疏矩阵。

一般来说，迭代法的求解速度高于直接法。但是，如果使用直接法时矩阵分解过程能够被很多后续计算重复使用，则后续的三角阵求解可以非常快速实现，此时直接法在性能上具有优势。典型例子是模拟电路瞬态仿真，这时需要多次以Newton-Raphson方法求解非线性方程，每一次求解均会在工作点附近展开为线性方程，而且所有线性方程的矩阵分解方式都是固定的，因此目前求解该类问题最好的方法是直接法。稀疏矩阵的矩阵分解在GPU上的实现是很困难的，主要难点在于现有算法的数据依赖性导致可利用的并行性不足。此外，矩阵元素的排列顺序对计算过程中间结果矩阵的非零元素个数有很大影响，同时矩阵分解后的非零元素的分布与原来矩阵可能很不相同。

迭代法的理论基础相对复杂，并且具有多种不同的具体算法，但其基本形式均为从一个猜测解出发，通过多次迭代逐渐收敛，当误差满足一定条件时迭代中止。共扼梯度法(CG)是迭代法的主流方法之一，特别适合于特征值为良态分布的对称正定方程组；其它迭代法包括Jacobi、逐次超松弛(SOR)、广义极小剩余(GMRES)、预条件共扼梯度(PCG)等。迭代法的核心算法是稀疏矩阵向量乘(SpMV)，因此实现SpMV的高效并行结构也是实现迭代法的基础。

直接法由高斯消元法发展向来，求解过程包括矩阵排序(matrix ordering)符号分解(symbolic factorization)、数值分解(numerical factorization)、三角方程组求解((triangular solves)四个步骤。其中，矩阵排序和符号分解属于预处理部分。矩阵排序通过启发算法置换稀疏矩阵的行列，试图在后续计算中维持矩阵的稀疏性或数值稳定性。符号分解则是预先对矩阵分解后的稀疏结构进行预测，预先分配存储空间并记录数据相关性。直接法的计算瓶颈在于数值分解部分和三角方程组求解部分，高效的直接法求解依赖于二者的高效实现。

对于一个稀疏线性方程组是选择直接法还是迭代法求解，一般有如下原则：对于低阶矩阵或大型带状矩阵所对应的线性方程组，用直接法求解；而对于大型(非带形)矩阵所对应的线性方程组，用迭代方法求解。实际上，选用何种方法还要看具体的应用背景，比如，对于线性规划和一些结构工程应用，只有直接法是切实可行的。对于精度要求很高的问题，还可以采用由直接法得到初始解再用迭代法进行迭代的方法求解，这种方法称为迭代精化法。

挑战面向算法加速器的稀疏矩阵算法并行化方法由于算法加速器采用了与通用处理器完全不同的体系结构和编程模型，因此需要专门设计针对算法加速器的并行程序。与通用处理器相比，算法加速器能够运行更多的线程，因此一方面需要将计算任务划分为更多的子任务，实现细粒度并行；另一方面需要考虑线程间的负载均衡和数据局部性，特别是线程间的通信和同步开销。与通用处理器相比，算法加速器更容易受到不规则访存和计算的影响，因此必须通过矩阵预处理等手段提高计算的规则性和数据的局部性。

通用处理器与算法加速器协同计算策略异构体系结构的计算单元包括通用处理器和算法加速器两部分。在学术研究中通常分别针对这两种计算单元设计并行算法，但在实际应用中，需要同时利用计算系统的所有计算资源以获取最高性能。然而，现有的针对通用处理器和算法加速器的并行算法分别采用了不同的编程模型和优化策略，无法直接互相迁移，并且在协同计算过程中，不但需要在两个平台上都能高效运行的并行程序，而且需要实现通信与计算的重叠以隐藏数据传输的开销，这些都要求对现有的算法针对异构平台进行重新设计。

面向稀疏结构特征的数据结构和算法映射不同应用领域的稀疏矩阵往往具有不同的稀疏结构特征，这些特征因问题而异；具有不同稀疏结构特征的稀疏矩阵往往对应不同的最优存储方案以及不同的并行算法实现。在计算过程中，对于稀疏矩阵算法的不同计算阶段，会因稀疏结构导致的计算特征的不同而适合在通用处理器或算法加速器上运行或是由二者协同计算。因此，研究矩阵的稀疏结构并设计相应的数据结构和算法映射策略，能够有效提高算法对于不同矩阵和不同计算平台的适应性。1