Pierre Mahé
email: pierre.mahe@xrce.xerox.com

Jean-Philippe Vert
email: jean-philippe.vert@ensmp.fr

Virtual screening with support vector machines and structure kernels

Support vector machines and kernel methods have recently gained considerable attention in chemoinformatics. They offer generally good performance for problems of supervised classification or regression, and provide a flexible and computationally efficient framework to include relevant information and prior knowledge about the data and problems to be handled. In particular, with kernel methods molecules do not need to be represented and stored explicitly as vectors or fingerprints, but only to be compared to each other through a comparison function technically called a kernel. While classical kernels can be used to compare vector or fingerprint representations of molecules, completely new kernels were developed in the recent years to directly compare the 2D or 3D structures of molecules, without the need for an explicit vectorization step through the extraction of molecular descriptors. While still in their infancy, these approaches have already demonstrated their relevance on several toxicity prediction and structure-activity relationship problems.

Support vector machines and kernels

. A pattern recognition problem occurs when one is given a finite set of objects that belong to two possible classes, and must learn from this training set a rule to automatically predict the class of objects with unknown class. This general and abstract formulation encompasses in fact a number of practical situations in chemoinformatics and beyond. We focus here in particular on situations where the objects available are small molecules, and the classes to be predicted represents various properties of interest such as toxic/non toxic, druggable/non-

Introduction

Computational approaches play an increasingly important role in modern drug discovery. In particular, accurate predictive models accounting for the biological activity and drug-likeliness of candidate molecules can help in the identification of promising molecules and screening for various side-effects, leading to substantial savings in terms of time and costs for the development of new drugs. Such predictive models aim at inferring a relationship between the structure of a molecule and its biological and chemical properties, including toxicity, pharmacocinetics and activity against a target. The development of high-throughput technologies to assay such properties for large numbers of candidate molecules, and the subsequent availability of increasing quantities of molecules with characterized properties, has triggered the use of statistical and machine learning approaches to automatically learn the structure-property relationship from these pools of characterized molecules.

Decades of research in machine learning and statistics have provided a profusion of methods for that purpose, ranging from classical least-square linear regression to artificial neural networks or decision trees [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. While each method has is specificities, strengths and weaknesses, a common issue when one wants to infer a structure-property relationship concerns the way molecules are represented. While small molecules are often represented as 2D or 3D structures in chemistry and chemoinformatics, most statistical methods, including linear models and nonlinear neural networks, require vectors an input. Molecules must therefore be first encapsulated as finite-dimensional vectors, using various molecular descriptors, before being presented as input to these algorithms. The construction of molecular descriptors is however a difficult task. Often a significant chemical expertise coupled with heuristic feature selection methods is needed to chose, among the plethora of possible molecular descriptors, the most relevant ones for a property to be predicted. The number of molecular descriptors must moreover be kept as small as possible to limit the complexity of the inference task.

An alternative to this issue has emerged recently with the advent of support vector machines (SVM) and related kernel methods in machine learning [START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF][START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF]. SVM is an algorithm for pattern recognition and regression that provides a useful framework to overcome the difficulty of data representations as vectors of low dimensions, both from a theoretical and a computational point of view. Theoretically, first, SVM are able to infer models in large or even infinite dimensions from a finite number of observations. Indeed the complexity of the learning task is not directly related to the dimension of the input vectors, but rather to some measure of complexity of the classification rules which are precisely controlled by SVM through the use of regularization [START_REF] Vapnik | Statistical Learning Theory[END_REF]. Practically, second, a computational trick known as the kernel trick allows the estimation of models with a complexity that does not depend on the dimension of the input, but only on the number of training points. Hence training a model with vectors of infinite dimension is no more computationally demanding than training a model for small fingerprints -as long as the so-called kernel function, which corresponds to the inner product of the vectors, can be computed efficiently. Combined together these properties give SVM the ability to work with molecules represented by vectors of large or even infinite dimension in a computationally efficient framework, leveraging the burden of feature selection and giving the modelers new opportunities to imagine large sets of molecular descriptors.

SVM often provide state-of-the-art performances on many classification and regression tasks, and enjoy therefore an increasing popularity in various application fields, including bioinformatics and chemoinformatics [START_REF] Schölkopf | Kernel Methods in Computational Biology[END_REF]. For example, SVM have been applied to the prediction of the activity of molecules on a number of target classes [START_REF] Burbidge | Drug design by machine learning: support vector machines for pharmaceutical data analysis[END_REF][START_REF] Weston | Feature selection and transduction for prediction of molecular bioactivity for drug design[END_REF][START_REF] Arimoto | Development of CYP3A4 inhibition models: comparisons of machine-learning techniques and molecular descriptors[END_REF][START_REF] Boser | Classifying "kinase inhibitor-likeness" by using machine-learning methods[END_REF][START_REF] Liu | QSAR and classification models of a novel series of COX-2 selective inhibitors: 1,5-diarylimidazoles based on support vector machines[END_REF][START_REF] Saeh | Lead hopping using SVM and 3D pharmacophore fingerprints[END_REF][START_REF] Tobita | A discriminant model constructed by the support vector machine method for HERG potassium channel inhibitors[END_REF], toxicological properties [START_REF] Kramer | Fragment generation and support vector machines for inducing SARs[END_REF][START_REF] Helma | Data mining and machine learning techniques for the identification of mutagenicity inducing substructures and structure activity relationships of noncongeneric compounds[END_REF][START_REF] Luan | Classification of the carcinogenicity of N-nitroso compounds based on support vector machines and linear discriminant analysis[END_REF], drug-likeliness [START_REF] Byvatov | Comparison of support vector machine and artificial neural network systems for drug/nondrug classification[END_REF][START_REF] Müller | Classifying 'drug-likeness' with Kernel-based learning methods[END_REF][START_REF] Takaoka | Development of a method for evaluating drug-likeness and ease of synthesis using a data set in which compounds are assigned scores based on chemists' intuition[END_REF], blood-brain barrier permeability [START_REF] Doniger | Predicting CNS permeability of drug molecules: comparison of neural network and support vector machine algorithms[END_REF], enantioselectivity (Aires-de Sousa and [START_REF] Vert | Prediction of enantiomeric excess in a combinatorial library of catalytic enantioselective reactions[END_REF], aqueous solubility [START_REF] Lind | Support vector machines for the estimation of aqueous solubility[END_REF], or isoelectric point (Liu et al., 2004a), to name just a few.

While most recent successful applications of SVM in chemoinformatics were obtained by just plugging classical molecular descriptors to the SVM, an increasing line of work seeks to investigate the unique opportunities offered by SVM to go beyond classical fingerprints and molecular descriptors, thanks to the kernel trick. This avenue was pioneered simultaneously and independently by [START_REF] Kashima | Marginalized Kernels between Labeled Graphs[END_REF] and [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF] who proposed to represent the 2D structure of a molecule by an infinite-dimensional vector of linear fragment counts and showed how SVM can handle this representation with the kernel trick. Later work quickly refined these 2D kernels [START_REF] Kashima | Kernels for graphs[END_REF][START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF][START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF] and proposed new infinite-dimensional representations of 3D structures [START_REF] Swamidass | Kernels for small molecules and the prediction of mutagenicity, toxicity and anti-cancer activity[END_REF]Mahé et al., 2006;[START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF].

These first attempts to enlarge the flexibility of molecular descriptor-based predictive models represent a promising direction for in silico modelling of structure-property relationship, because they illustrate the unique possibilities offered by SVM and more generally kernel methods in this context. We review them in this paper with the hope to offer a state-of-the-art description of the latest development in this field, and an invitation for the chemoinformatics community to further investigate these possibilities. For that purpose we first provide a quick introduction to SVM and kernels in Section 1, and illustrate the relevance of the kernel trick when working with 2D structures of molecules with a simple example of 2D kernel in Section 2. This example is further generalized and connected to recent work on 2D kernels in Section 3, and practical issues with these kernels are discussed in Section 4. In Section 5 we present another approach that focuses on the representation of 3D structures of molecules, and discuss practical issues for this approach in Section 6. We conclude by a discussion and suggestions for future work in Section 7.

druggable, or inhibitor/non-inhibitor of a given target. Hence a typical pattern recognition problem could be, given a list of toxic and non-toxic molecules, to learn a rule to predict whether a new candidate molecule is toxic or not.

More formally, we represent the training set available as a set of n objects x 1 , . . . , x n ∈ X , where X denotes the set of all possible objects, and associated binary labels y 1 , . . . , y n ∈ {-1; 1}. In our case each object x i represents a molecule, X denotes the set of all possible molecules, and the two classes 1 and -1 are arbitrary representations of two classes or interest, such as "toxic" and "non-toxic". Pattern recognition algorithms, such as SVM, use this training set to produce a classifier f : X → {-1; 1} that can be used to predict the class of any new data x ∈ X by the value f (x). When objects are d-dimensional vectors, that is, X = R d , the classifier output by SVM is based on the sign of a linear function:

f (x) = sign( w, x + b) , (1)
for some (w, b) ∈ X × R defined below. In this case the classifier has a geometric interpretation: the hyperplane w, x + b = 0 separates the input space X into two half-spaces, and the prediction of the class of a new point depends on its position on the one or on the other side of the hyperplane. The particular hyperplane selected by SVM is the one that solves the following optimization problem :

min w,b 1 2 ||w|| 2 + C n i=1 L (y i , w, x i + b) , (2)
where C is a parameter and L(y, t) is the hinge loss function equal to 0 if yt ≥ 1, and 1 -yt otherwise. For a given training example (x i , y i ), the hinge loss term L (y i , w, x i + b) quantifies how "good" the prediction w, x i + b of a candidate classifier (w, b) is, in the sense that the better the prediction, the smaller the loss. For example there is no loss when w and b are such that y i ( w, x i + b) ≥ 1, which means that w, x i + b has the sign of y i and is larger than 1 in absolute value. In other words the loss is zero when the prediction is correct and made with large confidence. Now the second term in the sum (2) is the average loss over the training set of the candidate classifier (w, b): it is small when the classifier fits well the training points, i.e., makes on average "good" predictions. On the other hand, the first term ||w|| 2 in (2) is small when the slope of the classifier is small. The two terms in (2) are often in conflicts, especially in large dimension, because it is often difficult to fit the training points well with linear functions of limited slopes. The rational behind the optimization problem (2) is indeed to find a linear classifier that reaches a trade-off between the goodness of fit on the training set (as quantified by the second term of this sum), and the smoothness of the classifier (as quantified by the first term).

The parameter C controls this trade-off, by balancing the importance of each term. In the extreme case when C = +∞ and the training points can be correctly separated by a hyperplane, then no error is allowed on the training set and the classifier with largest margin is found (Figure 1). It is often interesting to rewrite problem (2) in an equivalent way, using classical optimization theory. Indeed, this problem is equivalent to the following quadratic problem, called its dual:

max α∈R n    n i=1 α i - 1 4 n i,j=1 α i α j y i y j x i , x j    , subject to : n i=1 α i y i = 0 and 0 ≤ α i ≤ C , i ∈ [1 : n] . (3)
Both problems (2) and (3) are equivalent in the sense that the solution (w * , b * ) of the primal problem (2) can be deduced from the solution α * of the dual problem (3). In particular, it can be shown that w * = n i=1 α * i y i x i , and b * can also be deduced from α * . As a result, the decision function (1) can also be expressed in terms of the solution α * of the dual problem:

f (x) = sign n i=1 α * i x, x i + b * . (4)
Let us now consider the use SVM for pattern recognition with molecules, represented for example by their 2D or 3D structures. Such structures being not vectors, they can not be directly input to SVM. Instead we need to embed the set of 2D or 3D structures of molecules X to a vector space H through a mapping Φ : X → H. We can then apply the SVM algorithm to the training vectors Φ(x i ), i = 1, . . . , n, as illustrated in Figure 2. 2) is set to +∞, then the separating hyperplane selected by SVM is the one that maximizes the distance to the closest point on each side (γ on this picture). In general, some training points may be misclassified by the selected hyperplane to control overfitting.

Φ Φ

Figure 2: In order to use SVM with molecules, we need to define an embedding of the space of molecules to a vector space, i.e., a representation of each molecule x as a vector Φ(x). Note that, contrary to usual fingerprintbased approaches, the vector space might have a large or even infinite dimension.

An important point to notice is that in the dual formulation (3), the data are only present through dot-products: pairwise dot-products between the training points during the learning phase in (3), and dot-products between a new data and the training points during the prediction phase in (4). This means that instead of explicitly knowing Φ(x) for any x ∈ X , it suffices to be able to compute inner products of the form:

k(x, x ′ ) = Φ(x), Φ(x ′ ) , (5)
for any x, x ′ ∈ X . In that case the dual optimization problem (3) solved by SVM can be rewritten as follows:

max α∈R n    n i=1 α i - 1 4 n i,j=1 α i α j y i y j k(x i , x j )    , subject to : n i=1 α i y i = 0 and 0 ≤ α i ≤ C , i ∈ [1 : n] . (6)
Moreover the classification function (1) becomes:

f (x) = sign n i=1 α * i k(x, x i ) + b * . ( 7
)
Hence we see that for both the training of the SVM (6) and the prediction of the class of new points (7), the feature map Φ only appears through the function k, which is called a kernel. Importantly it is sometimes easier to compute directly the kernel k(x, x ′ ) between two points than their explicit representations as vectors in H.

In fact a classical result of [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] characterizes all functions k : X × X → R that are valid kernel, i.e., for which there exists a feature space H and a mapping Φ : X → H such that (5) holds (they constitute the so-called class of positive definite functions). Hence, with this characterization at hand, any kernel k can be used with a SVM as long as it satisfies the positive definiteness property. The formulation of SVM in terms of kernels (6-7) offers at least two major advantages over the formulation in terms of explicit vectors (3-4). First, it enables the straightforward extension of the linear SVM to non linear decision functions by using a nonlinear kernel, while keeping its nice properties intact (e.g., unicity of the solution, robustness to over-fitting, etc...). As an example, the Gaussian kernel k(x, x ′ ) = exp -||x -x ′ || 2 /2σ 2 is positive definite and can therefore be used as a kernel in the SVM algorithm (6). Plugging this kernel into (7) we see that the resulting discrimination function has the form:

f (x) = sign n i=1 α * i exp - ||x -x i || 2 2σ 2 ,
which is clearly a nonlinear function of x. Second, this formulation offers the possibility to directly apply SVM to non-vectorial data, such as 2D or 3D structures of molecules, provided a positive definite kernel function to compare these structures is defined. The definition of such structure kernels for molecules is explained in the following sections.

A simple kernel for 2D structures

It is common to describe the 2D structure of a molecule as a labeled undirected graph G = (V, E), with atoms as vertices V and covalent bounds as edges E. Here we assume that a label is assigned to each node and edge, typically to describe the type of atoms and bounds involved. In order to train linear models for structure-property relationship prediction, each labeled graph G representing a molecule must first be transformed into a vector Φ(G). In this section we describe a simple vector representation obtained by counting all walks of a given length n, and show the relevance of the kernel formulation in this case.

A walk of length n on a graph is a sequence of n adjacent vertices. We note that this definition allows a given vertex or edge to be present more than once in a walk. Clearly, the number of walks of length n on a graph G is finite, and we denote by W n (G) this set of walks in the following. By concatenating the labels of the vertices and of the edges of a walk w we obtain a sequence of labels which we denote by l(w), the label of the walk w. Moreover, we note L n the set of possible labels for walks of length n, i.e., all possible sequences alternating n vertex labels with n -1 edge labels. Figure 3 illustrates these definitions. Now a simple way to represent a graph G by a vector is to extract all walks of length n from its structure, sort them by label, and count in Φ(G) the number of walks with each possible label in L n . In other words the dimension of Φ(G) is equal to the size of L n , and for each possible walk label l ∈ L n we define the coordinate Φ l (G) as the number of walks in G having label l. More formally the feature Φ l (G) is defined by:

Φ l (G) = w∈Wn(G) 1(l(w) = l) . (8)
A direct approach to train a linear model with these vector representations would require the explicit computation and storage of Φ(G) for all graphs G in a dataset. This approach becomes problematic when n becomes large, because the number of walk labels increases exponentially with n. As an example, keeping only 6 types of atoms and 3 types of covalent bounds, the number of possible labels reaches 1, 944 for walks of length 3; 34, 992 for walks of length 4; 629, 856 for walks of length 5; and more than 3 billions for walks of length 8. This explosion in the dimension of Φ(G) suggests in practice either to restrict oneself to walks of length 2 or 3, or to compress the representation Φ(G). The later approach is widely used in chemoinformatics because fragments of length 5-10 are known to provide useful information in many structure-property relationship problems. The solution most often encountered is to use a hash table of limited size (typically 1024 or 2048) to map the vector Φ(G) onto a vector of smaller dimension, called a molecular fingerprint [START_REF] Gasteiger | Chemoinformatics : a Textbook[END_REF]. An obvious drawback of this solution is the danger of clashes, i.e., the mapping of different labels to the same position in the hashed vector.

An alternative solution for the use of large n values is to use kernels. As we now show, indeed, kernels allow the estimation of linear models for vectors Φ(G) without reducing their dimension nor requiring the computation and storage of the vectors. Indeed, remembering from Section 1 that SVM only need the definition of the inner product between vectors to estimate a linear problem, we only need to show how the inner product for the vector representation (8) can be computed efficiently. For that purpose, let us write this inner product more explicitly for any two graphs G and G ′ :

Φ(G), Φ(G ′ ) = l∈Ln Φ l (G)Φ l (G ′ ) = l∈Ln   w∈Wn(G) 1(l(w) = l)     w ′ ∈Wn(G ′ ) 1(l(w ′ ) = l)   = w∈Wn(G) w ′ ∈Wn(G ′ ) l∈Ln 1(l(w) = l)1(l(w ′ ) = l) = w∈Wn(G) w ′ ∈Wn(G ′ ) 1(l(w) = l(w ′ )) . (9)
In other words the inner product between Φ(G) and Φ(G ′ ) can be expressed exactly as the number of pairs of walks (w, w ′ ) of length n, respectively in G and G ′ , with the same label. In order to show how this number can be computed efficiently, it is useful to introduce the product graph G × G ′ which is a graph whose vertices are pairs of vertices of G and G ′ with the same label, and whose edges connect pairs of vertices which are connected both in G and G ′ (Figure 4). In other words the vertices of

G × G ′ are the pairs (v, v ′ ) ∈ V × V ′ with l(v) = l(v ′ )
, and there is an edge between (v 1 , v ′ 1 ) and (v 2 , v ′ 2 ) if and only if there is both an edge between v 1 and v 2 in G and an edge between v ′ 1 and v ′ 2 in G ′ , and if both edges have the same label. It is easy to see, then, that a walk in the product graph is a sequence of pairs of vertices (v, v ′ ), in G and G ′ , that are connected in G × G ′ and therefore in G and G ′ . Moreover both sequences of vertices in G and G ′ are made of pairs of vertices and pairs of edges with the same label, i.e., they form a pair of walks in G and G ′ with the same label. Conversely, given any walks w in G and w ′ in G ′ with same label l(w) = l(w ′ ), there is a walk in the product graph that corresponds to the pair of walks (w, w ′ ). In other words, there is a bijection between the pairs of walks in G and G ′ with the same label, on the one hand, and the walks on G × G ′ , on the other hand. Hence counting the number of pairs of walks of length n on G and G ′ with the same label is equivalent to simply counting the number of walks of length n on G × G ′ , as illustrated in Figure 4. It turns out that counting the number of walks of length n on a general graph (and in particular on a product graph for our purpose) can be easily computed by a recursion over n. Indeed, for a general graph, if we denote by A i (v) the number of walks 1'

2' 2,1' 1,2' 3,2' 1 2 3 X = Figure 4:
The product graph of two graphs (on the left) is obtained by considering all pairs of vertices with similar labels as vertices, and connecting two such vertices when the respective pairs of vertices in the initial graphs are connected (on the right). Each walk in the product graph (e.g., (2, 1 ′ ) -(3, 2 ′ )) is associated to a pair of walks in the initial graphs with same labels (e.g., 2 -3 and 1 ′ -2 ′ ), and vice-versa.

of length i starting at vertex v, then A 1 (u) = 1 for any vertex u and the following recursion formula holds:

A i+1 (v) = u∼v A i (u) , (10)
where the sum is over the neighbor vertices of v. A n (u) can therefore be computed for any u ∈ V by applying this formula recursively over i. The number of walks of length n on the graph is then simply obtained by summing A n (u) over the vertices u. We observe that if we denote by A the adjacency matrix of the graph and by 1 the vector whose entries are all equal to 1, then (10) simply expresses A i+1 1 as A × A i 1, and the count of walks of length n is equal to 1 ⊤ A n-1 1.

To summarize, we have shown that for the vector representation (8), the inner product between two graphs G and G ′ representing the 2D structures of two molecules can be computed by (i) constructing the adjacency matrix A of the product graph G × G ′ and (ii) computing 1 ⊤ A n-1 1 using the recursion (10). This computation is exact and efficient, although the dimension of the vectors can reach the billions. In particular, the complexity of the computation increases only linearly with n, while the number of features increases exponentially. Using this inner product with a kernel method for pattern recognition or regression allows to estimate a linear model in this space without ever computing nor storing any vector.

2D kernel extensions

The kernel for 2D structures presented in the previous section to illustrate the power of kernels can be used as such, but many extensions have been proposed to increase the flexibility and the expressiveness of the representation. In this section we review some of these extensions.

Walks of various lengths

In the computation of the kernel based on walks of length n, we note that kernels of length i < n are computed as intermediaries. The choice of n is arbitrary in practice and should depend on the targeted application and the data available. Alternatively we may decide not to choose a particular value of n, but to combine walks of different lengths in a joint feature model. The inner product being additive when new features are added, the kernel corresponding to the feature space ( 8) where all walks of length up to n are considered is the sum of the kernels corresponding walks of fixed length smaller than n. The complexity of the computation is barely increased for this extension: instead of performing the recursion (10) n times before summing the terms, one just need to increase a counter by the sum of the terms at each iteration.

When n increases the inner product in this "until-n" extension grows exponentially with n and diverges. A solution if one wishes to use large values for n, and even infinite n to be able to include all walks, is to weight the contributions of different walks by a factor λ(w) that will ensure convergence of the series, i.e., to consider the following kernel:

k(G, G ′ ) = ∞ n=1 w∈Wn(G) w ′ ∈Wn(G ′ ) λ(w)λ(w ′ )1(l(w) = l(w ′ )) .
(11)

As an example, [START_REF] Gärtner | Exponential and Geometric Kernels for Graphs[END_REF] proposed to weight the contribution of walks of length i in the inner product by a factor β i/2 , i.e., to consider the formula:

k(G, G ′ ) = ∞ n=1 w∈Wn(G) w ′ ∈Wn(G ′ ) β n 1(l(w) = l(w ′ )) = ∞ n=1 β n k n (G, G ′ ) , (12)
where k n denotes the kernel based on the count of walks of length exactly n. Remembering from the previous Section that

k n (G, G ′ ) is equal to 1 ⊤ A n-1 1,
where A is the adjacency matrix of G × G ′ , we can rewrite and factorize this kernel as follows for β small enough:

k(G, G ′ ) = ∞ n=1 β n 1 ⊤ A n-1 1 = β1 ⊤ ∞ n=0 β n A n 1 = β1 ⊤ (I -βA) -1 1 . ( 13
)
Hence the computation of the inner product in the infinite-dimensional space of all walk counts can be performed explicitly, at the cost of inverting the sparse matrix I -βA. In practice the first terms of the power series expansion provide a fast and good approximation to the complete kernel, and allow more flexibility in the weighting of the walks of different length.

Another weighting scheme for walks has been proposed independently by [START_REF] Kashima | Marginalized Kernels between Labeled Graphs[END_REF], who propose to define Markov random walks of each graph and weight the occurrence of each walk on a graph by its probability under the corresponding random walk model. As for the exponential decay, the random walk weighting scheme factorizes along the walks and can be computed with the same tricks as the exponential decay walk kernel.

Filtering tottering walks

In the previous section, we did not make any restriction on the definition of walks: they are simply defined as successions of connected graph vertices. Because molecular graphs are essentially undirected, this generic definition allows walks to have an erratic behaviour, which can lead in turn to a misleading information about the true structure of the graph in the kernel. Indeed, arbitrarily long walks can for instance be generated by simply alternating between two connected vertices. A natural way to increase the expressive power of walks with respect to the structure of the graphs is to prevent vertices from appearing more than once in a walk. In the terminology of graph theory, this corresponds to defining a kernel based on common paths instead of common

C C C C C C C C C
Figure 5: Tottering (red) and no-tottering (blue) walks. These two walks are labeled as a succession of 3 carbon atoms, but only the blue one involves 3 distinct atoms.

walks. Albeit very natural, this extension unfortunately renders the kernel computation untractable, as pointed out by [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF].

A computationally efficient alternative proposed by [START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF] is to disregard the tottering walks in the enumeration of walks. As illustrated in Figure 5, a tottering walk is a walk that comes back to a vertex it has just left. Although the notion of path is stronger than the notion of tottering walks for general graphs, they are equivalent on graphs without cycles. The relevance of the concept of tottering walks stems from computational advantages: as shown by [START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF] the set of tottering walks of a graph G corresponds to a set of walks of a transformed graph t(G), where the transformation t involves adding additional vertices and edges. As a result, the kernel for two graphs G and G ′ based on non-tottering walks only is easily computed as the standard walk kernel between the transformed graphs t(G) and t(G ′ ). More details about this transformation can be found in [START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF].

Increasing the expressiveness of walks

A second criticism that can be made to walk kernels is the fact that, because of their linearity, walks bear limited information about the structure of a graph. A principled way to address this issue, which is actually the topic of the next subsection, is to introduce subgraphs of a higher level of complexity in the kernel construction. In practice, however, this approach usually raises additional complexity issues that can be hard to circumvent. A simpler alternative is to keep a walk-based characterization of graphs and introduce some form of prior knowledge in the graph labeling function, in order to enrich the information brought by walks about the graph structure. This is in particular the approach taken in [START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF] where a new set of labels is defined for the vertices of a graph, based on the local environment of the atoms in the corresponding molecule. This method relies on a topological index, called the Morgan index, which is defined for each atom of the molecule according to the following iterative procedure. Initially, the index associated to every vertex is equal to 1. Then, at each step, the index of a vertex is defined as the sum of the indices associated to its neighbors at the previous iteration. This process is straightforward to implement in practice, since if we let M i be the vector of Morgan indices computed at the i-th iteration, it reads as M 0 = 1 and M i+1 = AM i , where 1 is the unity vector and A is the adjacency matrix of the graph.

As illustrated in Figure 6, Morgan indices make it possible to distinguish between atoms having the same type but different topological properties. When they are included in the labels of the vertices, these indices therefore define a walk as a sequence of atoms taken in a particular topological configuration. In practice, the advantage of this refinement is twofold. First, the introduction of topological information in walk labels enriches the information they bear with respect to the structure of the graphs to be compared. Second, because atoms are made more specific to the graph they belong to, as illustrated in Figure 6, the number of identically labeled atoms found in a pair of graphs automatically decreases, which has the effect of reducing the size of their product graph, hence the time of computing the kernel. Note that this computation advantage is surprisingly due to the increase in dimension of the feature space. We note however that while this Morgan process systematically reduces the cost of computing the kernel, performing too many iterations makes it impossible to detect common Initially, all atoms of the cycle are seen as identical. For increasing iterations, the presence of the NO2 branch is more and more reflected in the atoms of the cycle.

walks within a pair of graphs.

Subtree kernels

As mentioned in the previous subsection, the linear nature of walks limits their ability to properly encode the structure of a graph. This fact is emphasized by [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF] who show that graphs can be structurally different yet have the same walk content, which makes them indistinguishable by a kernel based on the count of common walks. Figure 7 illustrates this issue on a simple example. On the other hand, they also show that computing a perfect graph kernel, that is, a kernel mapping non-isomorphic graphs to distinct points in the feature space, is at least as hard as solving the graph isomorphism problem for which there is no known polynomial-time algorithm. This suggests that the expressiveness of graph kernels must be traded for their computational complexity. As a first step towards a refinement of the feature space used in walk-based graph kernels, [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF] introduce a kernel function comparing graphs on the basis of their common subtrees. As illustrated in Figure 8, this representation looks particularly promising for molecules, since it allows to capture in a principled way a wide range of functional features of molecules, that typically correspond to specific branching patterns on their associated graphs. On the practical side, this type of kernels can be computed by means of dynamic programming algorithms that recursively detect and extend identical neighborhood properties within the vertices of the graphs to be compared, in order to explicitly build their set of common subtrees. The relative contribution of subtrees of different sizes is typically controlled by means of a parameter playing a similar role to that of the parameter β in Equation ( 12). These algorithms have a prohibitive complexity in general, but they can be deployed for molecular graphs where, because of valence rules, the degree of the vertices is small in average. The relevance of this class of kernels, as well as its relationship with standard walk-based kernels, has been analyzed in details in Mahé and Vert (2006).

Figure 7: Two graphs having the same walk content, namely • : ×5 ; •→• : ×4 and •→•→• : ×2, and consequently mapped to the same point of the feature space corresponding a kernel based on the count of walks [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF].

2D kernels in practice

As a conclusion about kernels for 2D structures, we now discuss several issues related to their application in practice.

Implementation and complexity issues

As mentioned in Section 2 an elegant way to compute walk based kernels lies in the product graph formalism initially introduced by [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF]. The basic idea of the product graph construction is to merge the pair of graphs to be compared into a single graph, in such a way that a bijection is defined between the set of walks of the product graph and the set of common walks of the two initial graphs. It then follows that the number of walks of a given length occuring at the same time in the two graphs can be obtained by simple matrix products, which actually offers a closed form solution to the computation of kernels based on walks of infinite length for well chosen walk weighting schemes [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF][START_REF] Kashima | Kernels for graphs[END_REF]. As a result, even though the dimensionality associated to these kernels can be very large, and actually infinite, computing these kernels under the product graph formalism has a polynomial complexity with respect to the product of the size of the graphs to be compared1 . In practice, this type of product-graph implementations remains time consuming, even for relatively small graphs, which questions the suitability of these kernels for virtual screening applications, that typically involve large datasets of molecules. However, if only walks up to a given length are considered, which usually makes sense for real world applications, fast algorithms can be used to compute walk kernels, based for instance on trie tree structures and string kernel algorithms [START_REF] Leslie | The spectrum kernel: a string kernel for SVM protein classification[END_REF][START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF], or standard depth-first search procedures [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF]. Moreover, alternative implementations allowing to drastically reduce the time needed to compute such kernels in their general form have recently been proposed [START_REF] Vishwanathan | Fast Computation of Graph Kernels[END_REF].

Kernel normalization

A potential drawback of kernels comparing structured objects by means of their substructures lies in the fact that kernel values are highly dependent on the size of the objects to be compared. Indeed, big objects tend to be granted a higher degree of similarity than small objects for the only reason that they are made of a larger number of substructures. This fact can lead to a serious bias of the subsequent prediction model, and the classical way to tackle this issue it to apply a normalization operation in order to take into account the size of the objects in the value of the kernel function. In practice, the mainstream normalization scheme is given by the following expression:

k(x, y) = k(x, y) k(x, x)k(y, y) ,
where k is the original kernel, and k its normalized value. Note that this normalization operation has the effect of setting the diagonal of the kernel matrix to one, meaning that individual objects are given the same degree of self-similarity, whatever their size is. Geometrically, this amounts to scaling all vectors to unit norm before taking their inner product.

In the context of molecular graph kernels, alternative normalization schemes based on the Tanimoto similarity coefficient have recently been introduced [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF][START_REF] Swamidass | Kernels for small molecules and the prediction of mutagenicity, toxicity and anti-cancer activity[END_REF]. The Tanimoto coefficient is widely used in chemoinformatics to assess the similarity of molecular fingerprints. For a pair of fingerprints (A, B), it is defined as:

T AB = A ⊤ B A ⊤ A + B ⊤ B -A ⊤ B .
For binary fingerprints, it can be seen as the ratio between their intersection, that is, the number of bits set to one in both fingerprints, and their union. As pointed out by [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF], since it is based on inner product operations, this coefficient can be generalized to any kernel function, leading to the notion of Tanimoto kernel, defined for a kernel k as: , y) .

k(x, y) = k(x, y) k(x, x) + k(y, y) -k(x
This transformation provides an alternative way to normalize kernel functions in the sense that k(x, x) = 1 for all x. Several variations on this idea, that allow to generalize the classical Tanimoto coefficient in different ways, are proposed in [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF][START_REF] Swamidass | Kernels for small molecules and the prediction of mutagenicity, toxicity and anti-cancer activity[END_REF].

Kernel parameterization

Last but not least comes the issue of kernel parameterization. This question is of tremendous importance since a bad parameterization can seriously entail the success of the subsequent virtual screening application. First one must choose to consider the kernel based on walks or the kernel based on subtrees. As for now, this questions remains largely open since apart from the study of Mahé and Vert (2006), the relevance of subtrees in graph kernels has not been studied in details. While the preliminary results presented in this study suggest that subtree kernels may indeed improve over their walk-based counterparts, they also show that this class of substructures raises additional issues, related in particular to the computational complexity of the kernels as well as the explosion in the number of subtrees found in the graphs.

Concerning the parameterization of walk kernels, the main issue concerns the length(s) of the walks to consider: either walks of a precise length, up to a maximal (but finite) length, or even up to infinite length. In practice, this question is highly dependent on the problem considered. Optimally choosing this parameter can therefore hardly be made a priori but involves cross-validation procedures. Although focusing on walks of a precise length can be optimal in some cases2 , a safe default choice is to consider walks of length up to a limited value to be taken around 8 or 10. Actually, because kernels based on an infinite number of walks require to down-weight the contribution of walks depending on their length (as in Equation ( 12) for instance), long walks are in practice so penalized that their individual contribution is barely taken into account in the kernel. Explicitly limiting the length of the walks to be taken into account therefore makes sense in practice. Moreover, considering a finite number of walks provides a greater flexibility in the way to control their relative contribution in the kernel, and offers the practical advantage of paving the way to the deployment of computationally cheaper algorithms, as discussed in Section 4.1. A second important issue is related to kernel normalization. Although the impact of choosing the first or the second normalization scheme introduced in Section 4.2 has not been analysed in details, Tanimoto kernels led to good results in several validation studies [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF][START_REF] Swamidass | Kernels for small molecules and the prediction of mutagenicity, toxicity and anti-cancer activity[END_REF][START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF]. Finally, one may consider further refinements such as filtering tottering walks and introducing Morgan indices. As shown in [START_REF] Mahé | Graph kernels for molecular structure-activity relationship analysis with support vector machines[END_REF] while reducing their computational costs3 . Filtering tottering walks should be subject to caution however. Indeed, as shown in Mahé and Vert (2006), while this can indeed improve the models in some cases, it seems that the tottering phenomenon can also be helpful to detect similarity between structurally different compounds.

A 3D pharmacophore kernel

Motivated by the fact that the tridimensional structure of molecules have a central role in many biological mechanisms, including drug-target interactions for instance, recent attempts have been made to develop kernels for 3D structure of molecules. In this section, we introduce a class of kernels that relies on the notion of pharmacophore which is widely used in chemoinformatics. A pharmacophore is usually defined as a spatial arrangement of three to four atoms4 responsible for the biological activity of a drug molecule. In the following, we focus on three-points pharmacophores composed of three atoms, whose arrangement therefore forms a triangle in the 3D space (Figure 9), but similar ideas naturally apply to pharmacophores of different cardinalities5 . With a slight abuse we refer as pharmacophore below to any possible configuration of three atoms arranged as a triangle and present in a molecule, representing therefore a putative configuration responsible for the biological property of interest. More precisely, we consider a molecule m as a set of atoms in the 3D space, that is:

m = (x i , l i ) ∈ R 3 × L i=1,...,| m | ,
where | m | is the number of atoms that compose the molecule, and (x i , l i ) ∈ R 3 × L stands for its i-th atom, x i being its vector of (x,y,z) coordinates, and l i its label, such as its type for instance, but more generally taken from a set L of atom labels. With these notations at hand, the set of three-points pharmacophores that can be extracted from the molecule m can be formally defined as:

P(m) = (p 1 , p 2 , p 3 ) ∈ m 3 , p 1 = p 2 , p 1 = p 3 , p 2 = p 3 .
Following our discussion of Section 2, a simple way to represent a molecule m is to extract all its pharmacophores, sort them by type, and count in a vector Φ(m) the number of pharmacophores of each possible type. Clearly, the number of pharmacophores associated to a molecule is finite, but since their definition is based on the precise (x,y,z) coordinates of the atoms it is made of, or equivalently on continuous inter-atomic distances, the space of all possible pharmacophores is infinite. Defining such a vector representation therefore requires in practice to discretize the space of pharmacophores, which boils down to discretizing the range of inter-atomic distances into a pre-defined number of bins. Formally, if we consider n bins in the discretization, this operation defines a space of discrete pharmacophores T = L 3 × [1, n] 3 , where each pharmacophore corresponds to a triplet of atom labels, taken from the alphabet L, and a triplet of distance bin indices, taken in [1, n]. We can now define the vector representation Φ(m), in which each coordinate Φ t (m) is the number of pharmacophores extracted from the molecule m that correspond to the discrete pharmacophore t, that is,

Φ t (m) = p∈P(m) 1(disc(p) = t),
where the function 1(disc(p) = t) is one if the discretized version of the pharmacophore p is t, meaning that they are based on the same triplet of atom labels and the same triplet of distance bins, and zero otherwise. The number n of bins considered in the discretization specifies the resolution at which distinct pharmacophores are considered to be equivalent, and constitutes a critical parameterization issue. Indeed, small distance bins may prevent the detection of similar pharmacophores, while large distance bins can lead to a matching between unrelated pharmacophores. In practice, this parameter also defines the dimension of Φ(m). For example, considering 6 distinct types of atoms and 10 distance bins, which corresponds to bins of 2 angstroms if pairs of atoms are considered to lie within the 0-20 angstrom distance range, the cardinality of T , hence the dimension of Φ(m), is 216, 000. This number is raised up to 1, 728, 000 in order to reach a precision of 1 angstrom per interatomic distance bin. This explosion in the number of dimensions suggests again that in order to explicitly store the vector Φ(m), one should either consider a limited number of bins, thereby considering a poor resolution to characterize molecular structures, or rely on hashing algorithms to map the vector Φ(m) onto a vector of limited size, which as discussed previously, has the effect of inducing clashes between distinct pharmacophores. This representation highlights once again the benefit of using kernel functions since, following the lines of Equation 9, one can define the kernel:

k(m, m ′ ) = Φ(G) ⊤ Φ(G ′ ) = p∈P(m) p ′ ∈P(m ′ ) 1(disc(p) = disc(p ′ )), (14)
which, as will be discussed in Section 6.1, enables to map pairs of molecules and compute their inner product in feature spaces indexed by millions of pharmacophores, for a computational complexity that remains polynomial with respect to the product of their sizes. Of course, the idea of representing a molecule by means of its pharmacophoric content is not new, and the above approach bears strong similarity with well known pharmacophore fingerprint representations [START_REF] Brown | The information content of 2D and 3D structural descriptors relevant to ligandreceptor binding[END_REF][START_REF] Matter | Comparing 3D pharmacophore triplets and 2D fingerprints for selecting diverse compound subsets[END_REF][START_REF] Mcgregor | Pharmacophore fingerprinting. 1. Application to QSAR and focused library design[END_REF]. The above discussion nevertheless illustrates the interest of using kernel functions in this case, since they allow to exactly compute the inner products between very high-dimensional feature vectors without the need of computing nor storing them, which is not possible in general and comes at the price of an information loss. This is not, however, the major improvement made possible by kernel functions in this context. Indeed, as shown in Figure 10, the main drawback of this approach lies in the discretization of the pharmacophore space itself: not only the choice of the discretization step controls the precision required to match a pair of pharmacophores, but it also prevents pharmacophores falling on different sides of bins edges to be matched, although they can be very close, and actually even closer that two pharmacophores falling in the same bin. The kernel approach allows to circumvent this discretization issue by means of a simple generalization of Equation ( 14), where the binary function checking whether pairs of pharmacophores have the same discretized version or not is replaced by a general kernel between pharmacophores in order to continuously quantify their similarity. Letting k P be such a kernel, this leads to the general 3D kernel formulation:

k(m, m ′ ) = p∈P(m) p ′ ∈P(m ′ ) k P (p, p ′ ), (15)
which was introduced in Mahé et al. (2006). A meaningful kernel k P between pharmacophores should intuitively quantify at the same the similarity of the triplets of atoms the pair of pharmacophores to be compared are defined from, and the similarity of their spatial arrangement. A natural way to achieve this goal, which is at the same time compatible with the algorithm implementing the kernel (15) (see Section 6.1), consists in factorizing the kernel k P along the pairs of atoms and inter-atomic distances that define the pair of pharmacophores to be compared. Mahé et al. (2006) suggest for instance to introduce elementary kernel functions k At : L × L → R and k Dist : R × R → R comparing atoms and distances respectively, and to define the kernel k P as:

k P (p, p ′ ) = 3 i=1 k At (l i , l ′ i ) 3 i=1 k Dist (||x i -x i+1 ||, ||x ′ i -x ′ i+1 ||), (16) x1 x3 x2
Figure 10: Illustration of the discretization issue. x 1 , x 2 and x 3 correspond to pharmacophores living in a discretized bidimensional (Euclidean) space. x 1 is closer to x 3 than it is from x 2 , yet the discretization affects x 1 and x 2 to the same bin and x 3 to another bin. The kernel of Equation ( 15) allows to circumvent this issue.

where the pharmacophore p (resp. p ′ ) is defined as (l i , x i ) i=1:3 (resp. (l ′ i , x ′ i ) i=1:3 ), ||.|| denotes the Euclidean distance, and the index i + 1 is taken modulo 3. In this approach, the task of defining a kernel between 3D structures therefore boils down to defining a couple of kernels comparing atoms and inter-atomic distances. These kernels intuitively define the elementary notions of similarity involved in the pharmacophore comparison, which in turns define the overall similarity between molecules. A simple default choice for these kernels is to define the atom kernel k At as a binary kernel simply checking whether the pair of atoms to be compared have the same label or not, that is:

k At (l, l ′ ) = 1(l = l ′ ),
and to define the inter-atomic distance kernel k Dist as the following Gaussian radial basis function (RBF) kernel:

k Dist (x, y) = exp(- ||x -y|| 2 σ 2 ),
where σ is a bandwidth parameter. Under this parameterization, it is interesting to note that the continuous kernel of Equation ( 15) and its discretized counter part of Equation ( 14) share an important feature: because the atom kernel k At is binary, both kernels are based on pairs of pharmacophores defined by the same triplets of atom labels. The striking difference between the two formulations lies in the fact that in the kernel of Equation ( 15), the strength of the pharmacophore matching is continuously controlled by the parameter σ of the (RBF) kernel comparing inter-atomic distances. Choosing a small value of σ corresponds to imposing a strong constraint on the spatial similarity of pharmacophores, while a larger value of σ allows pairs of pharmacophores to be taken into account in the kernel although their spatial configurations may differ. We conclude this section by noting that the class of kernels defined by Equation ( 15) does not have an explicit inner product interpretation in general, and in particular using the above parameterization. Nevertheless, this construction is known to be valid as long as the kernel k P is a proper kernel function [START_REF] Haussler | Convolution Kernels on Discrete Structures[END_REF].

3D kernel in practice

In this Section, we discuss general considerations related to the application of 3D kernels in practice.

Implementation and complexity issues

Without going into technical details, it can be shown that the class of pharmacophore kernels introduced in Section 5 can be computed by algorithms derived from those used for the computation of 2D kernels. Indeed, while the 3D structure of a molecule was previously defined as a set of atoms in the 3D space, it can equivalently be seen as a fully connected labeled (and undirected) graph, with atoms as vertices and inter-atomic distances as edge labels. Under this representation, it is easy to see that computing the continuous kernel of Equation ( 15) can be interpreted as computing a walk kernel restricted to the walks that define cycles of length 3 on the graphs.

Moreover, provided the kernel k P factorizes along the pair of pharmacophore to be compared, which is the case of the kernel proposed in Equation ( 16), it is easy to show that this kernel can be computed by product-graph algorithms and simple matrix product operations, for a cubic complexity with respect to the product of the sizes of the molecules to be compared. While this complexity can be prohibitive for applications involving large datasets of molecules, the discretized version of the kernel can benefit from fast implementations derived, here also, from string kernel algorithms and trie-tree structures. We refer the interested reader to Mahé et al. (2006) for a detailed discussion about the implementation and the computational complexity of these kernels.

Kernel parameterization

In its discretized version, the only parameter entering the definition of the kernel is the number of bins n to discretize the inter-atomic distances. As already noted in Section 5, this parameter is of critical importance since it controls the precision up to which pharmacophores are considered to be identical or not. Unfortunately, this parameter can hardly be chosen a priori, and Mahé et al. (2006) suggest to optimize this parameter using crossvalidation procedures. In this study, when optimized over the grid {4, 6, 8, ..., 30} for a 0-20 angstrom interatomic distance range, this parameter was usually taken between 20 and 30, which suggests that the matching between a pair of pharmacophores should be subject to strong spatial constraints. On the other hand, such fine grained resolutions have the effect of increasing the impact of the discretization issue illustrated in Figure 10.

Under the parameterization proposed in Section 5, the only parameter entering the definition of the general kernel of Equation ( 15) is the bandwidth σ of the RBF kernel between inter-atomic distances. In the above study, small values of σ, which correspond to strong spatial constraints in the pharmacophore comparison, are usually selected by cross-validation procedures. While in these cases the discrete and continuous formulations of the kernel tend to coincide6 , the continuous formulation usually led to better performance in this study.

Molecule enrichment

Many mechanisms of interest in tridimensional virtual screening involve specific physicochemical properties of the molecules. In the case of drug-target interaction for instance, the molecular mechanisms responsible for the binding are known to depend on a precise 3D complementarity between the drug and the target, from both the steric and electrostatic perspectives. For this reason, standard pharmacophore based approaches, and in particular pharmacophore fingerprints, usually define pharmacophores from atoms or groups of atoms having particular properties. Typical molecular features of interest are positive and negative charges, high hydrophobicity, hydrogen donors and acceptors and aromatic rings [START_REF] Pickett | Diversity profiling and design using 3D pharmacophores : Pharmacophores-Derived Queries (PQD)[END_REF].

Similarly to the introduction of Morgan indices in 2D kernels discussed in Section 3.3, the atom-based kernel constructions presented in the previous section can naturally be extended to integrate this type of external information using specific label enrichment schemes. For instance, Mahé et al. (2006) use a simple scheme where the label of an atom is composed of its type and the sign of its partial charge. Positively-charged, neutral and negativaly-charged atoms of carbon are therefore labeled as {C + , C 0 , C -} in this approach. Alternative labeling schemes are considered in [START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF], based in particular on element hybridization, where for instance an sp3 carbon atom is labeled as C.3, and a typing of atoms according to conventional pharmacophoric features, such as polarity, hydrophobicity, and hydrogen-bond donors and acceptors. These studies show that, in general, such label enrichments have a positive influence on the subsequent structure-activity relationship models, while enabling to drastically reduce the computation cost of the kernels in some cases (Mahé et al., 2006).

Conformational analysis

For real-world applications, considering the tridimensional structure of molecules raises the additional issue of conformational analysis. Indeed, because of the presence of rotational bonds, molecules are not static in the 3D space, but can alternate between several spatial configurations of low-energy called conformations. The mainstream approach to conformational analysis is to represent a molecule as a set of structures, called conformers, sampled from its class of admissible conformations. On the methodological side, this operation casts the learning problem into the framework of multi-instance learning, that has been drawing a considerable interest in the machine learning community since its initial formulation [START_REF] Dietterich | Solving the Multiple Instance Problem with Axis-Parallel Rectangles[END_REF]. The SVM and kernel approaches lend themselves particularly well to this problem, due, on the one hand, to extensions of the SVM algorithm [START_REF] Andrews | Multiple Instance Learning with Generalized Support Vector Machines[END_REF], and, on the other hand, to the possibility to define kernels between sets of structures from a kernel between structures [START_REF] Gärtner | Exponential and Geometric Kernels for Graphs[END_REF]. A possible solution to the latter approach consists in averaging kernel values over all possible pairs of conformers. While more elaborated schemes can be adopted, such as [START_REF] Blaschko | Conformal Multi-Instance Kernels[END_REF] for instance, this simple configuration was already shown to be efficient in practice by [START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF].

Discussion

As a conclusion, it is probably fair to say that the empirical evaluations of the different kernel constructions introduced in this paper demonstrate the relevance of the approach based on structure kernels for virtual screening. Indeed, on the different tasks they have been tested on, including notably the prediction of high mutagenicity molecules and drug-target inhibitors, these kernels often compare favorably to state-of-the-art approaches. Moreover, because of the intrinsic modularity of kernel methods, this approach offers, to some extent, a unified approach to SAR and virtual screening, for two reasons. First, because they circumvent the need of selecting and extracting molecular descriptors, these kernels can straightforwardly be used to model different biological properties. Second, although we focused in this paper on classification applications, these kernels can be used in conjunction with the whole family of algorithms called kernel methods to solve a great variety of tasks which are relevant for virtual screening and chemoinformatics applications, such as, for instance, regression, clustering and similarity analysis. Concerning its practical use for the screening of large datasets however, it must be stressed that the approach based on kernel methods can be computationally demanding, even for relatively small datasets. Speeding up SVM and kernel methods for large datasets is currently a topic of interest in the machine learning community, and applications in virtual screening on large databases of molecules will certainly benefit from the advances in this field. The choice of a particular kernel, or even more importantly, of the 2D or 3D representation of molecular structures, should be dictated by the application considered. For example, while it is widely accepted that several drug-like properties, such as intestinal absorption [START_REF] Lipinski | Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings[END_REF] or mutagenicity [START_REF] King | Structure-activity relationships derived by machine learning: the use of atoms and their bond connectivities to predict mutagenicity by inductive logic programming[END_REF] for instance, can be efficiently deduced from the 2D structure of the molecule, target binding prediction is known to depend on a precise 3D complementarity between the structures of the drug and the target, from both the steric and electrostatic perspectives [START_REF] Böhm | Shortest-path kernels on graphs[END_REF]. Nevertheless, even in such problems that intrinsically depend on tridimensional mechanisms, it is not clear that models based on 3D kernels are more efficient than models based on 2D kernels. This fact is especially emphasized in [START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF] where 2D kernels are shown to outperform 3D kernels7 in general, which actually tallies previous fingerprint-based studies (Brown andMartin, 1996, 1997).

We see many potential extensions to the general kernel constructions presented above:

• First, the fact that the models could benefit from simple data enrichment schemes, based, for instance, on Morgan indices in the 2D case and partial charges in the 3D case, suggests that the introduction of a more thorough chemical knowledge could improve the expressive power of the kernels. In particular, several reduced representations of molecular structures exist, defined, for instance, by merging aromatic cycles and atoms that are part of the same functional groups in the 2D representation [START_REF] Gillet | Similarity searching using reduced graphs[END_REF], or by considering generic pharmacophoric features instead of isolated atoms in the 3D case [START_REF] Pickett | Diversity profiling and design using 3D pharmacophores : Pharmacophores-Derived Queries (PQD)[END_REF]. Applying such transformations in a pre-processing step is most likely to improve the characterization of the molecular structures in the kernels, while reducing their computational cost.

• Other important issues that, in our opinion, would be worth studying in more details are related to conformational analysis, and more precisely to the way the conformational space of a molecule is sampled and multi-instance kernels are defined. Although in their current form 3D kernels tend to be outperformed by their 2D counterparts [START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF], we believe that a proper handling of multi-conformers, together with a higher level of pharmacophoric characterization of molecules, can have a great impact for virtual screening applications.

• Another possible extension would be to adopt a global representation of molecules and to integrate the information derived from their 1D, 2D and 3D structures. A possible approach would be to consider a single kernel defined as a linear combination of kernels for 2D and 3D structures, together with a simple kernel based on global physicochemical properties. Several methods have been proposed to optimize such a kernel combination within the framework of support-vector machines based, for instance, on semidefinite programming [START_REF] Lanckriet | Kernel-based integration of genomic data using semidefinite programming[END_REF].

• Finally, in the case of drug-target prediction when additional information about the structure of the target is available, it would be interesting to combine the ligand-and the structure-based approaches to virtual screening, that would most likely benefit to each other in this context.

Last but not least, note that this gentle introduction to kernels for molecular structures and virtual screening applications only reflects our own view and experience, and was deliberately biased towards our own developments in this field. Indeed it must be stressed that, following the pioneer introduction of graph kernels by [START_REF] Kashima | Marginalized Kernels between Labeled Graphs[END_REF] and [START_REF] Gärtner | On graph kernels: hardness results and efficient alternatives[END_REF], several alternative kernel constructions have been proposed in recent years, among which:

• A graph kernel based on the detection of cyclic-and tree-patterns by [START_REF] Horváth | Cyclic pattern kernels for predictive graph mining[END_REF].

• A graph kernel based on the count of common paths by [START_REF] Böhm | Shortest-path kernels on graphs[END_REF]. However, because it is not possible to consider exhaustive sets of paths, as mentioned in Section 3.2, the kernel construction is restricted to the sets of shortest paths between pairs of vertices.

• An optimal assignment kernel, based on the idea of optimally assigning the atoms from one molecule to those of another, by [START_REF] Fröhlich | Optimal assignment kernels for attributed molecular graphs[END_REF]. This kernel formulates as the sum of a kernel between pairs of atoms, that has to be maximized over all possible assignment of the set of atoms of the smaller molecule to the set of atoms of the bigger one. Unfortunately, albeit very natural, this kernel is not positive definite and might require additional tricks to be used with kernel methods.

• Finally, borrowing techniques from computational geometry, standard walk-based graph kernels have recently been extended to kernels between tridimensional structures, based on graphs approximating molecular surfaces by [START_REF] Azencott | One-to four-dimensional kernels for virtual screening and the prediction of physical, chemical, and biological properties[END_REF].

Together with the references given in the above presentation, this list constitutes, to our knowledge, a comprehensive view of kernel for molecular structures with applications in virtual screening. As an ending remark, we would like to mention that open-source implementations of the family of kernels introduced in this paper can be found within the C++ ChemCpp toolbox, freely and publicly available at http://chemcpp.sourceforge.net. We hope that this introductory presentation, together with the availability of this software, will help and motivate the chemoinformatics community to further investigate SVMs and molecular kernels to model structure-activity relationship.

Figure 1 :

 1 Figure 1: SVM estimates a linear separation between the classes. When the training patterns are linearly separable and the trade-off parameter C in Equation (2) is set to +∞, then the separating hyperplane selected by SVM is the one that maximizes the distance to the closest point on each side (γ on this picture). In general, some training points may be misclassified by the selected hyperplane to control overfitting.

Figure 3 :

 3 Figure 3: The 2D structure of a molecule (on the left) can be represented by a labeled graph (on the right). Two walks on the graph are illustrated, together with their label and length.

Figure 6 :

 6 Figure6: Illustration of the Morgan process. Initially, all atoms of the cycle are seen as identical. For increasing iterations, the presence of the NO2 branch is more and more reflected in the atoms of the cycle.

Figure 8 :

 8 Figure8: Illustration of the tree-structured fragment representation: a graph G (left) and an extract of its feature space representation φ(G) (right). Note that the green tree corresponds to a walk structure.

Figure 9 :

 9 Figure 9: Left: the molecule of flavone. Right: a pharmacophore made of one hydrogen bond acceptor (topmost sphere) and two aromatic rings, with distances d 1 , d 2 and d 3 between the features that can be extracted from its structure, as shown in the middle.

More precisely, the worst case complexity is cubic with respect to the product of the size of the graphs.

For instance, walks of length 6 or 7 can be optimal to characterize molecules mainly made of aromatic cycles.

Actually, this is only true for product-graph implementations. For trie-tree implementations, Morgan indices have the opposite effect of increasing the cost of computing the kernel.

More generally, pharmacophore are defined as arrangements of groups of atoms having particular properties, such as positive or negative polarity, high hydrophobicity, and so on.

In particular, similar ideas were developed in[START_REF] Swamidass | Kernels for small molecules and the prediction of mutagenicity, toxicity and anti-cancer activity[END_REF] based on two-points pharmacophores, that is to say, distances between pairs of atoms.

Indeed, in the extreme case where σ tends to 0 and n to +∞, both formulations are equivalent.

including 3D kernels based on multi-conformers.

Acknowledgments

The work related to this paper was done while PM was PhD student in the Center for Computational Biology of the Ecole des Mines de Paris (France). PM would like to thank Xerox Research Center Europe for letting him spend some time writing this paper.


Florian Leprêtre
email: florian.lprt@gmail.com

Fabien Teytaud
email: fabien.teytaud@univ-littoral.fr

Julien Dehos
email: julien.dehos@univ-littoral.fr

Multi-armed bandit for stratified sampling: Application to numerical integration

Keywords: numerical integration, stratified sampling, machine learning, UCB

Contextual multi-armed bandits model decision problems, where the properties of the possible decisions are initially partially known, but may become better known as time passes. Such models have numerous applications and many algorithms have been proposed to provide approximate solutions.

In this paper, we propose an algorithm for computing multidimensional integration problems. Such problems are very common and can be solved using the Monte-Carlo method with the stratified sampling technique. This method consists in partitioning the integration domain then randomly sampling the partitions. Our algorithm considers the selection of the best partition to sample as a multi-armed bandit problem, which can be solved using the Upper Confidence Bound technique. We have experimented this approach for several integration problems and observed faster convergence rates.

I. INTRODUCTION

Computing an integral numerically is a common technique when the analytical result is not known. This oftenly occurs in some fields (computational physics, computational finance, image synthesis. . . ) where efficient numerical integration methods are needed. In one dimension, algorithms based on numerical quadrature are generally satisfactory and standard implementations (for instance QUADPACK [START_REF] Piessens | QUADPACK: A Subroutine Package for Automatic Integration[END_REF]) are available. However, the situation is less satisfactory for multidimensional integrals. Indeed, the number of function evaluations required by the quadrature methods, grows exponentially with the dimension leading to the so-called curse of dimensionality [START_REF] Lepage | A new algorithm for adaptive multidimensional integration[END_REF].

Several methods have been proposed for computing multidimensional integrals [START_REF] Hahn | Cuba -a library for multidimensional numerical integration[END_REF]. The cubature method is an extension of the quadrature method to high dimensions. It is often implemented with a deterministic algorithm using an adaptive subdivision strategy (DCUHRE algorithms [START_REF] Berntsen | An adaptive algorithm for the approximate calculation of multiple integrals[END_REF]). The sparse grid method [START_REF] Gerstner | Numerical integration using sparse grids[END_REF], [START_REF] Petras | Fast calculation of coefficients in the smolyak algorithm[END_REF] is also based on quadrature formulas but uses a linear combination of tensor products (the Smolyak algorithm [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain class of functions[END_REF]). These deterministic methods can give good results but are generally limited to a moderate number of dimensions or smooth integrands.

The Monte-Carlo method is a classic numerical integration method for high dimensions. Indeed, the method is based on function evaluations at random samples and does not depend on the dimension [START_REF] Metropolis | The monte carlo method[END_REF]. However, its convergence rate is slow so many variance-reduction techniques have been proposed to improve the method, for instance importance sampling and stratified sampling [START_REF] Fishman | Monte Carlo: Concepts, algorithms, and applications, ser. Springer Series in Operations Research[END_REF]. The idea of importance sampling is to choose a sampling distribution which encourages the interesting regions of the integrand. The sampling distribution greatly impacts the efficiency of the technique; it can be constructed iteratively, as in the VEGAS algorithm [START_REF] Lepage | A new algorithm for adaptive multidimensional integration[END_REF]. The idea of stratified sampling is to subdivide the integration domain and then to apply the Monte-Carlo method on each subdivision. Here the subdivision scheme is important; it can be constructed recursively as in the MISER algorithm [START_REF] Press | Recursive stratified sampling for multidimensional monte carlo integration[END_REF].

Finally, the Bayesian Quadrature method [START_REF] Rasmussen | Bayesian monte carlo[END_REF], [START_REF] Briol | Frankwolfe bayesian quadrature: Probabilistic integration with theoretical guarantees[END_REF] is a model-based approach for computing multidimensional integration. Here, the integrand is modeled as a Gaussian Process; a posterior distribution is inferred from prior samples and the integral is approximated as the mean of this distribution.

In this paper, we propose an improvement of the stratified sampling technique for the Monte-Carlo method. Our method uses a reinforcement learning technique to successively select the most interesting subdivision to sample (i.e. a multi-armed bandit problem). We compare our algorithm to other classic algorithms on numerous testing integrands, including the classic TESTPACK integrands [START_REF] Genz | A Package for Testing Multiple Integration Subroutines[END_REF].

II. BACKGROUND

A. Monte-Carlo integration

The Monte-Carlo method is commonly used for computing complex integration problems [START_REF] Fishman | Monte Carlo: Concepts, algorithms, and applications, ser. Springer Series in Operations Research[END_REF], [START_REF] Veach | Robust monte carlo methods for light transport simulation[END_REF]. Such problems can be defined by:

I = Ω f (x)dµ(x) (1)
where Ω is the integration domain (monodimensional or multidimensional), f : Ω → R is the function to integrate and µ is a measure function on Ω.

The Monte-Carlo method uses random sampling to estimate the value of the integral. In its simpliest form (where the domain Ω is the unit hypercube and the random samples X k are chosen uniformly and independently), the Monte-Carlo estimator can be written as:

F = 1 K K k=1 f (X k ). (2)
The estimator F converges to I at a rate of O(K -1/2 ); this can be shown by computing the expected value of F and the evolution of the standard deviation with K.

B. Stratified sampling

Stratified sampling is a variance reduction technique which can improve the convergence rate of the Monte-Carlo method. The basic idea is to partition the integration domain in strata Ω n :

Ω = N n=1 Ω n , (3)
and then to consider a given number of samples K n in each stratum Ω n , with:

K = N n=1 K n . (4)
The corresponding estimator is:

F = N n=1 v n 1 K n Kn k=1 f (X n,k ), (5)
where v n = µ(Ω n ) is the volume of the stratum Ω n . The variance of F is then:

Var(F ) = N n=1 v 2 n σ 2 n K n , (6)
where σ 2 n = Var(f (X n,k )). Assuming (for simplification) that K n = v n K, then the variance of the stratified estimator is:

Var(F ) = 1 K N n=1 v n σ 2 n . (7)
It can be shown that the variance of the unstratified estimator is [START_REF] Veach | Robust monte carlo methods for light transport simulation[END_REF]:

Var(F ) = 1 K N n=1 v n σ 2 n + N n=1 v n (µ n -I) 2 , (8)
where µ n is the mean value of f in the stratum Ω n . This implies that the variance of the stratified estimator is always less than or equal to the variance of the unstratified estimator, i.e. Var(F ) ≤ Var(F ). However, the effective reduction of variance depends on the partition of the integration domain and on the number of samples chosen in each stratum.

C. Multi-Armed Bandit

One historical paradigm for online learning is defined as the multi-armed bandit problem [START_REF] Lai | Asymptotically efficient adaptive allocation rules[END_REF]. This problem involves N bandit arms with unknown reward probabilities p n . At each step a player or a program selects an arm and receives a new reward r k where r k = 1 with probability p n and 0 otherwise. The goal is to maximize the cumulated reward gathered over all time steps K. This is equivalent to minimize the regret, defined as the loss incured compared to the best arm. The regret is defined as:

Kp * - K k=1 r k , (9)
where p * is the maximal reward probability among p 1 , . . . , p N . One recent way for dealing with this problem is the Upper Bound Confidence (UCB) algorithm [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF]. This algorithm provides an optimal asymptotic bound on the regret in O(ln(K)). At each step k, the algorithm selects an arm n which maximizes the formula:

pj,k + 2 ln N n=1 K n,k K j,k , (10)
where pj,k is the average reward for the arm j and K n,k is the number of times the arm n has already been selected. This formula is a trade-off between exploitation and exploration. The first part of the formula corresponds to the exploitation part as it tends to select the arm with the optimal average reward. The second part of the formula corresponds to the exploration part as it tends to select the arm which has been selected the most rarely.

III. PROPOSED METHOD

In this paper, we propose a multidimensional integration method, called UCBATURE, which combines the stratified sampling technique and the UCB algorithm, introduced in section II. The UCBATURE function is shown in Algorithm 1 and explained in this section.

A. Preliminary description

As the stratified sampling technique, our algorithm first divides the integration domain Ω into regions (the strata partitioning defined in equation 3). If the integrand is sufficiently "complex", some regions of the integration domain may be more interesting to sample. This can be seen as a multi-armed bandit problem: each region Ω n is an arm and we have to decide which arm to select (i.e. which region to sample).

Thus, rather than taking the same number of samples in every partitions (as stratified sampling does), we can implement an UCB strategy to distribute the samples in a cleverer way, using an online learning process. The algorithm will then tend to concentrate function evaluations in interesting regions, where the integrand varies the most or has the largest values.

B. Initializing with stratified sampling

The first step of the proposed method is the initialization of the regions (lines 1 to 10 of Algorithm 1). Any set of non-overlapping regions satisfying equation 3 can be used to partition the integration domain. To simplify implementation, we use a regular grid (i.e. all cells have the same volume), but this is not a restriction of the method. If we split the integration domain into N 0 parts for each of the D dimensions, we get N = N D 0 regions. Then, the algorithm computes K n,0 initial random samples in each region n to approximate the mean value and the variance of the regions. This can be seen as a stratified sampling step, giving a rough estimation of the relevance of every regions, in order to initialize the next step of the algorithm (online learning).

Algorithm 1: Ucbature Input: N : number of regions Ω n : region of index n v n : volume of the region Ω n K n,0 : number of initial samples for Ω n K: total number of samples f : integrand Output:

F : estimation of the integral of f 1 begin 2 {initialization using stratified sampling}

3 for n ← 1 to N do 4 for k ← 1 to K n,0 do 5 sample x in Ω n 6 y ← f (x) 7 S n ← S n + y 8 S n ← S n + y 2 9 k n ← K n,0 10 V n ← S n kn -Sn kn 2 11 {multi-armed bandit} 12 for k ← N n K n,0 to K do 13 n ← arg max n [Ucb(k, k n , V n )] 14 sample x in Ω n 15 y ← f (x) 16 S n ← S n + y 17 S n ← S n + y 2 18 k n ← k n + 1 19 V n ← S n kn -Sn kn 2 20 {final result} 21 F ← N n v n

Sn kn

In classic stratified sampling applications, it is generally recommended to compute at least 10 samples in each region [START_REF] Westfall | Post-stratified estimation: within-strata and total sample size recommendations[END_REF]. Here, the stratified sampling step only aims to initialize the following step of the algorithm so this size can be reduced (in our experiments, we used K n,0 = 2, ∀n). Moreover, if we have a total number of samples K, increasing the number of initial samples for the stratified sampling step implies decreasing the remaining number of samples for the online learning step (to K -

N n=1 K n,0 ).
With these initial samples, we can compute an initial estimation of the mean value and of the variance, in every regions. Since we have to update these values with new samples, in the online learning step, we compute the variances (without storing the samples) thanks to the König-Huygens formula [START_REF] Mood | Introduction to the Theory of Statistics[END_REF], defined by:

Var(X) = E[X 2 ] -E[X] 2 . ( 11
)
Indeed, if we store, for each region n, the sum of the integrand evaluations (S n ), the sum of the squared integrand evaluations (S n ) and the number of samples already taken in the region (k n ), then the current variance of the region is:

V n = S n k n - S n k n 2 . (12)

C. Sampling with UCB

Once regions are initialized (with their corresponding variances, sums and number of samples), we can easily apply an UCB strategy to iteratively select the most interesting region to sample (lines 11 to 19 of Algorithm 1).

To define the reward process of this multi-armed bandit problem (see equation 10), we specify the score of an arm (i.e. region) as:

Ucb(k, k n , V n ) = V n √ k n + R ln k k n , ( 13
)
where k is the total number of function evaluations, k n is the number of function evaluations for the arm n (which is roughly the number of times the arm has already been selected), V n is the variance of the respective region (interpreted as the reward of the arm), and R is an exploration parameter. The left term in equation 13 is the exploitation part: it tends to encourage regions that have a high reward (variance). The right term is the exploration part: it tends to encourage the arms that have not been selected many times. The parameter R defines the importance granted to exploration: if R has a high value the algorithm tends to sample equally in every region (similarly to the stratified sampling technique). Thus, equation 13 ensures the trade-off discussed in section II-C.

The main loop of our algorithm consists in determining the current best region to sample (according to UCB scores), computing a sample in this region, and then updating the region (variance, sums, number of samples). This process is repeated until we reach the number of function evaluations (K). In this way, we expect that the algorithm will preferably sample inside interesting regions, while not totally ignoring other regions that appear to be less attractive.

D. Computing final result

When all samples are computed, the last step consists in computing the final value of the integral. Since the algorithm is based on stratified sampling, equation 5 can be used. In Algorithm 1, this means computing:

F = N n=1 v n S n k n , ( 14
)
where v n is the volume of the region n.

IV. EXPERIMENTS

In this section, we present the integrand functions we use for evaluating the performance of the proposed method. We also detail the testing procedure and we present and discuss some of our results.

A. Integrand functions

We use the classic TESTPACK suite of multidimensional integrand functions [START_REF] Genz | A Package for Testing Multiple Integration Subroutines[END_REF]. We also use more complex functions, defined as follow:

f sphere (x) = 1 if x L2 < 1 0 otherwise , ( 15
)
f sinc (x) = 0.2 + 0.8 D D i=1 sin(πα i x i ) πα i x i (16)
and

f sinc2cos (x) = max [0 , f sinc2cos' (T(x))] ( 17
)
where T is a rotation group and

f sinc2cos' (y) = [0.5 + f sinc (0.5 × y)] × 0.5 D D i=1 cos(α i y i + β i ) × 1 + 0.1 D D i=1 cos(50 × y i ) .
In these definitions, D is the dimension of the integration domain, and α and β are two parameters that affect the shape of the integrand functions: α is a scaling vector (e.g. peak width) and β is a displacement vector (e.g. peak location). Thus, the previous equations define several families of functions and the two parameters enable us to define different functions and difficulties, from these families. Some of the test integrand functions used in our experiments are depicted Fig. 1.

B. Testing procedure

We use the classic procedure described in [START_REF] Genz | A Package for Testing Multiple Integration Subroutines[END_REF] for testing multidimensional integration algorithms. This procedure consists in applying the integration algorithm to various integrand functions using the function families described previously. To this end, we randomly choose two parameter vectors α and β and normalize α to a given difficulty (the greater the norm, the more difficult the problem is). Then, we estimate the groundtruth of the corresponding integral using a Monte-Carlo simulation performed with a very high number of samples (10 millions). Finally, we run the integration algorithms we want to compare on the integrand function and repeat the whole procedure, for another random integrand.

To evaluate the results of an algorithm, we compute the mean absolute error, defined by:

1 M M m=1 I m -Îm ( 18
)
where M is the number of random integrands (in our experiments: one thousand for each family), I m is the estimated ground-truth for the integrand i and Îm is the resulting integral returned by the algorithm.

C. Results and discussion

We compare the UCBATURE algorithm with two classic numerical integration algorithms: MONTE-CARLO and STRATIFIED-SAMPLING. For UCBATURE and STRATIFIED-SAMPLING, we partition each dimension into N 0 = 3 parts, which means the total number of regions is N = 3 D . For the exploration parameter R of UCBATURE, we tried several values and selected the best one, for each integrand family. These parameters may be tuned more finely, since the dimension, the number of regions and the bounds of the integration domain may affect the performance of the algorithms. We do not investigate deeply the tuning of the common parameters because this affect similarly all the algorithms, and our goal is to study their differences.

Our results show that UCBATURE generally converges faster than MONTE-CARLO and STRATIFIED-SAMPLING (see Table I and Fig. 2): for a given number of integrand evaluations, our algorithm has a lower average error than the two other algorithms. This benefit is clearly noticeable on single-peak integrands (e.g. gaussian and discontinuous), however it is less noticeable on more global integrands (e.g. oscillatory and sinc2cos). This can be explained by the fact that "global integrands" have a large but homogeneous variance across all regions, therefore no region is significantly more interesting than the other ones. Thus, UCBATURE mainly performs exploration hence performances are similar to STRATIFIED SAMPLING. A thinner partitioning should make the regions more heterogeneous and increase UCBATURE efficiency.

In higher dimensions, we observe a slightly slower learning process (see Table II and Fig. 3): with few integrand evaluations UCBATURE has generally greater errors than STRATIFIED-SAMPLING but with more integrand evaluations UCBATURE becomes better. This can be explained by the higher number of regions: since the algorithm privileges exploitation during the first integrand evaluations, it needs more samples to learn the interesting regions to sample.

V. CONCLUSION

In this paper, we propose an online learning approach to compute multidimensional numerical integration. Our method is based on the Monte-Carlo method with the stratified sampling technique, and uses the Upper Counfidence Bound technique to iteratively select the most interesting region to sample. Our experiments indicate that this approach brings noticeable benefits and, more generally, that machine learning methods seem promising for solving numerical integration problems.

Since we use a naive partitioning scheme, the proposed algorithm is limited to moderate dimensions [START_REF] Mitchell | Consequences of stratified sampling in graphics[END_REF]. In future works, we would investigate using an adaptive stratification scheme or a post-stratification method [START_REF] Holt | Post stratification[END_REF] to overcome this problem. Finally, we would simplify the use of the algorithm by automatically tuning the exploration parameter.

Fig. 1 :

 1 Fig.1: Some of the integrand functions used in our experiments (for a 2D integration domain).

Fig. 2 :

 2 Fig.2: Evolution of the integration error for the 2D gaussian integrand family, according to the number of integrand evaluations (R = 0.01, N 0 = 3).

Fig. 3 :

 3 Fig.3: Evolution of the integration error for the 5D gaussian integrand family, according to the number of integrand evaluations (R = 0.0125, N 0 = 3).

TABLE I :

 I Integration error of the three algorithms for the tested integrand families, in 2D.

	integrand	10 3 evaluations MONTE-CARLO STRATIFIED	UCBATURE	10 6 evaluations MONTE-CARLO STRATIFIED	UCBATURE
	oscillatory	0.070	0.068	0.052	0.0024	0.0023	0.0017
	product peak	39.89	31.83	19.04	1.27	1.04	0.70
	gaussian	0.0096	0.0086	0.0053	0.00030	0.00028	0.00019
	C 0 -continuous	0.0071	0.0062	0.0043	0.00024	0.00022	0.00016
	discontinuous	0.0039	0.0040	0.0030	0.00013	0.00012	0.00010
	sphere	0.041	0.037	0.03	0.0014	0.0011	0.0009
	sinusc	1.74	1.26	1.14	0.057	0.042	0.037
	sinc2cos	0.36	0.30	0.22	0.012	0.01	0.008

TABLE II :

 II Integration error of the three algorithms for the tested integrand families, in 5D.

	integrand	10 3 evaluations MONTE-CARLO STRATIFIED	UCBATURE	10 6 evaluations MONTE-CARLO STRATIFIED	UCBATURE
	oscillatory	0.54	0.32	0.30	0.017	0.011	0.011
	product peak	1342	1032	1058	41	34	28
	gaussian	0.020	0.017	0.019	0.00065	0.00057	0.00043
	C 0 -continuous	0.010	0.0084	0.0087	0.00034	0.00029	0.00024
	discontinuous	0.0020	0.0018	0.0013	6.40e-5	6.27e-5	5.18e-5
	sphere	0.31	0.25	0.31	0.01	0.008	0.007
	sinusc	8422	6545	8025	392	291	280
	sinc2cos	175	163	165	5.65	4.98	4.63