一种利用GPU优化大规模小方阵奇异值分解的新方法

李晓敏，鄢社锋，侯朝焕

(中国科学院 声学研究所，北京 100190)

传统意义的 GPU主要针对游戏加速和图形图像处理。NVIDIA公司于2007年发布了基于CUDA的GPU，这种意义上的GPU具有很强的并行计算能力和数据吞吐能力，在诸多领域里能使程序性能提高好几个数量级，对未来信息社会处理海量数据的需求具有很强的适应性。

基于CUDA的GPU自问世4年以来，由于其在性能、成本和开发周期上的显著优势，一经推出即在学术界和产业界引起了热烈反响，已经广泛应用于金融、石油、天文学、流体力学、信号处理、电磁仿真、模式识别、图像处理、视频压缩等众多领域[1]。

求解一般矩阵的 SVD最常用算法可以分为两大类：基于QR分解的算法和基于Jacobi旋转的分解[2]。前者首先通过一系列Householder变换实现双对角，然后通过QR变换实现对角。其中，QR分解是一种严格串行的过程，只能获得有限的并行度，因此不适合应用于GPU。后者使用平面旋转进行正交化。双边Jacobi每次正交化矩阵的两行和两列，调度时同时涉及到行和列，不适合 GPU的并行架构。而单边Jacobi每次旋转只需要正交化矩阵的两列，很容易实现并行访问。

目前关于在GPU上进行SVD的研究并不多。最初，文献[3]采用在GPU上实现双对角变换和在CPU上实现对角变换的混合编程法完成了SVD。2009年，文献[5]首次使用 GPU实现了双对角-对角算法。同年，文献[6]在 3D重建过程中将VL数据库里的SVD程序转换成CUDA语言，从而实现实时重建。2010年，文献[7]采用 Webb算法在GPU上实现了SVD。然而，这些方法相对于MATLAB和IntelMKL的优势仅仅体现在针对单个大矩阵SVD，并且当矩阵维数增大到成千上万时这种优势更加明显。到目前为止，仅有文献[8]实现了对单个512*512小型实矩阵 SVD的优化，而针对多个小型复矩阵SVD的优化几乎没有。

本文充分利用了 GPU内众多并行内核以及共享存储器和常数存储器隐藏访问时间的优势，借助Brent和Luk提出的Jacobi行列调度法[9]，采用单边 Jacobi旋转，迭代实现了大规模小方阵的 SVD优化。

1 单边Jacobi方法

对M*M维方阵A做SVD分解：

单边Jacobi方法通过一系列平面旋转(一般文献中在介绍Jacobi方法时为了方便也将平面旋转称为Jacobi旋转，即右乘正交阵V)使矩阵A正交化。首先，对矩阵A重复进行Jacobi旋转，直到旋转后的矩阵W各列两两正交，得到右奇异向量V；然后，对正交矩阵W归一化，得到奇异值向量和左奇异向量U。

1.1 旋转正交化过程

在矩阵 A 中任取一组列对(i，j)，计算第 i列内积存入a，计算第j列内积存入b，计算第i列和第j列的差积存入c；根据a，b，c计算旋转角cs和旋转因子sn；根据cs和sn对A中(i，j)两列的每个元素进行代数运算，从而实现这两列的正交，与此同时，对单位矩阵V中(i，j)两列也进行同样的代数运算。循环重复以上步骤，直到|c|/小于一个接近于0的正小数 tol，即表明A中任意两列均正交。

1.2 归一化和计算SVD

2 基于GPU的SVD实现

2.1 优化方法

CUDA是一种将GPU作为数据并行计算设备的软硬件体系，采用了比较容易掌握的类C语言进行开发。它是一个 SIMD(Single Instruction Multiple Data)系统，即一个程序编译一次以后，CUDA将计算任务映射为大量的可以并行执行的线程，并由拥有大量内核的硬件动态调度和执行这些线程，从而显著提高运算速度。如图1所示，将一个可以并行化执行的任务首先分配给若干个线程网格(Grid)，其次将每个Grid内的任务分配给若干个线程块(Block)，最后再将每个Block内的任务细分给若干个线程(Thread)。Grid中的所有Blocks并行执行，Block中的所有 Threads并行执行，这种两层并行模型是 CUDA最重要的创新之一。

图1 GPU内线程结构Fig.1 Configuration of threads in GPU

CUDA在GPU内定义了8种寄存器，图2所示展示了常用的4种。其中，共享存储器(Shared Memory)是GPU片内的高速存储器，通常用于存储中间结果或 Block的公共结果；常数存储器(ConstantMemory)是只读的地址空间，虽然位于显存，但是拥有缓存加速，在CUDA程序中用于存储需要频繁访问的只读参数。

本文的SVD优化方法将任务细分给多个Blocks内的多个 Threads，不同维数的矩阵，任务细分策略不同；同时，将需要频繁使用的矩阵数据和中间结果存入共享存储器，有效隐藏存储器访问延迟；另外，采用Brent和Luk提出的Jacobi行列调度法，将正交化过程中需要选取的所有列队序号制成表格，存入常数存储器，避免多次重复取数，使程序高效化执行。

图2 GPU内主要存储器结构Fig.2 Configuration of main memory in GPU

在每次任两列旋转正交之前，均执行算法1中的步骤。首先，初始化Lk=2k 1，Rk=2k。

算法1 Jacobi行列调度算法

1：如果 Lk

2：如果k=1，则outRk←Rk

3：否则如果k

4：如果k>1，则outLk← Rk，直到所有Rk赋值给 outLk；

5：如果 k>M/2，则 Rk←inRk，否则 Rk←Lk；6：如果k>1，则Lk←inLk

对于M维方阵，需要(M-1)轮正交化才能保证任意两列之间均正交化。以M=6为例：

第一步正交的列对：(1，2)，(3，4)，(5，6)；

第二步正交的列对：(1，4)，(2，6)，(3，5)；

第三步正交的列对：(1，6)，(4，5)，(2，3)；

第四步正交的列对：(1，5)，(6，3)，(4，2)；

第五步正交的列对：(1，3)，(5，2)，(6，4)；

在旋转变换之前，计算每一步需要正交的列对编号，并将其制成表格。

2.2 优化过程

本文显卡采用NVIDIAGTX460，拥有48KB共享存储器和16KB常数存储器。设矩阵维数为M，需要做SVD的矩阵个数为N。

由于单边 Jacobi方法每次需要正交化一个列对，所以矩阵维数最好为2的整数倍，本文分别选取16，32和64维的情况进行说明。执行正交化之前，文献[10]将所有复矩阵转化为具有等价奇异值和奇异向量的实矩阵，矩阵维数由原来的M增加到2*M。为了让描述简洁直观，之后的介绍中 M实际代表2*M。

由于本方法旨在充分利用共享存储器和常数存储器，因此要特别注意所用资源不应超过 GPU提供的上限，不同维数的矩阵，任务细分策略不同。当矩阵维数为16时，经过计算可以将A和V两个矩阵完全从全局存储器拷贝到共享存储器。使用N个Blocks，每个Block处理一个小方阵A，Block内有(M/2)*M个Threads，分别处理每次正交化过程中M/2个列对的M个方阵元素；当矩阵维数为32和64时候，经过计算 shared memory已经不能完全存储A和V。因此，使用N*(M/2)个Blocks，每个Block只读入一个列对到共享存储器中进行处理，内有2*M个Threads，正交化该列队的2*M个矩阵元素。

处理步骤如下：

1.将待处理的数据读入shared memory。

当矩阵维数为16时。每个Block将要处理的矩阵A从global memory中读入shared memory存为dA[M*M]，同时将单位阵V也读入shared memory存为dV[M*M]。每个Block里有(M/2)个列对，因此定义 cs[M/2]、ss[M/2]和 c[M/2]分别存放每个列对里第1列内积、第2列内积以及两列差积。

当矩阵维数为32和64时。根据当前Block的ID号查表得到所处理的是哪两列，每个Block将要处理的两列数据从矩阵A中读入shared memory存为DA[2*M]，同时将单位阵V的相应位置两列也读入 shared memory存为 DV[2*M]。每个Block里有1个列对，因此定义cs、ss和c分别存放第1列内积、第2列内积和两列差积。

2.求内积、差积和旋转因子。

当矩阵维数为16时。根据当前Thread的ID号查表得到所处理的是哪两列，求得内积和差积放入 shared memory中定义的临时区域 temp[M*(M/2)]，并对temp进行规约求和将结果放入cs、ss和c中。

当矩阵维数为32和64时。得内积和差积放入shared memory中定义的临时区域temp[2*M]，并对temp进行规约求和将结果放入cs、ss和c。

3.根据cs、ss和c求旋转角后对矩阵进行旋转，并将旋转后的 dA和 dV(M=16)及 DA和 DV(M=32，64)放回原矩阵A和V中。

4.重复2和3，直到 c小于一定的值，说明当前一系列列对的正交化已经基本完成。

5.在主机端重复调用1-4共(N-1)次，遍历完Jacobi行列调度表，保证经过旋转后的A任意两列之间均正交。

所有Blocks里的全部Threads同时执行以上1-4步，实现了最大并行化。所有4步均在sharedmemory中进行，其中 Jacobi行列调度表放在 constant memory中，大大减少了存储器访问时间，提高了数据吞吐量。

3 仿真结果

本节采用基于Intel Dual Core 2.10GHz CPU和NVIDIA GTX460显卡，分别测试了基于MATLAB 2008a、Intel MKL 10.0.0 LAPACK 和 NVIDIA GTX460上SVD的性能。

随机生成一批单精度32维方阵，图4显示采用MATLAB、MKL和GPU执行SVD的时间随矩阵个数的变化曲线。

图4 执行时间随矩阵个数的变化(M=32)

从图4可以看出，针对多个小规模方阵，GPU的执行时间明显小于MATLAB和MKL。同时，基于MATLAB和MKL的SVD执行时间随矩阵个数的增加成直线增长，而基于GPU的SVD执行时间随矩阵个数的增加明显趋于缓慢。

图5 GPU相对于Matlab和MKL的加速比随矩阵维数和矩阵个数的变化Fig.5 Speed up of GPU compared with Matlab and MKL mutative with matrix size and matrix number

图 5展示了针对不同数量不同维数小方阵SVD，GPU相对于MKL和MATLAB的加速比曲线图。

对于128个16维小方阵的SVD，GPU执行速度相比于 MATLAB提高了17.6倍，相比于 MKL提高了9.2倍。当矩阵维数大于16时，随着矩阵个数增加，GPU执行SVD相比于这两种方法的加速比均趋于稳定，分别约为5.1倍和3.4倍。

4 结束语

介绍了一种基于GPU针对大规模小方阵SVD的优化方法，适用于声纳和雷达信号处理算法。该方法充分利用了 GPU众多并行内核以及共享存储器和常数存储器隐藏访问时间的优势。在同等条件下，相比于MATLAB和IntelMKL都有一定加速。

[1]张舒，褚艳丽.GPU高性能计算之CUDA[M].北京：中国水利水电出版社，2009：1-3.

[2]张贤达.矩阵分析与应用[M].北京：清华大学出版社，2008：358-365.

[3]Bondhugula V，Govindaraju N，Manocha D.Fast Singular Value Decomposition on Graphics Processors[R].University of North Carolina at Chapel Hill，2006.

[4]Singular Value Decomposition Routines[DB/OL].http://www.culatools.com/features/lapack/.

[5]Sheetal L，Narayanan P J.Singular Value Decomposition on GPU using CUDA[J].IEEE Symposium on Parallel&Distributed Processing，2009：1-10.

[6]Kimikazu K，Tikara H.Singular Value Decomposition for Collaborative Filtering on a GPU[J].IOP Conference Series:Materials Science and Engineering，2010，10(1)：1-5.

[7]Real-time 3D registration of stereo-vision based range images using GPU[C].IEEE Workshop on Applications of Computer Vision(WACV)，2009：1-6.

[8]张舒，窦衡.基于CUDA的矩阵奇异值分解[J].计算机应用研究，2007，24(6)：104-107.

[9]Brent R P，Luk F.T.The solution of Singular Value and Symmetric Eigen problems on Multiprocessor Arrays[J].Sct Stat Comput，1985，6：69-84.

[10]李小波，薛王伟，孙志勇.一种求解复Hermite矩阵特征值的方法[J].数据采集与处理，2005，20(4)：403-406.