Julyan Arbel
email: julyan.arbel@inria.fr.

Olivier Marchal

Hien D Nguyen

On strict sub-Gaussianity, optimal proxy variance and symmetry for bounded random variables

We investigate the sub-Gaussian property for almost surely bounded random variables. If sub-Gaussianity per se is de facto ensured by the bounded support of said random variables, then exciting research avenues remain open. Among these questions is how to characterize the optimal sub-Gaussian proxy variance? Another question is how to characterize strict sub-Gaussianity, defined by a proxy variance equal to the (standard) variance? We address the questions in proposing conditions based on the study of functions variations. A particular focus is given to the relationship between strict sub-Gaussianity and symmetry of the distribution. In particular, we demonstrate that symmetry is neither sufficient nor necessary for strict sub-Gaussianity. In contrast, simple necessary conditions on the one hand, and simple sufficient conditions on the other hand, for strict sub-Gaussianity are provided. These results are illustrated via various applications to a number of bounded random variables, including Bernoulli, beta, binomial, uniform, Kumaraswamy, and triangular distributions.

X has a density, this means that the density is symmetric with respect to µ. A simple, and remarkable, equivalence holds for most of the standard bounded random variables.

Proposition 1.1. Let X be a Bernoulli, beta, binomial, Kumaraswamy or triangular random variable. Then, X is symmetric ⇐⇒ X is strictly sub-Gaussian.

. In this article, we provide proofs for the Bernoulli, binomial, Kumaraswamy and triangular distributions.

From Proposition 1.1, it may be tempting to conjecture that the equivalence holds true for any random variable having a bounded support. However, we establish that this is not the case. This was actually one of the starting points for the present work. More precisely, we shall provide a proof of the following result.

Introduction

Sub-Gaussian distributions are probability distributions that have tail probabilities that are upper bounded by Gaussian tails. More specifically, a random variable X with finite mean µ = E[X] is sub-Gaussian if there exists σ 2 > 0 such that:

E[exp(λ(X -µ))] ≤ exp λ 2 σ 2 2
, for all λ ∈ R.

(1)

The constant σ 2 is called a proxy variance and X is termed σ 2 -sub-Gaussian. For a sub-Gaussian random variable X, the smallest proxy variance is called the optimal proxy variance and is denoted σ 2 opt (X), or simply σ 2 opt . The variance always provides a lower bound on the optimal proxy variance: V[X] ≤ σ 2 opt (X). When σ 2 opt (X) = V[X], X is said to be strictly sub-Gaussian.

The sub-Gaussian property is increasingly studied and used in various fields of probability and statistics, primarily due to its intricate link with concentration inequalities (Boucheron et al., 2013, Raginsky and[START_REF] Raginsky | Concentration of measure inequalities in information theory, communications, and coding[END_REF], transportation inequalities [START_REF] Bobkov | Exponential integrability and transportation cost related to logarithmic Sobolev inequalities[END_REF]Götze, 1999, van Handel, 2014) and PAC-Bayes inequalities [START_REF] Catoni | PAC-Bayesian supervised classification: the thermodynamics of statistical learning[END_REF]. Applications include the missing mass problem [START_REF] Mcallester | Concentration inequalities for the missing mass and for histogram rule error[END_REF][START_REF] Berend | On the concentration of the missing mass[END_REF][START_REF] Ben-Hamou | Concentration inequalities in the infinite urn scheme for occupancy counts and the missing mass, with applications[END_REF], multi-armed bandit problems [START_REF] Bubeck | Regret analysis of stochastic and nonstochastic multi-armed bandit problems[END_REF] and singular values of random matrices [START_REF] Rudelson | Non-asymptotic theory of random matrices: extreme singular values[END_REF].

This paper focuses on the study of almost surely bounded random variables, where Bernoulli, beta, binomial, Kumaraswamy [START_REF] Jones | Kumaraswamy's distribution: A beta-type distribution with some tractability advantages[END_REF] or triangular [START_REF] Kotz | Beyond Beta: Other Continuous Families of Distributions with Bounded Support and Applications[END_REF] distributions are taken as standard and common examples. If sub-Gaussianity per se is de facto ensured because the support of said random variables is bounded, then exciting research avenues remain open in the area. Among these questions are (a) how to obtain the optimal sub-Gaussian proxy variance, and (b) how to characterize strict sub-Gaussianity?

Regarding question (a), we propose general conditions characterizing the optimal sub-Gaussian proxy variance, thus generalizing previous work [START_REF] Marchal | On the sub-Gaussianity of the beta and Dirichlet distributions[END_REF] that was tailored to the beta and Dirichlet distributions. Several techniques based on studying variations of functions are proposed. In illustrating our results with the Bernoulli distribution, we prove as a by-product of Proposition 4.1 the uniqueness of a global maximum of a function that was observed by [START_REF] Berend | On the concentration of the missing mass[END_REF] "as an intriguing open problem".

As for question (b), it turns out that the symmetry of the distribution plays a crucial role. By symmetry, we mean symmetry with respect to the mean µ = E[X]. That is, we say that X is symmetrically distributed if X and 2µ -X have the same distribution. Thus, if Proposition 1.2. Symmetry of X is neither (i) a sufficient condition, nor (ii) a necessary condition, for the strict sub-Gaussian property.

The proof of this result is presented in Section 3.2, where we demonstrate that (i) there exists simple symmetric mixtures of distributions (e.g., a two-components mixture of beta distribution and a three-components mixture distribution of Dirac masses) which are not strictly sub-Gaussian, and that (ii) there exists an asymmetric three-components mixture of Dirac masses which is strictly sub-Gaussian.

Before delving into detailing the strict sub-Gaussianity property in Section 3, we first investigate some conditions that characterize the optimal proxy variance σ 2 opt , in Section 2. The results of Sections 2 and 3 are then illustrated on a number of standard random variables on bounded supports, in Section 4. Technical results are presented in Appendix A.

2 Characterizations of the optimal proxy variance σ 2 opt Let X be an almost surely bounded random variable with mean µ = E[X]. Then, X is sub-Gaussian and satisfies Definition 1 for some σ 2 > 0.

An equivalent definition is that

∀ λ ∈ R : σ 2 ≥ 2 λ 2 K(λ),
where the function K, defined on R by:

K(λ) = ln E[exp(λ[X -µ])]
, corresponds to the cumulants generating function of X -µ. Thus the optimal proxy variance σ 2 opt can be defined as the supremum

σ 2 opt = sup λ∈R 2 λ 2 K(λ). ( 2
)
If X is almost surely bounded, then this supremum is attained, see Lemma A.1 for details. Note that the function h, defined on R by

h(λ) = 2 λ 2 K(λ), (3)
is continuous at λ = 0, since a standard series expansion demonstrates that:

h(λ) λ→0 = V[X] + o(1). (4)
Moreover, h may never vanish. In fact, since the logarithm function is strictly concave, Jensen's inequality implies that for any λ ∈ R,

h(λ) = 2 λ 2 ln E[e λ(X-µ) ] > 2 λ 2 E[ln e λ(X-µ) ] = 0. (5)
Equation ( 4) also explains directly why σ 2 opt ≥ V[X], since the variance is the value of the right-hand side (r.h.s.) function at λ = 0 and thus the maximum is always greater or equal to it. We therefore have the following result.

Proposition 2.1 (Characterization of σ 2 opt by h). The optimal proxy variance is given by:

σ 2 opt = max λ∈R h(λ) = max λ∈R 2 λ 2 K(λ). (6)
We may now present a necessary (but not always sufficient) system of equations for σ 2 opt . Indeed, since the maximum is achieved at a finite point, then this point must necessarily be a zero of the derivative of h, if h is differentiable (we will denote by D k the space of functions that are k times differentiable on R and by C k the space of functions that are k times differentiable on R and for which the k th derivative is continuous on R).

Thus, we obtain the following corollary.

Corollary 2.2 (Necessary condition for σ 2 opt , with respect to h). Let σ 2 opt be the optimal proxy variance, and assume that h and K are D 1 . Then there exists a finite λ 0 , such that

σ 2 opt = h(λ 0 ) and h (λ 0 ) = 0, (7)
which is equivalent to

σ 2 opt = 2 λ 2 0 K(λ 0 ) and λ 0 K (λ 0 ) = 2K(λ 0 ), (8)
using only the centered cumulants generating function K.

In practice, the previous set of equations has to used with caution, since there may be more than one solution to the second equation involving the derivative of h (or that of K), and a global maximizer is required to be picked among the stationary points, instead of a minimizer or a local maximizer. On a case-by-case basis, the following approach based on ordinary differential equations (ODEs), satisfied by h, can be used to demonstrate that it has a unique global maximum.

Proposition 2.3. If the function h is C 2 , then it is the unique solution of the ordinary differential equations:

h (λ) + 2 λ h(λ) = 2 λ 2 K (λ) with h(0) = V[X], (9)
or

h (λ) + 3 λ h (λ) = 2 λ K (λ) λ with h(0) = V[X] and h (0) = 1 3 E[(X -µ) 3 ]. ( 10
)
Proof. The result is directly obtained by differentiating h and via standard analysis theorems.

Remark. For cases such as the Bernoulli and uniform distributions, we may prove that the r.h.s. of (10) is strictly negative on R * := R \ {0}. This implies that if λ 0 is extremal (i.e., h (λ 0 ) = 0), then it satisfies h (λ 0 ) < 0 so that it is a local maximum. This implies that h has no local minimum and thus may only have one critical point which is necessarily the unique global maximum.

We conclude this section with another possible methodology for deriving a necessary and sufficient condition for σ 2 opt . To this end, the problem needs to be addressed from a different point of view, by studying the difference of the terms of Definition 1:

∆ : (σ 2 , λ) ∈ R * + × R → exp λ 2 σ 2 2 -E[exp(λ[X -µ])]. (11)
Proposition 2.4 (Characterization of σ 2 opt , with respect to ∆). If ∆ is C 1 , then the optimal proxy variance is characterized by:

λ → ∆(σ 2 opt , λ) ≥ 0 and ∃ λ 0 ∈ R, such that ∆(σ 2 opt , λ 0 ) = 0 and ∂ λ ∆(σ 2 opt , λ 0 ) = 0. ( 12
)
Proof. See Section A.2, in Appendix.

This proof technique was used by [START_REF] Marchal | On the sub-Gaussianity of the beta and Dirichlet distributions[END_REF] for obtaining the optimal proxy variance of the beta and Dirichlet distributions. However we find more convenient to use the conditions stated in Proposition 2.3 using the function h to address the issues presented in this article, except for the triangular distribution in Section 4.2 where this method is employed for a numerical evaluation of σ 2 opt . Remark. In general, we would like to remove the condition: λ → ∆(σ 2 , λ) ≥ 0 on the r.h.s. of Proposition 2.4, in order to have a simpler (and local) characterization of the optimal proxy variance, as a solution of (12). However, this is not possible, since we may not exclude that there exists a value σ 2 < σ 2 opt for which ∆(σ 2 , λ) presents a double zero λ 0 where locally it remains non-negative but at the same time a whole interval far from λ 0 where it would be strictly negative.

3 On strict sub-Gaussianity

Conditions based on the cumulants

Strict sub-Gaussianity is fulfilled when the optimal proxy variance equals the variance. In view of Equation ( 4), Proposition 2.1 can be rewritten as the following corollary in order to characterize the strict sub-Gaussianity property.

Corollary 3.1 (Corollary of Proposition 2.1). A distribution is strictly sub-Gaussian if and only if the maximum of function h, defined in (3), is attained in zero (and is automatically equal to V[X]). That is:

max λ∈R h(λ) = h(0) = V[X]. (13)
This characterization provides necessary conditions, based on cumulants, that are required for strict sub-Gaussianity to hold. Proposition 3.2 (Necessary conditions based on cumulants). If X is strictly sub-Gaussian, then the 3 rd and 4 th cumulants of X must satisfy

κ 3 = E[(X -E[X]) 3 ] = 0, and (14)
κ 4 = E[(X -E[X]) 4 ] -3 V[X] 2 ≤ 0. ( 15
)
Proof. By definition of the cumulant generating function K(λ) of X -µ,

K(λ) = ∞ i=1 κ i λ i i! , (16)
where κ i are the cumulants of X -µ. Since κ 1 = µ -µ = 0 and κ 2 = V[X], and using values for the third and fourth cumulants given in ( 14) and ( 15), we may write (locally around λ → 0):

h(λ) = V[X] + E[(X -µ) 3 ] λ 3 + E[(X -µ) 4 ] -3 V[X] 2 λ 2 12 + O(λ 3 ). ( 17
)
Therefore if E[(X -µ) 3 ] = 0, the maximum of h(λ) cannot be h(0) and thus strict sub-Gaussianity cannot be achieved. We conclude the proof by noting that if E[(X -µ) 3 ] = 0, we have the fact that λ = 0 can be a local maximum, only if

E[(X -µ) 4 ] ≤ 3 V[X] 2 .
Condition ( 14) requires that the third centered moment is zero and Condition (15) imposes a relation between the second and fourth centered moments. Note that the latter condition can be compactly formulated via an alternative condition on the kurtosis of X:

Kurt[X] = E[(X -E[X]) 4 ] E[(X -E[X]) 2 ] 2 ≤ 3.
More specifically, sub-Gaussianity requires that the random variable has kurtosis less than or equal to three, which is the kurtosis of a standard Gaussian random variable. Such distributions are referred to as platycurtic. The fourth cumulant defined in ( 15) is also termed excess kurtosis. Thus, strict sub-Gaussianity requires negative excess kurtosis.

When the above necessary conditions ( 14) and ( 15) hold, we are not able to obtain simple additional necessary conditions on the next cumulants. In particular, note that strict sub-Gaussianity does not imply symmetry (i.e., E[(X -E[X]) 2j+1 ] = 0, for any j ≥ 0), as will be discussed in the next section.

In contrast, more can be said when the distribution is symmetric. In fact, in the symmetric case, the moments of odd order are zero, and a simple sufficient condition can be readily obtained by comparing the Taylor expansions at λ = 0 of both terms of inequality (1), as stated in the following proposition.

Proposition 3.3 (Sufficient condition based on moments). If X is symmetric with respect to its mean µ = E[X], then a sufficient condition for X to be strictly sub-Gaussian can be stated in terms of all its even moments. That is, for X to be strictly sub-Gaussian, it is sufficient that

∀j ≥ 2, E[(X -µ) 2j ] (2j)! ≤ (V[X]) j 2 j j! (18)
holds.

Proof. The proof is based on series expansions at λ = 0 of both terms of inequality (1), when the proxy variance σ 2 is set to the variance V[X]. Namely:

E[exp(λX)] = ∞ j=0 E X 2j λ 2j (2j)! , and exp λ 2 V[X] 2 = ∞ j=0 (V[X]) j 2 j λ 2j j! , (19)
when compared term-by-term, leads to inequality (1), under assumption (18). Note that inequality (18) needs be checked only for j ≥ 2, as it trivially holds for j = 0, 1.

This technique was used by Marchal and Arbel (2017) (Section 2.2) for showing that a (symmetric) Beta(α, α) random variable is strictly sub-Gaussian. We also use it to address the cases of Bernoulli and binomial, and triangular distributions in Section 4.

Link with symmetry

The relationship between strict sub-Gaussianity and symmetry was discussed in the Introduction. Here, we provide a proof of Proposition 1.2, while the proof of Proposition 1.1 is deferred to Section 4.

3.2.1

Symmetry is neither a sufficient condition. . . Simple symmetric distributions which break the necessary condition of negative excess kurtosis can easily be constructed by hand. One such construction is by means of mixture of Dirac masses. First, consider the discrete random variable

X ∼ η 2 (δ -1 + δ 1 ) + (1 -η)δ 0 , (20)
which is a three-component mixture of Dirac masses at locations -1, 0 and 1, with η ∈ [0, 1]. It is symmetric, by construction, and its excess kurtosis equals

κ 4 = E[X 4 ] -3 V[X] 2 = η -3η 2 = η(1 -3η), (21)
which is strictly positive for all values η ∈ 0, 1 3 , hence X is not strictly sub-Gaussian for these values by virtue of Proposition 3.2. On the other hand when η → 1, the distribution of X degenerates to that of the so-called Rademacher random variable, which leads to the least possible excess kurtosis of -2.

Similar counter-examples to the sufficientness of symmetry can be built in the form of mixtures of two symmetric beta variables:

X ∼ η Beta(α, α) + (1 -η) Beta(β, β),
for η ∈ (0, 1) and α, β > 0. For any value of η ∈ (0, 1), values for α, β leading to positive excess kurtosis can be obtained. For instance, we may set (η, α, β) = (0.1, 1.5, 9), to obtain the excess kurtosis κ 4 ≈ 1.1 × 10 -4 .

. . . nor a necessary condition for strict sub-Gaussianity

Although most typical bounded random variables that are strictly sub-Gaussian are symmetric (see, e.g., Proposition 1.1), the symmetry of the distributions of such variables is not a necessary condition for strict sub-Gaussianity. Examples of such distributions include mixtures of Dirac masses. For example,

X ∼ 3 i=1 p i δ x i with 3 i=1 p i = 1 (22)
with (x 1 , x 2 , x 3 ) = -2, -1 2 , 5 4 and (p 1 , p 2 , p 3 ) = 1 13 , 4 7 , 32 91 . The function h for the random variable characterized by ( 22) is plotted in Figure 1b. Note that it attains its maximum in λ = 0. In (c), η varies from 0.25 (red curve, maximum not at zero) to 0.4 (blue curve, maximum at zero). In (d), we illustrate the function h with different zooming scales (around λ = 0): from [-1.5 5 , 1.5 5 ] (blue curve, maximum zoom-out) to [-1.5 -5 , 1.5 -5 ] (red curve, maximum zoom-in), with an adapted y-scale, showing that the maximum is attained in zero.

4 Results and applications to standard distributions

Bernoulli and binomial distributions

Consider a Bernoulli random variable, X ∼ Ber(µ) with µ ∈ (0, 1) and a binomial random variable, Y ∼ Bin(n, µ) which can be obtained as the sum of n independent Ber(µ) random variables, n a positive integer.

Proof of Proposition 1.1 for the Bernoulli and binomial distributions.

Starting with the Bernoulli: the third cumulant is equal to

E[(X -E[X]) 3 ] = 2µ(1 -µ) 1 2 -µ ,
thus, by virtue of Proposition 3.2, a non-degenerate Bernoulli random variable may only be strictly sub-Gaussian when µ = 1 2 . That is, when it is symmetric. Conversely, verifying the sufficient condition for the symmetric Bernoulli distribution Ber (1/2) is equivalent to assessing the condition for the Rademacher distribution instead. That is, the distribution of random variable X, where the events X = -1 and X = 1 have equal probability

P [X = -1] = P [X = 1] = 1 2 .
Since X 2 = 1, the variance of X and all of its even moments are

V[X] = E[X 2j ] = 1.
Therefore, to verify the sufficient condition of Proposition 3.3, we are required to demonstrate that (2j)! ≥ 2 j j!, for each j ≥ 2, which follows from the expansion

(2j)! = 2j × • • • × (j + 1) j × j! ≥ 2 j j!.
Thus, we have verification of the sufficient condition for the Rademacher distribution and hence the symmetric Bernoulli distribution, as a consequence.

Turning to the binomial distribution, we observe that the optimal proxy variance of a sum of i.i.d. (independent and identically distributed) variables is the sum of the optimal proxy variances. Thus, we immediately obtain the result that

σ 2 opt [Bin(n, µ)] = nσ 2 opt [Ber(µ)]. (23)
In particular, X ∼ Bin(n, µ) is strictly sub-Gaussian if and only if µ = 1 2 .

We now turn to the optimal proxy variance of a Bernoulli, which has the form

σ 2 opt = 1 2 -µ ln 1 µ -1 . ( 24
)
This fact is known via Theorem 2.1 and Theorem 3.1 of [START_REF] Buldygin | The sub-Gaussian norm of a binary random variable[END_REF]; see also the discussion in the introduction of [START_REF] Marchal | On the sub-Gaussianity of the beta and Dirichlet distributions[END_REF]. Here, we focus on a rather different approach, based on function h and Corollary 2.2, where

h(λ) = 2 λ 2 ln[µe λ + (1 -µ)] -µλ . (25)
Note that the study of the variations of h is observed by [START_REF] Berend | On the concentration of the missing mass[END_REF] (cf. their function g; Equation (2.1)). However, a formal proof that h has a single global maximum is left "as an intriguing open problem" by Berend and Kontorovich. This is stated in the next proposition, and formally proved, below. An illustration of this result is presented in Figure 2.

Proposition 4.1. If X ∼ Ber(µ), then the function

h : λ → 2 λ 2 ln(µe λ + 1 -µ) -µλ
admits a unique critical point which is a global maximum. The global maximizer is obtained at λ 0 = 2 ln 1-µ µ , which leads to the optimal proxy variance of form

σ 2 opt = h(λ 0 ) = 1 2 -µ ln 1 µ -1 .
Proof. Let us first prove that h admits a unique critical point, which is a global maximum, by using Proposition 2.3 and the remark that follows. ODEs ( 9) and ( 10) are respectively

h (λ) + 2 λ h(λ) := r(λ) = 2µ(1 -µ)(e λ -1) λ 2 (µe λ + 1 -µ) , h(0) = µ(1 -µ) (26)
and

h (λ) + 3 λ h (λ) = r (λ) + r(λ) λ = 2µ(1 -µ) λ 3 (µe λ + 1 -µ) (λ + 2µ -1)e λ -µe 2λ + 1 -µ , (27) with h (0) = µ(1-µ)(2µ-1) 3 . Let us denote g(λ) = (λ + 2µ -1)e λ -µe 2λ + 1 -µ, u = e λ > 0 and G(u) = u ln(u) + (2µ -1)u -µu 2 + 1 -µ. We have G (u) = 1
u -2µ so that G is positive on 0, 1 2µ and negative on [ 1 2µ , +∞). Since G 1 2µ = -ln(2µ) + 2µ(1 -2µ) < 0, for µ > 1 2 , we have the fact that G is always strictly negative on R + . Thus, G is a strictly decreasing function of u and hence g is also a strictly decreasing function of λ. Note that g(0) = 0, so that g is positive on R -and negative on R + . Thus, r (λ) + r (λ) λ is always strictly negative (there is a factor λ 3 in the denominator that changes sign at λ = 0). We conclude that a point λ 2 = 0, where h (λ 2 ) = 0, always satisfies h (λ 2 ) < 0 and therefore it is always a local maximum of h.

By differentiating h, we observe that the global maximizer of h is obtained as the unique solution (in λ 0 ) of the equation

2(µe λ 0 + 1 -µ) ln(µe λ 0 + 1 -µ) -µλ 0 (1 + µ)e λ 0 + 1 -µ = 0.
It is easy to verify that λ 0 = 2 ln 1-µ µ , and that this leads to the optimal proxy variance as stated. Figure 2: Illustration of the function h for Ber(µ) random variables with µ varying from 0.05 to 0.95 (gray curves), and of the particular case of Ber 1 2 (blue curve). The maximum of h is realized at zero, only in the Ber 1 2 case, which corresponds to a symmetric distribution. Maxima of h, corresponding to the optimal proxy variance, are depicted by the red curve. Varying λ on the x-axis.

Triangular distribution

We say that X ∼ Tri(a, b) is a triangular random variable on (-a, b), for any a, b > 0, if it is characterized by a density equal to

f (x; a, b) = 2 a(a+b) (x + a) if -a < x < 0, 2 b(a+b) (b -x) if 0 < x < b. (28)
See [START_REF] Kotz | Beyond Beta: Other Continuous Families of Distributions with Bounded Support and Applications[END_REF] for details and properties of such distributions. A recent review of research developments regarding the triangular distribution appears in [START_REF] Nguyen | Progress on a conjecture regarding the triangular distribution[END_REF].

Proof of Proposition 1.1 for the triangular distribution. The third cumulant is equal to

E[(X -E[X]) 3 ] = (b -a)(2a + b)(2b + a)
270 , so by virtue of Proposition 3.2, a triangular random variable may only be strictly sub-Gaussian when a = b. That is when it is symmetric.

Conversely, when the distribution is symmetric with a = b, we can easily express the moments of even order in the form

E[X 2j ] = 2a 2j (2j + 1)(2j + 2)
, so that the sufficient moment condition of Proposition 3.3 is equivalent to

E[X 2j ] (2j)! ≤ (V[X]) j 2 j j! ⇔ 2 (2j + 2)! ≤ 1 12 j j! .
In other words, the only remaining inequality is to show that

12 j j! ≤ (2j + 2)! 2 , ∀j ≥ 1.
The result is true for j ∈ {1, 2, 3, 4} by direct computation. We then make the decomposition:

(2j + 2)! 2 = (j + 1)2 j j i=1 (2i + 1) ≥ (j + 1)2 j 3 × 5 × j i=3 7 = 15 49 (j + 1)2 j 7 j ,
compute and we may write

h(λ) = 2 λ 2 ln 2 sh( λ 2 ) λ , (29)
which is a symmetric and is a C ∞ function. It remains to prove that it attains its global maximum at λ = 0.

To this end, we use Proposition 2.3 and the remark following it. ODEs ( 9) and ( 10) are respectively

h (λ) + 2 λ h(λ) = r(λ) = ch λ 2 λ 2 sh λ 2 - 2 λ 3 , with h(0) = V[X] = 1 12 , (30)
and

h (λ) + 3 λ h (λ) = r (λ) + r(λ) λ = - 2s(λ) λ 4 sh 2 ( λ 2 ) , (31)
with h(0) = 1 12 , h (λ) = 0, and s(λ) := λ 2 + λ sh(λ) -4ch(λ) + 4. We now apply the same method as for the Bernoulli case. That is, we need to show that the r.h.s. of ( 31) is negative, which amounts to showing that the function s is positive. Let us first observe that the result holds locally around λ = 0. Indeed, we have h(0) = 1 12 , h (0) = 0 and h (0) = -1 720 < 0. Then, observe that

s(λ) = λ 2 + λ sh(λ) -4 ch(λ) + 4, s (λ) = 2λ -3 sh(λ) + λ ch(λ), s (λ) = 2 -2 ch(λ) + λ sh(λ), s (3) (λ) = λ ch(λ) -sh(λ), and s (4) (λ) = λ sh(λ), (32)
from which we immediately conclude that

s(0) = 0, s (0) = 0, s (0) = 0, s (3) (0) = 0, and s (4) (0) = 0. ( 33
)
Obviously s (4) is strictly positive on R * , thus s (3) is a strictly increasing function on R. Since s (3) (0) = 0, we conclude that s (3) is strictly negative on R * -and strictly positive on R * + . Thus, s is strictly decreasing on R * -and strictly increasing on R * + . Finally, since s (0) = 0, we conclude that s is strictly positive on R * , therefore s is a strictly increasing function on R. Since s (0) = 0 then s is strictly negative on R * -and strictly positive on R * + so that s is strictly decreasing on R * -and strictly increasing on R * + . Since s(0) = 0, we conclude that s is strictly positive on R * . This proves that h has only one unique critical point, which is therefore the global maximizer.

In conclusion for X ∼ Unif([a, b]) with a < b, we have the celebrated result that

σ 2 opt = V[X] = 1 12 (b -a) 2 .

Sum of independent uniform variables

We may now consider the sum of independent (but not necessarily identically distributed) uniform random variables. Let (X 1 , . . . , X n ) be independent variables with

X i ∼ Unif([a i , b i ])
for i ∈ {1, . . . , n}, with a i < b i and denote S n = n i=1 X i . Since the family of uniform distributions is invariant under translation and multiplication by a constant, we have the standard result that

Z i = X i -a i b i -a i ∼ Unif([0, 1]).
Thus, since

X i = (b i -a i )Z i + a i , we have h X i (λ) = (b i -a i ) 2 h Unif([0,1]) ((b i -a i )λ),
and then, since the h function of a sum of independent random variables is the sum of the h functions of the variables and by independence of the variables (X i ) i≤n , we obtain

h Sn (λ) = n i=1 (b i -a i ) 2 h Unif([0,1]) ((b i -a i )λ).
The sum of the r.h.s. of the equation above is composed of functions that are all strictly increasing on R -and all strictly decreasing on R + . Thus, it too is strictly increasing on R - and strictly decreasing on R + . In particular, the global maximum is unique and obtained at λ = 0 for which we find:

σ 2 opt [S n ] = n i=1 V[X i ] = 1 12 n i=1 (b i -a i ) 2
Remark. Note in particular that the sum of two independent uniform variables is generically a trapezoid distribution, with symmetric triangular parts, or a symmetric (up to translation) triangular distribution. However the general asymmetric triangular case, considered in Section 4.2, cannot be expressed as a sum of independent uniform distributions.

Kumaraswamy distribution

Kumaraswamy distribution is characterized by the density on (0, 1):

f (x; α, β) = αβx α-1 (1 -x α ) β-1 ,
for α, β > 0, which yields the simple distribution function of form

F (x; α, β) = 1 -(1 -x α ) β .
The distribution was first studied in [START_REF] Kumaraswamy | A generalized probability density function for double-bounded random processes[END_REF] and was considered in details by [START_REF] Jones | Kumaraswamy's distribution: A beta-type distribution with some tractability advantages[END_REF].

Proof of Proposition 1.1 for the Kumaraswamy distribution. The Kumaraswamy distribution is symmetric if and only if α = β = 1 [START_REF] Jones | Kumaraswamy's distribution: A beta-type distribution with some tractability advantages[END_REF]. In this case, it reduces to the uniform distribution, which is strictly sub-Gaussian, as was proved in Section 4.3.

Conversely, let us now consider any potentially strictly sub-Gaussian Kumaraswamy distribution. It must then satisfy the necessary conditions of Proposition 3.2. The third cumulant κ 3 vanishes if and only if the parameters satisfy the relation

α = 1 β -(β -1)2 1 β . ( 34
)
In such a case, a numerical evaluation of the 4

th cumulant κ 4 = E[(X -E[X]) 4 ] -3 V[X] 2
demonstrates that it is negative, thus both necessary conditions of Proposition 3.2 hold. However, a numerical evaluation also shows that the maximizer of the h function is never located at zero (i.e., the condition of Corollary 3.1 is not satisfied), except for α = β = 1 (i.e., the uniform distribution). This is illustrated in Figure 4, where the function h is plotted for (α, β), satisfying relation (34), with β varying in the interval [10 -3 , 5]. The maximum of h is illustrated with the red curve, showing that the global maximizer always deviates from zero, except for the case of the uniform distribution (the black curve) and the degenerate symmetric Bernoulli distribution (the blue curve). This proves the necessity of symmetry and concludes the proof. . Maxima of h, corresponding to the optimal proxy variances, are depicted by the red curve. In particular, the maxima are located at zero only for the symmetric distributions (i.e., the uniform distribution and Ber 1 2 ). Varying λ on the x-axis. Log-scale on the y-axis.

Beta distribution

The optimal proxy variance for the Beta(α, β) distribution was derived in [START_REF] Marchal | On the sub-Gaussianity of the beta and Dirichlet distributions[END_REF], Theorem 2.1. In particular, this theorem states that the optimal proxy variance is equal to the variance if and only α = β. That is, if and only if the beta distribution is symmetric. This proves Proposition 1.1 for the beta distribution. The h function and the optimal proxy variance is illustrated on Figure 5.

A Technical results

A.1 A lemma regarding the supremum of h Lemma A.1. For variables with bounded support, the supremum in (2) is a maximum.

Proof. We have exp

(λ[X -µ]) ≤ exp(|λ|[M + |µ|]) ≤ exp(2|λ|M ),
where M is the maximum value of |X| (which is finite since X is almost surely bounded). Thus we obtain 2

λ 2 K(λ) ≤ 4M |λ| λ 2 = 4M |λ| → λ→∞ 0. (35)
Therefore, the supremum is not at infinity and since the function λ ∈ R → 2K(λ) λ 2 is continuous and positive, it must achieve its maximal value at finite values of λ.

A.2 Proof of Proposition 2.4

The proof of Proposition 2.4 is based on the study of the variations of the ∆ function, defined in Equation ( 11), which is the object of the next lemma.

We first observe that for any λ ∈ R * , the function σ 2 → ∆(σ 2 , λ) is strictly increasing. Moreover, at λ = 0 we have:

∆(σ 2 , λ) = σ 2 -V[X] λ 2 2 -E[(X -µ] 3 ) λ 3 3! + 3σ 4 -E[(X -µ) 4 ] λ 4 4! + O(λ 5 ). (36)
Therefore, for σ 2 > V[X], the function λ → ∆(σ 2 , λ) is strictly positive in a neighborhood of λ = 0, while for σ 2 < V[X], the function is strictly negative in a neighborhood of λ = 0. Thus we obtain the following lemma.

Lemma A.2. The variations of λ → ∆(σ 2 , λ) depend on the value of σ 2 , with respect to σ 2 opt , as follows:

1. for σ 2 > σ 2 opt , the function λ → ∆(σ 2 , λ) is strictly positive on R, 2. for σ 2 = σ 2 opt , the function λ → ∆(σ 2 opt , λ) is non-negative and there exists at least one point λ 0 ∈ R for which ∆(σ 2 opt , λ 0 ) = 0. In particular, since the function remains non negative, it implies that ∂ λ ∆(σ 2 opt , λ 0 ) = 0 and ∂ 2 λ ∆(σ 2 opt , λ 0 ) ≥ 0, and 3. for σ 2 < σ 2 opt , there exists an interval not reduced to a point on which λ → ∆(σ 2 , λ) < 0,

where the first and second derivatives of λ → ∆(σ 2 , λ) are denoted by ∂ λ ∆ and ∂ 2 λ ∆, respectively.

Proof. The proof is based on fact that σ 2 → ∆(σ 2 , λ) is strictly increasing, and the fact that (σ 2 , λ) → ∆(σ 2 , λ) is continuous.

1. Assume by contradiction that ∆ is not strictly positive. Then, since it is non-negative, there must exist at least one point λ 1 for which ∆(σ 2 , λ 1 ) = 0 with λ 1 = 0 (because for σ 2 > σ 2 opt ≥ V[X], we know that ∆ is strictly positive around λ = 0). Thus, for any σ2 < σ 2 we have ∆(σ 2 , λ 1 ) < ∆(σ 2 , λ 1 ) = 0, so that X is not σ2 -sub-Gaussian, hence σ2 < σ 2 opt and by taking the limit σ2 → σ 2 from below, we get σ 2 ≤ σ 2 opt , which is a contradiction. 2. At σ 2 = σ 2 opt , the function ∆ must vanish at least at one point λ 0 ∈ R, while remaining non-negative on R. Indeed, the function ∆ is non-negative by the sub-Gaussian definition, but if ∆ was strictly positive on R, then the continuity of ∆, relatively to (σ 2 , λ), would imply that we may lower σ 2 without ∆ vanishing. This would be in contradiction with the minimality of σ 2 opt . 3. For σ 2 < σ 2 opt , there exists at least a point λ 1 ∈ R for which ∆(σ 2 , λ 1 ) < 0. By continuity of the function λ → ∆(σ 2 , λ), this implies that there exists a neighborhood of λ 1 in which ∆ is strictly negative, which concludes the proof.

We are now ready to prove Proposition 2.4.

Proof of Proposition 2.4. Proposition A.2 indicates that if σ 2 = σ 2 opt , then ∆ ≥ 0, and there exists λ 0 ∈ R, such that ∆(σ 2 , λ 0 ) = 0 and ∂ λ ∆(σ 2 , λ 0 ) = 0. Conversely, let us assume that ∆ ≥ 0 and ∃λ 0 ∈ R, such that ∆(σ 2 , λ 0 ) and ∂ λ ∆(σ 2 , λ 0 ) = 0. Then, since ∆ ≥ 0, we have σ 2 ≥ σ 2 opt , and since for σ2 < σ 2 , ∆(σ 2 , λ 0 ) < 0, we also have σ 2 ≤ σ 2 opt . Thus σ 2 = σ 2 opt , which concludes the proof.

Figure 1 :

 1 Figure1: Illustration of the mixtures of Dirac masses, proving Proposition 1.2 described in Sections 3.2.1 (a,c) and 3.2.2 (b,d). In (c), η varies from 0.25 (red curve, maximum not at zero) to 0.4 (blue curve, maximum at zero). In (d), we illustrate the function h with different zooming scales (around λ = 0): from [-1.5 5 , 1.5 5 ] (blue curve, maximum zoom-out) to [-1.5 -5 , 1.5 -5 ] (red curve, maximum zoom-in), with an adapted y-scale, showing that the maximum is attained in zero.

Figure 4 :

 4 Figure 4: Function h for Kumaraswamy distributions for (α, β) satisfying relation (34), with β varying from 10 -3 (top gray curve) to 5 (bottom gray curve). The black curve represents the h function of the uniform distribution. The blue curve represents the h function of the symmetric Ber 12 . Maxima of h, corresponding to the optimal proxy variances, are depicted by the red curve. In particular, the maxima are located at zero only for the symmetric distributions (i.e., the uniform distribution and Ber 1 2 ). Varying λ on the x-axis. Log-scale on the y-axis.

Figure 5 :

 5 Figure 5: Illustrations of function h for Beta(a, a) (left), and Beta(a, 2 -a) (right) distributions, with a varying from 0.2 to 1.8. The black curve represents the h function of the uniform distribution. Maxima of h corresponding to the optimal proxy variance are depicted by the red curve. Varying λ on the x-axis.

Acknowledgements O.M. would like to thank Université Lyon 1, Université Jean Monnet and Institut Camille Jordan for material support. H.D.N. is funded by Australian Research Council grants: DE170101134 and DP180101192. This work was supported by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).

which yields 12 j j! ≤ (2j + 2)! 2 ⇔ j + 1 ≥ 15 49 ⇔ j ≥ 3, thus verifying the sufficient condition in the symmetric case.

For the general case, we first observe that

.

We further observe, numerically, that the difference ∆ introduced in (11) admits a unique minimum, and also that the h function (3) admits a unique (global) maximum; see Figure 3. The optimal proxy variance can be obtained numerically, by minimizing (11), as detailed in Proposition 2.4.

Uniform distribution

In this section, we prove that the uniform distribution is strictly sub-Gaussian using a similar proof as for obtaining the optimal proxy variance in the Bernoulli case (i.e., Proposition 4.1). First, we observe that after translation/dilatation, we may always reduce the problem to the case of X ∼ Unif([0, 1]). In this case, the moment generating function is straightforward to


Dexiong Chen
email: dexiong.chen@inria.fr

Laurent Jacob

Julien Mairal

Biological Sequence Modeling with Convolutional Kernel Networks

The growing number of annotated biological sequences available makes it possible to learn genotype-phenotype relationships from data with increasingly high accuracy. When large quantities of labeled samples are available for training a model, convolutional neural networks can be used to predict the phenotype of unannotated sequences with good accuracy. Unfortunately, their performance with medium-or small-scale datasets is mitigated, which requires inventing new data-efficient approaches. In this paper, we introduce a hybrid approach between convolutional neural networks and kernel methods to model biological sequences. Our method enjoys the ability of convolutional neural networks to learn data representations that are adapted to a specific task, while the kernel point of view yields algorithms that perform significantly better when the amount of training data is small. We illustrate these advantages for transcription factor binding prediction and protein homology detection, and we demonstrate that our model is also simple to interpret, which is crucial for discovering predictive motifs in sequences.

Introduction

Understanding the relationship between biological sequences and the associated phenotypes is a fundamental problem in molecular biology. Accordingly, machine learning techniques have been developed to exploit the growing number of phenotypic sequences in automatic annotation tools. Typical applications include classifying protein domains into superfamilies [START_REF] Leslie | Mismatch String Kernels for SVM Protein Classification[END_REF][START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF], predicting whether a DNA or RNA sequence binds to a protein [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF], its splicing outcome [START_REF] Jha | Integrative deep models for alternative splicing[END_REF], or its chromatin accessibility [START_REF] David R Kelley | Basset: learning the regulatory code of the accessible genome with deep convolutional neural networks[END_REF], predicting the resistance of a bacterial strain to a drug [START_REF] Drouin | Predictive computational phenotyping and biomarker discovery using reference-free genome comparisons[END_REF], or denoising a ChIP-seq signal [START_REF] Wei Koh | Denoising genome-wide histone chip-seq with convolutional neural networks[END_REF].

Choosing how to represent biological sequences is a critical part of methods that predict phenotypes from genotypes. Kernel-based methods [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF] have often been used for this task. Biological sequences are represented by a large set of descriptors, constructed for instance by Fisher score [START_REF] Jaakkola | A discriminative framework for detecting remote protein homologies[END_REF], k-mer spectrum up to some mismatches [START_REF] Leslie | Mismatch String Kernels for SVM Protein Classification[END_REF], or local alignment score [START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF]. By using the socalled kernel trick, these huge-dimensional descriptors never need to be explicitly computed as long as the inner-products between pairs of such vectors can be efficiently computed. A major limitation of traditional kernel methods is their use of fixed representations of data, as opposed to optimizing representations for a specific task. Another issue is their poor scalability since they require computing a n × n Gram matrix where n is the number of data points.

By contrast, methods based on convolutional neural networks (CNN) are more scalable and are able to optimize data representations for a specific prediction problem [START_REF] Lecun | Backpropagation applied to handwritten zip code recognition[END_REF].

Even though their predictive performance was first demonstrated for two-dimensional images, they have been recently successfully adopted for DNA sequence modeling [Alipanahi et al., 2015, Zhou and[START_REF] Zhou | Predicting effects of noncoding variants with deep learning-based sequence model[END_REF]. When sufficient annotated data is available, they can lead to good prediction accuracy, though they still suffer from some known limitations. An important one is their lack of interpretability: the set of functions described by the network is only characterized by its algorithmic construction, which makes both the subsequent analysis and interpretation difficult. CNNs for DNA sequences typically involve much fewer layers than CNNs for images, and lend themselves to some level of interpretation [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF][START_REF] Lanchantin | Deep motif dashboard: Visualizing and understanding genomic sequences using deep neural networks[END_REF], Shrikumar et al., 2017a]. However, a systematic approach is still lacking as existing methods rely on specific sequences to interpret trained filters [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF], Shrikumar et al., 2017a] or output a single feature per class [Lanchantin et al., 2017, (3.3)]. Correctly regularizing neural networks to avoid overfitting is another open issue and involves various heuristics such as dropout [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF], weight decay [START_REF] José | Comparing biases for minimal network construction with backpropagation[END_REF], and early stopping. Finally, training neural networks generally requires large amounts of labeled data. When few training samples are available, training CNNs is challenging, motivating us for proposing a more data-efficient approach.

In this paper we introduce CKN-seq, a strategy combining kernel methods and deep neural networks for sequence modeling, by adapting the convolutional kernel network (CKN) model originally developed for image data [START_REF] Mairal | End-to-end kernel learning with supervised convolutional kernel networks[END_REF]. CKN-seq relies on a continuous relaxation of the mismatch kernel [START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF]. The relaxation makes it possible to learn the kernel from data, and we provide an unsupervised and a supervised algorithm to do sothe latter being a special case of CNNs. On the datasets we consider, both approaches show better performance than DeepBind, another existing CNN [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF], especially when the amount of training data is small. On the other hand, the supervised algorithm produces task-specific and small-dimensional sequence representations while the unsupervised version dominates all other methods on small-scale problems but leads to higher dimensional representations. Consequently, we introduce a hybrid approach which enjoys the benefits of both supervised and unsupervised variants, namely the ability of learning low-dimensional models with good prediction performance in all data size regimes. Finally, the kernel point of view of our method provides us simple ways to visualize and interpret our models, and obtain sequence logos.

We investigate the performance of CKN-seq on a transcription factor binding prediction task as well as on a protein remote homology detection. We provide a free implementation of CKN-seq for learning from biological sequences, which can easily be adapted to other sequence prediction tasks.

Method

In this section, we introduce our approach to learning sequence representations. We first review CNNs and kernel methods over which our convolutional kernel network is built. Then, we present the construction of CKN followed by the learning method. We finish the section with discussions on the interpretation and visualization of a trained CKN.

Supervised learning problem

Let us consider n sequence samples x 1 , x 2 , . . . , x n in a set X of variable-length biological sequences. The sequences are assumed to be over an alphabet A. Each sequence x i is associated to a measurement y i in Y denoting some biological property of the sequence. For instance, Y may be binary labels {-1, 1} (e.g., whether the sequence is bound by a particular transcription factor or not) or R for continuous traits (e.g., the expression of a gene). The goal of supervised learning is to use these n examples {x i , y i } i=1,...,n to learn a function f : X → Y which accurately predicts the label of a new, unobserved sequence. Learning is typically achieved by minimizing the following objective:

min f ∈F 1 n n i=1 L(y i , f (x i )) + λΩ(f ), ( 1
)
where L is a loss function measuring how well the prediction f (x i ) fits the true label y i , and Ω measures the smoothness of f . F is a set of candidate functions over which the optimization is performed. Both CNNs and kernel methods can be thought of as manners to design this set.

Convolutional neural networks.

In neural networks, the functions in F perform a sequence of linear and nonlinear operations that are interleaved in a multilayer fashion. Specifically, the CNN DeepBind [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF] represents the four DNA characters respectively as the vectors (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), (0, 0, 0, 1), such that an input sequence x of length m is represented as a 4×m matrix. DeepBind then produces an intermediate representation obtained by one-dimensional convolution of the full sequence x with p convolution filters, followed by a pointwise non-linear function and a max pooling operation along each sequence, yielding a representation x in R p of the sequence. A final linear prediction layer is applied to x. The optimization in (1) acts on both the weights of this linear function and the convolution filters. Therefore, DeepBind simultaneously learns a representation x and a linear prediction function over this representation.

DeepBind additionally modifies the objective function (1) to enforce an invariance to reverse complementation of x. The loss term is replaced with L (y i , max (f (x i ), f (x i ))) where x denotes the reverse complement of x. Using this formulation is reported by [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF] to improve the prediction performance. Other versions have been then considered, by using a fully connected layer that allows mixing information from the two DNA strands [Shrikumar et al., 2017b], or by considering several hidden layers instead of a single one [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF]. Overall, across several versions, the performance of DeepBind with a single hidden layer turned out to be the best on average on ChIP-seq experiments from ENCODE [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF].

Kernel methods.

Like in CNNs, the main principle of kernel methods is to implicitly map each training point x i to a feature space in which simpler predictive functions are applied. For kernel methods, these feature spaces are generally high-(or even infinite-) dimensional vector spaces. This is achieved indirectly, by defining a kernel function K : X × X → R which acts as a similarity measure between input data. When the kernel function is symmetric and positive definite, a classical result [see [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF] states that there exists a Hilbert space F of functions from X to R, called reproducing kernel Hilbert space (RKHS), along with a mapping ϕ : X → F, such that ϕ(x), ϕ(x ′ ) F = K(x, x ′ ) for all (x, x ′ ) in X 2 , where ., . F is the Hilbertian inner-product associated with F. In other words, there exists a mapping of sequences into a Hilbert space, such that the kernel value between any sequence pairs is equal to the inner-product between their maps in the Hilbert space. Besides, any function f in F may be interpreted as a linear form f (x) = ϕ(x), f F for all x in X . A large number of kernels have been specifically designed for biological sequences [see Ben-Hur et al., 2008, and references therein].

In the context of supervised learning (1), training points x i can be mapped into ϕ(x i ) in F, and we look for a prediction function f in F. Interestingly, regularization is also convenient in the context of kernel methods, which is crucial for learning when few labeled samples are available. By choosing the regularization function Ω(f ) = f 2 F , it is indeed possible to control the regularity of the prediction function f : for any two points x, x ′ in X , the variations of the predictions are bounded by |f (x) -f (x ′ )| ≤ f F ϕ(x) -ϕ(x ′ ) F . Hence, a small norm f F implies that f (x) will be close to f (x ′ ) whenever x and x ′ are close to each other according to the geometry induced by the kernel.

Kernel methods have several assets: (i) they are generic and can be directly applied to any type of datae.g., sequences or graphs -as long as a relevant positive definite kernel is available; (ii) they are easy to regularize. However, as alluded earlier, naive implementations lack scalability. A typical workaround is the Nyström approximation [START_REF] Christopher | Using the nyström method to speed up kernel machines[END_REF], which builds an explicit q-dimensional mapping ψ : X → R q for a reasonably small q approximating the kernel, i.e., such that ψ(x), ψ(x ′ ) R q ≃ K(x, x ′ ). Then, solving the regularized problem (1) under this approximation amounts to learning a linear model with q dimensions. We will discuss how CKNs circumvent the scalability problem, while being capable to produce task-adapted data representations.

Convolutional kernel networks for sequences

We introduce convolutional kernel networks for sequences, and show their link with mismatch kernels [START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF].

Convolutional kernel for sequences

Given two sequences x and x ′ of respective lengths m and m ′ , we consider a window size k, and we define the following kernel, which compares pairwise subsequences of length k (k-mers) within x and x ′ :

K(x, x ′ ) = 1 mm ′ m i=1 m ′ j=1 K 0 (P i (x), P j (x ′ )), (2)
where P i (x) is a k-mer of x centered at position i, represented as a one-hot encoded vector of size p = |A|k and K 0 is a positive definite kernel used to compare k-mers. 1 We follow [START_REF] Mairal | End-to-end kernel learning with supervised convolutional kernel networks[END_REF] and use a homogeneous dot-product kernel such that for two vectors z and z ′ in R p ,

K 0 (z, z ′ ) = z z ′ κ z z , z ′ z ′ , ( 3
)
and κ : u → e 1 σ 2 (u-1) . Note that when z and z ′ are one-hot encoded vectors of subsequences, K 0 (z, z ′ ) = ke -1 2σ 2 k z-z ′ 2 (more details can be found in Supplementary Section A), and we recover a Gaussian kernel that involves the Hamming distance zz ′ 2 /2 between the two subsequences. Up to the normalization factors, this choice leads to the same kernel used by [START_REF] Morrow | Convolutional kitchen sinks for transcription factor binding site prediction[END_REF]. Yet, the algorithms we will present next are significantly different. While [START_REF] Morrow | Convolutional kitchen sinks for transcription factor binding site prediction[END_REF] use random features [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] to find a finite-dimensional mapping ψ : X → R q that approximates the kernel map, our approach relies on the Nyström approximation [START_REF] Christopher | Using the nyström method to speed up kernel machines[END_REF]. A major advantage of the Nyström method is that it may be extended to produce lower-dimensional task-dependent mappings [START_REF] Mairal | End-to-end kernel learning with supervised convolutional kernel networks[END_REF] and it admits a model interpretation in terms of sequence logos (see Section 3).

Learning sequence representation

The positive definite kernel K 0 defined in (3) implicitly defines a reproducing kernel Hilbert space F over k-mers, along with a mapping ϕ 0 : X → F. The convolutional kernel network model uses the Nyström method to approximate any point in F onto its projection on a finitedimensional subspace E defined as the span of some anchor points

E = Span (ϕ 0 (z 1 ), . . . , ϕ 0 (z p )) ,
where the z i 's are the anchor points in R |A|k . Subsequently, it is possible to define a coordinate system in E such that the orthogonal projection of ϕ 0 (z) onto E may be represented by a

x ∈ X x(u) ∈ A P i (x) k-mer ψ 0 (P i (x)) ∈ F kernel mapping approximation ψ 0 (P i (x)) = K -1 2 ZZ K Z (P i (x)) global pooling ψ(x) = 1 m ψ 0 (P i (x)) ψ(x) ∈ F y prediction layer w, ψ(x) P i (x) k-mer ψ 0 (P i (x)) ∈ F 0 x 1 ∈ F 0 pooling x 1 (w) ∈ F 0 P i (x 1 ) ψ 1 (P i (x 1 )) ∈ F 1 ψ(x) ∈ F prediction layer
Figure 1: Construction of single-layer (left) and multilayer (right) CKN-seq. For a single-layer model, each k-mer P i (x) is mapped to ϕ 0 (P i (x)) in F and projected to Πϕ 0 (P i (x)) parametrized by ψ 0 (P i (x)). Then, the final finite-dimensional sequence is obtained by the global pooling, ψ(x) = 1 m m i=0 ψ 0 (P i (x)). The multilayer construction is similar, but relies on intermediate maps, obtained by local pooling, see main text for details.

p-dimensional vector ψ 0 (z). Assume for now that the anchor points z i are given. Then, a finite-dimensional embedding [see Mairal, 2016, for details] is given by

ψ 0 (z) := K -1 2 ZZ K Z (z),
where K -1 2 ZZ is the inverse (or pseudo inverse) square root of the p × p Gram matrix [K 0 (z i , z j )] ij and K Z (z) = (K 0 (z 1 , z), . . . , K 0 (z p , z)) ⊤ . It is indeed possible to show that this vector preserves the Hilbertian inner-product in F after projection: Πϕ 0 (z), Πϕ 0 (z ′ ) F = ψ 0 (z), ψ 0 (z ′ ) R p for any z, z ′ in R |A|k , where Π denotes the orthogonal projection onto E. Assuming P i (x) and P j (x ′ ) map close enough to E, a reasonable approximation is therefore K 0 (P i (x), P j (x ′ )) ≈ ψ 0 (P i (x)), ψ 0 (P j (x ′ )) R p for all i, j in (2), and then

K(x, x ′ ) ≈ ψ(x), ψ(x ′ ) R p with ψ(x) = 1 m m i=1 ψ 0 (P i (x)).
Finally, the original optimization problem (1) can be approximated by min

w∈R p 1 n n i=1 L(y i , w, ψ(x i ) ) + λ w 2 . ( 4
)
We have assumed so far that the anchor points z i , i = 1 . . . , p were giveni.e., that the sequence representation ψ(x) was fixed in advance. We now present two methods to learn this representation. The overall approximation scheme is illustrated in the left panel of Figure 1.

Unsupervised learning of the anchor points. The first strategy consists in running a clustering algorithm such as K-means in order to find p centroids z i in R |A|k that "span" well the data. This is achieved by extracting a large number of k-mers from the training sequences and by clustering them. The method is simple, performs well in practice as shown in Section 3, and can also be used to initialize the training of the following supervised variant. However, the main drawback is that it generally requires a large number of anchor points (see Section 3) to achieve good prediction, which can be problematic for model interpretation.

Supervised learning of the anchor points.

The other strategy consists in jointly optimizing (4) with respect to the vector w in R p and to the anchor points that parametrize the representation ψ.

In practice, we adopt an optimization scheme that alternates between two steps: (a) we fix the anchor points (z i ) i=1,...,p , compute the finite-dimensional representations ψ(x 1 ), . . . , ψ(x n ) of all data points, and minimize function (4) with respect to w, which is convex if L is convex; (b) We fix w and update all the (z i ) i=1,...,p using one pass of a projected stochastic gradient descent (SGD) algorithm while fixing w, at a similar computational cost per iteration as a classical CNN. The optimization for the reverse-complement formulation can be done in the same way except that it is no more convex with respect to w, but we can still apply a fast optimization method such as L-BFGS [START_REF] Liu | On the limited memory bfgs method for large scale optimization[END_REF]]. We find this alternating scheme more efficient and more stable than using an SGD algorithm jointly on w and the anchor points.

Multilayer construction

We have presented CKNs with a single layer for simplicity, but the extension to multiple layers is straightforward. Instead of reducing the intermediate representation in the left panel of Figure 1 to a single point, the pooling operation may simply reduce the sequence length by a constant factor (right panel of Figure 1), in a similar way as pooling reduces image resolution in CNN. This leads to an intermediate sequence representation x 1 and we can define a valid kernel K 1 , the same as K 0 in (3), but on subsequences of x 1 . Then the same approximation described in Section 2.2.2 can be applied to K 1 . In this way, the previous process can be repeated and stacked several times, by defining a sequence of kernels K 1 , K 2 , . . . on subsequences from the previous respective layer representations, along with Hilbert spaces F 1 , F 2 , . . . and mapping functions ϕ 1 , ϕ 2 , . . . [see [START_REF] Mairal | End-to-end kernel learning with supervised convolutional kernel networks[END_REF]. Going up in the hierarchy, each point would carries information from a larger sequence neighborhood with more invariance due to the effect of pooling layers [START_REF] Bietti | Invariance and stability of deep convolutional representations[END_REF]. The training strategy is the same as for single-layer models.

Multilayer networks can potentially model larger motifs, with larger receptive fields, and possibly discover more interesting nonlinear relations between input variables than single-layer models. However, for the transcription factor binding prediction task under the setting of DeepBind or [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF], we have observed that increasing the number of convolutional layers for CKN-seq did not improve the predictive performance (Supplementary Figure 8), as also observed by [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] for CNNs. The use of multiple layers may be however important when processing very long sequences, as observed for instance by [START_REF] David R Kelley | Sequential regulatory activity prediction across chromosomes with convolutional neural networks[END_REF], who also use dilated convolutions to model even larger receptive fields than what regular CNNs can achieve.

Difference between supervised CKNs and CNNs

The main differences between CKN and CNN models are the choice of activation function (we used an exponential function in our experiments: κ(x) = e α(x-1) ) and the transformation by the inverse square root of the Gram matrix. From a kernel point of view, the inverse square root of the Gram matrix allows us to interpret the operation as a projection onto a finite-dimensional subspace of an RKHS. From a neural network point of view, this operation decorrelates the channel entries. This can be observed when using a linear activation function κ(u) = u. In such a case, the approximated mapping is then

ψ 0 (x) = (Z ⊤ Z) -1 2 Z ⊤ x = Z⊤ x, where Z⊤ = (Z ⊤ Z) -1
2 Z ⊤ is an orthogonal matrix. Encouraging orthogonality of the filters has been shown useful to regularize deep networks [START_REF] Cisse | Parseval networks: Improving robustness to adversarial examples[END_REF], and may provide intuition why our models perform better when small amounts of labeled data are available.

Data-augmented and hybrid CKN

As shown in our experiments, the unsupervised variant is sometimes more effective than the supervised one when there are only few training samples. In this section, we present a hybrid approach that can achieve similar performance as the unsupervised variant, while keeping a low-dimensional sequence representation that is easier to interpret. Before introducing this approach, we first present a classical data augmentation method for sequences, which consists in artificially generating additional training data, by perturbing the existing training samples. Formally, we consider random perturbations δ, such that given a sequence represented by a onehot encoded vector x, we denote by x + δ the one-hot encoding vector of a perturbed sequence obtained by randomly changing some characters. Each character is switched to a different one, randomly chosen from the alphabet, with some probability p. With such a data augmentation strategy, the objective (1) then becomes

min f ∈F 1 n n i=1 E δ∼∆ [L(y i , f (x i + δ))] + λΩ(f ), ( 5
)
where ∆ is a probability distribution of the variables δ corresponding to the perturbation process described above. The main assumption is that a perturbed sequence x i + δ should have the same phenotype y i when the perturbation δ is small enough. Whereas such an assumption may not be justified in general from a biological point of view, it led to significant improvements in terms of predictive accuracy. One possible explanation may be that for the tasks we consider, determining sequences may be short compared to the entire sequence: changing a few uniformly sample positions is therefore unlikely to perturb key bases.

As we show in Section 3, data-augmented CKN performs significantly better than its unaugmented counterpart when the amount of data is small. Yet, the unsupervised variant of CKN appears to be easier to regularize, and sometimes outperform all other approaches in such a lowdata regime. This observation motivates us to introduce the following hybrid variant. In a first step, we learn a prediction function f u based on the unsupervised variant of CKN, which leads to a high-dimensional sequence representation with good predictive performance. Then, we learn a low-dimensional model f s , whose purpose is to mimic the prediction of f u , by minimizing the cost function

min f ∈F 1 n n i=1 E δ∼∆ [L( ŷi (x i + δ), f (x i + δ))] + λΩ(f ), (6)
where ŷi (

x i + δ) = y i if δ = 0 and f u (x i + δ) otherwise.
Typically, the amount of perturbation that formulation (6) can afford is much larger than (5), as shown in our experiments, since it does not require to make the assumption that the sequence x i + δ should have exactly label y i , which is a wrong assumption when δ is large.

Model interpretation and visualization

As observed by [START_REF] Morrow | Convolutional kitchen sinks for transcription factor binding site prediction[END_REF], the mismatch kernel [START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF] for modeling sequences may be written as Eq. ( 2) when replacing K 0 with a discrete function I 0 that assesses whether the two k-mers are identical up to some mismatches. Thus, the convolutional kernel (2) can be viewed as a continuous relaxation of the mismatch kernel. Such a relaxation allows us to characterize the approximated convolutional kernel by the learned anchor points (the variables z 1 , . . . , z p in Section 2.2.2) that can be written as matrices in R |A|×k .

To transform these optimized anchor points z i into position weight matrices (PWMs) which can then be visualized as sequence logos, we identify the closest PWM to each z i : the kernel K 0 implicitly defines a distance between one-hot-encoded sequences of length k, which is approximated by the Euclidean norm after mapping with ψ 0 . Given an anchor point z i , the closest PWM µ according to the geometry induced by the kernel is therefore obtained by solving

min µ∈M ψ 0 (µ) -ψ 0 (z i ) 2 ,
where M is the set of matrices in R A×k whose columns sum to one. This projection problem can be solved using a projected gradient descent algorithm. The simplicial constraints induce some sparsity to the resulting PWM, yielding more informative logos. As opposed to the approach of [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF] which has relied on extracting k-mers sufficiently close to the filters in a validation set of sequences, the results obtained by our method do not depend on a particular dataset.

Application

We now study the effectiveness of CKN-seq on a transcription factor (TF) binding prediction and a protein homology detection problem.

Prediction of transcription factor binding sites

The problem of predicting TF binding sites has been extensively studied in the recent years with the continuously growing number of TF-binding datasets. This problem can be modeled as a classification task where the input is some short DNA sequence, and the label indicates whether the sequence can be bound by a TF of interest. It has recently been pointed out that incorporating non-sequence-based data modalities such as chromatin state can improve TF binding prediction [START_REF] Karimzadeh | Virtual chip-seq: Predicting transcription factor binding by learning from the transcriptome[END_REF]. However, since our method is focused on the modeling of biological sequences, our experiments are limited to sequence data only.

Datasets and evaluation metric

In our experiments, we consider the datasets used by [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF], consisting of fragment peaks in 506 different ENCODE ChIP-seq experiments. While negative sequences are originally generated by random dinucleotide shuffling, we also train our models with real negative sequences not bound by the TF, a task called motif occupancy by [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF]. Both datasets have a balanced number of positive and negative samples, and we therefore measure performances by the area under the ROC curve (auROC). As noted by [START_REF] Karimzadeh | Virtual chip-seq: Predicting transcription factor binding by learning from the transcriptome[END_REF], even though classical, this setting may lead to overoptimistic performance: the real detection problem is more difficult as it involves a few binding sites and a huge number of non-binding sites.

Hyperparameter tuning

We discuss here the choice of different hyperparameters used in CKN and DeepBind-based CNN models.

Hyperparameter tuning for CNNs.

In DeepBind [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF], the search for hyperparameters (learning rate, momentum, initialization, weight decay, DropOut) is computationally expensive. We observe that training with the initialization mechanism proposed by [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF] and the Adam optimization algorithm [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] leads to a set of canonical hyper-parameters that perform well across datasets, and to get rid of such an expensive dataset-specific calibration step. The results we obtain in such a setting are consistent with those reported by [START_REF] Alipanahi | Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[END_REF] (and produced by their software package) and by [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] (see Supplementary Figure 12 and13). Overall, this simplified strategy comes with great practical benefits in terms of speed.

Specifically, to choose the remaining parameters such as weight decay, we randomly select 100 datasets from DeepBind's datasets, and we use one quarter of the training samples as validation set, on which the error is used as a proxy of the generalization error. We observe that neither DropOut [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF], nor fully connected layers bring significant improvements, which leads to an even simpler model.

Hyperparameter tuning for CKNs.

The hyperparameters of CKNs are also fixed across datasets, and we select them using the same methodology described above for CNNs. Specifically, this strategy is used to select the bandwidth parameter σ and the regularization parameter λ (see Supplementary Figure 2 and 3), which is then fixed for all the versions of CKN and on either the DeepBind's or [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] datasets. For unsupervised CKN, the regularization parameter is dataset-specific and is obtained by a five-fold cross validation. To train CKN-seq, we initialize the supervised CKN-seq with the unsupervised method (which is parameter-free) and use the alternating optimization update presented in section 2.2. We use the Adam algorithm [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] to update the filters and the L-BFGS algorithm [START_REF] Zhu | Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization[END_REF] to optimize the prediction layer. The learning rate is fixed to 0.01 for both CNN and CKN. The logistic loss is chosen to be the loss function for both this and the next protein homology detection task. All the models only use one layer. The choice of filter size, number of filters, and number of layers are also discussed in Section 3.3.

Performance benchmark

We compare here the auROC scores on test datasets between different CKN and DeepBindbased CNN models.

Performance on entire datasets. Both supervised and unsupervised versions of CKN-seq show performance similar to DeepBind-based CNN models (Figure 2), on either the Deep-Bind [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] datasets.

Performance on small-scale datasets. When few labeled samples are available, unsupervised CKNs achieve better predictive performance than fully supervised approaches that are hard to regularize. Specifically, we have selected all the datasets with less than 5000 training samples and reevaluated the above models. The results are presented in the top part of Figure 3. As expected, we observe that the data-augmented version outperform the corresponding unaugmented version for all the models, while the supervised CKN is still dominated by the unsupervised CKN. Finally, the hybrid version of CKN-seq presented in section 2.3 performs nearly as well as the unsupervised one while only using 32 times fewer (only 128) filters. It is also more robust to the perturbation intensity used in augmentation than the data-augmented version (detailed choice and study of perturbation intensity can be found in Supplementary Figure 4 and5).

We obtain similar results on the [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] datasets as shown in the middle part of Figure 3, except that the data-augmented unsupervised CKN-seq does not improve performance over its unaugmented counterpart.

Protein homology detection

Protein homology detection is a fundamental problem in computational biology to understand and analyze the structure and function similarity between protein sequences. String kernels, see, e.g., [START_REF] Christina S Leslie | The spectrum kernel: A string kernel for svm protein classification[END_REF][START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF], [START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF], [START_REF] Rangwala | Profile-based direct kernels for remote homology detection and fold recognition[END_REF], have shown state-of-the-art prediction performance but are computationally expensive, which restricts their use to small-scale datasets. In comparison, CKN-seq and CNN models are much more computationally efficient and also turn out to achieve better performance, which we show in the rest of this section. Specifically, we consider the remote homology detection problem and benchmark different methods on the widely-used SCOP 1.67 dataset from [START_REF] Hochreiter | Fast model-based protein homology detection without alignment[END_REF], including 102 superfamily recognition tasks and extending the positive training samples with Uniref50. The number of training protein samples for each task is around 5000, whose length varies from tens to thousands of amino acids. Under our formulation, positive protein sequences are taken from one superfamily from which one family is withheld to serve as test samples, while negative sequences are chosen from outside the target family's fold.

Regarding the training of CNN and CKN-seq, we adopt the same setting as for the TF binding prediction task and the same methodology for the selection of hyper-parameters. A larger bandwidth parameter σ = 0.6 is selected (in contrast to σ = 0.3 in Section 3.1) due to the larger number of (20) characters in protein sequences. Further details about the validation scores obtained for various parameters are presented in Supplementary Figure 123. We also use max pooling in CKN-seq to aggregate feature vectors instead of mean pooling, which shows better performance in this problem. We fix the filter size to be 10 which seems computationally intractable for the exact algorithms, such as trie-based algorithm, for computing mismatch kernels [START_REF] Pavel P Kuksa | Scalable algorithms for string kernels with inexact matching[END_REF].

Profile-based methods [Kuang et al., 2005, Rangwala and[START_REF] Rangwala | Profile-based direct kernels for remote homology detection and fold recognition[END_REF] have shown very good performance on this task but suffer a few limitations as pointed out by [START_REF] Hochreiter | Fast model-based protein homology detection without alignment[END_REF], including computation time and interpretability. Nevertheless, we propose an approach which integrates profiles with CKN models. Specifically, we compute the position-specific probability matrix (PSPM) using PSI-BLAST for all the sequences in SCOP 1.67 dataset, following the same protocols as [START_REF] Rangwala | Profile-based direct kernels for remote homology detection and fold recognition[END_REF]. PSI-BLAST is performed against Uniparc2 filtering out all the sequences after 2015, which leads to a database similar to the NCBI nonredundant database used by [START_REF] Rangwala | Profile-based direct kernels for remote homology detection and fold recognition[END_REF]. We encode the sequences in Uniref50 using the BLOSUM62 position-independent probability matrix [START_REF] Henikoff | Amino acid substitution matrices from protein blocks[END_REF]] by replacing each character with its corresponding substitution probability in BLOSUM62. Finally, we train CKN models by replacing each sequence in our kernel (2) with the square root of its corresponding PSPM (or BLOSUM62). The training and hyperparameters remain unchanged. Performance on entire datasets. Besides auROC, we also use auROC50 (area under the ROC up to 50 false positives) as evaluation metric, which is extensively used in the litera-

Method auROC auROC50

GPkernel [START_REF] Håndstad | Motif kernel generated by genetic programming improves remote homology and fold detection[END_REF] 0.902 0.591 SVM-pairwise [START_REF] Liao | Combining pairwise sequence similarity and support vector machines for detecting remote protein evolutionary and structural relationships[END_REF] 0.849 0.555 Mismatch [START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF] 0.878 0.543 LA-kernel [START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF] 0.919 0.686 LSTM [START_REF] Hochreiter | Fast model-based protein homology detection without alignment[END_REF] [START_REF] Leslie | Fast string kernels using inexact matching for protein sequences[END_REF][START_REF] Saigo | Protein homology detection using string alignment kernels[END_REF]. Table 1 shows that unsupervised CKN-seq and CNN achieve similar performance and supervised CKN-seq achieves even better performance while they outperform all typical string kernels including local alignment kernel. They also outperform the LSTM model proposed by [START_REF] Hochreiter | Fast model-based protein homology detection without alignment[END_REF]. Finally, training CKN-seq is much faster than using string kernel-based methods. While training string kernel-based models requires hours or days [START_REF] Hochreiter | Fast model-based protein homology detection without alignment[END_REF], training CNN or CKN-seq are done in a few minutes. In our experiments, the average training time for CNN and supervised CKN-seq is less than 3 minutes on a single cluster with a GTX1080 TI GPU and 8 CPU cores of 2.4 GHz, while training an unsupervised CKN-seq with 16384 filters (which seems to be the maximal size that can be fit to GPU memory and gives 0.956 and 0.792 respectively for auROC and auROC50) needs 30 minutes in average. We also notice that using a random sampling instead of K-means in unsupervised CKN-seq reduces the training time to 6 minutes without loss of performance. By contrast, the training time for a local alignment kernel is about 4 hours. Profile-based CKN-seq models show substantial improvements over their non-profile counterparts, including the BLOSUM62-based CKN-seq which uses the position-independent BLO-SUM62 probability matrix instead of one-hot encoding to encode sequence characters. Supervised CKN-seq shows comparable results to the best performing methods. The performance may be further improved by computing the profiles for the extended sequences in Uniref50.

Performance on subsampled datasets. We simulate situations where few training samples are available by subsampling only 500 class-balanced training samples for each dataset. We reevaluate the above CNN and CKN models, the data-augmented versions and also the hybrid method. The results (bottom part of Figure 3) are similar to the ones obtained for the TF binding prediction problem except that supervised version of CKN-seq performs remarkably well in this task. We also notice that CKN-seq versions trained with only 500 samples outperform the best string kernel trained with all training samples.

Hyperparameter Study

We now study the effect of hyperparameters and focus on the supervised version of CKN, which is more interpretable than the unsupervised one. Table 2: Tomtom motif p-value comparison of CKN-seq and CNN for different distance functions, see [START_REF] Gupta | Quantifying similarity between motifs[END_REF].

Both CNN and CKN-seq with one layer achieve better performance with a filter size of 12 for every fixed number of filters (Supplementary Figure 7). Since this optimal value is only slightly larger than the typical length of the motifs for TFs [START_REF] Alexander | Why transcription factor binding sites are ten nucleotides long[END_REF], we deduce that the prediction mainly relies on a canonical motif while the nearby content has little contribution.

Increasing the number of filters improves the auROCs for both models regardless of the filter size, in line with the observation in [START_REF] Zeng | Convolutional neural network architectures for predicting DNA-protein binding[END_REF] for CNNs. This improvement saturates when more than 128 filters are deployed, sometimes leading to overfitting (Supplementary Figure 6). We observe the same behavior for the unsupervised version of CKN-seq (Supplementary Figure 6), but usually with much larger saturation bar (larger than 4096 for TF binding prediction and 32768 for protein homology detection). When using only 16 filters, CKN-seq shows better performance than DeepBind-based CNNs. This is an advantage as large numbers of filters make the model redundant and harder to interpret.

Model interpretation and visualization

In this section, we study the ability of a trained CKN-seq model to capture motifs and generate accurate and informative sequence logos. We use here simulated data since the true motifs are generally not known in practice. To simulate sequences containing some given motifs represented by a PWM, we follow the methodology adopted by Shrikumar et al. [2017a] and generate 500 training and 100 test samples. We train a 1-layer CKN-seq and CNN on two tasks of the respective motif FOXA1 and GATA1 [START_REF] Kheradpour | Systematic discovery and characterization of regulatory motifs in encode tf binding experiments[END_REF], using the same hyperparameter settings as previously. We fix the filter size and number of filters to 12 and 16 to avoid capturing too many redundant features. Both models achieve about 0.99 for the auROC on test set. The trained CNN is visualized by using the approach introduced by Alipanahi et al. [2015]. Specifically, all sequences from the test set are fed through the convolutional and rectification stages of the CNN, and only the k-mers that passed the activation threshold (which is 0 by default) were aligned to generate a PWM and the trained CKN is visualized by using the approach presented in section 2.4, i.e., solving min µ∈M ψ 0 (µ) -ψ 0 (z i ) 2 with a projected gradient descent method. The best recovered motifs (in the sense of information content) are compared to the true motifs using Tomtom [START_REF] Gupta | Quantifying similarity between motifs[END_REF].

Motifs recovered by CKN-seq and CNN are both aligned to the true motifs (Figure 16). The logos given by CKN-seq are more informative and match better with the ground truth in terms of any distance measures (Table 2). This suggests that CKN-seq may be able to find more accurate motifs. We also perform the same experiments with more training samples (see Supplementary Figure 16). We observe that CKN-seq achieves small p-values in both data regimes while p-values for CNN are larger when few training samples are available.

Discussion

We have introduced a convolutional kernel for sequences which combines advantages of CNNs and string kernels. The resulting CKN-seq is a special case of CNN which generalizes the mismatch kernel to motifs -instead of discrete k-mers -and makes it task-adaptive and scalable.

CKN-seq retains the ability of CNNs to learn sequence representations from large datasets, leading to slightly better performance than classical CNNs on a TF binding site prediction task and on a protein homology detection task. The unsupervised version of CKN-seq keeps the kernel formalism, which makes it easier to regularize and thus leads to good performance on small-scale datasets despite the use of a huge number of convolutional filters. A hybrid version of CKN-seq performs equally well as its unsupervised version but with much fewer filters. Finally, the kernel interpretation also makes the learned model more interpretable and thus recovers more accurate motifs.

The fact that CKNs retain the ability of CNNs to learn feature spaces from large training sets of data while enjoying a RKHS structure has other uncharted applications which we would like to explore in future work. First, it will allow us to leverage the existing literature on kernels for biological sequences to define the bottom kernel K 0 , possibly capturing other aspects than sequence motifs. More generally, it provides a straightforward way to build models for non-vector objects such as graphs, taking as input molecules or protein structures. Finally, it paves the way for making deep networks amenable to statistical analysis, in particular to hypothesis testing. This important step would be complementary to the interpretability aspect, and necessary to make deep networks a powerful tool for molecular biology beyond prediction.

C Hyperparameter study

We discuss here in more detail the effect of the number and size of convolutional filters and number of layers on CNN and CKN performances. We also present the discussions on the perturbation intensity in data-augmented and hybrid variants of CKN-seq. For some of the following comparisons, we also include the oracle model, which represents the best performance achievable by choosing the optimal parameter in comparison for each dataset (whereas parameters used in our experiments are fixed across datasets). The experiment shows that a dataset-dependent parameter calibration step could possibly improve the performance, but that the potential gain would be relatively small. Number of filters, filter size and number of layers. We show in Figure 6 that increasing the number of filters improved the performance for both supervised and unsupervised variants of CKN-seq. Furthermore, the improvement of prediction performance of the supervised one was saturated when more than 128 convolutional filters were deployed.

Both CNN and CKN-seq with one layer achieve better performance with a filter size of 12 for every fixed number of filters (Figure 7). Since this optimal value is only slightly larger than the typical length of the motifs for TFs, we deduce that the prediction mainly relies on a canonical motif while the nearby content has little contribution. However if one is interested in motif discovery only, running the algorithm with larger filter size may be of interest whenever one believes that some TF binding sites are explained by larger motifs. Increasing the number of convolutional layers in CNNs has been shown to decrease its performance. By contrast, it does not affect the performance of CKN-seq when using a sufficient number of convolutional filters (Figure 8). Multilayer architectures allow to learn richer or more complex descriptors such as co-motifs, but may require a larger amount of data. They would also make the interpretation of the trained models more difficult. When training with 2-layer CKN models, we also notice that increasing the number of filters from 64 to 128 at the first layer or that from 16 to 64 at the second layer does not improve performance (Figure 9). Perturbation intensity in data-augmented and hybrid CKN. We have shown that data augmentation improves both supervised and unsupervised CKN-seq. The hybrid approach has further improved data-augmented CKN-seq. We study here how the amount of perturbation used in augmenting training samples impacts performance. Specifically, we characterize the perturbation intensity by the percentage of changed characters in a sequence and show in Figure 11 the behavior of CKN-seq when increasing the amount of perturbation. By leveraging the best data-augmented unsupervised model on validation set, we train our hybrid variant and show its performance when increasing the amount of perturbation (Figure 10). We observe that the hybrid variant is more robust to larger amount of perturbation applied in the training samples than simply data-augmented one. Note that the results are consistent to those obtained on validation set (Section B).

D Effect of hyperparameter calibration in CNN

We study here how hyperparameter calibration as used in DeepBind could affect performance and training time for CNNs. For the calibrated variant of CNN, we used the same hyperparameter search scheme used in DeepBind for the CNN, with 30 randomly chosen calibration settings and 6 training trials across the data sets.

The calibrated variant slightly outperformed hyperparameter-fixed CNN and showed similar performance to CKN-seq in the TF binding prediction task while it didn't achieve better performance in the protein homology detection task (Figure 12).

On the other hand, training a calibrated CNN is much slower compared to hyperparameterfixed CNN or CKN-seq. To make a fair comparison, we reimplemented and evaluated both DeepBind and CKN-seq in Pytorch. Our reimplemented model achieved almost identical performance to the original DeepBind (left panel of Figure 13

Figure 4 :

 4 Figure 4: Motifs recovered by CKN-seq (middle row) and by CNN (bottom row) compared to the true motifs (top row)

Figure 1 :

 1 Figure 1: Validation of weight decay in CNNs for DeepBind's datasets (left) and SCOP 1.67 and its subsampled datasets (middle and right); k = 12 and 10 respectively for each task; p = 128 for both tasks.

Figure 2 :

 2 Figure 2: Validation of the bandwidth parameter σ for DeepBind's datasets (left) and SCOP 1.67 (right). The regularization parameter is fixed to 1e-6 and 1.0 and k = 12 and 10 respectively for each task; p = 128 for both tasks.

Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 34567 Figure 3: Validation of the regularization λ for DeepBind's datasets (left) and SCOP 1.67 (right). The bandwith parameter is fixed to 0.3 and 0.6 and k = 12 and 10 respectively for each task; p = 128 for both tasks.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Influence of the number of filters for 2-layer supervised CKN-seq on DeepBind's datasets.

Figure 12 :Figure 13 :

 1213 Figure 12: Comparison of calibrated CNN and universal models; left: DeepBind's dataset and right: SCOP 1.67 dataset

  ) in DeepBind's Datasets. In order to quantify the gain in training time for hyperparameter-fixed models, we measured the average training time on 50 different datasets for original DeepBind, our reimplemented DeepBind and CKN-seq on a Geforce GTX Titan Black GPU. The right panel of Figure 13 shows that training a CKN-seq model is about 25 times faster than training the original DeepBind model and 5 times faster than our reimplemented version.

Figure 15 :Figure 16 :

 1516 Figure 15: Pairwise comparison of CKN-seq and CNN on DeepBind, Zeng and SCOP 1.67 datasets. The metric is auROC for the two earlier datasets and auROC50 for the latter. The middle and bottom lines show performance of models trained on small-scale datasets.

  Figure 2: Performance comparison of CNN and CKN-seq on the DeepBind (left) and Zeng et al. [2016] (right) datasets. Number of filters for CNN and CKN-seq was set to 128 while it was 4096 for unsupervised CKN-seq. The average auROCs for CNN, CKN-seq and unsupervised CKN-seq are 0.931, 0.936, 0.937 on the DeepBind datasets and 0.803, 0.807, 0.804 on the Zeng et al. [2016] datasets. The pink and black lines respectively represent mean and median. Pvalues are from one-sided Wilcoxon signed-rank test. All the following figures are obtained in the same way.

		2.1e-11			2.6e-22
		1.00
		0.95			0.90
		0.90			0.85
	auROC	0.85		auROC	0.80
		0.80			0.75
		0.75			0.70
		0.70
					0.65
		CNN	CKN-seq Method	unsup CKN-seq	CNN	CKN-seq Method	unsup CKN-seq

Table 1 :

 1 Average auROC and auROC50 for SCOP 1.67 benchmark ture

	]	0.942	0.773
	CNN (128 filters)	0.960	0.799
	CKN-seq (128 filters)	0.965	0.819
	CKN-seq (128 filters) + BLOSUM62	0.973	0.835
	unsup CKN-seq (32768 filters)	0.958	0.806
	Profile-based methods
	Mismatch-profile on SCOP 1.53 [Kuang et al., 2005]	0.980	0.794
	SW-PSSM on SCOP 1.53 [Rangwala and Karypis, 2005]	0.982	0.904
	CKN-seq (128 filters) + profile	0.986	0.906
	unsup CKN-seq (4096 filters) + profile	0.968	0.863

It is also possible to introduce a concept of zero-padding for sequences, such that Pi(x) may contain characters outside of the original sequence, when i is close to the sequence boundary, see Section 3.

https://www.uniprot.org

Acknowledgements

We thank H. Zeng for sharing the experimental results of Zeng et al. [2016].

Funding

This work has been supported by the grants from ANR (MACARON project ANR-14-CE23-0003-01 and FAST-BIG project ANR-17-CE23-0011-01) and by the ERC grant number 714381 (SOLARIS).

Supplementary Material

In the supplementary material, we present details and additional experiments mentioned in the paper.

A Details about the convolutional kernel

In the definition of convolutional kernel, a bottom kernel K 0 was defined, for any z, z ′ in R d , as

where κ : u → e 1 σ 2 (u-1) . When z and z ′ are one-hot encoded vectors of k-mers, we have

which recovers a Gaussian kernel.

B Choice of model hyperparameters

We justify here the choice of the hyperparameters used in our experiments, including weight decay for CNNs, regularization parameter, bandwidth parameter in exponential kernel and perturbation intensity used in data-augmented CNN, CKN and hybrid model. We denote respectively by k the filter size and p the number of filters.

The scores for the following experiments are computed on a validation set, which is taken from one quarter of the training samples for each dataset and the models are trained on the rest of the training samples. For DeepBind's datasets, we only perform validation on 100 randomly sampled datasets, which save a lot of computation time and should give similar results when using all datasets.

Weight decay for CNN. The choice of weight decay is validated on the validation set as shown in Figure 1.

Bandwidth parameter in exponential kernel. The choice of the bandwidth parameter is only validated for supervised CKN-seq and the same value is used for the unsupervised variant.

Figure 2 shows the scores on the validation set when the other hyperparameters are fixed. The same choice as DeepBind's dataset is applied to Zeng's dataset.

Regularization parameter. The choice of the regularization parameter is validated following the same protocol as the bandwidth parameter. Figure 3 shows the scores on the validation set.

Perturbation intensity in data-augmented and hybrid model.

The perturbation amount used in the data-augmented CNN, CKN and the hybrid variant of CKN are also validated on the corresponding validation set. The scores are shown in Figure 5.

E Influence of fully connected layer in CNN

The authors of DeepBind have used a fully connected layer in their model. However, we found that there was no significant gain with this supplementary layer in our experiments, as shown in Figure 14.

F Pairwise comparison of CKN and CNN

We include here some scatter plots to illustrate the pairwise comparison on each individual dataset of DeepBind and Zeng. The results are shown in 15.

G Model interpretation and visualization

We perform the same experiments as in section 3.4 of the paper but on a larger datasets, with 9000 training samples and 1000 test samples. Motifs recovered by CKN-seq and CNN were aligned to the true motifs (Figure 16) while the logos given by CKN-seq are more informative and match better with the ground truth in terms of any distance measures (Table 3). The same conclusions can be drawn as in the small-scale case.