69、HyGrid：高效的CPU - GPU混合网格化算法

HyGrid：高效的CPU - GPU混合网格化算法

在处理不规则且大规模的天文数据时，传统的基于卷积的网格化方法在GPU上难以取得良好的性能。为解决这一问题，一种名为HyGrid的CPU - GPU混合网格化算法应运而生。

相关工作回顾

以往与基于卷积的网格化相关的研究主要集中在三个方面：
1. 基于GPU的散射（Scatters） ：由于存在数据竞争问题，提出的散射方法主要致力于避免同时更新设备内存。例如，Van Amesfoort等人通过为每个输入点分配一个小的私有网格来实现GPU上的散射，虽能从根本上避免数据竞争，但限制了输出网格的分辨率，不适用于最新的大型射电望远镜。Humphreys和Cornwell基于内存复制，为澳大利亚平方公里阵列探路者（ASKAP）在GPU上实现了散射。Romein设计了另一种基于GPU的散射方法，通过让线程尽可能长时间地在寄存器中累积加权和，减少设备内存访问。Merry进一步改进了该算法，使用线程粗化显著降低了设备内存的寻址开销。不过，这些方法严重依赖干涉仪数据的空间相干性，无法应用于单碟望远镜。
2. 预排序和收集（Pre - ordering and Gathers） ：提出的收集方法主要旨在减少非贡献输入点的不必要内存访问。Edgar等人使用Thrust对输入点进行预排序，虽提高了设备内存访问效率，但粗分箱策略仍导致大量不必要的内存寻址开销和卷积计算。Gai等人采用紧凑分箱方法对磁共振成像数据进行预排序，实现了GPU上的收集。Winkel等人提出基于HEALPix的紧凑分箱方法，在CPU上实现了名为cygrid的算法，但该算法在处理大规模数据时，数据布局和预排序策略会显著增加时间开销和内存消耗。
3. 优化策略 ：研究还关注提高GPU效率以及实现CPU和GPU之间的协同工作。Amesfoort的散射方法将卷积核作为一维纹理访问，并将公共数据仅加载到共享内存一次，有效利用了内存带宽，但不适用于大核或可变核。Merry采用了线程粗化和并行效率之间的权衡方法，Gai等人通过在GPU和CPU之间均匀划分网格化算子来改善负载平衡。

系统模型

预备知识

在应用FITS世界坐标系（WCS）时，HEALPix适合对球形采样空间进行分箱。给定网格分辨率（用$N_{side}$表示，$N_{side} \in {1, 2, 4, 8, \cdots}$），HEALPix将球面分层细分为$12N_{side}^2$个等面积像素。所有像素中心位于$4N_{side} - 1$个恒定纬度的环上，且在每个环上方位角等距。HEALPix软件库支持球面变换，可实现像素中心的WCS坐标与像素索引之间的相互转换。

问题描述

假设有$N$个输入点分布在采样平面上，每个输入点具有坐标和采样值。将平面划分为$I \times J$个规则网格单元后，基于卷积的网格化将为每个网格单元分配重采样值和总权重值。权重函数通常采用径向对称函数，如高斯核，用于集中每个贡献者的部分贡献到核的中心。

HyGrid算法

为限制搜索空间，HyGrid算法首先使用多CPU对输入点进行预排序，然后在GPU上为每个输出点分配一个专用线程来并行化卷积，并通过数据布局和线程粗化进一步优化网格化过程。

预排序

• 基于环的查找策略 ：为每个输入点找到其所属的像素和环，构建一个基于环的查找表$R_{start}$，记录每个环上的最低索引像素。通过该查找表，可以高效地找到给定像素在核半径内的贡献者。
• 多CPU上的预排序算法 ：HEALPix - 基于的预排序算法（HpxPreOrdering）采用并行块间接排序算法对输入点进行重排。该算法分为排序和生成查找表两个步骤，复杂度取决于块间接排序算法，最坏情况下为$O(N \log(N))$，最好情况下为$O(N)$。

Algorithm 1. HEALPix-based Pre-ordering Algorithm: HpxPreOrdering
Parameter:
Nside
input:
S, Spix, Sval
output:
S′, Spix, Rstart, r1, rl, rh
Initialize:
S′ ←∅, Rstart ←∅
▷sort input points
1: BlockIndirectSort(Spix, Sval, N)
2: for n = 1 to N do
3:
S′[n] ←S[Ssix[n]]
4: end for
▷generate lookup table
5: r1 = Pix2Ring(Spix[1], Nside)
6: rtotal = Pix2Ring(Spix[N], Nside)
7: Rstart[0] ←1, idx = 1
8: for cntRix = r1 + 1 to rtotal do
9:
pix = Ring2Start(cntRix, Nside)
10:
cntIdx = BiSearchLastPosLessThan(Spix, idx, N, pix)
11:
Rstart[cntRix −r1] = cntIdx + 1
12:
idx = cntIdx
13: end for

GPU上的收集算法

GPU - 基于的网格化收集算法（GPUGridGather）由单个线程执行每个输出点的计算。线程通过迭代相邻环，计算每个环上的最低和最高索引相邻像素，然后在$S’$中顺序扫描并应用卷积公式，最终得到输出点的重采样值和总权重值。

Algorithm 2. GPU-based Gridding Gather: GPUGridGather
Parameter:
Nside, rad
input:
S′, Spix, r1, Rstart, r1, rl, rh
output:
Wi,j, V [gi,j]
Initialize:
Wi,j = 0, V [gi,j] = 0
1: for cntRix = rl to rh do
▷iterate each neighboring ring
2:
pix = Ring2Start(cntRix, Nside)
3:
(α, β) ←Pix2Loc(pix, Nside)
4:
tanRad = Rad2Tan(βi,j, β, rad)
5:
pl = Loc2Pix(αi,j −tanRad, β, Nside)
▷lowest pixel
6:
ph = Loc2Pix(αi,j + tanRad, β, Nside)
▷highest pixel
7:
lIdx = Rstart[cntRix −r1]
8:
rIdx = Rstart[cntRix −rl + 1]
9:
n = BiSearchLastPosLessThan(Spix, lIdx, rIdx, pl) + 1
▷calculate Equation 1
10:
while Spix[n] ≤ph do
11:
if d(αi,j, βi,j; αn, βn) ≤rad then
12:
Wi,j = Wi,j + w(αi,j, βi,j; αn, βn)
13:
V [gi,j] = V [gi,j] + w(αi,j, βi,j; αn, βn) ∗V [s′
n]
14:
end if
15:
n = n + 1
16:
end while
17: end for

CUDA优化策略

在CUDA 8.0中实现HyGrid时，提出了两种优化策略：
1. 数据布局优化 ：将$S$和$G$以数组格式存储，将查找表作为一维纹理内存访问，将全局参数作为常量内存访问，并将线程的部分卷积结果保留在寄存器中，减少全局内存的刷新次数。
2. 线程粗化 ：在输出分辨率较大的情况下，由于GPU处理器的处理能力有限，无法同时执行所有线程。通过将位于同一环上的多个连续输出点的计算合并到一个线程中，减少线程数量。线程粗化的性能取决于粗化因子$\gamma$，即每个线程计算的输出点数。

实验环境

品牌	型号	核心数	SMs	Tflops	RAM (GB)	频率 (GHz)	缓存 (MB)	每核心线程数
GPU	NVIDIA Tesla K40	2880	15	1.4	12	-	-	-
CPU	Intel Xeon E5 - 2620	6	-	-	-	2.5	15	12

基准测试

通过一系列实验评估HyGrid在不同输入规模和输出分辨率下的性能，并与cygrid进行比较。
1. 性能与输入规模的关系 ：在采样场为$5^{\circ} \times 5^{\circ}$，输出分辨率$Or = 600’‘$的情况下，HyGrid在预排序、收集和整体网格化过程中均显著优于多线程的cygrid。随着输入点数$N$的增加，cygrid的预排序策略受内存墙限制，而HyGrid的GPUGridGather能获得更高的性能和数据吞吐量。
2. 性能与线程粗化的关系 ：选择固定输入点数和不同输出分辨率的基准模型进行测试。结果表明，线程粗化策略能显著提高性能，尤其是在输出分辨率较大时。但需要根据输出分辨率和硬件情况合理选择粗化因子$\gamma$，避免因卷积核形状问题导致输出图像质量下降。
3. 性能与数据有序程度的关系 ：实际的射电天文观测数据可能部分有序。实验表明，对于随机输入数据，预排序过程消耗的时间较多；而对于有序输入数据，最耗时的过程是收集阶段。因此，原始数据的有序程度对CPU和GPU之间的负载平衡有很大影响。
4. 实际应用性能 ：使用HyGrid（$\gamma = 3$）将1.7 GB的有序或随机数据网格化到$900 \times 900$的网格图像上，时间开销分别在8秒和32秒以内。处理1 TB的数据时，有序数据约需1.34小时，随机数据约需5.35小时。因此，对于小规模数据，单个工作站即可胜任；对于大规模数据，可使用小型集群同时运行多个HyGrid实例。为实现性能和节能的最佳平衡，应根据输入数据的有序程度合理选择CPU和GPU的数量。

综上所述，HyGrid算法通过引入基于HEALPix的高效查找表和合理利用数据布局及线程粗化，实现了GPU加速的收集过程，在处理不规则且大规模的天文数据时表现出显著的性能优势。

HyGrid：高效的CPU - GPU混合网格化算法

实验结果分析

为了更直观地展示HyGrid算法的性能，下面结合具体的实验数据和图表进行详细分析。

性能与输入规模的关系

在采样场为$5^{\circ} \times 5^{\circ}$，输出分辨率$Or = 600’‘$（即$I \times J = 90 \times 90$），粗化因子$\gamma = 1$的条件下，对HyGrid和cygrid进行性能对比。实验结果表明，HyGrid在各个方面都展现出了显著的优势。

算法	预排序加速比	收集加速比	整体网格化加速比
HyGrid vs cygrid	125.81倍	18.63倍	71.25倍

从数据中可以看出，随着输入点数$N$的增加，cygrid的预排序策略受到内存墙的限制，导致性能下降。而HyGrid的GPUGridGather算法能够充分利用GPU的并行计算能力，获得更高的性能和数据吞吐量。特别是当$N > 10^8$时，HyGrid的优势更加明显。下面是该实验结果的mermaid流程图：

graph LR
    A[输入点数N增加] --> B[cygrid预排序受内存墙限制]
    B --> C[cygrid性能下降]
    A --> D[HyGrid GPUGridGather发挥优势]
    D --> E[HyGrid性能提升]
    C --> F[HyGrid优势明显]
    E --> F

性能与线程粗化的关系

选择固定输入点数$N = 1.5 \times 10^8$，不同输出分辨率的基准模型进行测试。实验结果显示，线程粗化策略能够显著提高HyGrid的性能。

粗化因子$\gamma$	输出分辨率$Or$	$I \times J$	加速比
2	30’‘	600 × 600	2.28倍
3	20’‘	900 × 900	4.64倍

然而，粗化因子$\gamma$的选择需要根据输出分辨率和硬件情况进行合理调整。当$\gamma$过大时，卷积核的形状会变成矩形，可能导致输出图像质量下降。例如，在对$327×327$（$Or = 55’‘$）输出网格细胞进行网格化时，$\gamma = 2$是一个较好的选择；而对于$514×514$（$Or = 35’‘$）的输出网格细胞，$\gamma = 3$更为合适。下面是线程粗化性能影响的mermaid流程图：

graph LR
    A[选择粗化因子$\gamma$] --> B[输出分辨率和硬件情况合适]
    B --> C[性能提升]
    A --> D[粗化因子$\gamma$过大]
    D --> E[卷积核形状改变]
    E --> F[输出图像质量下降]

性能与数据有序程度的关系

实际的射电天文观测数据可能部分有序。为了研究数据有序程度对HyGrid性能的影响，分别对随机输入数据和有序输入数据进行实验。实验结果表明，原始数据的有序程度对CPU和GPU之间的负载平衡有很大影响。

数据类型	预排序与收集时间开销比
随机输入数据	预排序时间 > 收集时间
有序输入数据	预排序时间 < 收集时间

对于随机输入数据，预排序过程消耗的时间较多；而对于有序输入数据，最耗时的过程是收集阶段。这说明在实际应用中，需要根据数据的有序程度合理分配CPU和GPU的任务，以实现更好的负载平衡。下面是数据有序程度影响性能的mermaid流程图：

graph LR
    A[数据类型] --> B[随机输入数据]
    A --> C[有序输入数据]
    B --> D[预排序时间长]
    C --> E[收集时间长]
    D --> F[CPU和GPU负载平衡受影响]
    E --> F

实际应用性能

使用HyGrid（$\gamma = 3$）将1.7 GB的有序或随机数据网格化到$900 \times 900$的网格图像上，时间开销分别在8秒和32秒以内。处理1 TB的数据时，有序数据约需1.34小时，随机数据约需5.35小时。

数据规模	数据类型	时间开销
1.7 GB	有序数据	8秒以内
1.7 GB	随机数据	32秒以内
1 TB	有序数据	约1.34小时
1 TB	随机数据	约5.35小时

根据实验结果，对于小规模数据，单个工作站即可胜任；对于大规模数据，可使用小型集群同时运行多个HyGrid实例。为实现性能和节能的最佳平衡，应根据输入数据的有序程度合理选择CPU和GPU的数量。

总结与展望

HyGrid算法通过引入基于HEALPix的高效查找表，结合数据布局和线程粗化等优化策略，实现了GPU加速的收集过程。在处理不规则且大规模的天文数据时，HyGrid在性能上显著优于传统的cygrid算法。

未来的研究可以进一步探索以下几个方向：
1. 优化粗化因子选择 ：开发更智能的算法，根据不同的输入数据和硬件环境自动选择最优的粗化因子$\gamma$，以进一步提高性能。
2. 多GPU并行计算 ：研究如何将HyGrid算法扩展到多个GPU上，实现更高效的并行计算，处理更大规模的数据。
3. 应用领域拓展 ：尝试将HyGrid算法应用到其他领域，如医学成像、地理信息系统等，探索其在不同场景下的性能和适用性。

总之，HyGrid算法为处理不规则且大规模的数据提供了一种高效的解决方案，具有广阔的应用前景和研究价值。