Online Robust Low-Rank Tensor Learning

1 Introduction

The explosion of multidimensional data raises significant demand on efficient data analysis techniques. In recent years, low-rank tensor modeling methods have received increasing attention [Goldfarb and Qin, 2014; Lu et al., 2016] as they can discover intrinsic structure of tensors and provide more transparent and consolidated knowledge of the data. However, many realistic multi-dimensional data (e.g., surveillance videos) are usually generated in a streaming way from dynamic environments, which poses following new challenges: (1) how to develop scalable methods that can process large-scale dynamic tensors efficiently; (2) how to deal with data contamination such as noise or malicious outliers in the modeling process.

To handle noisy tensors, two batch robust tensor methods were developed: High-Order Robust Principal Component Analysis (HORPCA) [Goldfarb and Qin, 2014] and Tensor RPCA [Lu et al., 2016]. They both generalize RPCA [Cand`es et al., 2011] from matrix cases to tensors, i.e., solving the tensor RPCA problem minZ,E ‖Z‖∗+‖E‖1 with X = Z + E , to learn the inherent low-rank structure. Although they can handle noisy data, those methods are not scalable to largescale data due to their high memory cost which increases rapidly with the data size, and always seek a static lowrank component that cannot model dynamics of streaming data thus providing inferior performance. A proper method to analyze large-scale noisy tensor data is still absent.

In this work, we propose an Online Robust Low-rank Tensor Modeling (ORLTM) method for learning low-rank structures of tensors from streaming noisy tensor data. Different from batch approaches, ORLTM isolates the lowrank component into two factors in each mode m, i.e., Zm = WmRm where Wm models the low-rank subspace basis and Rm incorporates the corresponding tensor data representation coefficients w.r.t. Wm. Here, the basis and data representations are decoupled. Thus, the basis size is typically small and can be stored with low memory cost while data representation can be updated online. In this way, ORLTM saves the memory cost dramatically without performance drop.

Existing tensor RPCA approaches generally assume data reside in a single low-rank subspace. This may not hold for many realistic tensor data with complex inherent structures, e.g., those drawn from a union of multiple subspaces. Previous methods as well as the vanilla ORLTM are less competent in such circumstances. Hence, we further introduce a dictionary learning component to ORLTM. Dictionary provides a more flexible set of basis to represent the low-rank component of complex tensors. Instead of seeking a single set of basis, through learning proper dictionaries, ORLTM can capture sophisticated low-rank structures (e.g., mixture of multiple low-rank components) for better tensor modeling.

ORLTM has wide applications in large-scale and dynamic tensor data analysis. In this paper, we empirically examined its performance on both synthetic data and video background subtraction, the results of which have well justified the effectiveness of ORLTM.

2 RelatedWork

Low-rank models are useful for robust data recovery. A typical one is Robust Principal Component Analysis (RPCA) [Cand`es et al., 2011] that decomposes a given matrix into a low-rank component and a sparse component. Another is Low-Rank Representation (LRR) [Liu et al., 2013] which considers the data drawn from a union of multiple subspaces. To this end, many variants have emerged, e.g., [Liu and Yan, 2011] employed both observed and hidden data to enhance LRR; [Yin et al., 2015] developed a dual graph regularized LRR. However, these methods require tremendous memory to store all the samples. Hence, some online learning methods were developed, e.g., [Feng et al., 2013] proposed an Online Robust PCA (ORPCA) via stochastic optimization; [Zhan et al., 2016] studied online RPCA based on the recursive projected compressive sensing framework; [Shen et al., 2016] designed an online LRR implementation with convergence guarantees. Nevertheless, they fail to exploit the low-dimensional tensor structures.

Robust tensor models can utilize information across all the modes. Two principled forms are CANDECOMP/PARAFAC (CP) decomposition and Tucker decomposition [Sun et al., 2008]. For example, [Zhou et al., 2016] proposed an accelerated CP method for dynamic tensors. Based on the latter, Higher-Order RPCA (HORPCA) was presented in [Goldfarb and Qin, 2014]. Moreover, [Zhang et al., 2013] proposed a low-rank three-order tensor recovery method for rectification and alignment, while [Lu et al., 2016] proposed one tensor RPCA (TRPCA) method for image denoising using the framework in [Kilmer and Martin, 2011]. Nonetheless, these methods need large memory for batch optimization and cannot handle samples sequentially. In addition, there are rare works concerning robust online tensor analysis, and the ever-known one is Online Stochastic Tensor Decomposition (OSTD) [Sobral et al., 2015], which neglects the data with the mixture structure of multiple subspaces.

3 Notations and Preliminaries

Tensors are denoted by calligraphic letters, e.g., A, matrices by boldface uppercase letters, e.g., A, vectors by boldface lowercase letters, e.g., a. The order or mode of a tensor equals the number of dimensions. For a third-order tensor X ∈ Rn1×n2×n3 , its (i, j, k)-th entry is xijk, and its fiber is a column vector X(i, j, :); X(i, :, :) is the i-th horizontal slice; X(:, j, :) is the j-th lateral slice; X(:, :, k) or Xk is the k-th frontal slice. For the matrix, Aij denotes the (i, j)-th entry, Am without bracket denotes a traditional matrix.

Tensor Unfolding [Kolda and Bader, 2009] The modem unfolding of A ∈ Rn1×n2×...×nM , A(m) , is obtained by arranging the mode-m fibers in the matrix columns, i.e., unfoldm(A) = A(m) ∈ Rn¯m×nm , where n¯m = ∏Mj̸=m,j=1 nj and ∏Mj=1 nj = n1 × n2 × . . . × nM.

Tensor Folding [Kolda and Bader, 2009] Folding the mode-m unfolding of tensor A is defined as fold(A(m) ) = Am, which returns the corresponding tensor A in mode m.

Tensor Vectorization [Kolda and Bader, 2009] The vectorization of A is denoted by vec(A), which stacks the elements in a column vector. For mode-m unfolding, vec(Am) (a.k.a. vec(A(m) )) is done by aligning the columns of mode-m unfolding into a long column vector.

Mode-m Product [Goldfarb and Qin, 2014] The mode-m product of tensor A and matrix Um ∈ Rnm×nc on mode m is A×mU, which is defined as (A×mUm)(m) = A(m)Um. Tensor Norms [Goldfarb and Qin, 2014] The ℓ1-norm is ‖A‖1 = ‖vec(A)‖1 = ∑ijk |aijk|, and the Frobenius norm is ‖A‖F = √(vec(A)⊤vec(A)) = √(∑ijk a2ijk ). These norms can degenerate to the matrix or vector norms.

Tensor Rank [Goldfarb and Qin, 2014] The mode-m rank of a tensor A denoted by rankm(A) is the column rank of A(m) , and the set of M mode-m ranks is called Tuckerrank. However, it is an NP-hard problem and one often uses its convex surrogate CTrank(A) in practice [Hillar and Lim, 2013; Yu et al., 2015]. CTrank sums M nuclear norms of the mode-m unfoldings, i.e., CTrank(A) := ∑Mm=1 ‖A(m) ‖∗.

4 Online Robust Low-Rank Tensor Modeling 4.1 Formulation of ORLTM

We now introduce the objective formulation by starting with reviewing the batch-based High-Order RPCA (HORPCA) [Goldfarb and Qin, 2014]. This method factorizes the input tensor into a low-rank component plus a sparse noise component accounting for gross corruption or outliers, and its slice-wise model is formulated as

where λ1>0, L=fold(L(m) ), ∑Mm=1 ‖L(m) ‖∗ = CTrank(L), X,L, E ∈ Rn1×n2×...×nM are the observed noisy tensor, its low-rank component and the sparse component, respectively. On different mode-m unfoldings X(m) , the problem (1) leads to different low-rank components and sparse components. Thus we adopt the variable-splitting technique and introduce a set of auxiliary tensor variables, i.e., {Lm}Mm=1 , {Em}Mm=1 , for L and E respectively. So the problem (1) can be reformulated as

where X = fold(X(m) ), and each unfolding X(m) in all modes should return the common tensor X through the f old(·) operator. The tensors L and E can be approximated by average pooling, i.e., L = 1M ∑Mm=1 Lm, E = 1M ∑Mm=1 Em.

However, the constraints in (2) enforce the unfolding matrix to reside in a single low-rank subspace, failing to consider the situation where the data are with the mixture structure of several subspaces. Thus, this strict constraint hardly holds when data are drawn from a union of multiple subspaces, leading to degenerated performance. Therefore, to better discover the robust low-rank structure, we introduce a dictionary D(m) ∈ Rn ¯m×nm . The low-rank property is preserved by enforcing nuclear norm on Zm ∈ Rnm×nm , i.e., ‖Zm‖∗. As the inequality rank(D(m)Zm) ≤ rank(Zm) always holds, minimizing the nuclear norm of Zm amounts to bounding the rank of the clean data L(m) if it suffices to L(m) ≜ D(m)Zm, hence a low-rank structure. Thus, we can reach the following function

The product of D(m) and Zm yields the unfolding L(m) . To dynamically handle streaming data, we relax the constrains for online optimization by treating them as quadratic penalty terms, leading to

where λ2 > 0 governs the role of the nuclear norm term.

To process samples on the fly, we propose to employ the bi-factor form Zm = WmRm [Srebro et al., 2004], where Wm ∈ Rnm×p and Rm ∈ Rp×nm with p ≪ min(nm, n ¯m). Thus the rank of Zm is upper bounded by p. As indicated by [Fazel et al., 2001; Recht et al., 2010], minimizing ‖Zm‖∗ is equivalent to minimizing ‖Wm‖2F and ‖Rm‖2F simultaneously. So the unconstrained function in (4) can be reformulated as a non-convex optimization problem

which is the objective function of ORLTM, and updating the entries in Zm amounts to sequentially updating the rows of Wm and the columns of Rm.

The objective function (5) reveals a fact that the sizes of unfolding matrices increase with nM and the dictionary D(m) can only be partially accessed in each iteration. Besides, the row vectors of W(m) are coupled together as being left multiplied by the dictionary. To tackle these issues, we introduce another set of auxiliary variables B(m) = D(m)Wm ∈ Rn¯m×p , m = 1, 2, . . . ,M, to approximate the recovery part (X(m) -E(m) ) by B(m)Rm. This suggests that {B(m) }Mm=1 can be treated as Reinforced Basis Dictionaries while {Rm}Mm=1 are the coefficients or low-dimensional representation. Compared with D(m) , the dictionary B(m) is iteratively strengthened by respecting the two factored lowrank components of Zm. Therefore, we get an equivalent reformulation

where λ3 > 0 is used to control the reconstruction quality of the basis B(m) . This formulation is not only more appropriate for learning from streaming data, but also more informative since it explicitly models the basis of the union of multiple subspaces in all modes of tensor data. Thus, it achieves more accurate tensor subspace recovery and better tensor low-rank modeling.

4.2 Online Optimization for ORLTM

This section develops an efficient online optimization algorithm for optimizing the objective function in (6). This algorithm employs the stochastic optimization technique. First, two functions in mode m are defined as

where xm, em, dm, rm are column vectors of matrices X(m) , E(m) , D(m) ∈ Rn ¯m×nm , Rm ∈ Rp×nm , respectively, and wm is the row vector of Wm ∈ Rnm×p . These formulations allow us to rewrite (6) as

which amounts to minimizing the empirical loss function:

Now, we describe how to sequentially optimize the variables in (6), and we adopt the alternating method, i.e., solving one variable while holding the rest. The whole algorithmic framework is given in Algorithm 1, which consists of following variable updatings at each iteration t.

Computing rmt , emt . Given B(m) t-1 in the previous iteration for mode-m unfolding, we obtain the optimal {rmt , emt } by solving

The above problem has a closed-form solution rmt when em is fixed, i.e.,

and w.r.t. fixed rm, the local minimizer of em is given by the soft-thresholding operator [Hale et al., 2008]:

Both rmt and emt can be optimized by the coordinate descent algorithm [Wright, 2015] in an efficient way.

Computing wm t . Define an accumulation matrix G(m) t-1 = ∑t-1 i=1 dmi wm i ∈ Rn¯m×p , where Gm0 = 0, and then we can figure wm t out by minimizing

which leads to the solution in closed form:

While wm t depends on the entire dictionary D(m) , we only have access to the current atom dmt , which thus desires a proper estimation to approximately solve each wm t . The intuition behind the solution is that after observing the t-th atom dmt , we only update wm t by keeping others, and each wm t is sequentially updated only once after revealing all the atoms. The strategy is essentially the one-pass block coordinate descent algorithm, and can be flexibly refined as a multiple-pass version in demanding practical applications.

Computing B(m) t . The reinforced basis dictionary can be updated by minimizing the following surrogate function

One key to obtain efficient dictionary updates is that the surrogate ht(·) can be summarized by a few sufficient statistics updated iteratively, i.e., it is possible to describe ht(·) without explicitly storing previous samples [Mensch et al., 2016]. Actually, we can define two accumulation matrices:

Then, the solution to (15) is expressed by

where Φmt ∈ Rp×p and Θ(m) t ∈ Rn¯m×p are independent of the sample size N. This makes our approach potentially scalable to large data size.

In summary, the pipeline of our ORLTM method is described in Algorithm 1. For an M-th order tensor X, nM indicates the number of unit sub-tensor samples, i.e., one (M-1)-th order sub-tensor. If the number of unit sub-tensor samples is n, i.e., the entire tensor size, then the number of sub-tensor samples is N = n nM . At each iteration, ORLTM admits one sub-tensor, either single unit sub-tensor sample (nM = 1) or multiple unit sub-tensor samples.

Algorithm 1 Online Robust Low-Rank Tensor Modeling

Input: ObservedM-th order tensor X ∈ Rn1×n2×...×n and pre-given dictionary D, tradeoff parameters {λ1, λ2, λ3}, target rank p, sub-tensor size nM.

Output: L, B, E, {Wm,Rm}Mm=1 .

• 1: Initialize: Set all entries of {L, E} ∈ Rn1×n2×...×n , Wm ∈ Rnm×p , Rm ∈ Rp×nm , G(m) 0 , Θ(m) 0 , Φm0 to zero, initialized B ∈ Rn1×n2×...×p , N = n nM .
• 2: for t = 1 to N do
• 3: for m = 1 to M do
• 4: Access the t-th sample xmt by vec(X(m) t ).
• 5: Reveal the t-th atom dmt by vec(D(m) t ).
• 6: Obtain the coefficients rmt and sparse vectors emt by optimizing ˜ ℓ(xmt ,B(m) t-1 , rm, em) and utilize coordinate descent algorithm for Eqs. (11)(12) to derive the solutions.
• 7: Compute the coefficients wm t by minimizing ˜ ℓ2(dmt ,B(m) t-1 ,Gmt-1 ,wm), leading to Eq. (14).
• 8: Update the accumulation matrices as

• 9: Update the dictionary B(m) t by Eq. (16).
• 10: Update low-rank and sparse components by L(m) t ←B(m) t rmt and E(m) t ← emt , respectively.
• 11: end for
• 12: end for
• 13: Obtain low-rank tensor L, sparse tensor E, reinforced dictionary tensor B by average pooling on modem foldings, i.e., L = 1M ∑Mm=1 fold(L(m) ), E = 1M ∑Mm=1 fold(E(m) ), B = 1M ∑Mm=1 fold(B(m) ).

4.3 Complexity Analysis

Computational Complexity. There are four variables in each iteration. First, it is cheap to compute {rmt , emt } as one may utilize the block coordinate descent method in [Richt´arik and Tak´aˇc, 2014] which enjoys linear convergence. Second, wm t needs O(n¯mp) to conduct matrix-vector multiplication, where the rank p ≪ min(n¯m, nm). Besides, the sub-tensor size nM is usually quite small. Third, it requires O(n¯mp) to update accumulation matrices G(m) t , Θ(m) t , Φ(m) t . Fourth, computing Bmt costs O(n¯mp2 ). Therefore, the total computational cost is limited.

Memory Cost. Lower memory cost is a very promising advantage of ORLTM. On the one hand, we need O(n ¯mp) to load B(m) t and xmt to obtain {rmt , emt }, and also O(n¯mp) to exploit B(m) t , G(m) t for calculating wm t . On the other hand, because the history information of ht(B(m) t ) in (15) has been recorded by the accumulation matrices G(m) t , Θ(m) t , Φ(m) t , it only costs at most O(n¯mp). Hence, ORLTM only desires O(n¯mp) memory cost for every iteration.

5 Experiments

5.1 Synthetic Data

To generate synthetic data with low-rank structure, we construct an authentic tensor S of size 50×50×30 by rank3 factor matrices as in [Sobral et al., 2015], e.g., Sm ∈ R50×3 (here m denotes mode). Each factor matrix consists of three components [sin( 2πmim 50 ), cos( 2πmim 50 ), sgn(sin(0.5πim)], where im denotes the column element index in mode m, and the product of two factor matrices yields a frontal slice of mode-3. We randomly corrupt entries of the tensor by small noise from distribution N(0, 0.01) and outliers from uniform distribution U(-|H|, |H|) (H is interval bound magnitude). For evaluation, we use Root Relative Square Error (RRSE) as the measure criterion calculated by ‖ ˆ S-S‖F ‖S‖F , where ˆ S is the estimated lowrank tensor. In all, we evaluate four batch methods, i.e., RPCA [Cand`es et al., 2011], LRR [Liu et al., 2013; 2017], HOSVD [Goldfarb and Qin, 2014], and HORPCA [Goldfarb and Qin, 2014], and also four online methods (nM=1), i.e., ORPCA [Feng et al., 2013], OLRSC [Shen et al., 2016], OSTD [Sobral et al., 2015], and our ORLTM approach. Note that batch approaches have the hindsight knowledge of all the streaming data and provides performance upper bound for online counterparts. For compared methods, we use the default parameters in their papers, and the target rank p is set to the ground truth 3. For ORLTM, λ1 = 0.5/√log(n1n2), λ2 = 1, λ3 = λ1√log(t). The comparison results with a varying corruption percentage η from 0 to 80% and different outlier intervals 2|H| in [0 : 0.2 : 2] are plotted in Fig. 1. The figure shows that ORLTM indeed outperforms all online methods no matter how the corruption percentage and the outlier interval vary, and is comparable to batch approaches. Besides, our method is surprisingly better than LRR consistently and better than RPCA when the corruption exceeds 55%, which well justifies the superiority of ORLTM by exploiting the low-rank structure in all tensor modes.

In addition, we provide quantitative analysis with different sizes (n3) for online methods in Table 1. As seen from the table, the performance of compared methods is improved by increasing tensor size, which demonstrates that more robust tensor recovery could be achieved by taking in more samples sequentially. Meanwhile, our method still performs the best no matter how the tense size varies.

5.2 Video Background Subtraction

An input video is treated as a long third-order tensor, of which each frontal slice denotes a frame and multiple frames constitute a sub-tensor. We model the background (BG, i.e., static information) by the low-rank component L while the foreground (FG, i.e., moving objects) is pushed to the sparse component E. We test on the I2R database [Li et al., 2004] which involves various indoor and outdoor environments. We use eight video sequences containing over 15,000 frames with sizes ranging from 120×160 to 256×320. Each sequence provides 20 ground-truth Foreground Mask (FM) of frames. We compare ORLTM with online methods including ORPCA [Feng et al., 2013], OLRSC [Shen et al., 2016], OSTD [Sobral et al., 2015], and batch methods including RPCA [Cand`es et al., 2011], LRR [Liu et al., 2013], HORPCA [Goldfarb and Qin, 2014], for which parameters are set as indicated in their papers. For our method, parameters are fixed as λ1 = α/√log(n1n2) (α is 0.02 for Curtain, Lobby, Watersurface; and 0.1 for the rest), λ2 = 1, λ3 = λ1√log(t) (t is the iteration number), and each sub-tensor size is set to 1. For all methods, the target rank p is empirically defined as 10, and the foreground masks are obtained by thresholding the sparse component, i.e., FMij = 1 if E2ij /2 ≥ (std(vec(E)))2 , and zero otherwise, where std(·) denotes the standard deviation. For online methods, all frames are passed two epochs. For batch methods, they are applied to each small-batch (500 frames) per time considering limited memory resource. To initialize B(m) in ORLTM, we adopt the Bilateral Random Projections (BRP) [Zhou and Tao, 2012]. Concretely, given dense matrix Γ ∈ Rd×c , its low-rank approximation is built by ˜Γ = ΓY1(Y⊤ 2 ΓY1)-1Y⊤ 2 Γ, where Y1 ∈ Rc×r , Y2 ∈ Rd×r are random matrices, and the rank r ≤ min(d, c). In tests, we found LRR performs very poorly using the data set itself as the dictionary, so we use its mean matrix.

For evaluation, we employ F-score [Brutzer et al., 2011] computed by comparing each sequence with its available ground-truth frames. The results are recorded in Table 2. From the table, we observe that ORLTM has the most encouraging overall performance across all scenes, about 10% higher than second best RPCA (batch method). Especially, our method performs much better on Bootstrap, Curtain, Fountain, Lobby, and Hall, compared to the rest, and it also does well on the other sequences as ORLTM can achieve almost the same F-score to the best ones. The superior performance of ORLTM can be attributed to two facts: one is dictionary learning in ORLTM can better model video background despite in noisy and interrupted environments; the other is that learning the lowrank structures of tensor in all modes can sufficiently exploit the underlying information of video sequences. Note that batch methods process 500 frames per time, which might slightly degenerate its performance.

To examine the efficiency, we also report the Frames Per Second (fps) of our method and two batch methods in Table 3. As shown in the table, ORLTM is very efficient, e.g., it is over 12 times faster than RPCA on Bootstrap and up to 46 times faster than HORPCA on Fountain. This is because ORLTM as an online method only needs constantly small memory and its computational cost is low in each iteration. On the contrary, batch methods have to store all frames in the memory in advance while they require expensive singular value decompositions on batch frames. All tests were run with 3.06 GHz Core X5675 processor and 24GB RAM.

In addition, we display the foreground masks of some randomly selected frames in Fig. 2, showing ORLTM outperforms others in most scenes, e.g., curtain swinging (row 2), running fountain (row 3), illumination change

(row 4). This demonstrates our method can well model background and foreground in a rich variety of scenes.

6 Conclusion

We present an Online Robust Low-rank Tensor Modeling (ORLTM) method to learn low-rank structures of streaming noisy tensor data robustly. By developing the equivalent bifactor formulation of the nuclear norm, it can process samples sequentially, thus being scalable to large-scale tensor data. The objective function is solved by stochastic optimization, and in each iteration it saves memory cost by a factor of n compared to batch methods. Synthetic studies suggest ORLTM outperforms well-established online methods and its performance is comparable to batch ones when the data are corrupted by up to 80% gross noise with large magnitude. In video background subtraction, ORLTM gains the highest average F-score over 10% higher than the state of the arts, and is significantly faster than batch methods, up to 46 times faster than HORPCA.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grants 61502131, 61572169, 61472266, Zhejiang Provincial Natural Science Foundation of China under Grant LQ15F020012, and China Scholarship Council. The work of Jiashi Feng was supported in part by National University of Singapore Startup Grant R-263-000-C08-133 and Ministry of Education of Singapore AcRF Tier One Grant R-263-000-C21-112.