{'fraction_non_alphanumeric': 0.03903071752844762, 'fraction_numerical': 0.010789458168740051, 'mean_word_length': 4.534094616639478, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rapid Robust Principal Component Analysis: CUR Accelerated Inexact Low Rank Estimation

Robust principal component analysis (RPCA) is a widely used tool for dimension reduction. In this work, we propose a novel non-convex algorithm, coined Iterated Robust CUR (IRCUR), for solving RPCA problems, which dramatically improves the computational efficiency in comparison with the existing algorithms. IRCUR achieves this acceleration by employing CUR decomposition when updating the low rank component, which allows us to obtain an accurate low rank approximation via only three small submatrices. Consequently, IRCUR is able to process only the small submatrices and avoid expensive computing on the full matrix through the entire algorithm. Numerical experiments establish the computational advantage of IRCUR over the state-of-art algorithms on both synthetic and real-world datasets.


Introduction
Principal component analysis (PCA) is one of the fundamental tools for dimension reduction. A well-known drawback of the standard PCA approach, viz., singular value decomposition (SVD), is its high sensitivity to outliers. Robust principal component analysis (RPCA) is designed to overcome this shortcoming and raise the robustness of PCA to potentially corrupted data. RPCA has received a lot of attention in recent years and appears in a wide range of applications, e.g., image alignment and rectification [1,2], face recognition [3,4], sparse graph clustering [5], NMR spectroscopy signal recovery [6], and video background subtraction [7,8].
We consider the following problem setting for RPCA: we observe a sparsely corrupted data matrix D ∈ R n 1 ×n 2 , which is the sum of the underlying low rank matrix L and sparse outlier matrix S. Our goal is to recover L and S simultaneously from D = L + S. Intuitively, RPCA can be modeled as a non-convex optimization problem: where r is the rank of the underlying low rank matrix, the symbol · 0 denotes the 0-norm, and α is the sparsity level of the underlying sparse outlier matrix; for ease of notation, we use n1 = n2 =: n through the paper, but emphasize that all results can be easily extended to non-square matrices.
One can see that the RPCA model should handle outliers better than standard PCA since outliers can be separated into the sparse matrix S. However, the solution of (1) may not be unique if the low rank component is also sparse, or vice versa [9]. To ensure the uniqueness of the solution, the following assumptions are commonly made for RPCA: Assumption 2 (α-sparsity of S). S has no more than αn non-zero entries in each of its rows and columns.
Essentially, Assumption 1 ensures the low rank component is not too sparse and Assumption 2 ensures the sparse component is not locally dense. Note that some papers use random support assumption for S instead of Assumption 2; for example, Assumption 2 will be satisfied when the support of outliers is drawn from some commonly used stochastic process [10,Theorem 1.2]. With these two assumptions, the separation of L and S, viz., RPCA, becomes a well-posed problem.

Prior Art and Contribution
RPCA was raised and popularized by earlier works [9,11,12], wherein convex relaxed formulas for RPCA were proposed and studied. Unfortunately, these earlier approaches can only achieve sublinear convergence, and are therefore computationally intensive. Later, a number of non-convex approaches were studied to solve (1) directly. In particular, [13] proposes an alternating projections based non-convex algorithm and an accelerated version was studied in [14] later. Also, a gradient descent based method was proposed in [15], which was recently modified for accelerating with ill-conditioned problems [16]. All of the aforementioned non-convex methods provide linear convergence with at least O(n 2 r) flops per iteration.
In this paper, we propose a novel rapid non-convex algorithm for solving the non-convex RPCA model (1) directly. We introduce CUR decomposition [17] into the iterative RPCA framework as a substitute for the SVD, and dramatically reduce the complexity to O(nr 2 log 2 (n)) flops per iteration.

Background of CUR Decomposition
Consider a noiseless low rank data matrix. The standard dimension reduction tool, PCA, can achieve great approximation and compression for the data; however, the approximation may lose interpretability [17]. One method of resolving this difficulty is to utilize the self-expressiveness of the data, that is, the data is generally well-represented via linear combinations of the other data points rather than in some abstract basis like the singular vectors.
CUR decomposition is an efficient tool to fill the interpretability gap in dimension reduction. The classic CUR decomposition problem asks: given a matrix L ∈ R n×n with rank r, can we decompose it into terms involving only some of its columns and some of its rows? Particularly, if r columns of L which span the column space of L and r rows which span the row space of L are chosen, then we should be able to obtain L itself from these submatrices. The answer has been known to be yes for some time. Theorem 3 is essentially folklore; the reader may consult [18] for a history and proof. From the theorem statement, it is easy to see the key to the success of CUR decomposition is whether the rank of the mixing submatrix U is equal to the rank of L. In fact, with various sampling strategies, this condition can be guaranteed with high probability if we sample an appropriate number of rows and columns. For instance, Theorem 4 presents the sampling complexity for uniform sampling. The readers are directed to [20] for more information about the sampling complexity under various settings. For ease of presentation, we will focus on uniform sampling in this paper; however, similar results can be readily obtained for other sampling strategies.
Note that U is an O(r log(n)) × O(r log(n)) matrix under uniform sampling. By Theorem 3, the only computational cost is incurred by calculating the pseudo-inverse of U , which requires only O(r 3 log 2 (n)) flops. In contrast, computing the SVD requires O(n 2 r) flops. This confirms the computational efficiency of CUR decomposition with larger n and smaller r.

Proposed Algorithm
In this section, we introduce a novel rapid yet robust algorithm aiming to solve the non-convex RPCA problem (1) directly. Within the general alternating projections framework for RPCA [13,14], we propose a CUR-accelerated algorithm, dubbed Iterated Robust CUR (IRCUR). As summarized in Algorithm 1, there are two phases at the k + 1-st iteration of IRCUR: we first project D − L k to the set of sparse matrices via hard thresholding to update the estimate of S, then project D − S k+1 to the set of low rank matrices via CUR decomposition to update the estimate of L.
We shall first discuss Phase II of Algorithm 1, i.e. updating the estimate of L. In prior art, a common approach for updating L is to use the truncated SVD, which can be very costly when n is large. Inspired by Theorem 3, we instead employ CUR decomposition as an inexact low rank approximator here. More specifically, we let where the indices I, J are generated via uniform sampling, and Hr(·) denotes the best rank r approximation (i.e., truncated SVD) to the argument. Hence, the updated estimate of L, is rank r. Per Theorem 4, it costs O(r 3 log 2 (n)) flops to obtain U † k+1 . At first glance, it appears to cost O(n 2 r log(n)) flops to form L k+1 itself; however, we actually never need to form the entire L through 10: // Phase II: updating low rank component 11: 14: IRCUR, only the CUR components of L k+1 need to be saved and output. Thus, Phase II of IRCUR costs only O(r 3 log 2 (n)) flops. Note that other rank truncation methods have been studied for CUR decomposition [21,22,23], but we have found the one proposed here to be the most computationally efficient.
Next, we will discuss Phase I of IRCUR, i.e., updating the estimate of S. As shown in [6,13,14], the projection to the set of sparse matrices can be achieved via the hard thresholding operator T ζ defined as: At each iteration, with properly chosen thresholding value ζ k+1 , we can obtain a sparse S k+1 while keeping S −S k+1 ∞ under control. One strategy is to pick ζ k+1 ≥ L−L k ∞, which implies supp(S k+1 ) ⊆ supp(S) and S − S k+1 ∞ ≤ 2ζ k+1 . In practice, we found that iteratively decaying thresholding values achieves great success with proper parameter tuning. We now turn our attention to the computational complexity of Phase I. As discussed in Phase II, to form the CUR components for L k+1 , we only need the submatrices corresponding to the indices I and J . Therefore, there is no need to apply thresholding on the entire matrix D − L k , but only on the submatrices that we need to pass to Phase II. That is, computing and passing Overall, we only save and process the submatrices through the entire algorithm, thus IRCUR is also memory efficiency. Finally, at the output stage of IRCUR, a CUR decomposition of the estimated L is returned to the user, which allows for better interpretation of the low rank component. In the case of that the user is more interested in the traditional low rank expression (i.e., SVD), we also provide an efficient method for the conversion from CUR decomposition to SVD, which is summarized in Algorithm 2. This conversion involves two n × r QR decompositions and a r × r SVD, which lead its complexity to O(nr 2 ). Hence, this conversion between the two low rank representations does not increase the overall computational complexity.

Fixed vs. Resampled Indices
It is optional to resample the indices I and J in every iteration. This divides IRCUR into two variations: IRCUR-F for fixed indices and IRCUR-R for resampled indices. IRCUR-F has minimal required data access and therefore less runtime in reality, but we can get stuck with bad submatrices if we are unlucky, although the chance is very low. On the other hand, IRCUR-R uses more redundant data from different submatrices, and thus is able to correct any one-time bad-luck situations. One can expect IRCUR-R to tolerate outlier corruptions better than IRCUR-F. However, more data access may be forbidden in some application scenarios and also result a bit more computing (e.g., DI,: F in (4) will be re-computed).
Generally speaking, one should use IRCUR-F if the data access is restricted and extremely fast speed is desired; IRCUR-R should be used when data access is inexpensive and best corruption tolerance is needed. We will further study the empirical performance of these two variations in Section 4.

Parameter Tuning
There are a few parameters that need to be tuned in IRCUR. Firstly, for the hard thresholding operator to work properly, we need to take care of two parameters: the initial thresholding value ζ0, and the thresholding decay parameter γ. Intuitively, ζ0 should be a positive number that helps us to filter out the irregular values immediately at initialization. An ideal choice of ζ0 is L ∞, the maximum magnitude of the underlying low rank component, which implies supp(S0) ⊆ supp(S) while keeping the entrywise estimation error of S under control. We can easily estimate L ∞ by prior knowledge in many real-world problems. For example, in the image and video related applications, any pixels that fall out of the normal range of [0, 255] should be considered outliers immediately.
On the other hand, the decay parameter γ should be a value in (0, 1) that reflects the anticipated convergence rate. Generally speaking, harder problems (e.g., larger µ, r, α, etc.) will require larger γ for successful recovery, while easier problems can be more quickly recovered with smaller γ and still work with larger γ. Choosing larger γ will cause slower convergence, but it will also boost the robustness of IRCUR. Empirically, we recommend γ ∈ [0.5, 0.9].
The other parameters we need to tune are the numbers of rows and columns to be sampled. As stated in Theorem 4, we require |I| = cI r log(n) and |J | = cJ r log(n) for some constants cI , cJ > 0. Typically, in the presence of noise, more rows and columns need to be sampled to provide robustness (e.g., [24]). In our case, larger α increases the minimum requirement of cI and cJ . We will further study the empirical relationship between α and the sampling constants in Section 4.

Numerical Experiments
In this section, we compare the empirical performance of two versions of IRCUR, with fixed indices (IRCUR-F) and with resampled indices (IRCUR-R), against the state-of-the-art RPCA algorithms, AccAltProj [14] and GD [15], on both synthetic datasets and a real-world video background subtraction task. We use Matlab R2020a as our testing platform, and all the results are obtained on a laptop equipped with Intel i7-8750H and 32GB RAM. The codes for both AccAltProj and GD are downloaded from the authors' websites, and we manually tuned the parameters for their best performance. In particular, the actual µ, α and κ are provided to AccAltProj and GD for guiding the parameter tuning, while we simply set ζ0 = 2 L ∞ for IRCUR. To balance between robustness and convergence speed, we pick γ = 1.25 for GD, and γ = 0.65 for AccAltProj and IRCUR. All algorithms halt when e k ≤ 10 −5 is satisfied for fair comparison in all the following tests. Moreover, we provide a sample Matlab implementation for IRCUR at https://github.com/caesarcai/IRCUR.

Synthetic Datasets
Following the setup in [9,13,14,15], we form the underlying rank r matrix L = AB T ∈ R n×n via two random Gaussian matrices A, B ∈ R n×r . The underlying outlier matrix S has uniformly sampled support, and the values of its non-zero entries are i. square, we sample equal number of rows and columns in the following experiments, i.e., |I| = |J | = cI r log(n), but I may not equal to J .

Empirical Phase Transition
We investigate the recovery ability of IRCUR with different sampling constants cI and corruption rates α. Taking the problem dimension n = 1000, this experiment runs under 3 different rank settings: r = {5, 10, 20}. For each rank, we conduct 50 random tests for every given pair of (cI , α), and a recovery is considered success if the output of IRCUR satisfies The test results are summarized as Figure 1, whereas we see that the more rows/columns we sample, the more corruption/outliers IRCUR can handle. While the increment of sampling constant enhances the robustness rapidly when cI is small, the tolerance of corruption asymptotically reaches its limit later. This suggests that one should pick a medium value for cI to ensure robustness and efficient implementation; we find cI ∈ [3, 5] is a good balanced choice in practice. We also observe that IRCUR-R has slightly better recovery ability than IRCUR-F, which is expected. On the other hand, the maximum tolerance of corruption gets lower as the rank increases since higher rank creates harder problems.

Computational Efficiency
We evaluate the computational efficiency of the test algorithms. The experimental settings and results are summarized in Figure 2. The left subfigure shows that both versions of IRCUR have substantial speed advantage against AccAltProj and GD, especially when n is larger. The middle subfigure confirms that the computational complexity of IRCUR is indeed O(n log 2 (n)) and IRCUR-F is sightly more efficient than IRCUR-R. The right subfigure shows the linear convergence of all test algorithms, wherein IRCUR has lowest runtime per iteration.

Video Background Subtraction
We apply the test algorithms to the task of video background subtraction. Two popular videos, shoppingmall and restaurant 1 are used as benchmarks. We vectorize and stack the frames of a video to form the data matrix D, then apply RPCA to separate the background (i.e., the low rank component L) and the foreground (i.e., the sparse component S). For a stable static background, we set r = 2 and cI = 4 in this test. The  Figure 3: Video background subtraction results. The first row is corresponding to a selected frame from shoppingmall, and the second row is corresponding to a selected frame from restaurant. The first column is the original frames, the next 2 columns are the foreground and background outputted by IRCUR-F, and the last 2 columns are the foreground and background outputted by IRCUR-R.
video size information and runtime results for all test algorithms are recorded in Table 1. We again confirm that both versions of IRCUR are substantially faster than AccAltProj and GD in this real-world benchmark. Furthermore, all the test algorithms achieve visually desirable results under aforementioned parameter setting. For IRCUR, we present the separated background and foreground for a selected frame from each video in Figure 3. One can see that both versions of IRCUR enjoy crisp static backgrounds in both videos. Due to the page limit, we defer more and enlarged visual results to the supplementary document.