Bayesian Extreme Components Analysis

1 Introduction

Principal components analysis (PCA) is perhaps the most widely used algorithm in the fields of statistics, machine learning and data mining. It enjoys a number of desirable properties such as optimal reconstruction of the original signal in the L2 norm and retaining the maximal variance directions in the data. Despite this, there are many examples where it is not the principal components (PCs) that convey the important information, but rather the directions of small variance, or minor components (MCs). One can think of minor components as properties that are conserved in the data, i.e. constraints.

As a motivating example we consider the landmark measurements on the edge of a mosquito wing 1 , see Figure 1(a). In our experiments we have sub-sampled a number of these landmarks located on the edge of the wing and computed the relative distances between these landmarks and some of their nearest neighbors. We have used distances as features because they remove translational and rotational degrees of freedom (for more details see section 5.1.)

In Figure 1(b) we have plotted the deformation of the geometry of the landmark configuration as we change the coefficient for the principal component or the minor component of the data. One can observe that the first PC corresponds to a shift of the landmarks over the edge of the wing, but keeping the wing shape mostly invariant, indicating that the location where the veins of the wing intersect with the boundary is highly variable across mosquitos. Looking at the minor component of the data we see that it (when varied) would drastically change the shape of the wing, in particular the part where the wing is attached to the mosquito's body. Since minor components express variability which is absent in the data it implies that this type of shape change is highly unlikely in the mosquito population. We argue that these "conservation laws of biological evolution" are of more scientific interest than the directions of high variability.

1 Obtained from http://life.bio.sunysb.edu/morph/index.html

Figure 1: Mosquito Wing Landmarks

More generally, for an arbitrary dataset we would like to reliably determine which linear subspace constitutes an optimal description of the data. A statistical technique called eXtreme Components Analysis (XCA) was introduced in [Welling et al., 2003] to determine the optimal combination of principal and minor components automatically from data. In a statistical sense, it is very difficult to reliably estimate minor components from data if the number of data instances is relatively small compared to the number of attributes. This effect is illustrated in Figure 2 where we plot the sample spectrum computed from a multivariate normal distribution with unit variance in all directions. Even though this data should clearly have no preference for either principal or minor components, sample fluctuations always create artificial low variance directions which correspond to under-sampled directions in space. Since an eigenvalue decomposition searches for these directions it is highly prone to over-fitting to this type of sampling noise.

The contribution of this paper is to largely resolve this issue by introducing a conjugate prior to the parameters of the XCA model which in turn regularizes the eigenvalue decomposition. The effect of this is that minor components will only be incorporated in the model if they represent genuine constraints in the data and not merely under-sampled directions of space.

In our empirical evaluation we show that Bayesian XCA performs at least as good, and often better than 1) XCA based on frequentist shrinkage estimates of the sample covariance, 2) Bayesian PCA and 3) plain vanilla XCA. These results still hold true for an extension to mixtures of Bayesian XCA models that we discuss as well.

2 Extreme Components Analysis

XCA models a data cloud probabilistically by starting with an isotropic Gaussian pdf and either stretching in certain directions (the PCs) or contracting it in certain direction (the MCs). Denote these directions with wi (organized as rows of W ).

where z1 and z2 represent independent subspaces and the rows of V form an orthonormal basis in the orthogonal complement of the space spanned by W. The distribution of XCA

is normal with inverse covariance matrix given by

We can decompose W T using an SVD as follows,

where U ∈ RD×d,R ∈ Rd×d satisfying U T U = Id and RT R = Id. L is diagonal matrix with elements {li}. One can show that the variance of data in the direction ui (ith column of U) is given by li and in the remaining directions by v. Data is elongated in directions where li > v and contracted in directions where li < v. When all the li's are larger than v, XCA is equivalent to probabilistic PCA (PPCA) [Roweis, 1998; Tipping and Bishop, 1999b]. On the other hand, when all the li's are smaller than v, XCA is equivalent to probabilistic MCA (PMCA)[Williams and Agakov, 2002].

3 Bayesian XCA

Since XCA is more flexible than PCA and MCA, it will always have equal or higher likelihood on the training set. However, XCA can be sensitive to overfitting. The selection of PCs and MCs depends on the estimation of the sample spectrum, but both the variance as well as the bias of ordered eigenvalues of the sample covariance matrix increase as the ratio of the number of data instances to the number of attributes decreases. In the extreme case where the number of data instances is smaller than the number of attributes, the smallest eigenvalues are 0 and XCA will always pick these minor components. This results in a positive infinite log-likelihood of training data but a negative infinite loglikelihood of test data. Thus, regularization is necessary for XCA when the size of the training dataset is small compared to the number of attributes.

3.1 A Prior for XCA:

In this paper, we use a prior for the parameters of XCA to regularize the estimation of the eigenvalues. When the number of data cases increases, the effect of the prior will automatically diminish. We are aware of two approaches to Bayesian formulations of PPCA [Bishop, 1999; Minka, 2000]. Unlike [Bishop, 1999], the method of [Minka, 2000] can be extended to XCA as we will describe below in more detail 2 .

The probability of the dataset D given the covariance matrix C and mean m is under the XCA model is,

(4) where S = ∑n(xn -m)(xn -m)T , C-1 = 1 v (I -UU T )+ UL-1U T , N is the number of data points and D the dimensionality of data.

We use a conjugate prior for (U, L, R, v, m) controlled by two hyper-parameters α, β, similar to those used in [Minka,

2 We thank T. Minka for his prompt responses to our questions about his method. The prior in this paper is similar to that in an updated version of [Minka, 2000].

2000].

The prior for m is constant in an area large enough for the problem.This distribution factors into separate terms for (U, L, R, v):

The priors for li and v are χ-2 distributed with mode β while for other irrelevant parameters they are uninformative. The mean of li and v are approximately equal for large α. The width of the peak is controlled by α - 2 and α(D - d) - 2 respectively, and var(v) ≈ 1 D-d var(li). The fact that they have the same mode is consistent with our prior knowledge about the variance of both retained and discarded directions: Since li can be either larger or smaller than v, we don't place prior preference on this choice. The variance relationship between v and li is also consistent with the fact that at the maximum likelihood solution, v is the mean of D - d discarded eigenvalues.

3.2 TheMAP Estimator:

The full Bayesian approach is computationally expensive. In this paper, we use MAP estimation for the parameters, which is very fast and will be shown to perform well with a proper choice of hyper-parameters.

Multiplying the likelihood with the prior probability gives the posterior probability

where n = N +α. It's easy to show that the maximum of the posterior distribution is given by,

Plugging in det(C-1 ) = det(WW T )v-(D-d) and comparing the logarithm of the posterior distribution, equ 7 with equ 9 in [Welling et al., 2003],

we can see that these two equations are equivalent if we replace N in [Welling et al., 2003] by n, and S by (S+αβI)/n. Thus, the MAP estimator for U , L and v is obtained following the same derivation as in [Welling et al., 2003]:

where λi is the eigenvalue of the sample covariance matrix

belonging to the set of C. The ith column of U is the eigenvector of ˆ S/N corresponding to λi, and C and G are sets of respectively retained and discarded directions. As shown in [Welling et al., 2003], all the discarded eigenvalues are contiguous in the ordered eigenspectrum of the sample covariance matrix, and G is determined by comparing all the D - d + 1 candidate sets and choosing the set with the maximal posterior probability or equivalently with the minimal value of the following term:

The complete algorithm is summarized below:

Bayesian XCA MAP Solution

1. Compute the sample mean ˆm and covariance matrix ˆ S/N using eqn 8, 11.
2. Compute the eigenvalues of ˆ S/N and get their estimates, ˆli, using eqn 9
3. Find the optimal set G among D - d + 1 candidates, {j ≤ i ≤ j + d - 1}D-d+1 j=1 with the minimal value of K in eqn 12
4. Compute the estimate ˆ v using eqn 10

In the experiments of this paper, we remove the mean and normalize the variance of each attribute as a preprocessing step. The non-biased estimate of the mean eigenvalues of the covariance matrix is then always 1. Since the hyperparameter β acts as the best prior guess for the eigenvalues, we set it to 1 in this paper. The parameter α is chosen through cross validation. We can also apply this prior to PCA (MCA) to derive a Bayesian PCA (MCA) model by simply imposing a constraint li > v (li < v) on Bayesian XCA. The only difference is that no comparisons between eigenvalue sets are necessary because under the constraint Bayesian PCA (MCA) can only choose PCs (MCs).

The estimators ˆli and ˆ v satisfy some interesting properties. They are both a linear combination of the un-regularized estimator and a constant with a weight proportional to the number of data points N. When the dataset is large enough (N ≫ α), the effect of sampling noise can be neglected and the estimator ˆ li ≈ λi, ˆ v ≈ 1 D-d ∑i∈G λi, equivalent to the ML solution of XCA. When the dataset becomes smaller or α is larger, the constant term gets more weight. The eigenvalue estimator then trades bias for variance. Moreover, in terms of the ordered eigenvalues, as mentioned at the beginning of this section, the un-regularized estimator has a strong bias, which is possibly larger than the bias induced by the prior especially when N < D. Therefore, with properly chosen hyper-parameters, ˆli, ˆ v can reduce both the variance and bias of the ordered eigenvalue estimators.

3.3 Other Eigenvalue Estimators

Besides the MAP estimators, there are also other kinds of estimators for the eigenvalues of the covariance matrix. One is introduced by Lawley [Lawley, 1956] and cited in [Jackson, 2003]. It is able to correct the bias of ordered sample eigenvalues for Gaussian distributions up to O(1/N). We have implemented this method but did not include it in our experimental results because it was quite unstable and even produced negative variance estimates. Another estimator is a shrinkage method with an automatic selection of the shrinkage parameter [Sch¨afer and Strimmer, 2005] 3 . It's designed for the case where the number of data instances is close to or even smaller than the number of attributes. This estimator was included in our experiments described below.

4 Mixtures of Bayesian XCA

Mixtures of factor analyzers [Ghahramani and Beal, 2000] and mixtures of PPCAs [Tipping and Bishop, 1999a; Bishop and Tipping, 1998] are potentially powerful density estimators that combine a number of local dimensionally reduced models into a single joint model. It is not unreasonable to assume that these local models represent constraints, i.e. local patches of data shaped like high dimensional pancakes. Unfortunately, the issue of overfitting is exacerbated for local mixtures because the effective number of data items available to estimate each model component is much smaller. We propose that a mixture of Bayesian XCA models might resolve this.

We therefore introduce priors for each mixture component with shared hyper-parameters. The objective function we want to maximize is the log-likelihood including the regularizing prior given by

where θ = {W, v, m}, zn is the cluster index of xn, and πk = P (zn = k).

The EM algorithm is used to train this model. It has a regular E-step and in the M-step, we execute the BXCA algorithm for each mixture component. Iterating E- and M-steps until convergence is guaranteed to converge to a local maximum of the MAP regularized log-likelihood.

3 The code is available at http://strimmerlab.org/software.html

5 Experiments

5.1 Log-likelihood:

In this section we will compare the test log-likelihood on various datasets for five different models: XCA, PCA, Bayesian XCA, Bayesian PCA, and Shrinkage XCA.

Mosquito Wing Landmarks:

In section 1, we have shown that the first PC and MC of the mosquito wing landmark data represent respectively the large variability of the locations of veins and a constraint on the shape of wings. In this experiment, we want to compare the log-likelihood for different models on this dataset. A number of landmarks on the boundary of the wing are sub-sampled, and the distances of each landmark to its 4 nearest neighbors along the edge are used as features. For the 8 landmarks, there are a total of 13 distances as shown in the middle of two subfigures in Figure 1(b). This is the same number of degrees of freedom after removing translation and rotation information from the original coordinates.

Figure 3(a) shows the average log-likelihood of data points in the training (solid) and test (dashed) sets for Bayesian XCA, XCA and PCA with different numbers of retained directions. We use 50 data points in the training set and 77 in the test set. The inset shows the number of MCs picked by Bayesian XCA and XCA. The log-likelihood of Bayesian XCA and XCA are everywhere above PCA. Bayesian XCA and XCA always pick MCs when d < 11, suggesting that it's better to model the data with constraints. Overfitting is not very serious in this experiment and the plots of Bayesian XCA and XCA are close to each other. Figure 3(b) shows the comparison of Bayesian XCA, Bayesian PCA and Shrinkage XCA w.r.t. the average log-likelihood of test data. Clearly Bayesian XCA performs much better than the other two.

"Frey Faces" Image Data

We have repeated the experiment in [Welling et al., 2003] on the "Frey Faces" images 4 . This dataset contains 1965 images of size 20 × 28. On each pixel, the values are normalized across images to zero mean and unit variance. Figure 4(a),(b) show plots of the average log-likelihood per data point for the various methods. 1000 samples are used in the training set and the remaining 965 in the test set. The corresponding numbers of minor components picked by Bayesian XCA and XCA are plotted in the inset.

XCA and PCA are compared on the same dataset in [Welling et al., 2003] and it turns out that XCA overfits quickly because of the small size of the training set (see Figure 4(a)). With the prior, Bayesian XCA doesn't suffer from this problem. In fact, it remains to perform robustly even when we retain a large number of dimensions. In effect, it will not pick minor components unless sufficiently supported by the data. As a result it outperforms both XCA and PCA. From Figure 4(b) we find that the performance of Bayesian XCA is close to that of Bayesian PCA and Shrinkage XCA underscoring the fact that overfitting is the real issue here.

Figure 5 show the log-likelihood and corresponding numbers of MCs on different sizes of training set. Unlike XCA,

4 Obtained from http://www.cs.toronto.edu/˜roweis/data.html

(b)

Figure 3: (a),(b) show the average log-likelihood per data point of the Mosquito Wing Landmark data on the training set of 50 data points (solid) and test set of 77 data points (dashed) as a function of the number of retained directions. (a) compares Bayesian XCA with XCA and PCA, and (b) compares Bayesian XCA with Bayesian PCA and Shrinkage XCA. Hyper-parameter α = 0.02 is determined through cross validation with d = 5. The inset shows the number of minor components for different numbers of retained directions for Bayesian XCA (o) and XCA (x).

Bayesian XCA always chooses principal components on this dataset, and thus exhibits similar performance as Bayesian PCA. While XCA and PCA over-fit severely as the number of data points decreases the two Bayesian models show good performance across the board on the test set.

In all experiments we determined the value for α using cross validation using a single value of d (nr. of retained dimensions) and N (training size) and subsequently used this hyper-parameter setting for all the other values of d and N.

The experiments suggest that the performance of Bayesian XCA is not very sensitive to the choice of α. To confirm that we also estimated the optimal value of α for a range of d and N values. We found that α was relatively stable and that the test log-likelihood only marginally improved. Results are not presented due to space limitation.

5.2 Mixture Models

We have also tested the performance of Bayesian XCA on a dataset from the UCSD data-mining competition 5 which has

5 Obtained from http://mill.ucsd.edu/index.php?page=Datasets &subpage=Download in the standard classification task.

(b)

Figure 4: (a),(b) show the average log-likelihood per data point of the "Frey Faces" data on 1000 training data (solid) and 965 test data (dashed) as a function of the number of retained components. (a) compares Bayesian XCA with XCA and PCA, and (b) compares Bayesian XCA with Bayesian PCA and Shrinkage XCA. Hyperparameter α = 20 is determined through cross validation with d = 300. The inset shows the number of minor components for Bayesian XCA (o) and XCA (x).

19 attributes. These results are not included in this paper due to the space limitations. However, Bayesian XCA also shows its advantage to other methods on this dataset. Unlike the previous two experiments where either PCs or MCs were chosen by Bayesian XCA, on this dataset it finds a more balanced mix of the two.

We have compared the performance of a mixture of Bayesian XCA models with mixtures of respectively Bayesian PCA, XCA and PCA. We have tested these models on the UCSD dataset. We ran the proposed EM algorithm on 1000 samples with another 1000 samples for validation and a further 38000 samples as the test set. For each model, there are 20 mixture components and the number of retained directions is 10. The EM algorithm was terminated in two ways: it either ran until convergence on training set, or was stopped early by monitoring performance on the validation set. Results averaged over 100 runs are shown in Table 1.

The best log-likelihood is obtained by Bayesian XCA with early stopping, and the worst by XCA running until EM converges. Generally, Bayesian models perform better than nonBayesian models, and early stopping is better than running until convergence. However, the difference between these

α=20, #components=300

# components

200

100

0

10 1

10 2

10 3

# training data Bayesian XCA

Shrinkage XCA

Bayesian PCA

100

200

300

400

500

600

700

800

900

# training data

(b)

Figure 5: (a),(b) show the average log-likelihood per data point of the "Frey Faces" data on the training data (solid) and test data (dashed) as a function of the number of training data points. (a) compares Bayesian XCA with XCA and PCA, and (b) compares Bayesian XCA with Bayesian PCA and Shrinkage XCA. Hyperparameter α = 20. The inset shows the number of minor components for Bayesian XCA (o) and XCA (x).

two stopping criteria for Bayesian models is much less pronounced than for non-Bayesian methods. This suggests that Bayesian methods are able to reduce overfitting to a substantial extent. In a mixture model, Bayesian XCA is still able to combine PCs and MCs reliably and gives a better estimate of the probability density than Bayesian PCA/MCA.

6 Conclusions

We have argued that oftentimes constraints (minor components) better characterize the structure of data than the widely adopted principal components. In particular, scientific discovery seems to be mainly based on discovering conservation laws which directly correspond to minor components.

The primary goal of this paper is to introduce a practical method to reliably extract the optimal combination of principal and minor components from data. An earlier method, called "extreme components analysis" (XCA) [Welling et al., 2003] provided a first step towards this goal but was highly sensitive to sampling noise and as a result suffered from overfitting. The Bayesian extension of XCA we propose in this paper largely resolves this problem. This was verified empirically on various datasets. Moreover, Bayesian XCA pro-

L/N

0

-200

-400

-600

-800

-1000

-1200

-1400

vides a better density estimator than a number of alternative methods such PCA, Bayesian PCA, and XCA based on frequentist corrections to the spectrum. The proposed method is also highly scalable since it is based on a simple singular eigenvalue decomposition of the data matrix.

Code for Bayesian XCA will be released to the public soon. We hope that this will facilitate its widespread use in the scientific community.