Semi-Orthogonal Multilinear PCA with Relaxed Start

Introduction

Principal component analysis (PCA) is a classical unsupervised dimensionality reduction method [Jolliffe, 2002]. It transforms input data into a new feature space of lower dimension via orthogonal projections, while keeping most variance of the original data. PCA is widely used in areas such as data compression [Kusner et al., 2014], computer vision [Ke and Sukthankar, 2004], and pattern recognition [Anaraki and Hughes, 2014; Deng et al., 2014].

Many real-world data are multi-dimensional, in the form of tensors rather than vectors [Kolda and Bader, 2009]. The number of dimensions of a tensor is the order and each dimension is a mode of it. For example, gray images are second-order tensors (matrices) and video sequences are third-order tensors [Lu et al., 2013]. Tensor data are also common in applications such as data center monitoring, social network analysis, and network forensics [Faloutsos et al.,

2007]. However, PCA on multi-dimensional data requires reshaping tensors into vectors first. This vectorization often leads to breaking of original data structures, more complex model with lots of parameters, and high computational and memory demands [Lu et al., 2013]. Many researchers address this problem via multilinear extensions of PCA to deal with tensors directly, and there are two main approaches.

One approach is based on Tensor-to-Tensor Projection (TTP) that learns low-dimensional tensors from highdimensional tensors. The two-dimensional PCA (2DPCA) [Yang et al., 2004] is probably the first PCA extension to deal with images without vectorization. The generalized low rank approximation of matrices (GLRAM) [Ye, 2005] and the generalized PCA (GPCA) [Ye et al., 2004] further generalize 2DPCA from single-sided projections to two-sided projections via reconstruction error minimization and variance maximization, respectively. Concurrent subspace analysis (CSA) [Xu et al., 2005] and multilinear PCA (MPCA) [Lu et al., 2008] extend GLRAM and GPCA to general higher-order tensors, respectively.

Another approach is based on Tensor-to-Vector Projection (TVP) that learns low-dimensional vectors from highdimensional tensors in a successive way. The tensor rankone decomposition (TROD) [Shashua and Levin, 2001] minimizes reconstruction error via (greedy) successive residue calculation. The uncorrelated multilinear PCA (UMPCA) [Lu et al., 2009] maximizes variance with the zero-correlation constraint, following the successive derivation of PCA. However, the number of features that can be extracted by UMPCA is upper-bounded by the lowest mode dimension. For example, for a tensor of size 300×200×3, UMPCA can only extract three features, which have very limited usage.

Orthogonality constraint is popular in feature extraction [Hua et al., 2007; Kokiopoulou and Saad, 2007; Gao et al., 2013], tensor decomposition [Kolda, 2001], and low-rank tensor approximation [Edelman et al., 1998; Wang et al., 2015]. PCA also obtains orthogonal projections, and the TTP-based PCA methods produce orthogonal projection vectors in each mode. However, none of the existing TVP-based PCA methods derive orthogonal projections. Our study found that it is indeed ineffective to impose full orthogonality in all the modes for TVP-based PCA, due to low captured variance and limited number of extracted features.

In this paper, we present a new TVP-based multilinear PCA algorithm, Semi-Orthogonal Multilinear PCA (SOMPCA) with Relaxed Start, or SO-MPCA-RS, to be detailed in Sec. 3. There are two main contributions:

• We propose a novel SO-MPCA approach to maximize the captured variance via TVP with orthogonality constraint in only one mode, which is called semiorthogonality according to [Wang et al., 2015]. The semi-orthogonality results in more captured variance and more learned features than full-orthogonality. For the same tensor of size 300×200×3 discussed earlier, SO-MPCA can extract 300 features while fullorthogonal multilinear PCA can only extract three features (similar to UMPCA).
• We introduce a Relaxed Start (RS) strategy to get SOMPCA-RS by fixing the starting projection vectors for better generalization [Abu-Mostafa et al., 2012]. This strategy constrains the hypothesis space to a smaller set, leading to increased bias and reduced variance of the learning model. The experimental results in Sec. 4 show that SO-MPCA-RS outperforms other competing PCAbased methods on the whole. In addition, this is a new strategy for tensor-based algorithms and we show its effectiveness for other TVP-based PCA methods.

In the following, we cover the necessary background first.

Background

Notations and basic operations: We follow the notations in [Lathauwer et al., 2000] to denote vectors by lowercase boldface letters, e.g., x; matrices by uppercase boldface letters, e.g., X; and tensors by calligraphic letters, e.g., X. We denote their elements with indices in parentheses, and indices by lowercase letters spanning the range from 1 to the uppercase letter of the index, e.g., n = 1, · · · ,N. An Nthorder tensor A ∈ RI1×···×IN is addressed by N indices {in}. Each in addresses the n-mode of A. The n-mode product of an Nth-order tensor A by a vector u ∈ RIn , denoted by B = A×n uT , is a tensor with entries:

Tensor-to-vector projection: Elementary multilinear projections (EMPs) are the building blocks of a TVP. We denote an EMP as {u(1) , u(2) , · · · , u(N) }, consisting of one unit projection vector in each mode, i.e., ‖ u(n) ‖= 1 for n = 1, · · · , N, where ‖ · ‖ is the Euclidean norm for vectors. It projects a tensor X ∈ RI1×I2×···×IN to a scalar y through the N unit projection vectors as [Lu et al., 2013]:

The TVP of a tensor X to a vector y ∈ RP consists of P EMPs {u(1) p , · · · , u(N) p }, p = 1, · · · , P, which can bewritten concisely as {u(n) p , n = 1, · · · ,N}Pp=1 or {u(n) p }Pp=1 :

where the pth component of y is obtained from the pth EMP as:

SO-MPCA with Relaxed Start

This section presents the proposed SO-MPCA-RS by first formulating the SO-MPCA problem, then deriving the solutions with a successive and conditional approach, and finally introducing the relaxed start strategy for better generalization.

Formulation of Semi-Orthogonal MPCA

We define the SO-MPCA problem with orthogonality constraint in only one mode, i.e., semi-orthogonality [Wang et al., 2015], as follows:

The SO-MPCA problem: A set of M tensor data samples {X1, X2, · · · ,XM} are available for training. Each sample Xm ∈ RI1×I2×···×IN can be viewed a point in the tensor space RI1 ⊗ RI2 · · · ⊗ RIN , where In is the nmode dimension and ⊗ denotes the Kronecker product. SOMPCA considers a TVP, which consists of P EMPs {u(n) p ∈ RIn×1 , n = 1, · · · ,N}Pp=1 , that projects the input tensor space RI1 ⊗ RI2 · · · ⊗ RIN into a vector subspace RP , i.e.,

for m = 1, · · · ,M. The objective is to find a TVP to maximize the variance of the projected samples in each projection direction, subject to the orthogonality constraint in only one mode, denoted as the ν-mode. The variance is measured by the total scatter Sp defined as:

where ymp = Xm ×Nn=1 {u(n) p }, and ¯ yp = 1M ∑m ymp .

In other words, the objective of SO-MPCA is to obtain the P EMPs, with the pth EMP determined as:

where the orthogonality constraint (9) is imposed only in the ν-mode and there is no such constraint for the other modes (n = 1, · · · ,N, n ̸= ν). The normalization constraint (8) is imposed for all modes.

Bound on the number of features: Based on the proof of Corollary 1 in [Lu et al., 2009], we can derive that the number of features P that can be extracted by SO-MPCA is upper-bounded by the ν-mode dimension Iν: P ≤ Iν. Since we can choose any n as ν, we have the upper bound of P as P ≤ maxn In (i.e., the highest mode dimension).

Selection of mode ν: Although we are free to choose any mode n as ν to impose the orthogonality constraint (9), it is often good to have more features in practice. Thus, in this paper, we choose the mode with the highest dimension as ν:

such that P = maxn In = Iν. On the other hand, we can also obtain a total of ∑n In features by running SO-MPCA N times with ν = 1, · · · ,N. In this paper, we only focus on SO-MPCA with ν determined by (10).

Semi-orthogonality vs. full-orthogonality: If we impose the orthogonality constraint (9) in all modes, we can get Full-Orthogonal Multilinear PCA (FO-MPCA). However, our study found that FO-MPCA is not effective primarily due to two reasons:

• Due to the heavy constraints, the variance captured by FO-MPCA is quite low, even lower than UMPCA. In contrast, SO-MPCA can capture more variance than both FO-MPCA and UMPCA. This is illustrated in Fig. 1 in Sec. 4.
• Similar to UMPCA, the number of features that can be extracted by FO-MPCA is upper-bounded by the lowest mode dimension minn In, which can be quite limited. For instance, FO-MPCA can extract only three features for a tensor of size 300×200×3 while SO-MPCA can extract 300 features by choosing ν = 1 for the same tensor. This can be observed in Fig. 1 as well.

Successive Derivation of SO-MPCA

To solve the SO-MPCA problem, we follow the successive derivation in [Jolliffe, 2002; Lu et al., 2009] to determine EMPs one by one in P steps:

Step 1 (p = 1): Determine the first EMP {u(n) 1 , n = 1, · · · ,N} by maximizing S1 with the constraint (8).

Determine the pth EMP {u(n) p , n = 1, · · · ,N} by maximizing Sp with the constraints (8) and (9).

Conditional subproblem: In order to obtain the pth EMP {u(n) p , n = 1, · · · ,N}, we need to determine N vectors. We follow the approach of alternating least squares [Harshman, 1970]. Thus, we can only obtain locally optimal solutions as in many other tensor-based methods. For the pth EMP, the parameters of the n-mode projection vector u(n) p are estimated one mode by one mode separately conditioned on the projection vectors in all the other modes. Assuming the pth projection vectors in all but n-mode are given, we project the input tensor samples in these (N - 1) modes to obtain the partial multilinear projections as in [Lu et al., 2013]:

where ˜y(n) mp ∈ RIn . This conditional subproblem then becomes to determine u(n) p that projects the vector samples {˜y(n) mp ,m = 1, · · · ,M} onto a line to maximize the variance captured. Then the total scatter matrix ˜S(n) p corresponding to {˜y(n) mp ,m= 1, · · · ,M} becomes:

For p = 1 (step 1), the solution for u(n) 1 , where n = 1, · · · ,N, is obtained as the unit eigenvector of ˜S(n) 1 associated with the largest eigenvalue.

For p ≥ 2, we need to deal with the ν-mode and other modes differently. For modes other than ν, the solution for u(n) p , where n = 1, · · · ,N, n ̸= ν, is obtained as the unit eigenvector of ˜S(n) p associated with the largest eigenvalue.

Constrained optimization for ν-mode and p ≥ 2: When p ≥ 2, we need to determine u(ν) p by solving the following constrained optimization problem:

We solve this problem by the following theorem:

Theorem 1. The solution to the problem (13) is the (unitlength) eigenvector corresponding to the largest eigenvalue of the following eigenvalue problem:

where,

and IIn is an identity matrix of size In × In.

Proof. First, we use Lagrange multipliers to transform the problem (13) to include all the constraints as:

where λ and {µq, q = 1, · · · , p-1} are Lagrange multipliers.

Then we set the partial derivative of Lν with respect to u(ν) p to zero:

Premultiplying (17) by u(ν) p T , the third term vanishes and we get

Algorithm 1 Semi-Orthogonal Multilinear PCA with Relaxed Start (SO-MPCA-RS)

1: Input: A set of tensor samples {Xm ∈ RI1×···×IN ,m = 1, · · · ,M}, and the maximum number of iterations K. 2: Set ν = argmaxn In. 3: Set the first EMP: u(n) 1 = 1/ ‖ 1 ‖ for n = 1, · · · ,N . 4: for p = 2 to P do 5: Initialize u(n) p = 1/ ‖ 1 ‖ for n = 1, · · · ,N. 6: for k = 1 to K do 7: for n = 1 to N do 8: Calculate the partial multilinear projection {˜y(n) mp } for m= 1, · · · ,M according to (11). 9: if n == ν then 10: Calculate Γ(ν) p and ˜S(ν) p according to (15) and (12), respectively. Then, set u(ν) p to the eigenvector of Γ(ν) p ˜S(ν) p associated with the largest eigenvalue. 11: else 12: Calculate ˜S(n) p by (12). Set u(n) p to the eigenvector of ˜S(n) p associated with the largest eigenvalue. 13: end if 14: end for 15: end for 16: end for 17: Output The TVP {u(n) p , n = 1, · · · ,N}Pp=1 .

which indicates that λ is exactly the criterion to be maximized, with the orthogonality constraint.

Next, a set of (p - 1) equations are obtained by premultiplying (17) by u(ν) q T , q = 1, · · · , p - 1, respectively,

The second term vanishes and the summand in the third term is non-zero only for s = q. Thus, we get

Substituting (20) into (17), we get

Using the definition in (15), (22) can be rewritten as:

Since λ is the criterion to be maximized, this maximization is achieved by setting u(ν) p to the (unit) eigenvector of Γ(ν) p ˜S(ν) p associated with its corresponding largest eigenvalue .

Relaxed Start for Better Generalization

When we use SO-MPCA features for classification, we find the performance is limited. Therefore, we further introduce a simple relaxed start (RS) strategy to get SO-MPCA-RS by fixing the first EMP {u(n) 1 , n = 1, · · · ,N} (the starting projection vectors), without variance maximization. In this paper, we set this starting EMP u(n) 1 (for n = 1, · · · ,N) to the normalized uniform vector 1/ ‖ 1 ‖ for simplicity.

This idea is motivated by the theoretical studies in Chapter 4 of [Abu-Mostafa et al., 2012] showing that constraining a learning model could lead to better generalization. By fixing the first EMP as simple vectors, the following EMPs have less freedom due to the imposed semi-orthogonality, which increases the bias and reduces the variance of the learning model. Thus, the SO-MPCA-RS model has a smaller hypothesis set than the SO-MPCA model. The two algorithms differ only in how to determine the first (starting) EMP though the following EMPs will all be different due to their dependency on the first EMP.

This relaxed start strategy is not specific to SO-MPCA but generally applicable to any TVP-based subspace learning algorithm. We run controlled experiments in Sec. 4 to show that it can improve the performance of not only SO-MPCA but also TROD and UMPCA.

Algorithm 1 summarizes the SO-MPCA-RS algorithm. 1 The SO-MPCA algorithm can be obtained from Algorithm 1 by removing line 3, changing p = 2 in line 4 to p = 1 and setting Γ(ν) 1 (p = 1) in line 10 to an identity matrix.

Experiments

This section evaluates the proposed methods on both secondorder and third-order tensor data in terms of recognition rate, the number of extracted features, captured variance, and convergence. In addition, we also study the effectiveness of the relaxed start strategy on other TVP-based PCA algorithms.

Data: 2 For second-order tensors, we use the same subset of the FERET database [Phillips et al., 2000] as in [Lu et al., 2009], with 721 face images from 70 subjects. Each face image is normalized to 80 × 60 graylevel pixels. For thirdorder tensors, we use a subset of the USF HumanID "Gait Challenge" database [Sarkar et al., 2005]. We use the same gallery set (731 samples from 71 subjects) and probe A (727 samples from 71 subjects) as in [Lu et al., 2009], and we also test probe B (423 samples from 41 subjects) and probe C (420 samples from 41 subjects). Each gait sample is a (binary) silhouette sequence of size of 32×22×10.

Experiment setup: In face recognition experiments, we randomly select L = 1, 2, 3, 4, 5, 6, 7 samples from each subject as the training data and use the rest for testing. We repeat such random splits (repetitions) ten times and report the mean correct recognition rates. In gait recognition experiments, we follow the standard setting and use the gallery set as the training data and probes A, B, and C as the test data (so there is no random splits/repetitions), and report the rank 1 and rank 5 recognition rates [Sarkar et al., 2005].

1 Matlab code is available at: http://www.comp.hkbu.edu.hk/ ∼haiping/codedata.html

2 Both face and gait data are downloaded from: http://www.dsp. utoronto.ca/∼haiping/MSL.html

Algorithms and their settings: We first evaluate SOMPCA and SO-MPCA-RS against five existing PCA-based methods: PCA [Jolliffe, 2002], CSA [Xu et al., 2005], MPCA [Lu et al., 2008], TROD [Shashua and Levin, 2001], and UMPCA [Lu et al., 2009]. 3 CSA and MPCA produce tensorial features so they need to be vectorized. MPCA uses the full projection. For TROD and UMPCA, we use the uniform initialization [Lu et al., 2009]. For SO-MPCA and SOMPCA-RS, we set the selected mode ν = 1 for the maximum number of features. For iterative algorithms, we set the number of iterations to 20. All features are sorted according to the scatters (captured variance) in descending order for classification. We use the Nearest Neighbor Classifier with the Euclidean distance measure to classify the top P features. We test up to P = 80 features in face recognition and up to P = 32 features in gait recognition. The performance of FOMPCA is much worse than SO-MPCA so it is not included in the comparisons (except variance study) to save space.

Face recognition results: Table 1 shows the face recognition results for P = 1, 5, 10, 20, 50, 80 and L = 1, 2, 3, 4, 5, 6, 7, including both the mean and the standard deviation (std) over ten repetitions. We highlight the top two results in each row in bold fonts for easy comparison. Only SOMPCA-RS consistently achieves the top 2 results in all cases. Compared with existing methods (PCA, CSA, MPCA, TROD and UMPCA), SO-MPCA-RS outperforms the best performing existing algorithm (UMPCA) by 3.79% on average.

3 For second-order tensors, CSA and MPCA are equivalent to GLRAM [Ye, 2005] and GPCA [Ye et al., 2004], respectively.

Furthermore, for larger L = 5, 6, 7, SO-MPCA-RS outperforms the other five methods at least by 2.26% on average. For smaller L = 1, 2, 3, SO-MPCA-RS achieves a greater improvement of at least 5.28% over existing methods, indicating that SO-MPCA-RS is more superior in dealing with the small sample size (overfitting) problem.

Gait recognition results: Similarly, the gait recognition results are reported in Table 2 with the top two results highlighted. Again, only SO-MPCA-RS consistently achieves the top 2 results in all cases. In rank 1 rate, the best performing existing algorithm is PCA, which outperforms SO-MPCARS by 3.73% on average. While in rank 5 rate, SO-MPCARS outperforms the best performing existing algorithm (still PCA) by 6.76% on average.

Number of features: In the tables, we use '-' to indicate that there are not enough features. For PCA, there are at most 69 features for face data when L = 1 since there are only 70 samples for training. UMPCA can only extract 60 or 10 features for face and gait data, respectively. In contrast, SOMPCA and SO-MPCA-RS (with ν = 1) can learn 80 features for face data and 32 features for gait data.

Feature variance: We illustrate the variance captured by PCA, UMPCA, FO-MPCA, SO-MPCA, and SO-MPCA-RS in Fig. 1 for face data with L = 1 (not all methods are shown for clarity). Figure 1(a) shows the sorted variance. It is clear that semi-orthogonality captures more variance than full-orthogonality, as we discussed in Sec. 3.1. Moreover, both SO-MPCA and SO-MPCA-RS can capture more variance than UMPCA, but less than PCA (and also CSA, MPCA, and TROD, which are not shown). Though capturing less variance, SO-MPCA-RS achieves better overall classification performance than other PCA-based methods, with results consistently in the top two in all experiments.

We also show the unsorted captured variance in Fig. 1(b).

The variance captured by the first (fixed) EMP of SO-MPCARS is much less than other EMPs, which is not surprising since the variance is not maximized.

Convergence: We demonstrate the convergence of SOMPCA-RS in Fig. 2 for face data with L = 1. We can see that SO-MPCA-RS converges in just a few iterations. SO-MPCA has a similar convergence rate.

Effectiveness of relaxed start: To evaluate the proposed relaxed start strategy, we apply it to two other TVPbased methods, TROD and UMPCA, getting TROD-RS and UMPCA-RS, respectively. We summarize their performance in the last two columns of Tables 1 and 2 for face and gait recognition experiments, respectively.

Both tables show that relaxed start can help both TROD and UMPCA to achieve better recognition rates. From Table 1, TROD-RS improves over TROD by 0.32% and UMPCARS improves over UMPCA by 2.15% on average. From Table 2, TROD-RS achieves 3.03 % improvement over TROD and UMPCA-RS achieves 1.07% improvement over UMPCA for rank 1 rate. For rank 5 rate, TROD-RS improves 0.94% over TROD, and UMPCA-RS improves 5.28% over UMPCA. The relaxed start is most effective for our SO-MPCA. SO-MPCARS has an improvement of 9.97% on face data over SOMPCA. On gait data, SO-MPCA-RS outperforms SO-MPCA by 9.56% in rank 1 rate and 17.26% in rank 5 rate on average.

In addition, SO-MPCA-RS has better face recognition performance than TROD-RS and UMPCA-RS with an improvement of 8.08% and 1.64%, respectively. On gait data, SOMPCA-RS improves the rank 1 recognition rate by 3.03% over TROD-RS and 13.25% over UMPCA-RS on average, and SO-MPCA-RS improves the rank 5 recognition rate by 11.07% and 13.73% over TROD-RS and UMPCA-RS.

These controlled experiments show the effectiveness of relaxed start on SO-MPCA and other TVP-based multilinear PCA methods (UMPCA and TROD). One possible explanation is that RS increases the bias and reduces the variance of the learning model, while further investigation is needed.

Conclusion

This paper proposes a novel multilinear PCA algorithm under the TVP setting, named as semi-orthogonal multilinear PCA with relaxed start (SO-MPCA-RS). The proposed SO-MPCA approach learns features directly from tensors via TVP to maximize the captured variance with the orthogonality constraint imposed in only one mode. This semi-orthogonality can capture more variance and learn more features than fullorthogonality. Furthermore, the introduced relaxed start strategy can achieve better generalization by fixing the starting projection vectors to uniform vectors to increase the bias and reduce the variance of the learning model. Experiments on face (2D data) and gait (3D data) recognition show that SOMPCA-RS achieves the best overall performance compared with competing algorithms. In addition, relaxed start is also effective for other TVP-based PCA methods.

In this paper, we studied semi-orthogonality in only one mode. A possible future work is to learn SO-MPCA-RS features from each mode separately and then do a feature/scorelevel fusion.

Acknowledgments

This research was supported by Research Grants Council of the Hong Kong Special Administrative Region (Grant 22200014 and the Hong Kong PhD Fellowship Scheme).