Emilio Jorge
email: emilio.jorge@chalmers.se

Hannes Eriksson

Christos Dimitrakakis

Debabrota Basu

Divya Grover

Inferential Induction: A Novel Framework for Bayesian Reinforcement Learning

Bayesian Reinforcement Learning (BRL) offers a decision-theoretic solution to the reinforcement learning problem. While "model-based" BRL algorithms have focused either on maintaining a posterior distribution on models, BRL "modelfree" methods try to estimate value function distributions but make strong implicit assumptions or approximations. We describe a novel Bayesian framework, inferential induction, for correctly inferring value function distributions from data, which leads to a new family of BRL algorithms. We design an algorithm, Bayesian Backwards Induction (BBI), with this framework. We experimentally demonstrate that BBI is competitive with the state of the art. However, its advantage relative to existing BRL model-free methods is not as great as we have expected, particularly when the additional computational burden is taken into account.

Introduction

Many Reinforcement Learning (RL) algorithms are grounded on the application of dynamic programming to a Markov Decision Process (MDP) [START_REF] Richard | Reinforcement learning: An introduction[END_REF]. When the underlying MDP µ is known, efficient algorithms for finding an optimal policy exist that exploit the Markov property by calculating value functions. Such algorithms can be applied to RL, where the learning agent simultaneously acts in and learns about the MDP, through e.g. stochastic approximations, without explicitly reasoning about the underlying MDP. Hence, these algorithms are called modelfree.

In Bayesian RL (BRL) [START_REF] Ghavamzadeh | Bayesian reinforcement learning: A survey[END_REF], we represent our knowledge about the underlying MDP µ through some prior distribution β over a set M of possible MDPs. This explicit representation of uncertainty admits algorithms that can perform near-optimal exploration, in contrast to methods that only rely on a single empirical model of the MDP, which require exploration heuristics. [c.f. Vlassis et al., 2012] While model-based BRL is well-understood, many BRL algorithms try to become model-free by calculating distributions on value functions. Unfortunately, these methods typically make implicit assumptions or approximations about the underlying MDP distribution.

In the rest of this section, we provide background in terms of setting and related work. In Section 2, we explain our Inferential Induction framework and three different inference methods that emerge from it, before instantiating one of them into a concrete procedure. Based on this, Section 3 describes the BBI algorithm. In Section 4, we experimentally compare BBI with state-of-the art BRL algorithms. In the appendix additional experimental results are available in Appendix B and some implementation details in Appendix C.

Setting and Notation

In this paper, we generally use P and E to refer to probability (measures) and expectations while allowing some abuse of notation for compactness.

Reinforcement Learning (RL) is a sequential learning problem faced by agents acting in an unknown environment µ, typically modelled as a Markov decision process [c.f. [START_REF] Puterman | Markov Decision Processes : Discrete Stochastic Dynamic Programming[END_REF]. Definition 1.1 (Markov Decision Process (MDP)). An MDP µ with state space S and action space A is equipped with a reward distribution P µ (r | s) with corresponding expectation ρ µ (s) and a transition kernel P µ (s ′ |s, a) for states s, s ′ ∈ S and actions a ∈ A.

At time t, the agent observes1 the environment state s t , and then selects an action a t . Then, it receives and observes a reward r t and a next state s t+1 . The agent is interested in the utility U t ≜ T k=t γ k-t r k , i.e. the sum of future rewards r t . Here, γ ∈ (0, 1] is the discount factor and T ∈ [1, ∞] is the problem horizon. Typically, the agent wishes to maximise the expected utility, but other objectives are possible.

The agent acts in the environment using a policy π = (π 1 , . . . , π t , . . .) that takes an action a t at time t with probability π t (a t | s t , r t-1 , a t-1 , s t-1 , . . . , r 1 , a 1 , s 1 ). Dependence on the complete observation history is necessary, if the agent is learning from experience. However, when µ is known, the policy π * µ maximising expected utility over finite horizon is Markovian2 of the form π t (a t | s t ) and is computable using dynamic programming. A useful algorithmic tool for achieving this is the value function, i.e. the expected utility of a policy π from different starting states and action: Definition 1.2 (Value Function). The state value function of policy

π in MDP µ is V π µ,t (s) ≜ E π µ (U t | s t = s) and the corresponding state-action (or Q-)value function is Q π µ,t (s, a) ≜ E π µ (U t | s t = s, a t = a). P π
µ and E π µ denote probabilities and expectations under the process induced by π and µ.

Finally, the Bellman operator

B π µ V (s) ≜ ρ µ (s) + γ s ′ ∈S P π µ (s ′ | s)V (s ′ ) allows us to compute the value function recursively through V π µ,t = B π µ V π µ,t+1 . 3
Bayesian RL (BRL). In BRL, our subjective belief is represented as a probability measure β over possible MDPs. We refer to the initial belief β as the prior distribution. By interacting with the environment until time t, the agent obtains data D t = (s 1 , a 1 , r 1 , . . . , s t ). This data is used to calculate a posterior distribution β(µ | D t ) that represents agent's current knowledge about the MDP. 4 For a given belief and an adaptive policyπ β (s), we define the Bayesian value function to be:

V π β β,t (s) ≜ M V π β µ,t (s) dβ(µ). (1)
The Bayesian value function is the expected value function under the distribution β. The Bayesoptimal policy achieves the Bayes-optimal value function V

* β,t (s) = sup π V π β,t . Calculating V π β,t
involves integrating V π µ,t for all µ, while V * β,t typically requires exponential time. Information about the value function distribution can be a useful tool for constructing near-optimal policies, as well a way to compute risk-sensitive policies.

Related Work and Our Contribution

We arrange related work in the following areas. Firstly, "model-free" methods that do not explicitly try to take into account the uncertainty about the underlying MDP. Secondly, "model-based" methods, that explicitly maintain a distribution on MDP models. Finally, we discuss approximations to the Bayes-optimal solution, before we outline our contributions.

Model-free Bayesian Value Functions. Bayesian value function distributions have been considered extensively in model-free Bayesian Reinforcement Learning (BRL). One of the first methods was Bayesian Q-learning [START_REF] Dearden | Bayesian Q-learning[END_REF], which used a normal-gamma prior on the utility distribution. However, as i.i.d. utility samples cannot be obtained by bootstrapping from value function estimates, this idea had inherent flaws. [START_REF] Engel | Bayes meets Bellman: The Gaussian process approach to temporal difference learning[END_REF] developed a more sophisticated approach, the Gaussian Process Temporal Difference (GPTD) algorithm, which has a Gaussian process (GP) prior β(V ) on value functions. It then combines this with the likelihood function

P(D | V ) ∝ t i=1 exp{-|V (s i ) -r i -γV (s i+1 )| 2 }.
However, this makes the implicit assumption that the deterministic empirical MDP model is correct. [START_REF] Engel | Reinforcement learning with Gaussian process[END_REF] tried to relax this assumption by allowing for correlation between sequentially visited states. [START_REF] Deisenroth | Gaussian process dynamic programming[END_REF] developed a dynamic programming algorithm with a GP prior on value functions and an explicit GP model of the MDP. Finally, Tang and Agrawal [2018] introduced VDQN, generalising such methods to Bayesian neural networks. The assumptions that these model-free Bayesian methods implicitly make about the MDP are hard to interpret, and we find the use of an MDP model independently of the value function distribution unsatisfactory. We argue that explicitly reasoning about the joint value function and MDP distribution is necessary to obtain a coherent Bayesian procedure. Unlike the above methods, we calculate a value function posterior P(V |D) while simultaneously taking into account uncertainty about the MDP.

Model-based Bayesian Value Functions. If a posterior over MDPs is available, we can calculate a distribution over value functions in two steps: a) sample from the MDP posterior and b) calculate the value function of each MDP. Essentially, this amounts to performing posterior sampling (PSRL, [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF], Thompson [1933]) followed by policy evaluation. [START_REF] Dearden | Model based Bayesian exploration[END_REF] suggested an early version of this approach that obtained approximate upper bounds on the Bayesian value function and sketched a Bellman-style update for performing it online. Posterior sampling approach was later used to obtain value function distributions in the discrete case by [START_REF] Dimitrakakis | Robust Bayesian reinforcement learning through tight lower bounds[END_REF] and in the continuous case by [START_REF] Osband | Generalization and exploration via randomized value functions[END_REF]. Finally, O'Donoghue et al. [2018] derive bounds on the variance of the value function posterior. We instead focus on whether it is possible to compute complete value function distributions exactly or approximately through a backwards induction procedure. In particular, how can we obtain P(V i |D) from P(V i+1 |D)?

Our Contribution. We introduce Inferential Induction, a new Bayesian Reinforcement Learning (BRL) framework, which leads to a Bayesian form of backwards induction. Our framework allows Bayesian inference over value functions without any implicit assumption or approximation unlike its predecessors. The main idea is to calculate the conditional value function distribution at step i from the value function distribution at step i + 1 analogouslto backwards induction for the expectation (Eq. ( 2)). Following this, we propose a simple marginalisation techniques and design an appropriate Monte Carlo approximation for it. We can combine this procedure with a policy optimisation mechanism. We use a Bayesian adaptation of dynamic programming for this and propose the Bayesian backwards induction (BBI) algorithm. Our experimental evaluation shows that BBI is competitive to the current state of the art. Inferential Induction framework provides the opportunity to further design more efficient algorithms of this family.

Inferential Induction

Given a prior belief β, and data D t = (s 1 , a 1 , r 1 , . . . , s t ) obtained by interaction of the agent with the MDP, we wish to calculate the value function distribution from step t onwards, i.e. P π β (V t | D t ) for a policy π under the belief β, conditioned on the data D t . The main idea for doing this to inductively calculate

P π β (V i | D t ) from P π β (V i+1 | D t ) for i ∈ {t, .
. . , T } as follows:

P π β (V i | D t ) = V P π β (V i | V i+1 , D t ) d P π β (V i+1 | D t ). (2)
Let ψ i+1 be a (possibly approximate) representation of P π β (V i+1 | D t ). If we can calculate the above integral, then we can also obtain ψ i ≈ P π β (V i | D t ) recursively, from time T up to the current time step t. Then the estimation problem reduces to defining the term P π β (V i | V i+1 , D t ) appropriately. In this paper, we describe a simple Monte Carlo method for solving the estimation problem and an approximate dynamic programming algorithm for optimising the policy within this framework.

Integrating over P π β (µ | V i+1 , D t ).
A simple idea for dealing with the term linking the two value functions is to directly marginalise over the MDP as follows:

P π β (V i | V i+1 , D t ) = M P π µ (V i | V i+1 ) d P π β (µ | V i+1 , D t ). (3)
This equality holds because given µ, V i is uniquely determined by the policy π and V i+1 through the Bellman operator. However, it is crucial to note that

P π β (µ | V i+1 , D t ) ̸ = P β (µ | D t )
, as knowing the value function gives information about the MDP.5 This is a crucial difference with "model-free" Bayesian value function methods, which effectively hide strong assumptions about the MDP when they perform value function inference. We expect that this approach would give superior results. The remaining computations are rather straightforward Monte Carlo approximations.

The backwards induction step. To calculate the previous-step distribution from next-step distribution, let us now combine the induction step in (2) with the marginalisation in (3). We also substitute an approximate representation ψ i+1 for the next-step belief P(V i+1 | D t ), to obtain the following conditional probability measure on value functions:

ψ i (B) ≜ P π β (V i ∈ B|D t ) = V M 1 B π µ V i+1 ∈ B d P π β (µ|V i+1 , D t ) dψ i+1 (V i+1 )
Monte-Carlo approximation of the inner integral. We now have to leave generalities and commit to some hard choices for defining and calculating

P π β (µ | V i+1 , D t ).
Expanding this term we obtain, for subsets of MDPs A ⊆ M, the following measure:

P π β (µ ∈ A | V i+1 , D t ) = A P π µ (V i+1 ) dβ(µ | D t ) M P π µ (V i+1 ) dβ(µ | D t ) , ( 4
) since P π µ (V i+1 | D t ) = P π µ (V i+1
), as µ, π are sufficient for calculating V i+1 . To compute P π µ (V i+1 ), we can sample rollouts from an arbitrary starting state distribution q and the policy π for each sampled MDP µ and then marginalise over the resulting utility values U . In exact form, this is done through the integral:

P π µ (V i+1 ) = S dq(s) ∞ -∞ P π µ (V i+1 | U, s) P π µ (U |s) dU.
In order to understand the meaning of the term

P π µ (V i+1 | U, s), note that V π i+1,µ (s) = E π µ [U | s i+1 = s].
Thus, a rollout u from state s gives us partial information about the value function. Finally, the starting state distribution q is used to measure the goodness-of-fit, similarly to e.g. fitted Q-iteration6 .

As a design choice, we define the density of V i+1 given a rollout sample u m in MDP µ from state s m ∼ q to be a Gaussian with variance σ 2 . If we generate N µ number of MDPs µ (j) ∼ β(µ | D t ) and set:

w jk ≜ n m=1 e -| V (k) i+1 (sm )-u j m| 2 2σ 2 Nµ j ′ =1 n m=1 e - V (k) i+1 (sm )-u j ′ m 2 2σ 2 , ( 5
) we get E[w jk ] = P π β (µ ∈ M | V i+1 , D t ).
The weight w jk can be interpreted as a measure of how well the value function V (k) i+1 matches the rollouts obtained from µ (j) . This allows us to obtain value function samples for step i, V

(j,k) i ≜ B π µ (j) V (k) i+1 , (6)
each weighted by w jk , leading to the following Monte Carlo estimate of the value function distribution at step i

ψ i (B) = 1 N V N µ N V k=1 Nµ j=1 1 V (j,k) i ∈ B w jk . (7)
Here, N V is the number value function samples

V (k)
i+1 . This ends the general description of the Monte Carlo method. Detailed design of an algorithm depends on the representation that we use for ψ i and whether the MDP is discrete or continuous.

Section 3 gives algorithmic details, while specifications of hyperparameters and distributions are given in Section 4.

Algorithms

Algorithm 1 is a concise description of the Monte Carlo procedure that we develop. At each time step t, the algorithm is called with the prior and data D collected so far, and it looks ahead up to some lookahead factor H. 7 We instantiate it below for discrete and continuous state spaces.

Algorithm 1 Policy Evaluation with Method 1

1: Input: Prior β, data D, lookahead H, discount γ, policy π, N µ , N V . 2: Initialise ψ H . 3: Sample M ≜ µ (j) j ∈ [N µ ] from β(µ | D). 4: for i = H -1, . . . , 1 do 5: Sample V (k) ∼ ψ i+1 (v) for k ∈ [N V ]. 6:
Generate n utility samples u m 7:

Calculate w jk from (5) and V (j,k) i from (6).

8:

Calculate ψ i from (7). 9: end for 10: return {ψ i | i = 1, . . . , H} Discrete MDPs. When the MDPs are discrete, the algorithm is straightforward. Then the belief β(µ | D) admits a conjugate prior in the form of a Dirichlet-product for the transitions. In that case, it is also possible to use a histogram representation for ψ i , so that it can be calculated by simply adding weights to bins according to (7). However, as a histogram representation is not convenient for a large number of states, we model using a Gaussian ψ t . In order to do this, we use the sample mean and covariance of the weighted value function samples V (j,k) i :

m i = 1 N V N µ N V k=1 Nµ j=1 V (j,k) i w jk , Σ i = 1 N V N µ N V k=1 Nµ j=1 (V (j,k) i -m i )(V (j,k) i -m i ) ⊤ w jk . (8) such that ψ i = N (m i , Σ i ) is a multivariate normal distribution.
Continuous MDPs. In the continuous state case, we obtain ψ through fitted Q-iteration [c.f. [START_REF] Ernst | Tree-based batch mode reinforcement learning[END_REF]. For each action a in a finite set, we fit a weighted linear model

Q i (s i , a) = s T i ω a + ϵ a , where s i , ω a ∈ R d and ϵ a ∼ N (0, σ 2 a )
. Thus, we obtain the sample mean and covariance to be:

m a i = 1 N V N µ N V k=1 Nµ j=1 (s (j,k) i ) T ω a , Σ a i = σ 2 a N V N µ N V k=1 Nµ j=1 (Q (j,k) i (a) -m a i )(Q (j,k) i (a) -m a i ) ⊤ w jk .
(9) In practice, we often use a feature map ϕ : R d → R f for states.

Bayesian Backwards Induction

We now construct a policy optimisation component to use with the inferential induction based policy evaluation and the aforementioned two approximation techniques. We use a dynamic programming algorithm that looks ahead H steps, and at each step i calculates a policy maximising the Bayesian expected utility in the next i + 1 steps. We describe the corresponding pseudocode in Algorithm 2.

Algorithm 2 Line 8: Discrete MDPs. Just as in standard backwards induction, at each step, we can calculate π i by keeping π i+1 , . . . , π H fixed:

Q i (s, a) ≜ E β (U | s i = s, a i = a, D) = M ρ µ (s, a) + s ′ P (j) µ (s ′ |s, a)V πi+1,...,π H µ,i (s ′ ) ≈ j,k ρ µ (j) (s, a) + s ′ P (j) µ (s ′ |s, a)V (k) i+1 (s ′ ) w jk N µ N V . ( 10
)
Algorithm 2 Line 8: Continuous MDPs. As we are using fitted Q-iteration, we can directly use the state-action value estimates. So we simply set

Q i (s, a) = Qi (s, a).
The Q i estimate is then used to select actions for every state. We set π i (a|s) = 1 for a = arg max Q i (s, a) (Line 2.9) and calculate the value function distribution (Lines 2.10 and 2.11) for the partial policy (π i , π i+1 , . . . , π H ).

Algorithm 2 Bayesian Backwards Induction (BBI) with Method 1

1: Input: Prior β, data D, lookahead H, discount γ, N µ , N V . 2: Initialise ψ i . 3: Sample M ≜ µ (j) j ∈ [N µ ] from β(µ | D). 4: for i = H -1, . . . , 1 do 5: Sample V (k) ∼ ψ i+1 (v) for k ∈ [N V ]. 6:
Generate n utility samples u i 7:

Calculate w jk from (5).

8:

Calculate Q i from (10) or fitted Q-iteration. 9:

Set π i (a|s) = 1 for a ∈ arg max Q i (s, a).

10:

Calculate w jk from (5) and V (j,k) i from (6) with policy π i :

V (j,k) i ≜ B πi µ (j) V (k)
i+1 .

11:

Calculate ψ i from (8) or (9). 12: end for 13: return π = (π 1 , . . . , π H ).

Experimental Analysis

For performance evaluation, we compare Bayesian Backwards Induction (BBI, Alg. 2) with exploration by distributional reinforcement learning [VDQN, Tang and Agrawal, 2018]. We also compare BBI with posterior sampling [PSRL, [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF], Thompson, 1933] for the discrete MDPs and with Gaussian process temporal difference [GPTD, [START_REF] Engel | Bayes meets Bellman: The Gaussian process approach to temporal difference learning[END_REF] for the continuous MDPs. First, we describe the experimental setup and the priors used for implementation. In Section 4.1, we analyse the results obtained for different environments in terms of average reward obtained over time. Comparisons with Multi-MDP Backwards Induction [Dimitrakakis, 2011, MMBI], Bayesian Sparse Sampling [Wang et al., 2005, BSS] and Bayesian Q-Learning [Dearden et al., 1998, BQL] are available in the appendix.

Parameters. We run the algorithms for the infinite-horizon formulation of value function with discount factor γ = 0.99. Each algorithm updates its policy at steps t k = k(k+1)

2

, where k increases as 1, 2, . . .. This sequence of t k 's causes a total of 1413 updates in 10 6 steps. We set the lookahead H to 100 and 20 for discrete and continuous MDPs respectively. More implementation details can be found in the appendix.

Prior. For discrete MDPs, we use Dirichlet Dir(α) priors over each of the transition probabilities P(s ′ |s, a). The prior parameter α for each transition is set to 0.5. We use separate NormalGamma N G(µ, κ, α, β) priors for each of the reward distributions P(r|s, a). We set the prior parameters to

[µ 0 , κ 0 , α 0 , β 0 ] = [0, 1, 1, 1].
For continuous MDPs, we use factored Bayesian Multivariate Regression [START_REF] Thomas | Bayesian linear regression[END_REF] models as priors over transition kernels and reward functions for the continuous environments. This implies that the transition kernel P(s ′ |s, a) and reward kernel P(r|s, a) modelled as N (A Trans a s, Σ) and N (A Reward a s, σ 2 ). Σ is sampled from inverse Wishart distribution with corresponding d × d dimensional scale matrix, while σ is sampled from inverse Gamma with prior parameters ( 1 2 , 1 2 ). For transitions, we set the prior parameters to Ψ 0 = 0.001I and degrees of freedom ν 0 = rank(Ψ 0 ).

Feature Map. For InvertedPendulum we use the 10 dimensional feature map ϕ : R d → R f used by [START_REF] Lagoudakis | Least-squares policy iteration[END_REF] for BBI and GPTD such that the modelling becomes P(s ′ | ϕ(s), a) and P(r | ϕ(s), a). VDQN works directly on the underlying state s as it uses a neural network. We also add a regularizing term, λI, λ = 0.01 to the diagonal of the fitted Q-iteration fit. The choice of state distribution q(s) is of utmost importance in continuous environments. In the continuous environment, we experimented with a few options, trading of sampling states from our history, sampling from the starting configuration of the environment and sampling from the full support of the state space.

Description of Environments

We evaluate the algorithms on four discrete (NChain [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF], DoubleLoop [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF], LavaLake [START_REF] Leike | Ai safety gridworlds[END_REF] and Maze [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF]) and one continuous environment (InvertedPendulum [START_REF] Lagoudakis | Least-squares policy iteration[END_REF]). See the appendix for more details.

Experimental Results

The following experiments are intended to show that the general methodological idea is indeed sound, and can potentially lead to high performance algorithms. More results, comparative experiments, as well as tests of convergence, are provided in the appendix.

Figures 1a,1b,1c, 1d illustrate the evolution of average reward for BBI, PSRL and VDQN on the discrete MDPs.8 BBI is competitive with PSRL, which has good exploration properties, while VDQN generally performs worse. Figure 1e shows a comparison with state-of-the art algorithms, such as VDQN and GPTD. Our algorithm is competitive, and in particular performs much better than GPTD, while it performs similarly to VDQN, which is slightly worse initially and slightly better later in terms of average steps survived. This performance could partially be explained by the use of a linear value function Q(ϕ(s), a), in contrast to VDQN which uses a neural network.

Comparison of Computation Times

In Table 1 we can see the timing results for the different algorithms when run on Intel(R) Core(TM) i9-9900K CPU @ 3.60GHz. Here we can see that PSRL is faster and scales better due to the simplicity of using a single MDP. It is probable that BBI could obtain speedups through improve sampling

(V i | V i+1 , D
) using the data and to apply it inductively for computing the marginal value function distribution P π β (V i | D). We then designed a straightforward Monte Carlo approximation and combined it with a suitable policy optimisation mechanism and showed that it can be competitive with the state of the art.

In order to place it in context, standard backwards induction (i.e. value iteration) calculates

V i = B π µ V i+1 and distributional reinforcement learning methods calculate P π µ (U i | U i+1 )
, both for a given, underlying MDP. In an RL setting, this can be replaced either through stochastic approximation, such as Q-learning, or through an explicit empirical model. Inferential Induction differs from existing Bayesian value function methods, which essentially cast the problem into regression. For example, GPTD [START_REF] Engel | Bayes meets Bellman: The Gaussian process approach to temporal difference learning[END_REF]] can be written as Bayesian inference with a GP prior over value functions and a data likelihood that uses a deterministic empirical model of the MDP. While this can be relaxed by using temporal correlations as in [START_REF] Engel | Reinforcement learning with Gaussian process[END_REF], the fundamental problem remains. Even though such methods have practical value, we show that Bayesian estimation of value functions requires us to explicitly think about the MDP distribution as well.

Algorithm Design. We use specific approximations for discrete and continuous MDPs to propose the Bayesian Backwards Induction (BBI) algorithm. The algorithm we developed from this family appears promising, as we are able to outperform methods that implicitly use the empirical model of the MDP, such as GPTD. This shows that the inference is inherently sound. We see that BBI is also competitive with state-of-the art methods like PSRL, and it generally outperforms algorithms relying on approximate inference, such as VDQN.

Limitations. We have observed that the particular method, BBI, that we have developed requires a significant number of samples for it to obtain a good performance, which makes it inherently slow. Our timing analysis in Table 1 shows that methods like PSRL may work significantly faster as they have to keep track of only one MDP and are also simpler to implement. It is probable that BBI could obtain speedups from using sparse posteriors or sampling from state space in the importance sampling step of the algorithm.

For the continuous environments, we have used weighted linear model and linear priors for updating and tracking the value function distribution. Thus, we rely on pre-defined features (e.g. the features from [START_REF] Lagoudakis | Least-squares policy iteration[END_REF]) which are problem specific, not always efficient, and can be replaced more efficient functional approximators, such as neural networks (similar to Tang and Agrawal [2018]).

Since our aim has been to propose the framework, we have emphasised one policy evaluation and used simpler methods, such as Q-update and fitted Q-iteration, for policy update. It is possible to incorporate this method with policy gradient methods to design more efficient algorithms. This point shows the flexibility of our approach as a framework and weakness as a specifically designed algorithm, which is a topic of future investigation.

A Description of Environments

NChain. This is a discrete stochastic MDP with 5 states, 2 actions [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF]. Taking the first action returns a reward 2 for all states and transitioning to the first state. Taking the second action returns 0 reward in the first four states (and the state increases by one) but returns 10 for the fifth state and the state remains unchained. There is a probability of slipping of 0.2 with which its action has the opposite effect. This environment requires both exploration and planning to be solved effectively and thus acts as an evaluator of posterior estimation, efficient performance and effective exploration.

DoubleLoop. This is a slightly more complex discrete deterministic MDP with with two loops of states [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF]. Taking the first action yields traversal of the right loop and a reward 1 for every 5 state traversal. Taking the second action yields traversal of the left loop and a reward 2 for every 5 state traversal. This environment acts as an evaluator of efficient performance and effective exploration.

LavaLake. This is a stochastic grid world [START_REF] Leike | Ai safety gridworlds[END_REF] where every state gives a reward of -1, unless you reach the goal, in which case you get 50, or fall into lava, where you get -50. We tested on the 5 × 7 and a 10 × 10 versions of the environment. The agent moves in the direction of the action (up,down,left,right) with probability 0.8 and with probability 0.2 in a direction perpendicular to the action.

Maze. This is a grid world with four actions (ref. Fig. 3 in [START_REF] Strens | A Bayesian framework for reinforcement learning[END_REF]). The agent must obtain 3 flags and reach a goal.There are 3 flags throughout the maze and upon reaching the goal state the agent obtains a reward of 1 for each flag it has collected and the environment is reset. Similar to LavaLake, the agent moves with probability 0.9 in the desired direction and 0.1 in one of the perpendicular directions. The maze has 33 reachable locations and 8 combination of obtained flags for a total of 264 states.

InvertedPendulum. To extend our results for the continuous domain we evaluated our algorithm in a classical environment described in [START_REF] Lagoudakis | Least-squares policy iteration[END_REF]. The goal of the environment is to stabilize a pendulum and to keep it from falling. If the pendulum angle θ falls outside [ -π 2 , π 2 ] then the episode is terminated and the pendulum returned to its starting configuration. The state dimensionality is a tuple of the pendulum angle as well as its angular velocity, θ, s = (θ, θ). The environment is considered to be completed when the pendulum has been kept within the accepted range for 3000 steps. For further environment details, see [START_REF] Lagoudakis | Least-squares policy iteration[END_REF].

B Additional Results

Here we present some experiments that examine the performance of inferential induction in terms of value function estimation, inference and utility obtained. For the latter we present additional results against other algorithms, as well as one new continuous environment.

B.1 Bayesian Value Function Estimation

In this experiment, we evaluate the Bayesian (i.e. mean) value function of the proposed algorithm (BBI) with respect to the upper bound on the Bayes-optimal value function. The upper bound is calculated from M max π V π µ dβ(µ | D). We estimate this bound through 100 MDP samples for NChain. We plot the time evolution of our value function and the simulated Bayes bound in Figure 2 for 10 5 steps. We observe that this is becomes closer to the upper bound as we obtain more data.

B.2 Value Function Distribution Estimation

Here we evaluate whether inferential induction based policy evaluation (Alg. 1) results in a good approximation of the actual value function posterior. In order to evaluate the effectiveness of estimating the value function distribution using inferential induction (Alg. 1), we compare it with the Monte Carlo distribution and the mean MDP. We compare this for posteriors after 10, 100 and 1000 time steps, obtained with a fixed policy in NChain that visits all the states, in Figure 3 for 5 runs of Alg. 1. The fixed policy selects the first action with probability 0.8 and the second action with probability 0.2. The Monte Carlo estimate is done through 1000 samples of the value function vector (γ = 0.99). This shows that the estimate of Alg. 1 reasonably captures the uncertainty in the true distribution. For this data, we also compute the Wasserstein distance [START_REF] Fournier | On the rate of convergence in Wasserstein distance of the empirical measure[END_REF] between the true and the estimated distributions at the different time steps as can be found in Table 2. There we can see that the distance to the true distribution decreases over time. We have also run additional experiments with other Bayesian algorithms. In particular, here we show comparisons with MMBI [START_REF] Dimitrakakis | Robust Bayesian reinforcement learning through tight lower bounds[END_REF], BSS [Wang et al., 2005] and BQL [START_REF] Dearden | Bayesian Q-learning[END_REF].

As can be seen in Figures 4 to 6, BBI performs similarly to to MMBI and PSRL. This is to be expected, as the optimisation algorithm used in MMBI is close in spirit to BBI, with only the inference being different. In particular, this algorithm takes k MDP samples from the posterior, and then performs backward induction in all the MDP simultaneously to obtain a Markov policy. In turn, PSRL can be seen as a special case of MMBI with just one sample. This indicates that the BBI inference procedure is sound. The near-optimal Bayesian approximation performs slightly worse in this setting, perhaps because it was not feasible to increase the planning horizon sufficiently.9 Finally, the less principled approximations, like VDQN and BQL do not manage to have a satisfactory performance in these environments.

B.4 Analysis of Variation in Performance

In Figures 7 to 11, we illustrate the variation in performance of the different algorithms for each environment. The black lines illustrate the standard error and the 5th and 95th percentile performance is highlighted. The results indicate that BSS is the most stable algorithm, followed by BBI, MMBI and PSRL, which nevertheless have better mean performance. VDQN is quite unstable, however.

C Implementation Details

In this section we discuss some additional implementation details, in particular how exactly we performed the rollouts and the selection of some algorithm hyperparameters, as well as some sensitivity analysis.

C.1 Computational Details of Rollouts

To speed up the computation of rollouts, we have used three possible methods that essentially bootstrap previous rollouts or use value function samples:

u µ,πt t (s) = r(s, a) + γu µ,πt+1 t+1 (s ′ ) (11) u µ,πt t (s) = s ′ r(s, a) + γP (s ′ |s, a)u µ,πt+1 t+1 (s ′ ) (12) u µ,πt t (s) = s ′ r(s, a) + γP (s ′ |s, a)V t+1 (s ′ ) (13)
where V t+1 ∼ ψ t+1 . In experiments, we have found no significant difference between them. All results in the paper use the formulation in (13).

C.2 Hyperparameters

For the experiments, we use the following hyperparameters.

We use 10 MDP samples, a planning horizon T of 100, γ = 0.99 and we set the variance of the Gaussian to be σ 2 = V 2 span 10 -4 , where V span is the span of possible values for each environment (obtained assuming maximum and minimum reward). We use Eq. 13 for rollout computation with 10 samples from V t+1 and 50 samples from V t (20 for LavaLake 10 × 10 and Maze). If the weights obtained in ( 5) are numerically unstable we attempt to resample the value functions and then double σ until it works (but is reset to original value when new data is obtained). This is usually only a problem when very little data has been obtained.

Varying Horizon T . In order to check the sensitivity on the choice of horizon T , we perform a sensitivity analysis with T = 10, 20, 50, 100. In Figure 12, we can see that varying the horizon has a very small impact for NChain and Maze environments.

Figure 1 :

 1 Figure 1: Evolution of average reward for NChain, DoubleLoop, LavaLake 10 × 10, Maze and InvertedPendulum. The results are averaged over 30 runs of length 10 5 for the continuous Inverted Pendulum. The discrete environment runs are of length 10 6 (10 6 for Maze) with 50 runs for NChain and DoubleLoop and 30 runs for Maze and LavaLake 10 × 10. The runs are exponentially smoothed with a half-life 1000 and 2500 before averaging for the discrete and continuous runs respectively.

Figure 2 :

 2 Figure 2: Comparisons of the achieved value functions of BBI with the upper bound on Bayesoptimal value functions. Upper bound and BBI are calculated from 100 MDPs and plotted for 10 5 time steps.

Figure 3 :

 3 Figure 3: Comparison of value function posteriors obtained by inferential induction and Monte Carlo evaluations at different time steps for a fixed policy. We plot for five runs of inferential induction at each time step. The value of the mean MDP is shown by a vertical line.

Figure 4 :

 4 Figure4: Evolution of average reward for NChain and DoubleLoop environments, averaged over 50 runs of length 10 6 for each algorithm. For computational reasons BSS is only run for 10 4 steps. The runs are exponentially smoothened with a half-life 1000 before averaging.

Figure 5 :

 5 Figure5: Evolution of average reward for 5 × 7 and 10 × 10 LavaLake environments. The results are averaged over 20 and 30 runs respectively with a length of 10 6 for each algorithm. The runs are exponentially smoothened with a half-life 1000 before averaging.

Figure 6 :

 6 Figure 6: Evolution of average reward for the Maze environment. The results are averaged over 30 runs with a length of 10 6 for each algorithm. The runs are exponentially smoothened with a half-life 1000 before averaging.

Figure 7 :

 7 Figure7: Evolution of average reward for NChain environment with runs of length 10 6 for each algorithm. For computational reasons BSS is only run for 10 4 steps. The runs are exponentially smoothened with a half-life 1000. The mean as well as the 5th and 95th percentile performance is shown for each algorithm and the standard error is illustrated with black lines.

Figure 8 :

 8 Figure8: Evolution of average reward for DoubleLoop environment with runs of length 10 6 for each algorithm. For computational reasons BSS is only run for 10 4 steps. The runs are exponentially smoothened with a half-life 1000. The mean as well as the 5th and 95th percentile performance is shown for each algorithm and the standard error is illustrated with black lines.

Figure 11 :

 11 Figure11: Evolution of average reward for Maze environment with runs of length 10 7 for each algorithm. The runs are exponentially smoothened with a half-life 1000. The mean as well as the 5th and 95th percentile performance is shown for each algorithm and the standard error is illustrated with black lines.

Figure 12 :

 12 Figure12: Illustration of the impact of varying the horizon T in BBI. The results are averaged over 50 and 30 runs respectively with a length of 10 6 for each algorithm. The runs are exponentially smoothened with a half-life 1000 before averaging.

Table 1 :

 1 CPU-time (in seconds) used for each algorithm. Chain has 5 states and is run for 5000 steps (100 policy updates) while Maze has 264 states and is run for 20 steps (5 policy updates).

	BBI PSRL MMBI
	Chain 14	5	19
	Maze 921 6	256
	methods, alternative kernels, or a different marginalisation. We also compared it with MMBI, which
	samples of multiple MDP models from the posterior and performs approximate dynamic program-
	ming to obtain a policy, with the main difference being that MDPs are sampled only once, instead
	of at every backwards induction step.
	5 Discussion and Future Work

New Insights. We offered a new perspective on Bayesian value function estimation. The central idea is to calculate the conditional value function distribution P π β

  Yunhao Tang and Shipra Agrawal. Exploration by distributional reinforcement learning. In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 2710-2716. AAAI Press, 2018. W.R. Thompson. On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of two Samples. Biometrika, 25(3-4):285-294, 1933.

N. Vlassis, M. Ghavamzadeh, S. Mannor, and P. Poupart. Reinforcement Learning, chapter Bayesian Reinforcement Learning, pages 359-386. Springer, 2012. Tao Wang, Daniel Lizotte, Michael Bowling, and Dale Schuurmans. Bayesian sparse sampling for on-line reward optimization. In Proceedings of the 22nd international conference on Machine learning, pages 956-963, 2005.

Table 2 :

 2 Wasserstein distance to the true distribution of the value function, for Alg. 1 and the mean MDP model, for NChain. For Inferential Induction, the distances are averaged over 5 runs. The distances correspond to the plots in Figure3.

	Time steps Inf. Induction Mean MDP
	10	22.80	30.69
	100	16.41	17.90
	1000	4.18	4.27
	B.3 Further Comparisons on Discrete Environments: MMBI, BSS, BQL

In the partially-observable setting, the agent instead observes another variable dependent on the state.

For infinite horizon problems this policy is also stationary.

In the discounted setting, the value function converges to V π µ ≡ V π µ,1 as T → ∞.

This is expressible in closed form. When the MDP is discrete, a Dirichlet-product prior can be used, or when the MDP is continuous and the dynamics are assumed to be linear, a Gaussian-Wishart prior can be used[START_REF] Degroot | Optimal Statistical Decisions[END_REF]. Gaussian process inference can also be expressed in a closed-form but inference becomes approximate because the computational complexity scales quadratically with time.

Assuming otherwise results in a mean-field approximation. See Sec. 2.2. in the arXiv version of the paper.

As long as q has full support over the state space, any choice should be fine. For discrete MDPs, we use a uniform distribution q over states and sum over all of them, while we sample from q in the continuous case.

When the horizon T is small, we can set H = T -t.

More detailed figures with error bars and percentiles as well as comparisons with MMBI, BSS and BQL are available in the appendix.

For computational reasons we used a planning horizon of two with four next state samples and two reward samples in each branching step. We hope to be able to run further experiments with BSS at a later point.

Acknowledgments and Disclosure of Funding

Thank you to Nikolaos Tziortziotis for his useful discussions. This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The experiments were partly performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC).

for each algorithm. The runs are exponentially smoothened with a half-life 1000. The mean as well as the 5th and 95th percentile performance is shown for each algorithm and the standard error is illustrated with black lines.


Frédéric Dehais

Imad Rida

Raphaëlle N Roy

John Iversen
email: jiveren@ucsd.edu

Tim Mullen
email: jim.mullen@intheon.io

Daniel Callan
email: dcallan@nict.go.jp

A pBCI to Predict Attentional Error Before it Happens in Real Flight Conditions

Accident analyses have revealed that pilots can fail to process auditory stimuli such as alarms, a phenomenon known as inattentional deafness. The motivation of this research is to develop a passive brain computer interface that can predict the occurence of this critical phenomenon during real flight conditions. Ten volunteers, equipped with a dry-EEG system, had to fly a challenging flight scenario while responding to auditory alarms by button press. The behavioral results disclosed that the pilots missed 36% of the auditory alarms. ERP analyses confirm that this phenomenon affects auditory processing at an early (N100) and late (P300) stages as the consequence of a potential attentional bottleneck mechanism. Inter-subject classification was carried out over frequency features extracted three second epochs before the alarms' onset using sparse representation for classification (SRC), sparse and dense representation (SDR) and more conventional approach such as linear discriminant analysis (LDA), shrinkage LDA and nearest neighbor (1-NN). In the best case, SRC and SDR gave respectively a performance of 66.9% and 65.4% of correct mean classification rate to predict the occurrence of inattentional deafness, outperforming LDA (60.6%), sLDA (60%) and 1-NN (59.6%). These results open promising perspectives for the implementation of neuroadaptive automation with as ultimate goal to enhance alarm stimulation delivery so that it is perceived and acted upon.

I. INTRODUCTION

Operating aircraft is a demanding activity that involves the management of multiple visual (e.g. monitoring the flight parameter) and auditory tasks (e.g. radio communication) in a dynamic and uncertain environment [START_REF] Gateau | Real-time state estimation in a flight simulator using fnirs[END_REF], [START_REF] Causse | Cognitive aging and flight performances in general aviation pilots[END_REF], [START_REF] Dehais | Pilot flying and pilot monitoring's aircraft state awareness during go-around execution in aviation: A behavioral and eye tracking study[END_REF]. Distribution of attention is a key issue for piloting and relies on a tradeoff between focused attention (eg. performing a check-list) to avoid distraction, and diffused attention to detect unexpected changes (eg. failure). These top-down and bottom-up types of attention are respectively supported by the dorsal and ventral neural networks that are in close interaction [START_REF] Vossel | Dorsal and ventral attention systems: distinct neural circuits but collaborative roles[END_REF]. However, when task demand exceeds mental capacity, the homeostasis between these two neural pathways could be disrupted, leading to an impaired processing of unexpected stimuli [START_REF] Molloy | Inattentional deafness: visual load leads to time-specific suppression of auditory evoked responses[END_REF], [START_REF] Todd | Visual short-term memory load suppresses temporo-parietal junction activity and induces inattentional blindness[END_REF]. Although this shielding mechanism can prevent mental overload, missing critical information can jeopardize flight safety [START_REF] Dehais | Momentary lapse of control: A cognitive continuum approach to understanding and mitigating perseveration in human error[END_REF]. For instance, accident analyses [START_REF] Bliss | Investigation of alarm-related accidents and incidents in aviation[END_REF] and experiments conducted in flight simulators [START_REF] Dehais | Failure to detect critical auditory alerts in the cockpit: evidence for inattentional deafness[END_REF], [START_REF] Dehais | Missing critical auditory alarms in aeronautics: evidence for inattentional deafness?[END_REF], [START_REF] Dehais | The perseveration syndrome in the pilots activity: guidelines and cognitive countermeasures[END_REF] reported that auditory alarms could fail to reach awareness when engaged under demanding cognitive flying tasks. This phenomenon, known as inattentional deafness, has been shown to take place at an early stage of auditory processing [START_REF] Dehais | Inattentional deafness to auditory alarms: Inter-individual differences, electrophysiological signature and single trial classification[END_REF] through top-down inhibitory mechanisms implemented via the activation of cortical regions associated with an attentional bottleneck [START_REF] Durantin | Neural signature of inattentional deafness[END_REF]. Complementary evidence of this phenomenon was supported by an electrophysiological experiment conducted in real flight conditions which revealed that a reduction in phase resetting in alpha and theta band frequencies was a neural signature of inattentional deafness to auditory alarms [START_REF] Callan | Disruption in neural phase synchrony is related to identification of inattentional deafness in real-world setting[END_REF].

A relevant approach to improve flight safety is to implement passive brain computer interfaces (pBCI) or neuroadaptive technology [START_REF] Aricò | Passive bci in operational environments: Insights, recent advances, and future trends[END_REF], [START_REF] Zander | Neuroadaptive technology enables implicit cursor control based on medial prefrontal cortex activity[END_REF], [START_REF] Zander | Towards passive brain-computer interfaces: applying brain-computer interface technology to humanmachine systems in general[END_REF], [START_REF] Lotte | Brain-computer interface contributions to neuroergonomics[END_REF] to continuously monitor pilots' attentional state and to detect the possible occurrence of degraded states. Recently, [START_REF] Dehais | Inattentional deafness to auditory alarms: Inter-individual differences, electrophysiological signature and single trial classification[END_REF] implemented an offline pBCI to detect inattentional deafness to auditory alarms during a simulated flight. Such an approach opens promising perspectives for pilot-cockpit interaction, however a step further would be to predict rather than detect episodes of inattentional deafness [START_REF] Duprès | Pre-stimulus eeg engagement ratio predicts inattentional deafness to auditory alarms in realistic flight simulator[END_REF]. As a consequence, one could imagine the design of a smart cockpit that would implicitly adapt itself to the pilots' attentional state with a set of counter-measures.

To meet this challenging goal, supervised dictionary learning approach has been shown to be an efficient means to lead state-of-the-art results in many applications including signal classification [START_REF] Rida | Temporal signals classification[END_REF]. Indeed, in recent years there has been a growing interest in the use of techniques such as sparse representation for classification (SRC) or sparse and dense representation (SDR) in order to build discriminative representations by minimizing the intra-class homogeneity, maximizing class separability and promoting sparsity for more generalization ability [START_REF] Al Maadeed | Palmprint identification using sparse and dense hybrid representation[END_REF], [START_REF] Rida | Palmprint identification using an ensemble of sparse representations[END_REF], [START_REF] Rida | An efficient supervised dictionary learning method for audio signal recognition[END_REF]. This is done by learning a dictionary per class and making them dissimilar by boosting the pairwise orthogonality. When compared to conventional dictionary learning techniques [START_REF] Hamner | Learning dictionaries of spatial and temporal eeg primitives for brain-computer interfaces[END_REF] which they solely try to minimize the reconstruction error, supervised dictionary learning has the merit to be a very efficient way to classify EEG signals for BCI purposes [START_REF] Zhou | Discriminative dictionary learning for eeg signal classification in brain-computer interface[END_REF], [START_REF] Wen | Review of sparse representation-based classification methods on eeg signal processing for epilepsy detection, brain-computer interface and cognitive impairment[END_REF].

The objective of the present study was to develop a pBCI to predict auditory alarm misperception in the context of flying. Participants were asked to perform a demanding flying scenario along with an auditory alarm detection task in real flight conditions. In line with previous studies, we used dryelectrode EEG that has proven to be suitable for measuring This paper received the IEEE Brain Initiative award -IEEE SMC CONFERENCE -BMI 2019 Fig. 1. Left: ISAE-SUPAERO DR400 aircraft. Right: EEG dry electrode Enobio system used for the experiment brain response [START_REF] Callan | Disruption in neural phase synchrony is related to identification of inattentional deafness in real-world setting[END_REF], [START_REF] Dehais | Eegengagement index and auditory alarm misperception: an inattentional deafness study in actual flight condition[END_REF] and to perform single trial classification [START_REF] Dehais | Monitoring pilots mental workload using erps and spectral power with a six-dry-electrode eeg system in real flight conditions[END_REF], [START_REF] Gateau | Detecting pilot's engagement using fnirs connectivity features in an automated vs manual landing scenario[END_REF], [START_REF] Scholl | Classification of pilot-induced oscillations during in-flight piloting exercises using dry eeg sensor recordings[END_REF], [START_REF] Callan | Classification of singletrial auditory events using dry-wireless eeg during real and motion simulated flight[END_REF] in actual flight conditions. We proposed to use SRC as well as SDR techniques to predict inattentional deafness with frequency features computed over the EEG signal in a 3-second time-window preceding the onset of each auditory alarm. More conventional approaches including Nearest Neighbor (1-NN), Linear Discriminant Analysis (LDA) and shrinkage Linear Discriminant Analysis (sLDA) were also used as a benchmark.

II. METHODS

A. Experimental protocol

Ten pilots were recruited among the students of the ISAE-SUPAERO engineering school to participate in the study (10 males; 25-48 years old, with 40-230 flight hours experience). All had normal or corrected-to-normal vision and no history of neurological or psychiatric disorders. The study was approved by the European Aviation Safety Agency (EASA60049235) and all participants gave their informed written consent.

1) Experimental environment: DR400 aircraft:

The study was conducted using the ISAE-SUPAERO (Institut Supérieur de l'Aéronautique et de l'Espace -French Aeronautical University in Toulouse, France) experimental DR400 light aircraft (see Fig 1). The DR400 light aircraft was powered by a 180HP Lycoming engine and was equipped with classical gauges, radio navigation equipment, and actuators such as rudder, stick, thrust and switches to control the flight. A switch button was attached to the stick to collect the pilots response to auditory alarms.

2) Scenario:

The experiment lasted approximately one hour in duration and involved navigation tasks, diversions, simulated engine failure exercises, and series of touch-andgo (landings and takeoffs) involving low altitude circuits patterns. Along with the flying task, the participants had to perform an auditory oddball detection task with a total of 1000 auditory stimuli: 25% were alarm stimuli (250 normalized pure tone at 1100 Hz, 90 dB SPL) and 75% were non-alarm stimuli (750 normalized pure tone at 1000 Hz, 90 dB SPL). The inter-trial interval was set to 3000 ms with a 2000-ms jitter. The participants were instructed to respond quickly by pushing a button attached to the flight control stick after hearing the alarm stimuli. The auditory task started during takeoff and would end before or during landing.

3) Experimental protocol: The participants were first trained to perform the auditory task alone for five minutes while seated in a calm room at ISAE-SUPAERO hangar. After completing the auditory task, they were equiped with the 30 dry electrode Enobio EEG cap. The volunteers wer then left-seated on the aircraft and had to wear a Clarity Aloft headset that was used to trigger auditory stimuli from a PC via an audio cable. The participants could still communicate with the other crew members and the air traffic controllers when they received auditory alarms. The safety pilot was right seated and had the authority to stopping the experiment and taking over the control of the aircraft for any safety reason. The backseater was the experimenter: his role was to set the sensor, to trigger the experimental scenario and to supervise data collection.

4) Data recording: EEG data were recorded at 500 Hz using the 30 dry-electrode Enobio Neuroelectrics system positioned according to the 10-20 system. We used Lab Streaming Layer libraries (LSL, Swartz Center for Computational Neuroscience, UCSD, November 2018) to synchronize the oddball task in Matlab (Ver. 2017.b) and the response button with the Enobio acquisition software (NIC V2.0). The data were processed using EEGLab (Ver. 15).

5) Event Related Potential (ERP) and Frequency domain analysis with sLoreta:

Similarly to [START_REF] Callan | Disruption in neural phase synchrony is related to identification of inattentional deafness in real-world setting[END_REF], we computed ERP over Cz electrode. Our motivation was to check that we actually picked brain signals and that the N100 and P300 amplitudes were lower for misses that hits as demonstrated by [START_REF] Dehais | Inattentional deafness to auditory alarms: Inter-individual differences, electrophysiological signature and single trial classification[END_REF]. To do so, the continuous EEG data was filtered between 0.5-30 Hz (windowed-sinc FIR filter with an order of 250). Noisy portions of data (e.g., trials) were cleaned using the Artifact Subspace Reconstruction (ASR) [START_REF] Mullen | Real-time modeling and 3d visualization of source dynamics and connectivity using wearable eeg[END_REF]. Independent component analysis (ICA) was then ran over the filtered data and IClabel Label function was used to keep only brain components. The epochs for auditory misses and hits were extracted from the continuous data 0.2 s before and 1 s after stimuli onsets. The trials used for the ERP analyses were baseline normalized using data from 200 to 0 ms prior to the stimulus onset. We also ran descriptive frequency domain analyses using sLoreta over three second epochs extracted before the alarms' onset. The objective was to localize and identify potential neural mechanisms that could predict alarm misperception. To do so, we used the same pipeline described previously to the exception that the data were epoched starting 3000 ms before and ending 0 ms before the auditory misses and the auditory hits.

III. CLASSIFICATION

1) EEG processing pipeline for classification purpose:

The signal was first epoched starting 3000 ms before and ending 0 ms before the auditory alarms (see figure 2). Each epoch was then high-passed (0.5 Hz). Noisy portions of epoched data were cleaned using ASR. The ASR filter was calibrated using the first 30 s of EEG recording that were not used for the classification. Frequency based features were computed for each trial in the delta [1 4], theta [4 8] Hz, alpha [8 12] Hz, beta [12 30] Hz and gamma [30 80] Hz bands using a 250-order windowed sinc FIR-filter trial. Eventually, we computed an EEG engagement ratio β+γ α+θ adpated from [START_REF] Prinzel | A closed-loop system for examining psychophysiological measures for adaptive task allocation[END_REF].

A. Sparse representation for classification

Let us consider {(x j , y j )} n j=1 where x j ∈ R m be signal and y j ∈ {1, • • • , c} be its class label. We consider a class based dictionary {D i } c i=1 ∈ R m×ni the n

α α T = [α α α T 1 • • • α α α T i • • • α α α T c
] is given by:

min α α α 1 2 p -Dα α α 2 2 + λ α α α 1 (1)
where . 1 denotes the 1 -norm corresponding to the absolute sum of the vector α α α, λ is a parameter controlling the compromise between the reconstruction error and samplewise sparsity and α α α represents the sparse representation over the global dictionary D.

With the assumption that subspaces of distinct classes are independent to each other, the formulation (1) achieves a discriminative representation where significant nonzero coefficients are only associated to the correct subject. Thus the resulting sparse representation in (1) named in the literature Sparse Representation for Classification (SRC) [START_REF] Wright | Robust face recognition via sparse representation[END_REF] is suitable for classification. The classification is performed by computing residual reconstruction error of the test sample p using the training samples of each class i serving as a dictionary D i and their corresponding sparse coefficients α α α i as follows: The class label of the given test sample is assigned to class i that minimizes the reconstruction error using D i and α α α i .

e i = p -D i α α α i 2 2 i = 1, • • • , c (2)

B. Sparse and dense hybrid representation

Despite of the impressive results of SRC, a number of works put in doubt its effectiveness for classification [START_REF] Shi | Is face recognition really a compressive sensing problem?[END_REF] [47]. To solve this problem, it has been proposed to separate the class-specific information from others to allow the sparse representation of the query signal to coincide with the correct classification target specified by the class labels of the training data. Therefore, a query EEG signal p is decomposed into three main components as follows:

p = a + b + s ( 3
)
where a is the class-specific component, b the non-classspecific variations and s contains random sparse noise. Lets A a dictionary containing only class-specific component, the sparse representation α α α of the class-specific component of the query signal p can be computed using SRC as follows:

a = Aα α α + e a (4)
the sparse vector α α α of (4) will directly coincide with the class label of the query signal p as the both a and A only contain the class-specific information.

Unfortunately, separating the class-specific component a from a single unknown query signal p is a very challenging problem if not impossible. To address this problem, a nonclass-specific representation z is defined by the dictionary B containing non-class-specific information as follows:

b = Bz + e b ( 5
)
the component b does not contribute in the classification and hence z in ( 5) is defined as dense representation. The summation of ( 4), ( 5) and s yields the hybrid sparse-anddense representation (SDR) [START_REF] Jiang | Sparse and dense hybrid representation via dictionary decomposition for face recognition[END_REF] of the query signal as follows:

p = Aα α α + Bz + e (6)
where e = e a + e b + s is the combined representation error.

To represent a query signal p, every training sample only uses its class-specific component to compete against the others collaboratively with the non-class-specific component of all training samples. As α α α represents the class-specific information and contributes in the classification through a sparse minimization, it is chosen sparse. In the other hand, z stands for the non-class specific and does not contribute in the classification, as consequence is taken dense. The solution of the SDR, α α α, z and e is obtained by solving the following optimization problem [START_REF] Jiang | Sparse and dense hybrid representation via dictionary decomposition for face recognition[END_REF]:

       min α α α,z,e α α α 1 + β z 2 2 + γ e 1 s.t p = Aα α α + Bz + e (7)
the optimization problem ( 7) can be solved by the Augmented Lagrange Multiplier (ALM) scheme [START_REF] Bertsekas | Constrained optimization and Lagrange multiplier methods[END_REF].

The representation of a query signal p by the class-specific component of class i and the non-class-specific component of all classes collaboratively is given by:

p = AL i α α α + Bz + e i (8)
where L i ∈ R n×n is a selection operator given by

     L i (k, k) = 1 if k th training ∈ class i 0 otherwise (9)
The class-wise representation residual is defined by:

r i (p) = e -e i 2 = A(I -L i )α α α 2 ( 10
)
where I is an identity matrix. The query signal p is classified into the class that produces the minimum residual r i (p). It is very difficult if not impossible to decompose a single signal p into a class-specific component a and a non-classspecific component b. However, given a labeled training database D, it is possible to decompose it into a class-specific dictionary A, a non-class-specific dictionary B and a random sparse noise E based on machine learning from the labeled training database. The underlying principle is that the both dictionaries A and B must be low rank matrices. We can first initialize A by a very low rank matrix, for example, a matrix that contains the class means of all classes. Then we can gradually transfer more information from D to A so that after training, the class-specific dictionary A will contain much more class-specific information than just class means. This can be done by the iterative low rank matrix decomposition of A and B from D. It utilizes the low rank matrix recovery to transfer information from the supervised assigned dictionary B to A. More details can be found in [START_REF] Jiang | Sparse and dense hybrid representation via dictionary decomposition for face recognition[END_REF].

C. Classification pipeline

As explained in the section II-A.2, our participants faced 25% of auditory alarms and 75% of non-alarm stimuli. Our behavior results disclosed that our participants missed 36% of the 25% presented alarms stimuli. For classification purpose, we then used 1400 recorded EEG signals (i.e. 3s epoch before the auditory alarms): 700 miss alarms and 700 hit alarms. An equivalent number of misses and hits were selected for each participant. To avoid the dependency of the classification techniques to a specific training / testing set, we evaluate the robustness using different training and testing set, as a consequence we have divided our initial data into L = 10 different training and testing set. Supervised Dictionary techniques such as SRC as well as SDR were tested over the different extracted frequency features separately or all together. We then benchmarked their performance with more conventional approaches including Nearest Neighbor (1-NN), Linear Discriminant Analysis (LDA) and shrinkage Linear Discriminant Analysis (sLDA).

IV. RESULTS

A. Electrophysiological results

At the group level (see figure 3 left), statistical analyses disclosed lower N100 and P300 amplitude for the auditory misses than hits over Cz electrode (p = 0.01 bootstrap statistics with FDR for multiple comparisons). Descriptive analyses using sLoreta, ran over 3 second epochs preceding the alarms, disclosed that audio misses relative to hits had lower oscillatory activity in the low alpha band in the right inferior frontal gyrus (see figure 3 right) but also in the right middle frontal gyrus and the right insula.

B. Classification results

Inter-subject classification accuracy using different features as well as algorithms is depicted in Table I. It can be seen that the best results were obtained using all aggregated features.

V. DISCUSSION

The objective of this paper was to implement an EEG based pBCI to predict inattentional deafness to auditory alarms in aviation. This goal was challenging as the EEG data were collected in a real flight conditions. The behavioral results revealed that our participants missed 36% of auditory alarms on average, confirming that such a phenomenon could take place in the cockpit [START_REF] Dehais | Missing critical auditory alarms in aeronautics: evidence for inattentional deafness?[END_REF], [START_REF] Dehais | Failure to detect critical auditory alerts in the cockpit: evidence for inattentional deafness[END_REF], [START_REF] Callan | Disruption in neural phase synchrony is related to identification of inattentional deafness in real-world setting[END_REF], [START_REF] Durantin | Neural signature of inattentional deafness[END_REF]. ERPs analyses over Cz electrode disclosed that this phenomenon takes place at a perceptual (N100) and attentional (P300) level as previously demonstrated in flight simulator [START_REF] Dehais | Inattentional deafness to auditory alarms: Inter-individual differences, electrophysiological signature and single trial classification[END_REF]. Descriptive analyses using sLoreta indicated lower alpha oscillatory activity in brain regions generally involved in attentional bottleneck processing including the inferior frontal gyrus, the insula and the superior medial frontal cortex [START_REF] Tombu | A unified attentional bottleneck in the human brain[END_REF]. As higher alpha oscillation are thought to reflect inhibition mechanisms, lower alpha oscillation preceding misses may suggest, on the contrary, greater activation of the attentional bottleneck to prevent the processing of incoming alarms [START_REF] Callan | Disruption in neural phase synchrony is related to identification of inattentional deafness in real-world setting[END_REF]. Our classification results showed that the mean prediction accuracy rate reached almost 67% in the best case with SRC. This performance is close to the 70 % accuracy defined as a sufficient accuracy for BCI [START_REF] Kubler | Bci meeting 2005-workshop on clinical issues and applications[END_REF]. The use of supervised dictionary learning approach seems to outperform more conventional techniques such as 1-NN and LDA/sLDA (60%). These latter results are consisent of our previous study with shrinkage LDA [START_REF] Duprès | Pre-stimulus eeg engagement ratio predicts inattentional deafness to auditory alarms in realistic flight simulator[END_REF] which had 58% mean accuracy to predict the onset of inattentional deafness in simulated conditions. Moreover, this prior study considered intra-subject but not inter-subject classification as was used in the present study. Therefore, these findings show that our approach is a well-suited method for processing EEG signal as previously demonstrated by [START_REF] Zhou | Discriminative dictionary learning for eeg signal classification in brain-computer interface[END_REF], [START_REF] Wen | Review of sparse representation-based classification methods on eeg signal processing for epilepsy detection, brain-computer interface and cognitive impairment[END_REF]. Moreover, the performance obtained here with inter-subject classification and a dry-EEG system is altogether quite high in a general manner, but also very promising in its robustness to intersubject variability and to ecological settings. Taken together, these findings open good prospects for the implementation of a smart cockpit that would adapt to the user's state. For instance, one could imagine that the modality (e.g. tactile, visual) of the alarm could be optimized to increase the likelihood of reaching pilot's awareness. Another possibility would be to consider adaptive autonomy and to lower the pilots' engagement by giving more authority to the autoflight system to reduce their workload. Nonetheless, there is a need to improve the accuracy of this pBCI before it could be implemented into real cockpits. Future research should investigate other metrics such as connectivity features that have been shown to efficiently predict long term attentional performance [START_REF] Senoussi | Pre-stimulus antero-posterior eeg connectivity predicts performance in a uav monitoring task[END_REF]. Another possible approach would be to apply adaptive mixture independent component analysis to identify different brain network associated with auditory misses and hits, as it has previously been shown to identify brain dynamics underlying attention fluctuation in driving [START_REF] Hsu | Modeling brain dynamic state changes with adaptive mixture independent component analysis[END_REF].

  i training samples associated to class i. The global dictionary D = [D 1 • • • D c ] ∈ R m×n represents the concatenation of the class based dictionaries {D i } c i=1 . The sparse representation of a test sample p over the global dictionary D noted α

Fig. 2 .

 2 Fig. 2. Classification pipeline

Fig. 3 .

 3 Fig. 3. Left: Grand averaged waveforms of the ERPs for Fz with standard error. The black lines on the x axis specify the time range when the auditory misses and hits related ERPs amplitudes were significantly different (p = 0.01, corrected). Right: Frequency domain EEG source analysis using sLoreta. Lower alpha oscillatory activity in the right inferior frontal gyrus during misses than auditory hits.

  ± 3.29 57.29 ± 2.85 57.38 ± 4.06 58.21 ± 2.15 59.50 ± 3.01 58.04 ± 1.76 59.60 ± 2.70 LDA 60.20 ± 4.15 59.60 ± 2.79 58.71 ± 2.25 58.67 ± 3.06 58.50 ± 3.60 62.20 ± 2.50 60.60 ± 4.00 sLDA 60.75 ± 3.64 54.38 ± 3.45 53.38 ± 3.13 53.96 ± 3.20 56.25 ± 2.56 59.25 ± 3.33 60.00 ± 3.07 SDR 61.50 ± 3.50 62.60 ± 2.80 60.50 ± 1.80 60.40 ± 1.80 58.90 ± 1.60 62.50 ± 3.07 65.40 ± 2.80 SRC 65.60 ± 4.02 64.58 ± 2.25 63.83 ± 3.37 63.96 ± 3.42 64.08 ± 3.78 63.58 ± 2.94 66.90 ± 3.10

					TABLE I
				INTER-SUBJECT CLASSIFICATION ACCURACY.
	Methods Features
		Delta	Theta	Alpha	Beta	Gamma	Engagement	Fusion
	1-NN	59.08

ACKNOWLEDGMENT

This study was supported by the Agence Innovation Defense (AID -DGA), the AXA research fund and the Artificial and Natural Intelligence Toulouse Institute (ANITI). The authors wish to express their gratitude to Fabrice Bazelot, Stephane Juaneda and all the pilots who took part in the experiments.