I Introduction
The concept of multiway data analysis was introduced by Tucker in 1964 as an extension of standard twoway data analysis to analyze multidimensional data known as tensor. It is often used when traditional twoway data analysis methods such as Nonnegative Matrix Factorization (NMF), Principal Component Analysis (PCA) and Singular Value Decomposition (SVD) are not enough to capture the underlying structures inherited in multiway data. In the realm of multiway data, tensor decomposition methods such as
and (CP) [14] have been extensively studied and applied in various fields including signal processing, civil engineer and time series analysis. The CP decomposition has gained much popularity for analyzing multiway data due to its ease of interpretation. For example, given a tensor , CP method decomposes by loading matrices each represents one mode explicitly, where is the tensor order and each matrix represents one mode explicitly. In contrast to method where the three modes can interact with each other making it difficult to interpret the resultant matrices.The CP decomposition often uses the Alternating Least Squares (ALS) method to find the solution for a given tensor. The ALS method follows the batch mode training process which iteratively solves each component matrix by fixing all the other components, then it repeats the procedure until it converges [12]. However, ALS can lead to sensitive solutions [8][3]. Moreover, in the realm of big and nonstationary data, the ALS method raises many challenges in dealing with data that is continuously measured at high velocity from different sources/locations and dynamically changing over time. For example, structural health monitoring (SHM) data can be represented in a threeway form as which represents a large number of vibration responses measured over time by many sensors attached to a structure at different locations. This type of data can be found in many other application domains [1, 26, 15, 4]. In such applications, a naive approach would be to recompute the CP decomposition from scratch for each new incoming . Therefore, this would become impractical and computationally expensive as new incoming datum would have a minimal effect on the current tensor.
Zhou et al. [31] proposed a method called onlineCP to address the problem of online CP decomposition using ALS algorithm. The method was able to incrementally update the temporal mode in multiway data but failed for nontemporal modes [12]. In recent years, several studies have been proposed to solve the CP decomposition using stochastic gradient descent (SGD) algorithm which has the capability to deal with big data and online learning. However, such methods are inefficient and impractical due to slow convergence, numerical uncertainty and nonconvergence [9, 5, 18, 23]. To address the aforementioned problems, we propose an efficient solver method called NeCPD (Nesterov CP Decomposition) for analyzing largescale highorder data. The novelty of our proposed method is summarized in the following contributions:

Online CP Decomposition. Our method is capable to update in one single step. We realize this by employing the Stochastic Gradient Descent (SGD) algorithm to solve the CP decomposition as a baseline method.

Global convergence guarantee. We followed the perturbation approach which adds a little noise to the gradient update step to reinforce the next update step to start moving away from a saddle point toward the correct direction.

Optimal convergence rate. Our method employs Nesterov’s Accelerated Gradient (NAG) method into SGD algorithm to optimally accelerate the convergence rate [27]. It achieves a global convergence rate of comparing to for traditional SGD.

Empirical analysis on structural datasets. We conduct experimental analysis using laboratorybased and reallife datasets in the field of SHM. The experimental analysis shows that our method can achieve more stable and fast tensor decomposition compared to other known existing online and offline methods.
The remainder of this paper is organized as follows. Section II introduces background knowledge and review of the related work. Section III describes our novel NeCPD algorithm for CP decomposition based on SGD algorithm augmented with NAG method and perturbation approach. Section IV presents the motivation of this work. Section V shows the performance of NeCPD on structural datasets and presents our experimental results on both laboratorybased and reallife datasets. The conclusion and discussion of future research work are presented in Section VI.
Ii Background and Related work
Iia CP Decomposition
Given a threeway tensor , CP decomposes into three matrices , and , where is the latent factors. It can be written as follows:
(1) 
where ”
” is a vector outer product.
is the latent element, and are rth columns of component matrices , and . The main goal of CP decomposition is to decrease the sum square error between the model and a given tensor . Equation 2shows our loss function
needs to be optimized.(2) 
where is the sum squares of and the subscript is the Frobenius norm. The loss function presented in Equation 2 is a nonconvex problem with many local minima since it aims to optimize the sum squares of three matrices. Several algorithms have been proposed to solve CP decomposition [28, 17, 22]. Among these algorithms, ALS has been heavily employed which repeatedly solves each component matrix by locking all other components until it converges [21]. The rational idea of the least square algorithm is to set the partial derivative of the loss function to zero with respect to the parameter we need to minimize. Algorithm 1 presents the detailed steps of ALS.
IiB Stochastic Gradient Descent
Stochastic gradient descent algorithm is a key tool for optimization problems. Lets say our aim is to optimize a loss function , where is a data point drawn from a distribution and is a variable. The stochastic optimization problem can be defined as follows:
(3) 
The stochastic gradient descent method solves the above problem defined in Equation 3 by repeatedly updates to minimize . It starts with some initial value of and then repeatedly performs the update as follows:
(4) 
where is the learning rate and is a random sample drawn from the given distribution . This method guarantees the convergence of the loss function to the global minimum when it is convex. However, it can be susceptible to many local minima and saddle points when the loss function exists in a nonconvex setting. Thus it becomes an NPhard problem. In fact, the main bottleneck here is due to the existence of many saddle points and not to the local minima [9]. This is because the rational idea of gradient algorithm depends only on the gradient information which may have even though it is not at a minimum.
Previous studies have used SGD for matrix factorization and tensor decomposition with extensions to handle the aforementioned issues. Naiyang et al.[11] applies Nesterov’s optimal gradient method to SGD for nonnegative matrix factorization. This method accelerates the NMF process with less computational time. Similarly, Shuxin et al.[30] used an SGD algorithm for matrix factorization using Taylor expansion and Hessian information. They proposed a new asynchronous SGD algorithm to compensate for the delay resultant from a Hessian computation.
Recently, SGD has attracted several researchers working on tensor decomposition. For instance, Ge et al.[9] proposed a perturbed SGD (PSGD) algorithm for orthogonal tensor optimization. They presented several theoretical analysis that ensures convergence; however, the method is not applicable to nonorthogonal tensor. They also didn’t address the problem of slow convergence. Similarly, Maehara et al.[18] propose a new algorithm for CP decomposition based on a combination of SGD and ALS methods (SALS). The authors claimed the algorithm works well in terms of accuracy. Yet its theoretical properties have not been completely proven and the saddle point problem was not addressed. Rendle and Thieme [23] propose a pairwise interaction tensor factorization method based on Bayesian personalized rank. The algorithm was designed to work only on threeway tensor data. To the best of our knowledge, this is the first work applies SGD algorithm augmented with Nesterov’s optimal gradient and perturbation methods for optimal CP decomposition of multiway tensor data.
Iii Nesterov CP Decomposition (NeCPD)
Given an order tensor , NeCPD initially divides the CP problem into a convex subproblems since its loss function is nonconvex problem which may have many local minima. For simplicity, we present our method based on threeway tensor data. However, it can be naturally extended to handle a higherorder tensor.
In a threeway tensor , can be decomposed into three matrices , and , where is the latent factors. Following the SGD, we need to calculate the partial derivative of the loss function defined in Equation 2 with respect to the three modes and alternatively as follows:
(5)  
where is an unfolding matrix of tensor in mode . The gradient update step for and are as follows:
(6)  
The rational idea of SGD algorithm depends only on the gradient information of . In such nonconvex setting, this partial derivative may encounter data points with even though it is not at a global minimum. These data points are known as saddle points which may detente the optimization process to reach the desired local minimum if not escaped [9]. These saddle points can be identified by studying the secondorder derivative (aka Hessian) . Theoretically, when the , must be a local minimum; if , then we are at a local maximum; if
has both positive and negative eigenvalues, the point is a saddle point. The second order methods guarantee convergence, but the computing of Hessian matrix
is high, which makes the method infeasible for high dimensional data and online learning. Ge
et al.[9] show that saddle points are very unstable and can be escaped if we slightly perturb them with some noise. Based on this, we use the perturbation approach which adds Gaussian noise to the gradient. This reinforces the next update step to start moving away from that saddle point toward the correct direction. After a random perturbation, it is highly unlikely that the point remains in the same band and hence it can be efficiently escaped (i.e., no longer a saddle point). Since we interested in fast optimization process due to online settings, we further incorporate Nesterov’s method into the perturbedSGD algorithm to accelerate the convergence rate. Recently, Nesterov’s Accelerated Gradient (NAG) [19] has received much attention for solving convex optimization problems [11, 20, 10]. It introduces a smart variation of momentum that works slightly better than standard momentum. This technique modifies the traditional SGD by introducing velocity and friction , which tries to control the velocity and prevents overshooting the valley while allowing faster descent. Our idea behind Nesterov’s is to calculate the gradient at a position that we know our momentum is about to take us instead of calculating the gradient at the current position. In practice, it performs a simple step of gradient descent to go from to , and then it shifts slightly further than in the direction given by . In this setting, we model the gradient update step with NAG as follows:(7) 
where
(8) 
where is a Gaussian noise, is the step size, and is the regularization and penalization parameter into the norms to achieve smooth representations of the outcome and thus bypassing the perturbation surrounding the local minimum problem. The updates for and are similar to the aforementioned ones. With NAG, our method achieves a global convergence rate of comparing to for traditional gradient descent. Based on the above models, we present our NeCPD algorithm 2.
Iv Motivation
Numerous types of data are naturally structured as multiway data. For instance, structural health monitoring (SHM) data can be represented in a threeway form as . Arranging and analyzing the SHM data in multidimensional form would allow to capture the correlation between sensors at different locations and at the same time which was not possible using the standard twoway matrix
. Furthermore, in SHM only positive data instances i.e healthy state are available. Thus, the problem becomes an anomaly detection problem in higherorder datasets. Rytter
[24]affirms that damage identification requires also damage localization and severity assessment which are considered much more complex than damage detection since they require a supervised learning approach
[29].Given a positive threeway SHM data , NeCPD decomposes into three matrices and . The matrix represents the temporal mode where each row contains information about the vibration responses related to an event at time . The analysis of this component matrix can help to detect the damage of the monitored structure. Therefore, we use the matrix to build a oneclass anomaly detection model using only the positive training events. For each new incoming , we update the three matrices and incrementally as described in Algorithm 2
. Then the constructed model estimates the agreement between the new event
and the trained data.For damage localization, we analyze the data in the location matrix , where each row captures meaningful information for each sensor location. When the matrix is updated due to the arrival of a new event , we study the variation of the values in each row of matrix by computing the average distance from ’s row to nearest neighboring locations as an anomaly score for damage localization. For severity assessment in damage identification, we study the decision values returned from the oneclass model. This is because a structure with more severe damage will behave much differently from a normal one.
V Experimental Results
We conduct all our experiments using an Intel(R) Core(TM) i7 CPU 3.60GHz with 16GB memory. We use R language to implement our algorithms with the help of the two packages rTensor and e1071 for tensor tools and oneclass model.
Va Comparison on synthetic data
Our initial experiment was to compare our NeCPD method to SGD, PSGD, and SALS algorithms in terms of robustness and convergence. We generated a synthetic long time dimension tensor from 12 random loading matrices
in which entries were drawn from uniform distribution
. We evaluated the performance of each method by plotting the number of sample versus the root mean square error (RMSE). For all experiments we use the learning rate . It can be clearly seen from Figure 1 that NeCPD outperformed the SGD and PSGD algorithms in terms of convergence and robustnesses. Both SALS and NeCPD converged to a small RMSE but it was lower and faster in NeCPD.VB Comparison on SHM data
We evaluate our NeCPD algorithm on laboratorybased and reallife structural datasets in the area of SHM. The two datasets are naturally in multiway form. The first dataset is collected from a cablestayed bridge in Western Sydney, Australia. The second one is a specimen building structure obtained from Los Alamos National Laboratory (LANL) [16].
VB1 The CableStayed Bridge:
We used 24 uniaxial accelerometers and 28 strain gauges to measure the vibration and strain responses of the bridge at different locations. However, only accelerations data collected from sensors with were used in this study. Figure 2 illustrates the locations of the 24 sensors on the bridge deck. Since the bridge is in healthy condition and in order to evaluate the performance of our method in damage identification tasks, we placed two different stationary vehicles with different masses to emulate two different levels of damage severity [13, 7]. Three different categories of data were collected in this experiment: ”HealthyData” when the bridge is free of vehicles; ”CarDamage” when a light car vehicle is placed on the bridge close to location ; and ”BusDamage” when a heavy bus vehicle is located on the bridge at location . This experiment generates 262 samples (i.e., events) separated into three categories: ”HealthyData” (125 samples), ”CarDamage” data (107 samples) and ”BusDamage” data(30 samples). Each event consists of acceleration data for a period of 2 seconds sampled at a rate of 600Hz. The resultant event’s feature vector composed of 1200 frequency values.
VB2 Building Data:
This data is based on experiments conducted by LANL [16] using a specimen for a threestory building structure as shown in Figure 3. Each joint in the building was instrumented by two accelerometers. The excitation data was generated using a shaker placed at corner . Similarly, for the sake of damage detection evaluation, the damage was simulated by detaching or loosening the bolts at the joints to induce the aluminum floor plate moving freely relative to the Unistrut column. Three different categories of data were collected in this experiment: ”HealthyData” when all the bolts were firmly tightened; ”Damage3C” data when the bolt at location 3C was loosened; and ”Damage1A3C” data when the bolts at locations 1A and 3C were loosened simultaneously. This experiment generates 240 samples (i.e., events) which also were separated into three categories: HealthyData (150 samples), ”Damage3C” data (60 samples) and ”Damage1A3C” data(30 samples). The acceleration data was sampled at 1600 Hz. Each event was measured for a period of 5.12 seconds resulting in a vector of 8192 frequency values.
VB3 Feature Extraction:
The raw signals of the sensing data we collected in the aforementioned experiments exist in the time domain. In practice, time domainbased features may not capture the physical meaning of the structure. Thus, it is noteworthy to convert the generated data to a frequency domain. For all datasets, we initially normalize the timedomain features to have zero mean and one standard deviation. Then we use the fast Fourier transform method to convert them into the frequency domain. The resultant threeway data collected from the CableStayed Bridge now has a structure of 600 features
24 sensors 262 events. For Building data, we compute the difference between signals of two adjacent sensors which leads to 12 different joints in the three stories as in [16]. Then we select the first 150 frequencies as a feature vector resulting in a threeway data with a structure of 768 features 12 locations 240 events.VB4 Experimental Setup:
For all case studies we apply the following procedures:

The bootstrap technique to randomly select 80% of the healthy samples for training and the remaining 20% for testing in addition to the damage samples. All the reported accuracies are based on the average results over ten trials of the bootstrap experiment.

The core consistency diagnostic (CORCONDIA) technique described in [6] is used to decide the number of rankone tensors in the NeCPD.

We use oneclass support vector machine (OSVM)
[25] as our one class model for anomaly detection method. The Gaussian kernel parameter in OCSVM is tuned using the Edged Support Vector (ESV) algorithm [2], and the rate of anomalies was set to 0.05. 
The Fscore measure to compute the accuracy values. It is defined as where and (the number of true positive, false positive and false negative are abbreviated by TP, FP and FN, respectively).

The results of the competitive method SALS proposed in [18] are compared against our NeCPD method.
VB5 Results and Discussions:
VB6 The CableStayed Bridge:
Our NeCPD method with oneclass SVM was initially validated using vibration data collected from the cablestayed bridge described in Section VB1. The healthy training threeway tensor data (i.e. training set) was in the form of . The 137 examples related to the two damage cases were added to the remaining 20% of the healthy data to form a testing set, which was later used for the model evaluation. This experiment generates a damage detection accuracy Fscore of on the testing data. On the other hand, the Fscore accuracy of oneclass SVM using SALS is recorded at .
As we anticipated, tensor analysis with NeCPD is able to capture the underlying structure in multiway data with better convergence. This is demonstrated by plotting the decision values returned from oneclass SVM based NeCPD (as shown in Figure 4). We can clearly separate the two damage cases (”CarDamage” and ”BusDamage”) in this dataset where the decision values are further decreased for the samples related to the more severe damage cases (i.e., ”BusDamage”). These results suggest using the decision values obtained by NeCPD and oneclass SVM as structural health scores to identify the damage severity in a oneclass aspect. In contrast, the resultant decision values of oneclass SVM based on SALS are also able to track the progress of the damage severity in the structure but with a slight decreasing trend in decision values for ”BusDamage” (see Figure 4).
The last step was to analyze the location matrix obtained by NeCPD to locate the detected damage. Each row in this matrix captures meaningful information for each sensor location. Therefore, we calculate the average distance from each row in the matrix to nearest neighboring rows. Figure 5 shows the obtained nn score for each sensor. The first 25 events (depicted on the xaxis) represent healthy data, followed by 107 events related to ”CarDamage” and 30 events to ”BusDamage”. It can be clearly observed that NeCPD method successfully localized the damage in the structure. Whereas, sensors and related to the ”CarDamage” and ”BusDamage” respectively behave significantly different from all the other sensors apart from the position of the introduced damage. Also, we observe that the adjacent sensors to the damage location (e.g , , and ) react differently due to the arrival pattern of the damage events. The SALS method, however, is not able to successfully locate the damage since it fails to update the location matrix incrementally.
VB7 Building Data:
Our second experiment is conducted using the acceleration data acquired from 24 sensors instrumented on the threestory building as described in Section VB2. The healthy training threeway data (i.e.training set) was in the form of . The remaining 20% of the healthy data and the data obtained from the two damage cases were used for testing (i.e.testing set). The NeCPD with oneclass SVM achieve an Fscore of on the testing data compared to obtained from oneclass SVM based SALS.
Similar to the Bridge dataset, we further analyzed the resultant decision values which were also able to characterize damage severity. Figure 6 demonstrates that the more severe damage to the and location test data the more deviation from the training data with lower decision values.
The last experiment is to compute the nn score for each sensor based on the nearest neighboring of the average distance between each row of the matrix . Figure 7 shows the resultant nn score for each sensor. The first 30 events (depicted on the xaxis) represent the healthy data, followed by 60 events describing when the damage was introduced in location . The last 30 events represent the damage occurred in both locations and . It can be clearly observed that the NeCPD method successfully localized the structure’s damage where sensors and behave significantly different from all the other sensors apart from the position of the introduced damage. However, the SALS method is not able to successfully locate damage since it fails to update the location matrix incrementally.
Vi Conclusion
This paper investigated the CP decomposition with stochastic gradient descent algorithm for multiway data analysis. This leads to a new method named NeCPD for tensor decomposition. Our new method guarantees the convergence for a given nonconvex problem by modeling the second order derivative of the loss function and incorporating little noise to the gradient update. Furthermore, NeCPD employs Nesterov’s method to compensate for the delays of the optimization process and accelerate the convergence rate. Based on laboratory and real datasets from the area of SHM, our NeCPD, with oneclass SVM model for anomaly detection, achieve accurate results in damage detection, localization, and assessment in online and oneclass settings. Among the key future work is how to parallelize the tensor decomposition with NeCPD. Also, it would be useful to apply NeCPD with datasets from different domains.
Vii Acknowledgement
The authors wish to thank the Roads and Maritime Services (RMS) in New South Wales, New South Wales Government in Australia and Data61 (CSIRO) for provision of the support and testing facilities for this research work. Thanks are also extended to Western Sydney University for facilitating the experiments on the cablestayed bridge.
References
 [1] (2009) Unsupervised multiway data analysis: a literature survey. IEEE transactions on knowledge and data engineering 21 (1), pp. 6–20. Cited by: §I.

[2]
(2018)
Gaussian kernel parameter optimization in oneclass support vector machines.
In
2018 International Joint Conference on Neural Networks (IJCNN)
, pp. 1–8. Cited by: 3rd item.  [3] (2018) Regularized tensor learning with adaptive oneclass support vector machines. In International Conference on Neural Information Processing, pp. 612–624. Cited by: §I.
 [4] (2018) A tensorbased structural damage identification and severity assessment. Sensors 18 (1), pp. 111. Cited by: §I, Fig. 2.
 [5] (2016) Efficient approaches for escaping higher order saddle points in nonconvex optimization. In Conference on learning theory, pp. 81–102. Cited by: §I.
 [6] (2003) A new efficient method for determining the number of components in parafac models. Journal of chemometrics 17 (5), pp. 274–286. Cited by: 2nd item.
 [7] (2012) Indirect structural health monitoring in bridges: scale experiments. In Proc. Int. Conf. Bridge Maint., Safety Manag., Lago di Como, pp. 346–353. Cited by: §VB1.
 [8] (1980) Perturbation theory for the least squares problem with linear equality constraints. SIAM Journal on Numerical Analysis 17 (3), pp. 338–350. Cited by: §I.
 [9] (2015) Escaping from saddle points—online stochastic gradient for tensor decomposition. In Conference on Learning Theory, pp. 797–842. Cited by: §I, §IIB, §IIB, §III.
 [10] (2016) Accelerated gradient methods for nonconvex nonlinear and stochastic programming. Mathematical Programming 156 (12), pp. 59–99. Cited by: §III.
 [11] (2012) NeNMF: an optimal gradient method for nonnegative matrix factorization. IEEE Transactions on Signal Processing 60 (6), pp. 2882–2898. Cited by: §IIB, §III.
 [12] (2017) Smart infrastructure maintenance using incremental tensor analysis. In Proceedings of the 2017 ACM on Conference on Information and Knowledge Management, pp. 959–967. Cited by: §I, §I.
 [13] (2013) Identification of physically simulated damage on a footbridge based on ambient vibration data. In Structures Congress 2013: Bridging Your Passion with Your Profession, pp. 352–362. Cited by: §VB1.
 [14] (2009) Tensor decompositions and applications. SIAM review 51 (3), pp. 455–500. Cited by: §I.
 [15] (2008) Scalable tensor decompositions for multiaspect data mining. In 2008 Eighth IEEE international conference on data mining, pp. 363–372. Cited by: §I.
 [16] (1987) Los alamos national laboratory report no. Technical report LAUR86748. Cited by: Fig. 3, §VB2, §VB3, §VB.

[17]
(2014)
Speedingup convolutional neural networks using finetuned cpdecomposition
. arXiv preprint arXiv:1412.6553. Cited by: §IIA. 
[18]
(2016)
Expected tensor decomposition with stochastic gradient descent.
In
Thirtieth AAAI conference on artificial intelligence
, Cited by: §I, §IIB, 5th item.  [19] (2013) Introductory lectures on convex optimization: a basic course. Vol. 87, Springer Science & Business Media. Cited by: §III.
 [20] (2014) Stochastic proximal gradient descent with acceleration techniques. In Advances in Neural Information Processing Systems, pp. 1574–1582. Cited by: §III.
 [21] (2017) Tensors for data mining and data fusion: models, applications, and scalable algorithms. ACM Transactions on Intelligent Systems and Technology (TIST) 8 (2), pp. 16. Cited by: §IIA.
 [22] (2009) Learning optimal ranking with tensor factorization for tag recommendation. In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 727–736. Cited by: §IIA.
 [23] (2010) Pairwise interaction tensor factorization for personalized tag recommendation. In Proceedings of the third ACM international conference on Web search and data mining, pp. 81–90. Cited by: §I, §IIB.
 [24] (1993) Vibrational based inspection of civil engineering structures. Ph.D. Thesis, Dept. of Building Technology and Structural Engineering, Aalborg University. Cited by: §IV.

[25]
(2000)
Support vector method for novelty detection
. In Advances in neural information processing systems, pp. 582–588. Cited by: 3rd item.  [26] (2008) Incremental tensor analysis: theory and applications. ACM Transactions on Knowledge Discovery from Data (TKDD) 2 (3), pp. 11. Cited by: §I.

[27]
(2013)
On the importance of initialization and momentum in deep learning.
. ICML (3) 28 (11391147), pp. 5. Cited by: item 3.  [28] (2008) Tag recommendations based on tensor dimensionality reduction. In Proceedings of the 2008 ACM conference on Recommender systems, pp. 43–50. Cited by: §IIA.

[29]
(2006)
The application of machine learning to structural health monitoring
. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365 (1851), pp. 515–537. Cited by: §IV.  [30] (2017) Asynchronous stochastic gradient descent with delay compensation. In Proceedings of the 34th International Conference on Machine LearningVolume 70, pp. 4120–4129. Cited by: §IIB.
 [31] (2016) Accelerating online cp decompositions for higher order tensors. In Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pp. 1375–1384. Cited by: §I.