Multilinear Common Component Analysis via Kronecker Product Representation

Various statistical methodologies for extracting useful information from a large amount of data have been studied over the decades since the appearance of big data. In the present era, it is important to discover a common structure of multiple data sets. In an early study, Flury (1984) focused on the structure of the covariance matrices of multiple data sets and discussed the heterogeneity of the structure. The author reported that population covariance matrices differ among multiple data sets in practical applications. Many methodologies have been developed for treating the heterogeneity between covariance matrices of multiple data sets (see, Flury, 1986, 1988; Flury & Gautschi, 1986; Pourahmadi, Daniels, & Park, 2007; Wang, Banerjee, & Boley, 2011; Park & Konishi, 2020).

Among such methodologies, common component analysis (CCA; Wang et al., 2011) is an effective tool for statistics. The central idea of CCA is to reduce the number of dimensions of data while losing as little information of the multiple data sets as possible. To reduce the number of dimensions, CCA reconstructs the data with a few new variables that are linear combinations of the original variables. For considering the heterogeneity between covariance matrices of multiple data sets, CCA assumes that there is a different covariance matrix for each data set. There have been many papers on various statistical methodologies using multiple covariance matrices: discriminant analysis (Bensmail & Celeux, 1996), spectral decomposition (Boik, 2002), and a likelihood ratio test for multiple covariance matrices (Manly & Rayner, 1987). It should be noted that principal component analysis (PCA) (Pearson, 1901; Jolliffe, 2002) is a technique similar to CCA. In fact, CCA is a generalization of PCA; PCA can only be applied to one data set, whereas CCA can be applied to multiple data sets.

Meanwhile, in various fields of research, including machine learning and computer vision, the main interest has been in tensor data, which has a multidimensional array structure. In order to apply the conventional statistical methodologies, such as PCA, to tensor data, a simple approach is to first transform the tensor data into vector data and then apply the methodology. However, such an approach causes the following problems:

1. In losing the tensor structure of the data, the approach ignores the higher-order inherent relationships of the original tensor data.

2. Transforming tensor data to vector data substantially increases the number of features. It also has a high computational cost.

To overcome these problems, statistical methodologies for tensor data analyses have been proposed that take the tensor structure of the data into consideration. Such methods enable us to accurately extract higher-order inherent relationships in a tensor data set. In particular, many existing statistical methodologies have been extended for tensor data, for example, multilinear principal component analysis (MPCA) (Lu et al., 2008) and sparse PCA for tensor data analysis (Allen, 2012; Wang, Sun, Chen, Pang, & Zhou, 2012; Lai, Xu, Chen, Yang, & Zhang, 2014), as well as others (see Carroll & Chang, 1970; Harshman, 1970; Kiers, 2000; Badeau & Boyer, 2008; Kolda & Bader, 2009).

In this letter, we extend CCA to tensor data analysis, proposing multilinear common component analysis (MCCA). MCCA discovers the common structure of multiple data sets of tensor data while losing as little of the information of the data sets as possible. To identify the common structure, we estimate a common basis constructed as linear combinations of the original variables. For estimating the common basis, we develop a new estimation algorithm based on the idea of CCA. In developing the estimation algorithm, two issues must be addressed: the convergence properties of the algorithm and its computational cost. To determine the convergence properties, we investigate first the relationship between the initial values of the parameters and global optimal solution and then the monotonic convergence of the estimation algorithm. These analyses reveal that our proposed algorithm guarantees convergence of the mode-wise global optimal solution under some conditions. To analyze the computational efficacy, we calculate the computational cost of our proposed algorithm.

The rest of the letter is organized as follows. In section 2, we review the formulation and the minimization problem of CCA. In section 3, we formulate the MCCA model by constructing the covariance matrices of tensor data, based on a Kronecker product representation. Then we formulate the estimation algorithm for MCCA in section 4. In section 5, we present the theoretical properties for our proposed algorithm and analyze the computational cost. The efficacy of the MCCA is demonstrated through numerical experiments in section 6. Concluding remarks are presented in section 7. Technical proofs are provided in the appendixes. Our implementation of MCCA and supplementary materials are available at https://github.com/yoshikawa-kohei/MCCA.

where $tr(·)$ denotes the trace of a matrix. A crucial issue for solving the maximization problem 2.4 is the nonconvexity. Certainly the maximization problem is nonconvex since the problem is defined on a set of orthogonal matrices, which is a nonconvex set. Generally it is difficult to find the global optimal solution in nonconvex optimization problems. To overcome this drawback, Wang et al. (2011) proposed an estimation algorithm in which the estimated parameters are guaranteed to constitute the global optimal solution under some conditions.

In this section, we introduce a mathematical formulation of the MCCA, an extension of the CCA in terms of tensor data analysis. Moreover, we formulate an optimization problem of MCCA and investigate its convergence properties.

To solve the maximization problem efficiently and identify the inherent relationships, the maximization problem 3.5 can be decomposed into the mode-wise maximization problems represented in the following lemma.

Our estimation algorithm consists of two steps: initializing the parameters and iteratively updating the parameters. The initialization step gives us the initial values of the parameters near the global optimal solution for each mode. Next, by iteratively updating the parameters, we can monotonically increase the value of the objective function 3.7 until convergence.

Using the initial value of $w(k)$⁠, we can obtain the initial value of the parameter $V0(k)$ by maximizing $fk'(V(k))$ for each mode. The maximizer consists of $Rk$ eigenvectors, corresponding to the $Rk$ largest eigenvalues, obtained by eigenvalue decomposition of $MI(k)$⁠. The theoretical justification for this initialization is discussed in section 5.

Our estimation procedure comprises the above estimation steps. The procedure is summarized as algorithm 1.

This section presents the theoretical and computational analyses for algorithm 1. Theoretical analyses consist of two steps. First, we prove that the initial values of parameters obtained in section 4.1 are relatively close to the global optimal solution. If the initial values are close to the global maximum, then we can obtain the global optimal solution even if the maximization problem is nonconvex. Second, we prove that the iterative updates of the parameters in section 4.2 monotonically increase the value of objective function 3.7 by solving the surrogate problem 4.3. From the monotonically increasing property, the estimated parameters always converge at a stationary point. The combination of these two results enables us to obtain the mode-wise global optimal solution. In the computational analysis, we calculate computational cost for MCCA and then compare the cost with conventional methods. By comparing the costs, we investigate the computational efficacy of MCCA.

This section aims to provide the upper and lower bounds of the maximization problem 3.7. From the bounds, we find that the initial values in section 4.1 are relatively close to the global optimal solution. Before providing the bounds, we define a contraction ratio.

Note that a contraction ratio $α(k)$ satisfies $0≤α(k)≤1$ and $α(k)=1$ if and only if $Rk=Pk$⁠.

Using $fk'max$ and the contraction ratio $α(k)$⁠, we have the following theorem that reveals the upper and lower bounds of the global maximum in problem 3.7.

This theorem indicates that $fk'max→fkmax$ when $α(k)→1$⁠. Thus, it is important to obtain an $α(k)$ that is as close as possible to one. Since $α(k)$ depends on $V0(k)$ and $w(k)$⁠, $V0(k)$ depends on $w(k)$⁠. From this dependency, if we could set the initial value of $w(k)$ such that $α(k)$ is as large as possible, then we could obtain an initial value of $V0(k)$ that attains a value near $fkmax$⁠. The following theorem shows that we can compute the initial value of $w(k)$ such that $α(k)$ is maximized.

In fact, $α(k)$ is very close to one with the initial values given in theorem 2 even if $Rk$ is small. This resembles the cumulative contribution ratio in PCA.

We next verify that our proposed procedure for iteratively updating parameters maximizes the optimization problem 3.7. In algorithm 1, the parameter $Vs+1(k)$ can be obtained by solving the surrogate maximization problem 4.3. Theorem 3 shows that we can monotonically increase the value of the function $fk(V(k))$ in equation 3.7 by algorithm 1.

From theorem 1, we obtain initial values of the parameters that are near the global optimal solution. By combining theorems 1 and 3, the solution from algorithm 1 can be characterized by the following corollary.

Algorithm 1 does not guarantee the global solution due to the fundamental problem of nonconvexity, but it is enough for pragmatic purposes. We investigate the issue of convergence to global solution through numerical studies in section 6.3.

First, we analyze the computational cost. To simplify the analysis, we assume $P=argmaxjPj$ for $j=1,2,…,M$⁠. This implies that $P$ is the upper bound of $Rj$ for all $j$⁠. We then calculate the upper bound of the computational complexity.

The expensive computations of the each iteration in algorithm 1 consist of three parts: the formulation of $M(Vs(k))$⁠, the eigenvalue decomposition of $M(Vs(k))$⁠, and updating latent covariance matrices $Λg(k)$⁠. These steps are $O(GM2P3)$⁠, $O(P3)$⁠, and $O(GMP3)$⁠, respectively. The total computational complexity per iteration is then $O(GM2P3)$⁠.

Next, we analyze the memory requirement of algorithm 1. MCCA represents the original tensor data with fewer parameters by projecting the data onto a lower-dimensional space. This requires the $Pk×Rk$ projection matrices $V(k)$ for $k=1,2,…,M$⁠. MCCA projects the data with size of $N∏k=1MPk$ to $N∏k=1MRk$⁠, where $N=∑g=1GNg$⁠. Thus, the required size for the parameters is $∑k=1MPkRk+N∏k=1MRk$⁠. MPCA requires the same amount of memory as MCCA. Meanwhile, CCA and PCA need a projection matrix, which is size $R∏k=1MPk$⁠. The required size for the parameters is then $R∏k=1MPk+NR$⁠.

To compare the computational cost clearly, the upper bounds of computational complexity and the memory requirement are summarized in Table 1. Table 1 shows that the computational complexity of MCCA is superior to that of the other algorithms and the complexity of MCCA is not limited by sample size. In contrast, the MPCA algorithm is affected by the sample size (Lu, Plataniotis, & Venetsanopoulos, 2008). Additionally, MCCA and MPCA require a large amount of memory when the number of modes in a data set is large, but their memory requirements are much smaller than those of PCA and CCA.

To demonstrate the efficacy of MCCA, we applied MCCA, PCA, CCA, and MPCA to image compression tasks.

For the experiments, we prepared the following three image data sets:

• MNIST data set consists of data of handwritten digits $0,1,…,9$ at image sizes of $28×28$ pixels. The data set includes a training data set of 60,000 images and a test data set of 10,000 images. We used the first 10 training images of the data set for each group. The MNIST data set (Lecun, Bottou, Bengio, & Haffner, 1998) is available at http://yann.lecun.com/exdb/mnist/.

• AT&T (ORL) face data set contains gray-scale facial images of 40 people. The data set has 10 images sized $92×112$ pixels for each person. We used images resized by a factor of 0.5 to improve the efficiency of the experiment. The AT&T face data set is available at https://git-disl.github.io/GTDLBench/datasets/att_face_dataset/. All the credits of this data set go to AT&T Laboratories Cambridge.

• Cropped AR database has color facial images of 100 people. These images are cropped around the face. The size of images is $120×165×3$ pixels. The data set contains 26 images in each group, 12 of which are images of people wearing sunglasses or scarves. We used the cropped facial images of 50 males who were not wearing sunglasses or scarves. Due to memory limitations, we resized these images by a factor of 0.25. The AR database (Martinez & Benavente, 1998; Martinez & Kak, 2001) is available at http://www2.ece.ohio-state.edu/∼aleix/ARdatabase.html.

The data set characteristics are summarized in Table 2.

To compress these images, we performed dimensionality reductions by MCCA, PCA, CCA, and MPCA, as follows. We vectorized the tensor data set before performing PCA and CCA. In MCCA, the images were compressed and reconstructed according to the following steps:

1. Prepare the multiple image data sets $X(g)∈RP1×P2×⋯×PM×Ng$ for $g=1,2,…,G$⁠.

2. Compute the covariance matrix of $X(g)$ for $g=1,2,…,G$⁠.

3. From these covariance matrices, compute the linear transformation matrices $Vi∈RPi×Ri$ for $i=1,2,…,M$ for mapping to the $(R1,R2,…,RM)$-dimensional latent space.

4. Map the $i$th sample $X(g)i$ to $X(g)i×1V1×2V2⋯×MVM∈RR1×R2×⋯×RM$⁠, where the operator $×i$ is the $i$-mode product of tensor (Kolda & Bader, 2009).

5. Reconstruct $i$th sample $X˜(g)i$ = $X(g)i×1V1V1⊤×2V2V2⊤⋯×MVMVM⊤$⁠.

Meanwhile, PCA and MPCA each require a single data set. Thus, we aggregated the data sets as $X=[X(1),X(2),…,X(G)]∈RP1×P2×⋯×PM×∑g=1GNg$ and performed PCA and MPCA for data set $X$⁠.

For MCCA and MPCA, the reduced dimensions $R1$ and $R2$ were chosen as the same number, and then we fixed $R3$ as two. All computations were performed by the software R (version 3.6) (R Core Team, 2019). In the initialization of MCCA, solving the quadratic programming problem was carried out using the function ipop in the package kernlab. MPCA was implemented as the function mpca in the package rTensor. (The implementations of MCCA, PCA, and CCA are available at https://github.com/yoshikawa-kohei/MCCA.)

From Figure 1, we observe that the RER material MCCA is the smallest for any value of CR. This indicates that MCCA performs better than the other methods. In addition, note that CCA performs better than MPCA only for fairly small values of CR, even though it is a method for vector data, whereas MPCA performs better for larger values of CR. This implies the limitations of CCA for vector data.

Next we consider group size by comparing panels a, b, and c in Figure 1. The value of CR at the intersection of CCA and MPCA increases with increasing the group size. This indicates that MPCA has more trouble extracting an appropriate latent space as the group size increases. Since MPCA does not consider the group structure, it is not possible to properly estimate the covariance structure when the group size is large.

We examined the behavior of contraction ratio $α(k)$⁠. We performed MCCA on the AT&T(ORL) data set with the medium group size and computed $α(1)$ and $α(2)$ with the various pairs of reduced dimensions $(R1,R2)∈{1,2,…,25}×{1,2,…,25}$⁠.

We investigated the usefulness of determining the initial value of $w(k)$ by solving the quadratic programming problem 4.1. We applied MCCA to the AT&T(ORL) data set with the small, medium, and large number of groups. In addition, we used the smaller group size of three. For determining the initial value of $w(k)$⁠, we consider three methods: solving the quadratic programming problem 4.1 (MCCA:QP); setting all values of $w(k)$ to one (MCCA:FIX); and setting the values by random sampling according to the uniform distribution $U(0,1)$ (MCCA:RANDOM). We computed the $α(k)$ with the reduced dimensions $R1=R2(∈1,2,…,10)$ for each of these methods.

To evaluate the performance of these methods, we compared the values of $α(k)$ and the number of iterations in the estimation. The number of iterations in the estimation is the number of repetitions of lines 7 to 9 in algorithm 1. For MCCA(RANDOM), we performed 50 trials and calculated averages of each of these indices.

We have developed the multilinear common components analysis (MCCA) by introducing a covariance structure based on the Kronecker product. To efficiently solve the nonconvex optimization problem for MCCA, we have proposed an iteratively updating algorithm that exhibits some superior theoretical convergence properties. Numerical experiments have shown the usefulness of MCCA.

Specifically, MCCA was shown to be competitive among the initialization methods in terms of the number of iterations. As the number of groups increases, the overall number of samples increases. This may be the reason why all methods required almost the same number of iterations for small, medium, and large groups. Note that in this study, we used the Kronecker product representation to estimate the covariance matrix for tensor data sets. Greenewald, Zhou, and Hero (2019) used the Kronecker sum representation for estimating the covariance matrix, and it would be interesting to extend the MCCA to this and other covariance representations.

We provide two basic lemmas about Kronecker products before we prove lemma 1.

These lemmas are known as the mixed-product property and the spectrum property, respectively. See Harville (1998) for detailed proofs.

Theorem 1 can be easily shown from the following lemma.