# Mathematical and numerical assessment of Data-Driven Identification method applied to nonlinear elasticity

Nour Hachem, Adrien Leygue, Laurent Stainier

► To cite this version:

Nour Hachem, Adrien Leygue, Laurent Stainier. Mathematical and numerical assessment of Data-Driven Identification method applied to nonlinear elasticity. Computer Methods in Applied Mechanics and Engineering, 2025, 10.1016/j.cma.2025.118273. hal-04991913v3

HAL Id: hal-04991913

https://hal.science/hal-04991913v3

Submitted on 4 Aug 2025

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Mathematical and numerical assessment of Data-Driven Identification method applied to nonlinear elasticity

Nour Hachem, Adrien Leygue, Laurent Stainier

Nantes Université, École Centrale Nantes, CNRS, GeM, UMR 6183, 1 rue de la Noë, Nantes, F-44300, France

## Abstract

The Data-Driven Identification method estimates the mechanical response of material samples without postulating a predefined constitutive model. Leygue (2025) provided an initial mathematical analysis of this method, introducing a criterion to assess the well-posedness of the problem. This study builds on that work to deepen the mathematical understanding of DDI, with a focus on the selection of database size and algorithmic convergence. Three methods are presented to evaluate the well-posedness criterion, alongside the development of a new implementation process for the DDI algorithm based on one of these methods. The effectiveness of these approaches is demonstrated using synthetic data, and the convergence of the proposed algorithm is rigorously proven. This comprehensive theoretical and numerical analysis shows that the loss of uniqueness problem can be detected and corrected in this case. This work opens the door to further developments focusing on improving the method.

**Keywords:** Data-Driven Identification, Full-field measurement, Material identification, Alternating projection.

## 1. Introduction

Identifying the mechanical responses of materials is a fundamental challenge in solid mechanics. The constitutive equations close the boundary value problems of mechanics by establishing a relationship between kinematics (strain) and equilibrium (stress). Classically, the stress field is estimated using a constitutive model with parameters calibrated through identification
methods such as the finite element model updating (FEMU) [1], the virtual
fields method (VFM) [2, 3], the constitutive equation gap method (CEGM)
[4, 5], and the modified error in constitutive relation (MECR) [6, 7], among
others. For a global overview, the interested reader is referred to Avril et
al. [8] and Martins et al. [9]. These methods require the choice of a priori
constitutive models, and their parameters are adjusted by minimizing the dis-
tance between experimental measurements and numerical results. The initial
choice of the model strongly influences the results, potentially introducing
bias.

Other identification methods has been proposed, such as the unsupervised discovery of interpretable hyperelastic constitutive laws (EUCLID) [10] that proposed the use of sparse regression to identify relevant components from a library of (potential) models. In addition, machine learning based methods are being used to build surrogate models that learn material behavior directly from data [11, 12, 13]. However, these methods still rely on a models and require the adjustment of their parameters. Although adding more parameters to a model has been shown to improve the accuracy of identification [14], selecting the right parameters to include remains a nontrivial and open problem.

An alternative identification approach called Data-Driven Identification (DDI), proposed by Leygue et al. [15], aims to identify the material mechanical response of samples without postulating a constitutive model.

The main outputs of DDI are the statically admissible stress field on
the studied samples and a discrete representation of the constitutive relation
identified from full-field kinematic measurements [16, 17] and load cell mea-
surements. It relies on the availability of heterogeneous and rich data that
are clustered to build the strain-stress database. This database of admissible
strain-stress pairs samples the mechanical response of the material, and it is
inspired by the Data-Driven Computational Mechanics (DDCM) framework
[18]. Given the displacement and strain fields, DDI solves an inverse prob-
lem relying on the principle that clustered strain should yield similar stress,
which is achieved by minimizing the distance between the mechanical state
and its nearest representative in the database while ensuring equilibrium is
satisfied.

DDI has been tested and applied to synthetic and real data, including
studies on: small-strain elasticity [15], finite-strain elasticity [19, 20], vis-
coelasticity [21], and elastoplasticity [22, 23]. It has also been applied to
the case of bi-material samples [21]. Despite these promising applications,
a theoretical study of the method to validate the capacity of DDI to accu-
rately estimate and converge to the true mechanical stress is still missing.
Consequently, the parameters of the method have usually been selected em-
pirically, which may affect its robustness. Their influence was evaluated using
synthetic data as detailed in [24]. As these parameters are pushed to gain ac-
curacy, the DDI problem becomes ill-posed. Recently, Leygue [25] introduced
a new formulation that offers a criterion to determine the well-posedness of
the problem and an error estimate for the identified stress field under elastic
behavior.

In this paper, we investigate numerically and algorithmically the prelim-
inary work of Leygue [25]. First, the DDI method and the criterion for the
uniqueness of the obtained solution are recalled in Section 2. Second, we
present three methods to characterize the criterion, supported by a theoreti-
cal evaluation of their efficiency in Section 3. Then, in Section 4, we introduce
a new DDI algorithm inspired by the last method along with a proof of its
convergence. These methods are systematically tested on synthetic data,
with a detailed numerical analysis and discussion of their performance.

## 2. Data-Driven Identification

This section examines the recent DDI formulation presented in [25] and revisits the uniqueness criterion established by this formulation.
Throughout this work, all fields are assumed to be sufficiently regular. The equations are expressed after discretization utilizing a finite element mesh comprising *m* quadrature points and *n* nodes. The strain is described using a linearized kinematic measure, and the equilibrium equation is expressed in the reference configuration.

Index notations

The following index notation is used:

• $\cdot_e$ denotes a quantity evaluated at the integration point $e$, with $e \in \{1, \dots, m\}$.

• The index $a$ is also used to denote an integration point, occasionally re-
placing $e$ in expressions involving pairwise interactions between quadra-
ture points (e.g., double summations).

• $\cdot_k$ denote a quantity at node $k$, with $k \in \{1, \dots, n\}$.
• ·_i refers to a quantity associated with material state i, with i ∈ {1, ..., N_D}.

• Primed indices (e.g., e', a', k', i') are used as dummy summation index.

## 2.1. Formulation

The core idea of the method is to estimate the stress field at the quadra- ture points $\sigma_e$, with $e \in \{1, ..., m\}$, without postulating a constitutive model a priori. This is achieved by minimizing the distance between the current me- chanical state ($\varepsilon_e, \sigma_e$) and its nearest representative in a material database $\mathcal{D}$, under the equilibrium constraints.

The method relies on the following input data:

• the discrete divergence operator $\mathbf{D}_{ke}$, $\forall k \in \{1, ..., n\}$, $\forall e \in \{1, ..., m\}$,

• the external force vector $\mathbf{f}_k$, $\forall k \in \{1, ..., n\}$,

• the strain at all quadrature points, $\varepsilon_e$, $\forall e \in \{1, ..., m\}$.

The static admissibility is therefore expressed as:

$$ \sum_{e'=1}^{m} \mathbf{D}_{ke'} \cdot \boldsymbol{\sigma}_{e'} = \mathbf{f}_k \quad \forall k. \qquad (1) $$

The material database $\mathcal{D}$ samples the material behavior by a set of $N_D$ strain-stress pairs called *material states*,

$$ \mathcal{D} = \{(\mathcal{E}_i, S_i)\}_{i=1,...,N_D}. \qquad (2) $$

Each *mechanical state* ($\varepsilon_e, \sigma_e$) is connected to a material state ($\mathcal{E}_i, S_i$) through a $m \times N_D$ binary matrix pairing $P$. The matrix component $P_{ei}$ takes the value 1 if the mechanical state $e$ is paired with the material state $i$, or equivalently, if it belongs to cluster $i$, and 0 otherwise:

$$ P_{ei} \in \{0, 1\} \quad \text{for } 1 \le e \le m \text{ and } 1 \le i \le N_D, \qquad (3) $$

$$ \sum_{i'=1}^{N_D} P_{ei'} = 1 \quad \forall e. \qquad (4) $$
DDI is a minimization problem that estimates the mechanical stress field $\boldsymbol{\sigma}$, the material database $\mathcal{D}$, and the matrix $P$. The loss function $\mathcal{F}$ of the DDI problem is given by:

$$ \begin{aligned} \mathcal{F}(\boldsymbol{\sigma}_e, \boldsymbol{\varepsilon}_i, \boldsymbol{S}_i, P_{ei}) = & \frac{1}{2} \sum_{e'} w_{e'} \left( \boldsymbol{\varepsilon}_{e'} - \sum_{i'} P_{ei'} \boldsymbol{\varepsilon}_{i'} \right) : \mathbb{C} : \left( \boldsymbol{\varepsilon}_{e'} - \sum_{i'} P_{ei'} \boldsymbol{\varepsilon}_{i'} \right) \\ & + \frac{1}{2} \sum_{e'} w_{e'} \left( \boldsymbol{\sigma}_{e'} - \sum_{i'} P_{ei'} \boldsymbol{S}_{i'} \right) : \mathbb{C}^{-1} : \left( \boldsymbol{\sigma}_{e'} - \sum_{i'} P_{ei'} \boldsymbol{S}_{i'} \right) \\ = & \frac{1}{2} \sum_{e'} w_{e'} \left\| (\boldsymbol{\varepsilon}_{e'}, \boldsymbol{\sigma}_{e'}) - \sum_{i'} P_{ei'} (\boldsymbol{\varepsilon}_{i'}, \boldsymbol{S}_{i'}) \right\|_{\mathbb{C}}^2 . \end{aligned} \quad (5) $$

The functional $\mathcal{F}$ represents the distance between the mechanical states and the material states, measured under the norm $\|\cdot\|_{\mathbb{C}}$, where $\mathbb{C}$ is a fourth-order tensor and a parameter of the method. The norm is defined as:

$$ \|(A, B)\|_{\mathbb{C}}^2 = A : \mathbb{C} : A + B : \mathbb{C}^{-1} : B, $$

with $A$ and $B$ are two second-order tensor (see [18, 15] for details regarding the norm).

The problem is formulated as follows:

$$ (\boldsymbol{\sigma}_e, \boldsymbol{\varepsilon}_i, \boldsymbol{S}_i, P_{ei}) = \arg \min \mathcal{F}(\boldsymbol{\sigma}'_e, \boldsymbol{\varepsilon}'_i, \boldsymbol{S}'_i, P'_{ei}), \quad (6) $$

under the constraint of mechanical equilibrium given in Equation (1).

The only parameters of this method are: the scalar $N_{\mathcal{D}}$ referring to the number of clusters, the tensor $\mathbb{C}$ related to the norm, and the scalar weights $w_e$.

Recently, DDI has been formulated by a new loss function that no longer explicitly includes the material states, but implicitly through the matrix $\tilde{H}$.

The recent formulation [25] is given by:

$$ (\boldsymbol{\sigma}_e, P_{ei}) = \arg \min \left\{ \frac{1}{2} \sum_{e', a'} \tilde{H}_{e'a'} \boldsymbol{\varepsilon}_{e'} : \mathbb{C} : \boldsymbol{\varepsilon}_{a'} + \frac{1}{2} \sum_{e', a'} \tilde{H}_{e'a'} \boldsymbol{\sigma}_{e'} : \mathbb{C}^{-1} : \boldsymbol{\sigma}_{a'} \right\}, \quad (7) $$

$$ \text{subject to} \quad \sum_{e'=1}^{m} \mathbf{D}_{ke'} \cdot \boldsymbol{\sigma}_{e'} = \mathbf{f}_k \quad \forall k. $$
with the matrix $\tilde{H}$ defined as:

$$ \tilde{H}_{ea} = \sum_{e'} (I_{ee'} - \tilde{M}_{ee'}) w_{e'} (I_{e'a} - \tilde{M}_{e'a}) \quad \forall e, a. \qquad (8) $$

where

$$ \tilde{M}_{ea} = \sum_{i} \frac{1}{W_i} P_{ei} P_{ai} w_a \quad \forall e, a, \qquad (9) $$

and

$$ W_i = \sum_{e'} P_{e'i} w_{e'} \quad \forall i. \qquad (10) $$

The matrix $\tilde{H}$ in (7) depends implicitly on the matrix $P$, i.e. $\tilde{H} := \tilde{H}(P)$. As discussed in [25], the matrix $\tilde{M}$ computes the mean of each cluster and $\tilde{M}$ is a projection onto nullspace of $\tilde{H}$. The columns of $P$ form a basis of ker $H$ and $P$ is thus orthogonal by construction. For a more detailed discussion on the properties of $\tilde{M}$ and $\tilde{H}$, see [25].

The unknown material states of the DDI problem (6) can be found by this new formulation. They are simply the $w_e$-weighted average of the mechanical states over each cluster. Mathematically, this can be expressed as follows:

$$ \boldsymbol{\epsilon}_i = \frac{1}{W_i} \sum_{e'} P_{e'i} w_{e'} \boldsymbol{\epsilon}_{e'} \quad \forall i, \qquad (11) $$

and

$$ \boldsymbol{\sigma}_i = \frac{1}{W_i} \sum_{e'} P_{e'i} w_{e'} \boldsymbol{\sigma}_{e'} \quad \forall i. \qquad (12) $$

For algebraic convenience, the problem will be expressed in matrix notation instead of index notation. The problem (7) can be rewritten using Mandel notation [26] as follows:

$$ (\sigma, P) = \arg \min \left\{ \frac{1}{2} \varepsilon^{\top} H_{\varepsilon} \varepsilon + \frac{1}{2} \sigma^{\top} H_{\sigma} \sigma \right\} \quad \text{with} \quad D\sigma = F. \qquad (13) $$

where,

* $\varepsilon = [\varepsilon_e^M]$ denotes the collection vector corresponding to the second-order tensors $\varepsilon$, in Mandel notation,

* $\sigma = [\sigma_e^M]$ denotes the collection vector corresponding to the second-order tensors $\sigma$, in Mandel notation,
* $M = \tilde{M} \otimes I$ represents the Kronecker product of the matrix $\tilde{M}$ and identity matrix $I$,

* $H_{\varepsilon} = \tilde{H} \otimes C$ represents the Kronecker product of the matrix $\tilde{H}$ and $C$,

* $H_{\sigma} = \tilde{H} \otimes C^{-1}$ represents the Kronecker product of the matrix $\tilde{H}$ and $C^{-1}$,

* $C$ refers to the matrix representation of the fourth-order tensor $C$, in Mandel notation,

* $D$ and $F$ denote the discrete divergence operator $D$ and the external force vector $f$, in Mandel notation, respectively.

* $\cdot^{\top}$ denotes the conjugate transpose.

## 2.2. Uniqueness of the DDI solution

In [25], the author studies the well-posedness of the problem (13), which we will revisit in this section. Given that problem (13) is equivalent to the problem (6), it follows that they both exhibit the same properties of well-posedness.

**Theorem 2.1.** [25] The solution to the problem (13) is unique if and only if there does not exist a non-zero vector $\tilde{\sigma}$ such that $H_{\sigma}\tilde{\sigma} = 0$ and $D\tilde{\sigma} = 0$, that is, if and only if $\ker D \cap \ker H_{\sigma} = \{0\}$.

**Remark.** To relax the notation, we denote $H_{\sigma}$ as $H$, unless otherwise specified.

We consider a given clustering of deformations, defined through $P_{ei}$, i.e. the matrices $H_{\sigma}$ and $H_{\varepsilon}$ in Equation (13) are fixed. Then, the mechanical stresses $\sigma_{\varepsilon}$ are the solution of the following problem:

$$ \sigma = \arg \min \left\{ \frac{1}{2} \|\sigma\|_{H}^{2} \right\} \quad \text{with} \quad D\sigma = F. \qquad (14) $$

Here, $\| \cdot \|_{H}$ denotes a semi-norm, as $H$ is positive semi-definite [25], and defined as $\|\sigma\|_{H}^{2} := \sigma^{\top} H \sigma$.

*Proof.* We will prove only the implication of interest to us. (=$\Leftarrow$) Suppose that $\ker D \cap \ker H = \{0\}$. Our goal is to prove that the solution is unique under this assumption. Let $\sigma_{1}$ and $\sigma_{2}$ be two different solutions to the problem. Thus, we have:
1. Since $D\sigma_1 = F$ and $D\sigma_2 = F$, it follows that $\sigma_1 - \sigma_2 \in \ker D$.

2. As both $\sigma_1$ and $\sigma_2$ minimize the same quadratic functional, we have $\frac{1}{2}\sigma_1^\top H\sigma_1 = \frac{1}{2}\sigma_2^\top H\sigma_2$. This implies that $(\sigma_1 - \sigma_2)^\top H(\sigma_1 - \sigma_2) = 0$. Since $H$ is positive semi-definite, it follows that $\sigma_1 - \sigma_2 \in \ker H$.

From points 1 and 2, we can deduce that $\sigma_1 - \sigma_2 \in \ker D \cap \ker H = \{0\}$. As a result, we have $\sigma_1 = \sigma_2$. $\square$

### 3. How to characterize the intersection of the two kernels in practice?

This section presents three methods for characterizing the intersection of two kernels: (3.1) explicit calculation of the basis, (3.2) examination of the angles between the subspaces, and (3.3) alternating projection onto the two kernels.

The proofs of the theorems will not be detailed in the following subsections, with references provided for interested readers.

#### 3.1. Compute an intersection base explicitly

**Proposition 3.1.** [27] The intersection of the kernel of $H$ ($\ker H$) and the kernel of $D$ ($\ker D$) is the kernel of the larger matrix $G$, where $G = \begin{bmatrix} H \\ D \end{bmatrix}$.

In other words,

$$ \ker D \cap \ker H = \ker G. \qquad (15) $$

*Proof.* $G\sigma = 0$ implies both $H\sigma = 0$ and $D\sigma = 0$. So $\sigma$ is in the intersection of $\ker D$ and $\ker H$. $\square$

The matrices $H$ and $D$ are sparse, therefore, the matrix $G$ is also sparse. Numerically, the null space of a matrix is computed using singular value decomposition (SVD), which involves non-sparse linear algebra calculations. The computational complexity of SVD is $O(ab^2)$ where $a$ and $b$ are the dimensions of $G$, with $a \ge b$. For example, in the bidimensional case, we have $a = 3m + 2n$ and $b = 3m$, showing that this complexity is independent of the database size $N_D$. Nevertheless, it provides an exact basis for the intersection.
### 3.2. Angles between subspaces

Let $A \in \mathbb{R}^{d \times p}$ and $B \in \mathbb{R}^{d \times q}$ be two matrices containing basis of two given subspaces $A, B \subset \mathbb{R}^d$, respectively, and assume that $p = \dim A \ge \dim B = q \ge 1$.

**Definition 3.1.** [28] The principal angles $\theta_r \in [0, \pi/2]$ between $A$ and $B$ are defined for $r = 1, 2, \dots, q$ by:

$$ \cos \theta_r = \max_{v \in A} \max_{u \in B} u^\top v, \quad \|u\|_2 = 1, \|v\|_2 = 1, \qquad (16) $$

subject to

$$ u_j^\top u = 0, \quad v_j^\top v = 0, \quad j = 1, 2, \dots, r-1. \qquad (17) $$

The vectors $(u_1, u_2, \dots, u_q)$ and $(v_1, v_2, \dots, v_q)$ are called the principal vectors associated to $\{\theta_r\}_{r=1}^q$.

The following theorem illustrates the relation between singular values and the principal angles.

**Theorem 3.2.** [28] Let $Q_A$ and $Q_B$ two matrices such that their columns form unitary basis of $A$ and $B$ respectively. Let $N = Q_A^\top Q_B$ and let the singular values decomposition SVD of this $p \times q$ matrix be

$$ N = Y \Lambda Z^\top, \quad \Lambda = \operatorname{diag}(\lambda_1, \lambda_2, \dots, \lambda_q). \qquad (18) $$

Then, the principal angles and the principal vectors associated are given by,

$$ \cos(\theta_r) = \lambda_r(N), \quad U = Q_A Y, \quad V = Q_B Z, \qquad (19) $$

where $\lambda_r$ are the singular values of the matrix $N$.

*Proof.* The proof of this theorem is provided in Section 3 of [28]. □

If a principal angle $\theta_r$ equals zero, or equivalently, a singular value equals one, the solution is non-unique. The number of principal angles equal to zero corresponds to the dimension of the kernel intersection.

The theorem outlined above can also be used to compute an orthonormal basis for the intersection, as given by the following theorem.
**Theorem 3.3.** [29] Let {$\cos(\theta_r)}_{r=1}^q$ and {$u_r, v_r}_{r=1}^q$ be defined by Theorem 3.2. If the index $s$ is defined by $1 = \cos(\theta_1) = \cdots = \cos(\theta_s) > \cos(\theta_{s+1})$, then

$$ \operatorname{span}\{u_1, \dots, u_s\} = \operatorname{span}\{v_1, \dots, v_s\} = \mathcal{A} \cap \mathcal{B}. \qquad (20) $$

*Proof.* The proof is based on the observation that when $\cos(\theta_r) = 1$, it implies that $u_r = v_r$. $\square$

To practically determine the intersection dimension $s$, it is required to define a tolerance for when a computed singular value can be considered equal to 1.

We apply these theorems by setting $\mathcal{A} = \ker D$ and $\mathcal{B} = \ker H$. The procedure follows the following steps:

0. According to [25], the null space of $H$ is spanned by the columns of $P$,

1. Compute the kernel of $D$, numerically, using SVD decomposition,

2. Compute unit basis of $\ker D$ and $\ker H$, $Q_A$ and $Q_B$, respectively,

3. Perform SVD on matrix $N$, the singular values represent the principal angles between the subspaces according to Theorem 3.2. And from the principal vectors we can compute an intersection basis according to Theorem 3.3.

This method provides insight into the conditioning of our problem, as the presence of very small angles near zero indicates poor conditioning.

Moreover, this approach is computationally less expensive than the previous one. In the first step, we only need to compute the null space of $D$, which is less costly than computing the null space of matrix $G$. Although this method allows us to determine the basis for the intersection, it is based on a chosen tolerance. However, it still involves processing full matrices, which does not meet the sparsity criterion. To address this need, we propose an algorithm based on alternating projection onto the two subspaces and preserve the sparse matrix property of the matrices $H$ and $D$.

### 3.3. Sampling the intersection through alternating projection

The alternating projection is a simple algorithm to calculate a point at the intersection of two sets by using a sequence of projections on these sets. Suppose $\mathcal{X}$ and $\mathcal{Y}$ are two closed convex sets in $\mathbb{R}^d$, and let $P_{\mathcal{X}}$ and $P_{\mathcal{Y}}$ denote projection onto $\mathcal{X}$ and $\mathcal{Y}$, respectively. The algorithm begins with any point $x_0 \in \mathbb{R}^d$ and then subsequently alternates projections onto $\mathcal{X}$ and $\mathcal{Y}$,
$$y_j = P_Y(x_j), \quad x_{j+1} = P_X(y_j), \quad \text{for } j \ge 0. \qquad (21)$$

**Theorem 3.4.** [30]

1. If $X \cap Y \neq \emptyset$, then the two sequences $x_j$ and $y_j$ converge to $x^* \in X \cap Y$, meaning the alternating projection finds a point in the intersection.

2. If $X \cap Y = \emptyset$, then $x_j \to x^* \in X$ and $y_j \to y^* \in Y$, where $\|x^* - y^*\|_2 = \text{dist}(X, Y)$.

*Proof.* The convergence proof is given in Section 2 of [30]. □

The null space of a matrix is convex by definition, as it forms a subspace where all linear combinations are invariant. Moreover, this subspace is closed because the kernel is the inverse image of a closed set, the singleton zero $\{0\}$. As both kernels are subspaces, their intersection necessarily contains the zero vector. Therefore, our analysis focuses on the first part of the theorem, where we define $X = \text{ker } H$ and $Y = \text{ker } D$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Alternating projection onto $\text{ker } D$ and $\text{ker } H$.</figcaption></figure>

The algorithm is given by (see graphical illustration in Fig. 1):

$$\sigma_0 \xrightarrow{P_H} \sigma_1^H \xrightarrow{P_D} \sigma_1^D \xrightarrow{P_H} \sigma_2^H \xrightarrow{P_D} \sigma_2^D \xrightarrow{P_H} \dots \to \sigma^*, \qquad (22)$$

where the projections, for $j \ge 1$, are defined as follows:

* The projection $P_H$ onto $\text{ker } H$ is achieved by multiplying $\sigma$ by $M$ [25]:

$$P_H(\sigma_j) := \sigma_j^H = M\sigma_{j-1}, \qquad (23)$$
• The projection $P_D$ onto $ker D$ is computed by solving the following minimization problem:

$$ P_D(\sigma_j^H) := \sigma_j^D = \arg\min_{\sigma_j} \{(\sigma_j - \sigma_j^H)^\top w(\sigma_j - \sigma_j^H)\} \quad \text{with} \quad D\sigma_j = 0. \tag{24} $$

with $w$ is the weight.

To understand how this algorithm works, consider that the initial point is composed of $\sigma_0 = \sigma_0^I + \sigma_0^\perp$, where $\sigma_0^I$ lies in the intersection $ker D \cap ker H$ and $\sigma_0^\perp$ represents the remainder. The iterative process is defined by:

$$ (P_D P_H)^j(\sigma_0) = (P_D P_H)^j(\sigma_0^I) + (P_D P_H)^j(\sigma_0^\perp), \tag{25} $$

as $j \to \infty$, the convergence results hold:

• the sequence $(P_D P_H)^j \sigma_0$ converges to an intersection, denoted $\sigma^*$,

• since $\sigma_0^I$ is invariant under these projections, the sequence $(P_D P_H)^j \sigma_0^I$ remains unchanged and converges to $\sigma_0^I$,

• the term $(P_D P_H)^j \sigma_0^\perp$, representing the component orthogonal to the intersection, decreases over iterations and tends to zero.

In the case where $\sigma_0^I = 0$, or equivalently, $\sigma_0 \in (\ker H \cap \ker D)^\perp$, the algorithm can converge to $\sigma^* = 0$, even though the solution is not necessarily unique. Given that $\sigma_0$ in (22) is chosen randomly, it is highly improbable that $\sigma_0^I = 0$. This condition is similar to the requirement in the power method: the initial vector should not be orthogonal to the dominant eigenvector to ensure proper convergence. Bindel [31] and Demmel [32] discuss that the probability of randomly selecting an initial vector that satisfies this condition is zero. Based on this, we can conclude:

• if $\sigma^* \neq 0$, the solution is non-unique,

• if $\sigma^* = 0$, the solution is unique, with very high probability—assuming that $\sigma_0$ is chosen randomly.

In addition, this method allows us to compute a basis by performing multiple initializations and then using a method to check the linear independence of the resulting vectors, such as Gram-Schmidt.

The convergence rate of the alternating projection method is related to the angle between the two subspaces through the following theorem:
**Theorem 3.5.** [33] For all $j \in \mathbb{N}$,

$$ \|(P_X P_Y)^j - P_{X \cap Y}\| = \cos^{2j-1}(\theta_{s+1}), \qquad (26) $$

where $\theta_{s+1}$ is the smallest positive principal angle.

*Proof.* The proof of the rate of convergence theorem is in Chapter 9, Theorem 9.31, of [33]. □

Based on this theorem, the convergence is slow if the smallest positive angle $\theta_{s+1}$ is small.

This algorithm is computationally less expensive as it takes advantage of the sparsity of the matrix. As shown in Equation (24), there is no need to compute $ker D$ explicitly, thus, making the algorithm faster and easier to implement.

## 4. Solution

Following [15, 25], the DDI problem is solved using an alternating minimization approach regarding the discrete variable $P$ and the continuous variable $\sigma$.

The first minimization over $P$ with a fixed continuous variable $\sigma$ is given by applying $w_e$-weighted k-means clustering in the $\mathbb{C}$-norm on the mechanical states. For the second minimization with a fixed $P$, the problem (13) is simplified to:

$$ \sigma = \arg \min \left\{ \frac{1}{2} \|\sigma\|_H^2 \right\} \quad \text{with} \quad D\sigma = F. \qquad (27) $$

In the following section, we focus on this iteration where $P$ is fixed, and the goal is to solve the minimization problem given in Equation (27).

### 4.1. Limitations

The Lagrangian for the minimization problem (27) is given by:

$$ \begin{bmatrix} H & D^{\top} \\ D & 0 \end{bmatrix} \begin{bmatrix} \sigma \\ \eta \end{bmatrix} = \begin{bmatrix} 0 \\ F \end{bmatrix}, \qquad (28) $$

where $\eta$ denotes the Lagrange multiplier associated, and we define $T := \begin{bmatrix} H & D^{\top} \\ D & 0 \end{bmatrix}$.
<table><caption>Table 1: Comparison of connectivity patterns using the Reverse Cuthill-McKee permutation, with connectivity based on clustering (top row) and mesh (bottom row). The mesh used is the one shown in Figure 2, and for the clustering we take $N_D = 300$.</caption><thead><tr><th>Numbering</th><th><i>H</i></th><th><i>D</i><sup>T</sup><i>D</i></th></tr></thead><tbody><tr><td>Cluster-based</td><td></td><td></td></tr><tr><td>Mesh-based</td><td></td><td></td></tr></tbody></table>
The matrix $H$ encodes the connectivity between mechanical states based on clustering, while the matrix $D^\top D$ encodes the connectivity within the mesh. By applying the reverse Cuthill-McKee algorithm [29] to either of these matrices we can compute a renumbering of the degrees of freedom, either clustering based or mesh based. In table 1, we apply these renumbering to both matrices and observe that their sparsity patterns are incompatible and that the difference in the non-zero structures might result in significant fill-in during standard linear algebra operations.

In direct methods such as LU decomposition, this fill-in appears and increases memory usage and computational cost by introducing additional non-zero elements into the matrix during the factorization process. For example, in the example presented, the matrix $T$ has only 0.19% non-zero elements. However, after LU decomposition, the percentage of non-zero elements in $L+U$ rises to 13.14%.

To address these limitations, we propose an iterative approach based on the alternating projection described in Section 3.3.

## 4.2. New algorithm based on alternating projection

We recall from [25] that $M$ is a projection onto $ker H$, obtained by computing the mean of each cluster.

The solution of Equation (27) is obtained by alternately projecting $\sigma$ onto $ker H$ and $E := \{\sigma \text{ s.t. } D\sigma = F\}$. The algorithm for solving (27) is given by the following steps:

1. Initialization of $\sigma_0$,

2. Projection onto $ker H$:

$$ \sigma_j^H = M\sigma_{j-1} \quad \text{for } j \ge 1, \qquad (29) $$

3. Projection onto $E$ [18]:

$$ \sigma_j = \sigma_j^H + D^\top \eta \quad \text{with} \quad DD^\top \eta = F - D\sigma_j^H \quad \text{for } j \ge 1. \qquad (30) $$

*Proof of convergence.* Since both $ker H$ and $E$ are closed convex sets, we can safely apply Theorem 3.4, ensuring the convergence of the algorithm. According to this theorem, the algorithm will converge either to a distance or a point within their intersection. In both cases, the algorithm converges to the minimum of our loss function. $\square$
**Algorithm 1** DDI algorithm based on alternating projection

**Require:** $\varepsilon$, $D$ and $F$

1: Choose DDI parameters: $w_e$, $N_D$ and $\mathbb{C}$

2: Choose algorithm parameters: IterMax, tolerance

3: Initialize $\sigma_0$

4: Compute $P$ by performing a $w_e$-weighted k-means clustering of $\varepsilon$ in the $\mathbb{C}$-norm

5: **for** $j = 1$ to IterMax **do**

6: &nbsp;&nbsp;&nbsp;&nbsp; $\sigma_j^H = M\sigma_{j-1}$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Proj. onto ker $H$

7: &nbsp;&nbsp;&nbsp;&nbsp; Solve for $\eta$: $DD^\top\eta = F - D\sigma_j^H$

8: &nbsp;&nbsp;&nbsp;&nbsp; Update: $\sigma_j = \sigma_j^H + D^\top\eta$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ Proj. onto $E$

9: &nbsp;&nbsp;&nbsp;&nbsp; **if** $\|\sigma_{j-1} - \sigma_j\| \le \text{tolerance}$ **then**

10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **break** &nbsp;&nbsp;&nbsp;&nbsp; ▷ Stagnation reached

11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end if**

12: **end for**

**Ensure:** $\sigma$, $P$

We recommend initializing the algorithm with $\sigma_0 = 0$. In this case, the first projection $\sigma_1^H$ yields to:

$$ \sigma_1^H = M\sigma_0 = 0, $$

and the projection onto the set $E$ then gives:

$$ \sigma_1 = D^\top(DD^\top)^{-1}F = D^\dagger F, $$

where $D^\dagger$ is the pseudo-inverse of matrix $D$.

This means that the first iterate $\sigma_1$ is the pseudo-inverse solution (the minimum norm solution) of the equilibrium equation $D\sigma = F$. This provides a meaningful and stable starting point for the iterative projection process.

This algorithm can still be applied even if the matrix $T$ is singular, unlike direct methods, and it will converge to one of the possible solutions.

Similarly to the third method, the rate of convergence depends on the smallest positive principal angle according to Theorem 3.5. In the same way, the conditioning of the problem is related to the presence of very small angles close to zero, as discussed in Section 3.2.
## 5. Application and discussion

In this section, we apply the DDI method to synthetic data generated from a finite element simulation. The 2D finite element mesh includes 1340 nodes and 2416 quadrature points, as shown in Figure 2. The model is subjected to a 10% strain in the vertical direction, assuming plane stress. The material is assumed incompressible and is characterized by the following strain-stress relation:

$$ \hat{\sigma}(\varepsilon) = -p\mathbf{I} + 2G(1 + \alpha \mathbf{dev}(\varepsilon) : \mathbf{dev}(\varepsilon))\mathbf{dev}(\varepsilon), \quad (31) $$

where $\mathbf{dev}(\varepsilon)$ represents the deviatoric part of $\varepsilon$, and $p$ is determined by the plane stress condition. The parameter $\alpha$ is chosen to be $\alpha = 3$, the parameter $C$ is selected as:

$$ C = \left. \frac{\partial \hat{\sigma}}{\partial \varepsilon} \right|_{\varepsilon=0}, $$

and $w_e$ are set to the mesh quadrature weights. For all the results presented in this section, the DDI is carried out for different values of $N_D$ ranging from 10 to 630, more significantly between on average 4 to 242 elements by cluster. For each $N_D$, the pairing matrix $P_{ei}$ is computed from a $w_e$-weighted k-means clustering performed on the strain field only.

The stress field generated by the model $\hat{\sigma}$ is denoted as $\sigma_e^{FE} = \hat{\sigma}(\varepsilon)$. The mechanical stress resulting from the DDI will be denoted as $\sigma_e^{DDI}$. The DDI solution presented in this section is computed using the new algorithm described in Section 4.2. For each application of the DDI, the error is calculated as follows:

$$ \text{error} = \frac{1}{m} \sum_{e} \| \sigma_e^{DDI} - \sigma_e^{FE} \|_2^2. \quad (32) $$

Figure 3 illustrates the loss of well-posedness of DDI problem. The upper Figure 3a shows the cost function of DDI across different database size $N_D$ for two stress initialization $\sigma_0$ schemes: random initialization following a normal distribution and zero initialization. The results indicate that the cost function reaches the same minimum value regardless of the initialization, confirming that convergence is achieved. The lower Figure 3b shows the error as a function for the same $N_D$ and same two initialization schemes. The error decreases as $N_D$ increases, indicating that DDI converges toward the true stress field. However, above a certain $N_D$, the DDI converges to different results depending on the initialization. This highlights the loss of
<figure><img src="image_3.png" /><figcaption>Figure 2: Finite element mesh used to generate the synthetic data. The black line indicates the deformation contour.</figcaption></figure>

well-posedness of the problem: for this $N_D$, it is possible to find a non-trivial stress field that is constant within each cluster (kernel of $H$) and self-balanced (kernel of $D$). This stress field corresponds to the difference between the two computed solutions. This can be demonstrated through one of the methods described earlier in Section 3.

In Figure 4, we apply the first method, where the basis is explicitly computed. Figure 4 shows the dimension of the computed basis as a function of $N_D$. Instances with a zero basis dimension, indicating unique solutions, are marked in green circles, while red diamonds represent cases with a non-zero basis dimension, corresponding to a loss of uniqueness. Notably, the loss of uniqueness observed in Figure 3b aligns with these red points, as expected. Although the first method is used here, the same result would be obtained using any of the three methods detailed earlier, as supported by mathematical proof.

Table 2 compares the computational costs of the three methods outlined in Section 3 for determining the uniqueness of solutions and computing the intersection basis.

First method explicitly computes the intersection basis, requiring 5763 seconds for both calculations. It relies on dense matrix operations, which de-
<figure><img src="image_4.png" /><figcaption>(a) Cost function for different database sizes $N_D$ with two initializations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Comparison of the error using the finite element solution and DDI results for the same two initializations.</figcaption></figure>

<caption>Figure 3: Illustration of the loss of uniqueness.</caption>
<figure><img src="image_6.png" /><figcaption>Figure 4: Dimension of the intersection basis for various values of $N_D$.</figcaption></figure>

<table><caption>Table 2: Comparison of computational costs for verifying the well-posedness and computing the full intersection basis for different database sizes $N_D$. The CPU times were measured using the *tic/toc* functions in MATLAB® on an Intel Core i7-12800H processor (2.40 GHz) with 32 GB of RAM.</caption><thead><tr><th>Time (s)</th><th>Method 1</th><th>Method 2</th><th>Method 3</th></tr></thead><tbody><tr><td>Uniqueness</td><td>5763.00</td><td>456.51</td><td>41.85</td></tr><tr><td>Basis</td><td>5763.00</td><td>456.51</td><td>294.77</td></tr></tbody></table>
mand substantial time and memory. Second method reduces the computa-
tion time to 456.51 seconds for the two tasks using principal angles, with
a tolerance of 10⁻⁶ to identify zero angles. It also relies on dense matrices
and lacks scalability for large-scale problems. Third method uses an itera-
tive approach for computing a vector in the intersection of ker *D* and ker *H*
using alternating projection. If the vector is zero, the solution is unique.
Otherwise, there are multiple solutions. The time for uniqueness verification
corresponds to the computation of this single vector. To compute the full
intersection basis, the method is run repeatedly with different initializations,
and the Gram-Schmidt process is applied to ensure linear independence of
the resulting vectors. The total time to compute a basis of dimension *s* is the
time required to compute one vector (i.e. a single algorithm run), multiplied
by *s* + 1, for all *s* ≥ 0. By using sparse matrix structures, method three
significantly reduces computational time to 41.85 seconds for uniqueness ver-
ification and 294.77 seconds for basis computation.

All three methods produced the same results, although, for large $N_D$, slight
differences in the basis were observed with method 3. These differences come
from the choice of tolerance, the iteration limit, and the presence of very
small positive principal angles, which are nearly zero. Using Theorem 3.5,
we know:

$$
j = \frac{\ln(\text{tol})}{2 \ln(\cos(\theta_{s+1}))} + \frac{1}{2}, \quad (33)
$$

where tol = ||(PₓPᵧ)ⱼ - Pₓ∩ᵧ||. The convergence becomes slower as the
smallest positive angle decreases, requiring more iterations to achieve conver-
gence. For example, when the smallest positive angle is nearly zero θs+1 → 0,
the above equation shows that we need infinite iteration j → ∞ to achieve
convergence. In this analysis, we applied a tolerance of 10⁻⁵ for the Gram-
Schmidt process. For the alternating projection method, the stopping criteria
were either a maximum of 5000 iterations or stagnation indicating conver-
gence.

The last method (alternating projection) is the most effective and efficient
time and memory usage approach. For example, in larger meshes such as in
[20] and [34], methods 1 and 2 cannot be applied due to their high memory
requirements. The same discussion on the efficiency of the third method also
applies to the new algorithm given in Section 4.2.

Leygue et al. (2018) [15] conducted a preliminary study on noise influence
and showed that DDI is largely robust to displacement field noise due to the
averaging effect inherent in clustering measurement. However, it is important
<figure><img src="image_7.png" /><figcaption>Figure 5: Probability estimator of uniqueness (true/false) for different $N_D$ values across 600 different k-means initializations.</figcaption></figure>

to consider the uncertainty that arises from k-means clustering initialization. K-means is usually initialized with k-means++ [35], which provides a good heuristic for initialization. Since k-means++ relies on a random number generator to sample the data distribution, it introduces a natural variability in DDI clustering.

Figure 5 illustrates the influence of k-means initializations on the frequency of obtaining a unique solution for different numbers of clusters $N_D$. We use the last method, where the algorithm output is a vector in the intersection basis. A boolean value is defined to indicate "true" if uniqueness is achieved, equivalently when $\sigma^* = 0$, and "false" otherwise. This process is repeated across 600 different clustering for each $N_D$ to calculate the frequency of obtaining a unique solution for various $N_D$. The results reveal three distinct zones. In the low $N_D$ zone, the frequency of uniqueness is equal to one, representing a confidence zone where the problem is well-posed regardless of k-means initializations. In the high $N_D$ zone, the frequency is zero, indicating that the problem has multiple solutions. In between lies a long transition zone, where the chances of finding a unique solution are getting lower and depend on the clustering.

In Figure 6, the standard deviation of the basis dimension is analyzed as a
<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution of the average intersection basis dimension as a function of the database size. The shaded area corresponds to the standard deviation of said size.</figcaption></figure>

function of $N_D$, computed across 100 different k-means initializations. The second method, based on the angles between subspaces, is employed to determine the dimension of the intersection space. For small values of $N_D$, the standard deviation of the basis dimension remains close to zero, indicating minimal variability across initializations and confirming that the problem is well-posed in this range. After a certain threshold, as $N_D$ increases, the standard deviation rises slowly and, the computed basis dimension remains globally small.

As shown in Figure 3b, the error decreases as the size of the database increases, up to the point where uniqueness is lost. This leads to the question of whether it is beneficial to correct the solutions beyond this point to gain more accuracy. The goal of Figures 5 and 6 is to examine the consequences of losing the uniqueness of our solution, how severe the problem is, and determine whether it can be resolved. The result of both figures reveals that the transition zone is wide, and within this region, the dimension is generally small. Moreover, since the difference between two distinct solutions lies in the kernel intersection (as shown in the proof of Proposition 2.1), this suggests that the loss of well-posedness can be effectively corrected, for example, by adding a regularization term.
## 6. Conclusion

In this paper, we present a mathematical and numerical analysis of the robustness of the Data-Driven Identification (DDI) method based on the preliminary work of Leygue [25], which introduced a criterion for the uniqueness of solutions. The DDI solution is unique if and only if there does not exist a non-zero stress field that is constant per cluster (in ker $H$) and self-balanced (in ker $D$) simultaneously.

We propose three methods to numerically characterize this criterion: explicit calculation of the basis, examination of the angles between the subspaces, and alternating projection onto the two kernels. A comprehensive theoretical evaluation of their efficiency is provided. Furthermore, we address the limitations of the current DDI algorithm by introducing a new version based on alternating projection, with a proof of convergence and improved initialization. The numerical conditioning of the problem and the convergence rate have also been discussed.

These methods were tested on synthetic data, and their effectiveness was analyzed, including a comparison of the computational costs. Although all three methods gave similar results, the third method, based on alternating projection, demonstrated the highest computational efficiency in terms of time and memory usage. The DDI method can now be applied more effectively while determining whether the problem has a unique solution and assessing its convergence. This analysis provides a way for selecting the parameter $N_D$, representing the database size. We also examined the implications of losing the well-posedness of the problem, which can be identified and corrected using the proposed approaches.

Finally, this theoretical analysis of the DDI framework opens promising perspectives for optimizing experimental setups and refining clustering methods. It has been observed that while adjusting DDI parameters to gain accuracy, the problem becomes ill-posed. In such cases, stress estimates can be adjusted through additional non-invasive regularization techniques or by identifying the specific degrees of freedom responsible for the loss of uniqueness. The matrix $H$, associated with clustering, and the matrix $D$, associated with the mesh, have incompatible sparsity patterns. To overcome this incompatibility, alternative clustering methods can be applied to establish a stronger connection between clustering structures and mesh structures. Although the applications presented focused on nonlinear elastic materials and synthetic data, the theoretical results remain valid for more general
cases. Future work could explore applying the framework to real experi-
mental data, and more complex behaviors, such as bi-materials samples and
history-dependent responses.

Acknowledgments

The authors gratefully acknowledge the financial support of NExT ISite
program of Nantes Université through International Research Project iD-
DrEAM.

Declarations

A CC-BY 4.0 public copyright license has been applied by the authors to
the present document and will be applied to all subsequent versions up to
the Author Accepted Manuscript arising from this submission.

References

[1] K. T. Kavanagh, R. W. Clough, Finite element applications in the char-
acterization of elastic solids, International Journal of Solids and Struc-
tures 7 (1) (1971) 11–23.

[2] M. Grédiac, Principe des travaux virtuels et identification, Comptes ren-
dus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie,
Sciences de l'univers, Sciences de la Terre 309 (1) (1989) 1–5.

[3] E. Jones, J. D. Carroll, K. N. Karlson, S. L. B. Kramer, R. B. Lehoucq, P. L. Reu, D. Z. Turner, Parameter covariance and non-uniqueness in material model calibration using the virtual fields method, Computational Materials Science 152 (2018) 268–290.

[4] E. Florentin, G. Lubineau, Identification of the parameters of an elastic material model using the constitutive equation gap method, Computational Mechanics 46 (2010) 521–531.

[5] A. Moussawi, G. Lubineau, E. Florentin, B. Blaysat, The constitutive compatibility method for identification of material parameters based on full-field measurements, Computer methods in applied mechanics and engineering 265 (2013) 1–14.
[6] M. I. Diaz, W. Aquino, M. Bonnet, A modified error in constitutive equation approach for frequency-domain viscoelasticity imaging using interior data, Computer methods in applied mechanics and engineering 296 (2015) 129–149.

[7] M. Bonnet, P. Salasiya, B. B. Guzina, Modified error-in-constitutive-relation (mecr) framework for the characterization of linear viscoelastic solids, Journal of the Mechanics and Physics of Solids 190 (2024) 105746.

[8] S. Avril, M. Bonnet, A.-S. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Lautourte, D. Lemosse, S. Pagano, E. Pagnacco, et al., Overview of identification methods of mechanical parameters based on full-field measurements, Experimental Mechanics 48 (2008) 381–402.

[9] J. Martins, A. Andrade-Campos, S. Thuillier, Comparison of inverse identification strategies for constitutive mechanical models using full-field measurements, International Journal of Mechanical Sciences 145 (2018) 330–345.

[10] M. Flaschel, S. Kumar, L. De Lorenzis, Unsupervised discovery of interpretable hyperelastic constitutive laws, Computer Methods in Applied Mechanics and Engineering 381 (2021) 113852.

[11] F. As’ad, P. Avery, C. Farhat, A mechanics-informed artificial neural network approach in data-driven constitutive modeling, International Journal for Numerical Methods in Engineering 123 (12) (2022) 2738–2759.

[12] P. Thakolkaran, A. Joshi, Y. Zheng, M. Flaschel, L. De Lorenzis, S. Kumar, Nn-euclid: Deep-learning hyperelasticity without stress data, Journal of the Mechanics and Physics of Solids 169 (2022) 105076.

[13] A. Benady, E. Baranger, L. Chamoin, Nn-mcre: A modified constitutive relation error framework for unsupervised learning of nonlinear state laws with physics-augmented neural networks, International Journal for Numerical Methods in Engineering 125 (8) (2024) e7439.

[14] J. Neggers, F. Mathieu, F. Hild, S. Roux, N. Swiergiel, Improving full-field identification using progressive model enrichments, International Journal of Solids and Structures 118 (2017) 213–223.
[15] A. Leygue, M. Coret, J. Réthoré, L. Stainier, E. Verron, Data-based derivation of material response, Computer Methods in Applied Mechanics and Engineering 331 (2018) 184–196.

[16] M. Grédiac, F. Hild, Full-field measurements and identification in solid mechanics, John Wiley & Sons, 2012.

[17] N. McCormick, J. Lord, Digital image correlation, Materials today 13 (12) (2010) 52–54.

[18] T. Kirchdoerfer, M. Ortiz, Data-driven computational mechanics, Computer Methods in Applied Mechanics and Engineering 304 (2016) 81–101.

[19] M. Dalémat, M. Coret, A. Leygue, E. Verron, Measuring stress field without constitutive equation, Mechanics of Materials 136 (2019) 103087.

[20] L. Costecalde, A. Leygue, M. Coret, E. Verron, Data-driven identification of hyperelastic models by measuring the strain energy density field, Rubber Chemistry and Technology 96 (4) (2023) 443–454.

[21] G. Valdés-Alonzo, C. Binetruy, B. Eck, A. García-González, A. Leygue, Phase distribution and properties identification of heterogeneous materials: A data-driven approach, Computer Methods in Applied Mechanics and Engineering 390 (2022) 114354.

[22] R. Langlois, M. Coret, J. Réthoré, Non-parametric stress field estimation for history-dependent materials: Application to ductile material exhibiting piobert–lüders localization bands, Strain 58 (4) (2022) e12410.

[23] A. Vinel, R. Seghir, J. Berthe, G. Portemont, J. Réthoré, Experimental characterization of material strain-rate dependence based on full-field data-driven identification, International Journal of Impact Engineering 194 (2024) 105083.

[24] M. Dalémat, M. Coret, A. Leygue, E. Verron, Robustness of the data-driven identification algorithm with incomplete input data, Journal of Theoretical, Computational and Applied Mechanics (2023).
[25] A. Leygue, On the formulation and convergence of data-driven identification, Comptes Rendus. Mécanique 353 (G1) (2025) 761-773.

[26] J. Mandel, Généralisation de la théorie de plasticité de wt koiter, International Journal of Solids and Structures 1 (3) (1965) 273-295.

[27] G. Strang, Linear Algebra and Its Applications, 4th Edition, Cengage Learning, 2012.

[28] Å. Björck, G. H. Golub, Numerical methods for computing angles between linear subspaces, Mathematics of computation 27 (123) (1973) 579-594.

[29] G. H. Golub, C. F. Van Loan, Matrix computations, JHU press, 2013.

[30] S. Boyd, J. Dattorro, Alternating projections, EE392o, Stanford University (2003).

[31] D. Bindel, Lecture notes for CS 6210: Matrix computations, https://www.cs.cornell.edu/~bindel/class/cs6210-f16/lec/2016-10-17.pdf (2016).

[32] J. W. Demmel, Applied numerical linear algebra, SIAM, 1997.

[33] F. Deutsch, Best approximation in inner product spaces, Vol. 7, Springer, 2001.

[34] L. Stainier, A. Leygue, M. Ortiz, Model-free data-driven methods in mechanics: material data identification and solvers, Computational Mechanics 64 (2) (2019) 381-393.

[35] D. Arthur, S. Vassilvitskii, k-means++: The advantages of careful seeding, Tech. rep., Stanford (2006).