Stress representations for tensor basis neural networks: alternative
  formulations to Finger-Rivlin-Ericksen
Data-driven constitutive modeling frameworks based on neural networks and
classical representation theorems have recently gained considerable attention
due to their ability to easily incorporate constitutive constraints and their
excellent generalization performance. In these models, the stress prediction
follows from a linear combination of invariant-dependent coefficient functions
and known tensor basis generators. However, thus far the formulations have been
limited to stress representations based on the classical Rivlin and Ericksen
form, while the performance of alternative representations has yet to be
investigated. In this work, we survey a variety of tensor basis neural network
models for modeling hyperelastic materials in a finite deformation context,
including a number of so far unexplored formulations which use theoretically
equivalent invariants and generators to Finger-Rivlin-Ericksen. Furthermore, we
compare potential-based and coefficient-based approaches, as well as different
calibration techniques. Nine variants are tested against both noisy and
noiseless datasets for three different materials. Theoretical and practical
insights into the performance of each formulation are given.
Jan N. Fuhg, Nikolaos Bouklas, Reese E. Jones
Stress representations for tensor basis neural networks: alternative formulations to Finger-Rivlin-Ericksen

###### Abstract

Data-driven constitutive modeling frameworks based on neural networks and classical representation theorems have recently gained considerable attention due to their ability to easily incorporate constitutive constraints and their excellent generalization performance. In these models, the stress prediction follows from a linear combination of invariant-dependent coefficient functions and known tensor basis generators. However, thus far the formulations have been limited to stress representations based on the classical Rivlin and Ericksen form, while the performance of alternative representations has yet to be investigated. In this work, we survey a variety of tensor basis neural network models for modeling hyperelastic materials in a finite deformation context, including a number of so far unexplored formulations which use theoretically equivalent invariants and generators to Finger-Rivlin-Ericksen. Furthermore, we compare potential-based and coefficient-based approaches, as well as different calibration techniques. Nine variants are tested against both noisy and noiseless datasets for three different materials. Theoretical and practical insights into the performance of each formulation are given.

## 1 Introduction

Recently, there has been dramatically increased interest in machine learning (ML) in the computational sciences. This rise in popularity is due to: the ability of machine learning models to directly utilize experimental data in simulation environments, the potential speed up of ML models in comparison to traditional numerical models and methods, as well as the general utility and open-access ecosystem of ML tools. Nevertheless, many scientific ML (SciML) applications suffer from two interconnected bottlenecks: a lack of generalization capabilities due to poor extrapolations and a lack of trustworthiness due to the opaqueness of the trained models. The main premise in SciML is that the underlying data often comply with physical laws (known or yet to be discovered) or otherwise connect to known mathematical structure, which can help surmount the aforementioned bottlenecks via a physics-informed paradigm. The promise of SciML can lead to myriad benefits
such as: more accurate predictions, reduction of unnecessary human involvement, speed-up of the processing-performance-product development cycle, and minimization of the computational costs of detailed simulations. Particular to the focus of this work, an automated data-driven approach for constitutive modeling can have significant payoffs in material discovery, industrial engineering simulations and research. Many developments have been made in this arena for fluid closure models [1; 2; 3]; in this work we focus on constitutive models for solids.

A number of distinct approaches to forming constitutive models with ML have been investigated. ML tools have been utilized in parameter estimation of known constitutive models [4]. This is a task that becomes more complex as model parameters increase and experimental observations are limited. This is especially true for traditional optimization approaches due to the non-convex nature of the optimization problem at hand. Mixing traditional and ML approaches to representation and calibration via symbolic regression [5; 6; 7; 8; 9; 10] has been widely explored. This approach selects from a library of known models that directly enforce physical and mechanistic constraints (depending on the specific model choices) to distill parsimonious data-driven constitutive models. Notable developments include the approach of Wang _et al_. [11] who used reinforcement learning to turn model building into a competitive game. Also Schmidt _et al_. [12] used symbolic regression to distill constitutive laws from unlabeled data. Later, De Lorentzis and co-workers [13; 14] utilized sparse regression to discover interpretable constitutive laws for a wide array of material classes. An interesting extension to this work was the development of an unsupervised Bayesian framework for discovering hyperelasticity models which accounts for uncertainty [15]. Neural networks and Gaussian process models have been widely employed as replacements for human-selected, traditional model forms. In fact, the use of ML black-box constitutive models has been extensively studied for over 30 years. Starting from the influential works of Ghaboussi and collaborators [16; 17; 18], these tools have been employed for different material models with increasing complexity over the years [19; 20; 21; 22].

A significant current challenge is generating trustworthy models from low-data (constrained by experimental/computational cost) and limited-data (constrained by experimental design and observation). To this end efforts have been made to train data-driven constitutive models that do not only train with raw stress-strain data but incorporate additional physics-based restrictions to the trained model [23; 24; 25; 26; 27; 28]. These models, referred to as _physics-informed_ or _physics-guided data-driven_ constitutive models try to enforce a variety of physical principles and mechanics-informed assumptions. From enforcing objectivity, to material symmetries [1], thermodynamic constraints [29; 30] and polyconvexity [31; 32] there are approaches that enforce these condition weakly through the loss function [33; 34; 35] or strictly in the construction of the ML representation [1; 36; 37; 38]. A large majority of the proposed works in the literature for physics-guided constitutive models are based on neural networks [23; 39; 36; 24; 25; 26; 33] due to the flexibility of this paradigm.

Material frame indifference is a primary concern in developing constitutive models [40]. Ling _et al_. [1] introduced the tensor basis neural network (TBNN) to embed objectivity through an equivariant NN formulation. An anisotropic hyperelastic model was formed from the scalar invariants and tensor basis of the strain using atomistic crystal data, in addition to fluids applications. Later Frankel _et al_. [36] adapted the tensor basis representations to a Gaussian process formalism to represent general tensor functions and hyperelastic data. This was extended by Fuhg and Bouklas [37] to anisotropic materials, strictly enforcing known symmetries up to transverse isotropy; this work showed that this simplified learning approach led to significant generalization capabilities when the physics do not radically change outside of the training region. This approach was also utilized in Kalina _et al_. [41] integrated in a multiscale framework with automated data-mining. In
Fuhg _et al._[32], tensor basis NNs were utilized to discover the character of the anisotropy of the material through labeled data. Even though several works have focused on utilizing tensor basis representation theorems in learning of hyperelastic responses from labeled data pairs, there has not been an extensive study aimed at discovering the most efficient tensor basis representations for the learning tasks at hand in the context of finite deformation and hyperelasticity.

In the context of hyperelasticity, strict enforcement of polyconvexity requirements [42] for the strain energy density has also proven extremely useful towards generalization, discovery, and robustness. Input convex neural networks have been utilized for the enforcement of polyconvexity towards learning hyperelastic responses [31, 32], and in some cases even interpretability can be achieved [43] due to the non-parametric nature of the specific implementation. Alternately, neural ordinary differential equations have also been utilized towards strict enforcement of polyconvexity [44]. More recently Linden _et al._[38] presents a thorough review of techniques to enforce physical constraints and mechanistic assumptions towards learning hyperelasticity with NNs. Such approaches are crucial for the efficient utilization of the data and the development of robust material models that can efficiently generalize.

This work provides a limited survey of the wide variety of tensor basis techniques and contrasts their performance on representative data in the low-data regime (100 training points). We focus on stress representations for hyperelastic materials since they are the fundamental basis for finite deformation mechanics. The contributions of this work are: novel formulations based on the variety that the tensor basis framework affords, exploration of different methods of calibrating the models to data, and demonstration of the effects of noise and incompatible representations on physics-constrained formulations. To this end, we utilize well-known hyperelastic models as data generators.

In Sec. 2 we develop a multitude of equivariant _tensor basis_ neural networks (TBNNs) [1] formulations from classical representation theory. Then in Sec. 3 we give details of the data generation and training methodology. Sec. 4 presents the results of testing the models in and out of distribution and without and with additive noise. Finally in Sec. 5 we summarize the findings and conclude with avenues for future work.

## 2 Stress representations

In this work, we develop and compare a variety of tensor basis neural network (TBNN) formulations for stress representations. In this section, we introduce the fundamental differences between the representations and the neural network formulations that follow directly from the representations.

### Tensor basis models

Hyperelasticity is the prevailing theory for the description of finite deformation solid mechanics for continua in the absence of inelastic phenomena. The theory posits a potential \(\Psi\) from which the second Piola-Kirchhoff stress \(\mathbf{S}\) can be derived:

\[\mathbf{S}=2\boldsymbol{\partial}_{\mathbf{C}}\Psi(\mathbf{C})\, \tag{1}\]

as a function of the right Cauchy-Green stretch tensor \(\mathbf{C}=\mathbf{F}^{T}\mathbf{F}\). Here \(\mathbf{F}=\partial_{\mathbf{X}}\boldsymbol{\chi}\) is the deformation gradient of the spatial position \(\mathbf{x}=\boldsymbol{\chi}(\mathbf{X},t)\) at time \(t\) with respect to the corresponding reference position \(\mathbf{X}\) of the material. This potential ensures deformations are reversible, and is also utilized in some incremental formulations of large strain plasticity [45, 46].
In this work, we limit the discussion to isotropic hyperelasticity. In this case material frame invariance of the potential leads to the reduction of the inputs of \(\Psi\) to three scalar invariants \(I_{a}\) of \(\mathbf{C}\) and an equivariant stress function:

\[\mathbf{S}=2\,\boldsymbol{\partial}_{\mathbf{C}}\Psi(I_{1}(\mathbf{C}),I_{2}( \mathbf{C}),I_{3}(\mathbf{C})) \tag{2}\]

The chain rule results in the summation of material-specific, scalar derivative functions and an _a priori_ known tensor basis:

\[\mathbf{S}=2\boldsymbol{\partial}_{\mathbf{C}}\Psi(\mathbf{C})=2\,\sum_{a=1}^{ 3}\boldsymbol{\partial}_{I_{a}}\Psi\,\boldsymbol{\partial}_{\mathbf{C}}I_{a} \tag{3}\]

Typically the principal invariants,

\[I_{1}=\mathrm{tr}(\mathbf{C}),\quad I_{2}=\mathrm{tr}\big{(}\mathbf{C}^{-1} \big{)}\det(\mathbf{C}),\quad I_{3}=\det(\mathbf{C})\, \tag{4}\]

from the Cayley-Hamilton theorem

\[\mathbf{C}^{3}-I_{1}\mathbf{C}^{2}+I_{2}\mathbf{C}-I_{3}\mathbf{I}=\mathbf{0} \tag{5}\]

are employed. Note the second invariant is equivalently \(I_{2}=\frac{1}{2}(\mathrm{tr}\big{(}\mathbf{C})^{2}-\mathrm{tr}\big{(} \mathbf{C}^{2}\big{)})\). A three-term formula for the stress

\[\mathbf{S}=2\,[\partial_{I_{1}}\Psi+I_{1}\partial_{I_{2}}\Psi]\,\mathbf{I}-2 \partial_{I_{2}}\Psi\,\mathbf{C}+2I_{3}\,\partial_{I_{3}}\Psi\,\mathbf{C}^{-1} \tag{6}\]

comes from collecting terms with like powers of \(\mathbf{C}\). This is a well-known and arguably the most widely used stress representation for isotropic materials. It was first introduced by Finger [47] but was further popularized by Rivlin and Ericksen [48].

A generalization of this representation can be compactly written as a tensor basis expansion

\[\mathbf{S}=\sum_{a=1}^{3}c_{a}(\mathcal{I})\,\mathbf{B}_{a}\, \tag{7}\]

where the 3 coefficients \(c_{a}\) are functions of a set of 3 independent invariants \(\mathcal{I}\) and the basis \(\mathbf{B}_{a}\) must span \(\{\mathbf{C}^{a},a=-1,0,1\}\). For instance, Eq. (6) can be expressed as

\[\mathcal{I} = \{I_{a},\,a=1,2,3\} \tag{8}\] \[\mathcal{B} = \{\mathbf{C}^{a},\,a=0,1,-1\} \tag{9}\]

and

\[\begin{array}{rcl}c_{0}&=&2\,[\partial_{I_{1}}\Psi+I_{1}\partial_{I_{2}} \Psi]\\ c_{1}&=&-2\,\partial_{I_{2}}\Psi\\ c_{-1}&=&2I_{3}\,\,\partial_{I_{3}}\Psi\end{array} \tag{10}\]

Note that the Cayley-Hamilton theorem Eq. (5) allows the power basis to be shifted to higher or lower powers

\[\mathbf{C}^{3+k}=I_{1}\mathbf{C}^{2+k}-I_{2}\mathbf{C}^{1+k}+I_{3}\mathbf{C}^ {k}\text{ for }k\in\{\ldots,-2,-1,0,1,2,\ldots\}\, \tag{11}\]

for example

\[\mathcal{B}=\{\mathbf{C}^{a},\,a=0,1,2\}\, \tag{12}\]
via

\[I_{3}\mathbf{C}^{-1}=\mathbf{C}^{2}-I_{1}\mathbf{C}+I_{2}\mathbf{I}\;. \tag{13}\]

This basis together with the principal invariants (4) is another form of the Rivlin-Ericksen representation [48]. Also, the basis that results from the chain rule:

\[c_{a} = \{\partial_{I_{a}}\Psi,\,a=1,2,3\} \tag{14}\] \[\mathcal{B} = \{\boldsymbol{\partial}_{\mathbf{C}}I_{a},\,a=1,2,3\} \tag{15}\]

is part of an equally valid representation.

To calibrate Eq. (7), the model output can be regressed directly to stress data, or the coefficients for a given basis, e.g. \(\mathcal{B}=\{\mathbf{I},\mathbf{C},\mathbf{C}^{-1}\}\), can be determined at each data point \((\mathbf{C}_{i},\mathbf{S}_{i})\) via:

\[\begin{bmatrix}c_{1}\\ c_{2}\\ c_{3}\end{bmatrix}=\begin{bmatrix}1&\epsilon_{1}&\epsilon_{1}^{-1}\\ 1&\epsilon_{2}&\epsilon_{2}^{-1}\\ 1&\epsilon_{3}&\epsilon_{3}^{-1}\end{bmatrix}^{-1}\begin{bmatrix}\sigma_{1}\\ \sigma_{2}\\ \sigma_{3}\end{bmatrix} \tag{16}\]

using the fact that any power basis, such as Eq. (9), is collinear with \(\mathbf{S}\). Here \(\sigma_{a}\) and \(\epsilon_{a}\) are the eigenvalues of the stress and stretch tensors, herein \(\mathbf{S}\) and \(\mathbf{C}\), respectively. If the eigenvalues are distinct, Eq. (16) provides a unique solution for the coefficient values; however, multiplicity of strain eigenvalues requires special treatment, see Refs. [49, 50, 30] and App. A, which also outlines alternate solution procedures. Alternatively, we can use the Gram-Schmidt procedure

\[\mathbf{B}_{a}=\tilde{\mathbf{B}}_{a}-\sum_{b=1}^{a-1}\frac{\tilde{\mathbf{B} }_{a}:\mathbf{B}_{b}}{\mathbf{B}_{b}:\mathbf{B}_{b}} \tag{17}\]

to orthogonalize the basis \(\tilde{\mathbf{B}}_{a}\in\{\mathbf{I},\mathbf{C},\mathbf{C}^{-1}\}\), which results in

\[\mathcal{B} = \left\{\mathbf{I},\,\mathrm{dev}(\mathbf{C}),\,\mathrm{dev}( \mathbf{C}^{-1})-\left[\mathbf{C}^{-1}:\frac{\mathrm{dev}(\mathbf{C})}{\| \,\mathrm{dev}(\mathbf{C})\|}\right]\frac{\mathrm{dev}(\mathbf{C})}{\|\, \mathrm{dev}(\mathbf{C})\|}\right\} \tag{18}\]

if we keep the same scalar invariants. Herein \(\mathrm{dev}(\mathbf{C})=\mathbf{C}-1/3\,\mathrm{tr}(\mathbf{C})\mathbf{I}\). The fact that the Gram-Schmidt procedure starting with \(\mathbf{B}_{1}=\mathbf{I}\) and \(\mathbf{B}_{2}=\mathbf{C}\) leads to a spherical-deviatoric split is noteworthy. Orthogonality of the basis allows for direct determination of the coefficients:

\[c_{a}=\mathbf{S}:\mathbf{B}_{a} \tag{19}\]

Likewise, Gram-Schmidt applied to \(\{\mathbf{C}^{a},a=0,1,2\}\) gives the unnormalized basis

\[\mathcal{B} = \left\{\mathbf{I},\,\mathrm{dev}(\mathbf{C}),\,\|\,\mathrm{dev} (\mathbf{C})\|^{2}\,\mathrm{dev}(\mathbf{C}^{2})-(\mathbf{C}^{2}:\mathrm{dev} (\mathbf{C}))\,\mathrm{dev}(\mathbf{C})\right\} \tag{20}\]

Similarly, we can use a formulation inspired by the work Criscione _et al._[51] which effects an orthogonal spherical-deviatoric split of the basis via invariants:

\[\mathcal{I} = \left\{K_{1}=\mathrm{tr}(\mathbf{C}),\,K_{2}=\|\,\mathrm{dev}( \mathbf{C})\|,\,K_{3}=\det\left(\frac{\mathrm{dev}(\mathbf{C})}{\|\,\mathrm{ dev}(\mathbf{C})\|}\right)\right\} \tag{21}\] \[\mathcal{B} = \left\{\mathbf{I},\mathbf{A},-\frac{1}{3}\mathbf{I}-\mathrm{tr} \big{(}\mathbf{A}^{3}\big{)}\mathbf{A}+\mathbf{A}^{2}\right\} \tag{22}\]
where \({\bf A}={\rm dev}({\bf C})/\|\,{\rm dev}({\bf C})\|\). The resulting stress representation is

\[{\bf S} = \partial_{K_{1}}\Psi\,{\bf I}+\partial_{K_{2}}\Psi\,{\bf A}+ \partial_{K_{3}}\Psi\,\frac{1}{K_{2}}\left[-\frac{1}{3}{\bf I}-{\rm tr}\big{(}{ \bf A}^{3}\big{)}{\bf A}+{\bf A}^{2}\right]\] \[= \underbrace{\left[\partial_{K_{1}}\Psi-\frac{1}{3K_{2}}\, \partial_{K_{3}}\Psi\!\right]}_{c_{0}}\,{\bf I}+\underbrace{\left[\partial_{K _{2}}\Psi-\frac{3K_{3}}{K_{2}}\,\partial_{K_{3}}\Psi\right]}_{c_{1}}\,{\bf A} +\underbrace{\frac{1}{K_{2}}\,\partial_{K_{3}}\Psi}_{c_{2}}\,{\bf A}^{2}\]

Note Criscione _et al._[51] formulate the representation in terms of the spatial Hencky stretch, and here we apply the invariants of the same form to \({\bf C}\). This formulation combines derivative connection of the potential, invariants and the basis in the sense that the basis is a result of the choice of invariants, as in Eq. (14), and orthogonality of the basis. A better behaved set of related invariants

\[{\cal I}=\big{\{}{\rm tr}({\bf C}),\,{\rm tr}\,{\rm dev}({\bf C})^{2},\,{\rm det }\,({\rm dev}({\bf C}))\big{\}} \tag{24}\]

which eliminate the normalization in Eq. (21), leads to the basis

\[{\cal B}=\{{\bf I},2\,{\rm dev}\,{\bf C},({\rm det}({\rm dev}\,{\bf C}))\,{\rm dev }(({\rm dev}({\bf C}))^{-1})\} \tag{25}\]

See App. B for further details on the construction of an orthogonal basis.

With any of these representations, a densely connected neural network (NN) can be employed as a representation of the potential \(\Psi({\cal I})\) itself or the coefficient functions \(c_{a}({\cal I})\) directly. Summation of the coefficients with the known basis \({\cal B}\), as in Eq. (7), completes the formulation of a _tensor basis neural network_ (TBNN) [1]. Sec. 3.3, App. C and App. D provide details of the implementation of the TBNNs.

### Additional physical constraints

Other fundamental considerations, in addition to equivariance of the stress \({\bf S}({\bf C})\), constrain the form of the coefficient functions \(c_{a}({\cal I})\). Of the various constraints (rank-1 convexity, strong ellipticity, Hadamard stability, _etc._[40, 52, 53]), polyconvexity was proved by Ball [42] to ensure the existence of solutions in the context of hyperelasticity. For isotropic materials, polyconvexity requires that \(\Psi\) is convex in the triplet \(({\bf F},{\rm cof}\,{\bf F},{\rm det}\,{\bf F})\) which can be fulfilled when

\[\Psi=\Psi(I_{1},I_{2},I_{3})\mbox{ is convex in each of its arguments \@@cite[cite]{[\@@bibref{}{Ball}{}{}]}.} \tag{26}\]

An input convex neural network (ICNN) [54] satisfies these conditions and has been utilized for modeling hyperelastic materials in various recent studies [31, 38, 27]. Alternatively, if we assume that \(\Psi\) is polyconvex, we know that the derivatives of \(\Psi\) have to be non-decreasing, _i.e._\(\partial_{I_{a}}\Psi\geq 0\) with regards to \(I_{a}\). Assuming a representation of the form of Eq. (3)

\[{\bf S}=2\,\sum_{a=1}^{3}\mathbf{\partial}_{I_{i}}\Psi\,\mathbf{\partial}_{{\bf C}}I_{a}=\underbrace{\mathbf{\partial}_{I _{1}}\Psi}_{c_{0}}\,{\bf I}+\underbrace{\mathbf{\partial}_{I_{2}} \Psi}_{c_{1}}\,(I_{1}{\bf I}-{\bf C})+\underbrace{I_{3}\mathbf{\partial }_{I_{3}}\Psi}_{c_{-1}}\,{\bf C}^{-1}\,\,, \tag{27}\]

this implies that \(c_{1}(I_{1},I_{2}^{0},I_{3}^{0})\) is monotonically increasing in \(I_{1}\) for fixed \(I_{2}^{0}\) and \(I_{3}^{0}\). Note that the basis element \(I_{1}{\bf I}-{\bf C}=-\,{\rm dev}\,{\bf C}\) naturally arises from the Cayley-Hamilton/principal invariants, c.f. Eq. (4) and Eq. (10). We enforce this condition via an input monotone (or in fact monotonically
non-decreasing) neural network [27] which guarantees that the outputs of a neural network are monotonically non-decreasing in each of its inputs. Since to the best of our knowledge, no currently proposed neural network architecture enforces that each output individually is monotonically non-decreasing to only a subset of its inputs, and proposing a network of this kind is out of the scope of this work, we remark that this is an overconstrained way of enforcing the convexity condition.

Additional constitutive constraints resulting from mechanistic assumptions, include that the stress in the reference configuration is zero,

\[\mathbf{S}(\mathbf{C}=\mathbf{I})=\mathbf{0}\ \text{implies}\ \sum_{a}c_{a}(I_{1}^{0},I_{2}^{0},I_{3}^{0})=0. \tag{28}\]

with \(I_{1}^{0}=3,I_{2}^{0}=3,I_{3}^{0}=1\). One possible solution to enforcing this is to refactor the basis to form a Saint-Venant-like expansion:

\[\mathbf{S}=\sum_{a=1}^{3}c_{a}\mathbf{E}^{a} \tag{29}\]

where \(\mathbf{E}=1/2(\mathbf{C}-\mathbf{I})\). A \(\mathbf{C}\) based version is likewise:

\[\mathbf{S}=\sum_{a=1}^{3}c_{a}\mathbf{C}^{a}. \tag{30}\]

Note the coefficient functions \(c_{a}\) for these two representations are distinct but related, as are all the other representations introduced in this section. The requirements at the reference state \(\mathbf{C}=\mathbf{I}\) can be seen as a special case of the more general condition of symmetric loading where 2 or 3 of the eigenvalues are equal, examples include equibiaxial and hydrostatic/volumetric loadings.

The set of points where the eigenvalues are unique is dense in the invariant input space [55, 56], whereas highly symmetric cases are often used in testing and experiments since they are more easily understood and yet are sparse in the invariant input space. Since the unique case is dense there are continuous extensions for the coefficient functions to the case of eigenvalue multiplicity; however, the formula for the solution of the coefficients Eq. (16) does not provide them since the determinant of the system goes to zero. Although not well-cited, the important body of theoretical work starting with Serrin [55, 56, 57, 58] relates the smoothness of \(\mathbf{S}(\mathbf{C})\) or \(\mathbf{S}(\mathbf{E})\) to the smoothness of the coefficient functions with respect to the scalar invariants. Since most classical work treated only polynomial functions of the invariants, these developments have not been fully utilized; however, in the present context, we are forming general coefficient functions with neural networks. Man [56] proved that \(\mathbf{S}\) needs to be two degrees more continuous than the desired degree of smoothness of the coefficient functions, in particular \(\mathbf{S}(\mathbf{C})\) needs to be twice differentiable for \(c_{a}(\mathcal{I})\) to be continuous. Note that smooth solutions to the balance of momentum already require \(\mathbf{S}\) to be \(C^{1}\) and \(\Phi\) to be \(C^{2}\). Also, Scheidler [57] provided coefficient values from derivatives of the stress with respect to particular deformations, unlike Eq. (14).

Smoothness and growth considerations affect the choice of NN activations. For example, the St. Venant-like basis (29) incurs certain growth and asymptotic behavior. Refactoring the coefficients as

\[\tilde{c}_{a}=\|\mathbf{E}\|^{n}c_{a} \tag{31}\]

can enforce asymptotic behavior near \(\mathbf{C}\rightarrow\mathbf{I}\) as in Ref. [29]. The orthonormal basis formulations also need special consideration due to the normalization which creates non-smoothness in the coefficients, as in Eq. (18). An unnormalized basis, such as Eq. (20), avoids these issues.
### Summary of selected stress representations

In Sec. 4 we compare a number of distinct formulations of TBNNs for hyperelastic response listed in Table 1. Three are based on representing the strain energy potential \(\Psi(\mathcal{I})\) directly: (a) using the principal invariants \(\mathcal{I}\) as inputs to a standard feed-forward dense neural network (_Rivlin-Pot_), (b) using the principal invariants with an input convex neural network (_Convex-Pot_), and (c) using spherical-deviatoric split invariants in a standard dense neural network (_Crisc-Pot_). For these models the derivative of the potential with respect to these invariants through automatic differentiation provides the stress response. Six other models are based on coefficient-basis product formulations: (a) the customary power basis and coefficient functions in terms of the principal invariants (_Rivlin-Coeff_), (b) an input monotone neural network formulation of the coefficient functions with the power basis (_Mono-Coeff_) (c) the orthogonal basis with the Criscione invariants (_Crisc-Coeff_), (d) the orthogonal basis with the principal invariants (_Orthnorm-Coeff_), (e) an unnormalized orthogonal basis with more regular invariants (_Orth-Coeff_), and (f) a St.Venant-like basis with the principal invariants (_StV-Coeff_). For these both the coefficient functions \(c_{a}(\mathcal{I})\) and the basis \(\mathbf{B}_{a}\) are chosen. Table 1 summarizes the differences in the TBNN variants. In addition to these variations, we also explored how the method for calibration to stress data, e.g. via the coefficients found through regression or projection, or implicitly through direct calibration to stress, affects the model accuracy.

## 3 Data and training

For this study, we train the various NN models enumerated in Table 1 to stress data generated with classical hyperelastic models. In this section, we briefly discuss the classical data-generating models and give a detailed description of the data-generation process.

### Data and training

We remark that the complexity of the coefficient and potential functions is intrinsically connected to the stress measure, the basis function, and the invariants. To emphasize this consider the second Piola-Kirchhoff stress given by

\[\mathbf{S}=c_{0}^{C}\mathbf{I}+c_{1}^{C}\mathbf{C}+c_{-1}^{C}\mathbf{C}^{-1}. \tag{32}\]

\begin{table}
\begin{tabular}{|c|c c c|c c c|} \hline  & invariants & basis & coefficients & potential & convex & orthogonal basis \\ \hline Rivlin-Pot & Eq. (4) & Eq. (15) & Eq. (14) & \(\times\) & & \\ Convex-Pot & Eq. (4) & Eq. (15) & Eq. (14) & \(\times\) & \(\times\) & \\ Crisc-Pot & Eq. (21) & Eq. (22) & Eq. (23) & \(\times\) & & \(\times\) \\ Rivlin-Coeff & Eq. (4) & Eq. (12) & Eq. (7) & & & \\ Mono-Coeff & Eq. (4) & Eq. (27) & Eq. (27) & & \(\times\) & \\ Crisc-Coeff & Eq. (21) & Eq. (22) & Eq. (7) & & & \(\times\) \\ Orthnorm-Coeff & Eq. (21) & Eq. (18) & Eq. (7) & & & \(\times\) \\ Orth-Coeff & Eq. (24) & Eq. (20) & Eq. (7) & & & \(\times\) \\ StV-Coeff & Eq. (4) & Eq. (29) & Eq. (7) & & & \\ \hline \end{tabular}
\end{table}
Table 1: Tensor basis neural network variants.
Naively, one could presume that using the Kirchhoff stress tensor \(\mathbf{\tau}\) and the left Cauchy-Green tensor \(\mathbf{B}\) with an equivalent basis representation (\(\mathbf{I}\), \(\mathbf{B}\), \(\mathbf{B}^{-1}\)), i.e.

\[\mathbf{\tau}=c_{0}^{B}\mathbf{I}+c_{1}^{B}\mathbf{B}+c_{-1}^{B}\mathbf{B}^{-1} \tag{33}\]

the respective coefficients might be the same, e.g. \(c_{a}^{C}=c_{a}^{B}\) for \(a=0,1,-1\). However, recalling that the Kirchhoff stress can be expressed as \(\mathbf{\tau}=\mathbf{FSF}^{T}\), Eq. (32) can be rewritten as

\[\mathbf{\tau}=c_{0}^{C}\mathbf{B}+c_{1}^{C}\mathbf{B}^{2}+c_{-1}^{C}\mathbf{I}. \tag{34}\]

Hence, under the assumption that the eigenvalues are unique, we find that

\[c_{0}^{B}\mathbf{I}+c_{1}^{B}\mathbf{B}+c_{-1}^{B}\mathbf{B}^{-1}=c_{0}^{C} \mathbf{B}+c_{1}^{C}\mathbf{B}^{2}+c_{-1}^{C}\mathbf{I} \tag{35}\]

which yields

\[c_{0}^{B}=c_{-1}^{C}-c_{1}^{C}I_{2},\qquad c_{1}^{B}=c_{0}^{C}+c_{1}^{C}I_{1}, \qquad c_{-1}^{B}=c_{1}^{C}I_{3}. \tag{36}\]

via the Cayley-Hamilton theorem (5). The complexity of the two sets of coefficient functions \(\{c_{0}^{C},c_{1}^{C},c_{-1}^{C}\}\) and \(\{c_{0}^{B},c_{1}^{B},c_{-1}^{B}\}\) is therefore clearly different. Using the Cayley-Hamilton theorem (5) to transform the model representation would also alter the complexity of the coefficient functions. In order to make the following comparisons as fair as possible we have restricted ourselves to second Piola-Kirchhoff stress representations and data. Note that the Piola transform would also affect the orthogonality of the basis.

We furthermore remark that additively separable energies that are based on the Valanis-Landel hypothesis [59, 60, 61] lead to more trivial calibrations, i.e. if

\[\Psi(I_{1},I_{2},I_{3})=\Psi_{1}(I_{1})+\Psi_{2}(I_{2})+\Psi_{3}(I_{3}) \tag{37}\]

we can see from Eq. (10) that this would result in

\[\begin{split} c_{0}(I_{1},I_{2})&=2\left[\partial_{ I_{1}}\Psi_{1}(I_{1})+I_{1}\partial_{I_{2}}\Psi_{2}(I_{2})\right]\\ c_{1}(I_{2})&=-2\,\partial_{I_{2}}\Psi_{2}(I_{2}) \\ c_{-1}(I_{3})&=2I_{3}\partial_{I_{3}}\Psi_{3}(I_{3}).\end{split} \tag{38}\]

Hence, this leads to \(c_{1}\) and \(c_{-1}\) being functions of only one invariant and \(c_{0}\) reduced to a function of two invariants. In order to avoid these simplifications we use only hyperelastic models that are not additively decomposable with regards to their inputs.

Note the definition of the invariants can be engineered to reduce the complexity of the coefficient functions for a particular material dataset. In a limiting case, the coefficient functions are themselves invariants and hence present the simplest representation in some sense, albeit one that is hard to discover _a priori_ from the measured data. Representation complexity is particularly important in the low data regime which we explore.

### Data models

Three well-known compressible hyperelastic models were selected to generate training data: (a) Mooney-Rivlin [62, 63], (b) a modified version of Carroll's hyperelastic law [64, 65], and (c) a Gent-type model [66, 67]. Each is expressed in terms of the invariants \(I_{1}=\mathrm{tr}(\mathbf{C})\), \(I_{2}=\mathrm{tr}\big{(}\mathbf{C}^{-1}\big{)}\det(\mathbf{C})\), and \(J=\sqrt{\det\mathbf{C}}\)
The specific compressible Mooney-Rivlin model considered here has the strain energy function

\[\Psi=\theta_{1}\left(\frac{I_{1}}{J^{2/3}}-3\right)+\theta_{2}\left(\frac{I_{2}}{ J^{4/3}}-3\right)+\theta_{3}(J-1)^{2}\, \tag{39}\]

which yields a second Piola-Kirchhoff stress of the form:

\[\mathbf{S}=\underbrace{2\left(\frac{\theta_{1}}{J^{2/3}}+I_{1}\frac{\theta_{2 }}{J^{4/3}}\right)}_{c_{0}^{*}}\mathbf{I}\underbrace{-2\frac{\theta_{2}}{J^{4 /3}}}_{c_{1}^{*}}\mathbf{C}+\underbrace{J\left[-\frac{2}{3}\theta_{1}\frac{I_{ 1}}{J^{5/3}}-\frac{4}{3}\theta_{2}\frac{I_{2}}{J^{7/3}}+2\theta_{3}(J-1)\right] }_{c_{-1}^{*}}\mathbf{C}^{-1} \tag{40}\]

We use (scaled) material parameters ( \(\theta_{1}=0.92\) Pa, \(\theta_{2}=2.37\) Pa and \(\theta_{3}=10.001\) MPa) from fits to vulcanized rubber data, c.f. Ref. [68], for data generation.

Following Ref. [65], a modified Carroll model is defined by the strain energy function

\[\Psi=\theta_{1}\left(\frac{I_{1}}{J^{2/3}}-3\right)+\theta_{2}\left[\left( \frac{I_{1}}{J^{2/3}}\right)^{4}-81\right]+\theta_{3}\left(\sqrt{\frac{I_{2}} {J^{4/3}}}-\sqrt{3}\right)+\theta_{4}(J-1)^{2}. \tag{41}\]

This energy results in

\[\mathbf{S} =\underbrace{2\left[\frac{\theta_{1}}{J^{2/3}}+4\theta_{2}\left( \frac{I_{1}}{J^{2/3}}\right)^{3}+I_{1}\left(\theta_{2}\frac{1}{2J^{2/3}I_{2}} \right)\right]}_{c_{0}^{*}}\mathbf{I}\underbrace{-2\theta_{2}\frac{1}{2J^{2/3 }I_{2}}}_{c_{1}^{*}}\mathbf{C} \tag{42}\] \[+\underbrace{J\left(-\frac{2}{3}\theta_{1}\frac{I_{1}}{J^{4/3}}- \frac{8}{3}\theta_{2}\frac{I_{1}^{4}}{J^{11/3}}-\frac{2}{3}\theta_{3}\frac{ \sqrt{I_{2}}}{J^{5/3}}+\theta_{4}(J-1)\right)}_{c_{-1}^{*}}\mathbf{C}^{-1}\]

We use a scaled version of the material parameters reported in Ref. [65], in particular \(\theta_{1}=151.09387\) GPa, \(\theta_{2}=0.3028\) MPa, \(\theta_{3}=68.33070\) GPa, and \(\theta_{4}=500\) TPa.

Lastly, we also utilize the response of a compressible version of the Gent+Gent model, as named in Ref. [69], that is defined by the strain energy function

\[\Psi=-\frac{\theta_{1}}{2}J_{m}\log\left(1-\frac{I_{1}-3}{J_{m}}\right)-\theta _{2}\log\left(\frac{I_{2}}{J}\right)+\theta_{3}\left(\frac{1}{2}(J^{2}-1)-\log J \right)\, \tag{43}\]

where we choose \(\theta_{1}=2.4195\) MPa, \(\theta_{2}=1.8146\) MPa, \(\theta_{3}=1.2097\) MPa and \(J_{m}=77.931\). This strain energy yields a second Piola-Kirchhoff stress of the form:

\[\mathbf{S}=\underbrace{2\left[-\frac{\theta_{1}}{2}J_{m}\frac{1}{I_{1}-3-J_{m }}+\theta_{2}\frac{I_{1}}{I_{2}}\right]}_{c_{0}^{*}}\mathbf{I}\underbrace{-2 \theta_{2}\frac{1}{I_{2}}}_{c_{1}^{*}}\mathbf{C}+\underbrace{J\,\theta_{3}(J- \frac{1}{J})}_{c_{-1}^{*}}\mathbf{C}^{-1} \tag{44}\]

The Gent+Gent model is not polyconvex; however, it is convex over a limited range where \(I_{1}<J_{m}+3\). For simplicity, we refer to this model simply as Gent.

Note that hereafter \(c_{a}^{*}\) denote the true coefficients, which differ from the extracted coefficients \(c_{a}\) near ill-conditioned solves, and the fitted NN coefficients \(\hat{c}_{a}\).
### Training and validation

For sampling, we define a nine-dimensional space around the undeformed configuration of the deformation gradient as

\[\overline{F}_{ij}\in[F_{ij}^{L},F_{ij}^{U}]=\begin{cases}1-\delta\leq 1\leq 1+ \delta,&\text{when i=j},\\ -\delta\leq 0\leq\delta,&\text{otherwise}.\end{cases} \tag{45}\]

We then define a training region with \(\delta=0.2\) and a test region with \(\delta=0.3\) and use the space-filling sampling technique proposed in Ref. [37] to generate \(100\) training points and \(10,000\) test points that fill their respective spaces. Figure 1 shows the spread of these samples in invariant space \((I_{1},I_{2},J)\). Then given the triple \((I_{1},I_{2},I_{3}=J^{2})\) we can reconstruct the right Cauchy-Green tensor as

\[\mathbf{C}=\begin{bmatrix}\frac{1}{3}I_{1}-2\sqrt{H}\cos\left(\frac{\pi-\beta }{3}\right)&0&0\\ 0&\frac{1}{3}I_{1}-2\sqrt{H}\cos\left(\frac{\pi+\beta}{3}\right)&0\\ 0&0&\frac{1}{3}I_{1}-2\sqrt{H}\cos\left(\frac{\beta}{3}\right)\end{bmatrix} \tag{46}\]

where

\[H=\frac{1}{9}(I_{1}^{2}-3I_{2}),\qquad G=\frac{1}{3}I_{1}I_{2}-I_{3}-\frac{2}{ 27}I_{1}^{3},\qquad\beta=\arccos\left(-\frac{G}{2H^{3/2}}\right) \tag{47}\]

from which we can obtain the values of the invariants and the basis of all the investigated stress representations of Sec. 2. When training the model we use an \(80/20\) split to obtain \(20\) validation data points.

After obtaining a set of training and testing data we added noise to the resulting coefficient values to disrupt symmetries and analytic functional forms. In particular, we take the coefficient values \(c_{0},c_{1},c_{-1}\) corresponding to the representation

\[\mathbf{S}=c_{0}\mathbf{I}+c_{1}\mathbf{C}+c_{-1}\mathbf{C}^{-1} \tag{48}\]

for every data point, and define a noisy version \(\tilde{c}_{a}=c_{a}+\mathcal{N}(0,0.02\left|c_{a}\right|)\) for \(a=0,1,-1\) which then gives a noisy stress

\[\tilde{\mathbf{S}}=\tilde{c}_{0}\mathbf{I}+\tilde{c}_{1}\mathbf{C}+\tilde{c} _{-1}\mathbf{C}^{-1}. \tag{49}\]

This \(\tilde{\mathbf{S}}\) was then used as the target stress to obtain the coefficients for all other models, e.g. the Criscione model. Hence, we generate \(100\) noisy training samples that have the same invariants as the noiseless counterparts and use the same noiseless data for the test set. An example of a generated noisy test set is shown in Figure 2 for the Mooney-Rivlin model.

All the tensor basis neural network models [1] were implemented in _PyTorch_[70]. Potential models were formed from a multilayer, densely connected feedforward neural network (NN) with a single output

\[\Phi=N\!N(\mathcal{I}) \tag{50}\]

where the coefficients \(c_{a}=\partial_{I_{a}}\Phi\) are obtained through automatic differentiation. Summation with the known basis provides the stress

\[\mathbf{S}=\sum_{a}c_{a}\mathbf{B}_{a} \tag{51}\]
The coefficient-based models utilized a monolithic NN with 3 outputs

\[c_{a}=N\!\!N_{a}(\mathcal{I}) \tag{52}\]

and the same summation to form the stress. To be consistent all TBNN models consisted of 3 layers with 30 neurons per layer and a _Softplus_ activation function [71]. App. C and App. D provide additional details of the implementation of potential and coefficient-based TBNNs, respectively.

The training loss was formulated on the mean squared error of either the stress components or the coefficients

\[\text{MSE}=\lambda_{\mathbf{S}}\sum_{i=1}^{N_{\text{train}}}\|\mathbf{S}_{i}- \hat{\mathbf{S}}_{i}\|^{2}+\lambda_{c}\sum_{i=1}^{N_{\text{train}}}\sum_{a=1}^ {3}\|[c_{a}]_{i}-[\hat{c}_{a}]_{i}\|^{2}\, \tag{53}\]

since these are available from representative volume element (RVE)/experimental data, whereas the strain energy is less accessible. Although mixing both losses proved useful in preliminary studies, all reported data is from models trained with either \(\lambda_{\mathbf{S}}=1,\lambda_{c}=0\) or \(\lambda_{\mathbf{S}}=0,\lambda_{c}=1\). Note the coefficients scale as \(|c_{a}|\sim\|\mathbf{S}|\ \|\mathbf{B}_{a}\|^{-1}\) so the coefficient-based loss can suffer from numerical conditioning issues. All models were trained for \(50,000\) epochs using the _Adam_ optimizer [72] with a constant learning rate of \(0.001\).

We compared the performance of the models using the normalized mean squared error of the stress over the \(10,000\) testing data points

\[\text{RMSE}=\sqrt{\frac{\sum_{i=1}^{N_{\text{test}}}\|\mathbf{S}_{i}-\hat{ \mathbf{S}}(\mathbf{C}_{i})\|^{2}}{\sum_{i}\|\mathbf{S}_{i}\|^{2}}} \tag{54}\]

where \(\hat{\mathbf{S}}_{i}\) is the predicted stress and \(\mathbf{S}_{i}\) represents the ground truth at data point \(i\).

## 4 Results

First, we survey the test losses of the models enumerated in Table 1. Then we undertake detailed investigations of why the theoretically equivalent representations do or do not perform well by examining where the largest errors occur and how the predictions compare to held-out data.

Figure 1: Space-filling training/validation (colored by maximum stress component normalized by the mean for the Mooney-Rivlin data) and test samples in invariant space.
### Comparison of test losses

For each TBNN model described in Table 1 we assembled an ensemble of 30 parameterizations using random initialization of the NN parameters and shuffling the training/validation subsets of the 100 training points. Figure 3 shows the range of RMSE test errors for the six datasets from the 3 traditional models described in Sec.3, each with and without noise. Clearly, the various theoretically equivalent TBNN formulations perform differently and each of the datasets evoke different errors. Overall the polyconvex (Conv-Pot) and monotonic (Mono-Coeff) models appear to perform the best, although the standard potential-based (Rivlin-Pot) model has comparable performance to Mono-Coeff. The coefficient-based Rivlin-Coeff has considerably higher test errors than the potential-based Rivlin-Pot, despite Rivlin-Pot relying on automatic differentiation. The Criscione and other orthogonal models perform worse than the convex, monotonic and Rivlin models but they do better on the Gent data than on the other two datasets. We observe that the Gent model has a qualitatively different functional form than the selected Mooney-Rivlin and Carroll data-generating models, e.g. the presence of log terms in the energy Eq. (43). The Orth model with smoother invariants performs the best of the orthogonal basis models and is an anomaly in that it trains better indirectly to stress than to the extracted coefficients. This may be due to the conditioning issues with solving for coefficients of a power basis, mentioned in Sec. 3.3. The St. Venant model (both \(\mathbf{C}\) and \(\mathbf{E}\) based) is an outlier with large errors likely due to the mismatch in the growth of the tensor basis and the data, which necessitates more complex coefficient functions. Although small, the differences in calibration techniques can have an effect. Training to stress, instead of extracted coefficients, can regularize the trained coefficient functions, since stress is smoother than the coefficient functions as per the Man-Serrin theorem discussed in Sec. 2.2. Also training to stress can discourage a potential reinforcement of bias from individually trained coefficient functions that need to coordinate to form an accurate stress. Training to coefficients, on the other hand, removes the potentially ill-conditioned linear algebra implicit in training to stress. Projection of data onto expected bases may also remove discrepancies as with the noisy datasets.

The test errors seem to be largely dominated by testing the models in extrapolation, more insight will be given in the following sections.

Figure 2: Effective stress noise in percent of Frobenius norm error over the test data of a Mooney-Rivlin dataset.
### Locations of worst errors

The worst 1% errors of the 10,000 sample test set for each model are shown in Figure 4 for the Mooney-Rivlin data, in Figure 5 for the modified Carroll data, and in Figure 6 for the Gent data. For reference, the undeformed state is at \((I_{1}=3,I_{2}=3,I_{3}=1)\) which is inside the hull of sample points shown in these figures. Generally speaking, for most models, the worst errors are at the boundary of the test locus where they are forced into extrapolation. Note that \(I_{3}\) is associated with volumetric deformation, \(I_{1}\) can be interpreted at the linearization of \(I_{3}\), while \(I_{2}\) is sensitive to shear and deviatoric deformations.

Examining Figure 4, the Convex-Pot, Rivlin-Pot, and Orth-Coeff have the largest errors where the invariants (and eigenvalues of **C**) are large, while the Crisc-Coeff and the similar Orthnorm-Coeff have largest errors in the low region. The Crisc-Pot has relatively large errors at both extremes. Of the better-performing models, the Mono-Coeff formulation is an outlier since it has its worst errors in the midrange of the invariants, albeit still at the boundary. Likewise, the Orthnorm-Coeff has high error in the midrange, as well as the low range, while the St. Venant model performs particularly poorly in the low to midrange. The patterns are relatively unchanged with discrepancy

Figure 3: Comparison of TBNN variants summarized in Table 1 over test sets for the three different data models with and without noise. Note that coefficient-based variants with \({}^{*}\) were trained to stress data, and had to discover the coefficient values. Also the **E**-based variant of the St.Venant basis Eq. (29) had errors well above the displayed range so only the **C**-based variant Eq. (30) is shown.
added by noise, with the exception of the worst errors transitioning to the lower range for the most accurate model, Convex-Pot. The errors for the Carroll data shown in Figure 5 largely resemble those for the Mooney-Rivlin data, although for this data the Mono-Coeff and Orth-Coeff models seem more sensitive to noise. They, together with the Convex-Pot, shift their worst error locations with noise. The errors for the Gent data, shown in Figure 6, however, present different patterns. For this case, all models perform worst where the invariants are small, which can be ascribed to the Gent model being ill-defined where \(I_{1}\) becomes less than a value determined by the parameter \(J_{m}\). There are also scattered worst error locations in the midrange for the Rivlin-Coeff model.

Figure 7, Figure 8, and Figure 9 provide another view of the error patterns and corroborate the observations from the previous plots. These figures illustrate the correlation of the maximum errors with the difference in the largest and smallest (Rivlin-Ericksen) coefficient values for a particular data point in the test set. Large differences in the coefficient values are associated with less symmetric deformations. For each model, the figures show how the worst errors shift as a function of the difference in the coefficient values. For the Mooney-Rivlin data, only the worst-performing model, StV-C, has a single locus of maximal errors. Of the best-performing models, Convex-Pot, Rivlin-Pot, and Rivlin-Coeff have similar patterns that remain stable after the injection of noise. The Mono-Coeff model, on the other hand, changes the locations of where the worst errors occur relative to the difference in the coefficient values and also has the worst errors on par with Convex-Pot. Again the patterns for the Carroll data are similar to the Mooney-Rivlin data, while the Gent data present qualitatively different patterns. For the Gent data, all models, except the worst performing StV-C, have overlapping loci of worst errors that do not change appreciably with noise. This is perhaps due to discrepancies between the relatively simple TBNN models and the Gent data.

Next, we aim to test the performance of the TBNN variants to generalization. We conjecture that the differences in generalization performance can be attributed to how much of the complexity of the stress-strain mapping is intrinsically provided by the nonlinearity of the stress representation bases. Meaning, if the nonlinearity of the bases in the chosen representation is approaching the nonlinearity of the stress-strain mapping then the coefficient can be described by simpler lower-order functions. If the coefficients are lower-order functions, i.e. constant, linear, or quadratic, then accurately extrapolating this behavior with a TBNN is going to require less training data and a less complex NN to be accurate.

To check this hypothesis we offer the following approach. Consider the mapping \(\mathcal{M}_{IC}\) from the respective invariants to the respective coefficient values, e.g

\[\mathcal{M}_{IC}:\mathcal{I}\rightarrow[c_{1},c_{2},c_{3}]. \tag{55}\]

Let this mapping be approximated by a polynomial regressor of \(n\)-th order with interaction terms denoted by \(\hat{\mathcal{M}}^{n}_{IC}(\mathcal{I})\), e.g. for \(n=2\)

\[\hat{\mathcal{M}}^{2}_{IC}(\mathcal{I})=\hat{\mathbf{c}}^{2}_{IC}=\mathbf{\alpha}_{0} +\mathbf{\alpha}_{1}I_{1}+\mathbf{\alpha}_{2}I_{2}+\mathbf{\alpha}_{3}I_{3}+\mathbf{\alpha}_{4 }I_{1}^{2}+\mathbf{\alpha}_{5}I_{1}I_{2}+\mathbf{\alpha}_{6}I_{2}^{2}+\mathbf{\alpha}_{7}I _{2}I_{3}+\mathbf{\alpha}_{8}I_{3}^{2} \tag{56}\]

where \(\mathbf{\alpha}_{i}\in\mathbb{R}^{3}\). To check the potential complexity of this mapping we look at two scenarios. First, how polynomial regression fitted on the training data predicts the test data, and second, how well \(\hat{\mathcal{M}}^{n}_{IC}(\mathcal{I})\) predicts the test data coefficients if it was trained on the test data. The root-mean-squared error (RMSE) between the reference and predicted test data coefficients of the former is shown in Figure 10 for the three noiseless data cases. Surprisingly, polynomials of second order generalize the best for Rivlin-Coeff and Mono-Coeff. This leads us to the conjecture that the
complexity of the coefficient functions of Rivlin-Coeff and Mono-Coeff are generally lower than the other representation. This seems to correlate with the results of the TBNN generalization errors c.f. Figure 3. The coefficient error between the regressor of an increasing polynomial order, trained on the test data and evaluated on the test data is displayed in Figure 11. It can be seen that, compared to the previous case (Figure 3), Mono-Coeff and Rivlin-Coeff still have the lowest errors but more significantly that the complexity of the coefficient functions seems to have changed, i.e. while second order polynomials where best for models trained on the training data, now an increasing polynomial order seems to help to accurately fit the coefficient functions. Note that the upward trends after initial low discrepancy fits in some of the data in Figure 10 could possibly be attributed to overfitting of the polynomials to the 100-point training dataset.

Next, we aim to gauge what the contribution of the basis representation is on the accuracy. From the output of the polynomial regression of (56) an \(n\)-th order polynomial prediction of the stress can be obtained

\[\hat{\mathbf{S}}^{n}_{IC}=\sum_{a=1}^{3}\,\hat{c}^{n}_{a,IC}\mathbf{B}_{a}. \tag{57}\]

We furthermore assume an alternative mapping \(\mathcal{M}_{IS}\) from the invariants of the representation to

Figure 4: Mooney-Rivlin data, worst 1% errors for each model.
the symmetric components of the stress

\[\mathcal{M}_{IS}:\mathcal{I}\rightarrow[S_{11},S_{12},S_{13},S_{22},S_{23},S_{33}] \tag{58}\]

for which we build a similar \(n\)-th order polynomial referred to as \(\hat{\mathcal{M}}_{IS}^{n}\). Then, we can find the difference between the RMSEs of \(\hat{\mathbf{S}}_{IC}^{n}\) and \(\hat{\mathcal{M}}_{IS}^{n}\) evaluated on the test data, i.e.

\[\Delta\text{RMSE}^{n}(\mathcal{I}^{test})=\text{RMSE}(\hat{\mathcal{M}}_{IS}^ {n}(\mathcal{I}^{test}))-\text{RMSE}(\hat{\mathbf{S}}_{IC}^{n}(\mathcal{I}^{ test})).\]

Simplistically, this difference between the stress errors between these two regressors will help us judge the role and contribution of the basis generators, i.e. if \(\Delta\text{RMSE}^{n}>0\) then the bases have a positive contribution to the prediction which means that the RMSE of obtaining the stress from a linear combination of coefficients and bases \(\hat{\mathbf{S}}_{IC}^{n}=\sum_{a=1}^{3}\,\hat{c}_{a,IC}^{n}\mathbf{B}_{a}\) is lower than the mapping from invariants to stress directly. This would tell us that the basis components take complexity out of the system. On the other hand if \(\Delta\text{RMSE}^{n}<0\) the basis components make the mapping between invariants and stress more complex.

In Figure 12 we compare the RMSE-difference of models trained on the training data and evaluated on the test data while Figure 13 highlights the RMSE difference when the models were trained

Figure 5: Modified Carroll data, worst 1% errors for each model.
and tested on the test data. It can be seen that the basis components of the Mono-Coeff and Rivlin-Coeff representations generally help in reducing the complexity of the invariants-stress mapping, which is the reason why the invariants-stress mapping are not the same.

Figure 6: Gent data, worst 1% errors for each model.

Figure 7: Maximal errors for each model for the Mooney-Rivlin data. The horizontal lines indicate where the data has been clipped.
especially for lower polynomial order while the opposite is true for the remaining representations that were investigated.

Overall we believe that this investigation of the data through the lens of a polynomial regressor suggests that our hypothesis is valid. The TBNN has better generalization capabilities for some of the stress representations because the complexity of the mapping is reduced owing to nonlinearities introduced by the basis components.

Figure 8: Maximal errors for each model for the modified Carroll data. The horizontal lines indicate where the data has been clipped.

Figure 9: Maximal errors for each model for the Gent data. The horizontal lines indicate where the data has been clipped.
### Interpolation

The TBNN models have the ability to form smooth extensions to coefficient functions for high symmetry loading. Consider the parameterized invariants [73]

\[\begin{split} I_{1}(\gamma)&=3-2\,\gamma+\gamma^{2}\\ I_{2}(\gamma)&=3-4\gamma+2\,\gamma^{2}\\ I_{3}(\gamma)&=(1-\,\gamma)^{2}\end{split} \tag{59}\]

with \(\gamma\in[-0.2,0.2]\). A uniaxial extension can be observed for \(\gamma<0\) while \(\gamma>0\) yields an (equi)biaxial extension. This path is highlighted in Figure 14 in the projected invariant space and is inside the training region even though not explicitly part of the training data set. We specifically focus on this path due to the fact that it is characterized by two coalescent principal stretches \(\lambda_{2}=\lambda_{1}\) for \(\gamma<0\) and \(\lambda_{3}=\lambda_{2}\) for \(\gamma>0\) where \(\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}\). Here \(\lambda_{a}=\sqrt{\epsilon}_{a}\) are the principal stretches. As described earlier and as seen in App. A, this means that the solution matrix for the coefficients of some of the presented stress representations is ill-conditioned and special schemes are

Figure 11: RMSE between predicted coefficient values and reference coefficient values using a polynomial of \(n\)-th order as the regressor which was trained on all available data (test data) and tested on the same data.

Figure 10: RMSE between predicted coefficient values and reference coefficient values using a polynomial of \(n\)-th order as the regressor which was trained on the training data and tested on the test data set.
needed to be able to solve for the coefficients. We remark that:

1. Depending on the loading path, these schemes lead to discontinuities near the reference state \(\mathbf{C}=\mathbf{I}\). In particular, this is the case for the Rivlin-Coff (\(\mathbf{S}=c_{1}\mathbf{I}+c_{2}\mathbf{C}+c_{3}\mathbf{C}^{-1}\)) and St.V-Coeff (\(\mathbf{S}=c_{1}\mathbf{C}+c_{2}\mathbf{C}^{2}+c_{3}\mathbf{C}^{3}\)) representations for the path described in Eq. (59).
2. Due to the space-filling way the training data was generated the principal strains of all training data points are unique, apart from the undeformed configuration. This means that the trained models have not been trained on coefficients that came as a result of the special schemes described in App. A. Hence, examining the predicted coefficients in this high-symmetry case provides an interesting and interpretable test case.

Figure 15 shows the trained coefficients and stress predictions for the Rivlin-Coeff and St.V-Coeff models over \(\gamma\). Clearly, the extracted coefficient functions near the reference state \(\mathbf{C}=\mathbf{I}\) become discontinuous; however, the built-in continuity of the NN enables an approximate continuous

Figure 12: The RMSE difference between the stress prediction of an n-th order polynomial trained from the invariants to coefficients \(\hat{\mathbf{S}}^{n}_{IC}\) and an n-th order polynomial trained from the invariants directly to the stress \(\hat{\mathcal{M}}^{n}_{IS}\), i.e. \(\Delta\text{RMSE}^{n}(\mathcal{I}^{test})=\text{RMSE}(\hat{\mathcal{M}}^{n}_ {IS}(\mathcal{I}^{test}))-\text{RMSE}(\hat{\mathbf{S}}^{n}_{IC}(\mathcal{I}^{ test}))\). Both models were trained on the **training** data and tested on the **test** data.

Figure 13: The RMSE difference between the stress prediction of an n-th order polynomial trained from the invariants to coefficients \(\hat{\mathbf{S}}^{n}_{IC}\) and an n-th order polynomial trained from the invariants directly to the stress \(\hat{\mathcal{M}}^{n}_{IS}\), i.e. \(\Delta\text{RMSE}^{n}(\mathcal{I}^{test})=\text{RMSE}(\hat{\mathcal{M}}^{n}_ {IS}(\mathcal{I}^{test}))-\text{RMSE}(\hat{\mathbf{S}}^{n}_{IC}(\mathcal{I}^ {test}))\). Both models were trained on all available **test** data and tested on the same data.
extension. This approximation \(\hat{c}_{a}\) is different than the coefficients \(c_{a}\) extracted from the equation system, altered to accommodate multiplicity, but still yields an accurate stress representation for Rivlin-Coeff and a sufficient one for StV-Coeff. It seems that the smooth extension avoids large errors that would be incurred if the extracted coefficients were approximated.

Remarkably, the resulting predicted coefficients for Rivlin-Coeff are practically equivalent to the derived coefficients from the potential prediction of Convex-Pot, Figure 16. This is interesting since Convex-Pot has not been trained on the coefficients. We furthermore remark that (for this loading case) not all of the representations show a discontinuous coefficient behavior as a result of the coalescent principal strains. For example for Mono-Coff, whose coefficient matrix was also ill-conditioned and had to be adapted, the coefficient behavior is smooth over \(\gamma\), see Figure 17. In this case, the NN coefficients are basically identical to the extracted coefficients obtained from the altered equation system.

### Extrapolation

To highlight the predictive quality of selected representations consider a loading path given by the homogenous deformation [73]

\[\begin{split} I_{1}&=3-1.6\,\gamma+\gamma^{2}\\ I_{2}&=3-3.2\gamma+1.64\,\gamma^{2}\\ I_{3}&=(1-0.8\,\gamma)^{2}\end{split} \tag{60}\]

where \(\gamma\in[0,1]\). The loading path projected into invariant space is shown in Figure 18. We can see that it starts at the undeformed configuration and goes beyond even the range of the test data. The deformation leaves the hull of the training data points at \(\gamma\approx 0.25\). We highlight that this path is not explicitly part of the training data set. For five of the representations, Figure 19 shows the expected coefficients and trained coefficients respectively. The range of the training regime is highlighted. Surprisingly, all the selected models fit the reference coefficients sufficiently well, even St.V-Coeff which was by far the worst-performing approach in terms of generalization error. Comparing the fits of the models to the conclusions drawn from Figure 3, it becomes evident that the magnitude-wise largest coefficient predictions of the poorly performing representations tend to stagnate earlier. Similar reasoning can be followed when examining Figure 20 which plots three components of the expected stress and their predicted counterparts from the neural network models,

Figure 14: Uniaxial and Biaxial extension parameterization in invariant space given by \(\gamma\)
where we observe that the representations which have only one significant coefficient on the path are also the best performing.

## 5 Conclusion

In this work, we investigated the a wide variety of tensor basis neural network models, including previously unexplored alternatives to the classical Finger-Rivlin-Ericksen stress representation, such as TBNNs with an orthogonal basis and others with a St. Venant basis. We compared coefficient-based TBNNs against potential-based models and discussed and summarized techniques to obtain reference coefficient values from stress-strain data pairs.

In our cases studies involving six test datasets for three materials and two noise levels, the representations derived from the classical Finger-Rivlin-Ericksen formulation yield the best generalization performance. This was surprising to us, because we had initially believed that in particular, the

Figure 15: Reference and predicted coefficients around the undeformed configuration for uniaxial (\(\gamma<0\)) and biaxial (\(\gamma>0\)) extension for the Rivlin-Coeff and St.V-Coeff representations obtained from the Mooney-Rivlin model.
potential advantages of orthogonal bases, e.g. continuity of the coefficients and linear independence, would also translate into better extrapolations. We found that the generalization capabilities of the stress representations is largely dependent on the simplicity and lower complexity of the coefficient functions in the sense that the coefficient functions are smooth and monotonic, like low-order polynomials. This is the case when the stress generators (bases) already describe a majority of the stress-strain data complexities and hence the invariant-coefficient mappings are simpler. The introduction of physics-constrained extensions to the TBNN in particular convexity (potential) or monotonicity (coefficients) appears to be the most beneficial for accurate extrapolation and generalization. We also observed that the assumption of the existence of a potential helps the performance,

Figure 16: Reference and predicted coefficients around the undeformed configuration for uniaxial (\(\gamma<0\)) and biaxial (\(\gamma>0\)) extension for the Rivlin-Coeff and Convex-Pot representations obtained from the Mooney-Rivlin model.

Figure 17: Reference and predicted coefficients around the undeformed configuration for uniaxial (\(\gamma<0\)) and biaxial (\(\gamma>0\)) extension for the Mono-Coeff representation obtained from the Mooney-Rivlin model.
i.e. potential-based models generalize better than coefficient-based models for the same stress representation.

In future work the study will be extended to include anisotropic representations. We also aim to introduce monotonically increasing neural network formulations that eliminate the restrictiveness of current implementations, where a subset of the outputs is monotonically increasing with only a subset of the inputs, to potentially improve the performance of the monotonic Rivlin representation.

## Acknowledgments

REJ would like to thank Professor J.B. Estrada (University of Michigan) and NB would like to thank Professor K.T. Ramesh (Johns-Hopkins University) for independently pointing out Criscione's work on stress representations with an orthogonal basis.

JF and NB gratefully acknowledge support by the Air Force Office of Scientific Research under award number FA9550-22-1-0075.

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Figure 18: Parameterization in invariant space given by \(\gamma\)
Figure 19: Selected Mooney-Rivlin coefficient curves over parameter \(\gamma\). Dashed lines: extracted coefficients, solid lines: trained NN coefficients.
Figure 20: Selected Mooney-Rivlin stress-strain curves over parameter \(\gamma\). Dashed lines: data, solid lines: trained NN.