Paper deep dive
Benchmarking data-driven material models on the classic Treloar dataset
Hagen Holthusen, Moritz Flaschel, Denisa Martonová, Ellen Kuhl
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 8/17/2026, 4:55:35 AM
Summary
This paper benchmarks six data-driven machine learning frameworks for hyperelastic constitutive modeling against the classic Treloar dataset for natural rubber. The methods compared are EUCLID, Material Fingerprinting (MF), Adaptive Material Fingerprinting (AMF), Physics-Augmented Neural Networks (PANN), Constitutive Artificial Neural Networks (CANN), and Generalized Invariant-based CANN (GI-CANN). The study evaluates these methods based on fitting performance (R², MSE), computational cost (evaluation time), model complexity (number of parameters), and hyperparameter sensitivity. Results indicate that all methods reproduce the benchmark data remarkably well, with no single winner, highlighting trade-offs between predictive accuracy, complexity, and computational efficiency.
Entities (9)
Relation Signals (13)
EUCLID → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... EUCLID
Material Fingerprinting → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... Material Fingerprinting
Adaptive Material Fingerprinting → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... Adaptive Material Fingerprinting
Physics-Augmented Neural Networks → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... Physics-Augmented Neural Networks
Constitutive Artificial Neural Networks → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... Constitutive Artificial Neural Networks
Generalized Invariant-based CANN → benchmarkedon → Treloar dataset
confidence 95% · In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks... Generalized Invariant-based CANN
Generalized Invariant-based CANN → hasacronym → GI-CANN
confidence 95% · Generalized Invariant-based CANN GI-CANN
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Machine learning is rapidly reshaping constitutive modeling, offers new ways to learn material behavior directly from experimental data, and challenges long-established modeling paradigms. But with a growing number of machine-learning-based approaches available, how do they compare in practice? In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks for hyperelasticity: (Generalized-Invariant) Constitutive Artificial Neural Networks, Physics-Augmented Neural Networks, (Adaptive) Material Fingerprinting, and Efficient Unsupervised Constitutive Law Identification & Discovery. We compare their fitting performance, computational cost, hyperparameter sensitivity, and ease of implementation. Furthermore, we discuss the trade-offs between predictive accuracy and model complexity. The latter is assessed by quantifying both the number of material parameters in the discovered models and the computational time required to evaluate the constitutive model and its derivatives. The results show that all methods can reproduce the benchmark data remarkably well. Rather than identifying a single winner, we highlight the strengths and limitations of each approach and provide practical guidance for their use. The source code for all six methods, including the training and comparison scripts, as well as all results and data used in this study, is publicly available via this https URL.
Tags
Links
- Source: https://arxiv.org/abs/2608.14063v1
- Canonical: https://arxiv.org/abs/2608.14063v1
Trouble viewing inline? Open PDF directly →
Full Text
92,150 characters extracted from source content.
Expand or collapse full text
Benchmarking data-driven material models on the classic Treloar dataset Hagen Holthusen Corresponding author: Correspondence: hagen.holthusen@fau.de Address: Institute of Applied Mechanics, Egerlandstraße 5, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany Moritz Flaschel Address: Institute of Applied Mechanics, Egerlandstraße 5, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany Denisa Martonová Address: Institute of Applied Mechanics, Egerlandstraße 5, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany Ellen Kuhl Address: Institute of Applied Mechanics, Egerlandstraße 5, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany Address: Department of Mechanical Engineering, Stanford University, 440 Escondido Mall, California 94305, United States. Abstract Machine learning is rapidly reshaping constitutive modeling, offers new ways to learn material behavior directly from experimental data, and challenges long-established modeling paradigms. But with a growing number of machine-learning-based approaches available, how do they compare in practice? In this paper, we use the classic experimental data of Treloar to benchmark popular frameworks for hyperelasticity: (Generalized-Invariant) Constitutive Artificial Neural Networks, Physics-Augmented Neural Networks, (Adaptive) Material Fingerprinting, and Efficient Unsupervised Constitutive Law Identification & Discovery. We compare their fitting performance, computational cost, hyperparameter sensitivity, and ease of implementation. Furthermore, we discuss the trade-offs between predictive accuracy and model complexity. The latter is assessed by quantifying both the number of material parameters in the discovered models and the computational time required to evaluate the constitutive model and its derivatives. The results show that all methods can reproduce the benchmark data remarkably well. Rather than identifying a single winner, we highlight the strengths and limitations of each approach and provide practical guidance for their use. The source code for all six methods, including the training and comparison scripts, as well as all results and data used in this study, is publicly available via Zenodo. Keywords: machine learning, automated discovery, material modeling, hyperelasticity 1 Introduction Treloar’s experiments on natural rubber (37) constitute one of the most widely used benchmarks for the development and evaluation of constitutive models. A systematic classical comparison assessed fourteen phenomenological and micromechanically motivated hyperelastic models against these data and examined whether parameters identified from one deformation mode could predict the others (32). More than ten years later, the emergence of machine learning and automated model discovery methods has since changed the modeling landscape: rather than selecting exclusively from a limited set of prescribed constitutive forms, these approaches can learn or assemble the strain energy density function from data, with different balances of physical structure, flexibility, sparsity, and interpretability. This shift introduces new practical questions concerning generalization, optimization robustness, hyperparameter and initialization sensitivity, and computational cost. A renewed benchmark on the same established dataset is therefore timely. Here, we compare six representative data-driven approaches under common conditions. An overview of the considered methods is provided in Table 1. Fig. 1 shows their discovered constitutive functions, and Fig. 2 shows their model response in comparison with Treloar’s data. Table 1: Benchmark performance of data-driven methods (see Section 4 and Table 5 for details). Method Acronym Training R2R^2 Validation R2R^2 Core time (ms) # Parameters ‖0\| θ\|_0 Efficient Unsupervised Constitutive Law Identification & Discovery EUCLID 0.9973 0.9978 0.372 4 Material Fingerprinting MF 0.9781 0.9899 0.302 2 Adaptive Material Fingerprinting AMF 0.9991 0.9996 6.000 49 Physics-Augmented Neural Networks PANN 0.9973 0.9977 581.311 371 Constitutive Artificial Neural Networks CANN 0.9972 0.9975 320.947 12 Generalized Invariant-based CANN GI-CANN 0.9996 0.9993 664.988 8 To quantitatively assess the performance of the considered methods on Treloar’s dataset, we introduce several metrics that evaluate the fitting accuracy of the resulting models. In addition to fitting accuracy, we investigate the computational cost of the considered data-driven methods and evaluate the complexity of the resulting constitutive models using several complementary metrics. Constitutive model complexity is particularly relevant for computational simulations, where the repeated evaluation of the model and its derivatives at every integration point can have a significant impact on the overall computational cost. Table 1 summarizes some representative metrics used to quantify fitting accuracy, model complexity, and the computational cost of the inverse problem. Additional performance measures will be introduced and discussed throughout this work. Figure 1: Strain energy density in the principal stretch space predicted by the six selected best-performing models (see Section 4 for details). Figure 2: Comparison of Treloar’s experimental data and predictions from the best-performing discovered models. We note that the literature on data-driven constitutive modeling is extensive and continues to grow rapidly. Consequently, our comparison is restricted to a representative set of widely used approaches and does not claim to cover all state-of-the-art methods or identify the universally best-performing techniques. Specifically, we consider Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) (13), Material Fingerprinting (MF) (11; 29), Adaptive Material Fingerprinting (AMF) (9), Physics-Augmented Neural Networks (PANNs) (23; 25; 6), Constitutive Artificial Neural Networks (CANNs) (26; 27), and Generalized Invariant-based CANNs (GI-CANNs) (28). This list of methods is by no means exhaustive. Numerous other approaches have been proposed in the literature (17; 33), including, to name only a few, Gaussian processes (15; 18), Variational System Identification (VSI) (38), Neural Ordinary Differential Equations (NODE) (34), symbolic regression (1), spline-based approaches (39), non-smooth parameterizations of constitutive functions (4), and agentic artificial intelligence (35). In the following Section 2, we briefly review the classical Treloar dataset on natural rubber. In the subsequent Section 3, we introduce the considered data-driven methods for constitutive modeling. Next, we apply these methods to the benchmark dataset and analyze and discuss the results in Section 4. 2 Benchmark problem To benchmark different data-driven material modeling approaches with respect to predictive accuracy and computational efficiency, we apply several widely used methods to the same dataset under identical conditions. Before introducing the individual methods, we first describe the benchmark dataset and its preprocessing, the training and test data split, the training objective employed by the data-driven models, and the evaluation metrics used to quantify predictive performance. 2.1 Dataset and preprocessing In this work, we use the classical Treloar rubber dataset at 20∘20 C as a benchmark for data-driven constitutive model discovery (37). The dataset comprises stress–stretch measurements obtained from three canonical homogeneous deformation modes: uniaxial tension (UT), equibiaxial tension (ET), and pure shear (PS). The experimental data are given as stretch values λ and corresponding nominal first Piola-Kirchhoff stress values P11P_11. Since the constitutive models operate on the full three-dimensional deformation gradient, each scalar stretch value is first converted into a three-dimensional deformation gradient F. Assuming incompressibility, the deformation gradients are constructed as UT=diag(λ,λ−1/2,λ−1/2),ET=diag(λ,λ,λ−2),PS=diag(λ,1,λ−1), F^UT=diag (λ,λ^-1/2,λ^-1/2 ), F^ET=diag (λ,λ,λ^-2 ), F^PS=diag (λ,1,λ^-1 ), (1) for uniaxial tension, equibiaxial tension, and pure shear, respectively. For each experiment, a set of stretch–stress data tuples is available, see Fig. 2, UT=(λiUT,P11,iUT)|i=1,…,NUT,ET=(λiET,P11,iET)|i=1,…,NET,PS=(λiPS,P11,iPS)|i=1,…,NPS.D^UT= \( _i^UT,P_11,i^UT)\,|\,i=1,…,N_UT \, ^ET= \( _i^ET,P_11,i^ET)\,|\,i=1,…,N_ET \, ^PS= \( _i^PS,P_11,i^PS)\,|\,i=1,…,N_PS \. (2) For training, the uniaxial and equibiaxial datasets are concatenated into one calibration dataset, and the pure shear data are used for testing train=UT∪ET,test=PS.D_train=D^UT ^ET, _test=D^PS. (3) The pure shear data are not used during training and are retained as an independent validation case. This allows us to assess whether the discovered constitutive model generalizes to a deformation mode that was not part of the calibration data. To improve numerical stability during optimization, the nominal stress targets are normalized by the maximum absolute stress value in the training set P^11,i=P11,imaxj=1,…,Ntrain|P11,j|. P_11,i= P_11,i _j=1,…,N_train|P_11,j|. (4) The data-driven models are trained using the normalized stress values. During evaluation, the predicted stresses are transformed back to physical units using the corresponding normalization factor. No smoothing, filtering, or data augmentation is applied, such that the comparison remains directly based on the original experimental measurements. 2.2 Training objective All data-driven methods considered in this work seek to identify a parametric representation of the constitutive model. In the context of hyperelasticity, this corresponds to learning a strain energy density function W_ θ parameterized by the model parameters θ. The parameters are learned in a supervised manner by minimizing the discrepancy between the model predictions and the experimentally measured nominal stresses. For each deformation gradient i F_i in the training dataset, the model predicts the first Piola–Kirchhoff stress from the learned strain energy density (i)=∂W(i)∂i. P_ θ( F_i)= ∂ W_ θ( F_i)∂ F_i. (5) For incompressible materials, the stress is corrected by a pressure-like Lagrange multiplier so that the transverse stress condition is satisfied. In the present implementation, this pressure correction is chosen from the condition P33=0P_33=0 which removes the hydrostatic contribution from the predicted first Piola–Kirchhoff stress. The loss is therefore evaluated on the physically admissible incompressible stress response rather than on the unconstrained energy gradient. Since the experimental data provide the nominal stress component in the loading direction, the loss is formulated on the normalized P11P_11 component. The training objective is the mean squared error ℒdata()=1Ntrain∑i=1Ntrain(P^11,(i)−P^11,i)2,L_data( θ)= 1N_train _i=1^N_train ( P_11, θ( F_i)- P_11,i )^2, (6) where P^11,i P_11,i denotes the normalized experimental stress and P^11, P_11, θ denotes the corresponding model prediction normalized by the same factor as the experimental stress. The total optimization problem is given by ∗=argminℒdata(), θ = _ θL_data( θ), (7) possibly supplemented by architecture-dependent constraints or regularization terms that enforce constitutive structure, positivity, convexity, or sparsity. Finally, we note that not all data-driven methods considered in this work are formulated as continuous optimization problems of the form above. In particular, database-driven Material Fingerprinting methods instead rely on pattern-recognition algorithms operating on a pre-established material database. 2.3 Evaluation metrics In this work, we compare different data-driven methods based on several evaluation criteria. These include the computational time required to identify or discover them from the given data, the goodness of fit of the identified models, and the interpretability or simplicity of the resulting models. We introduce different metrics for quantifying the goodness of fit. For a given dataset D with a number of |||D| data points, we consider the mean squared error (MSE) MSE=1||∑i=1||(P11,(λi)−P11,i)2,MSE_D= 1|D| _i=1^|D| (P_11, θ( _i)-P_11,i )^2, (8) the root mean squared error (RMSE) RMSE=1||∑i=1||(P11,(λi)−P11,i)2,RMSE_D= 1|D| _i=1^|D| (P_11, θ( _i)-P_11,i )^2, (9) the range-normalized root mean squared error (NRMSE) NRMSE=RMSEmaxj=1,…,|||P11,j|−minj=1,…,|||P11,j|,NRMSE_D= RMSE_D _j=1,…,|D||P_11,j|- _j=1,…,|D||P_11,j|, (10) the mean absolute error MAE=1||∑i=1|||P11,(λi)−P11,i|,MAE_D= 1|D| _i=1^|D| |P_11, θ( _i)-P_11,i |, (11) and the coefficient of determination R2= 1−∑i=1||(P11,(λi)−P11,i)2∑i=1||(mean(P11,i)−P11,i)2,R^2_D\;=\;1- _i=1^|D| (P_11, θ( _i)-P_11,i )^2 _i=1^|D| (mean_D(P_11,i)-P_11,i )^2, (12) where mean(P11,i)mean_D(P_11,i) is the arithmetic mean of all experimental stress values in the dataset. In constitutive modeling, it is often the case that multiple models provide an equally accurate description of the available data, while differing substantially in their complexity, interpretability, and computational cost. For practical applications, particularly in large-scale computational simulations, simpler constitutive models are generally preferred, as they require less computational effort to evaluate the strain energy density function and its derivatives. Consequently, in addition to predictive accuracy, model complexity constitutes an important criterion for assessing and comparing constitutive models. We consider different metrics to quantify the model complexity of the identified models. The first measure of model complexity considered in this work is the effective number of material parameters (or trainable weights in the terminology of machine learning). Specifically, we count only those parameters that remain nonzero after the model identification procedure, as these determine the complexity of the resulting constitutive model. We denote this quantity by ‖0\| θ\|_0, in analogy with the L0L_0-pseudo-norm, which counts the number of nonzero entries in a vector. While the number of material parameters ‖0\| θ\|_0 provides a useful proxy for model complexity, it does not directly reflect the computational cost of evaluating a strain energy density function and its derivatives. Constitutive models with similar numbers of parameters may require substantially different numbers and types of operations, resulting in different computational costs when deployed in finite element simulations. To quantify this aspect, we implemented each identified strain energy density function in a common machine learning framework and used automatic differentiation to evaluate the energy W, its gradient P, and its Hessian ℍH. Each benchmark evaluated a batch of Neval=100N_eval=100 deformation gradients in double precision on the CPU. After two warm-up evaluations, the computation was repeated 100 times. The total runtime of each repetition was divided by NevalN_eval, yielding the runtime per deformation gradient. Compilation, parameter loading, sample generation, validation, and host result transfer were excluded. The resulting energy, gradient, and Hessian evaluation times are reported as the mean ± standard deviation over the 100 repetitions. 3 Methods In the following, we briefly introduce the data-driven material modeling approaches considered in this work. Detailed descriptions of the individual methods can be found in the corresponding original publications. An overview of the considered approaches is provided in Table 2, which summarizes the optimization algorithms and solvers employed for model identification from experimental data, together with the associated hyperparameters and modeling choices. These hyperparameters and modeling choices are discussed in detail for each method in the subsequent subsections. Table 2: Hyperparameters of machine learning methods for constitutive modeling. Method Solvers / Algorithms Hyperparameters Modeling Decisions EUCLID ∙ Coordinate Descent (CD) ∙ Least Angle Regression (LARS) ∙ Number of terms in the library ∙ Regularization parameter α ∙ Thresholding parameter θ→0 _ 0.5[1]$→$0 ∙ Choice of the library features ∙ Model input variables Material Fingerprinting ∙ Pattern-recognition ∙ Number of fingerprints in the database ∙ Parameter sampling during database generation ∙ Choice of the models in the database ∙ Model input variables Adaptive Material Fingerprinting ∙ Adaptive pattern-recognition ∙ Number of fingerprints in the database ∙ Parameter sampling during database generation ∙ Number of terms in the model NaN_a ∙ Step size s∈[0,1]s∈[0,1] ∙ Choice of the modeling features in the database ∙ Model input variables PANN Various (usually gradient-based) ∙ SGD ∙ Adam ∙ L-BFGS ∙ Learning rate η ∙ Number of epochs ∙ Number of layers ∙ Number of neurons per layer ∙ Optional ML hyperparameters (e.g., regularization, early stopping, dropout, scheduling) ∙ Choice of the activation functions ∙ Model input variables CANN and GI-CANN Various (usually gradient-based) ∙ SGD ∙ Adam ∙ L-BFGS ∙ Learning rate η ∙ Number of epochs ∙ Number of layers ∙ Number of neurons per layer ∙ Optional ML hyperparameters (e.g., regularization, early stopping, dropout, scheduling) ∙ Choice of the activation functions ∙ Model input variables For the benchmark problem considered in this work, all data-driven methods share the common objective of constructing a mathematical representation of the hyperelastic strain energy density function. In the present study, we restrict our attention to isotropic incompressible hyperelasticity, for which the strain energy density function can be expressed as W~()=W()−p[J−1], W( F)=W( F)-p[J-1], (13) where J=det()J= ( F) and p is a scalar Lagrange multiplier enforcing the incompressibility constraint J=1J=1. To ensure material isotropy, the constitutive contribution W()W( F) is commonly expressed in terms of invariants of the right Cauchy–Green deformation tensor =T C= F^T F. The most common choices are the first and second principal invariants, I1=tr()I_1=tr( C) and I2=12[(tr)2−tr(2)]I_2= 12 [(tr C)^2-tr( C^2) ], or the principal stretches λ1,λ2,λ3 _1, _2, _3, whose squares are the eigenvalues of C, i.e., λ12,λ22,λ32=eig()\ _1^2, _2^2, _3^2\=eig( C). These quantities satisfy I1=λ12+λ22+λ32I_1= _1^2+ _2^2+ _3^2, I2=λ12λ22+λ22λ32+λ32λ12I_2= _1^2 _2^2+ _2^2 _3^2+ _3^2 _1^2, and J=λ1λ2λ3=1J= _1 _2 _3=1. Alternatively, generalized invariant sets may be employed, as discussed below for the GI-CANN framework. By differentiating W~ W with respect to the deformation gradient and determining p from the constraint P33=0P_33=0 for the dataset introduced in Section 2, we obtain the P11P_11 component of the first Piola–Kirchhoff stress as P11=∂W~∂F11=∂W∂F11−∂W∂F33[cof()33]−1cof()11.P_11= ∂ W∂ F_11= ∂ W∂ F_11- ∂ W∂ F_33 [cof( F)_33 ]^-1cof( F)_11. (14) The following sections describe how EUCLID, Material Fingerprinting, Adaptive Material Fingerprinting, PANNs, CANNs, and GI-CANNs identify, discover, or learn the constitutive contribution W. 3.1 Efficient Unsupervised Constitutive Law Identification and Discovery (EUCLID) EUCLID is a method for the automated discovery of interpretable material models with a small number of modeling terms from data (10). The underlying idea is to construct a library or catalog of modeling terms and use sparse regression to identify the combination of terms that best describes the given data. In its original form (10), EUCLID is unsupervised, which means that it does not necessarily rely on stress-stretch data pairs but can instead be informed by full-field displacement measurements of complexly shaped specimens. However, the concepts of EUCLID can equally be applied in a supervised fashion when stress-stretch data pairs are available (13), as is the case for Treloar’s data. In the following, we demonstrate how EUCLID can be applied to Treloar’s data to efficiently discover material models. At the heart of EUCLID stands a set of candidate modeling terms that could potentially describe the material behavior. In the case of isotropic hyperelasticity, this can be, for example, the terms of the generalized Mooney-Rivlin model, which are motivated by a Taylor expansion of the invariant-based strain energy density about the undeformed configuration. We collect these modeling terms in the feature vector (I1,I2)=[[I1−3],[I2−3],[I1−3]2,[I1−3][I2−3],[I2−3]2,…]T. Q(I_1,I_2)= [[I_1-3],\,[I_2-3],\,[I_1-3]^2,\,[I_1-3][I_2-3],\,[I_2-3]^2,\,… ]^T. (15) In this work, we consider coupled polynomial features up to an order of five, resulting in Nfeat=20N_feat=20 modeling features. We note, however, that the library can be readily extended with more modeling features, such as those of the Ogden material model (13). A linear combination of the candidate modeling terms yields the strain energy density W(I1,I2)=(I1,I2)⋅,W(I_1,I_2)= Q(I_1,I_2)· θ, (16) where ∈ℝ≥0Nfeat θ _≥ 0^N_feat are material parameters, which we assume here to be non-negative. The objective of EUCLID is to identify the parameters θ such that model predictions are similar to the experimental measurements while, at the same time, removing candidate modeling terms that are irrelevant for describing the data. To this end, we consider the regularized optimization problem ∗=argmin∈ℝ≥0ℒdata()+α‖1, θ = _ θ _≥ 0L_data( θ)+α\| θ\|_1, (17) in which the L1L_1-norm penalizes solutions with many non-zero parameters and promotes sparsity in the solution (14; 36; 19; 5). The regularization parameter α≥0α≥ 0 determines how aggressively the parameters should be pushed to zero. Since our model is linear in the parameters θ (see Eq. 16), the loss term in the optimization problem above is quadratic in θ. The L1L_1-norm regularization term is non-smooth but convex in θ. The total objective function is convex and non-smooth and can be efficiently minimized using methods such as Coordinate Descent (CD) (16) or Least Angle Regression (LARS) (7). Such non-smooth optimization approaches have also been applied to constitutive model discovery and material modeling (8). Solving the problem above yields a parameter vector θ with many zero entries and, hence, a material model that is expressed by a small number of modeling terms. We note that some parameters may not be exactly zero but may instead attain values very close to zero. In such cases, a threshold parameter θ→0 _ 0.5[1]$→$0 can be introduced in a post-processing step, whereby parameters with absolute values below the threshold are set to zero (10). For the results presented in this paper, however, such a thresholding step was not necessary to obtain a sparse representation of the model and was therefore omitted. Furthermore, once EUCLID has identified a suitable sparse set of constitutive terms, it can be beneficial to perform a final regression with α=0α=0 to optimize only the remaining nonzero parameters. This post-processing step removes the bias introduced by the regularization and can further improve model accuracy, but it has not been employed in the present work. Hyperparameters As shown in Table 2, EUCLID depends on several hyperparameters, including the number of terms in the candidate library, the regularization parameter α, and the threshold parameter θ→0 _ 0.5[1]$→$0. In addition, its performance is influenced by several modeling choices made during the construction of the library. In particular, the variables on which the library features depend (here, the principal invariants), as well as the functional forms included in the library must be specified a priori. We finally note that previous works on EUCLID have also investigated the nonconvex LpL_p-regularization term with 0<p<10<p<1 (10; 30). However, since increasing the regularization parameter α allows the L1L_1-norm to achieve any desired level of sparsity, we restrict our attention here to the convex L1L_1-regularization. 3.2 Material Fingerprinting (MF) Material Fingerprinting is a database-driven method for constitutive model discovery that replaces the solution of a continuous optimization problem with a pattern-recognition task (11; 29). The key assumption is that every material exhibits a characteristic mechanical response when subjected to a standardized experimental protocol. This response defines the material’s fingerprint. Once a sufficiently rich database of fingerprints exists, material characterization reduces to identifying the closest fingerprint in the database. Material Fingerprinting consists of an offline and an online stage. During the offline stage, we generate a database of fingerprints by numerically simulating standardized experiments for different constitutive models (Blatz-Ko, Demiray, Gent, Holzapfel, Mooney-Rivlin, Neo-Hooke, Ogden) and parameter combinations. Each constitutive model is represented by its specific strain energy density function W and sampled over a prescribed range of material parameters (29). In the supervised setting considered here, the experiments correspond to homogeneous deformation modes, such as uniaxial tension and simple shear, and each fingerprint collects the stresses at prescribed deformation states. More generally, the method also applies to heterogeneous experiments that provide full-field displacement and reaction force measurements (12). Each database entry stores a fingerprint together with its constitutive model and material parameters ((i),(i)),i=1,…,Ndb, (f^(i), θ^(i) ), i=1,…,N_db, (18) where NdbN_db denotes the number of fingerprints in the database. To remove the dependence on the overall stiffness scale, we normalize every fingerprint ¯(i)=(i)‖(i)‖, f^(i)= f^(i)\|f^(i)\|, (19) which exploits the homogeneity of many hyperelastic constitutive models with respect to their material parameters. As a result, materials that differ only by a multiplicative scaling of their parameters share the same normalized fingerprint, which substantially reduces the required database size (11). During the online stage, we measure and normalize the fingerprint ∗f of an unknown material. We then compare the measured fingerprint with every database entry through the cosine similarity cs(i)=¯(i)⋅¯∗|¯(i)||¯∗|=¯(i)⋅¯∗,cs^(i)= f^(i)· f | f^(i) | | f |= f^(i)· f , (20) where the second equality follows because the fingerprints have unit norm. Equivalently, Material Fingerprinting identifies the fingerprint that forms the smallest angle with the experimental fingerprint. We therefore identify the constitutive model through i∗=argmaxics(i).i = _i\;cs^(i). (21) Unlike conventional parameter identification methods, Material Fingerprinting does not solve a continuous optimization problem during the online stage. Instead, it reduces constitutive model discovery to a matrix–vector multiplication followed by a maximum search. Moreover, because the algorithm searches the complete database, it identifies the discrete global optimum within the searchable model space instead of converging to a local optimum of a continuous optimization problem (11). Consequently, the method shifts the computational effort entirely to the offline stage, while the online characterization of new materials requires only milliseconds. In this work, we focus on supervised Material Fingerprinting for experiments with homogeneous deformation fields. We note that the computational speed-up is even more pronounced in the unsupervised case for experiments with heterogeneous deformation fields (12). Hyperparameters As shown in Table 2, the performance of Material Fingerprinting depends on several user-defined modeling choices and hyperparameters. For instance, the variables on which the constitutive models in the database are based (here, the principal invariants and principal stretches), as well as the admissible functional forms, must be specified a priori. In addition, the database generation requires selecting the sampling density of the material parameters, which determines the total number of fingerprints contained in the database. These hyperparameter choices affect only the database generation and therefore constitute a one-time effort. Once an optimized and sufficiently expressive database has been constructed, it can be reused repeatedly for efficient material model discovery. 3.3 Adaptive Material Fingerprinting (AMF) Adaptive Material Fingerprinting extends the original Material Fingerprinting framework by constructing the constitutive model incrementally instead of selecting a single predefined model from the database (9). Whereas the original method stores fingerprints of complete constitutive models, the adaptive method stores fingerprints of individual modeling features. It reconstructs the constitutive response step by step by successively adding the feature that best explains the discrepancy between the experimental data and the current model prediction. During the offline stage, we generate a database of feature fingerprints instead of fingerprints for complete constitutive models. We express the strain energy density as a linear combination of modeling features, W()=∑a=1NaθaQa(),W( F)= _a=1^N_a _a\,Q_a( F), (22) where NaN_a denotes the maximum number of feature terms included in the discovered constitutive model. We compute the stress response associated with every feature for different parameter choices. The feature library may contain isotropic and anisotropic terms based on principal stretches, principal area changes, or fiber stretches. This representation captures constitutive models such as the multi-term Ogden model or the Holzapfel-Gasser-Ogden model (9). Because we prescribe the admissible parameter ranges during database generation, we can also enforce physical constraints, such as polyconvexity, by restricting the database to admissible feature combinations. During the online stage, we initialize the predicted fingerprint with 0=.f_0=0. (23) At iteration k, we compute the residual fingerprint k=∗−k,r_k=f -f_k, (24) where ∗f denotes the experimental fingerprint and kf_k denotes the fingerprint predicted by the current constitutive model. We then compare the normalized residual with all feature fingerprints in the database through the cosine similarity ik∗=argmaxi¯(i)⋅¯k|¯(i)||¯k|.i_k = _i\; f^(i)· r_k | f^(i) | | r_k |. (25) After we identify the best matching feature, we add it to the constitutive model. Its contribution is scaled by a step size parameter s∈(0,1]s∈(0,1], which controls how much of the current residual is incorporated at each iteration. Thus, smaller values of s lead to more gradual updates, whereas s=1s=1 corresponds to incorporating the full residual contribution. We repeat the procedure for a prescribed maximum number of iterations NaN_a. Compared to the original Material Fingerprinting framework, the adaptive method increases the flexibility of the constitutive model because it combines multiple modeling features instead of selecting a single predefined model. At the same time, it retains the computational efficiency of database-based pattern recognition and, dependent on the choice of NaN_a, produces sparse, physically interpretable constitutive models without solving a continuous optimization problem (9). Hyperparameters As shown in Table 2, Adaptive Material Fingerprinting depends on a larger number of hyperparameters than the original Material Fingerprinting method. In particular, the number of terms NaN_a included in the model and the step size s must be specified a priori. Due to the computational efficiency of Adaptive Material Fingerprinting, suitable values for these hyperparameters can be determined using a simple grid search (9). 3.4 Physics-Augmented Neural Networks (PANNs) Physics-Augmented Neural Networks (23; 25) incorporate constitutive requirements directly into the neural network architecture. For the isotropic incompressible setting considered here, we represent the strain energy density by an input-convex neural network (ICNN) (2) that takes the first and second principal invariants as inputs, W(I1,I2)=ICNN(I1,I2)−ICNN(3,3),W(I_1,I_2)=ICNN_ θ(I_1,I_2)-ICNN_ θ(3,3), (26) where θ denotes the trainable network parameters. Since I1=I2=3I_1=I_2=3 in the undeformed reference configuration, the second term normalizes the strain energy density such that W(3,3)=0W(3,3)=0. This term is constant with respect to the deformation and therefore does not affect the predicted stresses or tangent operators. Moreover, subtracting a constant preserves the convexity and monotonicity properties of the ICNN. The additional stress normalization according to 25 is not required at this point, as the Lagrange multiplier p ensures a stress-free undeformed reference configuration. The choice of an ICNN is motivated by the construction of a polyconvex strain energy density. In three dimensions, the principal invariants can be written as I1=:,I2=cof():cof().I_1= F: F, I_2=cof( F):cof( F). (27) Thus, I1I_1 is a convex function of F, while I2I_2 is a convex function of cof()cof( F). If W(I1,I2)W(I_1,I_2) is convex and non-decreasing in each argument, the composition is jointly convex in F and cof()cof( F), regarded as independent arguments. Consequently, a strain energy density of the form W()=ICNN(:,cof():cof())−ICNN(3,3)W( F)=ICNN_ θ ( F: F,\,cof( F):cof( F) )-ICNN_ θ(3,3) (28) is polyconvex if the ICNN is convex and componentwise non-decreasing in its inputs. For incompressible deformations, I1≥3I_1≥ 3 and I2≥3I_2≥ 3. Componentwise monotonicity, together with the normalization at (I1,I2)=(3,3)(I_1,I_2)=(3,3), therefore additionally ensures W≥0W≥ 0 for all admissible deformation states. To describe how these properties are embedded into the network, we collect the invariants in the input vector =[I1,I2]T. x= bmatrixI_1,&I_2 bmatrix^T. (29) The first hidden representation is defined as 1=φ1(1+1), z_1= _1 ( U_1 x+ b_1 ), (30) and the subsequent representations are constructed recursively as ℓ=φℓ(ℓ−1+ℓ+ℓ),ℓ=2,…,L. z_ = _ ( W_ z_ -1+ U_ x+ b_ ), =2,…,L. (31) The final representation is scalar and defines the network output, ICNN()=L.ICNN_ θ( x)= z_L. (32) The ICNN can therefore be understood as a composition of the layer mappings gℓ(,)=φℓ(ℓ+ℓ+ℓ).g_ ( z; x)= _ ( W_ z+ U_ x+ b_ ). (33) Importantly, each layer depends not only on the preceding hidden representation but also directly on the original input through the term ℓ U_ x. These direct input paths constitute the skip connections of the ICNN architecture. Convexity of the ICNN with respect to x is ensured by choosing activation functions φℓ _ that are convex and non-decreasing and by imposing the componentwise constraints ℓ≥,ℓ=2,…,L. W_ ≥ 0, =2,…,L. (34) Indeed, the non-negative hidden-to-hidden weights preserve convexity when the representations from the preceding layer are combined, while the skip-connection term ℓ U_ x is affine in x. Input convexity alone would therefore not require a sign constraint on ℓ U_ . For the polyconvex construction above, however, the ICNN must additionally be non-decreasing in both I1I_1 and I2I_2. This is ensured by also constraining the input and skip-connection weights componentwise according to ℓ≥,ℓ=1,…,L. U_ ≥ 0, =1,…,L. (35) Together with the non-negative hidden-to-hidden weights and the non-decreasing activation functions, these constraints make every layer, and hence the complete network, componentwise non-decreasing in x. The biases ℓ b_ affect neither convexity nor monotonicity and remain unconstrained. Although the ICNN is convex with respect to its invariant inputs and the resulting strain energy density is polyconvex, the loss function is generally non-convex with respect to the network parameters. The identified model may therefore depend on their initialization. Hyperparameters As shown in Table 2, the design and training of PANNs depend on several user-defined modeling choices and hyperparameters. Prior to training, users must specify the network architecture, including the number of hidden layers, the number of neurons per layer, the activation functions, and the input variables on which the network depends. The training process is subsequently governed by hyperparameters such as the learning rate and the number of training epochs, as well as optional techniques including regularization, early stopping, dropout, and learning rate scheduling. 3.5 Constitutive Artificial Neural Networks (CANNs) Constitutive Artificial Neural Networks (26; 27) represent the strain energy density by a structured neural network whose architecture resembles a generalized constitutive model. Rather than learning an unrestricted mapping from deformation measures to energy, CANNs combine physically motivated invariant-based features through trainable activation functions and non-negative weights (21). For the isotropic incompressible setting considered here, the network takes the shifted first and second principal invariants I¯1=I1−3,I¯2=I2−3 I_1=I_1-3, I_2=I_2-3 (36) as inputs. Since I1≥3I_1≥ 3 and I2≥3I_2≥ 3 for incompressible deformations, both shifted invariants are non-negative over the admissible deformation domain. We construct the basic feature vector (I¯1,I¯2)=[I¯1,I¯2,I¯12,I¯22]T. q( I_1, I_2)= bmatrix I_1,& I_2,& I_1^2,& I_2^2 bmatrix^T. (37) Each feature enters the strain energy density both directly and through an exponential activation.11 1 For numerical stability, the arguments of the exponential activations are clipped to a finite interval in the implementation. This safeguard was introduced because the exponential branches acting on the quadratic features, of the form exp(w1,kexpI¯i2) (w_1,k^exp I_i^2), can otherwise overflow during training. The clipping is not part of the constitutive ansatz; the strict convexity argument applies over the range in which the clipping remains inactive. The resulting CANN can be written as W(I1,I2)=∑k=14w2,kidqk+∑k=14w2,kexp[exp(w1,kexpqk)−1],W(I_1,I_2)= _k=1^4w_2,k^id\,q_k+ _k=1^4w_2,k^exp [ (w_1,k^exp\,q_k )-1 ], (38) where w2,kidw_2,k^id are the output weights of the identity branches, w1,kexpw_1,k^exp control the exponential activations, and w2,kexpw_2,k^exp are the corresponding exponential output weights. Hence, the CANN comprises a total of 12 trainable weights. This construction ensures that W(3,3)=0W(3,3)=0 (39) is satisfied by construction, without requiring an additional energy normalization. All trainable weights are constrained componentwise according to w2,kid≥0,w1,kexp≥0,w2,kexp≥0,k=1,…,4.w_2,k^id≥ 0, w_1,k^exp≥ 0, w_2,k^exp≥ 0, k=1,…,4. (40) These constraints ensure that each constitutive branch contributes non-negatively to the strain energy density. Since the shifted invariants are non-negative, the linear and quadratic features qkq_k are convex and non-decreasing functions of their respective invariant. For w1,kexp≥0w_1,k^exp≥ 0, the exponential mapping qk↦exp(w1,kexpqk)−1q_k (w_1,k^exp\,q_k )-1 (41) is likewise convex and non-decreasing. A non-negative weighted sum of these terms therefore remains convex and componentwise non-decreasing in I1I_1 and I2I_2. As discussed for the PANN framework, these properties provide a sufficient construction for polyconvexity. In three dimensions, the invariants can be expressed in terms of F and cof()cof( F) according to Eq. 27 (20). Because the CANN is convex and non-decreasing in both invariants, its composition with these convex functions admits a jointly convex representation in F and cof()cof( F) (22). The resulting strain energy density is therefore polyconvex over the admissible incompressible deformation domain. Together with the invariant-based inputs, the architecture also ensures objectivity and material isotropy. Hyperparameters As shown in Table 2, the hyperparameters of CANNs are largely similar to those of the PANNs discussed earlier. The primary difference is that CANNs incorporate handcrafted and expert-guided functional forms directly into the network architecture. 3.6 Generalized-Invariant-based Constitutive Artificial Neural Networks (GI-CANNs) Generalized-Invariant-based Constitutive Artificial Neural Networks (28) extend the CANN framework by treating not only the functional form of the strain energy density, but also its invariant representation as trainable. Instead of restricting the model input to the classical invariants I1I_1 and I2I_2, GI-CANNs introduce the continuous family of generalized invariants (3) α=∑i=13λiα,J_α= _i=1^3 _i^α, (42) where λi _i are the principal stretches and α∈ℝα is an exponent that can be identified from the data. For incompressible materials, the classical invariants are recovered as the special cases 2=I1,−2=I2.J_2=I_1, _-2=I_2. (43) The trainable exponent therefore allows the network to continuously explore invariant representations beyond the fixed classical choices. In the general GI-CANN framework, the strain energy density depends on a set of generalized invariants, W=W(αkαk∈),W=W ( \J_ _k \_ _k ), (44) where ⊂ℝS denotes the set of invariant exponents. Fixing =2,−2S=\2,-2\ recovers the classical invariant representation underlying the standard-invariant-based CANN (27). In contrast, prescribing a discrete grid of exponents and combining the resulting generalized invariants through identity branches recovers the principal-stretch-based CANN (31). GI-CANN generalizes both approaches by treating the exponents as continuous trainable parameters and identifying the invariant representation and the strain energy density function simultaneously. For the benchmark implementation considered here, the exponent set contains one non-negative exponent and one non-positive exponent =α+,−α−,α+≥0,α−≥0.S= \ _+,- _- \, _+≥ 0, _-≥ 0. (45) The corresponding shifted generalized invariants are ¯+=∑i=13λiα+−3,¯−=∑i=13λi−α−3. J_+= _i=1^3 _i _+-3, J_-= _i=1^3 _i^- _--3. (46) In the numerical implementation, the principal stretches are computed as the singular values of the deformation gradient. Subtracting three ensures that both shifted generalized invariants vanish in the undeformed reference configuration. The general GI-CANN architecture may apply different powers and nonlinear transformations to the shifted generalized invariants. In the present implementation, we restrict the power set to the identity and use one linear and one exponential branch for each invariant. The strain energy density is therefore expressed as W=∑k∈+,−w2,kid¯k+w2,kexp[exp(w1,kexp¯k)−1].W= _k∈\+,-\ \w_2,k^id\, J_k+w_2,k^exp [ (w_1,k^exp J_k )-1 ] \. (47) Thus, the model simultaneously identifies two generalized invariant exponents, two weights controlling the exponential activations, and four output weights, resulting in a total of eight trainable parameters. All trainable parameters are constrained componentwise according to α+≥0,α−≥0,w1,kexp≥0,w2,kid≥0,w2,kexp≥0. _+≥ 0, _-≥ 0, w_1,k^exp≥ 0, w_2,k^id≥ 0, w_2,k^exp≥ 0. (48) For an incompressible deformation, the product of the principal stretches satisfies λ1λ2λ3=1 _1 _2 _3=1. Consequently, α≥3J_α≥ 3 for any real exponent α, and both shifted generalized invariants are non-negative. Together with the non-negative network weights, this ensures that every constitutive branch contributes non-negatively and that W≥0,W(1,1,1)=0.W≥ 0, W(1,1,1)=0. (49) Because the principal stretches enter through symmetric sums, the resulting strain energy density is objective and isotropic. Together with the incompressibility pressure correction, the undeformed configuration is stress-free. Noteworthy, the GI-CANN architecture admits a sufficient condition for polyconvexity in the cosidered incompressible case. For α+≥1 _+≥ 1, the generalized invariant α+=∑i=13σi()α+J_ _+= _i=1^3 _i( F) _+ (50) is a convex spectral function of F, where σi()=λi _i( F)= _i are its singular values. Under the incompressibility constraint, the inverse principal stretches coincide with the singular values of the cofactor. Consequently, −α−=∑i=13λi−α−=∑i=13σi(cof())α−J_- _-= _i=1^3 _i^- _-= _i=1^3 _i (cof( F) ) _- (51) is a convex spectral function of cof()cof( F) for α−≥1 _-≥ 1. Since the subsequent CANN mapping is convex and componentwise non-decreasing in both generalized invariants, the resulting strain energy density then admits a jointly convex representation in F and cof()cof( F), and is therefore polyconvex. In the present implementation, the exponent magnitudes are constrained only by α+,α−≥0 _+, _-≥ 0. Polyconvexity is therefore not guaranteed a priori if an active constitutive branch attains an exponent between zero and one. Imposing the stronger constraints α+≥1,α−≥1 _+≥ 1, _-≥ 1 (52) would provide a sufficient architectural guarantee of polyconvexity. These conditions only need to be satisfied by active constitutive branches; the exponent associated with a branch whose output weights vanish does not affect the resulting strain energy density. Hyperparameters As shown in Table 2, the hyperparameters of GI-CANNs are equivalent to those of the CANNs discussed earlier. 4 Results All computations use the same software implementation and parameter settings across the six methods. Detailed information on the execution command, number of repetitions and random seeds, runtime measurement procedure, and computational hardware is provided in Appendix A. 4.1 Predictive performance All six methods reproduce the experimental stress-stretch responses with high accuracy across uniaxial tension, equibiaxial tension, and pure shear, as depicted in Figs. 6 and 2. The predictions closely follow the experimental data over the complete stretch range, although the methods differ in their performance across deformation modes. The quantitative results in Table 4 and Figs. 4 and 5 show that all methods capture the dominant constitutive response. EUCLID provides the highest accuracy for uniaxial tension and pure shear, with R2=0.9985R^2=0.9985 and 0.99780.9978, respectively. GI-CANN performs particularly well for equibiaxial tension and reaches R2=0.9978R^2=0.9978, whereas EUCLID reaches R2=0.9930R^2=0.9930 for the same deformation mode. The differences between the methods become more apparent for equibiaxial tension, where MF, AMF, and PANN show larger deviations from the experimental response, see Fig. 6. In particular, PANN exhibits substantial variation across the different initializations. Overall, the results show that the methods differ less in their ability to reproduce the overall constitutive response than in their robustness and accuracy for individual deformation modes. EUCLID performs particularly well despite its sparse representation, while GI-CANN provides consistently high accuracy across all three loading modes. The results also demonstrate that good agreement in the training modes does not automatically lead to identical performance in the testing pure shear mode. 4.2 Best-performing models and model complexity We select one representative model for each method according to the minimum validation MSE in the testing pure shear mode. Table 5 summarizes the selected configurations, while Fig. 2 compares their predictions with the experimental data. AMF provides the highest validation accuracy among the selected models, with R2=0.9996R^2=0.9996. GI-CANN provides the highest training accuracy with R2=0.9996R^2=0.9996 and maintains a validation value of R2=0.9993R^2=0.9993. These results demonstrate that both methods can closely reproduce the experimental response, despite their fundamentally different model representations. The selected models also differ substantially in complexity. MF produces the sparsest model with only two nonzero parameters. EUCLID requires four nonzero parameters and still provides a validation R2R^2 of 0.99780.9978. GI-CANN uses eight nonzero parameters, whereas CANN, AMF, and PANN use 12, 49, and 371 nonzero parameters, respectively, see Table 5. PANN therefore requires substantially more parameters than the other methods without providing a corresponding increase in validation accuracy. These results reveal a clear accuracy-complexity trade-off. For Na=20N_a=20, AMF provides the highest selected-model validation accuracy but requires a substantially larger representation than EUCLID or GI-CANN. We note, however, that the sparsity of AMF can be optimized through a simple hyperparameter search (9). GI-CANN provides a particularly favorable compromise since it combines very high predictive accuracy with only eight nonzero parameters. EUCLID provides a different favorable compromise through its four-parameter representation and very low identification cost. 4.3 Variability and optimization behavior The repeated evaluations reveal clear differences in robustness between the methods. EUCLID and MF produce deterministic results for the selected settings, whereas AMF varies across the 20 investigated combinations of NaN_a and s. PANN, CANN, and GI-CANN vary across the 100 random initializations, see Tables 4, 4 and 5. For AMF, we note that we included Na=1N_a=1 for completeness, even though this reduces AMF to the classical MF with only one feature. PANN shows the largest variability, particularly for equibiaxial tension. Its mean R2R^2 reaches only 0.9039 for this deformation mode, with a standard deviation of 0.2450. In contrast, GI-CANN reaches R2=0.9978±0.0035R^2=0.9978± 0.0035 for the same deformation mode. CANN shows intermediate variability. The results, therefore, indicate substantially greater sensitivity to initialization for PANN than for GI-CANN. Fig. 6 illustrates this behavior through the prediction bands and training histories. PANN and CANN show pronounced differences between individual initializations, particularly for equibiaxial tension. GI-CANN also exhibits initialization-dependent convergence, but its prediction band remains comparatively narrow. This behavior suggests that the generalized invariant representation provides a more robust optimization landscape than the more flexible PANN architecture under the present benchmark conditions. EUCLID, MF, and AMF do not use iterative gradient-based training during the reported identification stage. Their variability, therefore, does not originate from random network initialization. Instead, AMF exhibits variation across its adaptive hyperparameters, while MF and EUCLID depend on their predefined model or feature spaces. 4.4 Computational cost of model identification The six methods show large differences in computational cost as shown in Tables 3 and 3. EUCLID and MF require less than 0.40.4 ms for their core identification procedures, with mean core times of 0.3370.337 ms and 0.3050.305 ms, respectively. AMF requires 1.0921.092 ms and therefore remains substantially faster than the neural-network methods. CANN, PANN, and GI-CANN require 324.539 ms, 588.495 ms, and 681.277 ms, respectively, for core optimization of 10,000 epochs. GI-CANN therefore requires more than three orders of magnitude more core computation than MF. The total identification time ranges from 0.4910.491 ms for EUCLID and 18.20118.201 ms for MF to approximately 0.40.4–0.80.8 s for the three neural-network methods, see Table 3. The results therefore separate two fundamentally different computational strategies. EUCLID, MF, and AMF perform sparse regression or database-based pattern recognition, whereas PANN, CANN, and GI-CANN solve nonlinear optimization problems. The latter approach provides greater flexibility in the constitutive representation but incurs a substantially higher identification cost. The MF timing requires particular consideration because the database generation occurs offline. The reported core time represents only the online pattern-recognition step. Once the database exists, the same database can support repeated material-characterization tasks without repeating the database-generation process. This feature makes MF particularly attractive for applications that require repeated characterization. The computational advantage can become even more pronounced in the unsupervised setting, where MF can process heterogeneous full-field measurements without solving a continuous optimization problem (12). The present benchmark therefore focuses on the supervised online cost and does not include the one-time database-generation effort. AMF follows the same offline database principle but adds an adaptive search over individual constitutive features. This additional search increases the online cost relative to MF, but AMF remains substantially faster than the neural-network methods while providing a more flexible constitutive representation. 4.5 Computational cost of constitutive model evaluation The cost of model identification is independent of the cost of constitutive evaluation during a finite element simulation. Such simulations repeatedly evaluate the strain energy density, stress, and consistent tangent at every integration point. We therefore benchmark the evaluation of W, P, and ℍH for the six selected models, as shown in Tables 6 and 8. For this purpose, we implemented all six methods within a consistent computational framework and employed JAX automatic differentiation to compute both stresses and tangent operators. We compare the evaluation times of the best-performing model for each method. MF provides the lowest evaluation time for all three quantities. It requires 0.151μs0.151\, per deformation gradient for W, 0.473μs0.473\, for P, and 3.263μs3.263\, for ℍH. The Hessian requires the largest computation for every method, with AMF reaching 7.750μs7.750\, compared with 3.263μs3.263\, for MF, see Table 6. The results show that the number of nonzero parameters alone does not determine the cost of constitutive evaluation. PANN uses 371 nonzero parameters, but its evaluation times remain comparable to those of the other neural-network models. MF combines the smallest model with the lowest evaluation cost. This result highlights the importance of considering both model structure and implementation when assessing constitutive-model efficiency. For EUCLID and the neural-network methods, we evaluate the complete selected model. We do not remove features or network components that become inactive during model identification. A reduced implementation could therefore further decrease the evaluation times, see Fig. 8. 4.6 Strain energy landscapes We finally compare the constitutive functions beyond the deformation paths used for calibration. Fig. 7 shows the strain energy density in principal-stretch space for the six selected models using identical contour levels and a common color scale. All models produce a smooth energy landscape with a minimum at the undeformed configuration (λ1,λ2)=(1,1)( _1, _2)=(1,1) and increasing energy with increasing deformation. Despite their different mathematical representations, the six models produce qualitatively similar energy landscapes. EUCLID, AMF, PANN, CANN, and GI-CANN show closely related contour shapes over the displayed deformation range. MF produces a more distinct landscape, with systematically shifted contours relative to the other models. This difference agrees with the larger deviations of MF in the equibiaxial tension response shown in Fig. 2. The direct contour overlay in Fig. 1 provides a more detailed comparison of the predicted energy levels. The contours remain closely aligned at low and moderate energy levels, whereas the differences increase toward larger stretches and higher energy levels. These regions extend beyond the central part of the experimental data and therefore require greater extrapolation from the identified constitutive models. Interestingly, Fig. 1 reveals that models based on principal invariants (EUCLID, CANN, PANN) yield similar strain energy densities, while models based on principal stretches (AMF, GI-CANN) likewise exhibit similar strain energy densities. This underscores the importance of selecting appropriate model input features. The individual landscapes in Fig. 7 further show that sparse representations such as EUCLID and GI-CANN can reproduce energy surfaces that closely resemble those of the more complex adaptive and neural-network representations. Despite their different mathematical representations, the six models produce qualitatively similar energy landscapes, see Figs. 1 and 7. EUCLID, PANN, and CANN exhibit nearly identical contour shapes across the displayed stretch range, while AMF and GI-CANN show similarly close agreement. MF produces shifted contours relative to the other models, which are consistent with its larger deviations in the equibiaxial tension response. Table 3: Computational performance of the six methods over 100 independent runs. Timing results are reported in milliseconds as mean ± standard deviation. The method-specific core time comprises the Coordinate Descent solve for EUCLID, pattern recognition for Material Fingerprinting, the adaptive loop for Adaptive Material Fingerprinting, and optimizer training after JAX warm-up for PANN, CANN, and GI-CANN. Setup time comprises the remaining preprocessing and initialization steps outside the core computation. Bold values indicate the best result within the corresponding column. Method Core time (ms) Setup time (ms) Total time (ms) EUCLID 0.337±0.1300.337± 0.130 0.137±0.0030.137± 0.003 0.491±0.1330.491± 0.133 MF 0.305±0.0060.305± 0.006 17.896±3.49417.896± 3.494 18.201±3.49318.201± 3.493 AMF 1.092±0.1401.092± 0.140 21.426±2.85421.426± 2.854 22.518±2.84022.518± 2.840 PANN 588.495±19.379588.495± 19.379 240.347±20.282240.347± 20.282 828.842±28.950828.842± 28.950 CANN 324.539±4.291324.539± 4.291 108.781±4.921108.781± 4.921 433.319±6.527433.319± 6.527 GI-CANN 681.277±14.110681.277± 14.110 117.587±12.135117.587± 12.135 798.863±18.259798.863± 18.259 Table 4: Predictive performance of the six methods for uniaxial tension (UT), equibiaxial tension (ET), and pure shear (PS). Results are reported as arithmetic mean ± standard deviation. The sample size n has a method-specific meaning: Adaptive Material Fingerprinting (AMF) comprises 20 combinations of the adaptive parameters NaN_a and s, and PANN, CANN, and GI-CANN each comprise 100 model initializations using different random seeds. MSE, RMSE, and MAE denote the mean squared, root mean squared, and mean absolute errors, respectively; NRMSE denotes the RMSE normalized by the range of the reference data. Bold values indicate the best mean performance within each loading mode. Method n MSE RMSE NRMSE MAE R2R^2 Uniaxial tension (UT) EUCLID 1 0.0025 ± 0 0.0501 ± 0 0.0112 ± 0 0.0389 ± 0 0.9985 ± 0 MF 1 0.0137 ± 0 0.1172 ± 0 0.0261 ± 0 0.0962 ± 0 0.9920 ± 0 AMF 20 0.0258 ± 0.0207 0.1443 ± 0.0722 0.0321 ± 0.0161 0.1033 ± 0.0588 0.9850 ± 0.0120 PANN 100 0.0103 ± 0.0156 0.0874 ± 0.0520 0.0195 ± 0.0116 0.0698 ± 0.0392 0.9940 ± 0.0091 CANN 100 0.0093 ± 0.0095 0.0887 ± 0.0385 0.0198 ± 0.0086 0.0725 ± 0.0295 0.9946 ± 0.0055 GI-CANN 100 0.0074 ± 0.0112 0.0619 ± 0.0601 0.0138 ± 0.0134 0.0507 ± 0.0502 0.9957 ± 0.0065 Equibiaxial tension (ET) EUCLID 1 0.0041 ± 0 0.0640 ± 0 0.0261 ± 0 0.0538 ± 0 0.9930 ± 0 MF 1 0.0407 ± 0 0.2016 ± 0 0.0823 ± 0 0.1725 ± 0 0.9310 ± 0 AMF 20 0.0355 ± 0.0407 0.1493 ± 0.1180 0.0610 ± 0.0482 0.0910 ± 0.0613 0.9397 ± 0.0690 PANN 100 0.0566 ± 0.1444 0.1682 ± 0.1692 0.0687 ± 0.0690 0.1219 ± 0.1008 0.9039 ± 0.2450 CANN 100 0.0150 ± 0.0099 0.1176 ± 0.0342 0.0480 ± 0.0140 0.0974 ± 0.0290 0.9746 ± 0.0168 GI-CANN 100 0.0013 ± 0.0020 0.0278 ± 0.0235 0.0114 ± 0.0096 0.0221 ± 0.0189 0.9978 ± 0.0035 Pure shear (PS) EUCLID 1 0.0008 ± 0 0.0277 ± 0 0.0152 ± 0 0.0221 ± 0 0.9978 ± 0 MF 1 0.0035 ± 0 0.0596 ± 0 0.0327 ± 0 0.0503 ± 0 0.9899 ± 0 AMF 20 0.0042 ± 0.0052 0.0504 ± 0.0420 0.0277 ± 0.0231 0.0378 ± 0.0308 0.9880 ± 0.0149 PANN 100 0.0037 ± 0.0055 0.0516 ± 0.0319 0.0283 ± 0.0175 0.0423 ± 0.0262 0.9895 ± 0.0158 CANN 100 0.0085 ± 0.0126 0.0805 ± 0.0450 0.0442 ± 0.0247 0.0671 ± 0.0409 0.9758 ± 0.0359 GI-CANN 100 0.0023 ± 0.0033 0.0374 ± 0.0308 0.0205 ± 0.0169 0.0293 ± 0.0272 0.9933 ± 0.0093 Table 5: Best-performing configuration identified for each method. Models were selected exclusively by the minimum validation mean squared error (MSE) in the pure-shear (PS) mode within the comparison run set; test data were not used for model selection. The table reports the principal hyperparameters, the random seed for neural network models, training and validation MSE, training and validation coefficient of determination R2R^2, the method-specific core computation time, and the number of nonzero parameters. Bold values denote the best predictive result or lowest computational cost in the corresponding column. Method Adaptive parameters Seed Training MSE Validation MSE Training R2R^2 Validation R2R^2 Core time (ms) # Parameters ‖0\| θ\|_0 EUCLID – – 3.304×10−33.304× 10^-3 7.663×10−47.663× 10^-4 0.9973 0.9978 0.372 4 MF – – 2.719×10−22.719× 10^-2 3.547×10−33.547× 10^-3 0.9781 0.9899 0.302 2 AMF Na=20N_a=20, s=0.25s=0.25 – 1.119×10−31.119× 10^-3 1.524×−1.524× 10^-4 0.9991 0.9996 6.000 49 PANN – 6 3.329×10−33.329× 10^-3 8.186×10−48.186× 10^-4 0.9973 0.9977 581.311 371 CANN – 93 3.535×10−33.535× 10^-3 8.925×10−48.925× 10^-4 0.9972 0.9975 320.947 12 GI-CANN – 39 4.436×−4.436× 10^-4 2.412×10−42.412× 10^-4 0.9996 0.9993 664.988 8 Table 6: Computational benchmark for evaluating the strain energy density function W, its gradient P, and its Hessian ℍH for the six methods. Times are reported in microseconds per sample as mean ± standard deviation over 100 repetitions after two warm-up runs. Each repetition evaluated 100 samples in a single chunk using double precision on the CPU. The reported times include one device synchronization per chunk and exclude compilation, parameter loading, sample generation, validation, and host result transfer. Bold values indicate the lowest mean evaluation time in the corresponding column. Method Energy W (μs /sample) Gradient P (μs /sample) Hessian ℍH (μs /sample) EUCLID 0.233±0.0440.233± 0.044 0.556±0.0390.556± 0.039 5.734±0.1365.734± 0.136 MF 0.151±0.0330.151± 0.033 0.473±0.0300.473± 0.030 3.263±0.1083.263± 0.108 AMF 0.430±0.0510.430± 0.051 1.262±0.0601.262± 0.060 7.750±0.3917.750± 0.391 PANN 0.428±0.0590.428± 0.059 0.827±0.0640.827± 0.064 5.496±0.2185.496± 0.218 CANN 0.175±0.0370.175± 0.037 0.548±0.0460.548± 0.046 4.578±0.0954.578± 0.095 GI-CANN 0.214±0.0360.214± 0.036 0.556±0.0470.556± 0.047 4.332±0.2724.332± 0.272 Figure 3: Method-specific core runtime of the six methods over 100 independent runs. Bars represent the arithmetic mean and error bars indicate one standard deviation. The core runtime comprises the Coordinate Descent solve for EUCLID, pattern recognition for Material Fingerprinting (MF), the adaptive loop for Adaptive Material Fingerprinting (AMF), and optimizer training after JAX warm-up for PANN, CANN, and GI-CANN. The runtime is reported in seconds on a logarithmic scale. Figure 4: Mean squared error (MSE) of the six methods for uniaxial tension (UT), equibiaxial tension (ET), and pure shear (PS). Bars represent the arithmetic mean and error bars indicate one standard deviation. EUCLID and Material Fingerprinting (MF) each contribute one deterministic result; Adaptive Material Fingerprinting (AMF) comprises 20 combinations of the adaptive parameters NaN_a and s; and PANN, CANN, and GI-CANN each comprise 100 model initializations using different random seeds. Lower values indicate better predictive performance. Figure 5: Coefficient of determination R2R^2 of the six methods for uniaxial tension (UT), equibiaxial tension (ET), and pure shear (PS). Bars represent the arithmetic mean and error bars indicate one standard deviation. EUCLID and Material Fingerprinting (MF) each contribute one deterministic result; Adaptive Material Fingerprinting (AMF) comprises 20 combinations of the adaptive parameters NaN_a and s; and PANN, CANN, and GI-CANN each comprise 100 model initializations using different random seeds. Values closer to one indicate better predictive performance. Figure 6: Predictive and training comparison of the six methods. Rows correspond to EUCLID, Material Fingerprinting (MF), Adaptive Material Fingerprinting (AMF), PANN, CANN, and GI-CANN. The first three columns show the first Piola–Kirchhoff stress component P11P_11 as a function of the stretch ratio λ for uniaxial tension (UT), equibiaxial tension (ET), and pure shear (PS), respectively. Black markers denote the experimental reference data, solid lines denote the mean predictions, and shaded regions indicate one standard deviation. For AMF, the variation reflects the investigated 20 hyperparameter configurations, whereas for PANN, CANN, and GI-CANN it reflects 100 model initializations using different random seeds. The fourth column shows the individual training loss histories on a logarithmic scale. No iterative loss history is available for EUCLID, MF, or AMF. Figure 7: Strain energy density W()W( F) predicted by the selected best-performing models: EUCLID, Material Fingerprinting (MF), Adaptive Material Fingerprinting (AMF), PANN, CANN, and GI-CANN. Each panel shows W as a function of the two independent principal stretches λ1 _1 and λ2 _2, with the third stretch determined by the incompressibility constraint λ3=(λ1λ2)−1 _3=( _1 _2)^-1. Identical contour levels and a common color scale are used in all six panels, enabling a direct comparison of the predicted energy landscapes. The marker at (λ1,λ2)=(1,1)( _1, _2)=(1,1) denotes the undeformed reference configuration. Figure 8: Runtime comparison for evaluating the strain energy density function W, its gradient P, and its Hessian ℍH for the six methods. Bars show the mean runtime per deformation gradient and error bars indicate one standard deviation over 100 repetitions after two warm-up runs. All evaluations were performed in double precision on the CPU using batches of 100 deformation gradients. The vertical axis uses a logarithmic scale. Compilation and other initialization steps are excluded from the reported runtimes. We note that for EUCLID and the neural network-based methods, we measured the evaluation times for the full model and did not exclude features that were identified as inactive during the inverse problem. Excluding such features could further reduce evauation times. 5 Conclusion We benchmark six data-driven constitutive modeling approaches on the classical Treloar rubber dataset under a common training and validation setup. AMF and GI-CANN provide the highest predictive accuracy among the selected models, while EUCLID and MF provide substantially lower identification costs and much sparser representations. GI-CANN combines high accuracy with a compact eight-parameter representation, whereas MF provides the lowest computational cost and the sparsest model. The benchmark therefore does not identify a single universally optimal method. Instead, the preferred approach depends on whether predictive accuracy, model sparsity, robustness, or computational efficiency receives the highest priority. The results further show that constitutive model evaluation cost can differ substantially from model identification cost and therefore represents an important consideration for computational deployment. Overall, the benchmark highlights complementary strengths across the considered approaches and provides a quantitative basis for selecting data-driven constitutive models according to the requirements of a specific application. Acknowledgments The authors acknowledge support from the European Research Council (ERC) Grant 101141626 DISCOVER funded by the European Union. Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Data and code availability All datasets, numerical results, and implementation files required to reproduce the analyses presented in this work have been deposited in a public Zenodo repository. The archive also contains the implementations of the six investigated approaches together with the corresponding training and evaluation scripts. The repository is available at https://doi.org/10.5281/zenodo.21915635. Disclosure of AI-assisted tools Generative artificial intelligence tools were used during the preparation of this manuscript to support language editing, improve clarity and readability, and assist with code development and debugging. All scientific content, methodological choices, numerical results, and conclusions were independently developed, verified, and approved by the authors, who take full responsibility for the content of the manuscript. References Abdusalamov et al. (2023) R. Abdusalamov, M. Hillgärtner, and M. Itskov Automatic generation of interpretable hyperelastic material models by symbolic regression. International Journal for Numerical Methods in Engineering, p. nme.7203 (en). External Links: ISSN 0029-5981, 1097-0207, Link, Document Cited by: §1. Amos et al. (2017) B. Amos, L. Xu, and J. Z. Kolter Input Convex Neural Networks. In Proceedings of the 34th International Conference on Machine Learning, Sydney, Australia, p. 10 (en). Cited by: §3.4. Anssari-Benam et al. (2024) A. Anssari-Benam, A. Goriely, and G. Saccomandi Generalised invariants and pseudo-universal relationships for hyperelastic materials: A new approach to constitutive modelling. Journal of the Mechanics and Physics of Solids 193, p. 105883 (en). External Links: ISSN 00225096, Link, Document Cited by: §3.6. Bleyer (2025) J. Bleyer Learning elastoplasticity with implicit layers. (en). External Links: Document Cited by: §1. Brunton et al. (2016) S. L. Brunton, J. L. Proctor, and J. N. Kutz Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences 113 (15), p. 3932–3937 (en). External Links: ISSN 0027-8424, 1091-6490, Link, Document Cited by: §3.1. Dammaß et al. (2025) F. Dammaß, K. A. Kalina, and M. Kästner When invariants matter: The role of I1 and I2 in neural network models of incompressible hyperelasticity. Mechanics of Materials, p. 105443. External Links: ISSN 0167-6636, Document Cited by: §1. Efron et al. (2004) B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani LEAST ANGLE REGRESSION. The Annals of Statistics 32 (2), p. 407–499 (en). Cited by: §3.1. Flaschel et al. (2026a) M. Flaschel, T. Hastie, and E. Kuhl Non-smooth optimization meets automated material model discovery. Journal of Computational Physics 563, p. 115043 (en). External Links: ISSN 00219991, Link, Document Cited by: §3.1. Flaschel et al. (2026b) M. Flaschel, H. Holthusen, D. Martonová, and E. Kuhl Adaptive Material Fingerprinting for the fast discovery of polyconvex feature combinations in isotropic and anisotropic hyperelasticity. International Journal of Engineering Science 227, p. 104594 (en). External Links: ISSN 00207225, Link, Document Cited by: §1, §3.3, §3.3, §3.3, §3.3, §4.2. Flaschel et al. (2021) M. Flaschel, S. Kumar, and L. De Lorenzis Unsupervised discovery of interpretable hyperelastic constitutive laws. Computer Methods in Applied Mechanics and Engineering 381, p. 113852 (en). External Links: Document Cited by: §3.1, §3.1, §3.1. Flaschel et al. (2026c) M. Flaschel, D. Martonová, C. Veil, and E. Kuhl Material Fingerprinting: A shortcut to material model discovery without solving optimization problems. Computer Methods in Applied Mechanics and Engineering 450, p. 118573 (en). External Links: ISSN 00457825, Link, Document Cited by: §1, §3.2, §3.2, §3.2. Flaschel et al. (2026d) M. Flaschel, M. A. Moreno-Mateos, S. Wiesheier, P. Steinmann, and E. Kuhl Unsupervised material fingerprinting: Ultra-fast hyperelastic model discovery from full-field experimental measurements. Computer Methods in Applied Mechanics and Engineering 461, p. 119256 (en). External Links: ISSN 00457825, Link, Document Cited by: §3.2, §3.2, §4.4. Flaschel et al. (2023) M. Flaschel, H. Yu, N. Reiter, J. Hinrichsen, S. Budday, P. Steinmann, S. Kumar, and L. De Lorenzis Automated discovery of interpretable hyperelastic material models for human brain tissue with EUCLID. Journal of the Mechanics and Physics of Solids 180, p. 105404 (en). External Links: ISSN 00225096, Link, Document Cited by: §1, §3.1, §3.1. Frank and Friedman (1993) I. E. Frank and J. H. Friedman A Statistical View of Some Chemometrics Regression Tools. Technometrics 35 (2), p. 109–135 (en). External Links: ISSN 0040-1706, 1537-2723, Link, Document Cited by: §3.1. Frankel et al. (2020) A. L. Frankel, R. E. Jones, and L. P. Swiler TENSOR BASIS GAUSSIAN PROCESS MODELS OF HYPERELASTIC MATERIALS. Journal of Machine Learning for Modeling and Computing 1 (1), p. 1–17 (en). External Links: ISSN 2689-3967, Link, Document Cited by: §1. Fu (1998) W. Fu Penalized Regressions: The Bridge versus the Lasso. (en). Cited by: §3.1. Fuhg et al. (2024) J. N. Fuhg, G. Anantha Padmanabha, N. Bouklas, B. Bahmani, W. Sun, N. N. Vlassis, M. Flaschel, P. Carrara, and L. De Lorenzis A Review on Data-Driven Constitutive Laws for Solids. Archives of Computational Methods in Engineering (en). External Links: ISSN 1134-3060, 1886-1784, Link, Document Cited by: §1. Fuhg et al. (2022) J. N. Fuhg, M. Marino, and N. Bouklas Local approximate Gaussian process regression for data-driven constitutive models: development and comparison with neural networks. Computer Methods in Applied Mechanics and Engineering 388, p. 114217 (en). External Links: ISSN 00457825, Link, Document Cited by: §1. Hastie et al. (2009) T. Hastie, R. Tibshirani, and J. Friedman The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Series in Statistics, Springer New York, New York, NY (en). External Links: ISBN 978-0-387-84857-0 978-0-387-84858-7, Link, Document Cited by: §3.1. Holthusen and Kuhl (2026) H. Holthusen and E. Kuhl A complement to neural networks for anisotropic inelasticity at finite strains. Computer Methods in Applied Mechanics and Engineering 450, p. 118612. External Links: ISSN 0045-7825, Document, Link Cited by: §3.5. Holthusen et al. (2024) H. Holthusen, L. Lamm, T. Brepols, S. Reese, and E. Kuhl Theory and implementation of inelastic Constitutive Artificial Neural Networks. Computer Methods in Applied Mechanics and Engineering 428, p. 117063 (en). External Links: ISSN 00457825, Link, Document Cited by: §3.5. Holthusen et al. (2026) H. Holthusen, K. Linka, E. Kuhl, and T. Brepols A generalized dual potential for inelastic constitutive artificial neural networks: a jax implementation at finite strains. Journal of the Mechanics and Physics of Solids 206, p. 106337. External Links: ISSN 0022-5096, Document, Link Cited by: §3.5. Klein et al. (2022) D. K. Klein, M. Fernández, R. J. Martin, P. Neff, and O. Weeger Polyconvex anisotropic hyperelasticity with neural networks. Journal of the Mechanics and Physics of Solids 159, p. 104703 (en). External Links: ISSN 00225096, Link, Document Cited by: §1, §3.4. Kopp (2008) J. Kopp Efficient numerical diagonalization of hermitian 3 × 3 matrices. International Journal of Modern Physics C 19 (03), p. 523–548. External Links: Document, Link Cited by: Appendix A. Linden et al. (2023) L. Linden, D. K. Klein, K. A. Kalina, J. Brummund, O. Weeger, and M. Kästner Neural networks meet hyperelasticity: A guide to enforcing physics. Journal of the Mechanics and Physics of Solids 179, p. 105363 (en). External Links: ISSN 00225096, Link, Document Cited by: §1, §3.4, §3.4. Linka et al. (2021) K. Linka, M. Hillgärtner, K. P. Abdolazizi, R. C. Aydin, M. Itskov, and C. J. Cyron Constitutive artificial neural networks: A fast and general approach to predictive data-driven constitutive modeling by deep learning. Journal of Computational Physics 429, p. 110010 (en). External Links: ISSN 00219991, Link, Document Cited by: §1, §3.5. Linka and Kuhl (2023) K. Linka and E. Kuhl A new family of Constitutive Artificial Neural Networks towards automated model discovery. Computer Methods in Applied Mechanics and Engineering 403, p. 115731 (en). External Links: ISSN 00457825, Link, Document Cited by: §1, §3.5, §3.6. Martonová et al. (2026a) D. Martonová, A. Goriely, and E. Kuhl Generalized invariants meet constitutive neural networks: A novel framework for hyperelastic materials. Journal of the Mechanics and Physics of Solids 206, p. 106352. External Links: ISSN 0022-5096, Document Cited by: §1, §3.6. Martonová et al. (2026b) D. Martonová, E. Kuhl, and M. Flaschel Material Fingerprinting for rapid discovery of hyperelastic models: First experimental validation. Journal of the Mechanics and Physics of Solids 208, p. 106463. External Links: ISSN 0022-5096, Document Cited by: §1, §3.2, §3.2. McCulloch et al. (2024) J. A. McCulloch, S. R. St. Pierre, K. Linka, and E. Kuhl On sparse regression, Lp_ p‐regularization, and automated model discovery. International Journal for Numerical Methods in Engineering 125 (14), p. e7481 (en). External Links: ISSN 0029-5981, 1097-0207, Link, Document Cited by: §3.1. Pierre et al. (2023) S. R. St. Pierre, K. Linka, and E. Kuhl Principal-stretch-based constitutive neural networks autonomously discover a subclass of Ogden models for human brain tissue. Brain Multiphysics 4, p. 100066. External Links: ISSN 2666-5220, Document Cited by: §3.6. Steinmann et al. (2012) P. Steinmann, M. Hossain, and G. Possart Hyperelastic models for rubber-like materials: consistent tangent operators and suitability for Treloar’s data. Archive of Applied Mechanics 82 (9), p. 1183–1217 (en). External Links: ISSN 0939-1533, 1432-0681, Link, Document Cited by: §1. Taç et al. (2024) V. Taç, K. Linka, F. Sahli-Costabal, E. Kuhl, and A. B. Tepole Benchmarking physics-informed frameworks for data-driven hyperelasticity. Computational Mechanics 73 (1), p. 49–65 (en). External Links: ISSN 0178-7675, 1432-0924, Link, Document Cited by: §1. Tac et al. (2022) V. Tac, F. Sahli Costabal, and A. B. Tepole Data-driven tissue mechanics with polyconvex neural ordinary differential equations. Computer Methods in Applied Mechanics and Engineering 398, p. 115248 (en). External Links: ISSN 00457825, Link, Document Cited by: §1. Tacke et al. (2025) M. Tacke, M. Busch, K. Bali, K. Abdolazizi, K. Linka, C. Cyron, and R. Aydin Constitutive scientific generative agent (CSGA): Leveraging large language models for automated constitutive model discovery. Machine Learning for Computational Science and Engineering 1 (1), p. 23 (en). External Links: ISSN 3005-1428, 3005-1436, Link, Document Cited by: §1. Tibshirani (1996) R. Tibshirani Regression Shrinkage and Selection via the Lasso. Journal of the Royal Statistical Society: Series B (Methodological) 58 (1), p. 267–288 (en). External Links: ISSN 00359246, Link, Document Cited by: §3.1. Treloar (1944) L. R. G. Treloar Stress-strain data for vulcanised rubber under various types of deformation. Transactions of the Faraday Society 40, p. 59 (en). External Links: ISSN 0014-7672, Link, Document Cited by: §1, §2.1. Wang et al. (2021) Z. Wang, J.B. Estrada, E.M. Arruda, and K. Garikipati Inference of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification. Journal of the Mechanics and Physics of Solids 153, p. 104474 (en). External Links: ISSN 00225096, Link, Document Cited by: §1. Wiesheier et al. (2024) S. Wiesheier, M. A. Moreno-Mateos, and P. Steinmann Versatile data-adaptive hyperelastic energy functions for soft materials. Computer Methods in Applied Mechanics and Engineering 430, p. 117208 (en). External Links: ISSN 00457825, Link, Document Cited by: §1. Appendix A Technical details For all results, we use the following command ⬇ python compare_all.py \ --n-timing 100 --n-seeds 100 --euclid-lambdas 1e-5 \ --study4-num-samples 100 --study4-repetitions 100 The option --n-timing 100 performs 100 independent runtime measurements for each of the six constitutive model identification methods. The option --n-seeds 100 evaluates the effect of neural-network initialization by training PANN, CANN, and GI-CANN with 100 different random seeds. The option --euclid-lambdas 1e-5 restricts the EUCLID hyperparameter study to the regularization parameter α=10−5α=10^-5. The option --study4-num-samples 100 generates one common set of 100 admissible incompressible deformation gradients and evaluates all six models for the strain energy density W, the first Piola stress =∂W∂ P= ∂ W∂ F, and the Hessian ℍ=∂2W∂∂H= ∂^2W∂ F∂ F. The option --study4-repetitions 100 repeats the evaluation of this complete set 100 times for each model–quantity combination. We divide the total runtime of each repetition by 100 to obtain the runtime per deformation gradient. We then calculate the mean, median, and sample standard deviation across the 100 repetitions. We exclude compilation and warm-up executions from these statistics and record the cold-start time separately. For the hyperparameters and method-specific configuration values, we refer to Table 7. We perform all computations on an Apple MacBook Pro with an Apple M4 processor that has a 10-core CPU with four performance cores and six efficiency cores and 24 GB of unified memory. The system runs macOS 26.6. We perform all reported runtime measurements on the CPU backend. Further, for the evaluation of the computational cost in Section 4.5, we compute the principal stretches as the singular values of the deformation gradient F. Instead of explicitly performing a singular-value decomposition, we exploit that the squared singular values of F are the eigenvalues of the right Cauchy–Green tensor =. C= F T F. (53) Denoting the eigenvalues of C by μi _i, the principal stretches follow as λi=μi,i=1,2,3. _i= _i, i=1,2,3. (54) Since C is symmetric, its eigenvalues are evaluated analytically using the closed-form trigonometric solution for a symmetric 3×33× 3 matrix (24). We first define q=13tr(),^=−q,q= 13tr( C), C= C-q I, (55) and p=16^:^.p= 16 C: C. (56) For p>0p>0, we introduce =^p,r=12det(),ϕ=13arccos(r), A= Cp, r= 12 ( A), φ= 13 (r), (57) where r is numerically restricted to the interval [−1,1][-1,1]. The largest and smallest eigenvalues of C are then obtained as μmax _ =q+2pcos(ϕ), =q+2p (φ), (58) μmin _ =q+2pcos(ϕ+2π3), =q+2p (φ+ 2π3 ), (59) while the remaining eigenvalue follows from the trace constraint as μmid=3q−μmax−μmin. _mid=3q- _ - _ . (60) The principal stretches are subsequently computed as λ1=μmax,λ2=μmid,λ3=μmin, _1= _ , _2= _mid, _3= _ , (61) such that λ1≥λ2≥λ3 _1≥ _2≥ _3. For nearly isotropic states, for which p falls below a machine-precision-dependent threshold, we instead use λ1=λ2=λ3=q, _1= _2= _3= q, (62) which avoids the numerical degeneracy of the trigonometric representation and ensures a well-defined evaluation of the undeformed reference configuration = F= I. Table 7: Fixed hyperparameters and method-specific configuration values. The configurations were selected by minimizing the validation mean squared error in the pure-shear (PS) mode without using test data. For PANN, CANN, and GI-CANN, the learning rate, AdamW optimizer settings, and number of epochs were fixed across the studies. Except for the method-specific learning rates, the default Optax settings of AdamW were used. The random seeds are therefore not listed as hyperparameters. For Adaptive Material Fingerprinting (AMF), the values shown correspond to the selected configuration, whereas the aggregate AMF performance statistics comprise 20 investigated combinations of the adaptive parameters NaN_a and s. Hyperparameter Value EUCLID Optimization method Coordinate Descent α 1×10−51× 10^-5 θth _th 0 Material Fingerprinting (MF) Material database Hyperelastic incompressible (HEI) Adaptive Material Fingerprinting (AMF) Material database Hyperelastic incompressible isotropic adaptive (HEIIA) NaN_a [1, 3, 5, 10, 20] s [0.25, 0.5, 0.75, 1.0] Physics-Augmented Neural Network (PANN) Optimizer AdamW (Optax) Weight constraints Projected Gradient Descent Learning rate 1×10−21× 10^-2 Remaining optimizer settings Optax defaults Number of epochs 10,000 Hidden layer widths [16,16][16,16] Output dimension 1 Hidden layer activation Softplus Output activation Softplus Constitutive Artificial Neural Network (CANN) Optimizer AdamW (Optax) Weight constraints Projected Gradient Descent Learning rate 1×10−41× 10^-4 Remaining optimizer settings Optax defaults Number of epochs 10,000 Generalized Invariant-based CANN (GI-CANN) Optimizer AdamW (Optax) Weight constraints Projected Gradient Descent Learning rate 1×10−31× 10^-3 Remaining optimizer settings Optax defaults Number of epochs 10,000