Benchmarking data-driven material models on the classic Treloar dataset

# Benchmarking data-driven material models on the classic Treloar dataset
Source: [https://arxiv.org/html/2608.14063](https://arxiv.org/html/2608.14063)
Hagen HolthusenCorresponding author:Correspondence: hagen\.holthusen@fau\.deAddress:Institute of Applied Mechanics, Egerlandstraße 5, Friedrich\-Alexander\-Universität Erlangen\-Nürnberg, 91058 Erlangen, GermanyMoritz FlaschelAddress:Institute of Applied Mechanics, Egerlandstraße 5, Friedrich\-Alexander\-Universität Erlangen\-Nürnberg, 91058 Erlangen, GermanyDenisa MartonováAddress:Institute of Applied Mechanics, Egerlandstraße 5, Friedrich\-Alexander\-Universität Erlangen\-Nürnberg, 91058 Erlangen, GermanyEllen KuhlAddress:Institute of Applied Mechanics, Egerlandstraße 5, Friedrich\-Alexander\-Universität Erlangen\-Nürnberg, 91058 Erlangen, GermanyAddress: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 ofTreloarto benchmark popular frameworks for hyperelasticity:\(Generalized\-Invariant\) Constitutive Artificial Neural Networks,Physics\-Augmented Neural Networks,\(Adaptive\) Material Fingerprinting, andEfficient 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](https://doi.org/10.5281/zenodo.21915635)\.

###### Keywords:

machine learning, automated discovery, material modeling, hyperelasticity

## 1Introduction

Treloar’s experiments on natural rubber\([37](https://arxiv.org/html/2608.14063#bib.bib9)\)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](https://arxiv.org/html/2608.14063#bib.bib23)\)\. 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[Table1](https://arxiv.org/html/2608.14063#S1.T1)\.[Fig\.1](https://arxiv.org/html/2608.14063#S1.F1)shows their discovered constitutive functions, and[Fig\.2](https://arxiv.org/html/2608.14063#S1.F2)shows their model response in comparison with Treloar’s data\.

Table 1:Benchmark performance of data\-driven methods \(see[Section4](https://arxiv.org/html/2608.14063#S4)and[Table5](https://arxiv.org/html/2608.14063#S4.T5)for details\)\.MethodAcronymTrainingR2R^\{2\}ValidationR2R^\{2\}Core time\(ms\)\# Parameters‖𝜽‖0\\\|\\bm\{\\theta\}\\\|\_\{0\}Efficient Unsupervised Constitutive Law Identification & DiscoveryEUCLID0\.99730\.99780\.3724Material FingerprintingMF0\.97810\.98990\.3022Adaptive Material FingerprintingAMF0\.99910\.99966\.00049Physics\-Augmented Neural NetworksPANN0\.99730\.9977581\.311371Constitutive Artificial Neural NetworksCANN0\.99720\.9975320\.94712Generalized Invariant\-based CANNGI\-CANN0\.99960\.9993664\.9888

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\.[Table1](https://arxiv.org/html/2608.14063#S1.T1)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[Section4](https://arxiv.org/html/2608.14063#S4)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](https://arxiv.org/html/2608.14063#bib.bib22)\), Material Fingerprinting \(MF\)\([11](https://arxiv.org/html/2608.14063#bib.bib28);[29](https://arxiv.org/html/2608.14063#bib.bib36)\), Adaptive Material Fingerprinting \(AMF\)\([9](https://arxiv.org/html/2608.14063#bib.bib31)\), Physics\-Augmented Neural Networks \(PANNs\)\([23](https://arxiv.org/html/2608.14063#bib.bib13);[25](https://arxiv.org/html/2608.14063#bib.bib15);[6](https://arxiv.org/html/2608.14063#bib.bib39)\), Constitutive Artificial Neural Networks \(CANNs\)\([26](https://arxiv.org/html/2608.14063#bib.bib7);[27](https://arxiv.org/html/2608.14063#bib.bib14)\), and Generalized Invariant\-based CANNs \(GI\-CANNs\)\([28](https://arxiv.org/html/2608.14063#bib.bib37)\)\. This list of methods is by no means exhaustive\. Numerous other approaches have been proposed in the literature\([17](https://arxiv.org/html/2608.14063#bib.bib16);[33](https://arxiv.org/html/2608.14063#bib.bib17)\), including, to name only a few, Gaussian processes\([15](https://arxiv.org/html/2608.14063#bib.bib12);[18](https://arxiv.org/html/2608.14063#bib.bib8)\), Variational System Identification \(VSI\)\([38](https://arxiv.org/html/2608.14063#bib.bib6)\), Neural Ordinary Differential Equations \(NODE\)\([34](https://arxiv.org/html/2608.14063#bib.bib4)\), symbolic regression\([1](https://arxiv.org/html/2608.14063#bib.bib24)\), spline\-based approaches\([39](https://arxiv.org/html/2608.14063#bib.bib18)\), non\-smooth parameterizations of constitutive functions\([4](https://arxiv.org/html/2608.14063#bib.bib27)\), and agentic artificial intelligence\([35](https://arxiv.org/html/2608.14063#bib.bib29)\)\.

In the following[Section2](https://arxiv.org/html/2608.14063#S2), we briefly review the classical Treloar dataset on natural rubber\. In the subsequent[Section3](https://arxiv.org/html/2608.14063#S3), 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[Section4](https://arxiv.org/html/2608.14063#S4)\.

## 2Benchmark 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\.1Dataset and preprocessing

In this work, we use the classical Treloar rubber dataset at20∘20^\{\\circ\}C as a benchmark for data\-driven constitutive model discovery\([37](https://arxiv.org/html/2608.14063#bib.bib9)\)\. 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λ\\lambdaand corresponding nominal first Piola\-Kirchhoff stress valuesP11P\_\{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𝑭\\bm\{F\}\. Assuming incompressibility, the deformation gradients are constructed as

𝑭UT=diag\(λ,λ−1/2,λ−1/2\),𝑭ET=diag\(λ,λ,λ−2\),𝑭PS=diag\(λ,1,λ−1\),\\bm\{F\}^\{\\mathrm\{UT\}\}=\\mathrm\{diag\}\\left\(\\lambda,\\lambda^\{\-1/2\},\\lambda^\{\-1/2\}\\right\),\\quad\\bm\{F\}^\{\\mathrm\{ET\}\}=\\mathrm\{diag\}\\left\(\\lambda,\\lambda,\\lambda^\{\-2\}\\right\),\\quad\\bm\{F\}^\{\\mathrm\{PS\}\}=\\mathrm\{diag\}\\left\(\\lambda,1,\\lambda^\{\-1\}\\right\),\(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](https://arxiv.org/html/2608.14063#S1.F2),

𝒟UT=\{\(λiUT,P11,iUT\)\|i=1,…,NUT\},𝒟ET=\{\(λiET,P11,iET\)\|i=1,…,NET\},𝒟PS=\{\(λiPS,P11,iPS\)\|i=1,…,NPS\}\.\\mathcal\{D\}^\{\\mathrm\{UT\}\}=\\left\\\{\(\\lambda\_\{i\}^\{\\mathrm\{UT\}\},P\_\{11,i\}^\{\\mathrm\{UT\}\}\)\\,\|\\,i=1,\\dots,N\_\{\\mathrm\{UT\}\}\\right\\\},\\quad\\mathcal\{D\}^\{\\mathrm\{ET\}\}=\\left\\\{\(\\lambda\_\{i\}^\{\\mathrm\{ET\}\},P\_\{11,i\}^\{\\mathrm\{ET\}\}\)\\,\|\\,i=1,\\dots,N\_\{\\mathrm\{ET\}\}\\right\\\},\\quad\\mathcal\{D\}^\{\\mathrm\{PS\}\}=\\left\\\{\(\\lambda\_\{i\}^\{\\mathrm\{PS\}\},P\_\{11,i\}^\{\\mathrm\{PS\}\}\)\\,\|\\,i=1,\\dots,N\_\{\\mathrm\{PS\}\}\\right\\\}\.\(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\.\\mathcal\{D\}\_\{\\mathrm\{train\}\}=\\mathcal\{D\}^\{\\mathrm\{UT\}\}\\cup\\mathcal\{D\}^\{\\mathrm\{ET\}\},\\quad\\mathcal\{D\}\_\{\\mathrm\{test\}\}=\\mathcal\{D\}^\{\\mathrm\{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\|\.\\widehat\{P\}\_\{11,i\}=\\frac\{P\_\{11,i\}\}\{\\max\\limits\_\{j=1,\\dots,N\_\{\\mathrm\{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\.2Training 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 functionW𝜽W\_\{\\bm\{\\theta\}\}parameterized by the model parameters𝜽\\bm\{\\theta\}\. 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\\bm\{F\}\_\{i\}in the training dataset, the model predicts the first Piola–Kirchhoff stress from the learned strain energy density

𝑷𝜽​\(𝑭i\)=∂W𝜽​\(𝑭i\)∂𝑭i\.\\bm\{P\}\_\{\\bm\{\\theta\}\}\(\\bm\{F\}\_\{i\}\)=\\frac\{\\partial W\_\{\\bm\{\\theta\}\}\(\\bm\{F\}\_\{i\}\)\}\{\\partial\\bm\{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 conditionP33=0P\_\{33\}=0which 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 normalizedP11P\_\{11\}component\. The training objective is the mean squared error

ℒdata​\(𝜽\)=1Ntrain​∑i=1Ntrain\(P^11,𝜽​\(𝑭i\)−P^11,i\)2,\\mathcal\{L\}\_\{\\mathrm\{data\}\}\(\\bm\{\\theta\}\)=\\frac\{1\}\{N\_\{\\mathrm\{train\}\}\}\\sum\_\{i=1\}^\{N\_\{\\mathrm\{train\}\}\}\\left\(\\widehat\{P\}\_\{11,\\bm\{\\theta\}\}\(\\bm\{F\}\_\{i\}\)\-\\widehat\{P\}\_\{11,i\}\\right\)^\{2\},\(6\)whereP^11,i\\widehat\{P\}\_\{11,i\}denotes the normalized experimental stress andP^11,𝜽\\widehat\{P\}\_\{11,\\bm\{\\theta\}\}denotes the corresponding model prediction normalized by the same factor as the experimental stress\.

The total optimization problem is given by

𝜽∗=arg⁡min𝜽​ℒdata​\(𝜽\),\\bm\{\\theta\}^\{\\ast\}=\\arg\\min\_\{\\bm\{\\theta\}\}\\mathcal\{L\}\_\{\\mathrm\{data\}\}\(\\bm\{\\theta\}\),\(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\.3Evaluation 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𝒟\\mathcal\{D\}with a number of\|𝒟\|\|\\mathcal\{D\}\|data points, we consider the mean squared error \(MSE\)

MSE𝒟=1\|𝒟\|​∑i=1\|𝒟\|\(P11,𝜽​\(λi\)−P11,i\)2,\\text\{MSE\}\_\{\\mathcal\{D\}\}=\\frac\{1\}\{\|\\mathcal\{D\}\|\}\\sum\_\{i=1\}^\{\|\\mathcal\{D\}\|\}\\left\(\{P\}\_\{11,\\bm\{\\theta\}\}\(\\lambda\_\{i\}\)\-\{P\}\_\{11,i\}\\right\)^\{2\},\(8\)the root mean squared error \(RMSE\)

RMSE𝒟=1\|𝒟\|​∑i=1\|𝒟\|\(P11,𝜽​\(λi\)−P11,i\)2,\\text\{RMSE\}\_\{\\mathcal\{D\}\}=\\sqrt\{\\frac\{1\}\{\|\\mathcal\{D\}\|\}\\sum\_\{i=1\}^\{\|\\mathcal\{D\}\|\}\\left\(\{P\}\_\{11,\\bm\{\\theta\}\}\(\\lambda\_\{i\}\)\-\{P\}\_\{11,i\}\\right\)^\{2\}\},\(9\)the range\-normalized root mean squared error \(NRMSE\)

NRMSE𝒟=RMSE𝒟maxj=1,…,\|𝒟\|⁡\|P11,j\|−minj=1,…,\|𝒟\|⁡\|P11,j\|,\\text\{NRMSE\}\_\{\\mathcal\{D\}\}=\\frac\{\\text\{RMSE\}\_\{\\mathcal\{D\}\}\}\{\\max\\limits\_\{j=1,\\dots,\|\\mathcal\{D\}\|\}\|P\_\{11,j\}\|\-\\min\\limits\_\{j=1,\\dots,\|\\mathcal\{D\}\|\}\|P\_\{11,j\}\|\},\(10\)the mean absolute error

MAE𝒟=1\|𝒟\|​∑i=1\|𝒟\|\|P11,𝜽​\(λi\)−P11,i\|,\\text\{MAE\}\_\{\\mathcal\{D\}\}=\\frac\{1\}\{\|\\mathcal\{D\}\|\}\\sum\_\{i=1\}^\{\|\\mathcal\{D\}\|\}\\left\|\{P\}\_\{11,\\bm\{\\theta\}\}\(\\lambda\_\{i\}\)\-\{P\}\_\{11,i\}\\right\|,\(11\)and the coefficient of determination

R𝒟2=1−∑i=1\|𝒟\|\(P11,𝜽​\(λi\)−P11,i\)2∑i=1\|𝒟\|\(mean𝒟​\(P11,i\)−P11,i\)2,R^\{2\}\_\{\\mathcal\{D\}\}\\;=\\;1\-\\frac\{\\sum\_\{i=1\}^\{\|\\mathcal\{D\}\|\}\\left\(\{P\}\_\{11,\\bm\{\\theta\}\}\(\\lambda\_\{i\}\)\-\{P\}\_\{11,i\}\\right\)^\{2\}\}\{\\sum\_\{i=1\}^\{\|\\mathcal\{D\}\|\}\\left\(\\text\{mean\}\_\{\\mathcal\{D\}\}\(\{P\}\_\{11,i\}\)\-\{P\}\_\{11,i\}\\right\)^\{2\}\},\(12\)wheremean𝒟​\(P11,i\)\\text\{mean\}\_\{\\mathcal\{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\\\|\\bm\{\\theta\}\\\|\_\{0\}, in analogy with theL0L\_\{0\}\-pseudo\-norm, which counts the number of nonzero entries in a vector\.

While the number of material parameters‖𝜽‖0\\\|\\bm\{\\theta\}\\\|\_\{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 energyWW, its gradient𝑷\\bm\{P\}, and its Hessianℍ\\mathbb\{H\}\. Each benchmark evaluated a batch ofNeval=100N\_\{\\mathrm\{eval\}\}=100deformation 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 byNevalN\_\{\\mathrm\{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±\\pmstandard deviation over the 100 repetitions\.

## 3Methods

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[Table2](https://arxiv.org/html/2608.14063#S3.T2), 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\.MethodSolvers / AlgorithmsHyperparametersModeling DecisionsEUCLID∙\\bulletCoordinate Descent \(CD\)
∙\\bulletLeast Angle Regression \(LARS\)∙\\bulletNumber of terms in the library
∙\\bulletRegularization parameterα\\alpha
∙\\bulletThresholding parameterθ→0\\theta\_\{\{\\mathrel\{\\scalebox\{0\.5\}\[1\]\{$\\rightarrow$\}\}\}\{0\}\}∙\\bulletChoice of the library features
∙\\bulletModel input variablesMaterial Fingerprinting∙\\bulletPattern\-recognition∙\\bulletNumber of fingerprints in the database
∙\\bulletParameter sampling during database generation∙\\bulletChoice of the models in the database
∙\\bulletModel input variablesAdaptive Material Fingerprinting∙\\bulletAdaptive pattern\-recognition∙\\bulletNumber of fingerprints in the database
∙\\bulletParameter sampling during database generation
∙\\bulletNumber of terms in the modelNaN\_\{a\}
∙\\bulletStep sizes∈\[0,1\]s\\in\[0,1\]∙\\bulletChoice of the modeling features in the database
∙\\bulletModel input variablesPANNVarious \(usually gradient\-based\)
∙\\bulletSGD
∙\\bulletAdam
∙\\bulletL\-BFGS∙\\bulletLearning rateη\\eta
∙\\bulletNumber of epochs
∙\\bulletNumber of layers
∙\\bulletNumber of neurons per layer
∙\\bulletOptional ML hyperparameters \(e\.g\., regularization, early stopping, dropout, scheduling\)∙\\bulletChoice of the activation functions
∙\\bulletModel input variablesCANN and GI\-CANNVarious \(usually gradient\-based\)
∙\\bulletSGD
∙\\bulletAdam
∙\\bulletL\-BFGS∙\\bulletLearning rateη\\eta
∙\\bulletNumber of epochs
∙\\bulletNumber of layers
∙\\bulletNumber of neurons per layer
∙\\bulletOptional ML hyperparameters \(e\.g\., regularization, early stopping, dropout, scheduling\)∙\\bulletChoice of the activation functions
∙\\bulletModel input variablesFor 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\],\\tilde\{W\}\(\\bm\{F\}\)=W\(\\bm\{F\}\)\-p\[J\-1\],\(13\)whereJ=det\(𝑭\)J=\\det\(\\bm\{F\}\)andppis a scalar Lagrange multiplier enforcing the incompressibility constraintJ=1J=1\. To ensure material isotropy, the constitutive contributionW⁡\(𝑭\)W\(\\bm\{F\}\)is commonly expressed in terms of invariants of the right Cauchy–Green deformation tensor𝑪=𝑭T​𝑭\\bm\{C\}=\\bm\{F\}^\{T\}\\bm\{F\}\. The most common choices are the first and second principal invariants,I1=tr⁡\(𝑪\)I\_\{1\}=\\operatorname\{tr\}\(\\bm\{C\}\)andI2=12​\[\(tr⁡𝑪\)2−tr⁡\(𝑪2\)\]I\_\{2\}=\\frac\{1\}\{2\}\\left\[\(\\operatorname\{tr\}\\bm\{C\}\)^\{2\}\-\\operatorname\{tr\}\(\\bm\{C\}^\{2\}\)\\right\], or the principal stretchesλ1,λ2,λ3\\lambda\_\{1\},\\lambda\_\{2\},\\lambda\_\{3\}, whose squares are the eigenvalues of𝑪\\bm\{C\}, i\.e\.,\{λ12,λ22,λ32\}=eig⁡\(𝑪\)\\\{\\lambda\_\{1\}^\{2\},\\lambda\_\{2\}^\{2\},\\lambda\_\{3\}^\{2\}\\\}=\\operatorname\{eig\}\(\\bm\{C\}\)\. These quantities satisfyI1=λ12\+λ22\+λ32I\_\{1\}=\\lambda\_\{1\}^\{2\}\+\\lambda\_\{2\}^\{2\}\+\\lambda\_\{3\}^\{2\},I2=λ12​λ22\+λ22​λ32\+λ32​λ12I\_\{2\}=\\lambda\_\{1\}^\{2\}\\lambda\_\{2\}^\{2\}\+\\lambda\_\{2\}^\{2\}\\lambda\_\{3\}^\{2\}\+\\lambda\_\{3\}^\{2\}\\lambda\_\{1\}^\{2\}, andJ=λ1​λ2​λ3=1J=\\lambda\_\{1\}\\lambda\_\{2\}\\lambda\_\{3\}=1\. Alternatively, generalized invariant sets may be employed, as discussed below for the GI\-CANN framework\.

By differentiatingW~\\tilde\{W\}with respect to the deformation gradient and determiningppfrom the constraintP33=0P\_\{33\}=0for the dataset introduced in[Section2](https://arxiv.org/html/2608.14063#S2), we obtain theP11P\_\{11\}component of the first Piola–Kirchhoff stress as

P11=∂W~∂F11=∂W∂F11−∂W∂F33​\[cof⁡\(𝑭\)33\]−1​cof⁡\(𝑭\)11\.P\_\{11\}=\\frac\{\\partial\\tilde\{W\}\}\{\\partial F\_\{11\}\}=\\frac\{\\partial W\}\{\\partial F\_\{11\}\}\-\\frac\{\\partial W\}\{\\partial F\_\{33\}\}\\left\[\\operatorname\{cof\}\(\\bm\{F\}\)\_\{33\}\\right\]^\{\-1\}\\operatorname\{cof\}\(\\bm\{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 contributionWW\.

### 3\.1Efficient 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](https://arxiv.org/html/2608.14063#bib.bib5)\)\. 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](https://arxiv.org/html/2608.14063#bib.bib5)\), 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](https://arxiv.org/html/2608.14063#bib.bib22)\), 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\.\\bm\{Q\}\(I\_\{1\},I\_\{2\}\)=\\left\[\[I\_\{1\}\-3\],\\,\[I\_\{2\}\-3\],\\,\[I\_\{1\}\-3\]^\{2\},\\,\[I\_\{1\}\-3\]\[I\_\{2\}\-3\],\\,\[I\_\{2\}\-3\]^\{2\},\\,\\dots\\right\]^\{T\}\.\(15\)In this work, we consider coupled polynomial features up to an order of five, resulting inNfeat=20N\_\{\\text\{feat\}\}=20modeling features\. We note, however, that the library can be readily extended with more modeling features, such as those of the Ogden material model\([13](https://arxiv.org/html/2608.14063#bib.bib22)\)\.

A linear combination of the candidate modeling terms yields the strain energy density

W⁡\(I1,I2\)=𝑸⁡\(I1,I2\)⋅𝜽,W\(I\_\{1\},I\_\{2\}\)=\\bm\{Q\}\(I\_\{1\},I\_\{2\}\)\\cdot\\bm\{\\theta\},\(16\)where𝜽∈ℝ≥0Nfeat\\bm\{\\theta\}\\in\\mathbb\{R\}\_\{\\geq 0\}^\{N\_\{\\text\{feat\}\}\}are material parameters, which we assume here to be non\-negative\.

The objective of EUCLID is to identify the parameters𝜽\\bm\{\\theta\}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

𝜽∗=arg⁡min𝜽∈ℝ≥0​ℒdata​\(𝜽\)\+α​‖𝜽‖1,\\bm\{\\theta\}^\{\\ast\}=\\arg\\min\_\{\\bm\{\\theta\}\\in\\mathbb\{R\}\_\{\\geq 0\}\}\\mathcal\{L\}\_\{\\mathrm\{data\}\}\(\\bm\{\\theta\}\)\+\\alpha\\\|\\bm\{\\theta\}\\\|\_\{1\},\(17\)in which theL1L\_\{1\}\-norm penalizes solutions with many non\-zero parameters and promotes sparsity in the solution\([14](https://arxiv.org/html/2608.14063#bib.bib3);[36](https://arxiv.org/html/2608.14063#bib.bib10);[19](https://arxiv.org/html/2608.14063#bib.bib21);[5](https://arxiv.org/html/2608.14063#bib.bib11)\)\. The regularization parameterα≥0\\alpha\\geq 0determines how aggressively the parameters should be pushed to zero\. Since our model is linear in the parameters𝜽\\bm\{\\theta\}\(see[Eq\.16](https://arxiv.org/html/2608.14063#S3.E16)\), the loss term in the optimization problem above is quadratic in𝜽\\bm\{\\theta\}\. TheL1L\_\{1\}\-norm regularization term is non\-smooth but convex in𝜽\\bm\{\\theta\}\. The total objective function is convex and non\-smooth and can be efficiently minimized using methods such as Coordinate Descent \(CD\)\([16](https://arxiv.org/html/2608.14063#bib.bib25)\)or Least Angle Regression \(LARS\)\([7](https://arxiv.org/html/2608.14063#bib.bib2)\)\. Such non\-smooth optimization approaches have also been applied to constitutive model discovery and material modeling\([8](https://arxiv.org/html/2608.14063#bib.bib30)\)\.

Solving the problem above yields a parameter vector𝜽\\bm\{\\theta\}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\\theta\_\{\{\\mathrel\{\\scalebox\{0\.5\}\[1\]\{$\\rightarrow$\}\}\}\{0\}\}can be introduced in a post\-processing step, whereby parameters with absolute values below the threshold are set to zero\([10](https://arxiv.org/html/2608.14063#bib.bib5)\)\. 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\\alpha=0to 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\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), EUCLID depends on several hyperparameters, including the number of terms in the candidate library, the regularization parameterα\\alpha, and the threshold parameterθ→0\\theta\_\{\{\\mathrel\{\\scalebox\{0\.5\}\[1\]\{$\\rightarrow$\}\}\}\{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 nonconvexLpL\_\{p\}\-regularization term with0<p<10<p<1\([10](https://arxiv.org/html/2608.14063#bib.bib5);[30](https://arxiv.org/html/2608.14063#bib.bib20)\)\. However, since increasing the regularization parameterα\\alphaallows theL1L\_\{1\}\-norm to achieve any desired level of sparsity, we restrict our attention here to the convexL1L\_\{1\}\-regularization\.

### 3\.2Material 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](https://arxiv.org/html/2608.14063#bib.bib28);[29](https://arxiv.org/html/2608.14063#bib.bib36)\)\. 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 functionWWand sampled over a prescribed range of material parameters\([29](https://arxiv.org/html/2608.14063#bib.bib36)\)\. 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](https://arxiv.org/html/2608.14063#bib.bib32)\)\. Each database entry stores a fingerprint together with its constitutive model and material parameters

\(𝐟\(i\),𝜽\(i\)\),i=1,…,Ndb,\\left\(\\mathbf\{f\}^\{\(i\)\},\\bm\{\\theta\}^\{\(i\)\}\\right\),\\qquad i=1,\\ldots,N\_\{\\mathrm\{db\}\},\(18\)whereNdbN\_\{\\mathrm\{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\)‖,\\bar\{\\mathbf\{f\}\}^\{\(i\)\}=\\frac\{\\mathbf\{f\}^\{\(i\)\}\}\{\\\|\\mathbf\{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](https://arxiv.org/html/2608.14063#bib.bib28)\)\.

During the online stage, we measure and normalize the fingerprint𝐟∗\\mathbf\{f\}^\{\\ast\}of an unknown material\. We then compare the measured fingerprint with every database entry through the cosine similarity

cs\(i\)=𝐟¯\(i\)⋅𝐟¯∗\|𝐟¯\(i\)\|​\|𝐟¯∗\|=𝐟¯\(i\)⋅𝐟¯∗,\\text\{cs\}^\{\(i\)\}=\\frac\{\\bar\{\\mathbf\{f\}\}^\{\(i\)\}\\cdot\\bar\{\\mathbf\{f\}\}^\{\\ast\}\}\{\\left\|\\bar\{\\mathbf\{f\}\}^\{\(i\)\}\\right\|\\left\|\\bar\{\\mathbf\{f\}\}^\{\\ast\}\\right\|\}=\\bar\{\\mathbf\{f\}\}^\{\(i\)\}\\cdot\\bar\{\\mathbf\{f\}\}^\{\\ast\},\(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∗=arg⁡maxi⁡cs\(i\)\.i^\{\\ast\}=\\arg\\max\_\{i\}\\;\\text\{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](https://arxiv.org/html/2608.14063#bib.bib28)\)\. 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](https://arxiv.org/html/2608.14063#bib.bib32)\)\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), 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\.3Adaptive 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](https://arxiv.org/html/2608.14063#bib.bib31)\)\. 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θa​Qa​\(𝑭\),W\(\\bm\{F\}\)=\\sum\_\{a=1\}^\{N\_\{a\}\}\\theta\_\{a\}\\,Q\_\{a\}\(\\bm\{F\}\),\(22\)whereNaN\_\{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](https://arxiv.org/html/2608.14063#bib.bib31)\)\. 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=𝟎\.\\mathbf\{f\}\_\{0\}=\\mathbf\{0\}\.\(23\)At iterationkk, we compute the residual fingerprint

𝐫k=𝐟∗−𝐟k,\\mathbf\{r\}\_\{k\}=\\mathbf\{f\}^\{\\ast\}\-\\mathbf\{f\}\_\{k\},\(24\)where𝐟∗\\mathbf\{f\}^\{\\ast\}denotes the experimental fingerprint and𝐟k\\mathbf\{f\}\_\{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∗=arg⁡maxi⁡𝐟¯\(i\)⋅𝐫¯k\|𝐟¯\(i\)\|​\|𝐫¯k\|\.i\_\{k\}^\{\\ast\}=\\arg\\max\_\{i\}\\;\\frac\{\\bar\{\\mathbf\{f\}\}^\{\(i\)\}\\cdot\\bar\{\\mathbf\{r\}\}\_\{k\}\}\{\\left\|\\bar\{\\mathbf\{f\}\}^\{\(i\)\}\\right\|\\left\|\\bar\{\\mathbf\{r\}\}\_\{k\}\\right\|\}\.\(25\)
After we identify the best matching feature, we add it to the constitutive model\. Its contribution is scaled by a step size parameters∈\(0,1\]s\\in\(0,1\], which controls how much of the current residual is incorporated at each iteration\. Thus, smaller values ofsslead to more gradual updates, whereass=1s=1corresponds to incorporating the full residual contribution\. We repeat the procedure for a prescribed maximum number of iterationsNaN\_\{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 ofNaN\_\{a\}, produces sparse, physically interpretable constitutive models without solving a continuous optimization problem\([9](https://arxiv.org/html/2608.14063#bib.bib31)\)\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), Adaptive Material Fingerprinting depends on a larger number of hyperparameters than the original Material Fingerprinting method\. In particular, the number of termsNaN\_\{a\}included in the model and the step sizessmust 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](https://arxiv.org/html/2608.14063#bib.bib31)\)\.

### 3\.4Physics\-Augmented Neural Networks \(PANNs\)

Physics\-Augmented Neural Networks\([23](https://arxiv.org/html/2608.14063#bib.bib13);[25](https://arxiv.org/html/2608.14063#bib.bib15)\)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](https://arxiv.org/html/2608.14063#bib.bib1)\)that takes the first and second principal invariants as inputs,

W⁡\(I1,I2\)=ICNN𝜽⁡\(I1,I2\)−ICNN𝜽⁡\(3,3\),W\(I\_\{1\},I\_\{2\}\)=\\operatorname\{ICNN\}\_\{\\bm\{\\theta\}\}\(I\_\{1\},I\_\{2\}\)\-\\operatorname\{ICNN\}\_\{\\bm\{\\theta\}\}\(3,3\),\(26\)where𝜽\\bm\{\\theta\}denotes the trainable network parameters\. SinceI1=I2=3I\_\{1\}=I\_\{2\}=3in the undeformed reference configuration, the second term normalizes the strain energy density such thatW⁡\(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](https://arxiv.org/html/2608.14063#bib.bib15)is not required at this point, as the Lagrange multiplierppensures 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\}=\\bm\{F\}:\\bm\{F\},\\qquad I\_\{2\}=\\operatorname\{cof\}\(\\bm\{F\}\):\\operatorname\{cof\}\(\\bm\{F\}\)\.\(27\)Thus,I1I\_\{1\}is a convex function of𝑭\\bm\{F\}, whileI2I\_\{2\}is a convex function ofcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\)\. IfW⁡\(I1,I2\)W\(I\_\{1\},I\_\{2\}\)is convex and non\-decreasing in each argument, the composition is jointly convex in𝑭\\bm\{F\}andcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\), regarded as independent arguments\. Consequently, a strain energy density of the form

W\(𝑭\)=ICNN𝜽\(𝑭:𝑭,cof\(𝑭\):cof\(𝑭\)\)−ICNN𝜽\(3,3\)W\(\\bm\{F\}\)=\\operatorname\{ICNN\}\_\{\\bm\{\\theta\}\}\\left\(\\bm\{F\}:\\bm\{F\},\\,\\operatorname\{cof\}\(\\bm\{F\}\):\\operatorname\{cof\}\(\\bm\{F\}\)\\right\)\-\\operatorname\{ICNN\}\_\{\\bm\{\\theta\}\}\(3,3\)\(28\)is polyconvex if the ICNN is convex and componentwise non\-decreasing in its inputs\. For incompressible deformations,I1≥3I\_\{1\}\\geq 3andI2≥3I\_\{2\}\\geq 3\. Componentwise monotonicity, together with the normalization at\(I1,I2\)=\(3,3\)\(I\_\{1\},I\_\{2\}\)=\(3,3\), therefore additionally ensuresW≥0W\\geq 0for 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\.\\bm\{x\}=\\begin\{bmatrix\}I\_\{1\},&I\_\{2\}\\end\{bmatrix\}^\{T\}\.\(29\)The first hidden representation is defined as

𝒛1=φ1​\(𝑼1​𝒙\+𝒃1\),\\bm\{z\}\_\{1\}=\\varphi\_\{1\}\\left\(\\bm\{U\}\_\{1\}\\bm\{x\}\+\\bm\{b\}\_\{1\}\\right\),\(30\)and the subsequent representations are constructed recursively as

𝒛ℓ=φℓ\(𝑾ℓ𝒛ℓ−1\+𝑼ℓ𝒙\+𝒃ℓ\),ℓ=2,…,L\.\\bm\{z\}\_\{\\ell\}=\\varphi\_\{\\ell\}\\left\(\\bm\{W\}\_\{\\ell\}\\bm\{z\}\_\{\\ell\-1\}\+\\bm\{U\}\_\{\\ell\}\\bm\{x\}\+\\bm\{b\}\_\{\\ell\}\\right\),\\qquad\\ell=2,\\ldots,L\.\(31\)The final representation is scalar and defines the network output,

ICNN𝜽⁡\(𝒙\)=𝒛L\.\\operatorname\{ICNN\}\_\{\\bm\{\\theta\}\}\(\\bm\{x\}\)=\\bm\{z\}\_\{L\}\.\(32\)The ICNN can therefore be understood as a composition of the layer mappings

gℓ​\(𝒛,𝒙\)=φℓ​\(𝑾ℓ​𝒛\+𝑼ℓ​𝒙\+𝒃ℓ\)\.g\_\{\\ell\}\(\\bm\{z\};\\bm\{x\}\)=\\varphi\_\{\\ell\}\\left\(\\bm\{W\}\_\{\\ell\}\\bm\{z\}\+\\bm\{U\}\_\{\\ell\}\\bm\{x\}\+\\bm\{b\}\_\{\\ell\}\\right\)\.\(33\)Importantly, each layer depends not only on the preceding hidden representation but also directly on the original input through the term𝑼ℓ​𝒙\\bm\{U\}\_\{\\ell\}\\bm\{x\}\. These direct input paths constitute the skip connections of the ICNN architecture\.

Convexity of the ICNN with respect to𝒙\\bm\{x\}is ensured by choosing activation functionsφℓ\\varphi\_\{\\ell\}that are convex and non\-decreasing and by imposing the componentwise constraints

𝑾ℓ≥𝟎,ℓ=2,…,L\.\\bm\{W\}\_\{\\ell\}\\geq\\bm\{0\},\\qquad\\ell=2,\\ldots,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𝑼ℓ​𝒙\\bm\{U\}\_\{\\ell\}\\bm\{x\}is affine in𝒙\\bm\{x\}\. Input convexity alone would therefore not require a sign constraint on𝑼ℓ\\bm\{U\}\_\{\\ell\}\. For the polyconvex construction above, however, the ICNN must additionally be non\-decreasing in bothI1I\_\{1\}andI2I\_\{2\}\. This is ensured by also constraining the input and skip\-connection weights componentwise according to

𝑼ℓ≥𝟎,ℓ=1,…,L\.\\bm\{U\}\_\{\\ell\}\\geq\\bm\{0\},\\qquad\\ell=1,\\ldots,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𝒙\\bm\{x\}\. The biases𝒃ℓ\\bm\{b\}\_\{\\ell\}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\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), 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\.5Constitutive Artificial Neural Networks \(CANNs\)

Constitutive Artificial Neural Networks\([26](https://arxiv.org/html/2608.14063#bib.bib7);[27](https://arxiv.org/html/2608.14063#bib.bib14)\)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](https://arxiv.org/html/2608.14063#bib.bib19)\)\.

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\\bar\{I\}\_\{1\}=I\_\{1\}\-3,\\qquad\\bar\{I\}\_\{2\}=I\_\{2\}\-3\(36\)as inputs\. SinceI1≥3I\_\{1\}\\geq 3andI2≥3I\_\{2\}\\geq 3for 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\.\\bm\{q\}\(\\bar\{I\}\_\{1\},\\bar\{I\}\_\{2\}\)=\\begin\{bmatrix\}\\bar\{I\}\_\{1\},&\\bar\{I\}\_\{2\},&\\bar\{I\}\_\{1\}^\{2\},&\\bar\{I\}\_\{2\}^\{2\}\\end\{bmatrix\}^\{T\}\.\(37\)Each feature enters the strain energy density both directly and through an exponential activation\.111For 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 formexp⁡\(w1,kexp​I¯i2\)\\exp\(w\_\{1,k\}^\{\\mathrm\{exp\}\}\\bar\{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,kid​qk\+∑k=14w2,kexp​\[exp⁡\(w1,kexp​qk\)−1\],W\(I\_\{1\},I\_\{2\}\)=\\sum\_\{k=1\}^\{4\}w\_\{2,k\}^\{\\mathrm\{id\}\}\\,q\_\{k\}\+\\sum\_\{k=1\}^\{4\}w\_\{2,k\}^\{\\mathrm\{exp\}\}\\left\[\\exp\\left\(w\_\{1,k\}^\{\\mathrm\{exp\}\}\\,q\_\{k\}\\right\)\-1\\right\],\(38\)wherew2,kidw\_\{2,k\}^\{\\mathrm\{id\}\}are the output weights of the identity branches,w1,kexpw\_\{1,k\}^\{\\mathrm\{exp\}\}control the exponential activations, andw2,kexpw\_\{2,k\}^\{\\mathrm\{exp\}\}are the corresponding exponential output weights\. Hence, the CANN comprises a total of 12 trainable weights\. This construction ensures that

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\}^\{\\mathrm\{id\}\}\\geq 0,\\qquad w\_\{1,k\}^\{\\mathrm\{exp\}\}\\geq 0,\\qquad w\_\{2,k\}^\{\\mathrm\{exp\}\}\\geq 0,\\qquad k=1,\\ldots,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 featuresqkq\_\{k\}are convex and non\-decreasing functions of their respective invariant\. Forw1,kexp≥0w\_\{1,k\}^\{\\mathrm\{exp\}\}\\geq 0, the exponential mapping

qk↦exp⁡\(w1,kexp​qk\)−1q\_\{k\}\\mapsto\\exp\\left\(w\_\{1,k\}^\{\\mathrm\{exp\}\}\\,q\_\{k\}\\right\)\-1\(41\)is likewise convex and non\-decreasing\. A non\-negative weighted sum of these terms therefore remains convex and componentwise non\-decreasing inI1I\_\{1\}andI2I\_\{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𝑭\\bm\{F\}andcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\)according to[Eq\.27](https://arxiv.org/html/2608.14063#S3.E27)\([20](https://arxiv.org/html/2608.14063#bib.bib35)\)\. Because the CANN is convex and non\-decreasing in both invariants, its composition with these convex functions admits a jointly convex representation in𝑭\\bm\{F\}andcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\)\([22](https://arxiv.org/html/2608.14063#bib.bib34)\)\. 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\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), 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\.6Generalized\-Invariant\-based Constitutive Artificial Neural Networks \(GI\-CANNs\)

Generalized\-Invariant\-based Constitutive Artificial Neural Networks\([28](https://arxiv.org/html/2608.14063#bib.bib37)\)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 invariantsI1I\_\{1\}andI2I\_\{2\}, GI\-CANNs introduce the continuous family of generalized invariants\([3](https://arxiv.org/html/2608.14063#bib.bib26)\)

𝒥α=∑i=13λiα,\\mathcal\{J\}\_\{\\alpha\}=\\sum\_\{i=1\}^\{3\}\\lambda\_\{i\}^\{\\alpha\},\(42\)whereλi\\lambda\_\{i\}are the principal stretches andα∈ℝ\\alpha\\in\\mathbb\{R\}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\.\\mathcal\{J\}\_\{2\}=I\_\{1\},\\qquad\\mathcal\{J\}\_\{\-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\\left\(\\left\\\{\\mathcal\{J\}\_\{\\alpha\_\{k\}\}\\right\\\}\_\{\\alpha\_\{k\}\\in\\mathcal\{S\}\}\\right\),\(44\)where𝒮⊂ℝ\\mathcal\{S\}\\subset\\mathbb\{R\}denotes the set of invariant exponents\. Fixing𝒮=\{2,−2\}\\mathcal\{S\}=\\\{2,\-2\\\}recovers the classical invariant representation underlying the standard\-invariant\-based CANN\([27](https://arxiv.org/html/2608.14063#bib.bib14)\)\. In contrast, prescribing a discrete grid of exponents and combining the resulting generalized invariants through identity branches recovers the principal\-stretch\-based CANN\([31](https://arxiv.org/html/2608.14063#bib.bib38)\)\. 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\.\\mathcal\{S\}=\\left\\\{\\alpha\_\{\+\},\-\\alpha\_\{\-\}\\right\\\},\\qquad\\alpha\_\{\+\}\\geq 0,\\quad\\alpha\_\{\-\}\\geq 0\.\(45\)The corresponding shifted generalized invariants are

𝒥¯\+=∑i=13λiα\+−3,𝒥¯−=∑i=13λi−α−−3\.\\bar\{\\mathcal\{J\}\}\_\{\+\}=\\sum\_\{i=1\}^\{3\}\\lambda\_\{i\}^\{\\alpha\_\{\+\}\}\-3,\\qquad\\bar\{\\mathcal\{J\}\}\_\{\-\}=\\sum\_\{i=1\}^\{3\}\\lambda\_\{i\}^\{\-\\alpha\_\{\-\}\}\-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=\\sum\_\{k\\in\\\{\+,\-\\\}\}\\left\\\{w\_\{2,k\}^\{\\mathrm\{id\}\}\\,\\bar\{\\mathcal\{J\}\}\_\{k\}\+w\_\{2,k\}^\{\\mathrm\{exp\}\}\\left\[\\exp\\left\(w\_\{1,k\}^\{\\mathrm\{exp\}\}\\bar\{\\mathcal\{J\}\}\_\{k\}\\right\)\-1\\right\]\\right\\\}\.\(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\.\\alpha\_\{\+\}\\geq 0,\\qquad\\alpha\_\{\-\}\\geq 0,\\qquad w\_\{1,k\}^\{\\mathrm\{exp\}\}\\geq 0,\\qquad w\_\{2,k\}^\{\\mathrm\{id\}\}\\geq 0,\\qquad w\_\{2,k\}^\{\\mathrm\{exp\}\}\\geq 0\.\(48\)For an incompressible deformation, the product of the principal stretches satisfiesλ1​λ2​λ3=1\\lambda\_\{1\}\\lambda\_\{2\}\\lambda\_\{3\}=1\. Consequently,𝒥α≥3\\mathcal\{J\}\_\{\\alpha\}\\geq 3for any real exponentα\\alpha, 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\\geq 0,\\qquad 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\\alpha\_\{\+\}\\geq 1, the generalized invariant

𝒥α\+=∑i=13σi​\(𝑭\)α\+\\mathcal\{J\}\_\{\\alpha\_\{\+\}\}=\\sum\_\{i=1\}^\{3\}\\sigma\_\{i\}\(\\bm\{F\}\)^\{\\alpha\_\{\+\}\}\(50\)is a convex spectral function of𝑭\\bm\{F\}, whereσi​\(𝑭\)=λi\\sigma\_\{i\}\(\\bm\{F\}\)=\\lambda\_\{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⁡\(𝑭\)\)α−\\mathcal\{J\}\_\{\-\\alpha\_\{\-\}\}=\\sum\_\{i=1\}^\{3\}\\lambda\_\{i\}^\{\-\\alpha\_\{\-\}\}=\\sum\_\{i=1\}^\{3\}\\sigma\_\{i\}\\left\(\\operatorname\{cof\}\(\\bm\{F\}\)\\right\)^\{\\alpha\_\{\-\}\}\(51\)is a convex spectral function ofcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\)forα−≥1\\alpha\_\{\-\}\\geq 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𝑭\\bm\{F\}andcof⁡\(𝑭\)\\operatorname\{cof\}\(\\bm\{F\}\), and is therefore polyconvex\.

In the present implementation, the exponent magnitudes are constrained only byα\+,α−≥0\\alpha\_\{\+\},\\alpha\_\{\-\}\\geq 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\\alpha\_\{\+\}\\geq 1,\\qquad\\alpha\_\{\-\}\\geq 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\.

HyperparametersAs shown in[Table2](https://arxiv.org/html/2608.14063#S3.T2), the hyperparameters of GI\-CANNs are equivalent to those of the CANNs discussed earlier\.

## 4Results

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[AppendixA](https://arxiv.org/html/2608.14063#A1)\.

### 4\.1Predictive 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](https://arxiv.org/html/2608.14063#S4.F6)and[2](https://arxiv.org/html/2608.14063#S1.F2)\. 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[Table4](https://arxiv.org/html/2608.14063#S4.T4)and[Figs\.4](https://arxiv.org/html/2608.14063#S4.F4)and[5](https://arxiv.org/html/2608.14063#S4.F5)show that all methods capture the dominant constitutive response\.

EUCLID provides the highest accuracy for uniaxial tension and pure shear, withR2=0\.9985R^\{2\}=0\.9985and0\.99780\.9978, respectively\. GI\-CANN performs particularly well for equibiaxial tension and reachesR2=0\.9978R^\{2\}=0\.9978, whereas EUCLID reachesR2=0\.9930R^\{2\}=0\.9930for 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](https://arxiv.org/html/2608.14063#S4.F6)\. 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\.2Best\-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\.[Table5](https://arxiv.org/html/2608.14063#S4.T5)summarizes the selected configurations, while[Fig\.2](https://arxiv.org/html/2608.14063#S1.F2)compares their predictions with the experimental data\. AMF provides the highest validation accuracy among the selected models, withR2=0\.9996R^\{2\}=0\.9996\. GI\-CANN provides the highest training accuracy withR2=0\.9996R^\{2\}=0\.9996and maintains a validation value ofR2=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 validationR2R^\{2\}of0\.99780\.9978\. GI\-CANN uses eight nonzero parameters, whereas CANN, AMF, and PANN use 12, 49, and 371 nonzero parameters, respectively, see[Table5](https://arxiv.org/html/2608.14063#S4.T5)\. 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\. ForNa=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](https://arxiv.org/html/2608.14063#bib.bib31)\)\. 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\.3Variability 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 ofNaN\_\{a\}andss\. PANN, CANN, and GI\-CANN vary across the 100 random initializations, see[Tables4](https://arxiv.org/html/2608.14063#S4.T4),[4](https://arxiv.org/html/2608.14063#S4.F4)and[5](https://arxiv.org/html/2608.14063#S4.F5)\. For AMF, we note that we includedNa=1N\_\{a\}=1for completeness, even though this reduces AMF to the classical MF with only one feature\.

PANN shows the largest variability, particularly for equibiaxial tension\. Its meanR2R^\{2\}reaches only 0\.9039 for this deformation mode, with a standard deviation of 0\.2450\. In contrast, GI\-CANN reachesR2=0\.9978±0\.0035R^\{2\}=0\.9978\\pm 0\.0035for 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](https://arxiv.org/html/2608.14063#S4.F6)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\.4Computational cost of model identification

The six methods show large differences in computational cost as shown in[Tables3](https://arxiv.org/html/2608.14063#S4.T3)and[3](https://arxiv.org/html/2608.14063#S4.F3)\. EUCLID and MF require less than0\.40\.4ms for their core identification procedures, with mean core times of0\.3370\.337ms and0\.3050\.305ms, respectively\. AMF requires1\.0921\.092ms 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 from0\.4910\.491ms for EUCLID and18\.20118\.201ms for MF to approximately0\.40\.4–0\.80\.8s for the three neural\-network methods, see[Table3](https://arxiv.org/html/2608.14063#S4.T3)\.

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](https://arxiv.org/html/2608.14063#bib.bib32)\)\. 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\.5Computational 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 ofWW,𝑷\\bm\{P\}, andℍ\\mathbb\{H\}for the six selected models, as shown in[Tables6](https://arxiv.org/html/2608.14063#S4.T6)and[8](https://arxiv.org/html/2608.14063#S4.F8)\.

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 requires0\.151​μ​s0\.151\\,\\mu\\mathrm\{s\}per deformation gradient forWW,0\.473​μ​s0\.473\\,\\mu\\mathrm\{s\}for𝑷\\bm\{P\}, and3\.263​μ​s3\.263\\,\\mu\\mathrm\{s\}forℍ\\mathbb\{H\}\. The Hessian requires the largest computation for every method, with AMF reaching7\.750​μ​s7\.750\\,\\mu\\mathrm\{s\}compared with3\.263​μ​s3\.263\\,\\mu\\mathrm\{s\}for MF, see[Table6](https://arxiv.org/html/2608.14063#S4.T6)\.

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](https://arxiv.org/html/2608.14063#S4.F8)\.

### 4\.6Strain energy landscapes

We finally compare the constitutive functions beyond the deformation paths used for calibration\.[Fig\.7](https://arxiv.org/html/2608.14063#S4.F7)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\)\(\\lambda\_\{1\},\\lambda\_\{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](https://arxiv.org/html/2608.14063#S1.F2)\.

The direct contour overlay in[Fig\.1](https://arxiv.org/html/2608.14063#S1.F1)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](https://arxiv.org/html/2608.14063#S1.F1)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](https://arxiv.org/html/2608.14063#S4.F7)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](https://arxiv.org/html/2608.14063#S1.F1)and[7](https://arxiv.org/html/2608.14063#S4.F7)\. 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±\\pmstandard 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\.MethodCore time\(ms\)Setup time\(ms\)Total time\(ms\)EUCLID0\.337±0\.1300\.337\\pm 0\.1300\.137±0\.003\\mathbf\{0\.137\\pm 0\.003\}0\.491±0\.133\\mathbf\{0\.491\\pm 0\.133\}MF0\.305±0\.006\\mathbf\{0\.305\\pm 0\.006\}17\.896±3\.49417\.896\\pm 3\.49418\.201±3\.49318\.201\\pm 3\.493AMF1\.092±0\.1401\.092\\pm 0\.14021\.426±2\.85421\.426\\pm 2\.85422\.518±2\.84022\.518\\pm 2\.840PANN588\.495±19\.379588\.495\\pm 19\.379240\.347±20\.282240\.347\\pm 20\.282828\.842±28\.950828\.842\\pm 28\.950CANN324\.539±4\.291324\.539\\pm 4\.291108\.781±4\.921108\.781\\pm 4\.921433\.319±6\.527433\.319\\pm 6\.527GI\-CANN681\.277±14\.110681\.277\\pm 14\.110117\.587±12\.135117\.587\\pm 12\.135798\.863±18\.259798\.863\\pm 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±\\pmstandard deviation\. The sample sizennhas a method\-specific meaning: Adaptive Material Fingerprinting \(AMF\) comprises 20 combinations of the adaptive parametersNaN\_\{a\}andss, 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\.MethodnnMSERMSENRMSEMAER2R^\{2\}Uniaxial tension \(UT\)EUCLID10\.0025±\\pm00\.0501±\\pm00\.0112±\\pm00\.0389±\\pm00\.9985±\\pm0MF10\.0137±\\pm00\.1172±\\pm00\.0261±\\pm00\.0962±\\pm00\.9920±\\pm0AMF200\.0258±\\pm0\.02070\.1443±\\pm0\.07220\.0321±\\pm0\.01610\.1033±\\pm0\.05880\.9850±\\pm0\.0120PANN1000\.0103±\\pm0\.01560\.0874±\\pm0\.05200\.0195±\\pm0\.01160\.0698±\\pm0\.03920\.9940±\\pm0\.0091CANN1000\.0093±\\pm0\.00950\.0887±\\pm0\.03850\.0198±\\pm0\.00860\.0725±\\pm0\.02950\.9946±\\pm0\.0055GI\-CANN1000\.0074±\\pm0\.01120\.0619±\\pm0\.06010\.0138±\\pm0\.01340\.0507±\\pm0\.05020\.9957±\\pm0\.0065Equibiaxial tension \(ET\)EUCLID10\.0041±\\pm00\.0640±\\pm00\.0261±\\pm00\.0538±\\pm00\.9930±\\pm0MF10\.0407±\\pm00\.2016±\\pm00\.0823±\\pm00\.1725±\\pm00\.9310±\\pm0AMF200\.0355±\\pm0\.04070\.1493±\\pm0\.11800\.0610±\\pm0\.04820\.0910±\\pm0\.06130\.9397±\\pm0\.0690PANN1000\.0566±\\pm0\.14440\.1682±\\pm0\.16920\.0687±\\pm0\.06900\.1219±\\pm0\.10080\.9039±\\pm0\.2450CANN1000\.0150±\\pm0\.00990\.1176±\\pm0\.03420\.0480±\\pm0\.01400\.0974±\\pm0\.02900\.9746±\\pm0\.0168GI\-CANN1000\.0013±\\pm0\.00200\.0278±\\pm0\.02350\.0114±\\pm0\.00960\.0221±\\pm0\.01890\.9978±\\pm0\.0035Pure shear \(PS\)EUCLID10\.0008±\\pm00\.0277±\\pm00\.0152±\\pm00\.0221±\\pm00\.9978±\\pm0MF10\.0035±\\pm00\.0596±\\pm00\.0327±\\pm00\.0503±\\pm00\.9899±\\pm0AMF200\.0042±\\pm0\.00520\.0504±\\pm0\.04200\.0277±\\pm0\.02310\.0378±\\pm0\.03080\.9880±\\pm0\.0149PANN1000\.0037±\\pm0\.00550\.0516±\\pm0\.03190\.0283±\\pm0\.01750\.0423±\\pm0\.02620\.9895±\\pm0\.0158CANN1000\.0085±\\pm0\.01260\.0805±\\pm0\.04500\.0442±\\pm0\.02470\.0671±\\pm0\.04090\.9758±\\pm0\.0359GI\-CANN1000\.0023±\\pm0\.00330\.0374±\\pm0\.03080\.0205±\\pm0\.01690\.0293±\\pm0\.02720\.9933±\\pm0\.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 determinationR2R^\{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\.MethodAdaptiveparametersSeedTrainingMSEValidationMSETrainingR2R^\{2\}ValidationR2R^\{2\}Core time\(ms\)\# Parameters‖𝜽‖0\\\|\\bm\{\\theta\}\\\|\_\{0\}EUCLID––3\.304×10−33\.304\\times 10^\{\-3\}7\.663×10−47\.663\\times 10^\{\-4\}0\.99730\.99780\.3724MF––2\.719×10−22\.719\\times 10^\{\-2\}3\.547×10−33\.547\\times 10^\{\-3\}0\.97810\.98990\.3022AMFNa=20N\_\{a\}=20,s=0\.25s=0\.25–1\.119×10−31\.119\\times 10^\{\-3\}1\.524×𝟏𝟎−𝟒\\mathbf\{1\.524\\times 10^\{\-4\}\}0\.99910\.99966\.00049PANN–63\.329×10−33\.329\\times 10^\{\-3\}8\.186×10−48\.186\\times 10^\{\-4\}0\.99730\.9977581\.311371CANN–933\.535×10−33\.535\\times 10^\{\-3\}8\.925×10−48\.925\\times 10^\{\-4\}0\.99720\.9975320\.94712GI\-CANN–394\.436×𝟏𝟎−𝟒\\mathbf\{4\.436\\times 10^\{\-4\}\}2\.412×10−42\.412\\times 10^\{\-4\}0\.99960\.9993664\.9888

Table 6:Computational benchmark for evaluating the strain energy density functionWW, its gradient𝑷\\bm\{P\}, and its Hessianℍ\\mathbb\{H\}for the six methods\. Times are reported in microseconds per sample as mean±\\pmstandard 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\.MethodEnergyWW\(μ​s\\mu\\mathrm\{s\}/sample\)Gradient𝑷\\bm\{P\}\(μ​s\\mu\\mathrm\{s\}/sample\)Hessianℍ\\mathbb\{H\}\(μ​s\\mu\\mathrm\{s\}/sample\)EUCLID0\.233±0\.0440\.233\\pm 0\.0440\.556±0\.0390\.556\\pm 0\.0395\.734±0\.1365\.734\\pm 0\.136MF0\.151±0\.033\\mathbf\{0\.151\\pm 0\.033\}0\.473±0\.030\\mathbf\{0\.473\\pm 0\.030\}3\.263±0\.108\\mathbf\{3\.263\\pm 0\.108\}AMF0\.430±0\.0510\.430\\pm 0\.0511\.262±0\.0601\.262\\pm 0\.0607\.750±0\.3917\.750\\pm 0\.391PANN0\.428±0\.0590\.428\\pm 0\.0590\.827±0\.0640\.827\\pm 0\.0645\.496±0\.2185\.496\\pm 0\.218CANN0\.175±0\.0370\.175\\pm 0\.0370\.548±0\.0460\.548\\pm 0\.0464\.578±0\.0954\.578\\pm 0\.095GI\-CANN0\.214±0\.0360\.214\\pm 0\.0360\.556±0\.0470\.556\\pm 0\.0474\.332±0\.2724\.332\\pm 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 parametersNaN\_\{a\}andss; 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 determinationR2R^\{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 parametersNaN\_\{a\}andss; and PANN, CANN, and GI\-CANN each comprise 100 model initializations using different random seeds\. Values closer to one indicate better predictive performance\.![Refer to caption](https://arxiv.org/html/2608.14063v1/prediction_method_grid_6x4.png)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 componentP11P\_\{11\}as a function of the stretch ratioλ\\lambdafor 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 densityW⁡\(𝑭\)W\(\\bm\{F\}\)predicted by the selected best\-performing models: EUCLID, Material Fingerprinting \(MF\), Adaptive Material Fingerprinting \(AMF\), PANN, CANN, and GI\-CANN\. Each panel showsWWas a function of the two independent principal stretchesλ1\\lambda\_\{1\}andλ2\\lambda\_\{2\}, with the third stretch determined by the incompressibility constraintλ3=\(λ1​λ2\)−1\\lambda\_\{3\}=\(\\lambda\_\{1\}\\lambda\_\{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\)\(\\lambda\_\{1\},\\lambda\_\{2\}\)=\(1,1\)denotes the undeformed reference configuration\.Figure 8:Runtime comparison for evaluating the strain energy density functionWW, its gradient𝑷\\bm\{P\}, and its Hessianℍ\\mathbb\{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\.

## 5Conclusion

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](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

- R\. Abdusalamov, M\. Hillgärtner, and M\. ItskovAutomatic generation of interpretable hyperelastic material models by symbolic regression\.International Journal for Numerical Methods in Engineering,pp\. nme\.7203\(en\)\.External Links:ISSN 0029\-5981, 1097\-0207,[Link](https://onlinelibrary.wiley.com/doi/10.1002/nme.7203),[Document](https://dx.doi.org/10.1002/nme.7203)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Amoset al\.\(2017\)B\. Amos, L\. Xu, and J\. Z\. KolterInput Convex Neural Networks\.InProceedings of the 34th International Conference on Machine Learning,Sydney, Australia,pp\. 10\(en\)\.Cited by:[§3\.4](https://arxiv.org/html/2608.14063#S3.SS4.p1.1)\.
- Anssari\-Benamet al\.\(2024\)A\. Anssari\-Benam, A\. Goriely, and G\. SaccomandiGeneralised invariants and pseudo\-universal relationships for hyperelastic materials: A new approach to constitutive modelling\.Journal of the Mechanics and Physics of Solids193,pp\. 105883\(en\)\.External Links:ISSN 00225096,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0022509624003491),[Document](https://dx.doi.org/10.1016/j.jmps.2024.105883)Cited by:[§3\.6](https://arxiv.org/html/2608.14063#S3.SS6.p1.1)\.
- Bleyer \(2025\)J\. BleyerLearning elastoplasticity with implicit layers\.\(en\)\.External Links:[Document](https://dx.doi.org/http%3A//dx.doi.org/10.2139/ssrn.5210734)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Bruntonet al\.\(2016\)S\. L\. Brunton, J\. L\. Proctor, and J\. N\. KutzDiscovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences113\(15\),pp\. 3932–3937\(en\)\.External Links:ISSN 0027\-8424, 1091\-6490,[Link](http://www.pnas.org/lookup/doi/10.1073/pnas.1517384113),[Document](https://dx.doi.org/10.1073/pnas.1517384113)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Dammaßet al\.\(2025\)F\. Dammaß, K\. A\. Kalina, and M\. KästnerWhen invariants matter: The role of I1 and I2 in neural network models of incompressible hyperelasticity\.Mechanics of Materials,pp\. 105443\.External Links:ISSN 0167\-6636,[Document](https://dx.doi.org/10.1016/j.mechmat.2025.105443)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Efronet al\.\(2004\)B\. Efron, T\. Hastie, I\. Johnstone, and R\. TibshiraniLEAST ANGLE REGRESSION\.The Annals of Statistics32\(2\),pp\. 407–499\(en\)\.Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Flaschelet al\.\(2026a\)M\. Flaschel, T\. Hastie, and E\. KuhlNon\-smooth optimization meets automated material model discovery\.Journal of Computational Physics563,pp\. 115043\(en\)\.External Links:ISSN 00219991,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0021999126003967),[Document](https://dx.doi.org/10.1016/j.jcp.2026.115043)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Flaschelet al\.\(2026b\)M\. Flaschel, H\. Holthusen, D\. Martonová, and E\. KuhlAdaptive Material Fingerprinting for the fast discovery of polyconvex feature combinations in isotropic and anisotropic hyperelasticity\.International Journal of Engineering Science227,pp\. 104594\(en\)\.External Links:ISSN 00207225,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0020722526001321),[Document](https://dx.doi.org/10.1016/j.ijengsci.2026.104594)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.3](https://arxiv.org/html/2608.14063#S3.SS3.p1.1),[§3\.3](https://arxiv.org/html/2608.14063#S3.SS3.p2.2),[§3\.3](https://arxiv.org/html/2608.14063#S3.SS3.p5.1),[§3\.3](https://arxiv.org/html/2608.14063#S3.SS3.p6.pic1.2.1.1),[§4\.2](https://arxiv.org/html/2608.14063#S4.SS2.p3.1)\.
- Flaschelet al\.\(2021\)M\. Flaschel, S\. Kumar, and L\. De LorenzisUnsupervised discovery of interpretable hyperelastic constitutive laws\.Computer Methods in Applied Mechanics and Engineering381,pp\. 113852\(en\)\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2021.113852)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p1.1),[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p5.1),[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p6.pic1.2.1.1)\.
- Flaschelet al\.\(2026c\)M\. Flaschel, D\. Martonová, C\. Veil, and E\. KuhlMaterial Fingerprinting: A shortcut to material model discovery without solving optimization problems\.Computer Methods in Applied Mechanics and Engineering450,pp\. 118573\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S004578252500845X),[Document](https://dx.doi.org/10.1016/j.cma.2025.118573)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p1.1),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p3.2),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p5.1)\.
- Flaschelet al\.\(2026d\)M\. Flaschel, M\. A\. Moreno\-Mateos, S\. Wiesheier, P\. Steinmann, and E\. KuhlUnsupervised material fingerprinting: Ultra\-fast hyperelastic model discovery from full\-field experimental measurements\.Computer Methods in Applied Mechanics and Engineering461,pp\. 119256\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0045782526005293),[Document](https://dx.doi.org/10.1016/j.cma.2026.119256)Cited by:[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p2.1),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p5.1),[§4\.4](https://arxiv.org/html/2608.14063#S4.SS4.p4.1)\.
- Flaschelet al\.\(2023\)M\. Flaschel, H\. Yu, N\. Reiter, J\. Hinrichsen, S\. Budday, P\. Steinmann, S\. Kumar, and L\. De LorenzisAutomated discovery of interpretable hyperelastic material models for human brain tissue with EUCLID\.Journal of the Mechanics and Physics of Solids180,pp\. 105404\(en\)\.External Links:ISSN 00225096,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0022509623002089),[Document](https://dx.doi.org/10.1016/j.jmps.2023.105404)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p1.1),[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p2.2)\.
- Frank and Friedman \(1993\)I\. E\. Frank and J\. H\. FriedmanA Statistical View of Some Chemometrics Regression Tools\.Technometrics35\(2\),pp\. 109–135\(en\)\.External Links:ISSN 0040\-1706, 1537\-2723,[Link](http://www.tandfonline.com/doi/abs/10.1080/00401706.1993.10485033),[Document](https://dx.doi.org/10.1080/00401706.1993.10485033)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Frankelet al\.\(2020\)A\. L\. Frankel, R\. E\. Jones, and L\. P\. SwilerTENSOR BASIS GAUSSIAN PROCESS MODELS OF HYPERELASTIC MATERIALS\.Journal of Machine Learning for Modeling and Computing1\(1\),pp\. 1–17\(en\)\.External Links:ISSN 2689\-3967,[Link](http://www.dl.begellhouse.com/journals/558048804a15188a,583c4e56625ba94e,651a2e6b0260f708.html),[Document](https://dx.doi.org/10.1615/JMachLearnModelComput.2020033325)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Fu \(1998\)W\. FuPenalized Regressions: The Bridge versus the Lasso\.\(en\)\.Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Fuhget al\.\(2024\)J\. N\. Fuhg, G\. Anantha Padmanabha, N\. Bouklas, B\. Bahmani, W\. Sun, N\. N\. Vlassis, M\. Flaschel, P\. Carrara, and L\. De LorenzisA Review on Data\-Driven Constitutive Laws for Solids\.Archives of Computational Methods in Engineering\(en\)\.External Links:ISSN 1134\-3060, 1886\-1784,[Link](https://link.springer.com/10.1007/s11831-024-10196-2),[Document](https://dx.doi.org/10.1007/s11831-024-10196-2)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Fuhget al\.\(2022\)J\. N\. Fuhg, M\. Marino, and N\. BouklasLocal approximate Gaussian process regression for data\-driven constitutive models: development and comparison with neural networks\.Computer Methods in Applied Mechanics and Engineering388,pp\. 114217\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S004578252100548X),[Document](https://dx.doi.org/10.1016/j.cma.2021.114217)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Hastieet al\.\(2009\)T\. Hastie, R\. Tibshirani, and J\. FriedmanThe 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](http://link.springer.com/10.1007/978-0-387-84858-7),[Document](https://dx.doi.org/10.1007/978-0-387-84858-7)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Holthusen and Kuhl \(2026\)H\. Holthusen and E\. KuhlA complement to neural networks for anisotropic inelasticity at finite strains\.Computer Methods in Applied Mechanics and Engineering450,pp\. 118612\.External Links:ISSN 0045\-7825,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.cma.2025.118612),[Link](https://www.sciencedirect.com/science/article/pii/S0045782525008849)Cited by:[§3\.5](https://arxiv.org/html/2608.14063#S3.SS5.p4.1)\.
- Holthusenet al\.\(2024\)H\. Holthusen, L\. Lamm, T\. Brepols, S\. Reese, and E\. KuhlTheory and implementation of inelastic Constitutive Artificial Neural Networks\.Computer Methods in Applied Mechanics and Engineering428,pp\. 117063\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0045782524003190),[Document](https://dx.doi.org/10.1016/j.cma.2024.117063)Cited by:[§3\.5](https://arxiv.org/html/2608.14063#S3.SS5.p1.1)\.
- Holthusenet al\.\(2026\)H\. Holthusen, K\. Linka, E\. Kuhl, and T\. BrepolsA generalized dual potential for inelastic constitutive artificial neural networks: a jax implementation at finite strains\.Journal of the Mechanics and Physics of Solids206,pp\. 106337\.External Links:ISSN 0022\-5096,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.jmps.2025.106337),[Link](https://www.sciencedirect.com/science/article/pii/S0022509625003084)Cited by:[§3\.5](https://arxiv.org/html/2608.14063#S3.SS5.p4.1)\.
- Kleinet al\.\(2022\)D\. K\. Klein, M\. Fernández, R\. J\. Martin, P\. Neff, and O\. WeegerPolyconvex anisotropic hyperelasticity with neural networks\.Journal of the Mechanics and Physics of Solids159,pp\. 104703\(en\)\.External Links:ISSN 00225096,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0022509621003215),[Document](https://dx.doi.org/10.1016/j.jmps.2021.104703)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.4](https://arxiv.org/html/2608.14063#S3.SS4.p1.1)\.
- Kopp \(2008\)J\. KoppEfficient numerical diagonalization of hermitian 3 × 3 matrices\.International Journal of Modern Physics C19\(03\),pp\. 523–548\.External Links:[Document](https://dx.doi.org/10.1142/S0129183108012303),[Link](https://doi.org/10.1142/S0129183108012303)Cited by:[Appendix A](https://arxiv.org/html/2608.14063#A1.p6.3)\.
- Lindenet al\.\(2023\)L\. Linden, D\. K\. Klein, K\. A\. Kalina, J\. Brummund, O\. Weeger, and M\. KästnerNeural networks meet hyperelasticity: A guide to enforcing physics\.Journal of the Mechanics and Physics of Solids179,pp\. 105363\(en\)\.External Links:ISSN 00225096,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0022509623001679),[Document](https://dx.doi.org/10.1016/j.jmps.2023.105363)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.4](https://arxiv.org/html/2608.14063#S3.SS4.p1.1),[§3\.4](https://arxiv.org/html/2608.14063#S3.SS4.p1.2)\.
- Linkaet al\.\(2021\)K\. Linka, M\. Hillgärtner, K\. P\. Abdolazizi, R\. C\. Aydin, M\. Itskov, and C\. J\. CyronConstitutive artificial neural networks: A fast and general approach to predictive data\-driven constitutive modeling by deep learning\.Journal of Computational Physics429,pp\. 110010\(en\)\.External Links:ISSN 00219991,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0021999120307841),[Document](https://dx.doi.org/10.1016/j.jcp.2020.110010)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.5](https://arxiv.org/html/2608.14063#S3.SS5.p1.1)\.
- Linka and Kuhl \(2023\)K\. Linka and E\. KuhlA new family of Constitutive Artificial Neural Networks towards automated model discovery\.Computer Methods in Applied Mechanics and Engineering403,pp\. 115731\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0045782522006867),[Document](https://dx.doi.org/10.1016/j.cma.2022.115731)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.5](https://arxiv.org/html/2608.14063#S3.SS5.p1.1),[§3\.6](https://arxiv.org/html/2608.14063#S3.SS6.p2.2)\.
- Martonováet al\.\(2026a\)D\. Martonová, A\. Goriely, and E\. KuhlGeneralized invariants meet constitutive neural networks: A novel framework for hyperelastic materials\.Journal of the Mechanics and Physics of Solids206,pp\. 106352\.External Links:ISSN 0022\-5096,[Document](https://dx.doi.org/10.1016/j.jmps.2025.106352)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.6](https://arxiv.org/html/2608.14063#S3.SS6.p1.1)\.
- Martonováet al\.\(2026b\)D\. Martonová, E\. Kuhl, and M\. FlaschelMaterial Fingerprinting for rapid discovery of hyperelastic models: First experimental validation\.Journal of the Mechanics and Physics of Solids208,pp\. 106463\.External Links:ISSN 0022\-5096,[Document](https://dx.doi.org/10.1016/j.jmps.2025.106463)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p1.1),[§3\.2](https://arxiv.org/html/2608.14063#S3.SS2.p2.1)\.
- McCullochet al\.\(2024\)J\. A\. McCulloch, S\. R\. St\. Pierre, K\. Linka, and E\. KuhlOn sparse regression,Lp\{\}\_\{\\textrm\{\{p\}\}\}‐regularization, and automated model discovery\.International Journal for Numerical Methods in Engineering125\(14\),pp\. e7481\(en\)\.External Links:ISSN 0029\-5981, 1097\-0207,[Link](https://onlinelibrary.wiley.com/doi/10.1002/nme.7481),[Document](https://dx.doi.org/10.1002/nme.7481)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p6.pic1.2.1.1)\.
- Pierreet al\.\(2023\)S\. R\. St\. Pierre, K\. Linka, and E\. KuhlPrincipal\-stretch\-based constitutive neural networks autonomously discover a subclass of Ogden models for human brain tissue\.Brain Multiphysics4,pp\. 100066\.External Links:ISSN 2666\-5220,[Document](https://dx.doi.org/10.1016/j.brain.2023.100066)Cited by:[§3\.6](https://arxiv.org/html/2608.14063#S3.SS6.p2.2)\.
- Steinmannet al\.\(2012\)P\. Steinmann, M\. Hossain, and G\. PossartHyperelastic models for rubber\-like materials: consistent tangent operators and suitability for Treloar’s data\.Archive of Applied Mechanics82\(9\),pp\. 1183–1217\(en\)\.External Links:ISSN 0939\-1533, 1432\-0681,[Link](http://link.springer.com/10.1007/s00419-012-0610-z),[Document](https://dx.doi.org/10.1007/s00419-012-0610-z)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p1.1)\.
- Taçet al\.\(2024\)V\. Taç, K\. Linka, F\. Sahli\-Costabal, E\. Kuhl, and A\. B\. TepoleBenchmarking physics\-informed frameworks for data\-driven hyperelasticity\.Computational Mechanics73\(1\),pp\. 49–65\(en\)\.External Links:ISSN 0178\-7675, 1432\-0924,[Link](https://link.springer.com/10.1007/s00466-023-02355-2),[Document](https://dx.doi.org/10.1007/s00466-023-02355-2)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Tacet al\.\(2022\)V\. Tac, F\. Sahli Costabal, and A\. B\. TepoleData\-driven tissue mechanics with polyconvex neural ordinary differential equations\.Computer Methods in Applied Mechanics and Engineering398,pp\. 115248\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0045782522003838),[Document](https://dx.doi.org/10.1016/j.cma.2022.115248)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Tackeet al\.\(2025\)M\. Tacke, M\. Busch, K\. Bali, K\. Abdolazizi, K\. Linka, C\. Cyron, and R\. AydinConstitutive scientific generative agent \(CSGA\): Leveraging large language models for automated constitutive model discovery\.Machine Learning for Computational Science and Engineering1\(1\),pp\. 23\(en\)\.External Links:ISSN 3005\-1428, 3005\-1436,[Link](https://link.springer.com/10.1007/s44379-025-00022-2),[Document](https://dx.doi.org/10.1007/s44379-025-00022-2)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Tibshirani \(1996\)R\. TibshiraniRegression Shrinkage and Selection via the Lasso\.Journal of the Royal Statistical Society: Series B \(Methodological\)58\(1\),pp\. 267–288\(en\)\.External Links:ISSN 00359246,[Link](http://doi.wiley.com/10.1111/j.2517-6161.1996.tb02080.x),[Document](https://dx.doi.org/10.1111/j.2517-6161.1996.tb02080.x)Cited by:[§3\.1](https://arxiv.org/html/2608.14063#S3.SS1.p4.2)\.
- Treloar \(1944\)L\. R\. G\. TreloarStress\-strain data for vulcanised rubber under various types of deformation\.Transactions of the Faraday Society40,pp\. 59\(en\)\.External Links:ISSN 0014\-7672,[Link](http://xlink.rsc.org/?DOI=tf9444000059),[Document](https://dx.doi.org/10.1039/tf9444000059)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p1.1),[§2\.1](https://arxiv.org/html/2608.14063#S2.SS1.p1.1)\.
- Wanget al\.\(2021\)Z\. Wang, J\.B\. Estrada, E\.M\. Arruda, and K\. GarikipatiInference 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 Solids153,pp\. 104474\(en\)\.External Links:ISSN 00225096,[Link](https://linkinghub.elsevier.com/retrieve/pii/S0022509621001459),[Document](https://dx.doi.org/10.1016/j.jmps.2021.104474)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.
- Wiesheieret al\.\(2024\)S\. Wiesheier, M\. A\. Moreno\-Mateos, and P\. SteinmannVersatile data\-adaptive hyperelastic energy functions for soft materials\.Computer Methods in Applied Mechanics and Engineering430,pp\. 117208\(en\)\.External Links:ISSN 00457825,[Link](https://linkinghub.elsevier.com/retrieve/pii/S004578252400464X),[Document](https://dx.doi.org/10.1016/j.cma.2024.117208)Cited by:[§1](https://arxiv.org/html/2608.14063#S1.p3.1)\.

## Appendix ATechnical details

For all results, we use the following command

pythoncompare\_all\.py\\

\-\-n\-timing100\-\-n\-seeds100\-\-euclid\-lambdas1e\-5\\

\-\-study4\-num\-samples100\-\-study4\-repetitions100

The option\-\-n\-timing 100performs 100 independent runtime measurements for each of the six constitutive model identification methods\. The option\-\-n\-seeds 100evaluates the effect of neural\-network initialization by training PANN, CANN, and GI\-CANN with 100 different random seeds\. The option\-\-euclid\-lambdas 1e\-5restricts the EUCLID hyperparameter study to the regularization parameterα=10−5\\alpha=10^\{\-5\}\.

The option\-\-study4\-num\-samples 100generates one common set of 100 admissible incompressible deformation gradients and evaluates all six models for the strain energy densityWW, the first Piola stress𝑷=∂W∂𝑭\\bm\{P\}=\\frac\{\\partial W\}\{\\partial\\bm\{F\}\}, and the Hessianℍ=∂2W∂𝑭​∂𝑭\\mathbb\{H\}=\\frac\{\\partial^\{2\}W\}\{\\partial\\bm\{F\}\\partial\\bm\{F\}\}\. The option\-\-study4\-repetitions 100repeats 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[Table7](https://arxiv.org/html/2608.14063#A1.T7)\.

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[Section4\.5](https://arxiv.org/html/2608.14063#S4.SS5), we compute the principal stretches as the singular values of the deformation gradient𝑭\\bm\{F\}\. Instead of explicitly performing a singular\-value decomposition, we exploit that the squared singular values of𝑭\\bm\{F\}are the eigenvalues of the right Cauchy–Green tensor

𝑪=𝑭𝖳​𝑭\.\\bm\{C\}=\\bm\{F\}^\{\\mathsf\{T\}\}\\bm\{F\}\.\(53\)Denoting the eigenvalues of𝑪\\bm\{C\}byμi\\mu\_\{i\}, the principal stretches follow as

λi=μi,i=1,2,3\.\\lambda\_\{i\}=\\sqrt\{\\mu\_\{i\}\},\\qquad i=1,2,3\.\(54\)Since𝑪\\bm\{C\}is symmetric, its eigenvalues are evaluated analytically using the closed\-form trigonometric solution for a symmetric3×33\\times 3matrix\([24](https://arxiv.org/html/2608.14063#bib.bib33)\)\. We first define

q=13​tr⁡\(𝑪\),𝑪^=𝑪−q​𝑰,q=\\frac\{1\}\{3\}\\operatorname\{tr\}\(\\bm\{C\}\),\\qquad\\widehat\{\\bm\{C\}\}=\\bm\{C\}\-q\\bm\{I\},\(55\)and

p=16​𝑪^:𝑪^\.p=\\sqrt\{\\frac\{1\}\{6\}\\widehat\{\\bm\{C\}\}:\\widehat\{\\bm\{C\}\}\}\.\(56\)Forp\>0p\>0, we introduce

𝑨=𝑪^p,r=12​det\(𝑨\),ϕ=13​arccos⁡\(r\),\\bm\{A\}=\\frac\{\\widehat\{\\bm\{C\}\}\}\{p\},\\qquad r=\\frac\{1\}\{2\}\\det\(\\bm\{A\}\),\\qquad\\phi=\\frac\{1\}\{3\}\\arccos\(r\),\(57\)whererris numerically restricted to the interval\[−1,1\]\[\-1,1\]\. The largest and smallest eigenvalues of𝑪\\bm\{C\}are then obtained as

μmax\\displaystyle\\mu\_\{\\max\}=q\+2​p​cos⁡\(ϕ\),\\displaystyle=q\+2p\\cos\(\\phi\),\(58\)μmin\\displaystyle\\mu\_\{\\min\}=q\+2​p​cos⁡\(ϕ\+2​π3\),\\displaystyle=q\+2p\\cos\\left\(\\phi\+\\frac\{2\\pi\}\{3\}\\right\),\(59\)while the remaining eigenvalue follows from the trace constraint as

μmid=3​q−μmax−μmin\.\\mu\_\{\\mathrm\{mid\}\}=3q\-\\mu\_\{\\max\}\-\\mu\_\{\\min\}\.\(60\)The principal stretches are subsequently computed as

λ1=μmax,λ2=μmid,λ3=μmin,\\lambda\_\{1\}=\\sqrt\{\\mu\_\{\\max\}\},\\qquad\\lambda\_\{2\}=\\sqrt\{\\mu\_\{\\mathrm\{mid\}\}\},\\qquad\\lambda\_\{3\}=\\sqrt\{\\mu\_\{\\min\}\},\(61\)such thatλ1≥λ2≥λ3\\lambda\_\{1\}\\geq\\lambda\_\{2\}\\geq\\lambda\_\{3\}\. For nearly isotropic states, for whichppfalls below a machine\-precision\-dependent threshold, we instead use

λ1=λ2=λ3=q,\\lambda\_\{1\}=\\lambda\_\{2\}=\\lambda\_\{3\}=\\sqrt\{q\},\(62\)which avoids the numerical degeneracy of the trigonometric representation and ensures a well\-defined evaluation of the undeformed reference configuration𝑭=𝑰\\bm\{F\}=\\bm\{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 parametersNaN\_\{a\}andss\.HyperparameterValueEUCLIDOptimization methodCoordinate Descentα\\alpha1×10−51\\times 10^\{\-5\}θth\\theta\_\{\\mathrm\{th\}\}0Material Fingerprinting \(MF\)Material databaseHyperelastic incompressible \(HEI\)Adaptive Material Fingerprinting \(AMF\)Material databaseHyperelastic incompressible isotropic adaptive \(HEIIA\)NaN\_\{a\}\[1, 3, 5, 10, 20\]ss\[0\.25, 0\.5, 0\.75, 1\.0\]Physics\-Augmented Neural Network \(PANN\)OptimizerAdamW \(Optax\)Weight constraintsProjected Gradient DescentLearning rate1×10−21\\times 10^\{\-2\}Remaining optimizer settingsOptax defaultsNumber of epochs10,000Hidden layer widths\[16,16\]\[16,16\]Output dimension1Hidden layer activationSoftplusOutput activationSoftplusConstitutive Artificial Neural Network \(CANN\)OptimizerAdamW \(Optax\)Weight constraintsProjected Gradient DescentLearning rate1×10−41\\times 10^\{\-4\}Remaining optimizer settingsOptax defaultsNumber of epochs10,000Generalized Invariant\-based CANN \(GI\-CANN\)OptimizerAdamW \(Optax\)Weight constraintsProjected Gradient DescentLearning rate1×10−31\\times 10^\{\-3\}Remaining optimizer settingsOptax defaultsNumber of epochs10,000