Iterative refinement for singular value decomposition based on matrix multiplication
基于矩阵乘法的奇异值分解迭代细化

Let $A$ be a real $m\times n$ matrix. In this paper, we propose a refinement algorithm for the singular value decomposition (SVD) of $A$ with $m\ge n$. If $m<n$, considering the SVD of $A^{\mathrm{T}}$ yields equivalent results. It is well known that the SVD has many applications in various fields, such as signal processing [1], [2], statistical analysis [3], [4], and so forth. Excellent overviews of the SVD can be found in [5], [6].
设 $A$ 为实 $m\times n$ 矩阵。在本文中，我们提出了一种 with $m\ge n$ 的 $A$ 奇异值分解 （SVD） 的细化算法。如果 $m<n$ ，考虑 的 SVD $A^{\mathrm{T}}$ 会产生等效的结果。众所周知，SVD 在各个领域都有许多应用，例如信号处理、统计分析等。SVD 的精彩概述可以在 、 中找到。

Throughout the paper, let $I_n$ and $O$ denote the $n\times n$ identity matrix and the zero matrix of appropriate size, respectively. Moreover, $\Vert \cdot \Vert$ denotes the spectral norm for matrices. If necessary, we distinguish between the approximate quantities and the computed results, e.g., for some quantity $\alpha$, we write $\tilde{\alpha }$ and $\hat{\alpha }$ as an approximation of $\alpha$ and a computed result for $\alpha$, respectively.
在整篇论文中，let $I_n$ 和 $O$ 分别表示适当大小的 $n\times n$ 单位矩阵和零矩阵。此外， $\Vert \cdot \Vert$ 表示矩阵的谱范数。如有必要，我们区分近似量和计算结果，例如，对于某个量 $\alpha$ ，我们分别将 和 $\hat{\alpha }$ 写 $\tilde{\alpha }$ 为 $\alpha$ 的近似 $\alpha$ 值和 计算结果 。

Let ${\sigma }_i\in \mathbb{R}$, $i=1,\dots ,n$, denote the singular values of $A$. We consider the (full size) SVD of $A$ such that $A=U\Sigma V^{\mathrm{T}},U\in {\mathbb{R}}^{m\times m},V\in {\mathbb{R}}^{n\times n},\Sigma \in {\mathbb{R}}^{m\times n},$where both $U$ and $V$ are orthogonal and $\Sigma$ is diagonal with ${\Sigma }_{ii}={\sigma }_i$. For simplicity, we assume that ${\sigma }_1>{\sigma }_2>\cdots >{\sigma }_n>0.$In other words, we consider the case that all singular values are simple, and $A$ has full column rank. If there are multiple or nearly multiple singular values, we need some special care as in [7], [8].
设 ${\sigma }_i\in \mathbb{R}$ ， $i=1,\dots ,n$ ， 表示 的 $A$ 奇异值。我们考虑 的（全尺寸）SVD 是 $A$ 这样的，其中 $U$ 和 $V$ 都是正交的，并且 $\Sigma$ 是与 ${\Sigma }_{ii}={\sigma }_i$ 的对角线。为简单起见，我们假设换句话说，我们考虑所有奇异值都是简单的，并且 $A$ 具有完整的列排名。如果存在多个或几乎多个奇异值，则需要特别小心，如 ， 。

Recently, the authors proposed refinement algorithms for symmetric eigenvalue decomposition in [7], [8]. In the same spirit of the previous papers, the use of higher-precision arithmetic in our proposed refinement algorithm for the SVD is primarily restricted to matrix multiplication, which accounts for most of the computational cost. There are several approaches for higher-precision matrix multiplication. For example, XBLAS (extra precise BLAS) [9] and fast and accurate algorithms for dot products [10] and matrix products [11] based on error-free transformations are available for efficient implementation.
最近，作者在 [7]、[8] 中提出了对称特征值分解的改进算法。本着与前几篇论文相同的精神，在我们提出的 SVD 细化算法中，更高精度算术的使用主要限于矩阵乘法，这占了大部分计算成本。有几种方法可以进行更高精度的矩阵乘法。例如，XBLAS （超精确 BLAS） [9] 和基于无差错变换的点积 [10] 和矩阵积 [11] 的快速准确算法可用于高效实现。

The idea of our algorithm is to use the following relations: (1)$U^{\mathrm{T}}U=I_m\text{(orthogonality of }U\text{)}$(2)$V^{\mathrm{T}}V=I_n\text{(orthogonality of }V\text{)}$(3)$U^{\mathrm{T}}AV=\Sigma \text{(diagonality of }A\text{ as the SVD)}$ Using these relations, we develop a refinement algorithm for the SVD in the same manner as Newton’s method. Thus, the proposed algorithm has quadratic convergence.
我们算法的思路是使用以下关系式：（1）$U^{\mathrm{T}}U=I_m\text{（U 的 正交性 }\text{）（}$2）$V^{\mathrm{T}}V=I_n\text{（V 的 正交性 }\text{）}$（3）$\text{Unexpected text node: '（'}$利用这些关系，我们以与牛顿方法相同的方式为 SVD 开发了一种细化算法。因此，所提出的算法具有二次收敛性。

There exist several refinement algorithms for SVD that are based on Newton’s method for nonlinear equations (cf. e.g., [12]). Since this sort of algorithm is designed to improve a triplet $(\sigma ,u,v)\in \mathbb{R}\times {\mathbb{R}}^n\times {\mathbb{R}}^n$ individually, where $\sigma$ is a singular value and $u$ and $v$ are corresponding left and right singular vectors, applying such an approach to all triplets requires $\mathcal{O}\left(n^4\right)$ arithmetic operations. In [13], Davies and Smith proposed an iterative refinement algorithm for updating the singular value decomposition in $\mathcal{O}\left(n^3\right)$ operations. However, similarly to the Davies–Modi algorithm [14] for the symmetric eigenvalue decomposition as mentioned in the previous paper [7], the Davies–Smith algorithm has the limitation of achievable accuracy of the results. The reason is as follows. The Davies–Smith algorithm assumes that a given real matrix $A$ is preconditioned to a nearly diagonal matrix such as ${\hat{U}}^{\mathrm{T}}A\hat{V}$, where $\hat{U}$ and $\hat{V}$ are computed SVD factors, i.e., $\hat{U}$ and $\hat{V}$ are approximately orthogonal matrices. Since $\hat{U}$ and $\hat{V}$ involve numerical errors, the matrix multiplications in ${\hat{U}}^{\mathrm{T}}A\hat{V}$ are generally not orthogonal transformations, and the singular values of ${\hat{U}}^{\mathrm{T}}A\hat{V}$ are slightly perturbed from the original matrix $A$. Then, the singular vectors are also perturbed. Therefore, even if the Davies–Smith algorithm provides accurate singular vectors of ${\hat{U}}^{\mathrm{T}}A\hat{V}$, they are not necessarily accurate ones of $A$. On the other hand, our proposed algorithm uses the original matrix $A$ for obtaining accurate singular vectors of $A$.
SVD 有几种基于牛顿非线性方程法的细化算法（参见 e.g.）。由于这种算法旨在单独改进三元组 $(\sigma ,u,v)\in \mathbb{R}\times {\mathbb{R}}^n\times {\mathbb{R}}^n$ ，其中 $\sigma$ 是奇异值 和 $u$ 和 $v$ 是相应的左和右奇异向量，因此将这种方法应用于所有三元组需要 $\mathcal{O}\left(n^4\right)$ 算术运算。在 中，Davies 和 Smith 提出了一种迭代细化算法，用于更新运算中的 $\mathcal{O}\left(n^3\right)$ 奇异值分解。然而，与上一篇文章中提到的对称特征值分解的 Davies-Modi 算法类似，Davies-Smith 算法在结果的可实现精度方面存在局限性。原因如下。Davies-Smith 算法假设给定的实数矩阵 $A$ 被预设为接近对角线的矩阵，例如 ${\hat{U}}^{\mathrm{T}}A\hat{V}$ ，其中 $\hat{U}$ 和 $\hat{V}$ 是计算的 SVD 因子，即 $\hat{U}$ 和 $\hat{V}$ 是近似正交矩阵。由于 $\hat{U}$ 和 $\hat{V}$ 涉及数值误差，因此 中的 ${\hat{U}}^{\mathrm{T}}A\hat{V}$ 矩阵乘法通常不是正交变换，并且 的 ${\hat{U}}^{\mathrm{T}}A\hat{V}$ 奇异值与原始矩阵 $A$ 略有不同。然后，奇异向量也会受到扰动。因此，即使 Davies-Smith 算法提供了 的 ${\hat{U}}^{\mathrm{T}}A\hat{V}$ 精确奇异向量 ，它们也不一定是 的 $A$ 准确向量。另一方面，我们提出的算法使用原始矩阵 $A$ 来获取 的 $A$ 精确奇异向量。

The rest of the paper is organized as follows. In Section 2, we present a refinement algorithm for the SVD. In Section 3, we provide a convergence analysis of the proposed algorithm. In Section 4, we present some numerical results showing the behavior and performance of the proposed algorithm.
本文的其余部分组织如下。在第 2 节中，我们提出了 SVD 的优化算法。在第 3 节中，我们提供了所提出的算法的收敛分析。在第 4 节中，我们提供了一些数值结果，显示了所提出的算法的行为和性能。

Let $\hat{U}\in {\mathbb{R}}^{m\times m}$ and $\hat{V}\in {\mathbb{R}}^{n\times n}$ be given approximation of $U$ and $V$, respectively. Let further $F\in {\mathbb{R}}^{m\times m}$ and $G\in {\mathbb{R}}^{n\times n}$ be correction matrices satisfying $U=\hat{U}(I_m+F)$ and $V=\hat{V}(I_n+G)$, respectively. Let $ϵ$ be defined as (4)$ϵ≔max({ϵ}_F,{ϵ}_G),{ϵ}_F≔\Vert F\Vert ,{ϵ}_G≔\Vert G\Vert .$We assume that $ϵ<1$. Then, both $I_m+F$ and $I_n+G$ are nonsingular, and ${(I_m+F)}^{-1}=I_m-F+{\Delta }_F,{\Delta }_F≔\sum _{k=2}^{\infty }{(-F)}^k,\Vert {\Delta }_F\Vert \le \frac{{ϵ}_F}{1-{ϵ}_F},{(I_n+G)}^{-1}=I_n-G+{\Delta }_G,{\Delta }_G≔\sum _{k=2}^{\infty }{(-G)}^k,\Vert {\Delta }_G\Vert \le \frac{{ϵ}_G}{1-{ϵ}_G}.$ Inserting $U=\hat{U}(I_m+F)$ into (1), we have $(I_m+F^{\mathrm{T}}){\hat{U}}^{\mathrm{T}}\hat{U}(I_m+F)=I_m$and ${\hat{U}}^{\mathrm{T}}\hat{U}={(I_m+F^{\mathrm{T}})}^{-1}{(I_m+F)}^{-1}=(I_m-F^{\mathrm{T}}+{\Delta }_F^{\mathrm{T}})(I_m-F+{\Delta }_F),$which yields (5)$F+F^{\mathrm{T}}=I_m-{\hat{U}}^{\mathrm{T}}\hat{U}+{\Delta }_1,{\Delta }_1≔{\Delta }_F+{\Delta }_F^{\mathrm{T}}+{(F-{\Delta }_F)}^{\mathrm{T}}(F-{\Delta }_F).$Similarly, inserting $V=\hat{V}(I_n+G)$ into (2), we have (6)$G+G^{\mathrm{T}}=I_n-{\hat{V}}^{\mathrm{T}}\hat{V}+{\Delta }_2,{\Delta }_2≔{\Delta }_G+{\Delta }_G^{\mathrm{T}}+{(G-{\Delta }_G)}^{\mathrm{T}}(G-{\Delta }_G).$Moreover, inserting $U=\hat{U}(I_m+F)$ and $V=\hat{V}(I_n+G)$ into (3), we have (7)$\Sigma -F^{\mathrm{T}}\Sigma -\Sigma G={\hat{U}}^{\mathrm{T}}A\hat{V}+{\Delta }_3,{\Delta }_3≔-\Sigma {\Delta }_G-{\Delta }_F^{\mathrm{T}}\Sigma -{(F-{\Delta }_F)}^{\mathrm{T}}\Sigma (G-{\Delta }_G).$Here, (8)$\Vert {\Delta }_1\Vert \le \frac{(3-2{ϵ}_F){ϵ}_F^2}{{(1-{ϵ}_F)}^2}\le \chi \left(ϵ\right){ϵ}^2,$(9)$\Vert {\Delta }_2\Vert \le \frac{(3-2{ϵ}_G){ϵ}_G^2}{{(1-{ϵ}_G)}^2}\le \chi \left(ϵ\right){ϵ}^2,$(10)$\Vert {\Delta }_3\Vert \le \frac{{ϵ}_F^2+{ϵ}_G^2+(1-{ϵ}_F-{ϵ}_G){ϵ}_F{ϵ}_G}{(1-{ϵ}_F)(1-{ϵ}_G)}\Vert \Sigma \Vert \le \chi \left(ϵ\right){ϵ}^2\Vert A\Vert ,$ where
设 $\hat{U}\in {\mathbb{R}}^{m\times m}$ 和 $\hat{V}\in {\mathbb{R}}^{n\times n}$ 分别得到 $U$ 和 $V$ 的近似值。设 further $F\in {\mathbb{R}}^{m\times m}$ 和 $G\in {\mathbb{R}}^{n\times n}$ 分别是满足 $U=\hat{U}(I_m+F)$ 和 $V=\hat{V}(I_n+G)$ 的校正矩阵。设 $ϵ$ 定义为 我们假设 $ϵ<1$ 。那么，和 都是 $I_m+F$ 奇异的，并且 插入 $U=\hat{U}(I_m+F)$ 到 ，我们有 和 ，得到 同样，插入 $V=\hat{V}(I_n+G)$ 到 ，我们有 此外，插入 $U=\hat{U}(I_m+F)$ 和 $V=\hat{V}(I_n+G)$ 进入 ，我们有 这里 $I_n+G$ ，其中(11)$\chi \left(ϵ\right)≔\frac{3-2ϵ}{{(1-ϵ)}^2}.$

Omitting the second-order terms ${\Delta }_1$, ${\Delta }_2$, and ${\Delta }_3$ from (5), (6), and (7) in a similar way to Newton’s method, we obtain a system of matrix equations for $\tilde{F}=\left({\tilde{f}}_{ij}\right)\in {\mathbb{R}}^{m\times m}$, $\tilde{G}=\left({\tilde{g}}_{ij}\right)\in {\mathbb{R}}^{n\times n}$, and $\tilde{\Sigma }=\mathrm{diag}\left({\tilde{\sigma }}_i\right)\in {\mathbb{R}}^{m\times n}$ as (12)$\left\{\begin{array}{c}\tilde{F}+{\tilde{F}}^{\mathrm{T}}=R,R≔I_m-{\hat{U}}^{\mathrm{T}}\hat{U}\hfill \\ \tilde{G}+{\tilde{G}}^{\mathrm{T}}=S,S≔I_n-{\hat{V}}^{\mathrm{T}}\hat{V}\hfill \\ \tilde{\Sigma }-{\tilde{F}}^{\mathrm{T}}\tilde{\Sigma }-\tilde{\Sigma }\tilde{G}=T,T≔{\hat{U}}^{\mathrm{T}}A\hat{V}\hfill \end{array}\right.$(13)$\iff \left\{\begin{array}{cc}{\tilde{f}}_{ij}+{\tilde{f}}_{ji}=r_{ij}\hfill & \text{for}1\le i,j\le m\hfill \\ {\tilde{g}}_{ij}+{\tilde{g}}_{ji}=s_{ij}\hfill & \text{for}1\le i,j\le n\hfill \\ {\tilde{\Sigma }}_{ij}-{\tilde{\sigma }}_j{\tilde{f}}_{ji}-{\tilde{\sigma }}_i{\tilde{g}}_{ij}=t_{ij}\hfill & \text{for}1\le i\le m,1\le j\le n\hfill \end{array}\right..$ All that remains is to solve (12) for $\tilde{F}$, $\tilde{G}$, and $\tilde{\Sigma }$.
以与牛顿方法类似的方式省略 （5）、（6） 和 （7） 中的二阶项 ${\Delta }_1$ 、 ${\Delta }_2$ 和 ${\Delta }_3$ ，我们得到一个矩阵 $\tilde{F}=\left({\tilde{f}}_{ij}\right)\in {\mathbb{R}}^{m\times m}$ 方程组 ， $\tilde{G}=\left({\tilde{g}}_{ij}\right)\in {\mathbb{R}}^{n\times n}$ $\tilde{\Sigma }=\mathrm{diag}\left({\tilde{\sigma }}_i\right)\in {\mathbb{R}}^{m\times n}$ 因为 (12)$\left\{\begin{array}{c}\tilde{F}+{\tilde{F}}^{\mathrm{T}}=R,R≔I_m-{\hat{U}}^{\mathrm{T}}\hat{U}\\ \tilde{G}+{\tilde{G}}^{\mathrm{T}}=S,S≔I_n-{\hat{V}}^{\mathrm{T}}\hat{V}\\ \tilde{\Sigma }-{\tilde{F}}^{\mathrm{T}}\tilde{\Sigma }-\tilde{\Sigma }\tilde{G}=T,T≔{\hat{U}}^{\mathrm{T}}A\hat{V}\end{array}\right.$(13)$\iff \left\{\begin{array}{cc}{\tilde{f}}_{ij}+{\tilde{f}}_{ji}=r_{ij}& \text{for}1\le i,j\le m\\ {\tilde{g}}_{ij}+{\tilde{g}}_{ji}=s_{ij}& \text{for}1\le i,j\le n\\ {\tilde{\Sigma }}_{ij}-{\tilde{\sigma }}_j{\tilde{f}}_{ji}-{\tilde{\sigma }}_i{\tilde{g}}_{ij}=t_{ij}& \text{for}1\le i\le m,1\le j\le n\end{array}\right..$ 剩下的就是求解 $\tilde{F}$ 、 和 $\tilde{G}$ $\tilde{\Sigma }$ 的 （12） 。