System and method for rapidly processing singular value decomposition of ultra-large scale low rank matrix

Abstract

A computer-implemented method for generating an approximate image of a matrix A includes constructing, in a processing unit, a decomposition of C, U, and R matrices of the matrix A, where the C matrix includes k-column vectors and the R matrix includes k-row vectors retrieved from a memory holding the matrix A; and performing singular value decomposition on C, U and R matrices in a processing unit to generate an approximate matrix # imgabs0 # such that # imgabs1 # for matrix A

Description

Technical Field

[0001] The present disclosure relates to a system and method for quickly processing singular value decomposition (SVD) of very large low-rank matrices. Background Art

[0002] Matrix factorization methods play an important role in computer applications such as computational mechanics, complex networks, machine learning, data analysis, image processing, and signal processing.

[0003] These computer applications that need to process large-scale low-rank matrices will always be performed in computer systems with limited resources. Under the growing demand for processing, analyzing and interpreting such high-dimensional data, methods for approximating matrices with low-rank counterparts have been used for practical applications. One of the most widely used methods for low-rank matrix analysis is singular value decomposition (SVD), which is a basic tool in linear algebra, statistics and machine learning. SVD used for dimensionality reduction techniques (e.g., principal component analysis (PCA)) maps data from high dimensions to low-dimensional space. Real matrix The SVD of A = UΣV T ,in, and are orthogonal matrices representing the left and right singular vectors, respectively, and Σ is an m×n diagonal matrix of singular values ​​in descending order.

[0004] Although SVD provides the best low-rank approximation, it requires high memory cost and large computation time, which is not possible for large matrices on ordinary desktop computers. The optimal SVD of a dense low-rank (k) matrix A of size m×n requires O(mn 2 +m 2 n+n 3 ) time. If optimality is not required, computationally efficient low-rank approximations based on randomized projection methods can be used to approximate the principal subspace of the data matrix. For example, the Lanczos method has a time complexity of O(mnk 2 ). However, to store the entire data matrix, the rank must be known in advance, and the process is still computationally expensive, especially for data projection.

[0005] Classical or randomized SVD methods for large data require reading the entire data at least once, storing large temporary matrices during computation, and performing large matrix multiplications during their projection steps.

[0006] Prior art matrix factorization methods are slow and limited by system specifications, making them non-scalable.

[0007] US Patent No. 10229092 discloses a method for p -norm space using low-rank matrix factorization, where p<2 (e.g., 1≤p<2). This method provides l 1 -Minimum absolute deviation of PCA technique (l 1 -norm) space. This method decomposes the residual matrix in the subspace of the observed data matrix into p -norm is minimized, thereby minimizing the l of the residual matrix 1 -norm, where p = 1. The alternating direction method of multipliers (ADMM) is applied to solve the subspace decomposition of the low-rank matrix factorization relative to the observed data matrix. The iteration of ADMM can include solving the subspace decomposition and calculating l 1 -norm. However, this approach is very time-consuming and does not scale well for large matrices.

[0008] U.S. Patent No. 9760537 discloses a computer-implemented method for obtaining a CUR decomposition. The method includes: constructing a matrix C based on a matrix A by a computer processor; constructing a matrix R based on the matrix A and the matrix C; constructing a matrix U based on the matrix A, C and R; the matrices C, U and R provide a CUR decomposition of the matrix A; the construction of the matrices C, U and R provides at least one of an input sparse temporal CUR and a deterministic CUR. This method is time-consuming and not well suited for large matrices. Summary of the invention

[0009] According to a first aspect of the present disclosure, there is provided a computer-implemented method for generating an approximate image of a matrix A. The computer-implemented method comprises:

[0010] Constructing, in a processing unit, a decomposition of C, U and R matrices of a matrix A, wherein the C matrix includes k column vectors and the R matrix includes k row vectors retrieved from a memory storing the matrix A;

[0011] In the processing unit, singular value decomposition is performed on the C, U and R matrices to generate an approximation matrix of the matrix A. Make

[0012] According to a first aspect, the C and R matrices comprise rows and columns of the matrix A given by indices I and J randomly generated from a uniform distribution.

[0013] According to a first aspect, the matrix U is The optimal solution of C is retrieved from the memory storing the matrix A using a deterministic sampling method. u , A u and R u rows and columns.

[0014] According to a first aspect, the deterministic sampling method is a discrete empirical interpolation method (DEIM).

[0015] According to the first aspect, the step of performing singular value decomposition on the C, U and R matrices comprises: generating a matrix [U c ,S c ,V c ] and [U r ,S r ,V r ], so that and

[0016] According to the first aspect, the step of performing singular value decomposition on the C, U and R matrices comprises: the step of generating a matrix Z such that

[0017] According to the first aspect, the approximation matrix By decomposing the matrix Z, we can generate

[0018] According to the first aspect, U k =U c U z ,∑ k =S z , and V k =V z V r .

[0019] According to a first aspect, k is greater than the rank of the matrix A.

[0020] According to the first aspect, k is an arbitrary integer.

[0021] According to a first aspect, the error tolerance ε is stored in a memory of the processing unit.

[0022] According to a first aspect, a distance function is generated for the singular value decomposition of the C and R matrices: and in, and

[0023] According to the first aspect, if one of the distance functions returns a value greater than the error tolerance ∈, k, the value of k will be increased to construct new C and R matrix decompositions with increased rows and columns in the processing unit until the distance function returns a value less than the error tolerance ∈.

[0024] According to the first aspect, the processing unit is used to calculate the approximate matrix using some elements of the matrix A.

[0025] According to a first aspect, the processing unit comprises a cache memory or a random access memory of a size smaller than that required to hold the matrix A.

[0026] According to the first aspect, the processing unit is configured to load only the partial elements of the matrix A into the cache memory or the random access memory for processing.

[0027] According to a first aspect, the matrix A comprises a plurality of digital images.

[0028] According to a first aspect, each of said digital images is represented as a column vector in a matrix A.

[0029] According to the first aspect, the matrix A includes time series data, wherein each of the column vectors represents snapshot data of a time interval. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Embodiments of the present disclosure will be described below by way of examples with reference to the accompanying drawings.

[0031] Figure 1 Schematic diagram of a computer server according to an embodiment of the present disclosure, wherein the computer server is used to implement a system for analyzing compressed video.

[0032] Figure 2 A block diagram of a computer-implemented method according to an embodiment of the present disclosure.

[0033] Figure 3a Included is a graph showing a comparison between fixed-rank CUR-SVD, SVD, SVD one-pass, and SVD randQB methods for computation time and storage usage when computing the first k singular values ​​and vectors for different matrix sizes and different values ​​of k.

[0034] Figure 3b Included is a graph showing a comparison between fixed-rank CUR-SVD, SVD, SVD one-pass, and SVD randQB methods for computation time and storage usage when computing the first k singular values ​​and vectors for different matrix sizes and different values ​​of k.

[0035] Figure 4aIncluded are plots showing a comparison of the fixed-rank CUR-SVD, SVD, SVD one-pass, and SVD randQB methods for approximation accuracy (relative error) when computing the first k singular values ​​and vectors for different matrix sizes and different values ​​of k.

[0036] Figure 4b Included are graphs showing the accuracy of the composite singular values ​​for three matrices with different singular value distributions (slow decay, fast decay, and sigmoid decay).

[0037] Figure 5a Included are graphs showing the impact of the parameters (block size, error tolerance, and matrix size) of the fixed precision CUR-SVD method of the present disclosure on the estimated rank, running time, relative error, and storage usage of the performance of SVD-based methods.

[0038] Figure 5b Included are graphs showing the impact of the parameters (block size, error tolerance, and matrix size) of the fixed precision CUR-SVD method of the present disclosure on the estimated rank, running time, relative error, and storage usage of the performance of SVD-based methods.

[0039] Figure 6a Included is a diagram showing the extended Yale B face dataset, CUR-SVD performance evaluation for PCA on the dataset, setting k to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0040] Figure 6b Included is a diagram showing the extended Yale B face dataset, CUR-SVD performance evaluation for PCA on the dataset, setting k to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0041] Figure 6c Included is a diagram showing the extended Yale B face dataset, CUR-SVD performance evaluation for PCA on the dataset, setting k to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0042] Figure 6d Included is a diagram showing the extended Yale B face dataset, CUR-SVD performance evaluation for PCA on the dataset, setting k to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0043] Figure 6e Included are graphs showing the CUR-SVD performance evaluation of PCA on the extended Yale B face dataset of Figures 6A, 6B, 6C, and 6D, with k set to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0044] Figure 6f Included are graphs showing the CUR-SVD performance evaluation of PCA on the extended Yale B face dataset of Figures 6A, 6B, 6C, and 6D, with k set to 300, containing 38 individuals, each with 64 images of size (168×192), which can be reshaped into a 32256×2432 matrix.

[0045] Figure 7a Included are graphs showing the performance of the proposed fixed-k (i.e., CUR-SVD) method under different constructions of the U matrix based on (i.e., CUR-SVD(I), CUR-SVD(O), and CUR-SVD(D)). For CUR-SVD(I) and CURSVD(O).

[0046] Figure 7b Included are graphs showing the performance of the proposed fixed-k (i.e., CUR-SVD) method under different constructions of the U matrix based on (i.e., CUR-SVD(I), CUR-SVD(O), and CUR-SVD(D)). For CUR-SVD(I) and CURSVD(O).

[0047] Figure 8 Included is a graph showing a comparison of singular value approximations between the proposed fixed-rank CUR-SVD and the ground truth, the matrix size is 1,000,000×1,000,000, where rank=100.

[0048] Fig. 9 Included are graphs showing verification of the accuracy of singular vectors computed by the fixed-rank-k algorithm compared to the SVD method under different singular value distributions. DETAILED DESCRIPTION

[0049] The present disclosure is designed to solve two problems in current computing systems that are physically limited in working memory and non-volatile memory: (1) representing large data matrices using a small number of rows and columns, thereby solving the memory problem; and (2) computing a satisfactory approximation of the SVD of the original matrix using these partial data sets.

[0050] The method of the present disclosure can be used in Figure 1The present invention can be implemented in the computing system shown or any embedded system to achieve significant improvements in speed and performance for processing applications that require manipulation of large-scale data matrices, such as AI matching engines, convolutional neural networks, etc.

[0051] In a preferred embodiment of the present disclosure, a computer-implemented method for generating an approximate image of a matrix A is provided, such as Figure 2 The computer-implemented method includes: constructing a decomposition of C, U and R matrices of matrix A in a processing unit, wherein the C matrix includes k column vectors and the R matrix includes k row vectors retrieved from a memory storing matrix A; and performing singular value decomposition on the C, U and R matrices in the processing unit to generate an approximate matrix of matrix A. Make In one embodiment, the processing unit includes a cache memory or a memory associated with a random access memory having a limited capacity. Typically, the capacity of the cache and RAM is less than the capacity required to hold the matrix A. During the execution of the method of the present disclosure, the matrix A will be stored in a hard disk drive. The processing unit is used to use the partial elements of the matrix A to calculate the approximate matrix

[0052] The C and R matrices consist of the rows and columns of matrix A given by indices I and J randomly generated from a uniform distribution. The matrix U is The optimal solution of , where C is retrieved from the memory holding the matrix A using a deterministic sampling method u , A u and R u The processing unit is used to load only part of the elements of the matrix A into the cache memory or the random access memory to process the method of the present disclosure.

[0053] The deterministic sampling method is the Discrete Empirical Interpolation Method (DEIM). The steps of performing singular value decomposition on the C and R matrices include: generating the matrix [U c ,S c ,V c ] and [U r ,S r ,V r ], so that and The steps of performing singular value decomposition on the C and R matrices include: generating a matrix Z such that Approximate Matrix By decomposing the matrix Z, we can generate Among them, U k =U c U z ,∑k =S z , and V k =V z V r .

[0054] The value of k is greater than the rank of the matrix A. In one embodiment, k is an arbitrary integer.

[0055] In one embodiment, when the rank of matrix A is unknown, multiple singular value decompositions of the C matrix and the R matrix are generated with increasing block sizes of the C matrix and the R matrix until the difference between 2 consecutive iterations is less than an error tolerance.

[0056] The method comprises the step of assigning a value to an error tolerance ε stored in a memory of the processing unit. Generating a distance function for the singular value decomposition of the C and R matrices and in and When one of the distance functions returns a value greater than the error tolerance ∈, the value of k will be increased to construct new C and R matrix decompositions with increased rows and columns in the processing unit until the distance function returns a value less than the error tolerance ∈.

[0057] The algorithm of an embodiment of the present disclosure has many applications, such as image processing and artificial intelligence. Generally, for image processing, matrix A includes multiple digital images, and each digital image is represented as a column vector in matrix A. In another application, matrix A includes time series data, wherein each column vector represents snapshot data of a time interval.

[0058] In order to provide solutions to these two problems, the present disclosure provides a novel and efficient method, namely CUR-SVD, to calculate the first k singular vectors and values ​​of a large-scale matrix by exploiting the low-rank structure of the large-scale matrix. Most large-scale data matrices in practical applications have a low-rank structure or can be well approximated by a low-rank matrix, allowing the large-scale data matrix to be represented by only a small part of the rows and columns. If the representation has a sufficiently high accuracy, efficient calculation of SVD can be performed for previously impossible applications.

[0059] Representing the data in a low-dimensional space with the actual data points as basis vectors based on interpretable and computationally efficient methods of selecting rows and columns from the data matrix to preserve sparsity or non-negativity, as in CUR matrix factorization, which represents the data matrix as the product of three small matrices. CUR is a low-rank approximation A≈CUR, where the matrix and is formed by extracting c columns and r rows of the matrix A, respectively, where c,r ≥ rank(A). The U matrix is ​​constructed so that the approximation error is very small. Randomized methods sample columns and rows based on the CUR decomposition from uniform or non-uniform distributions, while deterministic methods treat the approximation as an optimization problem and try to find the best decomposition that minimizes the approximation error. Low-rank approximation methods fall into two categories: (1) fixed rank, where the rank of the matrix A is known, and (2) fixed precision, where the matrix A is approximated by its approximation A. app The error between them is as small as possible to estimate the rank, ||AA app || F <∈, where is the precision tolerance given by ∈.

[0060] The system or method of the present disclosure is suitable for processing large-scale matrices that are too large to be stored in random access memory (RAM) or in a processing unit cache for calculating their SVD. Two algorithms (fixed-rank CUR-SVD and fixed-precision CUR-SVD) performed by a computing system or processing unit are proposed in the present disclosure to calculate the first k singular vectors and their corresponding singular values ​​of large-scale, low-rank data for cases where the rank is known or unknown, respectively.

[0061] In order to verify the accuracy and efficiency of the system or method of the present invention in practical situations, numerical experiments are carried out on desktop computers with data sets of different sizes and low-rank structures to demonstrate the effectiveness of CUR-SVD. The system or method of the present invention can provide approximate accuracy comparable to that of existing SVD algorithms for large-scale low-rank matrices, but is more efficient in terms of computational time and memory. For fixed-rank and fixed-precision algorithms, the first 200 singular vectors and their corresponding singular values ​​of a 1,000,000×1,000,000 matrix are calculated in approximately 3 minutes and 15 minutes, respectively, with a small approximation error and less than 3GB of storage usage. This makes CUR-SVD feasible on most current computing devices. However, traditional methods require 7.27TB of storage, which is impractical.

[0062] In this exemplary embodiment, the system and method of the present disclosure are implemented by a computer or processor. The computer can be implemented by any computing architecture, including: a portable computer, a tablet computer, a stand-alone personal computer (PC), a smart device, an Internet of Things (IoT) device, an edge computing device, a client / server architecture, a "dumb" terminal / host architecture, a cloud computing-based architecture, or any other appropriate architecture. The computing device can be appropriately programmed to implement the present disclosure.

[0063] like Figure 1, a schematic diagram of a computer system or computer server 100 is shown. The computer system or computer server 100 is used to implement an example embodiment of a system for analyzing compressed video. In this embodiment, the system includes the server 100. The server 100 includes appropriate components required to receive, store, and execute appropriate computer instructions. These components may include: a processing unit 102, which includes a central processing unit (CPU), a math co-processing unit (Math Processor), a graphics processing unit (GPU), or a tensor processing unit (TPU) for tensor or multi-dimensional array calculation or manipulation operations; a read-only memory (ROM) 104; a random access memory (RAM) 106; an input / output (I / O) device (e.g., a disk drive 108); an input device 110 (e.g., an Ethernet port, a USB port, etc.); a display 112, such as a liquid crystal display, a light-emitting display, or any other suitable display; and a communication link 114. The computing device 100 may include instructions that may be stored in the ROM 104, the RAM 106, or the disk drive 108 and executed by the processing unit 102. A plurality of communication links 114 may be provided that may be variously connected to one or more computing devices, such as servers, personal computers, terminals, wireless or handheld computing devices, IoT devices, smart devices, and edge computing devices. At least one of the plurality of communication links may be connected to an external computing network via a telephone line or other type of communication link. Without loss of generality, the processing unit 102, ROM 104, RAM 106, and disk drive 108 have limited physical limitations. For example, the processing 102 may have a 30MB cache memory and may support up to 192GB of RAM.

[0064] The server 100 may include a memory device such as a disk drive 108, which may include a solid-state drive, a hard disk drive, an optical drive, a tape drive, or a remote or cloud-based storage device. The server 100 may use a single disk drive or multiple disk drives, or a remote storage service 120. The server 100 may also have a suitable operating system 116 resident on the disk drive or in ROM of the server 100.

[0065] The computer or computing device may also provide the necessary computing power to operate or interface with a machine learning network such as a neural network to provide various functions and outputs. The neural network may be implemented locally, or may be accessed or partially accessed through a server or cloud-based service. The machine learning network may also be untrained, partially trained, or fully trained, and / or may also be retrained, modified, or updated over time.

[0066] The CUR-SVD algorithm executed in the processing unit of the computing system of the present invention includes two stages: constructing the CUR decomposition, and performing SVD on the C and R matrices of the small decomposition. Compared with existing methods, the SVD calculation of large-scale matrices is very efficient in terms of computing time and memory. Large-scale matrices are sometimes too large to be stored in computer memory. Therefore, the CUR-SVD algorithm of the present invention can be used in other valuable applications, such as the efficient calculation of the inverse matrix of extremely large-scale low-rank matrices.

[0067] Figure 2 The methods shown can be implemented on a computer system in hardware, software, or a combination thereof. Figure 2 An overview of the fixed-rank CUR-SVD method for computing the SVD of large-scale low-rank matrices is shown. CUR-SVD consists of two stages: constructing a CUR decomposition, and performing SVD on small decomposition matrices C and R. Compared with existing methods, CUR-SVD is highly computationally and memory efficient for SVD computation of large matrices.

[0068] Fixed (known) rank CUR-SVD method

[0069] For large matrices with low rank structure rank(A) = ρ < < min{m,n} Randomly select c columns of matrix A to form where c>ρ. Then, with high probability, the C matrix will capture the column span of the matrix A. Similarly, r rows of the matrix A are randomly selected to form where r>ρ. Then, with high probability, the R matrix captures the row span of the matrix A. The large-scale matrix A represented by the CUR decomposition is A≈CUR, where the U matrix is ​​the link between the columns and row spans of the matrix A, which makes the approximation error very small. Therefore, the left and right singular vectors and values ​​of the matrix A can be obtained from the small matrices C, U, and R without having to load the large-scale matrix A into memory, as Figure 2 The first stage is shown.

[0070] Below is the pseudo code of Algorithm 1 for generating an approximate decomposition matrix of matrix A, where

[0071] Algorithm 1 gives the pseudocode of the fixed-rank CUR-SVD method, which uses the small matrix of the CUR decomposition of the matrix A to calculate the first k singular vectors and values ​​of the large-scale matrix A. The first stage constructs the CUR decomposition of matrix A using only k (for simplicity, assume c=r=k) rows and columns of matrix A, and the second stage calculates the k-SVD of matrix A using the CUR decomposed matrix.

[0072] The input is a matrix A with k>>ρ, i.e. the number of rows and columns to be selected to form the C and R matrices, which is greater than the rank of the matrix A. The row and column indices I and J are randomly generated from a uniform distribution, and the R and C matrices are formed from the rows and columns of the matrix A given by I and J.

[0073] The role of the U matrix: The U matrix of CUR decomposition needs to be constructed so that the approximation error or ‖A-CUR‖ 2 As small as possible. The simplest method, the pseudo-inverse of the cross matrix between the C matrix and the R matrix, gives a high approximation error for high-rank or noisy data. The optimal solution is Among them, the U matrix is ​​formed by weighting all elements of the matrix A, thereby giving a small approximation error, but the computational overhead is high, especially for large-scale matrices. In order to balance the approximation accuracy and computational complexity, it is proposed to form the U matrix by selecting more rows from the C matrix and more columns from the R matrix, rather than the intersection between the c columns and the r rows, and then use the optimal form to calculate the U matrix. More rows and columns are obtained from the C matrix and the R matrix through a deterministic sampling method called discrete empirical interpolation method (DEIM). This method uses the information of singular vectors to capture the most informative rows and columns. DEIM is a point selection algorithm that was originally used to model order reduction in nonlinear dynamic systems. When forming the C matrix and the R matrix in the CUR decomposition, the point selection algorithm is used as a fully deterministic algorithm to improve the error bound.

[0074] Algorithm 3 shows the process of calculating the selected row (and similarly column) index I by the DEIM algorithm. The input is a full rank matrix U∈R^(m×k) containing the first k left singular vectors. The algorithm processes one singular vector at a time, and each vector provides an index, as shown in steps 2 to 9. Therefore, the maximum number of indexes generated is the number of available singular vectors, k.

[0075]

[0076] After constructing the C matrix and the R matrix, their SVD is calculated to give and As shown in lines 5 and 6 of Algorithm 1 and Figure 1Now, the left singular vectors (U c ) and the right singular vectors of the R matrix (V r ) is available, so DEIM is applied to obtain k indices for the rows of the C matrix and k indices for the columns of the R matrix. To avoid duplication, any duplicate indices are removed. Then C is constructed u , R u and the intersection between them (Au), as shown in lines 10 to 12 of Algorithm 1. The optimal form of SVD is used to compute the U matrix using the intersection between (≤2k) rows and (≤2k) columns, where (c=r=k).

[0077] The CUR-SVD of the U matrix based on intersection, optimal solution and DEIM method are denoted as CUR-SVD(I), CUR-SVD(O) and CUR-SVD(D), respectively. Figure 6 shows the impact of the U matrix on the running time and approximation error. All methods provide accurate results when k is greater than 20, where k is the rank of matrix A. The impact of the U matrix on the computation time and approximation error. It can be seen that all methods provide accurate approximations when k is greater than 20, where k is the rank of matrix A. For the case where k is less than the rank of matrix A, the approximation error is very high and unstable for CUR-SVD(I), but CUR-SVD(O) and CUR-SVD(D) methods have almost the same relative error values, which are comparable to SVD.

[0078] For the case where k is less than the rank of the matrix A, the approximation error is high and unstable for CUR-SVD(I), but similar to the SVD for CUR-SVD(O) and CUR-SVD(D). Figure 6b Their running times are shown. CUR-SVD(D) requires slightly higher running time than CUR-SVD(I), but both are faster than CUR-SVD(O). Subsequently, the U matrix is ​​constructed using the CUR-SVD(D) method.

[0079] Consider only the C matrix, U matrix and R matrix, and a total of (mc+nr+rc) elements, instead of all elements of the entire matrix A, and calculate the complexity separately and The complexity is replaced by the SVD of the C matrix, Z matrix and R matrix. Computation of SVD and complexity of matrix A.

[0080] In having All experiments were implemented using MATLAB R 2022a on a desktop computer with a Core TM17-10700 processor @ 2.90 GHz and 16 GB RAM. Figure 3aThe run time performance of the different methods for different matrix sizes and different values ​​of k is shown. The original SVD method provided in MATLAB is the most computationally expensive, and CUR-SVD outperforms the other methods. When computing the first k singular values ​​and vectors, CUR-SVD outperforms all other methods in terms of the memory required to compute k-SVD ( Figure 3b ). For a matrix of size 10,000×10,000, CUR-SVD requires only about 10MB of memory to compute the first 120 singular vectors and values.

[0081] The relative errors of different methods are Figure 3a As shown in Figure 2, CUR-SVD has an approximate accuracy comparable to that of SVD and SVD_rand QB, and is superior to the SVD single-pass method.

[0082] For the three different patterns, CUR-SVD gives more accurate singular values ​​than randomization and one-pass SVD, especially for the S-shaped pattern ( Figure 3b ). Figure 8 It shows that the orthogonality error for different singular value modes of CUR-SVD is very small, about 1e -14 , although it is larger than the error of the SVD method. The SVD method requires the entire matrix A, rather than just k rows and columns from the matrix A for CUR-SVD.

[0083] The timing, approximation error, and storage usage of CUR-SVD for computing the top-k-SVD for k = 100 and 200 from large matrices are given in Table 1. For the first 200 singular vectors and a 1,000,000 × 1,000,000 matrix, it takes about 82 seconds and only 3.052 GB of storage, achieving about 1.2e -12 This is superior to the prior art methods in terms of data memory and computation time requirements.

[0084] Table 1

[0085]

[0086] Figure 7 shows that for different low-rank singular value patterns, the singular values ​​from CUR-SVD are almost the same as the ground truth data. The original matrix occupies about 7.27TB of storage, which is impractical. However, CUR-SVD efficiently computes the SVD result of the matrix using 3.052GB of storage. In order to evaluate the PCACUR-SVD method relative to SVD, a face recognition benchmark dataset is used. When calculating the first 300 singular vectors and their corresponding singular values ​​of the Yale Extended B dataset, CUR-SVD is 12 times faster than SVD with high approximation accuracy.

[0087] By projecting one of the test images onto the UUTs of SVD and CUR - SVD while increasing the number of singular vectors of U from 20 to 300 in steps of 20, the Figure 6c and 6d reconstructed surfaces shown in Figure 6e were obtained, which demonstrated the accuracy of the CUR - SVD method.

[0088] Projecting 128 test images of two people by PCA using SVD and CUR - SVD methods showed their clear separation and classification ( Figure 6f , for the fifth and sixth PCA modes) and the accuracy of CUR - SVD.

[0089] Fixed - precision CUR - SVD method

[0090] The fixed - precision CUR - SVD method estimates the exact rank k and simultaneously computes the SVD without having to store the entire data set and the approximation matrix, making it efficient in terms of computational time and memory.

[0091] For a large matrix with a low - rank structure rank(A)=ρ << min{m,n} and C and R matrices from CUR decomposition If rank(C)=rank(A)=rank(R), then with high probability rank(U)=rank(A) and A≈CUR. Assuming that the number of columns chosen to form the C matrix is c >> ρ, then

[0092] rank(C)≈rank(A)=ρ (2.2.1)

[0093] By applying thin SVD on the C matrix, with singular values C, s 1 , …, s ρ , …, s c . For k < ρ < c, then

[0094]

[0095] If the rank k is set to ρ < k < c, then

[0096]

[0097] Similarly, for the R matrix, it is obtained by constructing an R matrix with r > ρ rows and truncating the R matrix with rank k > ρ:

[0098]

[0099] For the fixed - precision algorithm, the goal is to make ‖A - A k‖ 2 <∈, where ∈ is the error tolerance, find k that satisfies / / AΘk / / 2<<ε or find k that makes σi(A) approximately fixed for i=k+1,···,min{m,n}. This requires computing the SVD of the matrix A, increasing the rank k, and checking its singular values ​​(i.e., σi(A)) until it is almost fixed. Computing σi(A) is impractical for matrices of large size. Fixed-rank CUR-SVD can overcome this computational and memory complexity by selecting k columns and rows and computing the k-SVD using Algorithm 1. Then, the number of selected rows and columns can be increased and the k-SVD recomputed until σi(A) is almost fixed.

[0100] Recalculating k-SVD using fixed-rank CUR-SVD is the main computational cost in the entire algorithm, and as the number of selected rows and columns increases in each iteration, the SVD of the C matrix and the R matrix must be recalculated. Therefore, an embodiment of the present disclosure proposes to use the calculated singular vectors and corresponding singular values ​​from the previous iteration to calculate the k-SVD of the current iteration. Incremental SVD (given in Algorithm 4) is used to recalculate the SVD of the C matrix and the R matrix in each iteration, thereby reducing the computational complexity.

[0101] for c>>ρ, rank(C)≈rank(A)=ρ, and wr>>ρ, rank(R)≈rank(A)=ρ, then, rank(U)≈rank(A)=ρ and A=CUR. Therefore,

[0102] σ i (A)≈σ j (C)≈σ j (R)(2.2.5)

[0103] For i = ρ+1,…,min{m,n} and j = ρ+1,…,min{c,r}.

[0104] The stopping condition is the minimum value of j for which the singular values ​​of the C matrix and the R matrix become almost fixed. For the C matrix, let in is the i-th singular value of X and R, Therefore, δ c and δ r All need to be smaller than a tolerance parameter, ∈.

[0105] Algorithm 2 below gives the fixed precision CUR-SVD.

[0106]

[0107]

[0108] The input is a large matrix Tolerance parameter ε and block size b. The block size is the increment of the number of rows and columns to be incremented in each iteration. For simplicity, in this embodiment, it is assumed that c=r.

[0109] Initialization step: Generate column and row indices J and I using a uniform sampling strategy, and read b columns and rows of matrix A corresponding to the first b indices to construct C matrix and R matrix.

[0110] Procedure: Calculate the SVD of the C matrix and the R matrix and initialize δ c and δ r If necessary, new blocks of b columns and b rows of matrix A corresponding to the second subset of indices J and I are read into memory (from disk) to construct new matrices C′ and R′, and the SVD of the C and R matrices is updated using IncSVd (Incremental SVD Algorithm 4), where δ c and δ r is updated.

[0111] In one embodiment of the present disclosure, Brand's incremental SVD main algorithm is implemented in two variations. The main algorithm is modified to process a block of columns instead of a column in each iteration, and a modified Gram-Schmidt algorithm 4 is used to re-orthogonalize after each iteration. The symbols A(:, 1:l) and A l refers to the first l columns of matrix A, and the function Takes a vector as input and returns a square matrix with the elements of the input vector as its diagonal elements. Let matrix A be an m×n dense matrix, ai refers to the i-th column of matrix A, that is, A = [a 1 |a 2 |…|a n ]. Assume A l The SVD of A is calculated as follows: l =U∑V T ,in Is A l The diagonal matrix of the l(order) singular values ​​of ,

[0112] Incremental SVD aims to update the above U, Σ, and V matrices to find the SVD of A1+b=[A1 and B], where is the new b-column addition matrix.

[0113] Algorithm 4 shows the use of the calculated A l SVD of A and newly added column B to calculate l+b The SVD is the incremental SVD procedure.

[0114]

[0115]

[0116]

[0117] If the stopping condition ((δ c <∈)&&(δ r <∈)), then the estimated rank k is k=e×b, where e is the number of times the stopping condition is invalid. and The matrices and their SVDs [Uc, Sc, Vc] and [Ur, Sr, Vr] are the final results from IncSVD in rows 15 and 16. Thus, the rank k and SVD of the C and R matrices are computed simultaneously.

[0118] As with the fixed-rank CUR-SVD algorithm, after constructing the C matrix and the R matrix, the U matrix is ​​calculated. Finally, the intermediate matrix Z and its SVD [Uz, Sz, Vz] are calculated to give, and in

[0119] Estimated rank of the fixed precision CUR-SVD algorithm for all parameters (error tolerance ∈, block size b and matrix size n) for different tolerance parameters (row-wise) Running time, estimated error p n (E) and the impact of storage usage on Figure 4a The actual rank of the data matrix is ​​150, and the distribution of singular values ​​is around β = 0.5. CUR-SVD has better performance in terms of computational time and storage usage, and has similar performance to the truncated SVD method ( Figure 4b ) has a relatively similar approximation error.

[0120] For large error tolerance, ∈ = 1e -3 , estimated rank is quite different from the actual rank k. Especially for small block size b (first row, Figure 4a ), although the running time and memory usage are very small, it still produces a large approximation error p n (E). The error tolerance ∈ is 1e -8 ( Figure 4a ), the estimated rank Almost the actual rank when computing the first 160 singular values ​​and vectors of a 20,000×20,000 matrix, especially for small block sizes (with a runtime of 18 seconds, about 5e -12 Small P n(E) and about 50MB of storage usage). A trade-off must be made between ∈ and the accuracy of the approximation and between the block size b and the time and space complexity. Based on experiments, it is recommended to set ∈ to 1e -8 to ensure a small approximation error, and set b to 20 for efficient use of time and space.

[0121] Tables 2 and 3 show the performance of the fixed precision CUR-SVD algorithm for small and large matrices, respectively. The actual rank of the generated data is 150, the block size b is 20, and the error tolerance parameter ∈ is set to 1e -8 Under different values ​​of β, the estimated rank is 180. Therefore, CUR-SVD accurately estimates the columns and rows of , guaranteeing high approximation accuracy. For a 1,000,000×1,000,000 matrix, CUR-SVD requires about 2.71 GB of storage to compute the first 180 singular vectors and singular values ​​in about 15 minutes. For different matrix sizes and singular value patterns, the estimated error p i (E) is about 1e -13 .

[0122] Table 2

[0123]

[0124] Table 3

[0125]

[0126]

[0127] Embodiments of the present disclosure propose a novel method for calculating the highest k-SVD of very large-scale, low-rank matrices, which cannot be calculated by existing SVD algorithms. CUR-SVD is based on randomly observing a portion of the data using CUR decomposition. For different low-rank approximation problems, two algorithms, fixed k-rank and fixed precision CUR-SVD, are proposed, and the approximate accuracy comparable to the optimal simplified SVD algorithm is achieved. These methods do not analyze the entire matrix, but work through CUR decomposition and are more efficient than traditional SVD algorithms. When only a small portion of the matrix elements is observed, the method can perform SVD on extremely large matrices that cannot be loaded into computer memory, and perform SVD in much less time than the time required by standard methods.

[0128] method

[0129] The performance of CUR-SVD was evaluated by comparing with widely used methods for computing SVD. The MATLAB built-in command for truncated SVD, which is based on the Krylov subspace iteration method, was used and designated as SVD in this method. For the randomized SVD method, two widely used methods (SVD randQB and SVD one-pass) were used for performance comparison, which are the main randomized method and one-pass randomized SVD, respectively. The relative errors and running times for all experiments are reported as the average of ten trials. The relative error in the Frobenius norm was used as the accuracy measure of the approximation, which is defined as For low-rank matrices, they are generated directly based on SVD. and are two randomly generated real unitary matrices, and is a non-negative real diagonal matrix with artificially set diagonal values. Then, the low-rank matrix It will be expressed as A = V ∑ W T Since existing SVD algorithms cannot handle such huge matrices, the results are compared with the ground truth of experiments conducted on large-scale matrices using generated data.

[0130] Fixed-rank CUR-SVD

[0131] To evaluate the performance of CUR-SVD in processing U matrices (as shown in Figure 6), a random matrix A of 5,000×5,000 with rank 20 was generated, and k-SVD was calculated using SVD with different constructions of the U matrix of CUR and fixed-rank CUR-SVD. CUR-SVD with intersection-based U matrix, the optimal solution, and the DEIM-based method are denoted as CUR-SVD(I), CUR-SVD(O), and CUR-SVD(D), respectively. The value of k was varied from 1 to 25, and the relative error and running time for each value were calculated.

[0132] In order to evaluate the computational and spatial efficiency of the fixed-rank CUR-SVD method, randomized SVD and one-way SVD based on QB decomposition are used in one embodiment of the present disclosure. A set of matrices of different sizes with a rank of 100 are generated, and the singular values ​​are selected as ∑ k =diag(1 -β ,2 -β ,…,k -β ), as shown in Figure 3. For both randomized SVD and single-pass SVD, the oversampling parameter was set to 10. For each matrix size, the algorithm was run to compute the first k singular values ​​and singular vectors, where k varied from 20 to 120 in steps of 20 (for CUR-SVD, k represents the number of columns and rows selected to form the C and R matrices).

[0133] In order to check the accuracy of the obtained singular values ​​and vectors by CUR-SVD, three types of matrices are generated to represent the decay of different distribution shapes (slow, fast and S-shaped) of the singular values. The singular vectors U and V are random matrices with orthogonal columns, and the diagonal elements of Σ are σ for slow, fast and S-shaped decays, respectively. j =e pj / 2 , σ j =e -j / 7 and σ j =0.0001+(1+e j-30 ) -1 . For each mode of singular values, a matrix is ​​generated And calculate the first 200 singular vectors and singular vectors of A using different methods.

[0134] In one embodiment of the present disclosure, the orthogonality of the obtained singular vectors can be evaluated by the proposed method: three matrices are generated in the same way as above. For each matrix, SVD and the proposed method are used to calculate k-SVD, where the rank is varied from 10 to 200 with a step size of 10. For each k value, the orthogonality of the calculated left and right singular vectors is calculated as,

[0135]

[0136] Among them, ||x|| ∞ =max|x i |. refer to Fig. 9 , which shows a graph verifying the accuracy of singular vectors computed by the fixed-rank algorithm compared to the SVD method under different singular value distributions. Fig. 9 Included is a graph showing the accuracy verification of the singular vectors computed by the fixed rank k algorithm compared to the SVD method for different singular value distributions.

[0137] In order to estimate CUR-SVD with large-scale matrices, the classical power method is used to estimate the spectral error norm between the actual data matrix and the approximate matrix.

[0138] For large-scale data (Table 1), a set of matrices with sizes ranging from 100,000 × 100,000 to 1,000,000 × 1,000,000 were generated with rank 150 and slowly decaying singular value distributions. For each matrix, the rank k (the number of rows and columns chosen to construct the CUR) was set to 100 or 200. CUR-SVD was applied to large-scale matrices of 1,000,000 × 1,000,000 with four different singular value distributions: linear distribution with actual rank 100, multinomial distribution, cosine distribution, and trapezoidal distribution ( Figure 8 ).

[0139] In order to evaluate the PCACUR-SVD method against SVD, a face recognition benchmark dataset was used. The YaleExtended B dataset contains images of 38 individuals, each with 64 images under different lighting conditions, such as Figure 6a and 6b Each image has 192×168 pixels and is reshaped to 32,256 elements. Large column vector. In this example, images of 36 people are selected for training, and SVD and CUR-SVD methods are applied to the face data matrix. The k parameter of CUR-SVD is set to 300.

[0140] Fixed precision CUR-SVD

[0141] Matrix A k According to the three matrices U of its k-SVD k ,Σ k and V k To specify. The matrix U k and V k is randomly generated from a normal distribution, and then the orthogonal matrices of size m×k and k×n are obtained using QR decomposition. k The singular values ​​of are β∈{0.5, 1.0, 1.5, 2.0}. The small and large matrices are used to evaluate the performance of CUR-SVD (Tables 2 and 3 above, respectively). The small matrix with n The number of n is varied from 10,000 to 100,000 in steps of 10,000. The number of n is varied from 200,000 to 1,000,000 in steps of 100,000.

[0142] For each configuration, the computation time, storage usage, estimated rank k, and estimated error p are given. i (E). Given the matrix Provides a random column vector with independent and identically distributed elements Each element is distributed as a Gaussian random variable with zero mean and unit variance. Let The column vector is Then the spectral norm of the matrix A estimated by the power iteration method in the i-th step is given,

[0143]

[0144] For large values ​​of i, pi(A) estimates ||A|| with high probability. 2 In this exemplary embodiment, i=30 is used to ensure accurate estimation.

[0145] Assume U k ,∑k and V k Based on The actual matrix decomposed using k-SVD. Let and is the output matrix of the proposed method. Then, the reconstruction matrix is For large matrices, due to the calculation of A and The multiplication complexity and high storage cost of the original data matrix A and the approximate matrix The error between is calculated as However, in order to avoid computing A and , we can directly estimate the spectral norm of the error between them as follows The complexity of but

[0146]

[0147] pi(E) estimates with high probability The value of .

[0148] To evaluate the performance of the fixed-precision CUR-SVD algorithm compared to SVD using MATLAB, matrices of size 1,000 to 25,000 were generated with a step size of 1,000 using the parameter settings given above. The time required for each matrix size, the estimated error p, and the error p are reported for different decays of the singular value distribution β = 0.5. i (E) and storage use ( Figure 5b ).

[0149] Although not required, the embodiments described with reference to the accompanying drawings may be implemented as an application programming interface (API) or a series of libraries used by developers, or may be included in another software application, such as a terminal or personal computer operating system or a portable computing device operating system. Generally, since program modules include routines, programs, objects, components, and data files that help perform specific functions, those skilled in the art will understand that the functionality of a software application may be distributed among multiple routines, objects, or components to achieve the same functionality required herein.

[0150] It should also be understood that in the case where the methods and systems of the present disclosure are implemented in whole or in part by a computing system, any suitable computing system architecture may be used. Computing systems may include tablet computers, wearable devices, smart phones, IoT devices, edge computing devices, stand-alone computers, network computers, cloud-based computing devices, and dedicated hardware devices. Where the terms "computing system" and "computing device" are used, these terms are intended to cover any suitable configuration of computer hardware capable of implementing the described functionality.

[0151] It will be appreciated by those skilled in the art that various changes and / or modifications may be made to the present disclosure as shown in the specific embodiments without departing from the spirit or scope of the present disclosure as broadly described. Therefore, the present embodiments are to be considered in all aspects as illustrative and not restrictive.

[0152] Unless otherwise indicated, any reference to prior art contained herein is not to be taken as an admission that the information is common general knowledge.

Claims

1. A computer-implemented method for generating an approximate image of a matrix A, characterized in that The computer-implemented method comprises: Constructing, in a processing unit, a decomposition of C, U and R matrices of a matrix A, wherein the C matrix includes k column vectors and the R matrix includes k row vectors retrieved from a memory storing the matrix A; In the processing unit, singular value decomposition is performed on the C, U and R matrices to generate an approximation matrix of the matrix A. Make 2. The computer-implemented method of claim 1, wherein: in, The C and R matrices consist of the rows and columns of the matrix A given by indices I and J randomly generated from a uniform distribution.

3. The computer-implemented method of claim 2, wherein: in, The U matrix is The optimal solution of C is retrieved from the memory storing the matrix A using a deterministic sampling method. u , A u and R u rows and columns.

4. The computer-implemented method of claim 3, wherein: in, The deterministic sampling method is the discrete empirical interpolation method (DEIM).

5. The computer-implemented method of claim 4, wherein: in, The step of performing singular value decomposition on the C, U and R matrices includes: generating a matrix [U c ,S c ,V c ] and [U r ,S r ,V r ], so that and 6. The computer-implemented method of claim 5, wherein: in, The step of performing singular value decomposition on the C, U and R matrices includes: generating a matrix Z such that 7. The computer-implemented method of claim 6, wherein: in, The approximate matrix By decomposing the matrix Z, we can generate 8. The computer-implemented method of claim 7, wherein: in, U k = U c U z ,∑ k = S z and V k = V z V r 。 9. The computer-implemented method of claim 8, wherein: in, k is greater than the rank of the matrix A.

10. The computer-implemented method of claim 8, wherein: in, k is an arbitrary integer.

11. The computer-implemented method of claim 10, wherein: in, The error tolerance ε is stored in a memory of the processing unit.

12. The computer-implemented method of claim 11, wherein: in, Generate a distance function for the singular value decomposition of the C and R matrices and in and 13. The computer-implemented method of claim 12, wherein: in, If one of the distance functions returns a value greater than the error tolerance ∈, k, the value of k will be increased to construct new C, U and R matrix decompositions with increased rows and columns in the processing unit until the distance function returns a value less than the error tolerance ∈, k.

14. The computer-implemented method of claim 1, wherein: in, The processing unit is used to calculate the approximate matrix using some elements of the matrix A.

15. The computer-implemented method of claim 14, wherein: in, The processing unit includes a cache memory or random access memory of a size smaller than that required to hold the matrix A.

16. The computer-implemented method of claim 15, wherein: in, The processing unit is used to load only the partial elements of the matrix A into the cache memory or the random access memory for processing.

17. The computer-implemented method of claim 1, wherein: in, The matrix A includes a plurality of digital images.

18. The computer-implemented method of claim 17, wherein: in, Each of the digital images is represented as a column vector in the matrix A.

19. The computer-implemented method of claim 1, wherein: in, The matrix A includes time series data, wherein each column vector represents snapshot data of a time interval.

Images