Abdol Aziz

Ould Ismail

Drew Parker

Moises Hernandez-Fernandez

Steven Brem

Simon Alexander

Ofer Pasternak

Emmanuel Caruyer

Ragini Verma

Characterizing Peritumoral Tissue Using Free Water Elimination in Clinical DTI

Keywords: Diffusion tensor imaging (DTI), Fractional anisotropy (FA), Mean diffusivity

à 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

The characterization of tissue microstructure in the peritumoral region is crucial for accurate surgical planning (e.g., using tractography), to identify infiltration for targeted treatment, and to provide new radiomic features for distinguishing tumor subtypes. Accurate white matter (WM) microstructural modeling using diffusion magnetic resonance imaging (dMRI) has the potential to provide a characterization of the peritumoral tissue. Such microstructure modeling is expected to have a substantial impact on future clinical studies in pathologies.

A key aspect for exploring the microstructural properties is to estimate the free water component within a tissue of interest. Not accounting for free water can cause the fitted diffusion tensors to have an erroneously low fractional anisotropy (FA), causing some tracking algorithms to stop prematurely in the contaminated regions. Hence, a correct estimation of free water is crucial.

In one of the earliest attempts for estimating free water, a bi-tensor model was proposed [1], separating the derived diffusion signal into isotropic (free water) and tissue-based components. This approach acknowledged the fact that the bi-tensor model fitting has infinitely many solutions for single-shell data [2], and posited that the illposed nature of the problem could be addressed by appropriate initialization of the model parameters and by spatial regularization that stabilizes the fit. The regularized approach significantly contributed to the field of free water elimination (FWE), allowing a better reconstruction of healthy fornix tracts, as well as improving tractography in the peritumoral region. However, the initialization that was originally proposed is acquisition dependent [3]. In addition, there are various pathological conditions in which the initialization in [1] may lead to an inaccurate estimation of the free water compartment and diffusion indices (FA, MD) [START_REF] Ismail | Free watER iNvariant Estimation of Tensor (FERNET): Addressing the Issue of Edema in Clinically Feasible Acquisitions[END_REF].

In this work, we propose a novel initialization approach based on the estimation of free water maps from simulated ground truth data. This initialization relies on prior and heuristic constraints that were validated across a range of simulated free water fractions, anisotropy levels, and underlying diffusivities. Such an approach allows for optimizing the FWE in healthy tissue as well as regions affected by pathologies such as edema and infiltration. Therefore, the initialization method that we propose is applicable to standard clinical acquisitions, enabling retrospective investigations and providing novel insight into peritumoral regions. As the method is applicable on any clinically acquired data, the free water map is expected to become a crucial feature in radiomics.

METHODS

Free Water Elimination

The proposed FWE initialization is based on a bi-compartment model: a tensor for modeling the underlying tissue, and an isotropic free water compartment. Specifically, we fit [1]

𝐴 𝑖 = 𝑓 𝑒 -𝑏𝑞 𝑖 𝑇 𝐷𝑞 𝑖 + (1 -𝑓) 𝑒 -𝑏𝑑 [1]
where the first term models the tissue and the second term represents the free water compartment. f is the tissue volume fraction, Ai is the signal attenuation of the diffusion weighted image acquired along the i th gradient direction, b is the amount of diffusion weighting, 𝑞 𝑖 is i th gradient direction, D represents the diffusion tensor used for modeling the tissue compartment, and d is the diffusivity in the isotropic compartment, which is fixed at 3.0 x 10 -3 mm 2 /s. Fitting this model using a single-shell dMRI acquisition is a problem with infinitely many solutions. Previously, the fitting was stabilized by using a regularizer post-initialization [1]. There, Pasternak et al. proposed an initial estimate of the tissue compartment f based on scaling the unweighted images (𝑆 0 ) with respect to the mean signal of representative WM voxels (𝑆 𝑡 ) and CSF voxels (𝑆 𝑤 ) using:

𝑓 𝑡=0 = 1 - log ( 𝑆 0 𝑆 𝑡 ⁄ ) log ( 𝑆 𝑤 𝑆 𝑡 ⁄ ) [2]
In [1], f was additionally constrained within bounds that were enforced during the nonlinear fit stage via

𝐴 𝑚𝑖𝑛 -𝑒 -𝑏𝑑 𝑒 -𝑏𝜆 𝑚𝑎𝑥 -𝑒 -𝑏𝑑 < 𝑓 𝜆 𝑚𝑖𝑛,𝑚𝑎𝑥 < 𝐴 𝑚𝑎𝑥 -𝑒 -𝑏𝑑 𝑒 -𝑏𝜆 𝑚𝑖𝑛 -𝑒 -𝑏𝑑 [3]
where 𝜆 𝑚𝑎𝑥 and 𝜆 𝑚𝑖𝑛 are the expected maximal and minimal diffusivity of the tensor, and 𝐴 𝑚𝑎𝑥 and 𝐴 𝑚𝑖𝑛 are the maximum and minimum measured signal attenuations within a region of interest. This initialization approach suppressed the free water in healthy tissue, yet in the final estimation in this approach as well as other FWE approaches it is clear that free water content in healthy WM tissue is low but not zero [1,[START_REF] Hoy | Optimization of a free water elimination two-compartment model for diffusion tensor imaging[END_REF]. Moreover, the initialization proposed in eq. ( 2) is solely based on T2 contrast observed in b0 images. While this T2 contrast gives us information about the amount of free water, it may miss additional information that we can only infer using the diffusion-weighted attenuation.

An alternative approach is to constrain the corrected tissue signal attenuation by assuming that the mean diffusivity (MD) in the healthy tissue is around a fixed diffusivity, e.g., 𝑀𝐷 𝑡𝑖𝑠𝑠𝑢𝑒 = 0.60x10 -3 mm 2 /s (as suggested by FWE tensor maps using multishell acquisitions [START_REF] Hoy | Optimization of a free water elimination two-compartment model for diffusion tensor imaging[END_REF][START_REF] Kodiweera | Data of NODDI diffusion metrics in the brain and computer simulation of hybrid diffusion imaging (HYDI) acquisition scheme[END_REF]), and thus the free water compartment explains the increase in MD in voxels affected by partial voluming. As such:

𝑓 𝑀𝐷 = 𝑒 -𝑏𝑀𝐷 -𝑒 -𝑏𝑑 𝑒 -𝑏𝑀𝐷 𝑡𝑖𝑠𝑠𝑢𝑒 -𝑒 -𝑏𝑑 [4]
where MD is the mean diffusivity from the standard tensor fit in a voxel of interest. Although such an approach improves the estimation in the healthy tissue, it struggles to characterize regions that are highly contaminated by free water or have tumor cells that restrict the diffusivity further in the peritumoral area, i.e., where 𝑀𝐷 𝑡𝑖𝑠𝑠𝑢𝑒 is not a good initial value. The underlying problem in all these approaches, is that the initialization is designed to either address the healthy tissue or tissue with pathology, but not both.

To alleviate these issues, we propose an initialization to the FWE problem that is a hybrid between 𝑓 𝜆 𝑚𝑖𝑛,𝑚𝑎𝑥 and 𝑓 𝑀𝐷 . Aiming to obtain a better estimation of the free water compartment in both healthy tissue as well as voxels contaminated by edema or partial voluming, we introduce an initialization that is a log-linear interpolation between 𝑓 𝜆 𝑚𝑖𝑛,𝑚𝑎𝑥 and 𝑓 𝑀𝐷 . We propose:

𝑓 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 = 𝑓 𝑀𝐷 𝛼 × 𝑓 𝜆 𝑚𝑖𝑛,𝑚𝑎𝑥 1-𝛼 [5]
where 𝛼 is set to 𝛼 = 𝑓 𝑡=0 (see eq. ( 2)). That is, for tissue with healthy appearing T2, we constrain the free water using 𝑓 𝑀𝐷 (eq. ( 4)), and in regions that appear like CSF, or edematous in the T2 contrast, we estimate the free water compartment based on eq. (3).

Evaluation

We will evaluate the different FWE initialization approaches on simulated data (where the ground truth is known) as well as on brain tumor patients. On the patients, we will evaluate the fit in the healthy tissue, defined as white matter outside the peritumoral region, as well as in peritumoral region, and its effect on tracking for the purposes of surgical planning.

Simulated Data

To have a sample with ground truth, simulated data was generated with varying ground truth mean diffusivities, anisotropy, and free water volume fractions (Table 1). These simulated datasets follow a bi-tensor model where one of the compartments represents tissue and the second is isotropic with a fixed diffusivity (3.0 x 10 -3 mm 2 /s) [START_REF] Garyfallidis | Dipy, a library for the analysis of diffusion MRI data[END_REF]. For the unweighted images, we calculate the transverse magnetization as a linear sum of the contribution of each compartment:

𝑀 𝑇 = 𝑉𝐹 𝑤𝑚 𝜌 𝑊𝑀 (1 -𝑒 -𝑇𝑅 𝑇1 𝑊𝑀 ⁄ ) 𝑒 -𝑇𝐸 𝑇2 𝑊𝑀 ⁄ + 𝑉𝐹 𝐶𝑆𝐹 𝜌 𝐶𝑆𝐹 (1 -𝑒 -𝑇𝑅 𝑇1 𝐶𝑆𝐹 ⁄
) 𝑒

-𝑇𝐸 𝑇2 𝐶𝑆𝐹 ⁄ [START_REF] Kodiweera | Data of NODDI diffusion metrics in the brain and computer simulation of hybrid diffusion imaging (HYDI) acquisition scheme[END_REF] where VF is the volume fraction, 𝜌 is the proton density, and TR/TE are repetition and echo times. For every experiment, we perform 10 noise realizations by adding Rician noise equivalent to SNR=40.

Brain Tumor Datasets

138 brain tumor patients were included (88 glioblastoma / 50 metastasis). The data was acquired on Siemens 3T Verio scanner with TR/TE=5000/86ms, b=1000s/mm 2 , 3 b0, and 30 gradient directions. Automatic tumor and peritumoral region segmentation was performed on coregistered T1, T2, T2-FLAIR and T1-contrast enhanced images using GLISTR [START_REF] Gooya | GLISTR: glioma image segmentation and registration[END_REF].

Tractography

Tractography was performed, using Diffusion Toolkit [START_REF] Wang | Diffusion toolkit: a software package for diffusion imaging data processing and tractography[END_REF], on the FWE tensor field as well as the standard tensor field, in voxels with FA values exceeding 0.2, using the 2nd order Runge-Kutta algorithm with an angle threshold of 45 ο and a step size of 1 mm. Five bundles of interest (corticospinal tract, inferior frontal, inferior frontooccipital, uncinate and arcuate fasciculi) were extracted from each tractogram in each hemisphere, using the RecoBundles algorithm [START_REF] Garyfallidis | Recognition of white matter bundles using local and global streamline-based registration and clustering[END_REF], with a pruning parameter of 7 mm. Finally, a "coverage measure" was defined and calculated for each patient and each of the two tractograms. "Coverage Measure" is the percentage of voxels in the peritumoral edema region containing a streamline from any of the ten bundles of interest.

RESULTS

Evaluation on simulated data

We compare the proposed method to previous approaches [1,4] that fit a bicompartment model. Our results show that the FW compartment is estimated to be nearly 0 in healthy WM, in the previous approaches (Fig1. C.). In simulated edema scenarios (volume fraction (VF)=0.6), The proposed initialization of free water and FA is closer to the ground truth in simulated edematous/infiltrated regions (Fig. 1A-F), as compared to the other initialization approaches.

In-vivo Data: Application on Brain Tumor datasets

Fig. 2 shows the FW and FA maps, pre-and post-FWE, using the proposed approach, 𝑓 𝑀𝐷 initialization, as well as α =0 and α =1. Several things are observed: (i) FA maps from the corrected tensors show that various approaches were successful in increasing FA in the edema region (Fig. 2.A); (ii) FW is estimated to be nearly zero in many regions of healthy WM by previous approaches, as compared to our proposed method; (iii) a large number of physiologically implausible voxels (defined as MD<0.40x10 -3 mm 2 /s [START_REF] Hoy | Optimization of a free water elimination two-compartment model for diffusion tensor imaging[END_REF][START_REF] Kodiweera | Data of NODDI diffusion metrics in the brain and computer simulation of hybrid diffusion imaging (HYDI) acquisition scheme[END_REF]) were observed when using the method in [1] (Fig. 2

.C).

Although this was alleviated by use of α =0, this approach underestimated the free water [START_REF] Hoy | Optimization of a free water elimination two-compartment model for diffusion tensor imaging[END_REF][START_REF] Kodiweera | Data of NODDI diffusion metrics in the brain and computer simulation of hybrid diffusion imaging (HYDI) acquisition scheme[END_REF] in almost the entire healthy WM (Fig. 2.C). The proposed method, as well as the one that used α =1, were able to eliminate these implausible voxels, without underestimating the FW in the healthy tissue (Fig 2 .C).

Tractography

Fig. 3.A shows the comparative results of tracking in the peritumoral regions, pre-and post-FWE (using the proposed initialization) showing that the streamlines were able to travel through the peritumoral region more after FWE. Additionally, in Fig. 3.B., a quantitative evaluation using our coverage measure shows that the tracking coverage in the peritumoral region is significantly improved with FWE.

Fig. 3. A)

Comparison of tracking in the peri-tumoral region (gray blob), pre-and post FWE in the left arcuate and CST reconstructed in a tumor patient using the proposed initialization. In the overlap of the approaches, we see the significant increase of streamlines coverage in the edema. B) Tracking coverage in the edema was computed before and after FWE. A paired t-test analysis of tractography showed a significant increase in peritumoral region coverage using FWE tractography (t=6.91, p<1×10 -9 ).

DISCUSSION

In this manuscript, we proposed a novel initialization method for estimating the free water compartment from single-shell acquisitions. Our proposed method addresses initialization challenges of the two-compartment model, improves the fitting process and reduces implausible voxels in the FWE output.

Although there is a growing trend of dMRI Multi-shell protocols in the clinical environments, the single-shell acquisitions is still the most frequent sequence in clinical studies. Here we have shown that traditional approaches for characterizing the tissue microstructure with free-water elimination in single shell data, struggle in many scenarios due to the lack of a proper initialization that accounts properly for both healthy tissue and tissue contaminated by edema or cancer cells. The findings of our novel approach demonstrate that interpolating diffusivity and T2 information, in the initialization, significantly improved the estimation in the healthy and contaminated tissues when compared to previous initialization methods [1,[START_REF] Ismail | Free watER iNvariant Estimation of Tensor (FERNET): Addressing the Issue of Edema in Clinically Feasible Acquisitions[END_REF].

The fact that previous FWE initialization approaches [1,[START_REF] Ismail | Free watER iNvariant Estimation of Tensor (FERNET): Addressing the Issue of Edema in Clinically Feasible Acquisitions[END_REF] estimated free water to be zero in healthy white matter, is evident both in the simulated data (Fig. 1C), as well as the healthy tissue regions in Fig 2 .C. Initializations with 𝑓 𝑀𝐷 and our proposed method, were able to produce values that are closer to the ground truth in healthy tissue (Fig 2). However, simulations show that 𝑓 𝑀𝐷 approach is misleading when the underlying diffusivity is different than that of the healthy tissue (Fig. 1.A-B). This presents a clinical challenge in the case of the cytotoxic edema, where the underlying apparent diffusivity reduces significantly compared to that of the healthy tissue. On the other hand, the proposed approach yields a robust solution across different diffusivities in the simulated peritumoral region.

Furthermore, the proposed strategy for initializing the parameters of a bicompartment model leads to a better characterization of the underlying white matter in areas affected by edema. This affects the results of tractography, allowing the reconstruction of tracts impacted by edema (Fig. 3). These findings may play an important role in the pre-surgical planning, where the non-invasive reconstruction of white matter tracts is critical, especially in the regions with edema, so as to not miss eloquent tracts. This provides a safer resection margin for the surgeon.

An interesting question to pursue in future work is how different optimization routines affect the accuracy of free water estimation. In this study we have used the Gradient Descent optimization algorithm. In the future, we will compare this algorithm with Levenberg-Marquardt algorithm, and moreover, we will apply a Bayesian approach for fitting the model, Markov Chain Monte Carlo (MCMC). Additionally, the differences that are evident in the peritumoral tissue in the metastatic tumors and the GBMs, suggest the potential of free water measures in differentiating tumor subtypes. This will be undertaken as future work.

Fig. 1 .

 1 Fig.1. Three different ground truth (GT) MD values (GT-MD) (1 st , 2 nd , and 3 rd columns) were investigated using different GT FW and GT FA values. For details on the simulation parameters, see Table1. Four FWE initializations were performed (Pasternak et al. 2009, 𝑓 𝜆 𝑚𝑖𝑛,𝑚𝑎𝑥 (i.e., α =0), 𝑓 𝑀𝐷 , (i.e., α=1) and the proposed method with varying alpha). The free water (FW VF) and FA (FW FA) outputs of each method are compared to the corresponding ground truth values. The 1 st row represents the investigation of the FW maps in a simulated WM across different FW fractions. The 2 nd row is for FW FA maps for a simulated edema scenario (VF=0.60).

Fig. 2 .

 2 Fig. 2. A) FA maps from the various FWE initializations. B) Free water maps using different initialization strategies. C) Location of implausible and voxels with underestimated FW (FW volume fraction <0.02) in the different investigated methods.

  O, Sochen N, Gur Y, Intrator N, Assaf Y: Free water elimination and mapping from diffusion MRI. Magnetic resonance in medicine 2009, 62:717-730. 2. Scherrer B, Warfield SK. Why multiple b-values are required for multi-tensor models: evaluation with a constrained log-Euclidean model. InISBI 2010 (pp. 1389-1392). 3. Parker D, Ould Ismail A, Alexander S, Verma R. The Role of Bias Field Correction in the Free Water Elimination Problem. In ISMRM June, 2018

Table 1 .

 1 Parameters selected for simulating data to validate the proposed method

		Values	Number of experiments
	MD Values	0.4-1.0 x 10 -3 mm 2 /s 4
	FA Values	0-1	11
	Free Water Volume Fractions	0-1	11
	Rotation of the Tensor Directions -	10
	Noise Realizations	-	10
	Total		48,400


Nicolas Carrara
email: nicolas.carrara@orange.com

Romain Laroche
email: romain.laroche@microsoft.com

Jean-Léon Bouraoui
email: jeanleon.bouraoui@orange.com

Tanguy Urvoy
email: tanguy.urvoy@orange.com

Olivier Pietquin
email: olivier.pietquin@univ-lille1.fr

A Fitted-Q Algorithm for Budgeted MDPs

We address the problem of budgeted/constrained reinforcement learning in continuous state-space using a batch of transitions. For this purpose, we introduce a novel algorithm called Budgeted Fitted-Q (BFTQ). We carry out some preliminary benchmarks on a continuous 2-D world. They show that BFTQ performs as well as a penalized Fitted-Q algorithm while also allowing ones to adapt the trained policy on-the-fly for a given amount of budget and without the need of engineering the reward penalties. We believe that the general principles used to design BFTQ could be used to extend others classical reinforcement learning algorithms to budget-oriented applications.

Introduction

Classic reinforcement learning (RL) algorithms focus on the maximization of a unique expected reward signal. But many applications have multiple, possibly conflicting, objectives. For example, an autonomous driving car must optimize its travel time, its fuel consumption, as well as the safety of people. In this example, the first two objectives are of the same nature. Although conflicting on short term, they become coherent on the long term. If the fuel tank is depleted, then the travel time will be potentially infinite. In practice, by replacing this infinite value by the time needed to walk to the next fuel station, one gets a reasonable way to combine these two objectives into a single one. The safety objective is of different nature: it is quite difficult to assess the cost of a dramatic car accident and it's often more intuitive to think in term of probability for a given run to fail the constraint. *

*now with Google Brain, Paris

This problem is better formalized as a worst-case constrained or risk-sensitive RL problem: the car should optimize the expected travel time while guaranteeing a crash probability strictly lower than the one of a human driver. The natural formalism for constrained RL is given by Constrained Markov Decision Processes (CMDP [START_REF] Beutler | Optimal policies for controlled markov chains with a constraint[END_REF][START_REF] Altman | Constrained Markov Decision Processes[END_REF].

The most direct approach to these kind of problems is to relax the constraint by adding a reward penalty. That is to say, if the cost is C we replace the reward R by R -λC where λ > 0 is the strength of the penalty. But the calibration of this penalty requires the ability to simulate the policies.

Most of the existing solutions focus on the modelling of one policy for each budget. But in some practical situations, the budget can decrease with time. Let us consider the example of a military swarm of drones in a reconnaissance mission. In that case, the individual drones policies should become more risk-averse as the swarm population decreases. For these kind of problems, the budget is a parameter of the policy; they are formalized as Budgeted Markov Decision Process (BMDP [START_REF] Boutilier | Budget allocation using weakly coupled, constrained markov decision processes[END_REF]. Solving a BMDP is a harder problem than solving a CMDP.

The aim of our work is to handle budgeted RL problems:

• Directly from batch of transitions (i.e. without simulation or interaction with the environment during training);

• In a continuous-state space;

• With the budget as a parameter of the policy.

For this purpose we introduce Budgeted Fitted-Q (BFTQ), an extension of the Fitted-Q (FTQ) algorithm. It takes as input a batch of RL transitions and estimates the optimal parametric strategies under budget. We experiment BFTQ on 2-D worlds with forbidden areas. We show that BFTQ performs as well as a FTQ with reward penalties.

2 Background and Related Works

Markov Decision Process

A Markov Decision Process (MDP) is formalized as a tuple S, A, R, P, γ ; where S is the state set, A is the action set, R : S ×A×S → R is the (potentially stochastic) reward function, P : S × A × S → [0, 1] is the transition function, and γ is the reward discount factor. The behaviour of the agent is defined as a policy π : S → A, which can either be deterministic or stochastic. Solving an MDP consists in finding a policy π * that maximises the γ-discounted expected return

E π * [ t γ t R(s t , a t , s t+1 )].
The optimal policy π * satisfies Bellman's optimality equation [START_REF] Bellman | Dynamic programming and lagrange multipliers[END_REF]):

π * (s) = argmax a∈A Q * (s, a), Q * (s, a) =E s |s,a [R (s, a, s ) + γQ * (s , π * (s ))] If γ < 1,
{(s i , a i ), r i + γ * max a Q(s i , a )} i∈[0,N ]

Budgeted Markov Decision Processes

A Constrained Markov Decision Process (CMDP, [START_REF] Beutler | Optimal policies for controlled markov chains with a constraint[END_REF]; [START_REF] Altman | Constrained Markov Decision Processes[END_REF]) is an MDP augmented with a cost function C : S × A × S → R, a cost discount γ c , and a budget β. The cost function and its budget represent the "hard" constraint of the problem, while the reward function represents the task to complete. The problem of constrained RL consists in finding the policy which maximizes the expected return while keeping the discounted sum of the cost under the given budget β.

Formally:

π * = argmax π E π t γ t R t , (1)
subject to

E π t γ t c C t ≤ β (2)
where R t = R(s t , a t , s t+1 ) and C t = C(s t , a t , s t+1 ).

For the sake of simplicity we consider a weak form of risk in Eq. ( 2): the budget is only constrained on average. For sensitive applications the failure risk should be modelled in term of worst-deviation from the expected cumulative cost or in term of worst-percentile [START_REF] Chow | Risk-constrained reinforcement learning with percentile risk criteria[END_REF].

The Budgeted Markov Decision Process (BMDP, [START_REF] Boutilier | Budget allocation using weakly coupled, constrained markov decision processes[END_REF]) problem is a generalisation of the CMDP problem where the objective is to find a generic policy π * β which works for any CMDP of budget β > 0.

A remark on deterministic policies

A convenient property of MDPs, is that the optimal policy is unique, deterministic and greedy: π * (s) = argmax a Q(s, a). In a CMDP, and a fortiori in a BMDP, this is in general not the case. It has been shown indeed that the optimal policy under constraint is a random mixture of two deterministic policies Beutler and Ross (1985, Theorem 4.4).

To illustrate this fact, let us consider the trivial BMDP on the left of Fig. 1. On this example we have E π R = 10π 1 and E π C = π 1 . The deterministic policy consisting in always picking the safe action is feasible for any β ≥ 0.

But if β = 1/2, the most rewarding feasible policy is to randomly pick the safe and risky actions with equal probabilities. If we attempt to cast this BMDP into an MDP by replacing the costs by negative rewards, the policy we will obtain will be deterministic, hence suboptimal.

Related work

Optimization under constraint algorithms have proven their worth in many real life applications like operational efficiency for businesses or network design problems during the 90s [START_REF] Bertsekas | Constrained Optimization and Lagrange Multiplier Methods (Optimization and Neural Computation Series[END_REF]). But, these algorithms do not capture the temporal nature of a wide class of RL problems, such as in robotics, dialogue systems, control systems, or trading strategies.

Figure 1: On the left hand side, a simple risky-vs-safe BMDP. The probability of picking the risky action is π 1 . On the right hand side an attempt to relax the problem with negative rewards.

As a consequence, researchers got interested in resolving constrained sequential problems [START_REF] García | A Comprehensive Survey on Safe Reinforcement Learning[END_REF]), such as MDP under constraints [START_REF] Beutler | Optimal policies for controlled markov chains with a constraint[END_REF]; [START_REF] Altman | Constrained Markov Decision Processes[END_REF]). [START_REF] Undurti | Function Approximation for Continuous Constrained MDPs[END_REF] proposes an algorithm to solve CMDP with a continuous state-space using the same principles as FVI; but they assume that the environment (R, C and P ) is known.

Geibel and Wysotzki ( 2005 It is important to note, that these algorithms are designed to work for CMDP and require some adaptation to fit the general BMDP problem. On the contrary, Boutilier and Lu (2016) resolve a BMDP using a VI-like algorithm; however it only applies to finite and known environments.

Penalized Fitted-Q Iteration

As mentioned before, a natural way to relax the constraint of Eq. ( 2) is to penalize states with high cost. Assuming γ = γ = γ c , the MDP problem to solve becomes:

π * λ = argmax π E π t γ t (R t -λC t ) , (3)
As explained on Fig. 2, the optimal deterministic policy can be obtained by a line-search on the Lagrange multiplier values λ. Then, according to Beutler and Ross (1985, The optimal multiplier λ * avg is the smallest one to satisfy the budget constraint on average. Safer policies can also be selected according to the largest deviation form this mean cost. Theorem 4.4), the optimal policy is a randomized mixture of two deterministic policies: the safest deterministic policy that violates the constraint π λ-and the riskier of the feasible ones π λ+ .

Fitted-Q can be easily adapted continuous states CMDP and BMDP through this methodology, but, given the high variance, it requires a lot of simulations to get a proper estimate of the calibration curve. Our purpose is to avoid this calibration phase.

Budgeted Fitted-Q Iteration

In this section, we introduce Budgeted-Fitted-Q (BFTQ), a batch RL under constraints algorithm for BMDP, as an extension of the regular FTQ algorithm.

Intuition

We illustrate the principle of BFTQ on a simple BMDP described on Fig. 3. This example counts 7 states and 2 actions, its transition, reward, and constraint functions are fully deterministic. In order to maximise the cumulative reward while keeping the sum of the cost under β, there are only 3 parameters to optimise: π 0 = P(a 0 | start), π 1 = P(a 0 | s 2 ), and π 2 = P(a 0 | s 2 ).

First note that on the left hand side of the figure, on state s 1 , the action a 0 dominates a 1 for any budget: we can safely set π 1 = 1 and replace state s 1 by a terminal state of expected reward 10 and cost 5. We write: Q r (s 0 , a 0 ) = 10 and Q c (s 0 , a 0 ) = 2.

Assuming γ = γ c = 1, the optimization problem of Eq. 1 simplifies into:

max π0,π2 10π 0 + (1 -π 0 ) • (2 + 9π 2 ) s.t. 2π 0 + (1 -π 0 ) • (1 + 2π 2 ) ≤ β
The best policy depends on the budget: If β < 1 the problem is infeasible, if β ≥ 3 we maximize reward by playing a 1 , then a 0 . If β = 3/2, the best policy is to take π 0 = 1/2 and π 2 = 0, if β = 5/2, we rather take π 0 = 1/2 and π 2 = 1.

It is also possible to solve the problem locally on s 2 as long as the budget β ≥ 1 is left as a free parameter:

π * 2 (β) = argmax π2 9π 2 s.t. 2π 2 ≤ β -1 = min {1, (β -1)/2}
We then write:

Q r (s 0 , a 1 , β) = 2 + 9π * 2 (β) Q c (s 0 , a 1 , β) = 1 + 2π * 2 (β)
The final optimization problem becomes:

max π0,β1 π 0 Q π r (s 0 , a 0 ) + (1 -π 0 )Q π r (s 0 , a 1 , β 1 ) (4) s.t.: π 0 Q π c (s 0 , a 0 ) + (1 -π 0 )Q π c (s 0 , a 1 , β 1 ) ≤ β
If the values Q r and Q c are known for each state/action/budget triplet we can drive the agent locally according to this equation.

For tree-shaped BMDP like this one, the global optimization problem can always be separated in independent sub-problems (one for each sub-tree)1 . Intuitively speaking, the agent optimises its budget distribution over all possible actions.

Training Algorithm

The preceding example can be extended to continuousstates environments by the use of function approximation.

In BFTQ, two functions have to be approximated. The regular Q-function estimating the expected discounted sum of rewards, denoted here by Q r , is defined on a stateaction-budget triplet. The constraint Q-function for the sum of constraints, denoted here by Q c , is defined on a state-action-budget triplet as well. For both Q-functions, the budget variable contextualises with the current budget remaining when evaluating the function.

Notation Description n

Current iteration of BFTQ. i Index of the i th transition.

f (x i ) reg ← --y i Regression of f using database {(x i , y i )} i∈[1,N ] . ∆ A
Probability simplex over A.

Table 1: Notations used in BFTQ In Fig. 3, the environment transitions, rewards and constraints were assumed to be known. In most situations, the model of the environment is unknown. FTQ, like many others RL algorithms, overcomes this problem by using the agent interactions with the environment through a trials and errors process. Those interactions are collected as a base of transitions. To take the budget into account while constructing the optimal policy in BFTQ, this training base is enriched as a list: (s i , a i , r i , s i , c , β i ) i∈[0,N ] where r i and c i denote respectively the effective rewards and the costs of the transition and the β i are taken from a discretization of the admissible budgets2 .

The training phase of BFTQ relies on two policies: with the notations summarized on Table 1, the behaviouralpolicy π(s, a, β) : S ×R → ∆ A is the classical stochastic policy mapping states to distributions over actions, augmented with the budget dependency. The budget-policy b(s, a, β) : S × A × R → R allocates constraint budget to each triplet (action, state, budget).

The BFTQ algorithm consists in iterating over its two Q-functions: the reward action-value Q-function Q n r is evaluated with Eq. 5, and the constraint action-value Qfunction Q n c is evaluated with Eq. 6. Both can be trained with any supervised learning regression algorithm.

Q n+1 r (s i , a i , β i ) reg ← -- (5)
r i + γ a ∈A π n (s i , a , β i )Q n r (s i , a , b n (s i , a , β i )) Q n+1 c (s i , a i , β i ) reg ← -- (6)
c i + γ c a ∈A π n (s i , a , β i )Q n c (s i , a , b n (s i , a , β i ))
Contrary to regular FTQ where the policy is a simple greedy exploitation of the Q-function, BFTQ has to compute its two policies π n and b n at each iteration. Since taking some actions could exceed the budget, both budgetpolicy and behavioural-policy are jointly optimised under the constraint of belonging to set Ψ, described in Eq. 7. This equation denotes the search space for admissible (behavioural-policy, budget-policy) couples. Eq. 8 describes the optimisation process.

Ψ n =        π ∈ ∆ S×R A , b ∈ R S×A×R , such that, ∀s ∈ S, ∀β ∈ R, a∈A π(s, a, β)Q n c (s, a, b(s, a, β)) ≤ β        (7) (π n+1 ,b n+1 ) ← (8) argmax (π,b)∈Ψ n+1 a∈A π(s, a, β)Q n+1 r (s, a, b(s, a, β))
The stopping criterion of BFTQ is classically defined as a threshold η on Q-functions parameters differences.

The main result of this training phase is the couple of fitted Q-functions Q r and Q c . As a side result we also obtain approximations of the optimal action and budget policies π i and b i for each transition of the training base.

Theoretical Policy Execution

Let β t and s t be respectively the remaining budget and state of the agent at time step t. The agent can theoretically use the two Q-functions Q r and Q c according to Eq. ( 8) to determines the optimal distribution π(s t , •, β t ) over the actions and the corresponding budget repartition b(s t , •, β t ). It then samples an action according to the distribution and replaces its current budget according to the budget repartition and the action sampled.

a t ∼ π(s t , •, β t ) (9) β t+1 = b(s t , a t , β t ) (10) s t+1 ∼ P (s t , a t , •)
The advantage of this approach is that it only approximates on Q r and Q c . Its downside is that one has to solve a nonlinear problem under constraints at each time step of the policies execution: for time-critic applications, its use would be disastrous.

Policies approximations

Having to solve Eq. ( 8) is painful for both the training and deployment phases. To overcome this issue, we approximate the functions π and b using regressions on the training transitions. The optimisation described in Eq. 8 does not need to be computed exhaustively: it has to be performed only for the transitions in the learning base. For this purpose, we use the sets ψ n i , as defined in Eq. ( 11)), which are the projections of the set ψ n of Eq. ( 7) on the training transitions.

Ψ n i =    π i ∈ ∆ A , b i ∈ R A , such that: a∈A π i (a)Q n c (s i , a, b i (a)) ≤ β i    (11)
The projection of Eq. ( 8) on the training set is formalized in Eq. 12.

π n+1 (s i , •, β i ),b n+1 (s i , •, β i ) reg ← -- (12)
argmax

(πi,bi)∈Ψ n+1 i a∈A π i (a)Q n+1 r (s i , a, b i (a))
At the end of the training process, we can drop the Qfunctions and work directly with the fitted policies as in Eq. ( 9) and Eq. ( 10).

Solving the non-linear problem

The main remaining challenge is to solve the non-linear problem Eq 12 during the training phase. One may use some black-box optimization algorithm like random sampling, SLSQP [START_REF] Kraft | Slsqp, a nonlinear programming method with quadratic programming subproblems[END_REF]) but it would be inefficient. Hopefully we can exploit some specific properties of the Q-functions to speed-up the optimization. Let Q s,a (β) denote the following parametric function:

Q s,a (β) → (Q c (s, a, β), Q r (s, a, β)) (13)
This function is concave and strictly increasing in β along the axis of Q r . Indeed, given a budget β, a state s and an action a, the expected reward and actual budget consumption for a given action, Q r (s, a, β) and Q c (s, a, β) will be lower or equal to the ones yield by a larger budget β = β + ∀ ≥ 0. The intuition is as follows: given a budget β , in worst case scenarios (when extra budget is useless), one can always achieve at least the same performances as with the lower budget β; indeed, the space of admissible policies for β contains the one for β.

We need to solve Eq 8 for a state s. For the moment, we assume to have 2 possible actions a ⊥ and a . The functions Q s,a ⊥ and Q s,a are plotted on Fig 4 . The optimal Q r value is necessarily obtained by the probabilist combination of those actions. This combination of actions yield values (Q c , Q r ) on a straight line. Since Q functions are concave and strictly increasing, the line defining the optimal and admissible solutions is actually a tangent of both Q s,a ⊥ and Q s,a , lying on of both functions as shown on Fig 4 . We need to find both intersections points between the tangent and the Q functions to define the optimal solution. Let

(s, , ) Q c a ⊥ b ⊥ Q c (s, , ) Q c a ⊤ b ⊤ (s, , ) Q r a ⊥ b ⊥ (s, , ) Q r a ⊤ b ⊤ Q r Q s,a ⊤ Q s,a ⊥ β y β
(Q c (s, a ⊥ , b ⊥ ), Q r (s, a ⊥ , b ⊥ )) and (Q c (s, a , b ), Q r (s, a , b
)) denote the two intersections points. The optimal budget-policy in s is then defined as b(s, a ⊥ , β) = b ⊥ and b(s, a , β) = b . The optimal behavioural-policy in s is the following:

π(s, a , β) + π(s, a , β) = 1, if β ≤ Q c (s, a ⊥ , b ⊥ ), π(s, a , β) = 0, ( 14
) if β ≥ Q c (s, a , b ), π(s, a , β) = 1, otherwise, π(s, a , β) = β -Q c (s, a ⊥ , b ⊥ ) Q c (s, a , b ) -Q c (s, a ⊥ , b ⊥ )
.

One can extend this reasoning to cases with more than two actions by considering actions two by two. Finally, solving Eq 12 actually amounts to solve the following non linear program: where

π n+1 (s i , •, β i ),b n+1 (s i , •, β i ) reg ← -- (15
Ω n i =                                    π i , b i , such that: (a ⊥ , a ) ∈ A 2 , (b ⊥ , b ) ∈ R 2 , b i (a ) = b , b i (a ⊥ ) = b ⊥ , b i (a) = 0 ∀a ∈ A \ {a , a ⊥ } Q n c (s i , a ⊥ , b ⊥ ) ≤ β i < Q n c (s i , a , b ), π i (a) = 0 ∀a ∈ A \ {a , a ⊥ } π i (a ⊥ ) = βi-Q n c (si,a ⊥ ,b ⊥ ) Q n c (si,a ,b )-Q n c (si,a ⊥ ,b ⊥ ) , π i (a ) = 1 -π i (a ⊥ ), ∆ = ∆ ⊥ = ∆                                    with ∆ = Q n r (s i , a , b ) -Q n r (s i , a ⊥ , b ⊥ ) Q n c (s i , a , b ) -Q n c (s i , a ⊥ , b ⊥ ) ∆ ⊥ = ∂Q r (s i , a ⊥ , b ⊥ ) ∂β / ∂Q c (s i , a ⊥ , b ⊥ ) ∂β ∆ = ∂Q r (s i , a , b ) ∂β / ∂Q c (s i , a , b ) ∂β
One way to solve Eq 15 is to compute the convex hull of discrete version of Q functions. We propose an algorithm named SOLVE in Alg 1 to do so and show an example of its execution in Fig 5 . We assume that the function hull returns the convex hull H of a set in clockwise order. It starts with the lowest value in Qc, then the corresponding budgets B and finally the corresponding actions A P 3 . We also need to operate a reasonable discretization of Q-function and reduce the search space. For this reason we need to know the maximal cumulated constraint the environment could return. We call it β max . it more reward, since the path to the reward is not that risky. We can also observe a clear disparity between the λ values and their corresponding budgets. For instance it seems quite difficult to guess the right λ between 100 and 150 in order to stay tightly under a budget of value 0.25 without relaunching simulations. On the contrary, the BFTQ(β) policy allows one to fix any budget bound without worrying about the reward design.

On the second experiment, we compare the two agents on 30 randomly generated 10 × 10 worlds. Results are displayed on Fig. 7. Once again performances between algorithms are similar, although FTQ(λ) reaches slightly better reward return than BFTQ(β) for the same constraint return. This is mainly explained by the approximation of the behavioural/budget policies and the discretization of Q-functions involved in the computing of the convex hull. We can also notice that, on Fig. 7, BFTQ(β) agent's performances are less noisy than FTQ(λ)'s performances. The explication is as follows: the budget choosing phase is invariable for any world while the λ used in FTQ(λ) to keep the constraint return under a particular value may be different for two worlds. This is why FTQ(λ) needs a calibration phase -which is impossible when environment can not be simulated -while BFTQ(β) only needs to specify a budget.

Finally, some examples of BFTQ(β)'s policies executions are shown in Fig. 8. For a low budget, β = 0.1, the agent stops its episode after 1 or 2 actions. For β = 0.2, the agent tries to reach the goal with more insistence but still stops some episodes prematurely (selecting dont move action). For β = 0.6, the agent does not stop until it reaches the goal.

Conclusion

We have introduced BFTQ, a direct extension of FTQ to solve BMDP problem. To scale things up, we proposed to learn the policy using supervised learning. We have also shown that in continuous BMDP, Q-functions are concave and increasing. We exploited this property to approximately solve more efficiently non-linear program of BFTQ through a dedicated algorithm that we named SOLVE. Finally, we tested the algorithm on continuous 2-D worlds: these experiments confirm that BFTQ performs as well as a penalized FTQ in terms of reward/constraint raised couple without the need of a simulated calibration phase.

We believe that classic single reward oriented reinforcement learning algorithms could be extended to solve BMDP in the same way. It may also be of interest to find a close form solution for the program in Eq 15. This would speed up considerably the computation time of BFTQ and avoid the need of approximating policies.

  );[START_REF] Chow | Risk-constrained reinforcement learning with percentile risk criteria[END_REF];[START_REF] Achiam | Constrained policy optimization[END_REF] propose algorithms for CMDP, but they require interactions with the environment during the training (actor-critic, Q-learning and Policy Gradient, Constrained-Policy-Optimization).[START_REF] Abe | Optimizing debt collections using constrained reinforcement learning[END_REF] work with constraint-based RL but their solution is a Q-learning algorithm adapted to batch learning, which is known not to be a sample-efficient batch RL algorithm.[START_REF] Thomas | High confidence policy improvement[END_REF];[START_REF] Petrik | Safe policy improvement by minimizing robust baseline regret[END_REF];[START_REF] Laroche | Safe Policy Improvement with Baseline Bootstrapping[END_REF] introduced batch RL algorithms for risk sensitive problems but these algorithms were not adapted for continuous environments.

Figure 2 :

 2 Figure2: Calibration of a penalty multiplier according to the budget β. The optimal multiplier λ * avg is the smallest one to satisfy the budget constraint on average. Safer policies can also be selected according to the largest deviation form this mean cost.

Figure 3 :

 3 Figure 3: A simple deterministic finite BMDP.

Figure 4 :

 4 Figure 4: An example of Q-functions as defined in Eq 13. The dot lines are tangent to both curves. The tangent at the right defines the optimal policy when the budget β lies between Q c (s, a ⊥ , b ⊥ ) and Q c (s, a , b ).The value y β is then equal to π(s, a ⊥ , β)Q r (s, a ⊥ , β) + π(s, a , β)Q r (s, a , β).

Figure 5 :

 5 Figure 5: An execution of SOLVE with three actions. Each dots correspond to discretization w.r.t β. Red dots are considered to compute the convex hull. Blue dots are removed since they reach a worst value of Q r for greater Q c . Finally, pink cross are the dots retained to compute the optimal policy.

  Figure 6: BFTQ(β) and FTQ(λ) reward w.r.t the constraints. Values near the dots are the budgets β allowed for BFTQ and the value λ for FTQ. Each dot is the mean on 9 runs on a 5x5 2-D worlds with 1000 test episodes on each of them. The standard deviation plotted refer to the variability between all test episodes.

Figure 7 :Figure 8 :

 78 Figure 7: Each dot is the mean on 30 2-D worlds with 1000 test episodes on each of them. The standard deviation plotted refer to the variability between the random 2-D worlds.

To the best of our knowledge, this separability assumption is still a conjecture for general BMDPs (with loops).

For example, if the cumulative cost of a problem is in interval [1,

3], we can takeβ in {1 + n/100 | 0 ≤ n ≤ 200}.

Algorithm 1 SOLVE

In: s, Qr, Qc, K, β, βmax Out: π, b P = ∅ for a in A do for i = 0 to K do step = i/K * βmax p = (Qc(s, a, step), Qr(s, a, step)) P = P ∪ p end for end for H, B, AP = hull(P) f ound = F alse i = 0 while ¬ f ound do qc, qr = H the next section, we evaluate BFTQ on a set of 2-D worlds.

Experiments

We do experiments on N × N 2-D worlds with continuous state-space. The agent starts in top-left corner (0, 0). The goal is to reach the bottom-right corner, at coordinates (N -1, N -1), which ends the episode with a immediate reward of 1. On its way to the goal, the agent may fall into randomly placed holes. In that case, it also ends the episode with a constraint of magnitude 1. Consequently, constraints can be seen as failing probability. The same goes for rewards as success probability.There are three actions available: move right, move bottom and dont move. The two first actions yield stochastic transitions: a Gaussian noise of mean 0 and standard deviation σ is applied. Last action ensures that there always exists a solution for any starting state if β is zero.

In the following experiments, we use two different RL agents: BFTQ(β) uses BFTQ algorithm to learn its policy. Q r , Q c , π and b are approximated with regression trees. The algorithm used for computing the convex hull is the Graham scan [START_REF] Graham | An efficient algorithm for determining the convex hull of a finite planar set[END_REF]). β is the budget used when initiating a new episode. Note that all BFTQ(β) agents use Eq. 5 and 6. The second agent, FTQ(λ), uses FTQ with a regression tree. The penalty λ is used when rewarding a fall into a hole : r ← r -λ * c. The training dataset is uniformly distributed over the state-action sets.

On the first experiment, we compare BFTQ(β) with FTQ(λ) on a single 5 × 5 world. Results are displayed on Fig. 6. For a given constraint return, both agents achieve a similar reward return. In this setting, BFTQ(β) is slightly optimistic: until β = 0.3, there is an offset of 0.05 between the given budget and the constraint return. So BFTQ(β) may not fully respect the given budget. But the same phenomena occurs for FTQ(λ): when λ = 10 4 , we may think that the algorithm should return a constraint of value 0, but it is actually 0.07. This is because of the approximations involved in both algorithms. Note that when the constraint return exceeds 0.4, failing probabilities are capped as well as rewards. The explication is as follows: at some point, the agent had to take some risk to get the reward. But taking more risk won't necessary give