Frédéric Chyzak

Assia Mahboubi

Thomas Sibut-Pinote

Enrico Tassi

A Computer-Algebra-Based Formal Proof of the Irrationality of ζ(3)

This paper describes the formal verification of an irrationality proof of ζ(3), the evaluation of the Riemann zeta function, using the Coq proof assistant. This result was first proved by Apéry in 1978, and the proof we have formalized follows the path of his original presentation. The crux of this proof is to establish that some sequences satisfy a common recurrence. We formally prove this result by an a posteriori verification of calculations performed by computer algebra algorithms in a Maple session. The rest of the proof combines arithmetical ingredients and some asymptotic analysis that we conduct by extending the Mathematical Components libraries. The formalization of this proof is complete up to a weak corollary of the Prime Number Theorem.

Introduction

The irrationality status of the evaluations of the Riemann ζ-function at positive odd integers is a long-standing challenge of number theory. To date, ζ( 3) is the only one known to be irrational, although recent advances obtained by Rivoal [START_REF]Propriétés diophantiennes de la fonction zêta de Riemann aux entiers impairs[END_REF] and Zudilin [START_REF] Zudilin | One of the numbers ζ(5), ζ(7), ζ(9), ζ(11) is irrational[END_REF] showed that one at least of the numbers ζ(5), . . . , ζ [START_REF] Delahaye | Dealing with algebraic expressions over a field in Coq using Maple[END_REF] must be irrational. The number ζ(3) is sometimes referred to as the Apéry constant, after Roger Apéry who first proved that it is irrational [START_REF] Apéry | Irrationalité de ζ(2) et ζ(3)[END_REF]. As reported by van der Poorten [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF], Apéry announced this astonishing result by giving a rather obscure lecture that raised more skepticism than enthusiasm among the audience. His exposition indeed involved a number of suspicious assertions, proclaimed without a proof, among which was a mysterious common recurrence for two given sequences (see Lemma 2). After two months of work, however, Cohen, Lenstra, and van der Poorten completed, with the help of Zagier, a verification of Apéry's proof.

Theorem 1 (Apéry, 1978). The constant ζ(3) is irrational.

Almost at the same time, symbolic computation was emerging as a scientific area of its own, getting fame with the Risch algorithm [START_REF] Risch | The solution of the problem of integration in finite terms[END_REF] for indefinite integration. It gradually provided efficient computer implementations and got attention in experimental mathematics. Beside commutative algebra, differential and recurrence equations remained a central research topic of computer algebra over the years. In particular, the sequences used by Apéry in his proof belong to a class of objects well known to combinatorialists and computer-algebraists. Following seminal work of Zeilberger's [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF], algorithms have been designed and implemented in computer-algebra systems, which are able to obtain linear recurrences for these sequences. For instance the Maple packages gfun and Mgfun (both distributed as part of the Algolib [START_REF] Algolib | [END_REF] library) implement these algorithms, among other. Basing on this implementation, Salvy wrote a Maple worksheet [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF] that follows Apéry's original method but interlaces Maple calculations with humanwritten parts, illustrating how parts of this proof, including the discovery of Apéry's mysterious recurrence, can be performed by computations.

In the present paper, we describe a formal proof of Theorem 1, based on a Maple session, in the Coq proof assistant. The computer-algebra system is used in a skeptical way [START_REF] Harrison | A skeptic's approach to combining HOL and Maple[END_REF], to produce conjectures that are a posteriori proved formally. Alternative proofs are known for Theorem 1, as for instance the elegant one proposed by Beukers [START_REF] Beukers | A note on the irrationality of ζ(2) and ζ(3)[END_REF] shortly after Apéry. Our motivation however was to devise a protocol to obtain formal proofs of computer-algebra-generated recurrences in a systematic way. Interestingly, this work challenges the common belief in the computer-algebra community that such an a posteriori checking can be automatized. In addition to the formal verification of these computeralgebra-produced assertions, we have also machine-checked the rest of the proof of irrationality, which involves both elementary number theory and some asymptotic analysis. The latter part of the proof essentially consists in a formal study of the asymptotic behaviors of some sums and of tails of sums. Our formal proof is complete, up to a weak corollary of the repartition of prime numbers that we use as an assumption.

In Section 2, we outline a proof of Theorem 1. Section 3 presents the algorithms which are run in the Maple session we base on. Section 4 describes the formalization of the formal proof we obtain from the data produced by the computer-algebra system. Section 5 provides some concluding remarks and some perspectives for future work.

Our Maple and Coq scripts will be found at http://specfun.inria.fr/ zeta-of-3/.

From Apéry's recurrence to the irrationality of ζ(3)

In this section, we outline the path we have followed in our formalization, highlighting the places where we resorted to more elementary variants than Salvy or van der Poorten. In particular, Section 2.3 describes a simple argument we devised to simplify the proof of asymptotic considerations.

Overview

In all what follows, a Cauchy real (number) x is a sequence of rational numbers (x n ) n∈N for which there exists a function m x : Q → N, such that for any ǫ > 0 and any indices i and j,

having i ≥ m x (ǫ) and j ≥ m x (ǫ) implies |x i -x j | ≤ ǫ. Proposition 1. The sequence z n = n m=1 1 m 3 is a Cauchy real.
The Cauchy real of Proposition 1 is our definition for ζ [START_REF] Apéry | Irrationalité de ζ(2) et ζ(3)[END_REF]. Consider the two sequences a and b of rational numbers defined as:

a n = n k=0 n k 2 n+k k 2 , b n = a n z n + n k=1 k m=1 (-1) m+1 n k 2 n+k k 2 2m 3 n m n+m m . (1)
Introducing the auxiliary sequences of real numbers: 3) was a rational number, then for n large enough, every σ n would be a (positive) integer, preventing σ from tending to zero.

δ n = a n ζ(3) -b n , σ n = 2ℓ

Arithmetics, number theory

We extend the usual definition of binomial coefficients n k for n, k ∈ N to n, k ∈ Z by enforcing the Pascal triangle recurrence n+1 k+1 = n k+1 + n k for all n, k ∈ Z. Although this extension is not required by the present proof, it spares us some spurious considerations about subtraction over N. Binomial coefficients being integers, a n is also an integer for any nonnegative n ∈ N.

An important property of the sequence (b n ) n∈N is that for any n ∈ N, the product 2ℓ 3). We follow the argument described by Salvy in [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF], and show that each summand in the double sum defining b n has a denominator that divides 2ℓ 3 n : after a suitable re-organization in the expression of the summand, which uses standard properties of binomial coefficients, this follows easily from the following slightly less standard property of theirs: Lemma 1. For any integers i, j, n such that 1 ≤ j ≤ i ≤ n, j i j divides ℓ n . Lemma 1 is considered as folklore in number theory. Its proof consists in showing that for any prime p, the p-valuation of j i j is smaller than the one of ℓ n . Standard presentations of Apéry's proof make use of the asymptotic bound ℓ n = e n (1+o(1)) , which is a corollary of the distribution of the prime numbers. A bound 3 n is however tight enough for our purpose and has been proved by several independent and elementary proofs, for instance by Hanson [14] and Feng [START_REF] Feng | A simple elementary proof for the inequality dn < 3 n[END_REF]. However, we have not yet formalized any proof of this ingredient, which is completely independent from the rest of the irrationality proof. More precisely, our formal proof is parametrized by the following assumption: Proposition 2. There exists two positive rationals K and r, with r 3 < 33, such that for any large enough integer n, ℓ n < Kr n .

Consequences of Apéry's recurrence

The Cauchy sequence (b n /a n ) n∈N tends to ζ(3), thus δ n tends to zero. In this section, we prove that it does so fast enough to compensate for ℓ 3 n , while being positive. The starting point is Apéry's recurrence, (3) below: Lemma 2. For n ≥ 0, the sequences (a n ) n∈N and (b n ) n∈N satisfy the same second-order recurrence:

(n + 2) 3 y n+2 -(17n 2 + 51n + 39)(2n + 3)y n+1 + (n + 1) 3 y n = 0.

(3)

Salvy's worksheet [START_REF] Salvy | An Algolib-aided version of Apéry's proof of the irrationality of ζ(3)[END_REF] demonstrates in particular how to obtain this common recurrence by Maple calculations, performed by the Algolib library [START_REF] Algolib | [END_REF]. Following van der Poorten [START_REF] Van Der Poorten | A proof that Euler missed: Apéry's proof of the irrationality of ζ(3)[END_REF], we next use Lemma 2 (and initial conditions) to obtain a closed form of the Casoratian w n = b n+1 a n -b n a n+1 . Indeed, we prove w n = 6 (n+1) 3 for n ≥ 2. From this, we prove that δ n , and hence σ n , is positive for any n ≥ 2. We also use the closed form to estimate the growth of δ in terms of a.

The result is that there exists a positive rational number K such that δ n ≤ K an for large enough n. Finally, the zero limit of σ n follows from Proposition 2, the behaviour of δ, and Lemma 3 below, which quantifies that a grows fast enough.

Lemma 3. 33 n ∈ O(a n ).
Proof. Introduce the sequence ρ n = a n+1 /a n and observe that ρ 51 > 33. We now show that ρ is increasing. Define rational functions α and β so that the conclusion of Lemma 2 for a n rewrites to a n+2 -α(n)a n+1 +β(n)a n = 0 for n ≥ 0. Now, for any n ∈ N, introduce the homography

h n (x) = α(n) -β(n) x , so that ρ n+1 = h n (ρ n ). Let x n be the largest root of x 2 -α(n)x + β(n).
The result follows by induction on n from the fact that h( [1, x n ]) ⊂ [1, x n ] and from the observation that

ρ 2 ∈ [1, x 2 ].
⊓ ⊔

Algorithms on sequences in computer algebra

Lemma 2 is the bottleneck in Apéry's proof. Both sums a n and b n in there are instances of parametrised summation: they follow the pattern

F n = β(n) k=α(n) f n,
k in which the summand f n,k , potentially the bounds, and thus the sum, depend on a parameter n. This makes it appealing to resort to the algorithmic paradigm of creative telescoping, which was developed for this situation in computer algebra.

In order to operate on sequences, computer algebra substitutes implicit representations for explicit representations in terms of named sequences (factorial, binomial, etc). This is the topic of Section 3.1. A typical example of parametrised summation by this approach is provided by the identity n k=0 n k = 2 n : from an encoding of the summand n k by the recurrences

n + 1 k = n + 1 n + 1 -k n k , n k + 1 = n -k k + 1 n k , (4)
deriving the relation, with finite difference with respect to k in right-hand side,

n + 1 k -2 n k = n + 1 k + 1 - n k + 1 - n + 1 k - n k (5)
is sufficient to derive the explicit form 2 n , as will be explained below.

Recurrences as a data structure for sequences

The implicit representation fruitfully introduced by computer algebra to deal with sequences are systems of linear recurrences. In this spirit, ∂-finite sequences are algebraic objects that model mathematical sequences and enjoy nice algorithmic properties. Notably, the finiteness property of their definition makes algorithmic most operations under which the class of ∂-finite sequences is stable.

A ∂-finite sequence (see [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF] for a complete exposition of the subject) is an element of a module over the non-commutative ring A of skew polynomials in the indeterminates S n and S k , with coefficients in the rational-function field Q(n, k), and commutation rule

S i n S j k c(n, k) = c(n + i, k + j)S i n S j k . A skew polynomial P = (i,j)∈I p i,j (n, k)S i n S j k ∈ A acts on a "sequence" f by (P • f ) n,k = (i,j)∈I p i,j (n, k)f n+i,k+j
, where subscripts denote evaluation. For example for f n,k = n k , the recurrences (4) once rewritten as equalities to zero can be represented as P • f = 0 for P = S n -n+1 n+1-k and P = S k -n-k k+1 , respectively. To any ∂-finite sequence f , one associates the set of skew polynomials that annihilate it. This set, {P ∈ A : P • f = 0} is a left ideal of A, named the annihilating ideal of f , and denoted ann f . A non-commutative extension of the usual Gröbner-basis theory is available, together with algorithmic analogues. In this setting, a good representation of a ∂-finite sequence is obtained as a Gröbner basis of ann f for a suitable ordering on the monomials in S n and S k . For the example of f n,k = n k , a Gröbner basis consists of both already-mentioned skew polynomials encoding (4). In general, a Gröbner basis provides us with a (vectorial) basis of the quotient module A/ ann f . This basis can be explicitly written in the form B = {f n+i,k+j } (i,j)∈U , where the finite set U of indices is given as the part under the classical stair shape of the Gröbner-basis theory. Given a Gröbner basis GB for ann f , the normal form NF(p, GB) is unique for any p ∈ A. Again in the binomial example, the finite set is U = {(0, 0)}, and normal forms are rational functions. This is the basis of algorithms for a number of operations under which the ∂-finite class is stable, which all process by looking for enough dependencies between normal forms: application of an operator, addition, product. The case of summing a sequence (f n,k ) into a parametrised sum F n = n k=0 f n,k is more involved: it performs according to the method of creative telescoping [START_REF] Zeilberger | The method of creative telescoping[END_REF], in two stages. First, an algorithmic step determines pairs (P, Q) satisfying

P • f = (S k -1)Q • f ( 6
)
with P ∈ A ′ and Q ∈ A, where A ′ is the subalgebra Q(n) S n of A. To continue with our example f n,k = n k , Eq. ( 5) can be recast into this framework by choosing P = S n -2 and Q = S n -1. Second, a systematic but not fully algorithmic step follows: summing [START_REF] Chyzak | Do creative-telescoping algorithms provide complete proofs? a formal study of Apéry's theorem[END_REF] for k between 0 and n + deg Sn P yields

(P • F ) n = (Q • f ) k=n+deg Sn P +1 -(Q • f ) k=0 . ( 7
)
Continuing with our binomial example, summing (5) (or its equivalent form ( 6)) for k from 0 to n+1 (and taking special values into account) yields [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF] with right-hand side canceling to zero. The formula [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF] in fact assumes several hypotheses that hold not so often in practice; this will be formalized by Eq. ( 8) below.

n+1 k=0 n+1 k - 2 n k=0 n k = 0, a special form of

Apéry's sequences are ∂-finite constructions

The sequences a and b in (1) are ∂-finite: they have been announced to be solutions of (3). But more precisely, they can be viewed as constructed from

step explicit form GB operation input(s) 1 c n,k = n k 2 n+k k 2 C direct 2 an = n k=1 c n,k A creative telescoping C 3 dn,m = (-1) m+1 2m 3 ( n m )( n+m m ) D direct 4 s n,k = k m=1 dn,m S creative telescoping D 5 zn = n m=1 1 m 3 Z direct 6 u n,k = zn + s n,k U addition Z and S 7 v n,k = c n,k u n,k V product C and U 8 bn = n k=1 v n,k B creative telescoping V
Table 1: Construction of a n and b n : At each step, the Gröbner basis named in column GB, which annihilates the sequence given in explicit form, is obtained by the corresponding operation on ideals, with input(s) given on the last column.

"atomic" sequences by operations under which the class of ∂-finite sequences is stable. This is summarised in Table 1.

Both systems C and D are first-order systems obtained directly as easy consequences of (4); they consist respectively of expressions for c n+1,k and c n,k+1 in terms of c n,k and of expressions for d n+1,k and d n,k+1 in terms of d n,k . The case of Z is almost the same: it is not a parametrised summation but an indefinite summation. A (univariate, second-order) recurrence is easily obtained for it, without referring to any creative telescoping.

For each of C, D, and V , which undergo a summation operation, we obtain creative-telescoping pairs (P, Q): one for C for a set U = {(0, 0)}; one for V for a set U = {(0, 0), (1, 0), (0, 1)}; four for D for the same set. In all cases, we have had our computer-algebra program (informally) ensure that the corresponding P cancels the sum. For C and V , the single P thus obtained is (trivially) a Gröbner basis of the annihilating ideal of A or B, respectively. But for D, the four P have to be recombined, leading to three operators.

It should be observed that Gröbner bases is the only data used by computeralgebra algorithms in the program above, including when simplifying the righthand side in [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF] and in its generalization to come, Eq. ( 8) below. For instance, computer algebra computes the system V from the systems C and U alone, without resorting to any other knowledge of particular values of c and u in Table 1, and so does our formal proof to verify that the pointwise product xy is annihilated by V whenever C and U respectively annihilate sequences x and y. In other words, although a ∂-finite sequence is fully determined by a system of recurrences and sufficiently many initial conditions (that is, values of u i,j for small values of i and j), we do not maintain those values along our proofs.

Our formal proof as well models each sequence by a system of recurrences obtained solely from the operators of the Gröbner basis. We hence bet that the computer-algebra implementation of the algorithmic operations of addition, product, and summation, as well as the parts of our Maple script relying on less algorithmic operations do not take decisions based on private knowledge they could have on their input, viewed as a specific solution to the recurrence system used to encode it. Would the implementation do so without letting us know, then our a posteriori verification would have required guessing an appropriate description of this additional knowledge, like operators for specializations of the sequences. Fortunately, the Mgfun package we used has the wanted property.

Provisos and sound creative telescoping

Observe the denominators in (4): they prevent the rules to be used, respectively when k = n + 1 and k = -1. For example, one can "almost prove" Pascal's triangle rule by

n + 1 k + 1 - n k + 1 - n k = n + 1 n -k n -k k + 1 - n -k k + 1 -1 n k = 0 × n k = 0,
but this requires k = -1 and k = n. Therefore, this does not prove Pascal's rule for all n and k. The phenomenon is general: computer algebra is unable to take denominators into account. This incomplete modelling of sequences by algebraic objects may cast doubt on these computer-algebra proofs, in particular when it comes to the output of creative-telescoping algorithms. By contrast, in our formal proofs, we augmented the recurrences with provisos that restrict their applicability. In this setting, we validate a candidate identity like the Pascal triangle rule by a normalization modulo the elements of a Gröbner basis plus a verification that this normalization only involves legal instances of the recurrences. In the case of creative telescoping, Eq. ( 6) takes the form:

(n, k) / ∈ ∆ ⇒ (P • f ,k ) n = (Q • f ) n,k+1 -(Q • f ) n,k , (8)
where ∆ ⊂ Z 2 guards the relation and where f ,j denotes the univariate sequence obtained by specializing the second argument of f to j. Thus our formal analogue of Eq. ( 7) takes this restriction into account and has the shape

(P • F ) n = (Q • f ) n,n+β+1 -(Q • f ) n,α + r i=1 i j=1 p i (n) f n+i,n+β+j + α≤k≤n+β ∧ (n,k)∈∆ (P • f ,k ) n -(Q • f ) n,k+1 + (Q • f ) n,k , (9)
for F the sequence with general term F n = n+β k=α f n,k . The proof of identity ( 9) is a straightforward reordering of the terms of the left-hand side, (P • F ) n = r i=0 p(n) F n+i , after unfolding the definition of F and applying relation [START_REF] Cohen | Construction of real algebraic numbers in Coq[END_REF] everywhere allowed in the interval α ≤ k ≤ n + β. The first part of the righthand side is the usual difference of border terms, already present in Eq. [START_REF] Chyzak | Non-commutative elimination in Ore algebras proves multivariate identities[END_REF]. The middle part is a collection of terms that arise from the fact that the upper bound of the sum defining F n depends linearly on n and that we do not assume any nullity of the summand outside the summation domain. The last part, which we will call the singular part, witnesses the possible partial domain of validity of relation [START_REF] Cohen | Construction of real algebraic numbers in Coq[END_REF]. The operator P is a valid recurrence for the sequence F if the right-hand side of Eq. ( 9) normalizes to zero, at least outside of an algebraic locus that will guard the recurrence.

Formal proof of the common recurrence

This section describes the computer-algebra-aided formal proof of Lemma 2, based on a Maple session implementing the program described in Table 1.

Generated operators, hand-written provisos, and formal proofs

For each step in Table 1, we make use of the data computed by the Maple session in a systematic way. Figure 1 illustrates this pattern on the example of step 7. As mentioned in Section 3.3, we annotate each operator produced by the computer-algebra program with provisos (see below) and turn it this way into a conditional recurrence predicate on sequences. To each sequence in the program corresponds a file defining the corresponding conditional recurrences, for instance annotated_recs_c, annotated_recs_u, and annotated_recs_v for c, u, and v, respectively. More precisely these files contain all the operators obtained by the Maple script for a given sequence, not only the Gröbner basis. We use rounded boxes to depict the files that store the definitions of these predicates. These are generated by the Maple script which pretty-prints its output in Coq syntax, with the exception of the definition of provisos. Throughout this section, a maple leaf tags the files that are generated by our Maple script. Yet automating these annotations is currently out of reach.

In our formal proof, each step in Table 1 consists in proving that some conditional recurrences on a composed sequence can be proved from some conditional recurrences known for the arguments of the operation. We use square boxes to depict the files that store these formal proofs. The statement of the theorems proved in these files are composed from the predicates defined in the round boxes: a dashed line points to (predicates used to state) conclusions and a labelled solid line points to (predicates used to state) hypotheses. The constants Sn, Sk, and CT_premise are recurrences predicates, defined in terms of a bound variable c. Constants Sn and Sk are elements of the Gröbner basis. The definition of these recurrences is named to reflect the term it rewrites, e.g., the left-hand sides in (4): these names are the result of pretty-printing the (skew) monomial that encodes these left-hand sides, the prefix S standing for "shift". For example Sn is the name of a recurrence defining c n+1,k , while SnSk would be for c n+1,k+1 . Rewriting a given term with such an equation makes the term decrease for the order associated with the Gröbner basis. Another part of the file defines the recurrences obtained from a creative-telescoping pair (P, Q) generated for the purpose of the summation defining the sequence a. Observe that these generated definitions feature named provisos that are in fact placeholders. In the preamble of the file, displayed on Listing 1.2, we provide by a manual annotation a concrete definition for the proviso of each recurrence defined in the generated part. Observe however that part of these definitions can be inferred from the coefficients of the recurrences. For example the k = n + 1 condition in precond.Sn, the proviso of recurrence Sn, is due to the denominator (-n -1 + k) 2 of the coefficient (Sn00n k).

Module precond. Definition Sn n k :

= (k = n + 1) ∧ (n = -1). Definition Sk n k := (k + 1 = 0) ∧ (n = 0). Definition CT_premise n k := (n ≥ 0) ∧ (k ≥ 0) ∧ (k < n). End precond.

Listing 1.2: Hand-written provisos in annotated rec c

In the last part of the file, see Listing 1.3, a record collects the elements of the Gröbner basis C. Maple indeed often produces a larger set of annihilators for a given sequence, for instance CT_premise in Listing 1.1 is related to a creative telescoping pair but not to the Gröbner basis. Also, the Gröbner basis can be obtained by refining a first set of annihilators, which happens at step 4 of Table 1.

Formal proofs of a conditional recurrence

We take as a running example the file ops_for_a, which models step 2 in Table 1. This file proves theorems about an arbitrary sequence c satisfying the recurrences in the Gröbner basis displayed on Listing 1.3, and about the sequence a by definite summation over c. The formal proof of lemma P_eq_Delta_k_Q is an instance of Eq. 8. Using this property, we prove that the sequence a verifies a conditional recurrence associated to the operator P . As suggested in Section 3.1, this proof consists in applying the lemma punk.sound_telescoping, which formalizes a sound creative telescoping and in normalizing to zero the resulting right-hand side of Eq. 9. Listing 1.4 displays the first lines of the corresponding proof script, which select and name the three components of the right-hand side of Eq. 9, with self-explanatory names. The resulting proof context is displayed on Listing 1.5.

Theorem recApery_a n (nge2 : n ≥ 2) : P a n = 0. Proof. rewrite (punk.sound_telescoping P_eq_Delta_k_Q). set boundary_part := (X in X + _ + _). set singular_part := (X in _ + X + _). set overhead_part := (X in _ + _ + X).

Listing 1.4: Begining of a proof of sound creative telescoping

boundary_part := Q c n (n + 1) -Q c n 0 singular_part := \sum_(0 ≤ i < n + 1 | precond.CT_premise n i) P (c i) n -(Q c n (i + 1) -Q c n i) overhead_part := \sum_(0 ≤ i < degree P) \sum_(0 ≤ j < i) Pi n * c (n + i) (n + j + 1) ============================ boundary_part + singular_part + overhead_part = 0
Listing 1.5: Corresponding goal In Listing 1.5, (c i) denotes the expression (fun x => c x i), P i denotes the i-th coefficient of the polynomial P, and degree P the degree of P (two in this specific case). Note that we have access to degree P because in addition to the definition displayed on Listing 1.1, we also have at our disposal a list representation [:: P0; P1] of the same operator.

The proof of the goal of Listing 1.5, proceeds in three steps. The first step is to inspect the terms in singular_part and to chase ill-formed denominators, like n -n. These can arise from the specialisations, like k = n, induced when unrolling the definition of (the negation of) precond.CT_premise. In our formalization, a division by zero is represented by a conventional value: we check that these terms vanish by natural compensations, independently of the convention, and we keep only the terms in singular_part that represent genuine rational numbers. The second step consists in using the annihilator ann_c of the summand to reduce the resulting expression under the stairs of the Gröbner basis. In fact, this latter expression features several collections of terms, that will be reduced to as many independent copies of the stairs. In the present example, we observe two such collections: (i) terms that are around the lower bound (n, 0) of the sum, of the form c n,0 , . . . c n,s ; (ii) terms that are around the upper bound (n, n) of the summation, of the form c n,n , . . . , c n,n+s for a constant s. The border terms induce two such collections but there might be more, depending in particular on the shape of the precond.CT_premise proviso. For example, the sum n k=0 (-1) k n k 3k n = (-3) n leads to a proviso involving n = 3k + 1 and similar terms: an additional category of terms around (n, n/3) drifts away from both (n, 0) and (n, n) when n grows.

============================ P2 n * c (n + 2) n + P1 n * c (n + 1) n + P0 n * c n n + Q00 n n * c n n + P1 n * c (n + 1) (n + 1) + P2 n * c (n + 2) (n + 1) + P2 n * c (n + 2) (n + 2) = 0
Listing 1.6: Terms around the upper bound The collection of terms around the upper bound in our running example is displayed on Listing 1.6. The script of Listing 1.7 reduces this collection under the stairs of ann_c, producing the expression displayed on Listing 1.8. The premise of each rule in this basis being an integer linear arithmetic expression, we check its satisfiability using our front-end intlia to the lia proof command [START_REF] Besson | Fast reflexive arithmetic tactics the linear case and beyond[END_REF], which automates the formal proof of first-order formulae of linear arithmetics. The third and last step consists in checking that the rational-function coefficient of every remaining evaluation of c is zero. For this purpose, we start by unfolding the definitions of the coefficients P2, Sn0. Previous steps kept them carefully folded as these values play no role in the previous normalizations but can lead to pretty large expressions if expanded, significantly slowing down any other proof command. The resulting rational function is proved to be zero by a combination of the field [START_REF] Delahaye | Dealing with algebraic expressions over a field in Coq using Maple[END_REF] and lia [START_REF] Besson | Fast reflexive arithmetic tactics the linear case and beyond[END_REF] proof commands. The former reduces the rational equation into a polynomial one between the cross product of two rational functions. This equation is then solved by the ring proof command [START_REF] Mahboubi | Proving equalities in a commutative ring done right in Coq[END_REF]. The algebraic manipulations performed by field produce a set of non-nullity conditions for the denominators. These are solved by the lia proof command. To this end, our Maple script generates rational fractions with factored denominators, that happen to feature only linear factors in these examples.

Composing closures and reducing the order of B

Figure 2 describes the global dependencies of the files proving all the steps in Table 1. In order to complete the formal proof of Lemma 2, we verify formally in Fig. 2: Formal proofs of Table 1 Fig. 3: Formal proof of Lemma 2 file algo_closures that each sequence involved in the construction of a n and b n is a solution of the corresponding Gröbner system of annotated recurrence, starting from c n , d n , and z n and applying the lemmas proved in the ops_for_* files all the way to the the final conclusions of ops_for_a and ops_for_b. This proves that a n is a solution of the recurrence (3) but provides only a recurrence of order four for b n . In file reduce_order, we prove that b as well satisfies the recurrence (3) using four evaluations b 0 , b 1 , b 2 , b 3 that we compute in file initial_conds.

Conclusion

Formally proving the consequences of Apéry's recurrence

The present paper focuses on the computer-algebra-aided part of our formalization. Another significant part of it addresses the asymptotic study of some sequences of rational numbers. Formalizing this part of the proof requires introducing a type with more inhabitants than just rational numbers, for it deals with the properties of limits of sequences, notably the constant ζ(3) itself: we hence need to include at least computable real numbers. The construction of Cauchy reals proposed by Cohen in his PhD [START_REF] Cohen | Formalized algebraic numbers: construction and first-order theory[END_REF] turned out to perfectly suit our needs, although we had to enrich the existing library with a few basic properties. Indeed, we use the type of Cauchy real numbers mostly to state our theorems, like Theorem 1, but the proofs only involve reasoning with rational values of the sequences. We also benefited from the proof commands provided by Cohen [START_REF] Cohen | Formalized algebraic numbers: construction and first-order theory[END_REF] to postpone the definition of "large enough" bounds or moduli of convergence. The proof of Lemma 3 however involves non-rational real numbers that are not defined as converging sequences but as roots of polynomials. We resort in this case to Cohen's construction of real algebraic numbers [START_REF] Cohen | Construction of real algebraic numbers in Coq[END_REF] and benefit constructively from the decidability of comparison relations. Navigating between the different natures of numbers this irrationality proof features, integer, rational numbers, real algebraic numbers, and Cauchy reals was made possible by the genericity of the core Mathematical Component libraries [START_REF]Mathematical Components Libraries[END_REF].

The proof of the same Lemma 3 involves two computational steps, namely the observations that ρ 51 > 33 and that ρ 2 ∈ [1, x 2 ]. These numbers are rather large: for instance the normalized fraction ρ 51 features integers with more than 70 digits. The issue here is to obtain a formal proof of these facts that both fits with the rest of the proof of Lemma 3 and does not take ages to compute. We used the CoqEAL [START_REF] Cohen | Refinements for free! In CPP[END_REF] library and the framework it provides for writing parametric code which can be executed either symbolically, in abstract parts of a proof, or with appropriate data structures for efficient computations. The programs that implement the computation of ρ n and x n are very short and boil down to evaluating a (rather small) rational fraction defined recursively. At this place and when computing evaluations of the sequence (see Section 4.4), computations use (interpreted) binary integer arithmetics and are almost instantaneous.

Asymptotic behavior of lcm(1, . . . , n)

We plan to complete our irrationality proof by providing a formal proof of Proposition 2 following the proof of Hanson [14]. This involves considerations of elementary number theory as well as asymptotic-analysis studies that are rather similar to the ones we dealt with in the present work. The Prime Number Theorem has been machine-checked by Harrison [START_REF] Harrison | Formalizing an analytic proof of the Prime Number Theorem[END_REF] in HOL-Light. We considered importing this result into Coq by using the automated translation tool developed by Keller and Werner [START_REF] Keller | Importing HOL Light into Coq[END_REF]. Unfortunately, according to Keller, the import of such an involved proof is currently beyond reach because of its size.

Formal proofs on objects of size typical of computer algebra

The size of the mathematical objects in our formalization makes the interactive writing of proof scripts quite unusual and challenging. For example, the recurrence P • y = 0 of order four satisfied by the sequence b spans over 8,000 lines when pretty-printed. Both proofs of the fact P • b = 0, named recApery_b in our formalization, and of its premise

P • v = (S k -1)Q • v, named P_eq_Delta_k_Q,
end by normalizing rational functions to 0. Figure 4 reports on the size of the polynomials that occur in these proofs, together with the amount of time needed by Coq's field proof command to check that the rational functions are zero.

These objects can be easily manipulated by a computer-algebra system: Maple normalizes P • b to zero in less than 2 seconds. Coq, in its latest stable version (8.4pl3), takes roughly 4 minutes, but of course it produces and checks a formal proof of such normalization. To achieve this reasonable timing, we contributed to version 8.4pl3 a patch to substantially improve the performances of Navigating expressions of size typical of computer algebra Under the circumstances described above, where the goal to be proved cannot even be displayed, the need for uncommon proof commands arises. The pattern matching facilities offered by the set command [START_REF] Gonthier | A small scale reflection extension for the Coq system[END_REF] can be used to probe the goal for a subterm of the shape u (n + _) (k + _), or to identify all the occurrences of a closed terms.

Unfortunately not all uncommon needs are easily satisfiable with standard proof commands. A typical one is to invoke a computer-algebra system to rearrange an expression in order decide how to proceed with the proof. This was performed by hand, massaging text files to accommodate the small syntactic differences between the syntax of Coq and Maple. This task has been automated in the past in the Coq-Maple bridge developed by Delahaye and Mayero [START_REF] Delahaye | Dealing with algebraic expressions over a field in Coq using Maple[END_REF], for the purpose of enhancing proof automation in Coq. We foresee to generalise this work by allowing arbitrary Maple commands to let one explore algebraic expressions, and not necessarily prove the correctness of the results.

Proof-language implementation

The size of the mathematical expressions poses additional efficiency challenges when combined with the ability of the Coq logic to identify terms up to computation. For example a pattern like (3 * _) is matched by searching for the head symbol * verbatim, but its non-wildcard argument, 3 here, is compared taking computation into account. Given the size of the expressions, the fast filtering performed on the head symbol is not sufficient to discard enough false positives. Comparing two terms up to computation often triggers a full normalization before failing and this is an expensive operation.

The definitions of integer and rational numbers that are part of the Mathematical Components library are designed for the purpose of a comprehensive library of theorems, and not for fast computations: this makes the aforementioned pattern-matching operation often intractable. The commonly accepted idea of having two representations of the same concept, one good for reasoning and one good for computing, linked by a morphism does not work for expressions of this size, as switching from one representation to the other one relies again on pattern matching to push the morphism through the huge expression. For instance, we had to craft boilerplate code, which carefully controls computations and use locked casts for large constant numbers, in order to make the field and lia tactic work on these data-structures.

One possible work-around would be to have available, in the proof language, two distinct pattern-matching facilities: one that eases reasoning steps by identi-fying terms up to computation, and another one that performs a dumb syntactic comparison and is efficient even when applied to large, homogeneous expressions.

Theoretical and practical limitations of data structures

Algebraic rewriting with provisos An obstruction to turning our approach into a complete protocol is that we do not know how to determine the provisos other than by trial and error. This is connected to the fact that recurrences that have well understood rewriting properties when forming a Gröbner basis (a priori) lose all this nice structure when decorated with provisos. We do not understand what the new critical pairs have to be and how to ensure we give ourselves complete information with our selection of rules-with-provisos that simply lift a Gröbner basis. To the best of our knowledge, there is no literature on the topic yet.

We have been lucky with the few cases in the present formalization, in that we could just guide our reduction by the usual strategy of division by a Gröbner basis: all implied integer premises could be solved, without us understanding why they would.

Polynomial data structures Manipulating polynomials in computer algebra is made natural by the ease to keep them in collected form, whether automatically or by easy-to-call procedures. On the other hand, our treatment of recurrences amounts mostly to collecting polynomials and skew polynomials, through the encoding suggested in Section 3.1. A (skew) polynomial data structure, with computational product and division, and the corresponding theory, would help keeping goals tidy and more manageable, especially during the rewriting steps. In particular, the current strategy of reserving the call to field to the end is certainly a cause of blow-up of the goals and inefficiency of the computations.

In addition, making it generic with respect to the nature of skew-polynomial commutation would provide similar functionality for both recurrence and differential equations, paving the way to the formalization of integral identities by a known differential counterpart to creative telescoping.

Building on top of state-of-the-art libraries

It was a help to have at our disposal the broad scope of the Mathematical Components library. It was a cause for trouble, too. This research started as an interdisciplinary activity: half of the authors were newcomers to formal proofs. As a matter of fact understanding the sophisticated design of the library turned out to be a very arduous task for them.

We can identify in the combined use of implicit coercions, notations, and structure inference the major source of complexity. Their formalization activity was additionally hindered by the fact that various features of the Coq system do not interact well with the design patterns employed in the Mathematical Components library. The most notable example being search, that is almost ineffective in the context of the algebraic hierarchy.

We believe such problems could be alleviated by: introductory material specific to the Mathematical Components library and written with newcomers in mind; a "teaching/debug" mode in which type inference explicates steps at the desired level of detail; finally a proper searching facility.

Fig. 1 :

 1 Fig. 1: Proving that V is C × U

(Listing 1 . 1 :

 11 * Coefficients of every recurrence, P, and Q. *) Definition Sn00 n k := (n + 1 + k) 2 / (-n -1 + k) 2 . Definition Sk00 n k := (-n + k) 2 * (n + 1 + k) 2 / (k + 1) 4 . Definition P0 n := (n + 1) 3 . ... (* Conditional recurrences. *) Definition Sn c := ∀ n k, precond.Sn n k → c (n + 1) k = Sn00 n k * c n k Definition Sk c := ∀ n k, precond.Sk n k → c n (k + 1) = Sk00 n k * c n k (* Operators P and Q. *) Definition P c n := P0 n * c n + P1 n * c (n + 1) + P2 n * c (n + 2). ... (* Statement P = ∆ k Q. *) Definition CT_premise c := ∀ n k, precond.CT_premise n k → P (c k) n = Q c n (k + 1) -Q c n k. Generated part of annotated rec c

(

  * Choice of recurrences forming a Groebner basis. *) Record Annihilators c := { Sn : Sn c; Sk : Sk c }. Listing 1.3: Selection of a Gröbner basis

Require

  Import annotated_recs_c. Variables (c : int → int → rat) (ann_c : Annihilators c). Theorem P_eq_Delta_k_Q : CT_premise c. Proof. ... Qed. Let a n := \sum_(0 ≤ k < n + 1) c n k.

Listing 1 . 8 :

 18 rewrite (ann_c.Sk (n + 2) (n + 1)); last by intlia. rewrite (ann_c.Sk (n + 2) n); last by intlia. rewrite (ann_c.Sk (n + 1) n); last by intlia. rewrite (ann_c.Sn (n + 1) n); last by intlia. rewrite (ann_c.Sn n n); last by intlia. set cnn := c n n. Fail set no_more_c := c _ _. Listing 1.7: Reduction modulo the Gröbner basis of c ============================ P2 n * Sn00 (n + 1) n * Sn00 n n * cnn + P1 n * Sn00 n n * cnn + P0 n * cnn + Q00 n n * cnn + P1 n * Sk00 (n + 1) n * Sn00 n n * cnn + P2 n * Sk00 (n + 2) n * Sn00 (n + 1) n * Sn00 n n * cnn + P2 n * Sk00 (n + 2) (n + 1) * Sk00 (n + 2) n * Sn00 (n + 1) n * Sn00 n n * cnn = 0 Rational function with folded coefficients

Fig. 4 :

 4 Fig. 4: Statistics about the polynomials in ops_for_b and their normalization


Jonathan Grizou
email: jonathan.grizou@inria.fr

Iñaki Iturrate
email: inaki.iturrate@epfl.ch

Luis Montesano
email: montesano@unizar.es

Pierre-Yves Oudeyer
email: pierre-yves.oudeyer@inria.fr

Manuel Lopes
email: manuel.lopes@inria.fr

Calibration-Free BCI Based Control

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

Introduction

EEG-based brain-computer interfaces (BCI) have been used successfully to control different devices, such as robotic arms and simulated agents, using self-generated (e.g. motor imagery) and event-related potentials signals (see (Millán et al. 2010) for a review). Error-related potentials (ErrPs) are one kind of event-related potential appearing when the user's expectation diverges from the actual outcome [START_REF] Falkenstein | ERP components on reaction errors and their functional significance: A tutorial[END_REF]. Recently, they have been used as feedback instructions for devices to solve a user's intended task [START_REF] Chavarriaga | Learning from EEG error-related potentials in noninvasive brain-computer interfaces[END_REF]Iturrate, Montesano, and Minguez 2013a).

As in most BCI applications, ErrP-based BCI requires a calibration phase to learn a decoder (e.g. a classifier) that translates raw EEG signals from the brain of each user into meaningful instructions. This calibration is required due to specific characteristics of the EEG signals: non-stationary nature [START_REF] Vidaurre | Toward unsupervised adaptation of LDA for brain-computer interfaces[END_REF], large intra-and inter-subject variability [START_REF] Polich | On the relationship between EEG and P300: individual differences, aging, and ultradian rhythms[END_REF], and variations induced by the task (Iturrate, Montesano, and Minguez 2013b). The presence of an explicit calibration phase, whose length and frequency is hard to tune and is often tedious and impractical for users, hinders the deployments of BCI applications out of the lab. Thus, calibration free methods are an important step to apply this technology in real applications [START_REF] Millán | Combining braincomputer interfaces and assistive technologies: state-of-theart and challenges[END_REF].

Despite the importance of calibration-free BCI, there are only few related works. For long term operation using sensory-motor rhythms, it is possible to adapt the decoder online [START_REF] Vidaurre | Towards a cure for BCI illiteracy[END_REF]. In invasive BCI, [START_REF] Orsborn | Closed-Loop Decoder Adaptation on Intermediate Time-Scales Facilitates Rapid BMI Performance Improvements Independent of Decoder Initialization Conditions[END_REF] proposed a method to learn from scratch and in closed loop a decoder for known targets using predefined policies to each target. However, the approach needs a warm-up period of around 15 minutes. For P300 spellers, Kindermans et al. proposed a method to auto-calibrate the decoder by exploiting multiple stimulations and prior information (Kindermans, Verstraeten, and Schrauwen 2012;Kindermans et al. 2012;[START_REF] Tangermann | Zero training for BCI-reality for BCI systems based on event-related potentials[END_REF]. They exploit the particular fact that only one event out of six encodes a P300 potential in the speller paradigm.

Our main contribution is a calibration-free BCI method that infers simultaneously and seamlessly an EEG decoder of error-related potentials while controlling a device to achieve a sequential task. The core idea of the method is to assume a distribution of possible tasks, and infer the interpretation of EEG signals and the task by selecting the hypothesis which best explains the history of interaction. This inference can be continuously run and updated as new data comes in, which removes the need for an explicit calibration.

This method is inspired from our previous work [START_REF] Grizou | Robot learning simultaneously a task and how to interpret human instructions[END_REF] which considered a robotic setting and speech utterances as feedback signals. In the current work we improve the algorithm formalism, the robustness to noisy high-dimensional signals (e.g. EEG), and show that it is possible to use model-based planning relying on the uncertainty about the task and the feedback signals interpretation to explore the space efficiently while learning.

We also present an evaluation of this method with online experiments where four users control an agent in a virtual world. The results show that the proposed method allows to learn a good signal decoder and solve the task efficiently without any explicit calibration. Offline experiments show that our unsupervised trained decoder achieves similar performances than calibration based systems and illustrate the benefits of our planning strategy for speeding up learning.

Calibration-Free BCI Based Control BCI control based on feedback signals BCI control based on feedback signals differs from classical brain-computer interfaces in the sense that the user does not actively deliver commands to the device, but only delivers feedback about actions performed by the device [START_REF] Chavarriaga | Learning from EEG error-related potentials in noninvasive brain-computer interfaces[END_REF]Iturrate, Montesano, and Minguez 2013a). In this setting, the device needs to actively execute a sequence of several actions to solve the task and to be able to learn an intelligent behavior from the feedback. This idea can be seen as a shared control strategy [START_REF] Millán | Combining braincomputer interfaces and assistive technologies: state-of-theart and challenges[END_REF], where both the user and the device help each other to solve a task.

Essentially, this BCI control follows an iterative sequential process where the device performs an action which is in turn assessed by the user. This assessment will elicit potentials into the user's brain that can be recorded using EEG and will be different for "correct" and "wrong" assessments. The potentials elicited in the user's brain after performing assessments are called error-related potentials [START_REF] Ferrez | Error-related EEG potentials generated during simulated Brain-Computer interaction[END_REF]. After a calibration phase and once a usable decoder of these signals is available, user's assessments can be translated into (normally binary) feedback, which the device can use to adapt its behavior.

This control based on user's assessments decoded from brain signals can be exemplified for a reaching task, where the user wants to reach a target position unknown by the system. The device performs several discrete actions (e.g. moving left or right), and learns from the feedback given by the user. After several iterations, if the meanings of the EEG assessment signals provided by the users are known, the device can infer which is the user's desired position and how to reach it. The following section explains how we can achieve similar performances without knowing the brain signal decoder beforehand. We will use the term "virtual label" or "label" to denote the interpretation of a given EEG signal as a feedback instruction (e.g. "correct" or "wrong").

Simultaneous Estimation of Task and Signal Model

This section formalizes the problem of executing a task when the mapping between raw EEG signals to a discrete label among a set of pre-defined labels is unknown. The main idea is depicted in Figure 1 for a toy 1D example. The user wants the device to reach the right-most state. For each device's action, he provides a feedback signal which encode whether the action executed is "correct" or "wrong" according to the intended target. Such signals are generated from an underlying model which maps a binary label ("correct" or "wrong") to a continuous signal. However, neither the user's desired target nor the labels associated to the user's feedback signals are known.

Considering that we can define a finite set of task hypotheses (e.g. reaching one of a finite number of states), we can infer the labels that should be provided by the user with respect to each hypothesis. Then, given a particular interaction history, it is possible to compute a different signal decoder for each task hypothesis. The key point is that only the correct hypothesis will assign the correct labels to all feedback sig- On top is a 1D grid world with the hypothetic target state marked with a T letter. Notice that this hypothetic target is different for each of the three hypotheses. The arrows in each cell indicates what action should elicit a positive feedback, i.e. the optimal policy with respect to the hypothetic target. The user intended target is shown as the shaded blue state at the right extremity of the 1D world. Note that the system does not have access to this information. The correct hypothesis is the one on the Left [(a)] where the T state is the same as the shaded blue state. Below the 1D grid world, the signals received from the user are represented in a 2D feature space. They represent the user assessment signals of the past history of device's actions (e.g. moving randomly left and right). Our algorithm assigns virtual labels (green for "correct" and red for "wrong") to those signals with respect to their respective hypothetic target. Notice that the data points are the same for all three hypotheses, only their respective labels differ. With the virtual labels being assigned, we can compute the corresponding 2D Gaussian distributions estimates for each class (shown as colored ellipses) and each hypothesis. While for the correct hypothesis [(a)] the Gaussian distributions shows a large separability, the overlap increases as the hypothetic target (T) moves away from the real (blue shaded) one [(b), (c)]. This property can be exploited to estimate the correct hypothesis and the model generating the signals.

nals (Figure 1a), while the other hypotheses will gradually mix both classes as the hypothetic target gradually differs more from the correct one (Figure 1b and1c). Therefore, the hypothesis which provides the decoder with best accuracy and compactness can be selected as the most probable one. In the remainder of this section we show how this property can be exploited to estimate the target and the model generating the feedback signals.

Formally, we represent the problem as a discrete or continuous set of states s 2 S, a finite set of possible actions a 2 A, and a set of possible tasks, or targets, for which the system is able to plan the best sequence of actions. In the 1D example in Figure 1, the state space is composed of seven discrete states, and the agent can perform unitary directional actions, i.e. move one cell left or right. During an interactive learning session, the agent will proactively perform actions which will in turn be evaluated as "correct" or "wrong" by the user with respect to the desired target state on the grid. In our case the feedback signals are error-related potentials measured in the brain activity of the subject.

Let e i 2 R n denote the feature vector of the EEG measurements obtained at iteration i after the device performed action a i in state s i . The label z i 2 {c, w} of each feedback signal belongs to one of two classes ("correct" or "wrong"). Following [START_REF] Blankertz | Single-trial analysis and classification of ERP components: A tutorial[END_REF], we will model the EEG signals using independent multivariate normal distributions for each class (N (µ c , Σ c ) and N (µ w , Σ w )). We will denote by θ this set of parameters {µ c , Σ c ,µ w , Σ w }.

We assume the system has access to a set of task hypotheses ξ 1 ,...,ξ T which includes the task the user wants to solve. We do not make any particular assumption on how the task is represented but we assume that for each particular task ξ we are able to compute a policy π ξ which represents the probability of choosing a given action a in state s, π ξ (s, a)=p(a|s, ξ). As mentioned above, these are the policies that, conditioned on the task, provide labels to the feedback signals of a state-action pair (e.g. in a reaching task, progressing towards the goal will generate "correct" feedback while moving apart from it will generate "wrong"feedback).

The method aims to infer which task ξ the user wants to solve based on the user's feedback signals extracted from EEG measurements collected while the agent executes some actions. Following the analysis of Figure 1, a sensible option to estimate the task is to measure the coherence of the signal model for each possible task using the virtual labels provided by the target policy. In other words, the best (ξ, θ) pair provides the lowest predictive error on the observed signals p(e|s, a, ξ, θ), which themselves were collected in the recent history of interaction. One possible way of solving this problem is to maximize the expected classification rate:

ξ, θ = arg max ξ,θ E e [δ(z(s, a, ξ),z(e, θ))] (1)
where δ() is an indicator function, z(s, a, ξ) is the label (correct or wrong) corresponding to the execution of action a in state s under task ξ and z(e, θ) is the label provided by classifying the EEG signal e under the Gaussian classifier parameterized by θ. The expected classification rate (Ecr) can be explicitly written dependent on the task and decoder model:

Ecr (δ(z(s, a, ξ),z(e, θ))) = = X k∈{c,w} p(z = k|s, a, ξ)p(z = k|e, θ) (2)
where p(z = k|s, a, ξ) represents the probability of the user assigning label k when assessing action a in state s according to task ξ. We add a noise term to cope with those situations where the user assessment may be wrong. The model for probability of correct assessment is then:

p(z = c|s, a, ξ)= ⇢ 1 -α if a = arg max a π ξ (s, a) α otherwise
(3) with α modeling the assessment error rate of the user, which was set to 0.1 for our experiments. Finally, the term p(z = k|e, θ) is just the probability that the signal e belong to class k under the Gaussian model provided by θ and is given by:

p(z = k|e, θ)= p(e|z = k, θ)p(z = k) P l∈{c,w} p(e|z = l, θ)p(z = l) = N (e|µ k , Σ k )p(z = k) P l∈{c,w} N (e|µ l , Σ l )p(z = l) (4)
As we do not have a priori knowledge on the user intended meaning, we assume that it is equiprobable p(z = c)= p(z = w). We will factorize the optimization process using the fact that given a task ξ, the estimation of θ under the Gaussian model is straightforward. It basically requires to compute the maximum-likelihood estimate θ ML ξ using the labels associated to target ξ. We could also consider a prior distribution on the parameters and update it with new observations. Using the labels of target ξ, the estimation of θ under the Gaussian model described above simply returns to the computation of the posterior mean µ z and covariance Σ z for each class z 2 {c, w}. In order to avoid numerical problems when estimating the covariance for a low number of examples, a regularization term was applied to penalize very large and very small eigenvalues [START_REF] Friedman | Regularized discriminant analysis[END_REF]:

Σ z =(1-λ)Σ z + λ trace(Σ z ) n I n (5)
with n the feature dimension, I n the identity matrix of size n, and λ the regularization term which was set to 0.5 for our experiments. An automatic adaptation of the regularization could also be considered [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF].

Using equation 2 to estimate the expected classification rate is difficult because we ideally want to estimate it on future, never observed, data. A possible solution is to use cross-validation or bootstrapping methods using the available data. However, for small amounts of data, these methods result in estimates with high variance [START_REF] Bengio | No unbiased estimator of the variance of k-fold cross-validation[END_REF]) and computational cost.

Alternatively, we propose to use another approximation of the expected classification rate, the Bhattacharyya coefficient. This coefficient has been related to the classification error of Gaussian models [START_REF] Kailath | The divergence and Bhattacharyya distance measures in signal selection[END_REF] and is inversely proportional to the classification rate. Although there is no analytical relation between the coefficient and the classification rate, it is possible to derive bounds and good empirical approximations [START_REF] Lee | Bayes error evaluation of the Gaussian ML classifier[END_REF]. The Bhattacharyya coefficient ρ 2 [0, 1] between the Gaussian distributions associated to label "correct" (N (µ c , Σ c )) and "wrong" (N (µ w , Σ w )) is:

ρ = e -D B (θ) (6)
where

D B is the Bhattacharyya distance: D B (θ)= 1 8 (µ c -µ w ) T ( Σc+Σw 2 ) -1 (µ c -µ w )+ 1 2 ln ⇣ det( Σc +Σw 2 ) √ detΣcdetΣw
⌘ . Finally, we approximate the expected classification rate as:

Ecr / 1 -ρ (7)

Confidence on Target Estimation

Now that we have an estimation of the expected classification rate, we need to take a decision with respect to which task is the one intended by the user. To do so we should compare the expected classification rate of every task hypothesis ξ t with t 2 {1,...,T}. The hypothesis whose associated model has the highest expected classification rate, i.e. the lowest value of ρ, is expected to be the one intended by the user, however it is meaningless to define an absolute threshold on the value of the expected classification rate itself. Indeed, different people generate different signals which result in classifiers of different qualities. To bypass this problem we rely on a voting system where we attribute each hypothesis ξ t a weight that is updated at every iteration.

We rely on a pseudo-likelihood metric that for each hypothesis ξ t accumulate expected classification rate over time:

L(ξ t )= N Y i=1 1 -ρ ξt i (8)
with N the current number of iteration and ρ ξt i the Bhattacharyya coefficient associated to task ξ t using all data up to time i. By normalizing the pseudo-likelihood values between every hypothesis, we obtain what can be viewed as the probability of each target:

p(ξ t )= L(ξ t ) P u∈{1,...,T } L(ξ u ) (9)
Once a target reaches a probability threshold β we consider it as being the correct one, i.e. the one intended by the user. We used β =0.99.

Estimation of further Tasks and Online Re-Estimation of Signal Model

Once we have identified a first task, the user can change his desired target and our system has to identify the new target. However, the model of the feedback signals does not change and does not have to be re-learned from scratch. Indeed, once the system has correctly identified a task, it is possible to reuse acquired data by assigning to all the previously collected signals the labels associated to the previously estimated task.

The use of the Bhattacharyya coefficient provides a simple and efficient way to estimate a first target from scratch (when the number of examples is small) but does not allow a fast adaptation to new targets as the majority of collected signals now belong to the previous task. To avoid this problem we compute a classifier (e.g. a Gaussian Bayes classifier), and use it similarly to calibrated approaches in BCI. For this we factorize the joint distribution:

p(ξ t , θ | D ξt i ) / p(ξ t | D ξt i )p(θ | ξ t ,D ξt i ) (10)
where D ξt i contains all the quadruplet (e i ,s i ,a i ,z i ) up to time i, with the associated labels z i assigned with respect to task ξ t . The factorization makes explicit that given the task ξ t , the distribution p(θ | ξ t ,D ξt i ) can be easily evaluated using the labels of each target. We approximate this posterior using the maximum likelihood point estimate θ ML ξt per target. For the term p(ξ t | D ξt i ), we use a recursive Bayes filter:

p(ξ t | D ξt i ) / p(e i | ξ t , (s, a) i ,D ξt i-1 )p(ξ | D ξt i-1 ) ⇡ p(e i | θ ML ξt )p(ξ t | D i-1 ). (11)
Notice that we are keeping a different symbol model θ ML ξt for each possible target ξ t , the maximum likelihood estimation needs to be done in relation to a dataset D ξt i-1 which include the expected labels of target ξ t up to time i -1.

We use the same threshold mechanism as described in previous subsection to decide whether or not a task can be considered the correct one. Whenever a task is identified, its labels are transferred to the quadruplets D ξt i of all the other tasks to correct the prior for the next step with the right labels. This scheme performs a long term adaptation of θ to accommodate slight variations of EEG such as nonstationarities or variations induced by the task.

However, as tasks are identified, the prior becomes more informative and the adaptation may have problems to cope with drastic changes such as strongly modifying the signals.

To handle this problem, we chose to limit the size of the prior through the use of a sliding window on past iterations which allows to estimate a moving average of parameters θ.

In practice, we limited D ξt i to the last 250 elements.

Planning

Our algorithm is able to identify a task among a set of possible tasks. To do so it has to explore regions allowing to disambiguate among task hypotheses. There are several efficient model-based reinforcement learning exploration methods that add an exploration bonus for states that might provide more learning gains. Several theoretical results show that these approaches allow to learn tasks efficiently (Brafman and Tennenholtz 2003; [START_REF] Kolter | Near-Bayesian exploration in polynomial time[END_REF]. We define an uncertainty measure and use model-based planning to select sequences of actions that guide the agent towards states that better identify the desired task.

Our algorithm is based on comparing the expected classification rate between different task hypotheses. Therefore, to disambiguate between hypotheses, it has to collect data which are likely to affect differently the dataset quality of competing hypotheses. As we have access to the optimal policies of each task, we can choose for a given task an action in a given state and predict to receive a "correct" or "wrong" label. Such label is linked to a signal generation model (parameterized by θ ML ξt ) which differs from task to task, and to the optimal action at that particular state. A stateaction pair where the optimal actions and the signal models are the same for all hypotheses is unlikely to tell them apart and will be less informative that any other state-action where either optimal actions or signal models differ between hypotheses.

Thus, we define a measure of global uncertainty U (s, a) that is higher when, for a given state-action, there is a high incongruity between either optimal actions or signal models. For this we compute a similarity matrix S where each element S ij (s, a) corresponds to the similarity of the distributions of signals associated to the expected label from tasks i and j if action a is performed in state s. The final uncertainty value U (s, a) is computed as the opposite of the weighted sum of the similarity matrix elements:

U (s, a)=- T X i=1 T X j=1 S ij (s, a)p(ξ i )p(ξ j ) (12)
Computing the similarity between two Gaussian distributions for all state-action pairs was not feasible in real time. In order to improve computation efficiency we do not rely on a precise metric between Gaussian distributions and only consider the similarity between their means (empirical tests did show that this approximation does not impact the results).

Computed for each state and action, this measure is then used as an exploration bonus, whereas we switch to a pure exploitation of the task after reaching the desired confidence level β.

Experimental Protocol and Results

Control task

We consider a 5x5 grid world, where an agent can perform five different discrete actions: move up, down, left, right, or a target-reached action. The user goal is to teach the agent to reach one, yet unknown to the agent, of the 25 discrete states which represent the set of possible tasks (i.e. one task per possible target state). We thus consider that the agent has access to 25 different task hypotheses. We use Markov Decision Processes (MDP) to represent the problem [START_REF] Sutton | Reinforcement learning: An introduction[END_REF]. From a given task ξ, represented as a reward function, we can compute the corresponding policy π ξ using, for instance, Value Iteration [START_REF] Sutton | Reinforcement learning: An introduction[END_REF].

EEG-based feedback signals

EEG signals were recorded with a gTec system (2 gUS-Bamp amplifiers) with 32 electrodes distributed according to the 10/10 international system, with the ground on FPz and the reference on the left earlobe. The EEG signals were digitized with a sampling frequency of 256 Hz, commonaverage-reference (CAR) filtered and band-pass filtered at [0.5, 10] Hz.

During operation, the role of the users was to mentally assess the agent's actions as correct or wrong with respect to a selected target, obtaining this way error-related potentials. Previous studies have demonstrated that these signals can be detected online [START_REF] Ferrez | Error-related EEG potentials generated during simulated Brain-Computer interaction[END_REF]). Following these studies, features were extracted from two frontocentral channels (FCz and Cz) within a time window of [200,700] ms (0 ms being the action onset of the agent) and downsampled to 32 Hz. This leaded to a vector of 34 features that was the input for our system.

Calibration-Free Online BCI Control

This experiment evaluates if we can identify the task desired by the user even without an explicit calibration phase and without prior knowledge of the brain signals. The experiments were conducted with four subjects (aged between 25 and 28). Each subject was asked to mentally assess the agent's actions with respect to a given target. The system was not calibrated to decode the user EEG signals beforehand. Each subject performed 5 runs, for each run a new target was randomly selected and provided to the user. There was an action every three seconds. Each run lasted 200 actions, and the time between runs was around one minute.

The algorithm was able to identify the correct target for all runs of all the subjects, see Figure 2. There are strong variations among subjects but we note that our system identified each task in less iterations than a normal calibration phase requires (between 300 and 600 examples depending on the user performance [START_REF] Chavarriaga | Learning from EEG error-related potentials in noninvasive brain-computer interfaces[END_REF]Iturrate, Montesano, and Minguez 2010)).

Figure 2: Results from the online experiment: Evolution of the probability of the correct task for each subject and run. The algorithm was able to identify the correct target for each subjects and runs in less than 200 iterations.

Table 1 shows for each subject and run the number of iterations needed to reach the confidence threshold for the subject selected target. On average, the number of iterations needed to identify the target was of 85 ± 32. 1: Results from the online experiment: Number of iterations needed to identify the correct target for each subject and run. On average, the number of iterations needed to identify the target was of 85 ± 32.

Algorithm Offline Analysis

The objective of the offline analysis is to study the impact of our exploration method and evaluate if the classifier learned from scratch with our algorithm can be reused for learning new tasks. Finally we want to evaluate how robust the system is to abrupt changes in the signal properties. For these experiments, to ensure we have sufficient data to achieve statistically significant results, we rely on a large dataset of real EEG data. We used a dataset from (Iturrate, Montesano, and Minguez 2013b), which covers ten subjects that performed two different control problems (denoted T 1 and T 2). The role of the users was similar (assess the agent's actions), but the problems differs in the state-action space size and visual representations. For each subject, T 1 was composed of 1800 assessments, and T 2 of 1200. Despite the fact that both problems elicit error-related potentials, the EEG signals presented significant differences (Iturrate, Montesano, and Minguez 2013b).

For each subject, and each dataset (T 1 and T 2), we simulated 20 runs of 400 iterations following the control task. Each time the device performed an action, we sampled the dataset using the ground truth labels corresponding to the correct task and then removed the chosen signal from it. After a first task was identified, and following our approach, we continued running the system to identify new tasks.

We present most of the results in terms of the quality of the dataset, measured as the classification accuracy that a calibrated brain signal classifier would obtain. Results vary strongly between subjects and we will see that it is a direct consequence of the difficulty of finding a classifier with high accuracy.

Planning Methods We compared the average number of steps (with maximum values of 400 steps) needed to identify the first task when learning from scratch with different planning methods.

Uncertainty on task Greedy

Figure 3: Comparison of different exploration methods. Our proposed method, based on the uncertainty on the task and the signals interpretation, allows to lead the system to regions that improve disambiguation among hypotheses in a faster way. For the greedy method, all values were 400 which indicates it never allowed to identify any task.

Figure 3 shows the results averaged across subjects, runs and datasets. Values of 400 means the confidence threshold was not reached after 400 iterations. Our proposed method, based on the uncertainty on the task and the signals interpretation, allows to lead the system to regions that improve disambiguation among hypotheses in a faster way. Trying to follow the most probable task does not allow the system to explore sufficiently (Greedy), and at least some random exploration is necessary to allow a correct identification of the task (ε-greedy). Assessing uncertainty only on the task performs poorly as it does not take into account the signal interpretation ambiguity inherent to our problem. The large variability in the results is mainly due to the large variations in classification accuracy across subjects and datasets. Given these results, the remainder of this section will only consider our proposed planning method.

Online re-estimation of classifier After identifying the first task, and following our approach, we continued running the system and measured how many tasks were identified after 400 steps. The quality of our unsupervised method can be measured according to the percentage of labels correctly assigned (according to the ground truth label), see Figure 4. In general, having dataset with classification accuracies higher than 75% guaranteed that more than 90% of the labels were correctly assigned. This result shows that our algorithm can also be used to collect training data for calibrating any other state-of-the-art error-related potentials classifier, but has the important advantage of controlling the device at the same time. In general, having dataset with classification accuracies higher than 75% guaranteed that more than 90% of the labels were correctly assigned.

Figure 5 demonstrates the advantage of switching to a Bayes filter method after identification of a first target instead of keeping the estimation given by the Bhattacharyya coefficient. On the one hand, Bhattacharyya coefficient works very well for small amounts of data because it directly compares model parameters. On the other hand, when there is sufficient data, training a classifier allows for a faster identification since the classifier makes a much harder decision when evaluating a new EEG signal.

Figure 6 shows the number of tasks correctly and incorrectly identified in 400 iterations. For datasets of good qualities, we are able to identify more than 20 tasks in 400 iterations without the need for a calibration procedure (recap that previous works needed between 300 and 600 examples for the calibration phase [START_REF] Chavarriaga | Learning from EEG error-related potentials in noninvasive brain-computer interfaces[END_REF]Iturrate, Montesano, and Minguez 2010)). The number of correctly identified tasks is strongly correlated to the quality of the dataset. Figure 6: Number of targets correctly and incorrectly identified in 400 iterations (the markers show the median values and the error bars the 2.5th and 97.5th percentiles). For datasets of good qualities, we are able to identify more than 20 tasks in 400 iterations without the need for a calibration procedure.

Robustness to Abrupt Changes in the Signals' Properties We now want to determine whether the system is robust to changes in the signals' properties that occur when we change between different problem settings (Iturrate, Montesano, and Minguez 2013b). We modeled this by changing from Ti to Tj (i 6 = j) after 400 steps, and then executing 400 more steps from Tj. Both combinations (T 1 to T 2 and T 2 to T 1) were tested.

Figure 7 shows the number of tasks identified depending on the classification accuracy on Tj when training a classifier from dataset Ti. For abrupt changes which conserve a classification accuracy above 70% on the new signals, our method, based on a limited size prior, is able to recover and solve more than 20 tasks in 400 iterations. For those cases where the accuracy change is too drastic, starting from scratch may be a better solution than relying on the adaptation properties of our algorithm.

Conclusion

We introduced a novel method for calibration-free BCI based control of sequential tasks with feedback signals. The method provides an unsupervised way to train a decoder with almost the same performance as state-of-the-art supervised classifiers, while keeping the system operational and solving the task requested by the user since the beginning. The intuition for our method is that the classification of the brain signals is easier when they are interpreted according to the task desired by the user. The method assumes a distribution of possible tasks and relies on finding which pair of decoder-task has the highest expected classification rate on the brain signals.

The algorithm was tested with real online experiments, showing that the users were able to guide an agent to a desired position by mentally assessing the agent's actions and without any explicit calibration phase. Offline experiments show that we can identify an average of 20 tasks in 400 iterations without any calibration, while in previous works the calibration phase used between 300 and 600 examples. To improve the efficiency of the algorithm, we introduced a new planning method that uses the uncertainty in the decodertask estimation. Finally, we analyzed the performance of the system in the presence of abrupt changes in the EEG signals.

Our proposed method was able to adapt and reuse its learned models to the new signals. Furthermore, in those cases when the transfer is not possible, our method can still be used to recalibrate the system from scratch while solving the task.

A current limitation of the work is the need for a finite set of task hypotheses. This limitation could be solved by the use of a combination of particle filter and regularization on the task space. Additionally, our method can not dissociate fully symmetric hypotheses, e.g. right and left most state of our 1D grid world (Fig. 1), as the interpretation of feedback signals will also be symmetric and therefore as likely. This latter problem can be solved by redefining the set of hypotheses or the action set, for instance by adding a "stop" action valid only at the target state.

This work opens a new perspective regarding the global challenge of interacting with machines. It has application to many interaction problems which requires a machine to learn how to interpret unknown communicative signals. A promising avenue, outside the BCI field, lies in human robot interaction scenarios where robots must learn from, and interact with, many different users who have their own limitations and preferences.

Figure 1 :

 1 Figure 1: Representation of inferred signal labels for a 1D grid world in function of the hypothesized task. Three task hypotheses are displayed in column from left to right [(a), (b), (c)].On top is a 1D grid world with the hypothetic target state marked with a T letter. Notice that this hypothetic target is different for each of the three hypotheses. The arrows in each cell indicates what action should elicit a positive feedback, i.e. the optimal policy with respect to the hypothetic target. The user intended target is shown as the shaded blue state at the right extremity of the 1D world. Note that the system does not have access to this information. The correct hypothesis is the one on the Left [(a)] where the T state is the same as the shaded blue state. Below the 1D grid world, the signals received from the user are represented in a 2D feature space. They represent the user assessment signals of the past history of device's actions (e.g. moving randomly left and right). Our algorithm assigns virtual labels (green for "correct" and red for "wrong") to those signals with respect to their respective hypothetic target. Notice that the data points are the same for all three hypotheses, only their respective labels differ. With the virtual labels being assigned, we can compute the corresponding 2D Gaussian distributions estimates for each class (shown as colored ellipses) and each hypothesis. While for the correct hypothesis [(a)] the Gaussian distributions shows a large separability, the overlap increases as the hypothetic target (T) moves away from the real (blue shaded) one[(b), (c)]. This property can be exploited to estimate the correct hypothesis and the model generating the signals.

Figure 4 :

 4 Figure4: Percentage of labels correctly assigned according to the ground truth label (the markers show the median values and the error bars the 2.5th and 97.5th percentiles). In general, having dataset with classification accuracies higher than 75% guaranteed that more than 90% of the labels were correctly assigned.

Figure 5 :

 5 Figure5: Number of targets correctly identified in 400 iterations (the markers show the median values and the error bars the 2.5th and 97.5th percentiles). Comparison between switching to a Bayes filter method after identification of a first target instead of keeping the estimation given by the Bhattacharyya coefficient. The Bayes filter allows for a faster identification.

Figure 7 :

 7 Figure7: Number of targets correctly and incorrectly identified in 400 iterations after an abrupt change in the signals properties (the markers show the median values and the error bars the 2.5th and 97.5th percentiles). For abrupt changes which conserve a classification accuracy above 70% on the new signals, our method, based on a limited size prior, is able to recover and solve more than 20 tasks in 400 iterations.

Acknowledgments

The authors would like to thank the anonymous reviewers for their helpful comments which contributed in improving the quality of this paper. The authors from INRIA are with the Flowers Team, a join INRIA -Ensta ParisTech lab, France. Luis Montesano is with Instituto de Investigación en Ingeniería de Sistemas, Universidad de Zaragoza, Spain. Iñaki Iturrate is with the Chair in Non-invasive Brain-Machine Interface (CNBI) and Center for Neuroprosthetics and Institute of Bioengineering, EPFL, Switzerland. This work has been partially supported by INRIA, Conseil Régional d'Aquitaine and the ERC grant EXPLORERS 24007; and from the Spanish Ministry via DPI2011-25892 (L.M.), and DGA-FSE grupo T04.