Sylvain Takerkart
email: sylvain.takerkart[@]incm.cnrs-mrs.fr

Liva Ralaivola
email: liva.ralaivola[@]lif.univ-mrs.fr

MKPM: a Multiclass extension to the Kernel Projection Machine

We introduce Multiclass Kernel Projection Machines (MKPM), a new formalism that extends the Kernel Projection Machine framework to the multiclass case. Our formulation is based on the use of output codes and it implements a co-regularization scheme by simultaneously constraining the projection dimensions associated with the individual predictors that constitute the global classifier. In order to solve the optimization problem posed by our formulation, we propose an efficient dynamic programming approach. Numerical simulations conducted on a few pattern recognition problems illustrate the soundness of our approach.

Introduction

Many real world problems of pattern recognition pose the question of having at hand efficient and effective methods to tackle multiclass learning problems. In this paper, we propose a very simple yet effective new mutliclass kernel machine, meaning that the predictors f we are to consider are functions like

f (x) = f φ(x), (1)
defined on input space X , where φ : x → φ(x) = k(•, x) is the mapping from X to H associated with some positive definite kernel k : X × X → R [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]; here f φ(x) denotes the inner product between f ∈ H and φ(x) ∈ H.

Our approach rests on the idea of kernel methods regularized by finite-dimensional projections, methods that we broadly refer to Kernel Projection Machines1 (KPM) [START_REF] Blanchard | Finite-dimensional projection for classification and statistical learning[END_REF][START_REF] Vito | Learning from examples as an inverse problem[END_REF][START_REF] Zwald | Kernel projection machine: a new tool for pattern recognition[END_REF]. Given a training sample Z = {(x n , y n )} N n=1 where the data x n 's are from some space X and the targets y n are in Y, the strategy implemented by these machines is the following. First, compute a family Ψ = (ψ d ) d≥1 of functions in H (by, e.g., performing a Kernel principal component analysis of the x n 's) and its associated nested subspaces S 0 ⊂ S 1 ⊂ . . . ⊂ S d ⊂ . . ., where S d denotes the space spanned by the first d functions ψ 1 , . . . , ψ d -with S 0 = {0}. Then, for each d, solve the unregularized empirical minimization problem fd = argmin

f ∈S d 1 N N n=1 (f, x n , y n ), (2)
where : H × X × Y → R is a loss function, such that (f, x, y) accounts for the fact that f does not predicts the correct target y for x. Finally, select one of the fd 's according to some criterion: either based on theoretical arguments or based on an estimate of the generalization error computed on a validation set. The rationale that supports the use of unregularized empirical risk minimizer is that the process of considering nested subset of finite dimension precisely plays the role of regularization.

The approach that we propose for multiclass classification with KPMs builds on the original idea of KPM and the idea of output coding for multiclass classification [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF].

The formulation that we propose amounts to a minimization problem where several KPMs have to be learned in parallel and simultaneously regularized. This specifically means that 'optimal' projection dimensions should be found for each of the KPM learned and that the total number of dimensions retained has to be controlled, which gives rise to a combinatorial problem. In order to solve this problem, we propose a dynamic programming approach that allows us to come up with a very efficient learning procedure.

To summarize, our contribution is threefold: (i) a new model of multiclass prediction based on the idea of KPMs where the regularization parameter is a global dimension parameter, (ii) an efficient dynamic programming strategy to solve the combinatorial problem arising from our framework and (iii) numerical simulations conducted on pattern recognition problems that support the soundness of our approach. We may notice that, besides its conceptual simplicity, our approach makes it possible to compute the whole regularization path for free, which renders the validation procedure of the global regularization parameter quite inexpensive.

The paper is organized as follows. Section 2 introduces the notation that is used throughout the paper, briefly recalls the approach of output coding for multiclass problem and presents our new Multiclass Kernel Projection Machine (MKPM). Section 3 discusses several interesting features of our learning procedure, such as large-scale learning, model selection and generalization properties. In Section 4, results from numerical simulations conducted on benchmark datasets are reported as well as very positive results for a task of classifying images from a functional magnetic resonance imaging (fMRI) experiment; all the numerical results support the soundness of our proposed model.

Multiclass Kernel Projection Machines

From here on, the following notational conventions hold. For a positive integer N , [N ] denotes the set {1, . . . , N } and [N ] 0 the set {0, . . . , N }. For a vector x, x denotes its transpose, the same notation being used for matrices.

The input space X where the data to classify live is provided with a positive kernel function k : X × X → R (see, e.g., [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF] for a detailed account on kernel functions/kernel machines). The reproducing kernel Hilbert space H associated with k is such that, ∀f ∈ H (note that H is a functional space):

f (x) = f, k(•, x) H = f φ(x), ∀x ∈ X , (3)
where •, • H is the inner product of H and the notation f g is a shorthand for f, g H . The mapping φ :

X → H is such that φ(x) = k(•, x).
The set of labels is denoted by Y, and, to simplify notation, we further suppose that for a categorization problem with Q classes, the set of labels is Y = {1, . . . , Q}.

A training set Z = {(x n , y n )} n∈[N ] is made of N labeled pairs (x n , y n ) such that y n ∈ Y is the class of x n ∈ X ;
as commonly assumed in statistical learning, these pairs are independent realizations of a random variable (X, Y ) distributed according to an unknown and fixed distribution D on Z := X × Y.

Given Z, Φ denotes the (possibly infinite dimensional) matrix

Φ = [φ(x 1 ) • • • φ(x N )] (4)
such that the Gram matrix K = (k(x i , x j )) i,j∈[N ] of the data with respect to k may be written as

K = Φ Φ. (5)
Finally we will consider loss functions : H × X × Y → R + that are convex in their first argument. An example of such a loss function is the square loss square defined as:

square (f, x, y) = f φ(x) -y 2 = (f (x) -y) 2 , (6)
which we will subsequently make extensive use of.

Multiclass Learning with Output Codes

Multiclass prediction is a machine learning problem on it own, which raises questions related to modelling, as well as statistical and algorithmic issues. A very powerful and effective way to address this problem is to make use of output codes, which is thoroughly described in [START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF]. Output coding consists in associating a unique binary code c q ∈ {-1, 1} L of length L to each category q of Y. Given these codes, the traning set can be rewritten as

Z = {(x n , y n )} n∈[N ] , where y n ∈ {-1, 1} L is one of the Q codewords c q .
The initial learning problem is then converted to L learning problems defined with respect to the sequence of training sets

Z l = {(x n , y nl )} n∈[N ] , l ∈ [L], (7)
where, with a slight abuse of notation, we have reused the notation y and where y nl ∈ {-1, 1} is the l-th entry of y n . Given these L training sets, a usual strategy is to independently learn L predictors f 1 , . . . , f L such that each f l : X → R is inferred from Z l . The prediction ŷ(x) for a new instance x is then computed from the vector

f (x) = [f 1 (x) • • • f L (x)] as ŷ(x) = argmin q∈[Q] dist(f (x), c q ), (8)
where dist : R L × R L → R is some predefined divergence measure; for instance, this divergence measure may be the usual Euclidean distance, which we make use of later on.

The approach we propose builds on the framework of output codes and differs from the usual learning strategy in that the L classifiers inferred are simultaneously regularized or, in other words, co-regularized. This essentially means that their learning are no longer independent.

Proposed Learning Model

As briefly recalled in the introduction, Kernel Projection Machines work in the following way. Given an integer D and a family Ψ = (ψ d ) d≥1 of linearly independent functions of H, defining S d = span(ψ 1 , . . . , ψ d ) as the space spanned by the first d vectors of Ψ (with S 0 = {0}), KPMs implement a two-step strategy:

1. For each dimension d ∈ [D] 0 , solve the empirical risk minimization problem fd = argmin f ∈S d 1 N n∈[N ] (f, x n , y n ). (9)
2. Pick the optimal dimension d and corresponding classifier f d according to some model selection criterion.

As stated above, is a convex loss function. For classification as well as regression, it has been showed that constraining the search of a function f in a finite-dimensional subspace S d incurs some kind of regularization and that the pivotal parameter of regularization is precisely d [START_REF] Blanchard | Finite-dimensional projection for classification and statistical learning[END_REF][START_REF] Vito | Learning from examples as an inverse problem[END_REF]. This explains why there is no explicit regularization term in the empirical risk functional of ( 9).

The approach that we propose pushes the idea of KPM learning one step further. Namely, if we assume the same family Ψ as before and an output coding scheme with code length L for the multiclass learning problem. The twostep learning process based on the training sample

Z = {(x n , y n )} n∈[N ] is the following. 1. For each d ∈ [D] 0 , solve the empirical minimization problem δ(d) = argmin d 1 ,...,d L ∈[d] 0 l∈[L] d l =d R(d 1 , . . . , d L ) (10)
where the empirical risk R(d 1 , . . . , d L ) is defined as We now elaborate more precisely on how the first step of the procedure is performed: this comes down to studying how Equations ( 10) and ( 11) interplay with each other. As for the second step, we will assume that estimates of the generelization performances achieved for models associated with different values of d are computed on an held-out validation set or via a cross-validation procedure.

R(d 1 , . . . , d L ) := min f l ∈S d l 1 N n∈[N ] l∈[L] (f l , x n , y nl ). (11
Viewing (11) as a sum of precomputed terms. First note that the objective function of minimization problem [START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fmri data[END_REF] that defines R may be rewritten as:

1 N n∈[N ] l∈[L] (f l , x n , y nl ) = l∈[L] G l (f l ), (12)
where

G l (f ) := 1 N n∈[N ] (f, x n , y nl ). (13)
Hence, according to [START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fmri data[END_REF] and ( 12), the convex minimization problem that defines R has an objective that is the sum of L convex functions (the G l 's) defined on independent variables (the f l 's). In order to evaluate R(d 1 , . . . , d L ) for some d 1 , . . . , d L , it therefore suffices to independently solve the L convex minimization problems

G * l (d l ) := min f ∈S d l G l (f ), l ∈ [L], (14)
where the G * l (d l ) are uniquely defined since the G l 's are convex. In other words,

R(d 1 , . . . , d L ) = l∈[L] min f ∈S d l G l (f, Z l ) (15) = 1 L l∈[L] G * l (d l ), (16)
and, if all the G * l (d) for l ∈ [L] and d ∈ [D] are computed beforehand, i.e., each problem ( 14) is solved for all values of l and d l , R(d 1 , . . . , d L ) is merely a sum of L precomputed terms.

Problem [START_REF] Hussain | Theory of matching pursuit[END_REF] and dynamic programming. This view of R as a sum of L terms, makes it easier to observe that the problem posed by our learning strategy is a combinatorial minimization problem. Indeed, problem [START_REF] Hussain | Theory of matching pursuit[END_REF] rewrites as:

δ(d) = argmin d 1 ,...,d L ∈[d] 0 l∈[λ] d l =d l∈[L] G * l (d l ). (17)
A direct enumeration procedure to solve this problem would require to scan all the combinations d 1 , . . . , d L ∈ [d] 0 of L integers such that l d l = d, which is obviously prohibitive. Henceforth, a more clever way than enumeration should be devised to solve [START_REF] Williams | Using the Nystrm Method to Speed Up Kernel Machines[END_REF] and it turns out a dynamic programming [START_REF] Sniedovich | Dynamic Programming Foundations and Principles[END_REF] approach may be implemented to this end. Indeed, let us consider the following minimization problem

∆ λ (d) := min d 1 ,...,d λ ∈[d] 0 l∈[λ] d l =d l∈[λ] G * l (d l ). ( 18
)
When λ is set to L, then we observe that ∆ L (d) is directly linked to δ(d) of ( 10) and ( 17) by

∆ L (d) = R( δ(d)). (19)
Therefore, being capable of computing ∆ L (d) is equivalent to being able to solve [START_REF] Williams | Using the Nystrm Method to Speed Up Kernel Machines[END_REF], and we therefore focus on dealing with [START_REF] Zwald | Kernel projection machine: a new tool for pattern recognition[END_REF]. We have the following proposition, which gives us the dynamic programming scheme, i.e., a recursive formulation, to efficiently solve [START_REF] Zwald | Kernel projection machine: a new tool for pattern recognition[END_REF].

Proposition 1. If we define the series ∆λ (d)

λ∈[L] d∈[D] 0 such that ∆λ (d) :=    G * 1 (d) if λ = 1, min δ∈[d]0 ∆λ-1 (δ) + G * λ (d -δ) otherwise. (20) then, ∀λ ∈ [L], ∀d ∈ [D] 0 , ∆λ (d) = ∆ λ (d). (21)
Proof. Let us fix d ∈ [D] 0 and proceed by induction on λ.

For λ = 1, Equation ( 18) gives

∆ 1 (d) = min d 1 ∈[d] 0 l∈[1] d l =d l∈[1] G * l (d l ) = G * 1 (d). (22)
And, by definition of ( ∆λ (d)) in ( 20)

∆1 (d) = G * 1 (d). (23)
Hence for λ = 1 assertion (21) holds. Assume that ( 21) is valid for λ -1, that is

∀d ∈ [D] 0 , ∆λ-1 (d) = ∆ λ-1 (d). (24)
By construction, it is straightforward to see that ∆λ (d) writes as

∆λ (d) = l∈[λ] G * l ( dl ), (25)
for some d1 , . . . , dλ such that l∈[λ] dl = d. Hence, according to the definition of ∆ λ (d) (see [START_REF] Zwald | Kernel projection machine: a new tool for pattern recognition[END_REF]):

∆λ (d) ≥ ∆ λ (d). (26)
Conversely, for all d 1 , . . . ,

d λ ∈ [d] 0 such that l∈[λ] d l = d, l∈[λ] G * l (d l ) = l∈[λ-1] G * l (d l ) + G * λ (d λ ) (27)
≥ ∆ λ-1 (d -d λ ) + G * λ (d λ ) (28) = ∆λ-1 (d -d λ ) + G * λ (d λ ) (29)
≥ min

d λ ∈[d]0 ∆λ-1 (d -d λ ) + G * λ (d λ ) (30) = ∆λ (d), (31)
where (28) comes from ( 18) and (29) comes from induction assumption (24). Taking the minimum of the left-hand side of (27) with respect to Algorithm 1 summarizes the learning algorithm of MKPMs in the case of a generic convex loss function and a predefined family Ψ = (ψ d ) d≥1 of linearly independent functions from H.

d 1 , . . . , d λ such that l∈[λ] d l = d gives ∆ λ (d) = min d 1 ,...,d Λ ∈[d] 0 l∈[λ] d l =d l∈[λ] G * l (d l ) ≥ ∆λ (d), (32)

Algorithm 1 MKPM with output codes of length L

Require: Z = {(x n , y n )} n∈[N ] (training set), Z = {(x n , y n )} n∈[N ] (validation set), Ψ = (ψ d ) d≥1 , D ≥ 1 Ensure: d, d1 , . . . , dL ∈ [D] 0 , f1 , . . . , fL ∈ H for l = 1 to L do for d = 0 to D do G * l (d) = min f ∈S d 1 N N n=1 (f φ(x n ), y nl ) end for end for for d = 0 to D do ∆ 1 (d) = G * 1 (d) end for for λ = 2 to L do ∆ λ (d) = min δ∈[d]0 {∆ λ-1 (δ) + G * λ (d -δ)} end for
Choose and return d (and the corresponding d1 , . . . , dL ∈ [D] 0 , and f1 , . . . , fL ∈ H) according to an estimate of the generalization error computed on Z .

Least-Squares MKPM and KPCA

So far, we have described our new model with an abstract convex loss function and we have assumed that the family Ψ of vectors from H was predefined. In this section, we detail the computations that need be implemented for two particular choices for the setting presented in the previous section. Namely, we make use of the square loss function ( 6) and the family of vectors that we use corresponds to the principal components of the vectors φ(x 1 ), . . . φ(x N ).

Computing the family Ψ with Kernel PCA. Here, we just briefly recall how to perform a kernel PCA. Further details can be found in [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. The computation of the principal components ψ 1 , . . . , ψ D starts by observing that they are defined as the eigenvectors of the (empirical) covariance matrix C (we do not consider the problem of centering the data):

C = 1 N ΦΦ . (33)
A quick analysis shows that in order to compute these vectors, it suffices to extract the eigenvectors A = [a 1 • • • a D ] of the Gram matrix K, with corresponding strictly positive eigenvalues γ 1 , . . . γ D , to have

ψ d := 1 √ γ d Φa d , (34)
or, in matrix form

Ψ := ΦA √ Γ -1 , ( 35
)
where Γ is the diagonal matrix with diagonal elements γ 1 , . . . γ D .

Combining the square loss and the KPCA basis. We now recall how Kernel PCA and least-square regression gently combine to evaluate the values of the coefficients G * l (d). Suppose that we are interested in the lth subproblem Z l = {(x n , y nl )} n∈[N ] (among the L possible problems). Searching for f ∈ S d := span(ψ 1 , . . . , ψ d ) that solves the least square problem

min f ∈S d G l (f ) := N n=1 f φ(x n ) -y nl 2 , ( 36
)
may be achieved by parameterizing

f as f = Ψ [d] α, where α ∈ R d and Ψ [d]
is the matrix made of the first d columns of Ψ (cf. ( 35)), and by noting that G l can be written in matrix form as:

G l (f ) = f ΦΦ f -2f Φy 1 + y l 2 . ( 37
)
where

y l = [y 1l • • • y nl ] . Thus, G l (Ψ [d] α) = α Ψ [d] ΦΦ Ψ [d] α -2α Ψ [d] Φy l + y l 2 . (38)
According to the definition of Ψ (and thus

Ψ [d] ), Ψ [d] Φ = Γ [d] A [d] , (39)
where

Γ [d] is the upper left d × d block of Γ and A [d]
is the matrix made of the first d columns of A. This gives

G l (Ψ [d] α) = α Γ [d] α-2α Γ [d] A [d] y l + y l 2 . ( 40
)
Computing the gradient of G l (Ψ [d] α) with respect to α and having it be equal to 0 leads to:

α * = Γ -1 [d] A [d] y l . ( 41
)
This is then straightforward to compute

G * l (d) = G * l (Ψ [d] α * ) by G * l (d) = -A [d] y l 2 + y l 2 (42)
This means that when using the square loss, it suffices to extract the eigenvectors and eigenvalues of K to get a very simple expression for the f l 's and the G * l (d). This is the reason why we have chosen to implement this specific strategy in the numerical simulations presented below.

Discussion

Here we discuss a few interesting points of our proposed multiclass learning approach.

How does the learning procedure scale? One of the prominent questions related to kernel methods is to know whether they can adapt to large-scale datasets, as such types of datasets are plentiful in realworld applications. In order to answer this question, it is interesting to observe that the bulk of computations for MKPM is in (i) computing the family Ψ (for instance, KPCA has a O(N 3 ) time complexity) and (ii) evaluating the G * l (d)'s. As for the first point, many efficient methods exist to (help) perform approximate kernel principal component analysis: among the most popular, we may cite the Nyström method [START_REF] Drineas | On the nyström method for approximating a gram matrix for improved kernelbased learning[END_REF][START_REF] Williams | Using the Nystrm Method to Speed Up Kernel Machines[END_REF], and, the incomplete Gram-Schmidt decomposition [START_REF] Cristianini | Latent Semantic Kernels[END_REF][START_REF] Hussain | Theory of matching pursuit[END_REF]. These methods greatly reduce the cost of extracting a family Ψ as they run in O(v 2 N ) time where v N is a user defined parameter; meanwhile, the informativeness of the extracted family is almost as good as that of the family of principal components. As for the computations of the G * l (d)'s, which essentially depends on the largest dimension D to explore, we anticipate that this process is not this critical. Indeed, if the family Ψ is ordered in such a way that the first vectors explain most of the inertia of the φ(x i )'s, then a warm-start strategy that uses the solution computed for S d as a starting point for the optimization procedure when considering S d+1 should be very efficient.

What about the regularization path? We have mentioned in several occasions that d may be selected according to either theoretically-grounded criteria or empirically based estimates of the generalization error. We would like to elaborate on the empirical approach -the theoretical approach is related to the next question. We actually want to emphasize the fact that with our approach to multiclass Kernel Projection Machines, the regularization path is automatically computed at the time of learning. This means that the learning procedure is such that it produces all the multiclass predictors associated with dimensions 1 to D, even when only the predictor associated with D is sought. This comes from the dynamic programming approach described in Proposition 1 and it is a very nice feature of the learning process. It makes it possible to implement model selection procedures based on hold-out sets very easily without having to subsample the space of possible dimensions.

What are the generalization properties of MKPMs?

This is a question that goes beyond the scope of this paper. However, building on results presented in [START_REF] Blanchard | Finite-dimensional projection for classification and statistical learning[END_REF][START_REF] Vito | Learning from examples as an inverse problem[END_REF], we think that our model should inherit the nice adaptive properties of KPMs. On a closely related side, we would like to mention that the multiclass predictors learned with large values of d exhibited overfitting behaviours (not shown here). This suggests that d is precisely the complexity parameters associated with our model, just as it is the case for classical (i.e. not multiclass) KPMs.

Numerical results

For all numerical experiments described in this section, the model selection criterion used in the learning phase was to choose the value d that minimized the classification error rate on the training data. If a leave-one-out cross-validation was performed, this value was therefore different for each fold/iteration. Furthermore, we used output codes of length L = Q, where c q (l) = +1 if q = l and c q (l) = -1 if q = l. Finally, we used the Euclidean distance in the prediction rule defined in [START_REF] Haxby | Distributed and Overlapping Representations of Faces and Objects in Ventral Temporal Cortex[END_REF]. The criterion for measuring the performances of an algorithm was the rate of correct classification.

Benchmark datasets

We first applied our algorithm on three multiclass datasets issued from the Statlog database, the dna, segment and satimage datasets (downloaded from the libsvm website [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]), chosen because they are of medium size, thus limiting computing time, and because they were also included in the benchmarking study [START_REF] Hsu | A comparison of methods for multiclass support vector machines[END_REF]. In order to compare our results with the ones from this study, we limited ourselves to using the same RBF kernel K(x i , x j ) = e -γ||xi-xj || 2 , with different values of γ around the ones that were reported to be optimal γ o ; namely, we used 2 -v γ o , with v ∈ {-2, -1, 0, 1, 2}. We arbitrarily decided to use N = 1500 to further limit computing time.

The results of the simulations are reported in Table 1, where the last column represents the median performance of the five Multiclass SVM methods evaluated in [START_REF] Hsu | A comparison of methods for multiclass support vector machines[END_REF]. Overall, the classification rates obtained with our MKPM model were equivalent to the SVM performances, always for v = 0. Those performances might be improved by using more appropriate and possibly continuous output codes [START_REF] Crammer | On the learnability and design of output codes for multiclass problems[END_REF][START_REF] Dietterich | Solving multiclass learning problems via error-correcting output codes[END_REF] but the room for improvement on the segment and dna problems is potentially very small: i) as just noted, the results are comparable with those of SVM and ii) those SVM results have been obtained thanks to a very thorough validation procedure to choose the hyperparameters and they might be looked at as 'near optimal'. Finally, note that for the satimage dataset, the chosen d was very close to N , suggesting that considering higher values of N might entail better classification rates.

Functional MRI data

We now focus on one particular problem that can be formulated as a multiclass learning problem: the classification of functional Magnetic Resonance Imaging (fMRI) spatiotemporal data. This problem is sometimes refered to as "brain-reading" since it consists in trying to infer, from the spatial pattern of brain activity, what the subject was doing at a certain point in time, among a number of choices (the classes) defined by the experimental design. We prefer the more "down-to-earth" Multi-Voxel Pattern Analysis appelation (MVPA, see a review [START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fmri data[END_REF]); indeed it calls for using pattern recognition techniques to analyse fMRI data in a multivariate (multi-voxels) framework, thus possibly revealing effects that would not be detected with a standard univariate (in the sense that each voxel is processed independently) SPM-like approach [START_REF] Friston | Statistical parametric maps in functional imaging: A general linear approach[END_REF]. Despite the exponentially growing number of studies using this framework, KPMs have never been used to analyse fMRI data to our knowledge.

We used fMRI data from [START_REF] Haxby | Distributed and Overlapping Representations of Faces and Objects in Ventral Temporal Cortex[END_REF], which is the study that originally raised the interest of the neuroimaging community towards these approaches. The data of one of the six subjects that were scanned is available online at www.pni.princeton.edu/mvpa/downloads/ nifti_set.tar.gz. It contains ten acquisition sessions, each with examples of the eight categories of visual stimuli that were presented to the subject, which define eight classes. The inputs are points in a 577-dimensional space (the number of voxels in a cortical region defined anatomically, that was of particular interest for the authors; note that no feature selection was performed), whose coordinates are the raw activation values measured at single time-points (no temporal averaging was performed). We fed those examples both to our MKPM model and to a linear SVM (which is one of the algorithms giving the best performing according to the comparison study [START_REF] Misaki | Comparison of multivariate classifiers and response normalizations for pattern-information fmri[END_REF]), and performed a leave-one-session-out cross-validation to estimate the generalization capability of both methods.

In practice, we used the implementation of the linear SVM available in libsvm ( [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]), which uses a "one-againstone" method for the multiclass problem. We ran it with C = 10 -t (t integer so that -6 < t < 4) and obtained the best performances for t = -2 (classification rate = 0.681). To test our MKPM model, we used N = 577 (the number of features) and tried different standard kernels (linear, polynomial of different orders, and RBF with different γ). The best performances were obtained with the linear kernel (classification rate = 0.734), confirming that it is the most adapted for fMRI data; it is also to be noted that the order d of the selected model varied significantely across the iterations of the cross-validation (when using the linear kernel, it varied between 303 and 431, with a mean value of 390). We believe that our model outperforms the "one-against-one" SVM mostly because it relies on a true multiclass formulation; this should be tested on one hand with a multiclass SVM and on the other hand on different datasets, because if confirmed, the superiority of multiclass models would be an important result for the fMRI field.

Conclusion

In this paper, we have proposed MKPM, a multiclass extension of the kernel projection machines. The optimization problem posed by our approach is a combinatorial problem that can be nicely tackled by means of a dynamic programming procedure. We have conducted several numerical simulation on pattern recognition problems, from which MKPM appears to be very competitive with state-of-the-art methods. It seems that the co-regularization scheme implemented by our strategy is very relevant, but theoretical arguments need still be elaborated to support this feeling. This is going to be the topic of further research.

) 2 .

 2 Pick the (global) dimension d, the corresponding vector δ( d) = [ δ1 ( d) • • • , δL ( d)] of dimension L and the corresponding L predictors f1, d, . . . , fL, d that realize R( δ( d)) according to some model selection criterion.

  which, combined with (26), gives that ∆λ (d) = ∆ λ (d) and finishes the proof. Proposition 1 provides us with a way to compute ∆ λ (d) for all λ ∈ [L] and d ∈ [D] 0 . As a matter of fact, note that the computation of ∆ L (D) directly implies the computation of all the values of ∆ λ (d) for all λ ∈ [L] and d ∈ [D] 0 . If the G * λ (d) are precomputed, then the space complexity of evaluating ∆ L (D) is O(LD) and the time complexity is O(D 2 L).

Table 1 .

 1 Results obtained with our MKPM model and SVM.

	Dataset	N	d	Kernel	MKPM SVM
	segment 1500 451 γ = 2 -3	0.956	0.974
	dna	1500 690 γ = 2 -6	0.957	0.956
	satimage 1500 1492 γ = 2 -0	0.907	0.913
	fmri	577	390	linear	0.734	0.681

Here, we pervert the name Kernel Projection Machine as, when introduced in[START_REF] Blanchard | Finite-dimensional projection for classification and statistical learning[END_REF], it was specifically associated with the binary classification problem and the hinge loss. We provide the name with a broader meaning.


Dr Andrei V Nikitin
email: andreinikitin076@mail.ru

Gsma

New efficient algorithm for the calculation of energy levels of AB 3 type molecules

Keywords: symmetric-top molecules, variational calculation, irreducible tensor, finite basis representation, vibration energy levels symmetric-top molecules, variational calculation, irreducible tensor, finite basis representation, vibration energy levels

New efficient algorithm for the calculation of energy levels of AB3 type molecules

I. INTRODUCTION

The theory and computational modeling of molecular vibrations is on the top of interest in the molecular spectroscopy and related fields. Advances in computation of vibration-rotation spectra of three-atomic molecules (see ref. [1][2][3] and references therein) resulted in an interest to computations of spectra for polyatomic molecules . For more information on this subject, the reader is referred to the recent review [21]. An important part of polyatomic molecules are high symmetry molecules [4][5][6][7][8][9][10][11][12][13][14][18][19][20][21][22][23][24][25][26][27]. For the ABC 3 type molecules variational calculations using a finite basis representations and full C 3v symmetry have been applied to CH 3 F [10,11] and CH 3 Cl [27] molecules. In this paper C 3v and T d point groups are considered to be isomorphic to finite permutation groups (permutations of three or four atoms). In spite of a considerable number of publications devoted to computation of vibration-rotation spectra of the AB 4 and the AB 3 type molecules, they often do not completely use the molecular symmetry [4,5,8,9,19,25,26]. One variational method for calculating excited bending states of symmetric tetrahedral pentaatomic molecules based on the use of Radau coordinates and Jacobi polynomials as the basis functions has been presented in [30,31]. Only the technique based on the normal or symmetry coordinates does not result in special difficulties related to a construction of the symmetry basis [13,17,18,21]. See also rovibration calculations by Rey et al. [24] using rectilinear normal coordinates and the normal mode Hamiltonian expressed in terms of irreducible tensor operators.

When using the internal coordinates, such as three angles between the H-P-H bonds in PH 3 , in the general case, the basis functions constructed on the base of these coordinates are not orthogonal, even though those angles are independent. The construction of the symmetry basis using the q ij angles is a trivial task, but a problem is to work out an effective technique for calculations of multidimensional integrals using this basis. In the present work, a finite basis representations approach was applied for the AB 3 type molecules and simple formulas for calculations of the matrix elements of the KE operator and PES are reported . As a rule, we will consider only the angular part of the matrix elements, because the symmetry properties of radial part are obvious. Here we also assume that the molecule has only one deep minimum.

Tennyson et al. [32], Gatti et al. [33,34], Mladenovic´ [35][36] have presented derivations of kinetic energy operators for treating rotations and vibrations of polyatomic molecules using polyspherical coordinates. In these coordinates, the N nuclear position vectors are transformed to N-1 internal vectors and the nuclear center of mass vector, and the internal vectors are parametrized by spherical polar coordinates: Ri , θi , φi , i=1,...,N-1. Then, a body-fixed coordinate system is introduced, defined by aligning the body-fixed z axis along the first vector, and placing the second vector in the body-fixed xz plane. The remaining angles θi, φi are either referenced with respect to the body-fixed z axis, or other vectors.

At the same time, the symmetric form of vibration KE operators is known [30,31,37,38].

Instead of using the torsion angle, symmetric form use additional angles between molecular bonds. This form of KE operators is more tailored to our basis functions, and so we use it in the present work. Symmetric form is especially interesting because it has no sin(q) -2 singularity while the KE operator of Ref. [35] contains this singularity. When constructing the vibration basis set using the products of functions that depend on the bending q and torsion t angles, the angular basis f n (q) must be chosen to be vanishing at sin(q)=0. From a physical point of view, this choice of f n (q) is not always correct because in some vibration states, a molecule can be in a geometric configuration with the angle q=π with rather high probability.

Described in this work algorithm of calculation of vibrational energy levels uses nonorthogonal basis set of three interbond angles q ij . Of course, the non-orthogonal basis is not always easy-to-use, but this type of basis is used, for example, in the chemical reactions theory [39]. When calculating the matrix elements of PES, whenever possible, we try to keep the symmetrical form with respect to three angles. We transform to polar coordinates only at the final phase of the calculations. Similarly, when calculating the matrix elements of the KE operator in symmetric form, we first calculate the derivatives with respect to angles using matrix element symmetry, and only after this, we calculate the matrix element by using polar coordinates.

The present paper comprises seven sections. Section II describes an algorithm for constructing contracted angular basis set that uses three angles q ij between three identical atoms, and the norm calculation for this basis set. Sections III and IV are devoted to algorithm of Section V briefly outlines the general tree-like coupling scheme of PES and basis functions, a construction of the symmetry-adapted angular basis set, and algorithm for calculation of vibrational matrix elements for AB 3 molecules. In section VI we report the tests of the convergence of calculated vibrational PH 3 levels in the 0-7000 cm -1 range. Section VII is devoted to discussion and conclusion.

II. CONTRACTED ANGULAR BASIS FOR THREE IDENTICAL ATOMS AND NORM

CALCULATION

In the present work, the following coordinates were chosen as independent ones: three bond lengths {r 1 , r 2 , r 3 } and three interbond angles {q 12 , q 13 , q 23 }. For one-dimensional basis functions, the following functions were used:

∑ = l s l nl n q Y c q f ) 0 , ( 2 ) ( , π
, where Y are spherical harmonics written in the form of Ref. [START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF] that ensures that

π 4 1 ) 0 , ( 0 , 0 = q Y and 1 ) cos( ) 0 , ( 2 0 2 , = ∫ q d q Y s l π π . The normalized eigenfunctions 1 ) cos( ) ( 23 23 0 2 = ∫ q d q f n π
of one-dimensional angular problem can be used to obtain f n . Initial angular basis can be written in the form:

) ( ) ( ) ( 23 13 12 3 2 1 3 2 1 q f q f q f n n n n n n = . ( 1
)
This basis is not orthogonal one. The norm is calculated using the following three-dimensional integral:

∫ ∫ ∫ >= < π π π 0 0 2 0 23 13 12 23 23 13 13 12 12 3 2 1 3 2 1 ) cos( ) cos( ) ( ) ( ) ( ) ( ) ( ) ( | 3 3 2 2 1 1 dt q d q d q f q f q f q f q f q f k k k n n n k n k n k n
, where t 23 is the torsion angle.

To calculate the integral

) , ( ) , ( 2 ) ( ) ( 2 2 1 1 2 1 dt s q Y s q Y c c dt q f q f s l l l s l ml nl m n ∫ ∑ ∫ = π π
π the following expressions [START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF] were used:

∑ + + + = LM LM s l s l q Y LM s l s l C L l l C L l l q Y q Y ) 0 , ( ) , , ( ) 0 , 0 , 0 ( ) 1 2 ( 4 ) 1 2 )( 1 2 ( ) 0 , ( ) 0 , ( 23
∑ + = + m lm lm l t q Y q Y l t q q q q Y ) , ( ) 0 , ( 1
Therefore, from Eqs. (2, 3), we get:

∑ ∑ = L m Lm Lm L l l l l t q Y q Y K q Y q Y ) , ( ) 0 , ( ) 0 , ( ) 0 , ( 23 13 12 23 0 23 0 2 1 2 1 , (4)
where

2 2 1 2 1 ) 0 , 0 , 0 ( 1 2 ) 1 2 )( 1 2 ( 2 1 L l l C L l l K L l l + + + =
, and ) , , (

2 2 1 1 LM m l m l C is Clebsch-Gordan coefficient. Using Eq. ( 4
)
we get ( )

) cos( ) cos( ) 0 , ( ) 0 , ( 2
) ( ) ( ) ( ) ( | 13 12 13 0 12 0 2 0 0 13 13 12 12 3 2 1 3 2 1 2 1 2 1 2 3 1 3 2 2 1 1 q d q d q Y q Y K с с q f q f q f q f k k k n n n L L L L l l l l l k l n k n k n ∑ ∑ ∫ ∫ >= < π π π . ( 5
)
Let us denote the one-dimensional integral by

∫ ∑ ∫ = =
( 2 1 2 1 2 1 1 1 1 1 1 1 q d q Y q Y q Y c c q d q Y q f q f J L l l l l l k l n L k n L k n cos( ) 0 , ( ) 0 , ( ) 0 , ( 2 ) cos( ) 0 , ( ) ( )
Using Eq. ( 2) and taking into account that spherical harmonics Y lm with the same m are orthogonal we get:

( )

L l l l l l k l n l l l k l n L k n K c c L L l l C L l l c c J , , 2 2 1 2 1 0 , , , 2 1 2 1 2 1 1 1 2 1 2 1 1 1 1 1 1 2 4 1 ) 0 , 0 , 0 ( ) 1 2 ( 4 ) 1 2 )( 1 2 ( ∑ ∑ + = + + + = π π (6)
Using Eqs. (5,6), one can derive

( ) ( ) ∑ + >= < L L k n L k n L k n J J J L k k k n n n 0 , , 0 , , 0 , , 3 3 2 1 3 2 1 3 3 2 2 1 1 1 2 1 2 | π π (7)
For example, if q q q P P P . Due to simple recursive formulas [START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF] for

) ( ) cos( 0 , q Y q n
, this considerably facilitates the calculation of matrix elements. Therefore,

∑ > < = 3 2 3 2 ) 0 , ( | ) ( cos | 2 ) ( ) ( cos , 2 3 l l s l P kl k P q Y s l q s l c q f q π .
If s=0, similarly to Eqs. (5,6), we get the following expression for the matrix element:

( )

) cos( ) cos( ) 0 , ( ) 0 , ( 0 | cos | 0 ) ( ) ( ) ( ) ( 2 ) cos( 13 12 13 0 12 0 2 3 0 0 13 13 12 12 2 3 2 1 23 3 2 1 3 1 3 2 1 2 3 1 3 2 2 1 1 q d q d q Y q Y K l l с с q f q f q f q f k k k q n n n L L L L l l P l l l l k l n k n k n P ∑ ∑ ∫∫ > < >= < π π π
In comparison with Eq. ( 5), this expression contains an additional summation over the l 3 index and the <l 3 |cos p |l 2 > matrix element. Using the expression

π π 4 1 2 0 | cos | 0 2 3 1 3 2 1 2 3 1 3 3 3 2 3 0 , , , + > < = ∑ L K l l c c J L l l l l l P l k l n P L k n (8)
one can write the matrix element in the form:

( ) ( ) ( ) ∑ ∑ ∑ + = > < >= < L L k n L k n P L k n L k n L k n L L l l P l l l l k l n P J J J L J J K l l с с k k k q n n n 0 , , 0 , , 0 , , 3 0 , , 0 , , 2 3 2 3 2 1 23 3 2 1 3 3 2 2 1 1 2 2 1 1 3 1 3 2 1 2 3 1 3 1 2 1 2 0 | cos | 0 2 ) cos( π π π
Taking into account that q 23 is independent of q 12 and q 13 , it is easy to obtain the following general expression for the PES matrix elements:

( ) ( ) ∑ + >= < L P L k n P L k n P L k n P P P J J J L k k k q q q n n n 23 3 3 13 2 2 12 1 1 23 13 12 0 , , 0 , , 0 , , 3 3 2 1 23 13 12 3 2 1 1 2 1 2 ) cos( ) cos( ) cos( π π (9)

IV. CALCULATION OF MATRIX ELEMENTS OF KE

In this section, we use the mass-dependent orthogonal coordinates{

} ′ ′ ′ ′ ′ ′ 1 2 3 12 13 23
r ,r ,r ,q ,q ,q . The polar coordinates are defined in a standard way via three vectors { }

{ } ′ ′ ′ ′ ′ ′ 1 2 3 12 13 23
r ,r ,r ,q ,q ,q are defined in a similar way via three vectors { } ′ r i r :

∑ = + = 3 1 ' j j i i r d r r r r r
, where

c d µ 3 1 3 1 3 1 - + - = with P H H c m m m + = 3 µ , ( 10
)
The most convenient expression for the kinetic energy can be obtained in terms of these coordinates, while the analytical representation using internal polar coordinates includes additional terms in the kinetic energy operator [35]. Mass-dependent coordinates (10) keep the same symmetry properties as initial polar coordinates. Let us use the KE operator in the form Ref. [30,31,37] where the KE operator does not comprise the singular terms of the sin(q) -2 form.

Therefore one can use the ) 0 , ( 0 , q Y n spherical harmonics as basis functions that do not vanish at q=π. The kinetic energy of the AB 4 -and AB 3 -type molecules can be expressed in the form [37]:

j N j j k j N jk jk V h g h T ℜ + ℜ ℜ = - ∑ ∑ - - ∂ ∂ ∂ ∂ ∂ 6 3 2 6 3 2 ) 2 1 /( . (11)
For radial coefficients

i ij ij m g δ = .
The angular diagonal coefficients of the g matrix are written in the form:

) cos( 1 1 ) ( sin 2 2 2 2 2 ij j j i i ij q r m r m q ∂ ∂         +
. For a pair of angles q ij q lk , the angular off-diagonal elements of the g matrix vanish unless i=l or j=k:

) cos( ) cos( ) cos( ) cos( ) cos( 2 2 ) cos( ), cos( jk ij j j ik jk ij q q q q r m q q q g jk ij ∂ ∂ ∂ + - =
. Note, that in the g matrix, all off-diagonal radial-angular elements vanish. The radial and angular coefficients of the h matrix are correspondingly

i i i r m h 2 = and ) cos( 1 1 ) cos( 2 2 2 ) cos( ij j j i i ij q q r m r m q h ij ∂ ∂         + - =
. For AB 4 molecules, the summation in ( 11) is performed over four radial coordinates r i and six angles q ij , while for AB 3 molecules, only three radial coordinates r i and three angles q ij are used. In this

F
) 0 , ( ) 1 ( ) 0 , ( ) ( 0 , 0 , 2 2 ij l ij l ij ij ij q Y l l q Y q q ctg q + - =         ∂ ∂ + ∂ ∂ .
Note that only the full sum of all diagonal and off-diagonal angular terms results in symmetric matrix elements, while each of these sums, e.g. the summation only over the diagonal or offdiagonal terms, results in asymmetric matrix elements. In the present paper, to calculate matrix elements, we use the following property of the reduced matrix elements [START_REF] Zhilinskii | Method of Irreducible Tensorial Operators in the Theory of Molecular Spectra[END_REF][START_REF] Champion | Spherical top spectra in Spectroscopy of the Earth's Atmosphere and Interstellar Medium[END_REF]:

> >=< < R L R L PF O P PF F O F
ˆ where P -is an element of the molecular symmetry group. Due to this property, for AB 3 molecules, it is not mandatory to calculate all three diagonal and offdiagonal angular matrix elements. It is sufficient to calculate only one simplest diagonal and one off-diagonal angular matrix element. We calculate the matrix element for the sum of three diagonal operators using the expression:

( ) { }

) cos( ) cos( ) 0 , ( ) 0 , ( ) 0 , ( ) 0 , ( ) 0 , ( ) 0 , ( ) 1 ( 2 ) ( 3 ' 3 2 ' 2 1 ' 1 3 3 2 2 2 1 3 3 3 3 2 2 2 2 1 1 1 1 dt q d q d q Y q Y q Y q Y q Y q Y l l c c c c c c k k k q q ctg q n n n l l l l l l l l l l l l l k l n l k l n l k l n ∫∫ ∫ ∑ ′ ′ ′ ′ ′ ′ + ′ ′ = ∂ ∂ + ∂ ∂ π π π π Using the designation ( ) [ ] π 4 1 2 1 3 3 3 3 3 3 3 3 3 3 3 3 2 0 , , , + + ′ ′ = ′ ′ ′ ∑ L K l l c c J L l l l l l k l n D L k n , ( 12
)
we can write the diagonal matrix element it the following form:

( ) ( ) ∑ + = ∂ ∂ + ∂ ∂ L D L k n L k n L k n J J J L k k k q q ctg q n n n 2 0 , , 0 0 , , 0 0 , , 3 3 2 1 23 23 23 2 2 3 2 1 3 3 2 2 1 1 1 2 1 2 ) ( π π (13)
It is more convenient to calculate the off-diagonal matrix element using the polar coordinates.

Let us show, that the off-diagonal matrix element can be expressed in the following form:

( ) ∑ + = ∂ ∂ ∂ ∂ L L k n D L k n D L k n L J J J k k k q q t n n n ) 1 2 ( 1 2 ) cos( 0 , , , 1 1 , , , 1 1 , , , 3 3 2
1 13 12 23 3 2 1 3 3 2 2 1 1 π π , ( 14
)
where 1 1 ,

1 1 D L k n J
is a coefficient similar to the coefficients used in Eqs. ( 6), ( 8), (12). 4) for product of J functions depending on q 23 , integrate Eq. ( 14) over t 23 , and using the known formula for spherical harmonics [START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF] ) , ( )

- = - q Y q Y m l m m l 1 ( ) , ( , , ϕ ϕ -
, we get: ( )

) cos( ) cos( ) 0 , ( ) 0 , ( ) ( ) ( ) ( ) ( 2 ) cos( ) cos( ) ( ) ( ) ( ) cos( ) ( ) ( ) (
q d q d q Y q Y q f q f q q K q f q f c c dt q d q d q f q f q f q q t q f q f q f L L k k l l L L l l n n l k l n k k k n n n         ∂ ∂ ∂ ∂ × =         ∂ ∂ ∂ ∂ ∑ ∫ ∫ ∑ ∫ ∫ ∫ π π π π π π
If we denote the one-dimensional integrals in the right-hand member of this expression via

) cos( ) 0 , ( ) ( ) ( 12 12 1 , 12 12 12 1 1 , 1 1 1 1 q d q Y q f q q f J L k n D L k n ∫         ∂ ∂ =
, we get Eq. ( 14).

One can simplify the expression for the

) ( ) 1 ( ) 1 ( ) ( ) ( ) ( 1 , , , q Y m m l l q Y q mctg q Y q m l m l m l + + - + + = ∂ ∂ , we get: ) cos( ) 0 , ( ) 0 , ( ) 1 ( ) 0 , ( ) 2 ( 12 12 1 , 12 1 2 2 12 0 1 1 , 2 1 2 1 2 2 1 1 1 1 q d q Y q Y l l q Y c c J L l l l l l k l k D L k n ∫ ∑ + = π
. Using Eq. ( 2) with s 1 =0, s 2 =1, and taking into account the orthogonality of spherical harmonics, we get

) 1 , , 1 , , 0 , ( ) 0 , , 0 , , 0 , ( ) 1 2 ( 4 ) 1 2 )( 1 2 ( ) 1 ( 2 1 2 1 2 1 2 2 1 1 , 2 1 2 2 1 1 1 1 L l l С L l l С L l l l l c c J l l l k l k D L k n + + + + = ∑ π
For even values of l 1 +l 2 +L, it can be shown [START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF] that:

2 / 1 2 2 1 1 2 2 2 1 2 1 )] 1 ( ) 1 ( [ 2 ) 1 ( ) 1 ( ) 1 ( ) 0 , , 0 , , 0 , ( ) 1 , , 1 , , 0 , ( + +
+ - + + + = l l L L l l L L l l L l l С L l l С Since the ) 0 , , 0 , , 0 , ( 2 1
L l l С coefficients are not equal to zero only at even values of l 1 +l 2 +L , one can simplify the expression obtained for J in the following way:

2 / 1 2 / 1 1 1 2 2 1 1 , ] 1 2 [ )] 1 ( [ 4 ) 1 ( ) 1 ( ) 1 ( 2 1 2 1 2 2 1 1 1 1 + + + - + + + = ∑ L L L l l L L l l K c c J L l l l l l k l k D L k n π (15)

V. FULL SYMMETRY BASIS AND CALCULATION OF VIBRATIONAL LEVELS

The initial angular basis (1) can be easily symmetrized and tailored for calculations that use a symmetrical basis. Using the standard operators, one can construct the full symmetry basis set of functions:

S kG S n n n P kG n n n ... | ...
q q q f G BAB D and ) , , ( 3 2 1 , 3 r r r f G AB D
, respectively. The vibrational basis functions and the PES operators are expressed in terms of irreducible tensor operators [START_REF] Zhilinskii | Method of Irreducible Tensorial Operators in the Theory of Molecular Spectra[END_REF][START_REF] Champion | Spherical top spectra in Spectroscopy of the Earth's Atmosphere and Interstellar Medium[END_REF]. Each tensor is associated to a binary tree [START_REF] Nikitin | [END_REF]. The binary trees of vibrational coupling for AB 3 molecules are shown in Fig. 1. Each branch of the tree is characterized by its symmetry.

The binary tree of the PES parameters is shown in the left part of the figure while the binary tree of the vibrational basis coupling is shown in the right part of this figure. In order to determine the PES operators of the A 1 symmetry type, formed from the symmetrized coordinates [24,25] b a

E E A b E a E A SA SA SA SR SR SR , , , , , 1 1
we used the two step procedure. In first step, we construct the symmetrized power of the coordinates

∑ = i n i i n r r R K r r V ) ,... ( ) ,... ( 1 1
. The complete vibrational basis functions can be constructed by coupling the symmetrized basis functions of all coordinates into binary trees. For AB 3 and AB 4 molecules [26,27] (18) and the recoupling matrix elements formula:

σ σ σ ψ ψ ψ ψ σ σ σ C C C F V V C C C C C C
' ' 2 2 ' 2 ' ' 1 1 ' 1 ' ' 2 ' ' 1 2 1 ' 2 ' 1 2 2 2 1 1 1 ' ' 2 ' ' 1 2 1 ' 2 ' 1 ' ' 2 1 2 1 ' 2 1 ' ' ' ] ' ' ][ ][ ' [ ) ( ) ( ) ( C C C C C C C C C C V V C C C C C C C C V V ψ ψ ψ ψ Γ Γ Γ Γ = , ( 19
)
where the indices C and Γ denote point group irreducible representations (all other indices are omitted for the sake of simplicity) a multi-dimensional matrix element can be expressed as a sum of the one-dimensional matrix elements products.

The use of a nonorthogonal angular basis does not complicate the procedure for obtaining eigenvalues. There are two methods of solution eigenvalue vibrational problem. In first method, one can solve the eigenvalue angular problem, and then use the resulting orthogonal eigenvectors to construct the symmetrized orthogonal angular basis. In second method, one can use a nonorthogonal vibrational basis set. To solve the symmetric-definite generalized eigenvalue problem, standard computational programs, such as programs available in LAPACK library can be applied. Note that this technique does not consume too much additional main memory because

VI. EXAMPLE OF APPLICATION TO PH 3 MOLECULE

In the present work, we employ potential in the mass-dependent orthogonal coordinates that has been determined in the paper Ref. [25]. One-dimensional eigenfunctions were computed from the multi-dimensional PES by fixing all other coordinates to the equilibrium values. For the stretching coordinates, the following Morse type function was used:

)] ' ' ( exp[ 1 ) ; ' ( e i r r a a r f - - - = (20)
where a=1.9. This value of the a parameter ensures that the second order term of the potential provides a reliable representation for the one-dimensional stretching cut. The terms of higher orders result in relatively small corrections. For the interbond angular coordinates, the following functions were used: e ( ) cos( ) cos( ) gives a final set of the 6D expansion terms. Two types of localized one-dimensional wave functions were constructed: Q n (q ij ) -functions of bending angles q ij , and R n (r i ) -functions of radial coordinates r i . One can obtain the symmetrized vibrational basis set by coupling the

ij ij q q q φ = - (21)
A A G n n n σ , 23 13 12 and R R G n n n σ , 3 2 1
basis functions using Clebsch-Gordan coefficients of the C 3v symmetry group.

To ensure a rapid computation of the irreducible matrix elements, we applied the following technique. We first calculate all J-coefficients (6,8,12,15). After this, we calculate angular and radial three-dimensional matrix elements and store results into main memory. It is not necessary to use an optimization of calculations to get the radial matrix elements. At the same time, one needs to use an optimization of calculations to get the angular irreducible matrix elements.

Several convergence tests involving calculated vibration levels have been already performed using various analytical PES representations and C 2v symmetry [25]. A good agreement between vibration levels calculated using two PES fitted in orthogonal massdependent coordinates (10) and in internal mass-independent coordinates have been already reported in Ref. [25]. In this section, we check only the full symmetry nonorthogonal vibrational basis set convergence in orthogonal mass-dependent coordinates (10). Calculations with increasing dimensions for the basis cut-off suggest that our vibration energies up to 7000 cm -1 are converged in average to 0.01 cm -1 or better. The corresponding residuals increase gradually with energies. Similarly to Ref. [25], in the present paper, we use the nonsymmetrized vibrational basis dimension as the basis dimension. All calculations have been performed using a full symmetry-adapted set of basis functions. The used basis functions have been grouped into the A 1 , A 2 , and E symmetry classes. Let us consider four basis sets referred to as I, II. III, IV. The dimensions of those basis sets are 14643, 16047, 17889, and 19230 functions, respectively. The standard deviation between energy levels calculated using the basis sets I-IV, I-IV, III-IV is 0.034482, 0.002616, and 0.000571 cm -1 , respectively, up to 7000 cm -1 . Because all lower calculated energy levels show a slight difference with respect to those obtained in Ref. [25] only four observed band centers in the region 6700-7000 cm -1 were considered. Comparison of these observed band centres of PH 3 molecule and calculated from PES [25] in four basis sets I, II, III, IV given in Table 1. The deviations between vibration levels calculated using vibrational basis set IV and basis sets I, II, and III are shown in Fig. 2. Similarly to the results of calculations using internal coordinates [25], the main contribution to this standard deviation results from two levels above 6500 cm -1 shifted to 0.3 and 0.5 cm -1 (these levels are not shown in Fig. 2) and without these levels the standard deviation between all set is less 0.01 cm -1 . It is possible, that the considerable deviation for these two levels is resulted from wrong asymptotic behavior of the potential when the molecule is far from its equilibrium geometry. In general, one can come to the conclusion that calculations using the nonorthogonal C 3v basis and the basis reported in the Ref. [25] result in the similar convergence behavior. For example, according to results from Ref. [25], the standard deviation between two vibrational calculations using 12000 and 15000 nonsymmetrised basis functions is 0.048 cm -1 . At the same time, the use of the full symmetry provides that the computation speed of vibration levels is several times faster than for

VII. DISCUSSION AND CONCLUSION

The primary motivation for this study was to demonstrate a good convergence of new angular basis and to check algorithm of the matrix elements calculation. To do this, we dealt only with a solution of vibrational problem using orthogonal coordinates. At the same time, one can calculate the matrix elements of the vibration-rotation and nonorthogonal terms of the Hamiltonian in a similar way. For example, the matrix elements ( 14) have been calculated in polar coordinates. The functions of the torsion angle used to write vibration-rotation and nonorthogonal terms of the Hamiltonian can be expressed in terms of three angles q ij . This allows one to obtain a symmetrical form similar to (11). In addition, the dependence of the PES functions used in Ref. [17] on sin(q ij ) presents no essential difficulties.

The form of angular basis set based on cos(q ij ) product removes the sin(q) -2 singularity from the vibrational kinetic energy operator in polar coordinates [35] also for AB 4 , ABC 3 type molecules. However, the matrix elements calculations are complicated by the presence of a redundancy among the six angles q ij [19,10] for these molecules. On the other hand, symmetrisation of these basis functions can be achieved easily via projector operators technique.

For the CH 4 molecule, the study of spectra in the range higher than 5400 cm -1 is of considerable importance and it is difficult to obtain a good convergence of the calculated high energy levels without using the T d symmetry. The paper [8,9,26] uses the C 3v symmetry. This means that only 6 from 24 permutations in T d point group are in use. In paper [22], only 8 from 24 elements of T d point group were used. Thus, the cited papers do not employ the full symmetry of the molecule.

Of course, one can try to perform a numerical symmetrization of the angular basis, but it is a The calculation of the angular matrix elements for AB 4 molecules is a rather difficult task.

The calculation procedures for the AB 3 matrix elements employ only the (l 1 ,0,l 2 ,0|L,0) coefficients, while the calculations of the AB 4 matrix elements employ the (l 1 ,0,l 2 ,M|L,M) coefficients. For high values of L and M, the calculation of the (l 1 ,0,l 2 ,M|L,M) coefficients becomes a difficult problem [45,[START_REF] Heim | vibrational matrix elements calculation" should be "the calculation of vibrational matrix elements". Reply> OK, changed p3, line 24, should "numbers" be "quantum numbers" or "atoms"? Reply> OK, changed to atoms[END_REF]29]. For the AB 4 type molecules the same coupling scheme (Fig 1) as for the AB 3 type molecules could be applied. For ABC 3 molecules [27], the complete angular basis function could be received by coupling the three-dimensional symmetrized angular function ) , , ( of two torsion angles [26,27]. The symmetrized form of the PES used in the present paper contains much less terms than a nonsymmetrized form. For example, the 6 th order PES expansion for AB 3 , AB 4 and ABC 3 molecules contains only 196, 287 and 967 parameters, respectively. At the same time, the number of terms in nonsymmetrized forms [27] is 20-40 times greater than that in symmetrized forms. This proves that computations using symmetry-adapted forms of PES are very important.

In our future works, we will try to apply these computations to AB 4 and ABC 3 molecules. Binary tree of vibrational coupling for AB 3 molecules. In the left part: tree of the PES parameters.

In the right part: tree of the basis function coupling.

Figure 2

Deviations between vibration levels of PH 3 molecule calculated in various vibrational basis sets. In the Introduction (and possibly later) the reference list might include the papers of Kozin et al (2003-5) which also consider and apply some of the topics considered in this paper (orthogonal coordinates, contracted angular basis functions, efficient calculation of matrix elements ...) in the context of 4-atomic molecules. Reply> I added two references [28,29] Also in the first paragraph of the Intro, references [34,35] could be cited in the context of the symmetry of AB4 systems.

Reply> I added next sentence to the introduction: One variational method for calculating excited bending states of symmetric tetrahedral pentaatomic molecules based on the use of Radau coordinates and Jacobi polynomials as the basis functions has been presented in [30,31]. p4, line 43, which "New Algorithm"? The one to be described in This Work? Reply> Replaced on "Described in this work algorithm p5, line 54, has "t23" been defined? Reply> I added the definition p7, line 8, insert "angular" before "PES" Reply> OK p7, line 8, should "the power functions" read "a sum of product functions"? Reply> OK p8 (and elsewhere), I think the symbol q is being used in too many different and confusing ways: general coordinates with 1 subscript (as in equation 11) and also interbond angles with 2 subscripts (as in line 3) or no subscripts (as in line 36 , what are the "reduced matrix elements"? (They are defined later, but I think need some explanation here or earlier). Reply> I just added two references because there are explanation in next paragraph p11, line 3, "The not symmetry-adapted" should be "The initial" (and I think everything after "symmetrised" could be deleted in that sentence!).

Reply> OK, changed p11, line 21, there are 2 closing brackets missing in the first formula. Reply> OK, changed p11 & fig 1, the "binary tree" seems trivial -maybe all is explained in ref [START_REF] Zhilinskii | Method of Irreducible Tensorial Operators in the Theory of Molecular Spectra[END_REF] but perhaps a further sentence of explanation here would make the significance of this clearer? Reply> I agree that figure 1 is trivial. But readers very often don't understand coupling scheme. I added some details in the text. Figure 1 is also updated. p11, line 46, "the symmetrized coordinates Si" (as distinct from symmetrised basis functions) appear suddenly and without definition -again some explanation might help the reader. Reply> Si replaced and text is changed p13, line 3, "In the sequel" is not clear (does it mean "next" or "in practice" ...?) Reply> Text changed p15, lines 20-27, apparently a sentence is repeated but with different numbers! Reply> Deleted one sentence p15 lines 27-9 and Table 1, why only compare with 4 observed band centres? What makes these special? Are there no other observations? Also, how do the calculated vibrational levels compare with other variational calculations on this PES? Reply> I added one sentence to text. Calculated levels are very closed to [25]. The purpose of article just demonstrate that the algorithm work. Orthogonal and non orthogonal surfaces [25] are a little different so direct comparing is not informative. p15, the speculation about "some peculiarities of the potential" seems very vague -can the author be more specific? Can he say at least how the two large deviations vary with basis set?

  New efficient algorithm for the calculation of energy levels of AB3 type molecules

  angular matrix elements of PES and KE in contracted angular basis set.

r

  (i=1,2,3) linking the centre of mass with three identical atoms. Three stretching coordinates r i are the lengths of these vectors angular coordinates. The mass-dependent orthogonal coordinates[35]

  formula(14), we use Eq. (

  following formula[START_REF] Varshalovich | Quantum theory of Angular momentum[END_REF]:

  For example, for AB 3 molecules, we use the following six projection operators(17) to construct symmetrized three-dimensional angular

  step, we couple the symmetrized powers of different symmetrized coordinates into irreducible trees. A set of all possible trees of the A 1 symmetry type gives a final set of trees. The potential function was developed in power series of irreducible tensors

  , it is necessary to perform only one coupling of stretching basis function and bending basis function (

  is used only to transform the V(PES)+T(KE) matrix. The eigenvalues can be obtained then from the transformed V+T matrix. As a rule, the use of the full symmetry properties allows one to considerably decrease the dimension of the V+T matrices. The general technique used for the vibrational energy levels calculation includes several steps. Step 1: the solving of one-dimensional vibrational problems and calculation of the J coefficients like(6,8,12,15) for angular problem. Step 2: the symmetrization of wave functions and calculation of irreducible matrix elements for angular and radial trees in terms of a sum of products of onedimensional integrals. Step 3: Eq. (19) is recursively used for obtaining sub trees of basis functions until the subsequent sub tree becomes a tree representing only the same type coordinates or a leaf. The 9G symbols are saved in the memory. Vibrational reduced matrix elements are stored on disk separately for each type of symmetry. Step 4: the solution of the generalized eigenvalue problem.

  projection operators(17) for the group of permutations of three H-atoms. The subindices correspond to irreducible representations of the C 3V point group. In order to determine the six-dimensional (6D) expansion terms of the A 1 symmetry type, we applied the two step procedure described above. First step corresponds to the construction of the symmetrized powers, and second step to the coupling of the symmetrized powers of different symmetrized coordinates in irreducible trees. A set of all possible trees of the A 1 representation

  the C 2v basis. In addition, the computations using the full symmetry require less memory space that allows us to work with big basis sets.

  task, because an accurate solution of the vibration problem requires the use of the high excited angular basis functions.

  Figure captions

  Binary tree of vibrational coupling for AB3 molecules. In the left part: tree of the PES parameters.In the right part: tree of the basis function coupling. 215x279mm (600 x 600 DPI)Deviations between vibration levels of PH3 molecule calculated in various vibrational basis sets. 215x279mm (600 x 600 DPI)Reply to the Reviewer 1 (R1) of the manuscript TMPH-2010-0349 by A.V. Nikitin

Table I . Comparison of some observed band centers of PH 3 and calculated from PES[25] in four basis sets I, II, III, IV. All values are in cm -1 .

 I

	PH 3 band	Obs. †	I	II	III	IV
	3ν 3 (E)	6714.60 6714.769	6714.771	6714.766	6714.766
	2ν 1 +ν 3 (E)	6883.73 6886.345	6886.348	6886.344	6886.344
	ν 1 +2ν 3 (E)	6890.86 6890.122	6890.125	6890.121	6890.120
	3ν 3 (A 1 )	6971.16 6971.360	6971.362	6971.359	6971.358
	† See Ref [44].

1 basis sets:

ψ ψ ψ ψ Γ

    

Acknowledgement

The support of ANR(France) through Grant No. ANR-08-BLAN-0254-04 is acknowledged. The support of the Groupement de Recherche International SAMIA between CNRS (France), RFBR (Russia) and CAS (China) is acknowledged. We acknowledge the support from the computer centre Reims-Champagne-Ardenne and the computer centre SKIF Syberia of Tomsk. The support of RFBR(Russia) through Grants 09-05-92508 and 09-09-93103 is acknowledged. The support of CRDF(USA) through Grant RUG1-2954-TO-09 is also acknowledged.

in Fig 1). We used the following standard definition of the direct product of irreducible tensors [START_REF] Zhilinskii | Method of Irreducible Tensorial Operators in the Theory of Molecular Spectra[END_REF][START_REF] Champion | Spherical top spectra in Spectroscopy of the Earth's Atmosphere and Interstellar Medium[END_REF]:

where C,C′,C′′ are irreducible representations, σ, σ′, σ′′ are their rows, and ( )

symbols corresponding to Clebsch-Gordan coefficients of the C 3v symmetry group. Here [C] stands for the dimension of C. The recoupling scheme used for calculations of matrix elements is similar to one described in ref. [START_REF] Nikitin | [END_REF]. Using the standard definition of the reduced matrix elements [START_REF] Zhilinskii | Method of Irreducible Tensorial Operators in the Theory of Molecular Spectra[END_REF][START_REF] Champion | Spherical top spectra in Spectroscopy of the Earth's Atmosphere and Interstellar Medium[END_REF]: