Dimensionality Reduction with Generalized Linear Models

1 Introduction

Principal component analysis (PCA) is one of the most widely used methods for dimensionality reduction of multivariate data. Given a set of sample points yi ∈ Rd for i = 1, . . . , n, PCA finds a lower dimensional subspace that minimizes the sum of the square distances from the data points yi to their projections µi in the subspace. In other words, PCA finds a linear reconstruction of the data that is optimal in the least square sense. PCA also can be interpreted from a probabilistic viewpoint that each point yi is thought as a random draw from an isotropic Gaussian with mean µi and variance ψI where ψ is shared among all samples. In other words, each yi is considered as a noise-corrupted sample of the true data point µi which lies in a lower dimensional subspace.

However, the assumption of the Gaussian noise model and the linear reconstruction is not always appropriate for data of different domains. For example, binary data are more suitable to be modeled by Bernoulli with a logit or probit reconstruction (link) function. In order to easily obtain a proper representation for data using appropriate distributions and reconstruction functions, we want to have a general framework for data modeling, so that we can easily incorporate specific distribution and reconstruction function which are suitable for the data. The first step toward this goal was made by [Collins et al., 2002], in which the authors extended PCA to exponential family distributions (called exponential family PCA or EPCA). EPCA is based on the minimization of a generalized criterion in terms of the Bregman divergence, which has a duality relationship with exponential family distributions.

In this paper, we propose a probabilistic model called generalized linear principal component analysis (GLPCA), which is based on the generalized linear model to solve dimensionality reduction problems. Different from the algorithm for the EPCA model, which is based on Bregman divergence minimization, we propose a algorithm which directly maximizes the likelihood using the so called natural gradient ascent on the Riemannian manifold of the parameter space. The algorithm makes incorporating domain specific distributions and reconstruction functions a easy task which simply reduces to input the sufficient statistics to the algorithm. We demonstrate the logit and probit GLPCA algorithms for modeling documents.

It is worth mentioning that the paper [Gordon, 2003] also proposed a very general framework for both regression and factor analysis problems called (GL)2M. This framework in a sense, is more general than the generalized linear model. It uses three link functions instead of one in regular GLM. However, this formulation causes the problem highly non-convex. The subproblem of the optimization procedure is not only non-convex, but it can have exponentially many local minima [Kivinen and Warmuth, 2001]. The author proposed a Newton typed meta algorithm for their model. However, one has to make non-trivial effort to derive the required quantities (e.g. Hessian) in order to incorporate specific distributions and link functions.

In this paper, we denote scalars with lowercase letters (x), vectors with bold lowercase letters (x), matrices with bold uppercase letters (X) and high order tensors with calligraphy uppercase letters (T ).

2 Generalized Linear Principal Component Analysis Model

The classical PCA is a linear model which assumes that an observation yi is generated from a linear transformation of a latent low dimensional vector xi plus a bias termm and a Gaussian noise term ϵ,

The least square estimator is adapted to this assumption. However, the restriction of linearity is too strict for a variety of practical situations. For example, the linear model fails, when dealing with binary data (Bernoulli) or with count data (Poisson).

The generalized linear model (GLM) [Nelder and Wedderburn, 1972] was proposed to relax the assumption to a broad family of nonlinear models for regression problems. In [Rish et al., 2008], the authors utilized the GLM for supervised dimensionality reduction problems. In this paper, we propose the GLPCA model to extend PCA with GLM for dimensionality reduction problems. The essential feature of GLPCA is that the observation yi is a sample from certain distribution p(y). The expectation µi = E[yi] is a monotonic function of ηi such that

Here, we overload the notation f(η) as a vector function that each elements is obtained by applying f over each elements of vector η. The conditional distribution p(y|W, x,m) is a distribution of natural exponential family, which will be discussed later. The reconstruction function f is also called the link function in GLM literature. Given the i.i.d. samples {yi} ∈ Rd , our goal is to find the optimal W ∈ Rd×q , m ∈ Rd , and {xi} ∈ Rq in the sense that the likelihood of this model is maximized.

2.1 Exponential Family

We first introduce the natural exponential family, fromwhich we derive our GLPCA model. Some quantities of the natural exponential family are derived and will be used latter in the algorithm.

In the GLPCA model, we assume that the distribution of the observation yi is a member of the multivariate natural exponential family [Wainwright and Jordan, 2008]. The natural exponential family is a class of probability distributions where the first order statistics are of primary interest. It covers a broad range of distributions, such as Gaussian, Gamma, Bernoulli and Poisson distributions. The density functions of the natural exponential family distributions have a general form

y is a d dimensional observation. θ is the natural parameter. ψ is a scalar dispersion parameter and is a scale parameter of the variance. A(θ) is called the partition function which plays an important role in our latter derivation. The expectation and covariance of y are

For convenience, we denote

We can see that the expectation of y only depends on the natural parameter θ whereas the variance of y depends on both θ and the parameter ψ. Typically one assumes that the factor ψ is identical over all observations. One example is the isotropic Gaussian distribution which can be written in the standard form as

where θ = µ,ψ = σ2 ,A(θ) = θT θ/2, g(µ) = µ, Σ = I, and C(y, ψ) = -(yT y/ψ + ln(2πψ))/2),

If the degree of freedom of y is d, then the representation is said to be minimal; otherwise it is over-complete. If the representation is minimal, the multivariate density can be factorized as

Then it can be easily verified that A(θ) = ∑dj=1 A(θj ) where we overload the notation A(θj) to be the partition function of the univariate distribution p(yj|θj). As a result, the covariance matrix Σ becomes a diagonal matrix with diagonal elements are Σjj = ψA′′ (θj). This is the case for the isotropic Gaussian. For convenience, we define

For a over-complete representation, the degree of freedom of y is smaller than d. For example, if the sample vector represents a discrete distribution, i.e., ∑dj=1 yj = 1 and yj > 0, it is natural to model the data by the Dirichlet distribution. However, if we write the Dirichlet distribution in the form of natural exponential family (3), the representation is not minimal. Another example is that if the samples are normalized such that ‖y‖2 = 1, it is natural to model the data by the von Mises-Fisher distribution, of which the representation in the form of (3) again is not minimal. For certain distributions, transforming from the over-complete representations to the minimal ones is non-trivial. The previous EPCA algorithm is based on minimizing the Bregman divergence, and has an one-to-one map to the natural exponential family with the minimal representation [Wainwright and Jordan, 2008]. Therefore, EPCA only works with the natural exponential family with minimal representation. On the other hand, our GLPCA does not have this restriction. It works with both minimal and over-complete representations.

2.2 Link Function

Apart from the distribution of y, the link function is another important part of GLPCA. From previous section, we have

In the case that θ(η) = η, i.e., f = ∇A, the function f is called canonical link. For models with a canonical link, some theoretical and practical problems are easier to solve. If a canonical link is used and the data are modeled by the natural exponential family with the minimal representation, GLPCA is equivalent to EPCA [Collins et al., 2002] (see Section 4). Table 1 summarizes the characteristics for some exponential family distributions together with natural parameters and canonical link functions. GLPCA also can work with non-canonical link functions. For a specific distribution, a family of special link functions is valid. For example, the (canonical) logit link and cumulative Gaussian link can both be used with the Bernoulli distribution, called logistic and probit model respectively. A very flexible class of link functions is the class of power functions which are also called the Box-Cox transformations [Box and Cox, 1964]. They can be defined for all models for which we have observations with positive means. This family is usually specified as

We summarize some commonly used link functions in Table 2

3 MaximumLikelihood Estimation

In this section, we present the algorithm for maximum likelihood estimation based on the natural gradient ascent in the Riemannian manifold of distributions.

3.1 Natural Gradient in RiemannianManifold of Distributions

Consider a family S of parameterized distributions p(x|θ). Under some smoothness conditions, S = {pθ} can be considered as a differential manifold [Amari and Nagaoka, 2000]. The log-likelihood, denoted as ℓ(y; θ) = ln p(y|θ) is a scalar function in this manifold. The Fisher information matrix is defined as G = [gij ], where

By definition, the G is always symmetric and positive semidefinite. Under certain regularity condition over the distribution p(x|θ) (such as exponential family), the Fisher information matrix is positive definite. For positive definite Fisher information matrix, a Riemannian manifold is induced by defining the Riemannian metric over the differential manifold S as

which is also called Fisher metric or information metric in information geometry literature. It can be proved that the only Riemannian metric is Fisher metric that the geometry is invariant under coordinate transformations of θ and also under one-to-one transformations of random variable y. The technical detail can be found in [Amari and Nagaoka, 2000].

The steepest ascent direction of a function ℓ(θ) in a Riemannian manifold is given by

where ∇ℓ(θ) is the conventional gradient

˜∇ℓ(θ) is called natural gradient of ℓ in the Riemannian manifold [Amari, 1998]. When the manifold is Euclidean space and the coordinate system is orthonormal, we have ˜∇ℓ(θ) = ∇ℓ(θ). This argument suggests the natural gradient ascent algorithm of the form

where τ is the learning rate that determines the step size. Since the Riemannian manifold is invariant under transformation of coordinate θ and random variable y, the natural gradient ascent algorithm is also invariant under the transformation. The natural gradient method has been successfully applied on various machine learning models, including neural network and ICA [Amari and Nagaoka, 2000].

3.2 Maximizing Likelihood of GLPCA

For GLPCAmodel, the log-likelihood of the observations is

where

Our goal is to maximize L w.r.t. X = [xi]ni=1 ,W = [wj ]dj=1 , and m. Note that ψ and C(yi, ψ) have no influence on our objective.

To apply the natural gradient method to GLPCA model, we have to derive the gradients and the Fisher information matrices. However the parameter W and X are matrices. To proceed, in this paper we utilize the tensor notation [Kolda and Bader, 2009] for our derivation.

Taking the first derivative of (16) w.r.t. X, W and m yields the gradients

where

Σi = Σ(µi) and 1 = [1, . . . , 1]T . By definition, one can easily show that the Fisher information matrix is

For our problem, since the parameters are of the matrix form, the Fisher information are tensors, which are given by 1

1 Due to the restriction of the paper length, we omit the proof.

where T is a 4th-order tensor with entries

and all the other entries being zeros. The notation ×t is the tensor t-mode product, i.e., multiplying a tensor by a matrix (or a vector) in mode t.

To compute the natural gradient, we have to unfold the tensors G(·) into the matrices G(vec(·)). The unfolded Fisher information matrix G(vec(X)) is a block diagonal matrix with the diagonal blocks of size d × d given by

The Fisher information matrix G(vec(W)) comprises blocks of size d × d in which block jk is given by

The Fisher information matrix G(m) is simply

The gradients are easily unfolded as

Then the natural gradients w.r.t. to X,Wandm are

Our algorithm for solving the ML solution of (16) is to iterate through the following three steps

until converge. Note that in practice, the natural gradients in (27) need not be exactly computed. One can apply the preconditioned conjugate gradient descent method to solve the linear system

for ˜∇L in a few steps to produce approximated natural gradients with an acceptable precision which is sufficient for updating our parameters. This method is very efficient.

4 Exponential Family PCA

If a canonical link function is used, we recover the exponential family PCA. Following our natural gradient method for ML estimation, we obtain an efficient algorithm for exponential family PCA.

If the canonical link function f(·) = g(·) = ∇A(·) is used in GLPCA, the gradients in (18) become

The tensor T can be simplified with entries T·,i,·,i = Σi and all other entries being zeros. If the distribution is minimally represented in the form of natural exponential family (3), Σi is a diagonal matrix. Denote Vi = T·,i,·,i and Uj = Tj,·,j,·, we have

The algorithm (28) now becomes

In [Collins et al., 2002], the authors proposed the exponential family PCA (EPCA) model which minimizes the Bregman divergence

where BF is the Bregman divergence associated with function F. There is a duality relationship between the Bregman divergence BF (y‖g(θ)) and the natural exponential family p(y) = exp(yθ-A(θ)+C(y)). We can express the negative log-likelihood through Bregman divergence as [Azoury and Warmuth, 2001; Forster and Warmuth, 2002]

The model implicitly assumes that the distribution of the sample vector y satisfies the factorization condition (7). In other words, the distribution of the samples must be minimally represented by (3) and the degree of freedom of any sample yi must be d. This assumption excludes some commonly used distributions, such as von Mises-Fisher. The Bregman divergence based EPCA model is equivalent to a special case of our GLPCA model, where the distribution used is assumed satisfying the factorization condition (7) and the canonical link function corresponding to the distribution is used.

5 Example Models

5.1 Linear PCA

When the Gaussian noise model is assumed and the canonical link (identity function) is used, the GLPCA algorithm simply iterates through following steps

This algorithm is similar to the EM algorithm for PCA (EMPCA) which was introduced in [Roweis, 1998]. The difference is that in EMPCA the bias termm is set to be the mean of the data a priori and subtracted from the data before applying the EM algorithm. In our algorithm, the bias term is optimized together with other parameters. When the algorithm converges, the bias term obtained is not necessarily the mean of the data, since the solution is not unique. However the subspace spanned by W is identical to the one obtained by standard PCA.

5.2 Logistic PCA

Suppose we are given a data set of which each sample yi is a binary valued vector. It is convenient to model the data with the multivariate Bernoulli distribution, i.e.

The canonical link function of Bernoulli distribution is the logistic function. In GLPCA, if Bernoulli and logistic link function are used, we call this model logistic PCA. For logistic PCA, Following (31), we have

where ν(x) = x(1 - x) and µij is the j-th element of µi. Then the ML solution of logistic PCA can be obtained by iterating through (32) with the definition (37).

5.3 Probit PCA

Besides the logistic model, we can also use the probit model for binary data. Probit PCA has the same assumption of multivariate Bernoulli distribution as logistic PCA. The difference is that we use the cumulative Gaussian link µij = Φ(ηij), where

is the CDF of standard normal distribution. Following (22), we have

where φ(·) is the PDF of standard normal distribution and ηij is the j-th element of ηi. The ML solution of probit PCA can also be obtained with algorithm (32). Probit PCA provides a nice alternative to logistic PCA.

6 Experiments

In this section, we demonstrate our GLPCA method on some benchmark data sets. We compare the logistic and probit GLPCA with EPCA on binary data of varying size, dimensionality, and sparsity. Results of standard PCA based on singular value decomposition of mean-centered data is also presented. We use two datasets from the UCI machine learning repository. One is the anonymous Microsoft Web Database and we select the first 5000 instances with 285 binary variables in each instance for experiments. The proportion of nonzero variables in the data is ρ = 0.011. The other dataset is the Advertisement Data consisting of both continuous and binary variables. We only use the 1555 binary features for each instance and remove the corrupted ones. Finally we have 3264 instances with ρ = 0.0082.

We use the reconstruction accuracy as the evaluation criterion. The binary reconstructions are computed by thresholding the continuous values estimated from the lowdimensional representation. Two kinds of measures, a minimum error rate and a balanced error rate, are employed. The minimum error rate is obtained by choosing the best threshold to minimize the overall error rate, and the balanced error rate is obtained by choosing the threshold to equalize the false positive and false negative error rates. As the experimental data are very sparse, both error rates capture different notions of the performance.

The experimental results are shown in Table 3. We can see that both logistic and probit GLPCA outperform the standard PCA on all tasks. Our probit GLPCA outperforms the logistic EPCA algorithm [Schein et al., 2003] in some experiments. It shows different link functions may have strength for modeling different data sets. For binary data, the probit GLPCA is a good alternative of logistic model.

Table 3: Minimum and balanced error rates on two datasets for the task of binary data reconstruction. Each dataset is a d × n binary matrices with ρdn nonzero elements. The best results for each setting are shown in bold.

7 Conclusion

In this paper we proposed a general dimensionality reduction method based on a generalized linear model, called GLPCA. This general method admits a larger range of distributions and link functions. It further extends the previous exponential family PCA model to use non-canonical link functions. An unified meta algorithm based on natural gradient ascent in Riemannian manifold of distributions was proposed for maximizing the likelihood of the model. Deriving the algorithm for specific distribution and link function can be easily done by plugging in corresponding quantities. We also derived the logistic and probit GLPCA from our framework for binary data. Experimental results were shown to validate the effectiveness of the proposed method.