Outlier-Robust Multi-Aspect Streaming Tensor Completion and Factorization

1 Introduction

Tensors (or multi-way arrays) are higher order generalization of vectors and matrices. They are used to model multi-modal information in many real-world applications including recommendation systems [Karatzoglou et al., 2010], image processing [Lahat et al., 2015], social network analysis [Ermis¸ et al., 2015] and chemometrics [Mørup and Hansen, 2009]. A fundamental challenge is how to process, analyze and utilize such high-volume tensor data effectively and efficiently. In order to tackle this challenge, many algorithms such as CANDECOMP/PARAFAC (CP) factorization [Kolda and Bader,

∗Contact Author

2009] and Tucker factorization [Kolda and Bader, 2009] have been developed.

Due to the nature of data or cost of the data collection, the data tensor may be incomplete (i.e., tensor with missing values or partially observed tensor). The problem of filling the missing entries of incomplete tensor is called tensor completion, which has received much attention from researchers and practitioners in many domains such as recommender systems [Rendle et al., 2009] and image recovery [Liu et al., 2013].

(a) Single-Aspect Streaming Tensor

(b) Multi-Aspect Streaming Tensor

Figure 1: Illustration of Single-Aspect Streaming Tensors (traditional settings) and Multi-Aspect Streaming Tensors (t denotes the time-step) with missing data. The grey blocks represent the data at the previous time steps, while green ones represent the data at current time step. Black shapes with question marks are the missing entries.

For many modern applications, data is collected in a streaming fashion carrying time-varying information. For instance, in restaurant recommendation system, data is naturally structured as tensor with three modes (or dimensions), and number of users and restaurants may grow as time goes. In Facebook, users post 684,478 messages per minute. Because of increasingly massive amount of newly created data instances in streaming setting, using static algorithms for tensor completion and factorization is computationally expensive from time and space complexity perspectives. Learning latent factors in a streaming manner has benefit of avoiding irrelevant information and computational overhead of the longpast data. This raises a fundamental challenge for machine learning: how to obtain the complete tensor effectively and efficiently given incomplete streaming tensor.

Existing work on streaming tensor completion and factorization can be roughly classified into two categories: singleaspect approach and multi-aspect approach. Single-aspect approach [Kasai, 2016; Mardani et al., 2015] is usually based on the assumption that the streaming tensor evolves over one mode, while multi-aspect approach [Song et al., 2017; Nimishakavi et al., 2018] is built upon the fact that the tensor data may stream across multiple modes. Figure 1 shows an example of single-aspect and multi-aspect streaming tensors with missing values.

Following the first assumption, single-aspect streaming tensor completion and factorization have been extensively studied under the CP factorization. There are many work have been done for this track [Kasai, 2016; Mardani et al., 2015; Zhou et al., 2016]. These approaches iteratively update the latent factors of static modes with partial decomposition of the incremental mode.

Different from single-aspect approach, multi-aspect streaming tensor completion and factorization usually suffer from a high computational and space complexity as multiple modes evolve over time. Also, the incremental part may not be well-structured and therefore cannot be directly factorized. Song et al. [Song et al., 2017] proposed a method, named as MAST, that combines dynamic CP factorization and low-rank tensor completion to fill the missing entries of the incomplete streaming tensor and learn its latent factors without the whole tensor data reconstruction from scratch. MAST uses dynamic tensor factorization to update the model with the incremental part of multi-aspect streaming tensor sequence and track its low-rank subspace via nuclear norm for completion purpose. By leveraging the learned model in the previous time step, MAST accelerates filling the missing entries, while having comparable performance to traditional methods.

A challenging problem may arise when the incomplete tensor is contaminated by outliers. Outliers can be defined as significant deviations from real data values. Due to sensor failures, malicious tampering and system errors, realworld data can encounter outliers [Balcan and Zhang, 2016]. With outliers, although errors exist only on several parts of data, they can adversarially affect learning from data. Existing streaming tensor completion and factorization methods concentrate only on clean multi-aspect streaming tensors or outlier-corrupted single-aspect streaming tensors, which leads to the problem of outlier-contaminated streaming tensor completion and factorization. Hawkins et al. [Zhang and Hawkins, 2018] developed a method, named as BRST, for erroneous single-aspect streaming tensor via variational Bayesian inference. BRST models the whole single-aspect erroneous streaming tensor as the sum of a low-rank streaming tensor and a time-varying sparse component. To learn these two parts, bayesian statistical model is applied. This model imposes low-rank and sparsity using the hyperparameters and prior density functions. Then, the posterior density function of the latent factors is obtained by variational Bayesian method.

Nonetheless, BRST cannot be directly applied to multiaspect erroneous streaming tensors. To handle outliers in incomplete multi-aspect streaming tensors, we propose a novel Outlier-Robust Multi-Aspect Streaming Tensor Completion and Factorization (OR-MSTC) method based on CP. Different fromMAST, OR-MSTC decomposes the data tensor into a low-rank clean tensor and an error tensor, and imposes ℓ2,1 on the error tensor, aiming to learn better latent factors, while improving the robustness of completion and factorization against outliers. A clean tensor is the low-rank tensor without outlier for an outliers-corrupted tensor.

The error tensor captures the outliers in the tensor data, while the clean low-rank tensor encodes the underlying lowrank structure of the streaming tensor data. Without loss of generality, we assume that outliers are distributed on the 2nd mode of the tensor. We can hence use ℓ2,1 to characterize this sparsity property. The proposed OR-MSTC method imposes ℓ2,1 norm on error tensor. The decomposition into two parts and applying ℓ2,1 norm make the corresponding objective function of the proposed OR-MSTC method challenging to optimize. We present a new efficient optimization algorithm based on Alternating Direction Method of Multipliers (ADMM) [Boyd et al., 2011] to solve it. Our main contributions can be summarized as follows:

• To the best of our knowledge, OR-MSTC is the first work that recovers clean tensor under low-rank subspace, while capturing outliers for outlier-contaminated multi-aspect streaming tensors.
• OR-MSTC is the first method based on CP factorization that isolates the clean low-rank tensor from outliers for an outlier-corrupted multi-aspect streaming tensors.
• We propose a new efficient optimization algorithm based on ADMM to solve the objective function.
• Through extensive experiments on real-world datasets, we show that OR-MSTC is superior to several state-ofthe-art methods in tensor completion and factorization and robust against outliers.

2 Outlier-Robust Multi-Aspect Streaming Tensor Completion and Factorization

We begin with some necessary notations and concepts of tensor algebra. Table 1 lists basic symbols that will be used throughout the paper. Specifically, an element of a vector x, a matrix X, or a tensor X is denoted by xi, xi,j , xi,j,k etc., depending on the number of modes.

Given an N th order tensor X ∈ RI1×I2×···×IN and an integer R (rank), the CP factorization is defined by latent factor matrices Xn ∈ RIn×R for n ∈ [1 : N], respectively, such that X = ∑Rr=1 xr1 ◦ xr2 ◦ · · · ◦ xrN = JX1,X2, · · · ,XNK, where xrn ∈ RIn is the rth column of the latent factor matrix Xn, and J·K is used for shorthand notation of the sum of rank-one tensors. For convenience, in the following, we write x1 ◦ x2 ◦ · · · ◦ xN as ∏Nn=1 ◦xn.

We will first systematically propose a novel outlier-robust streaming tensor completion and factorization method, followed by a new efficient iterative algorithm to solve the formulated objective function.

Given a multi-aspect streaming N th order tensor X ∈ RI1×I2×···×IN with missing entries and possibly outliers, where X (t) is the snapshot of the tensor at time-step t and for any t ∈ Z+, X (t-1) is a sub-tensor of X (t) , the goal is to recover the missing and corrupted entries, and learn its latent factors.

2.1 Outlier-Robust Low-Rank Tensor Completion

The low-rank tensor completion aims at learning the low-rank complete tensor given an incomplete tensor, can be formulated as a rank-minimization problem as follows:

where X c ∈ RI1×I2×···×IN denotes the complete low-rank tensor, X indicates the incomplete tensor which is a partial observations of X c . Ω is a binary tensor with the same size as X c whose entries indicate whether each corresponding entry in X c is observed or not, where ωi1,··· ,iN = 1 if xci1,··· ,iN is observed and ωi1,··· ,iN = 0 otherwise. ∗ denotes the elementwise product. Since rank(X c ) is non-convex, the objective function in Eq. (1) is an NP-hard problem [H˚astad, 1990]. One way is to replace rank(X c ) with the summation of the nuclear norm of its factorized matrices and factorization error as follows [Liu et al., 2013]:

where αn are hyperparameters that balance impact of each mode.

X c might be also erroneous so that some of its entries are perturbed. As shown in Figure 2, we assume that the complete N th order tensor X c is a mixture of clean low-rank tensor L and error tensor E that captures sparse outliers in X c . For example, in network traffic, the clean low-rank tensor is the usual traffic flow and outliers are anomalies [Zhang and Hawkins, 2018]. This decomposition technique is called lowrank tensor decomposition and is formulated as X c = L+E.

Without loss of generality, our assumption is that outliers are distributed on 2nd dimension of tensor X c (and thus X) and they are very sparse in comparison to the data size. We can thus impose ℓ2,1 on the error tensor to characterize this sparsity property, which is defined as ‖E‖2,1 . Our proposed method can be easily extended to deal with multi-way outliers. Specially, similar to 2nd mode, we can make similar consideration and formulation for other modes. It can be seen that the sparsity-inducing property of ℓ2,1-norm pushes E to be sparse in each slice along 2nd mode. More specifically, 2nd mode shrinks to zero if the corresponding entries belongs to outliers. Using this decomposition, tensor completion is a rank minimization problem that can be approximated using tensor nuclear norm. The objective function for low-rank clean tensor decomposition can be stated as follows:

where ‖L‖∗ is the convex envelop of the tensor average rank within the unit ball [Lu et al., 2016] and λ is the hyperparameter to balance impact of the error tensor.

By replacing tensor nuclear norm minimization with Eq. (2), Eq. (3) can be written as follows:

We refer to the above loss as LORLTC. Due to outliers and missing entries in X, latent factors A1,A2, ..., and AN may not be exact. Thus, instead of using binary tensor Ω in Eq. (2), we define weighting tensor W with the same size as X c and follow an iterative approach to obtain better latent factors by taking incomplete nature of data into account as follows:

The weighting tensor W represents the reliability of the entries in tensor X . If a tensor entry is not missing, we set initial weight to not be 0. because that entry might be contaminated by outlier, which is set to 0.001 in our experiments. Otherwise, the weight of entry is set to be 0 as missing entry is not reliable. Having latent factors A1,A2, ..., and AN, the goal is to keep the most reliable entries in X, while iteratively updating those that have very low reliability. Inspired by [Shao et al., 2015], we use W to update X c as follows:

where 1 is a tensor with all ones. As the iteration goes up, more weight is assigned to less reliable entries. We then revise the weight tensor at the end of each iteration as follows:

Algorithm 1 Outlier-Robust Streaming Tensor Completion and Factorization (OR-MSTC)

Input: X, Ω

• 1: A(0) n = ˜ An, A(1) n = rand(dn,R), Zn =Yn = 0 where dn denotes the size of incremental part in mode n
• 2: repeat
• 3: η =min(η ∗ ρ, ηmax)
• 4: for n = 1 to N do
• 5: Update A(0) n and A(1) n using Eq. (13)
• 6: end for
• 7: Update X c using Eq. (6)
• 8: UpdateW using Eq. (7)
• 9: Update E = prox l21λ (JA1, ...,ANK -X c )
• 10: for n = 1 to N do
• 11: Yn =Yn + η(Zn -An)
• 12: Zn = SV Tαn η (An - Yn η )
• 13: end for

where β is a decay parameter, which we fix to 0.97 in our experiments. As the number of iterations increases, W converges to 1. The resulted objective function can be formulated as minX c ,E,A1,...,AN LORLTC, where X c is updated using Eqs. (6) and (7).

2.2 Dynamic CP Factorization

In streaming settings, a dynamic CP method obtains latent factors in an efficient and effective fashion. We use ˜ X and X to represent two consecutive snapshots X (t-1) (tensor at time step t - 1) and X (t) (tensor at time step t) such that ˜ X ⊊ X. The general CP factorization objective function for tensor X is defined as minA1,...,AN ‖X - JA1, A2, ..., AN K‖2F .

In CP factorization for ˜ X, we have ˜ X ≈ J ˜ A1, ˜ A2, ..., ˜ AN K ( ˜ A1, ˜ A2, ..., and ˜ AN denote latent factors for the former snapshot). We take advantage of former snapshot latent factors to learn latent factors for X by incorporating ˜ X into CP factorization of X. This in fact accelerates factorization of current snapshot. Dynamic CP factorization for X can be hence rewritten as summation of two terms: factorization related to former snapshot and factorization for the incremental part i.e., part in snapshot at time step t that didn't exist in snapshot at time step t - 1. The resulted objective function for X can be stated as follows (θ indicates entries of the incremental part):

whereA⊤n = [A(0) n ⊤, A(1) n ⊤]. Let dn denote increment along nth mode at time step t in comparison to time step t - 1. A(0) n ∈ R(In-dn)×R refers to partition of the latent factor w.r.t. tensor snapshot at time step t - 1 in the tensor snapshot at time step t, while A(1) n ∈ Rdn×R is for incremental part. Using CP factorization for ˜ X in Eq. (8) results in the following objective function:

where µ ∈ [0, 1] is a forgetting factor to control trade-off between previous and current factorizations for time step t - 1. As shown in Figure 3, second term can be defined as a summation of CP factorization for a set of sub-tensors [Song et al., 2017].

2.3 ProblemFormulation

We combine outlier-robust low-rank tensor completion and dynamic CP factorization as follows:

where X c is updated using Eqs. (6) and (7). In the second term, we use X c - E as substitution of X .

2.4 Optimization Procedure

Using ADMM, by inducing Z⊤n = [Z(0) n ⊤, Z(1) n ⊤] ∈ RIn×R and Lagrange multipliers Y⊤ n = [Y(0) n ⊤, Y(1) n ⊤] ∈ RIn×R, the augmented Lagrangian function for Eq. (10) is as follows:

where η is penalty parameter and X c is updated by Eqs. (6) and (7). Let P1, P2, P3 denote the following (L = X c - E):

The update rules for A(0) n and A(1) n are as follows:

where S0n = {(s1, ..., sN)| ∑Nk=1 sk ̸= 0, sk ∈ {0, 1}, sn = 0} and S1n = {(s1, ..., sN)|∀k ∈ [1, ..., N], sk ∈ {0, 1}, sn = 1} are sets of tuples to refer to subtensors. Zn can be updated as Zn = SV Tαn η (An - Yn η ), n = 1, 2, ...,N, where SV Tαn η is the singular value thresholding operator defined as SV Tω(A) = U(diag{σi - ω})+V⊤, where U(diag{σi})1<=i<=rV⊤ is the singular value decomposition of A. X+ = max{X, 0}, where max{., .} is an element-wise operator. The update rule for E can be stated as E = prox l21λ (JA1, ...,ANK -X c ) for incremental part. For 3rd order tensor R, prox l21λ (R) approximates Ri = ||Ri ||F -λ ||Ri ||F Ri , if ||Ri ||F > λ; 0, otherwise (Ri = R(:, i, :)). The framework of OR-MSTC is summarized in Algorithm 1. The convergence condition is that the relative changing of tensor X c in two consecutive iterations is smaller than the tolerance.

3 Experimental Evaluation

To evaluate the performance of the proposed OR-MSTC method, we conduct experiments on real-world datasets and compare with the following seven baselines: (1) Static NNCP is based on trace norm constrains and optimized by ADMM [Liu et al., 2015]. (2) Static CP-ALS employs CP factorization and optimized by Alternating Least Square (ALS) [Bro, 1998]. (3) Static Tucker-ALS is based on tucker factorization and optimized by ALS [Comon et al., 2009]. For all static baselines, the former-step factors is given as the initialization for the current step to improve learning. (4) OLSTEC is a single-aspect streaming tensor completion method optimized by ALS [Kasai, 2016]. (5) TeCPSGD that is pioneering CP factorization work for single-aspect streaming tensor completion and optimized by stochastic gradient descent [Mardani et al., 2015]. (6) MAST is a multi-aspect streaming tensor completion approach that assumes no error in the data [Song et al., 2017]. (7) BRST is a error-robust single-aspect streaming tensor completion approach based on bayesian inference [Zhang and Hawkins, 2018].

Similar to [Song et al., 2017], in order to apply singleaspect streaming tensor methods on multi-aspect streaming tensor, we divide the multi-aspect incremental parts into several single-aspect streaming ones and update them in random order in each iteration until convergence.

We apply grid search to identify optimal values for each hyperparameter from {10-9 , 10-8 , ..., 108 , 109 }. The tolerance rate is set to 10-4 , the maximum number of iterations to 500 for all the methods. The rank is tuned using 10 ranks varying from 5 to 40 based on relative error which we define later. In OR-MSTC and MAST, αn = 1 10N , n = 1, ...,N, η = 10-4 , ρ = 1.05 and ηmax = 106 . We tuned forgetting factor µ in our method and MAST based on the missing percentage. Following [Song et al., 2017], the initial completion and warm start matrices are calculated using NNCP method. We use Relative Error (RE) as performance evaluation metric, which is defined as follows:

where reconstructed tensor denotes the output complete tensor produced by the methods and real tensor indicates the true data tensor (i.e., ground truth). RE measures the deviation from the ground truth. The efficiency is measured by Average Running Time (ART) and defined as 1 T ∑Tt=1 RTt, where RTt is the running time at time step t. Each experiment is repeated for 5 times, and mean of the metric in each dataset is reported. The datasets used in the paper are Cardiac MRI (CMRI) [Sharif and Bresler, 2007], downsampled Yelp [Jeon et al., 2016] and Network Traffic (NT) [Lakhina et al., 2004], and summarized in Table 3. Following [Song et al., 2017], we randomly cover a percentage of data (missing percentage) ({10%, 20%}) and consider the remaining entries as observed information. We generate outliers on 2nd mode with mean of 1st slice of tensor data, magnitude of 5, variance of 0 and outlier percentage of {10%, 20%}.

Figure 4 and 5 show performance of the methods with various missing and outlier percentages. We don't report results for BRST on Yelp and NT as its RE is very large in comparison to the proposed OR-MSTC method and other baselines. On all datasets, OR-MSTC outperforms the dynamic baselines, especially BRST, error-robust streaming scheme, with varying percentages of missing entries and outliers. On CMRI and Yelp, OR-MSTC is significantly better than dynamic baselines (in average at least 0.4 less). In Yelp, number of businesses consistently increases as time evolves. This pattern reduces the quantity of substitution on corrupted and missing entries and results in noticeable improvement in performance of the proposed method. On NT, OR-MSTC always obtains nearly perfect result. The proposed OR-MSTC method shows strong stability in performance on all datasets as time goes. On Yelp, MAST, TeCPSGD and static methods suffer from high instability in performance. This is because incremental part at first 15 time steps contains high number of outlier and missing entries, and they account for most of relative errors in baselines. For various ratios of missing entries and outliers, OR-MSTC outperforms static baselines. In Table 2, OR-MSTC's ART is far better than previous error robust streaming scheme (BRST), but not necessarily better than the other streaming baselines.

On CMRI, performance of all methods is unaffected by varying ratio of missing data with the same outliers, while this observation is not true on varying ratios of outliers with the same missing data. Nevertheless, the proposed OR-MSTC method is more effective and robust with different ratios of missing data and outliers. Yelp is sparse in comparison to CMRI, and adding outliers and missing data does not therefore affect the dataset a lot. As ratio of missing data increases, performance of all methods decreases, but this drop is not high. Compared to baselines, the proposed OR-MSTC method shows stable performance. On NT, TeCPSGD and OLSTEC have large variation and drop in performance with increase in ratio of outliers and missing data compared to our method. They thus don't identify anomalies in network well.

4 Conclusion

We developed a method named OR-MSTC for multi-aspect streaming tensor completion and factorization. OR-MSTC has several advantages over existing streaming tensor completion and factorization methods. First, it effectively captures low-rank subspace of multi-aspect streaming tensor for completion and factorization tasks. Second, it handles outliers well. Third, iterative scalable optimization framework is proposed for OR-MSTC. Compared to existing streaming tensor completion and factorization methods, OR-MSTC showed better performance on three datasets.

Acknowledgements

This work is supported in part by NSF through grants IIS1526499, IIS-1763325, CNS-1626432, Guangdong Province through grant 2017A030313339 and NSFC 61672313.