Personalized Tucker Decomposition: Modeling Commonality and Peculiarity on Tensor Data

A tensor can be regarded as a data structure with more than 2 dimensions, also known as modes in tensor analysis (see Fig. 1). For example, in images, we use a vector of length $3$ to represent the RGB channel of the pixel. Thus, a picture can be represented by a 3-dimensional tensor with dimensions height$\times$width$\times$RGB. If we have multiple pictures of the same dimension, a dataset can be represented by a 4-dimensional tensor.

Suppose we have tensor data from $N$ sources. Each source has tensor data of order $K$. We use ${\bm{\mathcal{Y}}}_{n}$ to denote the data from source $n$, and assume that ${\bm{\mathcal{Y}}}_{n}$ has dimensions $I_{1}\times I_{2}\times\ldots\times I_{K}\times s_{n}$. Here, all dimensions across sources have the same length except for the last one. In particular, in practical applications, the last dimension $s_{n}$ denotes the number of samples from the source $n$, which often differs between sources.

Our approach relies on defining global and local components to model commonality and heterogeneity across different sources. To do so, we let the global components consist of shared global factor matrices ${\bm{U}}_{G,1},\ldots,{\bm{U}}_{G,K}$ and individual global core tensors $\bm{\mathcal{C}}_{G,1},\ldots,\bm{\mathcal{C}}_{G,N}$ for each source. The local components consist of individual core tensors $\bm{\mathcal{C}}_{L,1},\ldots,\bm{\mathcal{C}}_{L,N}$ and individual local factor matrices $\bm{V}_{n,1},\ldots,\bm{V}_{n,K}$. As such, the reconstructions of the global and local components for source $n$ are $\bm{\mathcal{C}}_{G,n}\times_{1}{\bm{U}}_{G,1}\ldots\times_{K}{\bm{U}}_{G,K}$ and $\bm{\mathcal{C}}_{L,n}\times_{1}\bm{V}_{n,1}\ldots\times_{K}\bm{V}_{n,K}$, respectively.

Based on the above definitions, we assume our data-generating process to be

where $\bm{\mathcal{E}}_{n}$ are tensors that represent additive noise.

Since global and local components should convey different information, they need to be distinguished so that each part can vary independently of each other. To do so, we require the orthogonality of the global and local tensors. Specifically, we assume that:

Given Proposition 1, we require local factor matrices to be orthogonal to global factor matrices for all sources in at least one mode. We define the set of such orthogonal modes by $\mathcal{K}$, $|\mathcal{K}|\geq 1$. Then our objective is to minimize the reconstruction loss of the data across all $N$ sources. This is written as

We assume that the dimension of the global core tensor for all sources is ${\bm{\mathcal{C}}}_{G,n}\in\mathbb{R}^{g_{1}\times\cdots\times g_{K}}$, and the dimension of the local core tensor for source $n$ is ${\bm{\mathcal{C}}}_{L,n}\in\mathbb{R}^{l_{n,1}\times\cdots\times l_{n,K}}$. This also defines the dimension of the global and local factor matrices.

A natural algorithm to solve the objective in (2) is block coordinate descent (BCD), where we iteratively optimize each variable. A general framework for BCD in our context is outlined in Algorithm 1.

In the rest of Sec. 3.3, we explain the update steps in Algorithm 1 in detail. We start with the update of global and local core tensors in Sec. 3.3.1 since the closed-form solution is a direct projection similar to the traditional Tucker decomposition owing to the orthogonality between global and local components. Then the solution consistently holds in the update of global and local factor matrices we introduced in Sec. 3.3.2 and Sec. 3.3.2. This simplifies the update of the factor matrices without the core tensors. Despite the challenges that pertain to the two distinct decomposition components within (2) and their orthogonality, a key result is that all updates can be done in closed form owing to the nice properties of the orthogonality constraint imposed on the model.

In this section, we provide the the convergence analysis of Algorithm 2. The special update rule in Proposition 4 brings challenges to the convergence analysis. In Algorithm 2, the update of local factors ${\bm{V}}_{n,k}$’s is different from the standard Tucker decomposition update, as ${\bm{V}}_{n,k}$ is required to be orthogonal to ${\bm{U}}_{G,k}$. As a result, updating the local factors does not necessarily decrease the objective value in (2). Thus, Algorithm 2 is not a strictly descent algorithm. Despite such subtleties, we can show that, when the proximal parameter $\rho$ is not too small, our algorithm can converge into stationary solutions.

We will present our theorem on global convergence in the following theorem. Recall that we use ${\bm{U}}_{G,k,t}$ to denote the $k$-th global factor ${\bm{U}}_{G,k}$ after iteration $t$, and ${\bm{V}}_{n,k,t}$ to denote the $k$-th local factor of source $n$ after iteration $t$.

Theorem 5 provides many key insights. First, it shows that the subspaces spanned by the column vectors of global and local factors all converge into fixed solutions. The result establishes the global convergence of factors, as it does not require careful initialization. Second, the convergence rates for global and local factors differ. Global factors converge at a rate of $O(\frac{1}{T})$, which is standard in non-convex optimization. However, since some local factors must be perpendicular to the global factors, they converge at a slightly slower rate of $O(\frac{1}{\sqrt{T}})$. Third, our result is based on having $|{\mathcal{K}}|\geq 2$. This requirement allows orthogonality to be maintained by each mode being updated.

To validate the convergence rates, we provide a proof-of-concept simulation study. Fig. 2 displays an example of the convergence of global and local factor matrices in this simulation. Both subspace errors of the local and global components go to 0. This slower rate of local components is primarily a result of the orthogonality requirement, and this result verifies Theorem 5. The detail of this simulation study is relegated to Appendix H

One simple approach is to initialize all components randomly. Alternatively, one may use a Tucker decomposition on all the data for initialization. To do so, let $s=\sum_{n=1}^{N}s_{n}$ be the total number of samples from all sources and recall that the data ${\bm{\mathcal{Y}}}_{n}$ from source $n$ has dimension $I_{1}\times\ldots\times I_{K}\times s_{n}$. Now the following steps can be taken:

1. 1.

   Construct a tensor ${\bm{\mathcal{Y}}}$ that encompasses all samples from every source, with dimensions $I_{1}\times\ldots\times I_{K}\times s$.

2. 2.

   Use Tucker decomposition ${\bm{\mathcal{Y}}}\approx{\bm{\mathcal{C}}}\times_{1}{\bm{U}}_{G,1}\ldots\times_{K}{\bm{U}}_{G,K}$ to initialize global factors.

3. 3.

   Use Proposition 2 to initialize the global core tensor for each source.

4. 4.

   For each source $n$, perform the Tucker decomposition on the local residual tensor ${\bm{\mathcal{R}}}_{L,n}={\bm{\mathcal{Y}}}_{n}-{\bm{\mathcal{C}}}^{\star}_{G,n}\times_{1}{\bm{U}}_{G,1}\ldots\times_{K}{\bm{U}}_{G,K}$ to initialize the local core tensor and the local factor matrices.

This initialization does not require orthogonality between global and local components. Therefore, we may observe an increase in the reconstruction error in the first iteration. But we have found that in practice, this method yields faster convergence.

In addition, in our model, one needs to choose the modes in which orthogonality is imposed. At the same time, our theorem suggests that $|{\mathcal{K}}|\geq 2$, and we found that the result is also true for $|{\mathcal{K}}|=1$ in various simulation studies. In practice, one can simply choose the mode with the largest dimension to impose the orthogonality constraint. Alternatively, cross-validation can be utilized to select the mode.

In this section, we introduce some practical applications of perTucker. Specifically, we shed light on its potential utility for improved classification, anomaly detection, and clustering. The key idea for all applications is to operate only on local components. This may allow for improved clustering, classification, and detection as differences become more explicit when shared knowledge is removed.

In this simulation work, each sample of the data is a grayscale image with dimensions 50 by 50. The construction of each sample is low-rank global component, heterogeneous local component, and i.i.d standard normal noise, as in Eq. (13)

We generate $N=3$ clients defined as $3$ patterns for the heterogeneous local component: Swiss pattern, oval pattern, and rectangle pattern, as in Fig. 3. The value of all the patterns is 5 while the rest part is 0. In each pattern, we generate $10$ sample images. There is some variability within each pattern, as shown in the left part of Fig. 3. The Swiss can be thin or thick; the oval can be vertical, horizontal, or circular; the rectangle can be wide, tall, or square.

For the global component, we randomly create orthonormal matrices ${\bm{U}}_{G,1}$ and ${\bm{U}}_{G,2}$ with dimension $50\times 5$ for the 3 clients. Then we randomly generate the global core tensor $\bm{\mathcal{C}}_{G,n}$ with dimension $5\times 5\times 10$ for each client. And each entry of the global core tensors follows i.i.d. $N(0,100)$. The global components are constructed by $\bm{\mathcal{Y}}_{G,n}=\bm{\mathcal{C}}_{G,n}\times_{1}{\bm{U}}_{G,1}\times_{2}{\bm{U}}_{G,2}$, $n=1,2,3$.

Therefore, the full data dimension is $3\times 50\times 50\times 10$. Some examples of the data generation structure are shown in the right part of Fig. 3. The three rows are for three patterns. The two columns show the global and local components, respectively, and the third column shows the sum of the global and local components along with the error term. With noise and global background, the local pattern can barely be recognized, which makes accurate identification of the local patterns challenging.

In the generated data, we apply perTucker to decouple the global and local components. For comparison, we also evaluate the performance of some benchmark algorithms.

1. 1.

   globalTucker: We first concatenate the samples of all clients into tensor ${\bm{\mathcal{Y}}}$ and then apply the Tucker decomposition on ${\bm{\mathcal{Y}}}$.

2. 2.

   localTucker: We apply a standard Tucker decomposition on each client ${\bm{\mathcal{Y}}}_{n}$ individually.

3. 3.

   robustTucker: We first concatenate samples from all clients in ${\bm{\mathcal{Y}}}$, then apply the method in Lu et al., (2019) to identify the low-rank and sparse components.

4. 4.

   perPCA: We apply perPCA (Shi and Kontar,, 2022) on the vectorized dataset where we vectorize $50\times 50$ images into vectors of length $2500$ and use perPCA to find global and local components. Although perPCA is designed for vector datasets, this comparison can highlight the need for personalized Tensor decompositions when data is in tensor form.

Fig. 4 depicts examples of global and local component reconstruction via perTucker. The first row represents the data, and the second row shows the reconstruction from perTucker. The columns represent global and local components. The three examples with different patterns indicate that perTucker can effectively reconstruct the shared and unique data patterns.

Furthermore, to numerically examine the performance and benchmark algorithms, we calculate a few performance metrics. $\hat{{\bm{U}}}_{G,k}$ and $\hat{{\bm{V}}}_{n,k}$ represent the estimated global and local factor matrices; $\bm{\mathcal{\hat{Y}}}_{G,n}$ and $\bm{\mathcal{\hat{Y}}}_{L,n}$ represent the estimated reconstruction of global and local components.

1. 1.

   Global subspace error: We compute the regularized global subspace error between the ground truth global factors $\{{\bm{U}}_{G,k}\}$ and the estimated global factors $\{\hat{{\bm{U}}}_{G,k}\}$ by $\varrho(\bigotimes_{k=1}^{K}{\bm{U}}_{G,k},\bigotimes_{k=1}^{K}\bm{\hat{U}}_{G,k})/\|\bigotimes_{k=1}^{K}{\bm{U}}_{G,k}\|_{F}^{2}$.

2. 2.

   Local subspace error: We first generate 100 images for each pattern, and use Tucker decomposition to estimate the local factors $\{{\bm{V}}_{n,k}\}$ for each pattern. Then the regularized local subspace error is calculated by $\varrho(\bigotimes_{k=1}^{K}{\bm{V}}_{n,k},\bigotimes_{k=1}^{K}\hat{{\bm{V}}}_{n,k})/\|\bigotimes_{k=1}^{K}{\bm{V}}_{n,k}\|_{F}^{2}$, and take the average of $3$ patterns.

3. 3.

   Global component error: defined by $\small\sum_{n}\|\bm{\mathcal{\hat{Y}}}_{G,n}-\bm{\mathcal{Y}}_{G,n}\|_{F}^{2}/\sum_{n}\|\bm{\mathcal{Y}}_{G,n}\|_{F}^{2}$.

4. 4.

   Local component error: defined by $\small\sum_{n}\|\bm{\mathcal{\hat{Y}}}_{L,n}-\bm{\mathcal{Y}}_{L,n}\|_{F}^{2}/\sum_{n}\|\bm{\mathcal{Y}}_{L,n}\|_{F}^{2}$.

5. 5.

   Denoised error: defined by $\small\sum_{n}\|(\bm{\mathcal{Y}}_{G,n}+\bm{\mathcal{Y}}_{L,n})-(\bm{\mathcal{\hat{Y}}}_{G,n}+\bm{\mathcal{\hat{Y}}}_{L,n})\|_{F}^{2})/\sum_{n}\|\bm{\mathcal{Y}}_{G,n}+\bm{\mathcal{Y}}_{L,n}\|_{F}^{2}$

We run each experiment 10 times from $10$ different random seeds and report their mean and standard deviation. The results of the measuring statistics for different methods are summarized in Table 1. From Table 1, we can conclude the following statements.

1) perTucker yields the best results in separating the global and local components due to the utilization of the low-rank tensor structure of both components, from the following observations: (a) Comparing the global component error of perTucker and perPCA, we can conclude that perTucker identifies the global component with better accuracy due to its use of low-rank tensor structures. (b) Comparing the local component error of perTucker and robustTucker, although robustTucker identifies the global subspace with decent accuracy, it yields a much larger error in terms of local component reconstruction due to the fact that the local component does not assume any low-rank structure.

2) The reconstruction error for local components is larger for all methods compared to the global components. Two factors contribute to this: (a) the reconstruction of the local component is generally harder compared to the global components since it is only shared within the same clients, resulting in lower accuracy and larger variance with the use of fewer datasets; (b) the local components are generated by the shape variations of Swiss, oval, and rectangle patterns, which are not exactly low-rank, thus, the true local rank subspace is also an approximation from the dataset.

In this section, we evaluate the classification performance of the perTucker algorithm. The training sample size ranges from 10 to 50 with a step of 10 for each pattern. Then we perform the perTucker decomposition and get the local factor matrices. Next, we create 50 new images for each pattern and use the method described in Sec. 3.6.1 to classify the 150 images. This procedure is repeated 100 times. Table 2 displays the mean and standard deviation of the classification accuracy for various parameters. From Table 2 we can see that perTucker exhibits excellent classification performance even when the sample size is small. Increasing the training sample size leads to improved prediction accuracy. In comparison, if we perform local Tucker on the three clients and classify the new figures, the accuracy will be around 33% regardless of the training sample size. This accuracy is close to “guessing” the class. This is because local Tucker will model the commonality and peculiarity simultaneously. Therefore, global features are also considered when making the decision, while the multiplication coefficients are randomly generated.

To visualize the classification process, we use the boxplot to show the summary of the test statistics in (12) in Fig. 5. The sample size is set to 50. The test statistics are centralized by the median statistics within each pattern. The test statistics for each specific true pattern are significantly higher than those for the other two patterns in the corresponding cases, implying that the classification procedure is effective and robust.

In this section, we study the clustering performance of perTucker under a different problem setting. When generating the data, the variability within each pattern is measured by a “ratio” variable. The performance of clustering is evaluated by the ability to cluster similar “ratio” within each pattern. We focus on the Swiss pattern and cut the range of “ratio” from $0.7$ to $1.4$ into the group of $7$ ratio intervals, each with a width of $0.1$. Then for each ratio interval, we generate $3$ clients, each with a sample size of $100$ figures. We make clustering plots for the $21$ clients and see the aggregation of the clients with similar ratio, as shown in Fig. 6.

The first example involves monitoring solar activities and detecting solar flares from a stream of solar images. A solar flare is a significant event in the sun that emits enormous amounts of energetic charged particles that could cause power grids to fail (Marusek,, 2007). Therefore, detecting solar flares promptly is critical to implementing preemptive and corrective measures. However, monitoring solar activities is challenging due to the high dimensionality of the solar thermal imaging data and the gradual changes in solar temperature over time. Existing detection methods, which rely on subtracting the functional mean (background) using the sample mean, are insufficient to detect small transient flares in the dynamic system (Zhao et al., 2022a, ). Other studies focus on the decomposition of solar images into their background and anomaly components (Yan et al.,, 2018).

This dataset, publicly available in (Zhao et al., 2022a, ), comprises a sequence of images of size 232 × 292 pixels captured by a satellite. We use a sample of 300 frames in this case study. To detect solar flares in real time, we begin by subtracting the mean from each sample and then preprocess the data using the method proposed by (Aharon et al.,, 2006). Following this, we use a sliding window of $8\times 8$ to divide each frame into small patches, resulting in a total of 1044 patches. The four right-most columns of pixels are discarded. Each patch or tile is then vectorized, yielding a data dimension of $300\times 1044\times 64$. After the preprocessing step, we apply the perTucker algorithm to the data to break them down into two components: global components representing the slowly changing background and local components indicating the detected solar flares.

Fig. 7 shows two frames where there is an abrupt change in the almost-stationary background. On the original images, such changes are not visible in the complicated background. However, after using perTucker to extract global and local components, one can clearly see the location and movement of small and rapid changes in local components. The experiments highlight perTucker’s ability to change the signal magnification.

Moreover, since the local component represents the anomaly, we can use the Frobenius norm of the local core tensor $\left\lVert{\bm{\mathcal{C}}}_{L,n}\right\rVert_{F}$ as monitoring statistics to detect the anomaly. The results are shown in Fig. 8, where we plot the logarithm of the Frobenius norm of the local core tensor $\left\lVert{\bm{\mathcal{C}}}_{L,n}\right\rVert_{F}$ for $300$ frames.

For the most time, the curve in Fig. 8 remains at a low level, suggesting that no significant changes in solar activities are detected. However, there are 2 clear peaks above the red control line, suggesting that there are two solar flares. Therefore, Fig. 8 provides an intuitive approach for detecting anomalous behaviors.

perTucker also applies to monitoring tonnage profiles in a multi-operation forging process that utilizes multiple strain gauge sensors. Specifically, the forging process employs four columns, each equipped with a strain gauge sensor that measures the tonnage force exerted by the press uprights, as illustrated in Fig. 9. Consequently, each operational cycle generates a four-channel tonnage profile. The dataset for this case study consists of 305 in-control profiles, collected under normal production conditions, and 69 out-of-control profiles for each of the four different fault classes. The length of each channel profile is 1201, resulting in a data dimension of ${\bm{\mathcal{Y}}}^{4\times 1201\times 305}$ for the “normal” class and ${\bm{\mathcal{Y}}}^{4\times 1201\times 69}$ for each fault class. Fig. 10 presents examples of profiles for both normal and the four different fault conditions.

In this case study, we define the five “clients” as normal and four fault conditions. We select the first 50 normal samples and the first 10 samples from each fault condition to form the training dataset. The test dataset comprises $255$ normal samples and $59$ samples for each fault condition. We further assume that only the global and local factor matrices along the signal length dimension are orthogonal to each other.

Applying the classification method in Sec. 3.6.1, we are able to achieve 100% accuracy for classifying all fault modes and normal cases. Fig. 11 shows the relative test statistics by computing the difference between the test statistics of the corresponding fault modes and the normal cases. Consequently, normal test statistics have been normalized as $y=0$ for each sample, as shown in a red line in Fig. 11.

We can see that in every fault mode, the range of the test statistics for the corresponding fault mode is much higher than those of the other three faults and the normal baseline. When the true label is normal, the relative test statistics for all fault modes are less than the normal baseline. Thus, perTucker provides informative statistics for anomaly detection and fault diagnostics.