An EM Based Probabilistic Two-Dimensional CCA with Application to Face Recognition
Mehran Safayani Seyed Hashem Ahmadi Homayun Afrabandpey Abdolreza Mirzaei



Probabilistic 2DCCA with EM 




Department of Electrical and Computer Engineering, Isfahan University of Technology, Isfahan 84156-83111, IRAN 

              Tel.: +98 31 33919063

              Fax: +98 31 33912451

              safayani@cc.iut.ac.ir(Corresponding Author)hashem.ahmadi@ec.iut.ac.irh.afraei@ec.iut.ac.ir mirzaei@cc.iut.ac.ir 






An EM Based Probabilistic Two-Dimensional CCA with Application to Face Recognition

    Mehran Safayani Seyed Hashem Ahmadi Homayun Afrabandpey Abdolreza Mirzaei 
    Received: date / Accepted: date
===================================================================================




Recently, two-dimensional canonical correlation analysis (2DCCA) has been successfully applied for image feature extraction. The method instead of concatenating the columns of the images to the one-dimensional vectors, directly works with two-dimensional image matrices. Although 2DCCA works well in different recognition tasks, it lacks a probabilistic interpretation. In this paper, we present a probabilistic framework for 2DCCA called probabilistic 2DCCA (P2DCCA) and an iterative EM based algorithm for optimizing the parameters. Experimental results on synthetic and real data demonstrate superior performance in loading factor estimation for P2DCCA compared to 2DCCA. For real data, three subsets of AR face database and also the UMIST face database confirm the robustness of the proposed algorithm in face recognition tasks with different illumination conditions, facial expressions, poses and occlusions.





§ INTRODUCTION


Although many real-world applications encounter high dimensional data, the most informative part of the data can be modeled in a low dimensional space. Moreover, processing high-dimensional data is a time consuming process and requires lots of resources. To tackle these problems, feature extraction has been used as a tool for finding a compact and meaningful data representation.

For single-mode source data, some subspace learning methods are conducted to learn more semantic description subspaces. Examples of these methods are principal component analysis (PCA) <cit.> and linear discriminant analysis (LDA). However, for observations from two sources that share some mutual information, canonical correlation analysis (CCA) <cit.> is a very popular approach for dimensionality reduction. CCA seeks a lower-dimensional space where two sets of variables are maximally correlated after projecting on it. This technique is widely used in different fields of pattern recognition, computer vision, bioinformatics, etc. <cit.>. In the CCA-based methods, it is necessary to vectorize 2D image matrices. Vectorization has three main drawbacks: (I) breaking the spatial structure of image data which may cause losing potentially useful structural information among column/rows <cit.>, (II) leading to a high-dimensional vector space and small sample size problem which in turn makes it difficult to calculate the covariance matrices <cit.> and (III) causing the covariance matrices to be very large which in turn makes the eigen-decomposition of such large matrices very time-consuming.

To overcome these drawbacks, in 2007 two-dimensional CCA (2DCCA) was introduced by Lee and Choi <cit.> which computes CCA directions based on 2D image matrices. The proposed 2DCCA overcomes the curse of dimensionality and significantly reduces the computational cost, by directly working with 2D images instead of reshaping them into 1D vectors. In <cit.>, higher recognition accuracies were reported using 2DCCA compared to CCA using two face databases and the time complexity has been improved.

However, an associated probabilistic model for observed data was notably absent form these feature extraction methods. A probabilistic feature extraction algorithm could be intuitively appealing for so many reasons <cit.>. To bridge the gap, in 1999, Tipping and Bishop proposed probabilistic PCA <cit.> based on a latent variable model known as factor analysis (FA) <cit.>. The proposed PPCA was then used as a framework for many other new formulations for PCA <cit.>.
Also, there have been some probabilistic models proposed for LDA <cit.>. In 2005, Bach and Jordan <cit.> also proposed a probabilistic interpretation of CCA and estimate the parameters of their proposed model using both maximum likelihood and expectation maximization. Recently, many inspiring research proceeded in the 1D CCA domain, including kernel based, semiparametric and nonparametric methods <cit.>, but in the 2D CCA domain, we feel that more work is required. To bridge the gap, a probabilistic model of 2DCCA was introduced by Safayani et al. in 2011 <cit.>. They showed that the maximum likelihood estimation of parameters, leads to the two dimensional canonical correlation directions. However, they didn't propose an EM based solution for their model. EM does not require the explicit eigen-decomposition of covariance matrices. Moreover, using EM it is possible to handle models with incomplete data such as mixture models where the cluster labels are the missing values <cit.>.

In this paper, we present a probabilistic interpretation of 2DCCA, referred to as P2DCCA, together with an EM based solution to estimate the parameters of the model. The proposed model can handle the small sample size problem effectively.

The rest of the paper is organized as follows: Section <ref> briefly reviews some related algorithms such CCA, PCCA and 2DCCA which are necessary to understand how the proposed algorithms work. The proposed P2DCCA model is introduced in Section <ref>. In Section <ref>, some experiments on synthetic data and several face databases are given to evaluate performance of the proposed algorithm; finally, the paper is concluded in Section <ref>.


§ BACKGROUND




 §.§ Canonical Correlation Analysis (CCA)

Imagine that we are given two sets of random vectors t_1 and t_2 where t_1,n∈R^D_1 and t_2,n∈R^D_2 for n∈1,2,...,N are realizations of the corresponding random vectors, respectively. CCA seeks transformation vectors w_1∈R^D_1 and w_2∈R^D_2 such that correlation between w_1^Tt_1 and w_2^Tt_2 are maximized. The correlation between w_1^Tt_1 and w_2^Tt_2 can be formulated as

    ρ = cov(w_1^Tt_1,w_2^Tt_2)/√(var(w_1^Tt_1)var(w_2^Tt_2))=w_1^TΣ_12w_2/√((w_1^TΣ_11w_1)(w_2^TΣ_22w_2))

where Σ_ij=1/N∑_n=1^N(t_i,n-μ_i)(t_j,n-μ_j)^T for i,j ∈1,2 is the cross-covariance matrix of t_1 and t_2 and μ_i=1/N∑_n=1^Nt_i,n for i ∈{1,2} denotes the mean vector of t_i.

Then, the objective function for CCA can be written as:

_
w_1,w_2
w_1^TΣ_12w_2   
s.t.   w_1^TΣ_11w_2=1
    w_2^TΣ_22w_2=1

Optimizing such a constrained maximization problem with respect to w_1 and w_2 leads to the following generalized eigenvalue problem:

    [    0 Σ_12; Σ_21    0 ][ w_1; w_2 ]
    =λ[ Σ_11    0;    0 Σ_22 ][ w_1; w_2 ]

By solving equation (<ref>), w_1 and w_2 that maximize the correlation between the projected data can be found.



 §.§ Probabilistic CCA


The generative model introduced by Bach and Jordan for CCA is as follows:

    t_i=W_iz+μ_i+ϵ_i     i∈{1,2}

In this model, W_i∈^D_i× d   i∈{1,2} are linear projections that map two sets of high dimensional observed random vectors t_i∈^D_i, i∈{1,2} to a set of lower dimensional latent vectors z∈^d. Here μ_i is the mean vector for x_i and ϵ_i is the error term which is assumed to follow a multivariate Gaussian distributions with zero mean and inverse covariance matrix Ψ_i. Bach and Jordan proved that the maximum likelihood estimation for the parameters of this model would lead to the canonical directions. The maximum likelihood estimates of the projection matrices are given by:

    W_i=Σ_iiU_idM_i,   i∈{1,2}

where Σ_ii is the sample covariance matrix, U_id are the first d canonical directions and M_i∈^d× d where i∈{1,2} are arbitrary matrices such that M_1M_2 = P_d, where P_d is the diagonal matrix of the first d canonical correlations. Figure <ref> shows a graphical representation of the model.





 §.§ Two-Dimensional CCA


Two-dimensional CCA (2DCCA) was proposed to tackle the problem of vectorizing data in CCA. For each random matrix T_1 and T_2, 2DCCA introduces left transforms u_i and right transforms v_i where i∈{1,2}. After the projection, data would have the form u_i^TT_iv_i. 2DCCA finds these left and right transforms in a way to maximize the correlation between projected data. Therefore, the objective function of 2DCCA can be formulated as:

   _
u_1,u_2,v_1,v_2
cov(u_1^TT_1v_1,u_2^TT_2v_2)   
   s.t.   var(u_1^TT_1v_1)=1
       var(u_2^TT_2v_2)=1,   

u_1 and u_2 can be obtained by solving the generalized eigenvalue problem (<ref>) with fixed v_1 and v_2:

[      0 Σ_12^r; Σ_21^r      0 ][ u_1; u_2 ]=λ[ Σ_11^r      0;      0 Σ_22^r ][ u_1; u_2 ]

In a similar way, given u_1 and u_2, right transforms v_1 and v_2 can be found by solving

[      0 Σ_12^l; Σ_21^l      0 ][ v_1; v_2 ]=λ[ Σ_11^l      0;      0 Σ_22^l ][ v_1; v_2 ]

where Σ_i,j^r, i,j∈{1,2} is the cross-covariance matrix between T_i and T_j and Σ_ii^r, i∈{1,2} is the auto-covariance matrix of T_i defined as follows, respectively:

Σ_ij^r=1/N∑_n=1^N(T_i,n-μ_i)v_iv_j^T(T_j,n-μ_j)^T,    i,j∈{1,2}


Σ_ij^l=1/N∑_n=1^N(T_i,n-μ_i)u_iu_j^T(T_j,n-μ_j)^T,    i,j∈{1,2}

where T_i,n for n∈{1,...,N} is the realization of the random matrix T_i and μ_i=1/N∑_n=1^NT_i,n is the corresponding matrix of mean values.

Left transforms (l_x and l_y) and right transforms (r_x and r_y) are obtained by iteratively solving equations (<ref>) and (<ref>), until convergence. The eigenvectors associated with d_1 largest eigenvalues in (<ref>) determine left transform matrices U_1 and U_2 and the eigenvectors associated with d_2 largest eigenvalues in (<ref>) determine right transform matrices V_1 and V_2. Using these transform matrices it is possible to project data from a high dimensional space to a new lower dimensional feature space.



§ PROBABILISTIC TWO DIMENSIONAL CCA (P2DCCA)


In this section, we propose probabilistic two-dimensional CCA and an EM-based solution for finding the parameters of the model. In our model, observed data are modeled as two-dimensional matrices as follows:

T_i=U_iZV_i^T+μ_i+Ξ_i   i∈{1,2}

where T_i∈^m_i× n_i for i∈{1,2} are observed matrices and Z∈^m× n is the latent matrix. U_i∈^m_i× m, V_i∈^n_i× n are projection matrices, μ_i is the mean matrix of the observed data and Ξ_i is the residual matrix. Based on this definition, parameters of the model are {U_i,V_i}_i=1^2 and the parameters of the distribution of Ξ_i.

Let D_i={T_i;n}_n=1^N where i∈{1,2} be a set containing N observed data matrices and {Z_n}_n=1^N be the corresponding latent variable set. Then the complete data would be (T_1,n,T_2,n,Z_n) and the log likelihood of the complete data can be written as

L_c=∑_n=1^Nlog   p(T_1,n,T_2,n,Z_n)

To estimate the parameters, first we must calculate expectation of the log-likelihood and then take the derivative of the expected log-likelihood with respect to each parameter. Unfortunately there is no closed-form solution for computing the projection matrices {U_i,V_i}_i=1^2 simultaneously. Inspired by <cit.>, a decoupled probabilistic model is employed to obtain projection matrices separately using an alternating optimization procedure. In such a model, we first assume that the value of one set of projection matrices, e.g. the right projections {V_i}_i=1^2, is known. Then observations are projected to the corresponding latent spaces. The projection procedure is a probabilistic one that introduced in section (<ref>). By doing this, the left probabilistic model is defined as

T_i^l=U_iZ^l+Ξ_i^l,   i=1,2

where Z^l is the left model latent matrix, μ_i^l is the mean matrix of the left projected observations and Ξ_i^l is the noise source for left probabilistic model where columns of the noise matrix follow a normal distribution with zero mean and covariance matrix Ψ_i^l. By such definition, parameter set for the left probabilistic model would be Θ^l={U_i,Ψ_i^l}_i=1^2 which can be estimated using expectation maximization procedure. The estimation procedure is explained later in this section.

In a similar procedure and parallel to the left probabilistic model, for the right probabilistic model, we assume that the left projection matrices, i.e. {U_i}_i=1^2, are known. Then observations are projected over the corresponding latent spaces,hence the right probabilistic model is defined as

T_i^r=V_iZ^r+Ξ_i^r,   i=1,2

Similar to the left probabilistic model, Z^r, μ_i^r and Ξ_i^r are defined for the right model where in this model the noise source have N(0,Ψ_i^r) distribution. The parameter set for the right probabilistic model would be Θ^r={V_i,Ψ_i^r}_i=1^2. Using these definitions, the decoupled predictive density p(T_1,T_2,Z) could be defined as

p(T_1,T_2,Z)∝p(T_1^l,T_2^l,Z^l)p(T_1^r,T_2^r,Z^r)

Now we can rewrite the log likelihood of equation (<ref>) as

   L_c=∑_n=1^Nlog(p(T_1,n^l,T_2,n^l,Z_n^l)p(T_1,n^r,T_2,n^r,Z_n^r))   
     =∑_n=1^Nlog  p(T_1,n^l,T_2,n^l,Z_n^l)+∑_n=1^Nlog  p(T_1,n^r,T_2,n^r,Z_n^r).   

To apply the EM algorithm to the decoupled probabilistic model, in E-step expectation of log likelihood for the left probabilistic model and the right probabilistic model is computed, separately. Then each of the expected log likelihood is maximized with respect to its parameters. In the following subsections we describe how to optimize left and right probabilistic model respectively.




 §.§ Optimizing the left probabilistic model


Let t_i,j^l be the j^th column vector of T_i^l∈^m_i× n. By assuming columns of T^l to be independent of each other, the distribution of T_i^l is defined as

p(T_i^l)=∏_j=1^np(t_i,j^l),     i ∈{1,2}

We also consider z_j^l∈^m× 1 as the j^th column vector of Z^l which has normal distribution of N(0,I) and also in the same way μ_i,j^l is the j^th column vector of μ_j^l. Based on equation (<ref>) and the distribution considered for Z^l and Ξ_i^l, it can be concluded that

p(t_i,j^l)∼N(μ_i,j^l,U_iU_i^T+Ψ_i^l),     i=1,2

Suppose τ_n,j^l= [(t_1,n,j^l)^T  (t_2,n,j^l)^T]^T∈^(m_1+m_2)× 1, V=[U_1^T  U_2^T]^T∈^(m_1+m_2)× m, m_j^l=[(μ_1,j^l)^T  (μ_2,j^l)^T]^T∈^(m_1+m_2)× 1 and Ψ^l = (
[ Ψ_1^l     0;     0 Ψ_2^l ])
for the left probabilistic model, where t_i,n,j^l refers to the j^th column vector of n^th image in the i^th observation set where i∈{1,2}. Therefore, distributions of p(τ_j^l) can be obtained as follows:

p(τ_j^l)∼N(m_j^l,Σ^l),    j∈[1,n]

where Σ^l=UU^T+Ψ^l and we assume Σ^l > 0.

Based on (<ref>) we can write:

p(T_1,n^l,T_2,n^l,Z_n^l)=∏_j=1^np(τ_n,j^l,z_n,j^l)=∏_j=1^np(τ_n,j^l|z_n,j^l)p(z_n,j^l)

To apply the EM algorithm to the decoupled probabilistic model, for each of the probabilistic models expectation of the log likelihood function is calculated in the E-step E(L_c^l) where the detail is given in the appendix, and then maximization step (M-step) is done by maximizing E(L_c^l) with respect to V and Ψ^l. By doing so, the values of the parameters are estimated as

U_t+1=Σ^l(Ψ_t^l)^-1U_t(M_t^l)^-1/(M_t^l)^-1+(M_t^l)^-1U_t^l(Ψ_t^l)^-1Σ^l(Ψ_t^l)^-1U_t(M_t^l)^-1


Ψ_t+1^l=(
[ (Σ^l-Σ^l(Ψ_t^l)^-1U_t(M_t^l)^-1U_t+1^T)_11                                          0;                                          0 (Σ^l-Σ^l(Ψ_t^l)^-1U_t(M_t^l)^-1U_t+1^T)_22 ]
)

where M^l=I+U^T(Ψ^l)^-1U and A_t shows the value of parameter A in iteration t and Σ^l is the sample covariance matrix of observed data for the left probabilistic model , i.e.

Σ^l=1/N∑_n=1^N[(T_1,n^l-μ_1^l)^T (T_2,n^l-μ_2^l)^T]^T[(T_1,n^l-μ_1^l)^T (T_2,n^l-μ_2^l)^T]



 §.§ Optimizing the right probabilistic model


In the manner similar to optimization of the left probabilistic model we have:

p(T_i^r)=∏_j=1^mp(t_i,j^r),     i ∈{1,2}

where t_i,j^r is the j^th column vector of T_i^r∈^n_i× m. Then p(t_i,j^r) is computed as:

p(t_i,j^r)∼N(μ_i,j^r,V_iV_i^T+Ψ_i^r),     i∈{1,2}

Let τ_n,j^r= [(t_1,n,j^r)^T  (t_2,n,j^r)^T]^T∈^(n_1+n_2)× 1, V=[V_1^T  V_2^T]^T∈^(n_1+n_2)× n, m_j^r=[(μ_1,j^r)^T  (μ_2,j^r)^T]^T∈^(n_1+n_2)× 1 and Ψ^r = (
[ Ψ_1^r     0;     0 Ψ_2^r ]), where t_i,n,j^r refers to the j^th column vector of n^th image in the i^th observation set. Then p(τ_j^r) and p(T_1,n^r,T_2,n^r,Z_n^r) are obtained as follows:

p(τ_j^r)∼N(m_j^r,Σ^r),    j∈[1,m]


p(T_1,n^r,T_2,n^r,Z_n^r)=∏_j=1^mp(τ_n,j^r,z_n,j^r)=∏_j=1^mp(τ_n,j^r|z_n,j^r)p(z_n,j^r)

where Σ^r=VV^T+Ψ^r>0. Given the details in appendix (A), after computing E(L_c^r) in the E-step, the parameters V and Ψ^r are computed by maximizing  the likelihood in the M-step. So, we have:

V_t+1=Σ^r(Ψ_t^r)^-1V_t(M_t^r)^-1/(M_t^r)^-1+(M_t^r)^-1V_t^r(Ψ_t^r)^-1Σ^r(Ψ_t^r)^-1V_t(M_t^r)^-1


Ψ_t+1^r=(
[ (Σ^r-Σ^r(Ψ_t^r)^-1V_t(M_t^r)^-1V_t+1^T)_11                                          0;                                          0 (Σ^r-Σ^r(Ψ_t^r)^-1V_t(M_t^r)^-1V_t+1^T)_22 ]
)

where M^r=I+V^T(Ψ^r)^-1V and Σ^r is computed as follows:

Σ^r=1/N∑_n=1^N[(T_1,n^r-μ_1^r)^T (T_2,n^r-μ_2^r)^T]^T[(T_1,n^r-μ_1^r)^T (T_2,n^r-μ_2^r)^T]





 §.§ Probabilistic projection and dimension reduction

We can project the observation matrices into the latent space using the standard projection matrices ,i.e., {U_1,U_2,V_1,V_2}. However, as described in <cit.>, it is more natural to use probabilistic projections. In this regard, we represent each projected observation matrix, T_i, by mean of distribution of corresponding latent space, i.e., E(Z|T_i). For the left model, it can be shown that

E(Z^l|T_1,T_2)=(M^r)^-1[
[ V_1^T V_2^T;       ]
][
[ Ψ_1   0;   0 Ψ_2 ]
]^-1[
[ T_1-μ_1; T_2-μ_2 ]
]


E(Z^l|T_1) and E(Z^l|T_2) are obtained by marginalizing (<ref>) over T_2 and T_1 respectively. So we have
   
E(Z^l|T_i)=(M^r)^-1V_i^T(Ψ_i^r)^-1(T_i-μ_i),     i∈{1,2}

 Similarly for the right model we have
   
E(Z^r|T_i)=(M^l)^-1U_i^T(Ψ_i^l)^-1(T_i-μ_i),     i∈{1,2}

The procedure for dimension reduction is given by sequential projection in left and right models as

E(Z|T_i)=(M^l)^-1U_i^T(Ψ_i^l)^-1(T_i-μ_i)((M^r)^-1V_i^T(Ψ_i^r)^-1)^T

The P2DCCA algorithm is summarized in Figure <ref>. The proposed P2DCCA model benefits the ability to extend to other methods such as Mixtures of P2DCCA, Bayesian P2DCCA and also to robust P2DCCA.





§ EXPERIMENTAL RESULTS


We evaluated our algorithm on both synthetic and real data. In synthetic data part we verified our implementation of P2DCCA algorithm by simple synthetic data and randomly generated projected matrices. we also compared our algorithm with 2DCCA method in projection matrices estimation . For real part evaluation, the proposed P2DCCA method was used for face recognition on two well-known face image databases (AR <cit.> and UMIST <cit.>). The AR database is divided into three subsets for evaluating the performance of the system in regard to different illumination, expression and occlusion conditions. The UMIST database is used to obtain the performance in dealing with pose variation.




 §.§ Experiments on synthetic data

 In this section, we aim to verify our implementation of the proposed method with the simplest possible scenario. So we generate some synthetic data and projection matrices. Then estimate the projection matrices using our method and compare them with the true ones. We know that P2DCCA estimations are up to rotation, hence to simplify comparison we assume Z to be 1 × 1 dimension. We set dimensions of T_1 and T_2 to 5 × 5. Then we generate 1000 samples of Z from a normal distribution with zero mean and unit variance, also we randomly generate the elements of U_i∈ ^5, V_i∈ ^5, i∈{1,2}  using a uniform distribution in [0, 1] interval and consider them as the ground truth projection matrices, then T_i will obtained using (<ref>). In this equation, each element of the residual matrices sampled from a gaussian distribution with 0 mean and σ_i^2 variance. Having the synthetic data we run P2DCCA algorithm as discussed in Figure <ref> and calculate U_i and V_i. We also run 2DCCA and obtain the corresponding projection matrices. Then we compare obtained matrices by these two algorithms with the ground truth projection matrices. To cancel the scale factors we divide each transform by its norm before comparison. Euclidean distance is utilized to compare the normalized transforms. Figure <ref> shows the results. It is obvious that in the worst-case the distance value becomes one and in the best-case it is zero. As it is depicted in this figure, P2DCCA estimation of U_i and V_i for i∈{1,2} are much closer to the ground truth compared to the those obtained by 2DCCA. In this experiment we used 1000 generated samples and set σ_i = 0.1. To examine the effect of these selections, we repeat our experiment with different values of sample numbers and also noise variances. Figure <ref> demonstrates the results. As it can be observed from this figure the P2DCCA estimation in all cases is much closer to the ground truth compared to the 2DCCA method.









 §.§ Experiments on the AR database


The AR face database contains over 4,000 color face images including frontal views of faces with different facial expressions, illumination conditions and occlusions. For most individuals, there are two sessions of images which were taken in two different time periods. Each session contains 13 images. In our experiments, we used the first session, because some individuals do not have the second session of images. We collected 1310 face images of 131 people (72 male and 59 female). For each person, there are 10 different face images in our collected images: three with different illumination conditions; three with different expressions; three with occlusions and the remaining images are those with neutral expression and no occlusion which are known as reference images in our experiments. To examine the performance of the proposed methods in different conditions, we partitioned the collected images into three subsets known as AR-1, AR-2 and AR-3. For each individual, AR-1 contains four images, three of which are images with different lighting conditions and the remaining one is the reference image. AR-2 is used to test performance of the algorithms when there exist expression variation. AR-2 involves four images per individuals: three images have different expressions and the last one is the reference image. AR-3 is prepared to test performance in the presence of occlusion. Again, this subset contains four images per individuals where three images were taken with glasses and the last one is the reference image. Figure <ref> shows exemplary face images of a man and a woman in AR-1, AR-2 and AR-3, respectively. Image are gray scaled, resized and then normalized to 50 × 50 pixels.

We compared the performance of the proposed method with a range of different supervised and unsupervised dimensionality reduction algorithms and different versions of them including PCA, LDA, CCA, PPCA, PCCA, 2DPCA, 2DLDA and 2DCCA. Both PCA and LDA based methods work with one set of data. Also LDA based methods are supervised while PCA and CCA based algorithms are unsupervised.
The task here is to investigate how well different algorithms can relate face images with varying illumination conditions, expressions and occlusion, in correspondence to the reference face images. The CCA, PCA, LDA, PPCA, PCCA, 2DCCA, 2DPCA, 2DLDA and P2DCCA are used to extract features from facial images and then a 1-NN classifier is employed for classification. Note that based on output type of each of the algorithms (vector or matrix), for 2DCCA, 2DPCA, 2DLDA and P2DCCA, Frobenius distance is used to calculate the distance between two feature matrices, while for CCA, PCA, LDA, PPCA and PCCA the common Euclidean distance measure is adopted. Furthermore, it should be noted here that since PCCA suffers the small sample size problem, implementing it using formulas introduced in <cit.> caused the covariance matrices to be singular. To solve the problem, we did dimension reduction using PCA  before implementing the algorithm.


To evaluate the recognition accuracy, we used three-fold cross-validation. As it is evident, CCA based algorithms need two sets of images for training, where in this paper the training sets are called left training set and Right training set. To form the training sets, e.g., for AR-1, neutral images (images with no illumination) are considered as Left training set, while to form the Right training set, one of the three images of each individual with different illumination conditions is selected randomly. The other two images are considered as test images. This procedure is repeated for three times, where each time a different image among the three images is selected for the Right training set, while neutral images are always used to form the Left training set. To test the performance of each algorithm, all images of both right and left training sets are projected on the new feature spaces using their corresponding transforms. Also, each of the test images is projected on both feature spaces, so we have two projection for every test image. Then we calculate the distance between each of the two projected test images and projected training images. The label of the training image with the nearest projection to any of the two test image projection determine the final class of the test image. This procedure iterates until we find the final class for all images in the test set. These final classes are compared to the real classes of the images and the recognition accuracy of each algorithm is calculated. Finally, the average recognition rate of the three round experiments is recorded as the final recognition accuracy.


  Since PCA and LDA based algorithms work with one set of data, to have a fair comparison we used two images as training and the other two images as the test data, where neutral images are always in the training data together with one of the other images with different illuminations in each iterations. Again the process is repeated three times and the final accuracy is the average of the three runs.  Figure <ref> shows the test process for P2DCCA.


Train and test procedure for AR2 and AR3 subsets are similar to AR1 and Table <ref> through Table <ref> demonstrate the recognition accuracy of evaluated algorithms for the experiments conducted on AR-1, AR-2 and AR-3, respectively. In these tables, d is the dimension of the reduced feature space. Note that output for two dimensional algorithms (2DCCA, 2DPCA, 2DLDA, P2DCCA and MP2DCCA) is of matrix type with dimension d×d, while for CCA, PCA, LDA, PPCA and PCCA methods output is a vector of dimension d. In these results we see that P2DCCA get the best performance among all tested methods.  We see about 10% improvement of recognition rate for P2DCCA over 2DCCA in AR-1 and AR-3, and about 3% improvement in AR-2.







Figure <ref> shows how the log-likelihoods of the left probabilistic model and the right probabilistic model of P2DCCA improve with each iteration. As it can be seen in the figure, both left and right models converge.



It should be noted that it is very common that in the algorithms for optimizing row and column projections only one iteration with the iterative algorithm is performed<cit.>. Therefore in all the experiments we use one iteration of the algorithm, i.e., T_max=1. This significantly reduces computational cost of the algorithm. We tried more iterations and got no significant improvement in the recognition rate. Figure <ref> shows the results for 1 to 5 iterations for AR1, AR2 and AR3. This results support the idea of choosing T_max=1.




 §.§ Experiments on the UMIST database

The UMIST face database also known as Sheffield face database <cit.> consists of 564 images of 20 subjects. Subjects have different races, sexes, and appearances. For each subject, there are images with different poses from profile to frontal view. Images have 256 grey levels with resolution of 220×220 pixels. In our experiment, 360 images with 18 samples per subject are used to examine performance of different algorithms when face orientation varies significantly. Figure <ref> shows 18 images of one subject.



We select frontal image as well as seven other randomly selected images for training set and the remaining images for the test set. In the training phase of CCA based methods, frontal image always is selected as the left training image and one of the seven other images as the right training image. 1-NN classifier is used for classification. This procedure is repeated for twenty times, and the average recognition rates of algorithms are reported. Table <ref> shows the recognition accuracy of evaluated algorithms for the experiments conducted on UMIST. In this test while P2DCCA achieved slightly better performance compared to 2DCCA, it is not the best. In fact LDA achieved the best performance. In the AR test, LDA had 2 images per class for training, but here it has 8 images per class which leaded to best performance for this supervised method. Ignoring LDA, we see that P2DCCA performance is higher than other methods.



Since the UMIST face dataset contains 20 subjects to be discriminated, the LDA features is limited to 19. To be able to compare the results of LDA with that of the other algorithms, we showed the results for d=5, d=10 and d=15 for LDA and larger values of d for other algorithms.



 §.§ Evaluation of the Experimental Results

The above experiments showed that the accuracy of P2DCCA is consistently better than other CCA based methods, i.e. CCA, PCCA and 2DCCA. But, a question sill remains: Are these differences statistically significant?. In this section we answered the question by evaluating the experimental results using independent-samples T-test (or independent t-test, for short). In this section and also the next section, we only considered CCA based algorithms including CCA, PCCA, 2DCCA and P2DCCA  since the goal of this paper is to compare the functionality of the newly proposed CCA based method with that of the other CCA based algorithms. The desired significance level is 0.05 and the null hypothesis is that there is no significant difference between recognition rates of P2DCCA and CCA, PCCA and 2DCCA, respectively. We reject the null hypothesis whenever the resulted ρ-value becomes lower than 0.05 and in this case the result can be considered statistically significant. It is necessary to note that to run the t-test in each dataset, for each algorithm we considered the highest recognition rate. Table <ref> shows the ρ-value of the test. As can be seen from this table, P2DCCA significantly outperforms other algorithms and the null-hypothesis has been rejected in all cases.








 §.§ Computational Complexity

This section compares the computational cost of the algorithms. To compare the time complexity of the algorithms, we consider input images of size m × m where we want to reduce their dimension to d × d in case of two-dimensional algorithms. For CCA and PCCA, vectorization caused the input data to have dimension m^2× 1 and the output to have dimension d × 1. However, for simplicity, d is considered to be equal to m, i.e. d=m in our analysis. Table <ref> shows the computational complexity of the algorithms. In this table, N is the number of random samples in the dataset. It should be noted that there are two types of iteration in the corresponding methods; one is the iteration necessary for the convergence of the EM part of the algorithm and the other is the iteration for alternating the optimization procedure between left and right model. We show the former by t and the latter by r in the table. However, r=1 in our experiments.







§ CONCLUSION

This paper proposed a probabilistic model for two dimensional CCA termed as P2DCCA together with an EM-based solution to estimate the parameters of the model. Experimental results demonstrated the functionality of the proposed method. The proposed P2DCCA has many advantages over 2DCCA where the most significant advantage is its ability to extend to a mixture of P2DCCA model. It may also be possible to develop a probabilistic Bayesian model for P2DCCA and gaining the benefits of a Bayesian model. These are our future works.





§ 

As it is mentioned in Section <ref>, each column of the latent matrix Z has the distribution N(0,I). Furthermore, based on (<ref>) and (<ref>) we can write:


   p(τ_n,j^l|z_n,j^l)∼N(Uz_n,j^l+m_j^l,Ψ^l)   





   p(τ_n,j^r|z_n,j^r)∼N(Vz_n,j^r+m_j^r,Ψ^r).   


Now from (<ref>) and (<ref>) and the distribution of columns of Z, we have


   p(T_1,n^l,T_2,n^l,Z_n^l)=∏_j=1^n[(2π)^-m_1+m_2/2|Ψ^l|^-1/2   
                   exp(-1/2(τ_n,j^l-Uz_n,j^l-m_j^l)^T(Ψ^l)^-1(τ_n,j^l-Uz_n,j^l-m_j^l))   
                   (2π)^-m/2exp(-1/2(z_n,j^l)^T(z_n,j^l))],   



   p(T_1,n^r,T_2,n^r,Z_n^r)=∏_j=1^m[(2π)^-n_1+n_2/2|Ψ^r|^-1/2   
                   exp(-1/2(τ_n,j^r-Vz_n,j^r-m_j^r)^T(Ψ^r)^-1(τ_n,j^r-Vz_n,j^r-m_j^r))   
                   (2π)^-n/2exp(-1/2(z_n,j^r)^T(z_n,j^r))],   


where |A| denotes the determinant of matrix A.

In the E-step, expectation of log likelihood for each of the probabilistic models is calculated as:


   E(L_c^l)=∑_n=1^N∑_j=1^n[-m_1+m_2/2log(2π)-1/2log|Ψ^l|   
   -1/2tr{(Ψ^l)^-1(τ_n,j^l)(τ_n,j^l)^T}+tr{(E[(z_n,j^l)]-m_j^l)^TU^T(Ψ^l)^-1τ_n,j^l}   
   -1/2tr{(E[z_n,j^l]-m_j^l)^TU^T(Ψ^l)^-1U(E[z_n,j^l]-m_j^l)}-m/2log(2π)   
   -1/2tr{E[(z_n,j^l)(z_n,j^l)^T]}]   



   E(L_c^r)=∑_n=1^N∑_j=1^m[-n_1+n_2/2log(2π)-1/2log|Ψ^r|   
   -1/2tr{(Ψ^r)^-1(τ_n,j^r)(τ_n,j^r)^T}+tr{(E[(z_n,j^r)]-m_j^r)^TV^T(Ψ^r)^-1τ_n,j^r}   
   -1/2tr{(E[z_n,j^r]-m_j^r)^TV^T(Ψ^r)^-1V(E[z_n,j^r]-m_j^r)}-n/2log(2π)   
   -1/2tr{E[(z_n,j^r)(z_n,j^r)^T]}]   



   E[z_n,j^l]=(M^l)^-1U^T(Ψ^l)^-1(τ_n,j^l-m_j^l)   


   E[(z_n,j^l)(z_n,j^l)^T]=(M^l)^-1+E[(z_n,j^l)]E[(z_n,j^l)]^T   


   E[z_n,j^r]=(M^r)^-1V^T(Ψ^r)^-1(τ_n,j^r-m_j^r)   


   E[(z_n,j^r)(z_n,j^r)^T]=(M^r)^-1+E[(z_n,j^r)]E[(z_n,j^r)]^T.   


By obtaining formulas of E(L_c^l) and E(L_c^r), the M-step is done by derivation of each expected log-likelihoods.

