Nicola Zaupa

Carlos Olalla

Isabelle Queinnec

Luis Martínez-Salamero

Luca Zaccarian

Hybrid Control of Self-Oscillating Resonant Converters with Three-Level Input

Keywords: Stability of hybrid systems, Lyapunov methods, Power electronics

We propose a hybrid feedback law inducing self-oscillating behavior in second-order resonant converters. With our controller, the converter switches at the resonant frequency of its tank, without the need of external oscillators. In addition, the output amplitude can be adjusted by a reference signal ranging from zero to π /2. The amplitude modulation is then performed while maintaining an approximately constant switching frequency. Theoretical results show uniqueness and almost global asymptotic stability of a nontrivial hybrid limit cycle. Experimental results show that a circuit implementing the new controller successfully matches the desirable simulated behavior.

I. INTRODUCTION

C ONTROL of resonant converters is a well-known prob- lem in power electronics associated with several solutions [START_REF] Erickson | Fundamentals of Power Electronics[END_REF]- [START_REF] Youssef | A review and performance evaluation of control techniques in resonant converters[END_REF]. Conventional approaches include frequency modulation, where the switching frequency is the control input [START_REF] Vorperian | Approximate small-signal analysis of the series and the parallel resonant converters[END_REF], [START_REF] Yin | Small-signal analysis of frequency-controlled electronic ballasts[END_REF], and amplitude modulation, where the control input is the phase-shift between the two legs of the switches [START_REF] Burdio | Comparison study of fixed-frequency control strategies for ZVS DC/DC series resonant converters[END_REF], [START_REF] Youssef | Self-sustained phase-shift modulated resonant converters: Modeling, design, and performance[END_REF]. Frequency modulation may induce desirable zero voltage switching (ZVS) but the regulation of the operating point may require wide frequency ranges. Conversely, amplitude modulation desirably operates at constant switching frequency, at the expense of losing ZVS.

This paper adopts a state-plane hybrid control approach inducing self-oscillation of the resonant converter without the need of external inputs. State-plane based methods that turn to self-oscillating schemes have been proposed in the past, either ensuring the oscillations regardless of the load condition [START_REF] Youssef | Self-sustained phase-shift modulated resonant converters: Modeling, design, and performance[END_REF], [START_REF] Pinheiro | Performance characterization of two self-oscillating controllers for parallel resonant converters operating with unity input power factor[END_REF], or proposing state-plane-induced input selection [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF], [START_REF] Bonache-Samaniego | Dynamic Modeling and Control of Self-Oscillating Parallel Resonant Converters Based on a Variable Structure Systems Approach[END_REF]. Our approach has several advantages with respect to these previous papers. First, we relax the conservative requirement imposed by [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF], [START_REF] Bonache-Samaniego | Dynamic Modeling and Control of Self-Oscillating Parallel Resonant Converters Based on a Variable Structure Systems Approach[END_REF] on the quality factor. Secondly, the proposed controller has useful links with the conventional control methods, as explained next.

This work was supported in part by: the Spanish Ministry of Science and Innovation under projects PID2019-111443RB-100 and PID2020-120151RB-I00; the ANR via grant HANDY, number ANR-18-CE40-0010; and the MUR under PRIN grant DOCEAT, number 2020RTWES4. (Corresponding author: Nicola Zaupa (e-mail: nzaupa@laas.fr).)

Nicola Zaupa, Isabelle Queinnec and Luca Zaccarian are with LAAS-CNRS, University of Toulouse, CNRS, Toulouse, France. Luca Zaccarian is also with Universit à di Trento, 38123 Trento, Italy.

Carlos Olalla and Luis Martínez-Salamero are with the Departament d'Enginyeria Electr ònica, El èctrica i Autom àtica, Escola T ècnica Superior d'Enginyeria, Universitat Rovira i Virgili, 43007 Tarragona, Spain. In our previous works [START_REF] Zaupa | Results on hybrid control of self-oscillating resonant converters[END_REF], [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF], we addressed the control of second-order parallel and series resonant converters (PRC and SRC, respectively) with the hybrid formulation of [START_REF] Goebel | Hybrid Dynamical Systems: modeling, stability, and robustness[END_REF] in a unified way. To represent the switching input voltage v s of the resonant tank we used there a discrete variable σ ∈ {-1, 1}. The hybrid control law allowed regulation of the converter by changing the operating frequency. Thus, the hybrid controller provided a self-oscillating frequency modulation behavior.

Given the limitations of using wide frequency ranges in resonant conversion, the objective of this paper is to propose and to validate a hybrid controller that allows modulation of the amplitude of the first harmonic of the input voltage v s . The hybrid controller produces, in a way, a self-oscillating amplitude modulation behavior. We introduce here a third value for σ ∈ {-1, 0, 1} corresponding to a zero voltage input, allowing us to operate at the resonant frequency of the tank across a range of output amplitudes.

The fact of having a three-level input voltage is beneficial in two ways: modulation of input voltage first harmonic is possible and the harmonic content is reduced. The last point is supported by [START_REF] Nabae | A New Neutral-Point-Clamped PWM Inverter[END_REF] since PRC and SRC are topologically similar to three-level neutral point clamped inverters (NPC) feeding a lossless LC two-port terminated at a resistive load. On the other hand, NPC requires an external sinusoidal reference tuned at the proper frequency. While, in our three-level input voltage SRC and PRC it is not required because the hybrid control induces an oscillation at the resonant frequency.

The paper is structured as follows. The system model and the proposed controller are discussed in Section II. Section III states and proves our main theoretical result and Section IV reports the experimental implementation and validation on a power converter prototype. Section V conclude the work.

Notation. R (R >0 ) [R ≥0 ] and Z (Z >0 ) [Z ≥0
] are the sets of (positive) [non-negative] real and integer numbers. R n denotes the n-dimensional Euclidean space. Given two vectors u ∈ R n and w ∈ R m , u ⊤ denotes the transpose of u, and (u, w) := u ⊤ w ⊤ ⊤ denotes their stacking. Given a (continuous, discrete, or hybrid) signal x, ẋ denotes its derivative with respect to continuous time t, while x + denotes

z 2 z 1 D1 D-10 z 2 z 1 1 1 D1 D10 φ z 2 z 1 φ D-1 D-10 φ z 2 z 1 φ 1 1 D-1 D10 φ φ 1 1 1 C 1 C 10 C-10 C -1 (σ, d) = (1, 1) (σ, = (0, 1) (σ, d) = (-1,-1) (σ, d) = (0,-1) 1 ξ ḡ(D1) λ(ξ) ξ λ(ξ)
Fig. 2. Jump sets (solid lines) and flow sets (sectors) for each selection of (σ, d). A sample trajectory is represented in black. Coordinate λ will be introduced and used in the proof.

its next value with respect to discrete time j.

II. SYSTEM MODEL AND PROPOSED CONTROLLER

We address parallel and series second-order resonant converters (PRC and SRC, respectively) as shown in Fig. 1. As in [START_REF] Zaupa | Results on hybrid control of self-oscillating resonant converters[END_REF], [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF], consider the following change of coordinates

z 1 := v C V g -σ, z 2 := 1 V g L C i C , (1)
based on the voltage and current on the capacitor (v C and i C ) and the discrete variable σ ∈ {-1, 0, 1} representing the switching input voltage v s of the resonant tank, so that v s = σV g , where V g is the external DC supply. Following the derivation in [11, eqs. (1) to (4)], the continuous-time dynamics of both converters is unified to

ż1 = ωz 2 , ż2 = -ωz 1 -βz 2 , (2)
where

ω := √ LC -1
is the natural frequency and β > 0 is the inverse of the time constant of the exponential decay associated with each resonant tank of Fig. 1:

β PRC := 1 RC , β SRC := R L . (3)
This novel approach to control (1) includes the possibility to apply a zero voltage input σ = 0. While this extra value does not change coordinates (1) nor their continuous dynamics (2), it radically changes the switching possibilities, which are here described by a hybrid model enforcing a precise periodic sequence for σ: 1, 0, -1, 0, . . . (in practice σ is a three-level square waveform):

    ż1 ż2 σḋ     = f (ξ) =     ωz 2 -ωz 1 -βz 2 0 0     ξ ∈ C(φ), (4a)
    z + 1 z + 2 σ + d +     = g(ξ) =     z 1 + d z 2 -d + σ -d + 2σ     ξ ∈ D(φ). (4b)
In our new model ( 4), ξ := (z 1 , z 2 , σ, d) collects the four states: two physical and two logical ones. The logical state σ denotes the switch position, as in [START_REF] Zaupa | Results on hybrid control of self-oscillating resonant converters[END_REF], [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF]. The new logical state d ∈ {-1, 1} plays an important role when σ = 0 as it stores the past value of σ, so that the next value will be assigned as the opposite one (i.e. -d), in order to generate the correct switching pattern. Clearly, any variation of σ and d must be instantaneous, so that σ = 0 and ḋ = 0.

The angle φ ∈ [0, π /2) is a reference input affecting the steady-state output voltage (minimum for φ = π /2 and maximum for φ = 0). Based on φ, we select the jump and flow sets as

C(φ) := C 1 (φ) ∪ C 10 (φ) ∪ C -1 (φ) ∪ C -10 (φ), D(φ) := D 1 (φ) ∪ D 10 (φ) ∪ D -1 (φ) ∪ D -10 (φ), (5a)
which correspond to the union of four jump and flow sets graphically illustrated in Fig. 2. Selection (5a) intuitively leads to a reduced phase mismatch between v s and i C , as shown by the sample solution of (4), ( 5) represented in black. More precisely, for each q ∈ {-1, 1}, define

C q (φ) := {(z, σ, d) : σ = q, d = q, d ((z 1 + σ) sin φ -z 2 cos φ) ≤ 0}, (5b) C q0 (φ) := {(z, σ, d) : σ = 0, d = q, (d (z 1 sin φ -z 2 cos φ) ≥ 0) ∨ (dz 2 ≥ 0)}, (5c) D q (φ) := {(z, σ, d) : σ = q, d = q, dz 2 ≥ 0, (z 1 + σ) sin φ -z 2 cos φ = 0}, (5d) D q0 (φ) := {(z, σ, d) : σ = 0, d = q, dz 2 ≤ 0, z 1 sin φ + z 2 cos φ = 0}. ( 5e
)
Note that D 1 (φ) and D -1 (φ) are horizontally translated with respect to the origin. Therefore, when σ = 0, solutions flow by spanning a constant cone centered at the origin, so that the flowing interval is constant regardless of the amplitude of z. Secondly, the zero crossing of the current i C (proportional to z 2 ) happens approximately at the middle of that flowing interval (the lower β the better this approximation, due to the higher quality factor). This results in keeping v s and i C mostly in phase and, for the SRC configuration, it ensures that the input voltage and the current of the resonant tank are almost synchronized.

Remark 1 Shifting horizontally the two jump sets D 1 and D -1 has a desirable effect on the relation between φ and the output amplitude. First, we can reach zero amplitude for φ = π /2. Secondly, the static relation between φ and •

III. STABILITY OF THE RESONANT BEHAVIOR

A. Main result

For our main stability result we require a mild underdamped assumption (see [12, Remark 1]) for the oscillating behavior, requiring that flowing solutions revolve in the phase-plane of Fig. 2, which ensures existence of a non trivial asymptotically stable hybrid limit cycle (as defined in [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF] and recalled below) for any reference input φ ∈ [0, π /2).

Assumption 1

The relation 2ω > β > 0 is satisfied, namely the resonant tank is underdamped. Equivalently, the roots of s 2 + βs + ω 2 = 0 are complex conjugate.

Definition 1 [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF] Given an hybrid system H = (C, f, D, g), a nontrivial hybrid periodic trajectory ζ is a complete solution (namely, a solution that evolves forever) that is not identically zero and for which there exists a pair (T, J) ∈ R ≥0 × Z ≥0 satisfying T + J > 0, such that (t, j) ∈ dom(ζ) implies (t + T, j + J) ∈ dom(ζ) and, moreover,

ζ(t, j) = ζ(t + T, j + J). ( 6
)
The image of ζ is a nontrivial hybrid periodic orbit.

We state next the main theoretical result of this work. It states the asymptotic stability of a nontrivial hybrid periodic orbit with basin of attraction corresponding to the whole space except for the four points (see the green dots in Fig. 4) in the set K := {(0, 0, 0, -1), (0, 0, 0, 1), (-1, 0, 1, 1), (1, 0, -1, -1)}. These points comprise weak equilibria1 and they do not belong to the basin of attraction: for example, solutions starting from (0, 0, 0, -1) and (0, 0, 0, 1) can flow forever staying there or may jump to (-1, 0, 1, 1) or (1, 0, -1, -1) and then: either flow and converge to the hybrid limit cycle, or continue to jump among the points in K.

Theorem 1 Under Assumption 1, for each selection of φ ∈ [0, π /2), the closed loop (4), ( 5) has a unique nontrivial hybrid periodic orbit O φ that is stable and almost globally attractive with basin of attraction Ξ = (R 2 × {-1, 0, 1} × {-1, 1}) \ K.

B. Preliminary observations to prove Theorem 1

Let us define the set Ω = Ω 1 + Ω -1 + Ω 10 + Ω -10 as the set of points visited by any solution flowing in C from g(D). More rigorously, denoting t → ψ ξ0 (t, 0) any flowing solution of (4) from ξ 0 , for each i ∈ {1, -1, 10, -10}, we define

Ω i = {ξ ∈ C : ξ = ψ g(ξ0) (t, 0) for some ξ ∈ D i , t ≥ 0}.
Sets Ω 1 and Ω -1 correspond to two cones centered at the origin. Set Ω 10 (resp. Ω -10 ) corresponds to a cone minus an area limited by what we may call boundary solution, which flows from the cone vertex z v = (1, 0) (resp. z v = (-1, 0)). About this boundary solution, it is well defined for any φ ∈ [0, π /2) since the vector field at the vertex is vertical, therefore the solution is forced to initially flow inside the cone and then, due to Assumption 1, it revolves clockwise while flowing, until it reaches D -1 (resp. D 1 ). We define the scalar λ m > 0 as the distance between the vertex (1, 0) (resp. (-1, 0)) and the point where the boundary solution intersects D -1 (resp. D 1 ) after flowing in the cone.

Fig. 4 depicts λ m , the set Ω, the set P := (D 1 ∪ D -1 ) ∩ Ω, and the coordinate λ ∈ [λ m , +∞) defined as

λ = Λ(ξ) := |[ 1 0 0 0 0 1 0 0 ] g(ξ)| , ∀ξ ∈ P, (7)
see also the first and third diagram in Fig. 2.

The following lemma is an immediate property of solutions and is stated without proof.

Lemma 1 Set Ω is forward invariant and all solutions from Ω are unique. Moreover, all solutions from (C ∪ D) \ K converge to Ω in uniform finite time.

Introduce now the function ϕ P (ξ 0 ) : Ω → P returning the value ψ ξ0 (t P , j P ) ∈ P, where ψ ξ0 is the solution to (4) starting at ξ 0 (unique by Lemma 1), and (t P , j P ) is

(t P , j P ) = argmin (t,j)∈dom ψ ξ 0 (t + j) s.t. ψ ξ0 (t, j) ∈ P. (8)
Intuitively speaking, for any ξ 0 ∈ Ω, ϕ P (ξ 0 ) denotes the value of the unique solution starting at ξ 0 the first time that it belongs to P.

C. Proof of Theorem 1

We first prove existence and global asymptotic stability (GAS) from Ω of a hybrid periodic orbit O φ contained in the interior of Ω, based on a Poincaré analysis. Then by the uniform convergence established in Lemma 1, the result follows.

After the first jump, all solutions must transit through set P. Therefore, we may define a Poincaré map characterizing a solution half-cycle as ξ + = ϕ P (g(ξ)) ∈ P for any ξ ∈ P. Due to the central symmetry of the dynamics, we can also parameterize the Poincaré map in terms of λ, as defined in [START_REF] Youssef | Self-sustained phase-shift modulated resonant converters: Modeling, design, and performance[END_REF], through description

λ + = Φ(λ), λ ∈ [λ m , +∞), (9)
with Φ(λ) := Λ(ϕ P (ξ λ )), and Λ -1 (λ) denoting a selection of ξ λ ∈ P ensuring Λ(ξ λ ) = λ. The discrete dynamics (9) plays the role of a Poincaré map for (4) restricted to Ω. We show below that there exists a unique globally asymptotically stable equilibrium λ ⋆ = Φ(λ ⋆ ) for [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF]. Energy function. Due to Lemma 1, we consider next only points ξ ∈ P. Then, the proof of existence and uniqueness focuses on the behavior of the z components of the unique solution to (4) inducing Φ, in terms of the energy function

E(z) := 1 2 |z| 2 . ( 10
)
Due to this fact, to simplify our notation, in the sequel, we use g to denote the z component of the jump map g. Also, we will employ overlines on the jump and flow sets, to denote their projection on the z plane. With a slight abuse of notation, we will only refer to the z component of the solutions, whenever the σ and d components are clear from the (flow or jump) set under consideration.

Let φ ∈ [0, π /2) be fixed and, for each λ ≥ λ m , consider any pair of points z F (λ) ∈ g(D 1 ) and z B (λ) ∈ g(D -1 ) both of them at the same distance λ from the origin, where "F " stands for "Forward" and "B" stands for "Backwards". A pictorial representation of these points and related solutions is shown in Fig. 4. Due to the symmetry of the jump sets, we have Forward solution. About z F (λ) and ζ F (λ), the forward solution flows in C 10 and then jumps after it reaches D 10 (which necessarily happens after a finite time, due to the revolving nature of flowing solutions, as induced by Assumption 1). In particular, since C 10 is a cone, regardless of λ, the solution flows for a fixed time T and the energy dissipation during such a flowing interval corresponds to Σ F (λ) = 1 2 λ 2 (1 -µ(φ) 2 ), where µ(φ) (the dependence on φ is dropped in the following) is a constant that can be easily computed from the flowing solution of ż = Az, by solving the equality e AT λ cos φ sin φ = λµ

E(z F (λ)) = E(z B (λ)) = 1 2 λ 2 . Consider the points ζ F (λ) ∈ g(D 10
δ F (λ) = 1 2 |(e AT z F ) + | 2 -|e AT z F | 2 = 1 2 + λµ cos φ. (11
) Overall, we may characterize the energy at point ζ F (λ) as

E F (λ) := E(ζ F (λ)) = 1 2 λ 2 -Σ F (λ) + δ F (λ). ( 12
)
Backward solution. Consider now z B (λ) and ζ B (λ), and also here notice that the backward solution of interest jumps backwards to a point in D -1 denoted as z B (λ

) -= λ -cos φ -sin φ +[ 1 0 ] = 1-λ cos φ -λ sin φ .
The variation of energy across this backward time jump corresponds to

δ B (λ) = 1 2 |z B (λ) -| 2 -|z B (λ)| 2 = 1 2 -λ cos φ. ( 13
)
Continuing in the backwards motion of this solution, for any λ ≥ λ m , we see a flowing interval in the set C -1 where the solution reaches (in finite backwards time) the set g(D 10 ), due to the revolving properties induced by Assumption 1.

Applying [12, Lemma 2] such a flowing interval is associated with an energy increase (in backwards time) corresponding to Σ B (λ) := β ω Π(z B (λ) -), where Π(z B (λ) -) is the (unsigned) area hatched between the graph of the backwards trajectory and the axis z 2 = 0. Note that Σ B is a positive and quadratically increasing function of λ.

Overall, we may characterize the energy at point ζ B (λ) as

E B (λ) := E(ζ B (λ)) = 1 2 λ 2 + Σ B (λ) + δ B (λ). (14)
Energy balance. Consider the mismatch function

E(λ) := E F (λ) -E B (λ) = δ(λ) -Σ(λ), (15)
where Σ is a quadratically increasing function of λ (it is strongly convex):

Σ(λ) := Σ B (λ) + Σ F (λ) = Σ B (λ) + 1 2 λ 2 (1 -µ 2 ), (16)
and δ is a linearly increasing function of λ:

δ(λ) := δ F (λ) -δ B (λ) = λ(1 + µ) cos φ. (17)
Lemma 2 There exist a unique value λ ⋆ > λ m such that E(λ ⋆ ) = 0. Moreover, λ ⋆ is a fixed point for Φ in [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF].

Proof: Let us show that E(λ m ) = δ(λ m ) -Σ(λ m ) > 0, which from (15) is equivalent to showing that E F (λ m ) - E B (λ m ) > 0. For λ = λ m , a) the forward solution reaches ζ F (λ m ) with energy E F (λ m ) = 1 2 + 1 2 λ 2 m µ 2 + λ m µ cos φ; b) by definition of λ m , the backward solution reaches the vertex ζ B = (1, 0) with energy E B (λ m ) = 1 2 . Therefore E F (λ m ) -E B (λ m ) = 1 2 λ 2 m µ 2 + λ m µ cos φ,
which is the sum of two positive scalars. Since function λ → ψ(λ) := -E(λ) is negative at λ m and strongly convex, then, from [12, Lemma 1], there exists a unique positive value λ ⋆ > λ m such that E(λ ⋆ ) = 0. Finally, we can express the energy difference, in terms of λ, between the final and initial points Φ(λ) and λ, as

E(λ) = 1 2 (Φ(λ) 2 -λ 2 )
. From which we conclude that λ ⋆ is a fixed point for [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF].

The corresponding points ξ ⋆ ∈ P that satisfy Λ(ξ ⋆ ) = λ ⋆ , and the image of the unique solution originating from those points is a unique nontrivial hybrid periodic orbit in the interior of Ω (because λ ⋆ > λ m ), characterized as

O φ = {ξ ∈ Ω : Λ(ϕ P (ξ)) = λ ⋆ }. ( 18
)
Lyapunov function.

In order to prove GAS of the hybrid periodic orbit we reduce the stability analysis to the Poincaré map [START_REF] Bonache-Samaniego | Design of self-oscillating resonant converters based on a variable structure systems approach[END_REF]. We have already shown that there exists a unique equilibrium point λ ⋆ = Φ(λ ⋆ ). Then, we introduce the quadratic Poincaré-like Lyapunov function,

V (λ) := (λ -λ ⋆ ) 2 . ( 19
)
Similar to the final part of the proof of [START_REF] Bisoffi | Relay-based hybrid control of minimal-order mechanical systems with applications[END_REF]Lemma 2], to prove global asymptotic stability, we show that:

V (Φ(λ)) -V (λ) < 0 ∀λ ≥ λ m , λ ̸ = λ ⋆ , or, equivalently, (Φ(λ) -λ)(Φ(λ) + λ -2λ ⋆ ) < 0 ∀λ ≥ λ m , λ ̸ = λ ⋆ . (20)
To show this, let us divide the study for λ > λ ⋆ and λ < λ ⋆ . From the proof of Lemma 2, for λ > λ ⋆ , we have that E(λ) < 0, which implies Φ(λ) -λ < 0, in other terms λ > Φ(λ) > λ ⋆ where Φ(λ) > λ ⋆ comes from the uniqueness of flowing solutions so that they cannot cross. Next, we assess the positivity of the second term in (20) as Φ(λ) + λ -2λ ⋆ > 2(Φ(λ)-λ ⋆ ) > 0. For the case λ < λ ⋆ , similar considerations hold, leading to (20). Finally, continuity of Φ follows from similar reasonings to [15, Lemma 2] and the stability of λ ⋆ for (9) proved above, imply global asymptotic stability of O φ in (18) from Ω, following the steps in [15, Proposition 1], which are here omitted due to length restrictions. Finally, GAS of O φ from Ξ follows from the uniform convergence established in Lemma 1.

IV. EXPERIMENTAL VALIDATION

The results of Theorem 1 are validated on an experimental SRC prototype and via simulations. The experimental setup consists in an SRC designed for a 50 kHz resonant frequency. The DC power supply is V g = 24 V, the resonant tank components are L = 94.5 µH and C = 100 nF. The control law is implemented in an FPGA where we use ADCs to sense the current and the voltage on the capacitor (needed to compute (z 1 , z 2 ) as illustrated in [START_REF] Zaupa | Results on hybrid control of self-oscillating resonant converters[END_REF]). The control law is implemented in the FPGA as an automaton with four states, corresponding to the four parts in which one period can be decomposed, based on the variables σ and d: (σ, d) ∈ {(1, 1), (0, 1), (-1, -1), (0, -1)}. (see Fig. 2). The jump sets are modeled as half-planes to ensure robustness and simplicity in the implementation. The automaton is allowed to go to the next sequential state whenever 1) jump signal is active and 2) the automaton is in the preceding state (e.g (σ, d) = (1, 1) and D 1 active: (z 1 + 1) sin φ -z 2 cos φ ≥ 0). Differently from [START_REF] Zaupa | Hybrid control of self-oscillating resonant converters[END_REF], to evaluate the jump conditions, a generic scaled version of (z 1 , z 2 ) is not enough: we need to know the scaling factor in the acquisition chain and compensate for it. This is due to the sets D 1 and D -1 being shifted with respect to the origin.

Multiple experiments have been run for two load selections: R = 10.1 Ω and R = 21.8 Ω. Fig. 5 compares the simulated and experimental steady-state responses for different values of φ. We emphasize that no startup controller is needed due to the (almost) global attractivity established in Theorem 1, thus the device starting at rest effectively spirals out to reach the pre-scribed oscillating behavior. The corresponding experimental time waveforms are reported in Fig. 6. Note that, the input current (i.e. the capacitor current, due to the SRC configuration) is approximately in phase with the input voltage. Having the first harmonics aligned implies improved performance in terms of efficiency, since more power is transmitted to the load and less power is dissipated in the resonant tank. Finally, Figs. 7 and 8 compare the simulated and experimental frequencies and amplitudes. The black horizontal dotted lines show the resonant frequencies evaluated by characterizing the resonant tank, with the load, at low power. The physical controller presents some limits for low amplitudes (large φ) due to the resolution of the ADCs and the amplitude of the noise. Therefore, only a subset of values for φ has been validated experimentally.

The experiments confirm the simulation results. The frequency mismatch visible in both figures amounts to an acceptable 4%. Minor amplitude discrepancies are also visible, probably due to model uncertainties (e.g. nonlinear dependence of the inductor on the current, parasitic effect on the PCB, sampling effects) but the reliable closed-loop behavior illustrates the robustness of the proposed control.

V. CONCLUSIONS AND FUTURE WORK

We have proposed a novel state-plane control approach modulating the first harmonic of the input voltage as an emergent self-oscillating behavior: an amplitude modulation approach. Through a reference input φ ∈ [0, π /2) we can control the amplitude of the converter output and we can always generate a self-oscillating hybrid limit cycle in underdamped conditions. The expected theoretical and simulated response has been confirmed by experiments on an SRC prototype.

The proposed hybrid controller opens interesting possibilities for future research. First, it can be anticipated that using the reference input for regulation of the converter will present advantageous dynamic properties. Secondly, novel conditions on the hybrid control law ensuring zero voltage switching would reduce the switching losses and open new possibilities

Fig. 1 .

 1 Fig. 1. Parallel resonant converter and series resonant converter.

Fig. 3 .

 3 Fig.3. Amplitude of the first harmonic of z 2 with respect to φ ∈ [0, π /2). Interpolating cosine function (red, solid) vs simulations results (blue * ).

Fig. 4 .

 4 Fig.4. Projection on the z-plane of the sets: Ω (shaded areas), P (solid red lines), K (green dots) and representation of the boundary solutions defining λm (black dotted arcs). The dashed-dotted lines represent sample solutions used in the proof.

  ) and ζ B (λ) ∈ g(D 10 ) obtained, respectively, from the forward flowing solution (unique from Lemma 1) starting at z F (λ) and the backward solution starting at z B (λ), jumping back to the unique point z B (λ) -∈ D -1 and then (uniquely) flowing backwards in time in C -1 until they both reach g(D 10 ).

-

  sin φ , resulting in µ = cos φ -sin φ ⊤ e AT cos φ sin φ ∈ (0, 1]. After this flowing interval, the solution belongs to D 10 and must jump (from e AT z F (λ) = λµ cos φ sin φ to ζ F (λ)). This jump is associated with an increase δ F of the energy E, amounting to

Fig. 5 .Fig. 6 .

 56 Fig. 5. Phase portraits in the (z 1 , z 2 ) coordinates: simulations (blue circles) and experiments (orange solid line) with different values of φ. The φ-tilted dotted black lines represents the jumps sets. The load is R = 10.1 Ω.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Simulated and experimental responses with R = 10.1 Ω. The black dotted line is the estimated resonant frequency f0 = 50 kHz.

A weak equilibrium is a point from which there exist constant solutions but also non-constant ones.


Brian Staber
email: <brian.staber@safrangroup.com>

Sébastien Da Veiga
email: <sebastien.da-veiga@ensai.fr>.

Benchmarking Bayesian neural networks and evaluation metrics for regression tasks

Due to the growing adoption of deep neural networks in many fields of science and engineering, modeling and estimating their uncertainties has become of primary importance. Despite the growing literature about uncertainty quantification in deep learning, the quality of the uncertainty estimates remains an open question. In this work, we assess for the first time the performance of several approximation methods for Bayesian neural networks on regression tasks by evaluating the quality of the confidence regions with several coverage metrics. The selected algorithms are also compared in terms of predictivity, kernelized Stein discrepancy and maximum mean discrepancy with respect to a reference posterior in both weight and function space. Our findings show that (i) some algorithms have excellent predictive performance but tend to largely over or underestimate uncertainties (ii) it is possible to achieve good accuracy and a given target coverage with finely tuned hyperparameters and (iii) the promising kernel Stein discrepancy cannot be exclusively relied on to assess the posterior approximation. As a by-product of this benchmark, we also compute and visualize the similarity of all algorithms and corresponding hyperparameters: interestingly we identify a few clusters of algorithms with similar behavior in weight space, giving new insights on how they explore the posterior distribution.

Introduction

Due to their recent achievements in the last decade, deep neural networks have been widely adopted in many research fields and industries. However, they still suffer from several shortcomings that prevent their deployement in fields where decisions involve high stakes. These limitations are mainly due to their inability to provide uncertainty estimates which are crucial for many real-world applications. Uncertainty quantification in deep learning has attracted a lot of attention and several approaches have been investigated. Thorough overviews are provided by the recent review papers of [START_REF] Abdar | A review of uncertainty quantification in deep learning: Techniques, applications and challenges[END_REF] and [START_REF] Gawlikowski | A survey of uncertainty in deep neural networks[END_REF]. Amongst the available approaches, Bayesian neural networks offer a simple and flexible formulation but raise several challenges. The high-dimensional posterior distribution of the network parameters is intractable, possibly multimodal, and the influence of the prior distribution remains an open question. The posterior distribution is typically approximated using sampling methods, variational inference, or Gaussian approximations. When dealing with complex posterior distributions, gradient-based sampling methods such as Markov Chain Monte Carlo (MCMC) methods are often adopted. In particular, Hamiltonian Monte Carlo (HMC) is usually considered as a gold standard sampling algorithm but is unfortunately extremely computationally demanding [START_REF] Izmailov | What are bayesian neural network posteriors really like[END_REF][START_REF] Cobb | Scaling hamiltonian monte carlo inference for bayesian neural networks with symmetric splitting[END_REF]. Since the work of [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF], stochastic gradient MCMC methods have been extensively studied (see, e.g., [START_REF] Ma | A complete recipe for stochastic gradient mcmc[END_REF] and [START_REF] Nemeth | Stochastic gradient markov chain monte carlo[END_REF]), where the Metropolis-Hastings correction is omitted and a mini-batched stochastic gradient is used. Variational approaches approximate the posterior distribution by minimizing the Kullback-Leibler divergence over a family of tractable distributions [START_REF] Hoffman | Stochastic variational inference[END_REF]. Several scalable methods have been proposed such as the Bayes by Backprop (BBB) method of [START_REF] Graves | Practical variational inference for neural networks[END_REF] and [START_REF] Blundell | Weight uncertainty in neural network[END_REF], the multiplicative normalizing flows proposed by [START_REF] Louizos | Multiplicative normalizing flows for variational bayesian neural networks[END_REF], and the Monte Carlo dropout method [START_REF] Gal | Dropout as a bayesian approximation: Representing model uncertainty in deep learning[END_REF]. Although ensemble methods are generally used to improve the generalization capabilities of neural networks, they can also be interpreted as a Bayesian approach [START_REF] Lakshminarayanan | Simple and scalable predictive uncertainty estimation using deep ensembles[END_REF][START_REF] Fort | Deep ensembles: A loss landscape perspective[END_REF], where uncertainties are predicted thanks to random initialization and data shuffling. Gaussian approximations have also been investigated by several authors, such as the Laplace approximation [START_REF] Ritter | A scalable laplace approximation for neural networks[END_REF][START_REF] Daxberger | Laplace redux-effortless Bayesian deep learning[END_REF] and the highly scalable stochastic weight averaging Gaussian method [START_REF] Maddox | A simple baseline for bayesian uncertainty in deep learning[END_REF].

Regardless of the approximation method, evaluating the quality of the predictive uncertainties remains an open issue. In contrast to previous works, we propose for the first time to estimate sensible coverage probabilities by taking into account the variability induced by the training dataset. We also compare the selected algorithms with the help of discrepancy measures, namely, the maximum mean discrepancy [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF] and the kernelized Stein discrepancy [START_REF] Liu | A kernelized stein discrepancy for goodness-of-fit tests[END_REF]. [START_REF] Izmailov | What are bayesian neural network posteriors really like[END_REF] have studied the performance of Bayesian neural networks (BNNs) by relying on exhaustive full-batch Hamiltonian Monte Carlo. The performance is assessed in terms of RMSE for regression problems and accuracy for classification tasks. Most notably, [START_REF] Izmailov | What are bayesian neural network posteriors really like[END_REF] have applied BNNs to practical deep architectures and large datasets thanks to impressive computational ressources. [START_REF] Wenzel | How good is the Bayes posterior in deep neural networks really?[END_REF] have studied the influence of posterior temperature on the performance of stochastic gradient sampling methods. The performance is reported in terms of accuracy for classification tasks only. More closely related to this work, [START_REF] Yao | Quality of uncertainty quantification for bayesian neural network inference[END_REF] has compared several approximation methods, for regression and classification tasks, and reported the prediction interval coverage probability in order to evaluate the quality of the approximated prediction intervals.

Related works

The objective of this work is to assess the performance of Bayesian neural networks on simple regression tasks by comparing several metrics. In contrast to previous works, we assess the validity of the confidence intervals with marginal and conditional coverage probabilities instead of predictive interval coverage probability, which are defined in section 4. It is worth emphasizing that the aim of this work is not to suggest to use marginal and conditional coverage for practical problems, but to investigate which performance BNNs can achieve. In particular, we investigate if there are correlations between accuracy, coverage metrics and kernel Stein discrepancy, as well as if some algorithms have the same exploration behavior in weight and function spaces.

Bayesian neural networks

Let D = {(X i , Y i )} N i=1 denote a dataset made of N ob- servations of input X i ∈ R D and output Y i ∈ R M pairs.
The observations take the form Y i = f (X i ) + ε i where f : R D → R M represents the unknown latent regression function, and ε i ∼ N (0, σ 2 (X i )I M ) is an additive noise modeling aleatoric uncertainties. The latent function f is approximated by a neural network f (•; w) with parameters w ∈ R d . The observations Y 1 , . . . , Y N are assumed to be i.d.d. and the likelihood function takes the form p(Y i |X i , w) = N (Y i ; f (X i ; w), σ 2 I M ) for any i ∈ {1, . . . , N }. Given a prior distribution p(w) over the network parameters, the posterior distribution p(w|X, Y) can be deduced using Baye's formula: p(w|X, Y) ∝ p(Y|X, w)p(w). The prediction of the model for a new test input x is then given by

p(y|x, X, Y) = R d p(y|x, w)p(w|X, Y)dw.
Unfortunately, this integral is intractable and we must have recourse to approximate inference. The true posterior distribution is approximated by some q(w) ≈ p(w|X, Y) obtained via for instance MCMC methods, variational inference, or Gaussian approximations. In the rest of this section, we briefly describe the approximation methods that we consider in this paper, additional details are given in Appendix A.

Monte Carlo Markov Chain. Monte Carlo Markov chain (MCMC) methods generate a Markov chain whose stationary distribution is the target posterior distribution. In this work, we consider the Hamiltonian Monte Carlo (HMC) algorithm [START_REF] Neal | Mcmc using hamiltonian dynamics[END_REF] which requires the knowledge of the unnormalized posterior p(Y|X, w)p(w) and the gradient of the potential function U (w) = -log(p(Y|X, w)p(w)), given by

∇U (w) = - N i=1 ∇U i (w) -∇ log(p(w)) , where U i (w) = log(p(Y i |X i , w)).
Stochastic gradient Monte Carlo Markov Chain. In order to alleviate the computational cost of classical MCMC methods, many efforts have been dedicated to the development of stochastic gradient MCMC methods [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF][START_REF] Ahn | Bayesian posterior sampling via stochastic gradient fisher scoring[END_REF][START_REF] Ding | Bayesian sampling using stochastic gradient thermostats[END_REF][START_REF] Chen | Stochastic gradient hamiltonian monte carlo[END_REF][START_REF] Ma | A complete recipe for stochastic gradient mcmc[END_REF][START_REF] Li | Preconditioned stochastic gradient langevin dynamics for deep neural networks[END_REF][START_REF] Zhang | Cyclical stochastic gradient mcmc for bayesian deep learning[END_REF]. Here, the Metropolis-Hastings correction step is omitted and the gradient of the potential function is approximated by a stochastic, mini-batched, gradient:

∇U (w) = - N |B| i∈B ∇U i (w) -∇ log(p(w)) .
Despite being computationally more efficient than traditional MCMC, stochastic gradient MCMC methods introduce asymptotic bias. In this work, we consider the stochastic gradient Lagenvin dynamics (SGLD) and Hamiltonian Monte Carlo (SGHMC). We also consider variants of SGM-CMC methods that include variance reduction techniques [START_REF] Baker | Control variates for stochastic gradient mcmc[END_REF][START_REF] Dubey | Variance reduction in stochastic gradient langevin dynamics[END_REF], where the stochastic gradient ∇U is replaced by the following estimation:

∇U (w) = ∇U (w)| w=η + ∇U (w) -∇U (w)| w=η .
We then consider the variants SGLD-CV and SGHMC-CV where η is set to a MAP estimate w MAP , together with the variants SGLD-SVRG and SGHMC-SVRG where η is first set to a MAP estimate, and then updated to the current value of w every m iterations. Finally, we also evaluate the preconditioned method pSGLD of [START_REF] Li | Preconditioned stochastic gradient langevin dynamics for deep neural networks[END_REF], and the cyclical algorithms C-SGLD and C-SGHMC that rely on a cyclical step size [START_REF] Zhang | Cyclical stochastic gradient mcmc for bayesian deep learning[END_REF].

Gaussian approximations. We consider two Gaussian approximations: SWAG and the Laplace approximation with a kronecker factored log likelihood Hessian approximation (LA-KFAC), which both require a pre-trained neural network with parameters w MAP . The SWAG approximation is constructed by collecting values of the parameters along a SGD trajectory with a possibly high step size.

Evaluation metrics

Validity of the confidence intervals. We assess the validity of the confidence intervals produced by an approximation method with coverage probabilities. There are several related notions of coverage probabilities such as the prediction interval coverage probability, marginal coverage probability, and the conditional coverage probability [START_REF] Lin | Locally valid and discriminative prediction intervals for deep learning models[END_REF].

Given a target confidence level 1 -α, the confidence interval ĈD α is said to have a valid prediction interval coverage probability (PICP) if

P{Y ⋆ ∈ ĈD α (X ⋆ )|D} ≥ 1 -α,
where the probability is taken only over a test data (X ⋆ , Y ⋆ ). This coverage probability does not take into account the variability induced by the training dataset D. A more suited property would be the marginal coverage probability (MCP):

P{Y ⋆ ∈ ĈD α (X ⋆ )} ≥ 1 -α,
where the probability is taken over to both the training data D and the test data (X ⋆ , Y ⋆ ). An even stronger property can be obtained by conditioning on the test input data X ⋆ , leading to the conditional coverage probability (CCP):

P{Y ⋆ ∈ ĈD α (X ⋆ )|X ⋆ = x} ≥ 1 -α,
for almost all x ∈ X , where the probability is taken over the training dataset D.

Distance to the HMC reference (weight and function space). We compute the maximum mean discrepancy (MMD) between each approximation and the approximation obtained via exhaustive HMC. Let k : R d × R d → R be a kernel function and let H(k) be the reproducing kernel Hilbert space with kernel k. The MMD between two probability measures Q and Q ′ is defined as [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]:

MMD(Q, Q ′ ) = ∥µ Q -µ Q ′ ∥ H(k) , µ Q = k(•, w) dQ , where µ Q is called the kernel mean embedding of Q in H(k).
Here k is chosen as the distance-based kernel function k(w, w ′ ) = ∥w∥ 2 + ∥w ′ ∥ 2 -∥w -w ′ ∥ 2 proposed by [START_REF] Sejdinovic | Equivalence of distance-based and rkhs-based statistics in hypothesis testing[END_REF]. In this setting, MMD(Q, Q ′ ) = 0 implies that Q = Q ′ . For a given training dataset D and a set of hyperparameters, we compute the MMD between each approximation Q in weight space obtained by SGCM-CMC, SWAG, LA-KFAC, MC-Dropout, or Deep ensembles, and the approximation Q HMC obtained via exhaustive HMC that is considered as the reference. In addition, given a sample {w i } m i=1 from an approximation Q, we can deduce the predictions of the neural network x → f (x; w) for each algorithm and compute the MMD distance to the HMC reference in function space.

Distance to the target posterior (weight space). We also assess the performance of an approximation method by measuring a distance to the target posterior distribution. We rely on the kernelized Stein discrepancy which, in contrast to other distances, only requires the knowledge of the gradient of the log-posterior distribution. The squared KSD between the target posterior measure P and any other probability measure Q is defined as

KSD 2 (P, Q) = E w∼Q E w ′ ∼Q k p (w, w ′ ) ,
where k p denotes the Stein kernel given by

k p (w, w ′ ) = ⟨∇ w , ∇ w ′ k(w, w ′ )⟩ + ⟨s p (w), ∇ w ′ k(w, w ′ )⟩ + ⟨s p (w ′ ), ∇ w k(w, w ′ )⟩ + ⟨s p (w), s p (w ′ )⟩k(w, w ′ ) .
Here, s p denotes the score function of the target posterior distribution, namely, s p (w) = ∇ w log(p(Y|X, w)p(w)). The Stein kernel also depends on an additional kernel function, k, which has to be carefully chosen. Based on the theoretical results of [START_REF] Gorham | Measuring sample quality with kernels[END_REF], the kernel k is chosen as the inverse multi-quadratic (IMQ) kernel, which is defined as k(w, w ′ ) = (1 + ∥w -w ′ ∥ Γ ) -1/2 with Γ = ℓ 2 I. In this setting, the KSD defines a distance between two probability measures such that KSD(P, Q) = 0 implies that P = Q. The lengthscale ℓ of the IMQ kernel is chosen as the median of the pairwise distances, estimated with a subsample of the entire collection of samples generated by all the approximation methods.

Similarities between the algorithms (weight and function space). We also use the MMD in order to establish possible similarities between the algorithms. More precisely, the MMD is computed between each pair ( Q, Q′ ) of approximations obtained with the considered approximation methods and set of hyperparameters. The resulting matrix of pairwise MMD distances is visualized thanks to multidimensional scaling [START_REF] Torgerson | Multidimensional scaling: I. theory and method[END_REF] and analyzed to highlight any similarities between the approximation methods. We also proceed the same in function space.

Experimental setup

The performances of the selected algorithms are studied for several regression tasks. In this section, we summarize the considered datasets and neural network architectures. For each algorithm, we also study the influence of some hyperparameters, which are summarized in the sequel of this section. Additional details about the experiment setup have been reported in the Supplementary Material.

Regression tasks. We consider 4 synthetic regression problems where the output is one-dimensional but the input may be one or multi-dimensional. For each synthetic regression problem, we generate N D = 500 independent datasets D 1 , . . . , D N D used for MCP and CCP computations. We also consider an additional test dataset

D ⋆ = {(X ⋆ i , Y ⋆ i )} N ⋆ i=1
for accuracy, PICP, MCP and CCP. For some problems it may have out-of-distribution (OOD) samples. Each element

(X i , Y i ) of a training dataset D j is such that Y i = f (X i ) + ε i , where ϵ i ∼ N (0, σ 2 ), i = 1, . . . , N .
The underlying latent regression function f and the variance of the noise σ are both known. The four regression problems are constructed with analytical functions, which are summarized in Table 1 below. More details about the regression problems and their generative process can be found in Appendix B. w),w ∼ N (0,I) 0.02 2 120 120 ✓ for every regression problem. The number of parameters ranges from 2, 651 to 20, 501 (see Appendix C for more details). We use a centered normalized Gaussian prior distribution for the weights in all the experiments. Note that for each training data D j , j = 1, . . . , N D , a MAP estimate is computed by training the neural network with an Adam optimizer and an exponentially decaying learning rate. These N D MAP estimates are subsequently used in LA-KFAC, SWAG, SGMCMC-CV, and SGMCMC-SVRG.

Neural network architectures. A feed-forward neural network with ReLU activations is used as a surrogate model

Hyperparameters. Besides the LA-KFAC approximation, all the algorithms depend on one or several hyperparameters. For all the remaining methods (SGMCMC, Deep ensembles, MC-Dropout, and SWAG), we study the influence of the step size ϵ > 0 (i.e., the learning rate) by considering 10 equally log-spaced values ϵ 1 < • • • < ϵ 10 . The ranges are carefully chosen for each type of algorithm and are reported in Appendix A. We also investigate the influence of two additional hyperparameters. Five dropout rates in the MC-Dropout method are considered (0.1, 0.2, . . . , 0.5) and three cycle lengths in the C-SGLD and C-SGHMC methods are investigated (10, 100, and 1000). Other hyperparameters are held fixed, see Appendix A.3 for more details.

HMC reference. For a given training dataset D, a reference sample is generated by Hamiltonian Monte Carlo and subsequently used to evaluate the performance of the selected algorithms (see section 4). We run 3 HMC chains of 200 iterations, discard the first 100 iterations as burn-in, and perform 10, 000 leapfrogs steps. The step size is selected such that the Metropolis-Hastings accept rates are at least above 80%.

Regression without OOD testing

We first present the results for the first experiment (AF#1 in Table 1). Here, the training and test datasets are sampled from the same distribution so that any test input lies within the range of the input training data points.

Coverage probabilities. Figure 1 gathers the graphs of the marginal coverage probabilities and the mean absolute error in conditional coverage with respect to the step size ϵ for a target coverage of 0.95. In addition, the coefficient of determination on the testing set, denoted by

Q 2 = 1 - n i=1 Y ⋆ i -f (X ⋆ i ; w) 2 /Var Y ⋆
, is also reported. It can be observed that (a) The SGLD and SGHMC methods yield similar performances.

(b) The performance of the variants with control variates (SGLD-CV and SGHMC-CV) drops significantly for high step sizes. As pointed out by [START_REF] Nemeth | Stochastic gradient markov chain monte carlo[END_REF], if the state of the Markov chain gets far from the control variate (here a MAP estimate), then the variance of the stochastic gradient can increase instead of being reduced, hence explaining the observed behavior.

(c) Updating the control variates seems to address this issue as illustrated by the performances of SGLD-SVRG, but SGHMC-SVRG is still affected by high step sizes. Given that SGHMC mixes better than SGLD, we believe that this behavior can be corrected by increasing the frequency at which the control variates are updated in SGHMC-SVRG.

(d) The preconditioned SGLD method easily reaches the target coverage probability but has lower coefficients of determination than the other methods.

(e) Deep ensembles and SWAG both yield high coefficient of determinations but only deep ensembles is able to reach the target coverage probability.

(f) MC-dropout has the lowest performances in this example. In particular, the high variability of the absolute error in conditional coverage suggests that the confidence intervals are not smooth with respect to the input space.

(g) Finally, for the LA-KFAC method we obtain Q 2 = 0.99 ± 0.004, MCP = 0.99 ± 0.01, and MAE = 0.048 ± 0.003. As a result, the mean prediction of LA-KFAC is very accurate but the high coverage probability suggests that it overestimates the uncertainties.

We also investigate the marginal coverage probability with respect to various target levels (1 -α), α ∈]0, 1[, which are shown in Figure 2 for 6 selected algorithms (see also Appendix C) . Here, we see that pSGLD and LA-KFAC easily overestimate the target level regardless of the step size. Finally, we report the best marginal coverage probabilities (closest to 0.95) and the best Q 2 coefficients (closest to 1) for each algorithm in Table 2. For each best coverage (resp.

Q 2 ), we also report the associated Q 2 (resp. coverage). Interestingly, we observe that the best coverages lying in 0.95±0.1 do not especially correspond to the best regression Deep ensembles is the only method for which the MMD distance in weight space decreases below 5. In contrast, the MMD distances in function space may decrease towards 0. In particular, the SGLD-SVRG, deep ensembles, and SWAG methods yield the lowest MMD distances in function space. With the LA-KFAC method, we obtained MMD distances of 5.76 ± 0.13 in weight space, and of 1.27 ± 0.21 in function space. Overall and perhaps surprisingly, the behavior with respect to the step size is not the same in weight and function space.

KSD distance to the target posterior. The kernelized Stein discrepancies between the target posterior P and the approximations Q are shown in Figure 4. Here, the lowest discrepancies are obtained with deep ensembles and SWAG, while SGMCMC methods yield higher KSDs. Surprisingly, the LA-KFAC method yields the highest KSDs, with values between 10 5 and 10 7 . If for some algorithms the conclusions from the MMD and the KSD coincide, this is not the case in general: we discuss below some potential explanations related to KSD limitations.

Similarities between the algorithms. The similarities between the algorithms is depicted in Figure 5. At least four groups of approximations can be distinguished in weight space (left panel in Figure 5 Discussion. Based on the results presented in this section, a few observations can be made. First, we see that low KSD values do not especially correspond to good coverage probabilities and Q 2 coefficients. This can be seen in the case of SWAG, which achieves the lowest KSD values but poor coverage probabilities. Amongst the SGMCMC methods, it can also be seen that the lowest KSD values are often obtained for the lowest step sizes. As a result, the KSD seems often unable to identify approximation methods that generate valid confidence intervals. Several studies have shown that the KSD suffers from strong pathologies in sim- ple experiments [START_REF] Wenliang | Blindness of score-based methods to isolated components and mixing proportions[END_REF][START_REF] Korba | Kernel stein discrepancy descent[END_REF]. In particular, the KSD is known to be insensitive to weight proportions in multimodal distributions, which may explain our observations. A similar behavior is observed with the MMD distance to the HMC reference in weight space. For instance, this distance increases with the step size for SGLD but the best predictive performances are obtained for the highest step sizes. In contrast, the MMD distance to the HMC reference in function space is correlated to the predictive performance. A low MMD distance in function space typically corresponds to the best coverage probabilities and Q 2 coefficients obtained by a given approximation method.

Finally, it can be observed that in this experiment, deep ensembles, SGLD-SVRG, and SGHMC-SVRG provide similar approximations as shown by the similarity measures in Figure 5, the MMD distances in function space, and the performances in terms of coverage probabilities and mean predictions. While the LA-KFAC method yields the highest KSD values and higher MMD distances than other methods, it offers a high precision on the test dataset, and a high coverage probability way above the target level.

Regression with OOD testing

In this section, we present the results obtained for the second regression problem (AF#2 in Table 1). Here, the distribution of test dataset differs from the distribution of the training dataset, which is more challenging in terms of coverage probabilities and predictive performances. We find that the behavior of the KSD with respect to the step size and approximation method is similar to the one observed in the previous experiment (see Figure 4). Furthermore, we also observe that the similarities between the approximation methods is close to the one observed in Figure 5. All the associated figures are reported in Appendix C for brevity.

Coverage probabilities. The coverage probabilities metrics and the regression coefficients Q 2 are reported in Figure 6. Compared to the previous experiment of section 6, we see that achieving the target coverage level is more challenging for most of the approximation methods due to the test inputs that are out of the distribution of the training data. Several SGMCMC methods (such as SGLD, SGHMC, and their cyclical variants) easily yield high regression coefficients but struggle to reach the target confidence level. A similar behavior is observed for deep ensembles. This suggests that the width of the confidence intervals are not big enough to properly cover the test predictions. In this experiment, we see that MC-dropout has slightly better performances than in the previous section but remains less effective than the other approximation methods. Finally, we find that the performances of pSGLD and SWAG remain similar to AF#1.

MMD distances to the HMC reference. In terms of MMD distances, we find that deep ensembles and SGLD-SVRG yield the closest approximations to the HMC reference in function space. For brevity, the graphs of the MMD distances have been reported in Appendix C (see Figure 17). Here again, we find that the MMD distance in weight space increases with the step size while it mainly decreases in function space.

Discussion about coverage probabilities

Previous works that assess the quality of the confidence intervals in BNNs rely on the PICP (see, e.g., [START_REF] Yao | Quality of uncertainty quantification for bayesian neural network inference[END_REF]). In contrast to MCP and CCP, PICP only integrates over the test set and does not take into account the variability of the confidence interval with respect to the training set. This is convenient in practice given that we generally do not have access to the distribution of the training data. However, PICP may lead to wrong conclusions about the validity of the confidence intervals. To support our claim, in Figure 7 we compare the three coverage metrics: the PICP histogram obtained by considering all the PICPs (mean over test set) for each training set, while the CPP histogram gathers all the CCPs (mean over training set) for each sample from the test set. We also represent the MCP (mean over both training and test set) and the target confidence level. In both cases, the MCP and CPP metrics clearly indicate that the target coverage is far from being attained. On the contrary, the PICP exhibits large variability over the training sets, with a significant probability of being close to the target coverage for SGLD on the left panel. This implies that, depending on the random sampling of the training set, the PICP may lead to consider SGLD as an efficient approximation of the posterior, while in reality the MCP and CPP are way more pessimistic and tend to conclude the opposite. Consequently, we advocate to use the PICP very cautiously.

Since in practice the MCP and CPP cannot be computed by sampling many independent training sets as we did in our experiments, we suggest instead to estimate them by resampling methods such as the bootstrap, in order to compensate for the PICP limitations illustrated above.

Summary

This work compares several approximation methods for Bayesian neural networks with 14 algorithms including 9 SGMCMC methods, deep ensembles, LA-FKAC, MC-Dropout, SWAG, and HMC which is taken as a reference. In contract to previous works, we assess the performance of the selected methods by evaluating the validity of the generated confidence intervals, the quality of the approximation in weight and function spaces, and the precision of the mean prediction on test data. We find that SGM-CMC methods and deep ensembles are able to achieve good coverage probabilities and predictive performances. Most notably, in our experiments, deep ensembles provide the best approximation to the HMC reference in both weight and function spaces. While SWAG has excellent performances in terms of regression, it tends to underestimate the uncertainties. We also propose a novel similarity analysis between the algorithms which shows that LA-KFAC and pSGLD with high step sizes provide close approximations in weight space. SGMCMC-SVRG and deep ensembles are also very similar in weight space. We finally consider the KSD but unfortunately we observe that it cannot exclusively be relied on to assess the quality of an approximation to the target posterior. In summary, our experiments show that several approximations methods can be efficient, but selecting appropriate hyperparameters yielding valid confidence intervals and regression accuracy is still challenging. A promising way would be to rely on resampling methods on the training set to estimate the MCP and the CCP or use the rencetly introduced conformal prediction [START_REF] Barber | Predictive inference with the jackknife+[END_REF]. We also hope that the proposed new comparisons in weight and function spaces will give insights on how the different algorithms explore the posterior distribution.

We briefly synthesize the considered SGMCMC methods that we consider in our study, with details on software implementation and practical considerations in our numerical experiments.

A.1. Stochastic gradient Markov chain Monte Carlo algorithms

Stochastic gradient Markov chain Monte Carlo (SGMCMC) methods are a family of gradient-based algorithms that rely on the resolution of an Itô stochastic differential equation [START_REF] Ma | A complete recipe for stochastic gradient mcmc[END_REF] (see also [START_REF] Nemeth | Stochastic gradient markov chain monte carlo[END_REF] for a review). Amongst this family of methods, we consider the stochastic gradient Langevin dynamics (SGLD) [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF] and the stochastic gradient Hamiltonian Monte Carlo (SGHMC) [START_REF] Chen | Stochastic gradient hamiltonian monte carlo[END_REF] methods, together with other variants that include variance reduction [START_REF] Dubey | Variance reduction in stochastic gradient langevin dynamics[END_REF][START_REF] Baker | Control variates for stochastic gradient mcmc[END_REF] or preconditioning techniques [START_REF] Li | Preconditioned stochastic gradient langevin dynamics for deep neural networks[END_REF][START_REF] Springenberg | Bayesian optimization with robust bayesian neural networks[END_REF].

Stochastic gradient Langevin dynamics. The SGLD algorithm proposed by [START_REF] Welling | Bayesian learning via stochastic gradient langevin dynamics[END_REF] takes the form

k ≥ 0 , w k+1 = w k -ϵ k ∇U (w k ) + √ 2ϵ k ∆W k+1 ,
where ∆W k+1 is centered and normalized Gaussian random variable. Herein, SGLD is implemented with a constant step size ϵ k = ϵ for all k ≥ 0.

Stochastic gradient Hamiltonian Monte Carlo. The SGHMC algorithm has been proposed by [START_REF] Chen | Stochastic gradient hamiltonian monte carlo[END_REF] where one leapfrog step reads as:

k ≥ 0 : w k+1 = w k + η k M -1 r k , r k+1 = r k -η k ∇U (w k ) -ϵ k CM -1 r k + 2Cη k ∆W k+1 ,
where r denotes the momentum and ∆W k+1 ∼ N (0, I). By introducing the change of variable v k = ϵ k M -1 r k , the above equations can be rewritten as k ≥ 0 :

w k+1 = w k + v k , v k+1 = (1 -α)v k -ϵ k ∇U (w k ) + √ 2αϵ k ∆W k+1 ,
where ϵ k = η 2 k M -1 and α = ηM -1 C. The step size ϵ k is constant in our experiments with SGHMC, and the momentum v is resampled every 10 leapfrog steps following a centered and normalized Gaussian distribution.

Preconditioned SGLD. [START_REF] Ma | A complete recipe for stochastic gradient mcmc[END_REF] proposed the Riemannian SGLD by building upon the Riemannian manifold Langevin and Hamitonian Monte Carlo methods [START_REF] Girolami | Riemann manifold langevin and hamiltonian monte carlo methods[END_REF]. The Riemannian SGLD takes the form:

k ≥ 0 , w k+1 = w k -ϵ k D(w k ) ∇U (w k ) + Γ(w k ) + 2ϵ k D(w k )∆W k+1 , (1)
where D ∈ R d×d can be seen as a preconditioning matrix and the vector Γ ∈ R d is given by

Γ i (w) = j ∂ ∂ς j D ij (w) . ( 2
)
Several forms of the preconditioner D can be used, the most common choice being a constant diagonal matrix D diag of the form D diag = N -1 I d , where N denotes the size of the training dataset. The expected Fisher information matrix is also a natural choice but it is usually intractable, and a few alternatives have been proposed in the literature. Herein we consider the pSGLD method of [START_REF] Li | Preconditioned stochastic gradient langevin dynamics for deep neural networks[END_REF] that neglects the vector Γ and uses a preconditioner inspired by the RMSprop optimization algorithm. The gradient is scaled using a moving average of its norm at each iteration:

D(w k ) = diag λI + L(w k ) -1 , L(w k ) = αL(w k-1 ) + (1 -α) ∇U (w k-1 ) • ∇U (w k-1 ) ,
where α is a parameter with values in [0, 1], λ is a regularization constant, and • denotes the Hadamard (element-wise) product.

SGMCMC-CV and SGMCMC-SVRG. As described in section 3, we consider two variants of the SGLD and SGHMC algorithms that rely on control variates to reduce the variance of stochastic approximation ∇U of ∇U . Following [START_REF] Baker | Control variates for stochastic gradient mcmc[END_REF], SGLD-CV and SGHMC-CV use an estimation of the maximum a posteriori (MAP) for the control variable η.

Hence, the neural network is first trained in a classical fashion using an optimization algorithm, and the resulting values of the network parameters are stored. The remaining variants SGLD-SVRG and SGHMC-SVRG use the strategy proposed by [START_REF] Dubey | Variance reduction in stochastic gradient langevin dynamics[END_REF] which consists in setting η to a MAP estimate as well, but then updating η every m iterations as follows: η = w ℓ if mod(ℓ, m) = 0. The update frequency m has been fixed to 100 in every experiment.

Cyclical SGMCMC. The cyclical variants of SGLD and SGHMC use a step size of the form [START_REF] Zhang | Cyclical stochastic gradient mcmc for bayesian deep learning[END_REF])

ϵ k = ϵ 0 2 cos πmod(k -1, ⌈K/M ⌉) ⌈K/M ⌉ + 1 ,
where ϵ 0 is the initial step size, K denotes total number of iterations, and M is the cycle length. The performance of the cyclical variants is studied with respect to the initial step size and the cycle length, as described in section 4.

A.2. Implementation

Our experiments are implemented with JAX [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] and PyTorch (Paszke et al., 2019). For the LA-FKAC approximation, we use the Laplace package that implements various types of Laplace approximations [START_REF] Daxberger | Laplace redux-effortless Bayesian deep learning[END_REF] and relies on PyTorch. The remaining methods are implemented with Google's JAX framework. Regarding the SGMCMC methods, we rely on BlackJAX [START_REF] Lao | A sampling library for JAX[END_REF], where SGLD and SGHMC are already implemented. We use our own implementations of variants with variance reduction, integrated within BlackJAX. The remaining approximation methods are implemented from scratch with JAX.

A.3. Hyperparameters for numerical experiments

The selected step sizes for each SGMCMC algorithm are given in Table 3. Any other hyperparameter, such as the number of iterations or batch size, is held fixed. For all SGMCMC methods, the number of iterations is set to 10 5 and the first 50 4 iterations are discarded as burn-in. The remaining iterations are automatically thinned by selecting 2000 particles which are obtained by minimizing the maximum mean discrepancy [START_REF] Gretton | A kernel method for the two-sample-problem[END_REF]. This procedure, referred to as MMD thinning, is applied to every SGMCMC output and works as follows. Given T samples w 1 , . . . , w T of a network parameters w ∈ R d , we find a sequence of indices π ∈ {1, . . . , T } m such that w π(1) , . . . , w π(m) represent the selected samples. The sequence of indices π is obtained thanks to a greedy quantization algorithm [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF] that minimizes the maximum mean discrepancy. Let then P T be the empirical distribution of the T samples w 1 , . . . , w T , P T = T i=1 δ(w i ). At the (i + 1)-th iteration of the quantization algorithm, the index π(i + 1) is obtained by minimizing the MMD as follows:

π(i + 1) = arg min j∈{1,...,T } MMD 2 (P T , Q i+1 (j)) ,
where Q i+1 (j) denotes the empirical measure of the already i selected samples w 1 , . . . , w i , and an additional sample w j to be determined:

Q i+1 (j) = 1 i + 1 i ℓ=1 δ(w π(ℓ) ) + 1 i + 1 δ(w j ) .
The underlying kernel function k is chosen as the following characteristic distance-based kernel [START_REF] Sejdinovic | Equivalence of distance-based and rkhs-based statistics in hypothesis testing[END_REF])

k(w, w ′ ) = ∥w∥ 2 + ∥w ′ ∥ 2 -∥w -w ′ ∥ 2 .
In our all experiments, we select m = 2000 samples in order to represent each SGMCMC output. In this case of the variance reduction technique SVRG, the control variates are updated every 100 iterations. We use 10 leapfrog steps in SGHMC and each of its variants. For deep ensembles, 200 neural networks are independently trained and their weights are initialized from a centered normalized Gaussian distribution. For each of the remaining approximation methods (LA-KFAC, MC-Dropout, and SWAG), samples of sizes 2000 are saved.

B. Description of the training and test datasets

The synthetic regression problems are described in more details below.

Regression problem AF#1. This first regression test case is a homoscedastic regression of a one-dimensional function.

The dataset -3, 3]), and Y i = cos(2X i ) + sin(X i ) + ϵ i , where ϵ i ∼ N (0, σ), with σ = 0.2. The test dataset D ⋆ has the same distribution.

D = {X i , Y i } N i=1 is such that X i ∼ U([
Regression problem AF#2. This test case is taken from [START_REF] Yao | Quality of uncertainty quantification for bayesian neural network inference[END_REF] which consists in a homoscedastic regression of a one-dimensional function as well. The inputs X 1 , . . . , X N of a training dataset D are uniformly sampled from [-4, -1]∪[1, 4], while test inputs are uniformly distributed between [-4, 4]. The observations are defined as Y i = 0.1X 3 i +ϵ i where ϵ i ∼ N (0, 0.25).

Regression problem AF#3. This test case is taken from [START_REF] Yao | Quality of uncertainty quantification for bayesian neural network inference[END_REF] which consists in a homoscedastic regression of a one-dimensional function. The first 80 inputs X 1 , . . . , X 80 of a training dataset D are uniformly samples in [-6, -2] ∪ [2, 6], and the last two samples X 81 , X 82 are uniformly sampled in [-2, 2]. In contrast, the inputs of a test dataset D ⋆ are all uniformly distributed in [-6, 6]. The output observations are given by Y i = -(1 + X i ) sin(1.2X i ) + ϵ i where ϵ ∼ N (0, 0.25).

Regression problem AF#4. This last regression problem is taken from the work of [START_REF] Izmailov | What are bayesian neural network posteriors really like[END_REF]. The inputs of a training dataset 14,18]. The observations Y i are defined as Y i = f (X i ; w 0 ) + ϵ i , with ϵ i ∼ N (0, 0.02), and where f (•; w 0 ) denotes a feed-forward neural network with three hidden layers of sizes 100 (and thus 20, 501 parameters in total). The weight parameters w 0 are sampled once from a centered normalized Gaussian distribution N (0, I d ). The inputs of the test dataset are uniformly distributed in [-12, 22].

D = {(X i , Y i )} N i=1 are uniformly distributed in [-10, -6] ∪ [6, 10] ∪ [

C. Additional results

In this section, we gather additional results about the empirical experiments, and a few additional details about the networks architecture. For the first synthetic problem (AF #1), the architecture is made of two hidden layers with 100 hidden features, leading to a network with d = 10, 401 parameters. Following [START_REF] Yao | Quality of uncertainty quantification for bayesian neural network inference[END_REF], we consider two hidden layers with 50 features for the second and third regression problems (AF #2 and AF #3), such that the network is made of d = 2, 651 parameters. Finally, following [START_REF] Izmailov | What are bayesian neural network posteriors really like[END_REF], a network made of three hidden layers with 100 features each is use for the last synthetic problem (AF #4), leading to a network made of d = 20, 501 parameters.

C.1. Best marginal coverage probabilities and regression coefficients

We first gather below in Tables 4-6 the best marginal coverages and best Q 2 coefficients obtained for the experiments AF#2, AF#3, and AF#4.

C.2. Similarities between the algorithms

The similarities between the algorithms are gathered for all the experiments in Figure 8. The similarities in weight space seem to exhibit the same structure in our four experiments. However, the similarities in function space do not exhibit any particular structure.

C.3. Additional results for regression problem AF#1

Results obtained for the cyclical variants of SGMCMC have been reported herein. Figure 9 shows the Q 2 regression cofficient and coverage probabilities obtained with CSGMCMC. We also report here graphs of the marginal coverage probability with respect to the target confidence level (1 -α), for some α ∈ [0.05, 0.95], and all the considered values for the underlying hyperparameters, which are gathered in Figures 10111213.

We find that most SGMCMC methods (without variance reduction) give rise to marginal coverage probabilities that do not increase monotonically with the target level, and that are easily much higher than the target level. The cyclical step size reduces the marginal coverage probabilities. The pSGLD methods yield very different results as the target level is easily achieved, even for the lowest step sizes. LA-FKAC shows a similar behavior, where the marginal coverage probability is always above the target level.

The MC-Dropout method is able to reach the target level whenever α is small enough, but struggles to reach an appropriate marginal coverage as α increases. Finally, we find that the marginal coverage obtained with deep ensembles varies monotonically with the target level for sufficiently high step sizes.

C.4. Additional results for regression problem AF#2

The coverage metrics and regression coefficient obtained with CSGLD and CSGHMC are shown in Figure 13. The performances of the cyclical variants are similar to those of SGLD and SGHMC (see Figure 6). The graphs of the kernelized Stein discrepancies and maximum mean discrepancies are shown in Figures 14 and 17

Figure 1 .

 1 Figure 1. Problem AF#1. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95. Results obtained for the cyclical SGMCMC variants are reported in the Appendix.

Figure 2 .

 2 Figure 2. Problem AF#1. Graphs of the marginal coverage probability with respect to the target level and for 10 step sizes.

  ): SGLD, SGLD-CV and CS-GLD methods, SGHMC, SGHMC-CV and CSGHMC, pSGLD and LA-KFAC, and the remaining methods form the last group. Amongst the SGLD and SGHMC methods, it is clearly observed that CV variants yield distinct approximations. Surprisingly, SGLD-SVRG and SGHMC-SVRG generate similar approximations to deep ensembles, in weight space. The similarities in function space do not exihibit any particular structure and have been reported in Appendix C.2 (see Figure C.2).

Figure 3 .

 3 Figure 3. Problem AF#1. Maximum mean discrepancies MMD(Q, QHMC) between the approximated posterior distributions Q and the reference HMC sample, in weight and function spaces.

Figure 4 .

 4 Figure 4. Problem AF#1. Kernelized Stein discrepancies KSD(P, Q) between the target posterior measure P and the approximated posteriors Q. KSD values for the cyclical SGMCMC methods are reported in the Appendix.

Figure 5 .

 5 Figure 5. Dataset AF#1. Similarities between the algorithms as measured by the MMD in weight space and represented in a twodimensional space build with multidimensional scaling. The markers sizes are proportional to the value of the step size ϵ.

Figure 6 .

 6 Figure 6. Problem AF#2. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95.

Figure 7 .

 7 Figure 7. Problem AF#1. Comparison of the three coverage probabilities obtained with SGLD and SGLD-SVRG.

Figure 8 .

 8 Figure 8. Similarities between the algorithms as measured by the MMD and represented in a two-dimensional space build with multidimensional scaling. Left panels: weight space, right panels: function space. The markers sizes are proportional to the value of the underlying step size ϵ.

  -18.C.5. Additional results for regression problems AF#3 and AF#4All the results obtained for the third and fourth regression problems are shown in Figures[19][20][21][22][23][24]. The distributions of the test data differ significantly from the training data in problem AF#3, and even more in problem AF#4. It can be seen in the associated results that these problems are challenging for all the methods, which are not able to reach the target coverage probability, nor obtain good regression coeffcients.

Figure 9 .

 9 Figure 9. Problem AF#1. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95.

Figure 10 .

 10 Figure 10. Problem AF#1. Graphs of the marginal coverage probability with respect to the target level for SGMCMC methods. The curves are colored by the value of the underlying step size: dark blue corresponds to the lowest, while dark red corresponds to the highest step size.

Figure 11 .

 11 Figure 11. Problem AF#1. Graphs of the marginal coverage probability with respect to the target level for MC-Dropout, SWAG, LA-KFAC, deep ensembles, and pSGLD.

Figure 12 .

 12 Figure 12. Problem AF#1. Kernelized Stein discrepancies KSD(P, Q) between the target posterior measure P and the approximated posteriors Q.

Figure 13 .

 13 Figure 13. Problem AF#2. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95.

Figure 14 .

 14 Figure 14. Problem AF#2. Kernelized Stein discrepancies KSD(P, Q) between the target posterior measure P and the approximated posteriors Q.

Figure 15 .

 15 Figure 15. Problem AF#2. Graphs of the marginal coverage probability with respect to the target level for SGMCMC methods. The curves are colored by the value of the underlying step size: dark blue corresponds to the lowest, while dark red corresponds to the highest step size.

Figure 16 .

 16 Figure 16. Problem AF#2. Graphs of the marginal coverage probability with respect to the target level for MC-Dropout, SWAG, LA-KFAC, deep ensembles, and pSGLD.

Figure 17 .

 17 Figure 17. Problem AF#2. Maximum mean discrepancies MMD(Q, QHMC) between the approximated posterior distributions Q and the reference HMC sample, in weight and function spaces.

Figure 18 .

 18 Figure 18. Problem AF#2. Maximum mean discrepancies between the approximated posterior distributions and the reference HMC sample, in weight and function spaces.

Figure 19 .

 19 Figure 19. Problem AF#3. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95.

Figure 20 .

 20 Figure 20. Problem AF#3. Maximum mean discrepancies MMD(Q, QHMC) between the approximated posterior distributions Q and the reference HMC sample, in weight and function spaces.

Figure 22 .

 22 Figure 22. Problem AF#3. Kernelized Stein discrepancies KSD(P, Q) between the target posterior measure P and the approximated posteriors Q.

Figure 23 .

 23 Figure 23. Problem AF#3. Graphs of the marginal coverage probability with respect to the target level for SGMCMC methods. The curves are colored by the value of the underlying step size: dark blue corresponds to the lowest, while dark red corresponds to the highest step size.

Figure 24 .

 24 Figure 24. Problem AF#3. Graphs of the marginal coverage probability with respect to the target level for MC-Dropout, SWAG, LA-KFAC, deep ensembles, and pSGLD.

Figure 25 .

 25 Figure 25. Problem AF#4. Coverage metrics and Q 2 coefficient with respect to the step size ϵ. The target coverage is set to 0.95.

Figure 27 .

 27 Figure 27. Problem AF#4. Maximum mean discrepancies between the approximated posterior distributions and the reference HMC sample, in weight and function spaces.

Figure 28 .

 28 Figure 28. Problem AF#4. Kernelized Stein discrepancies KSD(P, Q) between the target posterior measure P and the approximated posteriors Q.

Figure 29 .

 29 Figure 29. Problem AF#4. Graphs of the marginal coverage probability with respect to the target level for SGMCMC methods. The curves are colored by the value of the underlying step size: dark blue corresponds to the lowest, while dark red corresponds to the highest step size.

Figure 30 .

 30 Figure 30. Problem AF#4. Graphs of the marginal coverage probability with respect to the target level for MC-Dropout, SWAG, LA-KFAC, deep ensembles, and pSGLD.

  The SWAG approximation reads as p(w|X, Y) ≈ N (w|w SWAG , Σ SWAG ) where w SWAG is given by the running mean of the SGD iterates, and Σ SWAG is a diagonal plus low rank approximation. The LA-KFAC approximation is given by p(w|X, Y) ≈ N (w|w MAP , Σ KFAC ) where Σ KFAC denotes the Kronecker-factored approximate curvature approximation (see, e.g.,[START_REF] Daxberger | Laplace redux-effortless Bayesian deep learning[END_REF]).

	Variational methods. The Monte Carlo dropout method
	(MC-Dropout) proposed by Gal & Ghahramani (2016) is
	considered herein. The neural network is augmented with
	dropout layers which are activated during both training and
	inference stages. The dropout mechanism is here applied to
	the output features of each layer.
	Deep ensembles. We also investigate the deep ensem-
	bles method of Lakshminarayanan et al. (2017); Fort et al.
	(2019), which consists in training several neural networks in-
	dependently with random initializations, and gathering their
	predictions to obtain a mean prediction and uncertainties.

Table 1 .

 1 Description of the synthetic datasets where AF stands for analytical functions.

	Task	Latent function	σ D N N ⋆ OOD
	AF#1	cos(2x) + sin(x)	0.2 1 100 200 ✗
	AF#2	0.1x 2	0.25 1 100 200 ✓
	AF

Table 2 .

 2 Best marginal coverage probabilities (MCP) and best Q 2 coefficients obtained with each algorithm.

	2 MCP

Method

Best MCP

Q 2 Best Q 2 MCP Method Best MCP Q 2 Best Q

Table 3 .

 3 Selected step sizes for each algorithm, equally log-spaced.

	Step size SGLD variants SGHMC variants pSGLD, Ensemble, MC Dropout LA-KFAC	SWAG
	ϵ 1	1e-08	1e-08	0.0001	0.000001 1.000000e-07
	ϵ 2	2.154435e-08	1.668101e-08	0.000215	0.000004 1.668101e-07
	ϵ 3	4.641589e-08	2.782559e-08	0.000464	0.000013 2.782559e-07
	ϵ 4	1.000000e-07	4.641589e-08	0.001000	0.000046 4.641589e-07
	ϵ 5	2.154435e-07	7.742637e-08	0.002154	0.000167 7.742637e-07
	ϵ 6	4.641589e-07	1.291550e-07	0.004642	0.000599 1.291550e-06
	ϵ 7	1.000000e-06	2.154435e-07	0.010000	0.002154 2.154435e-06
	ϵ 8	2.154435e-06	3.593814e-07	0.021544	0.007743 3.593814e-06
	ϵ 9	4.641589e-06	5.994843e-07	0.046416	0.027826 5.994843e-06
	ϵ 10	1e-05	1e-06	0.1	0.1	1e-05

Table 4 .

 4 AF#2: Best marginal coverage probabilities (MCP) and best Q 2 coefficients obtained with each algorithm.

	Method	Best MCP Q 2 Best Q 2 MCP Method	Best MCP Q 2 Best Q 2 MCP
	SGLD	0.94	0.98 0.989 0.696 CSGHMC(10)	0.947 0.988 0.988 0.947
	SGLD-CV	0.809	-4e8 0.954 0.497 CSGHMC(100)	0.946 0.990 0.99 0.946
	SGLD-SVRG	0.954 0.995 0.997 0.908 CSGHMC(1000) 0.928 0.99 0.99 0.928
	SGHMC	0.96	0.99 0.99 0.960 pSGLD	0.962 0.961 0.993 0.996
	SGHMC-CV	0.825 -13.05 0.944 0.583 Deep ensemble	0.964 0.997 0.997 0.736
	SGHMC-SVRG 0.942 0.992 0.995 0.867 SWAG	0.273 0.997 0.997 0.273
	CSGLD(10)	0.949 0.978 0.989 0.655 MC Drop.(0.1)	0.890 0.994 0.994 0.89
	CSGLD(100)	0.92	0.976 0.99 0.652 MC Drop.(0.2)	0.917 0.989 0.989 0.917
	CSGLD(100)	0.91	0.980 0.989 0.648 MC Drop.(0.3)	0.959 0.979 0.984 0.909

Table 5 .

 5 AF#3: Best marginal coverage probabilities (MCP) and best Q 2 coefficients obtained with each algorithm.

	Method	Best MCP Q 2 Best Q 2 MCP Method	Best MCP Q 2 Best Q 2 MCP
	SGLD	0.921	0.803 0.823 0.885 CSGHMC(10)	0.834 0.825 0.867 0.781
	SGLD-CV	0.551 -29.357 -0.091 0.034 CSGHMC(100)	0.803 0.794 0.850 0.753
	SGLD-SVRG	0.919	0.859 0.879 0.877 CSGHMC(1000) 0.805 0.817 0.851 0.740
	SGHMC	0.837	0.788 0.847 0.801 pSGLD	0.919 0.851 0.879 0.907
	SGHMC-CV	0.387 -28.251 -0.165 0.131 Deep ensemble	0.94 0.933 0.944 0.821
	SGHMC-SVRG 0.736	-1e9 0.655 0.614 SWAG	0.751 0.538 nan 0.654
	CSGLD(10)	0.902	0.829 0.833 0.856 MC Drop.(0.1)	0.624 0.792 0.799 0.588
	CSGLD(100)	0.893	0.739 0.803 0.840 MC Drop.(0.2)	0.554 0.662 0.662 0.554
	CSGLD(1000)	0.900	0.818 0.859 0.851 MC Drop.(0.3)	0.415 0.436 0.437 0.397

Table 6 .

 6 AF#4: Best marginal coverage probabilities (MCP) and best Q 2 coefficients obtained with each algorithm.

	Method	Best MCP	Q 2	Best Q 2 MCP Method	Best MCP	Q 2	Best Q 2 MCP
	SGLD	0.646	-2.136 0.785 0.188 CSGHMC(10)	0.589	-2.316	0.814 0.133
	SGLD-CV	0.696	-57.602 0.781 0.130 CSGHMC(100)	0.836	0.446	0.780 0.123
	SGLD-SVRG	0.529	0.488	0.869 0.390 CSGHMC(1000) 0.475	0.189	0.830 0.266
	SGHMC	0.661	-1.575 0.792 0.098 pSGLD	0.957 -20858.904 -324.041 1.0
	SGHMC-CV	0.616 -1284.229 0.786 0.148 Deep ensemble	0.413	0.420	0.925 0.203
	SGHMC-SVRG 0.607	-0.367 0.829 0.176 swag	0.537	0.862	0.862 0.537
	CSGLD(10)	0.639	-40.265 0.769 0.087 MC Drop.(0.1)	0.232	0.802	0.851 0.227
	CSGLD(100)	0.654 -5.929760 0.785 0.139 MC Drop.(0.2)	0.400	0.796	0.860 0.270
	CSGLD(1000)	0.563	-0.639 0.769 0.126 MC Drop.(0.3)	0.453	0.549	0.821 0.314

Supplementary Material

A. Considered algorithms