Pass-efficient randomized LU algorithms for computing low-rank matrix approximation

A Linear algebra basics

We briefly review the definitions and properties of the orthogonal projection and the pseudoinverse [10].

B SVD for low-rank approximation problem

Suppose the SVD decomposition for matrix 𝐀 ∈ ℝ m × n , m ≥ n is the following:

where 𝐔 ∈ ℝ m × n has orthonormal columns, 𝐕 ∈ ℝ n × n is an orthogonal matrix, and 𝚺 ∈ ℝ n × n is a diagonal matrix whose diagonal entries are the singular values of 𝐀 in descending order. Then the optimal solution to the fixed rank problem (1.1) for a given rank k is 𝐀 k = 𝐁𝐂 , where we set 𝐁 = 𝐔 ( : , 1 : k ) 𝚺 ( 1 : k , 1 : k ) and 𝐂 = 𝐕 ( : , 1 : k ) T .^[3] This gives us the explicit matrix norms of our rank k approximation errors as

C Randomized SVD Algorithm

Formally, for a given matrix 𝐀 ∈ ℝ m × n , we multiply it with a random Gaussian matrix 𝛀 ∈ ℝ n × l , where l = k + q , k is the target rank and q is the oversampling parameter. Then 𝐘 = 𝐀 ⁢ 𝛀 will capture the most action approximately (i.e., capture the space spanned by the first l singular vector approximately) of 𝐀 and more importantly the dimension of the space spanned by the columns of 𝐘 is much smaller than the column space of 𝐀 . Suppose 𝐐 ∈ ℝ m × l is an orthonormal basis for the column space of 𝐘 . Then we have the following random QB low-rank approximation

Here 𝐁 = 𝐐 T ⁢ 𝐀 ∈ ℝ l × n , which is a smaller matrix than 𝐀 on which we could perform some standard factorization such as the SVD.

For a matrix whose singular values decay slowly, the above procedure may not provide good results. However, the power iteration ( 𝐀𝐀 T ) p ⁢ 𝐀 can be applied to solve this problem, where p is the exponent of the power iteration. Suppose 𝐀 has the SVD decomposition (B.1). Then

Equation (C.1) implies that the singular values of the matrix ( 𝐀𝐀 T ) p ⁢ 𝐀 decay faster than those of 𝐀 . However, both of them have the same left and right singular vectors. A basic randomized SVD algorithm (RandSVD) with subspace iteration is shown here as Algorithm 9 (see [13]). In practice, to remove the rounding errors introduced by (C.1), we use Algorithm 10 in practice to compute it [13, 36]. An optimized version of reorthogonalization was proposed in [17], where the authors used LU factorization to replace QR every time except the last iteration. As we cans see that for the RandSVD Algorithm 9, we need to access the matrix an even number times. In some cases, the input matrix is enormous, and accessing the matrix is very expensive. Randomized SVD algorithms that minimize accesses of the matrix have been proposed in [1], where the authors evaluated the co-range sketch of ( 𝐀 T ⁢ A ) p ⁢ 𝛀 for p ≥ 1 .

The idea for the generalized subspace iteration is very simple. Let p = ⌊ v - 1 2 ⌋ when v is odd in Algorithm 9. We modify steps 1–3 of Algorithm 9 when v ≥ 2 is even as follows:

1. Generate a random Gaussian matrix 𝛀 ∈ ℝ n × l .

2. Compute ( 𝐀𝐀 T ) p ⁢ 𝐀 ⁢ 𝛀 and its orthonormal basis 𝐕 ;

Combine above discussion, we introduce the generalized subspace iteration here as Algorithm 11.

Algorithm 9 (A Basic Randomized SVD Algorithm.).

Algorithm 10 (Power Iteration Reorthogonalization for Algorithm 9.).

Algorithm 11 (Generalized Subspace Iteration for Algorithm 9.).