Randomized SVD 算法介绍与实现

导语

对于大型矩阵的分解，我们往往会想到用SVD算法。然而当矩阵的维数与奇异值个数k上升到一定程度时，SVD分解法往往因为内存溢出而失败。因此，文本介绍一种Randomized SVD算法，相比于SVD，它更能适应大型矩阵分解的要求，且速度更快。实验发现，在Tesla平台相同的资源配置下，同样的矩阵分解任务，SVD算法运行失败，而Randomized SVD算法在复杂的迭代计算过程下也仅耗时1.6h。

1. Randomized SVD算法原理介绍

1.1 Randomized SVD基本原理

简单的说，Randomized SVD算法也是一种矩阵分解算法。之前的文章《矩阵奇异值分解法SVD介绍》中详细介绍了SVD分解算法，本文的Randomized SVD分解算法是在SVD算法基础上实现的，下面将详细介绍该算法的原理。

Randomized SVD算法主要是在文章[1]中提出来的，它的主要计算过程分为两步：

构建一个能够捕捉到原始矩阵”行为”的低维子矩阵 将原始矩阵限制在低维子空间，并计算其标准的矩阵分解，例如SVD 上述两步的详细步骤为：

• Stage 1 从原始输入矩阵A的列空间中获得一个近似基Q，并满足如下条件：

其中，矩阵Q的列向量是正交的，指定A的行数与列数分别为m和n 上述过程中最关键的就是求取满足要求的矩阵Q，文章中给出了具体的方法，即使用随机采样方式构建矩阵Q： 1.通过设置需要获得的奇异值个数k以及过采样参数p，构建一个由k + p个n维的随机列向量组成的矩阵 Ω，要求(k + p) <= min{m,n}，每一个列向量中的值均采样自标准正态分布，因此，这些采样的列向量线性独立 2.进行矩阵乘积运算Y=AΩ，由于向量集合Y也是线性独立的，因此，Y形成了矩阵A的列向量空间 3.通过求取向量集合Y的正交基，从而得到A的近似基Q

• Stage 2

1.构建低维矩阵B，满足：

2.计算低维矩阵B的SVD分解，使得

从1中的公式我们可以看到，B是一个k+p行n列的矩阵，相比初始矩阵A(mn)，B的行数非常小，更易于进行SVD分解

3.计算矩阵A实际的左奇异向量U

通过以上的步骤就实现了初始矩阵A的SVD分解，可以看到，与之前的SVD分解过程相比，Randomized SVD算法主要多了一个构建随机向量的过程。

1.2 Randomized SVD算法的一般过程

根据上述的算法原理，我们知道Randomized SVD的计算过程为：

• 算法一：Randomized SVD算法的一般过程 Input: m行n列的初始矩阵A，奇异值个数k，过采样参数p，要求满足k+p<=min(m，n) Output: A的SVD分解结果，分别为左奇异向量UA ，奇异值矩阵ΣA，以及右奇异向量VA

1.构建一个n∗(k+p)维的高斯随机矩阵Ω 2.进行矩阵乘积运算Y=AΩ 3.利用QR分解获得Y的正交基Q=qr(Y) 4.构建低维矩阵B=QT A 5.矩阵B的SVD分解，[UB,ΣB,VB]=SVD(B) 6.用Q更新左奇异向量，UB=QUB

7.得到最终结果UA=UB(:,1:k),ΣA=ΣB(1:k,1:k),VA=VB(:,1:k)

1.3 Randomized SVD算法的改进

当矩阵的奇异值衰减较快的时候，算法一的表现非常好。然而，当我们遇到较大的矩阵，或者矩阵的奇异值下降缓慢时，算法一的结果往往不准确。这主要是由于较小的奇异值对应的奇异值向量干扰了计算的结果[2]。为此，Randomized SVD算法提出了改进算法，主要是通过矩阵的多次(the power of matrix)乘积运算来减小这些奇异值向量相对于主要奇异值向量的权重。对照算法原理，主要修改了Stage 1阶段的过程。下面主要阐述Stage 1阶段的过程。

• 算法二：Randomized SVD算法的Power 迭代过程

Input: m行n列的初始矩阵A，奇异值个数k，过采样参数p，要求满足k+p<=min(m，n)，power指数q Output: A的近似基Q

1.构建一个n∗(k+p)维的高斯随机矩阵Ω 2.交替使用A与AT 构建q轮的迭代过程Y=(AAT )qAΩY 3.利用QR分解获得Y的正交基Q=qr(Y)

另外，为了避免上述的Power迭代过程中数值较小的奇异值所携带的信息在计算过程中丢失，文章[1]中进一步提供了改进策略，即每次进行A或者A^T的乘积运算时，都对采样矩阵的列向量进行正交化。具体的过程见算法三。

• 算法三：Randomized SVD算法的子空间迭代过程

Input: m行n列的初始矩阵A，奇异值个数k，过采样参数p，要求满足k+p<=min(m，n)，power指数q Output: A的近似基Q

1.构建一个n∗(k+p)维的高斯随机矩阵Ω 2.矩阵乘积运算Y=AΩ，并通过QR分解获得其正交向量Y0 =Q0 R0

3.进行q轮的迭代过程，for j = 1,2,…q

4.Q=Qq

2. Randomized SVD算法的分布式实现

以上就是Randomized SVD算法的原理，接下来我们主要探讨Randomized SVD算法的分布式实现。当矩阵的维数非常大时，我们通常都会想到将这个矩阵进行分布式存储，并且采用以spark为平台实现的SVD算法来对矩阵以分布式的方式进行分解，目前这个算法已经发布在tesla平台，然而这种分解方式不光会占用大量的时长，同时还有内存溢出的可能，特别当矩阵的列数n非常大，或者要求奇异值的个数k非常大时，这种分解任务往往失败。因此，下面我们重点分析Randomized SVD在spark上的实现原理。

在上文的原理介绍中我们知道，利用k+p个随机采样的向量可以将原始矩阵的维数缩减至k+p维。当原始矩阵的维数m非常大时，k+p将远小于m，这时矩阵被缩减成一个非常小的矩阵，甚至不需要像原始矩阵那样采用分布式的方式存储，而是可以直接存储在本地。这种缩减方式将极大的降低算法的空间复杂度，同时由于SVD分解过程在维数较低的矩阵上进行，因此也节约了整个算法的运行时间。

2.1 QR分解的分布式实现

通过以上分析我们总结Randomized SVD算法在spark上的实现过程。首先当原始矩阵很大时，我们采用分布式的方式存储。其次Randomized SVD算法的关键就在于对原始矩阵的降维，从算法一、二到三我们可以看出，这就需要原始分布式矩阵右乘一个本地矩阵Ω，这在spark上是比较容易实现的。乘积的结果是一个分布式矩阵，所以接下来要对分布式矩阵进行QR分解，注意这里要分解的矩阵是一个m行(k+p)列的，由于k+p远小于m和n，因此QR分解的分布式方式通常可以满足要求。QR分解的分布式实现主要参考了文章[3]，这里简要介绍一下原理。

我们可以把一个分布式矩阵看做如下的形式：

那么矩阵A就可以分解为如下形式：

继续对等式右边的第二个矩阵进行QR分解：

上述等式右边的矩阵中最右边的矩阵R就是QR分解矩阵中的R，那么矩阵Q就可以通过矩阵A右乘以R的逆得到

根据以上公式我们可以看到，当把分布式的矩阵A划分成多个本地矩阵，并对每个本地矩阵进行QR分解，以及整合他们的R矩阵再进行QR分解就可以并行的获得最终的R矩阵。这也是QR分解的分布式实现的主要思想。

当然，如果整合的多个R矩阵依然比较大时，我们还可以继续借用这种思想。如下：

这里对整合的R矩阵进行分布式的QR分解。因此，按照上述思想继续分解，则A矩阵的QR分解最终将转化为如下形式：

下面是m行n列矩阵的分布式QR分解示意图

2.2 两个大型矩阵乘积的实现

解决了分布式矩阵的QR分解问题，接下来我们继续分析。由于QR分解的Q矩阵仍然是一个分布式矩阵，接下来这个矩阵（m行k+p列）将与原始矩阵（m行n列）进行乘积运算。由于这两个矩阵都非常大，这个过程将非常占内存，对于算法二来说，这种矩阵乘积方式将执行q次，比较耗时。对于算法三来说这个乘积过程不光要执行q次，同时每次还需要进行QR分解，这就会占用更多的时长。因此对于算法二和三来说，这个过程是实现的重点。考虑到算法二和三比算法一的准确率更高，使用场景更广，因此，tesla上的Randomized SVD算法主要根据算法二和三来实现。接下来我们重点阐述这里的实现过程。

算法二和三的过程基本上是相似的，都需要进行两个大型矩阵的乘积，在实际问题中，大型矩阵通常分为稀疏型和稠密型，稠密型的较多。由于稀疏与稠密矩阵的存储方式有较大差异，实现时也分为两种方式。先说稠密型矩阵，AT Q是两个大型矩阵乘积的表达形式（对于算法二，Q=AΩ）。实现时将A与Q都按行进行分布式存储，并根据矩阵乘积原理，将两个矩阵通过每行的索引采用join连接起来，再按A矩阵的列计算乘积结果的每一行。示意图如下：

对于稀疏型矩阵A，由于仅存储了存在的数值，占用内存较小，因此可以将整个矩阵A存放至本地，再根据等式AT Ω=(QT A)T ，求取Q矩阵的逆并以分布式的方式存储，然后将本地矩阵A广播至每个节点，与Q矩阵的逆的每一行进行乘积运算，在相乘时，仅需要遍历稀疏矩阵A上存在的数值，这极大地较小了时间复杂度。乘积的结果是一个本地矩阵，对本地矩阵进行转置即可获得结果。示意图如下：

根据稠密型与稀疏型矩阵的不同实现原理，我们可以看出，与稀疏型的计算方式相比稠密型的仅适用于行数与列数相对较小的矩阵，过大的行数与列数很容易造成内存溢出，这也是使用时要注意的地方。

2.3 本地矩阵B的SVD分解

以上步骤得到了矩阵A的近似基Q后，Stage1的任务已经完成。接下来我们只需要采用同样的分布式矩阵乘积方式计算B=QT A，得到本地矩阵B即可，然后对这个矩阵B采用常规的SVD分解等过程，得到最终的左奇异向量UA ，奇异值矩阵ΣA ，以及右奇异向量VA 。然而实际应用中往往没有这么顺利，这主要是由于矩阵A的列向量过大导致的。经过推导我们知道B是一个（k+p）行n列的矩阵，而k+p远远小于n，这样矩阵B将是一个行数远小于列数的矩阵，同时当n很大，例如5万维左右，直接对矩阵B进行SVD分解会因为计算格莱姆矩阵BT B导致内存溢出。直接计算是不行的，这里考虑将矩阵B进行转置，这样计算的格莱姆矩阵是（k+p）* （k+p），维数将大大减小，非常有利于计算接下来的特征值与特征向量。然而，矩阵B转置后的SVD分解不能直接用来计算最终的结果，我们还需要对其进行转化。推导如下，

如果A的SVD分解表达为： A=UΣVT ，则 AT =(UΣVT )T =VΣUT 可以看出，转置后的左、右奇异值向量将发生互换。

根据以上原理，我们得到矩阵B实际的左、右奇异值向量，再根据算法一中的6、7步得到最终矩阵A的各个分解结果。

3. Tesla上的Randomized SVD算法

3.1 算法使用介绍

Tesla上现已发布了Randomized SVD算法，见如下位置：

在使用方面，我们提供了样例，具体参考如下模块

该模块可以自动识别输入的矩阵是稀疏的还是稠密的，并据此采用2.2节中的方式进行计算。下面介绍使用中需要特殊注意的地方。

• 矩阵乘积的迭代类型 同上文所说，为了达到较为精确的结果，Randomized SVD模块根据算法二和三进行实现。因此，在tesla中，我们提供了两种迭代方式：QR与none，通过矩阵乘积的迭代类型来选择，其中，QR代表每轮的矩阵乘积过程中都要采用QR分解，即算法三的过程；none代表每轮矩阵乘积无需进行QR分解，是算法二的过程。如下图：

• 矩阵乘积的迭代轮数 同时，模块提供了迭代指数：矩阵乘积的迭代轮数，默认情况下选择“auto”，表示模块将根据奇异值个数k与Min(矩阵行数， 矩阵列数) * 0.1进行比较，k值较小，则迭代轮数为7，否则为4。上述规则是根据经验给定，后期可做相应修改。当然除了默认的情况外，用户也可以自己给定迭代轮数。

3.2 运行情况比较

在对比Randomized SVD算法与SVD算法的运行情况时，使用了两种类型的数据：稠密型与稀疏型。各配置如下：

1. 8001行1850列的稠密型矩阵，进行k值为800的矩阵分解，其中Randomized SVD算法的迭代类型选择none，过采样参数为5，迭代轮数为2。其他参数同SVD算法；
2. 760万行6万列的稀疏型矩阵，进行k值为2000的矩阵分解，其中Randomized SVD算法的迭代类型选择QR，过采样参数为10，迭代轮数为2。其他参数同SVD算法。在相同的资源配置下，运行时长的结果对比如下：

除了时长对比，在计算结果的数值对比上，稠密型矩阵的两种算法的计算结果相同。根据实验结果可以看到，与SVD算法相比，Randomized SVD算法在矩阵较大以及k值较大的情况下有很大优势。

致谢

感谢海量计算组kevinzwyou(游遵文)同学以及SNG业务安全组的sixiangpeng(彭思翔)、quintyqian(钱淑钗)同学为算法提出的宝贵建议以及在实现上带来的帮助。同时感谢海量计算组全体工作人员对本工作的支持。

参考文献

Nathan Halko, Per-Gunnar Martinsson, Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev., Survey and Review section, Vol. 53, num. 2, pp. 217-288, June 2011 V. Rokhlin, A. Szlam, and M. Tygert, A randomized algorithm for principal component analysis, SIAM J. Matrix Anal. Appl., 31 (2009), pp. 1100–1124. Paul G. Constantine,David F. Gleich.Tall and skinny QR factorizations in MapReduce architectures.ACM,June 08 - 08, 2011,pp.43-50