A Convengent Solution to Tensor Subspace Learning

1 Introduction

Subspace learning algorithms [Brand, 2003] such as Principal Component Analysis (PCA) [Turk and Pentland, 1991] and Linear Discriminant Analysis (LDA) [Belhumeur et al., 1997] traditionally express the input data as vectors and often in a high-dimensional feature space. In real applications, the extracted features are usually in the form of a multidimensional union, i.e. a tensor, and the vectorization process destroys this intrinsic structure of the original tensor form. Another drawback brought by the vectorization process is the curse of dimensionality which may greatly degrade the algorithmic learnability especially in the small sample size cases.

Recently substantial efforts have been devoted to the employment of tensor representation for improving algorithmic learnability [Vasilescu and Terzopoulos, 2003]. Among them, 2DLDA [Ye et al., 2004] and DATER [Yan et al., 2005] are tensorized from the popular vector-based LDA algorithm. Although the initial objectives of these algorithms are different, they all end up with solving a higherorder optimization problem, and commonly iterative procedures were used to search for the solution. A collective problem encountered by their solution procedures is that the iterative procedures are not guaranteed to converge, since in each iteration, the optimization problem is approximately simplified from the Trace Ratio form arg maxU k T r(U kT S p k U k )/T r(U kT S k U k ) 1 to the Ratio Trace form argmaxU k Tr[ (U kT S kU k )-1 (U kT S p k U k ) ] in order to obtain a closed-form solution for each iteration. Consequently, the derived projection matrices are unnecessary to converge, which greatly limits the application of these algorithms since it is unclear how to select the iteration number and the solution is not optimal even in the local sense.

In this work, by following the graph embedding formulation for general dimensionality reduction proposed by [Yan et al., 2007], we present a new solution procedure for subspace learning based on tensor representation. In each iteration, instead of transforming the objective function into the ratio trace form, we transform the trace ratio optimization problem into a trace difference optimization problem maxU k T r[U kT (S p k -λS k )U k ] where λ is the objective function value computed from the solution (U k |nk=1 ) of the previous iteration. Then, each iteration is efficiently solved with the eigenvalue decomposition method [Fukunaga, 1991]. A detailed proof is presented to justify that λ, namely the value of the objective function, will increase monotonously, and also we prove that the projection matrix U k will converge to a fixed point based on the point-to-set map theories [Hogan, 1973].

It is worthwhile to highlight some aspects of our solution procedure to general subspace learning based on tensor representation here:

1. The value of the objective function is guaranteed to monotonously increase; and the multiple projection matrices are proved to converge. These two properties ensure the algorithmic effectiveness and applicability.
2. Only eigenvalue decomposition method is applied for iterative optimization, which makes the algorithm extremely efficient; and the whole algorithm does not suffer from the singularity problem that is often encountered by the traditional generalized eigenvalue decomposition method used to solve the ratio trace optimization problem.
3. The consequent advantage brought by the sound theoretical foundation is the enhanced potential classification

1 Matrices Sp k and Sk are both positive semidefinite and more detailed definitions are described afterward.

capability of the derived low-dimensional representation from the subspace learning algorithms.

The rest of this paper is organized as follows. Section II reviews the general subspace learning based on tensor representation, and then we introduce our new solution procedure along with the theoretical convergency proof in section III. By taking the Marginal Fisher Analysis (MFA) algorithm proposed in [Yan et al., 2007] as an example, we verify the convergency properties of the new proposed solution procedure and the classification capability of the derived lowdimensional representation is examined with a set of experiments on the real-world databases in Section IV.

2 Subspace Learning with Tensor Data

In this section, we present a general subspace learning framework by encoding data as tensors of arbitrary order, extended from the one proposed by [Yan et al., 2007] and taking the data inputs as vectors. The concepts of tensor inner production, mode-k production with matrix, and mode-k unfolding are referred to the work of [Yan et al., 2005].

2.1 Graph Embedding with Tensor Representation

Denote the sample set as X = [X1,X2, . . . ,XN],Xi ∈ Rm1×m2×...×mn , i = 1, . . . N, with N as the total number of samples. Let G = {X, S} be an undirected similarity graph, called an intrinsic graph, with vertex set X and similarity matrix S ∈ RN×N . The corresponding diagonal matrix D and the Laplacian matrix L of the graph G are defined as

The task of graph embedding is to determine a lowdimensional representation of the vertex set X that preserves the similarities between pairs of data in the original highdimensional feature space. Denote the low-dimensional embedding of the vertices as Y = [Y1,Y2, . . . ,YN], where Yi ∈ Rm′1×m′2×...×m′n is the embedding for the vertex Xi, with the assumption that Yi is the mode-k production of Xi with a series of column orthogonal matrices U k ∈ Rmk×m′k ,

where Im′k is anm′k -by-m′k identitymatrix. Tomaintain similarities among vertex pairs according to the graph preserving criterion [Yan et al., 2007], we have

where f(U k |nk=1 ) is a function that poses extra constraint for the graph similarity preserving criterion. Here U k |nk=1 means the sequence U 1 , U 2 to U n and so for the other similar representations in the following parts of this work. Commonly, f (Uk|nk=1 ) may have two kinds of definitions. One is for scale normalization, that is,

where B is a diagonal matrix with non-negative elements. The other is a more general constraint which relies on a new graph, referred to as penalty graph with similarity matrix S p , and is defined as

Without losing generality, we assume that the constraint function is defined with penalty matrix for simplicity; and for scale normalization constraint, we can easily have the similar deduction for our new solution procedure. Then, the general formulation of the tensor-based subspace learning is expressed as

Recent studies [Shashua and Levin, 2001] [Ye, 2005] [Ye et al., 2004] [Yan et al., 2005] have shown that dimensionality reduction algorithms with data encoded as high-order tensors usually outperform those with data represented as vectors, especially when the number of training samples is small. Representing images as 2D matrices instead of vectors allows correlations between both rows and columns to be exploited for subspace learning.

Generally, no closed-form solution exists for (7). Previous works [Ye et al., 2004] [Yan et al., 2005] utilized iterative procedures to search for approximate solutions. First, the projection matrices U 1 , . . . , U n are initialized arbitrarily; then each projection matrix U k is refined by fixing the other projection matrices U 1 , . . . , U k-1 , U k+1 , . . . , U n and solving the optimization problem:

where Y k i is the mode-k unfoldingmatrix of the tensor ˜Yi = Xi ×1 U 1 . . . ×k-1 U k-1 ×k+1 U k+1 . . . ×n U n and S k = ∑i̸=j Sij(Y k i - Y k j )(Y k i - Y k j )T , Sp k = ∑i̸=j S p ij (Y k i - Y k j )(Y k i -Y k j )T .

The optimization problem in (9) is still intractable, and traditionally its solution is approximated by transforming the objective function in (9) into a more tractable approximate form, namely, Ratio Trace form,

which can be directly solved with the generalized eigenvalue decomposition method. However, this distortion of the objective function leads to the computational issues as detailed in the following subsection.

2.2 Computational Issues

As the objective function in each iteration is changed from the trace ratio form (9) to the ratio trace form (10), the deduced solution can satisfy neither of the two aspects: 1) the objective function value in (7) can monotonously increase; and 2) the solution (U 1, U 2 , . . . , U n) can converge to a fixed point. In this work, we present a convergent solution procedure to the optimization problem defined in (7).

3 Solution Procedure and Convergency Proof

In this section, we first introduce our new solution procedure to the tensor-based subspace learning problems, and then give the convergency proof to the two aspects mentioned above.

As described above, there does not exist closed-form solution for the optimization problem (7), and we solve the optimization problem also in an iterative manner. For each iteration, we refine one projection matrix by fixing the others and an efficient method is proposed for this refinement. Instead of solving a ratio trace optimization problem (10) for an approximate solution, we transform the trace ratio optimization problem (9) into a trace difference optimization problem defined as

where λ is the value of objective function (7) computed from the projection matrices of the previous iteration.

Though the iterative procedure may converge to a local optimum for the optimization problem (7), it can monotonously increase the objective function value as proved later, which directly leads to its superiority over the ratio trace based optimization procedure, since the step-wise solution of the latter is unnecessarily optimal for (9).

We iteratively refine the projection matrices, and the detailed solution procedure to solve the tensor-based general subspace learning problem is listed in Algorithm 1.

3.1 Analysis of Monotonous Increase Property

Rewrite the objective function of (7) as

and then we have the theory as below:

Theorem-1. By following the terms in Algorithm-1 and Eqn. (13), we have

Proof. Denote g(U) = T r(U T (S p k - λS k )U) where

then we have

Algorithm 1 . Procedure to Tensor-based Subspace Learning

• 1: Initialization. InitializeU 1 0 , U 2 0 , . . . , U n 0 as arbitrary column orthogonal matrices.
• 2: Iterative optimization.

1. ComputeS k andS p k as in (9) based on the projection matrices U 1 t , . . . , U k-1 t and U k+1 t-1 , . . . , U n t-1 .
2. Conduct Eigenvalue Decomposition:

where vj is the eigenvector corresponding to the jth largest eigenvalue λj .

1. Reshape the projection directions for the sake of orthogonal transformation invariance:

• (b) Let S v = V V T (∑i X k i X k i T )V V T , where X k i is the mode-k unfolding of the tensor Xi;
• (c) Conduct Eigenvalue Decomposition as

1. Set the column vectors of matrix U k t as the leading eigenvectors, namely, U k t = [u1, u2, . . . , um′k ].

End If ‖U k t - U k t-1‖ < √mkm′k ε, k = 1, 2, . . . , n (ε is set to 10-4 in this work), then break.

End

• 3: Output the projectionmatrices U k =U k t , k=1, 2, . . . , n.

Moreover, fromU T U = Im′k , it is easy to prove that

From Algorithm 1, we have g(U k t ) = ∑m′k j=1 λj , and hence

Then, Tr(U k t T (S p k - λS k )U k t ) ≥ 0. As matrix S k is positive semidefinite 2 , we have

2 Though Sk may have zero eigenvalues, T r(U k t T SkU k t ) will be positive when m′k is larger than the number of the zero eigenvalues.

From theorem-1, we can conclude that the value of the objective function monotonously increases.

3.2 Proof of Convergency

To prove the convergency of the projection matrices U 1 , U 2 , . . . , U n, we need the concept of point-to-set map. The power set ℘(χ) of a set χ is the collection of all subsets of χ. A point-to-set map Ω is a function: χ → ℘(χ). In our solution procedure to tensor-based subspace learning, the map from (U k t-1|nk=1 ) to (U k t |nk=1 ) can be considered as a point-to-set map, since each U k t is invariant under any orthogonal transformation.

Strict Monotony. An algorithm is a point-to-set map Ω:χ → ℘(χ). Given an initial point x0, an algorithm generates a sequence of points via the rule that xt ∈ Ω(xt-1). Suppose J : χ → R+ is continuous, non-negative function, an algorithm is called strict monotony if 1) y ∈ Ω(x) implies that J(y) ≥ J(x), and 2) y ∈ Ω(x) and J(y) = J(x) imply that y = x.

Let set χ be the direct sum of the orthogonal matrix space Omk×m′k, that is, the data space χ = Om1×m′1 ⊕ Om2×m′2 ⊕ . . . ⊕ Omn×m′n, then the Algorithm 1 produces a point-to-set algorithm with respect to J(x) = G(U k |nk=1 ), and it can be proved to be strictly monotonic as follows.

Theorem-2. The point-to-set map from Algorithm 1 is strictly monotonic.

Proof. From theorem-1, we have G(U k t-1 |nk=1 ) ≤ G(U k t |nk=1 ), and hence the first condition for strict monotony is satisfied. For the second condition, we take U 1 as an example to prove that this condition is also satisfied. If G(U k t-1 |nk=1 ) = G(U k t |nk=1 ), then from the proof of theorem1, we have g(U 1 t-1) = g(U 1 t ) with λ = G(U k t-1|nk=1 ) and S k , Sp k computed from (U k t-1 |nk=1 ). From the proof of theorem-1, we can have that there only exists one orthogonal transformation 3 between U 1 t-1 and U 1 t . As shown in Algorithm 1, this kind of orthogonal transformation has been normalized by the reshaping step, hence we have U 1 t-1=U 1 t . Similarly, we can prove that U k t = U k t-1 for k = 1, 2, . . . , n, hence the second condition is also satisfied and the Algorithm 1 is strictly monotonic.

3 This claim is based on the assumption that there do not exist duplicated eigenvalues in (11).

Theorem-3 [Meyer, 1976]. Assume that the algorithm Ω is strictly monotonic with respect to J and it generates a sequence {xt} which lies in a compact set. If χ is normed, then ‖xt - xt-1‖ → 0.

From theorem-3, we can have the conclusion that the obtained (U k t |nk=1 ) will converge to a local optimum, since the χ is compact and with norm definition.

4 Experiments

In this section, we systematically examine the convergency properties of our proposed solution procedure to tensor-based subspace learning. We take the Marginal Fisher Analysis (MFA) as an instance of general subspace learning, since MFA has shown to be superior to many traditional subspace learning algorithms such as Linear Discriminant Analysis (LDA); more details on the MFA algorithm is referred to [Yan et al., 2007]. Then, we evaluate the classification capability of the derived low-dimensional representation from our solution procedure compared with the traditional procedure proposed in [Ye et al., 2004] and [Yan et al., 2005]. For tensor-based algorithm, the image matrix, 2r d tensor, is used as input, and the image matrix is transformed into the corresponding vector as the input of vector-based algorithms.

4.1 Data Sets

Three real-world data sets are used. One is the USPS handwritten dataset 4 of 16-by-16 images of handwritten digits with pixel values ranging between -1 and 1. The other two are the benchmark face databases, ORL and CMU PIE 5 . For the face databases, affine transform is performed on all the samples to fix the positions of the two eyes and the mouth center. The ORL database contains 400 images of 40 persons, where each image is normalized to the size of 56-by-46 pixels. The CMU PIE (Pose, Illumination, and Expression) database contains more than 40,000 facial images of 68 people. In our experiment, a subset of five near frontal poses

4 Available at: http://www-stat-class.stanford.edu/ tibs/ElemStatLearn/data.html

5 Available at http://www.face-rec.org/databases/.

(C27, C05, C29, C09 and C07) and illuminations indexed as 08 and 11 are used and normalized to the size of 32-by-32.

4.2 Monotony of Objective Function Value

In this subsection, we examine the monotony property of the objective function value from our solution procedure compared with the optimization procedure that step-wisely transforms the objective function into the ratio trace form. The USPS, ORL and PIE databases are used for this evaluation. The detailed results are shown in Figure 1. It is observed that the traditional ratio trace based procedure does not converge, while our new solution procedure guarantees the monotonous increase of the objective function value and commonly our new procedure will converge after about 4-10 iterations. Moreover, the final converged value of the objective function from our new procedure is much larger than the value of the objective function for any iteration of the ratio trace based procedure.

4.3 Convergency of the Projection Matrices

To evaluate the solution convergency property compared with the traditional ratio trace based optimization procedure, we calculate the difference norm of the projection matrices from two successive iterations and the detailed results are displayed in Figure 2. It demonstrates that the projection matrices converge after 4-10 iterations for our new solution procedure; while for the traditional procedure, heavy oscillations exist and the solution does not converge. As shown in Figure 3, the recognition rate is sensitive to the oscillations caused by the unconvergent projection matrices and the classification accuracy is degraded dramatically.

4.4 Face Recognition

In this subsection, we conduct classification experiments on the benchmark face databases. The Tensor Marginal Fisher Analysis algorithm based on our new solution procedure (TMFA TR) is compared with the traditional ratio trace based Tensor Marginal Fisher Analysis (TMFA RT), LDA, Ratio Trace based MFA (MFA RT) and Trace Ratio based MFA (MFA TR), where MFA TR means to conduct tensor-based MFA by assuming n=1. To speed up model training, PCA is conducted as a preprocess step for vector-based algorithms. The PCA dimension is set as N-Nc (N is the sample number and Nc is the class number), which is equivalent to the case for Fisherface algorithm [Belhumeur et al., 1997]. The same graph configuration with nearest neighbor k = 3 for the intrinsic graph and kp = 40 for the penalty graph is adopted for all the MFA based algorithms. Since the traditional tensor subspace learning algorithms do not converge, we terminate the process after 3 iterations.

For comparison, the classification result on the original gray-level features without dimensionality reduction is also reported as the baseline, denoted as 'w/o DR.' in the result tables. In all the experiments, the Nearest Neighbor method is used for final classification. All possible dimensions of the final low-dimensional representation are evaluated, and the best results are reported. For each database, we test various configurations of training and testing sets for the sake of statistical confidence, denoted as ′GxPy′ for which x images of each subject are randomly selected for model training and the remaining y images of each subject are used for testing. The detailed results are listed in Table 1 and 2. From these results, we can have the following observations:

1. TMFA TR mostly outperforms all the other methods concerned in this work, with only one exception for the case G5P 5 on the CMU PIE database.
2. For vector-based algorithms, the trace ratio based formulation (MFA TR) is consistently superior to the ratio trace based one (MFA RT) for subspace learning.
3. Tensor representation has the potential to improve the classification performance for both trace ratio and ratio trace formulations of subspace learning.

5 Conclusions

In this paper, a novel iterative procedure was proposed to directly optimize the objective function of general subspace learning based on tensor representation. The convergence of the projection matrices and the monotony property of the objective function value were proven. To the best of our knowledge, it is the first work to give a convergent solution for general tensor-based subspace learning.

6 Acknowledgement

The work described in this paper was supported by grants from the Research Grants Council of the Hong Kong Special Administrative Region and DTO Contract NBCHC060160 of USA.