An Accelerated Accuracy-enhanced Randomized Singular Value Decomposition for Factorizing Matrices with Low-rank Structure

Big data has in recent years gained ground in many scientific and engineering problems. It seems to some extent prohibitive for traditional matrix decomposition methods (i.e. QR, SVD, EVD, etc.) to handle such large-scale problems involving data matrix. Many researchers have developed several algorithms to decompose such big data matrices. An accuracy-enhanced randomized singular value decomposition method (referred to as AE-RSVDM) with orthonormalization recently becomes the state-of-the-art to factorize large data matrices with satisfactory speed and accuracy. In our paper, low-rank matrix approximations based on randomization are studied, with emphasis on accelerating the computational efficiency on large data matrices. By this, we accelerate the AE-RSVDM with modified normalized power iteration to result in an accelerated version. The accelerated version is grounded on a two-stage scheme. The first stage seeks to find the range of a sketch matrix which involves a Gaussian random matrix. A low-dimensional space is then created from the high-dimensional data matrix via power iteration. Numerical experiments on matrices of different sizes demonstrate that our accelerated variant achieves speedups while attaining the same reconstruction error as the AE-RSVDM with orthonormalization. And with data from Google art project, we have made known the computational speed-up of the accelerated variant over the AE-RSVDM algorithm for decomposing large data matrices with low-rank form.


Introduction
Matrices in many applications appear with a low-rank structure and matrix factorization approaches are often employed to develop compacted and informative observations to make data calculation and interpretation easier. A long-established method is the almighty singular value decomposition (SVD) which discovers the finest low-rank approximation of a data matrix. Various SVD algorithms in scientific computing depend on low-rank approximation via matrix decomposition of the form This illustration permits data analysts to analyze or work with the matrix D using the factor matrices E and F as an alternative of the full matrix, which is more efficient both computationally and memory usage. Besides, these smaller factor matrices can offer definite structures to achieve the desired result by analyzing a data matrix [1]. SVD provides a numerically robust matrix factorization which can facilitate to calculate the pseudo-inverses of general matrices [2], get low-rank approximations [3] and evaluate the least-squares and minimum-norm outcomes for a linear model [4]. Besides, the SVD is the algorithm behind several machine learning concepts, including sparse coding, matrix completion, dictionary learning, PCA and robust PCA. For a Factorizing Matrices with Low-rank Structure comprehensive technical outline of SVD, one can refer to the work of Trefethen et. al. [5].

Some Applications of Low-rank Approximation
Many applications in scientific computing encompass matrices with low-rank structure. Singular value decomposition is employed in the following applications: A basic technique in statistics is to calculate directions of highest variance in vector-valued data by performing PCA on an input matrix [6].
Singular value decomposition also plays an essential part in tensor decomposition [7].
Another on the list of applications includes image compression [8].

Early Works of Randomized Algorithm
Recent algorithms for performing singular value decomposition are grounded on Krylov approaches, like Lanczos algorithm [9]. These approaches are precise and are particularly prevailing for approximating structured, and sparse matrices.
In the last 20 years, random matrix factorization algorithms have become more common to measure low-rank matrix approximation. Introducing the "Monte Carlo" SVD, a strenuous approach was used in the paper of Frieze, Kannan and Vempala to competently measure the low-rank SVD based on a uniform column sample row [10]. Several works used a random projection based on a robust method [11,12]. Specifically, the features of random vectors are dominated so that the range of a matrix can be captured effectively. Woolfe et al. further improved computational performance by using the characteristics of very structured matrices, that allow matrix multiplication speed [13]. The significant effort by Halko et al. united and extended previous work on randomized SVD and presented advanced prototype algorithms to calculate low-rank singular value decomposition which is closer to the optimum value [3]. The algorithmic work of Halko et al. is summarized in algorithm 1 below: ; U QU = To improve on accuracy, Martinsson introduced the "Accuracy-enhanced randomized SVD with orthonormalization (AE-RSVDM)" [14] which is displayed in algorithm 2. (2) ( ); ; U QU =

Main Contribution
Inspired by the work of Martinsson [14], we found it attractive to improve upon the speed of the AE-RSVDM (with orthonormalization). This manuscript focuses on a conceptually workable algorithm for accelerating the computational performance of the numerical low-rank approximation of singular value decomposition (SVD) via random projection. The algorithms presented in this manuscript are obtained by a significant disparity on the normalized power iteration scheme [15] and are easy to implement. We present the so-called "accelerated AE-RSVDM with modified normalized power iteration" algorithm which is an extension of Martinsson's work. Section 4 of this work shows several numerical experiments on both 'fat' and 'thin' matrices and comparing the performance with the AE-RSVDM (with orthonormalization) [14]. With the accelerated version of the AE-RSVDM, efficiency is improved by the use of LU factorization (with partial pivoting) in step 2 of the modified normalized power iteration instead of QR decomposition to re-normalize after each step of the first q applications of Y , the sketch matrix. The precision is also observed because the reconstruction errors are not significantly different.

Notation
Throughout this manuscript, all vectors are measured in n R using their 2 l norm (Euclidean norm), Matrices are denoted with capital letters such as , .
A Q For simplicity, we prefer to work with only real matrices in this manuscript e.g., respectively [14].

Orthonormalization
Given an m × l matrix B , with m ≥ l , we introduce a function ( ) , to signify orthonormalization of the columns of B . In other words, Q will be an m × l orthonormal matrix whose columns make a basis for the range of B [14]. In practice, the keyword 'qr' means to use a householder transformation to perform a QR factorization; e.g., in MATLAB, we would write ( )

Conceptual Overview of Singular Value Decomposition
Every real matrix of arbitrary rectangular dimensions can be factored into the form where and , U V both have orthonormal columns, and whose columns hold the left and right singular vectors, respectively. The matrix ∑ is diagonal and contains singular values which match the singular vectors. By convention, we denote the th i singular value as , and arrange them and their corresponding singular vectors into

Geometrical Interpretation of SVD
Proof: We know that We know that .
Geometrically, U and V can be assumed as rotation matrices

Randomized Algorithm Enhanced
The basis matrix Q usually fails to deliver a good estimation for the range of the input matrix. The reason being that, most real-world data matrices do not have an exact rankk , but instead have a range of singular values that gradually decay [14]. Efficiency can be significantly improved with the concept of over-sampling [3] and power repetition as suggested in Martinsson's paper [14].

Oversampling
Majority of data matrices do not have an exact rank, which is normally enough to achieve a strong basis similar to the best possible one [3].

The Idea of Power Iteration
The idea of power sample iterations is another technique for enhancing the value of the basis Q . The method outlined in algorithm 3 works well for matrices with some kind of decay in singular values but they can generate a weak basis for a flat singular spectrum or a very big input matrix. The main clue about power iteration is for enhancing the exactness of random algorithms, which were originally proposed in the manuscript of Rokhlin et. al. [15]. To review, presumably, we have calculated a rankk approach to an m n The concept demonstrates that the spectral norm error is restricted to a factor that measures ( )

The Modified Normalized Power Iteration Algorithm
Martinsson's algorithm tends to be accurate but time-wasting to compute the approximated rankk SVD. We therefore, found it necessary to improve upon speed and also maintain the accuracy of the RSVD via the modified normalized power iteration. We therefore, introduce the modified normalized power iteration algorithm.

Algorithm 3: Modified Normalized Power Iteration
Inputs: An input matrix m n A × ∈ R and the sketch matrix

The Accelerated AE-RSVDM Algorithm
The prototype randomized scheme [3] gives accurate results for matrices which rapidly decay with respect to its singular value but tends to produce sub-optimal results when they do not [14]. We therefore, introduce the accelerated AE-RSVDM with modified normalized power iteration algorithm. We give a quick overview of the accelerated variant algorithm. As usual, the two-stage scheme was used: Stage A-The range finder: For us to obtain a low-rank SVD approximation of and m n ≥ , we require a random test matrix . By this, p denotes a parameter for oversampling, and we will see that a small amount of oversampling is often sufficient in practice.
The definition of a random projection is widely used to sample the range of the data matrix A to construct such an orthogonal projection scheme efficiently. Random projections are orthodox data and generated from the standard normal distribution Equation (10) can be efficiently executed in parallel.
Therefore, let us generate the random test matrix n× Ω ∈ l R , which is indeed generated from the standard normal distribution, whose columns are the vectors { } 1 , is then acquired by the random test matrix after the input matrix has been multiplied : .
Once Y is found, it only remains to use the modified normalized power iteration in algorithm 3 to renormalize and orthonormalize the sketch (sample) matrix to update the sketch matrix such that We now aim to find a smaller matrix . This means we are therefore constructing a low-dimensional space from a high-dimensional data matrix such that : .
We now compute the SVD of the smaller matrix B to obtain two orthonormal matrices The left singular vectors are then updated by multiplying the orthonormal basis Q and the truncated matrix U ; ; U QU = Efficiency is enhanced by the use of LU factorization (with partial pivoting) in step 2 of algorithm 3 instead of QR decomposition to re-normalize after each step of the first q applications of , Y the sketch matrix. This normalization ensures speedup and maintains accuracy as compared to the AE-RSVDM (with orthonormalization) as proposed by Martinsson [14]. This is summarized in algorithm 3.

Analysis of Computational Complexity
We analyze the significance of each step to contextualize the computational complexity of both the accelerated AE-RSVDM with modified normalized power iteration and the AE-RSVDM (with orthonormalization) in floating-point operation. The computational complexities are summarized in Table 1

Theoretical Error Bound of the Accelerated AE-RSVDM with Modified Normalized Power Iteration
To state the theorem of the error bound of algorithm 4, we make reference to two most important propositions about random measurement matrix [17][18][19].
Here, the operator E refers to the expectation in relation to the Gaussian measurement matrix Ω , and e is known to be Euler's number. Furthermore, p the over-sampling parameter is presumed to be larger than or equal to 2. From the error bound in (18), it follows that the approximation error can be controlled both by the over-sampling parameter, p and the power iteration scheme parameter, q . The power iteration scheme parameter q speeds up the decay rate in the singular values of the sampled matrix by keeping the same eigenvectors [14].

Proof:
We proceed with the proof. Partition the SVD of A as follows: From theorem (2.9.1) of [3], assuming 1 Ω is not singular, it holds that Here, ⋅ represents either the 2 l -operator norm, or the Frobenius norm. Because Ω is obtained from a Gaussian distribution, the matrices 1 Ω and 2 Ω also exhibits the features of a Gaussian distribution. As a consequence, the matrices U and V are not involved in the analysis and we can simply assume that A is diagonal, is typically large, and very unstable [3].
First observe that the expectation of the third term of (21) and making use of proposition 1, † † † 2 2 2 2 1 1 1

Numerical Experiments
In this section, we include numerical instances to demonstrate the computational efficiency of the accelerated AE-RSVDM with modified normalized power iteration scheme and to compare its efficiency with the AE-RSVDM with orthonormalization found in the manuscript of Martinsson [14]. Both algorithms are executed in MATLAB R2018a software. All computations were done on a machine with the following specifications: A desktop computer with a 4-core (8 logical processors) Intel Core i7-7700 CPU (3.60GHz), and 32GB DDR4 RAM. The relative reconstruction error is calculated using the Frobenius norm as where k is the target rank and M k is the approximated matrix. The spectral norm was also used.

Comparison of Execution Speed
We compare the run time and the reconstruction errors of the accelerated variant to that of the AE-RSVDM with orthonormalization. A synthetic 'thin' matrix whose singular values do not decay easily of size 10000 8500 × was used to test the efficiency of the two algorithms. We controlled the efficiency of the algorithms with varying q . A target rank 3000 = l was used. The results are shown in Table 2.     From Table 2 and Table 3, the accelerated AE-RSVDM with modified normalized power iteration achieves speed-ups over the AE-RSVDM (with orthonormalization) over different , q while attaining nearly the same reconstruction errors on both matrices used in the numerical computations. Refer to Figure 4 and Figure 5.

Image Compression
SVD is a common method for image compression. An image can be regenerated from a subset of the principal singular values and singular vectors. We use Benvenuto Tisi's well-known painting "Garofalo Annunciation" which was obtained from the Google Art Project. This image has a high resolution of 30000 22882.

×
We seek to compute the low-rank approximation of this high-resolution image via the two randomized algorithms to ascertain the time each algorithm takes.
To begin, we reduced the resolution to 5120 3072 × pixels in order for the computer to bear the computational load. We seek to show the reconstruction of the data image with different . k p = + l Figure 6 depicts the original input of the "Garofalo Annunciation" image.

Conclusion
Reduced dimensionality and the associated theory of low-rank matrix approximations are crucial algorithmic techniques in many fields of study. Big data are, however, increasingly difficult to compute in classical matrix algorithms. The randomized SVD is undoubtedly the greatest influential and ubiquitous method for randomized algorithms. This work has examined both the runtime and the relative reconstruction errors of the accelerated AE-RSVDM with modified normalized power iteration against the AE-RSVDM (with orthonormalization) as formulated by [15]. Inspired by randomization, the computation of the approximate low-rank SVD of a big data matrix from a small compressed matrix is based on an efficient two-pass algorithm. The results show that the accelerated AE-RSVDM variant achieves a faster runtime than the AE-RSVDM (with orthonormalization) but does not achieve any significant reconstruction errors. Future research focuses on improving on the relative approximate error as well as implementing the accelerated variant in facial recognition.