GOCPT: Generalized Online Canonical Polyadic Tensor Factorization and Completion

1 Introduction

Streaming tensor data becomes increasingly popular in areas such as spatio-temporal outlier detection [Najafi et al., 2019], social media [Song et al., 2017], sensor monitoring [Mardani et al., 2015], video analysis [Kasai, 2019] and hyper-order time series [Cai et al., 2015]. The factorization/decomposition of such multidimensional structural data is challenging since they are usually sparse, delayed, and sometimes incomplete. There is an increasing need to maintain the low-rank Tucker [Sun et al., 2006; Xiao et al., 2018; Nimishakavi et al., 2018; Fang et al., 2021; Gilman and Balzano, 2020] or CP [Du et al., 2018; Phipps et al., 2021] structure of the evolving tensors in such dynamics, considering model efficiency and scalability.

Several online (we also use "streaming" interchangeably) settings have been discussed before. The most popular two are online tensor decomposition [Zhou et al., 2016; Song et al., 2017] and online tensor completion [Kasai, 2019], where the temporal mode grows with new incoming slices. Some pioneer works have been proposed for these two particular settings. For the factorization problem, [Nion and Sidiropoulos, 2009] incrementally tracked the singular value decomposition (SVD) of the unfolded third-order tensors to maintain the CP factorization results. Accelerated methods have been proposed for the evolving dense [Zhou et al., 2016] or sparse [Zhou et al., 2018] tensors by reusing intermediate quantities, e.g., matricized tensor times Khatri-Rao product (MTTKRP). For the completion problem, recursive least squares [Kasai, 2016] and stochastic gradient descent (SGD) [Mardani et al., 2015] were studied to track the evolving subspace of incomplete data.

However, most existing methods are designed for the growing patterns. They cannot support other possible patterns (e.g., missing entries in the previous tensor can be refilled or existing values can be updated) or more complex scenarios where the tensors evolve with a combination of patterns.

Motivating application: Let us consider a public health surveillance application where data is modeled as a fourthorder tensor indexed by location (e.g., zip code), disease (e.g., diagnosis code), data generation date (GD) (i.e., the date when the clinical events actually occur) and data loading date (LD) (i.e., the time when the events are reported in the database) and each tensor element stands for the count of medical claims with particular location-disease-date tuples. The tensor grows over time by adding new locations, diseases, GD's, or LD's. For every new LD, missing entries before that date may be filled in, and existing entries can be revised for data correction purposes [Qian et al., 2021]. Dealing with such a dynamic tensor is challenging, and very few works have studied this online tensor factorization/completion problem. The most recent work in [Qian et al., 2021] can only handle a special case with the GD dimension growing but not with data correction or two more dimensions growing simultaneously.

To fill the gap, we propose GOCPT - a general framework that deals with the above-mentioned online tensor updating scenarios. The major contributions of the paper are summarized as follows:

• We propose a unified framework GOCPT for online tensor

factorization/completion with complex evolving patterns such as mode growth, missing filling, and value update, while previous models cannot handle such general settings.

• We propose GOCPTE, i.e., a more memory and computational efficient version of GOCPT, to deal with cases where historical tensor elements are unavailable due to limited storage or privacy-related issues.
• We experimentally show that both GOCPT and GOCPTE work well under a combination of online tensor challenges. The GOCPT improves the fitness by up to 2.8% and 9.2% on real-world Covid and medical claim datasets, respectively. In comparison, GOCPTE provides comparable fitness scores as GOCPT with 20%+ complexity reduction compared to the baseline methods.

The first version of GOCPT package has been released in PyPI 1 and open-sourced in GitHub 2 . Supplementary of this paper can be found in the same GitHub repository.

2 ProblemFormulation

Notations. We use plain letters for scalars, e.g., x or X, boldface uppercase letters for matrices, e.g., X, boldface lowercase letters for vectors, e.g., x, and Euler script letters for tensors or sets, e.g., X. Tensors are multidimensional arrays indexed by three or more modes. For example, an Nth-order tensor X is an N-dimensional array of size I1×I2×· · ·×IN , where xi1i2···iN is the element at the (i1, i2, · · · , iN) location. For a matrix X, the r-th row is denoted by xr. We use �∗ for Hadamard product (element-wise product), ⊙ for Khatri-Rao product, and �·� for Kruskal product (which inputs matrices and outputs a tensor).

For an incomplete observation of tensor X , we use a mask tensor Ω to indicate the observed entries: if xi1i2···iN is observed, then Ωi1i2···iN = 1, otherwise 0. Thus, X �∗ Ω is the actual observed data. In this paper, Ω can also be viewed as an index set of the observed entries. We define ∥(·)Ω∥2F as the sum of element-wise squares restricted on Ω, e.g.,

where X and Y may not necessarily be of the same size, but Ω must index within the bounds of both tensors. We describe basic matrix/tensor algebra in appendix A, where we also list a table to summarize all the notations used in the paper.

2.1 ProblemDefinition

Modeling Streaming Tensors. Real-world streaming data comes with indexing features and quantities, for example, we may receive a set of disease count tuples on a daily basis, e.g., (location, disease, date; count), where the first three features can be used to locate in a third-order tensor and the counts can be viewed as tensor elements.

Formally, a streaming of index-element tuples, e.g., represented by (i1, · · · , iN; xi1···iN ), can be modeled as an evolving tensor structure. This paper considers three typical types of evolution, shown in Fig. 1.

1 https://pypi.org/project/GOCPT/

2 https://github.com/ycq091044/GOCPT

• (i)Mode growth. New (incomplete) slices are added along one or more modes. Refer to the blue parts.
• (ii) Missing filling. Some missing values in the old tensor is received. Refer to the green entries in the figure.
• (iii) Value update. Previously observed entries may change due to new information. Refer to yellow entries.

To track the evolution process, this paper proposes a general framework for solving the following problem on the fly.

Problem 1 (Online Tensor Factorization and Completion). Suppose tensor X t-1 ∈ RI t-1 1 ×···×I t-1 N admits a lowrank approximation Y t-1 at time (t - 1), parametrized by rank-R factors {An,t-1 ∈ RI t-1 n ×R}Nn=1 , i.e., Y t-1 = �A1,t-1 , · · · ,AN,t-1 �. Given the mask Ωt-1 , it satisfies

Following the evolution patterns, given the new underlying tensor X t ∈ RI t1×···×I tN and the mask Ωt at time t, our target is to find N new factor matrices An,t ∈ RI tn×R, I tn ≥ I t-1 n , ∀n and define Y t = �A1,t , · · · ,AN,t �, so as to minimize,

where l(·) is an element-wise loss function and the data and parameters are separated by semicolon.

2.2 Regularization and Approximation

This section expands the objective defined in Eqn. (2) with two regularization objectives.

• Regularization on Individual Factors. We add a generic term Lreg({An,t }) to Eqn. (2), and it can be later customized based on application background.
• Regularization on the Reconstruction. Usually, from time (t - 1) to t, the historical part of the tensor will not change much, and thus we assume that the new reconstruction Y t , restricted on the bound of previous tensor, will not change significantly from the previous Y t-1 .

Considering these two regularizations, the generalized objective can be written as,

Approximation. In Eqn. (2), the current tensor data masked by Ωt consists of two parts: (i) old unchanged data (indicating dark elements in Fig. 1), we denote it by Ωt,old , which is a subset of Ωt-1 ; (ii) newly added data (blue part, green and yellow entries in Fig. 1), denoted by set subtraction ˜Ωt = Ωt \ Ωt,old .

The old unchanged data can be in large size. Sometime, this part of data may not be entirely preserved due to (i) limited memory footprint or (ii) privacy related issues. By replacing the old tensor with its reconstruction [Song et al., 2017], we can avoid the access to the old data. Thus, we consider an approximation for the first term of Eqn. (3),

In the above derivation, the rationale of the approximation is the result in Eqn. (1). In Eqn. (4), we find that the term ∑Ωt,old l(yt-1 i1···iN ; yt i1···iN ) is part of the reconstruction regularization Lrec(Y t-1 ; Y t ) and thus can be absorbed. Therefore, we can use L˜Ωt(X t ; Y t ) to replace the full quantity LΩt(X t ; Y t ) in Eqn. (3), which results in a more efficient objective for streaming data processing,

For this new objective, the unchanged part Ωt,old is not counted in the first term (only the new elements ˜Ωt counted), however, it is captured implicitly in the second term.

In sum, L and LE are two objectives in our framework. They have a similar expression but with different access to the tensor data (i.e., the former with mask Ωt and the latter with ˜Ωt ). Generally, L is more accurate while LE is more efficient and can be applied to more challenging scenarios.

2.3 Generalized Optimization Algorithm

Structure of Parameters. For two general objectives defined in Eqn. (3) and (5), the parameters {An,t ∈ RI tn×R} are constructed by the upper blocks {Un,t ∈ RI t-1 n ×R} and the lower blocks {Ln,t ∈ R(I tn-I t-1 n )×R}, i.e.,

where Un,t corresponds to the old dimensions along tensor mode-n and Ln,t is for the newly added dimensions of moden due to mode growth. To reduce clutter, we use "factors" to refer to An,t , and "blocks" for Un,t , Ln,t .

Parameter Initialization. We use the previous factors at time (t - 1) to initialize the upper blocks, i.e.,

Note that, we use hat ˆ to denote the initialized block/factor.

Each lower block is initialized by solving the following objective, individually, ∀ n,

Here, Y t is constructed from the parameters {An,t }, and ik is the row index of the mode-k factor Ak,t . The summation is over the indices bounded by the intersection of Ωt and an Ndim cube, where other N - 1 modes use the old dimensions and mode-n uses the new dimensions. Thus, this objective essentially uses { ˆUk,t }k̸=n to initialize Ln,t .

ParameterUpdating. Generally, for any differentiable loss l(·), e.g, Frobenius norm [Yang et al., 2021] or KL divergence [Hong et al., 2020], we can apply gradient based methods, to update the factor matrices. The choices of loss function, regularization term, optimization method can be customized based on the applications.

3 GOCPT Framework

3.1 GOCPT Objectives

To instantiate a specific instance from the general algorithm formulation, we present GOCPT with

• Squared residual loss: l(x; y) = (x- y)2 ;
• L2 regularization: Lreg({An,t }) = ∑Nn=1 ∥An,t ∥2F ;
• Alternating least squares optimizer (ALS): updating factor matrices in a sequence by fixing other factors.

Then, the general objectives in Eqn. (3) (5) becomes

while the form of L is identical but with a different mask Ωt . Here, LE has only access to the new data {xti1···iN }˜Ωt but L has full access to the entire tensor {xti1···iN }Ωt up to time t.

In this objective, the second term (regularization on the reconstruction defined in Sec. 2.2) is restricted on {(i1, · · · , iN) : 1 ≤ in ≤ I t-1 n , ∀n}, so only the upper blocks Un,t ∈ RI t-1 n ×R of factors An,t ∈ RI tn×R, ∀n, are involved. Note that, L and LE are optimized following the same procedures.

3.2 GOCPT Optimization Algorithm

For our objectives in Eqn. (8), we outline the optimization flow: we first initialize the factor blocks; next, we update the upper or lower blocks (by fixing other blocks) following this order: U1,t ,L1,t , · · · ,UN,t ,LN,t .

Factor Initialization. As mentioned before, the upper blocks, {Un,t }Nn=1 , are initialized in Eqn. (6). In this specific framework, the minimization problem for initializing lower blocks {Ln,t }Nn=1 (in Eqn. (7)) can be reformed as

where Ωn,t denotes the intersection space in Eqn. (7). The initialization problem in Eqn. (9) and the targeted objectives in Eqn. (8) can be consistently handled by the following.

Factor Updating Strategies. Our optimization strategy depends on the density of tensor mask, e.g., ˜Ωt in Eqn. (8). For an evolving tensor with sparse new updates (for example, covid disease count tensor, where the whole tensor is sparse and new disease data or value correction comes irregularly at random locations), we operate only on the new elements by sparse strategy, while for tensors with dense updates (for example, multiangle imagery tensors are collected real-time and update slice-by-slice while each slice may have some missing or distortion values), we first impute the tensor to a full tensor, then apply dense operations by our dense strategy. For example, we solve Eqn. (8) as follows:

• Sparse strategy. If the ˜Ωt (the index set of newly added data at time t) is sparse, then we extend the CP completion alternating least squares (CPC-ALS) [Karlsson et al., 2016] and solve for each row of the factor.

Let us focus on an,t in ∈ R1×R, which is the in-th row of factor An,t . To calculate its derivative, we define the intermediate variables,

Here, ˜Ωt,in (we slightly abuse the notation) indicates the in-th slice of mask ˜Ωt along the n-th mode, and

The same convention works for �∗ . Here, Pn,t in ∈ RR×R is a positive definite matrix and qn,t in ∈ R1×R.

Then, the derivative w.r.t. each row can be expressed as,

Next, given that the objective is a quadratic function w.r.t. an,t in , we set the above derivative to zero and use the global minimizer to update each row,

Note that, this row-wise updating can apply in parallel. If ˜Ωt,in is empty, then we do not update an,t in .

• Dense strategy. If ˜Ωt is dense, then we extend the EMALS [Acar et al., 2011] method, which applies standard ALS algorithm on the imputed tensor. To be more specific,

we first impute the full tensor by interpolating/estimating the unobserved elements,

where we use ˆ X t to denote the estimated full tensor, and { ˆAn,t }Nn=1 are the initialized factors. Then, the first term of the objective in Eqn. (8) is approximated by a quadratic formw.r.t. each factor/block,

To calculate the derivative, let us define

where ˆXtn is the mode-n unfolding of ˆ X t , and ( ˆXtn )U ∈ RI t-1 n ×∏n̸=k I k,t n is the first I t-1 n rows of ˆXtn , while ( ˆXtn )L ∈ R(I tn-I t-1 n )×∏n̸=k I k,t n is the the remaining (I tn - I t-1 n ) rows of ˆXtn . Here, Pn,t U ,Pn,t L ∈ RR×R are positive definite, Qn,t U ∈ RI t-1 n ×R and Qn,t L ∈ R(I tn-I t-1 n )×R.

We express the derivative of the upper and lower blocks by the intermediate variables defined above,

Here, the overall objective LE is a quadratic function w.r.t. Un,t or Ln,t . By setting the derivative to zero, we can obtain the global minimizer for updating,

To sum up, the optimization flow for LE is summarized in Algorithm 1. For optimizing L, we simply modify the algorithm by replacing the input {xti1···iN }˜Ωt with {xti1···iN }Ωt . We analyze the complexity of our GOCPT in appendix B.

Algorithm 1: Factor Updates for LE at time t

• 1 Input: {An,t-1 }Nn=1, {xti1···iN }˜Ωt , α, β;;
• 2 Parameters: {Un,t }Nn=1, {Ln,t }Nn=1;;
• 3 Initialize parameters by Eqn. (6) and (9);;
• 4 for n ∈ [1, · · · ,N] do
• 5 update Un,t using Eqn. (10) or (11);;
• 6 update Ln,t using Eqn. (10) or (11);;
• 7 end

• 8 Output: new factors {An,t = [Un,t Ln,t ]}N n=1 .

4 Unifying Previous Online Settings

Several popular online tensor factorization/completion special cases can be unified in our framework. Among them, we focus on the online tensor factorization and online tensor completion. We show that those objectives in literature can be achieved by GOCPT. Our experiments confirm GOCPT can obtain comparable or better performance over previous methods in those special cases.

4.1 Online Tensor Factorization

Online tensor factorization tries to maintain the factorization results (e.g., CP factors) while the tensor is growing, shown in Fig. 2. Zhou et al. [Zhou et al., 2016] proposed an accelerated algorithm for the common setting where only the temporal mode grows. [Song et al., 2017] considered the multi-aspect setting where multiple mode grows simultaneously. We discuss the multi-aspect setting below.

Problem 2 (Multi-aspect Online Tensor Factorization). Suppose tensor X t-1 ∈ RI t-1 1 ×I t-1 2 ×I t-1 3 at time (t - 1) admits a low-rank approximation, induced by Ut-1 ∈ RI t-1 1 ×R, Vt-1 ∈ RI t-1 2 ×R,Wt-1 ∈ RI t-1 3 ×R, such that

At time t, we want to learn a new set of factors, Ut ∈ RI t1×R,Vt ∈ RI t2×R,Wt ∈ RI t3×R to approximate the growing new tensor X t ∈ RI t1×I t2×I t3 , which satisfies,

Unification. To achieve the existing objective in the literature, we simply reduce our LE by changing the L2 regularization into nuclear norm (i.e., the sum of singular values),

where γn, ∀ n ∈ [1, 2, 3] are the hyperparameters and they sum up to one. This regularization is the convex relaxation of CP rank [Song et al., 2017]. Then, the objective becomes

where Ωt and Ωt-1 are the bounds of the new and previous tensors and ˜Ωt = Ωt \ Ωt-1 (in this case, Ωt,old = Ωt-1 ) indicates the newly added elements at time t. This form is identical to the objective proposed in [Song et al., 2017].

In experiment Sec. 5.3, we evaluate on the temporal mode growth setting. Our methods use the L2 regularization as listed in Sec. 3 and follow Algorithm 1.

4.2 Online Tensor Completion

Online tensor completion studies incomplete tensors, where new incomplete slices will add to the temporal mode. The goal is to effectively compute a new tensor completion result (i.e., CP factors) for the augmented tensor. This problem is studied in [Mardani et al., 2015; Kasai, 2016; Kasai, 2019]. We show the formulation below.

Problem 3 (Online Tensor Completion). Suppose X t-1 is the underlying tensor at time (t-1), which admits a low-rank approximation, Ut-1 ∈ RI1×R, Vt-1 ∈ RI2×R, Wt-1 ∈ RI t-1 3 ×R, such that given the mask Ωt-1 ,

At the time t, given a new slice with incomplete entries, i.e., Ωt∆ �∗ Xt∆ ∈ RI1×I2 , the new data is concatenated along the temporal mode (the third-mode) by Ωt = [Ωt-1 ; Ωt∆] and X t = [X t-1 ;Xt∆]. We want to update the approximation with Ut ∈ RI1×R, Vt ∈ RI2×R, Wt ∈ RI t3×R, where I t3 = I t-1 3 +1, such that

Unification. To handle this setting, we remove the reconstruction regularizer (second term) from our L and the reduced version can be transformed into,

where Ωt,k and Xt,k indicate the k-th slices of the tensor along the temporal mode, and wtk is the k-th row of Wt . From Eqn. (12), we could make α exponential decaying w.r.t. time steps (i.e., αk = γ I t3-k ), the β within summation be time-variant (i.e., βk = λk × γ I t3-k ) and the β outside be λI t3 , where γ, {λk} are hyperparameters. We then obtain the identical objective as proposed in [Mardani et al., 2015]. The experiments are shown in Sec. 5.4.

5 Experiments

In the experiment, we evaluate our model on various settings. Our models are named GOCPT (with the objective L) and GOCPTE (with objective LE). Dataset statistics are listed in Table 1. More details can be found in appendix C.

5.1 Experimental Setups

Metrics. The main metric is percentage of fitness (PoF) [Acar et al., 2011], which is defined for the factorization or completion problem respectively by (the higher, the better)

where Ω and X are the mask and underlying tensor, {A1, · · · ,AN} are the low-rank CP factors. We also report the total time consumption as an efficiency indicator.

Baselines. We simulate three different practical scenarios for performance comparison. Since not all existing methods can deal with the three cases, we select representative stateof-the-art algorithms in each scenario for comparison:

• For the general case in Sec. 5.2, most previous models cannot support this setting. We adopt EM-ALS [Acar et al., 2011; Walczak and Massart, 2001] and CPC-ALS [Karlsson et al., 2016] as the compared methods, which follow the similar initialization procedure in Sec. 3.2.
• For the online tensor factorization in Sec. 5.3, we employ OnlineCPD [Zhou et al., 2016]; MAST [Song et al., 2017] and CPStream [Smith et al., 2018], which use ADMM and require multiple iterations; RLST [Nion and Sidiropoulos, 2009], which is designed only for third-order tensors.
• For the online tensor completion in Sec. 5.4, we implement EM-ALS and its variant, called EM-ALS (decay), which assigns exponential decaying weights for historical slices; OnlineSGD [Mardani et al., 2015]; OLSTEC [Kasai, 2016; Kasai, 2019] for comparison.

We compare the space and time complexity of each model in appendix B. All experiments are conduct with 5 random seeds. The mean and standard deviations are reported.

5.2 General Case with Three Evolving Patterns

Datasets and Settings. We use (i) JHU Covid data [Dong et al., 2020] and (ii) proprietary Patient Claim data to conduct the evaluation. The JHU Covid data was collected from Apr. 6, 2020, to Oct. 31, 2020 and the Patient Claim dataset collected weekly data from 2018 to 2019. To mimic tensor value updates on two datasets, we later add random perturbation to randomly selected 2% existing data with value changes uniformly of [-5%, 5%] at each time step. The leading 50% slices are used as preparation data to obtain the initial factors with rank R = 5. The results are in Table 2.

Result Analysis. Overall, our models show the top fitness performance compared to the baselines, and the variant GOCPTE shows 20% efficiency improvement over the best model with comparable fitness on both datasets. CPC-ALS and EM-ALS performs similarly and they are inferior to our models in both fitness and efficiency.

5.3 Special Case: Online Tensor Factorization

Dataset andSetups. We present the evaluation on low-rank synthetic data. In particular, we generate three CP factors from uniform [0, 1] distribution and then construct a low-rank tensor (I1, I2, I3,R) = (50, 50, 500, 5). We use the leading 10% slices along the (third) temporal mode as preparation; then, we add one slice at each time step to simulate mode growth. We report the mean values in Figure 4.

Also, we show that our GOCPT can provide better fitness than all baselines and GOCPTE outputs comparable fitness with state of the art efficiency on two real-world datasets: (i) ORL Database of Faces (FACE-3D) and (ii) Google Covid Search Symptom data (GCSS). The result tables are moved to appendix C.3 due to space limitation.

5.4 Special Case: Online Tensor Completion

Datasets and Setups. Using the same synthetic data described in Sec. 5.3, we randomly mask out 98% of the entries and follow the same data preparation and mode growth settings. The results of mean curves are shown in Fig. 4.

We also evaluate on two real-world datasets: (i) Indian Pines hyperspectral image dataset and (ii) a proprietary Covid disease counts data: location by disease by date, we call it the health tensor (CovidHT), and the results (refer to appendix C.4) show that GOCPT has the better fitness and GOCPTE has decent fitness with good efficiency.

6 Conclusion

This paper proposes a generalized online tensor factorization and completion framework, called GOCPT, which can support various tensor evolving patterns and unifies existing online tensor formulations. Experiments confirmed that our GOCPT can show decent performance on the general setting, where most previous methods cannot support. Also, our GOCPT provides comparable or better performance in each special setting. The GOCPT package is currently open-sourced.

Acknowledgements

This work was in part supported by the National Science Foundation award SCH-2014438, PPoSS 2028839, IIS1838042, the National Institute of Health award NIH R01 1R01NS107291-01 and OSF Healthcare. We thank Navjot Singh and Prof. Edgar Solomonik for valuable discussions.