Michel Duprez
email: mduprez@math.cnrs.fr

Stéphane Pierre

Alain Bordas

Marek Bucki
email: marek.bucki@texisense.com

Huu Phuoc Bui

Franz Chouly
email: franz.chouly@univ-fcomte.fr

Vanessa Lleras
email: vanessa.lleras@umontpellier.fr

Claudio Lobos
email: clobos@inf.utfsm.cl

Alexei Lozinski
email: alexei.lozinski@univ-fcomte.fr

Pierre-Yves Rohan
email: pierre-yves.rohan@ensam.eu

Satyendra Tomar
email: tomar.sk@iitkalumni.org

Huu Phuoc Bui
email: huu-phuoc.bui@math.cnrs.fr

Quantifying discretization errors for soft-tissue simulation in computer assisted surgery: a preliminary study

Keywords: computer assisted surgery, computational biomechanics, goal oriented error estimates, adaptive finite elements

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

Introduction

Patient-specific finite element models of soft tissue and organs have received an increasing amount of interest in the last decades. Such finite element models are widely employed to investigate both, the underlying mechanisms that drive normal physiology of biological soft tissues [START_REF] Avril | Anisotropic and hyperelastic identification of in vitro human arteries from full-field optical measurements[END_REF][START_REF] Franquet | Identification of the in vivo elastic properties of common carotid arteries from MRI: A study on subjects with and without atherosclerosis[END_REF][START_REF] Gras | Hyper-elastic properties of the human sternocleidomastoideus muscle in tension[END_REF], and the mechanical factors that contribute to the onset and development of diseases such as pressure ulcers [START_REF] Loerakker | The effects of deformation, ischemia, and reperfusion on the development of muscle damage during prolonged loading[END_REF][START_REF] Shilo | Identification of capillary blood pressure levels at which capillary collapse is likely in a tissue subjected to large compressive and shear deformations[END_REF], atherosclerosis or aneurysms [START_REF] Floc'h | Vulnerable Atherosclerotic Plaque Elasticity Reconstruction Based on a Segmentation-Driven Optimization Procedure Using Strain Measurements: Theoretical Framework[END_REF][START_REF] Badel | Finite Element simulation of buckling-induced vein tortuosity and influence of the wall constitutive properties[END_REF][START_REF] Romo | In vitro analysis of localized aneurysm rupture[END_REF], or multilevel lumbar degenerative disc diseases [START_REF] Park | Effects of degenerated intervertebral discs on intersegmental rotations, intradiscal pressures, and facet joint forces of the whole lumbar spine[END_REF][START_REF] Schmidt | Effect of multilevel lumbar disc arthroplasty on spine kinematics and facet joint loads in flexion and extension: a finite element analysis[END_REF], to name a few. Finite element models are also valuable tools that contribute to the development of medical devices, see e.g. vascular stent-grafts [START_REF] Perrin | Patient-specific numerical simulation of stent-graft deployment: Validation on three clinical cases[END_REF], artificial facet systems for spinal fusion [START_REF] Goel | Anatomic Facet Replacement System (AFRS) Restoration of Lumbar Segment Mechanics to Intact: A Finite Element Study and In Vitro Cadaver Investigation[END_REF] or knee braces [START_REF] Pierrat | Evaluation of the mechanical efficiency of knee orthoses: A combined experimentalnumerical approach[END_REF]. They have the potential to improve prevention strategies [START_REF] Luboz | Personalized modeling for real-time pressure ulcer prevention in sitting posture[END_REF][START_REF] Trabelsi | Patient specific stress and rupture analysis of ascending thoracic aneurysms[END_REF], surgical planning [START_REF] Buchaillard | Simulations of the consequences of tongue surgery on tongue mobility: implications for speech production in post-surgery conditions[END_REF] and pedagogical simulators for medical training [START_REF] Bro-Nielsen | Real-time Volumetric Deformable Models for Surgery Simulation using Finite Elements and Condensation[END_REF][START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF][START_REF] Buttin | Biomechanical simulation of the fetal descent without imposed theoretical trajectory[END_REF][START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF] . In this context, one major issue is meshing, since the reliability of the predicted mechanical response arising from computer simulation heavily relies on the quality of the underlying finite element mesh: if some elements of the mesh are too distorted or if the mesh is too coarse in some regions, the numerical solution may deteriorate significantly [START_REF] Grätsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. Whilst generating meshes of complex shapes with tetrahedral elements is generally possible thanks to advanced meshing algorithms, e.g. [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF][START_REF] Remacle | A two-level multithreaded delaunay kernel[END_REF], low-order Lagrange tetrahedral elements are unsuitable for most biomechanical problems due to volumetric locking. To circumvent such issues, strain smoothing approaches were developed [START_REF] Ong | On stability, convergence and accuracy of bes-fem and bfs-fem for nearly incompressible elasticity[END_REF][START_REF] Lee | Strain smoothing for compressible and nearly-incompressible finite elasticity[END_REF], which have the drawback of leading to larger system matrix bandwidth, but the advantage of being easily parallelized on graphical processing units thanks to nodal integration. Hexahedral low-order Lagrange elements alleviate locking issues, but research on the automatic generation of hexahedral meshes is still ongoing [START_REF] Bommes | Quad-mesh generation and processing: A survey[END_REF][START_REF] Pellerin | Combining tetrahedra into hexahedra: a vertex based strategy[END_REF][START_REF] Ray | Hexahedral meshing: Mind the gap!Preprint[END_REF], spurred by a recrudescent surge in research on polyhedral mesh generation [START_REF] Hu | Centroidal voronoi tessellation based polycube construction for adaptive all-hexahedral mesh generation[END_REF][START_REF] Mouly | Minimizing finite-volume discretization errors on polyhedral meshes[END_REF][START_REF] Budninskiy | Optimal voronoi tessellations with hessian-based anisotropy[END_REF][START_REF] Chen | An interface-fitted mesh generator and virtual element methods for elliptic interface problems[END_REF][START_REF] Walton | Advances in co-volume mesh generation and mesh optimisation techniques[END_REF] and approximations such as the virtual finite element method [START_REF] Beirão Da Veiga | Virtual element method for general second-order elliptic problems on polygonal meshes[END_REF], hybrid high-order (HHO) methods [START_REF] Di Pietro | A hybrid high-order locking-free method for linear elasticity on general meshes[END_REF], polyhedral finite elements [START_REF] Talebi | Stress analysis of 3d complex geometries using the scaled boundary polyhedral finite elements[END_REF][START_REF] Kraus | Finite element formulations for 3d convex polyhedra in nonlinear continuum mechanics[END_REF][START_REF] Natarajan | A new one point quadrature rule over arbitrary star convex polygon/polyhedron[END_REF][START_REF] Francis | Linear smoothed polygonal and polyhedral finite elements[END_REF], smoothed finite elements [START_REF] Sellam | Smoothed polygonal finite element method for generalized elastic solids subjected to torsion[END_REF], scaled finite elements [START_REF] Talebi | Stress analysis of 3d complex geometries using the scaled boundary polyhedral finite elements[END_REF][START_REF] Natarajan | A scaled boundary finite element formulation over arbitrary faceted star convex polyhedra[END_REF], for various applications, including large strain hyper-elasticity [START_REF] Rajagopal | ???? Hyperelastic analysis based on a polygonal finite element method[END_REF][START_REF] Abbas | Hybrid High-Order methods for finite deformations of hyperelastic materials[END_REF] and optimization [START_REF] Chau | A polytreebased adaptive polygonal finite element method for multi-material topology optimization[END_REF][START_REF] Chau | A polytreebased adaptive polygonal finite element method for multi-material topology optimization[END_REF][START_REF] Nguyen-Hoang | A new polyhedral element for the analysis of hexahedral-dominant finite element models and its application to nonlinear solid mechanics problems[END_REF]. The patient-specific mesh has to be built from segmented medical images (CT, MRI, ultrasound), and has to conform to anatomical details with potentially complex topologies and geometries [START_REF] Telfer | Simplified versus geometrically accurate models of forefoot anatomy to predict plantar pressures: A finite element study[END_REF][START_REF] Bijar | Atlas-based automatic generation of subject-specific finite element tongue meshes[END_REF][START_REF] Al-Dirini | Development and Validation of a High Anatomical Fidelity FE Model for the Buttock and Thigh of a Seated Individual[END_REF], which led to the design of algorithms that aim to optimize the quality of the generated mesh by reducing the distortion of the elements [START_REF] Zhang | Patient-Specific Vascular NURBS Modeling for Isogeometric Analysis of Blood Flow[END_REF][START_REF] Laville | Parametric and subject-specific finite element modelling of the lower cervical spine. Influence of geometrical parameters on the motion patterns[END_REF][START_REF] Shang | Hexahedral mesh generation via constrained quadrilateralization[END_REF]. These algorithms may also have to satisfy a number of additional constraints such as minimizing human intervention (automation), preservation of certain important anatomical details or robustness with respect to data [START_REF] Couteau | The Mesh-Matching algorithm: an automatic 3D mesh generator for finite element structures[END_REF][START_REF] Alliez | Variational tetrahedral meshing[END_REF][START_REF] Ito | Octree-based reasonable-quality hexahedral mesh generation using a new set of refinement templates[END_REF]116,[START_REF] Lobos | Techniques for the generation of 3D Finite Element Meshes of human organs[END_REF]. In general the quality of a given mesh can be assessed through purely geometrical criteria, that allow in some way to quantify the distortion of the geometry of the elements and how far they are from their ideal shape [START_REF] Frey | Mesh generation, 2nd Edition[END_REF][START_REF] Bucki | Jacobian-based repair method for finite element meshes after registration[END_REF][START_REF] Burkhart | Finite element modeling mesh quality, energy balance and validation methods: A review with recommendations associated with the modeling of bone tissue[END_REF][START_REF] George | Geometric validity (positive jacobian) of high-order lagrange finite elements, theory and practical guidance[END_REF]. To circumvent or simplify the mesh generation issue, implicit/immersed boundary approaches have been proposed, where the mesh does not conform to the geometry, which is treated by implicit functions such as level sets and enriched finite element methods. This idea was proposed in [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] and later generalized in [START_REF] Moumnassi | Finite element analysis on implicitly defined domains: An accurate representation based on arbitrary parametric surfaces[END_REF][START_REF] Moumnassi | Analysis using higher-order xfem: implicit representation of geometrical features from a given parametric representation[END_REF], [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF][START_REF] Burman | CutFEM: discretizing geometry and partial differential equations[END_REF], [START_REF] Liu | Stress analysis of stl models by octree mesh and scaled boundary finite element method[END_REF][START_REF] Kalameh | Semi-implicit representation of sharp features with level sets[END_REF] and [START_REF] Ródenas | Creation of patient specific finite element models of bone-prosthesissimulation of the effect of future implants[END_REF][START_REF] Tur | Stabilized method of imposing dirichlet boundary conditions using a recovered stress field[END_REF][START_REF] Fries | Higher-order meshing of implicit geometries. Part I: Integration and interpolation in cut elements[END_REF], in combination with (goal-oriented) error estimates. Although promising, applications of such approaches to patient-specific geometries remain in their infancy, see a review of related methods in [START_REF] Bordas | Recent advances towards reducing the meshing and re-meshing burden in computational sciences[END_REF]. Yet another approach consists in directly using the image as a model which could enable simulations without any mesh generation [START_REF] Li | Biomechanical model for computing deformations for whole-body image registration: A meshless approach[END_REF]. In "image as a model", the boundary of the domain is smeared, which significantly complicates imposition of boundary conditions, particularly contact. Finally, meshless methods [START_REF] Nguyen | Meshless methods: a review and computer implementation aspects[END_REF] are possible alternatives which may simplify biomechanics simulations by relaxing the constraints posed on mesh quality and simplifying local mesh adaptation. Comparatively to Galerkin meshfree methods [START_REF] Lim | On the use of meshfree methods and a geometry based surgical cutting algorithm in multimodal medical simulations[END_REF][START_REF] De | The point collocation-based method of finite spheres (PCMFS) for real time surgery simulation[END_REF][START_REF] Lim | Real time simulation of nonlinear tissue response in virtual surgery using the point collocation-based method of finite spheres[END_REF], point collocation methods [START_REF] Banihani | A comparison of some model order reduction methods for fast simulation of soft tissue response using the point collocation-based method of finite spheres[END_REF], also known as generalized finite difference methods, are potentially competitive as they do not require any numerical integration. Quality control of such collocation methods remains an open problem, as well as conditioning of the system matrix, which is strongly related to the choice of the stencil [START_REF] Davydov | Optimal stencils in Sobolev spaces[END_REF]. On the other hand, collocation methods are easily parallelized, for instance on graphical processing units. Beyond mesh quality, mesh density is another, related, parameter which must be controlled during biomechanics simulations. Solutions must be obtained on commodity hardware within clinical time scales: milliseconds (for surgical training); minutes (for surgical assistance); hours (for surgical planning). Therefore, and although this would lead to the most accurate solution, it is impractical to use a uniformly fine mesh over the whole domain. This remark begs the question: "given a tolerable error level, what is the coarsest possible mesh which will provide the required accuracy." This leads to the notion of "mesh optimality," which is achieved for an optimal balance between the accuracy in a given quantity of interest to the user and the associated computational cost. It is probably intuitively understood that this "optimality" criterion, and the resulting optimized mesh both depend on the quantity of interest and that, in general, the optimal mesh will be non-uniform, displaying local refinement around specific regions. A possible criterion for mesh adaptation can be any a priori knowledge of the problem or its solution such as geometry, material properties or boundary layers e.g. localized loads, contacts, sharp features, material interfaces. Similarly, knowledge of the quantity of interest can help guide local mesh refinement. Nevertheless, such mesh refinement guidelines are generally ad hoc and cannot guarantee the resulting mesh will be optimal. To summarize, the, standard, piecewise linear or quadratic Lagrange-based, finite element method is still the most popular technology to solve biomechanics problems. The choice of an optimal mesh, in particular its local refinement level for given problems and quantities of interest remains an open issue. Moreover, without knowing the finite element solution itself, it is practically impossible to quantify the adequacy of a given mesh only from heuristics or other ad hoc criteria derived from a priori knowledge of the problem or its exact solution.

As a result, we aim at addressing the following two questions in this paper:

1. For a patient-specific finite element computation, can we provide some information to the user about the accuracy of the numerical solution, namely can we compute an approximate discretization error caused by the choice of the mesh? By discretization error, we mean the difference between the finite element solution and the exact solution of the same boundary value problem on the same geometry. 2. Can the numerical solution be used to optimize the mesh in the critical regions only, to achieve maximum accuracy for a given computational cost, or, conversely, to achieve a given accuracy with a minimum computational cost?

For the sake of simplicity we do not consider 1. modeling errors, which arise due to the approximation of the geometry, physical assumptions, and uncertainty on material parameters, 2. numerical errors, which arise due to linearization, iterative solvers, and machine precision.

In this paper, we investigate the capability of a posteriori error estimates [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Verfürth | A posteriori error estimation techniques for finite element methods[END_REF] to provide useful information about the discretization error. A posteriori error estimates are quantities computed from the numerical solution, that indicate the magnitude of the local error. These estimates are at the core of mesh adaptive techniques [START_REF] Nochetto | Theory of adaptive finite element methods: an introduction[END_REF]. Many a posteriori error estimation methods have been developed in the numerical analysis community. These methods have different theoretical and practical properties. However, despite their great potential, error estimates, to the best of our knowledge, have rarely been considered for patient-specific finite element simulations in the biomechanical community. The only reference known to us which addresses discretization error estimation in biomechanics is the very recent paper [START_REF] Bui | Real-time error control for surgical simulation[END_REF] who consider simple but real-time error estimation approaches for needle insertion. We limit our study to a simplified setting in order to gain preliminary insights into the behaviour of such a posteriori error estimates and to address with the first technical difficulties. We focus on two-dimensional linear elasticity (plane strain) problems, with simple boundary conditions (prescribed displacements and tractions), and we assume triangular meshes. This is somehow restrictive in comparison to current practice in softtissue simulation. Among the existing a posteriori error estimates, we focus on Dual Weighted Residuals (DWR), as presented in, e.g., [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] (see also [START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF][START_REF] Paraschivoiu | A hierarchical duality approach to bounds for the outputs of partial differential equations[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF][START_REF] Maday | Numerical analysis of a posteriori finite element bounds for linear functional outputs[END_REF][START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF][START_REF] Bangerth | Adaptive finite element methods for differential equations[END_REF]). Indeed this method allows to estimate the error for a given quantity of interest. As a matter of fact, for the majority of applications, controlling the error in the energy norm is not relevant, and the error must be controlled for a specific quantity of interest to the user (e.g., stress intensity factors, shear stress or strain intensity at specific locations). The DWR method is conveniently implemented in the standard finite element library FEniCS [START_REF] Logg | Automated Solution of Differential Equations by the Finite Element Method[END_REF] and we make use of the implementation described in detail in the paper of Rognes and Logg [START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF], with minor modifications. This paper is organized as follows. In Section 2, we present the linear elastic problem, the corresponding finite element method, the a posteriori error estimates as well as the algorithm for mesh refinement. In Section 3, we consider two simplified test-cases, inspired by patient-specific biomechanics, where the current methodology is applied. The results are discussed in Section 4.

Material and methods

We first present the general problem considered in this contribution, that represents a simplified setting for contractile soft-tissue simulation. We then describe in details the computation of the a posteriori error estimate: a global estimator that provides an estimation of the discretization error and a local estimator that drives the mesh refinement.

We end this section with the description of a simple algorithm for mesh refinement. We first introduce some useful notations. In what follows, bold letters such as u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K) represent functional sets involving vector fields. As usual, we denote by (H

s (•)) d , s ∈ R, d = 1, 2, 3,
the Sobolev spaces in one, two or three space dimensions [START_REF] Adams | Sobolev spaces[END_REF]. In the sequel the symbol | • | will either denote the Euclidean norm in R d , or the measure of a domain in R d .

Setting: a "toy" boundary value problem in linear elasticity

We consider an elastic body whose reference configuration is represented by the domain Ω in R 2 . We consider the plane strain formulation, and allow only small deformations. We suppose that ∂Ω consists of two disjoint parts Γ D and Γ N , with meas(Γ D ) > 0. The unit outward normal vector on ∂Ω is denoted by n. A displacement u D = 0 is applied on Γ D , and the body is subjected to volume forces f ∈ (L 2 (Ω)) 2 and surface loads

F ∈ (L 2 (Γ N )) 2 .
We introduce the bilinear form

a(v, w) := Ω σ(v) : ε(w) dx,
which represents the (internal) virtual work associated to passive elastic properties. The notation ε(v) = 1 2 (∇v + ∇v T )/2 represents the linearized strain tensor field, and σ = (σ ij ), 1 ≤ i, j ≤ 2, stands for the stress tensor field, assumed to be given by Hooke's law. The linear form l E (w) :

= Ω f • w dx + Γ N F • w ds
stands for the virtual work of external loads in the body and on its surface. Finally we represent in a very simplified manner the active properties of soft-tissue as a linear anisotropic pre-stress

l A (w) := -βT ω A (ε(w)e A ) • e A dx,
where ω A is the part of the body where muscle fibers are supposed to act, T ≥ 0 is a scalar which stands for the tension of the fibers, e A is a field of unitary vectors that stands for muscle fibers orientation, and β ∈ [0, 1] is the activation parameter. When β = 0 there is no activation of the muscle fibers, and the value β = 1 corresponds to the maximum activation. This modeling can be viewed as a linearization of some more sophisticated active stress models of contractile tissues (see, e.g., [START_REF] Cowin | Cardiovascular soft tissue mechanics[END_REF][START_REF] Payan | Biomechanics of living organs : hyperelastic constitutive laws for finite element modeling[END_REF]). We want to solve the following weak problem

Find a displacement u ∈ V such that a(u, v) = l(v), ∀ v ∈ V, (1)
where

l(•) = l E (•) + l A (•)
, and where u and v lie in the space of admissible displacements

V := v ∈ H 1 (Ω) 2 | v = 0 on Γ D .
From the displacement field, we are interested in computing a linear quantity

J : V u → J(u) ∈ R, (2)
which can be defined according to a specific application and the interest of each practitioner. Thereby, the quantity J will be aptly called quantity of interest (QoI). We will provide its expression(s) for each test case.

Finite element method

Consider a family of meshes (K h ) h>0 constituted of triangles and assumed to be subordinated to the decomposition of the boundary ∂Ω into Γ D and Γ N . For a mesh K h , we denote by E h the set of edges, by E int h := {E ∈ E h : E ⊂ Ω} the set of interior edges, and by E N h := {E ∈ E h : E ⊂ Γ N } the set of boundary edges that correspond to Neumann conditions (we assume that any boundary edge is either inside Γ N or inside Γ D ). For an element

K of K h , we set E K the set of edges of K, E int K := E K ∩E int h and E N K := E K ∩E N h .
We also assume that each element K is either completely inside ω A or completely outside it. Let σ be a second-order tensorial field in Ω, which is assumed to be piecewise continuous. We define the jump of σ across an interior edge E of an element K, at a point y ∈ E, as follows

[[σ]] E,K (y) := lim α→0 + (σ(y + αn E,K ) -σ(y -αn E,K )) n E,K ,
where n E,K is the unit normal vector to E, pointing out of K. The finite element space V h ⊂ V is built upon continuous Lagrange finite elements of degree k = 1, 2 (see, e.g., [START_REF] Ern | Theory and practice of finite elements[END_REF]), i.e.

V h := v h ∈ (C 0 (Ω)) d : v h|K ∈ (P k (K)) d , ∀K ∈ K h , v h = 0 on Γ D .

Problem (1) is approximated by

Find u h ∈ V h such that a(u h , v h ) = l(v h ), ∀v h ∈ V h . (3)

Goal-oriented error estimates

We compute goal-oriented error estimates using the Dual Weighted Residual (DWR) technique [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] (see also [START_REF] Bangerth | Adaptive finite element methods for differential equations[END_REF][START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF][START_REF] Maday | Numerical analysis of a posteriori finite element bounds for linear functional outputs[END_REF][START_REF] Paraschivoiu | A hierarchical duality approach to bounds for the outputs of partial differential equations[END_REF][START_REF] Paraschivoiu | A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems: application to the control of pointwise errors[END_REF]), which is inspired from duality arguments (see, e.g., [START_REF] Eriksson | Introduction to adaptive methods for differential equations[END_REF]). We follow the framework described in [START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF], with some minor changes and adaptations. Let us consider u h the solution to Problem (3). The weak residual is defined for all

v ∈ V by r(v) := l(v) -a(u h , v).
Let z denote the solution to the dual problem:

Find z ∈ V such that a(v, z) = J(v), ∀ v ∈ V. (4)
The DWR method, in a linear setting, relies on the fundamental observation that

J(u) -J(u h ) = a(u, z) -a(u h , z) = l(z) -a(u h , z) = r(z). (5)
From this, we design an error estimator of J(u) -J(u h ) as an approximation of the residual r(z). We detail the different steps below.

Numerical approximation of the dual problem and global estimator

The exact solution z to the dual system (4) is unknown in most of the practical situations, and thus needs to be approximated. Let us consider a finite element space V h ⊂ V. This space is assumed to be finer than V h , for instance, made of continuous piecewise polynomials of order k + 1. The approximation z h of the solution to the dual problem z is obtained by solving the following approximate dual problem

Find z h ∈ V h such that a( v h , z h ) = J( v h ), ∀ v h ∈ V h . (6)
We define

η h := |r( z h )| (7)
as the global estimator that approximates the residual r(z).

Derivation of local estimators

Following [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF], we provide a local estimator of the error |J(u) -J(u h )|, that can be written in a general form

K∈K h η K , η K := K R K • ( z h -i h z h )dx + E∈E K E R E,K • ( z i h -i h z h )ds , ∀K ∈ K h , (8)
where the notation i h stands for the Lagrange interpolant onto V h . The local element-wise and edge-wise residuals are given explicitly by

R K := f K + div σ A (u h ) and R E,K :=    - 1 2 [[σ A (u h )]] E,K if E ∈ E int K , F E -σ A (u h )n E,K if E ∈ E N K , where σ A := σ(u h ) + βT (e A ⊗ e A )χ A .
The notation χ A stands for the indicator function of ω A , i.e. χ A = 1 in ω A and χ A = 0 elsewhere. The quantity σ A represents the sum of passive and active contributions within the stress field. The quantity f K (resp. F E ) is a computable approximation of f (resp. F).

The following bound always holds

η h ≤ K∈K h η K ,
since compensation effects (balance between positive and negative local contributions) can occur for η h , see, e.g., [START_REF] Nochetto | A safeguarded dual weighted residual method[END_REF]. Thus η h is expected to be sharper than

K∈K h η K . In practice, K∈K h
η K aims at quantifying the local errors for mesh refinement.

Remark 2.1. Each local estimator η K is made up of two contributions. On one hand, the residuals R K and R E,K represent the local error in the natural norm. On the other hand, the contribution ( z i h -i h z h ) coming from the dual problem can be interpreted as a weight (or a sensitivity factor) that measures the local impact on the quantity of interest J(•), see, e.g., [14, Remark 3.1].

Remark 2.2. In [START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF] the local residuals R K and R E,K are computed implicitly through local problems, in a generic fashion. No significant difference has been observed numerically between their technique and an explicit computation.

Remark 2.3. We have chosen to compute z h through the approximate dual system computed in V h ⊂ V ( i.e. the space made of continuous piecewise polynomials of order k +1). Other strategies are possible: see, e.g., [14, Section 5.1] for a discussion. For example, the authors of [START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF] use extrapolation of the approximate dual system computed in V h . We can also mention [START_REF] Becker | Weighted marking for goal-oriented adaptive finite element methods[END_REF], where the weight is estimated using a residual a posteriori error estimate for the dual system, approximated in V h . The aforementioned techniques are cheaper since the same space is used for the primal and dual solutions, but they can be less accurate.

Algorithm for goal-oriented mesh refinement

In the last sections, we have described the different steps to construct the global and local error estimators. Using the Dörfler marking strategy [START_REF] Dörfler | A convergent adaptive algorithm for Poisson's equation[END_REF], we now describe, in Algorithm 1, a simple algorithm to refine the mesh by taking into account these quantities. In this algorithm, there are two independent numerical parameters: first a parameter 0 < α ≤ 1 that controls the level of refinement in Dörfler marking, and then a tolerance threshold ε > 0 for the global estimator, that serves as a stopping criterion.

Results

We present numerical results for two different test cases: the biomechanical response of both a human tongue and an artery predicted using finite element analysis, inspired from studies [START_REF] Bijar | Atlas-based automatic generation of subject-specific finite element tongue meshes[END_REF] and [START_REF] Floc'h | Vulnerable Atherosclerotic Plaque Elasticity Reconstruction Based on a Segmentation-Driven Optimization Procedure Using Strain Measurements: Theoretical Framework[END_REF], respectively. We propose to assess the discretization error for the two quantities of interest

J 1 (u) := ω (u x + u y ) dx and J 2 (u) := ω div u dx, (9)
Algorithm 1: Refinement algorithm

Initialization :

Select an initial triangulation K h of the domain Ω.

Build the finite elements spaces V h and V h .

While η h > do 1. Compute u h ∈ V h : a(u h , v h ) = l(v h ), ∀v h ∈ V h . 2. Compute z h ∈ V h : a( v h , z h ) = J( v h ), ∀ v h ∈ V h . 3.
η K := K R K •( z h -i h z h )dx + E∈E K E R E,K •( z i h -i h z h )ds , ∀K ∈ K h .
6. Sort the cells {K 1 , ..., K N } by decreasing order of η K . 7. Dörfler marking: mark the first M * cells for refinement where

M * := min M ∈ N M i=1 η K i ≥ α K∈K h η K .
8. Refine all cells marked for refinement (and propagate refinement to avoid hanging nodes). 9. Update correspondingly the finite element spaces V h and V h .

where u x and u y are the two components of u in a Cartesian basis. The first quantity J 1 (u) is physically related to the displacement in the region of interest ω ⊂ Ω. This corresponds to a quantity that can easily be measured experimentally and that is therefore of practical interest. The second quantity J 2 (u) physically corresponds to the internal strain I 1 = tr(ε(u)). This is also of practical interest because many of the mechanisms driving the onset of pathologies are related to shear strains or principal strains. The region of interest ω will be specified in both situations. All the simulations of this section are performed with Lagrange finite elements of degree k = 2, and the space V h in which z h is computed is built from Lagrange finite elements of degree k = 3. In Algorithm 1, the parameter α for Dörfler marking is fixed at 0.8, and the stopping criterion ε will be specified for each application. In the following, the exact value of J(u) is unknown but is estimated using computations on a very fine uniform mesh.

Human tongue with fiber activation

In the first example, we focus on the case study for the activation of the posterior genioglossus (GGp), that is a lingual muscle located at the root of the tongue and inserted in the front to the mandible. The activation of this muscle compresses the tongue in the lower part and generates a forward and upward movement of the tongue body, because of the incompressibility of tongue tissues, for example during the production of the phonemes /i/ or /s/. The 2D mesh used in this example has been derived from the generic 3D mesh presented in [START_REF] Bijar | Atlas-based automatic generation of subject-specific finite element tongue meshes[END_REF] where the authors developed a process to generate subject-specific meshes. More precisely an automatic atlas-based method was proposed that generates subject-specific meshes via a registration guided by Magnetic Resonance Imaging. The domain Ω is depicted in Figure 1 (left). The width and height of the tongue are respectively equal to 73.8 mm and 53.7 mm. For the passive tissue material properties, we use the values reported in [START_REF] Gerard | Non-linear elastic properties of the lingual and facial tissues assessed by indentation technique: Application to the biomechanics of speech production[END_REF] based on indentation experiments on a cadavers tongue. The authors initially proposed an incompressible two parameter Yeoh hyperelastic material model and fitted the material constants to the data. In this work, a linear elastic material model is assumed. According to [START_REF] Tracqui | Transmission of mechanical stresses within the cytoskeleton of adherent cells: a theoretical analysis based on a multi-component cell model[END_REF], linearisation of the model proposed in [START_REF] Gerard | Non-linear elastic properties of the lingual and facial tissues assessed by indentation technique: Application to the biomechanics of speech production[END_REF] yields E 6c10 = 0.6 MPa. For the sake of simplicity Poisson ratio is assumed to be ν = 0.4. No volumic force field is applied: f = 0. The direction of the fibers e A is depicted in Figure 1 (center) and corresponds approximately to the posterior genioglossus muscle [START_REF] Bijar | Atlas-based automatic generation of subject-specific finite element tongue meshes[END_REF]. Other parameters for fiber activation have been chosen as T = 2e-5 MPa and β = 1. The tongue is attached to the hyoid bone and to the mandible, which are supposed to be fixed. This leads to a homogeneous Dirichlet boundary condition such as depicted in Figure 1 (right). On the remaining part of the boundary a homogeneous Neumann condition (F = 0) is applied. The orange part depicts the region ω A where fibers are supposed to be located. The green part depicts the region of interest ω for the computation of J 1 and J 2 . The resulting displacement is depicted in Figure 2 (left). We computed the relative displacement and the strain intensity, which maximal values are of 5.7 % and 4.8 %, respectively: thus the small displacement and small strain assumptions are both verified in this case. The parameter T has been chosen accordingly in order to respect these assumptions. In Figure 2, the dual solutions for the quantities of interest J 1 (center) and J 2 (right) are represented. As mentioned in Remark 2.1, the dual solution z is used as a weight in the computation of the estimators, and influences the local refinement. We present the final mesh after 2 and 8 iterations of Algorithm 1 for both quantities of interest J 1 and J 2 , in Figure 3 and Figure 4, respectively. We first remark that the refinement occurs in some specific regions such as those near Dirichlet-Neumann transitions and concavities on the boundary. Note as well that the refinement is stronger for J 2 at the boundary of the region of interest ω,. versus N , the number of cells of the mesh, both for uniform refinement (blue) and adaptive refinement (red). The stopping criterion ε has been fixed to 2e-4 and 1e-6, respectively. In each situation, we observe that, as expected, adaptive refinement performs better: not only it leads to a lower error but it also converges much faster when the number of cells N is increased. 10 -4 10 -3

10 -2 N |J 1 (u) -J 1 (u h )|/|J 1 (u)|
Uniform refinement Adaptive refinement 10 3 10 4

10 -4

10 -3

10 -2 N |J 2 (u) -J 2 (u h )|/|J 2 (u)|
Uniform refinement Adaptive refinement Finally in Figure 6 we depict the efficiency indices for the global estimator η h and the sum of local estimators K η K . For both quantities J 1 and J 2 , the two estimators provide an estimation of the discretization error with an efficiency index around 1. In the case of J 2 , we observe a slight overestimation for K η K and a slight underestimation for η h .

10 0 N η h /|J 1 (u) -J 1 (u h )| K η K /|J 1 (u) -J 1 (u h )| 10 3 10 4 10 -1 10 0 N η h /|J 2 (u) -J 2 (u h )| K η K /|J 2 (u) -J 2 (u h )|

Human artery with fiber activation

As a second example we showcase the performance of the proposed algorithm for the analysis of the mechanical response of an artery with vulnerable coronary plaque to internal loading. Rupture of the cap induces the formation of a thrombus which may obstruct the coronary artery, cause an acute syndrome and the patient death. The geometry (see Figure 7 (left)) comes from [START_REF] Floc'h | Vulnerable Atherosclerotic Plaque Elasticity Reconstruction Based on a Segmentation-Driven Optimization Procedure Using Strain Measurements: Theoretical Framework[END_REF] where the authors develop a methodology to reconstruct the thickness of the necrotic core area and the calcium area as well as the Young's moduli of the calcium, the necrotic core and the fibrosis. Their objective is the prediction of the vulnerable coronary plaque rupture. As represented in Figure 7 (left), the diameter of the Fibrosis is equal to 5 mm. Following [START_REF] Floc'h | Vulnerable Atherosclerotic Plaque Elasticity Reconstruction Based on a Segmentation-Driven Optimization Procedure Using Strain Measurements: Theoretical Framework[END_REF], we set different elastic parameters in each region: E = 0.011 MPa, ν = 0.4 in the necrotic core and E = 0.6 MPa, ν = 0.4 in the surrounding tissue. No volumetric force field is applied: f = 0. We consider muscle fibers only in the media layer, where smooth muscle cells are supposed to be perfectly oriented in the circumferential direction e A = e θ , where (e r , e θ ) is the basis for polar coordinates, see Figure 7 (center). Other parameters for fiber activation have been chosen as T = 0.01 MPa and β = 1. As depicted in Figure 7 (right), the artery is fixed on the red portion of external boundary Γ D . Elsewhere, on the remaining part of the boundary, a homogeneous Neumann condition is applied: F = 0. In the same figure, the green part represents the region of interest ω. This choice is relevant in the study of vulnerable coronary plaque rupture. As in the previous example, we computed the relative displacement and the strain intensity, which maximal values are of 6.15 % and 0.3 %, respectively. This ensures that small displacement and small strain assumptions are verified. Figure 8 depicts the magnitude of the solution in terms of displacements (left) and the dual solutions associated to J 1 (center) and J 2 (right). In Figure 9, we present the final mesh after 2 and 6 iterations of Algorithm 1 for the quantity of interest J 1 . As in the previous example, the refinement occurs in some specific regions, such as those near Dirichlet-Neumann transitions and concavities on the boundary. Our results also show that the proposed method leads to the strong refinement near the interface between the necrotic core and the fibrosis, where stresses are localized because of the material heterogeneity. Conversely to the previous example, the refined meshes obtained for J 2 (not depicted) are very similar to those obtained for J 1 . Figure 10 (left) depicts the relative goal-oriented error |J 1 (u) -J 1 (u h )|/|J 1 (u)| versus the number N of cells in the mesh, both for uniform refinement (blue) and adaptive refinement (red). The stopping criterion ε has been fixed at 5e-6. In Figure 10 (right), we depict the efficiency indices for the global estimator η h and the sum of local estimators K η K . The same observations as in the previous example can be stated, and the estimators provide acceptable value of the discretization error. Moreover η h performs better though it still underestimates slightly the error. Results we obtained for the quantity J 2 are very similar.

N |J 1 (u) -J 1 (u h )|/|J 1 (u)|
Uniform refinement Adaptive refinement

10 3 10 4 10 -1 10 0 N η h /|J 1 (u) -J 1 (u h )| K η K /|J 1 (u) -J 1 (u h )|

Discussion

In the first part, we discuss about the ability of the proposed methodology to assess and reduce the discretization error. In a second part, we comment on some further issues to improve and guarantee the accuracy of the error estimator, and to optimize the mesh refinement algorithm. Finally we address the issue of tackling more complex problems that arise in current practice for clinical biomechanics, and point out the main limitations of the current study as well as some perspectives.

Towards quantification of the discretization error: first achievements

The numerical results obtained in the last section show the ability of the proposed framework to provide relevant information about the discretization error : though the global estimator η h provides only an approximation of the error in the quantity of interest |J(u) -J(u h )|, this is often sufficient in practice. Moreover, the local estimators η K provide a means to evaluate "relative" errors and thereby drive mesh refinement (Algorithm 1). Both the local and global errors can be significantly reduced without much computational effort. For instance, in the first test-case 3.1, and for J 1 , the error is reduced by a factor of almost 4, after two successive refinements, and with only 20 % of extra cells. In order to quantify more precisely the computational gains provided by the adaptive procedure, the computational time required to compute the error estimator and to regenerate or adapt the mesh should be thoroughly computed and analyzed, as was done for three-dimensional fracture treated by enriched finite element methods [START_REF] Jin | Error-controlled adaptive extended finite element method for 3d linear elastic crack propagation[END_REF]. Let us emphasize the well-known fact that sources of discretization errors are local, and concentrated mostly in regions where the solution is not smooth, e.g. subjected to strong variations, discontinuities or singularities. As a consequence, uniform refinement is highly suboptimal, while adaptive refinement performs much better by optimizing the number of elements, their size and location within the domain. Moreover, the proposed adaptive procedure is fully automatic, and no a priori knowledge of the critical regions is needed. For goal-oriented error estimation, the refined mesh obtained by the algorithm can in fact be counter-intuitive, because it is driven by the sensitivity of the quantity of interest with respect to the local error. This sensitivity is obtained by solving the dual problem (see for instance Figure 2 in Section 3.1) whose solution is, indeed, often not intuitive and difficult to interpret from a physical viewpoint. In comparison to widespread error techniques implemented in most of commercial finite element software, the DWR technique allows to estimate and to improve the error for an arbitrary quantity of interest J. Each practitioner can choose the relevant quantity of interest J and obtain an approximation of the error on this quantity of interest |J(u) -J(u h )|, as well as a map of the local error. The authors emphasize that the results obtained in the current study also demonstrate that the optimal refinement strategy depends significantly on the choice of the quantity of interest J. In general, such a goaloriented refinement strategy leads to meshes which may differ significantly from those obtained by minimizing the error in energy. Remark that such goal-oriented approaches were also developed for the Zienkiewicz-Zhu error estimators [132] in [START_REF] González-Estrada | Mesh adaptivity driven by goal-oriented locally equilibrated superconvergent patch recovery[END_REF][START_REF] González-Estrada | Locally equilibrated stress recovery for goal oriented error estimation in the extended finite element method[END_REF] and for explicit residual based estimates in [START_REF] Rüter | Goal-oriented explicit residual-type error estimates in xfem[END_REF] and [START_REF] Wick | Goal functional evaluations for phase-field fracture using PU-based DWR mesh adaptivity[END_REF].

Some further mathematical / computational issues

It is desired that the global estimator η h provide reliable information on the error in the quantity of interest |J(u) -J(u h )|, providing quality measures to the user. In theory, this error is a guaranteed upper bound, with an explicit constant equal to 1. Yet, the theory assumes that the dual solution z is exactly known. This is never the case in practice as the dual problem is also solved using finite elements. Our numerical experiments show, however, that |J(u)-J(u h )| is estimated with reasonable accuracy and that the effectivity indices are close to 1, meaning that the approximate error on the quantity of interest is close to the (unknown) exact error on this quantity. The numerical experiments provided in this paper confirm those of the literature on DWR technique, e.g., [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF][START_REF] Rognes | Automated goal-oriented error control I: Stationary variational problems[END_REF] showing that the DWR is, in most situations, a reliable approach to compute goal-oriented error estimates. However, in certain situations, the DWR estimator is not as reliable as desired, since the effect of approximating the dual solution is difficult to control. This issue has been already pointed in the literature: see e.g. [START_REF] Carstensen | Estimation of higher Sobolev norm from lower order approximation[END_REF][START_REF] Nochetto | A safeguarded dual weighted residual method[END_REF][START_REF] Ainsworth | Guaranteed computable bounds on quantities of interest in finite element computations[END_REF] and earlier considerations in, e.g., [START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF][START_REF] Bangerth | Adaptive finite element methods for differential equations[END_REF]. Especially, in [START_REF] Nochetto | A safeguarded dual weighted residual method[END_REF] a simple situation where η h provides a poor estimation on a coarse mesh is detailed. There is up to now no simple, cheap and general technique to address this issue, but first solutions have been suggested in [START_REF] Carstensen | Estimation of higher Sobolev norm from lower order approximation[END_REF][START_REF] Nochetto | A safeguarded dual weighted residual method[END_REF][START_REF] Ainsworth | Guaranteed computable bounds on quantities of interest in finite element computations[END_REF]. They consist in modifying the DWR estimator so as to take into account the approximation of z. This is a stimulating perspective for further research. Moreover, the issue of computing a cheaper approximation of z, without compromising the reliability and efficiency of the estimator still needs to be addressed in depth. Concerning mesh refinement, though the local estimator η K combined with Algorithm 1 provides acceptable results, no effort has been spent on finding the value of parameter α in the Dörfler marking that yields improved refined meshes. On this topic, our global strategy for error estimation and mesh refinement is only a first attempt, and can be improved. For instance, in [START_REF] Becker | Weighted marking for goal-oriented adaptive finite element methods[END_REF], an adaptive method based on specific weighting of the residuals of the primal and dual problems has been designed, and leads to quasi-optimal adapted meshes. Such a method could be tested and compared to the current one.

Applicability for patient-specific biomechanics?

Though the preliminary results presented in this paper demonstrate the relevance and practicability of a posteriori error estimators for providing quality control in quantities of interest to the biomechanics practitioner, and to drive mesh adaptation, much effort is still needed for the approaches developed here to address practical, personalized, clinicallyrelevant Finite Element (FE) simulations for biomechanical applications. First, the compressible linear framework considered here is inadequate in practice and must be replaced by a fully non-linear, incompressible, time and history dependent model [START_REF] Payan | Biomechanics of living organs : hyperelastic constitutive laws for finite element modeling[END_REF]. Non-linearities also occur due to boundary conditions, when, for instance, contact is present [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF]. Moreover, most of the quantities of interest in biomechanics are non-linear (norm of the displacements, local shear stress, maximum admissible stress and strain, etc). It is important to point out here that the DWR method for goal-oriented error estimation is already capable of tackling non-linearities: see, e.g., [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] for the general framework, and, e.g., [START_REF] Larsson | Strategies for computing goal-oriented a posteriori error measures in non-linear elasticity[END_REF][START_REF] Whiteley | Error estimation and adaptivity for incompressible hyperelasticity[END_REF] for first applications in non-linear elasticity and [START_REF] Rüter | Goal-oriented explicit residual-type error estimates in xfem[END_REF] for fracture mechanics. Nevertheless, this non-linear framework needs to be adapted and tested in the specific case of hyperelastic soft-tissue. The major limitation of our work is that it assumes that the mathematical model used to describe the biomechanics problem is able to reproduce the physical reality. Unfortunately, in general, selecting the proper mathematical model for a given biomechanics problem is probably the most challenging part of the simulation process. The large, and increasing, number of papers dealing with the choice of constitutive model, for example, testifies for this difficulty. For a wide range of problems, indeed, modeling errors are the most significant. Estimating rigorously and systematically the impact of these errors is extremely challenging, in particular when dealing with patient-specific simulations. Dealing with this issue is the focus of ongoing research in our teams but is far beyond the scope of this paper. We would nonetheless like to make the following remarks. The first problem which must be addressed is the choice of a model (hyperelastic, viscous, porous, single/multi-scale...). The chosen model has parameters which must be estimated through inverse analysis. Once estimates, or probability distributions for these parameters are available, their importance on quantities of interest must be evaluated, through sensitivity analysis and uncertainty quantification. The major difficulty is, therefore, to select the proper model, and its parameters for a given patient. As in vivo experiments are in general not possible, data must be extracted as the patient is being treated, e.g. during an operation. This can be done using Bayesian methods, which provide a reconciliation between expert knowledge on patient cohorts (prior) and actual properties of a given patient [START_REF] Rappel | Bayesian inference for the stochastic identification of elastoplastic material parameters: Introduction, misconceptions and insights[END_REF][START_REF] Rappel | Bayesian inference to identify parameters in viscoelasticity[END_REF]. Real-time machine-learning-like methods such as Kalman filters demonstrated as well promising results [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: perspectives in medical imaging[END_REF][START_REF] Haouchine | Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery[END_REF]. To evaluate the effects of uncertainties on such material parameters, accelerated Monte-Carlo methods are possible avenues of investigation [START_REF] Hauseux | Accelerating monte carlo estimation with derivatives of high-level finite element models[END_REF]. An exciting question is the comparative usefulness and combination of physical models (potentially learnt during medical treatment) and machine-learning algorithms, mostly based on data acquired during the intervention. Last but not least, note that the DWR method is based on optimal control principles, that makes it suitable for extensions to parameter calibration (viewed as an optimal control problem). In such a setting, it allows to combine sensitivity analysis with goal-oriented a posteriori error estimation, see [START_REF] Becker | Mesh refinement and numerical sensitivity analysis for parameter calibration of partial differential equations[END_REF]. In the same spirit, the interplay between a posteriori error estimation and uncertainty quantification has been object of recent research interests [START_REF] Eigel | An adaptive multilevel Monte Carlo method with stochastic bounds for quantities of interest with uncertain data[END_REF][START_REF] Guignard | A posteriori error estimation for elliptic partial differential equations with small uncertainties. Numerical Methods for Partial Differential Equations[END_REF]. We also note that if users can obtain some estimate, even rough, of modeling errors, they will also be able to compare discretization and model errors. This enables the coarsening of the mesh if the discretization error is unnecessarily small in comparison to the modeling error as is done, e.g. in [START_REF] Akbari Rahimabadi | Scale selection in nonlinear fracture mechanics of heterogeneous materials[END_REF] and [START_REF] Talebi | Molecular dynamics/xfem coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture[END_REF][START_REF] Budarapu | Efficient coarse graining in multiscale modeling of fracture[END_REF][START_REF] Talebi | A computational library for multiscale modeling of material failure[END_REF][START_REF] Silani | A semi-concurrent multiscale approach for modeling damage in nanocomposites[END_REF][START_REF] Talebi | Concurrent multiscale modeling of three dimensional crack and dislocation propagation[END_REF] for adaptive scale selection. Conversely, for specific applications where modeling errors are small or moderate, the mesh can be refined efficiently to increase the precision. With our methodology, practitioners spending a large amount of time and effort in patientspecific mesh generation can obtain useful information on the impact of the quality of the mesh on quantities of interest to them. This information goes well beyond purely geometrical criteria for the regularity of the elements which are typically provided in commercial software. This information can be used directly to optimize the choice of the discretization/mesh in view of minimizing the error on a specific quantity of interest. Fast/real-time numerical methods which provide real-time predictions have been intensively researched since the beginning of the 1990's. Those approaches are critical to build surgical planning and guidance tools, for example. Reliable error estimation is critical in these situations to guarantee the accuracy, but has been extremely scarcely addressed in the literature. As a first step in this direction, the recent work of [START_REF] Bui | Real-time error control for surgical simulation[END_REF] provides a real-time mesh refinement algorithm for needle insertion. Mesh refinement is driven by a ZZ error estimate, for the global norm. It would be interesting to extend such a method for goal-oriented error estimation, e.g. on the motion of a target, or reaction/friction force along the needle shaft. We should also mention alternative approaches to (implicit, standard) finite elements for fast nonlinear finite element analysis: for instance the solution of total lagrangian formulation of the equilibrium equations on graphics processing unit for neurosurgical simulation [START_REF] Joldes | Suite of finite element algorithms for accurate computation of soft tissue deformation for surgical simulation[END_REF][START_REF] Joldes | Non-locking tetrahedral finite element for surgical simulation[END_REF][START_REF] Joldes | Real-time nonlinear finite element computations on GPU Application to neurosurgical simulation[END_REF], or model order reduction techniques for the real-time, interactive simulation of tissue tearing during laparoscopic surgery [START_REF] Niroomandi | Real-time simulation of surgery by reduced-order modeling and X-FEM techniques[END_REF][START_REF] Niroomandi | Real-time simulation of biological soft tissues: a PGD approach[END_REF][START_REF] Quesada | ???? Haptic simulation of tissue tearing during surgery[END_REF]. A perspective consists in extending the current framework to such numerical methods where error control is particularly challenging. For explicit approaches, the interplay between the choice of the time-step and that of the mesh size is a difficult topic, especially for domains with significant stiffness differences where adaptive and multi-time-step schemes should be investigated.

Figure 1 :

 1 Figure 1: Tongue model: initial geometry (left), fiber orientation (center) and region of interest (right).

Figure 2 :

 2 Figure 2: Tongue model: displacement (left), dual solutions for J 1 (center) and for J 2 (right).

Figure 3 :

 3 Figure 3: Tongue mesh: refinement driven by the QoI J 1 . Initial mesh (left) with 426 cells and a relative error of 1.07e-2, adapted meshes after 2 iterations (center) with 523 cells and a relative error of 2.82e-3 and after 8 iterations (right) with 5143 cells and a relative error of 3.82e-05.

Figure 5

 5 Figure 5 depicts the relative goal-oriented errors |J 1(u)-J 1 (u h )|/|J 1 (u)| (left) and |J 2 (u)-J 2 (u h )|/|J 2 (u)| (right)versus N , the number of cells of the mesh, both for uniform refinement (blue) and adaptive refinement (red). The stopping criterion ε has been fixed to 2e-4 and 1e-6, respectively. In each situation, we observe that, as expected, adaptive refinement performs better: not only it leads to a lower error but it also converges much faster when the number of cells N is increased.

Figure 4 :

 4 Figure 4: Tongue mesh: refinement driven by the QoI J 2 . Initial mesh (left) with 426 cells and a relative error of 2.51e-2, adapted meshes after 2 iterations (center) with 766 cells and a relative error of 2.21e-3 and after 8 iterations (right) with 13513 cells and a relative error of 2.44e-5.

Figure 5 :

 5 Figure 5: Tongue model: relative error for the QoI J 1 (left) and J 2 (right) vs. the number N of cells in the case of uniform (blue) and adaptive (red) refinement.

Figure 6 :

 6 Figure 6: Tongue model: efficiency indexes for η h (blue) and K η K vs. the number N of cells for the QoI J 1 (left) and J 2 (right).

Figure 7 :

 7 Figure 7: Artery model: geometry (left), fiber orientation (center) and region of interest (right).

Figure 8 :

 8 Figure 8: Artery model: displacement (left), dual solution for J 1 (center) and for J 2 (right).

Figure 9 :

 9 Figure 9: Artery mesh: refinement driven by the QoI J 1 . Initial mesh (left) with 1242 cells and a relative error of 3.83e-1, adapted meshes after 2 iterations (center) with 2079 cells and a relative error of 5.25e-2 and after 6 iterations (right) with 15028 cells and a relative error of 3.37e-3.

Figure 10 :

 10 Figure 10: Artery model. Left: relative error for the QoI J 1 vs. the number N of cells in the case of uniform (blue) and adaptive (red) refinement. Right: efficiency indexes of η h (blue) and K η K vs. the number of cells N for the QoI J 1 .

  Evaluate the global error estimator η h = |r( z h )|. 4. If η h ≤ ε, then stop. 5. Evaluate the local estimators

Acknowledgements

The authors thank warmly Mathias Brieu for the organization of the EuroMech 595 workshop. They also thank Roland Becker, Jean-Louis Martiel, Jacques Ohayon and Yohan Payan for their support and helpful comments that allowed to improve the paper, as well as the following people for inspiring discussions: Jack Hale, Florence Hubert and Pascal Perrier. For funding, F. Chouly thanks Région Bourgogne Franche-Comté ("Convention Région 2015C-4991. Modèles mathématiques et méthodes numériques pour l'élasticité non-linéaire"), the Centre National de la Recherche Scientifique ("Convention 20 232789 DEFI InFIniTI 2017 -Projet MEFASIM") and the Agence Maths Entreprises (AMIES) ("Projet Exploratoire PEPS2 MethASim").


Luc Le Magoarou

Stéphane Paquelet

Bias-variance tradeoff in MIMO channel estimation

Keywords: Channel estimation, physical model, MIMO

come

Bias-variance tradeoff in MIMO channel estimation

I. INTRODUCTION

Multiple-input multiple-output (MIMO) communication systems allow for a dramatic increase in channel capacity, adding space to the classical time and frequency dimensions [1], [START_REF] Tse | Fundamentals of wireless communication[END_REF]. This is done by using several antennas at the transmitter (N t ) and at the receiver (N r ). The capacity of MIMO systems is maximized if the channel state is perfectly known at both ends of the link.

Channel estimation is deeply impacted by the transition from single antenna to MIMO systems. Indeed, it amounts to determine a complex gain for each transmit/receive antenna pair, the narrowband (single carrier) MIMO channel as a whole being usually represented as a matrix H ∈ C Nr×Nt of such gains. The number of real parameters to estimate is thus 2N r N t , which may be very large for massive MIMO systems, i.e. systems with up to several hundreds of antennas [START_REF] Rusek | Scaling up mimo: Opportunities and challenges with very large arrays[END_REF], [START_REF] Erik G Larsson | Massive mimo for next generation wireless systems[END_REF].

In a massive MIMO context, the classical least squares (LS) estimator is not adapted because the high dimensionality of the parameter space leads to an ill-posed problem. Therefore, in order to add prior information, it has been proposed classically to use Bayesian estimation and thus model the channel matrix as random, giving rise to estimators such as the linear minimum mean square error (LMMSE) [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF], [START_REF] Bazzi | Downlink training sequence design for fdd multiuser massive mimo systems[END_REF]. Another possibility is to use a parametric model based on the physics of wave propagation, in which the channel is expressed as a sum of p paths [START_REF] Akbar | Deconstructing multiantenna fading channels[END_REF]. Whereas LS and LMMSE estimators have been studied extensively in terms of optimal training sequences and performance [START_REF] Biguesh | Training-based mimo channel estimation: a study of estimator tradeoffs and optimal training signals[END_REF], a similar study is still lacking for channel estimators based on a physical model. Contributions. In this paper, the performance of MIMO channel estimators based on a physical model is studied, highlighting a bias-variance tradeoff driven by the number of paths p considered by the model. Moreover, physical channel estimators are compared to classical ones in a massive MIMO context. Finally, the impact of channel estimation error on the system data rate is theoretically assessed. This work has been performed in the framework of the Horizon 2020 project ONE5G (ICT-760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.

II. PROBLEM FORMULATION

Notations. Matrices and vectors are denoted by bold upper-case and lower-case letters: A and a (except 3D "spatial" vectors that are denoted -→ a ); the ith column of a matrix A by a i ; its entry at the ith line and jth column by a ij . A matrix transpose, conjugate and transconjugate is denoted by A

A. Channel estimation

Consider a narrowband block fading channel between a transmitter and a receiver with respectively N t and N r antennas. It is represented by the matrix H ∈ C Nr×Nt , in which h ij corresponds to the channel between the jth transmit and ith receive antennas.

In order to carry out channel estimation, N s known pilot symbols are sent through the channel by each transmit antenna. The corresponding training matrix is denoted X ∈ C Nt×Ns . It obeys a power constraint of the form X 2 F = P where P is the total training energy. In addition, the transmit power is assumed constant during the training phase, leading to x i 2 2 = P e ,i = 1,...,N s , which implies P = N s P e . The signal at the receive antennas is thus expressed as HX+N, where N is a noise matrix (usually vec(N) ∼ CN (0,σ 2 Id)). One useful quantity is the potential signal to noise ratio (pSNR) defined as

pSNR P e H 2 F σ 2 . (1)
It is an upper bound on the signal to noise ratio (SNR), which is used to quantify the performance of channel estimators in section III. Note that pSNR and SNR coincide for a rank one channel combined with an optimal precoding at the transmitter, and that with no precoding at all, SNR = pSNR-log 10 (N t ) (in dB). Due to the high cost and power consumption of millimeter wave Radio Frequency (RF) chains, it has been proposed to have less RF chains than antennas in both the transmitter and receiver [START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF]- [START_REF] Akbar | Millimeter-Wave MIMO Transceivers: Theory, Design and Implementation[END_REF]. Such systems are often referred to as hybrid architectures. Mathematically speaking, this translates into sensing the channel through analog precoders v i ∈ C Nt , i = 1,...,N RF , with N RF ≤ N t , as well as observing the signal at the receiver through analog combiners denoted w j ∈ C Nr , j = 1,...,N c . The observed data is thus expressed in all generality as

Y = W H HX+W H N, (2)
where W (w 1 ,...,w Nc ) and the training matrix is constrained to be of the form X = VZ, where Z ∈ C NRF×Ns is the digital training matrix, and V ∈ C Nt×NRF is the analog precoding matrix. For more convenience and to lighten notations, one can introduce the equivalent vectorized observation model:

y = Mh+n, (3)
where y vec(Y) ∈ C NcNs , h vec(H) ∈ C NrNt , n vec(W H N) ∈ C NcNs and M X T ⊗W H ∈ C NcNs×NrNt is called the observation matrix. In the remaining of the paper, it is assumed that the analog combiners are mutually orthogonal and of unit norm: w H i w j = δ ij . This yields W 2 F = N c and n ∼ CN (0,σ2 Id). Objective. Channel estimation aims at retrieving h (or equivalently H) from the observation of y, knowing the observation matrix M and the distribution of the noise vector n. The channel estimator is denoted ĥ (or equivalently Ĥ in matrix form). At first sight, the search space of channel estimation is thus of dimension 2N r N t , which may be very large in massive MIMO systems (up to several thousands). For this reason, classical estimation methods such as the least squares (LS) may not be appropriate. In order to overcome this limitation, some information about the channel has to be used to regularize the problem. For example, the channel can be considered as a random vector whose distribution is known, yielding Bayesian channel estimation [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF], [START_REF] Bazzi | Downlink training sequence design for fdd multiuser massive mimo systems[END_REF]. Another possibility is to use a physical channel model, as presented in the next subsection.

B. Physical channel model

Inspired by the physics of propagation under the plane waves assumption, it has been proposed to express the channel as a sum of rank one matrices, each corresponding to a single physical path between transmitter and receiver. This approach is used by most MIMO channel simulators [START_REF] Sun | A novel millimeter-wave channel simulator and applications for 5g wireless communications[END_REF]-[14] and has been validated by propagation measurements. Adopting this kind of modeling and generalizing it to take into account any three-dimensional antenna array geometry, channel matrices take the form

H = P i=1 c i e r ( -→ u r,i )e t ( -→ u t,i ) H , (4)
where P is the total number of considered physical paths (up to a few hundreds in classical simulators), c i ρ i e jφi is the complex gain of the ith path, -→ u t,i is the unit vector corresponding to its direction of departure (DoD) and -→ u r,i the unit vector corresponding to its direction of arrival (DoA)(both being described in spherical coordinates by an azimuth angle η and an elevation angle ψ). The steering vectors e r ( -→ u ) ∈ C Nr and

e t ( - → u ) ∈ C Nt are defined as (e x ( - → u )) j = 1 √ Nx e -j 2π λ --→ ax,j . - → u for
x ∈ {r,t}. The set A x { --→ a x,1 ,..., ---→ a x,Nx } gathers the positions of the antennas with respect to the centroid of the array (transmit if x = t, receive if x = r). The same model is used in [START_REF] Le | Parametric channel estimation for massive MIMO[END_REF]. Estimation model. It makes sense to have the channel estimator Ĥ take a form similar to that of the physical channel H given in (4). It has indeed been proposed in several contexts [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]- [START_REF] Venugopal | Channel estimation for hybrid architecture-based wideband millimeter wave systems[END_REF], leading to estimators expressed as

Ĥ = p k=1 d k e r ( --→ v r,k )e t ( -→ v t,k ) H , ( 5
)
where p is the number of paths considered for estimation, called virtual paths. In practice, one takes p P (with p at most a few dozens), because several physical paths can be merged into a single virtual path without harming a lot description accuracy. Indeed, steering vectors associated to close enough directions are almost collinear, yielding a limited angular resolution for the system [START_REF] Akbar | Deconstructing multiantenna fading channels[END_REF]. Moreover, a smaller p leads to a better conditioned problem, as will be evidenced in the next section. Let us define the model set

M p A ∈ C Nr×Nt ,A = p k=1 d k e r ( --→ v r,k )e t ( -→ v t,k
) H to which the estimator made of p virtual paths belongs. By abuse of notation, the vectorized version of this set will be denoted the same way, such that

M p = a ∈ C NrNt , a = p k=1 d k e t ( -→ v t,k ) * ⊗ e r ( --→ v r,k ) . Such sets obey the inclusion relation M q ⊂ M q+1 .
The estimation model can be seen as a parametric model with

parameter vector θ = {θ (k) (ρ k ,φ k ,η r,k ,ψ r,k ,η t,k ,ψ t,k ), k = 1,...,p}.
There are thus 6p real parameters in this model (the complex gain, DoD and DoA of every path are described with two parameters each). Of course, the model is most useful for estimation in the case where 6p 2N r N t , since the number of parameters is thus greatly reduced. Conversely, taking p too small may harm the descriptive power of the model. The number of virtual paths p thus drives a tradeoff between complexity and expressiveness, whose study is the object of the next sections.

III. PERFORMANCE OF ESTIMATORS

Let us consider an estimator ĥ of a channel h. Its performance is assessed using the relative mean square error (rMSE) defined as

rMSE( ĥ) E h-ĥ 2 2 h 2 2 . (6)
For any estimator, this error can be decomposed [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF] as

rMSE( ĥ) = h-E[ ĥ] 2 2 h 2 2 + Tr[cov( ĥ)] h 2 2 , (7)
where the first term called bias represents the error due to the deviation of the average estimate from the true channel and the second called variance represents the error due to fluctuations around the average. In this section and the following, the rMSE of various channel estimators is studied. The emphasis is put on physical channel estimators, for which a bias-variance tradeoff driven by the number of virtual paths p is exhibited.

The study of this tradeoff is made possible by the definition of an appropriate oracle estimator.

A. Classical channel estimators

Let us now review the performance of the most classical MIMO channel estimators that are the least squares (LS) and the linear minimum mean square error (LMMSE). This subsection essentially summarizes the main results of [START_REF] Biguesh | Training-based mimo channel estimation: a study of estimator tradeoffs and optimal training signals[END_REF], in a way that takes into account hybrid systems. LS. The least squares estimator is an unbiased estimator defined (using the vectorized notation of (3)) as

ĥLS argmin g y-Mg
It requires to send N s ≥ N t pilot symbols to exist (for the matrix M H M to be invertible). In that case, the variance (and thus the rMSE) is minimized if the observation matrix obeys the condition

M H M opt = P e N c N s N r N t Id, (9)
leading to the optimal performance

rMSE opt ( ĥLS ) = N r N t pSNR N r N t N c N s . ( 10
)
Note that in the case where N s = N t (which corresponds to the minimum number of pilot symbols) and N c = N r (which corresponds to a non-hybrid receiver), the optimal rMSE is proportional to N r N t , showing that such an estimator is not adapted to massive MIMO systems in which this quantity is large. LMMSE. The linear minimum mean square error estimator is a Bayesian estimator, meaning that the channel is assumed to be random and to follow a known distribution h ∼ CN (0,R).

A channel realization can then be seen as drawing at random a user location in a given region in a fixed environment: the smaller the region the closer to singular R. Note that in practice the covariance matrix R has to be estimated using previous channel estimates.If R is perfectly known, the LMMSE is the linear estimator minimizing the MSE, it takes the form

ĥLMMSE RM H MRM H +σ 2 Id -1 y, (11)
with

RM H MRM H +σ 2 Id -1 = argmin A E h-Ay 2 2 ,
where the expectation is taken over both the noise and the channel distributions. In that case, the rMSE is minimized if the observation matrix obeys the condition (for a high enough SNR, see [8, section V] for more details)

M H M opt = P e N c N s N r N t + σ 2 N r N t Tr R -1 Id-σ 2 R -1 .
(12) This leads to the optimal performance rMSE opt ( ĥLMMSE ) = 1

E[pSNR] NrNt NcNs NrNt + Tr[R -1 ]Tr[R] (NrNt) 2 . ( 13
)
In this expression, the first term at the denominator

E[pSNR] NrNt

NcNs

NrNt corresponds to the expected inverse of the optimal LS performance, indicating that the LMMSE should exhibit the same behaviour for massive MIMO systems. However, the second term

Tr[R -1 ]Tr[R] (NrNt) 2
may compensate for the first one. This term is lower-bounded by one (if all the eigenvalues of R are equal), and grows as the eigenvalues distribution becomes more uneven (if R is closer to be singular). Simply put, the distribution of the eigenvalues of R determines to which extent knowing it is informative for channel estimation. The covariance being estimated in practice with previous channel estimates, the smaller the used time interval, the closer to singular R (since the user stayed in a smaller region), so that it gives more information and the optimal rMSE is reduced.

B. Physical channel estimators

As seen in the previous subsection, the LS and LMMSE estimators have been theoretically studied, and their performance is well understood. In contrast, estimators based on a physical model taking the form of ( 5) have not been investigated as deeply. This may be due to the difficulty to study jointly the bias and variance terms of the rMSE. To overcome this difficulty, an oracle estimator is introduced here which allows to separate the analysis of the bias from that of the variance. This oracle should be seen as a tool to study theoretically estimators based on a physical model, and its relevance with respect to practical estimators is assessed in section VI. Oracle definition. First, one can notice that for an estimator ĥ ∈ M p , rMSE( ĥ) ≥ h-proj Mp (h)

2 2 h 2 2
.

where proj Mp (u) argmin a∈Mp u-a . Then, the oracle estimator ĥoracle is defined as the unbiased efficient estimator of proj Mp (h). In other words, the oracle property amounts to consider that Mproj Mp (h)+n is observed instead of Mh+n, and thus ignore the part of the channel that is orthogonal to the model set M p . Mathematically speaking, using the oracle allows to study separately the bias and variance terms of the rMSE, as done below. Intuitively, the variance should increase with p whereas the bias should decrease, since p represents the complexity/expressiveness of the model.

Variance. The oracle being efficient means that

Tr[cov( ĥoracle )] = CRB
where CRB is the Cramér-Rao lower bound [START_REF] Calyampudi | Information and the accuracy attainable in the estimation of statistical parameters[END_REF], [21]. This bound has been derived in [START_REF] Le | Parametric channel estimation for massive MIMO[END_REF], leading to the following bound for the variance term of [START_REF] Akbar | Deconstructing multiantenna fading channels[END_REF].

Theorem 1. The variance of the oracle estimator is given by

Tr[cov( ĥoracle )] h 2 2 ≥ 3p pSNR . (14)
The important feature of this result is the fact that the bound on the variance is proportional to the number of virtual paths p (see [START_REF] Le | Parametric channel estimation for massive MIMO[END_REF] for the proof of the theorem). Moreover, denoting

proj Mp (h) = p i=1 b i e t ( --→ w t,i ) * ⊗e r ( --→ w r,i ), ( 15
)
the inequality is replaced by an equality in the theorem if the observation matrix M = X T ⊗W H obeys the two conditions

span p i=1 e t ( --→ w t,i ), ∂et( --→ wt,i) ∂ηt,i , ∂et( --→ wt,i) ∂ψt,i ⊂ im(X), span p i=1 e r ( --→ w r,i ), ∂er( --→ wr,i) ∂ηr,i , ∂er( --→ wr,i) ∂ψr,i ⊂ im(W).
Note that the variance is in that case directly proportional to the number of parameters to estimate. Moreover the result is pretty general. Indeed, for example if the receiver (or transmitter) has only one antenna, the theorem remains valid with the right-hand side of (14) becoming 2p pSNR , since the two parameters corresponding to the DoA (or DoD) of each virtual path disappear, leading to only four parameters per path.

Bias. The oracle being unbiased with respect to proj Mp (h) means that E[ ĥoracle ] = proj Mp (h). In that case, the bias term of (7) becomes

h-E[ ĥoracle ] 2 2 h 2 2 = h-proj Mp (h) 2 2 h 2 2 . ( 16
)
The study of this bias term is the main contribution of this paper. It is carried out in the next section.

IV. BIAS OF PHYSICAL CHANNEL ESTIMATORS Computing the bias of the oracle estimator ĥoracle defined above amounts to compute the projection proj Mp (h). Unfortunately, even considering a discretized set of candidates DoAs and DoDs, this problem (which then becomes a sparse approximation problem) is NP-hard [START_REF] Tropp | Computational methods for sparse solution of linear inverse problems[END_REF]. The projection can be approximated numerically using sparse recovery methods, as will be done in section VI. However the objective of this section is to study theoretically the bias and give an interpretable upper bound. Assumptions. Let us consider in this section the simpler case in which the transmitter or receiver has only one antenna (N t = 1 or N r = 1). In that case, the physical channel h and its estimator ĥ ∈ M p are expressed

h = P i=1 c i e( - → u i ), ĥ = p k=1 d k e( - → v k ),
where the indices denoting transmitter or receiver have been dropped since they are useless (they will be dropped in each notation defined in section II-B used in this section). This assumption allows to lighten considerably the notations in the following development. However, the general method used to bound the bias remains valid in the general model given by ( 4) and [START_REF] Björnson | Massive mimo networks: Spectral, energy, and hardware efficiency[END_REF].

Bound on the bias. Starting from [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF], the inequality

h-E[ ĥoracle ] 2 ≤ h-ĥ 2
is true for any ĥ ∈ M p (by definition of the projection). In particular, it is true when using the following suboptimal strategy to define ĥ: 1) Partition the set of P physical paths into p subsets R k ,k = 1,...,p, of paths with close directions. 2) Assign one virtual path of direction -→ v k to each subset R k .

3) Set d k as the optimal coefficient to approximate the paths in each subset R k by

d k e( - → v k ), given the direction - → v k .
The two first steps of this strategy amount to write

h-ĥ = p k=1 i∈R k c i e( - → u i )-d k e( - → v k ) ,
and the last one yields

d k = i∈R k c i e( - → v k ) H e( - → u i ) (this is the coefficient of the orthogonal projection onto e( - → v k )). In that case, h-ĥ 2 = p k=1 i∈R k c i e( - → u i )-e( - → v k ) H e( - → u i )e( - → v k ) 2 ≤ p k=1 i∈R k c i e( - → u i )-e( - → v k ) H e( - → u i )e( - → v k ) 2 = p k=1 i∈R k |c i | 1-|e( - → v k ) H e( - → u i )| 2 ,
where the second line is obtained by applying the triangle inequality twice and the third by simply developing the norm in the second one. At this point, the bias is bounded by a sum over all physical paths, each term being a quantity 1-|e( -→ v k ) H e( -→ u i )| 2 measuring the non-collinearity of the steering vectors corresponding to the physical path and to the associated virtual path, weighted by the modulus of its physical coefficient |c i |. To go further, let us analyze in details the scalar product between steering vectors appearing in the previous inequality. Its expression is given by

e( - → v k ) H e( - → u i ) = 1 N N n=1 e -j 2π λ -→ an.( - → ui--→ v k ) .
Using the infinite series representation of the exponential function leads to the following result.

Lemma 1. provided - → u i -- → v k 2 < 1 √ 2π -→ an λ 2
,∀n,

1-|e( - → v k ) H e( - → u i )| 2 ≤ 2π - → u i -- → v k 2 1 N N n=1 -→ a n λ 2 2 cos 2 -→ a n ,( - → u i -- → v k ) .
The proof is given in appendix A. If the condition of the lemma is fulfilled by each physical path (if there are enough virtual paths with appropriate directions), this yields

h-E[ ĥoracle ] 2 ≤ p k=1 i∈R k 2π|c i | - → u i -- → v k 2 1 N N n=1 -→ a n λ 2 2 cos 2 -→ a n ,( - → u i -- → v k ) .
Interpretations. First, notice that the part under the square root represents the angular sensitivity of the antenna array: it is bounded by the quantity κ(A)

1 N N n=1
-→ an λ 2 2 which depends only on the antenna array geometry (larger antenna arrays exhibit more angular sensitivity). Second, let us assess the dependency of the bound to p. To do so, assume that the assignment of physical paths to virtual paths is done by partitioning a given portion of the sphere in p regions and assigning all the physical paths whose direction falls into the kth region to the kth virtual path, whose direction is at the center of the region. That way, the bound is obviously decreasing when p increases, since adding a virtual path reduces the distance -→ u i --→ v k 2 for paths of the new region. The rate at which it decreases depends on the physical directions distribution. For example, in the worst-case scenario of uniformly distributed paths over a given portion of the sphere and regions of equal size, each of the p regions gathers of the order of 1 p physical paths and the typical distance to the center of each region is of the order of 1 √ p . This leads to a bound of the order of p× 1 p × 1 √ p = 1 √ p and consequently to a bias in [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF] of the order of 1 p at worst. The bias is evaluated numerically for realistic paths distributions in section VI. Finally, this bound lends itself to nice interpretations, but may be loose. A way of improving the bound is discussed in appendix B.

V. IMPACT ON THE DATA RATE

In order to assess the influence of channel estimation on the system performance, the rMSE in itself is of little use. Instead, the data rate loss caused by the estimation error (which leads to imperfect precoding) should be evaluated. To do so, assuming a transmission model of the form z = Hx + n after channel estimation has been performed and optimal nearest neighbour decoding [START_REF] Ayach | Spatially sparse precoding in millimeter wave mimo systems[END_REF], the channel capacity

C = log 2 det Id+E nn H -1 HE xx H H H
should be studied. Given the noise covariance E nn H = σ 2 Id and the channel singular value decomposition (SVD) H = UΛV H , the optimal signal covariance is given by E xx H opt = VDV H , where the diagonal matrix D corresponds to the optimal power allocation (computed by water-filling). Let us assume here that the channel is fixed but unknown to the system, which uses an estimate Ĥ instead (see [START_REF] Chul | Design and analysis of mimo spatial multiplexing systems with quantized feedback[END_REF] for a similar model). Moreover, given the SVD Ĥ = Û Λ VH , assume that the signal covariance is taken according to the imperfect estimate as E xx H = V D VH , leading to a mismatched precoder.

In that case, C = C opt -C loss , with [START_REF] Chul | Design and analysis of mimo spatial multiplexing systems with quantized feedback[END_REF] for more details). This expression is difficult to study in the general case. However, it simplifies if the receiver has only one antenna (N r = 1), to

C opt = log 2 det Id+ 1 σ 2 ΛDΛ H and C loss = -log 2 det Id+A Λ V H V D VH V-D Λ with A = 1 σ 2 Id+ 1 σ 2 ΛDΛ H -1 (see
C opt = log 2 (1+pSNR) and C loss = C opt -log 2 1+ | h, ĥ | 2 h 2 2 ĥ 2 2 pSNR . (17)
This capacity loss can be linked to the rMSE of channel estimation by the following bound:

C loss ≤ C opt -log 2 1+pSNR 1-rMSE( ĥ) 2-rMSE( ĥ) , (18)
proven in appendix C. Note that the bound is valid as soon as rMSE ≤ 1. This bound gives a worst-case data rate loss depending only on the rMSE and the pSNR. It is used in section VI to assess this loss for realistic channels.

VI. EXPERIMENTS

The objective of this section is to assess empirically the mathematical developments of the previous sections by: determining if using the oracle estimator makes sense, comparing physical channel estimators to others and quantifying the data rate loss caused by physical channel estimation. All experiments performed here are done using realistic channels generated with help of the NYUSIM channel simulator [START_REF] Sun | A novel millimeter-wave channel simulator and applications for 5g wireless communications[END_REF], in a millimeter wave massive MIMO downlink context. In particular, the frequency is set to f = 28 GHz and the distance between transmitter and receiver to d = 30 m to obtain the DoDs, DoAs, gains and phases of each path. The channel matrix is then obtained from (4) (with the total number of physical paths P between fifty and a hundred), considering a square uniform planar array (UPA) with half-wavelength separated antennas with N t = 64 (unless otherwise stated) at the transmitter and a single antenna receiver (N r = 1). All results shown here are averages over one hundred channel realizations. Relevance of the oracle. First, the goal is to determine whether or not the estimator ĥoracle has a behavior close to practical estimators. To do so, the bias of the oracle ( 16) is numerically evaluated by computing proj Mp (h) with help of the orthogonal matching pursuit algorithm (OMP) [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF], called directly on h (without noise). The variance of ĥoracle is taken directly as ( 14) with an equality (assuming optimal observations), and its rMSE is obtained by summing bias and variance. The oracle is compared to a practical estimator taking as input noisy observations of h, as in (3) with an optimal observation matrix, using OMP to get channel estimates. Results in terms of rMSE with respect to the number of virtual paths p for various pSNRs are shown on figure 1, in which the bias and variances of the oracle are also plotted (in black) to enhance clarity. First, one can notice that the oracle and the practical estimator exhibit a similar bias variance tradeoff, especially at high pSNR (note that for 64 antennas and no precoding, pSNR = SNR+18 dB). This allows to confirm that the oracle estimator is relevant, because it allows a separate study of the bias and variance while providing accurate performance predictions. Note that for a high p, the OMP estimator leads to a lower rMSE than the oracle, this is because it yields biased estimates with respect to the projection onto the model (thus the bound on the variance does not apply). Finally, note that the optimal number of virtual paths p yielding the smallest rMSE is very small compared to the number of physical paths P , (no more than a dozen against more than fifty), showing empirically that physical paths can indeed be merged into fewer virtual paths with little accuracy loss. Comparison with LMMSE. Let us now compare the oracle to other estimators presented in section III-A. Actually, the comparison is performed here with the LMMSE only, the comparison with the LS estimator not being shown for readability and brevity reasons. This is justified by the fact that the LS is always worse than the LMMSE in terms of optimal rMSE, but the same trends are observed, due to the similar expressions of the optimal performance (10) and ( 13), leading to the same conclusions. The performance of LMMSE is taken according to [START_REF] Jaeckel | Quadriga: A 3-d multi-cell channel model with time evolution for enabling virtual field trials[END_REF], assuming the covariance matrix R is perfectly known, and taking Tr

[R -1 ]Tr[R] (NrNt) 2
= 2 (this choice is rather arbitrary, but the point here is to show the relative behavior of the methods, assessing this quantity in practice would of course be valuable but is out of the scope of the paper). Results are shown on figure 2 for various number of transmit antennas. One important thing to notice is that especially at high pSNR, having a lot of transmit antennas degrades the performance of LMMSE more than it degrades the performance of the oracle estimator. This indicates that for massive MIMO with a lot of antennas, physical channel estimation methods may be better suited than LMMSE. Second, one can see that the optimal number of virtual paths to consider for the oracle estimator grows with the number of antennas. This behavior was expected from the bound on the bias given in section IV, since the quantity κ(A) grows with N t (the angular resolution increases). Finally, note that even though the curves indicate that the estimation task is harder for a greater number of antennas at constant pSNR, high pSNRs are easier to attain with a lot of antennas since according to (1), pSNR is proportional to N t . Impact on data rate. Finally, let us quantify the data rate loss caused by the estimation error when using a physical model, with respect to the number of considered virtual paths. To do so, the bound ( 17) is used with the rMSE of the oracle estimator, the results being shown on figure 3, where the relative capacity is defined as

Copt-Closs Copt

. These curves show it is sufficient to estimate very few virtual paths to guarantee a data rate close to the optimal one (obtained if the channel is perfectly known) in the tested millimeter wave context. For example, with a pSNR of 20 dB, only four virtual paths are sufficient to attain 94% of the optimal capacity. The higher the pSNR, the higher the optimal (in terms of capacity) number of virtual paths to estimate, which makes sense, since more physical paths are then above the noise level.

VII. CONCLUSION

In this paper, the performance of MIMO channel estimators using a physical model was theoretically studied. To do so, an appropriate oracle estimator was defined which allowed to study separately the bias and variance terms of the rMSE. The variance term having been already studied in a previous work [START_REF] Le | Parametric channel estimation for massive MIMO[END_REF], the present paper focused on the bias term, which was bounded by an interpretable quantity depending both on the propagation properties of the channel and on the antenna arrays geometries. Moreover, the data rate loss caused by the channel estimation error was bounded by a quantity depending only on the rMSE and the pSNR. Then, the defined oracle, whose relevance was empirically confirmed on realistic channels in a millimeter wave massive MIMO context, was compared to the classical LMMSE estimator, showing that channel estimators based on a physical model may be better suited to systems with a large number of antennas. Finally, it was shown empirically that estimating very few virtual paths is sufficient to guarantee a data rate close to the optimum in the aforementioned context.

In the future, it would be of great interest to perform extensive comparisons between physical channel estimators and other types of estimators such as the LMMSE, with the objective to see emerge different characteristic settings in which one or the other type of estimator is more adapted. Such a study would be very useful for system design. On a more technical side, improving the bound on the bias proposed here should be possible, as mentioned in appendix B.

Fig. 1 :

 1 Fig. 1: Oracle compared to the OMP algorithm for various pSNRs, with respect to p.

Fig. 2 :

 2 Fig. 2: Comparison of the oracle with the LMMSE for various numbers of transmit antennas and pSNRs, with respect to p.

Fig. 3 :

 3 Fig. 3: Relative capacity with respect to p for various pSNRs.

APPENDIX A. Proof of lemma 1

 1 Using the series representation of the exponential e x = +∞ l=0x l l! , the scalar product is expressed e( -→ v k ) H e( -→ u n = 0 since the antennas are located with respect to the centroid of the array, the term l = 1 is null, leading toe( -→ v k ) H e( -→ u i ) = 1+ +∞ l=2 n .( -→ u i --→ v k ) l .

  T , A * and A H respectively. The trace of a linear transformation represented by A is denoted Tr(A). For matrices A and B, A ≥ B means that A -B is positive semidefinite. The linear span of a set of vectors A is denoted: span(A). The Kronecker product and vectorization operators are denoted by ⊗ and vec(•) respectively. The identity matrix is denoted by Id.

CN (µ,Σ) denotes the standard complex gaussian distribution with mean µ and covariance Σ. E(•) denotes expectation and cov(•) the covariance of its argument.

= (M H M) -1 M H y. (8)

Both the real and imaginary parts of this series are alternating series. The terms are decreasing in magnitude if

,∀n.

This condition is fulfilled as soon as

,∀n, (taking l = 0 and applying the Cauchy-Schwarz inequality). In that case, Re e( -

Developing the scalar product proves the lemma.

B. Improving the bound

The bound of this paper is obtained applying the triangle inequality twice in a row, which may lead to a loose bound. Instead, one could apply it only once to get

At this point, one can notice that each term of the sum, when squared, corresponds to a quadratic form

Studying the properties of the matrix Q may lead to a tighter bound, but is not guaranteed to lead to interpretable results.

C. Proof of (18) Starting from [START_REF] Alkhateeb | Channel estimation and hybrid precoding for millimeter wave cellular systems[END_REF], one can first notice that