Arthur Tenenhaus
email: arthur.tenenhaus@imed.jussieu.fr

Alain Giron

Emmanuel Viennet

Michel Béra

Gilbert Saporta

Bernard Fertil

Kernel logistic PLS: a tool for supervised nonlinear dimensionality reduction and binary classification

Keywords: Classification, Kernel, PLS regression, Logistic regression, Dimensionality reduction

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In this paper, we revisit the classification problem by focusing on kernel-based methods such as support vector machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF][START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF], kernel logistic regression (KLR) [START_REF] Zhu | Kernel logistic regression and the import vector machine[END_REF] and others such as kernel partial least squares (KPLS) regression [START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF][START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF][START_REF] Bennett | An optimization perspective on kernel partial least squares regression[END_REF]. In kernel-based methods, a complexity control step is usually performed. Vapnik-Chervonenkis (VC) theory provides a general framework for complexity control called structural risk minimization (SRM). VC theory relies on a nested set of models of increasing complexity (called a structure). In this context, the choice of a suitable structure is a decisive step. This paper proposes a novel way to build such a nested structure based on the partial least squares (PLS) regression. The PLS algorithm is especially convenient when handling high-dimensional and highly 1 correlated data [START_REF] Wold | The multivariate calibration problem in chemistry solved by the PLS method[END_REF][START_REF] Tenenhaus | La Régression PLS[END_REF], which is the very issue addressed when using kernels. This dual aspect (complexity control procedure and explicit linear model designed for high-dimensional and highly correlated data) makes the association of kernel theory with PLS regression (PLS-R) a clearly fruitful approach that has already given rise to many developments. For instance, [START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF] have proposed a nonlinear PLS algorithm via kernel: the KPLS regression.

PLS-R was not originally designed for classification. However, based on the algorithmic structure of PLS-R, two new linear approaches have been proposed: the PLS logistic regression [START_REF] Tenenhaus | La Régression Logistique PLS validée par bootstrap[END_REF][START_REF] Bastien | PLS generalized linear regression[END_REF][START_REF] Tenenhaus | Kernel logistic PLS: a new tool for complex classification[END_REF] and the PLS for discrimination [START_REF] Barker | Partial least squares for discrimination[END_REF]. Rosipal et al. have also suggested a kernelized version of the PLS for discrimination [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF].

Based on [START_REF] Bennett | An optimization perspective on kernel partial least squares regression[END_REF] work, we present a nonlinear extension of the PLS logistic regression based on the factorization of the kernel matrix: the Kernel Logistic PLS (KL-PLS).

The paper is organized as follows: in Section 2, we briefly review some principles of SRM and PLS-R. In Section 3, we focus on PLS for classification. In Section 4, we review some nonlinear extensions of PLS-R. In Section 5, we propose the KL-PLS algorithm and present numerical and graphical results in Section 6.

Partial least squares regression and structural risk minimization

Structural risk minimization (SRM)

Given a set of examples (x 1 ,y 1 ),...,(x n ,y n ), x i ∈ R p , y i ∈ R, drawn from an unknown distribution P(x, y), expected risk minimization is achieved by finding a function f opt ∈ F ={f , ∈ }, which provides the smallest possible expected risk: R( ) = V(y,f (x))P (x, y) dx dy, (1) where each function f in F is uniquely labeled by the adjustable parameter and V defines a cost function measuring the discrepancy between the y and f (x).

Usually, the probability distribution P(x, y) is unknown, preventing minimizing expected risk. However, a sampling of P(x, y) is available from which the empirical risk can be minimized:

R emp ( ) = 1 n n i=1
V(y i ,f (x i )).

(

) 2
Vapnik and Chervonenkis have provided an inequality connecting expected risk with empirical risk, with a probability of 1 -: R( ) R emp ( ) + h(log(2n/ h) + 1)log( /4) n ∀ ∈ ,

where h is the VC-dimension of F [START_REF] Vapnik | Statistical Learning Theory[END_REF]. The right-hand side of the inequality ( 3) is called the "structural risk". It is clear from inequality (3) that a small value of the empirical risk does not necessarily imply a small value of the expected risk. In order to make the expected risk small, both the empirical risk and the VC-dimension should be minimized at the same time. In the context of SRM, a set of possible models F ={f (x) : ∈ } is decomposed into a nested structure of subsets F k ={f (x) : ∈ k } of increasing complexity (or flexibility to fit the data), so that

F 1 ⊂ F 2 ⊂ ••• ⊂ F k ⊂ ..., and the (finite) VC-dimension h k of F k verifies h 1 h 2 ••• h k ...
. Thus, a model can be chosen among the nested sequence based on the smallest value of the upper bound of (3). Application of SRM, in practice, depends on a chosen structure. An example of generic structure is a dictionary representation where the set of approximating functions is

f m (x,c,w)= m i=1 c i b(x, w i ), (4)
where b(x, w i ) is a set of basis functions with adjustable parameters w i and c i are linear coefficients [START_REF] Cherkassky | Model complexity control for regression using VC generalization bounds[END_REF]. Both w i and c i are estimated from the training data. Representation (4) defines a structure, since

f 1 ⊂ f 2 ⊂ ••• ⊂ f k ⊂ ....
Hence the number of terms m in (4) specifies a particular f i of the structure. We formulate, in the next section, a practical way to build such a nested structure based on the PLS-R.

Partial least squares regression (PLS-R)

PLS-R [START_REF] Wold | The multivariate calibration problem in chemistry solved by the PLS method[END_REF]) is a technique for summarizing two data sets X and Y by latent variables (or PLS components), taking into account that the block Y is a set of responses and the block X is a set of predictors. PLS components, denoted by t 1 ,...,t m are related to X and constrained to be orthogonal, and PLS components, denoted by u 1 ,...,u m are related to Y and are not constrained to be orthogonal. Let X 0 and Y 0 be the centralized versions of X and Y. [START_REF] Höskuldsson | PLS regression methods[END_REF] has shown that the hth PLS components t h = X h-1 w h and u h = Y h-1 c h , are obtained by maximizing the [START_REF] Tucker | An inter-battery method of factor analysis[END_REF]

criterion: cov 2 (X h-1 w h ,Y h-1 c h ) = var(X h-1 w h )corr 2 (X h-1 w h ,Y h-1 c h ) var(Y h-1 c h ), (5)
subject to the constraints w h = c h =1. X h-1 is the residual of the regression of X 0 on t 1 ,...,t h-1 and Y h-1 is the residual of the regression of Y 0 on t 1 ,...,t h-1 . The solution w h (respectively, c h ) of this problem is obtained from the normalized eigenvector of the matrix

X t h-1 Y h-1 Y t h-1 X h-1 (respectively, Y t h-1 X h-1 X t h-1 Y h-1
) associated with the largest eigenvalue.

Therefore, PLS creates orthogonal latent variables t 1 ,t 2 ,...,t m (PLS components), linear combinations of the p original variables x 1 ,...,x p using Y for their determination. Let us focus on the one-dimensional output case. The PLS model with m components fits expansion (4) as follows:

f m (x,c,w)= m i=1 c i b(x, w i ) = m h=1 c h t h = m h=1 c h p j =1 x h-1,j w hj = m h=1 c h p j =1 x j w * hj , (6)
where w h = arg w h max cov(X h-1 w h ,y), c h is the regression coefficient of t h in the regression of y on t 1 ,...,t m and w * h is defined in Section 5.5. It is worth noting that y is integrated in the construction of the nested structure. Subsequently, the nested structure of f i is built in a supervised way and complexity control is achieved here by supervised dimensionality reduction.

Linear PLS methods for classification

Although PLS was not inherently designed for classification, it is routinely used for this purpose. PLS discriminant analysis (PLS-DA) is usually defined as a PLS-R of the dummy matrix derived from y on the set of predictors X. Thus PLS-DA is a technique for modeling linear relation between a set of binary variables Y on the set of predictors X.In addition, two linear PLS methods have been designed especially for classification: PLS for discrimination and PLS logistic regression.

PLS for discrimination (PLS-D)

In this section, we assume that Y is the dummy matrix. [START_REF] Barker | Partial least squares for discrimination[END_REF] suggested to remove the not meaningful Y-space penalty, var(Y h-1 c h ), in criterion (5). They have noted that this penalty is, in fact not adequate for classification and proposed instead to obtain the PLS-D components t h = X h-1 w h and u h = Yc h , h = 1,...,m,by maximizing criterion:

var(X h-1 w h )corr 2 (X h-1 w h ,Yc h ), (7)
subject to the constraints w h = c h =1. X h-1 is the residual of the regression of X 0 on t 1 ,...,t h-1 . The solution w h is obtained from the normalized eigenvector of the matrix X t h-1 Y(Y t Y) -1 Y t X h-1 associated with the largest eigenvalue. We may note that t h = X h-1 w h is the first principal component of the PCA of the X h-1 variables projected onto the Y variable space.

As far as one-dimensional output is concerned, the PLS components extracted from PLS-D and PLS-DA are identical. This follows from the fact that Y t Y is a scalar. PLS-DA and PLS-D only extract PLS components and a supervised method (such as the logistic regression) applied on selected components, provides a classification rule.

PLS logistic regression (PLS-LR)

When the set of responses Y (=y) is limited to one variable, the algorithm can be simplified because the ranks of the matrices X t yy t X and y t XX t y are equal to one. The following results are obtained: [START_REF] Garthwaite | An interpretation of partial least squares[END_REF] has noticed that each coordinate of the weight vector w 1 can be written as

w 1 = X t y/ X t y , t 1 = Xw 1 = XX t y/ X t y .
w 1j = x t j y/ X t y = x t j x j x t j y x t j x j • X t y ∝ s 2 j a j j s 2 j a j ,
where s 2 j is the variance of x j and a j the regression coefficient of x j in the regression of y on x j . Therefore, the weight vector w 1 is obtained by a succession of simple regressions of the dependent variable y on each explanatory variable x j . This result remains valid for the next components (w h is obtained by a succession of multiple regressions of y on t 1 ,...,t h-1 and each residual variable x h-1,j ).

Using the [START_REF] Garthwaite | An interpretation of partial least squares[END_REF] approach and replacing the successions of simple regressions by simple generalized linear regressions, [START_REF] Bastien | PLS generalized linear regression[END_REF] extended the PLS-R to the PLS generalized linear regression (PLS-GLR). The PLS-GLR algorithm provides a natural tool for classification through the logistic regression (since binary logistic regression belongs to the generalized linear model family) giving rise to the PLS logistic regression (PLS-LR) [START_REF] Tenenhaus | La Régression Logistique PLS validée par bootstrap[END_REF][START_REF] Tenenhaus | La régression logistique PLS[END_REF][START_REF] Bastien | PLS generalized linear regression[END_REF]. The PLS-GLR algorithm is described in Appendix A.

The PLS-GLR of a response y on a set of predictors x 1 ,...,x p with m components is written as

g( ) = m i=1 c h ⎛ ⎝ p j =1 w * hj x j ⎞ ⎠ , ( 8
)
where the parameter is the mean of a numerical variable. The link function, g, is chosen by the user according to the probability distribution of y.

Below are quoted four important properties of PLS-LR:

(i) Contrary to many data analysis methods, it needs only to compute the inverse of small-dimensional matrices (size equal to the number of PLS-LR components). The number of components is generally small and usually chosen by cross-validation. (ii) It allows handling data with more variables than observations. (iii) It allows handling highly correlated variables. (iv) Eqs. ( 8) and ( 4) have a similar expansion so that PLS-LR (as well as PLS-R) fits the SRM paradigm. PLS-DA, PLS-D and PLS-LR are designed to operate with high-dimensional and highly correlated data. The ability of these methods to do inference in high-dimensional space makes them ideal candidates for a nonlinear analysis of data based on an increase of their dimension.

Nonlinear PLS methods

Several methods are devoted to the extension of the linear PLS methods to the nonlinear case. Focusing here on binary classification we briefly review pros and cons of some of them, in order to precisely substantiate our contribution.

The simplest approach

One simple way of introducing nonlinearity in the model is to increase the dimension of the input matrix X by using higher order and cross terms of the original input variables, and then to use this new matrix X in the PLS algorithm [START_REF] Wold | Non-linear PLS modelling[END_REF]. However, such an increase is workable only if p is originally very small. It is worth noting that the interpretability of the final model (coefficients and variables) is preserved.

The Wold's approach

Without discarding PLS as a tool for dimensionality reduction, the linear inner relation

ŷ = m h=1 c h p j =1 w * hj x j = m h=1 c h t h , (9)
may be replaced with a nonlinear relation yielding the following general form of ŷ:

ŷ = m h=1 g h ⎛ ⎝ p j =1 w * hj x j ⎞ ⎠ = m h=1 g h (t h ), (10)
where g h 's are smooth nonlinear functions. This was first proposed by Wold using quadratic polynomials without cross terms [START_REF] Wold | Non-linear PLS modelling[END_REF]) and further developed with smooth functions approximated by splines [START_REF] Wold | Non-linear partial least squares modelling. II Spline inner function[END_REF] or neural networks [START_REF] Baffi | Non-linear projection to latent structures revisited (the neural networks PLS algorithm)[END_REF]. The initial PLS dimension reduction is particularly useful when a set of variables with no discriminant information jeopardizes data segmentation. However, if the discriminant information relies on nonlinear relationships between variables and target, the initial PLS dimension reduction may have a disastrous effect.

KPLS approach

Kernel-based learning algorithms [START_REF] Cortes | Support vector network[END_REF][START_REF] Schölkopf | Learning With Kernel-Support Vector Machines Regularization Optimization and Beyond[END_REF][START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF] proceed by mapping original data X into a (possibly high-dimensional) Hilbert space F, usually called feature space, and then by using linear pattern analysis to detect relations in the feature space. Mapping is performed by specifying the inner product between each pair of data rather than by explicitly calculating their coordinates (the so-called "kernel trick"). This approach has several advantages, the most important being that the inner product in the feature space is often much more easily computed than the coordinates of the points (notably when the dimensionality of the feature space is high). Given an input space X and a feature space F, we consider a function : X → F. Given two points u ∈ X and v ∈ X, the kernel function returns the inner product between their images mapped by . The choice of the kernel function (and associated parameter(s)) defines the relative position of the data in the feature space. The linear kernel (k(u, v) = u, v ) and the Gaussian kernel

(k(u, v) = exp(-u -v 2 /2 2
)) are two of the most popular kernel functions.

The nonlinear KPLS regression method is based on the mapping of the original input data X into a high-dimensional feature space F. The objective is to construct a linear PLS-R in F or equivalently a nonlinear PLS-R model in X. Rosipal and Trejo developed a nonlinear form of the PLS-R based on the kernel trick, the so-called KPLS regression [START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF]. In the case of one-dimensional output, the first PLS component t PLS 1 is proportional to XX t y and thus, PLS can be rigorously adapted to the kernel approach, giving rise to the KPLS.

Denote the matrix of mapped input data into the feature space F.

t KPLS 1 ∝ t y. (11)
Instead of an explicit nonlinear mapping, the kernel trick can be used.

t KPLS 1 ∝ t y = Ky, ( 12
)
where K is the n × n kernel matrix (see Section 5.1).

In Section 2.2, we assumed a centralized PLS-R problem. To centralize the mapped data in a feature space F we can apply (13) [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF][START_REF] Rosipal | Kernel partial least squares regression in reproducing kernel hilbert space[END_REF]:

K 0 = I n - 1 n 1 n 1 t n K I n - 1 n 1 n 1 t n , ( 13
)
where I n is an n × n identity matrix and 1 n represents a n × 1 vector of ones. Like PLS, KPLS is an iterative process. After extraction of the unit norm t KPLS h-1 , the (h -1)th KPLS component, the algorithm is reiterated using the following deflated matrix:

h 2,K h-1 = (I -t KPLS h-1 t KPLS t h-1 )K h-2 (I -t KPLS h-1 t KPLS t h-1 ). (14)
Eq. ( 14) is based on the fact that h-1 is deflated as

h 2, h-1 = h-2 -t KPLS h-1 t KPLS t h-1 h-2 , ( 15
)
and the hth unit norm KPLS component is defined by

t KPLS h = K h-1 y. ( 16
)
Like PLS, KPLS is not designed for classification purposes. Therefore, [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF] have suggested a kernelized version of the PLS for discrimination (KPLS-D). As it is for the linear framework, KPLS-D and KPLS-DA for binary classification are equivalent. However, KPLS-D provides an appealing solution for the multiclass classification [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF] and should be preferred to KPLS-DA.

It is worth pointing out that KPLS allows avoiding drawbacks of the two first nonlinear PLS methods described above. Indeed, PLS component in the feature space allows capturing nonlinear relation between X and y. Moreover, it is "only" necessary to handle a n × n kernel matrix and KPLS limitations are much more related to the number of observations (n) than to the number of variables (p). However, the interpretability of the model is lost.

Direct kernel PLS (DK-PLS)

Bennett and Embrechts (2003) introduced the DK-PLS regression. They suggested applying PLS-R of y on K 0 in order to obtain a low rank approximation of the kernel matrix (suitable for regression), used to obtain the final linear model. DK-PLS setting is based on the so-called "empirical kernel map" (also called "direct kernel mapping") [START_REF] Schölkopf | Learning With Kernel-Support Vector Machines Regularization Optimization and Beyond[END_REF][START_REF] Webb | Nonlinear feature extraction with radial basis functions using a weighted multidimensional scaling stress measure[END_REF]. Columns of the kernel matrix are the new input variables. It is worth noting that DK-PLS is a particular case of the first approach (Section 4.1). In addition, when using the kernel matrix, DK-PLS allows avoiding drawbacks of the two first approaches (cf. Sections 4.1 and 4.2). However, limitation of the third approach (cf. Section 4.3) related to the number of observations is not avoided, although combining empirical kernel map with sampling of the columns of K may help escaping it. We highlight here the flexibility of the empirical kernel map. Finally, like KPLS, interpretability of the model is lost.

In the case of a one-dimensional output, the first DK-PLS component t DKPLS , is proportional to K 0 K t 0 y and K 0 being symmetric:

t DKPLS 1 = K 2 0 y. ( 17
)
In a similar way, the hth DKPLS component is defined by

t DKPLS h = K 2 h-1 y, (18)
where K h-1 is the deflated matrix obtained by the following procedure:

K h = K h-1 -t DKPLS h t DKPLS t h K h-1 . ( 19
)
KPLS and DK-PLS are strongly related. In fact, for one-dimensional output response, it is easy to show that the PLS-R of y on (KPLS) is the PLS-R of y on K 1/2 (DK-PLS).

Kernel logistic PLS ( KL-PLS )

Noting that KPLS and DKPLS are competitive with other kernel-based classifiers, we propose a nonlinear classification algorithm that combines the flexibility of the empirical kernel map with the supervised properties of the PLS logistic regression: the KL-PLS. KL-PLS can be viewed as a low-rank approximation of the kernel matrix oriented for classification. The main idea of KL-PLS is to look for a discriminant space spanned by the KL-PLS components (t 1 ,...,t m ), where a simple model, such as the logistic regression may become efficient for classification. We use the kernel matrix to map the data into a higher-dimensional space where the probability of finding the hyperplane increases with the dimension of the space.

Therefore, KL-PLS computes orthogonal latent variables in the space induced by the kernel matrix before performing logistic regression in the derived feature space.

KL-PLS is consequently a 3-step algorithm:

Step 1: Computation of the kernel matrix.

Step 2: Computation of the KL-PLS components.

Step 3: Logistic regression of Y on the m retained KL-PLS components.

Construction of the kernel matrix

Let X be an n × p matrix representing n observations on p explanatory variables x k , k = 1,...,p and y ∈{0, 1} a binary variable (the target) observed on the n units. Let K be a kernel matrix associated with X. K may be generated by the kernel function given in Section 4.3. The dimension of the kernel matrix is n × n. Each cell k ij is an inner product between the individuals i and j (holds a measure of similarity) in the feature space F. Each column k j ,j = 1,...,nof K, measures similarities between individual j and the whole data set.

Construction of the first KL-PLS component

The first KL-PLS component t 1 provides the first discriminant axis. This component is built by performing a succession of n simple logistic regressions:

Step 1: Computation of the coefficients a 1j of k j in the binary logistic regression of y on each k j , j = 1,...,n.

Step 2: Normalization of the column vector a 1 made by a 1j 's:

w 1 = a 1 / a 1 .
Step 3: Computation of the first KL-PLS component as

t 1 = Kw 1 .
Consequently, there is no matrix inversion during the computation of the first KL-PLS component.

Deflation procedure

Let us assume that the KL-PLS components t 1 ,...,t h-1 are already computed. This step is designed to obtain components orthogonal to t 1 ,...,t h-1 , holding residual information on y. The hth KL-PLS component is subsequently computed from the residual of the regression of each k j , j = 1,...,non t 1 ,...,t h-1 .

Step 1: Computation of the residual k h-1,1 ,...,k h-1,n from the multiple regression of each k j , j = 1,...,n on t

1 ,...,t h-1 . Let K h-1 =[k h-1,1 ,...,k h-1,
n ] be the residual matrix.

The residual matrix (or the deflated matrix) is performed by computing a succession of n regressions on a small number of uncorrelated variables.

Computation of the hth KL-PLS component

The hth KL-PLS component has to capture the discriminant information not available in the h -1 previous ones. We compute the hth KL-PLS component by performing n logistic regressions on h uncorrelated variables.

Step 1: Computation of the coefficients a hj of k h-1,j in the binary logistic regression of y on t 1 ,...,t h-1 and each k h-1,j , j = 1,...,n.

Step 2: Normalization of the column vector a h made by a hj 's: w h = a h / a h .

Step 3: Computation of the hth KL-PLS component as t h = K h-1 w h .

Step 4: Expression of the component t h in terms of K as t h = Kw * h . Consequently, during the construction of the hth KL-PLS component, we only compute the inverse of n h-dimensional matrices.

Expression of KL-PLS components in terms of original variables

Expression of KL-PLS components, in terms of original variables, is a fundamental step to analyze new data. Let Ktest be the new data set. The product Ttest = Ktest × W * allows to compute the values of the KL-PLS components for the new data set.

Computation of w * h :

(a) The first KL-PLS component is already expressed in terms of original variables: t 1 = Kw 1 and w * 1 = w 1 . (b) The second KL-PLS component is expressed in terms of residuals in the regression of the original variables on t 1 . Let p 1j be the regression coefficient of t 1 in the regression of k j on t 1 . From K = t 1 p t 1 + K 1 and t 2 = K 1 w 2 we get

t 2 = K 1 w 2 = (K -t 1 p t 1 )w 2 = (K -Kw 1 p t 1 )w 2 = K(I -w 1 p t 1 )w 2 = Kw * 2 .
(c) In a similar way, it can be shown that t h is expressed in terms of the original variables as

t h = K h-1 w h = K - h-1 i=1 t i p t i w h = K - h-1 i=1 Kw * i p t i w h = K I - h-1 i=1 w * i p t i w h = Kw * h .

Final model

We perform a binary logistic regression of y on the m retained KL-PLS components:

P(y = 1/X = x) = e 0 + m h=1 h t h 1 + e 0 + m h=1 h t h , ( 20
)
where the choice of the dimension m is guided by model selection techniques such as cross validation 1 . As mentioned in Section 5.5, the KL-PLS components t 1 ,...,t h are linear combinations of the columns k j of K. Thus, Eq. ( 20) can also be rewritten as

P(y = 1/X = x) = e 0 + m h=1 h t h 1 + e 0 + m h=1 h t h = e 0 + n j =1 j k j 1 + e 0 + n j =1 j k j . ( 21
)
Consequently, in order to compute P(y = 1/X = x new ), for a new observation x new , it is worth noting that there is no need to compute the coordinates of x new in the space spanned by the KL-PLS components.

Kernel storage properties

Strategies for improving kernel methods through factorization have been receiving increasing attention [START_REF] Bach | Predictive low-rank decomposition for kernel methods[END_REF][START_REF] Smola | Sparse greedy matrix approximation for machine learning[END_REF][START_REF] Williams | Effect on the input density distribution on kernel-based classifiers[END_REF]. In fact, for example, KPLS requires storing in memory the full kernel matrix, which is a limiting factor for large training data set. Factorization appears as a very fruitful line of research. Below are some facts that impact kernel storage for KL-PLS:

(i) During the construction of the latent variables, KL-PLS only requires an individual column of the kernel matrix. (ii) As suggested by Bennett and Embrechts, the introduction of an intercept when constructing the latent variables (step 1) avoids the kernel centering method proposed by [START_REF] Wu | The kernel PCA algorithms for wide data. Part II: fast cross validation and application in classification of NIR data[END_REF]. (iii) The succession of simple regressions used within KL-PLS algorithm provides an efficient deflation of the kernel matrix.

Consequently, it is not necessary to store the entire kernel matrix in memory but only columns one at a time. Such an implementation, however, significantly slows down the algorithm because it is necessary to reevaluate the Gram Matrix for each KL-PLS component. However, the KL-PLS components may be calculated from a subset of variables resulting from a sampling of the columns of K when the analysis of large dataset is of concern. Sampling scheme could be based on the Wald test or the likelihood ratio test: only variables that are significantly related to y in the simple logistic regression of y on each k j , j = 1,...,n, contribute to the construction of K. This approach is not explored in the current paper.

Complete separation

Let us investigate three data configurations: complete separation, quasi-complete separation and overlap. According to [START_REF] Allison | Logistic regression using the SAS system: theory and application[END_REF], there is a complete separation of data points if a vector w such as w t x j > 0 for y j = 1 w t x j < 0 for y j = 0 exists and quasi-complete separation of data points if a vector w such as w t x j 0 for y j = 1 w t x j 0 for y j = 0 exists. The maximum likelihood estimate exists only in the overlap configuration. Therefore, when applying logistic regression on the KL-PLS components, it is important to care about the data structure. For a complete or quasicomplete separation, Linear Discriminant Analysis or the Heinze and Schemper approach [START_REF] Heinze | A solution to the problem of separation in logistic regression[END_REF] consisting of reducing the bias of maximum likelihood using [START_REF] Firth | Bias reduction of maximum likelihood estimates[END_REF] technique, could be used. The reader is referred to [START_REF] Allison | Logistic regression using the SAS system: theory and application[END_REF] and [START_REF] Albert | On the existence of maximum likelihood estimates in logistic regression models[END_REF] for details on problems arising from (quasi-) complete data separation in the logistic regression framework.

Validation

In this section, we present experimental results for 11 benchmark data sets and three medical classification problems.

Benchmarks

KL-PLS efficiency is tested on 11 benchmark data sets (two-class discrimination) used in [START_REF] Rätsch | Soft margin for adaboost[END_REF] and available at http://ida.first.fraunhofer.de/projects/bench/benchmarks.htm. Details can be found in [START_REF] Rätsch | Soft margin for adaboost[END_REF]. Each data set consists of 100 different training and testing partitions. Dimension of data, training sample size and test sample size of each benchmark are presented in Table 1.

Different methods have been performed for each data set and each partition. Mean and standard deviation of the test set error rate of each method have been computed using results of the 100 runs (for KLR, Zhu and Hastie computed mean and standard deviation using only 20 partitions). These values are reported in Table 2. [START_REF] Rätsch | Soft margin for adaboost[END_REF] performed SVMs. [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF] evaluated Kernel PLS-SVC (KPLS-SVC). KPLS-SVC is based on the KPLS for discrimination as dimensionality reduction method followed by SVM on retained KPLS components for the final classification. [START_REF] Zhu | Kernel logistic regression and the import vector machine[END_REF] proposed the KLR. [START_REF] Rätsch | Soft margin for adaboost[END_REF], KPLS-SVC [START_REF] Rosipal | Kernel PLS-SVC for linear and nonlinear classification[END_REF], KLR [START_REF] Zhu | Kernel logistic regression and the import vector machine[END_REF] Mean rank 3.12 .91 .91 .9

The last column provides the width of the Gaussian kernel and the number of retained KL-PLS components. For each data set, the methods were ordered from rank 1 to rank 4 according to their performances; the last row of the table provides the mean rank for SVM, KLR, KPLS-SVC and KL-PLS. The star indicates that results of the concerned method and KL-PLS are statistically different (p < 0.05) (pairwise T-test for all comparisons except for KLR (two-sample T-test)).

The methods presented above are compared to KL-PLS. Gaussian kernel is used throughout this study because of its flexibility. Thus, KL-PLS efficiency relies on the value of width of the Gaussian and the number of retained KL-PLS components. These values have been selected based on the minimum classification error observed after cross validation on the first five training sets.

A pairwise T-test is used to compare KL-PLS to the other classifiers for each data set, except for KLR. A less powerful two-sample T-test is used instead, since individual classification errors are not available for that classifier. The significance level for all tests has been set to = 0.05.

KL-PLS vs. SVM:

The classification error rate of KL-PLS is lower for 10 data sets out of 11. The null hypothesis is rejected for five of these data sets (banana, diabetes, ringnorm, twonorm, and waveform). Conversely, the classification errors rate of SVM is lower for one data set (flare-solar) and the null hypothesis is rejected in this case.

KL-PLS vs. KPLS-SVC:

The classification error rate of KL-PLS is lower for four data sets out of 11. However, the null hypothesis is never rejected. Conversely, the classification errors rate of KPLS-SVC is lower for six data sets and the null hypothesis is rejected for banana, breast-cancer, flare-solar, twonorm, waveform.

The good behavior of KPLS-SVC may be related to the use of the linear SVM model, which provides a boundary following the maximal margin principle. However, the quality of results relies on a larger set of free parameters (one additional parameter is required for the SVM step), which may be bothersome in practice.

KL-PLS vs. KLR: The classification error rate of KL-PLS is lower for eight data sets out of 10 and the null hypothesis is rejected for Ringnorm and Waveform. Conversely, the classification errors rate of KLR is lower for two data sets and the null hypothesis is rejected for Banana.

Ranks of performances over all data sets offer another way to compare models (Table 2). It seems that models fall into 2 groups for which the complexity control differs: KPLS-SVC and KL-PLS-which are supervised dimensionality reduction-based models-most often provide the best results as compared to SVM and KLR for which the complexity control is achieved by Tikhonov regularization [START_REF] Tikhonov | Solution of Ill-Posed Problem[END_REF][START_REF] Evgeniou | Regularization networks and support vector machines[END_REF][START_REF] Wahba | Support vector machines, reproducing kernel Hilbert spaces and randomized GACV[END_REF].

It is worth noting that both KLR and KL-PLS are based on the logistic regression model and that the main difference concerns the complexity control scheme. To what extent difference in performances relies on the complexity control method remains to be studied.

From these experiments, we can conclude that KL-PLS works at least as well as the other methods. It is remarkable that best methods achieved very comparable results on the benchmark data sets. In fact, there may be an irreducible level of error (linked to the data structure) already reached by some of them.

To conclude, the enlarged feature space generated by the use of kernel-based methods gives flexibility, and complexity control allows good generalization. Results presented on these 11 benchmark datasets highlight the interest of the supervised dimensionality reduction approach.

Banana data set

Banana is an artificial 2-dimensional data set with two classes (400 training set, 4600 testing set) allowing a visual inspection of the data separation in the original space. According to the 2-dimensional representation, it is essentially a nonlinear classification problem. Subsequently, it seems interesting to examine classification results in more details. A nice linear separation of the two classes can be observed in the derived feature space showing that a linear logistic regression is adequate to achieve an efficient classification.

Banana data projection onto the two first components found by KL-PLS

KL-PLS and probability

KL-PLS offers a natural estimate of the conditional probability, p(x) = P(y = 1/X = x) by using logistic regression as final function. Fig. 2 depicts banana data onto the original space (left) and the two first KL-PLS components (right). Colors symbolize conditional probability values p(x) and lines symbolize isoprobability contours.

Linear decision boundaries in the space induced by the KL-PLS components correspond to complex nonlinear boundaries in the original space.

Banana test error as a function of the number of KL-PLS components

Fig. 3 shows performances of the KL-PLS model as function of the number of retained KL-PLS components. For banana data set, we note that even when we consider an important number of KL-PLS components, there is no overfitting. Moreover, the tuning of KL-PLS is quite easy. In addition to the kernel parameters, there is only one intuitive parameter (with a small number of discrete values) to adjust as compared to the continuous parameter (C) for SVM for example.

Kernel logistic PLS and robustness

Robustness is defined here as the quality of the model to tolerate noise. In this experiment, the training variable y is degraded by flipping the labels of randomly selected observations while preserving class balance. The test set classification accuracy is then observed as a function of the percentage of noise injected in the target. Fig. 4 shows that KL-PLS performs as well as SVM light [START_REF] Joachims | Making large-scale SVM learning practical[END_REF] in this example.

The test set classification accuracy curves for SVM and KL-PLS have a similar shape. It is remarkable that the accuracy is almost unchanged up to 30% degradation.

KL-PLS and high-dimensional data

The ovarian cancer and the lung cancer datasets from [START_REF] Li | Kent Ridge Biomedical Dataset Repository[END_REF] were chosen to test the efficiency of KL-PLS for high-dimensional data. Gene expression profiles obtained from DNA micro-arrays are used to classify different types of tumors. In these situations, there is generally a huge number of variables and a small number of observations: ovarian (respectively, lung) is a 253 × 15154, (respectively, 181 × 12 533) data set with two classes. The evaluation As it frequently happens for high-dimensional problems (e.g. text mining, genomic), a linear kernel is effective in performing an efficient discrimination. Therefore, the number of KL-PLS components is the unique parameter to adjust. Here, its value is based on the minimum classification error observed by cross validation on the first five training sets. It is set to 19 for ovarian and six for lung. Linear SVM classifier results were obtained by [START_REF] Shen | Dimension reduction-based penalized logistic regression for cancer classification using microarray data[END_REF]. Results of KL-PLS and SVM are reported in Table 3.

To test whether KL-PLS is statistically more accurate than SVM on ovarian and lung dataset, T-test between the means of the testing errors of the two classifiers is used. From these experiments, it appears that KL-PLS outperforms SVM. Nevertheless, we can note that the error rates observed are low for both methods and are always due to a few observations only.

We may emphasize in these examples that the n>p configuration is particularly convenient for KL-PLS. The construction of the KL-PLS components is computationally efficient because it is obtained from a "small" n × n kernel matrix.

Fig. 5 shows the test error rate and the CPU time of KL-PLS as a function of the number of retained KL-PLS components for ovarian cancer and lung cancer.

Horizontal lines provide the CPU time of SVM light . For lung cancer, the CPU time of KL-PLS with six components (Matlab 6.5 implementation) is about five times lower than the CPU time of SVM light (C implementation). The CPU time of KL-PLS with 19 components for ovarian cancer is greater than the CPU time of SVM light . Nevertheless, at equal error rate (obtained with 14 components), the CPU time of KL-PLS is 30% lower than the CPU time of SVM light .

Melanoma detection from digital medical imaging

In this section, we applied KL-PLS to a medical problem based on skin tumor detection (benign lesion vs. melanoma). From 227 images of pigmented skin lesion provided by dermatology departments of British Hertfort Hospital and Louis Mourier Hospital, we have extracted five geometric and 38 colorimetric features. Each lesion is thus characterized by 43 variables. Thirty-two melanoma and 195 benign lesions make up the database. The objective of this application is to provide a diagnosis based on KL-PLS (melanoma detection) that can be compared with the diagnoses of five senior dermatologists. The histology (melanoma lesion or benign lesion) is the variable to predict (gold standard).

In medical diagnosis applications, accuracy is not the main issue. Sensitivity and specificity (or predictive positive and negative values) are far more relevant because of different misclassification costs. Fig. 6 provides, in addition to dermatologist's performances, the ROC curve of KL-PLS. See [START_REF] Hanley | Receiver operating characteristic methodology: the state of the art[END_REF] for details on the ROC curve.

Each dermatologist is characterized by two points: the first point (respectively, the second point) M i (respectively, E i ) corresponds to the accuracy of diagnosis (respectively, accuracy of excision decision) for dermatologist i (in terms of sensitivity and 1-specificity). We may observe that KL-PLS reaches performances close to the dermatologist's performances in terms of diagnosis and excision decision. Results are obtained by leave-one-out cross validation technique. We have retained three components with = 100. KL-PLS provides, for each tested lesion, a probability to be a melanoma and the ROC curve is built from these probabilities.

We highlight that the output probability is useful here because it gives rise to the possibility to reproduce the dermatologist's behavior. Indeed, diagnostic and therapeutic decisions correspond to a modulation of the sensitivity level. This decision flexibility is very important for the melanoma detection problem because of the unbalanced misclassification cost: sensitivity is far more important than specificity since disregarding a melanoma is a severe error.

Conclusion

This work demonstrates competitiveness of KL-PLS with other state-of-the art classification methods for 11 benchmark discrimination data sets as well as for three medical classification problems. The algorithm is simple to implement since it is mainly composed of ordinary least squares and logistic regressions. The complexity control achieved by the supervised dimensionality reduction performed by KL-PLS offers a potentially fruitful alternative to Tikhonov regularization based methods. Indeed, the projection of the high-dimensional data onto a small number of KL-PLS components allows visual inspection of data layout. It provides clues about data structure and selection of kernel parameter(s). Visualization of the spatial relationships between classes helps evaluating quality of classification and feature space; Are classes linearly separables? How well-separated and coherent are they? Finally, KL-PLS offers an estimate of the conditional probability, often of interest itself, in medical field for example.

We would like to point out that the direct factorization of the kernel matrix in the framework of KL-PLS presents two significant advantages:

(i) Rectangular kernels can be factorized by KL-PLS. Derived from the fact that KL-PLS has only to handle a n × n matrix, the algorithmic limitation is much more connected to the number of observations than to the number of variables. Therefore, this approach allows management of very high-dimensional data. Large training set size can also be managed, thanks to the sampling that can be performed within KL-PLS with no modification of the algorithm.

(ii) The requirement of semi-definiteness may rule out the most natural pairwise similarity functions for a given problem. In the KL-PLS setting, only similarity measures are required and there is no need defining a dot product in the feature space; the kernel matrix need not verify the Mercer's conditions (positive definite).

Extension of KL-PLS to the multi-class problem is easy for ordinal output y (see Appendix B) and more complicated for categorical output (work in progress). Validation of the Kernel Generalized PLS, a nonlinear version of the PLS-GLR that could be derived from a generalization of DKPLS, is an interesting perspective.

1

 1

Fig. 1

 1 depicts projection of the original training (left) and testing data (right) onto the two first KL-PLS components (training data).

Fig. 1 .

 1 Fig. 1. Banana data depicted onto the two first KL-PLS components.

Fig. 2 .

 2 Fig. 2. Banana depicted onto the original space and the two first KL-PLS components. Colors symbolize conditional probability P(y = 1/X = x).

Fig. 3 .

 3 Fig. 3. Test error as a function of the number of KL-PLS components retained for Banana data. Boxplots show the distribution of the 100 individual test error rates. The black points represent the means and the horizontal red lines the median.

Fig. 4 .

 4 Fig. 4. Comparison between KL-PLS (m = 10 and sigma = 0.9) and SVM (C = 1 and sigma = 0.9) test set accuracy rate as a function of percentage of degraded y.

Fig. 5 .

 5 Fig. 5. Test error rate (respectively, CPU time) of KL-PLS as the function of the number of retained components for ovarian data and lung data (left) (respectively, (right)). Horizontal lines show the CPU time of SMV light . Same legend for the boxplots as in Fig. 3.

Fig. 6 .

 6 Fig. 6. Roc curve of KL-PLS. M i : Accuracy of dermatologist i in terms of sensitivity and (1-specificity) for the melanoma diagnosis. E i : Accuracy of dermatologist i in terms of sensitivity and (1-specificity) for the excision decision.

Table 1

 1 Basic information about the benchmark data sets used in the study

	Data set	Dimension	Training sample size	Test sample size
	Banana	2	400	4600
	B. cancer	9	200	77
	Diabetes	8	468	300
	German	20	700	300
	Heart	13	170	100
	Ringnorm	20	400	7000
	F. solar	9	666	400
	Thyroid	5	140	75
	Titanic	3	150	2051
	Twonorm	20	400	7000
	Waveform	21	400	4600

Table 2 Comparison

 2

of classification errors (mean ± standard deviation) (test sets) for SVM

According to[START_REF] Bastien | PLS generalized linear regression[END_REF] the computation of the PLS-LR component t h may be simplified by setting non-significant regression coefficients a hj to 0 according to the Wald test. Consequently, only variables that are significantly related to y, contribute to the computation of t h . So, the number m of KL-PLS components to be retained may be chosen by observing that the component t m+1 is not significant because none of the coefficients a m+1,j is significantly different from 0.

Acknowledgments

The Ph.D. fellowship of Arthur Tenenhaus is granted by KXEN (Knowledge eXtraction ENgines), a global analytic software company and the National Association of Technical Research. We thank Fanny Aziza, Michel Tenenhaus, Léon Bottou and Yann Lecun for helpful and fruitful discussions.


Farrusseng

Cécile Daniel

Yves Schuurman

David Farrusseng

Discovery of very active catalysts for methanol carboxylation into DMC by screening of a large and diverse catalyst library

Keywords: ceria, zirconia, flame spray pyrolysis, DMC, CO 2

published or not. The documents may come

Introduction

Among organic carbonates, dimethyl carbonate (DMC) is a promising green chemical because of its low toxicity and rapid biodegradability 1 2 • Its versatile chemical properties account for its wide use as an agent in carbonylation and methylation reactions, replacing DMSO, COCl2 and CH 3 X 34• It is also largely used as a precursor for the production of polycarbonates and polyurethanes. DMC offers excellent solving properties and is therefore used as a safe solvent in organic synthesis 5• It has a high oxygen content (53%) and is a promising oxygenate additive to gasoline, compared to other chemicals such as MTBE, ethanol or methanol 6. Due to its high dielectric constant, DMC is used as an electrolyte in lithium-ion batteries 7• DMC was historically synthesized by the harmful methanolysis of phosgene. Current industrial processes include the hazardous, corrosive oxidative carbonylation in gas or in liquid phase 8• Other possible routes are transesterification of propylene carbonate and methanolysis of urea 9 1 0 •

The direct synthesis of DMC from methanol and C02 (reaction 1) appears particularly attractive: it is a greener alternative to the other routes that simultaneously allows C02 valorization 1 0 • However, the methanol conversion, XMeoH, is strongly limited by thermodynamic equilibrium. The standard Gibbs enthalpy was evaluated at 26 kJ.mo1-1 11 • 1 2 • As a consequence, the DMC yield, assuming 100% selectivity, can hardly exceed 2% even at favorably high pressure 13 and at temperatures below 473K.

Dehydration of reactive media can be used to increase the DMC yield. Several researchers have used in situ reactive water traps such as trimethyl phosphate, butylene oxide 14, acetonitrile 15 and 2cyanopyridine 16, reaching per-pass yields of up to 6%. The use of an inorganic water trap such as zeolite 3A in a multi-pass process yields up to 40% DMC 17 . Recent studies demonstrate that a direct synthesis process can technically achieve high DMC yield, when an efficient dehydration process is applied. With this in mind, the application of catalytic membrane reactors allowing continuous water withdrawal during reaction has also been proposed, and a successful proof of concept was recently reported 18 19 .

A large diversity of catalysts have already been tested for the direct synthesis of DMC 20 21 , including organic metal-alkoxy compounds, metal oxides, metal-supported catalysts and ionic liquids. It is generally proposed that well-balanced acidic and basic properties play a key role in DMC synthesis for methanol activation 22 2324 . The majority of metal oxide studies deal with ceria or zirconia 25 26 27 . The effects of crystal structure and morphology on activity are still under debate. Ceria catalysts exhibiting nanorods, nanocubes, octahedrons and spindle-like morphologies show different activities 28 29 30 .

Generally, the observed tendency is that mixed oxides often perform better than pure oxides. Hence, authors have explored modified ceria and zirconia such as H 3 PO 4 -functionalized ZrO 2 23 31 , ceria doped with Al 2 O 3 , Fe 2 O 3 32 33 , as well as mixed oxides of ceria -zirconia 24 34 22 35 36 . Although many catalyst formulations have been evaluated in the literature, one can hardly identify which formulations are the most active. One of the reasons is that because no standard testing protocol has been established, the processes developed have been very diverse. Several studies have been carried out in the gas phase in a fixed-bed reactor at low pressure, enabling continuous withdrawal of water, and therefore conversion beyond thermodynamic equilibrium. Most of the measurements are carried out under pressure in the liquid phase in batch reactors, resulting in very low conversion (<2%) limited by thermodynamic equilibrium. Remarkable differences in yields for both processes are reviewed in 37 .

Unfortunately, most of the studies do not report catalytic activities, i.e., moles converted per unit of time and catalyst mass, but only conversion data. Moreover, data measured at conditions close to thermodynamic equilibrium cannot be used to calculate the catalytic activity. In addition, the equilibrium conversion is difficult to calculate due to the lack of accurate thermodynamic data under reaction conditions that are often non-ideal and the fact that traces of water, notably those contained in the solid catalyst, have a strong impact on the equilibrium. Unfortunately, this last issue has been neglected in most studies 38 39 . Finally, water formed during the reaction may cause the deactivation of some catalysts 40 . As a conclusion, a ranking of catalytic activities across the literature can hardly be established on a rational basis. This is true, for example, for pure cerium oxides, one of the most studied catalysts 28 29 30 41 .

This study aims to rank catalyst activities and to identify very highly active catalysts using a systematic process allowing the measurements of "true" catalytic activities. Based on an assessment of critical protocol parameters, a protocol for the measurement of catalytic activity far from equilibrium was developed, with special attention paid to critical experimental parameters. Following this protocol, a screening of different pure and mixed oxides was carried out. Finally, an optimization study was carried out on Zr-Ce mixed oxides obtained from different synthesis methods, leading to the discovery of the catalyst which, to the best of our knowledge, is the most active reported so far.

Experimental

Catalyst sourcing and preparation

In total 61 catalysts were tested. The list of catalyst formulation and catalytic activities is available in Supporting Information.

Pure oxides were purchased from Sigma Aldrich, Fluka, Acros and Alfa-Aesar. Other samples were synthesized using three different synthesis procedures described below.

Co-precipitated mixed oxides referenced as CP-were prepared from nitrate salt precursors.

Predetermined amounts of salts were mixed in aqueous solutions and added dropwise to a solution of NaOH heated at 333K. The precipitate was washed until conductivity dropped below 50 µS. The final sample was calcined under air at 773K for 5 h. Pure zirconia and pure ceria used to prepare impregnated samples, referenced as P-ZrO 2 and P-CeO 2 , respectively, were synthesized according to the same single precipitation method. Mixed co-precipitated oxides were labeled A x B 1-x O 2 , where A= Ce or Zr, B=Al, Gd, Nb, Ti, Fe, La or Zr, and x stands for the molar fraction of A.

Pure zirconia, pure ceria and zirconia-ceria mixed oxides were also prepared by flame spray pyrolysis, referenced as FSP. The typical equipment and synthesis parameters are described elsewhere 42 43 . The powders obtained were calcined under air at 773K for 5 h. The molar composition is noted

A x B 1-x O 2 .
Finally, a third series of catalysts was prepared by impregnation of ceria P-CeO 2 with a second metal M in solution and calcined at 773K under air for 5 h. They were labeled as z % M/CeO 2 , where M is the impregnated metal and z the targeted weight percentage.

Catalyst Characterization

The surface area was calculated by the Brunauer-Emmett-Teller (BET) method from N 2 physisorption isotherms using a Belsorp Mini instrument (Bel Japan). Before N 2 physisorption, the samples were heated at 473K for 3 h in a vacuum ramp (Belprep, Bel Japan).

Powder X-ray diffraction patterns (XRD) were recorded to assess solid crystallinity. Diffractograms were collected between 4 and 90° + L. with steps of 0.02° and 1 s per step with a Bruker D8 Advance A25 diffractometer using Cu =M radiation at N = 1.5418 Å. The average crystallite size was calculated using the Scherrer method.

Elemental analysis of the catalysts was performed using an ICP-OES (Inductively Coupled Plasma-Optical Emission Spectrometry) ACTIVA from HORIBA Jobin Yvon equipped with a CCD detector for the determination of metal loadings.

Method of carboxylation of methanol and CO 2

Anhydrous methanol was purchased from Carlos Erba with a water content below 300 ppm. Anhydrous toluene was purchased from Sigma-Aldrich with a water content below 10 ppm. Liquid CO 2 was supplied by Messer with a purity of 99.995% and a water content below 5 ppm.

The three-phase catalytic carboxylation of methanol was carried out in a high-pressure Parr Instrument vessel of 50 ml volume equipped with a pressure indicator, used as a batch reactor. Typically, an adjusted amount of dried catalyst in powder form (ranging from 20 mg to 500 mg) and 12 g of methanol with 0.1 wt % anhydrous toluene as internal standard were introduced into the vessel. The reactor was purged 3 times with 10 bar of dry N 2 . A volume of 18 ml of liquid CO 2 kept at 281K and 50 bar was fed to the autoclave with an Isco 250D syringe pump. The reactor was heated in an oil bath at 413 K for 4 h, under magnetic stirring. We observed that the rate at which CO 2 was introduced had an effect on the final pressure. The dissolution of CO 2 into methanol is rather slow as it takes several minutes. As a result, the final amount of dissolved CO 2 depends on the feeding rate. Consequently, in order to obtain a reproducible final pressure, the same CO 2 feeding rate and feed duration were applied by means of a syringe pump. In our protocol the vessel reached a stable pressure of 120 bar after 1 hour of reaction (Table 1S).

Under such conditions, the reaction takes place in the liquid phase with the suspended catalyst, which is in equilibrium with the vapor phase. The reaction was stopped after 4 hours by cooling the vessel in an ice bath. An aliquot from the liquid phase was sampled through a 0.45 µm filter (Millipore) in order to remove powder catalyst. We have analyzed by GCMS the composition of the liquid phase after reaction (Shimadzu QP 2010 GCMS equipped with a ZB Wax column). The only product of the liquid phase is DMC. In addition no DME could be noticeable in the gas phase after reaction carried out with P-CeO 2 and P-ZrO 2 , the later being acid catalyst. As consequence, we do not expect production of DME with other oxides screened in the test conditions applied here and gas phase was not systematically analyzed.

Reaction products were quantified by an Agilent GC 6850 equipped with a capillary column DB WAXetr (30 m) and a flame ionization detector. The GC detection limit of DMC was evaluated at 0.015 wt %.

Each value refering to conversion is evaluated on methanol basis.

The catalyst activity was calculated on the basis of the amount of DMC produced after 4 hours of reaction with an initial molar ratio MeOH/CO 2 =1 at T=413 ± 2 K, and P=120 ± 2 bar and expressed as the number of moles DMC produced per gram of catalyst per second. In contrast to most studies in which activities are calculated with methanol-based conversions approaching thermodynamic equilibrium, the amount of catalyst in this study was adjusted to always keep the conversion below 0.6%. For example, the catalyst loading was 500 mg and 20 mg for the least and most active samples, respectively, while all other testing parameters were kept constant.

The water content was checked by coulometric Karl Fischer titration with a Titrando 905 Methrom apparatus. Hydranal Coulomat (Sigma Aldrich) was used for titration.

3 Results and discussion

Determination of equilibrium conversion and repeatability

The presence of water in the reaction mixtures interferes with the catalytic measurements in different ways. First of all, the reaction produces water and the equilibrium conversion under typical reaction conditions is very low, i.e., 0.7% on a dry basis (see details in E.S.I.). When measuring the catalytic activity, the conversion needs to be sufficiently below the equilibrium conversion to neglect the effect of the reverse reaction on the catalytic activity. An accurate estimation of the equilibrium conversion is thus necessary and requires knowledge of the initial water content. Water can also inhibit the reaction rate by chemisorption on the active sites in competition with methanol chemisorption on the same sites. A variation of the initial water content might lead to a variation of the catalytic activity, so it is important that all samples are evaluated under the same conditions. Further detailed kinetic studies will be necessary to reveal this inhibition effect. Finally, water might lead to catalyst deactivation. Therefore, the different origins and quantities of the introduction of water into the reaction mixture were evaluated.

Assuming the presence of 300 ppm of water in methanol and 5 ppm in CO 2 , the initial water content represents less than 10% (0.002 g) of the expected water production by reaction assuming a conversion of 0.7% on a dry basis (0.0236 g) (details in SI). Another possible source of water contamination can be physisorbed water on the catalyst if it is not properly dried prior to catalytic measurements; this may represent up to 30% of the expected water production (see E.S.I. Table 2S). In light of these preliminary estimations, all catalysts were desorbed at 473K under primary vacuum in our protocol. A fresh bottle of methanol was used at each catalytic run. The water concentration was always below 300 ppm (Karl Fischer measurements).

As already stated, because of the lack of accurate thermodynamic data for the given operating conditions and the fact that water content has a strong impact on the equilibrium conversion, the conversion at thermodynamic equilibrium was experimentally measured. Conversion measurements were carried out with four different catalysts at 413K and 120 bar for an initial MeOH:CO 2 ratio of 1 and an estimated initial water content of 300 ppm. Conversion was measured after 24 hours of reaction, which is six times longer than the time period used for activity measurements. The conversion at thermodynamic equilibrium (X eq ) was found to be near 0.8% ±0.02 (sd) for the four catalysts and is therefore catalyst independent as expected (Fig. 1). The value of 0.8% was taken as the thermodynamic equilibrium, which is in good agreement with literature data between 413 and 443K 39 26 . The reaction conversion as a function of time was followed over 48h over CP-Ce 0.63 Zr 0.37 O 2 (see E.S.I.

Fig 4S)

. The initial rate of reaction was found to be linear until 4h. According to this result, the reaction time was set at 4 hours for the screening protocol which was used to measure the catalytic activity of the catalyst library. The stirring rate was not found to influence catalytic activity measurements. Increasing the stirring rate from 1500 rpm to 2000 rpm did not change results in the conditions described in Fig. 2.

A repeatability study of the activity measurement protocol was carried out through 10 identical runs over the same catalyst. We established an average conversion of 0.29% ± 0.02 (sd), corresponding to a relative error of ± 7% that has been applied for the calculation of activities.

Library of pure oxides

Commercial oxides were selected to complement literature results covering different acid-base scales according to 44 45 . These commercial oxides were revealed to have low specific surface areas below 10 m 2 g -1 . Among all oxides tested, only ceria and zirconia were active. CeO 2 and ZrO 2 were fully selective in DMC synthesis. No other by-products, such as DME, were found. 46 . Different morphologies of ceria such as nanorods, octahedrons and nanocubes were also used to produce DMC 30 . The authors concluded that the crystal planes as well as the number of acid and basic sites are decisive factors for active DMC catalysts. Hence, the difference in surface properties and crystal defects may explain the difference in activity observed in Table 2 for the pure ceria samples.

Ceria mixed oxides and zirconia mixed oxides

Impregnation method A series of ceria mixed oxides was prepared by impregnation of a second metal followed by a thermal treatment. The initial pristine ceria (noted P-CeO 2 ), was prepared by precipitation and provided by Johnson Matthey. The addition of a second metal (Al, Zn, Fe, La, Y, Gd, Sm, Zr, Nd) was carried out by a wetness impregnation process to reach a loading between 1 -5 wt.%. After heating the sample at 773K in air, the surface area of the serie ranged from 85 to 97 m 2 g -1 were obtained, which was 15 to 25% lower than the surface of the pristine ceria material. Two additional catalysts were prepared by adding a second thermal treatment at 1023K in air on 1% Al/CeO 2 and 5 % La/CeO 2 . Taking the activity of precipitated pristine ceria (P-CeO 2 ) as a benchmark, we can define metal additives that do not change the activity, that improve the activity or that on the contrary are detrimental.

Obviously, ceria impregnated with Fe, Y, Gd, Sm or Zr oxides and treated at 773K did not show any significant effect on the activity with respect to the pristine ceria. We want to stress here that we do not conclude that the addition of one of these metal oxides using a different synthesis procedure may not lead to higher activities. Indeed, as we will show below, Ce-Zr mixed oxides obtained by other means are much more active. Impregnation with Al made the activity drop dramatically, despite a high measured specific surface area of 85 m 2 .g -1 . This might be due to coverage of active ceria sites by a layer of amorphous alumina. Nevertheless, even a five-hour thermal treatment at 1023K for 1 wt.% Al/CeO 2 , which could have led to a better alumina dispersion, did not improve the performance, although its specific surface area remained high (52 m 2 g -1 ).

The observed trend indicates a significantly enhanced activity in the following order: La > Nd = Zn. The increase of the La loading from 1 wt.% to 5 wt.% results, however, in a loss of the activity. On the other hand, the activity of 5 wt.% La/CeO 2 -750 is surprisingly boosted after a thermal treatment at 1023K despite a reduced specific surface area of 63 m 2 g -1 .

Ceria and zirconia mixed oxides prepared by co-precipitation with a second metal oxide All mixed oxides prepared by co-precipitation were nanocrystalline with a crystallite size between 4 and 8 nm (estimated by the Debye-Scherrer method) and specific surface areas (BET method -N 2 )

ranging from 90 to 130 m 2 g -1 (see Table 5S , E.S.I). XRD analysis showed that cerium-based oxide was crystallized in a cubic structure, whereas zirconia-based samples had a tetragonal structure (Fig. 2S and Fig. 3S). The co-precipitation with niobium oxide inhibited the activity, suggesting that large amount of Nb is not favorable to DMC synthesis 47 . A 10 wt.% concentration of titanium oxide decreases the catalytic activity, whereas a low incorporation (2 wt.%) has no significant effect.

Third generation: ceria-zirconia mixed oxides by co-precipitation and flame spray pyrolysis

A number of literature studies have highlighted high activities of mixed ceria-zirconia samples 22 31 34 35 36 . However, the effect of the mixing is not obvious with respect to the pure oxides, as both are active. In addition, we show here that the synthesis of mixed oxides by post-impregnation of zirconium oxide precursor on cerium oxide does not improve the catalytic activity (Fig 3) whereas co-precipitation leads to a synergistic effect (Fig 4). Apparently, the preparation method is a key parameter for obtaining very active ceria-zirconia catalysts. Therefore, we decided to complete the screening of ceriazirconia catalysts prepared by co-precipitation (CP) by investigating the effect of the composition range from pure ceria to pure zirconia and the effect of the heat-treatment temperature from 423 to 1123K.

In addition, we have explored the catalytic activities of ceria-zirconia catalysts prepared by flame spray pyrolysis. In order to compare catalytic results between CP and FSP synthesis methods, we have targeted the same molar composition (noted x hereafter) and the same treatment temperatures were applied.

Co-precipitated mixed oxides

Co-precipitated mixed oxides are true nanocrystalline solid solutions as far it can be determined by XRD (see Fig 2S). The surface areas of the solids obtained ranged from 87 to 126 m 2 g -1 , for crystallite sizes between 3.6 and 8 nm (see E.S.I Table 6S). The activities of CexZr 1 _x02 mixed oxides prepared by co-precipitation and calcined at 773K are presented in Fig. 5. Interestingly, all intermediate compositions were more active than the pure oxides Ce02 and Zr02, i.e., there is a clear synergistic phenomenon. The maximum of activity is reached for the composition Ce0_77Zr0.2302 (x=0.77). These results are in good agreement with previous studies showing that the catalytic activity strongly depends on the Ce/Zr ratio. In 22 the authors prepared CexZr 1 _x02 by a sol-gel method and calcined the catalysts at 773K. The amount of DMC produced showed a "volcano-shaped" curve with respect to the cerium content, with a maximum at x=0.6 corresponding to Ce0_6Zr00. 4 02• Another study found the same composition-activity profile with the highest activity for Ce0_5Zr 0_502 prepared by a citric acid sol-gel method 4 8 • However, for reasons that are unknown to us, the activity of the latter is about one order of magnitude lower.

The co-precipitated Ce0_88Z0. 1 202 was selected for optimization by post-thermal treatment. We observe, as expected, a decrease of the surface area upon increasing the temperature that is due to sintering (Fig. 6). Unexpectedly, the activity follows the opposite trend, i.e., it increases when the

Flame spray pyrolysis ceria-zirconia

In good agreement with previous studies 49 so the characterization of flame spray pyrolysis ceriazirconia samples revealed solid solutions with specific surface areas between 88 and 99 m 2 g-1 and crystallite sizes between 6 and 10 nm (E.S.I Table 75).

The flame spray pyrolysis ceria-zirconia were about one order of magnitude more active than the corresponding co-precipitated solids of the same composition. Whereas the studies of co-precipitated ceria-zirconia were carried out using 0.1 g of solid, the amount of catalyst was adjusted down to 0.02 g for the flame spray pyrolysis samples, otherwise the equilibrium conversion would have been reached well before the sampling time (see section 3.1.3).

The activities of flame spray pyrolysis ceria-zirconia as a function of the composition are compared to the co-precipitated (CP) solids (Fig. 7). The activity profile follows a "volcano" curve as reported above and is similar in trend to the co-precipitated samples. However, we cannot define here which composition corresponds to the highest activity, because the region with 0.5 < x < 1 was not experimentally investigated. Nevertheless, we can assume from the shape of the curve and the trend of the co-precipitated ceria-zirconia oxide that the maximum activity should be found in the region of equimolar Zr:Ce composition (x=0.5).

Astonishingly, the activity of the FSP mixed oxide Ce0_5Zr0_50 2 (x=0.5) reached an activity of 6 µmol g-1 s-1, which is 20 times higher than that of the CP mixed oxide with the closest molar composition (x=0.45). FSP zirconia was six times more active than precipitated zirconia, whereas the difference in activity was negligible for pure ceria. The activity scale is presented as a log scale. FSP: The catalytic data obtained in this study for ceria-zirconia mixed oxides with an equimolar composition or nearly equimolar composition is compared with literature data (Fig. 8). We have restricted our comparison to data obtained under three-phase conditions without a water removal system. It is very difficult to compare catalytic activities under the same experimental conditions, as no standard testing protocol exists. Nevertheless, for each set of published data, the highest conversion reached versus ceria-zirconia composition was selected. Activities were calculated from the parameters of the corresponding study (mass of catalyst and reaction duration). Obviously the FSP catalyst is far more active than any other catalyst reported (6 µmol .s -1 .g -1 ), while the co-precipitated catalyst showed activities in the same order of magnitude as other good ceria-zirconia mixed oxides reported in the literature.

The unexpectedly high catalytic activity of the flame spray pyrolysis ceria-zirconia led us to investigate other ceria-based oxides prepared by flame spray pyrolysis. The second metal was selected among the best found by impregnation or co-precipitation methods, namely Zn, La and Nd. Ceria mixed oxides were obtained using the same process as for flame spray pyrolysis ceria-zirconia. The specific surface areas for these solids were around 95 (±8) m 2 g -1 (E.S.I, Table 8S). The catalytic results are benchmarked against the flame spray pyrolysis ceria (Fig. 9). The addition of La, Nd and Zn significantly increased the activity compared to pure FSP ceria. This result confirmed the value of introducing a small amount of Zr, La or Nd into the ceria as previously reported for impregnated samples. Nevertheless, the activity gap compared to the co-precipitated samples was not as remarkable as it was for ceria-zirconia.

Conclusions

The present work describes a robust method for measuring catalytic activity for DMC synthesis from methanol and CO 2 far from equilibrium conditions. This method made it possible to rank the activity of more than 60 catalysts which, to the best of our knowledge, is the largest consistent dataset yet reported. Different synthesis methods and compositions were tested. By selecting the most active composition of ceria-zirconia samples and further exploring different synthesis approaches, we found that the synthesis of mixed ceria-zirconia oxides by flame spray pyrolysis yielded the most active catalysts reported so far. We do not claim that other catalysts with higher activity could not exist, as the screening was of moderate size. On the contrary, we believe that better catalysts may exist with different Ce-Zr ratios, in particular Ce-rich compositions or solids that have undergone an optimized thermal treatment.

It is common to say that the structure of mixed oxides, especially the surface composition, strongly depends on the synthesis method and on the thermal treatment. We have observed that the addition of alumina yields better catalysts when it was added by co-precipitation, whereas a drop in activity was obtained by impregnation. We have also observed that for the latter method, high catalytic activities is obtained after a post-treatment at high temperature. Surprisingly, the catalytic activity of precipitated Ce 0.88 Zr 0.12 O 2 treated at different temperature is inversely correlated to its surface area.

This suggests major modifications of the surface upon thermal treatment.

The activity of ceria-zirconia mixed oxides is about one order of magnitude higher when they are prepared by flame spray pyrolysis as opposed to by co-precipitation. The characterization of the structure by XRD or the texture by N 2 physisorption did not reveal significant differences that could account for such a difference in catalytic activities. More detailed investigations are required for the elucidation of structural features that are responsible for the high catalytic activity of mixed ceriazirconia obtained by flame spray pyrolysis. This study will be published separately.

Fig. 1 :

 1 Fig. 1: Evaluation of conversion at thermodynamic equilibrium for the set of conditions: 413K, 120 bar, ratioMeOH:CO 2 :1:1, t=24 h, Wcat=0.5 g, batch operation.

Fig. 2 :

 2 Fig. 2 : Repeatability of conversion over CP-Ce 0.63 Zr 0.37 O 2 at 413K for 4 h, molar ratio MeOH/catalyst =577 mol/mol. P=120 bar, MeOH:CO2=1:1,Wcat= 0.1 g.

1 )Fig. 3 :

 13 Fig.3: Activity of ceria mixed oxides obtained by impregnation-heat treatment. Except for those noted with the label -750, all samples were heat-treated at 773K in air. P-CeO 2 is pristine ceria before impregnation. The horizontal dashed line is a guide for the eyes.

Fig. 5 :

 5 Fig. 5: Activity of Ce x Zr 1-x O 2 prepared by co-precipitation and calcined at 773K.

Fig. 6 :

 6 Fig. 6: Evolution of the activity(•) and surface area (o) of Ce0_88Zr0.1202 prepared by co-precipitation as a function of the temperature treatment.

Fig. 7 :

 7 Fig. 7: Comparison of catalytic activities of flame spray pyrolysis CexZr1_x02 compared to co-precipitated solids. The activity scale is presented as a log scale. FSP: •, CP: o

1 )Fig. 9 :

 19 Fig. 9: Activities of flame-sprayed ceria-based oxides. The horizontal dashed line is a guide for the eyes. Dopants content are expressed in wt%.

Table 1 :

 1 Activity of pure oxides for direct synthesis of DMCVarious ceria sources with different specific surface areas and particle size distributions were evaluated from commercial suppliers. Synthesis methods of commercial oxides are unknown. All ceria samples were active in DMC synthesis and fully selective.

			Activity @413K	S BET
	Oxide	Properties
			µmol DMC.g -1 . s -1	m 2 .g -1
	TiO 2	Amphoteric	0	7.5.
	La 2 O 3	Lanthanoids, Basic	0	5.2
	Pr 6 O 11	Lanthanoids, Basic	0	4
	Nb 2 O 5	Acidic	0	0.5
	SrO	Basic	0	2.6
	Gd 2 O 3	Lanthanoids, Basic	0	2
	Ga 2 O 3	Amphoteric	0	11
	ZnO	Amphoteric	0	7.5

Table 2 :

 2 Activity and specific surface areas of various ceria.

	Entry		Specific surface
		Ceria sample	area	Activity
			(m 2. g -1 )	(µmol.g -1 .s -1 )
	A	Ceria 1 (Johnson Matthey)	100	0.51
	B	Ceria 2 (Johnson Matthey)	135	0.25
	C	Ceria 3 flame spray pyrolysis (Johnson Matthey)	99	0.31
	D	Ceria precipitated	115	0.38
	E	Ceria HSA (Solvay)	237	0.2
	F	CeO 2 Aldrich nanopowder (particle size < 25 nm)	30	0.16
	G	CeO 2 Aldrich microsized (particle size < 5 µm)	20	0.03
	H	CeO 2 fused (3-6 mm pellets)	2	0.01

Although only ceria is concerned, no clear relation can be established between the activity data and the specific surface area (see Figure

1S

in E.S.I). On one hand, low specific surface area samples below 50 m 2 .g -1 generally exhibited poor activity. But on the other hand, high specific surface area alone is not sufficient for high activity, as demonstrated by the sample with a surface area of 237 m 2 .g -1 (entry E), but an activity of only 0.2 µmol.g -1 .s -1 . Yoshida et al. have concluded that high crystallinity of ceria calcined below 873K is a favorable property for catalytic activity

  For zirconia-based catalysts obtained by co-precipitation, namely Zr 0.98 Ti 0.02 O 2 , Zr 0.91 Ti 0.09 O 2 , Zr 0.99 Al 0.01 O 2 , Zr 0.97 Al 0.03 O 2 , we can hardly identify which metal oxides have a positive or negative effect due to the very low activity. We can observe that the Al-or Fe-containing oxides increased the activity of about 20%, in agreement with Aresta et al. who have reported higher activity for mixed oxides, namely 3-10 wt.% Al/CeO 2 32 33 and 7 wt.% Fe/CeO 2 prepared by co-precipitation.

			0.6
	Activity	(µmol.g -1 .s -1	0.2 0.4
			0.0

)

Fig. 4: Catalytic activity of co-precipitated Ce-and Zr-mixed oxides (composition in molar composition). The horizontal dashed line that corresponds to the activity of P-CeO 2 is a guide for the eyes.

At a first glance, ceria-based catalysts prepared by co-precipitation are more active than zirconia-based samples.

  The X-ray diffraction patterns of co-precipitated and flame spray pyrolysis ceria-zirconia are compared in the SI(Fig 2S and 3S). For each set of catalysts, the continuous shift of the patterns versus the molar composition indicates that both CP and FSP series are solid solutions of ceria-zirconia. Although the XRD is not accurate enough to draw a conclusion regarding the homogeneity of the solid solution at nanoscale level, we cannot see a significant structural difference in the bulk phase51 . Hence, the crystallite sizes of CP and FSP (Debye-Scherrer) are in the same range (E.S.I Table6S, 7S). The specific surface area of whole CP and FSP mixed oxides vary relatively little (from 93 m 2 g -1 up to 126 m 2 g -1 ). In contrast, the catalytic activity varies by more than an order of magnitude across the CP and FSP mixed oxides. Such a difference in activity cannot be the result of a different number of the same active sites and should originate from a different kind of active sites. A more detailed investigation of the crystalline structure and the surface properties by complementary techniques is required to postulate a possible relation between structural features and the catalytic activity. This study goes well beyond the scope of this paper and will be published elsewhere. Comparison of catalytic activities in log scale of ceria-zirconia mixed oxides with an equimolar composition or nearly equimolar composition from this study and compared with literature data A: 48 ; B:24 ; C:

	Activity
	µmol.g -1 .s -1
	10	140°C	413K
	1
	0.1		140°C 413K	110°C 493K	170°C 443K	170°C E 443K
			100°C 373K		C	D 393K 120°C
	0.01
	Fig. 8:

•, CP: o 22 ; D:

34

; E:

35

. For each condition, the temperature is mentioned above the bar.

Acknowledgements

The work leading to these results received funding from the European Union Seventh Framework Program FP7-NMP-2010-Large-4, under grant agreement no 263007 (acronym CARENA). We thank Solvay and Johnson-Matthey for providing samples.