<figure><img src="image_1.png" /></figure>

HAL
open science

# A Fully Bayesian Extension to FEMU for Identification of Spatially Varying Elastic Properties from Digital Image and Volume Correlation Measurements

Armand Touminet, Sabine Cantournet, Victor Fabre, Pierre Kerfriden

► To cite this version:

Armand Touminet, Sabine Cantournet, Victor Fabre, Pierre Kerfriden. A Fully Bayesian Extension to FEMU for Identification of Spatially Varying Elastic Properties from Digital Image and Volume Correlation Measurements. 2024. [hal-04587210](https://hal.science/hal-04587210)

HAL Id: hal-04587210

https://hal.science/hal-04587210

Preprint submitted on 27 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à 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.

Public Domain
# A Fully Bayesian Extension to FEMU for Identification of Spatially Varying Elastic Properties from Digital Image and Volume Correlation Measurements

Armand Touminet¹,², Sabine Cantournet¹, Victor Fabre², and Pierre Kerfriden¹

¹Mines Paris, Centre des Matériaux, 91100, France

²Hutchinson, Centre de Recherche et Innovation, Chalette-sur-Loing, 45120, France

## Abstract

We present a fully Bayesian framework for identifying spatially varying elastic parameters and their covariance properties using noisy displacement observations obtained with DIC or DVC trials. Our method is a generalization of identification procedures such as FEMU or I-DIC to materials with spatially varying properties and stochastic mesostructures. The identified variables can be used to generate new samples with similar covariance properties. We formulate the approach as a hierarchical Bayesian PDE-constrained inverse problem and MAP estimates are obtained through gradient based optimization. We resort to an adjoint based formulation and leverage automatic differentiation to derive the parameter sensitivities. We show how modelling unknown parameters with Gaussian Random Fields leads to a natural Bayesian regularization and develop the use of Whittle-Matérn priors. Hyperparameter estimation is discussed, and we propose an empirical Bayes approach to avoid numerical shortcomings related to a standard hierarchical model. A set of numerical examples is presented to assess the performance of the proposed method, based on synthetic data generated through Matérn Random fields. In particular, we show how data noise is naturally modelled by the Bayesian formulation and impacts spatial covariance of identified parameters.

## 1 Introduction

Full field identification procedures from in-situ observations such as Digital Image Correlation (DIC), stereo-DIC, or Digital Volume Correlation (DVC) have widely been developed in the past decades to calibrate constitutive laws of a wide range of materials. Examples of such identification procedures include the Finite Elements Updating Method (FEMU), Integrated DIC (I-DIC) [Leclerc et al., 2009] or Mechanical Image Correlation (MIC) [Réthoré, 2010]. In these methods, a digital twin of a mechanical trial is developed, and the simulation to experiment gap is iteratively minimized by an optimization algorithm [Avril et al., 2008]. Usually, the parameters of interest are a set of scalar variables from a complex and non-linear constitutive law [Réthoré, 2010, Leclerc et al., 2009, Neggers et al., 2017]. Typically, the identified parameters are assumed to be constant across the specimen of interest, which is suitable to study homogeneous materials, or heterogeneous materials which mechanical properties present no significant variability at the macroscopic scale. In this work, we aim at extending those identification procedures to materials with a significant spatial variability of mechanical properties or with stochastic micro- or mesostructures. For instance, discontinuous composites include a large range of materials with such properties, including compression moulded Discontinuous Long Fiber composites (DLFs) with prepreg based Sheet Moulding Compounds (SMCs) or thermoplastic Tow Based Discontinuous Composites (TBDCs), and injection moulded short fiber composites. In these materials, a strong spatial variability is observed due to random fiber bundle layout, and fiber realignment during moulding [Martulli et al., 2019], which are of great importance for modelling such materials [Alves et al., 2023]. Other examples include concrete, where strength has been shown to be spatially variable [Geyer et al., 2023], or functionally graded materials [Watanabe et al., 2008], where mechanical parameters are purposely continuously varying across a given specimen. Although the authors aim at applying the present developments to DLF composites, this paper remains purely
# A Fully Bayesian Extension to FEMU for Identification of Spatially Varying Elastic Properties from Digital Image and Volume Correlation Measurements

Armand Touminet¹,², Sabine Cantournet¹, Victor Fabre², and Pierre Kerfriden¹

¹Mines Paris, Centre des Matériaux, 91100, France

²Hutchinson, Centre de Recherche et Innovation, Chalette-sur-Loing, 45120, France

## Abstract

We present a fully Bayesian framework for identifying spatially varying elastic parameters and their covariance properties using noisy displacement observations obtained with DIC or DVC trials. Our method is a generalization of identification procedures such as FEMU or I-DIC to materials with spatially varying properties and stochastic mesostructures. The identified variables can be used to generate new samples with similar covariance properties. We formulate the approach as a hierarchical Bayesian PDE-constrained inverse problem and MAP estimates are obtained through gradient based optimization. We resort to an adjoint based formulation and leverage automatic differentiation to derive the parameter sensitivities. We show how modelling unknown parameters with Gaussian Random Fields leads to a natural Bayesian regularization and develop the use of Whittle-Matérn priors. Hyperparameter estimation is discussed, and we propose an empirical Bayes approach to avoid numerical shortcomings related to a standard hierarchical model. A set of numerical examples is presented to assess the performance of the proposed method, based on synthetic data generated through Matérn Random fields. In particular, we show how data noise is naturally modelled by the Bayesian formulation and impacts spatial covariance of identified parameters.

## 1 Introduction

Full field identification procedures from in-situ observations such as Digital Image Correlation (DIC), stereo-DIC, or Digital Volume Correlation (DVC) have widely been developed in the past decades to calibrate constitutive laws of a wide range of materials. Examples of such identification procedures include the Finite Elements Updating Method (FEMU), Integrated DIC (I-DIC) [Leclerc et al., 2009] or Mechanical Image Correlation (MIC) [Réthoré, 2010]. In these methods, a digital twin of a mechanical trial is developed, and the simulation to experiment gap is iteratively minimized by an optimization algorithm [Avril et al., 2008]. Usually, the parameters of interest are a set of scalar variables from a complex and non-linear constitutive law [Réthoré, 2010, Leclerc et al., 2009, Neggers et al., 2017]. Typically, the identified parameters are assumed to be constant across the specimen of interest, which is suitable to study homogeneous materials, or heterogeneous materials which mechanical properties present no significant variability at the macroscopic scale. In this work, we aim at extending those identification procedures to materials with a significant spatial variability of mechanical properties or with stochastic micro- or mesostructures. For instance, discontinuous composites include a large range of materials with such properties, including compression moulded Discontinuous Long Fiber composites (DLFs) with prepreg based Sheet Moulding Compounds (SMCs) or thermoplastic Tow Based Discontinuous Composites (TBDCs), and injection moulded short fiber composites. In these materials, a strong spatial variability is observed due to random fiber bundle layout, and fiber realignment during moulding [Martulli et al., 2019], which are of great importance for modelling such materials [Alves et al., 2023]. Other examples include concrete, where strength has been shown to be spatially variable [Geyer et al., 2023], or functionally graded materials [Watanabe et al., 2008], where mechanical parameters are purposely continuously varying across a given specimen. Although the authors aim at applying the present developments to DLF composites, this paper remains purely
numerical and material agnostic, and application to a specific material is deferred to subsequent works.

Extending identification procedures to spatially varying parameters introduces numerous numerical challenges to overcome. First, a proper spatial discretization of the parameters of interests is introduced, resulting in a high number of variables to optimize. In typical FEMU or I-DIC methods, the simulation to experiment gap is minimized with a quasi-Newton algorithm, and the gradient vector is usually approximated with finite differences, and requires solving multiple linear systems for each evaluation of the objective function. With a fine spatial discretization for the variable to optimize, this process becomes numerically intractable, and in this paper we leverage an adjoint based approach to compute the gradient vector efficiently, independently of the number of variables to optimize. Secondly, parameter identification from displacement observations is known to be an ill-posed inverse problem, and requires spatial regularization to avoid local minima of the objective function and to avoid overfitting observation noise. The use of spatial regularization, such as the well known Tikhonov regularizer family, introduces so-called regularization parameters that are often cumbersome to select, although some selection methods exist, such as those based on the Morozov principle [Nair, 2009], or the L-curve approach [Hansen, 1992]. The present framework is formulated as a Bayesian inverse problem, meaning that both parameters of interest and observation noise are modelled as random variables, and we derive Maximum A Posteriori (MAP) estimates for parameter identification. We use Gaussian priors to model the sought parameters and leverage covariance matrices derived from Matérn class kernels, which have the advantage of having statistically and geometrically interpretable parameters, namely the prior variance and length scale. More specifically, we use stationary Whittle-Matérn priors, which have been shown to be very numerically efficient to implement by following the SPDE approach [Lindgren et al., 2011], which leads to a sparse approximation of the precision matrix on unstructured grids. For hyperparameter selection, we propose a Bayesian approach by introducing a hierarchical model, which is in turn approximated following a Parametric Empirical Bayes (PEB) approach. Three choices of estimators are investigated for hyperparameter update, namely a maximum likelihood estimator, a maximum a posteriori estimator, and an empirical estimator based on variogram regression. We discuss the numerical issues and shortcomings of each approach and link numerical results to the existing literature. In particular, we show that a maximum likelihood update strategy alongside a Whittle-Matérn prior leads to an unstable numerical scheme, and propose a numerically robust alternative. We assess the performance of the PEB approach on synthetic samples and show the impact of data noise on estimated parameters and hyperparameters.

Although identification from surface (DIC) and volume (DVC) measurements typically constitute separate fields of research, we abstract the underlying method used to obtain displacement fields from image correlation, and treat noisy surface or volume displacement fields as input to our framework, that is presented hereafter to handle both types of measurements.

For implementation, we leverage automatic differentiation to derive the adjoint operator automatically and independently of the chosen stiffness model and show examples of parameter identification on synthetic data with a selected set of stiffness models. In the present paper, we develop the proposed method in a stiffness model agnostic way. However, for simplicity, we assume that the stiffness model is linear, though it could be adapted to non-linear constitutive laws with minor modifications only.

## 2 Related Works and Main Contributions

Following the developments of Digital Image Correlation (see, e.g., Sutton et al. [2009] for details about common image correlation techniques), many kinds of parameters identification methods were developed, among which we distinguish non-integrated procedures, such as FEMU, in which the simulation to experiment gap is formulated in terms of a displacement misfit, and integrated procedures, such as I-DIC, in which an image residual is considered to perform the identification, removing the need to obtain correlated displacement fields beforehand. Many other methods exist, and for a comprehensive review, we refer to [Avril et al., 2008]. Although the differences between all of these methods are not important for the present paper, it was shown that data noise can impact the identification results significantly, especially on setups with a low signal-to-noise ratio [Ruybalid et al.]. An interesting variation is the regularized FEMU (FEMU-R) method, in which the displacement gap
numerical and material agnostic, and application to a specific material is deferred to subsequent works.

Extending identification procedures to spatially varying parameters introduces numerous numerical challenges to overcome. First, a proper spatial discretization of the parameters of interests is introduced, resulting in a high number of variables to optimize. In typical FEMU or I-DIC methods, the simulation to experiment gap is minimized with a quasi-Newton algorithm, and the gradient vector is usually approximated with finite differences, and requires solving multiple linear systems for each evaluation of the objective function. With a fine spatial discretization for the variable to optimize, this process becomes numerically intractable, and in this paper we leverage an adjoint based approach to compute the gradient vector efficiently, independently of the number of variables to optimize. Secondly, parameter identification from displacement observations is known to be an ill-posed inverse problem, and requires spatial regularization to avoid local minima of the objective function and to avoid overfitting observation noise. The use of spatial regularization, such as the well known Tikhonov regularizer family, introduces so-called regularization parameters that are often cumbersome to select, although some selection methods exist, such as those based on the Morozov principle [Nair, 2009], or the L-curve approach [Hansen, 1992]. The present framework is formulated as a Bayesian inverse problem, meaning that both parameters of interest and observation noise are modelled as random variables, and we derive Maximum A Posteriori (MAP) estimates for parameter identification. We use Gaussian priors to model the sought parameters and leverage covariance matrices derived from Matérn class kernels, which have the advantage of having statistically and geometrically interpretable parameters, namely the prior variance and length scale. More specifically, we use stationary Whittle-Matérn priors, which have been shown to be very numerically efficient to implement by following the SPDE approach [Lindgren et al., 2011], which leads to a sparse approximation of the precision matrix on unstructured grids. For hyperparameter selection, we propose a Bayesian approach by introducing a hierarchical model, which is in turn approximated following a Parametric Empirical Bayes (PEB) approach. Three choices of estimators are investigated for hyperparameter update, namely a maximum likelihood estimator, a maximum a posteriori estimator, and an empirical estimator based on variogram regression. We discuss the numerical issues and shortcomings of each approach and link numerical results to the existing literature. In particular, we show that a maximum likelihood update strategy alongside a Whittle-Matérn prior leads to an unstable numerical scheme, and propose a numerically robust alternative. We assess the performance of the PEB approach on synthetic samples and show the impact of data noise on estimated parameters and hyperparameters.

Although identification from surface (DIC) and volume (DVC) measurements typically constitute separate fields of research, we abstract the underlying method used to obtain displacement fields from image correlation, and treat noisy surface or volume displacement fields as input to our framework, that is presented hereafter to handle both types of measurements.

For implementation, we leverage automatic differentiation to derive the adjoint operator automatically and independently of the chosen stiffness model and show examples of parameter identification on synthetic data with a selected set of stiffness models. In the present paper, we develop the proposed method in a stiffness model agnostic way. However, for simplicity, we assume that the stiffness model is linear, though it could be adapted to non-linear constitutive laws with minor modifications only.

## 2 Related Works and Main Contributions

Following the developments of Digital Image Correlation (see, e.g., Sutton et al. [2009] for details about common image correlation techniques), many kinds of parameters identification methods were developed, among which we distinguish non-integrated procedures, such as FEMU, in which the simulation to experiment gap is formulated in terms of a displacement misfit, and integrated procedures, such as I-DIC, in which an image residual is considered to perform the identification, removing the need to obtain correlated displacement fields beforehand. Many other methods exist, and for a comprehensive review, we refer to [Avril et al., 2008]. Although the differences between all of these methods are not important for the present paper, it was shown that data noise can impact the identification results significantly, especially on setups with a low signal-to-noise ratio [Ruybalid et al.]. An interesting variation is the regularized FEMU (FEMU-R) method, in which the displacement gap
is weighted by a displacement noise correlation matrix [Leclerc et al., 2009, Réthoré, 2010], leading to the same data fidelity term as typically used in MAP estimates, hinting towards a Bayesian approach, although it was never formulated as such. In the aforementioned works, the cost function is commonly minimized by the means of a quasi-Newton algorithm with a Gauss-Newton approximation for the hessian, and the gradient vector is usually approximated by finite differences.

The adjoint method has been used to solve inverse problems in linear elasticity, in a context unrelated to DIC or DVC, in e.g. [Seidl et al., 2019], to identify a piecewise constant shear modulus, and also resorts to total variation (TV) regularization to enforce piecewise constant solutions. In [Bonnet and Constantinescu, 2005], the adjoint method is briefly presented to identify an unknown Poisson ratio distribution, however the questions of regularization and data noise are not investigated thoroughly. Also, in these last references, the spatial covariance structure of the unknown parameter fields are not investigated.

The use of Gaussian random fields for modelling spatial data has been developed and applied to a wide range of data, including soil properties [Ricketts et al., 2023], concrete failure parameters [Hai and Lyu, 2023], or even laminate composite tensile modulus [Sriramula and Chryssanthopoulos, 2009]. In such cases, random fields models can be used to generate new realisations of the underlying spatial process, and are a good alternative to Monte-Carlo simulations that can be cumbersome. For this purpose, a proper selection of covariance parameters is key, as it can highly influence the simulation results.

In the Bayesian inversion literature, identification of a spatially varying unknown parameter is a common problem, and the use of so-called Bayesian priors enforce statistical information about the solution. For a review of common methods and algorithms, we refer to [Dashti and Stuart, 2015]. A specific challenge arising when the unknown is discretized on an unstructured mesh, which is common when dealing with PDE-constrained inverse problems, is that evaluating the prior term typically involves assembling a dense matrix, which has a quadratic complexity in time and memory and is often intractable in large scale problems. This was partially solved in [Lindgren et al., 2011] by proposing a sparse approximation of precision matrices arising from some cases of Matérn class correlation kernels. In later works [Roininen et al., 2014], such priors were used to identify unknown PDE parameters through MAP estimates. In more recent works [Roininen et al., 2019], the case where Matérn prior hyperparameters are spatially varying was studied, and a hierarchical model based on Cauchy hyper-priors was proposed to estimate these hyperparameters through Markov Chain Monte Carlo (MCMC) sampling. However, papers tackling approximation of Bayesian priors on unstructured meshes are largely unrelated to parameter estimation from DIC or DVC measurements.

In the present paper, we bridge the gap between all the previously mentioned topics to construct a comprehensive framework to model spatially varying elastic properties with Gaussian random fields, and perform parameter and hyperparameter identification from DIC and DVC measurements:

* We extend the FEMU method to a fully Bayesian formulation that takes into account all sources of measurement noise,

* An optimization strategy based on the adjoint method and automatic differentiation is presented to derive maximum a posteriori estimates,

* We leverage Whittle-Matérn priors to model the covariance structure of unknown parameters and obtain an efficient discretization scheme,

* A hierarchical model is presented for hyperparameter estimation, and a numerical method based on the PEB approach is presented and thoroughly discussed.

* The results can be used to generate new realisations of the spatially varying parameter, with estimated covariance parameters.

## 3 Notations

Throughout this article, and unless specified otherwise, zeroth, first, second and fourth order tensors and tensor valued functions are denoted by $a$, $\underline{a}$, $\underline{\underline{a}}$ and $\bar{a}$ respectively. The simple and double contraction are written $\cdot$ and $:$ respectively, and we employ the Einstein notation to express them in a Cartesian basis:

$$ \underline{a} \cdot \underline{b} = a_i b_j, \quad \underline{\underline{a}} : \underline{\underline{b}} = a_{ij} b_{ij}, \quad \bar{a} : \underline{\underline{b}} = a_{ijkl} b_{kl}. $$
is weighted by a displacement noise correlation matrix [Leclerc et al., 2009, Réthoré, 2010], leading to the same data fidelity term as typically used in MAP estimates, hinting towards a Bayesian approach, although it was never formulated as such. In the aforementioned works, the cost function is commonly minimized by the means of a quasi-Newton algorithm with a Gauss-Newton approximation for the hessian, and the gradient vector is usually approximated by finite differences.

The adjoint method has been used to solve inverse problems in linear elasticity, in a context unrelated to DIC or DVC, in e.g. [Seidl et al., 2019], to identify a piecewise constant shear modulus, and also resorts to total variation (TV) regularization to enforce piecewise constant solutions. In [Bonnet and Constantinescu, 2005], the adjoint method is briefly presented to identify an unknown Poisson ratio distribution, however the questions of regularization and data noise are not investigated thoroughly. Also, in these last references, the spatial covariance structure of the unknown parameter fields are not investigated.

The use of Gaussian random fields for modelling spatial data has been developed and applied to a wide range of data, including soil properties [Ricketts et al., 2023], concrete failure parameters [Hai and Lyu, 2023], or even laminate composite tensile modulus [Sriramula and Chryssanthopoulos, 2009]. In such cases, random fields models can be used to generate new realisations of the underlying spatial process, and are a good alternative to Monte-Carlo simulations that can be cumbersome. For this purpose, a proper selection of covariance parameters is key, as it can highly influence the simulation results.

In the Bayesian inversion literature, identification of a spatially varying unknown parameter is a common problem, and the use of so-called Bayesian priors enforce statistical information about the solution. For a review of common methods and algorithms, we refer to [Dashti and Stuart, 2015]. A specific challenge arising when the unknown is discretized on an unstructured mesh, which is common when dealing with PDE-constrained inverse problems, is that evaluating the prior term typically involves assembling a dense matrix, which has a quadratic complexity in time and memory and is often intractable in large scale problems. This was partially solved in [Lindgren et al., 2011] by proposing a sparse approximation of precision matrices arising from some cases of Matérn class correlation kernels. In later works [Roininen et al., 2014], such priors were used to identify unknown PDE parameters through MAP estimates. In more recent works [Roininen et al., 2019], the case where Matérn prior hyperparameters are spatially varying was studied, and a hierarchical model based on Cauchy hyper-priors was proposed to estimate these hyperparameters through Markov Chain Monte Carlo (MCMC) sampling. However, papers tackling approximation of Bayesian priors on unstructured meshes are largely unrelated to parameter estimation from DIC or DVC measurements.

In the present paper, we bridge the gap between all the previously mentioned topics to construct a comprehensive framework to model spatially varying elastic properties with Gaussian random fields, and perform parameter and hyperparameter identification from DIC and DVC measurements:

* We extend the FEMU method to a fully Bayesian formulation that takes into account all sources of measurement noise,

* An optimization strategy based on the adjoint method and automatic differentiation is presented to derive maximum a posteriori estimates,

* We leverage Whittle-Matérn priors to model the covariance structure of unknown parameters and obtain an efficient discretization scheme,

* A hierarchical model is presented for hyperparameter estimation, and a numerical method based on the PEB approach is presented and thoroughly discussed.

* The results can be used to generate new realisations of the spatially varying parameter, with estimated covariance parameters.

## 3 Notations

Throughout this article, and unless specified otherwise, zeroth, first, second and fourth order tensors and tensor valued functions are denoted by $a$, $\underline{a}$, $\underline{\underline{a}}$ and $\bar{a}$ respectively. The simple and double contraction are written $\cdot$ and $:$ respectively, and we employ the Einstein notation to express them in a Cartesian basis:

$$ \underline{a} \cdot \underline{b} = a_i b_j, \quad \underline{a} : \underline{b} = a_{ij} b_{ij}, \quad \bar{a} : \underline{b} = a_{ijkl} b_{kl}. $$
If $x$ is a continuous random variable, we note $D(x)$ the density of $x$. If $x$ follows a normal or uniform law, we respectively note $x \sim \mathcal{N}(x^0, \sigma_x^2)$ where $x^0$ is the mean and $\sigma_x^2$ the variance, and $x \sim U(I)$ where $I \subset \mathbb{R}$ is a bounded interval.

If $U$ and $V$ are two Banach spaces, and $F: U \to V$ is a differentiable function, the Gateaux derivative of $F$ at point $x \in U$ along the direction $\delta x \in U$ is denoted by $dF(x; \delta x)$. If $U_1, ..., U_n$ are Banach spaces, and $G: U_1 \times ... \times U_n \to V$ is differentiable function, we note $\partial_i G(x_1, ..., x_n; \delta x_i)$ its Gateaux derivative with the respect to its $i$-th variable, defined as

$$ \partial_i G(x_1, ..., x_n; \delta x_i) = \lim_{h \to 0} \frac{G(x_1, ..., x_i + h\delta x_i, ..., x_n) - G(x_1, ..., x_n)}{h} $$

If $U_1, U_2$ are of finite dimension, bilinear forms $a: U_1 \times U_2 \to \mathbb{R}$ are identified with their matrix representations, denoted with a bold upper case letter $\mathbf{A}$, and linear forms $a: U_1 \to \mathbb{R}$ with their vector representation, by a bold lower case letter $\mathbf{a}$. Given bases for $U_1$ and $U_2$, $\mathbf{a} \in \mathbb{R}^{\dim(U_1)}$ and $\mathbf{A} \in \mathbb{R}^{\dim(U_1) \times \dim(U_2)}$. Given an open and bounded domain $\Omega \subset \mathbb{R}^n$, $\partial \Omega = \bar{\Omega} \setminus \Omega$ denotes the boundary of $\Omega$, and we introduce the common Hilbert and Sobolev spaces used in variational analysis, namely $L^\infty(\Omega)$ the space of bounded functions over $\Omega$, $L^2(\Omega)$ the space of square integrable functions over $\Omega$ and $H^1(\Omega)$.

# 4 Stiffness Modelling with Gaussian Random Fields

Throughout this paper, we consider on open and bounded domain $\Omega \subset \mathbb{R}^3$ representing a specimen made of an abstract material with spatially varying elastic properties, in the sense that the stiffness tensor $\bar{C}(\underline{x})$ is spatially varying. We also assume that the stiffness tensor field can be parameterized with a set of $N_p$ scalar valued parameter fields $(p_1(\underline{x}), ..., p_{N_p}(\underline{x})) \in \mathcal{D} \subset \mathbb{R}^{N_p}$ (e.g. isotropy, transverse isotropy, etc.), where $\mathcal{D}$ is called the parameter space. The mapping $(p_1, ..., p_{N_p}) \mapsto \bar{C}(p_1, ..., p_{N_p})$ is hereafter referred to as the stiffness model, and is assumed to be well-defined, in the sense that $\bar{C}(\mathcal{D})$ is a set containing tensors that satisfy the usual symmetry assumptions in linear elasticity. We further assume that spatial variations are due to some unknown spatial random process, and model the stiffness parameters as Gaussian random fields:

$$ \forall 1 \le i \le N_p, \quad p_i \sim \mathcal{N}(p_i^0, c_i), \qquad (1) $$

where $p_i^0$ is known and $c_i$ is an unknown covariance function. By combining the stiffness model with a specified covariance model for the stiffness parameters, it is possible to sample realizations of a synthetic specimen that accounts for spatial variability and randomness of the underlying material. The goal of this paper is to infer estimates for stiffness parameters and covariance models from displacement measurements only. The inferred model can then be used to generate numerous specimens for simulations. The method is formulated as a Bayesian inverse problem. Section 5 presents how the inverse problem is formulated and derives MAP estimates for the unknown stiffness parameters and covariance hyperparameters, calibration of the covariance model is detailed in section 6, section 8 presents how the adjoint method is used to solve the inverse problem, and section 9 illustrates the proposed method with two numerical examples.

# 5 Problem Description and Bayesian Modelling

In this section, we present the steps to extend the FEMU method to identify a spatially varying material parameter, with a fully Bayesian approach. The developments presented throughout this paper can be adapted to the I-DIC method in the same way, and for completeness, details are given in appendix A.

## 5.1 PDE Modelling and Discretization

We consider the spatial domain $\Omega$ to represent a trial specimen subject to an (unknown) boundary displacement $u_D$ over $\Omega_D \subset \partial \Omega$, modelling some mechanical trial. We also consider an observation zone $\Omega^{\text{obs}} \subset \Omega$ modelling a domain where in-situ imaging is performed. For DIC, $\Omega^{\text{obs}} \subset \partial \Omega$
If $x$ is a continuous random variable, we note $D(x)$ the density of $x$. If $x$ follows a normal or uniform law, we respectively note $x \sim \mathcal{N}(x^0, \sigma_x^2)$ where $x^0$ is the mean and $\sigma_x^2$ the variance, and $x \sim U(I)$ where $I \subset \mathbb{R}$ is a bounded interval.

If $U$ and $V$ are two Banach spaces, and $F: U \to V$ is a differentiable function, the Gateaux derivative of $F$ at point $x \in U$ along the direction $\delta x \in U$ is denoted by $dF(x; \delta x)$. If $U_1, ..., U_n$ are Banach spaces, and $G: U_1 \times ... \times U_n \to V$ is differentiable function, we note $\partial_i G(x_1, ..., x_n; \delta x_i)$ its Gateaux derivative with the respect to its $i$-th variable, defined as

$$ \partial_i G(x_1, ..., x_n; \delta x_i) = \lim_{h \to 0} \frac{G(x_1, ..., x_i + h\delta x_i, ..., x_n) - G(x_1, ..., x_n)}{h} $$

If $U_1, U_2$ are of finite dimension, bilinear forms $a: U_1 \times U_2 \to \mathbb{R}$ are identified with their matrix representations, denoted with a bold upper case letter $\mathbf{A}$, and linear forms $a: U_1 \to \mathbb{R}$ with their vector representation, by a bold lower case letter $\mathbf{a}$. Given bases for $U_1$ and $U_2$, $\mathbf{a} \in \mathbb{R}^{\dim(U_1)}$ and $\mathbf{A} \in \mathbb{R}^{\dim(U_1) \times \dim(U_2)}$. Given an open and bounded domain $\Omega \subset \mathbb{R}^n$, $\partial \Omega = \bar{\Omega} \setminus \Omega$ denotes the boundary of $\Omega$, and we introduce the common Hilbert and Sobolev spaces used in variational analysis, namely $L^\infty(\Omega)$ the space of bounded functions over $\Omega$, $L^2(\Omega)$ the space of square integrable functions over $\Omega$ and $H^1(\Omega)$.

# 4 Stiffness Modelling with Gaussian Random Fields

Throughout this paper, we consider on open and bounded domain $\Omega \subset \mathbb{R}^3$ representing a specimen made of an abstract material with spatially varying elastic properties, in the sense that the stiffness tensor $\bar{C}(\underline{x})$ is spatially varying. We also assume that the stiffness tensor field can be parameterized with a set of $N_p$ scalar valued parameter fields $(p_1(\underline{x}), ..., p_{N_p}(\underline{x})) \in \mathcal{D} \subset \mathbb{R}^{N_p}$ (e.g. isotropy, transverse isotropy, etc.), where $\mathcal{D}$ is called the parameter space. The mapping $(p_1, ..., p_{N_p}) \mapsto \bar{C}(p_1, ..., p_{N_p})$ is hereafter referred to as the stiffness model, and is assumed to be well-defined, in the sense that $\bar{C}(\mathcal{D})$ is a set containing tensors that satisfy the usual symmetry assumptions in linear elasticity. We further assume that spatial variations are due to some unknown spatial random process, and model the stiffness parameters as Gaussian random fields:

$$ \forall 1 \le i \le N_p, \quad p_i \sim \mathcal{N}(p_i^0, c_i), \qquad (1) $$

where $p_i^0$ is known and $c_i$ is an unknown covariance function. By combining the stiffness model with a specified covariance model for the stiffness parameters, it is possible to sample realizations of a synthetic specimen that accounts for spatial variability and randomness of the underlying material. The goal of this paper is to infer estimates for stiffness parameters and covariance models from displacement measurements only. The inferred model can then be used to generate numerous specimens for simulations. The method is formulated as a Bayesian inverse problem. Section 5 presents how the inverse problem is formulated and derives MAP estimates for the unknown stiffness parameters and covariance hyperparameters, calibration of the covariance model is detailed in section 6, section 8 presents how the adjoint method is used to solve the inverse problem, and section 9 illustrates the proposed method with two numerical examples.

# 5 Problem Description and Bayesian Modelling

In this section, we present the steps to extend the FEMU method to identify a spatially varying material parameter, with a fully Bayesian approach. The developments presented throughout this paper can be adapted to the I-DIC method in the same way, and for completeness, details are given in appendix A.

## 5.1 PDE Modelling and Discretization

We consider the spatial domain $\Omega$ to represent a trial specimen subject to an (unknown) boundary displacement $u_D$ over $\Omega_D \subset \partial\Omega$, modelling some mechanical trial. We also consider an observation zone $\Omega^{\text{obs}} \subset \Omega$ modelling a domain where in-situ imaging is performed. For DIC, $\Omega^{\text{obs}} \subset \partial\Omega$
<figure><img src="image_1.png" /><figcaption>Figure 1: The two considered specimen: a thin dog bone geometry for uniaxial tensile test, and a thick biaxial geometry.</figcaption></figure>

represents the specimen surface captured by the camera, as for DVC, $\Omega^{\text{obs}}$ represents a volume subregion imaged through, e.g., x-ray tomography. We represent in figure 1 and 11 two specimen geometries where the different regions are annotated, illustrating both surface and volume observation zones. We model the results of image correlation by a set of pointwise, noisy displacement observations

$$ \underline{u}^{\text{obs}} = (\underline{u}^{\text{obs}}(\underline{x}_1), \dots, \underline{u}^{\text{obs}}(\underline{x}_{N^{\text{obs}}})) , \quad (2) $$

$$ \underline{u}^{\text{obs}}(\underline{x}_i) = \underline{u}(\underline{x}_i) + \underline{\xi}, \quad 1 \le i \le N^{\text{obs}}, \quad (3) $$

where $N^{\text{obs}}$ is the number of pointwise observations, and $\underline{u}^{\text{obs}}$ is a vector of size $3 \times N^{\text{obs}}$. Here, $\underline{u}$ denotes the true displacement field, and $\underline{\xi}$ some observation noise. The observed macroscopic force $F^{\text{obs}} \in \mathbb{R}$ models the measured force $F$, typically by a force sensor, that is also corrupt by some noise $\eta$, such that

$$ F^{\text{obs}} = F + \eta. \quad (4) $$

Without loss of generality, we assume that there is only a single force sensor. In other cases (e.g. biaxial traction), $F$ should be taken as vector valued. Let $\underline{y}^{\text{obs}}$ and $\underline{\epsilon}$ be vectors of size $3 \times N^{\text{obs}} + 1$ collecting all measured values and noise values respectively:

$$ \underline{y}^{\text{obs}} = (\underline{u}^{\text{obs}}(\underline{x}_1), \dots, \underline{u}^{\text{obs}}(\underline{x}_{N^{\text{obs}}}), F^{\text{obs}}), \quad (5) $$

$$ \underline{\epsilon} = (\underline{\xi}(\underline{x}_1), \dots, \underline{\xi}(\underline{x}_{N^{\text{obs}}}), \underline{\eta}) . \quad (6) $$

We assume that $\underline{y}^{\text{obs}}$ depends on a scalar valued material parameter field $p$, defined on some Hilbert space $V$, through some non-linear operator $\mathcal{A}$ defined on $V$ as

$$ \underline{y}^{\text{obs}} = \mathcal{A}(p) + \underline{\epsilon}. \quad (7) $$

Before going into the Bayesian formulation details, we introduce a digital twin of the experiment, which is modelled as the static equilibrium boundary value problem over $\Omega^{\text{obs}}$ Find $\underline{u} \in H^1(\Omega^{\text{obs}})$ such that

$$ \nabla \cdot \underline{\underline{\sigma}} = \underline{0} \quad \text{over } \Omega^{\text{obs}}, \quad (8) $$

$$ \underline{\underline{\sigma}} = \bar{C}(p) : \underline{\underline{\epsilon}}, \quad (9) $$

$$ \underline{\underline{\epsilon}} = \frac{\nabla \underline{u} + (\nabla \underline{u})^T}{2}, \quad (10) $$

$$ \underline{u} = \underline{u}_D^{\text{obs}} \text{ over } \partial \Omega_D^{\text{obs}}, \quad (11) $$

$$ \underline{\underline{\sigma}} \cdot \underline{n} = \underline{0} \quad \text{over } \partial \Omega_N^{\text{obs}}, \quad (12) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: The two considered specimen: a thin dog bone geometry for uniaxial tensile test, and a thick biaxial geometry.</figcaption></figure>

represents the specimen surface captured by the camera, as for DVC, $\Omega^{\text{obs}}$ represents a volume subregion imaged through, e.g., x-ray tomography. We represent in figure 1 and 11 two specimen geometries where the different regions are annotated, illustrating both surface and volume observation zones. We model the results of image correlation by a set of pointwise, noisy displacement observations

$$ \underline{u}^{\text{obs}} = (\underline{u}^{\text{obs}}(\underline{x}_1), \dots, \underline{u}^{\text{obs}}(\underline{x}_{N^{\text{obs}}})) , \quad (2) $$

$$ \underline{u}^{\text{obs}}(\underline{x}_i) = \underline{u}(\underline{x}_i) + \underline{\xi}, \quad 1 \le i \le N^{\text{obs}}, \quad (3) $$

where $N^{\text{obs}}$ is the number of pointwise observations, and $\underline{u}^{\text{obs}}$ is a vector of size $3 \times N^{\text{obs}}$. Here, $\underline{u}$ denotes the true displacement field, and $\underline{\xi}$ some observation noise. The observed macroscopic force $F^{\text{obs}} \in \mathbb{R}$ models the measured force $F$, typically by a force sensor, that is also corrupt by some noise $\eta$, such that

$$ F^{\text{obs}} = F + \eta. \quad (4) $$

Without loss of generality, we assume that there is only a single force sensor. In other cases (e.g. biaxial traction), $F$ should be taken as vector valued. Let $\underline{y}^{\text{obs}}$ and $\underline{\epsilon}$ be vectors of size $3 \times N^{\text{obs}} + 1$ collecting all measured values and noise values respectively:

$$ \underline{y}^{\text{obs}} = (\underline{u}^{\text{obs}}(\underline{x}_1), \dots, \underline{u}^{\text{obs}}(\underline{x}_{N^{\text{obs}}}), F^{\text{obs}}), \quad (5) $$

$$ \underline{\epsilon} = (\underline{\xi}(\underline{x}_1), \dots, \underline{\xi}(\underline{x}_{N^{\text{obs}}}), \underline{\eta}) . \quad (6) $$

We assume that $\underline{y}^{\text{obs}}$ depends on a scalar valued material parameter field $p$, defined on some Hilbert space $V$, through some non-linear operator $\mathcal{A}$ defined on $V$ as

$$ \underline{y}^{\text{obs}} = \mathcal{A}(p) + \underline{\epsilon}. \quad (7) $$

Before going into the Bayesian formulation details, we introduce a digital twin of the experiment, which is modelled as the static equilibrium boundary value problem over $\Omega^{\text{obs}}$ Find $\underline{u} \in H^1(\Omega^{\text{obs}})$ such that

$$ \nabla \cdot \underline{\underline{\sigma}} = \underline{0} \quad \text{over } \Omega^{\text{obs}}, \quad (8) $$

$$ \underline{\underline{\sigma}} = \bar{C}(p) : \underline{\underline{\epsilon}}, \quad (9) $$

$$ \underline{\underline{\epsilon}} = \frac{\nabla \underline{u} + (\nabla \underline{u})^T}{2}, \quad (10) $$

$$ \underline{u} = \underline{u}_D^{\text{obs}} \text{ over } \partial \Omega_D^{\text{obs}}, \quad (11) $$

$$ \underline{\underline{\sigma}} \cdot \underline{n} = \underline{0} \quad \text{over } \partial \Omega_N^{\text{obs}}, \quad (12) $$
where $\underline{u}$ is the displacement field, $\underline{\varepsilon}$ and $\underline{\sigma}$ are respectively the strain and stress tensor fields, $\overline{C}(p) \in L^\infty(\Omega^{\text{obs}})$ is the stiffness tensor, and $\underline{u}_D^{\text{obs}}$ is a prescribed displacement boundary condition over $\partial\Omega_D^{\text{obs}} \subset \partial\Omega^{\text{obs}}$, and $\underline{n}$ is a outward boundary normal unit vector. Here, we emphasize on the fact that $\overline{C}(p)$ is not assumed to be constant over $\Omega^{\text{obs}}$. To keep notations concise, and without loss of generality, the stiffness model is assumed to depend on a single parameter $p$, and is assumed to be differentiable. The stiffness model is typically material specific, and to remain as general as possible, no further assumption is made in this section. The considered constitutive relation (9) is linear as we're only interested in identification of elastic properties, and extensions to other cases are out of scope of this paper. The measured force is then modelled as the average stress over the traction direction, which can be expressed as a volume integral as

$$F = \frac{|\partial \Omega_D^{\text{obs}}|}{|\Omega^{\text{obs}}|} \int_{\Omega^{\text{obs}}} \underline{n}_F \cdot \underline{\underline{\sigma}}(\underline{u}(p)) \cdot \underline{n}_F \, d\underline{x}, \quad (13)$$

where

$$|\partial \Omega_D^{\text{obs}}| = \int_{\partial \Omega_D^{\text{obs}}} d\underline{S}, \quad (14)$$

$$|\Omega^{\text{obs}}| = \int_{\Omega^{\text{obs}}} d\underline{x}, \quad (15)$$

and $\underline{n}_F$ is a unit vector aligned with the loading direction. As usual, we introduce the following Sobolev function spaces

$$H_D^1(\Omega^{\text{obs}}) = \{ \underline{v} \in H^1(\Omega^{\text{obs}}) \mid \underline{v}|_{\partial \Omega_D^{\text{obs}}} = \underline{u}_D \}, \quad (16)$$

$$H_0^1(\Omega^{\text{obs}}) = \{ \underline{v} \in H^1(\Omega^{\text{obs}}) \mid \underline{v}|_{\partial \Omega_D^{\text{obs}}} = \underline{0} \}, \quad (17)$$

and reformulate equations 8-11 in their weak form: Find $\underline{u} \in H_D^1(\Omega^{\text{obs}})$ and $\underline{v} \in H_0^1(\Omega^{\text{obs}})$ such that

$$a(\underline{u}, \underline{v}) = l(\underline{v}), \quad (18)$$

$$a(\underline{u}, \underline{v}) = \int_{\Omega} \underline{\varepsilon}(\underline{u}) : \overline{C}(p) : \underline{\underline{\varepsilon}}(\underline{v}) \, d\underline{x}, \quad (19)$$

$$l(\underline{v}) = 0. \quad (20)$$

The computational domain $\Omega$ is discretized using a standard finite element mesh, which we enforce to be conforming with the boundary of the observation zone $\partial\Omega_{\text{obs}}$. We consider a mesh of $\Omega^{\text{obs}}$ made of $N_T$ cells ($T_k)_{0 \le k < N_T$}, and introduce the first-order polynomial discretization space,

$$V(\Omega^{\text{obs}}) = \{ u \in \mathcal{C}^0(\Omega^{\text{obs}}) \mid \forall 0 \le k < N_T, u|_{T_k} \in P^1 \}, \quad (21)$$

where $P^1$ denotes the space of first order polynomials. To keep notations simple, the $P^1$ discretization space for vector valued functions is noted $\underline{V}(\Omega^{\text{obs}})$. We introduce finite elements bases for $V(\Omega^{\text{obs}})$ and $\underline{V}(\Omega^{\text{obs}})$:

$$V(\Omega^{\text{obs}}) = \operatorname{span}(\varphi_0, \dots, \varphi_{N_V}), \quad (22)$$

$$\underline{V}(\Omega^{\text{obs}}) = \operatorname{span} \left( \underline{\varphi}_0, \dots, \underline{\varphi}_{N_V} \right). \quad (23)$$

With these notations, the simulated displacement and unknown parameter are discretized in $\underline{V}(\Omega^{\text{obs}})$ and $V(\Omega^{\text{obs}})$ respectively:

$$\underline{u} = \sum_{k=0}^{N_V} u_k \underline{\varphi}_k, \quad p = \sum_{k=0}^{N_V} p_k \varphi_k. \quad (24)$$

We note $\underline{u} = (u_0, \dots, u_{N_V})$ and $\underline{p} = (p_0, \dots, p_{N_V})$. For simplicity, we assume that observation points coincide with mesh nodes, i.e. degrees of freedom coordinates of space $\underline{V}(\Omega^{\text{obs}})$ are $\underline{x}_1, \dots, \underline{x}_{N_{\text{obs}}}$, such that the notation $\underline{u}^{\text{obs}} - \underline{u}$ makes sense. With regard to the previously introduced operator, we pose $\underline{y} = (u_0, \dots, u_{N_V}, F) = \mathcal{A}(p)$.
where $\underline{u}$ is the displacement field, $\underline{\varepsilon}$ and $\underline{\sigma}$ are respectively the strain and stress tensor fields, $\overline{C}(p) \in L^\infty(\Omega^{\text{obs}})$ is the stiffness tensor, and $\underline{u}_D^{\text{obs}}$ is a prescribed displacement boundary condition over $\partial\Omega_D^{\text{obs}} \subset \partial\Omega^{\text{obs}}$, and $\underline{n}$ is a outward boundary normal unit vector. Here, we emphasize on the fact that $\overline{C}(p)$ is not assumed to be constant over $\Omega^{\text{obs}}$. To keep notations concise, and without loss of generality, the stiffness model is assumed to depend on a single parameter $p$, and is assumed to be differentiable. The stiffness model is typically material specific, and to remain as general as possible, no further assumption is made in this section. The considered constitutive relation (9) is linear as we're only interested in identification of elastic properties, and extensions to other cases are out of scope of this paper. The measured force is then modelled as the average stress over the traction direction, which can be expressed as a volume integral as

$$F = \frac{|\partial \Omega_D^{\text{obs}}|}{|\Omega^{\text{obs}}|} \int_{\Omega^{\text{obs}}} \underline{n}_F \cdot \underline{\underline{\sigma}}(\underline{u}(p)) \cdot \underline{n}_F \, d\underline{x}, \quad (13)$$

where

$$|\partial \Omega_D^{\text{obs}}| = \int_{\partial \Omega_D^{\text{obs}}} d\underline{S}, \quad (14)$$

$$|\Omega^{\text{obs}}| = \int_{\Omega^{\text{obs}}} d\underline{x}, \quad (15)$$

and $\underline{n}_F$ is a unit vector aligned with the loading direction. As usual, we introduce the following Sobolev function spaces

$$H_D^1(\Omega^{\text{obs}}) = \{ \underline{v} \in H^1(\Omega^{\text{obs}}) \mid \underline{v}|_{\partial \Omega_D^{\text{obs}}} = \underline{u}_D \}, \quad (16)$$

$$H_0^1(\Omega^{\text{obs}}) = \{ \underline{v} \in H^1(\Omega^{\text{obs}}) \mid \underline{v}|_{\partial \Omega_D^{\text{obs}}} = \underline{0} \}, \quad (17)$$

and reformulate equations 8-11 in their weak form: Find $\underline{u} \in H_D^1(\Omega^{\text{obs}})$ and $\underline{v} \in H_0^1(\Omega^{\text{obs}})$ such that

$$a(\underline{u}, \underline{v}) = l(\underline{v}), \quad (18)$$

$$a(\underline{u}, \underline{v}) = \int_{\Omega} \underline{\varepsilon}(\underline{u}) : \overline{C}(p) : \underline{\varepsilon}(\underline{v}) \, d\underline{x}, \quad (19)$$

$$l(\underline{v}) = 0. \quad (20)$$

The computational domain $\Omega$ is discretized using a standard finite element mesh, which we enforce to be conforming with the boundary of the observation zone $\partial\Omega_{\text{obs}}$. We consider a mesh of $\Omega^{\text{obs}}$ made of $N_T$ cells ($T_k)_{0 \le k < N_T$}, and introduce the first-order polynomial discretization space,

$$V(\Omega^{\text{obs}}) = \{ u \in \mathcal{C}^0(\Omega^{\text{obs}}) \mid \forall 0 \le k < N_T, u|_{T_k} \in P^1 \}, \quad (21)$$

where $P^1$ denotes the space of first order polynomials. To keep notations simple, the $P^1$ discretization space for vector valued functions is noted $\underline{V}(\Omega^{\text{obs}})$. We introduce finite elements bases for $V(\Omega^{\text{obs}})$ and $\underline{V}(\Omega^{\text{obs}})$:

$$V(\Omega^{\text{obs}}) = \operatorname{span}(\varphi_0, \dots, \varphi_{N_V}), \quad (22)$$

$$\underline{V}(\Omega^{\text{obs}}) = \operatorname{span}(\underline{\varphi}_0, \dots, \underline{\varphi}_{N_V}). \quad (23)$$

With these notations, the simulated displacement and unknown parameter are discretized in $\underline{V}(\Omega^{\text{obs}})$ and $V(\Omega^{\text{obs}})$ respectively:

$$\underline{u} = \sum_{k=0}^{N_V} u_k \underline{\varphi}_k, \quad p = \sum_{k=0}^{N_V} p_k \varphi_k. \quad (24)$$

We note $\underline{u} = (u_0, \dots, u_{N_V})$ and $\underline{p} = (p_0, \dots, p_{N_V})$. For simplicity, we assume that observation points coincide with mesh nodes, i.e. degrees of freedom coordinates of space $\underline{V}(\Omega^{\text{obs}})$ are $\underline{x}_1, \dots, \underline{x}_{N_{\text{obs}}}$, such that the notation $\underline{u}^{\text{obs}} - \underline{u}$ makes sense. With regard to the previously introduced operator, we pose $\underline{y} = (u_0, \dots, u_{N_V}, F) = \mathcal{A}(p)$.
## 5.2 Bayesian Formulation

In a Bayesian inference setting, the problem of parameter identification is to estimate $\boldsymbol{p}$ given observations $\boldsymbol{y}^{\text{obs}}$. Following the Bayesian approach, we model both $\boldsymbol{p}$ and $\boldsymbol{\epsilon}$ as random variables. First, we assume that the observation noise $\boldsymbol{\xi}$ and $\boldsymbol{\eta}$ are Gaussian:

$$ \boldsymbol{\xi} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{\boldsymbol{u}}), \qquad (25) $$

$$ \boldsymbol{\eta} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\sigma}_{\boldsymbol{F}}^{2}), \qquad (26) $$

and we note

$$ \boldsymbol{\epsilon} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{\boldsymbol{y}}), \qquad (27) $$

where $\boldsymbol{\Sigma}_{\boldsymbol{u}}$ is the displacement noise covariance matrix, $\boldsymbol{\sigma}_{\boldsymbol{F}}^2$ is the measured force uncertainty variance, and $\boldsymbol{\Sigma}_{\boldsymbol{y}}$ is the global noise covariance matrix, defined blockwise as

$$ \boldsymbol{\Sigma}_{\boldsymbol{y}} = \begin{pmatrix} \boldsymbol{\Sigma}_{\boldsymbol{u}} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\sigma}_{\boldsymbol{F}}^2 \end{pmatrix}. \qquad (28) $$

We assume $\boldsymbol{\Sigma}_{\boldsymbol{y}}$ to be fully known, and its exact expression depends on the experimental setup and to the algorithm used to perform the correlation. For instance, when performing global DIC, the displacement noise is not a white noise, but instead is spatially correlated, and closed form expressions can be obtained in some cases. We give details and point to the relevant literature in appendix B. Following the Bayesian approach, we are interested in the posterior distribution of parameter $\boldsymbol{p}$, given observations $\boldsymbol{y}^{\text{obs}}$, which, using the Bayes theorem, can be expressed as follows:

$$ D(\boldsymbol{p}|\boldsymbol{y}^{\text{obs}}) \propto D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p})D(\boldsymbol{p}), \qquad (29) $$

where $D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p})$ is the likelihood, and $D(\boldsymbol{p})$ is the prior density, which remains to be specified. The posterior density can be sampled from using a Markov chain Monte Carlo (MCMC) algorithm, which is useful to perform uncertainty quantification. However, in the present paper we are only interested in point estimates of the unknown parameter, which are usually obtained by the means of a Maximum A Posteriori (MAP) estimator, which maximises the posterior density, and can be computed by the following optimization problem

$$ \boldsymbol{p}_{\text{MAP}} = \underset{\boldsymbol{p}}{\operatorname{argmin}} - \log D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p}) - \log D(\boldsymbol{p}). \qquad (30) $$

Using assumptions (7) and (27), the likelihood can be expressed as

$$ D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p}) \propto \exp \left[ -\frac{1}{2} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y})^{\mathrm{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y}) \right]. \qquad (31) $$

As explained in the introduction, in order to model spatial variability of elastic properties, we model $\boldsymbol{p}$ as a stationary Gaussian random field, which translates to a multivariate Gaussian prior on $\boldsymbol{p}$:

$$ \boldsymbol{p} \sim \mathcal{N}(\boldsymbol{p}^{0}, \boldsymbol{\Sigma}_{\boldsymbol{p}}(\boldsymbol{\mu})), \qquad (32) $$

where the mean $\boldsymbol{p}^0$ is assumed to be known, $\boldsymbol{\Sigma}_{\boldsymbol{p}}$ is the prior covariance matrix, and $\boldsymbol{\mu}$ is a set of covariance parameters, referred to as hyperparameters for the rest of this paper. In this section, $\boldsymbol{\mu}$ is assumed to be known. Formula (30) now writes

$$ \boldsymbol{p}_{\text{MAP}} = \underset{\boldsymbol{p}}{\operatorname{argmin}} \frac{1}{2} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y})^{\mathrm{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y}) + \frac{1}{2} (\boldsymbol{p} - \boldsymbol{p}^{0})^{\mathrm{T}} \boldsymbol{\Sigma}_{\boldsymbol{p}}^{-1} (\boldsymbol{p} - \boldsymbol{p}^{0}). \qquad (33) $$

This formulation for the MAP estimate leads to a scalar valued function $J$ to minimize. For the subsequent explanations, we split $J$ into three terms:

$$ J(\boldsymbol{p}) = J_{\boldsymbol{u}}(\boldsymbol{p}) + J_{\boldsymbol{F}}(\boldsymbol{p}) + J_{\boldsymbol{p}}(\boldsymbol{p}), \qquad (34) $$

$$ J_{\boldsymbol{u}}(\boldsymbol{p}) = \frac{1}{2}(\boldsymbol{u}^{\text{obs}} - \boldsymbol{u})^{\mathrm{T}}\boldsymbol{\Sigma}_{\boldsymbol{u}}^{-1}(\boldsymbol{u}^{\text{obs}} - \boldsymbol{u}), \qquad (35) $$

$$ J_{\boldsymbol{F}}(\boldsymbol{p}) = \frac{(\boldsymbol{F} - \boldsymbol{F}^{\text{obs}})^2}{2\boldsymbol{\sigma}_{\boldsymbol{F}}^2}, \qquad (36) $$

$$ J_{\boldsymbol{p}}(\boldsymbol{p}) = \frac{1}{2}(\boldsymbol{p} - \boldsymbol{p}^{0})^{\mathrm{T}}\boldsymbol{\Sigma}_{\boldsymbol{p}}^{-1}(\boldsymbol{p} - \boldsymbol{p}^{0}). \qquad (37) $$

Using a gradient based optimization algorithm to minimize $J$ is not straightforward for two reasons:
## 5.2 Bayesian Formulation

In a Bayesian inference setting, the problem of parameter identification is to estimate $\boldsymbol{p}$ given observations $\boldsymbol{y}^{\text{obs}}$. Following the Bayesian approach, we model both $\boldsymbol{p}$ and $\boldsymbol{\epsilon}$ as random variables. First, we assume that the observation noise $\boldsymbol{\xi}$ and $\boldsymbol{\eta}$ are Gaussian:

$$ \boldsymbol{\xi} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{\boldsymbol{u}}), \qquad (25) $$

$$ \boldsymbol{\eta} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\sigma}_{\boldsymbol{F}}^{2}), \qquad (26) $$

and we note

$$ \boldsymbol{\epsilon} \sim \mathcal{N}(\boldsymbol{0}, \boldsymbol{\Sigma}_{\boldsymbol{y}}), \qquad (27) $$

where $\boldsymbol{\Sigma}_{\boldsymbol{u}}$ is the displacement noise covariance matrix, $\boldsymbol{\sigma}_{\boldsymbol{F}}^2$ is the measured force uncertainty variance, and $\boldsymbol{\Sigma}_{\boldsymbol{y}}$ is the global noise covariance matrix, defined blockwise as

$$ \boldsymbol{\Sigma}_{\boldsymbol{y}} = \begin{pmatrix} \boldsymbol{\Sigma}_{\boldsymbol{u}} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{\sigma}_{\boldsymbol{F}}^2 \end{pmatrix}. \qquad (28) $$

We assume $\boldsymbol{\Sigma}_{\boldsymbol{y}}$ to be fully known, and its exact expression depends on the experimental setup and to the algorithm used to perform the correlation. For instance, when performing global DIC, the displacement noise is not a white noise, but instead is spatially correlated, and closed form expressions can be obtained in some cases. We give details and point to the relevant literature in appendix B. Following the Bayesian approach, we are interested in the posterior distribution of parameter $\boldsymbol{p}$, given observations $\boldsymbol{y}^{\text{obs}}$, which, using the Bayes theorem, can be expressed as follows:

$$ D(\boldsymbol{p}|\boldsymbol{y}^{\text{obs}}) \propto D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p})D(\boldsymbol{p}), \qquad (29) $$

where $D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p})$ is the likelihood, and $D(\boldsymbol{p})$ is the prior density, which remains to be specified. The posterior density can be sampled from using a Markov chain Monte Carlo (MCMC) algorithm, which is useful to perform uncertainty quantification. However, in the present paper we are only interested in point estimates of the unknown parameter, which are usually obtained by the means of a Maximum A Posteriori (MAP) estimator, which maximises the posterior density, and can be computed by the following optimization problem

$$ \boldsymbol{p}_{\text{MAP}} = \underset{\boldsymbol{p}}{\operatorname{argmin}} - \log D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p}) - \log D(\boldsymbol{p}). \qquad (30) $$

Using assumptions (7) and (27), the likelihood can be expressed as

$$ D(\boldsymbol{y}^{\text{obs}}|\boldsymbol{p}) \propto \exp \left[ -\frac{1}{2} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y})^{\text{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y}) \right]. \qquad (31) $$

As explained in the introduction, in order to model spatial variability of elastic properties, we model $\boldsymbol{p}$ as a stationary Gaussian random field, which translates to a multivariate Gaussian prior on $\boldsymbol{p}$:

$$ \boldsymbol{p} \sim \mathcal{N}(\boldsymbol{p}^{0}, \boldsymbol{\Sigma}_{\boldsymbol{p}}(\boldsymbol{\mu})), \qquad (32) $$

where the mean $\boldsymbol{p}^0$ is assumed to be known, $\boldsymbol{\Sigma}_{\boldsymbol{p}}$ is the prior covariance matrix, and $\boldsymbol{\mu}$ is a set of covariance parameters, referred to as hyperparameters for the rest of this paper. In this section, $\boldsymbol{\mu}$ is assumed to be known. Formula (30) now writes

$$ \boldsymbol{p}_{\text{MAP}} = \underset{\boldsymbol{p}}{\operatorname{argmin}} \frac{1}{2} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y})^{\text{T}} \boldsymbol{\Sigma}_{\boldsymbol{y}}^{-1} (\boldsymbol{y}^{\text{obs}} - \boldsymbol{y}) + \frac{1}{2} (\boldsymbol{p} - \boldsymbol{p}^{0})^{\text{T}} \boldsymbol{\Sigma}_{\boldsymbol{p}}^{-1} (\boldsymbol{p} - \boldsymbol{p}^{0}). \qquad (33) $$

This formulation for the MAP estimate leads to a scalar valued function $J$ to minimize. For the subsequent explanations, we split $J$ into three terms:

$$ J(\boldsymbol{p}) = J_{\boldsymbol{u}}(\boldsymbol{p}) + J_{\boldsymbol{F}}(\boldsymbol{p}) + J_{\boldsymbol{p}}(\boldsymbol{p}), \qquad (34) $$

$$ J_{\boldsymbol{u}}(\boldsymbol{p}) = \frac{1}{2}(\boldsymbol{u}^{\text{obs}} - \boldsymbol{u})^{\text{T}}\boldsymbol{\Sigma}_{\boldsymbol{u}}^{-1}(\boldsymbol{u}^{\text{obs}} - \boldsymbol{u}), \qquad (35) $$

$$ J_{\boldsymbol{F}}(\boldsymbol{p}) = \frac{(\boldsymbol{F} - \boldsymbol{F}^{\text{obs}})^2}{2\boldsymbol{\sigma}_{\boldsymbol{F}}^2}, \qquad (36) $$

$$ J_{\boldsymbol{p}}(\boldsymbol{p}) = \frac{1}{2}(\boldsymbol{p} - \boldsymbol{p}^{0})^{\text{T}}\boldsymbol{\Sigma}_{\boldsymbol{p}}^{-1}(\boldsymbol{p} - \boldsymbol{p}^{0}). \qquad (37) $$

Using a gradient based optimization algorithm to minimize $J$ is not straightforward for two reasons:
1. The prior precision matrix $\Sigma_p^{-1}$ is dense and thus intractable on larger meshes,

2. Computing $\nabla J$ involves computing $\frac{du}{dp}$ and a finite difference approximation (which is the common approach in the FEMU literature) is intractable on larger meshes.

We address both points in sections 6.1 and 8.1 respectively.

## 5.3 Parametric Empirical Bayes Formulation

In this section, we no longer assume the hyperparameters $\mu = (\mu_k)_{0 \le k < N_\mu}$ to be known and propose a Bayesian approach for hyperparameter estimation. More precisely, we model $\mu$ as a random vector, and formulate a log-normal hyperprior on each component of $\mu$:

$$ \forall 0 \le i < N_\mu, \log \mu_i \sim \mathcal{N}(\tau_i^0, \sigma_{\tau_i}^2). \quad (38) $$

The choice of the log-normal hyperprior is due to the fact that covariance hyperparameters are often positive. In the rest of this section, we pose $\tau = (\log \mu_k)_{i \le k < N_\mu}$. The posterior distribution (29) is then updated as follows:

$$ D(p, \tau|y^{\text{obs}}) \propto D(y^{\text{obs}}|p)D(p|\tau)D(\tau). \quad (39) $$

The expression for the prior $D(p|\tau)$ in (37) must be updated to take into account the normalization term that depends on $\mu$:

$$ J_p(p) = \frac{1}{2}(p - p^0)^T \Sigma_p^{-1} (p - p^0) - \frac{1}{2} \log \det \Sigma_p^{-1} + \frac{N_V}{2} \log(2\pi). \quad (40) $$

In the resulting hierarchical Bayesian formulation, the posterior distribution can be approximated with MCMC sampling (e.g. Gibbs sampling, see section 11.3 of [Bishop, 2006]). As we are only interested in the MAP estimate of the unknown parameter, like in section 5.2, it is natural to approximate the hierarchical model with a *Parametric Empirical Bayes* approach (PEB), where the hyperparameters are iteratively updated from the parameter MAP estimator, following the Bayes formula at the hypermodel level:

$$ D(\tau|p_{\text{MAP}}) \propto D(p_{\text{MAP}}|\tau)D(\tau). \quad (41) $$

The formulated hypermodel 38 suggests updating the hyperparameters with a MAP estimator $\tau_{\text{MAP}}$ of $\tau$:

$$ \tau_{\text{MAP}} = \underset{\tau}{\operatorname{argmin}} J_p(p_{\text{MAP}}, \tau) + J_{\tau}(\tau), \quad (42) $$

where $J_\tau(\tau)$ is a new term corresponding to the log-normal hyperprior:

$$ J_{\tau}(\tau) = \sum_{i=0}^{N_{\mu}-1} \frac{(\tau_i - \tau_i^0)^2}{2\sigma_{\tau_i}}. \quad (43) $$

The resulting formulation in quite inconvenient as we now have $2N_\mu$ hyper-hyperparameters to specify, which are typically also unknown. A common way to tackle this issue is to specify an improper hyperprior as opposed to a log-normal distribution like in equation (38):

$$ \forall 0 \le i < N_\mu, \log \mu_i \sim \mathcal{U}(\mathbb{R}), \quad (44) $$

which leads to $J_\tau(\tau) = 0$. From an empirical Bayes point of view, this is equivalent to updating the hyperparameters with a maximum likelihood (ML) estimator as opposed to a MAP estimator:

$$ \tau_{\text{ML}} = \underset{\tau}{\operatorname{argmin}} J_p(p_{\text{MAP}}, \tau). \quad (45) $$

In section 6.2, we show the numerical shortcomings of using the ML estimator for updating the hyperparameters, and in particular we explain why the resulting optimization problem is unstable when the Matérn covariance is used for prior modelling. To address this issue, we propose in section 6.2 an alternative updating scheme, where the hyperparameter are iteratively identified through variogram regression, another variation of the PEB approach.
1. The prior precision matrix $\Sigma_p^{-1}$ is dense and thus intractable on larger meshes,

2. Computing $\nabla J$ involves computing $\frac{du}{dp}$ and a finite difference approximation (which is the common approach in the FEMU literature) is intractable on larger meshes.

We address both points in sections 6.1 and 8.1 respectively.

## 5.3 Parametric Empirical Bayes Formulation

In this section, we no longer assume the hyperparameters $\mu = (\mu_k)_{0 \le k < N_\mu}$ to be known and propose a Bayesian approach for hyperparameter estimation. More precisely, we model $\mu$ as a random vector, and formulate a log-normal hyperprior on each component of $\mu$:

$$ \forall 0 \le i < N_\mu, \log \mu_i \sim \mathcal{N}(\tau_i^0, \sigma_{\tau_i}^2). \qquad (38) $$

The choice of the log-normal hyperprior is due to the fact that covariance hyperparameters are often positive. In the rest of this section, we pose $\tau = (\log \mu_k)_{i \le k < N_\mu}$. The posterior distribution (29) is then updated as follows:

$$ D(p, \tau|y^{\text{obs}}) \propto D(y^{\text{obs}}|p)D(p|\tau)D(\tau). \qquad (39) $$

The expression for the prior $D(p|\tau)$ in (37) must be updated to take into account the normalization term that depends on $\mu$:

$$ J_p(p) = \frac{1}{2}(p - p^0)^T \Sigma_p^{-1} (p - p^0) - \frac{1}{2} \log \det \Sigma_p^{-1} + \frac{N_V}{2} \log(2\pi). \qquad (40) $$

In the resulting hierarchical Bayesian formulation, the posterior distribution can be approximated with MCMC sampling (e.g. Gibbs sampling, see section 11.3 of [Bishop, 2006]). As we are only interested in the MAP estimate of the unknown parameter, like in section 5.2, it is natural to approximate the hierarchical model with a *Parametric Empirical Bayes* approach (PEB), where the hyperparameters are iteratively updated from the parameter MAP estimator, following the Bayes formula at the hypermodel level:

$$ D(\tau|p_{\text{MAP}}) \propto D(p_{\text{MAP}}|\tau)D(\tau). \qquad (41) $$

The formulated hypermodel 38 suggests updating the hyperparameters with a MAP estimator $\tau_{\text{MAP}}$ of $\tau$:

$$ \tau_{\text{MAP}} = \underset{\tau}{\operatorname{argmin}} J_p(p_{\text{MAP}}, \tau) + J_{\tau}(\tau), \qquad (42) $$

where $J_\tau(\tau)$ is a new term corresponding to the log-normal hyperprior:

$$ J_{\tau}(\tau) = \sum_{i=0}^{N_{\mu}-1} \frac{(\tau_i - \tau_i^0)^2}{2\sigma_{\tau_i}}. \qquad (43) $$

The resulting formulation in quite inconvenient as we now have $2N_\mu$ hyper-hyperparameters to specify, which are typically also unknown. A common way to tackle this issue is to specify an improper hyperprior as opposed to a log-normal distribution like in equation (38):

$$ \forall 0 \le i < N_\mu, \log \mu_i \sim \mathcal{U}(\mathbb{R}), \qquad (44) $$

which leads to $J_\tau(\tau) = 0$. From an empirical Bayes point of view, this is equivalent to updating the hyperparameters with a maximum likelihood (ML) estimator as opposed to a MAP estimator:

$$ \tau_{\text{ML}} = \underset{\tau}{\operatorname{argmin}} J_p(p_{\text{MAP}}, \tau). \qquad (45) $$

In section 6.2, we show the numerical shortcomings of using the ML estimator for updating the hyperparameters, and in particular we explain why the resulting optimization problem is unstable when the Matérn covariance is used for prior modelling. To address this issue, we propose in section 6.2 an alternative updating scheme, where the hyperparameter are iteratively identified through variogram regression, another variation of the PEB approach.
# 6 Spatial Covariance Modelling Through Whittle-Matérn Priors

In the previous section, a multivariate normal prior was formulated on $\mathbf{p}$, however a covariance model remains to be specified. In this section, we present and discuss the use of the Matérn covariance model.

## 6.1 Covariance Precision Discretization

Our choice is motivated by the requirement to obtain an efficient discretization on unstructured grids. In this regard, the Matérn covariance is a suitable choice, as it can lead to a sparse approximation of the precision matrix [Lindgren et al., 2011]. More precisely, the Matérn covariance is defined as

$$M(r; \sigma^2, \ell, \nu) = \frac{\sigma^2 2^{1-\nu}}{\Gamma(\nu)} \left(\frac{r}{\ell}\right)^{\nu} K_{\nu}\left(\frac{r}{\ell}\right), \quad (46)$$

where $\sigma^2$ is the marginal variance, $\ell$ is the length scale and $\nu$ is the smoothness parameter. With regard to the notations introduced in section 5.2, $\mu = (\sigma^2, \ell, \nu)$ is the hyperparameter vector. We also introduce the correlation length $\rho = \ell\sqrt{8\nu}$, which is the distance at which the correlations are near 0.1 [Lindgren et al., 2011]. It was shown in [Lindgren et al., 2011, Roininen et al., 2014] that Gaussian random fields with Matérn covariance are the solutions of the following stochastic PDE:

$$(\mathrm{Id} - \ell^2 \Delta^2)^{(\nu+d/2)/2} x = \sigma \sqrt{\frac{\ell^d}{\beta}} W, \quad (47)$$

where

$$\beta = \frac{\Gamma(\nu)}{2^d \pi^{d/2} \Gamma(\nu + d/2)}, \quad (48)$$

and $d = 2, 3$ is the spatial dimension. We also pose $\alpha = \nu + d/2$. For a broad overview of the different applications related to this link between the Matérn covariance and the SPDE (47), we refer to [Lindgren et al., 2022]. When considering an unstructured mesh discretization, the precision matrix of $\mathbf{x}$ can have a sparse approximation. We refer to [Lindgren et al., 2011] for a proof. More precisely, let $\mathbf{x} = (x_i)_{0 \le i < N_v}$ be the degrees of freedom of a solution of equation (47) discretized in $V(\Omega^{\mathrm{obs}})$, and let $\mathbf{Q}_{\mathbf{x}}(\alpha, \ell)$ be the precision matrix of $\mathbf{x}$. We introduce the usual FEM mass and stiffness matrices for a scalar valued $P^1$ discretization:

$$M = \left( \int_{\Omega^{\mathrm{obs}}} \varphi_i \varphi_j \, d\mathbf{x} \right)_{0 \le i,j < N_v}, \quad (49)$$

$$K = \left( \int_{\Omega^{\mathrm{obs}}} \nabla \varphi_i \cdot \nabla \varphi_j \, d\mathbf{x} \right)_{0 \le i,j < N_v}. \quad (50)$$

We have the following sparse approximations:

$$\mathbf{Q}_{\mathbf{x}}(1, \ell) = \frac{\beta}{\sigma^2 \ell^d} (\mathbf{M} + \ell^2 \mathbf{K}), \quad (51)$$

$$\mathbf{Q}_{\mathbf{x}}(2, \ell) = \frac{\beta}{\sigma^2 \ell^d} (\mathbf{M} + \ell^2 \mathbf{K}) \mathbf{M}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}), \quad (52)$$

$$\forall n \in \mathbb{N}, n \ge 3, \quad \mathbf{Q}_{\mathbf{x}}(n, \ell) = (\mathbf{M} + \ell^2 \mathbf{K}) \mathbf{M}^{-1} \mathbf{Q}_{\mathbf{x}}(n - 2, \ell) \mathbf{M}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}). \quad (53)$$

Formulas (51)-(53) correspond to $\nu = 0, 1, 2, ...$ for $d = 2$ and $\nu = 0.5, 1.5, 2.5, ...$ for $d = 3$. For other values $\alpha > 0$ that are not integers, more recent papers give another approximation of $\mathbf{Q}_{\mathbf{x}}(\alpha, \ell)$ based on, e.g., a truncated Taylor expansion of the underlying pseudodifferential operator or a quadrature of its Brochner integral, and we refer to, e.g. [Bolin et al., 2018, Roininen et al., 2018] for details. Note that formulas (51)-(52) are slightly different to those in [Lindgren et al., 2011], as the base SPDE (47) was taken from [Roininen et al., 2014] which defines different SPDE coefficients. In equation (47), the marginal variance and length scale of $\mathbf{x}$ are controlled by separate hyperparameters, as opposed
# 6 Spatial Covariance Modelling Through Whittle-Matérn Priors

In the previous section, a multivariate normal prior was formulated on $\mathbf{p}$, however a covariance model remains to be specified. In this section, we present and discuss the use of the Matérn covariance model.

## 6.1 Covariance Precision Discretization

Our choice is motivated by the requirement to obtain an efficient discretization on unstructured grids. In this regard, the Matérn covariance is a suitable choice, as it can lead to a sparse approximation of the precision matrix [Lindgren et al., 2011]. More precisely, the Matérn covariance is defined as

$$M(r; \sigma^2, \ell, \nu) = \frac{\sigma^2 2^{1-\nu}}{\Gamma(\nu)} \left(\frac{r}{\ell}\right)^{\nu} K_{\nu}\left(\frac{r}{\ell}\right), \quad (46)$$

where $\sigma^2$ is the marginal variance, $\ell$ is the length scale and $\nu$ is the smoothness parameter. With regard to the notations introduced in section 5.2, $\mu = (\sigma^2, \ell, \nu)$ is the hyperparameter vector. We also introduce the correlation length $\rho = \ell\sqrt{8\nu}$, which is the distance at which the correlations are near 0.1 [Lindgren et al., 2011]. It was shown in [Lindgren et al., 2011, Roininen et al., 2014] that Gaussian random fields with Matérn covariance are the solutions of the following stochastic PDE:

$$(\mathrm{Id} - \ell^2 \Delta^2)^{(\nu+d/2)/2} x = \sigma \sqrt{\frac{\ell^d}{\beta}} W, \quad (47)$$

where

$$\beta = \frac{\Gamma(\nu)}{2^d \pi^{d/2} \Gamma(\nu + d/2)}, \quad (48)$$

and $d = 2, 3$ is the spatial dimension. We also pose $\alpha = \nu + d/2$. For a broad overview of the different applications related to this link between the Matérn covariance and the SPDE (47), we refer to [Lindgren et al., 2022]. When considering an unstructured mesh discretization, the precision matrix of $\mathbf{x}$ can have a sparse approximation. We refer to [Lindgren et al., 2011] for a proof. More precisely, let $\mathbf{x} = (x_i)_{0 \le i < N_V}$ be the degrees of freedom of a solution of equation (47) discretized in $V(\Omega^{\text{obs}})$, and let $\mathbf{Q}_{\mathbf{x}}(\alpha, \ell)$ be the precision matrix of $\mathbf{x}$. We introduce the usual FEM mass and stiffness matrices for a scalar valued $P^1$ discretization:

$$M = \left( \int_{\Omega^{\text{obs}}} \varphi_i \varphi_j \, d\mathbf{x} \right)_{0 \le i,j < N_V}, \quad (49)$$

$$K = \left( \int_{\Omega^{\text{obs}}} \nabla \varphi_i \cdot \nabla \varphi_j \, d\mathbf{x} \right)_{0 \le i,j < N_V}. \quad (50)$$

We have the following sparse approximations:

$$\mathbf{Q}_{\mathbf{x}}(1, \ell) = \frac{\beta}{\sigma^2 \ell^d} (\mathbf{M} + \ell^2 \mathbf{K}), \quad (51)$$

$$\mathbf{Q}_{\mathbf{x}}(2, \ell) = \frac{\beta}{\sigma^2 \ell^d} (\mathbf{M} + \ell^2 \mathbf{K}) \mathbf{M}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}), \quad (52)$$

$$\forall n \in \mathbb{N}, n \ge 3, \quad \mathbf{Q}_{\mathbf{x}}(n, \ell) = (\mathbf{M} + \ell^2 \mathbf{K}) \mathbf{M}^{-1} \mathbf{Q}_{\mathbf{x}}(n - 2, \ell) \mathbf{M}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}). \quad (53)$$

Formulas (51)-(53) correspond to $\nu = 0, 1, 2, \dots$ for $d = 2$ and $\nu = 0.5, 1.5, 2.5, \dots$ for $d = 3$. For other values $\alpha > 0$ that are not integers, more recent papers give another approximation of $\mathbf{Q}_{\mathbf{x}}(\alpha, \ell)$ based on, e.g., a truncated Taylor expansion of the underlying pseudodifferential operator or a quadrature of its Brochner integral, and we refer to, e.g. [Bolin et al., 2018, Roininen et al., 2018] for details. Note that formulas (51)-(52) are slightly different to those in [Lindgren et al., 2011], as the base SPDE (47) was taken from [Roininen et al., 2014] which defines different SPDE coefficients. In equation (47), the marginal variance and length scale of $\mathbf{x}$ are controlled by separate hyperparameters, as opposed
to the range parameter in [Lindgren et al., 2011] which affects both the correlation length and the marginal variance. In (52)-(53), the inverse of the mass matrix can be approximated by the inverse of its lumped approximation, defined as

$$ \forall 0 \le i < N_V, \quad \tilde{M}_{ii} = \sum_{j=0}^{N_V-1} M_{ij}, \qquad (54) $$

$$ \forall i \neq j, \quad \tilde{M}_{ij} = 0. \qquad (55) $$

In regard to the previous Bayesian formulation for the material parameter $p$, fixing $\alpha \ge 2$ to an integer value leads to a sparse expression for the prior (37), that can be implemented efficiently. For instance, with $\alpha = 2$, we can plug (52) into the MAP estimate expression (37), resulting in a so-called Whittle-Matérn prior:

$$ J_p(p) = \frac{\beta}{2\sigma^2\ell^2} (p - p^0)^T (M + \ell^2 K) \tilde{M}^{-1} (M + \ell^2 K) (p - p^0). \qquad (56) $$

## 6.2 Whittle-Matérn Hyperparameters Estimation

As discussed in section 5.3, hyperparameter estimation using the PEB method essentially consists in iteratively estimating $\mu = (\sigma^2, \ell, \nu)$ from parameter estimates $p_{MAP}$. In this section, we discuss several issues related to the use of the Matérn covariance, and in particular, the use of a sparse approximation with a fixed smoothness parameter $\nu$.

First, given the sparse approximation $Q_x$ of $\Sigma_p^{-1}$ in formulas (51)-(53), the mapping $\nu \mapsto Q_x(\alpha(\nu), \ell)$ is not defined for all $\nu > 0$, and thus not differentiable. A gradient based optimization procedure for computing $\mu_{MAP}$ is thus not possible. A possible workaround would be to set $\alpha = 2$ and estimate only $(\sigma^2, \ell)$. Surprisingly, many authors employ this solution, claiming that the data does not contain much information about $\nu$, which is often poorly identified [Diggle and Ribiero, 2007, Bose et al., 2018]. However, this claim is questionable, and for a detailed discussion on this matter, we refer to [Oliveira and Han, 2022].

More precisely, fixing $\nu$ to an arbitrary leads to estimability issues of $\sigma^2$ and $\ell$: there does not exist an estimator for the pair $(\sigma^2, \ell)$ that is weakly consistent under infill asymptotics, i.e. when considering a sequence of increasing number of pointwise observations $((p_k)_{0 \le k < n})_{n \ge 1}$, and a sequence of estimators $(\hat{\mu}_n)_{n \ge 1}$, $\hat{\mu}_k = (\hat{\sigma}_k^2, \hat{\ell}_k)$, then $\hat{\mu}_n$ does not converge when $n \to \infty$. This result is proven in [Zhang, 2004] and thoroughly discussed in [Zhang, 2004, Tang et al., 2021]. Another result of [Zhang, 2004] is that the so-called microergodic parameter $c = \sigma^2/\ell$ is consistently estimable. In the light of these results, in the hierarchical Bayesian model of section 5.3, when setting $\nu$ to an arbitrary value, and considering an improper hyperprior for $(\sigma^2, \ell)$, the minimization problem for the maximum likelihood estimator $\mu_{ML}$ used in the PEB hyperparameters update is unstable, in the sense that $\mu_{ML}$ does not converge when the computational mesh is refined, and that the result of the optimization process is highly sensitive to initialization values due to ridges in the likelihood profile (see, e.g., figure 1 of [Lalchand et al., 2022]). In order to stabilize the hierarchical model, a proper hyperprior (or a fixed value), e.g. a log-normal hyperprior as presented in equation (38), must be set on either $\sigma^2$ or $\ell$, which results in a hybrid approach between ML and MAP estimator for PEB hyperparameter updating, where is referred to as *hybrid-MAP* for the rest of this paper.

Another shortcoming of this approach is that, assuming $\sigma^2$ is known, $\nu$ set to an arbitrary value and $\ell$ estimated using an improper hyperprior, then there is no guarantee that the resulting correlation length $\rho$ matches the statistics of the underlying random process. To illustrate this problem, let's consider a set of one dimensional Matérn field realizations $f_\alpha$ with fixed hyperparameters $\sigma_f^2, \rho_f$ and increasing values for $\alpha$ (i.e. increasing $\nu$ and decreasing $\ell = \rho_f/\sqrt{\nu}$). We wish to identify the correlation length $\rho_f$ using a maximum likelihood estimator, computed by using formula (56) with a smoothness parameter $\alpha_0$ set to an arbitrary value ($\alpha_0 = 2$ in this case). The estimator is then
to the range parameter in [Lindgren et al., 2011] which affects both the correlation length and the marginal variance. In (52)-(53), the inverse of the mass matrix can be approximated by the inverse of its lumped approximation, defined as

$$ \forall 0 \le i < N_V, \quad \tilde{M}_{ii} = \sum_{j=0}^{N_V-1} M_{ij}, \qquad (54) $$

$$ \forall i \neq j, \quad \tilde{M}_{ij} = 0. \qquad (55) $$

In regard to the previous Bayesian formulation for the material parameter $p$, fixing $\alpha \ge 2$ to an integer value leads to a sparse expression for the prior (37), that can be implemented efficiently. For instance, with $\alpha = 2$, we can plug (52) into the MAP estimate expression (37), resulting in a so-called Whittle-Matérn prior:

$$ J_p(p) = \frac{\beta}{2\sigma^2\ell^2} (p - p^0)^T (M + \ell^2 K) \tilde{M}^{-1} (M + \ell^2 K) (p - p^0). \qquad (56) $$

## 6.2 Whittle-Matérn Hyperparameters Estimation

As discussed in section 5.3, hyperparameter estimation using the PEB method essentially consists in iteratively estimating $\mu = (\sigma^2, \ell, \nu)$ from parameter estimates $p_{MAP}$. In this section, we discuss several issues related to the use of the Matérn covariance, and in particular, the use of a sparse approximation with a fixed smoothness parameter $\nu$.

First, given the sparse approximation $Q_x$ of $\Sigma_p^{-1}$ in formulas (51)-(53), the mapping $\nu \mapsto Q_x(\alpha(\nu), \ell)$ is not defined for all $\nu > 0$, and thus not differentiable. A gradient based optimization procedure for computing $\mu_{MAP}$ is thus not possible. A possible workaround would be to set $\alpha = 2$ and estimate only $(\sigma^2, \ell)$. Surprisingly, many authors employ this solution, claiming that the data does not contain much information about $\nu$, which is often poorly identified [Diggle and Ribiero, 2007, Bose et al., 2018]. However, this claim is questionable, and for a detailed discussion on this matter, we refer to [Oliveira and Han, 2022].

More precisely, fixing $\nu$ to an arbitrary leads to estimability issues of $\sigma^2$ and $\ell$: there does not exist an estimator for the pair $(\sigma^2, \ell)$ that is weakly consistent under infill asymptotics, i.e. when considering a sequence of increasing number of pointwise observations $((p_k)_{0 \le k < n})_{n \ge 1}$, and a sequence of estimators $(\hat{\mu}_n)_{n \ge 1}$, $\hat{\mu}_k = (\hat{\sigma}_k^2, \hat{\ell}_k)$, then $\hat{\mu}_n$ does not converge when $n \to \infty$. This result is proven in [Zhang, 2004] and thoroughly discussed in [Zhang, 2004, Tang et al., 2021]. Another result of [Zhang, 2004] is that the so-called microergodic parameter $c = \sigma^2/\ell$ is consistently estimable. In the light of these results, in the hierarchical Bayesian model of section 5.3, when setting $\nu$ to an arbitrary value, and considering an improper hyperprior for $(\sigma^2, \ell)$, the minimization problem for the maximum likelihood estimator $\mu_{ML}$ used in the PEB hyperparameters update is unstable, in the sense that $\mu_{ML}$ does not converge when the computational mesh is refined, and that the result of the optimization process is highly sensitive to initialization values due to ridges in the likelihood profile (see, e.g., figure 1 of [Lalchand et al., 2022]). In order to stabilize the hierarchical model, a proper hyperprior (or a fixed value), e.g. a log-normal hyperprior as presented in equation (38), must be set on either $\sigma^2$ or $\ell$, which results in a hybrid approach between ML and MAP estimator for PEB hyperparameter updating, where is referred to as *hybrid-MAP* for the rest of this paper.

Another shortcoming of this approach is that, assuming $\sigma^2$ is known, $\nu$ set to an arbitrary value and $\ell$ estimated using an improper hyperprior, then there is no guarantee that the resulting correlation length $\rho$ matches the statistics of the underlying random process. To illustrate this problem, let's consider a set of one dimensional Matérn field realizations $f_\alpha$ with fixed hyperparameters $\sigma_f^2, \rho_f$ and increasing values for $\alpha$ (i.e. increasing $\nu$ and decreasing $\ell = \rho_f/\sqrt{\nu}$). We wish to identify the correlation length $\rho_f$ using a maximum likelihood estimator, computed by using formula (56) with a smoothness parameter $\alpha_0$ set to an arbitrary value ($\alpha_0 = 2$ in this case). The estimator is then
computed with

$$ \rho_{\alpha}^{\text{ML}} = \sqrt{1.5} \operatorname{argmin}_{\ell} [-\log L_{\ell}(\mathbf{f})] = \frac{\beta}{2\sigma^2\ell} (\mathbf{f}_{\alpha} - \mathbf{f}_{\alpha}^0)^T (\mathbf{M} + \ell^2 \mathbf{K}) \tilde{\mathbf{M}}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}) (\mathbf{f}_{\alpha} - \mathbf{f}_{\alpha}^0) - \frac{1}{2} \log \det \left[ \frac{\beta}{\sigma^2\ell} (\mathbf{M} + \ell^2 \mathbf{K}) \tilde{\mathbf{M}}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}) \right]. \quad (57) $$

We represent in figure 2 realisations $\mathbf{f}_2, \mathbf{f}_4, \mathbf{f}_{16}$ and the profile of the objective function $-\log L_{\ell}(\mathbf{f}_{\alpha})$

<figure><img src="image_2.png" /><figcaption>Figure 2: Realisations of a one dimensional Matérn field $f$ with several smoothness parameters $\alpha$, and fixed correlation length $\rho_f = \ell\sqrt{\nu}$. On the right, the negative log-likelihood profile of each realisation w.r.t. $\ell$ for a fixed smoothness parameter $\alpha_0 = 2$ for discretization of the Matérn precision.</figcaption></figure>

with the values of the estimated correlation length. The right correlation length is recovered only for $\alpha = 2$, which corresponds to the case where the precision is discretized with the correct smoothness parameter. For all other cases, the correlation length is overestimated by an order of magnitude. This simplified example shows that the hybrid-MAP approach highly overestimates the correlation length when the posterior mean $\mathbf{p}_{\text{MAP}}$ is smoother than the arbitrarily chosen value set for $\nu$.

In regard to this last issue, we also explore an alternative hyperparameter updating scheme that uses variogram regression. This idea closely follows the developments of [Brown et al., 2020]. More precisely, if $\mathbf{z} = (z(\underline{x}_k))_{0 \le k < N_{\text{obs}}}$ are pointwise observations of an underlying Gaussian process $z$ with Matérn covariance $\tilde{M}(\cdot; \sigma^2, \ell, \nu)$, then the empirical variogram is defined as

$$ \hat{\gamma}(r \pm \delta; \mathbf{z}) = \frac{1}{2|\mathcal{N}(r \pm \delta)|} \sum_{(i,j) \in \mathcal{N}(r \pm \delta)} |z_i - z_j|^2, \quad (58) $$

$$ \mathcal{N}(r \pm \delta) = \{0 \le i, j < N_{\text{obs}} : |\underline{x}_i - \underline{x}_j| = r \pm \delta\}, \quad (59) $$

where $\delta$ is the separation distance which is used for binning. The Matérn variogram [Brown et al., 2020] is defined as

$$ \gamma_M(r; a_0, \sigma^2, \ell, \nu) = \begin{cases} 0 & \text{if } r = 0 \\ a_0 + (\sigma^2 - a_0) \left(1 - \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{r}{\ell}\right)^{\nu} K_{\nu}\left(\frac{r}{\ell}\right)\right) & \text{if } r > 0, \end{cases} \quad (60) $$

where the nugget $a_0 \ge 0$ is the variance value at a distance just greater than 0, and the sill $\sigma^2$ is also the marginal variance can be fitted to the empirical variogram $\hat{\gamma}$ through various methods, e.g. weighted least-squares [Cressie, 1985]:

$$ \psi_{\text{opt}} = \operatorname{argmin}_{\psi} J_{\text{vario}}(\psi, \mathbf{z}) := \sum_{r} \frac{|\mathcal{N}(r \pm \delta)|}{\gamma_M(r; \psi)^2} (\hat{\gamma}(r \pm \delta; \mathbf{z}) - \gamma_M(r; \psi))^2, \quad (61) $$

where $\psi = (a_0, \sigma^2, \ell, \nu)$. Although the nugget $a_0$ is estimated by variogram fitting, it is not useful for the considered problem, as it does not appear in the prior expression. A typical fitting method will yield non-integer values for $\nu$. In order to simplify numerical implementation, it is useful to approximate it with the nearest integer values (in 2D) or the nearest half-integer (in 3D). In order to preserve the correlation length, it is then necessary to update the identified length scale accordingly. More precisely, if $\ell_{\text{opt}}, \nu_{\text{opt}}$ are identified hyperparameters, and $\lfloor \nu \rfloor_d$, $d = 2, 3$ denotes the nearest integer or half-integer, then choosing

$$ \tilde{\ell}_{\text{opt}} = \ell_{\text{opt}} \sqrt{\frac{\nu_{\text{opt}}}{\lfloor \nu \rfloor_d}} \quad (62) $$
computed with

$$ \rho_{\alpha}^{\text{ML}} = \sqrt{1.5} \operatorname{argmin}_{\ell} [-\log L_{\ell}(\mathbf{f})] = \frac{\beta}{2\sigma^2\ell} (\mathbf{f}_{\alpha} - \mathbf{f}_{\alpha}^0)^T (\mathbf{M} + \ell^2 \mathbf{K}) \tilde{\mathbf{M}}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}) (\mathbf{f}_{\alpha} - \mathbf{f}_{\alpha}^0) - \frac{1}{2} \log \det \left[ \frac{\beta}{\sigma^2\ell} (\mathbf{M} + \ell^2 \mathbf{K}) \tilde{\mathbf{M}}^{-1} (\mathbf{M} + \ell^2 \mathbf{K}) \right]. \quad (57) $$

We represent in figure 2 realisations $\mathbf{f}_2, \mathbf{f}_4, \mathbf{f}_{16}$ and the profile of the objective function $-\log L_{\ell}(\mathbf{f}_{\alpha})$

<figure><img src="image_3.png" /><figcaption>Figure 2: Realisations of a one dimensional Matérn field $f$ with several smoothness parameters $\alpha$, and fixed correlation length $\rho_f = \ell\sqrt{\nu}$. On the right, the negative log-likelihood profile of each realisation w.r.t. $\ell$ for a fixed smoothness parameter $\alpha_0 = 2$ for discretization of the Matérn precision.</figcaption></figure>

with the values of the estimated correlation length. The right correlation length is recovered only for $\alpha = 2$, which corresponds to the case where the precision is discretized with the correct smoothness parameter. For all other cases, the correlation length is overestimated by an order of magnitude. This simplified example shows that the hybrid-MAP approach highly overestimates the correlation length when the posterior mean $\mathbf{p}_{\text{MAP}}$ is smoother than the arbitrarily chosen value set for $\nu$.

In regard to this last issue, we also explore an alternative hyperparameter updating scheme that uses variogram regression. This idea closely follows the developments of [Brown et al., 2020]. More precisely, if $\mathbf{z} = (z(\underline{x}_k))_{0 \le k < N_{\text{obs}}}$ are pointwise observations of an underlying Gaussian process $z$ with Matérn covariance $\tilde{M}(\cdot; \sigma^2, \ell, \nu)$, then the empirical variogram is defined as

$$ \hat{\gamma}(r \pm \delta; \mathbf{z}) = \frac{1}{2|\mathcal{N}(r \pm \delta)|} \sum_{(i,j) \in \mathcal{N}(r \pm \delta)} |z_i - z_j|^2, \quad (58) $$

$$ \mathcal{N}(r \pm \delta) = \{0 \le i, j < N_{\text{obs}} : |\underline{x}_i - \underline{x}_j| = r \pm \delta\}, \quad (59) $$

where $\delta$ is the separation distance which is used for binning. The Matérn variogram [Brown et al., 2020] is defined as

$$ \gamma_M(r; a_0, \sigma^2, \ell, \nu) = \begin{cases} 0 & \text{if } r = 0 \\ a_0 + (\sigma^2 - a_0) \left(1 - \frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{r}{\ell}\right)^{\nu} K_{\nu}\left(\frac{r}{\ell}\right)\right) & \text{if } r > 0, \end{cases} \quad (60) $$

where the nugget $a_0 \ge 0$ is the variance value at a distance just greater than 0, and the sill $\sigma^2$ is also the marginal variance can be fitted to the empirical variogram $\hat{\gamma}$ through various methods, e.g. weighted least-squares [Cressie, 1985]:

$$ \psi_{\text{opt}} = \operatorname{argmin}_{\psi} J_{\text{vario}}(\psi, \mathbf{z}) := \sum_{r} \frac{|\mathcal{N}(r \pm \delta)|}{\gamma_M(r; \psi)^2} (\hat{\gamma}(r \pm \delta; \mathbf{z}) - \gamma_M(r; \psi))^2, \quad (61) $$

where $\psi = (a_0, \sigma^2, \ell, \nu)$. Although the nugget $a_0$ is estimated by variogram fitting, it is not useful for the considered problem, as it does not appear in the prior expression. A typical fitting method will yield non-integer values for $\nu$. In order to simplify numerical implementation, it is useful to approximate it with the nearest integer values (in 2D) or the nearest half-integer (in 3D). In order to preserve the correlation length, it is then necessary to update the identified length scale accordingly. More precisely, if $\ell_{\text{opt}}, \nu_{\text{opt}}$ are identified hyperparameters, and $\lfloor \nu \rfloor_d$, $d = 2, 3$ denotes the nearest integer or half-integer, then choosing

$$ \tilde{\ell}_{\text{opt}} = \ell_{\text{opt}} \sqrt{\frac{\nu_{\text{opt}}}{\lfloor \nu \rfloor_d}} \quad (62) $$
will preserve the correlation length. Using the pair $(\tilde{\ell}_{\text{opt}}, [\nu]_d)$, it is then possible to use formula (53) for an efficient prior precision discretization. With the variogram fitting approach, the posterior for the hierarchical problem is then approximated with an alternate optimization approach described in algorithm 1.

<caption>**Algorithm 1:** Hyperparameter Identification Procedure.</caption>

1 **Function IdentifyHyperparameters**($\sigma_0^2$, $\ell_0$, $\mathbf{y}^{\text{obs}}$, $\mathbf{p}_0$, $\varepsilon_{\text{tol}}$, $\eta_{\text{tol}}$):
2 &nbsp;&nbsp;&nbsp;**$\mu_{\text{opt}} \leftarrow (\sigma_0^2, \ell_0)$; $\mu \leftarrow \mu_{\text{opt}}$** // Initial hyperparameters guess
3 &nbsp;&nbsp;&nbsp;**$\mathbf{p}_{\text{opt}} \leftarrow \mathbf{p}_0$**
4 &nbsp;&nbsp;&nbsp;**$r \leftarrow 2 * \varepsilon_{\text{tol}}$**
5 &nbsp;&nbsp;&nbsp;**while** $r \ge \varepsilon_{\text{tol}}$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mathbf{p}_0 \leftarrow \mathbf{p}_{\text{opt}}$**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mathbf{p}_{\text{opt}} \leftarrow \text{IdentifyParameter}(\mathbf{y}^{\text{obs}}, \mathbf{p}_0, \eta_{\text{tol}})$** // Algorithm 2
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**($a_0, \sigma^2, \ell, \nu$) $\leftarrow \text{argmin}_{\psi} J_{\text{vario}}(\psi, \mathbf{p}_k)$** // Variogram regression
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\tilde{\ell} \leftarrow \ell \sqrt{\nu/[\nu]_d}$**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mu_{\text{opt}} \leftarrow (\sigma^2, \tilde{\ell}, [\nu]_d)$** // Rounding of smoothness parameter
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$r \leftarrow |\mu_{\text{opt}} - \mu|$**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mu \leftarrow \mu_{\text{opt}}$**
13 &nbsp;&nbsp;&nbsp;**return** $\mu_{\text{opt}}$**

For hyperparameter initialization, the authors of [Brown et al., 2020] recommend fitting the Matérn variogram on the observed field. However, in the present case, due to the non-linear relation between the stiffness parameter $\mathbf{p}$ and the displacement $\mathbf{u}$, the observed displacement field has a very different spatial correlation structure than the material parameter. We therefore recommend initializing $\mu_0$ to an arbitrary value. Note that with the variogram regression approach, the issues related to the identifiability of $(\sigma^2, \ell)$ remain present, as the marginal variance estimated by fitting the empirical variogram is simply the empirical variance of the observed sample and does not converge to the true marginal variance of the underlying Gaussian process when refining the computational mesh. The only difference when compared to the hybrid-MAP approach with the empirical variance as the hyperprior mean is that the smoothness parameter $\nu$ is estimated from the data as opposed to set to an arbitrary value, and the discrete prior precision matrix is approximated for the nearest integer (or half integer) value of $\nu$.

# 7 Comparison with a Classical Regularization Approach

Instead of a Bayesian formulation as developed in section 5.2, a classical regularization approach could have been used instead. In this case, a similar objective function could be defined:

$$J^c(\mathbf{p}; \gamma_R) = (\mathbf{u}^{\text{obs}} - \mathbf{u})^T (\mathbf{u}^{\text{obs}} - \mathbf{u}) + \gamma_F (F^{\text{obs}} - F)^2 + \gamma_R \|\mathbf{p}\|_L^2, \quad (63)$$

where $\|\cdot\|_L$ is a regularization norm, $\gamma_F$ a weight between the two data terms, and $\gamma_R$ a regularization parameter. The weight $\gamma_F$ is usually fixed, and a common choice is to set $\gamma_F = \sigma_u^2/\sigma_F^2$, where $\sigma_u^2$ is the average pointwise measurement displacement noise variance, so that each data term is weighted with the corresponding measurement noise. A common choice for $\|\cdot\|$ is the so-called Tikhonov regularization which is simply the $L^2$ norm, that writes in a discrete setting:

$$\|\mathbf{p}\|_{L^2}^2 = \mathbf{p}^T M \mathbf{p}. \quad (64)$$

However, in order to enforce a smooth solution, a more relevant choice is the $H^1(\Omega)$ semi-norm:

$$\|\mathbf{p}\|_{H^1(\Omega)}^2 = \int_{\Omega} (\nabla \mathbf{p})^2 \, dx, \quad (65)$$

which takes the following form in a discrete setting:

$$\|\mathbf{p}\|_{H^1(\Omega)}^2 = \mathbf{p}^T K \mathbf{p}. \quad (66)$$
will preserve the correlation length. Using the pair $(\tilde{\ell}_{\text{opt}}, \lfloor \nu \rfloor_d)$, it is then possible to use formula (53) for an efficient prior precision discretization. With the variogram fitting approach, the posterior for the hierarchical problem is then approximated with an alternate optimization approach described in algorithm 1.

<caption>**Algorithm 1:** Hyperparameter Identification Procedure.</caption>

1 **Function IdentifyHyperparameters**($\sigma_0^2$, $\ell_0$, $\mathbf{y}^{\text{obs}}$, $\mathbf{p}_0$, $\varepsilon_{\text{tol}}$, $\eta_{\text{tol}}$):
2 &nbsp;&nbsp;&nbsp;**$\mu_{\text{opt}} \leftarrow (\sigma_0^2, \ell_0)$; $\mu \leftarrow \mu_{\text{opt}}$** // Initial hyperparameters guess
3 &nbsp;&nbsp;&nbsp;**$\mathbf{p}_{\text{opt}} \leftarrow \mathbf{p}_0$**
4 &nbsp;&nbsp;&nbsp;**$r \leftarrow 2 * \varepsilon_{\text{tol}}$**
5 &nbsp;&nbsp;&nbsp;**while** $r \ge \varepsilon_{\text{tol}}$ **do**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mathbf{p}_0 \leftarrow \mathbf{p}_{\text{opt}}$**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mathbf{p}_{\text{opt}} \leftarrow \text{IdentifyParameter}(\mathbf{y}^{\text{obs}}, \mathbf{p}_0, \eta_{\text{tol}})$** // Algorithm 2
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$(a_0, \sigma^2, \ell, \nu) \leftarrow \text{argmin}_{\psi} J_{\text{vario}}(\psi, \mathbf{p}_k)$** // Variogram regression
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\tilde{\ell} \leftarrow \ell \sqrt{\nu / \lfloor \nu \rfloor_d}$**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mu_{\text{opt}} \leftarrow (\sigma^2, \tilde{\ell}, \lfloor \nu \rfloor_d)$** // Rounding of smoothness parameter
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$r \leftarrow |\mu_{\text{opt}} - \mu|$**
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**$\mu \leftarrow \mu_{\text{opt}}$**
13 &nbsp;&nbsp;&nbsp;**return** $\mu_{\text{opt}}$**

For hyperparameter initialization, the authors of [Brown et al., 2020] recommend fitting the Matérn variogram on the observed field. However, in the present case, due to the non-linear relation between the stiffness parameter $\mathbf{p}$ and the displacement $\mathbf{u}$, the observed displacement field has a very different spatial correlation structure than the material parameter. We therefore recommend initializing $\mu_0$ to an arbitrary value. Note that with the variogram regression approach, the issues related to the identifiability of $(\sigma^2, \ell)$ remain present, as the marginal variance estimated by fitting the empirical variogram is simply the empirical variance of the observed sample and does not converge to the true marginal variance of the underlying Gaussian process when refining the computational mesh. The only difference when compared to the hybrid-MAP approach with the empirical variance as the hyperprior mean is that the smoothness parameter $\nu$ is estimated from the data as opposed to set to an arbitrary value, and the discrete prior precision matrix is approximated for the nearest integer (or half integer) value of $\nu$.

# 7 Comparison with a Classical Regularization Approach

Instead of a Bayesian formulation as developed in section 5.2, a classical regularization approach could have been used instead. In this case, a similar objective function could be defined:

$$J^c(\mathbf{p}; \gamma_R) = (\mathbf{u}^{\text{obs}} - \mathbf{u})^T (\mathbf{u}^{\text{obs}} - \mathbf{u}) + \gamma_F (F^{\text{obs}} - F)^2 + \gamma_R \|\mathbf{p}\|_L^2, \quad (63)$$

where $\|\cdot\|_L$ is a regularization norm, $\gamma_F$ a weight between the two data terms, and $\gamma_R$ a regularization parameter. The weight $\gamma_F$ is usually fixed, and a common choice is to set $\gamma_F = \sigma_u^2/\sigma_F^2$, where $\sigma_u^2$ is the average pointwise measurement displacement noise variance, so that each data term is weighted with the corresponding measurement noise. A common choice for $\|\cdot\|$ is the so-called Tikhonov regularization which is simply the $L^2$ norm, that writes in a discrete setting:

$$\|\mathbf{p}\|_{L^2}^2 = \mathbf{p}^T M \mathbf{p}. \quad (64)$$

However, in order to enforce a smooth solution, a more relevant choice is the $H^1(\Omega)$ semi-norm:

$$\|\mathbf{p}\|_{H^1(\Omega)}^2 = \int_{\Omega} (\nabla \mathbf{p})^2 \, dx, \quad (65)$$

which takes the following form in a discrete setting:

$$\|\mathbf{p}\|_{H^1(\Omega)}^2 = \mathbf{p}^T K \mathbf{p}. \quad (66)$$
Compared to the Bayesian formulation, there is only a single regularization parameter, which must be tuned. A common approach is to select $\gamma_R$ with an approach based on the L-curve, a parametric curve defined as:

$$ \gamma_R \mapsto (J_{\text{data}}^c(\gamma_R), J_{\text{reg}}^c(\gamma_R)), \qquad (67) $$

where

$$ J_{\text{data}}^{c}(\gamma_{R}) = (\mathbf{u}^{\text{obs}} - \mathbf{u}(\mathbf{p}_{\gamma_{R}}))^{T}(\mathbf{u}^{\text{obs}} - \mathbf{u}(\mathbf{p}_{\gamma_{R}})) + \gamma_{F}(F^{\text{obs}} - F(\mathbf{p}_{\gamma_{R}}))^{2}, \qquad (68) $$

$$ J_{\text{reg}}^{c}(\gamma_{R}) = \|\mathbf{p}_{\gamma_{R}}\|_{L}, \qquad (69) $$

$$ \mathbf{p}_{\gamma_R} = \underset{\mathbf{p}}{\operatorname{argmin}} \quad J^c(\mathbf{p}; \gamma_R). \qquad (70) $$

The L-curve represents the tradeoff between the regularized solution norm and the residual norm. For more background and details, we refer to [Hansen, 1992] (and [Gulliksson, 1998] for application to nonlinear problems). A common method to select the regularization parameter is to choose the value that corresponds to the L-curve point of maximum curvature, which yields the value satisfying the Morozov discrepancy principle [Engl and Grever, 1994]. In practice, the point of maximum curvature can be obtained graphically, or by the means of more sophisticated algorithms that aim at reducing the number of L-curve points to construct (see, e.g., [Cultrera and Callegaro, 2020]).

We note that even though an analogy between classical and Bayesian inverse problems can be formulated [Polson and Sokolov, 2019], the Whittle-Matérn prior formulated in section 6 is not equivalent to a standard Tikhonov regularization. In particular, the single regularization parameter $\gamma_R$ does not allow controlling the amplitude and length scale of the solution separately, as it is the case with the Matérn covariance parameters ($\sigma^2, \rho$).

In section 9, we show an example of identified parameter and regularization parameter through the classical approach. Although it leads to satisfactory results for parameter estimation, we lose the statistical interpretation of the regularization parameter and the sampling capabilities from the tuned covariance matrix. This is why the presented framework based on Gaussian random fields modelling and the Bayesian inverse problem formulation is considered more relevant for modelling spatially varying material parameters.

# 8 Numerical Optimization : Adjoint Formulation and Numerical Scheme

In this section, we discuss the numerical approximation of the MAP estimator of the unknown parameter defined in section 5. We first present how to efficiently compute the loss function gradient through the adjoint method, and then discuss the optimization algorithm used to derive MAP estimates.

## 8.1 Efficient Gradient Computation with the Adjoint Method

In typical parameter identification procedures such as FEMU or I-DIC, the parameter sensitivities are approximated with a finite difference scheme. This is tractable in the case where there are few parameters to identify. In the present case, as we introduced a finite elements discretization of the unknown parameter $p$, there is a total number of $N_V$ degrees of freedom to identify, which scales with mesh size. In this case a finite differences approximation of $\nabla J$ is usually intractable for larger meshes. The adjoint method introduces a so-called *adjoint problem* which can be solved using the finite elements method and provides an efficient way to compute the gradient. The adjoint method has been used in a wide variety of settings to compute gradients in PDE-constrained optimization problems [Richter, 2021, Giannakoglou and Papadimitriou, 2008, Asch et al., 2016], however the method is often presented in a problem specific way and is not always straightforward to adapt to another problem. For clarity, we derive the equations in the specific case of problem (8)-(11). Reusing the notations introduced in section 5.1, given a displacement field $\underline{u} \in H_D^1(\Omega^{\text{obs}})$ and a stiffness parameter field $p \in L^\infty(\Omega^{\text{obs}})$, we define the weak residual as the linear form

$$ \underline{v} \mapsto R(\underline{u}, \underline{v}, p) = l(\underline{v}) - a(\underline{u}, \underline{v}, p), \qquad (71) $$
Compared to the Bayesian formulation, there is only a single regularization parameter, which must be tuned. A common approach is to select $\gamma_R$ with an approach based on the L-curve, a parametric curve defined as:

$$ \gamma_R \mapsto (J_{\text{data}}^c(\gamma_R), J_{\text{reg}}^c(\gamma_R)), \qquad (67) $$

where

$$ J_{\text{data}}^{c}(\gamma_{R}) = (\boldsymbol{u}^{\text{obs}} - \boldsymbol{u}(\boldsymbol{p}_{\gamma_{R}}))^{T}(\boldsymbol{u}^{\text{obs}} - \boldsymbol{u}(\boldsymbol{p}_{\gamma_{R}})) + \gamma_{F}(F^{\text{obs}} - F(\boldsymbol{p}_{\gamma_{R}}))^{2}, \qquad (68) $$

$$ J_{\text{reg}}^{c}(\gamma_{R}) = \|\boldsymbol{p}_{\gamma_{R}}\|_{L}, \qquad (69) $$

$$ \boldsymbol{p}_{\gamma_R} = \underset{\boldsymbol{p}}{\operatorname{argmin}} \quad J^c(\boldsymbol{p}; \gamma_R). \qquad (70) $$

The L-curve represents the tradeoff between the regularized solution norm and the residual norm. For more background and details, we refer to [Hansen, 1992] (and [Gulliksson, 1998] for application to nonlinear problems). A common method to select the regularization parameter is to choose the value that corresponds to the L-curve point of maximum curvature, which yields the value satisfying the Morozov discrepancy principle [Engl and Grever, 1994]. In practice, the point of maximum curvature can be obtained graphically, or by the means of more sophisticated algorithms that aim at reducing the number of L-curve points to construct (see, e.g., [Cultrera and Callegaro, 2020]).

We note that even though an analogy between classical and Bayesian inverse problems can be formulated [Polson and Sokolov, 2019], the Whittle-Matérn prior formulated in section 6 is not equivalent to a standard Tikhonov regularization. In particular, the single regularization parameter $\gamma_R$ does not allow controlling the amplitude and length scale of the solution separately, as it is the case with the Matérn covariance parameters ($\sigma^2, \rho$).

In section 9, we show an example of identified parameter and regularization parameter through the classical approach. Although it leads to satisfactory results for parameter estimation, we lose the statistical interpretation of the regularization parameter and the sampling capabilities from the tuned covariance matrix. This is why the presented framework based on Gaussian random fields modelling and the Bayesian inverse problem formulation is considered more relevant for modelling spatially varying material parameters.

# 8 Numerical Optimization : Adjoint Formulation and Numerical Scheme

In this section, we discuss the numerical approximation of the MAP estimator of the unknown parameter defined in section 5. We first present how to efficiently compute the loss function gradient through the adjoint method, and then discuss the optimization algorithm used to derive MAP estimates.

## 8.1 Efficient Gradient Computation with the Adjoint Method

In typical parameter identification procedures such as FEMU or I-DIC, the parameter sensitivities are approximated with a finite difference scheme. This is tractable in the case where there are few parameters to identify. In the present case, as we introduced a finite elements discretization of the unknown parameter $p$, there is a total number of $N_V$ degrees of freedom to identify, which scales with mesh size. In this case a finite differences approximation of $\nabla J$ is usually intractable for larger meshes. The adjoint method introduces a so-called *adjoint problem* which can be solved using the finite elements method and provides an efficient way to compute the gradient. The adjoint method has been used in a wide variety of settings to compute gradients in PDE-constrained optimization problems [Richter, 2021, Giannakoglou and Papadimitriou, 2008, Asch et al., 2016], however the method is often presented in a problem specific way and is not always straightforward to adapt to another problem. For clarity, we derive the equations in the specific case of problem (8)-(11). Reusing the notations introduced in section 5.1, given a displacement field $\underline{u} \in H_D^1(\Omega^{\text{obs}})$ and a stiffness parameter field $p \in L^\infty(\Omega^{\text{obs}})$, we define the *weak residual* as the linear form

$$ \underline{v} \mapsto R(\underline{u}, \underline{v}, p) = l(\underline{v}) - a(\underline{u}, \underline{v}, p), \qquad (71) $$
where we made the dependency with respect to $p$ appear explicitly in the notation. We consider an arbitrary scalar valued loss function $J(\underline{u}, p) : H_D^1(\Omega^{\text{obs}}) \times L^\infty(\Omega^{\text{obs}}) \to \mathbb{R}$ that is assumed to be Gateaux differentiable. With these notations, the parameter identification problem can be written as follows:

$$ (\underline{u}^{\text{opt}}, p^{\text{opt}}) = \underset{(\underline{u}, p) \in H^1(\Omega^{\text{obs}}) \times L^\infty(\Omega^{\text{obs}})}{\operatorname{argmin}} \quad J(\underline{u}, p) \qquad (72) $$

$$ \text{such that } \forall \underline{v} \in H_0^1(\Omega^{\text{obs}}), R(\underline{u}^{\text{opt}}, \underline{v}, p^{\text{opt}}) = 0. \qquad (73) $$

The PDE constraint (73) introduces a dependency between $\underline{u}$ and $p$ and $J$ can thus be seen as a composed function $p \mapsto J(\underline{u}(p), p)$. We now consider $p, \delta p \in L^\infty(\Omega^{\text{obs}})$ and $\underline{u} \in H_D^1(\Omega^{\text{obs}})$, and expand the Gateaux derivative of $J$ and $R$ at point $p$ in the direction $\delta p$, using the chain rule:

$$ dJ(\underline{u}(p), p; \delta p) = \partial_1 J(\underline{u}(p), p; d\underline{u}(p; \delta p)) + \partial_2 J(\underline{u}(p), p; \delta p), \qquad (74) $$

$$ \forall \underline{v} \in H_0^1(\Omega^{\text{obs}}), dR(\underline{u}(p), \underline{v}, p; \delta p) = \partial_1 R(\underline{u}(p), \underline{v}, p; d\underline{u}(p; \delta p)) + \partial_3 R(\underline{u}(p), \underline{v}, p; \delta p). \qquad (75) $$

Using the constraint (73), we can write $dR(\underline{u}(p), \underline{v}, p, \delta p) = 0$ for all $\underline{v} \in H_0^1(\Omega^{\text{obs}})$, and then

$$ \partial_3 R(\underline{u}(p), \underline{v}, p; \delta p) = -\partial_1 R(\underline{u}(p), \underline{v}, p; d\underline{u}(p; \delta p)). \qquad (76) $$

Now, let's consider $\underline{\lambda}$, the unique solution to the variational adjoint problem: Find $\underline{\lambda} \in H_0^1(\Omega^{\text{obs}})$ such that

$$ \forall \delta \underline{u} \in H^1(\Omega^{\text{obs}}), \partial_1 R(\underline{u}(p), \underline{\lambda}, p; \delta u) = \partial_1 J(\underline{u}(p), p; \delta u). \qquad (77) $$

Combining (77) with $\delta \underline{u} = d\underline{u}(p; \delta p)$ and (76) with $\underline{v} = \underline{\lambda}$ inside (74) yields

$$ dJ(\underline{u}(p), p; \delta p) = \partial_2 J(\underline{u}(p), p; \delta p) - \partial_3 R(\underline{u}, \underline{\lambda}, p; \delta p). \qquad (78) $$

Typically, both derivatives in (78) can be computed analytically, provided the adjoint solution $\underline{\lambda}$ is available. We also note that in the present case, as the operator $(\underline{u}, \underline{v}) \mapsto a(\underline{u}, \underline{v}, p)$ is bilinear and self-adjoint, the adjoint equation (77) reduces to

$$ a(\delta \underline{u}, \underline{\lambda}, p) = \partial_1 J(\underline{u}(p), p; \delta u). \qquad (79) $$

Now, let's apply the parameter sensitivity formula (78) to the discretized Bayesian formulation introduced in section 5.2. We recall that the Bayesian loss is composed of three terms $J(\underline{p}) = J_u(\underline{p}) + J_F(\underline{p}) + J_p(\underline{p})$. For the third term, computing the gradient is straightforward:

$$ \nabla J_p(\underline{p}) = (\underline{p} - \underline{p}^0)^T \underline{\Sigma}_p^{-1}. \qquad (80) $$

For $\nabla J_u$ and $\nabla J_F$, we apply the adjoint method. First, we note $\underline{A}$ the matrix of operator $a$ in basis $(\underline{\varphi}_k)_{0 \le k < N_{\underline{V}}}$:

$$ \underline{A} = (a(\underline{\varphi}_i, \underline{\varphi}_j))_{0 \le i, j < N_{\underline{V}}}. \qquad (81) $$

To keep notations simple, we assume that basis functions are ordered in such a way that $\underline{A}$ can be written in the following blocked form

$$ \underline{A} = \begin{bmatrix} \underline{A}_{DD} & \underline{A}_{DU} \\ \underline{A}_{DU} & \underline{A}_{UU} \end{bmatrix}, \qquad (82) $$

where $D$ is the set of DOF indices located inside $\partial \Omega_D^{\text{obs}}$ and $U$ is the set of DOF indices located inside $\Omega^{\text{obs}} \setminus \partial \Omega_D^{\text{obs}}$. The linear system satisfied by $\underline{u}$ given boundary conditions (11) writes

$$ \underline{A}_0 \underline{u} = \underline{b}, \qquad (83) $$

where

$$ \underline{A}_0 = \begin{bmatrix} I_{DD} & \underline{0} \\ \underline{0} & \underline{A}_{UU} \end{bmatrix}, \quad \underline{b} = \begin{bmatrix} \underline{u}_D^{\text{obs}} \\ -\underline{A}_{UD}\underline{u}_D^{\text{obs}} \end{bmatrix}, \qquad (84) $$
where we made the dependency with respect to $p$ appear explicitly in the notation. We consider an arbitrary scalar valued loss function $J(\underline{u}, p) : H_D^1(\Omega^{\text{obs}}) \times L^\infty(\Omega^{\text{obs}}) \to \mathbb{R}$ that is assumed to be Gateaux differentiable. With these notations, the parameter identification problem can be written as follows:

$$ (\underline{u}^{\text{opt}}, p^{\text{opt}}) = \underset{(\underline{u}, p) \in H^1(\Omega^{\text{obs}}) \times L^\infty(\Omega^{\text{obs}})}{\operatorname{argmin}} \quad J(\underline{u}, p) \quad (72) $$

$$ \text{such that } \forall \underline{v} \in H_0^1(\Omega^{\text{obs}}), R(\underline{u}^{\text{opt}}, \underline{v}, p^{\text{opt}}) = 0. \quad (73) $$

The PDE constraint (73) introduces a dependency between $\underline{u}$ and $p$ and $J$ can thus be seen as a composed function $p \mapsto J(\underline{u}(p), p)$. We now consider $p, \delta p \in L^\infty(\Omega^{\text{obs}})$ and $\underline{u} \in H_D^1(\Omega^{\text{obs}})$, and expand the Gateaux derivative of $J$ and $R$ at point $p$ in the direction $\delta p$, using the chain rule:

$$ dJ(\underline{u}(p), p; \delta p) = \partial_1 J(\underline{u}(p), p; d\underline{u}(p; \delta p)) + \partial_2 J(\underline{u}(p), p; \delta p), \quad (74) $$

$$ \forall \underline{v} \in H_0^1(\Omega^{\text{obs}}), dR(\underline{u}(p), \underline{v}, p; \delta p) = \partial_1 R(\underline{u}(p), \underline{v}, p; d\underline{u}(p; \delta p)) + \partial_3 R(\underline{u}(p), \underline{v}, p; \delta p). \quad (75) $$

Using the constraint (73), we can write $dR(\underline{u}(p), \underline{v}, p, \delta p) = 0$ for all $\underline{v} \in H_0^1(\Omega^{\text{obs}})$, and then

$$ \partial_3 R(\underline{u}(p), \underline{v}, p; \delta p) = -\partial_1 R(\underline{u}(p), \underline{v}, p; d\underline{u}(p; \delta p)). \quad (76) $$

Now, let's consider $\underline{\lambda}$, the unique solution to the variational adjoint problem: Find $\underline{\lambda} \in H_0^1(\Omega^{\text{obs}})$ such that

$$ \forall \delta \underline{u} \in H^1(\Omega^{\text{obs}}), \partial_1 R(\underline{u}(p), \underline{\lambda}, p; \delta u) = \partial_1 J(\underline{u}(p), p; \delta u). \quad (77) $$

Combining (77) with $\delta \underline{u} = d\underline{u}(p; \delta p)$ and (76) with $\underline{v} = \underline{\lambda}$ inside (74) yields

$$ dJ(\underline{u}(p), p; \delta p) = \partial_2 J(\underline{u}(p), p; \delta p) - \partial_3 R(\underline{u}, \underline{\lambda}, p; \delta p). \quad (78) $$

Typically, both derivatives in (78) can be computed analytically, provided the adjoint solution $\underline{\lambda}$ is available. We also note that in the present case, as the operator $(\underline{u}, \underline{v}) \mapsto a(\underline{u}, \underline{v}, p)$ is bilinear and self-adjoint, the adjoint equation (77) reduces to

$$ a(\delta \underline{u}, \underline{\lambda}, p) = \partial_1 J(\underline{u}(p), p; \delta u). \quad (79) $$

Now, let's apply the parameter sensitivity formula (78) to the discretized Bayesian formulation introduced in section 5.2. We recall that the Bayesian loss is composed of three terms $J(\underline{p}) = J_u(\underline{p}) + J_F(\underline{p}) + J_p(\underline{p})$. For the third term, computing the gradient is straightforward:

$$ \nabla J_p(\underline{p}) = (\underline{p} - \underline{p}^0)^T \underline{\Sigma}_p^{-1}. \quad (80) $$

For $\nabla J_u$ and $\nabla J_F$, we apply the adjoint method. First, we note $\underline{A}$ the matrix of operator $a$ in basis $(\underline{\varphi}_k)_{0 \le k < N_{\underline{V}}}$:

$$ \underline{A} = \left( a(\underline{\varphi}_i, \underline{\varphi}_j) \right)_{0 \le i, j < N_{\underline{V}}}. \quad (81) $$

To keep notations simple, we assume that basis functions are ordered in such a way that $\underline{A}$ can be written in the following blocked form

$$ \underline{A} = \begin{bmatrix} \underline{A}_{DD} & \underline{A}_{DU} \\ \underline{A}_{DU} & \underline{A}_{UU} \end{bmatrix}, \quad (82) $$

where $D$ is the set of DOF indices located inside $\partial \Omega_D^{\text{obs}}$ and $U$ is the set of DOF indices located inside $\Omega^{\text{obs}} \setminus \partial \Omega_D^{\text{obs}}$. The linear system satisfied by $\underline{u}$ given boundary conditions (11) writes

$$ \underline{A}_0 \underline{u} = \underline{b}, \quad (83) $$

where

$$ \underline{A}_0 = \begin{bmatrix} I_{DD} & \underline{0} \\ \underline{0} & \underline{A}_{UU} \end{bmatrix}, \quad \underline{b} = \begin{bmatrix} \underline{u}_D^{\text{obs}} \\ -\underline{A}_{UD}\underline{u}_D^{\text{obs}} \end{bmatrix}, \quad (84) $$
and $u_D^{\text{obs}}$ is a vector collecting degrees of freedom of $\underline{u}^{\text{obs}}$ over $\partial\Omega_D^{\text{obs}}$. Constructing $A_0$ and $b$ is a common way to apply non-homogeneous Dirichlet boundary conditions to a linear system. As specified in problem (77), the adjoint solution $\lambda$ is sought in $H_0^1(\Omega^{\text{obs}})$ and therefore satisfies homogeneous Dirichlet boundary conditions. The linear system for the discretized adjoint solution $\lambda$ then writes

$$A_0\lambda = -b_{\text{adj}}, \qquad (85)$$

where

$$b_{\text{adj}} = \left[ \left( \frac{\mathbf{0}_D}{\partial \mathbf{u}} \right)_U \right], \quad \frac{\partial J}{\partial \mathbf{u}} = -(\mathbf{u}^{\text{obs}} - \mathbf{u})^T \mathbf{\Sigma}_{\mathbf{u}}^{-1} + \frac{\partial J_F}{\partial \mathbf{u}}. \qquad (86)$$

The expression for the gradient is then derived from (78), provided $\mathbf{u}$ and $\lambda$ have already been computed:

$$\nabla J(\mathbf{p}) = \frac{\partial J}{\partial \mathbf{p}} - \frac{\partial}{\partial \mathbf{p}} (\lambda^T \mathbf{R} \mathbf{u}) \qquad (87)$$

$$= (\mathbf{p} - \mathbf{p}^0)^T \mathbf{\Sigma}_{\mathbf{p}}^{-1} - \frac{\partial}{\partial \mathbf{p}} (\lambda^T \mathbf{R} \mathbf{u}), \qquad (88)$$

where

$$\mathbf{R} = \left( -a(\underline{\varphi}_i, \underline{\varphi}_j, p) \right)_{0 \le i,j < N_v}. \qquad (89)$$

Note that the second term of (88) is not straightforward to compute and depends on the stiffness model $p \mapsto \bar{C}(p)$. We show in section 8.2 how to leverage automatic differentiation to compute this term easily.

**Remark.** In the case where multiple stiffness parameter fields $p_1, \dots, p_{N_p}$ are sought, as derivatives with respect to $p_i$ are not required to construct the adjoint system, only two linear systems must be solved to compute $\nabla J(p_1, \dots, p_{N_p})$, regardless of the number of stiffness parameters.

## 8.2 Implementation

The FEMU loss function is minimized with a quasi-Newton method. As opposed to most FEMU implementations, e.g. [Leclerc et al., 2009, Neggers et al., 2017], the L-BFGS method [Liu and Nocedal, 1989] is used instead of Gauss-Newton, as the latter requires computing the Jacobian matrix of the residual vector to approximate the hessian, and in our case storing the full Jacobian has a memory complexity of $O(N_V^2)$, which can make it intractable for larger meshes. Instead, L-BFGS requires only the loss gradient, and stores only a few columns of the hessian approximation, which keeps the overall memory complexity to $O(N_V)$. The iterative optimization procedure for parameter identification is summarized in algorithm 2, and the global hyperparameter identification procedure in algorithm 1.

Our implementation of the proposed method is written in Python using free and open source software only. For finite elements operators assembly, we use the DOLFINx library [Habera et al., 2020]. Using the FEiCSx ecosystem, both variational problems (18) and (77), as well as the stiffness model $p \mapsto \bar{C}(p)$ and its gradient can be modelled using the UFL language [Alnæs et al., 2014]. We provide in figure 3 a simplified version of the code to compute the PDE residual gradient and the gradient of the macroscopic force misfit $\nabla J_F$, which are both dependent on the stiffness model and thus not straightforward to differentiate analytically. Instead, by leveraging the automatic differentiation capabilities of UFL and FFCX, the required derivatives do not need to be expressed analytically. As the optimization loop involves solving two linear systems, we leverage the PETSc library [Balay et al., 2022] for all solves. A stabilized conjugate gradient solver is used alongside an algebraic multi-grid preconditioner. For hyperparameter estimation, the GSTools library is used [Müller et al., 2022, Müller and Schüler, 2023], which provides routines for empirical variogram estimation and variogram fitting.
and $u_D^{\text{obs}}$ is a vector collecting degrees of freedom of $\underline{u}^{\text{obs}}$ over $\partial\Omega_D^{\text{obs}}$. Constructing $A_0$ and $b$ is a common way to apply non-homogeneous Dirichlet boundary conditions to a linear system. As specified in problem (77), the adjoint solution $\lambda$ is sought in $H_0^1(\Omega^{\text{obs}})$ and therefore satisfies homogeneous Dirichlet boundary conditions. The linear system for the discretized adjoint solution $\lambda$ then writes

$$A_0\lambda = -b_{\text{adj}}, \qquad (85)$$

where

$$b_{\text{adj}} = \left[ \left( \frac{\mathbf{0}_D}{\partial \mathbf{u}} \right)_U \right], \quad \frac{\partial J}{\partial \mathbf{u}} = -(\mathbf{u}^{\text{obs}} - \mathbf{u})^T \Sigma_{\mathbf{u}}^{-1} + \frac{\partial J_F}{\partial \mathbf{u}}. \qquad (86)$$

The expression for the gradient is then derived from (78), provided $\mathbf{u}$ and $\lambda$ have already been computed:

$$\nabla J(\mathbf{p}) = \frac{\partial J}{\partial \mathbf{p}} - \frac{\partial}{\partial \mathbf{p}} (\lambda^T \mathbf{R} \mathbf{u}) \qquad (87)$$

$$= (\mathbf{p} - \mathbf{p}^0)^T \Sigma_{\mathbf{p}}^{-1} - \frac{\partial}{\partial \mathbf{p}} (\lambda^T \mathbf{R} \mathbf{u}), \qquad (88)$$

where

$$\mathbf{R} = \left( -a(\underline{\varphi}_i, \underline{\varphi}_j, p) \right)_{0 \le i,j < N_v}. \qquad (89)$$

Note that the second term of (88) is not straightforward to compute and depends on the stiffness model $p \mapsto \bar{C}(p)$. We show in section 8.2 how to leverage automatic differentiation to compute this term easily.

**Remark.** In the case where multiple stiffness parameter fields $p_1, \dots, p_{N_p}$ are sought, as derivatives with respect to $p_i$ are not required to construct the adjoint system, only two linear systems must be solved to compute $\nabla J(p_1, \dots, p_{N_p})$, regardless of the number of stiffness parameters.

## 8.2 Implementation

The FEMU loss function is minimized with a quasi-Newton method. As opposed to most FEMU implementations, e.g. [Leclerc et al., 2009, Neggers et al., 2017], the L-BFGS method [Liu and Nocedal, 1989] is used instead of Gauss-Newton, as the latter requires computing the Jacobian matrix of the residual vector to approximate the hessian, and in our case storing the full Jacobian has a memory complexity of $O(N_V^2)$, which can make it intractable for larger meshes. Instead, L-BFGS requires only the loss gradient, and stores only a few columns of the hessian approximation, which keeps the overall memory complexity to $O(N_V)$. The iterative optimization procedure for parameter identification is summarized in algorithm 2, and the global hyperparameter identification procedure in algorithm 1.

Our implementation of the proposed method is written in Python using free and open source software only. For finite elements operators assembly, we use the DOLFINx library [Habera et al., 2020]. Using the FEiCSx ecosystem, both variational problems (18) and (77), as well as the stiffness model $p \mapsto \bar{C}(p)$ and its gradient can be modelled using the UFL language [Alnæs et al., 2014]. We provide in figure 3 a simplified version of the code to compute the PDE residual gradient and the gradient of the macroscopic force misfit $\nabla J_F$, which are both dependent on the stiffness model and thus not straightforward to differentiate analytically. Instead, by leveraging the automatic differentiation capabilities of UFL and FFCX, the required derivatives do not need to be expressed analytically. As the optimization loop involves solving two linear systems, we leverage the PETSc library [Balay et al., 2022] for all solves. A stabilized conjugate gradient solver is used alongside an algebraic multi-grid preconditioner. For hyperparameter estimation, the GSTools library is used [Müller et al., 2022, Müller and Schüler, 2023], which provides routines for empirical variogram estimation and variogram fitting.
<caption>**Algorithm 2: Parameter Identification Procedure**</caption>

1 **Function IdentifyParameter(y<sup>obs</sup>, p<sub>0</sub>, η<sub>tol</sub>):**
2 &nbsp;&nbsp;&nbsp; p ← p<sub>0</sub> // Initial stiffness parameter guess
3 &nbsp;&nbsp;&nbsp; while J(p, y<sup>obs</sup>) > η<sub>tol</sub> do
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Assemble A<sub>0</sub>, b;
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Solve A<sub>0</sub>u = b // Forward problem
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Assemble dJ/du and b<sub>adj</sub>;
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Solve A<sub>0</sub>λ = b<sub>adj</sub> // Adjoint problem
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Assemble ∂J/∂p and ∂/∂p(λ<sup>T</sup> Ru);
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dJ/dp ← ∂J/∂p - ∂/∂p(λ<sup>T</sup> Ru);
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; δp ← BGFS (dJ/dp) // L-BFGS update formula
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; p ← p + δp;
12 &nbsp;&nbsp;&nbsp; return p

1 from typing import Callable
2 from ufl import *
3 from dolfinx.fem import Function, form, assemble_vector, assemble_scalar
4
5 def epsilon(u):
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; """strain"""
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return sym grad(u)
8
9 def sigma(u, C):
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; """stress"""
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; return inner(C, epsilon(u))
12
13 def dRdp(u_h: Function, lmbda_h: Function, p: Function, stiffness_model: Callable):
14 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; """
15 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Evaluate PDE residual sensitivity
16 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; u_h: forward problem solution
17 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; lmbda_h: adjoint problem solution
18 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; p: current model parameter values
19 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; stiffness_model: parameter to stiffness mapping in UFL
20 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; """
21 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; R = -dot(sigma(epsilon(u_h), stiffness_model(p)), epsilon(lmbda_h)) * dx
22 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dRdp = derivative(R, p)
23
24 # call FFCX JIT compiler and generate optimized C code
25 dRdp_form = form(dRdp)
26
27 # assemble DOF array
28 return assemble_vector(dRdp_form)
29
30 def dJFdp(u_h: Function, p: Function, Fobs: float, sigma_F: float, n_F, stiffness_model: Callable):
31 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; """Evaluate macroscopic force missfit"""
32 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; F = dot(n_F, dot(sigma(epsilon(u_h), stiffness_model(p)), n_F)) * dx
33 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dFdp = derivative(F, p)
34
35 F_eval = assemble_scalar(form(F))
36 dFdp_eval = assemble_vector(form(dFdp))
37
38 return (F_eval - Fobs) / sigma_F**2

<caption>Figure 3: Simplified implementation of residual and loss sensitivities using UFL and DOLFINx assemblers</caption>
<caption>Algorithm 2: Parameter Identification Procedure</caption>

1 **Function IdentifyParameter(y<sup>obs</sup>, p<sub>0</sub>, η<sub>tol</sub>):**
2 &nbsp;&nbsp;&nbsp;&nbsp;`p ← p<sub>0</sub>` // Initial stiffness parameter guess
3 &nbsp;&nbsp;&nbsp;&nbsp;**while** J(`p`, y<sup>obs</sup>) > η<sub>tol</sub> **do**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Assemble** A<sub>0</sub>, **b**;
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Solve** A<sub>0</sub>u = **b** // Forward problem
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Assemble** d/J/du and b<sub>adj</sub>;
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Solve** A<sub>0</sub>λ = b<sub>adj</sub> // Adjoint problem
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Assemble** ∂J/∂p and ∂/∂p(λ<sup>T</sup>Ru);
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;dJ/dp ← ∂J/∂p - ∂/∂p(λ<sup>T</sup>Ru);
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;δp ← BGFS(dJ/dp) // L-BFGS update formula
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;`p ← p + δp`;
12 &nbsp;&nbsp;&nbsp;&nbsp;**return** p

1 from typing import Callable
2 from ufl import *
3 from dolfinx.fem import Function, form, assemble_vector, assemble_scalar
4
5 def epsilon(u):
6 """strain"""
7 return sym(grad(u))
8
9 def sigma(u, C):
10 """stress"""
11 return inner(C, epsilon(u))
12
13 def dRdp(u_h: Function, lmbda_h: Function, p: Function, stiffness_model: Callable):
14 """
15 Evaluate PDE residual sensitivity
16 u_h: forward problem solution
17 lmbda_h: adjoint problem solution
18 p: current model parameter values
19 stiffness_model: parameter to stiffness mapping in UFL
20 """
21 R = -dot(sigma(epsilon(u_h), stiffness_model(p)), epsilon(lmbda_h)) * dx
22 dRdp = derivative(R, p)
23
24 # call FFCX JIT compiler and generate optimized C code
25 dRdp_form = form(dRdp)
26
27 # assemble DOF array
28 return assemble_vector(dRdp_form)
29
30 def dJFdp(u_h: Function, p: Function, Fobs: float, sigma_F: float, n_F, stiffness_model: Callable):
31 """Evaluate macroscopic force missfit"""
32 F = dot(n_F, dot(sigma(epsilon(u_h), stiffness_model(p)), n_F)) * dx
33 dFdp = derivative(F, p)
34
35 F_eval = assemble_scalar(form(F))
36 dFdp_eval = assemble_vector(form(dFdp))
37
38 return (F_eval - Fobs) / sigma_F**2

<caption>Figure 3: Simplified implementation of residual and loss sensitivities using UFL and DOLFINx assemblers</caption>
# 9 Results and Discussion

To showcase our method, we consider two numerical experiments, one modelling standard DIC with surface observations, and one with volume observations.

## 9.1 Identification of a Young Modulus from Surface Displacement

The first model problem is posed on the thin dog-bone specimen represented in figure 1. We consider an isotropic stiffness model, with a spatially varying Young modulus $E$, and a constant Poisson ratio $\nu_r$ (we employ the $r$ subscript to avoid confusion with $\nu$ that denotes the Matérn smoothness hyperparameter). The constitutive law 9 writes:

$$ \underline{\underline{\sigma}} = \frac{E\nu_r}{(1 + \nu_r)(1 - 2\nu_r)} \text{tr}(\underline{\underline{\varepsilon}})I_3 + \frac{E}{1 + \nu_r}\underline{\underline{\varepsilon}}. \quad (90) $$

The considered Young modulus is a realisation of a log-Gaussian random field with a Matérn covariance

$$ \mathbf{p} := \ln \mathbf{E} \sim \mathcal{N}(\ln E^0, \mathbf{\Sigma}(\sigma^2, \ell, \nu)), \quad (91) $$

with fixed smoothness parameter $\nu = 1$. A tensile test is simulated on $\Omega$, where the traction direction is represented by the unit vector $\underline{n}_F$. A uniform stress boundary condition $\underline{\underline{\sigma}} \cdot \underline{n} = s_0 \underline{n}_F$ is applied on $\partial\Omega_D^{\text{obs}}$ (in this case, $\underline{n} = \underline{n}_F$), and the surface displacement $\underline{u}$ is simulated with the finite element method. The observed DIC displacement is simulated by adding some noise $\underline{\underline{\xi}}$ to the simulated displacement degrees of freedom $\mathbf{u}^{\text{obs}} = \mathbf{u} + \underline{\underline{\xi}}$. For simplicity, we consider an uncorrelated Gaussian white noise

$$ \underline{\underline{\xi}} \sim \mathcal{N}(0, \sigma_u^2 \mathbf{I}), \quad (92) $$

and the measured force is calculated from the simulated stress field with formula (13) and corrupt with Gaussian white noise $\eta \sim \mathcal{N}(0, \sigma_F^2)$. The generated data is represented in figure 4. In order to study the impact of measurement noise on the identified parameter and hyperparameters, we define the observed signal-to-noise ratio (SNR) for the displacement and macroscopic force as

$$ \text{SNR}_{\underline{u}}^{\text{dB}}(\underline{u}, \sigma_u) = 20 \log_{10} \left( \frac{1}{N_{\text{obs}} \sigma_u} \sum_{k=1}^{N_{\text{obs}}} m_k \right), \quad (93) $$

$$ m_k = \sqrt{(u_{3k}^{\text{obs}})^2 + (u_{3k+1}^{\text{obs}})^2 + (u_{3k+2}^{\text{obs}})^2}, \quad (94) $$

$$ \text{SNR}_{F}^{\text{dB}}(F, \sigma_F) = 20 \log_{10} \left( \frac{F^{\text{obs}}}{\sigma_F} \right). \quad (95) $$

In the following results, we choose $\sigma_u$ and $\sigma_F$ such that $\text{SNR}_{\underline{u}}^{\text{dB}} \approx \text{SNR}_{F}^{\text{dB}}$. For a wide range of noise levels, we identify the MAP estimate for the posterior density, and the Matérn hyperparameters with the various PEB approaches presented in section 6.2: the ML estimator approach, the hybrid-MAP approach and the variogram regression approach. We recall that for the ML estimator and hybrid-MAP approaches, considering the approximations formulated in section 6.2, the discrete Whittle-Matérn prior is constructed by fixing the smoothness parameter to $\nu = 1$ (so that $\alpha = 2$). For the variogram regression approach, formula (53) is applied, using the identified smoothness parameter truncated to the nearest integer, as reported in algorithm 1. The MAP estimates $E_{\text{MAP}}$ and associated kinematic fields are showcased in figure 5 for decreasing noise levels. For high SNR values, the identified Young modulus field approximates very well the true field. For decreasing SNR values, $\mathbf{p}_{\text{MAP}}$ appears more smoothed out than the true field. Note that even though the true and simulated strain fields are represented in figures 4 and 5, no strain data is used for parameter or hyperparameter estimation. The reason is that the observed strain is computed from the noisy observed displacement field:

$$ \underline{\underline{\varepsilon}}^{\text{obs}} = \nabla (\underline{u} + \underline{\underline{\xi}}). \quad (96) $$
# 9 Results and Discussion

To showcase our method, we consider two numerical experiments, one modelling standard DIC with surface observations, and one with volume observations.

## 9.1 Identification of a Young Modulus from Surface Displacement

The first model problem is posed on the thin dog-bone specimen represented in figure 1. We consider an isotropic stiffness model, with a spatially varying Young modulus $E$, and a constant Poisson ratio $\nu_r$ (we employ the $r$ subscript to avoid confusion with $\nu$ that denotes the Matérn smoothness hyperparameter). The constitutive law 9 writes:

$$ \underline{\underline{\sigma}} = \frac{E\nu_r}{(1 + \nu_r)(1 - 2\nu_r)} \text{tr}(\underline{\underline{\varepsilon}})I_3 + \frac{E}{1 + \nu_r}\underline{\underline{\varepsilon}}. \quad (90) $$

The considered Young modulus is a realisation of a log-Gaussian random field with a Matérn covariance

$$ \mathbf{p} := \ln \mathbf{E} \sim \mathcal{N}(\ln E^0, \mathbf{\Sigma}(\sigma^2, \ell, \nu)), \quad (91) $$

with fixed smoothness parameter $\nu = 1$. A tensile test is simulated on $\Omega$, where the traction direction is represented by the unit vector $\underline{n}_F$. A uniform stress boundary condition $\underline{\underline{\sigma}} \cdot \underline{n} = s_0 \underline{n}_F$ is applied on $\partial\Omega_D^{\text{obs}}$ (in this case, $\underline{n} = \underline{n}_F$), and the surface displacement $\underline{u}$ is simulated with the finite element method. The observed DIC displacement is simulated by adding some noise $\underline{\underline{\xi}}$ to the simulated displacement degrees of freedom $\mathbf{u}^{\text{obs}} = \mathbf{u} + \underline{\underline{\xi}}$. For simplicity, we consider an uncorrelated Gaussian white noise

$$ \underline{\underline{\xi}} \sim \mathcal{N}(0, \sigma_u^2 \mathbf{I}), \quad (92) $$

and the measured force is calculated from the simulated stress field with formula (13) and corrupt with Gaussian white noise $\eta \sim \mathcal{N}(0, \sigma_F^2)$. The generated data is represented in figure 4. In order to study the impact of measurement noise on the identified parameter and hyperparameters, we define the observed signal-to-noise ratio (SNR) for the displacement and macroscopic force as

$$ \text{SNR}_{\underline{u}}^{\text{dB}}(\underline{u}, \sigma_u) = 20 \log_{10} \left( \frac{1}{N_{\text{obs}} \sigma_u} \sum_{k=1}^{N_{\text{obs}}} m_k \right), \quad (93) $$

$$ m_k = \sqrt{(u_{3k}^{\text{obs}})^2 + (u_{3k+1}^{\text{obs}})^2 + (u_{3k+2}^{\text{obs}})^2}, \quad (94) $$

$$ \text{SNR}_{F}^{\text{dB}}(F, \sigma_F) = 20 \log_{10} \left( \frac{F^{\text{obs}}}{\sigma_F} \right). \quad (95) $$

In the following results, we choose $\sigma_u$ and $\sigma_F$ such that $\text{SNR}_{\underline{u}}^{\text{dB}} \approx \text{SNR}_{F}^{\text{dB}}$. For a wide range of noise levels, we identify the MAP estimate for the posterior density, and the Matérn hyperparameters with the various PEB approaches presented in section 6.2: the ML estimator approach, the hybrid-MAP approach and the variogram regression approach. We recall that for the ML estimator and hybrid-MAP approaches, considering the approximations formulated in section 6.2, the discrete Whittle-Matérn prior is constructed by fixing the smoothness parameter to $\nu = 1$ (so that $\alpha = 2$). For the variogram regression approach, formula (53) is applied, using the identified smoothness parameter truncated to the nearest integer, as reported in algorithm 1. The MAP estimates $E_{\text{MAP}}$ and associated kinematic fields are showcased in figure 5 for decreasing noise levels. For high SNR values, the identified Young modulus field approximates very well the true field. For decreasing SNR values, $\mathbf{p}_{\text{MAP}}$ appears more smoothed out than the true field. Note that even though the true and simulated strain fields are represented in figures 4 and 5, no strain data is used for parameter or hyperparameter estimation. The reason is that the observed strain is computed from the noisy observed displacement field:

$$ \underline{\underline{\varepsilon}}^{\text{obs}} = \nabla (\underline{u} + \underline{\underline{\xi}}). \quad (96) $$
<figure><img src="image_3.png" /><figcaption>Figure 4: Generated data for model problem. Simulated displacement (a), simulated strain field (b), generated Matérn field $p = \ln E$ (c), corresponding Young modulus field $E = \exp(p)$ (d).</figcaption></figure>

Differentiating noisy data is an ill-posed problem [Knowles and Renka, 2014] that involves some kind of underlying regularization (e.g. Tikhonov regularization). In typical DIC software, a regularization parameter is usually set to an arbitrary value in order to compute and visualize the strain field. Performing covariance estimation from the observed strain fields would yield more information about that regularization parameter than the covariance properties of the underlying material parameter. Even in a low noise setting, as the strain field is fully computed from the displacement field, including strain data in the Bayesian modelling would not be more informative than using the displacement field only.

<figure><img src="image_4.png" /><figcaption>Figure 5: Plots of observed displacement (a), identified displacement (b), identified stiffness parameter $\ln E$ (c), identified strain (d) for decreasing noise levels.</figcaption></figure>

In figure 10, we compare the identified hyperparameters with respect to noise level for the three proposed PEB updating approaches. For the variogram regression and hybrid-MAP approaches, the identified marginal variance and correlation length stabilize close to the true values, whereas for the ML approach, the identified hyperparameters are erroneous even with a high SNR, despite the identified Young modulus field being close to the other PEB methods. These results can be interpreted as an illustration of the results of [Zhang, 2004]: when fixing $\nu$ (here, $\nu$ is set to 1), it is not possible to differentiate two hyperparameter pairs ($\sigma_0^2, \ell_0$) and ($\sigma_1^2, \ell_1$) if $\sigma_0^2/\ell_0 = \sigma_1^2/\ell_1$ with a ML estimator. In the hybrid-MAP approach, a hyperprior is specified on $\sigma^2$ only. Numerically, this is equivalent to estimate the microergodic hyperparameter $\sigma^2/\ell$, which is shown to be consistently estimable in [Zhang, 2004]. Numerically, the instabilities with the ML estimator approach translate to ridges
<figure><img src="image_4.png" /><figcaption>Figure 4: Generated data for model problem. Simulated displacement (a), simulated strain field (b), generated Matérn field $p = \ln E$ (c), corresponding Young modulus field $E = \exp(p)$ (d).</figcaption></figure>

Differentiating noisy data is an ill-posed problem [Knowles and Renka, 2014] that involves some kind of underlying regularization (e.g. Tikhonov regularization). In typical DIC software, a regularization parameter is usually set to an arbitrary value in order to compute and visualize the strain field. Performing covariance estimation from the observed strain fields would yield more information about that regularization parameter than the covariance properties of the underlying material parameter. Even in a low noise setting, as the strain field is fully computed from the displacement field, including strain data in the Bayesian modelling would not be more informative than using the displacement field only.

<figure><img src="image_5.png" /><figcaption>Figure 5: Plots of observed displacement (a), identified displacement (b), identified stiffness parameter $\ln E$ (c), identified strain (d) for decreasing noise levels.</figcaption></figure>

In figure 10, we compare the identified hyperparameters with respect to noise level for the three proposed PEB updating approaches. For the variogram regression and hybrid-MAP approaches, the identified marginal variance and correlation length stabilize close to the true values, whereas for the ML approach, the identified hyperparameters are erroneous even with a high SNR, despite the identified Young modulus field being close to the other PEB methods. These results can be interpreted as an illustration of the results of [Zhang, 2004]: when fixing $\nu$ (here, $\nu$ is set to 1), it is not possible to differentiate two hyperparameter pairs ($\sigma_0^2, \ell_0$) and ($\sigma_1^2, \ell_1$) if $\sigma_0^2/\ell_0 = \sigma_1^2/\ell_1$ with a ML estimator. In the hybrid-MAP approach, a hyperprior is specified on $\sigma^2$ only. Numerically, this is equivalent to estimate the microergodic hyperparameter $\sigma^2/\ell$, which is shown to be consistently estimable in [Zhang, 2004]. Numerically, the instabilities with the ML estimator approach translate to ridges
in the log-likelihood profile $\tau \mapsto J_p(p_{\text{MAP}}, \tau)$ (see, e.g., figure 1 of [Lalchand et al., 2022]), which means that the identified hyperparameters are also highly dependent on initialization values. In this example, as the true smoothness parameter is set to $\nu = 1$, the discretized prior (56) constructed with $\alpha = 2$ models the correct smoothness. In that case, given the results of figure 10, the hybrid-MAP and variogram regression approaches are equivalent.

As an illustration of the classical approach detailed in section 7, the Young modulus field is also identified through gradient penalty regularization, with a regularization parameter selected from the constructed *L*-curve. For a given SNR level, identified parameters for various values of $\gamma_R$ are represented in figure 6, alongside with the associated empirical *L*-curve. This figure illustrates that the solution is smoother as the value of $\gamma_R$ increases. Figure 7 illustrates selection of $\gamma_R$ through the

<figure><img src="image_5.png" /><figcaption>Figure 6: Representation of the *L*-curve for a fixed SNR level, and identified parameters for different values of $\gamma_R$.</figcaption></figure>

*L*-curve and associated identified parameter, for two different noise levels. As expected, the magnitude of regularization is increased when the noise level is increased. We can also compare the identified parameters of figure 5 and figure 7 to conclude that both the classical and Bayesian approach yield similar results for parameter identification, and that the main difference between the two approaches lies in the statistical interpretation of covariance parameters, and the sampling capabilities of the Bayesian approach.

<figure><img src="image_6.png" /><figcaption>Figure 7: Identified classical regularization parameters from two *L*-curves associated with two different SNR levels, and Young modulus field with optimal values for $\gamma_R$.</figcaption></figure>

From the identified hyperparameters, Matérn random fields with similar covariance properties can then be sampled. We illustrate this in figure 8, where new Matérn fields are sampled from the prior distribution with tuned covariance. The samples were generated using the Python library GSTools [Müller et al., 2022, Müller and Schüler, 2023]. These samples illustrate the impact of measurement noise on the spatial statistics of the resulting generative model.
in the log-likelihood profile $\tau \mapsto J_p(p_{\text{MAP}}, \tau)$ (see, e.g., figure 1 of [Lalchand et al., 2022]), which means that the identified hyperparameters are also highly dependent on initialization values. In this example, as the true smoothness parameter is set to $\nu = 1$, the discretized prior (56) constructed with $\alpha = 2$ models the correct smoothness. In that case, given the results of figure 10, the hybrid-MAP and variogram regression approaches are equivalent.

As an illustration of the classical approach detailed in section 7, the Young modulus field is also identified through gradient penalty regularization, with a regularization parameter selected from the constructed *L*-curve. For a given SNR level, identified parameters for various values of $\gamma_R$ are represented in figure 6, alongside with the associated empirical *L*-curve. This figure illustrates that the solution is smoother as the value of $\gamma_R$ increases. Figure 7 illustrates selection of $\gamma_R$ through the

<figure><img src="image_6.png" /><figcaption>Figure 6: Representation of the *L*-curve for a fixed SNR level, and identified parameters for different values of $\gamma_R$.</figcaption></figure>

*L*-curve and associated identified parameter, for two different noise levels. As expected, the magnitude of regularization is increased when the noise level is increased. We can also compare the identified parameters of figure 5 and figure 7 to conclude that both the classical and Bayesian approach yield similar results for parameter identification, and that the main difference between the two approaches lies in the statistical interpretation of covariance parameters, and the sampling capabilities of the Bayesian approach.

<figure><img src="image_7.png" /><figcaption>Figure 7: Identified classical regularization parameters from two *L*-curves associated with two different SNR levels, and Young modulus field with optimal values for $\gamma_R$.</figcaption></figure>

From the identified hyperparameters, Matérn random fields with similar covariance properties can then be sampled. We illustrate this in figure 8, where new Matérn fields are sampled from the prior distribution with tuned covariance. The samples were generated using the Python library GSTools [Müller et al., 2022, Müller and Schüler, 2023]. These samples illustrate the impact of measurement noise on the spatial statistics of the resulting generative model.
<figure><img src="image_7.png" /><figcaption>Figure 8: Identified stiffness parameter field $p^{\text{opt}} = \ln E^{\text{opt}}$ and random samples from the prior distribution with tuned Matérn covariance through the variogram regression approach, for increasing SNR levels.</figcaption></figure>

A more interesting result is presented in figure 9, where a tensile test was simulated for the same specimen, using the stiffness tensor distribution resulting from the sampled Matérn fields with tuned hyperparameters using the variogram regression approach. For each simulation, the boundary conditions presented at the beginning of section 9.1 were used. The resulting simulated displacement and strain fields illustrate the expressiveness of the generative model obtained when combining the tuned covariance model with the selected stiffness model: our method was able to learn the spatial statistics of kinematic fields that cannot simply be modelled as Gaussian random fields.

<figure><img src="image_8.png" /><figcaption>Figure 9: Identified fields and generated fields from the tuned hyperparameters and stiffness model: log Young modulus (a), tensile displacement (b), tensile strain (c). The identified fields were obtained using the variogram regression approach with a SNR of 55 dB.</figcaption></figure>

**Note.** In figure 10, even though the identified parameter field and hyperparameters appear to stabilize with increasing SNR (for two of the three approaches), there is no proper convergence due to the ill-posedness nature of the inverse elasticity problem. More precisely, for any fixed set of hyperparameters $\mu$, when $\text{SNR}^{\text{dB}} \to \infty$, then $\sigma_F^2 \to 0$ and $\|\Sigma_u^{-1}\| \to \infty$ for any norm, and consequently the loss function (34) is not well-defined. For very high SNR levels, the prior term of the loss becomes negligible and $J(p) \approx J_u(p) + J_F(p)$. This is essentially equivalent to a classical noise free inverse problem without regularization, which is numerically unstable due to the non-invertible nature of the forward operator $\mathcal{A}$. A Bayesian way to modify the modelling to ensure convergence with respect to SNR would be to model the ill-posedness of operator $\mathcal{A}$ as artificial operator noise $\zeta_A \sim \mathcal{N}(0, \sigma_A^2)$ and update the likelihood accordingly:

$$D(y^{\text{obs}}|p) \propto \exp \left[ -\frac{1}{2} (y^{\text{obs}} - y)^T (\Sigma_y + \sigma_A^2 I)^{-1} (y^{\text{obs}} - y) \right]. \quad (97)$$
<figure><img src="image_8.png" /><figcaption>Figure 8: Identified stiffness parameter field $p^{\text{opt}} = \ln E^{\text{opt}}$ and random samples from the prior distribution with tuned Matérn covariance through the variogram regression approach, for increasing SNR levels.</figcaption></figure>

A more interesting result is presented in figure 9, where a tensile test was simulated for the same specimen, using the stiffness tensor distribution resulting from the sampled Matérn fields with tuned hyperparameters using the variogram regression approach. For each simulation, the boundary conditions presented at the beginning of section 9.1 were used. The resulting simulated displacement and strain fields illustrate the expressiveness of the generative model obtained when combining the tuned covariance model with the selected stiffness model: our method was able to learn the spatial statistics of kinematic fields that cannot simply be modelled as Gaussian random fields.

<figure><img src="image_9.png" /><figcaption>Figure 9: Identified fields and generated fields from the tuned hyperparameters and stiffness model: log Young modulus (a), tensile displacement (b), tensile strain (c). The identified fields were obtained using the variogram regression approach with a SNR of 55 dB.</figcaption></figure>

**Note.** In figure 10, even though the identified parameter field and hyperparameters appear to stabilize with increasing SNR (for two of the three approaches), there is no proper convergence due to the ill-posedness nature of the inverse elasticity problem. More precisely, for any fixed set of hyperparameters $\mu$, when $\text{SNR}^{\text{dB}} \to \infty$, then $\sigma_F^2 \to 0$ and $\|\Sigma_u^{-1}\| \to \infty$ for any norm, and consequently the loss function (34) is not well-defined. For very high SNR levels, the prior term of the loss becomes negligible and $J(\mathbf{p}) \approx J_u(\mathbf{p})+J_F(\mathbf{p})$. This is essentially equivalent to a classical noise free inverse problem without regularization, which is numerically unstable due to the non-invertible nature of the forward operator $\mathcal{A}$. A Bayesian way to modify the modelling to ensure convergence with respect to SNR would be to model the ill-posedness of operator $\mathcal{A}$ as artificial operator noise $\zeta_{\mathcal{A}} \sim \mathcal{N}(0, \sigma_{\mathcal{A}}^2)$ and update the likelihood accordingly:

$$D(\mathbf{y}^{\text{obs}}|\mathbf{p}) \propto \exp \left[ -\frac{1}{2}(\mathbf{y}^{\text{obs}} - \mathbf{y})^T (\mathbf{\Sigma}_y + \sigma_{\mathcal{A}}^2 \mathbf{I})^{-1} (\mathbf{y}^{\text{obs}} - \mathbf{y}) \right]. \quad (97)$$
The added term acts as a classical regularization, and is equivalent to artificially increasing the prior variance and thus the marginal variance hyperparameter no longer has a statistical interpretation. The noise free case can also be treated following the classical approach developed in section 7 with $\gamma_F$ set to an arbitrary value. However, this asymptotic noise free case is not very important for actual DIC parameter identification, as it can only occur in numerical experiments.

## 9.2 Identification of a Spatial Orientation from Volume Observations

In this second example, we showcase our method on volume displacement observations, modelling a DVC experiment, in a case where the results for the hybrid-MAP and variogram fitting approaches differ. We consider the thick dog-bone specimen represented in figure 11, in which the observation zone is a volume subdomain $\Omega^{\text{obs}} \subset \mathbb{R}^3$. The considered stiffness model represents a transversely isotropic material, where the symmetry plane is orthogonal to the loading direction represented by the unit vector $\underline{n}_F$. The local orientation is defined in spherical coordinates $(\theta, \varphi)$ and represents the direction orthogonal to the symmetry plane in which the material is assumed to be stiffer. The local orientation is assumed to be spatially varying, which can be seen as a simplified model for waviness of fibers in a unidirectional composite. The transversely isotropic stiffness tensor is denoted by $\bar{C}_{\text{ref}}(E_1, E_2, G_{12}, \nu_{12}, \nu_{23})$ where $E_1, E_2, G_{12}, \nu_{12}, \nu_{23}$ are known constants. The stiffness of the wavy material is defined as

$$ \bar{C}(\theta, \varphi) = \underline{\underline{R}}(\theta, \varphi) : \underline{\underline{R}}(\theta, \varphi) : \bar{C}_{\text{ref}}(E_1, E_2, G_{12}, \nu_{12}, \nu_{23}) : \underline{\underline{R}}(\theta, \varphi) : \underline{\underline{R}}(\theta, \varphi), \quad (98) $$

where $\underline{\underline{R}}(\theta, \varphi)$ is a rotation matrix. For the generated data, we give the following closed form for $\theta^{\text{true}}$ and $\varphi^{\text{true}}$:

$$ \theta^{\text{true}}(x, y, z) = a \left( \cos(x/b) \cos(y/b) \sin(z/b) - 1/2 \right), \quad (99) $$

$$ \varphi^{\text{true}}(x, y, z) = a \left( \sin(x/b) \sin(y/b) \cos(y/b) - 1/2 \right), \quad (100) $$

where $a, b > 0$ are fixed values. As with the previous example, the noisy displacement is observed in the volume subdomain $\Omega^{\text{obs}}$, and the unknown discrete stiffness parameters $(\mathbf{p}_0, \mathbf{p}_1) = (\theta, \varphi)$ are modelled by a stationary Gaussian random field with Matérn covariance. The unknown hyperparameters are assumed to be identical for both stiffness parameters $(\mu_0 = \mu_1)$. The generated parameters $\phi_{\text{true}}$ and $\theta_{\text{true}}$, and MAP estimates $\phi_{\text{MAP}}$ and $\theta_{\text{MAP}}$ are represented in figure 12, where the variogram regression approach was used for hyperparameter identification. The results show that both MAP estimates are close to the true parameters, even though the parameters were not initially sampled from the prior distribution. What's interesting is that, given the analytical expressions (99)-(100), $\phi_{\text{true}}$ and $\theta_{\text{true}}$ are of class $C^\infty$ and thus cannot be modelled accurately with a rough Matérn covariance model (with $\nu = 1/2$). As discussed in section 6.2, the hybrid-MAP approach will overestimate the correlation length by an order of magnitude. This is represented in figure 13, where samples from the prior distribution are represented. In the first case, the hyperparameters are identified with the variogram regression approach, and in the other case the hybrid-MAP approach is used. A high SNR (60 dB) is used in both cases. The results show that the correlation length of the new samples is highly overestimated when the hybrid-MAP approach is used, whereas the variogram regression approach leads to the right correlation length. In this case, the variogram regression approach is thus better suited for hyperparameter identification than the hybrid-MAP approach.

## 9.3 Discussion

The numerical examples presented in sections 9.1 and 9.2 show that our method is able to estimate stiffness parameters from noisy surface or volume displacement observations. One of the main advantages of the proposed Bayesian formulation is that it accounts for measurement noise explicitly. The example shown in section 9.1 shows that the covariance structure of the identified stiffness parameters is highly dependent on the measurement noise level. More precisely, in a low SNR setting, the variance of the parameter of interest tends to be underestimated, and its correlation length overestimated, which is a common result in the Gaussian process regression literature (see, e.g. section 5.4.1 of Rasmussen [2004]). Unlike the classical regularization approach, the use of Whittle-Matérn priors models the variance and correlation length of the parameters of interest explicitly, which makes them
The added term acts as a classical regularization, and is equivalent to artificially increasing the prior variance and thus the marginal variance hyperparameter no longer has a statistical interpretation. The noise free case can also be treated following the classical approach developed in section 7 with $\gamma_F$ set to an arbitrary value. However, this asymptotic noise free case is not very important for actual DIC parameter identification, as it can only occur in numerical experiments.

## 9.2 Identification of a Spatial Orientation from Volume Observations

In this second example, we showcase our method on volume displacement observations, modelling a DVC experiment, in a case where the results for the hybrid-MAP and variogram fitting approaches differ. We consider the thick dog-bone specimen represented in figure 11, in which the observation zone is a volume subdomain $\Omega^{\text{obs}} \subset \mathbb{R}^3$. The considered stiffness model represents a transversely isotropic material, where the symmetry plane is orthogonal to the loading direction represented by the unit vector $\underline{n}_F$. The local orientation is defined in spherical coordinates $(\theta, \varphi)$ and represents the direction orthogonal to the symmetry plane in which the material is assumed to be stiffer. The local orientation is assumed to be spatially varying, which can be seen as a simplified model for waviness of fibers in a unidirectional composite. The transversely isotropic stiffness tensor is denoted by $\bar{C}_{\text{ref}}(E_1, E_2, G_{12}, \nu_{12}, \nu_{23})$ where $E_1, E_2, G_{12}, \nu_{12}, \nu_{23}$ are known constants. The stiffness of the wavy material is defined as

$$ \bar{C}(\theta, \varphi) = \underline{\underline{R}}(\theta, \varphi) : \underline{\underline{R}}(\theta, \varphi) : \bar{C}_{\text{ref}}(E_1, E_2, G_{12}, \nu_{12}, \nu_{23}) : \underline{\underline{R}}(\theta, \varphi) : \underline{\underline{R}}(\theta, \varphi), \quad (98) $$

where $\underline{\underline{R}}(\theta, \varphi)$ is a rotation matrix. For the generated data, we give the following closed form for $\theta^{\text{true}}$ and $\varphi^{\text{true}}$:

$$ \theta^{\text{true}}(x, y, z) = a \left( \cos(x/b) \cos(y/b) \sin(z/b) - 1/2 \right), \quad (99) $$

$$ \varphi^{\text{true}}(x, y, z) = a \left( \sin(x/b) \sin(y/b) \cos(y/b) - 1/2 \right), \quad (100) $$

where $a, b > 0$ are fixed values. As with the previous example, the noisy displacement is observed in the volume subdomain $\Omega^{\text{obs}}$, and the unknown discrete stiffness parameters $(\mathbf{p}_0, \mathbf{p}_1) = (\theta, \varphi)$ are modelled by a stationary Gaussian random field with Matérn covariance. The unknown hyperparameters are assumed to be identical for both stiffness parameters $(\mu_0 = \mu_1)$. The generated parameters $\phi_{\text{true}}$ and $\theta_{\text{true}}$, and MAP estimates $\phi_{\text{MAP}}$ and $\theta_{\text{MAP}}$ are represented in figure 12, where the variogram regression approach was used for hyperparameter identification. The results show that both MAP estimates are close to the true parameters, even though the parameters were not initially sampled from the prior distribution. What's interesting is that, given the analytical expressions (99)-(100), $\phi_{\text{true}}$ and $\theta_{\text{true}}$ are of class $C^\infty$ and thus cannot be modelled accurately with a rough Matérn covariance model (with $\nu = 1/2$). As discussed in section 6.2, the hybrid-MAP approach will overestimate the correlation length by an order of magnitude. This is represented in figure 13, where samples from the prior distribution are represented. In the first case, the hyperparameters are identified with the variogram regression approach, and in the other case the hybrid-MAP approach is used. A high SNR (60 dB) is used in both cases. The results show that the correlation length of the new samples is highly overestimated when the hybrid-MAP approach is used, whereas the variogram regression approach leads to the right correlation length. In this case, the variogram regression approach is thus better suited for hyperparameter identification than the hybrid-MAP approach.

## 9.3 Discussion

The numerical examples presented in sections 9.1 and 9.2 show that our method is able to estimate stiffness parameters from noisy surface or volume displacement observations. One of the main advantages of the proposed Bayesian formulation is that it accounts for measurement noise explicitly. The example shown in section 9.1 shows that the covariance structure of the identified stiffness parameters is highly dependent on the measurement noise level. More precisely, in a low SNR setting, the variance of the parameter of interest tends to be underestimated, and its correlation length overestimated, which is a common result in the Gaussian process regression literature (see, e.g. section 5.4.1 of Rasmussen [2004]). Unlike the classical regularization approach, the use of Whittle-Matérn priors models the variance and correlation length of the parameters of interest explicitly, which makes them
<figure><img src="image_9.png" /></figure>

(a)

<figure><img src="image_10.png" /></figure>

(b)

<figure><img src="image_11.png" /></figure>

(c)

<caption>Figure 10: Identified Matérn hyperparameters with respect to observation SNR, using the three PEB approaches: variogram regression (10a), hybrid-MAP estimator with a hyperprior on $\sigma^2$ (10b), ML estimator (10c).</caption>
<figure><img src="image_10.png" /><figcaption>Figure 10: Identified Matérn hyperparameters with respect to observation SNR, using the three PEB approaches: variogram regression (10a), hybrid-MAP estimator with a hyperprior on $\sigma^2$ (10b), ML estimator (10c).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11: Thick dog-bone geometry for example with volume observations</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Generated displacement (a), generated parameters $\phi^{\text{true}}$ and $\theta^{\text{true}}$ (b), and identified parameters $\phi^{\text{opt}}$ and $\theta^{\text{opt}}$ with a 60 dB SNR (c). Hyperparameter selection was performed with the variogram regression approach.</figcaption></figure>

interpretable with regard to microstructural features. As an ill-posed inverse problem, the Whittle-Matérn prior also acts as a regularizer, and we have shown the similarities with using a classical inverse problem formulation using a gradient penalty regularization. The main advantage of our approach is that formulating a Matérn random field prior on the unknown material parameter leads to a generative model after hyperparameter tuning: combining the identified correlation length $\rho$, the identified marginal variance $\sigma^2$ with a carefully chosen stiffness model $p \mapsto \bar{C}(p)$ enables modelling a wide range of random spatial elastic data $\bar{C}(x)$.

In this paper, we investigated various approaches for hyperparameters estimation, based on a Parametric Empirical Bayes approach. The various issues related to estimability of the Matérn covariance parameter pointed out in section 6.2 are due to several numerical approximations of the underlying hierarchical Bayesian model, namely iteratively estimating hyperparameters from the data, and the use of sparse approximation of the precision matrix with fixed smoothness parameters. A possible extension of this work would be to consider the full hierarchical model instead of the PEB approximation, and sample from the posterior distribution for $p$ and $\mu$ using MCMC sampling. Considering the recent papers that present numerical methods for a sparse approximation of the Matérn precision matrix for any $\nu > 0$ (see, e.g. [Roininen et al., 2018, Bolin et al., 2018]), MCMC sampling for all Matérn hyperparameters could be investigated, although the computational cost would be much higher than the presented framework. Under the PEB approximation, we've shown that the estimator for hyperparameters must be chosen carefully to avoid numerical instabilities. We've discussed the issues related to the use of a maximum likelihood estimator, and linked the numerical results to the existing literature. The conclusions of our investigations are that the variogram regression approach is the best choice for all cases. Using this approach, we were able to identify the parameter of interest and their covariance structure accurately in both numerical experiments.
<figure><img src="image_11.png" /><figcaption>Figure 11: Thick dog-bone geometry for example with volume observations</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Generated displacement (a), generated parameters $\phi^{\text{true}}$ and $\theta^{\text{true}}$ (b), and identified parameters $\phi^{\text{opt}}$ and $\theta^{\text{opt}}$ with a 60 dB SNR (c). Hyperparameter selection was performed with the variogram regression approach.</figcaption></figure>

interpretable with regard to microstructural features. As an ill-posed inverse problem, the Whittle-Matérn prior also acts as a regularizer, and we have shown the similarities with using a classical inverse problem formulation using a gradient penalty regularization. The main advantage of our approach is that formulating a Matérn random field prior on the unknown material parameter leads to a generative model after hyperparameter tuning: combining the identified correlation length $\rho$, the identified marginal variance $\sigma^2$ with a carefully chosen stiffness model $p \mapsto \bar{C}(p)$ enables modelling a wide range of random spatial elastic data $\bar{C}(x)$.

In this paper, we investigated various approaches for hyperparameters estimation, based on a Parametric Empirical Bayes approach. The various issues related to estimability of the Matérn covariance parameter pointed out in section 6.2 are due to several numerical approximations of the underlying hierarchical Bayesian model, namely iteratively estimating hyperparameters from the data, and the use of sparse approximation of the precision matrix with fixed smoothness parameters. A possible extension of this work would be to consider the full hierarchical model instead of the PEB approximation, and sample from the posterior distribution for $p$ and $\mu$ using MCMC sampling. Considering the recent papers that present numerical methods for a sparse approximation of the Matérn precision matrix for any $\nu > 0$ (see, e.g. [Roininen et al., 2018, Bolin et al., 2018]), MCMC sampling for all Matérn hyperparameters could be investigated, although the computational cost would be much higher than the presented framework. Under the PEB approximation, we've shown that the estimator for hyperparameters must be chosen carefully to avoid numerical instabilities. We've discussed the issues related to the use of a maximum likelihood estimator, and linked the numerical results to the existing literature. The conclusions of our investigations are that the variogram regression approach is the best choice for all cases. Using this approach, we were able to identify the parameter of interest and their covariance structure accurately in both numerical experiments.
<figure><img src="image_14.png" /><figcaption>Figure 13: Random samples from the prior distribution with tuned Matérn covariance, for the hybrid-MAP approach (a) and the variogram regression approach (b).</figcaption></figure>

In this work, we made the assumption that the covariance model structure was stationary. This may not always be relevant to model materials with non-stationary microstructural features, such as when edge effects are present. In this case, spatially varying hyperparameters could be considered. Some preliminary work is presented in [Roininen et al., 2019] where Whittle-Matérn priors with a spatially varying length scale is considered. The main challenge with this extension is that a correlation structure for the hyperparameters must be specified and tuned. Another limitation of our work is that we did not investigate the case where the stiffness model depends on multiple parameters, each having a different covariance structure. In this case, the number of hyperparameters to identify is increased, and the estimability of all hyperparameters should be investigated. As this paper was focused on identification of linear elastic properties, only a single displacement field observation was required to perform the identification. The Bayesian framework could be extended to time series observations in order to identify spatially varying parameters from a non-linear constitutive law. In order to achieve this, the loss function should account for the observed displacement for all time steps $u_t^{\text{obs}}$, $t = 0, ..., N_t$, and the forward and adjoint problem solvers should be modified accordingly. However, the Bayesian modelling would remain unchanged.

# 10 Conclusion

The literature on numerical methods for materials with spatially varying, and random elastic properties is scarce, perhaps due to the fact that such materials are used in very specific contexts only. In this work, a Bayesian extension to parameter identification procedures from noisy displacement observations was presented, in which the unknown parameter is modelled as a Matérn random field. Our approach naturally leads to a generative model in which the identified covariance parameters can be used to sample new realisations from the random field prior, and generate spatially varying elasticity tensor fields. Also, the identified covariance hyperparameters are informative about the spatial statistics of underlying microstructural features. The numerical aspects of the proposed framework were presented and discussed. As shown, extending existing procedures to a spatially varying context involves quite a lot of specific challenges, and spatial covariance properties must be accounted for. The approach presented in this paper leverages numerical methods from fields of research largely unrelated to structural mechanics, which we believe to be very relevant for modelling such materials. As explained in the introduction, the authors aim at applying the presented framework to DLF composites, in which spatial covariance of elastic properties is of major importance, and is going to be the main topic of upcoming works.

# Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
<figure><img src="image_13.png" /><figcaption>Figure 13: Random samples from the prior distribution with tuned Matérn covariance, for the hybrid-MAP approach (a) and the variogram regression approach (b).</figcaption></figure>

In this work, we made the assumption that the covariance model structure was stationary. This may not always be relevant to model materials with non-stationary microstructural features, such as when edge effects are present. In this case, spatially varying hyperparameters could be considered. Some preliminary work is presented in [Roininen et al., 2019] where Whittle-Matérn priors with a spatially varying length scale is considered. The main challenge with this extension is that a correlation structure for the hyperparameters must be specified and tuned. Another limitation of our work is that we did not investigate the case where the stiffness model depends on multiple parameters, each having a different covariance structure. In this case, the number of hyperparameters to identify is increased, and the estimability of all hyperparameters should be investigated. As this paper was focused on identification of linear elastic properties, only a single displacement field observation was required to perform the identification. The Bayesian framework could be extended to time series observations in order to identify spatially varying parameters from a non-linear constitutive law. In order to achieve this, the loss function should account for the observed displacement for all time steps $u_t^{\text{obs}}$, $t = 0, ..., N_t$, and the forward and adjoint problem solvers should be modified accordingly. However, the Bayesian modelling would remain unchanged.

# 10 Conclusion

The literature on numerical methods for materials with spatially varying, and random elastic properties is scarce, perhaps due to the fact that such materials are used in very specific contexts only. In this work, a Bayesian extension to parameter identification procedures from noisy displacement observations was presented, in which the unknown parameter is modelled as a Matérn random field. Our approach naturally leads to a generative model in which the identified covariance parameters can be used to sample new realisations from the random field prior, and generate spatially varying elasticity tensor fields. Also, the identified covariance hyperparameters are informative about the spatial statistics of underlying microstructural features. The numerical aspects of the proposed framework were presented and discussed. As shown, extending existing procedures to a spatially varying context involves quite a lot of specific challenges, and spatial covariance properties must be accounted for. The approach presented in this paper leverages numerical methods from fields of research largely unrelated to structural mechanics, which we believe to be very relevant for modelling such materials. As explained in the introduction, the authors aim at applying the presented framework to DLF composites, in which spatial covariance of elastic properties is of major importance, and is going to be the main topic of upcoming works.

# Declaration of Competing Interest

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

Hutchinson SA, a company of TotalEnergies, is acknowledged for the financial support of this study. This work was also supported by the French Ministry of Research (ANRT).

## A Extension to I-DIC

Adapting the I-DIC method to our framework is similar to the development for FEMU. Following [Leclerc et al., 2009, Réthoré, 2010, Neggers et al., 2017], we introduce the following notations:

* $f$ is the reference image,

* $g$ is the deformed image,

* $\tilde{g}$ is the back transformed image, defined for a given displacement field $\underline{u}$ as

$$ \tilde{g}(\underline{x}) = g(\underline{x} + \underline{u}(\underline{x})), $$

* $N$ is the number of pixels of all considered images.

In all this section, an image $I$ of shape $(w, h) \in \mathbb{N}^2$ is represented by a smooth function $I: [0; w] \times [0; h] \to [0; 1]$ where $I(\underline{x})$ is defined by the bicubic interpolation of neighbour pixel values. When written with a bold letter, an image $I$ denotes a $w \times h$ vector collecting pixel values.

The I-DIC loss is formulated in terms of a gray level misfit between the back transformed image and the reference image. In the Bayesian formulation, the gray level noise is assumed to be Gaussian and spatially uncorrelated. We note $\sigma_f$ the variance of the white Gaussian noise affecting the captured gray level, and as commonly done in, e.g. [Réthoré, 2010], to simplify notations, only the reference image is assumed to be noisy with a variance $2\sigma_f^2$. Let $\underline{u}^{\text{true}}$ denote the true displacement field that maps $f$ to $g$. Then

$$ \forall 1 \le i \le N, \quad g(\underline{x}_i + \underline{u}^{\text{true}}(\underline{x}_i)) - f(\underline{x}_i) \sim \mathcal{N}(0, 2\sigma_f^2). \qquad (101) $$

The I-DIC loss can be written as the following likelihood:

$$ J_I(\underline{u}(\underline{p})) = \frac{1}{2\sigma_f^2} \sum_{i=1}^{N} (g(\underline{x}_i + \underline{u}(\underline{p})(\underline{x}_i)) - f(\underline{x}_i))^2. \qquad (102) $$

Using (102) instead of (35) inside (34), the loss function for the MAP estimate using the I-DIC approach now writes

$$ J_{\text{IDIC}}(\underline{p}) = J_I(\underline{p}) + J_F(\underline{p}) + J_p(\underline{p}), \qquad (103) $$

where $J_F$ and $J_p$ have the same expressions as for the FEMU method. Before deriving the parameter sensitivities for the I-DIC loss, it is important to note that the displacement field $\underline{u}$ present in $J_I$ is not exactly the solution to the forward problem, but the simulated displacement mapped to image coordinates, as the simulated displacement $\underline{u}$ is defined on $\Omega$ and $\underline{u}$ is defined on $[0; w] \times [0; h]$. Hence, we introduce an invertible affine transformation $\underline{\phi}$ satisfying $\underline{\phi}(\Omega^{\text{obs}}) \subset [0; w] \times [0; h]$ that maps the computational domain to the reference image, and we pose $\underline{\psi} = \underline{\phi}^{-1}$. Also, after discretization, $\underline{u}$ is a displacement image that is not represented using finite elements shape functions. To avoid any confusion, we hereafter use the following notations:

* $\underline{u}: \Omega \to \mathbb{R}^3$ is the discrete forward problem solution, expressed in a finite element basis $(\underline{\varphi}_0, \dots, \underline{\varphi}_{N_V})$ as

$$ \underline{u} = \sum_{i=0}^{N_V} u_i \underline{\varphi}_i, $$

* $\tilde{u}: [0; w] \times [0; h] \to \mathbb{R}^3$ is the forward problem solution in image coordinates, defined as

$$ \tilde{u}(\underline{x}) = d\underline{\psi} \circ \underline{u} \circ \underline{\psi}, $$
# Acknowledgments

Hutchinson SA, a company of TotalEnergies, is acknowledged for the financial support of this study. This work was also supported by the French Ministry of Research (ANRT).

## A Extension to I-DIC

Adapting the I-DIC method to our framework is similar to the development for FEMU. Following [Leclerc et al., 2009, Réthoré, 2010, Neggers et al., 2017], we introduce the following notations:

* $f$ is the reference image,

* $g$ is the deformed image,

* $\tilde{g}$ is the back transformed image, defined for a given displacement field $\underline{u}$ as

$$ \tilde{g}(\underline{x}) = g(\underline{x} + \underline{u}(\underline{x})), $$

* $N$ is the number of pixels of all considered images.

In all this section, an image $I$ of shape $(w, h) \in \mathbb{N}^2$ is represented by a smooth function $I: [0; w] \times [0; h] \to [0; 1]$ where $I(\underline{x})$ is defined by the bicubic interpolation of neighbour pixel values. When written with a bold letter, an image $I$ denotes a $w \times h$ vector collecting pixel values.

The I-DIC loss is formulated in terms of a gray level misfit between the back transformed image and the reference image. In the Bayesian formulation, the gray level noise is assumed to be Gaussian and spatially uncorrelated. We note $\sigma_f$ the variance of the white Gaussian noise affecting the captured gray level, and as commonly done in, e.g. [Réthoré, 2010], to simplify notations, only the reference image is assumed to be noisy with a variance $2\sigma_f^2$. Let $\underline{u}^{\text{true}}$ denote the true displacement field that maps $f$ to $g$. Then

$$ \forall 1 \le i \le N, \quad g(\underline{x}_i + \underline{u}^{\text{true}}(\underline{x}_i)) - f(\underline{x}_i) \sim \mathcal{N}(0, 2\sigma_f^2). \qquad (101) $$

The I-DIC loss can be written as the following likelihood:

$$ J_I(\underline{u}(\underline{p})) = \frac{1}{2\sigma_f^2} \sum_{i=1}^{N} (g(\underline{x}_i + \underline{u}(\underline{p})(\underline{x}_i)) - f(\underline{x}_i))^2. \qquad (102) $$

Using (102) instead of (35) inside (34), the loss function for the MAP estimate using the I-DIC approach now writes

$$ J_{\text{IDIC}}(\underline{p}) = J_I(\underline{p}) + J_F(\underline{p}) + J_p(\underline{p}), \qquad (103) $$

where $J_F$ and $J_p$ have the same expressions as for the FEMU method. Before deriving the parameter sensitivities for the I-DIC loss, it is important to note that the displacement field $\underline{u}$ present in $J_I$ is not exactly the solution to the forward problem, but the simulated displacement mapped to image coordinates, as the simulated displacement $\underline{u}$ is defined on $\Omega$ and $\underline{u}$ is defined on $[0; w] \times [0; h]$. Hence, we introduce an invertible affine transformation $\underline{\phi}$ satisfying $\underline{\phi}(\Omega^{\text{obs}}) \subset [0; w] \times [0; h]$ that maps the computational domain to the reference image, and we pose $\underline{\psi} = \underline{\phi}^{-1}$. Also, after discretization, $\underline{u}$ is a displacement image that is not represented using finite elements shape functions. To avoid any confusion, we hereafter use the following notations:

* $\underline{u}: \Omega \to \mathbb{R}^3$ is the discrete forward problem solution, expressed in a finite element basis $(\underline{\varphi}_0, \dots, \underline{\varphi}_{N_V})$ as

$$ \underline{u} = \sum_{i=0}^{N_V} u_i \underline{\varphi}_i, $$

* $\tilde{u}: [0; w] \times [0; h] \to \mathbb{R}^3$ is the forward problem solution in image coordinates, defined as

$$ \tilde{u}(\underline{x}) = d\underline{\psi} \circ \underline{u} \circ \underline{\psi}, $$
• $\underline{U}: [0; w] \times [0; h] \to \mathbb{R}^2$ is a discrete displacement image, that is linked to $\tilde{u}$ through some image interpolation operator $\underline{\Pi}$ such that $\underline{U}(\underline{x}_i) = (\underline{\Pi} \circ \tilde{u})(\underline{x}_i)$.

Also, in regard to the previously introduced discrete notations, $\underline{u} = (u_0, ..., u_n)$ denotes the coordinates of $\underline{u}$ in basis $\mathcal{B} = (\underline{\varphi}_0, ..., \underline{\varphi}_n)$, $\tilde{u} = (\tilde{u}_0, ..., \tilde{u}_n)$ denotes the coordinates of $\underline{u}$ in basis $\tilde{\mathcal{B}} = (\tilde{\underline{\varphi}}_0, ..., \tilde{\underline{\varphi}}_n)$ where

$$\forall 0 \le i \le N_V, \quad \tilde{\underline{\varphi}}_i = \underline{d\psi} \circ \underline{\varphi}_i \circ \underline{\psi},$$

and $\underline{U}$ is a $2 \times w \times h$ tensor such that $\underline{U}_{ikl} = \underline{U}(k, l)_i$. With these notations, equation (102) now writes after discretization

$$J_I = \frac{1}{2\sigma_f^2} \sum_{i=1}^{N} \left( g(\underline{x}_i + (\underline{\Pi} \circ \underline{d\psi} \circ \underline{u} \circ \underline{\psi}(\underline{x}_i))) - f(\underline{x}_i) \right)^2. \quad (104)$$

The adjoint variational problem is derived in a same manner as for FEMU, the main difference is the expression of the source term:

$$\frac{\mathrm{d}J_{\mathrm{IDIC}}}{\mathrm{d}\underline{u}} = \frac{\mathrm{d}J_I}{\mathrm{d}\underline{u}} + \frac{\mathrm{d}J_F}{\mathrm{d}\underline{u}}. \quad (105)$$

The derivative of $J_F$ remains unchanged, and the derivative of $J_I$ is expanded using the chain rule

$$\frac{\mathrm{d}J_I}{\mathrm{d}\underline{u}} = \frac{1}{\sigma_f^2}(\tilde{g} - f) \frac{\mathrm{d}\tilde{g}}{\mathrm{d}\underline{U}} \frac{\mathrm{d}\tilde{u}}{\mathrm{d}\underline{u}}. \quad (106)$$

In this last expression, each term of the chain rule expansion is explained below:

* $\frac{\mathrm{d}\tilde{u}}{\mathrm{d}\underline{u}}$ is a block-diagonal sparse matrix of size $(3(N_V + 1))^2$ that is the discrete representation of $\underline{d\psi}$ in basis $\mathcal{B}$. If we note $\underline{\psi}(\underline{x}) = \underline{\underline{A}} \cdot \underline{x} + \underline{b}$, then

$$\frac{\mathrm{d}\tilde{u}}{\mathrm{d}\underline{u}} = \begin{pmatrix} \underline{\underline{A}} & 0 \\ \vdots & \vdots \\ 0 & \underline{\underline{A}} \end{pmatrix}, \quad (107)$$

* $\frac{\mathrm{d}U}{\mathrm{d}\tilde{u}}$ is a sparse matrix of size $(2N) \times (3N_V)$ that stores the finite element basis values at pixel coordinates. In terms of interpolation, if $\underline{U}$ is seen as a vector valued image, then that image can be interpreted as a structured mesh made of quadrilateral cells, on which can be constructed a $P^1$ finite element basis $\mathcal{B}_{\mathrm{pix}}$. That being said, $\frac{\mathrm{d}U}{\mathrm{d}\tilde{u}}$ is simply the discrete interpolator from $\tilde{\mathcal{B}}$ to $\mathcal{B}_{\mathrm{pix}}$. The matrix entries are given by

$$\frac{\mathrm{d}U}{\mathrm{d}\tilde{u}} = \left( (\underline{\varphi}_k(\underline{x}_l))_i \right)_{\substack{i=0,1 \\ 0 \le k \le N_V \\ 0 \le l < N}} \quad (108)$$

where $\underline{x}_l$ denotes pixel coordinates,

* $\frac{\mathrm{d}\tilde{g}}{\mathrm{d}\tilde{u}}$ is a block-diagonal sparse matrix of size $N \times 2N$ whose entries are made of the back-transformed image gradient $\nabla\tilde{g}$. There are multiple approaches to compute the discrete gradient of an image, but in order to remain consistent with the bicubic gray level interpolation, we use a centred finite difference scheme:

$$(\nabla_x \tilde{g})(\underline{x}) = \lim_{\epsilon \to 0} \frac{\tilde{g}(\underline{x} + \epsilon \underline{e}_x) - \tilde{g}(\underline{x} - \epsilon \underline{e}_x)}{2\epsilon}, \quad (109)$$

$$(\nabla_y \tilde{g})(\underline{x}) = \lim_{\epsilon \to 0} \frac{\tilde{g}(\underline{x} + \epsilon \underline{e}_y) - \tilde{g}(\underline{x} - \epsilon \underline{e}_y)}{2\epsilon}, \quad (110)$$

where $\tilde{g}(\underline{x} \pm \epsilon \underline{e}_x)$ is defined through bicubic interpolation, $\underline{e}_x = (1, 0)$ and $\underline{e}_y = (0, 1)$. From there,

$$\frac{\mathrm{d}\tilde{g}}{\mathrm{d}\tilde{u}} = \begin{pmatrix} \nabla\tilde{g}(\underline{x}_0) & 0 \\ \vdots & \vdots \\ 0 & \nabla\tilde{g}(\underline{x}_N) \end{pmatrix} \quad (111)$$

where

$$\nabla\tilde{g}(\underline{x}_l) = ((\nabla_x\tilde{g})(\underline{x}_l), (\nabla_y\tilde{g})(\underline{x}_l)).$$
• $\underline{U}: [0; w] \times [0; h] \to \mathbb{R}^2$ is a discrete displacement image, that is linked to $\tilde{u}$ through some image interpolation operator $\underline{\Pi}$ such that $\underline{U}(\underline{x}_i) = (\underline{\Pi} \circ \tilde{u})(\underline{x}_i)$.

Also, in regard to the previously introduced discrete notations, $\underline{u} = (u_0, ..., u_n)$ denotes the co-ordinates of $\underline{u}$ in basis $\mathcal{B} = (\underline{\varphi}_0, ..., \underline{\varphi}_n)$, $\tilde{\underline{u}} = (\tilde{u}_0, ..., \tilde{u}_n)$ denotes the coordinates of $\underline{u}$ in basis $\tilde{\mathcal{B}} = (\tilde{\underline{\varphi}}_0, ..., \tilde{\underline{\varphi}}_n)$ where

$$ \forall 0 \le i \le N_V, \quad \tilde{\underline{\varphi}}_i = \underline{d\psi} \circ \underline{\varphi}_i \circ \underline{\psi}, $$

and $\underline{U}$ is a $2 \times w \times h$ tensor such that $\underline{U}_{ikl} = \underline{U}(k, l)_i$. With these notations, equation (102) now writes after discretization

$$ J_I = \frac{1}{2\sigma_f^2} \sum_{i=1}^{N} \left( g(\underline{x}_i + (\underline{\Pi} \circ \underline{d\psi} \circ \underline{u} \circ \underline{\psi}(\underline{x}_i))) - f(\underline{x}_i) \right)^2. \quad (104) $$

The adjoint variational problem is derived in a same manner as for FEMU, the main difference is the expression of the source term:

$$ \frac{dJ_{\text{IDIC}}}{d\underline{u}} = \frac{dJ_I}{d\underline{u}} + \frac{dJ_F}{d\underline{u}}. \quad (105) $$

The derivative of $J_F$ remains unchanged, and the derivative of $J_I$ is expanded using the chain rule

$$ \frac{dJ_I}{d\underline{u}} = \frac{1}{\sigma_f^2}(\tilde{\underline{g}} - \underline{f}) \frac{d\tilde{\underline{g}}}{d\underline{U}} \frac{d\underline{U}}{d\tilde{\underline{u}}} \frac{d\tilde{\underline{u}}}{d\underline{u}}. \quad (106) $$

In this last expression, each term of the chain rule expansion is explained below:

* $\frac{d\tilde{\underline{u}}}{d\underline{u}}$ is a block-diagonal sparse matrix of size $(3(N_V + 1))^2$ that is the discrete representation of $\underline{d\psi}$ in basis $\mathcal{B}$. If we note $\underline{\psi}(\underline{x}) = \underline{A} \cdot \underline{x} + \underline{b}$, then

$$ \frac{d\tilde{\underline{u}}}{d\underline{u}} = \begin{pmatrix} \underline{A} & 0 \\ \vdots & \vdots \\ 0 & \underline{A} \end{pmatrix}, \quad (107) $$

* $\frac{d\underline{U}}{d\tilde{\underline{u}}}$ is a sparse matrix of size $(2N) \times (3N_V)$ that stores the finite element basis values at pixel coordinates. In terms of interpolation, if $\underline{U}$ is seen as a vector valued image, then that image can be interpreted as a structured mesh made of quadrilateral cells, on which can be constructed a $P^1$ finite element basis $\mathcal{B}_{\text{pix}}$. That being said, $\frac{d\underline{U}}{d\tilde{\underline{u}}}$ is simply the discrete interpolator from $\tilde{\mathcal{B}}$ to $\mathcal{B}_{\text{pix}}$. The matrix entries are given by

$$ \frac{d\underline{U}}{d\tilde{\underline{u}}} = \left( (\underline{\varphi}_k(\underline{x}_l))_i \right)_{\substack{i=0,1 \\ 0 \le k \le N_V \\ 0 \le l < N}} \quad (108) $$

where $\underline{x}_l$ denotes pixel coordinates,

* $\frac{d\tilde{\underline{g}}}{d\underline{u}}$ is a block-diagonal sparse matrix of size $N \times 2N$ whose entries are made of the back-transformed image gradient $\nabla \tilde{\underline{g}}$. There are multiple approaches to compute the discrete gradient of an image, but in order to remain consistent with the bicubic gray level interpolation, we use a centred finite difference scheme:

$$ (\nabla_x \tilde{\underline{g}})(\underline{x}) = \lim_{\epsilon \to 0} \frac{\tilde{\underline{g}}(\underline{x} + \epsilon \underline{e}_x) - \tilde{\underline{g}}(\underline{x} - \epsilon \underline{e}_x)}{2\epsilon}, \quad (109) $$

$$ (\nabla_y \tilde{\underline{g}})(\underline{x}) = \lim_{\epsilon \to 0} \frac{\tilde{\underline{g}}(\underline{x} + \epsilon \underline{e}_y) - \tilde{\underline{g}}(\underline{x} - \epsilon \underline{e}_y)}{2\epsilon}, \quad (110) $$

where $\tilde{\underline{g}}(\underline{x} \pm \epsilon \underline{e}_x)$ is defined through bicubic interpolation, $\underline{e}_x = (1, 0)$ and $\underline{e}_y = (0, 1)$. From there,

$$ \frac{d\tilde{\underline{g}}}{d\tilde{\underline{u}}} = \begin{pmatrix} \nabla \tilde{\underline{g}}(\underline{x}_0) & 0 \\ \vdots & \vdots \\ 0 & \nabla \tilde{\underline{g}}(\underline{x}_N) \end{pmatrix} \quad (111) $$

where

$$ \nabla \tilde{\underline{g}}(\underline{x}_l) = ((\nabla_x \tilde{\underline{g}})(\underline{x}_l), (\nabla_y \tilde{\underline{g}})(\underline{x}_l)). $$