Renye Cai

Frédéric Holweck

Zhi-Qiang Feng

François Peyraut
email: francois.peyraut@utbm.fr

Integrity basis of polyconvex invariants for modeling hyperelastic orthotropic materials -Application to the mechanical response of passive ventricular myocardium

The present paper proposes a new Strain Energy Function (SEF) for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. In order to build our SEF, we have followed a classical strategy based on exponential functions, but we have chosen to work with polyconvex invariants instead of the standard ones. Actually, in the context of hyperelastic problems, the polyconvexity of the strain energy density is considered as a prerequisite for ensuring the existence of solutions. By selecting a set of polyconvex invariants, we demonstrate that our model can predict the experimental data with 6 different shear modes applied to passive ventricular myocardium.

Introduction

Understanding the behavior of anisotropic hyperelastic materials is of major importance for scientists because their modeling has a wide range of applications in engineering biosciences such as in health therapeutic, medical prosthesis, ergonomics or virtual surgery. The mechanical study of the shear deformation of myocardial layers is for example useful because these deformations are considered to play an important role in the mechanical behavior of the heart [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF].

These past ten years, many works have been performed to investigate the structurally based model originally proposed by Holzapfel and Ogden in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] in relation to the tests carried out in [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF] where the orthotropic nature of the ventricular myocardium has been proven. Erikson et al. [START_REF] Eriksson | Modeling the dispersion in electromechanically coupled myocardium[END_REF] have for example improved this model with structure parameters allowing the quantification of the degree of dispersion based on measured fiber and sheet angle data. Melnik et al. [START_REF] Melnik | A generalised structure tensor model for the mixed invariant 𝐼 8[END_REF] have modified the hyperelastic Holzapfel-Ogden model by using Generalized Structure Tensors (GSTs) in order to account for fiber dispersion in every term of the strain-energy function. McEvoy et al. [START_REF] Mcevoy | Compressibility and anisotropy of the ventricular myocardium: Experimental analysis and microstructural modeling[END_REF] have studied the effect of compressibility by combining a nonlinear volumetric hyperelastic component with the Holzapfel-Ogden anisotropic hyperelastic component for myocardium fibers. Palit et al. [START_REF] Palit | In vivo estimation of passive biomechanical properties of human myocardium[END_REF] have implemented the Holzapfel-Ogden constitutive law in the MSCMarc finite element software in order to inversely estimate the constitutive parameters of the model. Very recently, Li et al. [START_REF] Li | Insights into the passive mechanical behavior of left ventricular myocardium using a robust constitutive model based on full 3D kinematics[END_REF] have extended the Holzapfel-Ogden model by accounting for the mixed-invariants in the fiber-normal and sheet-normal directions.

All these models, like in most of the papers published in the literature, separate the SEF into an isotropic part and an anisotropic part. The first part is used to model the low strain behavior of the ground matrix and the second part accounts for the behavior of the fibers at higher strain. A first alternative to this standard approach is to introduce free-invariants SEF [START_REF] O'shea | Hyperelastic constitutive modelling for transversely isotropic composites and orthotropic biological tissues[END_REF]. Another alternative consists in using an integrity basis of invariants as proposed in [START_REF] Ta | A new invariant-based method for building biomechanical behavior laws-application to an anisotropic hyperelastic material with two fiber families[END_REF] and mixing them in a single SEF [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]. This second alternative, mathematically justified by the theory of invariant polynomials and by Noether's theorem, is inspired by the pioneering work of Thionnet et al. [START_REF] Thionnet | A new constructive method using the theory of invariants to obtain material behavior laws[END_REF]. One advantage of this second alternative is to work with polyconvex invariants, the polyconvexity of the strain energy density being considered as a prerequisite for ensuring the existence of solutions in compatibility with physical requirements [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. A wide survey with many proofs on polyconvexity of isotropic and transversely isotropic functions can be viewed in [START_REF] Schröder | Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions[END_REF].

The present paper proposes a new SEF for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. In order to build our SEF, we have followed a strategy based on exponential functions, as proposed by Holzapfel et al. [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF], but we have selected polyconvex invariants instead of the standard ones generally used in the literature [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF][START_REF] Eriksson | Modeling the dispersion in electromechanically coupled myocardium[END_REF][START_REF] Melnik | A generalised structure tensor model for the mixed invariant 𝐼 8[END_REF][START_REF] Mcevoy | Compressibility and anisotropy of the ventricular myocardium: Experimental analysis and microstructural modeling[END_REF][START_REF] Palit | In vivo estimation of passive biomechanical properties of human myocardium[END_REF][START_REF] Li | Insights into the passive mechanical behavior of left ventricular myocardium using a robust constitutive model based on full 3D kinematics[END_REF]. Working with the set of polyconvex invariants exhibited in [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF] allows to replace the classical mixed invariant 𝐼 8 , which is non polyconvex (the proof is in Section 4), by the polyconvex invariant 𝐿 4 defined by Eq. [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]. This replacement provides consistent numerical results with [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] and our model can perfectly match the experimental data obtained by Dokos et al. [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF] with 6 different shear modes applied to passive ventricular myocardium. In addition, sufficient conditions ensuring the convexity of the SEF are given in Section 7. Finally, our model can represent stretch effect in the 3 specific directions of passive ventricular myocardium, the myocyte axis direction, the direction lying within the muscle layer and transverse to the myocyte axis, and the direction normal to the muscle layer.

Notations

A bold-face Latin lowercase letter, say a, and a bold-face Latin capital letter, say A, will denote a vector and second-order tensor, respectively. The standard Euclidean inner product ⟨., .⟩ in a n vector space dimension, its related norm ‖.‖, and the product .⊗. between two vectors a and b, are respectively defined by:

⟨𝑨𝒂, 𝒂⟩ = 𝑛 ∑ 𝑖=1 𝐴 𝑖𝑗 𝑎 𝑗 𝑎 𝑖 ; ‖𝒂‖ = √ √ √ √ 𝑛 ∑ 𝑖=1 ( 𝑎 𝑖 ) 2 ; (𝒂 ⊗ 𝒃) 𝑖𝑗 = 𝑎 𝑖 𝑏 𝑗

Kinematics and basic continuum mechanics

Consider a continuum body V with particle P ∈ V which is embedded in the three-dimensional Euclidean space at a given instant t of time. As the continuum body V moves in space from one instant of time to another it occupies a continuous sequence of geometrical regions denoted by 𝛺 0 … 𝛺 𝑖 . 𝛺 0 is referred to as the fixed reference (undeformed or Lagrangian) configuration (𝑡 = 0) of the body V while the configuration t is called the current (deformed or Eulerian) configuration. The position 𝒙 of particle P in the current (or deformed) configuration can be deduced by the motion 𝜑 that depends on its initial position vector 𝑿 and the time 𝑡:

𝒙 = 𝜑(𝑿, 𝑡) (1)
The deformation gradient matrix 𝑭 is defined by:

𝑭 = 𝜕𝒙 𝜕𝑿 = 𝑰 + 𝜕𝑼 (𝑿, 𝑡) 𝜕𝑿 (2)
𝑰 is the unity tensor and 𝑼 (𝑿, 𝑡) = 𝒙 -𝑿 is the displacement vector of the particle. The tensors 𝑪 and B are the so-called right and left Cauchy-Green strain tensors:

𝑪 = 𝑭 𝑇 𝑭 ; 𝑩 = 𝑭 𝑭 𝑇 (3)
In terms of stress tensor, we remind that the second Piola-Kirchhoff stress tensor 𝑺 and the corresponding Cauchy stress tensor 𝝈 are obtained by differentiating a SEF W with respect to C:

𝑺 = 2 𝜕𝑊 𝜕𝑪 -𝑝𝑪 -1 ; 𝝈 = 𝐽 -1 𝑭 𝑺𝑭 𝑇 (4)
where the extra pressure 𝑝 is included in the formulation to account for the incompressibility condition 𝐽 = det (𝑭 ) = 1. Combining the two equalities of Eq. ( 4) yields to:

𝝈 = 2𝐽 -1 𝑭 𝜕𝑊 𝜕𝑪 𝑭 𝑇 -𝑝𝑰 (5)
Note that, in the following, the time 𝑡 will be omitted because we are only interested in static problems.

Integrity basis of polyconvex invariants for modeling orthotropic materials

The mechanical response of passive ventricular myocardium to simple shear loading is sensitive to the shear mode as described in [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF]. In this context, the fiber direction 𝒇 𝟎 , the direction 𝒔 𝟎 transverse to the fiber direction in plane of sheets, and the normal 𝒏 𝟎 to sheets play a key role to define the material symmetry group (Fig. 1). In particular, the 3 planes perpendicular to the three directions 𝒇 𝟎 , 𝒔 𝟎 and 𝒏 𝟎 are planes of symmetry for the material. The material properties remain therefore invariant under the action of [START_REF] Eriksson | Modeling the dispersion in electromechanically coupled myocardium[END_REF] 6)-( 7) is stable under multiplication, contains the neutral element I and each element has an inverse in 𝑆 8 . This confers a group structure to 𝑆 8 . Using a mathematical argument based on the Reynolds operator and on the Noether's theorem, Ta et al. [START_REF] Ta | A new invariant-based method for building biomechanical behavior laws-application to an anisotropic hyperelastic material with two fiber families[END_REF] have presented a systematic method to find a set of invariants associated with 𝑆 8 and demonstrated that the following 7 polynomial invariants 𝐾 𝑖 form an integrity basis of the ring of invariant polynomials under the material symmetry group 𝑆 8 :

𝐾 1 = 𝜌 1 ; 𝐾 2 = 𝜌 2 ; 𝐾 3 = 𝜌 3 ; 𝐾 4 = 𝜌 2 4 ; 𝐾 5 = 𝜌 2 5 ; 𝐾 6 = 𝜌 2 6 ; 𝐾 7 = 𝜌 4 𝜌 5 𝜌 6 (8) { 𝜌 1 = ⟨𝑪𝒇 𝟎 , 𝒇 𝟎 ⟩ ; 𝜌 2 = ⟨𝑪𝒔 𝟎 , 𝒔 𝟎 ⟩ ; 𝜌 3 = ⟨𝑪𝒏 𝟎 , 𝒏 𝟎 ⟩ 𝜌 4 = ⟨𝑪𝒇 𝟎 , 𝒔 𝟎 ⟩ ; 𝜌 5 = ⟨𝑪𝒇 𝟎 , 𝒏 𝟎 ⟩ ; 𝜌 6 = ⟨𝑪𝒔 𝟎 , 𝒏 𝟎 ⟩ (9)
Later, Cai et al. [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF] combined the invariants 𝐾 𝑖 in order to build an integrity basis of 7 polyconvex invariants 𝐿 𝑖 :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝐿 1 = 𝐾 1 ; 𝐿 2 = 𝐾 2 ; 𝐿 3 = 𝐾 3 ; 𝐿 4 = ( 𝐾 1 + 𝐾 2 ) 2 + 4𝐾 4 ; 𝐿 5 = ( 𝐾 1 + 𝐾 3 ) 2 + 4𝐾 5 𝐿 6 = ( 𝐾 2 + 𝐾 3 ) 2 + 4𝐾 6 ; 𝐿 7 = 𝐾 7 + 1 2 ( 𝐾 1 𝐾 2 𝐾 3 -𝐾 1 𝐾 6 -𝐾 2 𝐾 5 -𝐾 3 𝐾 4 ) (10)
It is also demonstrated in [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF] that 𝐿 7 contributes in Eq. ( 5) under the form of a redundant term with the extra pressure p. 𝐿 7 can therefore be removed from the set of invariants. Consequently, in order to benefit from the property of polyconvexity, we assume that the SEF W depends on the six first polyconvex invariants 𝐿 𝑖 . Therefore, Eq. ( 5) yields to:

𝝈 = 2𝐽 -1 𝑭 ( 6 ∑ 𝑖=1 𝜕𝑊 𝜕𝐿 𝑖 𝜕𝐿 𝑖 𝜕𝑪 ) 𝑭 𝑇 -𝑝𝑰 (11)
The derivatives

𝜕𝐿 𝑖
𝜕𝑪 embedded in Eq. ( 11) are calculated straightforwardly from Eqs. ( 8)- [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]:

𝜕𝐿 1 𝜕𝑪 = 𝒇 𝟎 ⊗ 𝒇 𝟎 ; 𝜕𝐿 2 𝜕𝑪 = 𝒔 𝟎 ⊗ 𝒔 𝟎 ; 𝜕𝐿 3 𝜕𝑪 = 𝒏 𝟎 ⊗ 𝒏 𝟎 ( 12
)
𝜕𝐿 4 𝜕𝑪 = 2 { ( 𝐿 1 + 𝐿 2 ) ( 𝒇 𝟎 ⊗ 𝒇 𝟎 + 𝒔 𝟎 ⊗ 𝒔 𝟎 ) + √ 𝐿 4 - ( 𝐿 1 + 𝐿 2 ) 2 ( 𝒇 𝟎 ⊗ 𝒔 𝟎 + 𝒔 𝟎 ⊗ 𝒇 𝟎 ) } ( 13
)
𝜕𝐿 5 𝜕𝑪 = 2 { ( 𝐿 1 + 𝐿 3 ) ( 𝒇 𝟎 ⊗ 𝒇 𝟎 + 𝒏 𝟎 ⊗ 𝒏 𝟎 ) + √ 𝐿 5 - ( 𝐿 1 + 𝐿 3 ) 2 ( 𝒇 𝟎 ⊗ 𝒏 𝟎 + 𝒏 𝟎 ⊗ 𝒇 𝟎 ) } ( 14
)
𝜕𝐿 6 𝜕𝑪 = 2 { ( 𝐿 2 + 𝐿 3 ) ( 𝒔 𝟎 ⊗ 𝒔 𝟎 + 𝒏 𝟎 ⊗ 𝒏 𝟎 ) + √ 𝐿 6 - ( 𝐿 2 + 𝐿 3 ) 2 ( 𝒔 𝟎 ⊗ 𝒏 𝟎 + 𝒏 𝟎 ⊗ 𝒔 𝟎 ) } (15)
In order to calculate the derivatives 𝜕𝑊 𝜕𝐿 𝑖 also included in Eq. ( 11), it is necessary to define how the SEF W depends on the invariants 𝐿 i . This issue is discussed in the next section.

Strain energy function

Dokos et al

. [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF] have demonstrated that passive ventricular myocardium behaves differently depending on the plane of shear. The 6 different modes of shear are shown in Fig. 1 where the myocyte axis representing the fiber direction is depicted with green full lines. Building a consistent model thus requires accounting for the myocyte axis direction 𝒇 𝟎 , the direction 𝒔 𝟎 lying within the muscle layer and transverse to 𝒇 𝟎 , and the direction 𝒏 𝟎 normal to the muscle layer.

Due to the exponential trends of the stress observed experimentally, Holzapfel et al. [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] have proposed a SEF made of exponential functions:

𝑊 = 𝑎 2𝑏 𝑒𝑥𝑝 [ 𝑏 ( 𝐼 1 -3 )] + ∑ 𝑖=𝑓 ,𝑠 𝑎 𝑖 2𝑏 𝑖 { 𝑒𝑥𝑝 [ 𝑏 𝑖 ( 𝐼 4𝑖 -1 ) 2 ] -1 } + 𝑎 𝑓 𝑠 2𝑏 𝑓 𝑠 { 𝑒𝑥𝑝 [ 𝑏 𝑓 𝑠 𝐼 8𝑓 𝑠 2 ] -1 } (16)
where a, b, 𝑎 𝑓 , 𝑎 𝑠 , 𝑏 𝑓 , 𝑏 𝑠 , 𝑎 𝑓 𝑠 and 𝑏 𝑓 𝑠 are 8 positive material constants, 𝐼 1 is the first principal invariant of C, and 𝐼 4f , 𝐼 4s and 𝐼 8fs are the classical mixed invariants combining C with 𝒇 𝟎 and 𝒔 𝟎 :

𝐼 1 = 𝑇 𝑟 (𝑪) ; 𝐼 4𝑓 = ⟨𝑪𝒇 𝟎 , 𝒇 𝟎 ⟩ ; 𝐼 4𝑠 = ⟨𝑪𝒔 𝟎 , 𝒔 𝟎 ⟩ ; 𝐼 8𝑓 𝑠 = ⟨𝑪𝒇 𝟎 , 𝒔 𝟎 ⟩ (17)
It is easy to check that the standard invariants (17) are connected to the polyconvex invariants 𝐿 𝑖 introduced by Eqs. ( 8)-( 10) by:

𝐼 1 = 𝐿 1 + 𝐿 2 + 𝐿 3 ; 𝐼 4𝑓 = 𝐿 1 ; 𝐼 4𝑠 = 𝐿 2 ; 𝐼 8𝑓 𝑠 = 1 2 √ 𝐿 4 - ( 𝐿 1 + 𝐿 2 ) 2 (18)
And conversely:

𝐿 1 = 𝐼 4𝑓 ; 𝐿 2 = 𝐼 4𝑠 ; 𝐿 3 = 𝐼 1 -𝐼 4𝑓 -𝐼 4𝑠 ; 𝐿 4 = 4𝐼 8𝑓 𝑠 2 + ( 𝐼 4𝑓 + 𝐼 4𝑠 ) 2 (19)
However, it should be also noted that, if the classical invariants 𝐼 1 , 𝐼 4f and 𝐼 4s are polyconvex, as convex functions of 𝑭 , it is not the case for 𝐼 8fs . Actually, from Eq. ( 17), the second derivative of 𝐼 8fs with respect to 𝑭 is:

𝑑 2 𝐼 8𝑓 𝑠 𝑑𝑭 2 (𝑭 ) (𝒅𝑭 ) (𝒅𝑮) = ⟨𝒅𝑭 𝒇 𝟎 , 𝒅𝑮𝒔 𝟎 ⟩ + ⟨𝒅𝑭 𝒔 𝟎 , 𝒅𝑮𝒇 𝟎 ⟩ (20) ⇒ 𝑑 2 𝐼 8𝑓 𝑠 𝑑𝑭 2 (𝑭 ) (𝒅𝑭 ) (𝒅𝑭 ) = 2 ⟨𝒅𝑭 𝒇 𝟎 , 𝒅𝑭 𝒔 𝟎 ⟩ (21)
In order to demonstrate that 𝐼 8fs is non polyconvex, we just have to find a counterexample with a matrix 𝒅𝑭 giving a negative value in Eq. ( 21). We choose this 𝒅𝑭 as follows:

𝒅𝑭 = -𝒇 𝟎 ⊗ 𝒇 𝟎 + 𝒇 𝟎 ⊗ 𝒔 𝟎 + 𝒔 𝟎 ⊗ 𝒇 𝟎 (22)
One reports Eq. ( 22) in the inner product of Eq. ( 21):

⟨𝒅𝑭 𝒇 𝟎 , 𝒅𝑭 𝒔 𝟎 ⟩ = ⟨-𝒇 𝟎 + 𝒔 𝟎 , 𝒇 𝟎 ⟩ = ⟨-𝒇 𝟎 , 𝒇 𝟎 ⟩ = -1 < 0 (23)
The proof of the non-polyconvexity of the classical mixed invariant 𝐼 8fs is complete. Besides, it is underlined in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] that 𝐼 8fs generally has a destabilizing influence. It turns out that if polyconvexity is a prerequisite for building SEF in order to ensure existence of solutions, it would be better to choose the set of invariants defined by Eqs. ( 8)- [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF] rather than the ones defined by Eq. ( 17). In order to make further comparison between the model introduced in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] and our model, it is useful at this stage to remind the Cauchy stress calculated in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF]:

𝝈 = 𝑎 𝑒𝑥𝑝 [ 𝑏 ( 𝐼 1 -3 )] 𝑩 -𝑝𝑰 + 2𝑎 𝑓 ( 𝐼 4𝑓 -1 ) 𝑒𝑥𝑝 [ 𝑏 𝑓 ( 𝐼 4𝑓 -1 ) 2 ] 𝒇 ⊗ 𝒇 + 2𝑎 𝑠 ( 𝐼 4𝑠 -1 ) 𝑒𝑥𝑝 [ 𝑏 𝑠 ( 𝐼 4𝑠 -1 ) 2 ] 𝒔 ⊗ 𝒔 + 𝑎 𝑓 𝑠 𝐼 8𝑓 𝑠 𝑒𝑥𝑝 [ 𝑏 𝑓 𝑠 𝐼 8𝑓 𝑠 2 ] (𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 ) (24)
where 𝒇 and 𝒔 represent the deformed directions, respectively for the fiber and for the direction transverse to the fiber direction in plane of sheets:

𝒇 = 𝑭 𝒇 𝟎 ; 𝒔 = 𝑭 𝒔 𝟎 ; 𝒏 = 𝑭 𝒏 𝟎 ( 25
)
To build our model, we proceed at the beginning in the same manner as in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF], by expressing the SEF in the form of exponential functions such as in Eq. ( 16). However, instead of using the classical invariants introduced by Eq. ( 17), we select the set of polyconvex invariants defined by Eqs. ( 8)- [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]:

𝑊 = 3 ∑ 𝑖=1 𝛼 𝑖 2𝛽 𝑖 { 𝑒𝑥𝑝 [ 𝛽 𝑖 ( 𝐿 𝑖 -1 ) 2 ] - 1
} + 𝛼 4 2𝛽 4 { 𝑒𝑥𝑝 [ 𝛽 4 ( 𝐿 4 -4 ) 2 ] -1 } (26)
The 8 coefficients 𝛼 1 , 𝛼 2 , 𝛼 3 , 𝛼 4 , 𝛽 1 , 𝛽 2 , 𝛽 3 and 𝛽 4 are the material parameters of our model, the 𝛼 𝑖 parameters having the unit of stress while the parameters 𝛽 𝑖 are dimensionless. Note that we have omitted in Eq. ( 26) the dependence of W with respect to the invariants 𝐿 5 and 𝐿 6 in order to keep the same number of invariants as in the model proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF].

It is also remarked that the SEF defined by Eq. ( 26) is null if the material is at rest. Actually, if the displacement is equal to zero, it comes from Eqs. ( 2)-( 3) that both tensors F and C are equal to the identity matrix I . Therefore, it follows from Eqs. ( 8)-( 10) that: 𝐿 1 = 𝐿 2 = 𝐿 3 = 1; 𝐿 4 = 4. Consequently, 𝑊 = 0. We now calculate the Cauchy stress corresponding to our model by reporting Eqs. ( 12)-(15) in Eq. ( 11) and by assuming the incompressible nature of the material (𝐽 = 1):

𝝈 = 2𝑭 {[ 𝜕𝑊 𝜕𝐿 1 + 2 𝜕𝑊 𝜕𝐿 4 ( 𝐿 1 + 𝐿 2 ) ] 𝒇 𝟎 ⊗ 𝒇 𝟎 + [ 𝜕𝑊 𝜕𝐿 2 + 2 𝜕𝑊 𝜕𝐿 4 ( 𝐿 1 + 𝐿 2 ) ] 𝒔 𝟎 ⊗ 𝒔 𝟎 + 𝜕𝑊 𝜕𝐿 3 𝒏 𝟎 ⊗ 𝒏 𝟎 +2 𝜕𝑊 𝜕𝐿 4 √ 𝐿 4 - ( 𝐿 1 + 𝐿 2 ) 2 ( 𝒇 𝟎 ⊗ 𝒔 𝟎 + 𝒔 𝟎 ⊗ 𝒇 𝟎 ) } 𝑭 𝑇 -𝑝𝑰 (27)
We next report the well-known property 𝒏 𝟎 ⊗ 𝒏 𝟎 = 𝑰 -𝒇 𝟎 ⊗ 𝒇 𝟎 -𝒔 𝟎 ⊗ 𝒔 𝟎 in Eq. ( 27) and we use Eq. ( 25):

𝝈 = 2 𝜕𝑊 𝜕𝐿 3 𝑩 -𝑝𝑰 + 2 [ 𝜕𝑊 𝜕𝐿 1 - 𝜕𝑊 𝜕𝐿 3 + 2 𝜕𝑊 𝜕𝐿 4 ( 𝐿 1 + 𝐿 2 ) ] 𝒇 ⊗ 𝒇 + 2 [ 𝜕𝑊 𝜕𝐿 2 - 𝜕𝑊 𝜕𝐿 3 + 2 𝜕𝑊 𝜕𝐿 4 ( 𝐿 1 + 𝐿 2 ) ] 𝒔 ⊗ 𝒔 + 4 𝜕𝑊 𝜕𝐿 4 √ 𝐿 4 - ( 𝐿 1 + 𝐿 2 ) 2 (𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 ) (28)
We finally calculate the derivatives 𝜕𝑊 𝜕𝐿 𝑖 from Eq. ( 26):

𝜕𝑊 𝜕𝐿 𝑖 = 𝛼 𝑖 ( 𝐿 𝑖 -1 ) 𝑒𝑥𝑝 [ 𝛽 𝑖 ( 𝐿 𝑖 -1 ) 2 ] , 𝑖 = 1, 2, 3; 𝜕𝑊 𝜕𝐿 4 = 𝛼 4 ( 𝐿 4 -4 ) 𝑒𝑥𝑝 [ 𝛽 4 ( 𝐿 4 -4 ) 2 ] (29)
In order to compare the Cauchy stress of Eq. ( 28) with Eq. ( 24) extracted from [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF], one reports Eq. (29) in Eq. ( 28):

𝝈 = 𝐴𝑩 -𝑝𝑰 + [𝐵 -𝐴 + 𝐶] 𝒇 ⊗ 𝒇 + [𝐷 -𝐴 + 𝐶] 𝒔 ⊗ 𝒔 + 𝐸 (𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 ) (30)
𝐴 = 2𝛼 3 ( 𝐿 3 -1 ) 𝑒𝑥𝑝 [ 𝛽 3 ( 𝐿 3 -1 ) 2 ] ; 𝐵 = 2𝛼 1 ( 𝐿 1 -1 ) 𝑒𝑥𝑝 [ 𝛽 1 ( 𝐿 1 -1 ) 2 ]
(31)

𝐶 = 4𝛼 4 ( 𝐿 4 -4 ) ( 𝐿 1 + 𝐿 2 ) 𝑒𝑥𝑝 [ 𝛽 4 ( 𝐿 4 -4 ) 2 ] ; 𝐷 = 2𝛼 2 ( 𝐿 2 -1 ) 𝑒𝑥𝑝 [ 𝛽 2 ( 𝐿 2 -1 ) 2 ] ( 32
) 𝐸 = 4𝛼 4 ( 𝐿 4 -4 ) √ 𝐿 4 - ( 𝐿 1 + 𝐿 2 ) 2 𝑒𝑥𝑝 [ 𝛽 4 ( 𝐿 4 -4 ) 2 ] (33)
It is observed that the same 5 tensors 𝑩, 𝑰, 𝒇 ⊗ 𝒇 , 𝒔 ⊗ 𝒔 and 𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 are involved in both equations ( 24) and (30). But there exist fundamental differences since Eq. ( 24) is a sum with a single exponential associated to each term while several exponentials coming from the polyconvex invariants are mixed in the 3 last terms of Eq. ( 30), according to Eqs. ( 31)-(33). Beyond working with polyconvex invariants, the interest of Eq. ( 30) is to account for the invariant 𝐼 4𝑛 = ⟨𝑪𝒏 𝟎 , 𝒏 𝟎 ⟩ and therefore to distinguish explicitly the effect of the stretch in the myocyte axis direction 𝒇 𝟎 , in the direction 𝒔 𝟎 lying within the muscle layer and transverse to 𝒇 𝟎 , and in the direction 𝒏 𝟎 normal to the muscle layer. To prove this claim, Eq. ( 30) is reorganized by factorizing the common terms with 𝐴 and 𝐶:

𝝈 = 𝐴 (𝑩 -𝒇 ⊗ 𝒇 -𝒔 ⊗ 𝒔) + 𝐵𝒇 ⊗ 𝒇 + 𝐷𝒔 ⊗ 𝒔 -𝑝𝑰 + 𝐶 (𝒇 ⊗ 𝒇 + 𝒔 ⊗ 𝒔) + 𝐸 (𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 ) (34)
The first parenthesis in Eq. ( 34) is simplified by reporting inside the following result which holds in the basis (𝒇 𝟎 , 𝒔 𝟎 , 𝒏 𝟎 ):

𝑩 = 𝒇 ⊗ 𝒇 + 𝒔 ⊗ 𝒔 + 𝒏 ⊗ 𝒏 (35) ⇒ 𝝈 = 𝐴𝒏 ⊗ 𝒏 + 𝐵𝒇 ⊗ 𝒇 + 𝐷𝒔 ⊗ 𝒔 -𝑝𝑰 + 𝐶 (𝒇 ⊗ 𝒇 + 𝒔 ⊗ 𝒔) + 𝐸 (𝒇 ⊗ 𝒔 + 𝒔 ⊗ 𝒇 ) (36)
It is noted from Eqs. ( 19), ( 31) and (32) that the coefficients 𝐵 and 𝐷 are linked to the invariants 𝐼 4𝑓 and 𝐼 4𝑠 :

𝐵 = 2𝛼 1 ( 𝐼 4𝑓 -1 ) 𝑒𝑥𝑝 [ 𝛽 1 ( 𝐼 4𝑓 -1 ) 2 ] ; 𝐷 = 2𝛼 2 ( 𝐼 4𝑠 -1 ) 𝑒𝑥𝑝 [ 𝛽 2 ( 𝐼 4𝑠 -1 ) 2 ] (37)
That means that the second and third terms of Eq. (36) represent the effect of the stretch in the directions 𝒇 𝟎 and 𝒔 𝟎 . Regarding the first term of Eq. ( 36), it follows from Eqs. ( 17), ( 19) and (31) that 𝐿 3 is equal to the mixed invariant 𝐼 4𝑛 :

𝐿 3 = 𝑇 𝑟 ( 𝑪 [ 𝑰 -𝒇 𝟎 ⊗ 𝒇 𝟎 -𝒔 𝟎 ⊗ 𝒔 𝟎 ]) = 𝑇 𝑟 ( 𝑪𝒏 𝟎 ⊗ 𝒏 𝟎 ) = 𝐼 4𝑛 (38)
The coefficient 𝐴 is therefore linked to 𝐼 4𝑛 via Eqs. ( 31) and ( 38):

𝐴 = 2𝛼 3 ( 𝐼 4𝑛 -1 ) 𝑒𝑥𝑝 [ 𝛽 3 ( 𝐼 4𝑛 -1 ) 2 ] ( 39
)
That means that the first term of Eq. ( 36) represents the effect of the stretch in the direction 𝒏 𝟎 . It is thus possible to give a physical meaning to each of the 3 first terms of Eq. ( 36) with respect to the specific geometry of the material. This would be useful if one needed to focus specifically on the effect of stretching in the 𝒏 𝟎 direction. This is a significant difference with the model proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] which uses the 3 invariants 𝐼 1 , 𝐼 4𝑓 and 𝐼 4𝑠 instead of 𝐼 4𝑓 , 𝐼 4𝑠 and 𝐼 4𝑛 .

Simple shear test

A shear experiment is driven by the shear deformation k as shown in Fig. 2. The displacement U , related for example to a simple shear deformation in the (𝒏 𝟎 , 𝒇 𝟎 ) plane, is expressed with respect to k as follows:

𝑼 = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑘𝑛 0 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ (40)
It follows from Eqs. ( 2), ( 3) and (40) that:

𝑭 = ⎛ ⎜ ⎜ ⎝ 1 0 𝑘 0 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ; 𝑪 = ⎛ ⎜ ⎜ ⎝ 1 0 𝑘 0 1 0 𝑘 0 𝑘 2 + 1 ⎞ ⎟ ⎟ ⎠ ; 𝑩 = ⎛ ⎜ ⎜ ⎝ 𝑘 2 + 1 0 𝑘 0 1 0 𝑘 0 1 ⎞ ⎟ ⎟ ⎠ (41)
Applying to the 6 different modes of shear the same kind of calculation as the one used for establishing Eqs. ( 40) and (41), it is easy to obtain the displacement 𝑼 , the gradient deformation matrix 𝑭 and the strain tensors 𝑪 and 𝑩 related to each of these 6 modes (Table 1). By using the data contained in Table 1, the expressions of the polyconvex invariants 𝐿 1 , 𝐿 2 , 𝐿 3 and 𝐿 4 with respect to k (Table 2) come from Eqs. ( 8)- [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]. Finally, by using Tables 1 and2, Eq. ( 28) gives the Cauchy stress expressions corresponding to each of the 6 shear modes (Table 3). Since the (nf) and (ns) modes have the same invariant 𝐿 3 , they also have the same shear stresses. This is consistent with the experimental observations since the responses for the (nf) and (ns) modes are indistinguishable as reported in [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF].

Comparison between the model of Eq. (26) and the experimental data of Dokos et al. [1]

In order to identify the 8 material parameters 𝛼 1 , 𝛼 2 , 𝛼 3 , 𝛼 4 , 𝛽 1 , 𝛽 2 , 𝛽 3 and 𝛽 4 of the SEF W defined by Eq. ( 26), we have used the classical coefficient of determination 𝑅 2 :

𝑅 2 = 1 - 𝑆 𝑟𝑒𝑠 𝑆 𝑡𝑜𝑡 ( 42
)
where 𝑆 𝑟𝑒𝑠 and 𝑆 𝑡𝑜𝑡 are respectively the residual sum and the total sum of squares over the number of experimental data n:

𝑆 𝑟𝑒𝑠 = ‖𝒚 -𝒇 ‖ 2 = 𝑛 ∑ 𝑖=1 (𝑦 𝑖 -𝑓 𝑖 ) 2 ; 𝑆 𝑡𝑜𝑡 = ‖ ‖ 𝒚 -𝒚 ‖ ‖ 2 = 𝑛 ∑ 𝑖=1 (𝑦 𝑖 -𝑦) 2 (43)

Table 1

Displacement and strain tensors for the 6 shear modes.

Shear modes

𝑈 𝐹 𝐶 𝐵 (nf) ⎧ ⎪ ⎨ ⎪ ⎩ 𝑘𝑛 0 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 0 𝑘 0 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 𝑘 0 1 0 𝑘 0 𝑘 2 + 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝑘 2 + 1 0 𝑘 0 1 0 𝑘 0 1 ⎞ ⎟ ⎟ ⎠ (ns) ⎧ ⎪ ⎨ ⎪ ⎩ 0 𝑘𝑛 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 𝑘 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 𝑘 0 𝑘 𝑘 2 + 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 𝑘 2 + 1 𝑘 0 𝑘 1 ⎞ ⎟ ⎟ ⎠ (fn) ⎧ ⎪ ⎨ ⎪ ⎩ 0 0 𝑘𝑓 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 0 𝑘 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝑘 2 + 1 0 𝑘 0 1 0 𝑘 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 𝑘 0 1 0 𝑘 0 𝑘 2 + 1 ⎞ ⎟ ⎟ ⎠ (fs)
⎧ ⎪ ⎨ ⎪ ⎩ 0 𝑘𝑓 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 0 0 𝑘 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝑘 2 + 1 𝑘 0 𝑘 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 𝑘 0 𝑘 𝑘 2 + 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ (sn) ⎧ ⎪ ⎨ ⎪ ⎩ 0 0 𝑘𝑠 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 0 0 𝑘 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 𝑘 2 + 1 𝑘 0 𝑘 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 0 0 0 1 𝑘 0 𝑘 𝑘 2 + 1 ⎞ ⎟ ⎟ ⎠ (sf) ⎧ ⎪ ⎨ ⎪ ⎩ 𝑘𝑠 0 0 0 ⎫ ⎪ ⎬ ⎪ ⎭ ⎛ ⎜ ⎜ ⎝ 1 𝑘 0 0 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 1 𝑘 0 𝑘 𝑘 2 + 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ 𝑘 2 + 1 𝑘 0 𝑘 1 0 0 0 1 ⎞ ⎟ ⎟ ⎠

Table 2

Polyconvex invariants for the 6 shear modes.

Shear modes Polyconvex invariants

𝐿 1 𝐿 2 𝐿 3 𝐿 4 (nf) 1 1 𝑘 2 + 1 4 (ns) 1 1 𝑘 2 + 1 4 (fn) 𝑘 2 + 1 1 1 ( 𝑘 2 + 2 ) 2 (fs) 𝑘 2 + 1 1 1 ( 𝑘 2 + 2 ) 2 + 4𝑘 2 (sn) 1 𝑘 2 + 1 1 ( 𝑘 2 + 2 ) 2 (sf) 1 𝑘 2 + 1 1 ( 𝑘 2 + 2 ) 2 + 4𝑘 2

Table 3

Shear component of the Cauchy stresses for the 6 shear modes.

Shear modes

Cauchy stress -shear component

(nf) 2𝑘 𝜕𝑊 𝜕𝐿 3 (ns) 2𝑘 𝜕𝑊 𝜕𝐿 3 (fn) 2𝑘 [ 𝜕𝑊 𝜕𝐿 1 + 2 ( 𝑘 2 + 2 ) 𝜕𝑊 𝜕𝐿 4 ] (fs) 2𝑘 [ 𝜕𝑊 𝜕𝐿 1 + 2 ( 𝑘 2 + 4 ) 𝜕𝑊 𝜕𝐿 4 ] (sn) 2𝑘 [ 𝜕𝑊 𝜕𝐿 2 + 2 ( 𝑘 2 + 2 ) 𝜕𝑊 𝜕𝐿 4 ] (sf) 2𝑘 [ 𝜕𝑊 𝜕𝐿 2 + 2 ( 𝑘 2 + 4 ) 𝜕𝑊 𝜕𝐿 4 ]
𝑦 𝑖 stands for the experimental data, 𝑓 𝑖 for the theoretical data and 𝑦 for the mean of the experimental data: . The experimental data were extracted from [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF] and the theoretical values were computed by using the stress expressions of Table 3 in combination with Eq. ( 29):

𝑦 = 1 𝑛 𝑛 ∑ 𝑖=1 𝑦 𝑖 (44)
𝜎 𝑛𝑓 = 𝜎 𝑛𝑠 = 2𝑘 3 𝛼 3 𝑒𝑥𝑝 ( 𝛽 3 𝑘 4 ) (45)
𝜎 𝑓 𝑛 = 2𝑘 (47)

𝜎 𝑠𝑛 = 2𝑘 3 [ 𝛼 2 𝑒𝑥𝑝 ( 𝛽 2 𝑘 4 ) + 2𝛼 4 ( 𝑘 2 + 2 ) ( 𝑘 2 + 4 ) 𝑒𝑥𝑝 ( 𝛽 4 𝑘 4 ( 𝑘 2 + 4 ) 2 )] (48)
𝜎 𝑠𝑓 = 2𝑘 3 [ 𝛼 2 𝑒𝑥𝑝 ( 𝛽 2 𝑘 4 ) + 2𝛼 4 ( 𝑘 2 + 4 ) ( 𝑘 2 + 8 ) 𝑒𝑥𝑝 ( 𝛽 4 𝑘 4 ( 𝑘 2 + 8 ) 2 )] (49)
We have implemented equations Eqs. ( 45)-(49) in the Octave free software [START_REF] Eaton | GNU Octave version 4.2.0 manual: a high-level interactive language for numerical computations[END_REF] and performed the material identification by using the optimization routine fminunc based on gradient search. The identified material parameters are shown in Table 4 and the corresponding coefficients of determination 𝑅 2 in Table 5. It is observed an excellent level of correlation, confirmed by Fig. 3 which shows a very good fit between the theoretical curves and the experimental data. That proves that a model using polyconvex invariants is as efficient as a model using standard invariants for predicting the mechanical response of passive ventricular myocardium. The worst prediction, with a coefficient of determination equal to 0.982, is obtained with the (nf) = (ns) modes, but it is not surprising because these modes concern the lowest level of stress and it is generally difficult for a model to match perfectly experimental low values. Finally, we have observed that the theoretical curves drawn in Fig. 3 are identical to the ones provided in Figure 6 of [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF], proving the consistency of our model with the literature.

Convexity of the SEF defined by Eq. (26)

The convexity of a SEF is an important issue regarding computational efficiency, and particularly for ensuring the uniqueness of the solution of a hyperelastic problem. In order to present sufficient conditions providing this desired property of convexity, we proceeded at the beginning in the same way as in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF]. Since the SEF defined by Eq. ( 26) is expressed as a sum over exponential functions, a sufficient condition for ensuring convexity is to study separately the convexity for each term. The contribution in the SEF of the first invariant 𝐿 1 is for example written as follows from Eq. (26):

𝑊 1 (𝑪) = 𝑓 1 ( 𝐿 1 (𝑪) ) (50)
With the real-valued function 𝑓 1 defined by:

𝑓 1 (𝑥) = 𝛼 1 2𝛽 1 { 𝑒𝑥𝑝 [ 𝛽 1 (𝑥 -1) 2 ] -1 } (51)
From Eq. ( 8)-( 10), the invariant 𝐿 1 (𝑪) can be expressed linearly with respect to 𝑪 thanks to the trace operator 𝑇 𝑟:

𝐿 1 (𝑪) = ⟨𝑪𝒇 𝟎 , 𝒇 𝟎 ⟩ = 𝑇 𝑟 ( 𝒇 𝟎 ⊗ 𝒇 𝟎 𝑪 ) (52)
Deriving twice Eq. (50), one obtains from Eq. (52):

𝑊 1 ′′ (𝑪) (𝒅𝑪) (𝒅𝑪) = 𝑓 1 ′′ ( 𝐿 1 (𝑪) ) [ 𝑇 𝑟 ( 𝒇 𝟎 ⊗ 𝒇 𝟎 𝒅𝑪 )] 2 (53)
where the second derivative of 𝑓 1 is obtained from Eq. (51):

𝑓 1 ′′ ( 𝐿 1 (𝑪) ) = 𝛼 1 { 1 + 2𝛽 1 ( 𝐿 1 (𝑪) -1 ) 2 } 𝑒𝑥𝑝 [ 𝛽 1 ( 𝐿 1 (𝑪) -1 ) 2 ] (54)
It is deduced from Eqs. ( 53) and (54) that a sufficient condition for 𝑊 1 to be strictly convex is:

𝛼 1 > 0; 𝛽 1 > 0 (55)
Following the same logic with the second and third polyconvex invariants 𝐿 2 and 𝐿 3 , it is obtained the additional sufficient conditions:

𝛼 2 > 0; 𝛼 3 > 0; 𝛽 2 > 0; 𝛽 3 > 0 (
Noting that 𝐿 4 ′′ (𝑪) (𝒅𝑪) (𝒅𝑪) ≥ 0, as a sum over positive quantities, and that 𝐿 4 -4 ≥ 0 in each of the 6 shear loadings studied in this work (last row of Table 2), it results from Eqs. ( 60)-( 63) that a sufficient condition for 𝑊 4 to be convex is:

𝛼 4 > 0; 𝛽 4 > 0 (64)
Note that the material parameters identified in Table 4 satisfy the convexity conditions given by Eqs. ( 55), ( 56) and (64).

Conclusions

In this study, a new strain energy function (SEF) was proposed for modeling incompressible orthotropic hyperelastic materials with a specific application to the mechanical response of passive ventricular myocardium. The SEF was built by combining exponential functions as proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] but, instead of using standard mixed invariants, we have selected the integrity basis of polyconvex invariants introduced in [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF].

Using this new set of invariants allows to replace the classical mixed invariant 𝐼 8 , which is demonstrated to be non-polyconvex, by the polyconvex invariant 𝐿 4 defined by Eqs. ( 8)- [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF]. With a SEF built on this basis, according to Eqs. [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF] and (26), we demonstrated that our model is capable to accurately predict the experimental data obtained in [START_REF] Dokos | Shear properties of passive ventricular myocardium[END_REF] with 6 different shear modes applied to passive ventricular myocardium. In addition, our model is also consistent with results previously obtained in the literature [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF]. Besides, we established sufficient conditions ensuring the convexity of the SEF and expressed them as inequalities involving the material parameters (Eqs. (55), ( 56) and ( 64)). Finally, our model can explicitly represent the stretch effect in the 3 specific directions of a passive ventricular myocardium: the myocyte axis direction 𝒇 𝟎 , the direction 𝒔 𝟎 lying within the muscle layer and transverse to 𝒇 𝟎 , and the direction 𝒏 𝟎 normal to the muscle layer. We demonstrated that 3 of the polyconvex invariants are equal to the squared stretches 𝐼 4𝑓 , 𝐼 4𝑠 and 𝐼 4𝑛 . Consequently, three terms directly connected to the 3 possible stretch effects are embedded in the stress tensor (Eq. ( 34)). This is a significant difference with the model proposed in [START_REF] Holzapfel | Constitutive modelling of passive myocardium: a structurally based framework for material characterization[END_REF] which uses the 3 invariants 𝐼 1 , 𝐼 4𝑓 and 𝐼 4𝑠 instead of 𝐼 4𝑓 , 𝐼 4𝑠 and 𝐼 4𝑛 .

Because polyconvexity is often considered as a prerequisite for ensuring the existence of solutions, the next step of our research work will be to implement the new SEF introduced in this paper in a university finite element code.

CRediT authorship contribution statement

Fig. 1 .

 1 Fig. 1. 6 possible modes of simple shear test in different planes.

Fig. 2 .

 2 Fig. 2. Simple shear test in the (𝒏 𝟎 , 𝒇 𝟎 ) plane.

Fig. 3 .

 3 Fig. 3. Comparison between the model and the experimental data.

  Cai: Methodology, Software, Formal Analysis, Investigation, Writing -review & editing. Frédéric Holweck: Conceptualization, Writing -review & editing, Supervision. Zhi-Qiang Feng: Conceptualization, Writing -review & editing, Supervision. François Peyraut: Conceptualization, Methodology, Software, Formal analysis, Writingoriginal draft, Writing -review & editing, Visualization, Supervision, Project administration.

around 𝒇 𝟎 , 𝒔 𝟎 and 𝒏 𝟎 ,

  and under 𝑰 and -𝑰:

	reflections 𝑹 related to these 3 planes; 3 rotations 𝑸 𝝅 ( 𝒇 𝟎 ) , 𝑸 𝝅 ( 𝒇 𝟎 ) , 𝑹 ( 𝒔 𝟎 ) and 𝑸 𝝅 ( 𝒔 𝟎 ) and 𝑹 ( 𝒏 𝟎	( 𝒏 𝟎 ) by )
	an angle 𝜋 𝑹 ( 𝒇 𝟎 ) = 𝑰 -2𝒇 𝟎 ⊗ 𝒇 𝟎 = -𝑸 𝝅	( 𝒇 𝟎	)	; 𝑹	( 𝒔 𝟎	)	= 𝑰 -2𝒔 𝟎 ⊗ 𝒔 𝟎 = -𝑸 𝝅	( 𝒔 𝟎	)
	𝑹 The set 𝑆 8 = ( 𝒏 𝟎 ) = 𝑰 -2𝒏 𝟎 ⊗ 𝒏 𝟎 = -𝑸 𝝅 { 𝑹 ( 𝒇 𝟎 ) , 𝑹	( 𝒏 𝟎 ( 𝒔 𝟎	) )	, 𝑹	( 𝒏 𝟎	)	, 𝑸 𝝅	(	𝒇 𝟎	)	, 𝑸 𝝅	( 𝒔 𝟎	)	, 𝑸 𝝅	(6) (7) 𝒏 𝟎 ( ) ,
	𝑰, -𝑰} built with the 8 invariant matrix operators defined by Eqs. (

Table 4

 4 Identified material parameters.

	𝛼 1 (kPa)	𝛼 2 (kPa)	𝛼 3 (kPa)	𝛼 4 (kPa)
	18.877	2.495	3.184	0.168
	𝛽 1 (-)	𝛽 2 (-)	𝛽 3 (-)	𝛽 4 (-)
	19.39	20.113	11.543	0.107

Table 5

 5 Coefficient of determination 𝑅 2 . 𝑅 2 is, the best the fit of the experimental data by the theoretical data will be. The aim is thus to find the set of material parameters( 𝛼 1 , 𝛼 2 , 𝛼 3 , 𝛼 4 , 𝛽 1 , 𝛽 2 , 𝛽 3 , 𝛽 4) minimizing the ratio

	Shear modes	(nf)=(ns)	(fn)	(fs)	(sn)	(sf)
	𝑅 2	0.982	0.998	0.997	0.993	0.998
	The closest to 1 𝑆 𝑟𝑒𝑠
						𝑆 𝑡𝑜𝑡

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Valentina Zantedeschi
email: vzantedeschi@gmail.com

Paul Viallard

Emilie Morvant

Rémi Emonet

Amaury Habrard

Pascal Germain
email: pascal.germain@ift.ulaval.ca

Benjamin Guedj
email: benjamin.guedj@inria.fr

Remi Emonet

Learning Stochastic Majority Votes by Minimizing a PAC-Bayes Generalization Bound

come   L'archive ouverte pluridisciplinaire

Introduction

By combining the outcomes of several predictors, ensemble methods [START_REF] Thomas | Ensemble methods in machine learning[END_REF] have been shown to provide models that are more accurate and more robust than each predictor taken singularly. The key to their success lies in harnessing the diversity of the set of predictors [START_REF] Ludmila | Combining Pattern Classifiers: Methods and Algorithms[END_REF]. Among ensemble methods, weighted Majority Votes (MV) classifiers assign a score to each base classifier (a.k.a. voter) and output the most common prediction, given by the weighted majority. When voters have known probabilities of making an error and make independent predictions, the optimal weighting is given by the so-called Naive Bayes rule [START_REF] Berend | A finite sample analysis of the naive Bayes classifier[END_REF]. However, in most situations these assumptions are not satisfied, giving rise to the need for techniques that estimate the optimal combination of voter predictions from the data.

Among them, PAC-Bayesian based methods are well-grounded approaches for optimizing the voter weighting. Indeed, PAC-Bayes theory (introduced by [START_REF] Shawe | A PAC analysis of a bayesian estimator[END_REF][START_REF] Mcallester | Some PAC-Bayesian theorems[END_REF]-we refer to Guedj, 2019 for a recent survey and references therein) provides not only bounds on the true error of a MV through Probably Approximately Correct (PAC) generalization bounds (see e.g. [START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF], [START_REF] Seeger | PAC-Bayesian generalisation error bounds for gaussian process classification[END_REF], [START_REF] Maurer | A note on the PAC Bayesian theorem[END_REF], [START_REF] Langford | PAC-Bayes & margins[END_REF], [START_REF] Germain | Risk bounds for the majority vote: from a PAC-Bayesian analysis to a learning algorithm[END_REF]), but is also suited to derive theoretically grounded learning algorithms (see e.g. [START_REF] Germain | Pac-bayesian learning of linear classifiers[END_REF], [START_REF] Roy | From PAC-Bayes bounds to quadratic programs for majority votes[END_REF], [START_REF] Parrado-Hernández | Pac-bayes bounds with data dependent priors[END_REF], [START_REF] Alquier | On the properties of variational approximations of Gibbs posteriors[END_REF]). Contrary to the most classical PAC bounds [START_REF] Leslie | A theory of the learnable[END_REF], as VC-dimension [Vapnik, 2000] or Rademacher-based bounds [START_REF] Mohri | Foundations of Machine Learning[END_REF], PAC-Bayesian guarantees do not stand for all hypotheses (i.e. are not expressed as a worst-case analysis) but stand in expectation over the hypothesis set. They involve a hypothesis space (formed by the base predictors), a prior distribution on it (i.e. an a priori weighting) and a posterior distribution (i.e. an a posteriori weighting) evaluated on a learning sample. The prior brings some prior knowledge on the combination of predictors, and the posterior distribution is learned (adjusted) to lead to good generalization guarantees; the deviation between the prior and the posterior distributions plays a role in generalization guarantee and is usually captured by the Kullback-Leibler (KL) divergence. In their essence, PAC-Bayesian results do not bound directly the risk of the deterministic MV, but bound the expected risk of one (or several) base voters randomly drawn according to the weight distribution of the MV [START_REF] Langford | PAC-Bayes & margins[END_REF], Lacasse et al., 2006[START_REF] Lacasse | Learning with randomized majority votes[END_REF][START_REF] Germain | Risk bounds for the majority vote: from a PAC-Bayesian analysis to a learning algorithm[END_REF][START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF]. This randomization scheme leads to upper bounds on the true risk of the MV that are then used as a proxy to derive PAC-Bayesian generalization bounds. However, the obtained risk certificates are generally not tight, as they depend on irreducible constant factors, and when optimized they can lead to sub-optimal weightings. Indeed, by considering a random subset of base predictors, state-of-the-art methods do not fully leverage the diversity of the whole set of voters. This is especially a problem when the voters are weak, and learning to combine their predictions is critical for good performance.

Our contributions. In this paper, we propose a new randomization scheme. We consider the voter weighting associated to a MV as a realization of a distribution of voter weightings. More precisely, we analyze with the PAC-Bayesian framework the expected risk of a MV drawn from the posterior distribution of MVs. The main difference with the literature is that we propose a stochastic MV, while previous works aim at studying randomized evaluations of the true risk of the deterministic MV. Doing so, we are able to derive tight empirical PAC-Bayesian bounds for our model directly on its expected risk, in Section 4. We further propose, in Section 3 two approaches for optimizing the generalization bounds, hence learning the optimal posterior: the first optimizes an analytical and differentiable form of the empirical risk that can be derived when considering Dirichlet distributions; the second optimizes a Monte Carlo approximation of the expected risk and can be employed with any form of posterior. In our experiments, reported in Section 5, we first compare these two approaches, highlighting in which regimes one is preferable to the other. Finally, we assess our method's performance on real benchmarks w.r.t. the performance of PAC-Bayesian approaches also learning MV classifiers. These results indicate that our models enjoy generalization bounds that are consistently tight and non-vacuous both when studying ensembles of data-independent predictors and when studying ensembles of data-dependent ones.

Societal impact. Our work abides by the ethical guidelines enforced in contemporary research in machine learning. Given the theoretical nature of our contributions we do not foresee immediate potentially negative societal impact.

Notation and background

In this section, we formally define weighted Majority Vote (MV) classifiers and review the principal PAC-Bayesian approaches for learning them.

Weighted majority vote classifiers

Consider the data random variable pX, Y q, taking values in X ˆY with X Ď R d a d-dimensional representation space and Y the set of labels. We denote P the (unknown) data distribution of pX, Y q. We define a set (dictionary) of base classifiers D"th j : X Ñ Yu M j"1 . The weighted majority vote classifier is a convex combination of the base classifiers from D. Formally, a MV is parameterized by a weight vector θ P r0, 1s M , such that ř M j"1 θ j " 1, as follows:

f θ pxq " argmax yPY M ÿ j"1 θ j 1ph j pxq " yq, ( 1
)
where 1p¨q is the indicator function. Let W θ pX, Y q be the random variable corresponding to the total weight assigned to base classifiers that predict an incorrect label on pX, Y q, that is

W θ pX, Y q " M ÿ j"1 θ j 1ph j pXq ‰ Y q. (2)
In binary classification with |Y|"2, the MV errs whenever W θ pX, Y q ě 0.5 [START_REF] Lacasse | Learning with randomized majority votes[END_REF][START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF]. Hence the true risk (w.r.t. 01-loss) of the MV classifier can be expressed as Rpf θ q " E P 1pW θ pX, Y q ě 0.5q fi PpW θ ě 0.5q.

(3)

Similarly, the empirical risk of f θ on a n-sample S"tpx i , y i q"Pu n i"1 is given by Rpf θ q " 1 n ř n i"1 1pW θ px i , y i q ě 0.5q. Note that the results we introduce in the following are stated for binary classification, but are valid also in the multi-class setting (|Y|ą2). Indeed, in this context Equation (3) becomes a surrogate of the risk: we have Rpf θ q ď PpW θ ě 0.5q (see the notion of ω-margin with ω " 0.5 proposed by [START_REF] Laviolette | Risk upper bounds for general ensemble methods with an application to multiclass classification[END_REF]).

A PAC-Bayesian perspective on the majority vote

The current paper stands in the context of PAC-Bayesian MV learning. The PAC-Bayesian framework has been employed to study the generalization guarantees of randomized classifiers. It is known to provide tight risk certificates that can be used to derive self-bounding learning algorithms [e.g. [START_REF] Germain | Pac-bayesian learning of linear classifiers[END_REF][START_REF] Gintare | Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data[END_REF][START_REF] Pérez-Ortiz | Tighter risk certificates for neural networks[END_REF]. In the PAC-Bayes literature of the analysis of MV, θ is interpreted as the parameter of a categorical distribution Cpθq over the set of base classifiers D [e.g. [START_REF] Germain | Risk bounds for the majority vote: from a PAC-Bayesian analysis to a learning algorithm[END_REF][START_REF] Sloth Lorenzen | On PAC-Bayesian bounds for random forests[END_REF][START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF]. In this sense, f θ corresponds to the MV predictor f θ pxq" argmax yPY E h"Cpθq 1phpxq"yq and W θ corresponds to the expected ratio of wrong predictors W θ pX, Y q " E h"Cpθq 1phpXq ‰ Y q. The PAC-Bayesian analysis provides a sensible way to find such a categorical distribution, called posterior distribution, that leads to a model with low true risk Rpf θ q. However, an important caveat is that the PAC-Bayesian generalization bounds cannot be derived directly on the true risk Rpf θ q, without making assumptions on the distribution of W θ and raising fidelity problems. Thus, a common approach is to consider oracle1 upper bounds on the true risk in terms of quantities from which PAC-Bayesian bounds can be derived, that are typically related to the statistical moments of W θ . By doing so, oracle bounds act as a proxy for estimating the cumulative density function for W θ ě 0.5 [START_REF] Langford | PAC-Bayes & margins[END_REF], Germain et al., 2006[START_REF] Lacasse | Learning with randomized majority votes[END_REF][START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF]. Generalization guarantees for Rpf θ q are hence derived, involving the empirical counterpart of the oracle bound, the KLdivergence between the posterior categorical distribution Cpθq and a prior one. An overview of the existing oracle bounds described below is represented in Figure 1.

First order bound. The most classical "factor two" oracle bound [START_REF] Langford | PAC-Bayes & margins[END_REF] is derived considering the relation between MV's risk and the first moment of W θ , a.k.a. the expected risk of the randomized classifier: R 1 pθqfiE h"Cpθq Rphq"E P W θ . By applying Markov's inequality, we have Rpf θ q ď 2R 1 pθq. This "factor two" bound is close to zero only if the expectation of the risk of a single base classifiers drawn according to Cpθq is itself close to zero. Therefore, it does not take into account the correlation between base predictors, which is key to characterize how a MV classifier can achieve strong predictions even when its base classifiers are individually weak (as observed when performing e.g. Boosting [START_REF] Schapire | Improved boosting algorithms using confidence-rated predictions[END_REF]). This explains why R 1 pθq can be a very loose estimate of Rpf θ q when the base classifiers are adequately decorrelated.

Randomized MV. A generalization of the first order approach was proposed in [START_REF] Lacasse | Learning with randomized majority votes[END_REF], where the true risk of the MV is estimated by drawing multiple (N ) base hypotheses and computing the probability that at least N 2 make an error (which is given by the binomial with parameter W θ :

W θ N pX, Y q " ř N k" N 2 `N k ˘W k θ p1
´Wθ q pN ´kq . The higher N , the better W θ N pX, Y q approximates the true risk, but the looser the bound, as the KL term worsens by a factor of N . Moreover, with this approach it is not possible to derive generalization bounds directly on the true risk, and the corresponding oracle bound still presents a factor two:

Rpf θ q ď 2 E P W θ N .
Second order bound. A parallel line of works focuses on improving the bounds by accounting for voter correlations, i.e. considering the agreement and/or disagreement of two random voters. [START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF] recently proposed a new upper bound for the true risk depending on the second moment of W θ , a.k.a. tandem loss or joint error, the risk that two predictors make a mistake on the same point: R 2 pθq"E P E h"Cpθq,h 1 "Cpθq 1phpXq‰Y ^h1 pXq‰Y q"E P W 2 θ . By applying the second order Markov's inequality, we have Rpf θ qď4R 2 pθq. Masegosa et al. [2020] show that in the worst case (i.e. when the base classifiers are identical) the second order bound could be twice worse than the first order bound (R 2 pθq«2R 1 pθq), while in the best case (i.e. the M base classifiers are perfectly decorrelated), the second order bound could be an order of magnitude tighter (R 2 pθq« 1 M R 1 pθq).

C-bound. Originally proposed by [START_REF] Breiman | Random forests[END_REF] in the context of Random Forest, a bound known as the C-Bound in the PAC-Bayesian literature [START_REF] Lacasse | PAC-Bayes bounds for the risk of the majority vote and the variance of the gibbs classifier[END_REF][START_REF] Roy | From PAC-Bayes bounds to quadratic programs for majority votes[END_REF][START_REF] Germain | Risk bounds for the majority vote: from a PAC-Bayesian analysis to a learning algorithm[END_REF][START_REF] Laviolette | Risk upper bounds for general ensemble methods with an application to multiclass classification[END_REF] is derived by considering explicitly the joint error and disagreement between two base predictors. Using Chebyshev-Cantelli's inequality, the C-Bound can be written in terms of the first and second moments of

W θ : Rpf θ q ď R2pθq´E h"Cpθq pEP 1phpXq‰Y qq 2 R2pθq´R1pθq`1 4
. This bound is tighter than the second order one if the disagreement E h,h 1 E P p1phpxq ‰ h 1 pXqqq is greater than R 2 pθq [START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF].

From a practical point of view, one could minimize the generalization bounds of one of the above methods to learn a weight distribution over an ensemble of predictors. However, this could lead to sub-optimal MV classifiers. To illustrate this behavior we plot in Figure 3 the decision surfaces learnt by the minimization of a PAC-Bayesian bound on each of the aforementioned oracle bounds. These plots provide evidence that, when the base classifiers are weak, state-of-the-art PAC-Bayesian methods do not necessarily build powerful ensembles (failing to improve upon a Naive Bayes approach [START_REF] Berend | A finite sample analysis of the naive Bayes classifier[END_REF] [2020], while Second Order and C-Bound fail to leverage the diversity in the whole set of classifiers. On the contrary, Randomized is able to fit the problem by drawing more than 2 voters, but it provides generalization bounds that are loose even when the learned model exhibits good generalization capabilities, as in this case.

3 Our approach: stochastic weighted majority vote

In our framework for deriving generalization bound for Majority Vote classifiers, we consider f θ as a realization of a distribution of majority vote classifiers, with parameter θ P Θ Ď R M and probability measure ρ, as represented in Figure 2. The main advantage of considering a stochastic majority vote is that it allows to derive and optimize PAC-Bayesian generalization bounds directly on its true risk and to fully leverage the set of base classifiers. The true risk of the proposed stochastic weighted majority vote takes into account the whole distribution of MVs ρpθq, as follows:

ż Θ Rpf θ qρpdθq " ż Θ E P 1pW θ pX, Y q ě 0.5qρpdθq " E P ż Θ 1pW θ pX, Y q ě 0.5qρpdθq. (4)
Given its stochastic nature, in order to evaluate and/or minimize Equation (4) we can either (i) compute its closed form or (ii) approximate it (e.g. through Monte Carlo methods). In both cases, assumptions have to be made on the form of ρ. As the components of θ are constrained to sum to one, θ lies in the (M -1)-simplex: Θ " ∆ pM -1q , hence natural choices for its probability measure are e.g. the Dirichlet or the Logit Normal distributions. In the following, we show that under Dirichlet assumptions, we can derive an analytical and differentiable form of the risk of a stochastic MV.

Exact risk under Dirichlet assumptions

First of all, we recall that the probability density function of the Dirichlet distribution is defined by:

θ " Dpα 1 , . . . , α M q, ρpθq " 1 Bpαq M ź j"1 pθ j q αj ´1, (5)
with α " rα j P R `sM j"1 the vector of concentration parameters and Bpαq a normalization factor (see App A for its definition). Notice that by taking α as the unit vector, the distribution corresponds to a uniform distribution over the (M -1)-simplex.

Under these assumptions for the MV distribution, an closed form can be derived for the expected risk:

Lemma 1. For a given px, yq " P, let w " tj|h j pxq ‰ yu be the set of indices of the base classifiers that misclassify px, yq and c " tj|h j pxq " yu be the set of indices of the base classifiers that correctly classify px, yq. The expected error (or 01-loss) for px, yq of the stochastic majority vote under θ " Dpα 1 , . . . , α M q is equal to

ż Θ 1pW θ px, yq ě 0.5qρpdθq " I 0.5 ˜ÿ jPc α j , ÿ jPw α j ¸, (6)
with I 0.5 p¨q the regularized incomplete beta function evaluated at 0.5.

Proof. We rewrite W θ as W θ px, yq " ř M j"1 θ j 1ph j pxq ‰ yq " ř jPw θ j , and use the aggregation property of Dirichlet distributions (Lemma 2 in Appendix) to show that W θ follows a bivariate Dirichlet distribution (a.k.a. Beta distribution) over the two sets of wrong and correct base classifiers:

PrW θ px, yq " ωs " P « ÿ jPw θ j " ω ff P « ÿ jPc θ j " 1 ´ωff " D ˜ÿ jPw α j , ÿ jPc α j ¸by aggregation. (7)
Finally, notice that the expected error is related to the cumulative probability function of W θ , the incomplete beta function I p : R `ˆR `Ñ r0, 1s:

ż Θ 1pW θ px, yq ě 0.5qρpdθq " ż 1 0.5 PrdW θ px, yqs (8)
" 1 ´I0.5 ˜ÿ jPw α j , ÿ jPc α j ¸" I 0.5 ˜ÿ jPc α j , ÿ jPw α j ¸. (9)
Equation ( 9) follows by symmetry of the incomplete beta function: I p pa, bq " 1 ´I1´p pb, aq.

The expected risk can be then expressed as follows:

Rpρq

" E P ż Θ 1pW θ px, yq ě 0.5qρpdθq " E P I 0.5 ˜ÿ jPc α j , ÿ jPw α j ¸. (10)
Importantly, this exact form of the risk is differentiable, hence can be directly optimized by gradient-based methods.

Monte Carlo approximated risk

We now propose a relaxed Monte Carlo (MC) optimization scheme for those distributions that do not admit an analytical form of the expected risk, unlike the Dirichlet one. This second strategy is also suited to speed up training, in some cases, as the derivatives of the exact risk depend on the hyper-geometric function and can be slow to evaluate (see App. A.4). With the approximated scheme, in order to update α by gradient descent we need to relax the true risk as the gradients of the 01-loss are always null for discrete W θ . In practice, we make use of a tempered sigmoid loss σ c pxq "

2. Compute the relaxed empirical risk

ř T t"1 Rσc pθq " ř T t"1 ř n i"1 σ c pW θ px i , y i q ´0.5q;
3. Update α by gradient descent.

Notice that when considering Dirichlet distributions for the posterior and the prior, at inference time the empirical PAC-Bayesian bounds can still be evaluated using the exact form of Lemma 1.

A drawback of the approximated scheme is that it has a complexity linear in the number of MC draws T , but also linear in the number of predictors M , as sampling over the simplex requires sampling from OpM q distributions, one per base classifier. As an example, sampling from a Dirichlet over the (M -1)-simplex is usually implemented as sampling from M Gamma distributions and normalizing the samples so that they lie on the simplex. In contrast, the exact formulation's complexity is constant in M as it depends only on the sets of wrong and of correct predictors, hence on a constant number of variables (2) no matter the number of predictors M . In Section 5.1 we empirically study the trade-off between training time and accuracy, showing in which regimes it is more convenient to optimize the relaxed MC risk than optimizing the exact one, and viceversa.

PAC-Bayesian generalization guarantees

We now derive PAC-Bayesian generalization upper bounds for the proposed stochastic MV. In our context, upper bounds can be derived for studying the gap between true and empirical risk considering a prior distribution π over the hypothesis space Θ. In this paper, we make use of one of the tightest classical PAC-Bayesian bound [START_REF] Seeger | PAC-Bayesian generalisation error bounds for gaussian process classification[END_REF][START_REF] Maurer | A note on the PAC Bayesian theorem[END_REF]:

Theorem 1 (Seeger's bound). For any π over Θ and δ P p0, 1q with probability at least 1´δ over samples S " tpx i , y i q"Pu n i"1 of size n we have simultaneously for any posterior ρ over Θ:

ż Θ Rpf θ qρpdθq ď kl ´1 ¨żΘ Rpf θ qρpdθq , KLpρ, πq `ln ´2? n δ n ', (11)
with Rpf θ q " 1 n ř n i"1 1pW θ px i , y i q ě 0.5q the empirical risk on sample S, KLpρ, πq " ş Θ ρpθq log ρpθq πpθq dθ the KL divergence and kl ´1pq, q the inverse of the binary KL divergence defined as kl ´1pq, q " maxtp P r0, 1s | klpq, pq ď u.

A proof of Theorem 1 can be found in [START_REF] Seeger | PAC-Bayesian generalisation error bounds for gaussian process classification[END_REF]. The kl ´1 function can be evaluated via the bisection method and optimized by gradient descent, as proposed in [START_REF] Reeb | Learning Gaussian Processes by Minimizing PAC-Bayesian Generalization Bounds[END_REF]. Note that our contributions do not restrict the choice of generalization bound.

Importantly, Theorem 1 is valid when the prior π is independent from the data. Thus it cannot be evaluated with base classifiers learned from the training sample. However, it is known that considering a data-dependent prior can lead to tighter PAC-Bayes bounds [START_REF] Gintare | On the role of data in pac-bayes[END_REF]. Following recent works on PAC-Bayesian bounds wit data-dependent priors [START_REF] Thiemann | A strongly quasiconvex PAC-Bayesian bound[END_REF][START_REF] Mhammedi | PAC-Bayes un-expected bernstein inequality[END_REF], we derive a cross-bounding certificate that allows us to learn and evaluate the set of base classifiers without held-out data. More precisely, we split the training data S into two subsets (S ďm "tpx i , y i q P Su m i"1 and S ąm "tpx i , y i q P Su n i"m`1 ) and we learn a set of base classifiers on each data split independently (determining the hypothesis spaces Θ ďm and Θ ąm ). We refer to the prior distribution over D ďm as π ďm and to the prior distribution over D ąm as π ąm . In the same way, we can then define a posterior distribution per hypothesis space: ρ ďm and ρ ąm . The following theorem shows that we can bound the expected risk of any convex combination of the two posteriors, as long as their empirical risks are evaluated on the data split that was not used for learning their respective priors.

Theorem 2 (Seeger's bound with informed priors). Let π ďm and ρ ďm be the prior and posterior distributions on Θ ďm , and π ąm and ρ ąm the prior and posterior distributions on Θ ąm . For any pPp0, 1q and δ P p0, 1q with probability at least 1´δ over samples S"tpx i , y i q"Pu n i"1 we have kl ´pRpρ ąm q `p1 ´pqRpρ ďm q › › ›p Rďm pρ ąm q `p1 ´pq Rąm pρ ďm q ď p KLpρ ąm , π ąm q m `p1 ´pq KLpρ ďm , π ďm q n ´m `ln

4 ? mpn´mq δ n ,
with Rpρ ąm q " ş Θąm Rpf θ qρpdθq, and Rpρ ďm q " ş Θďm Rpf θ qρpdθq, and Rďm pρ ąm q" ş Θąm Rpf θ qρpdθq, and Rďm pρ ąm q" ş Θďm Rpf θ qρpdθq. The result follows by kl's convexity. The complete proof is reported in App. B. In practice, following [START_REF] Mhammedi | PAC-Bayes un-expected bernstein inequality[END_REF] we set m " n 2 and p " m n , and we learn the base classifiers by empirical risk minimization.

Comparison with existing bounds Until now, we considered bounds for the expected risk over the space of MVs but not for a single realization f θ , unlike state-of-the-art methods. Under Dirichlet assumptions, we can bound the risk of the expected MV f θ by twice the expected risk (see App A.5): Rpf θ q " RpE ρpθq f θ q ď 2E ρpθq Rpf θ q. The obtained oracle bound would comprise an irreducible factor, but its KL term would not degrade, unlike state-of-the-art bounds that account for voter correlation [START_REF] Lacasse | PAC-Bayes bounds for the risk of the majority vote and the variance of the gibbs classifier[END_REF][START_REF] Lacasse | Learning with randomized majority votes[END_REF][START_REF] Andrés | Second order PAC-Bayesian bounds for the weighted majority vote[END_REF]. The empirical bound on the deterministic MV could be also refined leveraging works on the disintegration of PAC-Bayesian guarantees [START_REF] Blanchard | Occam's Hammer[END_REF][START_REF] Catoni | PAC-Bayesian Supervised Classification: The Thermodynamics of Statistical Learning[END_REF], Rivasplata et al., 2020[START_REF] Viallard | A general framework for the derandomization of PAC-Bayesian bounds[END_REF].

A downside of our method is that the complexity of the KL term grows with the number of base classifiers M , unlike the KL on categoricals that tends to 0 in its limit. This results in making the generalization bounds increasingly looser and conservative with growing M , even for low empirical risks. From a generalization perspective, our guarantees are hence able to reflect the complexity of the model, expressed as the size of the hypothesis space M . However, our risk certificates do not account for redundancy in the voter set. For instance, they are not able to distinguish scenarios where base classifiers are highly correlated (hence less complex hypothesis space) from scenarios where base classifiers are independent (more complex). An expedient for ensuring that generalization bounds are tight consists in learning the hypothesis space: the number of base predictors can then be limited without degrading the performance of the majority vote.

Experiments

In this section, we empirically evaluate StocMV, and we compare its generalization bounds and test errors to those obtained with PAC-Bayesian methods learning majority votes. We show that our method allows to derive generalization bounds that are consistently tight (i.e. close to the test errors) and non-vacuous (i.e. smaller than 1) both when studying ensembles of weak predictors and when studying ensembles of strong ones.

In the following, we consider Dirichlet distributions for the prior and the posterior of our method, and refer to the model obtained by optimizing the exact risk as exact and the one obtained by optimizing the approximated one as MC. We consider as baselines the PAC-Bayesian methods described in Section 2.1: We refer to the methods optimizing the First Order, Second Order and Randomized empirical bound as FO, SO and Rnd respectively. We do not compare with the C-Bound, as it is hard to optimize on large-scale datasets and existing algorithms are suited only for binary classification. All generalization bounds are evaluated with a probability 1´δ"0.95 and all prior distributions are set to the uniform (we provide a study for different priors in App. C.3). The posterior parameters (α for our method, θ for the others) are initialized uniformly in r0.01, 2s (and normalized to sum to 1 for SO, FO and Rnd ). Finally, for MC the sigmoid's slope parameter c is set to 100 and for Rnd the number of voters drawn at each iteration is set to N "100. Code was implemented in pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] and all experiments were run on a virtual machine with 8 vCPUs and 128Gb of RAM.

Comparison of exact and MC variants

For this set of experiments, we optimize Seeger's Bound (Equation ( 1)) by (batch) Gradient Descent, for 1, 000 iterations and with learning rate equal to 0.1. We study the performance of our method on the binary classification two-moons dataset, with 2 features, 2 classes and N p0, 0.05q Gaussian noise. We draw 1, 000 points for training, and 1, 000 points for testing. Figure 4 reports a comparison of exact and MC variants in terms of error, generalization bound and training time (in seconds). Increasing the number of MC draws T unsurprisingly allows to recover exact's performance, but at a higher computational cost for M ą 10 3 in this case. At the same time, when M ą 102 , MC achieves exact's errors even with T " 1, probably because of increased redundancy in the base predictor set. We also observe that the error rates and bounds gradually degrade for higher values of M for both methods. This is due to the KL term increasing with M , as highlighted in Section 4, becoming a too strong regularization during training and making the bound looser.

Experiments on real benchmarks

We now compare the considered methods on real datasets and on two different scenarios, depending on the type of PAC-Bayesian bounds that are evaluated: When making use of data-independent priors, we chose as voters axis-aligned decision stumps, with thresholds evenly spread over the input space (10 per feature); When making use of data-dependent priors, we build Random Forests [START_REF] Breiman | Random forests[END_REF] as set of voters, with M "100 trees learned bagging n 2 points and sampling ? d random features to ensure voter diversity, optimizing Gini impurity score and, unless stated otherwise, without bounding their maximal depth.

We consider several classification datasets from UCI [START_REF] Dua | UCI machine learning repository[END_REF], LIBSVM 2 and Zalando [START_REF] Xiao | Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms[END_REF], of different number of features and of instances. Their descriptions and details on any preprocessing are provided in App. C.3. We train the models by Stochastic Gradient Descent (SGD) using Adam [Kingma and Ba, 2015] with p0.9, 0.999q running average coefficients, batch size equal to 1024 and learning rate equal to 0.1 with a scheduler reducing this parameter of a factor of 10 with 2 epochs patience. We fix the maximal number of epochs to 100 and patience equal to 25 for early stopping, and for MC we fix T " 10 to increase randomness.

We report the test errors and generalization bounds in Figure 5 (additional results are reported in the appendix, in Tables 1 and2 and Figure 15): We compare the different methods on binary datasets and with data-independent priors in Figure 5a, and on multi-class datasets and with data-dependent priors in Figure 5b. First we notice that the bounds obtained by our method are consistently non vacuous and tighter than those obtained by the baselines on all datasets. Regarding the error rates, our method's performance is generally aligned with the baselines, while it achieved error rates significantly lower than FO and SO on the perfectly separable two-moons dataset. Sensitivity to noise could explain why our method does not outperform the baselines on the studied real problems, as these usually present label and input noise. Indeed our learning algorithm optimizes the 01-loss, which does not distinguish points with margins close or far from 0.5 because of its discontinuity in W θ " 0.5. Preliminary results reported in App C.2 seem to confirm this supposition. Finally, to gain a better understanding of the relation between base classifier strength and performance of the models obtained with the different methods, we further study their performance and generalization guarantees with varying voter strength. As hypothesis set, we learn a Random Forest with 100 decision trees, as before, but for this experiment we bound their maximal depth between 1 and 10. Constraining the tree depth allows to indirectly control how well the voters fit the training set (as shown in Appendix, Figure 14). Following [START_REF] Sloth Lorenzen | On PAC-Bayesian bounds for random forests[END_REF], we assess the voter strength by computing the expected accuracy of a random voter. All methods' results improve overall when voters get stronger, even though FO at a slower pace. Notice that on the considered datasets SO is the most sensitive method, particularly suffering from weak base predictors. Our method generally provides test errors comparable with the best baselines and consistently tighter bounds.

Future work

We propose a stochastic version of the classical majority vote classifier, and we directly analyze and optimize its expected risk through the PAC-Bayesian framework. The benefits on the model accuracy of this direct optimization are however reduced in presence of input noise, and fostering robustness in noisy contexts is the subject of future work. Another potential improvement would consist in tackling the discussed looseness of our generalization bounds with increasing number of base predictors, by accounting for redundancy in the hypothesis space.

A Derivation details under Dirichlet assumptions

A.1 Dirichlet distribution

The probability density function of the Dirichlet distribution is defined as follows:

θ " Dpα 1 , . . . , α M q, ρpθq " 1 Bpαq M ź j"1 pθ j q αj ´1, (12)
with α " rα j P R `sM j"1 the vector of the distribution parameters and Bpαq the normalized multi-variate beta function, defined using the gamma function Γ:

Bpαq " ś M j"1 Γpα j q Γp ř M j"1 α j q , Γpaq " ż 8 0 t a´1 e ´tdt. (13)
A.2 Aggregation property of Dirichlet distributions Lemma 2. If θ " Dpα 1 , . . . , α M q, then for any j P r1, M s and j 1 P r1, M sztju the variable θ 1 formed by dropping θ j and θ j 1 and adding their sum also follows a Dirichlet distribution pθ 1 , . . . , θ M , θ j `θj 1 q " Dpα 1 , . . . , α M , α j `αj 1 q.

A.3 KL divergence between Dirichlet distributions

Let ρpθq " Dpαq and πpθq " Dpβq, with α 0 " ř M j"1 α j and β 0 " ř M j"1 β j . The KL divergence between ρ and π is equal to:

KLpρ, πq " lnpΓpα 0 qq ´M ÿ j"1 lnpΓpα j qq ´lnpΓpβ 0 qq `M ÿ j"1 lnpΓpβ j qq `M ÿ j"1 pα j ´βj qpψpα j q ´ψpα 0 qq (14)
with ψpaq the digamma function: the first derivative of lnpΓpaqq (see Figure 7).

Proof. We have

KLpρ, πq " ż Θ ρpθq ln ˆρpθq πpθq ˙dθ (15) " ż Θ ln ˜Bpβq Bpαq ś M j"1 θ αj ´1 j ś M j"1 θ βj ´1 j ¸ρpθqdθ (16) " ż Θ ˜ln Bpβq ´ln Bpαq `ln ˜M ź j"1 θ αj ´βj j ¸¸ρpθqdθ (17) " ż Θ ˜ln Bpβq ´ln Bpαq `M ÿ j"1 pα j ´βj q ln θ j ¸ρpθqdθ (18)
" ln Bpβq ´ln Bpαq `M ÿ j"1 Equation ( 20) follows by definition of Dirichlet's geometric mean: ż Θ ln θ j ρpθqdθ " ψpα j q ´ψpα 0 q.

pα j ´βj q ż Θ ln θ j ρpθqdθ (19)

A.4 Partial derivatives for gradient-based optimization

When optimizing our objective functions by gradient descent, we make use of the following partial derivatives w.r.t. the Dirichlet parameters tα j u M j"1 . The partial derivatives of the KL divergence are ∇ αj KLpρ, πq " ψpα 0 q ´ψpα j q `ψpα j q ´ψpα 0 q `pα j ´βj qpψ 1 pα j q ´ψ1 pα 0 qq (21)

" pα j ´βj qpψ 1 pα j q ´ψ1 pα 0 qq.

For a given data point px, yq " P, recall that w " tj|h j pxq ‰ yu is the set of indices of the base classifiers that misclassify px, yq and c " tj|h j pxq " yu is the set of indices of the base classifiers that correctly classify px, yq. Let us define wpαq " ř jPw α j and cpαq " ř jPc α j " α 0 ´wpαq. For any j P c we have ∇ αj I 0.5 pcpαq, wpαqq " pln 0.5 ´ψpcpαqq `ψpα 0 qqI 0.5 pcpαq, wpαqq ´0.5 cpαq ΓpcpαqqΓpα 0 q Γpwpαqq 3 F 2 pcpαq, cpαq, 1 ´wpαq; cpαq `1, cpαq `1; 0.5q

and for any j P w we have ∇ αj I 0.5 pcpαq, wpαqq " p´ln 0.5 `ψpwpαqq ´ψpα 0 qqI 0.5 pcpαq, wpαqq `0.5 wpαq ΓpwpαqqΓpα 0 q Γpcpαqq 3 F 2 pwpαq, wpαq, 1 ´cpαq; wpαq `1, wpαq `1; 0.5q

where 3 F 2 pa, b, c; d, e; zq is the generalized hyper-geometric function:

3 F 2 pa, b, c; d, e; zq " T ÿ t"1 paq t pbq t pcq t z t pdq t peq t t! with p.q t the rising factorial.

The hyper-geometric function can be slow to evaluate, as its convergence rate varies depending on cpαq and wpαq. A possible strategy for speeding up this evaluation, apart from parallelizing computations, would be dynamic programming, i.e. storing the gradients of the incomplete beta function, as it is likely to be evaluated several times for the same cpαq and wpαq.

A.5 Oracle bound on expected Majority Vote

We now derive an oracle bound for the risk of the expected majority vote, given a distribution ρpθq having a Dirichlet form and concentration vector α " tα j u M j"1 . The oracle bound can then be used to derive empirical PAC-Bayesian bounds on the true risk of the expected MV. The expected MV is parameterized by the mean weighting vector θ, which for a Dirichlet distribution is given by

θ " E ρ θ " α α 0 with α 0 " M ÿ j"1 α j .
For a given px, yq " P, let us define wpx, yq " tj|h j pxq ‰ yu the set of indices of the base classifiers that misclassify px, yq and cpx, yq " tj|h j pxq " yu be the set of indices of the base classifiers that correctly classify px, yq. The risk of the expected MV can then be measured as follows:

Rpf θ q " P " W θ ě 0.5 ‰ ( 25
) " P rE W θ ě 0.5s (
Equation ( 26) follows by linearity of W θ . We recall that the expected risk is given by E px,yq"P I 0.5 ¨ÿ jPcpx,yq α j , ÿ jPwpx,yq α j '" E px,yq"P I 0.5 `α0 ´αwpx,yq , α wpx,yq with α wpx,yq " ř jPwpx,yq α j . As I 0.5 pq is monotonically increasing in its second argument and I 0.5 `α0 2 , α0 2 ˘" 0.5, we can then relate the two risks:

Rpf θ q ď 2E ρ Rpf θ q. ( 29
)
See Figure 8 for an illustration.

B Proof of Theorem 2

Let us define the following empirical risks

Rďm pρ ąm q " ż Θąm Rďm pf θ qρ ąm pdθq and Rąm pρ ďm q " ż Θďm Rąm pf θ qρ ďm pdθq and true risks Rpρ ąm q " ż Θąm Rpf θ qρ ąm pdθq and Rpρ ďm q " ż Θďm Rpf θ qρ ďm pdθq.

Theorem 2 (Seeger's bound with informed priors). Let π ďm and ρ ďm be the prior and posterior distributions on Θ ďm and π ąm and ρ ąm the prior and posterior distributions on Θ ąm . For any p P p0, 1q and δ P p0, 1q with probability at least 1´δ over samples S " tpx i , y i q"Pu n i"1 of size n we have simultaneously:

kl ´pRpρ ąm q `p1 ´pqRpρ ďm q › › ›p Rďm pρ ąm q `p1 ´pq Rąm pρ ďm q ď p KLpρ ąm , π ąm q m `p1 ´pq KLpρ ďm , π ďm q n ´m `ln 4 ? mpn´mq δ n with Rpρ ąm q " ş Θąm Rpf θ qρpdθq, and Rpρ ďm q " ş Θďm Rpf θ qρpdθq, and Rďm pρ ąm q " ş Θąm Rpf θ qρpdθq, and Rďm pρ ąm q " ş Θďm Rpf θ qρpdθq.

Proof. From the joint convexity of the binary kl divergence, we have for any p P r0, 1s kl ´pRpρ ąm q `p1 ´pqRpρ ďm q › › › p Rďm pρ ąm q `p1 ´pq Rąm pρ ďm q

ď pkl ´Rpρ ąm q › › › Rďm pρ ąm q ¯`p1 ´pqkl ´Rpρ ďm q › › › Rąm pρ ďm q ¯.
For @ρ ąm defined on Θ ąm :

P Sďm"P m « p kl ´Rpρ ąm q › › › Rďm pρ ąm q ¯ď p KLpρ ąm , π ąm q `ln 4 ? m δ m ff ě 1´δ 2 ñ P S"P n « p kl ´Rpρ ąm q › › › Rďm pρ ąm q ¯ď p KLpρ ąm , π ąm q `ln 4 ? m δ m ff ě 1´δ 2 ; (30)
and for @ρ ďm defined on Θ ďm :

P Sąm"P n´m « p1 ´pqkl ´Rpρ ďm q › › › Rąm pρ ďm q ¯ď p1 ´pq
KLpρ ďm , π ďm q `ln 4 ? n´m δ n ´m ff ě 1´δ 2 ñ P S"P n « p1 ´pq kl ´Rpρ ďm q › › › Rąm pρ ďm q ¯ď p1 ´pq KLpρ ďm , π ďm q `ln 4

? n´m δ n ´m ff ě 1´δ 2 . ( 31
)
Combining Equation ( 30) and Equation (31) using the union bound, we obtain the desired result with 1 ´δ probability.

C Additional experimental results

C.1 Analysis of the role of bound regularization on the posterior

Synthetic dataset We study the behavior of the method during training on a simple toyset, built from two normal distributions N 1 pr´1, 0s, diagpr0.1, 1sqq and N 2 pr1, 0s, diagpr0.1, 1sqq, one per class, so that the two classes are almost perfectly separable on the first dimension. from each class distribution for training, n points are drawn, and 1, 000 points for testing.

Base predictors and prior

We fix the set of base voters to M "4 axis-aligned decision stumps: 2 per class and centered in 0 on each dimension, so that the problem is well specified as the optimal classifier is in the predictor set. Then, we set the parameters β of the prior distribution all to 0.1 to encourage sparse solutions for the posterior, and we initialize the posterior's parameters by drawing α j " Up0.01, 2q. Study of the impact of optimizing the PAC-Bayesian bound In Figure 9, we compare the evolution of α during training for the exact method and n"50. In each plot we compare the posterior parameters α obtained when optimizing Bound (Equation ( 1)) with the one obtained when optimizing only the empirical risk (hence, without any regularization). Namely, we study the evolution during training of }α} 2 , of α's gradients and of the KL divergence w.r.t. prior: KLpρ, πq. We notice that without regularization, α diverges, while it tends to a constant and smaller value when optimizing the bound. This behaviour is not due to optimization instability, as shown by the gradient smoothly tending to 0 for both methods. Instead, it can be explained considering that without regularization the stochastic MV tends to concentrate around the optimal MV, i.e. a single θ. This behavior results in reducing the variance of the distribution, thus increasing the concentration parameters α. In contrast, when optimizing the bound the model fits the empirical risk and finds optimal solutions that are also close to the prior.

C.2 Results on two-moons dataset

Comparison of exact and MC variants In Figure 10 we report a second comparison of exact and MC variants in terms of error, generalization bound and training time (in seconds). This time we fix the number of base classifiers M and vary the size of the training set n. We first remark that T " 100 samples for MC are consistently enough to recover exact's performance, while allowing for significantly faster optimization. In general, as MC is easily parallelizable, its training times have a better dependence on n.

Full comparison with baselines

We report the full comparison of FO, SO, Rnd and our method on the binary classification two-moons dataset. We report in Figure 11 the results obtained after optimizing the respective PAC-Bayesian bounds. Notice that our bounds are tighter than the baselines, for any value of M : this is principally due to the fact that our method consistently obtains lower error rates, being able to Figure 12: Comparison in terms of error rates and PAC-Bayesian bounds, depending on the magnitude of input noise N p0, σ 2 q. We report the means and standard deviations over 10 different runs.

better fit the training set (as shown in Figure 3), but also to the fact that our bound is not based on an oracle upper bound on the true risk, unlike the baselines.

Sensitivity to noise

We observe that on two-moons our method outperforms the baselines also in terms of error rates. However, this does not seem to be the case on real benchmarks where its test errors are generally aligned with those of the baselines. We conduct an additional experiment to study whether this phenomenon is due to sensitivity to noise, such as input noise. Indeed our learning algorithm optimizes the 01-loss, which does not distinguish points with margins close or far from 0.5 because of its discontinuity in W θ " 0.5. In Figure 12, we assess training and test errors, and Seeger's bound values with increasing input noise. For this experiment, we generated training and test sets as before, and with Gaussian noise N p0, σ 2 q added to the inputs. As expected, all methods degrade with increasing noise. In particular, the test errors of exact and Rnd worsen the fastest and the benefits of using them vanish from σ 2 ą 0.35.

C.3 Additional results on real benchmarks

Dataset descriptions We consider several classification datasets from UCI [START_REF] Dua | UCI machine learning repository[END_REF], LIB-SVM https://www.csie.ntu.edu.tw/~cjlin/libsvm/ and Zalando [START_REF] Xiao | Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms[END_REF], of different number of features and of points:

• Haberman (UCI): prediction of survival of n " 306 patients who had undergone surgery from d " 3 anonymized features;

• TicTacToe (UCI): determination of a win for player x at TicTacToe game of any of the n " 958 board configurations (d " 9 categorical states);

• Svmguide1 (LIBSVM): d " 4 features, n " 7, 089 instances and 2 classes (no description available);

• Mushrooms (UCI): prediction of edibility of n " 8, 124 mushroom sample, given their d " 22 categorical features describing their aspect;

• Phishing (LIBSVM): prediction of phishing websites (n " 2456 websites and d " 68 binary encoded features);

• Adult (LIBSVM a1a): determining whether a person earns more than 50K a year (n " 32, 561 people and d " 123 binary features);

• CodRNA (LIBSVM): detection of non-coding RNAs among n " 59, 535 instances and from d " features;

• Pendigits (UCI): recognition of hand-written digits (10 classes, d " 9 features and n " 12, 992);

• Protein (LIBSVM): d " 357 features, n " 24, 387 instances and 3 classes;

• Shuttle (UCI): d " 9 features, n " 58, 000 and 7 classes;

• Sensorless (LIBSVM): prediction of motor condition (n " 58, 509 instances and 11 classes), with intact and defective components, from d " 48 features extracted from electric current drive signals;

• MNIST (LIBSVM): prediction of hand-written digits (n " 70, 000 instances and 10 classes) from d " 28 ˆ28 gray-scale images;

• Fashion-MNIST (Zalando): prediction of cloth articles (n " 70, 000 instances and 10 classes) from d " 28 ˆ28 gray-scale images.

At each run of an algorithm, we randomly split a dataset in training and test sets of sizes 80% ´20% respectively, apart from Svmguide1, Pendigits, Protein, Shuttle, MNIST and Fashion-MNIST for which we keep the original train-test splits. Note that we do not make use of a validation set, as we use the PAC-Bayesian bounds as estimate of the test error for model selection. Finally, we convert all categorical features to numerical using an ordinal encoder and z-score all features using the statistics of the training set.

Choice of prior In all the other experiments, we fixed the prior distribution (parameterized by rβ j s M j"1 ) to the uniform, i.e. β j " 1, @j. This choice was to make the comparison with the baselines as fair as possible, as their prior was also fixed to the uniform (categorical). However, we can bias the sparsity of the posterior, or conversely its concentration, by choosing a different value for the prior distribution parameters. In some cases, tuning the prior parameters allows to obtain better performance, as reported in Figure 13. In particular, on Protein encouraging sparser solutions generally provides better results, confirmed by the fact that the best baseline on this dataset, FO, is known to output sparse solutions. On the contrary, on datasets where methods accounting for voter correlation outperform FO, such as on MNIST, encouraging solutions to be concentrated and close to the simplex mean yields better performance. In general, these results suggest that the choice of prior distribution has a high impact on the learned model's performance and tuning its concentration parameters would be a viable option for improving results.

C.4 Impact of voter strength

We report the complete study on the impact of voter strength on the learned models. More precisely we provide results for additional datasets as well as the study of the expected strength of a voter as a function of the tree maximal depth. Recall that as hypothesis set, we learn a Random Forest with 100 decision trees for which we bound the maximal depth between 1 and 10. In Figure 14, we can see that limiting the maximal depth is an effective way for controlling the strength of the voters, measured as the expected accuracy of a random voter. Apart from Protein, where decision trees do not seem to be a good choice of base predictor, increasing the strength of the voters generally yields more powerful ensembles for all methods. Our method has error rates comparable with the best baselines and enjoys tight and non-vacuous generalization guarantees for any tree depth. Finally, by comparing SO's training and test errors we notice that this method tends to overfit the dataset especially when the base classifiers are weaker (tree depth close to 1). We fix all M prior's parameters to β, on the x axis: the smaller β, the sparser the posterior is encouraged to be. We plot average and standard deviations over 4 trials.

C.5 Summary of main results

We finally report the detailed comparison on real benchmarks in Figure 15 and Tables 1 and2.

Figure 1 :

 1 Figure 1: Oracle upper bounds for the true risk. The risks are the areas under the respective curves, for an arbitrary distribution of W θ (typically different from the uniform).

Figure 2 :

 2 Figure 2: Visualization of the density measure ρ : R 3 Ñ r0, 1s taking the form of a Dirichlet distribution, with concentration parameters α. The darker the color, the higher ρpθq. Each θ on the simplex corresponds to a majority vote classifier f θ and has an associated probability depending on α.

Figure 4 :

 4 Figure4: Average performance for 10 trials of exact and MC variants of our method as a function of the number of decision stumps M (represented in log scale), for n " 1, 000 and different number of MC draws T . The same analysis but for fixed M and increasing n is reported in Appendix C.2.

  Figure 5: Comparison in terms of test error rates and PAC-Bayesian bound values. We report the means (bars) and standard deviations (vertical, magenta lines) over 10 different runs.

Figure 6 :

 6 Figure6: Test error rates and PAC-Bayesian bound values as a function of the maximal depth of the voters (decision trees here). We report the means and standard deviations over 4 different runs, and mark with a red horizontal line the threshold above which Seeger's bounds are vacuous. Additional results are available in the appendix.

Figure 8 :

 8 Figure 7: Functions of the gamma family.

Figure 9 :Figure 10 :

 910 Figure9: Evolution of }α} 2 (left), its gradient (center) and KLpα, βq (right) during training using the exact method. In each plot we compare the values obtained when optimizing Seeger's Bound (Equation (1)) with those obtained when optimizing only the empirical risk.

Figure 11 :

 11 Figure 11: Comparison in terms of test error rates and PAC-Bayesian bound values. Each block corresponds to a different number of predictors M . We report the means (bars) and standard deviations (vertical, magenta lines) over 10 different runs.

Figure 13 :

 13 Figure13: Study of impact of prior on posterior's performance. We fix all M prior's parameters to β, on the x axis: the smaller β, the sparser the posterior is encouraged to be. We plot average and standard deviations over 4 trials.

Figure 14 :

 14 Figure 14: Comparison voter strength (1st column), training error (2nd column), test error (3rd column) and Seeger's bound (4th column) as a function of the tree maximal depth. We mark with a red horizontal line the threshold above which the bounds are vacuous. Results are averaged over 4 trials.

Figure 15 :

 15 Figure 15: Comparison of First Order (FO), Second Order (SO) and Randomized MV (Rnd ) and our methods (exact, MC ) in terms of training and test error rates and PAC-Bayesian bound values. Each row of subfigures corresponds to a dataset, where we marked its number of features d and number of training instances n. The dashed horizontal line in the rightmost column plots marks the threshold above which the bounds are vacuous. We report the means (bars) and standard deviations (vertical, black lines) over 10 different runs.

Table 1 :

 1 Summary of the main results on binary datasets with data agnostic prior. For the PAC-Bayesian methods, we report the performance of the model obtained after optimizing Seeger's Bound (Theorem (1)). Bayesian Naive Bayes (BayesianNB) corresponds to the weighting strategy proposed in[START_REF] Berend | A finite sample analysis of the naive Bayes classifier[END_REF]. For each method we report the average (˘standard deviations) of the test errors and generalization bounds (for the PAC-Bayesian methods) over 10 runs. We bold the results (test error and bound value) that are significantly better (smaller) than the other results for a dataset. Notice that the PAC-Bayesian methods generally improve upon BayesianNB and that our method consistently provides the tightest and non-vacuous generalization bounds over all the datasets.

	PHISHING CODRNA ADULT	44.49 ˘0.05 66.67 ˘0.00 24.07 ˘0.00	11.26 ˘0.58 23.52 ˘0.35 22.14 ˘0.36	10.10 ˘0.73 22.93 ˘1.43 17.02 ˘0.29	6.84 ˘0.41 12.41 ˘0.23 16.24 ˘0.37	8.12 ˘0.35 11.77 ˘0.16 23.49 ˘1.99	7.74 ˘0.42 13.13 ˘3.43 24.07 ˘0.00	25.75 ˘0.32 49.45 ˘0.47 46.93 ˘0.19	35.67 ˘0.21 60.59 ˘0.20 53.33 ˘0.09	24.64 ˘0.20 36.24 ˘0.13 41.06 ˘0.13	13.50 ˘0.19 15.16 ˘0.17 25.42 ˘1.07	14.62 ˘0.28 16.48 ˘2.86 26.46 ˘0.09
	HABERMAN TICTACTOE SVMGUIDE MUSHROOMS	25.65 ˘2.10 30.21 ˘3.06 34.51 ˘0.00 11.65 ˘0.59	25.48 ˘3.37 29.22 ˘2.77 6.83 ˘0.00 4.60 ˘0.36	26.45 ˘2.99 30.52 ˘3.55 6.83 ˘0.00 4.60 ˘0.36	26.45 ˘2.19 27.76 ˘3.71 6.04 ˘0.32 1.08 ˘0.22	29.84 ˘7.34 30.91 ˘2.10 6.11 ˘0.25 1.30 ˘0.35	30.10 ˘7.21 30.84 ˘2.38 6.27 ˘0.40 1.23 ˘0.23	75.24 ˘1.78 76.48 ˘1.48 18.62 ˘0.00 12.31 ˘0.22	110.60 ˘1.65 100.25 ˘0.99 31.06 ˘0.36 21.27 ˘0.31	89.19 ˘10.59 79.55 ˘0.90 26.37 ˘1.34 11.49 ˘0.27	46.97 ˘17.69 42.49 ˘0.53 11.38 ˘0.61 4.77 ˘0.26	47.34 ˘17.57 43.07 ˘0.72 11.79 ˘0.68 5.42 ˘0.20
	Method	BayesianNB	Fist Order	Second Order	RandomizedMV	ours-exact	ours-MC	Fist Order	Second Order	RandomizedMV	ours-exact	ours-MC
			Test error		Seeger's b.

Table 2 :

 2 Summary of the main results on multi-class datasets with informed priors. For the PAC-Bayesian methods, we report the performance of the model obtained after optimizing Seeger's Bound with informed priors (Theorem (2)). Bayesian Naive Bayes (BayesianNB) corresponds to the weighting strategy proposed in[START_REF] Berend | A finite sample analysis of the naive Bayes classifier[END_REF]. For each method we report the average (˘standard deviations) of the test errors and generalization bounds (for the PAC-Bayesian methods) over 10 runs. We bold the results (test error and bound value) that are significantly better (smaller) than the other results for a dataset. Notice that the PAC-Bayesian methods generally improve upon BayesianNB and that our method consistently provides the tightest and non-vacuous generalization bounds over all the datasets.

	FASHION-MNIST	89.14 ˘0.05	25.94 ˘0.22	21.67 ˘0.14	21.80 ˘0.07	21.90 ˘0.05	21.90 ˘0.05	52.12 ˘0.16	57.09 ˘0.08	42.93 ˘0.12	21.83 ˘0.06	21.83 ˘0.06
	MNIST	86.09 ˘0.09	21.45 ˘0.17	13.50 ˘0.08	13.62 ˘0.11	13.84 ˘0.08	13.84 ˘0.09	46.18 ˘0.37	44.24 ˘0.12	30.46 ˘0.14	15.57 ˘0.07	15.49 ˘0.17
	PENDIGITS PROTEIN SHUTTLE SENSORLESS	82.41 ˘0.11 82.09 ˘0.00 20.90 ˘0.03 86.27 ˘0.21	13.01 ˘1.12 54.51 ˘0.69 0.08 ˘0.01 2.94 ˘0.12	7.59 ˘0.28 64.15 ˘0.39 0.05 ˘0.01 0.41 ˘0.06	7.86 ˘0.19 55.24 ˘0.47 0.07 ˘0.00 0.43 ˘0.05	8.04 ˘0.23 55.00 ˘0.78 0.07 ˘0.01 0.46 ˘0.05	8.05 ˘0.23 54.74 ˘0.80 0.07 ˘0.01 0.46 ˘0.05	20.65 ˘0.45 112.29 ˘0.46 0.32 ˘0.02 6.80 ˘0.26	16.38 ˘0.20 139.81 ˘0.23 0.42 ˘0.01 3.82 ˘0.06	8.32 ˘0.16 122.92 ˘0.41 0.28 ˘0.01 1.29 ˘0.05	4.34 ˘0.11 59.28 ˘0.60 0.15 ˘0.01 0.67 ˘0.03	4.35 ˘0.11 59.25 ˘0.63 0.15 ˘0.01 0.67 ˘0.03
	Method	BayesianNB	Fist Order	Second Order	RandomizedMV	ours-exact	ours-MC	Fist Order	Second Order	RandomizedMV	ours-exact	ours-MC
			Test error		Seeger's b.

Oracle bounds are expressed in terms of the unknown data distribution; their exact value cannot be computed.

1`expp´cxq with slope parameter c P R `. De facto this corresponds to solving a relaxation of the problem and not its exact form[Nesterov,

2005]. At each iteration of the MC optimization algorithm we perform:1. Draw a sample tθ t " ρpαqu T t"1 using the implicit reparameterization trick[Figurnov et al., 2018, Jankowiak and[START_REF] Jankowiak | Pathwise Derivatives Beyond the Reparameterization Trick[END_REF];

https://www.csie.ntu.edu.tw/~cjlin/libsvm/

Acknowledgement

This work was supported by the French Project APRIORI ANR-18-CE23-0015. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr). P. Germain is supported by the Canada CIFAR AI Chair Program, and the NSERC Discovery grant RGPIN-2020-07223.