Cholesky分解的单Pass随机算法在秩亏最小二乘中的应用

The Application of Single-Pass Randomized Algorithm for Cholesky Decomposition in Rank Deficit Least Squares

摘要: 通常针对大规模的对称半正定矩阵求解最小二乘问题时，往往得不到精确解，而传统的矩阵分解方法在规模增大的情况下大大地增加了求解时间和复杂度，Cholesky分解的单pass随机算法只需对输入矩阵访问一次，因此本文将Cholesky分解的单pass随机算法与两种常见的求解最小二乘问题的方法结合，提出了两种新的求解最小二乘问题的算法，并将两者进行了比较，最后通过数值实验证明了这两种算法的有效性及可行性。

Abstract: Usually for large-scale symmetric semi-positive definite matrices to solve the least squares problem, often do not get the exact solution, and the traditional matrix factorization method greatly increases the solution time and complexity in the case of increasing scale, the Cholesky decomposition of the single-pass randomized algorithm only needs to access the input matrix once, so this paper combines the single-pass random algorithm of Cholesky decomposition with two common methods to solve the least squares problem, proposes two new algorithms for solving the least squares problem, compares the two, and finally proves the effectiveness and feasibility of the two algorithms through numerical experiments.

1. 问题的提出

随着计算机技术的不断发展，最小二乘问题的背后是一种优化方法，旨在找到一条直线(或更一般的一个模型)，以最小化观测数据点与该直线(或模型)之间的垂直距离的平方和[1]。在实际应用中，最小二乘问题通常涉及到大量数据的分析和模型拟合，以便从观测数据中提取出有用的信息并进行预测[2]。其背后的数学原理和算法在现代科学和工程中扮演着重要的角色[3]。

近年来，随机算法的不断完善，越来越多的人将随机算法应用到最小二乘问题的求解上面，并取得了不错的结果，在文献[4]中，Christos等人将随机投影的思想运用到求解非负最小二乘问题上去，在文献[5]中，Avron等人提出了Blendenpik算法，是一种基于随机算法求解最小二乘问题的算法。在文献[6]中，Drineas和Mahoney等人提出了两种随机算法求解最小二乘问题，一种是基于随机采样的思想，另一种是基于随机投影的思想。紧接着在文献[7]中，董银丽提出了一种基于随机采样方法求解最小二乘问题的新算法。

针对大规模的对称非正定矩阵分解来求解最小二乘问题，传统的矩阵分解算法至少需要读取输入矩阵两次。而面对存储在内存外部的大规模数据和流型数据时，往往只有一次读取数据的机会。在文献[8]中，提出的Cholesky分解的单pass随机算法是从随机矩阵投影的角度构造的一种新算法。本文尝试利用这种算法去求解对称半正定矩阵的最小二乘问题。根据求解最小二乘的方法，我们在本文第三章提出了两种新的算法用来求解秩亏最小二乘问题。

2. Cholesky分解的单Pass随机算法

任意给定一个对称半正定矩阵$A\in R^{n\times n}$ ，则我们可以对A进行Cholesky分解，将A表示成一个下三角矩阵L与L的转置的乘积，即$A=L{L}^{\text{T}}$ 。在大数据分析和机器学习中，Cholesky分解已经成为一种关键的分析工具。随着大数据时代的到来，传统的LU和QR算法在面对超大规模的对称半正定矩阵分解时，高内存消耗和高计算复杂度已无法满足时代需求，而单pass随机Cholesky分解只需对输入矩阵访问一次，而且有很高的误差精度，具体内容见算法1。

3. 用Cholesky分解的单Pass随机算法求解秩亏最小二乘问题

设矩阵A是一个$n\times n$ 阶矩阵，且b是一个n元列向量。秩亏最小二乘问题为：求解向量$x_{n\times 1}$ ，使得$\Vert Ax-b\Vert$ 达到最小。在本章中我们提出了两种求解最小二乘的方法，具体思想如下。

算法1. Cholesky分解的单pass随机算法(SPRCH)

输入：半正定矩阵$A\in {ℝ}^{n\times n}$ ，目标值r，过采样参数$p\ge 2$ ，$l=r+p$ ；

输出：置换矩阵$P\in {ℝ}^{n\times n}$ ，下三角矩阵$L\in {ℝ}^{n\times \left(r+p\right)}$ ，使得$\Vert PA{P}^{\text{T}}-L{L}^{\text{T}}\Vert \le \epsilon$ ，$\epsilon$ 是误差界。

生成一个$n\times \left(r+p\right)$ 阶的高斯采样矩阵${\Omega }_1$ ；

计算$Y_1=A{\Omega }_1$ ，$Y_2={\Omega }_1^{\text{T}}{Y}_1$ ；

对$Y_1$ 进行列主元LU分解，得到$P{Y}_1=L_y{U}_y$ ；

令$B={\left({\Omega }_1^{\text{T}}{P}^{\text{T}}{L}_y\right)}^{†}{Y}_2{\left(L_y^{\text{T}}P{\Omega }_1\right)}^{†}$ ；

对B作Cholesky分解，得$B=L_b{L}_b^{\text{T}}$ ；

令$L=L_y{L}_b$ ；

返回L，P。

第一种：在算法1中，我们通过Cholesky分解的随机算法得到$PA{P}^{\text{T}}\approx L{L}^{\text{T}}$ ，其中L是一个$n\times l$ 阶列满秩矩阵，U是一个$l\times n$ 阶矩阵。所以秩亏最小二乘问题可以转化为：

最小化$\Vert Ax-b\Vert \approx \Vert P^{\text{T}}L{L}^{\text{T}}Px-b\Vert =\Vert L{L}^{\text{T}}Px-Pb\Vert$ 。

我们可以令$y=L^{\text{T}}Px$ ，$c=Pb$ 。又因为矩阵L列满秩，所以$y=L^{†}c$ 。接下来再令$z=Px$ ，矩阵$L^{\text{T}}$ 不是列满秩矩阵，所以对$L^{\text{T}}$ 与z进行分块：

$L^{\text{T}}=\left(L_1,L_2\right)$ ，$z=\left(\begin{array}{l}{z}_1\hfill \\ z_2\hfill \end{array}\right)$

其中$L_1$ 是一个$l\times l$ 阶矩阵，$L_2$ 是一个$l\times \left(n-l\right)$ 阶矩阵。$z_1$ 是一个l元列向量，$z_2$ 是一个$\left(n-l\right)$ 元列向量。由此可知这个方程有无穷个解，为了简便起见，令$z_2=0$ ，$z_1=L_1^{-1}y$ 。最后通过$x=P^{\text{T}}z$ 可求出最小二乘解x。

第二种：秩亏最小二乘问题还可以转化为：

最小化$\Vert Ax-b\Vert \approx \Vert P^{\text{T}}L{L}^{\text{T}}Px-b\Vert =\Vert B{B}^{\text{T}}x-b\Vert$ ，我们可以令$B=P^{\text{T}}L$ ，$y=B^{\text{T}}x$ ，$c=Pb$ 。因为矩阵L列满秩则$B=P^{\text{T}}L$ 仍为列满秩矩阵，所以$y=B^{†}b$ 。接下来直接求解三角方程组$y=B^{\text{T}}x$ 求出最小二乘解x。以上两种算法具体如算法2和算法3：

算法2. 第一种方法求解秩亏最小二乘问题(SPRCH-LS1)

输入：分解矩阵$A\in {ℝ}^{n\times n}$ ，目标秩$r<n$ ，随机矩阵的列数$l=r+p$ ，其中p是过采样参数，$b\in {ℝ}^{n\times 1}$ ；

输出：正交置换矩阵$P\in {ℝ}^{n\times n}$ ，$x\in {ℝ}^{n\times 1}$ ，使得$\Vert Ax-b\Vert \le \epsilon$ ，$\epsilon$ 表示一个比较小的误差界；

1) 运用算法1对A做分解得到L, P；

2) 计算$c=Pb$ ；

3) 计算$y=L^{†}c$ ；

4) 计算$L^{\text{T}}$ ，并对$L^{\text{T}}$ 作截断得到$L_1$ ；

5) 计算$z_1=L_1^{-1}y$ ，令$z_2=0$ ；

6) 计算$x=P^{\text{T}}z$ ；

7) 返回x。

算法3. 第二种方法求解秩亏最小二乘问题(SPRCH-LS2)

输入：分解矩阵$A\in {ℝ}^{n\times n}$ ，目标秩$r<n$ ，随机矩阵的列数$l=r+p$ ，其中p是过采样参数，$b\in {ℝ}^{n\times 1}$ ；

输出：正交置换矩阵$P\in {ℝ}^{n\times n}$ ，$x\in {ℝ}^{n\times 1}$ ，使得$\Vert Ax-b\Vert \le \epsilon$ ，$\epsilon$ 表示一个比较小的误差界；

1) 运用算法1对A做分解得到L，P；

2) 计算$B=P^{\text{T}}L$ ；

3) 计算$y=B^{†}b$ ；

4) 计算$x=B^{\text{T}}\backslash y$

5) 返回x。

4. 数值实验

4.1. 对奇异值慢速衰减的固定输入矩阵

本节所处理的矩阵是随机生成的奇异值慢速衰减的半正定矩阵。生成方式如下：

1) 生成$\Sigma =\mathrm{diag}\left(1,\cdots ,1,{\text{e}}^{-s},{\text{e}}^{-2s},\cdots ,{\text{e}}^{-\left(n-r\right)s}\right)$ ，这里的r和s都是大于0的实数，r用来控制输入矩阵的秩，s用来控制衰减速度，令$s=190/\left(n-r-1\right)$ ，让输入矩阵A的奇异值从1到${\text{e}}^{-200}$ 呈指数型慢速下降；

2) 生成一个$n\times n$ 的随机正交矩阵U；

3) 令待分解矩阵$A=U\Sigma U^{\text{T}}$ 。

(a) (b)

Figure 1. Error and time comparison of Algorithm 2 and Algorithm 3 under a 2000 × 2000 dimensional matrix with slowly decaying singular values. (a) Error comparison, (b) Time comparison

图1. 算法2和算法3在2000 × 2000维且奇异值缓慢衰减的矩阵下的误差和时间比较。(a) 误差比较，(b) 时间比较

从图1中我们可以看出，随着秩的增加，这两种算法的运行时间都在增大，但是误差的变换趋势相同，误差精度都保持在10⁻⁶左右。

4.2. 对奇异值快速衰减的固定输入矩阵

与4.1节类似，这里我们令$s=900/\left(n-r-1\right)$ ，让输入矩阵A的奇异值从1到${\text{e}}^{-1000}$ 呈指数型快速下降。

(a) (b)

Figure 2. Comparison of error and time between Algorithm 2 and Algorithm 3 under a 2000 × 2000 dimensional matrix with fast decay of singular values. (a) Error comparison, (b) Time comparison

图2. 算法2和算法3在2000 × 2000维且奇异值快速衰减的矩阵下的误差和时间比较。(a) 误差比较，(b) 时间比较

从图2可以看出在奇异值快速衰减的条件下，随着秩的增加，这两种算法的误差和运行时间基本相似。

5. 总结与展望

5.1. 全文总结

本文所提出的这两种新算法在对处理大规模对称半正定矩阵的最小二乘问题具有很大的优势。这两种算法的误差精度和运行时间成本都很好，其次我们发现当输入矩阵的奇异值慢速衰减时，算法2略有优势，而当输入矩阵的奇异值快速衰减时这两种算法的处理结果很接近。

5.2. 本文展望

本文主要提出了两种关于Cholesky分解的单pass随机算法在秩亏最小二乘中应用的新算法，这两种新算法为求解大规模对称半正定矩阵的最小二乘问题提供了合适的方案，但仍存在不足，比如除了以上两种求解最小二乘的思想，是否还能将Cholesky分解的单pass随机算法在运用到其他的求解二乘方法中，希望在以后的研究中可以给出更多的结果。