3010在TensorCore上通过WY表示进行快速对称特征值分解0张少帅，鲁琪{szhang36,rkshah5}@uh.edu休斯顿大学计算机科学系美国德克萨斯州0大和田浩之，横田理央{ootomo.h@rio.,rioyokota@}gsic.titech.ac.jp 东京工业大学全球科学信息与计算中心日本东京0吴盼若 pwu7@uh.edu休斯顿大学计算机科学系美国德克萨斯州0摘要0对称特征值分解（EVD）是许多科学领域中使用的一种基本的分析和数值工具。从性能上看，最先进的算法通常是两阶段三对角化方法。两阶段三对角化的第一阶段称为连续带约化（SBR），它将对称矩阵约化为带状形式，其计算成本通常占主导地位。当使用TensorCore（专用矩阵计算加速器）来加速昂贵的EVD时，传统的基于ZY表示的方法由于矩阵计算的不利形状而导致性能不佳。在本文中，我们提出了一种新的方法，使用WY表示而不是ZY表示（详见第3.2节），它可以更好地结合局部性和并行性，从而在TensorCore上表现更好。实验结果表明，与最先进的实现相比，所提出的方法在SBR中可以带来高达3.7倍的加速，在整个EVD中可以带来2.3倍的加速。0CCS概念 • 计算数学 → 数学软件性能； • 计算理论 →并行算法； • 计算方法学 → 并行算法。0关键词 特征值分解，GPGPU，数值线性代数，TensorCore，HPC，混合精度计算，矩阵计算，奇异值分解，低秩逼近0未经ACM许可，可以免费制作本作品的全部或部分数字或印刷副本，但不得为了牟利或商业优势而制作或分发副本，并且副本必须带有本通知和第一页的完整引用。必须尊重ACM以外其他人拥有的本作品组成部分的版权。允许带有信用的摘要。要复制其他内容，或重新发布，或发布到服务器上，或重新分发到列表上，需要事先获得特定的许可和/或支付费用。请向permissions@acm.org申请权限。PPoPP'23，2023年2月25日至3月1日，加拿大蒙特利尔，QC，©2023年计算机协会。ACM ISBN 979-8-4007-0015-6/23/02...$15.00https://doi.org/10.1145/3572848.357751601 引言0随着处理器架构的进步，NvidiaVolta、Turing和Ampere架构的专用单元TensorCore被设计用于执行高效的半精度（16位浮点格式）矩阵乘法，这些乘法在单精度中累积。最近，TensorCore已被广泛用于通过加速GEMM（通用矩阵乘法）来训练大规模深度神经网络。研究人员还尝试在TensorCore上部署线性代数算法的可能性。这些算法包括LU分解[20]，QR分解[41]和一些BLAS3操作[42]。特征值分解是数值线性代数中另一个重要的矩阵计算，因为它在各种应用中非常有用，包括量子化学[33]、量子力学[15]、量子物理学[32]、数据驱动分析以及众多机器学习和信号处理任务。其中一些应用可能需要高精度算术，如双精度，但在许多新兴的数据驱动方法中，单精度甚至更低精度已经足够。例如，主成分分析[1]、低秩逼近[29]、深度学习中的二阶优化[17]以及相关应用[35,43]。TensorCore通常可以用于加速这些算法和应用。对称特征值问题可以定义为：0� = � × Λ × � − 10其中 � 是对称矩阵，� 是包含特征向量的正交矩阵，Λ是一个对角矩阵，其对角线上有特征值。计算对称稠密矩阵的特征值最关键的步骤是三对角化。三对角化步骤可以通过著名的Householder三对角化[13]可靠地计算。另外，LAPACK中的块变体涉及面板分解和双侧矩阵尾部更新[2]，可以显著改善局部性。并行实现可以在[18]和[24]中找到。Bischof和Sun[7]提出了更详细的方法，进一步改进了分层存储系统的性能，提出了两阶段连续带状缩减。1）第一阶段涉及正交相似变换，将对称稠密矩阵分解为带状形式：3020PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔，魏绍帅，鲁奇∙沙，大友浩之，横田理央，吴攀若0该阶段涉及正交相似变换，将对称稠密矩阵分解为带状形式：0� = � × � × � − 10其中 �是一个对称带状矩阵；2）在第二阶段，进行凸包追踪变换，将带状形式减少为三对角形式。一些相关研究表明，使用两阶段方法可以减少通信和计算瓶颈，从而获得更好的吞吐量[24]。然而，目前还没有关于基于TensorCore的特征值分解或连续带状缩减的研究。存在两个挑战。首先，TensorCore中的低精度算术需要仔细设计算法。例如，混合精度矩阵分解算法通常采用近似迭代方案，先从低精度分解得到近似解（预条件器），然后迭代地将解精确到更高精度。这种近似迭代方法在特征值分解中更加困难。根据以前的文献，只有SICE算法[38]部分满足要求，因为该算法只能处理部分所需的特征值和特征向量。其次，传统算法的矩阵形状无法充分加速TensorCore，无论是传统的三对角化还是两阶段三对角化都无法避免大量的高而瘦的GEMM运算，在TensorCore上运行速度较慢。特别是，在SBR中使用ZY表示法可以利用尾部矩阵更新中的对称性质（使用syr2k例程），但是TensorCore不支持这种类型的GEMM运算。因此，为了提高性能，我们专注于SBR，并证明在本文中使用修改后的WY表示法而不是传统的ZY表示法可以显著提高性能，通过“挤压”GEMM运算。此外，对于精度问题，如果需要单精度，我们还引入了一种名为误差相关Tensor CoreGEMM（EC-TCGEMM）的新技术，从一开始就将精度提高到单精度。因此，我们认为我们的贡献是：0•我们设计、实现和评估了一个快速稳定的高而瘦QR面板，与MAGMA和cuSOLVER面板分解相比，速度提高了约0•我们分析和研究了在SBR中使用WY表示法[5]和传统的ZY表示法[12]的性能，这解决了通常不会被考虑的算法（基于WY的）在新架构（TensorCore）上可以比传统算法（基于ZY的）表现更好的问题。0•我们提供了几种适用于不同情况的SBR实现。与最先进的MAGMA相对例程相比，加速比为3.7倍。0（半精度）和1.8倍（单精度），如果不需要特征向量。0•我们检查了将MAGMA的bulge-chasing过程与分而治之和我们的SBR实现相结合的可行性，整体EVD性能提升高达2.3倍。本文的其余部分组织如下：第二节介绍了相关工作，包括带状化、三对角化和TensorCore技术；第三节介绍了关于三对角化、Householder变换和对称带状化的一些背景知识；第四节分析了基于WY的算法性能更高的原因，并展示了关于内部GEMM的一些实验结果；第五节介绍了我们的实现和优化的细节，并展示了一些相关的实验结果；第六节将性能和精度结果与最先进的软件和工具进行了比较，第七节得出了结论并描述了我们的未来工作。02 相关工作对对称特征值分解问题的研究有着悠久的历史。此外，该领域的算法和方法已经得到了广泛的研究。02.1 三对角化和三对角求解器0最广泛使用和知名的算法是QR算法[39]，它重复调用QR分解和GEMM，并最终收敛到一个包含其对角线上特征值的对角矩阵。然而，在QR和GEMM迭代之前，通常会进行一个三对角化步骤作为“预处理器”，以减少QR迭代的次数。通常，三对角化是通过Householder变换[14]来执行的，为了提高在现代高性能架构上的执行效率，通常会应用WY表示技术[5,34]来进行变换过程。另一种三对角化方法称为两阶段三对角化[19]，它首先将矩阵降低到带状形式（第一阶段），然后将带状形式降低到三对角矩阵（第二阶段）。这种方法在多核架构上被证明非常高效[24,25]。另一种流行的方法是分而治之，这种方法在大多数线性代数包中都有实现，包括LAPACK[2]，MAGMA[37]和CuSOLVER[1]。另一种灵活的方法称为二分法[10]，它旨在找到特征值的一个子集，例如最大/最小的100个或在区间[�,�]内的所有特征值。此外，在2004年，提出了MRRR[11]方法。它寻求计算精确的正交特征向量，而无需昂贵的（�(�^3)最坏情况）重新正交化。01 https://docs.nvidia.com/cuda/cusolver/index.html3030通过WY表示在Tensor Core上进行快速对称特征值分解 PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔02.2 其他特征值分解方法极分解[22]将矩阵分解为一个正交矩阵和一个正定对称矩阵的乘积。它与特征值分解和奇异值分解有关，因此最近提出了一些基于极分解的新算法。QDWH-eig（基于QR的动态加权Halley特征值分解）[30]使用QR分解来计算极分解，然后使用迭代子空间方法对导出的正交矩阵进行因式分解。2016年，提出了QDWH-eig和QDWH-SVD的GPU实现[36]，但该工作将QDWH-eig替换为两阶段特征值分解。另一种计算极分解的方法称为缩放牛顿法[8]，它的数学运算比QDWH少。然而，它高度依赖于矩阵的反稳定逆。最近，随机线性代数[28]引起了越来越多的关注，特别是用于计算低秩近似特征值/奇异值分解的随机子空间迭代。其中两种算法是随机子空间迭代[16]和随机分块Lanczos[40]。这两种算法在现实应用中被证明是高效的，特别是在现代高性能架构上[35,43]。然而，这些算法通常基于乘以一个随机生成的矩阵，这意味着它们只能应用于对精度不敏感的应用。02.3基于张量核心的线性代数算法Nvidia在其Volta架构[27]中于2017年引入了张量核心加速器技术。后来在2020年，发布了更强大的Tensor Core单元的A100GPU[9]，其半精度GEMMs执行速率可以提升到300TFLOPS。张量核心单元旨在加速深度神经网络中的GEMMs。一些其他应用程序寻求利用张量核心加速线性代数算法，包括BLAS3操作[42]、LU分解[20, 23]和QR分解[41, 44]。03 背景0本节将介绍一些背景知识和相关技术和算法来进行带状化降阶。03.1三对角化三对角化过程通常是特征值分解的预处理步骤。三对角化的目标可以表达如下：0�=�−1×�×�0�是正交矩阵，�是三对角矩阵。使用Householder反射的传统方法消除除三对角部分以外的元素。然而，只有50%的计算可以完全块化（BLAS3）。换句话说，50%的计算是BLAS2操作，这些操作0图1. 2阶段三对角化0在现代硬件上的受益较少。事实上，根据MAGMAssyevd例程的实验结果，未块化的计算占据了三对角化（ssytrd例程）执行时间的90%以上。增加块化百分比的另一种方法称为2阶段三对角化[19, 24,25]。该方法在三对角化之前添加了连续的带状降阶。从数学上讲，连续的带状降阶可以用以下方程表示：�=�−1×�×�其中�是带状矩阵，带宽为�。此外，这是第一阶段。第二阶段将带状形式减少为三对角形式，称为bulgechasing[6]，第二阶段的性能取决于带宽�。图1显示了2阶段三对角化的一般步骤。03.2 Householder变换Householder反射器是通过反射来构造的正交投影，用于将给定的向量正交地转换为一个轴（从而消除除一个之外的所有分量）。具体来说，给定一个向量�，正交矩阵�(�)=�−2���/(���)将将�映射到第一个轴：�(�)�=[||�||, 0, ...,0]�。上述变换是一个秩-1更新，换句话说，是一个BLAS2操作。幸运的是，我们可以将多个Householder变换累积到WY表示[5]中的一个块中，以丰富BLAS3操作。假设我们有�个Householder矩阵[�1, �2, ..., ��]，则WY表示将为：0����−1...�2�1=�−�����如果(�+1)块被因式分解，那么我们有：0��+1=�−��+1���+1��+1=[��|��+1]��+1=[��|��+1−�������+1]03.3 从完整矩阵到带状矩阵的算法0对称矩阵Householder更新，也称为2边更新，非常简单。参见图2。第一个块包括�1和名为Panel的灰色区域。因为我们希望将矩阵减少为带状形式，所以只需要对Panel进行因式分解。第一步完成后，我们有��(�����)=(�−���)�，矩阵�将覆盖Panel的上三角部分，同时下三角部分𝐴2 = 𝐼𝑊𝑌𝑇 −1𝐴2 𝐼𝑊𝑌𝑇(1)𝑍 = 𝐴𝑊 − 12𝑌𝑊 𝑇𝐴𝑊(2)𝐴2 = 𝐴2𝑌𝑍𝑇𝑍𝑌𝑇(3)𝑘TC-GEMMSGEMMTC-GEMMSGEMM326.289.3620.029.316411.699.6533.309.8512824.4410.2249.8310.0225642.6510.3397.4110.2351266.5710.36122.8910.33102485.7310.40138.8210.372048112.0812.91121.5513.134096133.1715.31140.8514.333040PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔，魏沙帅，鲁奇∙沙，大友浩之，横田理央，吴攀若0图2.对具有块大小�的完整矩阵进行2边Householder更新的单步骤0面板的其余部分将被设置为0。进一步的尾部矩阵更新是双边的：0实际上，我们通常使用数学上等价的ZY表示[ 12 ]对尾部矩阵� 2 进行秩为2k的更新：0其余的更新过程类似，我们可以将更新后的 � 2视为一个新的完整矩阵，并进行迭代地分解矩阵。04 ZY-based和WY-based算法之间的性能分析0本节将主要讨论ZY-based和WY-based算法之间的差异，并解释为什么我们的新WY-based算法在TensorCore的环境中表现更好。我们还将提供一些实验结果来证明WY-based算法在Tensor Core上具有更好的GEMM性能。04.1高瘦GEMM在TensorCore上的一般性能表1显示了SBR中的两种类型的GEMM。第一种是一个方形矩阵乘以一个高瘦矩阵，出现在方程1和2中（� × � ）；第二种是一个外积，出现在方程1和3中（ � × � � ， × � � 和 � × � � ）。注意，Tensor Core尚不支持 ��� 2 �例程，因此方程3必须涉及两个外积。根据表1，GEMM的性能与Tensor Core上的 � 密切相关，而SGEMM的性能在 �增加时更加稳定。这可能是因为高瘦GEMM在 �较小时变为内存限制而不是计算限制，并且在TCGEMM中启动内核的时间成本并不微不足道。显然，要改善SBR的整体性能0表1. 在A100GPU上的TCGEMM和SGEMM性能，以TFLOPS为单位，当 �从32变化到4096，且固定 � = 32768。在第2-3列中， � ∈ R �× � ， � ∈ R � × � 。在第4-5列中， � ∈ R � × � ， � ∈ R � × � 。0在Tensor Core上，我们希望 �尽可能大。然而，不幸的是，在传统的WY和ZY表示算法中，�被固定为通常小于256的带宽。由于突出追踪的计算复杂度为� ( �� 2 ) ，将块大小设置得太大会有一定的代价。04.2一种改变GEMM形状的替代方法受基于TensorCore的QR分解[ 41]的启发，从迭代方法到递归方法的简单修改可以改变GEMM的形状并显著加速。然而，递归策略似乎对SBR不起作用。与QR分解不同，SBR是一种双边分解；尾部矩阵不能简单地分为左半部分和右半部分。04.2.1尝试使用ZY表示改变GEMM形状0假设矩阵 � 和 � 是SBR的第 � 次迭代中的ZY表示。那么整个 �和第 � 次迭代后的 � 矩阵将分别为 � = [ � 1 | � 2 | ... | � � ] 和 � =[ � 1 | � 2 | ... | � � ]。为了“压缩”GEMM，我们可以尝试将几个Z和Y向量合并在一起，并一次更新尾部矩阵。一种解决方案是找到一个类似于基于Tensor Core的QR分解[ 41]的递归策略。然而，与QR分解不同，SBR中的尾部矩阵更新是双向的。当我们尝试更新下一个面板时，我们需要先前的尾部矩阵完全更新。将高瘦GEMM转换为递归QR分解算法[ 41]中的方形GEMM的关键思想是在需要分解尾部矩阵时更新尾部矩阵，而不是在每次迭代中更新整个尾部矩阵。这种方法在递归QR分解中效果很好，因为它是单边的3050通过WY表示的快速对称特征值分解在Tensor Core上的PPoPP '23，2023年2月25日至3月1日，加拿大蒙特利尔0图3.使用ZY表示的SBR的一个假设，绿色区域是要更新的面板，蓝色轮廓内的矩阵是尾部矩阵0因式分解。在两侧因式分解中，情况变得不同。图3展示了我们在SBR中使用ZY表示的假设。假设我们有8个块，块大小为�，并且已经使用ZY表示因式分解了第一个块。现在我们只想通过 �� = �� − �� ( 1 : �, : ) � − �� ( 1 : �, : ) �来更新矩阵中的绿色区域（我们可以称之为 �� ），而棕色区域（ ���）则从未被触及。下一步（图3中的第2张图片）是因式分解矩阵中的红色区域并更新绿色区域（ �� ）。下一步是更新第2张图片中的 �� ，如果我们有 � = [ � 1 | � 2 ] ; � = [ � 1 | � 2] ，那么可以通过 �� = − �� ( 1 : 2 � �, : ) � − �� ( 1 : 2 � �, : ) � 来轻松更新 ��。然而，问题是如何获得 � 2 。正如我们之前提到的，矩阵 �与尾部矩阵有关，实际上，如果我们已经有了 � 2 和 � 2，并且假设带有蓝色轮廓的区域为 �� ，我们可以通过 � 2 = �� × � 2 − 1 来计算 � 2 。02 � 2 � � 2 × �� × � 2 。然而，要推导出正确的 �� ，子矩阵 需要通过 �� = �� − � 1 � 1 ( � + 1 : �, : ) � − � 1 � 1 ( � 进行更新，这也包括了高瘦矩阵的乘法。我们没有找到绕过使用ZY表示更新整个尾部矩阵的方法。04.2.2使用WY表示改变GEMM形状的替代策略。WY表示与ZY表示在尾部矩阵更新方面的区别在于，WY表示通过乘法来更新矩阵，而ZY表示则使用连续的减法。此外，WY表示中更新 �矩阵不需要更新整个尾部矩阵，而只需要更新面板。WY表示的这些特性使我们能够累积Householder矩阵并一次性更新尾部矩阵，但代价是需要更多的计算和内存使用，因为它必须在每次迭代中更新 � 矩阵。0图4.基于WY表示的SBR，红色区域是要因式分解的面板，绿色区域是尾部矩阵。蓝色轮廓内的矩阵表示原始矩阵0图4描述了在第一个大迭代之前WY-basedSBR的几个步骤。显然，第一步是在(a)中对红色面板进行QR分解并得到 � 1 和 � 1 。为了只更新(b)中的绿色区域（ ��），我们将让 � = � 1 ( 1 : �, : ) ，并在右乘法中使用 �而不是整个 � 1 。然后 �� = ( � − � 1 � � 1 ) × �� × ( � − � 1 � � )。下一步将是对(c)中的红色区域进行面板因式分解，并形成新的 � = [ � 1 | � 2 ] 和 � = [ � 1 | � 2 − � 1 � � 1 � 2 ]。在递归之前，我们需要更新整个尾部矩阵。因为我们以前从未触及过(d)中的绿色区域，所以我们仍然需要使用原始矩阵的带有蓝色轮廓的子矩阵（ �� ）。然后我们将 � 替换为 � ( 2 � � + 1 : �, : ) ，最后更新 �� = ( � − �� � ) × �� × ( � − �� � )。算法1给出了基于WY表示的SBR的Matlab原型。参数 ��应该作为一个更大的块大小给出，其中只有面板将被更新；在大块外部，整个尾部矩阵将被更新。04.3 GEMM性能评估04.3.1算术运算：传统的ZY表示和算法1之间的一个区别是算术运算的数量。算法1中的运算增量包括在内循环中构建�矩阵以及算法1中的第9行和第13行的较大的GEMM，因为我们总是需要使用内循环中的原始矩阵的大子矩阵（第3行）。dominates the overall performance improvement. Neverthe-less, when 𝑛𝑏 is larger than 1024, the rise of computationsprevents the performance from improving.Fix 𝑛𝑏 to be 1024, and Figure 6 compares the total execu-tion time of the GEMMs between the two different strategies.When the matrix is as small as 4096, the WY-based algorithmdoes not have an advantage over the ZY-based algorithm.This is obvious because the GEMM execution rate in the twoalgorithms is similar when the matrix size is small. Thereforethe change of GEMM shapes does not help. In contrast, theincrement of the mathematical operations cannot be over-looked. Furthermore, that is why when the size is 4096 and8192, the GEMM execution time in the WY-based algorithmis longer than in the ZY-based algorithm. However, when thesize is quite large, the better shapes of the GEMMs improvethe performance significantly. For instance, when the matrixsize is 32768, the inside GEMMs in the WY-based algorithmcan reach 240 TFLOPs, while the highest GEMM rate in theZY-based algorithm is only around 50 TLFOPs. This meansthat even though we are doing more computations in theWY-based algorithm, the overall performance of GEMMs canstill benefit a lot from the relatively square GEMMs (around1.5x speedup in terms of GEMMs). We also perform anotherexperiment on testing the SGEMMs and show the results inFigure 7. According to Table 1, the SGEMM cannot benefitfrom the square shapes. As a result, we can find in Figure 7that if TCGEMMs are replaced with SGEMMs, the ZY-basedalgorithm can have better performance, which means theWY-based algorithm only brings speedup with Tensor Coresupport.3060PPoPP ’23, 2023年2月25日至3月1日，加拿大蒙特利尔，张少帅，Ruchi Shah，Hiroyuki Ootomo，Rio Yokota和Panruo Wu0算法1 递归的基于WY的对称带约减法，块大小为nb01 function [A] = sy2sb (A, oriA)02 n = length(A);03 OA=oriA(b+1:n, b+1:n);04 for i=1:b:nb05 [w,y]= PanelQR(A(i+b+1:n, i:i+b)) ;06 Y=[Y y];07 W=[W w-W*Y'*w];08 GA=A(i+b+1:n, i+b+1:i+2*b);09 GA=(I-W*Y')'*OA*(I-W*Y(i:i+nb,:)');010 A(i+b+1:n, i+b+1:i+2*b)=GA;011 结束012 GA=A(nb+1:n, nb+1:n);013 GA=(I-W*Y(nb+1:n,:)')'*OA*(I-W*Y(nb+1:n,:)');014 A(nb+1:n,nb+1:n)= sy2sb(GA,GA) ;015 结束0ZY WY0块大小 128 128 256 512 1024 2048 40960FLOPS 0.70 0.93 1.05 1.12 1.17 1.22 1.310表2.基于ZY的SBR（带宽为128）和基于WY的SBR（块大小从128到4096不同）的实际算术运算数量；数字的指数为10的14次方0在内循环中，ZY算法和WY算法之间的一个区别是算术运算的数量增加。为了直观地显示算术运算的增量，我们分别计算了基于ZY和基于WY的算法在矩阵大小为32768×32768时的操作数量，结果如表2所示。随着块大小的增加，操作数量显著增加，这说明为什么ZY算法在以前的硬件上更受欢迎。04.3.2 GEMM性能的实验结果0WY算法确实将一些高而瘦的矩阵乘法转换为相对较方的矩阵乘法。然而，问题在于它增加了总的数学运算量，包括形成分块的�矩阵和更大的矩阵大小，因为原始矩阵的大子矩阵（算法1中的第3行）在大块内始终被重复使用。因此，WY算法能否带来加速仍然未知。让我们从块大小��变化时的性能评估开始。较小的��导致较少的计算量但更多的高而瘦的GEMM，而较大的��导致更多的计算量和方形的GEMM。我们的直觉是存在一个��可以给我们最佳性能，图5支持了我们的猜测。当��小于1024时，尽管浮点运算增加，但方形GEMM的增长主导了整体性能改进。然而，当��大于1024时，计算量的增加阻碍了性能的提升。将��固定为1024，图6比较了两种不同策略下GEMM的总执行时间。当矩阵大小只有4096时，WY算法在性能上并没有优势。这是显而易见的，因为当矩阵大小较小时，两种算法中的GEMM执行速率是相似的。因此，GEMM形状的变化并没有帮助。相反，数学运算量的增加不能被忽视。此外，这就是为什么当大小为4096和8192时，WY算法中的GEMM执行时间比ZY算法中的GEMM执行时间长的原因。然而，当大小非常大时，更好的GEMM形状显著提高了性能。例如，当矩阵大小为32768时，WY算法中的内部GEMM可以达到240TFLOPs，而ZY算法中的最高GEMM速率仅约为50TLFOPs。这意味着即使我们在WY算法中进行了更多的计算，GEMM的整体性能仍然可以从相对方形的GEMM中受益（GEMM方面的加速约为1.5倍）。我们还对SGEMM进行了另一个实验，并在图7中显示了结果。根据表1，SGEMM无法从方形形状中受益。因此，我们可以从图7中发现，如果将TCGEMM替换为SGEMM，ZY算法的性能会更好，这意味着WY算法只有在有Tensor Core支持时才能带来加速。0图5.当块大小从128变化到4096时，TCGEMM的总耗时。点上的数字是TCGEMM的TFLOPs。04.4 形成特征向量0当需要特征向量时，基于 WY 的算法中的 GEMM更加高效。因为在内积中形成 �矩阵（反变换）根本不会浪费。由于在每个内循环中完全获得� 矩阵，所以最终矩阵 � 将具有以下形式：� = [�1 | �2 | ... |��]，其中 � = � / ��。然后形成整个 �将非常容易，并且可以以递归方式形成。3070通过张量核心的 WY 表示进行快速对称特征值分解 PPoPP '23，2023 年 2 月 25 日至 3 月 1 日，加拿大蒙特利尔0图 6. WY-based 算法和 ZY-based 算法中 TCGEMM的总耗时比较，WY-based 算法中的块大小固定为 10240图 7. WY-based 算法和 ZY-based 算法中 SGEMM的总耗时比较，WY-based 算法中的块大小也固定为 10240算法 2 递归 � 构造，块大小为 nb01 function [W] = FormW (W,Y)02 n=length(W);03 if n<=2*nb04 W1=W(:,1:nb);05 W2=W(:,nb+1:2*nb);06 Y1=Y(:,1:nb);07 W=[W1|W2-W1*Y1'*W2];08 end09 %左递归010 W1= FormW (W(:,1:n/2), Y(:,1:n/2));011 %右递归012 W2= FormW (W(:,n/2+1:n), Y(:,n/2+1:n));013 W=[W1|W2-W1*Y1'*W2];014 end0基于我们的实验，给定一个大小为 32768 × 32768 的矩阵，0图 8. MAGMA、cuSOLVER 和我们的 TSQR 实现之间面板QR 分解耗时的比较，0在 WY-based 和 ZY-based算法中，反变换所需的时间分别为 320ms 和420ms，这提供了约 10% 的加速。05 实现和优化 在本节中，我们将介绍我们的 GPU实现和优化策略，以提高性能和准确性。05.1 高瘦面板 QR 分解 图 2显示了当带宽较小时，面板是高瘦的。由于低数据局部性和并行性，对这种矩阵形状进行因式分解比对一个方阵进行因式分解要慢得多，特别是在 GPU 上。避免通信的QR（CAQR）[3]因式分解算法可以高效地执行 QR 分解。高瘦 QR（TSQR）是 CAQR的一个特例，只处理高瘦矩阵。关于一般思想和实现可以在基于张量核心的 QR论文[41]中找到。但仍然有一些修改：1）我们在每个小块上应用 Householder变换而不是修改的 Gram-Schmidt过程以保持稳定性；2）每个线程束负责一列而不是一行，以获得更好的性能。我们通过将矩阵大小从 4096 变化到 32768 来评估带状约简中的面板因式分解，通过比较 MAGMA中面板因式分解的执行时间（ssytrd_sy2sb 20例程），cuSOLVER（sgeqrf 和 sorgqr 例程）和我们的TSQR；参见图 8 。请注意，图 8 中的耗时不是纯粹的 QR分解；它还包括在 cuSOLVER 和 TSQR 实现中重构 � 和 �矩阵，这将在下一节中讨论。02 对称矩阵转换为对称带状形式3080PPoPP '23，2023 年 2 月 25 日至 3 月 1 日，加拿大蒙特利尔，Shaoshuai Zhang、Ruchi Shah、Hiroyuki Ootomo、Rio Yokota 和 Panruo Wu05.2 重构 Householder 向量 尽管 TSQR 可以加速面板与cuSOLVER SGEQRF() 相比，但仍存在一个问题，即与cuSOLVER SGEQRF() 不同，我们从 TSQR 中获得的是显式的� 。相比之下，cuSOLVER SGEQRF() 提供了 Householder向量。此外，使用显式的 � 而不是 Householder向量，进一步的尾部矩阵更新将导致不稳定的结果。因此，有必要开发一种算法，可以从 TSQR 形成的显式的 � 中输出Householder 向量。一种解决方案是通过显式的 � 重构Householder 向量[4]。其思想如下。给定一个正交矩阵，�可以表示为 � = � - � × � × � �（内存高效的 WY表示），并且这个方程也可以修改为 � - � = � × � × � � 。�是一个下三角矩阵，�是一个上三角矩阵。因此，它可以被视为具有 � - � = � × � = (�) (� × � �) 的 LU 分解。文献[4]还报告了 LU 分解提供了唯一的 � �，并且部分主元选取是不必要的。在我们的算法中，我们对矩阵 � - �的上部分进行因式分解，然后进行三角求解（STRSM）以获得整个 � 矩阵。在通过 LU 分解获得 Householder 向量 �后，我们进行另一个三角求解以构造 � ；参见算法 3 。0算法3 通过显式由TSQR生成的�重构WY表示01 function [W,Y] = 重构WY(Q)02 [m,n] = size(Q);03 I = 单位矩阵(m,n);04 A = I-Q;05 [L1,U] = 非主元LU(A(1:n,:));06 L2 = A(n+1:m,:)/U;07 Y=[L1;L2];08 W=A/Y';09 结束0Ballard等人在�−�[4]上执行LU分解，其中�是与HouseholderQR算法内部的符号选择相对应的对角符号矩阵。请注意，这一步也是避免LU分解中秩缺失的关键。将此算法与我们的TSQR实现相结合，我们可以得到生成�和�的面板分解。05.3 基于张量核心的GEMM的误差校正优化的TCGEMM[31]是一种基于张量核心的GEMM，它提供与CUBLASSGEMM使用FP32相同的精度，但超过了FP32的理论峰值性能，从而实现更高的吞吐量。它们是Markidis[26]提出的误差校正方法的改进版本，用于在张量核心上执行矩阵乘法。实现误差校正方法可以0不能提供FP32SIMT核心上矩阵乘法的准确结果。实现高精度的关键是处理张量核心内的舍入误差。优化的TCGEMM保存尾数损失并使用它来纠正矩阵乘法的精度。这种高精度、高性能和低功耗的实现是在NVIDIA的CUTLASS库中开发的，使用NVIDIAA100GPU上的TF32张量核心，对于有限指数范围可达到51TFLOPS/s，对于完整指数范围可达到33TFLOPS/s。该实现超过了19.5TFLOPS/s的理论FP32SIMT核心峰值性能。具体来说，考虑一个GEMM�=�×�。要在张量核心上计算这个GEMM，需要提前将�和�从FP32截断为FP16。然后我们有�=˜�+Δ�和�=˜�+Δ�，其中˜�和˜�以FP16精度存储。带有FP32累加的TCGEMM实际上执行的是�=˜�ט�，准确的结果应该是�=(˜�+Δ�)×(˜�+Δ�)=˜�ט�+˜�×Δ�+Δ�ט�+Δ�×Δ�。在实践中，Δ�×Δ�非常微小，通常可以忽略。因此，执行两个额外的TCGEMM似乎可以恢复精度损失。但是张量核心内的舍入误差和下溢问题仍然阻止了朴素的恢复解决方案获得准确的结果。因此，作者还应用了另外两种方法来解决这个问题：1）尝试避免张量核心内累加的舍入以恢复FP32的精度；2）尝试对矩阵进行缩放以减少下溢。最后，优化的TCGEMM可以超越SGEMM，并同时将精度恢复到单精度。我们还将在图10中展示使用误差相关TCGEMM对我们的SBR的性能影响。06 实验评估0我们在所有实验中使用了5.4.0-9
我的内容管理
展开
