%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2025-11-28 12:25:09
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Jailin_20250751.bib}

\usepackage{siunitx}
\usepackage{subcaption}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\newcommand{\ie}{i.e.}
\newcommand{\eg}{e.g.}

\let\bm\boldsymbol
\let\mathbf\boldsymbol

\newcommand{\BC}{\mathrm{BC}}
\newcommand{\GT}{\mathrm{GT}}
\newcommand{\Imag}{\mathrm{Im}}
\newcommand{\ICNN}{\mathrm{ICNN}}
\newcommand{\PANN}{\mathrm{PANN}}
\newcommand{\valid}{\mathrm{valid}}
\newcommand{\intern}{\mathrm{int}}

\DeclareMathOperator{\bias}{bias}
\DeclareMathOperator{\corr}{corr}

\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\abs}{abs}
\DeclareMathOperator{\cof}{cof}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\trace}{trace}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}
\DeclarePairedDelimiter{\angles}{\langle}{\rangle}

%%%---------------------------------------------------------------------------

%%% FIGURES

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

%%%---------------------------------------------------------------------------


\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}

\title{Noise-bias compensation for the unsupervised learning of constitutive laws}
\alttitle{Compensation du biais induit par le bruit pour l’apprentissage non supervisé des lois de comportement}

\author{\firstname{Clément} \lastname{Jailin}}
\address{Universit\'e Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Laboratoire de M\'ecanique Paris-Saclay (LMPS), 91190 Gif-sur-Yvette, France}
\email{clement.jailin@ens-paris-saclay.fr}

\author{\firstname{Stéphane} \lastname{Roux}}
\address[1]{Universit\'e Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Laboratoire de M\'ecanique Paris-Saclay (LMPS), 91190 Gif-sur-Yvette, France}

\author{\firstname{Antoine} \lastname{Benady}}
\address[1]{Universit\'e Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Laboratoire de M\'ecanique Paris-Saclay (LMPS), 91190 Gif-sur-Yvette, France}
\address{Department of Mechanical and Process Engineering, ETH Zürich, Zürich, Switzerland}

\author{\firstname{Emmanuel} \lastname{Baranger}}
\address[1]{Universit\'e Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Laboratoire de M\'ecanique Paris-Saclay (LMPS), 91190 Gif-sur-Yvette, France}

\thanks{This project was made possible by the Agence Nationale de la Recherche (ANR) grant No.~R-22-CPJ2-0046-01}
\CDRGrant[ANR]{R-22-CPJ2-0046-01}

\thanks{Part of this project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No.~101002857)}
\CDRGrant[ERC]{101002857}

\keywords{\kwd{Artificial intelligence} \kwd{PANN} \kwd{EUCLID} \kwd{experimental noise} \kwd{digital image correlation} \kwd{constitutive modeling}}
\altkeywords{\kwd{Intelligence artificielle} \kwd{réseaux de neurones augmentés par la physique (PANN)} \kwd{EUCLID} \kwd{bruit expérimental} \kwd{corrélation d'images numériques} \kwd{lois de comportement}}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\begin{abstract}
The Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) framework allows the non-supervised learning of constitutive laws from full-field displacement data and global reaction forces. Nonetheless, its accuracy is adversely affected by measurement noise, resulting in biased material parameter identification due to uniform nodal weighting and mesh dependencies. To mitigate these issues, covariance-based weighting and systematic bias compensation strategies are proposed, which account for measurement uncertainties and counteract noise-induced errors. Additionally, Gaussian smoothing is introduced as a low-cost alternative to reduce noise impact by averaging nodal force residuals. These methods were evaluated using both a simplified toy model and a complex numerical test case with realistic noise levels. Results demonstrate that the proposed compensation strategies significantly enhance EUCLID’s robustness and accuracy, achieving up to 93\% improvement in validation metrics under high-noise conditions. Furthermore, mesh dependency issues are addressed, enabling mesh-independent learning. These advancements substantially improve the reliability of constitutive law identification in noisy experimental environments.
\end{abstract}

\begin{altabstract}
Le cadre EUCLID (Efficient Unsupervised Constitutive Law Identification and Discovery) permet l’apprentissage non supervisé de lois de comportement à partir de champs de déplacements et de forces globales. Toutefois, son exactitude est fortement dégradée en présence de bruit de mesure, ce qui induit une identification biaisée des paramètres matériaux en raison d’un poids nodal uniforme et d’une dépendance au maillage. Pour remédier à ces limitations, nous proposons une pondération fondée sur les covariances ainsi que des stratégies de compensation du biais, capables de tenir compte des incertitudes de mesure et de corriger les erreurs induites par le bruit. En complément, un lissage gaussien est introduit comme solution peu coûteuse pour réduire l’impact du bruit en moyennant les résidus nodaux de forces. Ces approches sont évaluées à la fois sur un modèle simplifié et sur un cas test numérique complexe intégrant des niveaux de bruit réalistes. Les résultats montrent que les stratégies de compensation proposées améliorent significativement la robustesse et la précision d’EUCLID, avec jusqu’à 93\% d’amélioration des métriques de validation en conditions fortement bruitées. Par ailleurs, les problèmes de dépendance au maillage sont corrigés, permettant un apprentissage indépendant de la discrétisation. Ces avancées renforcent considérablement la fiabilité de l’identification de lois de comportement dans des environnements expérimentaux bruités.
\end{altabstract}

\dateposted{2026-01-05}
\begin{document}
%\input{CR-pagedemetas}
%\end{document}
\maketitle

\section{Introduction}

Artificial Intelligence (AI) holds the potential to transform mechanical science, particularly when modeling constitutive laws~\cite{fuhg2024review}. By leveraging its flexibility and capacity to identify complex patterns, AI has enabled the accurate representation of intricate material behaviors that has the capacity to represent traditional physical methods. This shift in paradigm has accelerated advancements in material science and empowered industries to achieve more ambitious goals in efficiency and innovation~\cite{dornheim2024neural,herrmann2024deep, jin2023recent}.

Several strategies have been developed for integrating machine learning into material behavior modeling:
\begin{description}
\item[Data-driven techniques] These techniques rely almost exclusively on data interpolation to identify constitutive models~\cite{leygue2018data, leygue2019non, stainier2019model}. They involve collecting large datasets of material responses under various conditions and using AI to interpolate and predict the material response. However, they often struggle with data quality and relevance issues, mainly when applied to real-world scenarios.
\item[Physics-informed methods] Physics-informed methods embed physical knowledge into the machine learning training process~\cite{raissi2019physics, abueidda2023enhanced, masi2022multiscale}. This integration addresses common issues such as physical inconsistency and limited generalization capabilities. By embedding known physical laws into the learning algorithm, these methods enhance the robustness and reliability of the resulting models.
\item[Physics-Augmented Neural Networks (PANNs)] PANNs take the integration of machine learning and physical knowledge one step further by embedding physical principles directly into the model architecture~\cite{linden2023neural, fernandez2021anisotropic, kalina2023fe}. These techniques generally use Input-Convex Neural Networks (ICNN)~\cite{amos2017input}, which enforce the convexity of the free energy to obey thermodynamic principles. Other thermodynamically consistent architectures exist~\cite{hernandez2021structure} in which a network is constructed to preserve the metriplectic structure of dissipative systems. This approach ensures that the models adhere to fundamental physical laws while leveraging the flexibility and learning capability of neural networks. PANNs have shown promise in learning complex and subtle behaviors from measurable data, such as anisotropic elasticity~\cite{klein2022polyconvex}, viscoelasticity~\cite{rosenkranz2024viscoelasticty}, or elastoplasticity~\cite{fuhg2023modular}.
\end{description}
On both numerical and theoretical fronts, research teams are continually enhancing AI constitutive model architectures, comparing methods, and conducting benchmarks~\cite{tacc2023benchmarking, rosenkranz2023comparative}. These AI-based strategies have achieved remarkable quality and diversity, demonstrating their potential in material science applications~\cite{fuhg2024review}.

Measurement techniques like Digital Image Correlation (DIC)~\cite{lucas1981iterative, sutton1983determination,hild2025future} and Digital Volume Correlation (DVC)~\cite{bay2008methods, buljac2018digital} measure displacement fields between a reference image, typically taken at the beginning of a test, and an image in a deformed configuration. By using global DIC~\cite{hild2012comparison, roux2024comprehensive}, the displacement field is expressed through a finite element mesh, which not only regularizes the problem but also facilitates a smooth integration with a finite-element software for identifying mechanical properties, \eg~through Finite Element Model Updating~\cite{mathieu2015estimation,roux2020optimal}. The extensive data obtained from these measurements provide the kinematic information essential for feeding data-driven models. It should be noted that the characteristics of measurement noise when processing through image correlation can be estimated~\cite{roux2024comprehensive}.

A significant limitation of data-based methods is their dependence on extensive stress measurements, which are inaccessible experimentally and must be estimated via ancillary measurements like force sensors. Typically, experiments are equipped with limited and simple sensors (uniaxial or biaxial). A supervised learning procedure cannot be employed. The Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) framework~\cite{flaschel2021unsupervised} addresses these challenges by requiring only global reaction forces and full-field displacement data, which are obtainable through modern measurement techniques~\cite{li2022equilibrium, jailin2024,bourdyot5358124learning}. EUCLID maps displacement measurements to residual balanced forces using a local form of equilibrium (internal and boundary conditions). The loss is written from the residual forces, enabling unsupervised identification of internal stresses and learning of a constitutive law.

Current non-supervised learning methods like EUCLID are highly sensitive to measurement noise, leading to significant biases in the identified material parameters~\cite{flaschel2021unsupervised, benady2024physics,jailin2025material}. This sensitivity arises because these frameworks typically assign equal weighting to all nodal contributions, disregarding the levels of uncertainty inherent in different measurement regions. In the same order, the mesh dependency introduced by uniformly weighted loss functions can result in suboptimal performance, for example, in complex geometries, with non-uniform meshes or when dealing with multiple experiments meshed differently. In their paper, Flaschel \etal~introduced white Gaussian noise on the deformation fields, resulting in significant model errors~\cite{flaschel2021unsupervised}. Given that measurement noise is unavoidable, it is crucial to control its impact through correction and compensation methods. In addition, a correct weighting of the input data related to their signal to noise ratio is important.

Recent advancements have focused on designing loss functions that are less sensitive to noise. A modified constitutive relation error framework for learning PANNs, more robust to noise impact, has been developed~\cite{benady2024unsupervised,benadynnmcre}. This framework utilizes a global noise model but does not account for uncertainty at each node. In addition to EUCLID, other frameworks also employ the local form of equilibrium with internal and boundary force loss terms, such as Equilibrium-based Convolution Neural Network (ECNN) methods~\cite{li2022equilibrium}. Furthermore, some frameworks have introduced virtual field losses to compensate for noise bias~\cite{conti2020data,platzer2021finite}. Nevertheless, the design of those fields is generally not related to the noise information. These approaches aim to mitigate the effects of noise, enhancing the robustness and accuracy of constitutive law identification in noisy environments.

Existing approaches often lack mechanisms to effectively propagate and account for measurement uncertainties through the non-linear transformations inherent in neural network models, limiting their robustness and generalization capabilities. These challenges necessitate the development of training specific noise management and uncertainty quantification strategies.

This work aims to evaluate the impact of measurement noise on models within the EUCLID (or NN-EUCLID) framework and proposes effective compensation strategies. By leveraging knowledge of noise characteristics, its influence on the constitutive law identification process can be compensated. A new functional is proposed, weighting all measurement contributions with its uncertainty, leading to a mesh independent learning with drastically damped noise bias.

The paper is structured as follows: Section~\ref{sec:METHOD} details the EUCLID formalism and measurement noise characteristics. Section~\ref{sec:APPLI} introduces a toy model of springs and a numerical test case contaminated with realistic noise. Section~\ref{sec:STRATEGY} outlines various noise compensation methods. Section~\ref{sec:RESULTS} compares the performance of these methods using metrics on the energy density function and predicted stress states. Finally, Section~\ref{sec:DISCUSS} discusses the advantages and limitations of the proposed approaches.

\section{Unsupervised constitutive learning}\label{sec:METHOD}

Consider a displacement field $\bm{u}(\bm{x})$ defined at spatial positions $\bm{x}$. This displacement field, obtained from Digital Image Correlation (DIC) measurement procedures, is assumed to be known on the sample surface. From global DIC~\cite{hild2012comparison}, the displacement field is decomposed over a finite element mesh with $N_n$ nodes, employing finite element shape functions $\Phi_i(\bm{x})$ such that:
\begin{equation}
\bm{u}(\bm{x}) = \sum_{i=1}^{N_n} \bm{u}_i \Phi_i(\bm{x}),
\end{equation}
where $\bm{u}_i$ are the nodal displacements. The deformation gradient field is approximated as follows:
\begin{equation}
\bm{F}(\bm{x}) = \bm{I} + \sum_{i=1}^{N_n} \bm{u}_i \bm \nabla \bm \Phi_i(\bm{x}),
\end{equation}
where $\bm{I}$ is the identity matrix and $\bm \nabla$ is the gradient operator defined in Lagrangian coordinates.

A hyperelastic constitutive model $\mathcal{M}$ maps deformation tensors $\bm{F}$ onto the first Piola--Kirchhoff stresses $\bm{P}(\bm{F}) = \mathcal{M}[\bm{F}]$. The sample is loaded with $N_R$ reaction forces, $\bm R_i$, for $i \in [1, N_R]$, evaluated from external force sensors. The EUCLID framework involves minimizing a loss function based on the force balance residuals, consisting of two terms: $\mathcal{L}_{\intern}$ for internal nodes and $\mathcal{L}_{\BC}$ for controlled boundary conditions. The $N_n$ degrees of freedom of the finite element mesh are similarly partitioned into an internal part $N_{\intern}$ (two times the number of internal nodes in 2D) and a controlled (measured forces) scalar boundary condition part $N_{\BC}$, with $N_n = N_{\intern} + N_{\BC}$.

In the absence of body forces, the nodal residual is expressed by integrating the stress over the reference domain $\Omega$ (in the following expression, implicit summation on repeated indices $j$ has been used):
\begin{equation}
f_i = \int_{\Omega} P_{ij} \nabla_j \Phi_i \dif V = \int_{\Omega} \mathcal{M}[\bm{F}]_{ij} \nabla_j \Phi_i \dif V.
\end{equation}

Let us consider the vector forces $\bm f_{\BC}$ being the sum over the concerned nodes of each nodal force balancing the imposed, $\bm f_{\intern}$ being the internal nodal forces, and measured, reaction forces~$\bm R$. The force residual can be written
\begin{equation}
\{ \bm f \}=
\begin{Bmatrix}
\bm f_{\intern} \\ \bm f_{\BC}
\end{Bmatrix} -
\begin{Bmatrix}
\bm 0 \\ \bm R
\end{Bmatrix}
\end{equation}
The entire process, including the derivation of displacement, inference in the constitutive model, and nodal force computation, is represented by $\Gamma$ such that: $\Gamma[\mathcal{M}, \bm u,\bm R] = \{ \bm f \}$.

In the EUCLID framework, the loss minimizes the $L_2$ norm of nodal force residuals, ensuring equilibrium. A general form of the loss function with a metric matrix $[\bm{A}]$ is defined as:
\begin{equation}\label{eq:InternalLoss}
\mathcal{L}_{[\bm A]}(t) = \braces[\big]{\bm{f}(t)}^\top [\bm{A}]^{-1} \braces[\big]{\bm{f}(t)}.
\end{equation}
Typically, all nodal contributions are weighted equally by choosing an identity metric, $[\bm A]=[\bm 1]$, or normalized by the number of degrees of freedom, $[\bm A]=(1/N) [\bm 1]$, where $N$ is the number of degrees of freedom. Although this choice appears natural, it introduces mesh dependencies and suboptimal weighting, which are addressed in subsequent sections. Note that a diagonal term $\omega$ as a weighting factor can also be added on the boundary condition terms, as proposed in~\cite{flaschel2021unsupervised}.

For the later noise compensation, a separation of the internal and boundary condition contribution is done. The framework can be split as a first loss term minimizing norm on the internal nodal force residuals (ensuring a balanced state), with a metric matrix $[\bm{A}']$:
\begin{gather}
[\hat {\bm A}] =
\begin{bmatrix}
[\bm{A}']& \bm{0}\\
\bm{0} & [\bm{A}''] \\
\end{bmatrix},
\\	\mathcal{L}_{\intern[\bm A']}(t)=\braces[\big]{\bm{f}_{\intern}(t)}^\top [\bm{A}']^{-1} \braces[\big]{\bm{f}_{\intern}(t)},
\end{gather}
and a second term expects to balance the reaction forces with a metric matrix $[\bm{A}'']$:
\begin{equation}
\mathcal{L}_{\BC[\bm A'']}(t)=\braces[\big]{\bm{f}_{\BC}(t) - \bm R}^\top [\bm{A}'']^{-1} \bracks[\big]{\bm{f}_{\BC}(t)- \bm R}.
\end{equation}
The total loss $\mathcal{L}_{l}$ is a sum of these two terms over all $N_t$ time steps
\begin{equation}
\mathcal{L}_{l[\bm A]} = \dfrac{1}{N_t} \sum_{t=1}^{N_t} \mathcal{L}_{[\bm A]}(t).
\end{equation}
The EUCLID method has been extended to incorporate neural network-based constitutive models, referred to as NN-EUCLID~\cite{thakolkaran2022nn}, which utilize the same loss function.

\section{Mechanical data}\label{sec:APPLI}

\subsection{Acquisition noise}

\paragraph{Image noise}
Displacement field measurements obtained from Digital Image Correlation (DIC) are subject to acquisition noise. Various noise models exist, adapted to different types of images (\eg~2D visible images, X-ray projections, CT reconstruction, MEB). Noise affecting pixel gray level values can be characterized, for instance, by evaluating its characteristics on repeated acquisitions.

\paragraph{Displacement uncertainty}
In DIC, not all degrees of freedom exhibit the same signal-to-noise ratio. Nodes near high gray-level image gradients in the image can be tracked more accurately than nodes in less contrasted areas. An extensive description of DIC noise can be found in~\cite{roux2024comprehensive}.

Considering image noise in the form of a centered white Gaussian term, with a standard deviation $\sigma_{\Imag}$, the nodal displacement noise is conveniently computed from the sensitivity fields, $\bm{S}(\bm{x})$, equal to the scalar product of the image gradient $\bm{\nabla} \Imag(\bm{x})$ by the DIC finite element shape functions $\Phi_i(\bm{x})$:
\begin{equation}
\bm S_i(\bm{x}) = \bm{\nabla} \Imag(\bm{x}) \Phi_i(\bm{x}).
\end{equation}
The covariance matrix of the nodal displacement $[\bm C_{u}]$ can be written as
\begin{equation}
[\bm C_{u}] = 2\sigma_{\Imag}^2 \parens[\big]{[\bm S]^\top[\bm S]}^{-1}.
\end{equation}

\begin{rema*}
A common procedure in full field measurement consists in performing DIC on multiple images with known displacement fields (typically with no motion) to empirically estimate the nodal uncertainty. This method involves capturing a series of static images and analyzing the variations in the measured displacements to estimate the noise.
\end{rema*}

Synthetic applications are used in the following. Noise is directly considered at the displacement level, chosen to be white and Gaussian distributed with a standard deviation $\sigma_u$, denoted $\mathcal{N}(0,\sigma_u^2)$, hence
\begin{equation}
[\bm C_{u}]=\sigma_{u}^2 [\bm 1].
\end{equation}
This assumption can be seen as a mean-field approximation (\ie~considering a constant image gradient), equally integrated across all nodal values, leading to $[\bm S]^\top[\bm S] \propto [\bm 1]$.

In this work, the covariance matrix $[C_u]$ is used only to generate synthetic noisy displacement fields with prescribed statistical properties. The identity assumption corresponds to uncorrelated and spatially uniform noise. Using a different covariance would simply modify the structure of the generated noise but would not affect the formulation of the proposed method. In experimental conditions, a set of repeated acquisitions naturally provides realistic noise realizations, making the explicit definition of $[\bm C_{u}]$ unnecessary.

For context, a standard deviation of 0.01 pixels corresponds to the uncertainty commonly achieved by efficient modern DIC techniques, generally considered very precise. In the applications, four noise amplitudes are considered in the $x$ and $y$ directions: null, low, medium, and intense, corresponding respectively to standard deviations $\sigma_u$ of 0, 0.01, 0.1, or 0.5 pixels.

\paragraph{Deformation uncertainty}
The noise at the deformation level is obtained by derivation of the noise at the displacement level~\cite{roux2024comprehensive}. Four components of the deformation tensor for each triangular linear element of the finite element mesh (subjected to 6 nodal displacement noise components) are obtained.

\paragraph{Force uncertainty}
The uncertainty associated with the force sensor can also be measured. This additional uncertainty can be combined with the nodal force uncertainty of the boundary condition. Standard force sensors typically have very low uncertainty~\cite{roux2024comprehensive}. The covariance $ [\bm{C}_R] $ of the reaction forces defined from the force sensor variance $\sigma_R^2$ and assumed independent measurements is written with $\bm{I}_{N_{\BC}}$ the identity matrix of size $N_{\BC}$:
\begin{equation}
{[\bm{C}_{R}^k] = \sigma_R^2
\begin{bmatrix}
\bm{0} & \bm{0}\\
\bm{0} & \bm{I}_{N_{BC}} \\
\end{bmatrix}}.
\end{equation}
Note that a specific variance associated to each sensor could be easily included, as well as extra dependencies such as a varying noise variance with the magnitude of the measurand.

\subsection{The effect of noise in EUCLID illustrated with a toy example}

This toy model development gives an overview of the noise impact and compensation strategy. Specifically, it demonstrates how measurement noise affects constitutive law identification within the EUCLID framework, and provides insight into how the proposed compensation strategies address these issues. This simplified setting clarifies the underlying mechanisms of noise-induced biases and the effectiveness of the noise mitigation approaches introduced in this work.

Consider a beam with a varying section $S(x)$ and length $L$ loaded with a force $R$, where the displacement at position $x$ is denoted by $u(x)$, with $u(0)=0$. The constitutive law (stress $\sigma$, strain~$\epsilon$) is expressed as:
\begin{equation}
\varepsilon = a\sigma + b\sigma^2,
\end{equation}
where $a$ and $b$ are material parameters. This constitutive law includes two parameters and a non-linear behavior. The stress and tangent stiffness $k=\frac{\dif \sigma}{\dif \varepsilon}$ can be extracted:
\begin{equation}
\sigma(\varepsilon,a,b) = \frac{\sqrt{a^2 + 4b\varepsilon} - a}{2b}, \qquad k(\varepsilon,a,b) = \frac{1}{\sqrt{a^2 + 4b\varepsilon}}.
\end{equation}
Considering a slowly varying beam section~($\frac{\dif S}{\dif x} \approx 0$), the balance equation is
\begin{equation}
\frac{\dif \parens[\big]{\sigma S(x)}}{\dif x} = S(x) \frac{\dif \sigma}{\dif x} =0.
\end{equation}
Let's consider noise in the displacement:
\begin{equation}
v(x) = \eta(x) + u(x), \qquad w(x)=\frac{\dif v}{\dif x},
\end{equation}
with $\eta(x)$ a white and Gaussian noise term, centered with a variance $\gamma_\eta$. The balance equation now is
\begin{equation}
S(x) \frac{\dif \sigma}{\dif x} = S(x) \frac{\dif \sigma}{\dif \varepsilon} \frac{\dif \varepsilon}{\dif x} = S(x) k(w,a,b) \frac{\dif w}{\dif x} = S(x) k(w,a,b)\frac{\dif^2\eta}{\dif x^2}.
\end{equation}
Finally, the boundary condition loss and the expectation of the internal one is obtained from the integration. The internal equilibrium involves the second derivative of the noisy displacement (in a discrete noise formulation), $\dif^2\eta/\dif x^2$. When taking the expectation, the squared term $E\bracks[\big]{(\dif^2\eta/\dif x^2)^2}$ yields a constant variance $\gamma_\eta$, while the local section $S(x)$ remains within the integral:
\begin{gather}
	\mathcal{L}_\BC = \bracks[\big]{S(0) \sigma\parens[\big]{w(0), a, b}-R}^2 + \bracks[\big]{S(L) \sigma\parens[\big]{w(L), a, b}-R}^2,
\\	E[\mathcal{L}_\intern] = \gamma_\eta \int S(x)^2 k(w, a, b)^2 \dif x.
\end{gather}
At first order in the zero-mean noise on $w$, the expectation of the stiffness $k(w,a,b)$ remains unchanged.

While the boundary condition term tends to be minimal when the correct solution is found (\eg~equilibrated with the applied forces), $\mathcal{L}_{\intern}$ shows a minimum when internal forces are balanced. This highlights a degenerated solution when stiffness is zero. This degenerate solution seeks to reduce noise impact when mapped into noise-induced forces. Thus, the two losses do not share a similar optimum, leading to a systematic bias.

Concerning the internal loss term, it is composed of three factors: the noise variance $\gamma_\eta$ shows a dependency on the noise level, $S(x)^2$ suggests a geometry impact, and $k(w, a, b)$ highlights that the optimal solution tends to lower the stiffness. On the opposite, the boundary condition term does not lead to a noise bias as the solution is driven towards the force measurement $R$.

For the numerical applications, default dimensionless parameters are chosen, the external force is $R=1$, $L=2000$, the varying section is $S(x)=\abs\bracks[\big]{4({2x}/{L}-1)}+1$, the noise standard deviation is $\gamma_\eta=0.005$. The ground truth (GT) material parameters are chosen to be $a=1$, and $b=0.75$.

As it stands, the (NN-)EUCLID loss presents several noise-related drawbacks:
\begin{itemize}
\item noise induces a bias in $\mathcal{L}_{\intern}$, resulting in an underestimation of the true stiffness;
\item the current procedure is mesh-dependent as it aggregates contributions from degrees of freedom; finer meshed areas will have a larger impact than coarse meshed areas; combining different training experiments with different meshes will not result in a consistent loss of weight;
\item the current loss does not account for measurement noise characteristics and weights all nodal contributions equally.
\end{itemize}

\subsection{Test-cases: training and validation}

The application of this paper is based on an isotropic hyperelastic behavior following a neo-Hookean law. Considering the deformation invariants defined as $I_1 = \tr(\bm C)$ and $J = \det(\bm C)^{1/2}$, with $\bm C = \bm F^\top \cdot \bm F$, the strain energy density function can be written as:
\begin{equation}
W = \dfrac{\mu}{2} (I_1 - 3) - \mu \log(J) + \dfrac{\lambda}{2} \log(J)^2.
\end{equation}
The material parameters are chosen arbitrarily, $E=100$~MPa, and $\nu=0.3$ giving the Lamé coefficients $\mu=E/\parens[\big]{2(1+\nu)}$ and $\lambda=\nu E/\parens[\big]{(1+\nu)(1-2\nu)}$.

Two different structures are used, one for the training and one for the validation/test. The structures are shown in Figure~\ref{fig:samples}. The training structure is a rectangular sample with an inner hole loaded with an axial force. This simple test could be easily obtained experimentally with a simple uni-axial testing machine. The validation structure is a V-shape with multiple holes loaded with an oblique force.
All simulations are performed in two dimensions under the plane-stress assumption, consistent with surface DIC measurements where the out-of-plane strain component cannot be directly measured.

Figure~\ref{fig:samples} illustrates the training and validation structures with the geometrical characteristics. The thickness of the two samples is 10~mm. The mesh characteristics are presented in Table~\ref{tab:meshes}. All simulations are based on FEniCS~\cite{alnaes2015fenics} software. The training structure is loaded continuously up to a tensile force of $R_{t}=(0,2000)$~N. 10 loading steps, regularly distributed from 0 up to the maximum load, are acquired. The validation structure is loaded proportionally both in tension and shear with a unique force: $R_v=(500,1000)$~N.

\begin{figure}[h]
{\includegraphics[width=0.6\textwidth]{samples.png}}
\caption{Application samples. Left: Training sample. Right: Validation sample.}\label{fig:samples}
\end{figure}

\begin{table}[ht]
\caption{Mesh characteristics of the two simulated structures.}\label{tab:meshes}
\begin{tabular}{lcccc}
\hline
\\[-2.7ex]
\hline
& \textbf{Number of nodes} & \boldmath{$N_{\textbf{BC}}$} & \textbf{Number of elements} & \boldmath{$l_e$} \textbf{[pixels]} \\
\hline
\\[-2.7ex]
\hline
\textbf{Training} & 934 & 21 & 1708 & 25 \\
\textbf{Validation} & 1678 & -- & 3050 & 25 \\
\hline
\end{tabular}
\end{table} 

\paragraph{Realistic noise amplitude}
Selecting an appropriate noise amplitude for displacement measurements necessitates careful consideration of the specific application scenario. This involves choosing a representative sample and loading condition that generates measurable displacement amplitudes, as well as selecting a measurement setup with a resolution that mimics typical experimental conditions. In the current application, it is considered that the imaged sample is imaged with a length of 1000 pixels, a scale easily achievable with standard laboratory cameras.

In this study, synthetic noise is introduced into the displacement field of the training structure to assess its impact on the learned behavior. Four distinct noise levels are evaluated: no noise, low noise, medium noise, and intense noise~\cite{hild2025future}, corresponding to uncertainties in the nodal displacement field of $\sigma_u = 0$, $0.01$, $0.1$, $0.5$ pixels, respectively. While no additional force uncertainty is incorporated in this example, it is important to recognize that the boundary conditions are inherently influenced by the displacement uncertainty.

This approach ensures that the noise amplitudes considered are both realistic and relevant to practical experimental setups, thereby providing meaningful insights into the robustness and accuracy of the constitutive law identification under varying noise conditions.

Figure~\ref{fig:F} shows the four deformation components polluted by a 0.5~pixel displacement noise on the training structure.
\begin{figure}[t]
{\includegraphics[width=0.7\textwidth]{F_train_0.5.png}}
\caption{Training data at maximum load: deformation components polluted by a $\sigma_u=0.5$~pixel noise.}\label{fig:F}
\end{figure}

\subsubsection*{PANN architecture model}\label{sec:PANN}

The neural network model used in this application is a PANN, $\mathcal{M}_{\PANN}$, designed for isotropic hyperelastic applications. The strain energy density is written as a function of the Green--Lagrange strain tensor $\bm E = (\bm C-\bm 1)/2$, which represents a sufficient condition for objectivity ($\bm C = \bm F^T \cdot \bm F$ being the right Cauchy--Green deformation tensor). Following~\cite{as2022mechanics,linden2023neural}, the strain energy density is written as:
\begin{equation}
W(\bm F) = W_{\ICNN}\parens[\big]{\bm E(\bm F)} + W_{\bm F=\bm 1} + \bm H:\bm E.
\end{equation}
The first term corresponds to the energy output of an input convex neural network $W_{\ICNN}\parens[\big]{\bm E(\bm F)}$)~\cite{amos2017input}, and the last two terms correspond to energy and stress normalization. Various physical constraints are implemented following reference papers in the literature~\cite{klein2022polyconvex,linden2023neural}. We follow the architecture proposed in~\cite{as2022mechanics}:
\begin{itemize}
\item Manual definition of the deformation invariants before the first layer to enforce isotropy. Note that the invariant $I$ is convex in $\bm F$, the invariant $I_2 = \tr\parens[\big]{\cof(\bm C)}$ is convex in cof$(\bm F)$, and the invariant $I_3=\det(\bm F)^2$ is convex in $\det(\bm F)$. The invariants layer receives the deformation gradient $\bm F$ from the input layer and computes strain invariants.
\item Thermodynamic consistency: the first Piola--Kirchhoff stress is obtained by derivation of the strain energy density function $W$: $\bm P(\bm F)=\frac{\partial W(\bm F)}{\partial \bm F}$. An intermediate output of the neural networks is hence a scalar value $W$.
\item Stability of the potential is obtained by the polyconvexity with respect to $\bm F$, $\det(\bm F)$ and $\cof(\bm F)$ by using the convex and non-decreasing \emph{softplus} activation functions for the hidden layers~\cite{amos2017input}: $\mathcal{A}(y)=\log(1+e^y)$ and positive model weights $\bm \omega>0$. In order to satisfy the polyconvexity constraints, all the invariants $I$, $I_2$, and $J$ are typically modeled as convex and monotonically non-decreasing. However, due to this overconstraint (the monotonically increasing constrain is not required in $J$), these models are unable to represent negative stresses~\cite{klein2022polyconvex,linden2023neural}, which leads to the addition of an extra invariant $-J$ to the input set. The relationship between two hidden layers, composed of $N_n$ neurons can hence be written, with $\bm \omega^i$ and $b^i$ the weights and bias of layer $i$, and $\bm z^{i-1}$ and $\bm z^{i}$ the neuron values of layers $i-1$ and $i$
\begin{equation}
\bm z^{(i)}=\mathcal{A}\circ \bracks[\Bigg]{\sum_{k=1}^{N_n}\bm \omega_{jk}^{(i)} z_k^{(i-1)}+ b_j^{(i)}}.
\end{equation}
\item Normalization of the energy and stresses: $\bm P(\bm F=\bm 1)=\bm 0$, and $W(\bm F=\bm 1)= 0$. This normalization is obtained with the two last terms in the energy formulation, with
\begin{equation}
W_{\bm F=\bm 1}=-W_{\ICNN}\parens[\big]{\bm E(\bm 1)}, \qquad \bm H = -\left. \dfrac{\partial W_{\ICNN}}{\partial \bm{F}}\right|_{\bm F = \mathbf{I}}.
\end{equation}
Note that this stress normalization term does not guarantee polyconvexity as $E$ is not convex on every of its component. This normalization choice follow~\cite{as2022mechanics} and has been enhanced in recent literature~\cite{linden2023neural}.
\end{itemize}

The core architecture is composed of 5 dense layers with 16 neurons, leading to 1,233 trainable weights. The optimization is performed with an ADAM optimizer~\cite{kingma2014adam} with a learning rate following a predefined schedule. Weights are initialized randomly. The training on each structure takes approximately 10 minutes on a laptop with an 8Gb GPU unit (Nvidia RTX-A2000). The neural network architecture is developed using the TensorFlow library in Python.

\section{Noise compensation strategy}~\label{sec:STRATEGY}

To address the previous limitations, several complementary strategies are introduced. First, a reformulation of the loss is proposed to reduce its sensitivity to high-frequency noise, notably by applying spatial filtering to the residual force fields. Second, a statistical characterization of the measurement noise is leveraged to inform a covariance-based weighting of the residuals, enabling both noise-aware and mesh-independent learning. A systematic bias compensation term is introduced to explicitly counteract the noise-induced distortion in the internal loss. Together, these approaches provide a coherent framework for reducing noise effects while preserving the physical interpretability and generalization capacity of the learned constitutive laws.

\goodbreak

The two complementary noise management strategies are hence presented:
\begin{description}
\item[reducing noise impact] reformulating a new internal loss norm with lower sensitivity to the force noise, thereby reducing the bias term impact (Section~\ref{ssec:gauss});
\item[compensating noise impact] tracing the noise throughout the computation, it is proposed either to account carefully for the uncertainty or to compensate for the noise-induced bias (Section~\ref{ssec:cov}):
\begin{itemize}
\item accounting for measurement uncertainty considering the uncertainty of the measurement and an appropriate node weighting scheme (allowing mesh independence and managing the weighting ratio between internal and boundary condition terms); this requires propagating the measurement uncertainty into the non-linear behavior model;
\item removing the systematic noise bias: since the (internal) functional is affected by the noise effect, an opposite potential is built to mitigate the noise effect.
\end{itemize}
\end{description}
A simple representation of the second approach is proposed in Figure~\ref{fig:scheme}.

\begin{figure}[h!]
\includegraphics[width=0.85\textwidth]{image_global.png}
\caption{Representation of the noise compensation approach.}\label{fig:scheme}
\end{figure}

\subsection{Reducing noise impact}\label{ssec:gauss}

The first solution involves focusing on low-frequency components to reduce noise impact. In a similar spirit to the virtual field method~\cite{pierron2012virtual}, a low-frequency mechanically consistent extractor could be used to compute the loss term. Emphasis is placed on low-frequency filtering (less affected by high-frequency noise) by convolving the internal nodal forces residual with a Gaussian kernel $G_\lambda$, characterized by a scale $\lambda$. This filtering is applied to the residual field and not to the original displacement quantity in order to play a role on the unbalanced nodal forces. Note that only the internal nodes are filtered in this process. This prevents the boundary condition error from spreading on the sample. This first filtering compensation is the only place in this paper where an internal/boundary condition split is performed.

Note that even for a linear elastic model, the stiffness operator and the Gaussian smoothing operator do not commute, as they act on different quantities (derivation versus spatial averaging). For instance, an infinite smoothing would yield a constant displacement field and thus no internal forces after applying the tangent operator, while smoothing the resulting unbalanced forces would not cancel them.

Let $\bm{x}_i$ denote the location of node $i$ in the mesh. The filtered residual nodal forces $\widehat{ \bm{f}}_i$ is computed as
\begin{equation}
\widehat{ \bm{f}}_i=
\frac{\sum_j\exp\parens[\big]{-(\bm x_i -\bm x_j)^2/\lambda^2} \bm{f}_j} {\sum_k\exp\parens[\big]{-(\bm x_i -\bm x_k)^2/\lambda^2}}
\equiv [\bm G_\lambda]\{ \bm f\}.
\end{equation}
The modified internal loss term is hence given by
\begin{equation}
\hat{\mathcal{L}}_{\intern[\bm A]'}(t)= \{ \bm{f}_\intern\}^\top [\bm{G}_\lambda]^\top [\bm{A}']^{-1} [\bm{G}_\lambda] \{ \bm{f}_\intern\}.
\end{equation}
Then, with
\begin{equation}
[\hat {\bm A}]^{-1} =
\begin{bmatrix}
[\bm{G}_\lambda]^\top [\bm{A}']^{-1} [\bm{G}_\lambda] & \bm{0}\\
\bm{0} & [\bm{A}'']^{-1} \\
\end{bmatrix},
\end{equation}
the whole loss can be written:
\begin{equation}
\hat{\mathcal{L}}_{[\hat{\bm A}]}(t)= \{ \bm{f}\}^\top [\hat{\bm{A}}]^{-1} \{ \bm{f}\}.
\end{equation}

Using large Gaussian kernel scales allows unbalanced nodal forces to be evaluated over extensive regions. In the following application, different lengths for the Gaussian kernels, $\lambda=\alpha l_e$, are considered, where $l_e$ represents the average element length (computed as the average distance between nodes in an element), and $\alpha$ is a dimensionless scalar value chosen in $\{0.5, 1, 2, 3\}$.

It is noteworthy that pre-processing the input Digital Image Correlation (DIC) fields to mitigate noise effects is a viable strategy, as discussed in Flaschel et~al.~(2021)~\cite{flaschel2021unsupervised}. However, this approach was not adopted in the current study to preserve the integrity of the measured field information.

\subsection{Correcting noise impact}\label{ssec:cov}

\subsubsection*{Monte Carlo estimation of the force covariance}

As the PANN model $\mathcal{M}$ accounts for a non-linear transformation, propagating image (or displacement) noise to internal force uncertainty is not straightforward. To address this, the covariance matrix is estimated using Monte Carlo (MC) simulations and updated every $N_{MC}$ epochs.

To formalize the MC procedure with synthetic noise generation, we follow the following steps at training epoch $k$. When computing the force covariance, that of the stroke sensor can also be added. Note that the covariance now depends on the loading step as it is related to the material states. For a complete experiment with multiple loading steps, the covariance should be estimated at each step $t$. From the generation of synthetic noise on samples, that can be obtained from a noise model (\eg~Gaussian), or repeated acquisitions:
\begin{equation}
{b}_i \sim \mathcal{N}(0, \sigma_u), \quad i = 1, 2, \ldots, M.
\end{equation}
The inference of the noisy signal through the model at epoch $k$, $\mathcal{M}^k$ gives:
\begin{equation}
\bm{z}^k(\bm x,\bm{u},\bm{b},\mathcal{M}^k,\bm R,t) = \Gamma \bracks[\big]{\mathcal{M}^k, \bm u(\bm x,t) + \bm{b}(\bm x), \bm R(t)} - \Gamma \bracks[\big]{\mathcal{M}^k, \bm u(\bm x,t), \bm R(t)},
\end{equation}
and allows to compute the final covariance matrix, averaged over noise realizations, considering a centered force noise term:
\begin{equation}
\bracks[\big]{\bm C_{\bm{f}}^k(t)} = [\bm{C}_{R}^k] +\dfrac{1}{N_n}\angles[\big]{(\bm{z}_i^k)(\bm{z}_i^k)^\top}_{\bm b}.
\end{equation}

In summary, the process involves generating Gaussian noise samples, perturbing the initial signal with each noise sample, passing the perturbed signals through the model, and calculating the empirical covariance matrix of the model outputs. Note that this empirical procedure may result in a covariance matrix that is difficult to invert (large matrix with 4\textsuperscript{th} order coupling), necessitating the introduction of a Tikhonov regularization term.

At the beginning of the learning phase, the behavior model is unknown. Therefore, the model covariance is initialized with the identity matrix: $C_{\mathbf{f}}^{k=0}(t)=[\bm{1}]$, or by computing it without MC sampling and updating it at selected epochs to reduce computation time per epoch.

It is now proposed to use the inverse covariance matrix as the metric $[\bm{A}]$ in $\mathcal{L}_{\intern}$ (see equation~\ref{eq:InternalLoss}):
\begin{equation}
\mathcal{L}_{[\bm{C}^k_{\bm{f}}]}(t) = { \bm{f}(t) }^\top [\bm{C}_{\bm{f}}^k]^{-1} { \bm{f}(t)}.
\end{equation}
This covariance term allows for the weighting of the nodes based on their uncertainty, derived from measurements (gray level or displacement) to force. A node located in an area with low image texture will exhibit high uncertainty in the DIC analysis. The resulting inverse covariance adequately weights the kinematic measurements.

\begin{rema*}
This remark addresses the coupling of an appropriate covariance weighting strategy with the previously discussed filtering approach. First, as highlighted in this section, weighting nodal information by its associated uncertainty provides an optimal metric. Consequently, any linear operation applied to the nodal quantity is directly transferred in the covariance computation and, hence vanishes. Therefore, there is no need to use both techniques at the same time. With $\{\widehat{ \bm{z}} \}$ being the filtered force noise,
\begin{equation}
\{\widehat{ \bm{z}} \} = [\bm G_\lambda]\cdot\{ \bm z\},
\end{equation}
the noise filtered covariance reads $\bracks[\big]{\widehat{\bm C_{\bm{f}}^k(t)}} = [\bm G_\lambda] \bracks[\big]{\bm C_{\bm{f}}^k(t)} [\bm G_\lambda]^\top$, hence:
\begin{equation}
\begin{split}
\hat{\mathcal{L}}_{[\widehat{\bm C_{\bm{f}}^k(t)}]}(t)
	& = \{ \bm{f}\}^\top [\bm{G}_\lambda]^\top \parens[\big]{[\bm G_\lambda] \bracks[\big]{\bm C_{\bm{f}}^k(t)} [\bm G_\lambda]^\top}^{-1} [\bm{G}_\lambda] \{ \bm{f}\},
\\	& = \braces[\big]{\bm{f}(t)}^\top [\bm{C}_{\bm f}^k]^{-1} \braces[\big]{\bm{f}(t)},
\\	& = {\mathcal{L}}_{\bracks{\bm C_{\bm{f}}^k(t)}}(t)
\end{split}
\end{equation}
\end{rema*}

\subsubsection*{Bias correction due to noise}\label{ssec:Nbias}

A noise bias compensation scheme is proposed to cancel the bias computed from a known noise. The noise characteristics can be extracted experimentally (by performing multiple acquisitions of a similar state) or from a general noise model. Since the mechanical model is non-linear, the impact of noise depends on the material state. This bias can be isolated and written considering noise terms, $\bm{b}$, applied to the original displacement, and averaged over noise realization
\begin{equation}
\mathcal{L}_{\bias}\parens[\big]{\Gamma(\mathcal{M}, \bm R)} =
\angles[\big]{\mathcal{L}\bracks[\big]{\bm{z}^k(\bm x,t)}}_{\bm b}.
\end{equation}
Additionally, if the chosen metric $[\bm A]$ is the usual identity matrix, then
\begin{equation}
\mathcal{L}_{[\bm 1]\bias}= \trace\parens[\big]{[\bm C^k_{\bm f}](t)}.
\end{equation}
In contrast, if it is chosen to be the inverse covariance of $ f$, as above suggested, then the simple following result is obtained
\begin{equation}
\mathcal{L}_{[\bm C^k_{\bm f}]\bias}=1,
\end{equation}
irrespective of the mesh and load history. Let us stress that the target is to reach unbalanced forces such that their variance reaches precisely the value expected from noise.

The general corrected functional $\mathcal{L}_{\corr}$ is hence constructed with the bias compensation term, and using a positive (L1) norm to ensure the stability, such that:
\begin{equation}
\mathcal{L}_{[\bm A]\corr} = \abs
\parens[\big]{\mathcal{L}_{[\bm A]}\bracks[\big]{\Gamma(\mathcal{M}, \bm u, \bm R)} - \mathcal{L}_{[\bm A]\bias}\bracks[\big]{\Gamma(\mathcal{M}, \bm R)}},
\end{equation}
or using the inverse covariance matrix metric
\begin{equation}
\mathcal{L}_{[\bm C_{f}]\corr} = \abs
\parens[\big]{\mathcal{L}_{[\bm C_{f}]}\bracks[\big]{\Gamma(\mathcal{M}, \bm u, \bm R)} - 1}.
\end{equation}
With this corrected functional, the model is trained to ensure equilibrium (weak form) while the noise bias is opposed by its evaluation performed using synthetic noise.

\begin{rema*}
To ensure numerical stability and optimize computational resources, an approximated version of the proposed loss is employed. It consists of approximating the covariance function by the diagonal terms only:
\begin{equation}
\mathcal{L}_{[\diag\bm C_{f}]\corr} = \abs\parens[\big]{\{ \bm{f}\}^\top \bracks[\big]{\diag\ \bm C^k_{\bm f} (t)}^{-1} [\bm{f}\} -1}.
\end{equation}
By considering only the diagonal elements of the covariance matrix, the computational complexity is significantly reduced, and complex matrix inversion is avoided. With noise information derived from experimental procedures --- such as extracting fields from a hundred repeated images --- this approach becomes feasible without necessitating thousands of computations. 200 or 50 repeated acquisitions, which are experimentally accessible, allow respectively estimating for each node the expected value of a Gaussian noise distribution with a standard deviation error of 5\% and 10\%. Although this method disregards nodal coupling, it remains a fair approximation under the assumption that the image speckle is homogeneous.
\end{rema*}

\subsection{Proposed losses}

The different loss expressions are summarized in Table~\ref{tab:losses}.

\begin{table}[ht]
\caption{Internal loss choices.}\label{tab:losses}
\begin{tabular}{l|c}
\hline
\\[-2.7ex]
\hline
\textbf{Name} & \textbf{Expression} \\
\hline
\\[-2.7ex]
\hline
EUCLID~\cite{flaschel2021unsupervised} & $\{\bm{f}\}^\top [\bm{1}]^{-1} \{\bm{f}\}\vphantom{\parens[\big]{C^k_{\bm f}}}$ \\
Reducing noise impact & $\{ \bm{f}\}^\top [\bm{G}_\lambda]^\top [\bm{1}]^{-1} [\bm{G}_\lambda] \{ \bm{f} \}\vphantom{\parens[\big]{C^k_{\bm f}}}$\\
Correcting noise impact & $\abs\parens[\big]{\{ \bm{f}\}^\top \bracks[\big]{[\diag\ \bm C^k_{\bm f} (t)}^{-1} [\bm{f}\} -1}$\\
\hline
\end{tabular}
\end{table}

To prevent the boundary condition error from propagating through the sample and slowing down convergence, it is decided to separate it from $\mathcal{L}_{\intern}$. The filtering approach is thus only applied to the internal nodes. Noise has a lower impact on the boundary condition as it is already written as a sum of BC nodes.

\subsection{Metrics}

To compare the different approaches, the model is trained five times for each configuration. Since it is a synthetic test case, the ground-truth free energy function is known, $W_{\GT}$, and can also be estimated with the constitutive model: $W_{\PANN}$. To evaluate performance, an energy error metric derived from the $L_2$ norm is proposed, computed on the validation structure, $\bm{u}_{\valid}$:
\begin{equation}
m_{\valid} = \norm[\big]{W_{\PANN}(\bm{u}_{\valid}) - W_{\GT}(\bm{u}_{\valid})}^2.
\end{equation}
Note that the ground truth is the one obtained from the theoretical model. An un-noisy computation modeled with a PANN will also serve as a reference for the best achievable approximation with the chosen machine learning architecture.

A key point to be highlighted is that the validation structure is not affected by noise. Consequently, this metric provides a direct and unbiased measure of the model ability to replicate the ground-truth behavior. A relative improvement quantity, in percent, is also computed from the reference, non-noisy case, to appreciate the improvement. This metric captures two critical aspects of model performance:
\begin{itemize}
\item accuracy of the learned behavior: to evaluate how well the model reproduces the ground-truth free energy behavior, both for the training and the validation data;
\item generalization capability: by assessing the model on a validation structure that differs from the training structure, this metric highlights the model's ability to generalize to unseen data.
\end{itemize}

\begin{rema*}
For the simple toy model, a visual evaluation of the learned stiffness over all the loadings will be proposed to better appreciate the two parameters, $a$ and $b$, impacts.
\end{rema*}

\section{Results}\label{sec:RESULTS}

\subsection{Toy model results}

The toy model application is evaluated 50 times to appreciate the variability due to noise pollution and compensation. The scenario with the 1D Gaussian filter is performed using a characteristic length of 1 neighbor.

\begin{figure}[!h]
\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{toy_ini.png}
\caption{EUCLID}\label{fig:toy_a}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{toy_filtered.png}
\caption{Filtered}\label{fig:toy_b}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{\textwidth}
\includegraphics[width=\textwidth]{toy_corr.png}
\caption{Corrected}\label{fig:toy_c}
\end{subfigure}
\caption{Toy model: comparison of loss functions without noise and with Gaussian noise. \eqref{fig:toy_a}~Uncorrected loss functions show significant bias with noise. \eqref{fig:toy_b}~Filtered loss functions effectively reduce bias. \eqref{fig:toy_c}~The proposed loss compensation model.}\label{fig:toy}
\end{figure}

Figure~\ref{fig:toy} presents the internal (left), boundary condition (center) and total loss (right) for the three approaches: EUCLID (a), Gaussian filtering (b), and proposed loss (c). When noise is introduced, the boundary condition loss exhibits minor variations but still maintains a minimum close to the ground truth parameter values, only slightly influenced by measurement noise. In contrast, the internal loss term $\mathcal{L}_{\intern}$ displays an intense gradient that tends to increase both $a$ and $b$, thus a reduction in stiffness. Applying Gaussian filtering with a characteristic scale set to the first neighbor mean distance, or $\alpha = 1$ neighbor, significantly improves the accuracy. With the proposed noised-compensated loss, the obtained model parameters are close to the GT. In this toy model, the filtering and the noise bias compensation appear to display comparable improvements. Since the boundary condition loss term is less sensitive to noise, increasing its weight in the computation could be a potential solution (as suggested in the original EUCLID paper). However, this approach would lower the benefit of the rich full-field information measured within the specimen.

Because the proposed law couples $a$ and $b$ and it may be difficult to interpret them as stiffnesses, all the learned constitutive laws are presented in Figure~\ref{fig:toy_laws}. As anticipated, the stiffness derived from the noise-polluted functional is significantly lower than the ground truth (GT) value. Both the filtered and corrected approaches exhibit a marked improvement, closely aligning with the GT behavior. Specifically, the filtered approach effectively dampens the noise effects, resulting in stiffness measurements that are substantially closer to the GT, although a slight reduction in stiffness remains. In contrast, the corrected approach achieves results that are tightly centered around the GT values, demonstrating a more accurate compensation of noise-induced biases. This indicates that while Gaussian filtering primarily attenuates the noise impact, the corrected method not only reduces the noise influence but also effectively offsets the resulting bias, ensuring a more faithful reproduction of the true constitutive behavior.

\begin{figure}[!h]
\includegraphics[width=.7\textwidth]{toymod.png}
\caption{Toy model: reconstruction of the identified constitutive models. The black dashed curve corresponds to the ground truth. The large drawn width represent one standard deviation computed from the 50 runs.}\label{fig:toy_laws}
\end{figure}

\subsection{Noise impact on the test case}

Figure~\ref{fig:stresses} shows the relative error (expressed in percent) between the computed and predicted fields of the Piola--Kirchhoff stress. The noise directly influences the learned behavior, as the predicted stress states show. The color dynamics is chosen at $\pm 5\%$ of the initial stress dynamics, showing important stress errors with a noise of 0.1\%, which is already considered an acceptable correlation procedure.

\begin{figure}[t!]
\begin{subfigure}{.48\textwidth}
{\includegraphics[width=\textwidth]{stress_0_0.png}}
\caption{$\sigma_u=0$~pixel}\label{fig:stresses_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
{\includegraphics[width=\textwidth]{stress_0.01_0.png}}
\caption{$\sigma_u=0.01$~pixel}\label{fig:stresses_b}
\end{subfigure}

\vskip .5\baselineskip

\begin{subfigure}{.48\textwidth}
{\includegraphics[width=\textwidth]{stress_0.1_0.png}}
\caption{$\sigma_u=0.1$~pixel}\label{fig:stresses_c}
\end{subfigure}
\hfill
\begin{subfigure}{.48\textwidth}
{\includegraphics[width=\textwidth]{stress_0.5_0.png}}
\caption{$\sigma_u=0.5$~pixel}\label{fig:stresses_d}
\end{subfigure}
\caption{First Piola--Kirchhoff stresses error with respect to the ground truth. The model training is performed with different noise levels: \eqref{fig:stresses_a}~no noise: $\sigma_u=0$~pixel, \eqref{fig:stresses_b}~low noise: $\sigma_u=0.01$~pixel, \eqref{fig:stresses_c}~medium noise: $\sigma_u=0.1$~pixel, and \eqref{fig:stresses_d}~intense noise: $\sigma_u=0.5$~pixel. The color bar is set to $\pm 5\%$ of the ground truth stress dynamics in all figures.}\label{fig:stresses}
\end{figure}

Concerning the energy error metric, Figure~\ref{fig:errors_ini} demonstrates the sensitivity of the standard EUCLID method to varying noise levels, presented through a log-log plot of metrics across low (0.01 pixel) to intense (0.5 pixel) noise amplitudes. A geometric series of twelve noise amplitudes (regularly distributed on a log scale) is chosen. The results reveal a sharp escalation in the metric as noise intensity increases. This performance degradation underscores the limitations of the EUCLID method in preserving accuracy as noise levels escalate, thereby emphasizing the imperative for enhanced noise management strategies.

\begin{figure}[t]
{\includegraphics[width=300pt]{results_EUCLID.png}}
\caption{Impact of the training noise, evaluated on the validation structure using the energy error metric.}\label{fig:errors_ini}
\end{figure}

\begin{figure}[h!]
{\includegraphics[width=360pt]{results_EUCLID_proposed.png}}
\caption{Comparison of EUCLID and proposed methods under varying noise levels, evaluated using the energy error metric. The proposed method consistently achieves lower error metrics across all noise levels, demonstrating enhanced robustness against noise.}\label{fig:results_EUCLID_proposed}
\end{figure}

\subsection{Test case results}

Figure~\ref{fig:results_EUCLID_proposed} compares the standard EUCLID method with the proposed approach across different noise scenarios (0, 0.01, 0.1, and 0.5 pixels). All computations are performed 10 times to evaluate the variability when considering different random aspects for the noise generation and the compensation strategy.

The reference shows that the model architecture can effectively catch the prescribed hyperelastic behavior up to an energy error metric of 1.2. The proposed method consistently outperforms EUCLID as noise levels increase. For a small noise of 0.01, the EUCLID and noise-corrected versions show similar results. At a noise level of 0.1 pixels, EUCLID results begin to deviate significantly while the proposed method keeps much lower values. For higher noise levels, such as 0.5 pixels, the proposed method maintains relatively low metrics, whereas EUCLID exhibits substantial error.

\begin{figure}[!h]
{\includegraphics[width=360pt]{gaussian_results.png}}
\caption{Performance comparison of EUCLID, Gaussian smoothing, and the proposed method under medium and intense noise levels, respectively 0.01, 0.1, and 0.5 pixels. The bars, in a vertical log plot, represent the mean values (over 10 runs) and are associated with the standard deviation.}\label{fig:gaussian_results}
\end{figure}

Figure~\ref{fig:gaussian_results} further illustrates the performance comparison, showing that Gaussian smoothing with kernel lengths of 0.05 and 0.1 (one or two neighboring nodes) highly reduces noise impact.
At a noise level of 0.1 pixels, the new Table~\ref{tab:noise_model_performance} reveals that EUCLID’s error metric increases to 729, while the proposed method is much lower at 67, representing a 90.9\% improvement. Still, Gaussian smoothing with small kernel lengths can achieve even lower metrics (38 and 36) and correspondingly higher relative improvements (95\%). For higher noise levels, such as 0.5 pixels, the proposed method remains much better than EUCLID (186 vs.\ 2584).

Table~\ref{tab:noise_model_performance} summarizes the noise impact on model performance metrics and the relative improvements achieved by the proposed method compared to EUCLID and Gaussian smoothing.
At low noise levels~(0.01 pixels), the proposed method shows marginal improvement. However, at medium (0.1 pixels) and intense noise levels~(0.5 pixels), the proposed method exhibits substantial relative improvements of up to 93.2\%, highlighting its effectiveness in mitigating noise-induced biases.
The table also indicates that at 0.1 pixels of noise, Gaussian smoothing methods reach up to 95.1\% improvement, whereas the proposed method is at 90.9\%. For a 0.5 pixels noise, the Gaussian-0.05 approach reaches a 94.9\% improvement, slightly exceeding the proposed method. Hence, it shows that all three strategies (Gaussian smoothing and the proposed approach) greatly outperform EUCLID in moderate-to-intense noise, with Gaussian smoothing sometimes yielding the very best result.

\begin{table}[h!]
\caption{Noise impact on model performance metrics and relative improvements}\label{tab:noise_model_performance}
\begin{tabular}{c|c|c|c}
\hline
\\[-2.7ex]
\hline
\begin{tabular}{@{}c@{}}
	\textbf{Noise}
\\	\textbf{amplitude}
\end{tabular}
& \textbf{Model} &
\begin{tabular}{@{}c@{}}
	\textbf{Validation}
\\	\textbf{metric}
\end{tabular}
&
\begin{tabular}{@{}c@{}}
	\textbf{Relative}
\\	\textbf{improvement}
\end{tabular}\\
\hline
\\[-2.7ex]
\hline
0 & Reference & 1.2 & -- \\
\hline
0.01 & EUCLID & 10.0 & -- \\
& Gaussian-0.05 & 2.6 & +84.4\% \\
& Gaussian-0.1 & 3.5 & +73.6\% \\
& Proposed & 9.8 & +2.2\% \\
\hline
0.1 & EUCLID & 729 & --\\
& Gaussian-0.05 & 38 & +94.8\% \\
& Gaussian-0.1 & 36 & +95.1\% \\
& Proposed & 67 & +90.9\% \\
\hline
0.5 & EUCLID & 2584 & --\\
& Gaussian-0.05 & 130 & +94.9\% \\
& Gaussian-0.1 & 156 & +93.9\% \\
& Proposed & 186 & +93.2\% \\
\hline
\end{tabular}
\end{table}

\section{Discussion}\label{sec:DISCUSS}

The influence of measurement noise on the EUCLID framework, a local form of equilibrium for constitutive law identification, has been thoroughly examined in this study. Effective strategies to mitigate its impact have been proposed and evaluated through both a simplified toy model and a more complex numerical test case with realistic noise levels.

The results demonstrate that the proposed noise compensation method significantly outperforms traditional EUCLID when noise amplitude increases, although the results table reveals that Gaussian smoothing settings can surpass the proposed approach at medium and high noise levels. Specifically, the new table shows that at a noise of 0.1 pixels, Gaussian smoothing can yield up to a 95.1\% improvement, whereas the proposed method achieves 90.9\%. At a high noise of 0.5 pixels, Gaussian smoothing can surpass 94\%, compared to 93.2\% for the proposed method. One reason is that the training structure does not exhibit strong heterogeneity and localized states. Filtering the residual forces hence does not betray local features (\eg~cracks, plastic zone, singularities).
Despite not always producing the very lowest metric, the proposed method still offers large improvements over EUCLID across all tested noise scenarios, highlighting its robustness and ability to effectively counteract noise-induced biases.

The toy model results further corroborate these findings, illustrating how the internal loss term becomes highly sensitive to noise, leading to significant biases in stiffness identification. The proposed compensation strategies effectively mitigate these biases, ensuring accurate and reliable identification of material properties even under challenging noise conditions.

A comprehensive analysis of the noise impact on model performance is presented in Table~\ref{tab:noise_model_performance}, which summarizes the validation metrics and relative improvements across different noise amplitudes. At a low noise level of 0.01 pixels, the proposed method achieves a marginal improvement of 2.2\% over the standard EUCLID approach, demonstrating its initial effectiveness in mitigating noise-induced biases. As noise amplitude increases to 0.1 pixels, the proposed method exhibits a substantial relative improvement of 90.1\%, significantly outperforming both the EUCLID and Gaussian smoothing methods. This trend continues at an intense noise level of 0.5 pixels, where the proposed method achieves an impressive 93.2\% relative improvement, highlighting its robustness and superior capability in handling high-noise scenarios. These results affirm that the proposed noise compensation strategies not only reduce noise impact but also effectively compensate for it, ensuring that the learned constitutive behavior remains accurate despite significant measurement uncertainties.

These updated findings confirm that both approaches, Gaussian smoothing, and noise-bias compensation method, can reduce noise impact significantly. Choosing between them may depend on the ease of tuning smoothing parameters, the desired level of automation, and the particular geometry or local phenomena of interest.

Key factors contributing to the enhanced effectiveness of the proposed approach include:
\vspace{-.1\baselineskip}
\begin{description}
\item[Covariance-based weighting] By weighting the internal loss term with the inverse covariance matrix, the method ensures that nodes with higher uncertainty (\ie, those more affected by noise) are appropriately down-weighted in the loss function. This node-specific weighting allows for a more nuanced and accurate representation of the true material behavior, as it accounts for the varying levels of confidence in different parts of the displacement field. This integration eliminates the need for arbitrary weighting factors, which are often required in traditional approaches to balance different loss terms.
\item[Systematic bias compensation] Unlike smoothing, which passively dampens noise, the proposed method actively compensates for the noise-induced bias in the internal loss term. By introducing a bias compensation term, the method ensures that the loss function aligns more closely with the true equilibrium state of the system. This compensation is crucial for maintaining the integrity of the learned constitutive model, as it directly addresses the systematic errors introduced by noise.
\item[Mesh independence] The proposed approach inherently addresses mesh dependency issues by incorporating covariance-based weighting, which adjusts for the local noise characteristics. This feature ensures that the method effectiveness is not tied to the mesh density or distribution.
\end{description}

The noise introduced in this study reflects realistic measurement uncertainties typically encountered in experimental settings. By simulating noise directly at the displacement level and considering various amplitudes, the study ensures that the findings are relevant and applicable to practical scenarios where measurement noise is an inevitable factor. It is important to acknowledge that the effectiveness of the proposed noise correction methods may vary depending on the specific application and data representation employed. Factors such as the complexity of the material behavior, the geometry of the structure, and the fidelity of the finite element mesh can influence the extent to which noise impacts the learning process and how well the compensation strategies perform. For instance, in applications involving highly localized phenomena like crack propagation or in cases with intricate geometrical features, the current smoothing approach might require further adaptation to maintain its effectiveness.

\begingroup
A Bayesian version of EUCLID~\cite{joshi2022bayesian}, providing a probabilistic framework for uncertainty quantification, has recently been introduced in the literature and could complement the presented noise compensation strategies by offering a robust statistical framework for uncertainty propagation and quantification.
\endgroup

\goodbreak

\subsection{Limitations}

While the proposed methods substantially improve noise robustness, they also present certain limitations.

\paragraph{Limits on the Gaussian smoothing}
The current smoothing strategy is geometry-agnostic, which does not account for material state or structural features. This limitation could lead to inaccuracies in localized phenomena, such as cracks, where smoothing may blur critical discontinuities. Future work could explore adaptive smoothing techniques that adjust to local material states or employ specialized kernels tailored to specific structural features, such as cracks or stress concentrations~\cite{williams1957stress}. The selection of appropriate kernel types and sizes (rather than using a Gaussian one) presents a clear opportunity to increase accuracy and convergence speed. For example, in the case of localization (\eg, crack propagation), specific crack-related kernels, based on the Williams series, usually used in integrated DIC~\cite{rethore2009extended}, could be implemented.

\paragraph{Limits on the proposed approach}
The success of the covariance-based corrections depends on accurate noise characterization, which may be challenging in real-world applications with complex noise patterns. A major consideration is computational cost. Large 3D meshes with tens of thousands of nodes can make the many required iterative Monte Carlo sampling costly. In the current application, with only the diagonal values of the covariance being computed, the training time was increased approximately 4 times. In opposite, the filtered approach has no additional cost compare to standard EUCLID.

\section{Conclusion}\label{sec:CONCLUSION}

In this study, the impact of measurement noise on the Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) framework was thoroughly examined. Recognizing the inherent sensitivity of EUCLID to noise, which introduces significant biases in the identified material parameters, robust noise-bias correction strategies were proposed to enhance the framework reliability and accuracy.

The proposed methods, including covariance-based weighting and systematic bias compensation, were found to effectively mitigate the adverse effects of noise. By incorporating a covariance matrix that encapsulates the uncertainty associated with each nodal force measurement, nodes with higher uncertainty were appropriately down-weighted in the loss function. Additionally, the introduction of a bias compensation term derived from Monte Carlo simulations ensured that the loss function was aligned more closely with the true equilibrium state, thereby counteracting the systematic errors introduced by noise.

In practice, the tested approach simply consists of performing the inference on noisy signals (obtained from repeated acquisitions) to noise-induced forces and subtracting it from the original loss.

The efficacy of these strategies was demonstrated through both a simplified toy model and a more complex numerical test case with realistic noise levels. The results unequivocally showed that the proposed methods significantly outperformed traditional approaches such as Gaussian smoothing, especially at large noise amplitudes. In such a case, notably for a noise standard deviation of 0.5 pixel, relative improvements of up to 93.2\% in validation metrics were achieved, underscoring the robustness of the method in high-noise environments.
Furthermore, mesh dependency issues were inherently addressed by adjusting for local noise characteristics through covariance-based weighting, ensuring versatility across different finite element mesh configurations. This mesh independence, coupled with effective noise bias compensation, provided a more accurate and reliable framework for constitutive behavior identification.

The dependence on accurate noise characterization and the increased computational overhead due to covariance estimation and Monte Carlo simulations were identified as areas for future enhancement. Future work will focus on optimizing these computational aspects and exploring adaptive smoothing techniques tailored to specific structural features, such as crack propagation.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\setlength\bibitemsep{.1\baselineskip}

\printbibliography

\end{document}