Serge Iovleff
email: serge.iovleff@inria.fr

Probabilistic Auto-Associative Models and Semi-Linear PCA

Keywords: Auto-Associative Models, Non-Linear PCA, Probabilistic Non-Linear PCA

Auto-Associative models cover a large class of methods used in data analysis, among them are for example the famous PCA and the auto-associative neural networks. In this paper, we describe the general properties of these models when the projection component is linear and we propose and test an easy to implement Probabilistic Semi-Linear Auto-Associative model in a Gaussian setting. We show that it is a generalization of the PCA model to the semi-linear case. Numerical experiments on simulated datasets and a real astronomical application highlight the interest of this approach.

Introduction

Principal component analysis (PCA) [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF]; [START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF]; [START_REF] Jolliffe | Principal Component Analysis[END_REF] is a well established tool for dimension reduction in multivariate data analysis. It benefits from a simple geometrical interpretation. Given a set of n points Y = (y 1 , . . . , y n ) with y i ∈ R p and an integer 0 ≤ d ≤ p, PCA builds the d-dimensional affine subspace minimizing the Euclidean distance to the scatter-plot [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF]. The application of principal component analysis postulates implicitly some form of linearity. More precisely, one assumes that the data cloud is directed, and that the data points can be well approximated by their projections to the affine hyperplane corresponding to the first d principal components.

Starting from this point of view, many authors have proposed nonlinear extensions of this technique. Principal curves or principal surfaces methods [START_REF] Hastie | Principal curves[END_REF]; [START_REF] Hastie | The elements of statistical learning[END_REF]; [START_REF] Delicado | Another look at Principal curves and surfaces[END_REF] as well as non-linear transformation of the original data set [START_REF] Durand | Generalized principal component analysis with respect to instrumental variables via univariate spline transformations[END_REF]; [START_REF] Besse | A fixed effect curvilinear model[END_REF] belong to this family of approaches. The auto-associative neural networks can also be viewed as a nonlinear PCA model [START_REF] Baldi | Neural networks and principal component analysis: Learning from examples without local minima[END_REF]; [START_REF] Lu | Quasi-objective nonlinear principal component analysis[END_REF]; [START_REF] Bishop | Pattern recognition and machine learning, information science and statistics[END_REF]; [START_REF] Hinton | Modeling the manifolds of images of handwritten digits[END_REF]. In [START_REF] Girard | Auto-Associative models and generalized principal component analysis[END_REF] we propose the auto-associative models (AAM) as candidates to the generalization of PCA using a projection pursuit regression algorithm [START_REF] Friedman | Projection pursuit regression[END_REF]; [START_REF] Klinke | Projection pursuit regression[END_REF] adapted to the auto-associative case. A common point of these approaches is that they have the intent to estimate an auto-associative model whose definition is given hereafter.

Definition 1 A function g is an auto-associative function of dimension d if it is a map from R p to R p that can be written as g = R • P where P (the "Projection") is a map from R p to R d (generally d < p) and R (the "Restoration" or the "Regression") is a map from R d to R p .

An auto-associative model (AAM) of dimension d is a manifold M g of the form

M g = {y ∈ R p , y -g(y) = 0}
where g is an auto-associative function of dimension d.

For example PCA constructs an auto-associative model using as auto-associative function an orthogonal projector on an affine subspace of dimension d. More precisely we have

g(y) = m + d ∑ i=1 a i , y -m a i , y ∈ R p y, m, a i ∈ R p .
m is a position parameter and the vectors a i are chosen in order to maximize the projected variance. g can be written g = R • P with P(y) = a 1 , ym , . . . , a d , ym and R(x) = m + x 1 a 1 + . . . + x d a d with x = (x 1 , . . . , x d ) . The AAM is then the affine subspace given by the following equation

M g = y ∈ R p ; y -m - d ∑ i=1 a i , y -m a i = 0
Interested reader can check that principal curves, principal surfaces, auto-associative neural networks, kernel PCA Schölkopf et al (1999), ISOMAP Tenenbaum (2000) and local linear embedding [START_REF] Prautzsch | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF] have also the intent to estimate an AAM.

In the PCA approach the projection and the restoration function are both linear. It is thus natural to say that the PCA model is a Linear Auto-Associative Model. In the general case, the manifold M g set can be empty (i.e. the auto-associative function g has no fixed point) or very complicated to describe. Our aim in this paper is to study from a theoretical and practical point of view the properties of some Auto-Associative models in an intermediary situation between the PCA model and the general case: we will assume that the projection function is linear and let the regression function be arbitrary. We call the resulting AAM the Semi-Linear Auto-Associative Models (SLAAM).

Having restricted our study to the SLAAM, we have to give us some criterion to maximize. As we said previously, PCA tries to maximize the projected variance or, equivalently, to minimize the residual variance. Common AAM approaches also used the squared reconstruction error as criterion, or more recently a penalized criterion [START_REF] Hastie | The elements of statistical learning[END_REF]. However as pointed out by M. E. Tipping and C. M. [START_REF] Tipping | Probabilistic principal component analysis[END_REF], one limiting disadvantage of this approach is the absence of a probability density model and associated likelihood measure. The presence of a probabilistic model is desirable as the definition of a likelihood measure permits comparison between concurrent models and facilitates statistical testing, -A single AAM may be extended to a mixture of such models, if a probabilistic AAM is used to model the class conditional densities in a classification problem, the posterior probabilities of class membership may be computed.

We propose thus a Gaussian generative model for the SLAAM and try to estimate it using a maximum likelihood approach. In the general case we are faced with a difficult optimization problem and we cannot go further without additional assumptions.

It will appear clearly that if the projection P, identified hereafter with a matrix P, is known then the estimation problem of a SLAAM is very close to an estimation problem in a regression context. However, there are some differences we will underline. In particular it will appear that in order to get tractable maximum-likelihood estimates, we have to impose some restrictions to the noise. We call the resulting model of all these assumptions/simplifications a Semi-Linear Principal Component Analysis. It does not seem possible to add non-linearity to the PCA model and get tractable likelihood estimate for P. But clearly, the assumption that P is known is too strong in practice. We thus propose to estimate it in a separate step using either PCA or a contiguity analysis [START_REF] Lebart | Contiguity analysis and classification[END_REF]. Finally, even if P is assumed to be known, it remains to estimate the regression function R which is a non-linear function from R d to R p . If d > 1 and p is moderately high the task becomes very complicated. Thus the implementation of the model proposed in this paper simplifies it once more. In this implementation, we assume that R is additive inspired by the Generalized Additive Model (GAM) approach [START_REF] Hastie | Generalized Additive Models[END_REF].

In view of the experiments we have performed and we present here, it seems we obtain a practical and simple model which generalizes the PCA model to the nonlinear case in an understandable way. The main advantages of the proposed model are:

1. the projection step being linear, it is easy to interpret the projected values X = (x 1 , . . . , x n ) in term of the original variables Y and to represent them in a (kind of) bi-plot graph (see [START_REF] Gower | Biplots[END_REF]), 2. the regression functions can be drawn individually, i.e. we can fix p -1 coordinates to a value, say for example x k = 0 if k = j, and, using the additivity, represent the maps x j → R j (x j ) for j = 1, . . . , d, in order to get a better understanding of the (eventually) non-linearity present in the data in term of x j , 3. the probabilistic model allows to choose the intrinsic dimension of the projection and the regression function using a model selection criterion.

Some illustrative tables and graphs can be found in the examples given in this paper (section 5) and illustrate these facts. The main disadvantages of the method are 1. its use of intensive numerical computations, so that it cannot be applied to very big data sets (either in the number of individuals or number of variables), 2. the use of linear projection functions doesn't allow to modelize data presenting too strong non-linearity (see proposition 3).

The paper is organized as follows. Section 2 introduces the Probabilistic Semi-Linear Auto-Associative Models and relate them to the PCA and Probabilistic PCA models. In section 3 we present the Semi-Linear PCA models and the estimation of theirs parameters conditionally to the knowledge of the projection matrix P. Section 4 focus on implementation details, in particular section 4.1 is devoted to the determination of the projection matrix P using either PCA or contiguity analysis, and section 4.2 presents the additive B-Spline regression. Data sets and experiments are detailed in Section 5 with a real astronomical data set. Some comparisons are done with the usual PCA method. Finally, some concluding remarks are proposed in Section 6.

2 Semi-Linear Auto-Associatif Models (SLAAM)

Geometrical properties of the SLAAM

Let us first consider a general auto-associative model M g with g = R • P as given in the definition 1. We have the following evident property

Proposition 1 Let H = {P(y); y ∈ M g } ⊂ R d .

On H the projection function and the regression function verify

P • R = Id d (1)
where Id d denote the identity function of R d .

Proof Let y ∈ M g and let x = P(y), then

x = P(y) = P(g(y)) = P(R(P(y))) = P(R(x)).
As a consequence, we have the following "orthogonality" property verified by an AAM when P is an additive function Proposition 2 Let V = {P(y); y ∈ R p } and assume that the property (1) extends on V , let y ∈ R p , ȳ = R(P(y)) and ε = yȳ. If P is additive, i.e. P(y+y ) = P(y)+P(y ), then P(ε) = 0.

Proof Using the property (1), we have on one hand P(ȳ) = P(R(P(y))) = P(y).

While on the other hand P(ȳ) = P(y -ε) = P(y) -P(ε) giving the announced result.

Clearly we have H ⊂ V and the assumption given in this proposition seems quite natural. We focus now on the semi-linear case and we assume that

P(y) = a 1 , y , . . . , a d , y = Py. (2)
with P = (a 1 , . . . , a d ) a matrix of size (d, p).

Proposition 3 Let g = R • P be an auto-associative function, with P given in (2) and R verifying the property (1). Let B = a 1 , . . . , a d , a d+1 , . . . , a p be an orthonormal basis of R p with (a d+1 , . . . , a p ) chosen arbitrarily. Let y ∈ M g , and let ỹ and R denote respectively the vector y and the auto-associative function R in the basis B, then

          ỹ1 . . . ỹd ỹd+1 . . . ỹp           =           ỹ1 . . . ỹd Rd+1 ( ỹ1 , . . . , ỹd ) . . . Rp ( ỹ1 , . . . , ỹd )           . ( 3
)
Proof It is sufficient to notice that the change of basis matrix Q is given by

Q = a 1 , . . . , a d , a (d+1) , . . . , a p ,
thus the left multiplication of y and R by Q, using (1), will give (3).

From this last proposition we can see that the Semi-Linear Auto-Associative models have a relatively simple geometrical structure and that we cannot expect to model highly non-linear models with them.

Probabilistic Semi-Linear Auto-Associative Models

In the sequel, we will denote by V the subspace spanned by the set of vectors (a 1 , . . . , a d ), and give us an arbitrary orthonormal basis of V ⊥ denoted by (a d+1 , . . . , a p ). We will denote by P the matrix (a 1 , . . . , a d ) and by P the matrix (a d+1 , . . . , a p ) . As in proposition 3, Q represents the unitary matrix (P| P) = (a 1 , . . . , a p ) .

General Gaussian Setting

Definition 2 Let x be a d-dimensional Gaussian random vector:

x ∼ N (µ x , Σ x ) (4)
and let ε be a p-dimensional centered Gaussian random vector, independent of x, with a diagonal covariance matrix Σ ε = Diag(σ 1 , . . . , σ p ). The p-dimensional vector y is a Probabilistic Semi-Linear Auto-Associative Model (semi-linear PCA) if it can be written as

y = Q                     x 1 . . . x d rd+1 (x) . . . rp (x)           + ε          = R(x) + ε, (5)
where the r j (x), d + 1 ≤ j ≤ p, are arbitrary real functions from R d to R.

Link with the Principal Component Analysis

Assume that:

1. r j (x) = μ j for all j ∈ {d + 1, . . . p}, 2. the covariance matrix of x,

Σ x = Diag(σ 2 1 , . . . , σ 2 d ) is diagonal with σ 1 ≥ σ 2 ≥ . . . ≥ σ d , 3
. the Gaussian noise ε has the following covariance matrix Σ ε = Diag(0, . . . , 0, σ 2 , . . . , σ 2 ) (d zeros and pd sigmas on the diagonal) with σ < σ d then the vector y is a Gaussian random vector

y ∼ N (µ, Σ ) with µ = Q           μ1 . . . μd μd+1 . . . μp           and Σ = Q           σ 1 . . . 0 σ d σ 0 . . . σ           Q
and a 1 , . . . , a d are the d first eigenvectors given by PCA.

Link with the Probabilistic Principal Component Analysis

The probabilistic PCA [START_REF] Tipping | Probabilistic principal component analysis[END_REF] is a model of the form

y = µ + Wx + ε, (6)
with W a (p, d) matrix, x a d-dimensional isotropic Gaussian vector, i.e. x ∼ N (0, I d ), and ε a p-dimensional centered Gaussian random vector with covariance matrix σ 2 I p . The law of y is not modified if W is right multiplied by a (d, d) unitary matrix, it is thus possible to impose to the rows of W to be orthogonal (assuming that W is of full rank).

The following proposition is then straightforward Proposition 4 Assume that ε (and thus ε) is an isotropic Gaussian noise, i.e. Σ ε = σ 2 I p , take r j = μ j for all d + 1 ≤ j ≤ p and set

W = P       σ 1 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 σ d      
.

The resulting Probabilistic Semi-Linear Auto-Associative Model is a Probabilistic Principal Component Analysis.

For this simple model there exists a close form of the posterior probability of y and for the maximum likelihood of the parameters of the model. In particular, the matrix W can be estimated up to a rotation and spans the principal subset of the data.

3 Semi-Linear PCA Our aim is now to generalize the PCA model we present in part 2.2.2 to the semilinear case and we assume in the rest of this article that [N] the Gaussian noise ε have the following covariance matrix

Σ ε = Diag(0, . . . , 0, σ 2 , . . . , σ 2 ).
Σ ε is a diagonal matrix with d zeros and pd sigmas on the diagonal.

As we said in the introduction, adding non-linearity make everything more difficult and a practicable/tractable implementation for a real usage forces us to make some simplifications. The model is estimated using a two stages strategy: we first determine the projection matrix P using either the first components of a PCA or using a contiguity analysis (section 4.1), and in a second step we estimate the parameters and the regression function R (section 4.2).

Using [N] and expressing y in the basis B (definition 2) by rotation we get the following expression for ỹ:

          ỹ1 . . . ỹd ỹd+1 . . . ỹp           =           x 1 . . . x d rd+1 (x) . . . rp (x)           +           0 . . . 0 εd+1 . . . εp           . (7)
In other words, the coordinates of ỹ can be split in two sets. The first d coordinates are the Gaussian random vector x, while the remaining pd coordinates are a random vector z which is conditionally to x a Gaussian random vector N r(x), σ 2 I p-d .

Observe that the regression functions are dependents of the choice of the vectors a d+1 , . . . , a p and that, as the noise ε lives in the orthogonal of V , we have x = Py.

The parameters we have to estimate are the position and correlation parameters µ x and Σ x for the x part and the noise and the regression functions (σ 2 , r) for the non-linear part. Given a set of n points Y = (y 1 , . . . , y n ) in R p , we get by projection two sets of n points X = (x 1 , . . . ,

x n ) = YP in R d , and Z = (z 1 , . . . , z n ) = Y P in R p-d .
Standard calculation give the Gaussian maximum likelihood for µ x and

Σ x μx = 1 n n ∑ i=1 x i . ( 8
) and Σx = 1 n n ∑ i=1 (x i -μx )(x i -μx ) . (9)
The maximum likelihood of σ 2 is given by

σ 2 = 1 n(p -d) n ∑ i=1 z i -r(x i ) 2 . ( 10
)
It remains to estimate r. This will be obtained by minimizing σ 2 . Many parametric or non-parametric regression tools can be used for this purpose. A practical method based on B-Spline is presented in the next section.

Computational aspects

This section is devoted to the description of two main computational steps: the determination of a "good" projection matrix (section 4.1), and the estimation of the regression function (section 4.2). We also present briefly the model selection criteria AIC and BIC (section 4.3) and we terminate with some practical implementation details (4.4). All the methods described are implemented in a C++ library called aam.

Computing the Projection matrix : PCA and Contiguity Analysis

Given (a 1 , . . . , a d ) an orthonormal set of vector in R p , an index I: R p×d → R + is a functional measuring the interest of the projection of the vector y on Vec(a 1 , . . . , a d ) with a non negative real number. The choice of the index I is crucial in order to find "good" parametrization directions for the manifold to be estimated. We refer to [START_REF] Huber | Projection Pursuit[END_REF] and [START_REF] Jones | What is projection pursuit?[END_REF] for a review on this topic in a projection pursuit setting. The meaning of the word "good" depends on the considered data analysis problem. For instance, Friedman et al [START_REF] Friedman | A Projection Pursuit algorithm for exploratory data analysis[END_REF]; [START_REF] Friedman | Exploratory projection pursuit[END_REF], and more recently Hall [START_REF] Hall | On polynomial-based projection indices for exploratory projection pursuit[END_REF], have proposed an index which measure the deviation from the normality in order to reveal more complex structures of the scatter plot. An alternative approach can be found in [START_REF] Caussinus | Metrics for finding typical structures by means of Principal Component Analysis[END_REF] where a particular metric is introduced in PCA in order to detect clusters. We can also mention the index dedicated to outliers detection [START_REF] Pan | Multiple outlier detection in multivariate data using projection pursuit techniques[END_REF].

A widely used choice of I is I( a 1 , y , . . . , a d , y ) = tr(Var [Py]), the projected variance. This is the criterion maximized in the usual PCA method [START_REF] Jolliffe | Principal Component Analysis[END_REF]. It can be computed using the total variance matrix

V = 1 n n ∑ i=1 (y i -μ)(y i -μ) . ( 11
)
and by maximizing

I(a 1 , . . . , a d ) = d ∑ i=1 a i Va i , a i , a j = δ i j , (12)
where δ i j represent the Kronecker's delta. The resulting vectors are then the eigenvectors associated to the d largest eigenvalues of V.

An other approach generalizes the one presented in [START_REF] Girard | Auto-Associative models and generalized principal component analysis[END_REF] and consists in defining a contiguity coefficient similar to Labart one's [START_REF] Lebart | Contiguity analysis and classification[END_REF] whose maximization allows to unfold nonlinear structures. A contiguity matrix is a n × n boolean matrix M whose entry is m i j = 1 if data points i and j are "neighbors" and m i j = 0 otherwise. Lebart proposes to use a threshold r 0 to the set of n(n -1) distances in order to construct this matrix (but the choice of r 0 could be delicate if the data scale is not homogeneous) or to use a k-contiguity matrix, i.e. to set m i j = 1 iff j is one of the k-nearest neighbor of i. This is the approach we have implemented.

The contiguity matrix being determinate, we compute the local covariance matrix

V * = 1 2kn n ∑ i=1 n ∑ j=1 m i j (y i -y j )(y i -y j ) . (13)
The axis of projection are then estimated by maximizing the contiguity index

I(a 1 , . . . , a d ) = d ∑ i=1 a i V * a i a i Va i . ( 14
)
which can be expected to unfold the scatter-plot since distant observations are ignored in computing the local covariance matrix. Using standard optimization techniques, it can be shown that the resulting axis are the d eigenvectors associated with the largest eigenvalues of the matrix V * -1 V . The behavior of the PCA and contiguity indexes is illustrated on the next two figures (figure 1). Three hundred data have been simulated uniformly on the intervals [-1, 1] 2 and [-2, 2] 2 and the third coordinate have been computed on the parabola z = x 2 +y 2 . In the first case PCA and contiguity methods give almost the same results, while in the second case the contiguity method preserve the local structure of the scatter-plot. Remark: In all the examples k has been fixed to 3. It seems to be a good compromise between the computational time (the computation of V * increase with k) and the stabilization of the index.

Estimating the Regression function: Linear regression and additive B-Spline regression

Linear Regression

In the linear case, we are looking for a vector µ and a (d,

p -d) matrix R minimizing n ∑ i=1 z i -µ -R x i 2 .
It is easily verified that

μ = 1 n n ∑ i=1 z i -R x i = μz -R μx .
Setting X = X -1 μ x and Z = X -1 μ z , where 1 represent a vector of size n with 1 on every coordinates. Assuming that the matrix X X is invertible, standard calculus show that R = ( X X) -1 X Z.

Finally, using the decomposition in eigenvalues of the covariance matrix of Y, it is straightforward to verify the following theorem Theorem 1 If the d orthonormal vectors a 1 , . . . , a d are the eigenvectors associated with the first d eigenvalues of the covariance matrix of Y then the estimated autoassociative model is the on obtain by PCA.

Additive B-Spline regression

B-Spline regression allows to approximate a curve

s(x) defined from the interval [a, b] into R p using a set of n observations (x i , y i ) n i=1 ranging over [a, b] × R p . Given a set of m + k + 1 real values t i , called knots, such that t 0 ≤ t 1 ≤ . . . ≤ a = t k-1 ≤ . . . ≤ t m+1 = b ≤ . . . ≤ t m+k .
A B-Spline curve of degree k is defined by

s(x) = m ∑ l=0 α l s l (x), x ∈ [t k-1 ,t m+1 ].
Points (α l ) m l=0 are called the control points and the functions (s l ) m l=0 are the m + 1 B-Spline basis of order k computed using the Cox-De Boor recursion formula (see de Boor (1978) for details). m and k are numbers fixed by the user.

Given d sets of m+k +1 knots t j i , an additive B-spline curve of degree k is defined by

s(x) = d ∑ j=1 m ∑ l=0 α l s jl (x j ), x j ∈ [t j k-1 ,t j m+1 ].
For simplicity, the number of control point have been set equal for all j.

In order to estimate the regression function r, we express it as an additive linear combination of m + 1 B-Spline basis functions r jl . We have thus to estimate the set of coefficients (α jl ), j = 1, . . . , d, l = 0, . . . m by minimizing

n ∑ i=1 z i -α 0 - d ∑ j=1 m ∑ l=0 α jl r jl (x i j ) 2 .
Standard regression techniques give then the estimates α = ((S S) -1 S Z) where S is the design matrix which depends on the knots position, the degree of the B-Spline and the number of control points chosen by the user [START_REF] Prautzsch | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF]. The regression function r is then estimated by the formula

r(x) = α0 + d ∑ j=1 m ∑ l=0 α jl r jl (x j ) = α0 + d ∑ j=1 r j (x j ).

Remark:

In the examples presented below we use cubic B-Spline basis but the aam library allow the user to chose its own degree. The knots t j are regularly spaced in the range of (x i j ) n i=1 .

Model Selection

Since a Semi-linear PCA model depends highly of the projection matrix P and the regression function, model selection allows to select among various candidates the best projection and the number of control points. Several criteria for model selection have been proposed in the literature and the widely used are penalized likelihood criteria. Classical tools for model selection includes Akaike [START_REF] Akaike | A new look at the statistical mode identification[END_REF] and Bayesian [START_REF] Schwarz | Estimating the dimension of a model[END_REF] information criteria. The Akaike Information Criterion (AIC) is a measure of the relative goodness of fit of a statistical model given by -2 ln(L) + 2γ(M ) and the Bayesian Information Criterion (BIC) consists in selecting the model which penalizes the log-likelihood by γ(M ) 2 log(n) where γ(M ) is the number of parameters of the model M and n is the number of observations. Remark: In this paper we use BIC because of its popularity but the aam library allows to use AIC.

Estimation in practice

The drawback of the previous maximum likelihood equations ( 7) and ( 10) is that we have to perform a rotation of the original data set in order to estimate the regression function and next to perform and inverse rotation of the estimated model. In practice, having fixed the dimension d of the model and the number of control point m, we avoid such computations by estimating the model using the following steps:

-Center and (optionally) standardize the data set Y: obtain Ȳ, -Compute the projected data set X = ȲP (X is centered), -Compute the regression Ȳ ∼ X (without intercept), -Compute the residual variance σ 2 , the log-likelihood and the BIC.

As we can see the main difference is in the regression step: we estimate directly a function from R d to R p . In practice, as the non-linear part of the model is in V ⊥ , the regression function we obtain numerically give the identity function in the V space.

These estimations are repeated for d (the dimension of the auto-associative model) varying from 1 to a fixed value d max , and for m (the number of control point) varying from 4 to a fixed value m max (for a cubic B-Spline interpolation we cannot have less than 4 control points). The model with smallest BIC is selected.

Examples

We first present two illustrations of the estimation principle of semi-linear PCA on low dimensional data (Section 5.1 and 5.2). Second, semi-linear PCA is applied to an astronomical analysis problem in section 5.3. In all the examples, we use an additive B-Spline regression model for the estimation of the regression function R (section 4.2.2). The B-Spline are of degree 3 and we select the number of control points using BIC. All the results are compared with PCA results.

First example on simulated data

The data are simulated using a one-dimensional regression function in R 3 . The equation of the AA model is given by u → (u, sin u, cos u),

and thus P(u, v, w) = u. The first coordinate of the random vector is sampled 1000 times from a centered Gaussian distribution with standard deviation σ x = 3. An independent noise with standard deviation σ = 1 is added to the v and w coordinates. The axis of projection have been computed by PCA and by contiguity analysis (section 4.1) using the 3 nearest neighbors for the proximity graph. The correlations between the projected data set x and the original data set are y 1 y 2 y 3 x (contiguity) 0.999968 -0.000579458 0.0089 x (PCA) 0.999997 0.008813556 0.0118 which show that the first axis given either by contiguity analysis or PCA is very close from the x-axis as it was expected. The projected variance on the first axis find by contiguity analysis is 9.20496 which is also very close to 9. We use BIC in order to select the dimension of the model and the number of control points using an additive B-Spline regression. A summary of the tested model is given in the table 1 and the estimated AAM using B-Spline regression or PCA is drawn in the figures 2 and 3.

Second example on simulated data

In our second example the AAM is given by

(u, v) → u, v, u 2 + v 2 (16)
and thus P(u, v, w) = (u, v). The first two coordinates of the random vector are sampled from a centered Gaussian distribution with covariance matrix Σ x = 0.5 2 0 0 0.55 2 and n = 1000 points are simulated. An independent noise with standard deviation σ = 0.2 has then been added to the z coordinate.

The result of the contiguity analysis with k = 3 neighbors is displayed in the figure 4 (but only the 2-neighbors graph is drawn). The correlation circle shows that the first component obtains by the contiguity analysis is essentially correlated with the second variable (y-axis) and that the second component is very correlated with the first variable (x-axis) and a bit with the third variable (z-axis). The last graphic shows that in contrast the first component obtain by PCA is essentially correlated with the z-axis.

The model is estimated using an additive B-Spline regression. BIC selects 6 control points (40 parameters) and the residual standard deviation is 0.309132 which overestimate a bit the true level of noise. The true model and the estimated model obtained with an additive B-Spline regression are given in the figure 5 and compared with the usual PCA. The figure 6 draw the estimation of the regression functions (R j k ) for j = 1, . . . , 2 and k = 1, 2, 3.

Example in spectrometry analysis

Finally we illustrate the performance of the semi-linear PCA on a real data set. The data consists of 19-dimensional spectral information of 487 stars [START_REF] Stock | Quantitative stellar spectral classification[END_REF]; [START_REF] Garcia | Quantitative Stellar Spectral Classification. III. Spectral Resolution[END_REF]; [START_REF] Stock | Quantitative Stellar Spectral Classification. II. Early Type Stars[END_REF]; [START_REF] Garcia | Quantitative Stellar Spectral Classification. IV. Application to the Open Cluster IC 2391[END_REF] and they have 1. The axis of projection given by the PCA outperform largely the results we obtain with the contiguity analysis, for any choice of the number of control points. The scatter plot in the main PCA space with the correlation circle and the individuals regression functions in the plane (6,7) of the PCA. are given in the figure 7. A summary of some tested model are given in the table 2. We present a visualization of some of the estimated regression functions in figure 8. There is 5 × 19 = 95 regression functions to plot but we show only the first five variables (y 1 , . . . , y 5 ) in term of the projected variables x = (x 1 , . . . , x 5 ). The first projected variable x 1 fits very well the variables y j . Table 2 Values of BIC for various number of control points (given in the first column). BIC selects the model of dimension 5 with 9 control points using as projection matrix the 5 axis given by PCA. The total variance (inertia) of the data set was 26.59832.

Fig. 8 The estimated regression functions from R 5 to R 19 . Only the first five dimensions are given.

Conclusion

We have presented a class of auto-associative models for data modeling and visualization called semi-linear auto-associative models. We provided theoretical groundings for these models by proving that the principal component analysis and the probabilistic principal component analysis are special cases. This model allows to model data set with a simple non-linear component and is truly generative with an underlying probabilistic interpretation. However it does not allow to model data with a strong non-linear component and it depends highly on the choice of the projection matrix. The Semi-Linear PCA have been implemented in C++ using the stk++ library [START_REF] Iovleff | The Statitiscal ToolKit[END_REF] and is available at: https://sourcesup.renater.fr/projects/ aam/. The program is accompanied with a set of R scripts which allows to simulate data sets and display the results of the aam program.

Fig. 1

 1 Fig.1Planes of projection using the PCA index (dark gray) and the contiguity index (light gray) with 3 neighbors. The planes have been shifted from the origin for display purpose.

Fig. 2

 2 Fig.2(a) The simulated scatter-plot, the estimated AAM (dark grey) and the true AAM (thick grey) using a contiguity analysis and a B-Spline regression. (b) Results using an usual PCA (the straight line). This graphic is obtained with R using the draw3d command of the aam library.

Fig. 3

 3 Fig. 3 Estimated regression functions x → R k (x) and the scatter plots (x i , y ik ) n i=1 for k = 1, 2, 3. The true model is given in the equation 15. This graphic is obtained with R using the draw.regression command of the aam library.

Fig. 4

 4 Fig.4The AAM components and the 2-neighbors graph, the correlation circle of the AAM components with the variables and representation of the scatter-plot in the PCA plan. These pictures have been obtained using the plot command of the aam library.

Fig. 6

 6 Fig.6Estimated regression functions x j → R j k (x j ) and the scatter plots (x i j , y ik ) n i=1 for k = 1, 2, 3 and j = 1, 2. The true model is given in the equation 16. This graphic is obtained with R using the draw.regression command of the aam library.

Fig. 7

 7 Fig. 7 Scatter plot in the PCA space, correlation circle of the first two components with the variables and plane (6,7) of the PCA with the estimated regression functions.

Table 1

 1 Values of BIC

			Contiguity Analysis		PCA
		dim	BIC	Residual Variance	BIC	Residual Variance
	linear	1	9987.13(5)	1.439	11495.8	1.43864
		2	10609.7(10)	1.40815
	9	1	11247.6(29)	1.06557	11058.4	1.06407
		2	11621(58)	1.1086
	10	1	11060.7(32)	0.986316	10926.3	0.985777
		2	11469(64)	0.913605
	11	1	10853.7(35)	0.941064	10855	0.92711
		2	11503(70)	0.90682
	12	1	10844.6(38)	0.92717	10845(38)	0.941661
		2	11525(76)	0.892669
	13	1	10871.4(41)	0.929965	10871.6	0.927976
		2	11568.4(82)	0.891119
	14	1	108887.5(44)	0.927834	10888.3	0.927976
		2	11604(88)	0.885939

Acknowledgements The author thanks two anonymous referees for comments that improved the readability of this article.


Parmeet Bhatia

Serge Iovleff

Gérard Govaert

blockcluster: An R Package for Model Based Co-Clustering

Keywords: model based clustering, block mixture model, EM and CEM algorithms, simultaneous clustering, co-clustering, R, blockcluster

published or not. The documents may come

Introduction

Cluster analysis is an important tool in a variety of scientific areas such as pattern recognition, information retrieval, micro-array, data mining, and so forth. Although many clustering procedures such as hierarchical clustering, k-means or self-organizing maps, aim to construct an optimal partition of objects or, sometimes, of variables, there are other methods, called block clustering methods, which consider simultaneously the two sets and organize the data into homogeneous blocks. Let x denotes a n × d data matrix defined by x = {(x ij ); i ∈ I and j ∈ J}, where I is a set of n objects (rows, observations, cases etc) and J is a set of d variables (columns, attributes etc). The basic idea of these methods consists in making permutations of objects and variables in order to draw a correspondence structure on I × J. For illustration, consider Figure 1 where a binary data set defined on set of n = 10 individuals I = A, B, C, D, E, F, G, H, I, J and set of d = 7 binary variables J = 1, 2, 3, 4, 5, 6, 7 is reorganized into a set of 3 × 3 clusters by permuting the rows and columns. In recent years, co-clustering have found numerous applications in the fields ranging from data mining, information retrieval, biology, computer vision and so forth. [START_REF] Dhillon | Co-Clustering Documents and Words Using Bipartite Spectral Graph Partitioning[END_REF]] published an article on text data mining by simultaneously clustering the documents and content (words) using bipartite spectral graph partitioning. This is quite useful technique for instance to manage huge corpus of unlabeled documents. [START_REF] Xu | Co-Clustering Analysis of Weblogs Using Bipartite Spectral Projection Approach[END_REF]] presents another co-clustering application (again using bipartite spectral graph) to understand subset aggregates of web users by simultaneously clustering the users (sessions) and the page view information. [START_REF] Giannakidou | Co-Clustering Tags and Social Data Sources[END_REF]] employed a similarity metric based co-clustering technique for social tagging system. In field of Bio-informatics, co-clustering is mainly used to find structures in gene expression data. This is useful for instance to find set of genes correspond to particular kind of disease. Some of the pioneer material in this context can be found in [START_REF] Kluger | Spectral Biclustering of Microarray Data: Coclustering Genes and Conditions[END_REF]]. Recently many co-clustering based algorithms have also been developed to target computer vision applications. For instance, [START_REF] Qiu | Image and Feature Co-Clustering[END_REF]] demonstrated the utility of co-clustering in Image grouping by simultaneously clustering images with their low-level visual features. [START_REF] Guan | Spectral Images and Features Co-Clustering with Application to Content-Based Image Retrieval[END_REF]] extended this work and presented opportunity to develop a novel content based image retrieval system. Similarly, [START_REF] Rasiwasia | Holistic Context Modeling Using Semantic Co-Occurrences[END_REF]] used co-clustering to model scenes. We have presented here with a very brief survey demonstrating the diversity of co-clustering applications. Unfortunately, many of the algorithms developed above for co-clustering are suitable in that particular context only and cannot be considered as generic framework. Here, we take the opportunity to emphasize on the fact that the algorithms used in our software can be more or less adopted in most or all the above context.

In practice, several softwares exist for block clustering such as Bicat [START_REF] Barkow | BicAT: A Biclustering Analysis Toolbox[END_REF]], biclust [START_REF] Kaiser | A Toolbox for Bicluster Analysis in R[END_REF]], BiGGEsTS [START_REF] Goncalves | BiGGEsTS: Integrated Environment for Biclustering Analysis of Time Series Gene Expression Data[END_REF]], or refer to [START_REF] Madeira | Biclustering Algorithms for Biological Data Analysis: A Survey[END_REF]] for a review. In this paper, we introduce our R [R Development Core Team (2008)] package blockcluster, which provides a bridge between a newly developed C++ library and the R statistical computing environment. The R package blockcluster allows to estimate the parameters of the co-clustering models [START_REF] Govaert | Clustering with Block Mixture Models[END_REF]] for binary, contingency and continuous data. This package is unique from the point of view of generative models it implements (latent blocks), the used algorithms (BEM, BCEM) and, apart from that, special attention has been given to design the library for handling very huge data sets in reasonable time. The R package is already available on CRAN at http://cran.r-project.org/web/packages/blockcluster/index.html.

The rest of the article is organized as follows. In Section 2, we give the mathematical foundation (without going into detailed proofs) of latent block models for binary, contingency and continuous data-sets. In Section 3, we elaborate various details of blockcluster package and illustrate its usage with the help of various examples. This part can also act as first-hand comprehensive tutorial for the audience who are interested in running the blockcluster package directly. Finally we illustrate our package on real data-sets in Section 4 and terminate with some concluding remarks in Section 5.

Block mixture models

Mixture models: Classic formulation

In model-based clustering, it is assumed that the data arises from a mixture of underlying probability distributions, where each component k of the mixture represents a cluster. Hence the matrix data x=(x 1 , . . . , x n ) is assumed to be i.i.d and generated from a probability distribution with the density

f (x i ; θ) = k π k f k (x i ; α),
where f k is the density function for kth component. Generally, these densities belong to the same parametric family (say for example Gaussian). The parameters (π k ) g k=1 give the probabilities that an observation belongs to the kth component and the number of components in the mixture is g which is assumed to be known. The parameter θ of this model is the vector (π 1 , . . . , π g , α) containing the parameter α and the mixing proportions (π 1 , ..., π g ). Using this, the density of the observed data x can be expressed as

f (x; θ) = i k π k f k (x i ; α). (1)
An important property of the pdf f (x; θ) shown by [START_REF] Govaert | Clustering with Block Mixture Models[END_REF] is that the density (1) can be rewritten as

f (x; θ) = z∈Z p(z)f (x|z; α) (2)
where Z denotes the set of all possible partitions of I in g clusters, f (x|z; α) = i f z i (x i ; α) and p(z) = i π z i . With this formulation, the data matrix is assumed to be a sample of size 1 from a random (n, d) matrix.

Latent block model: General formulation

Let the data-set x be defined on a set I with n elements (individuals) and a set J with d elements (variables). We represent a partition of I into g clusters by z = (z 11 , . . . , z ng ) with z ik = 1 if i belongs to cluster k and z ik = 0 otherwise, z i = k if z ik = 1 and we denote by z .k = i z ik the cardinality of row cluster k. Similarly, we represent a partition of J into m clusters by w = (w 11 , . . . , w dm ) with w j = 1 if j belongs to cluster and w j = 0 otherwise, w j = if w j = 1 and we denote w . = j w j the cardinality of column cluster . Furthermore, note that for any variable V ij , we denote V .k to represent i V ik and similarly

V i. shall represent k V ik .
To study the block clustering problem, the previous formulation ( 2) is extended to propose block mixture model defined by the following probability density function

f (x; θ) = u∈U p(u; θ)f (x|u; θ)
where U denotes the set of all possible labellings of I × J and θ contains all the unknown parameters of this model. By restricting this model to a set of labellings of I × J defined by a product of labellings of I and J, and further assuming that the labellings of I and J are independent of each other, we obtain the decomposition

f (x; θ) = (z,w)∈Z×W p(z; θ)p(w; θ)f (x|z, w; θ) (3)
where Z and W denote the sets of all possible labellings z of I and w of J.

By extending the latent class principle of local independence to our block model, each data point x ij will be independent once z i and w j are fixed. Hence we have

f (x|z, w; θ) = i,j f z i w j (x ij ; α)
where f z i w j (x, α) is either a probability density function (pdf) defined on the real set R, or a discrete probability function. Denoting θ = (π, ρ, α), where π = (π 1 , . . . , π g ) and ρ = (ρ 1 , . . . , ρ m ) are the vectors of probabilities π k and ρ that a row and a column belong to the kth row component and to the th column component respectively, we obtain the latent block mixture model with pdf

f (x; θ) = (z,w)∈Z×W i,j π z i ρ w j f z i w j (x ij ; α). (4)
Using above formulation, the randomized data generation process can be described by the threes steps row labellings (R), column labellings(C), and data generation (X) as follows:

(R) Generate the labellings z = (z 1 , . . . , z n ) according to the distribution π = (π 1 , . . . , π g ).

(C) Generate the labellings w = (w 1 , . . . , w d ) according to the distribution ρ = (ρ 1 , . . . , ρ m ).

(X) Generate for i = 1, ..., n and j = 1, ..., d a value x ij according to the (density) distribution f z i w j (.; α).

Model parameters estimation

EM based algorithms are the standard approach to estimate model parameters by maximizing observed data log-likelihood. In our case, the complete data is represented as a vector (x, z, w)

where unobservable vectors z and w are the labels. The complete data log-likelihood can then be written

L C (z, w, θ) = k z .k log π k + w . log ρ + i,j,k, z ik w j log f k (x ij ; α). (5)
The EM algorithm maximizes the log-likelihood L(θ) iteratively by maximizing the conditional expectation Q(θ, θ (c) ) of the complete data log-likelihood given a previous current estimate θ (c) and x:

Q(θ, θ (c) ) = i,k t (c) ik log π k + j, r (c)
j log ρ + i,j,k, e (c) ikj log f k (x ij ; α)
where

t (c) ik = P (z ik = 1|x, θ (c) ), r (c)
j = P (w j = 1|x, θ (c) ), e (c)
ikj = P (z ik w j = 1|x, θ (c) )
Unfortunately, difficulties arise owing to the dependence structure in the model, and more precisely, to determinate e (c) ikj . To solve this problem a new approximate solution is proposed [START_REF] Govaert | Clustering with Block Mixture Models[END_REF]] using the [START_REF] Hathaway | Another Interpretation of the EM Algorithm for Mixture Distributions[END_REF]] and [START_REF] Neal | A View of the EM Algorithm that Justifies Incremental, Sparse, and Other Variants[END_REF]] interpretation of the EM algorithm. By constraining the probability p to verify the relation p(z i , w j ) = p(z i )p(w j ) ∀i, j, it can be shown that the fuzzy clustering criterion proposed by [START_REF] Neal | A View of the EM Algorithm that Justifies Incremental, Sparse, and Other Variants[END_REF]] takes the form

FC (t, r; θ) = L C (t, r, θ) + H(t) + H(r) (6)
in case of latent block model. where

H(t) = - ik t ik log t ik and H(r) = - j r j log r j
and L c is the fuzzy complete-data log-likelihood associated to the block latent model given by

L C (t, r; θ) = k t .k log π k + . r log ρ + i,j,k, t ik r j log f k (x ij ; α). (7)
This is a variational approach and its objective is to maximize the function FC . Doing that, it has replace the maximization of the likelihood by the maximization of an approximation of this likelihood defined by

L(θ) = argmax t,r FC (t, r, θ).
To simplify some notations, we will also use the restricted complete data log-likelihood

L CR (z, w, θ) = i,j,k, z ik w j log f k (x ij ; α) (8)
and the restricted fuzzy complete-data log-likelihood

L CR (t, r; θ) = i,j,k, t ik r j log f k (x ij ; α). (9)

Algorithms

Here we give brief outline of two algorithms used in our package.

Block expectation maximization (BEM) algorithm

The fuzzy clustering criterion given in ( 6) can be maximized using EM algorithm. We here only outline the various expressions evaluated during E and M steps.

E-

Step: Here we compute conditional row and column class probabilities which are respectively given by:

t ik = log π k + j r j log f k (x ij ; α). ( 10
)
and d . Also the estimate of model parameters α will be obtained by maximizing i,j,k,

r jl = log ρ + i,k t ik log f k (x ij ; α) (11)
t ik r j log f k (x ij ; α) (12)
which will depend on the pdf's f k and shall be treated later separately in subsections discussing block mixture models for binary, contingency and continuous case.

Using the E and M steps defined above, BEM algorithm can be enumerated as follows:

1. Initialize t (0) , r (0) and θ (0) = (π (0) , ρ (0) , α (0) ).

2. Compute t (c+1) , π (c+1) , α (c+(1/2)) by using EM algorithm for the data matrix u i = j r j x ij and starting from t (c) , π (c) , α (c) .

3. Compute r (c+1) , ρ (c+1) , α (c+1) by using EM algorithm for the data matrix v jk = i t ik x ij and starting from r (c) , ρ (c) , α (c+(1/2)) .

4. Iterate steps ( 2) and (3) until convergence.

Block classification expectation maximization (BCEM) Algorithm

To apply the complete maximum likelihood (CML) approach for parameter estimation of the latent block model ( 4), the partitions z and w are added to the parameters to be estimated, and the objective is to maximize complete data log-likelihood associated with the latent block model (4). Unlike the BEM situation, this maximization does not require an approximation and can be done,for instance, by maximizing with fixed z and ρ and then with fixed w and π. This algorithm, can be seen as a variant of BEM algorithm described above and to perform this algorithm, it is sufficient to add a C-step in each of the EM algorithms, which converts the t ik 's and r jl 's to a discrete classification before performing the M-step by assigning each individual and variable to the cluster having the highest posterior probability.

Block mixture models for binary, contingency and continuous datasets

Binary data

The parameters α of the underlying distribution of a binary data set is given by the matrix p = (p k ) where p k ∈]0, 1[ ∀ k = 1, . . . , g and = 1, . . . , m and the probability distribution

f k (x ij ; p) = f (x ij ; p k ) is the Bernoulli distribution f (x ij ; p k ) = (p k ) x ij (1 -p k ) 1-x ij .
we re-parameterize the model density as follows:

f k (x ij ; α) = (ε kj ) |x ij -a k | (1 -ε kj ) 1-|x ij -a k |
where

a k = 0, ε k = p k if p k < 0.5 a k = 1, ε k = 1 -p k if p k > 0.5.
Hence the parameters p k of the Bernoulli mixture model are replaced by the following parameters:

The binary value a k , which acts as the center of the block k, and which gives, for each block, the most frequent binary value,

The value ε k belonging to the set ]0, 1/2[ that characterizes the dispersion of the block k, and which is, for each block, represents the probability of having a different value than the center.

For this model, using Equation 10 and 11 the row and column class conditional probabilities are respectively given by:

t ik = log π k + |u i -r . a k | log ε k 1 -ε k + r . log (1 -ε k ) and r j = log ρ + k |v jk -t .k a k | log ε k 1 -ε k + k t .k log (1 -ε k )
where

u i = j r j x ij and v jk = i t ik x ij .
Furthermore, we can write the restrained data log-likelihood for this model in terms of previously defined explicit parameters as

L CR (t, r, α) = k |y k -t .k r . a k | log ε k 1 -ε k + t .k r . log (1 -ε k ) ( 13
)
where y k = i,j t ik r j x ij . It can be easily shown that the values of model parameters maximizing Equation 13 is given by

a k = 0, ε k = y k t .k r . if y k t .k r . < 0.5 a k = 1, ε k = 1 -y k t .k r .

otherwise.

blockcluster: An R Package for Model Based Co-Clustering For detailed discussion on Bernoulli latent block models we would kindly refer our readers to [START_REF] Govaert | Block Clustering with Bernoulli Mixture Models: Comparison of Different Approaches[END_REF]].

Contingency data

To apply the block latent mixture model on contingency data, it is assumed that for each block k, , the values x ij are distributed according to Poisson distribution P(µ i ν j γ k ) where the Poisson parameter is split into µ i and ν j the effects of the row i and the column j respectively and γ k the effect of the block k [START_REF] Govaert | Latent Block Model for Contingency Table[END_REF]]. Then, we have

f k (x ij ; α) = e -µ i ν j γ k (µ i ν j γ k ) x ij x ij !
where α = (µ, ν, γ) with µ = (µ 1 , . . . , µ n ), ν = (ν 1 , . . . , ν d ) and γ = (γ 11 , . . . , γ gm ).

Unfortunately, this parameterization is not identifiable. It is therefore not possible to estimate simultaneously µ i , the ν j and γ k . The following two solutions are considered to tackle this problem:

µ i and ν j are known: In this case, it only remains to estimate the matrix γ k .

Constraints are imposed on the parameter θ = (π, ρ, µ, ν, γ) : In this situation, following Constraints

k π k γ kl = ρ γ kl = 1, i µ i = 1, j ν j = 1 (14)
appear to make the model identifiable, but it is yet remains to prove it.

The transformation of the parameter θ on a parameter θ defined as follows

π k = π k and ρ = ρ µ i = µ i N and ν j = ν j N γ k = γ k N ,
shows that it is possible to parameterize the model with the constraints

k π k γ kl = ρ γ kl = 1/N, i µ i = j ν j = N ( 15
)
where N is any positive real.

Case-1 : The row and column effects are known In this situation, α = γ = (γ 11 , . . . , γ gm ) and f k (x ij ; γ) can be written as

f k (x ij ; γ) = K ij .g k (x ij ; γ)
where

g k (x ij ; γ) = e -µ i ν j γ k (γ k ) x ij and K ij = (µ i ν j ) x ij x ij ! .
Then,we can replace the function f k by the function g k in the log-likelihoods and in the expressions t ik and r j . Hence, using Equation 10 and 11 the row and column class conditional probabilities are respectively given by:

t ik = log π k + j r j x ij log γ k -µ i γ k ν r j = log ρ + i,k t ik x ij log γ k -ν j k γ k µ k .
Also the restrained complete data log-likelihood is given by:

L CR (t, r, α) = k, y k log γ k -µ k ν γ k . ( 16
)
where y k = ij t ik r j x ij . The maximization of Equation 16gives:

γ k = y k µ k ν ∀k,
Case-2 : The row and column effects are unknown Using the constraints given by Equation 15, we obtain E(x i. ) = µ i and E(x .j ) = ν j . We can propose as an estimator of µ i and ν j margins x i. and x .j and take the results of case-1 to calculate row and column class conditional probabilities by replacing µ i and ν j by x i. and x .j respectively. Hence the row and column class conditional probabilities are respectively given by:

t ik = log π k + j r j x ij log γ k -x i. γ k y . r j = log ρ + i,k t ik x ij log γ k -x .j k γ k y k.
Furthermore the derivation of γ k is exactly same as in previous case. Using this particular form of γ k , the calculation of t ik and r j is simplified and becomes

t ik = log π k + j r j x ij log γ k -x i. r j = log ρ + i,k t ik x ij log γ k -x .j

Continuous data

In this case, the continuous data is modeled using unidimensional normal distribution. Hence the density for each block is given by:

f k (x ij ; α) = 1 2πσ 2 k exp -{ 1 2σ 2 k (x ij -µ k ) 2 }
The parameters of the model is α = (α 11 , . . . , α gm ) where α k = (µ k , σ 2 k ). For this model, using Equation 10 and 11 the row and column class conditional probabilities are respectively given by:

t ik = log π k - 1 2 (r . log σ 2 k + 1 σ 2 k a ik )
where

a ik = j r j (x ij -u i /r . ) 2 + r . ( u i r . -µ k ). 2 and r j = log ρ - 1 2 k (t .k log σ 2 k + 1 σ 2 k b jk ) where b jk = i t ik (x ij -v jk /t .k ) 2 + t .k k (v jk -µ k ). 2
The computation of parameters α is given by maximization of restrained log-likelihood which in this case is given by:

L CR (t, r, α) = - 1 2 k, (t .k r . log σ 2 k + 1 σ 2 k i,j t ik r j (x ij -µ k ) 2 ). ( 17
)
It can be deduced easily that the value of µ k and σ 2 k is given by minimization

of ij t ik r j (x ij - µ k ) 2 and t .k r . log σ 2 k + 1 σ 2 k ij t ik r j (x ij -µ k ) 2 ∀k
, respectively. Hence we obtain:

µ k = j r j v jk t .k r . , σ 2 k = ij t ik r j x 2 ij t .k r . -µ 2 k

Block clustering with package blockcluster

blockcluster is a newly developed R package for co-clustering of binary, contingency and continuous data. The core library is written in C++ and blockcluster API acts as a bridge between C++ core library and R statistical computing environment. A brief introduction of the C++ library is given in Appendix A. In the subsequent sections, we shall provide detailed explanation of the various functionalities and utilities of the package with help of examples wherever possible.

Strategy for parameters estimation

Correct parameters estimation in EM based algorithms has always been a challenging problem due to various local optima. In case of co-clustering algorithms, this problem could be even more apparent and hence the final results shall be directly affected by the choice of a strategy. As being explained in previous sections, we have two block clustering algorithms namely BEM and BCEM and the way we run these algorithms in our library is named respectively as 'XEMStrategy' and 'XCEMStrategy'. Both the strategies are equivalent except that they run different algorithms. The various steps involved in 'XEMStrategy' (equivalently in 'XCEMStrategy') are:

1. Set the pseudo log-likelihood L XEM to minimum possible value, 2. Run the BEM algorithm for some pre-defined number of times (we denote this number by 'nbxem' in our package) as follows: 3. Starting with the parameters α start , run the algorithm with a low value of epsilon (low tolerance) and a high number of iterations and calculate the current pseudo loglikelihood L current .

4. If L current is greater than L XEM , copy the currently estimated parameters to α f inal and set L XEM to L current .

5. Repeat step 2 through 3 for some pre-set times (we denote this number by 'nbtry' in our package) in order to obtain a good estimation of the parameters.

Although the above strategy does not explore all the possible solutions and hence cannot ensure to obtain global optimum, it is still a very promising approach in practice. The tuning of the values of 'nbxem' and 'nbtry' need to be done intuitively, and could have a substantial effect on final results. A good way to set these values is to run co-clustering few number of times and check if final log-likelihood is stable. If not, one may need to increase these values.

In practice, it is better to increment 'nbxem' as it could lead to better (stable) results without compromising too much the running time.

In the package blockcluster, we have a function called cocluststrategy 1 which allow to set the values of various input parameters including 'nbxem' and 'nbtry'. In the following example, we get an object defaultstrategy of S4 class strategy by calling the function cocluststrategy without any arguments and then we called the overloaded function summary to see various default values. 1 This function is also use to set all other input parameters like number of iterations and tolerance values. They are not described here explicitly for the sake of brevity. We refer our readers to the package documentation for details. One thing which is worth noting in the summary output (above) is that there is no default value for initialization method. It will be set automatically depending on the type of input model as explained in Section 3.2. Any of these input parameters can be set by passing appropriate arguments to function cocluststrategy as shown in example below where we set nbtry, nbxem and algo parameters. R > newstrategy <-cocluststrategy(nbtry=5, nbxem=10, algo='XCEMStrategy')

Model initialization

Initializing model parameters is arguably the most important step in getting the good final estimate of model parameters. The three most important initialization strategies for EM algorithms are: Random: In this method, the model parameters are initialized randomly using the input data. One may go with several random initializations and chose the one which gives the best log-likelihood.

CEM:

The model parameters are first initialized using random values by making use of input data, and then CEM algorithm [START_REF] Celeux | A Classification EM Algorithm for Clustering and Two Stochastic Versions[END_REF]] is run for small number of iterations with high tolerance. The CEM algorithm optimizes directly the completed log-likelihood by replacing the M step by a classification step based on MAP. In general, the CEM algorithm converges very quickly (not necessarily to good optimum) which is an advantage as we spend few time in initialization compared to the overall run time.

FuzzyCEM (Stochastic EM):

This method is equivalent to CEM except for the fact that instead of classifying samples into various clusters based on MAP, we stochastically select the class for each sample using the corresponding class probability distribution.

Unfortunately, in blockcluster, we do not have the freedom to initialize any model with all of the methods above. Instead, for every model we only have one possible initialization method at present (mainly CEM) as shown in Table 1. Hence the initialization strategy will be selected automatically by the library based on the input model. This is the reason why we have no default initialization strategy. This is indeed on the priority list of our future works to provide as many initialization strategies as possible for various kinds of models.

Examples with simulated datasets

Binary Data

We have simulated binary data using the data generation process explained in Section 2 and with parameters given in Table 2. The data consist of 1000 rows (samples) and 100 columns (variables) with two clusters on rows and three clusters on columns.

a, 0, 0.1 0, 0.3 1, 0.1 1, 0.3 1, 0.2 0, 0.1 π .6 .4 ρ .3 .3 .4

Table 2: Parameters for simulation of binary data.

The following R commands shows how to load the library, process the data and visualize/summarize results using blockcluster. abstracts) making a total of 2431 documents. Furthermore, we used all the words (excluding stop words) as features making a total of 9275 unique words. The data matrix consist of words on the rows and documents on the columns with each entry giving the term frequency, that is the number of occurrences of corresponding word in corresponding document. To run our package, we assume that the term frequency follows Poisson distribution. Hence we can apply the model 'pik_rhol_unknown' available in our package for contingency (Poisson Family) datasets with unknown row and column effects. Table 5 shows the confusion matrix and compare our results with classical bipartite spectral graph partitioning algorithm of [START_REF] Dhillon | Co-Clustering Documents and Words Using Bipartite Spectral Graph Partitioning[END_REF]] where we have obtained 100 percent correct classification. [START_REF] Dhillon | Co-Clustering Documents and Words Using Bipartite Spectral Graph Partitioning[END_REF]] (a), and Results using blockcluster (b). The difference in number of Cranfield documents is because we made use of the readymade data extracted from the documents and there are two less documents data in the same.

Conclusion and future directions

This article presents our newly developed R package for co-clustering algorithms based on latent block models. In the first part of article, we have presented brief mathematical foundation of co-clustering algorithms based on latent block models. In the second part, we have given details on the strategy adopted to run these algorithms. Furthermore, we have completed and providing his support in proposing the project Plan to INRIA-Lille. Furthermore we are thankful to Alhousseiny Pam who worked partly on the graphical part of the R package.

Finally we are grateful to INRIA-Lille for funding this project.

Appendices

A. The core library: Structure and usage

In this paper we have given detailed analysis on the usage of R package blockcluster to estimate parameters of latent block models, but behind the scenes the kernel (computational) part lies in the C++ core library2 . Our purpose to put this section is only to briefly introduce the kernel part and partly for developers who might be interested in integrating the library with their own codes, but we highly recommend to use only the R package blockcluster for research purposes as lot of utilities are absent in the C++ codes for instance the graphics part.

A.1. Library structure

The C++ library is based on famous Strategy design pattern [START_REF] Freeman | Head First Design Patterns[END_REF]] and (not so complete) architecture is revealed in Figure 7. The choice of the architecture is motivated from the fact that it allows easy extension to the library (accommodating any future research) with as little changes as possible in the existing codes. Hence it follows famous design principle "Open for extension but closed for modification" [START_REF] Freeman | Head First Design Patterns[END_REF]].

A.2. Usage

In our library, we have provided a simple main function (to illustrate its usage) that can be run from console as follows.

$ coclustmain -f path/to/optionfile.txt coclustmain is binary executable file which take path to a text file as command line argument. A sample input text file is shown below.

Figure 1 :

 1 Figure 1: Binary data set (a), data reorganized by a partition on I (b), by partitions on I and J simultaneously (c) and summary matrix (d).

  (a) Set the pseudo log-likelihood L xem to minimum possible value, (b) Initialize the model parameters, (c) Run the algorithm with a high tolerance (epsilon) and few iterations using the initialized parameters, and calculate the current pseudo log-likelihood L current , (d) If L current is greater than L xem , copy the currently estimated parameters to α start and set L xem to L current .

  *************************************************************** Algorithm: XEMStrategy Initialization method(There is no default value): Stopping Criteria: Parameter Various Iterations

*

  ***************** Number of global iterations while running initialization: 10 Number of iterations for internal E-step: 5 Number of EM iterations used during xem: 50 Number of EM iterations used during XEM: 500 Number of xem iterations: 5 Number of tries: 2 Various epsilons **************** Tolerance value used while initialization: 0.01 Tolerance value for internal E-step: 0.01 Tolerance value used during xem: 1e-04 Tolerance value used during XEM: 1e-10 ******************************************************************

  *************************************************************** Model Family : Bernoulli Latent block model Model Name : pik_rhol_epsilonkl Model Parameters..

Figure 5 :

 5 Figure 5: Original and co-clustered (segmented) image.

Figure 6 depicts the 2 ×

 2 Figure 6 depicts the 2 × 2 checkerboard pattern in the data matrix, hence confirming the more frequent occurrence of particular set of words in one document and vice-versa. Please note that the data matrix images are extremely sparse (data points almost invisible) and have been processed using simple image processing tools for visualization purpose only. Medline Cranfield Medline 1026 0 Cranfield 7 1400 (a)

  _ rhol _ epsilonklThis file contains all the input options needed to run the library. The only mandatory inputs are the data-type, number of row/column clusters and the path to data file. As the reader may have already guessed, this file is more or less similar to blockcluster package function cocluststrategy() where the optional parameters in the sample text file are equivalent to arguments of the function. By looking at the coclustmain source file, one can easily understand how to integrate the core library with their own C++ codes.

Figure 7 :

 7 Figure 7: C++ core library software architecture based on strategy design pattern

Table 1 :

 1 Various models available in package blockcluster.

	Model	Datatype	Proportions Dispersion/Variance Initialization
	pik rhol epsilonkl	binary	unequal	unequal	CEM
	pik rhol epsilon	binary	unequal	equal	CEM
	pi rho epsilonkl	binary	equal	unequal	CEM
	pi rho epsilon	binary	equal	equal	CEM
	pik rhol sigma2kl continuous	unequal	unequal	CEM
	pik rhol sigma	continuous	unequal	equal	CEM
	pi rho sigma2kl	continuous	equal	unequal	CEM
	pi rho sigma2	continuous	equal	equal	CEM
	pik rhol unknown contingency	unequal	N.A	CEM
	pi rho unknown	contingency	equal	N.A	CEM
	pik rhol known	contingency	unequal	N.A	Random
	pi rho known	contingency	equal	N.A	Random

Table 5 :

 5 Confusion Matrix: Results reported in

The source of the C++ library can be checkout anonymously from the repository svn://scm.gforge. inria.fr/svnroot/cocluster

The following R command is used to plot the original and coclustered data (Figure 2(a)) with aspect ratio set to false (it is true by default).

R > plot(out, asp = 0)

To Plot the various block distributions (Figure 2(b)), the following R command is used with type set to 'distribution' (type is 'cocluster' by default). R > plot(out, type = 'distribution')

Contingency data

In this case, we simulated contingency data as explained earlier in Section 2 and using the parameters given in Table 3. Again the data dimensions are 1000 × 100 with two clusters on rows and three clusters on columns.

Continuous data

In this final example, we generated continuous data using parameters given in Table 4. The data-set again contains 1000 rows and 100 columns with two clusters on rows and three clusters on columns. the analysis of blockcluster package by providing the examples on simulated datasets for each data-type. Finally we have shown applicability of our software on real datasets. As the package facilitates co-clustering for all types of data, it seems to provide a generic framework and a one stop solution for wide variety of co-clustering applications.

In the future release(s) of blockcluster package, we have planned to provide parallel implementations of the software by exploiting various data and task parallelism in latent block model estimation. Furthermore, we are also working on new initialization strategies as well as model selection criterion that should be part of upcoming releases of the package. Finally we have plans to extend the library by providing Stochastic EM algorithms for parameter estimation as well as to implement semi-supervised co-clustering in which some row and/or column classes are provided by the user.

Acknowledement

The authors are grateful to Christophe Biernacki for his leading support and constant motivation throughout the project. We are also thankful to Gilles Celeux for his contributions