¹

Key words. Polyconvexity, Input Convex Neural Network, relaxation, isotropic damage model, parameter-dependent
AMS subject classifications. 49J45, 49J10, 74G65, 74B20, 68T07, 74A45

1 Introduction↩︎

Many relevant problems in solid mechanics seek to identify the minimising deformations \(u\colon \Omega \to \mathbb{R}^d\) of a body \(\Omega\) in spatial dimension \(d \in \{2,3\}\) with respect to the energy functional \[\label{eq:energyfunctional} I(u) = \int_{\Omega} W(\nabla u; \zeta) \,\mathrm{d}x,\tag{1}\] where the energy density \(W \colon \mathbb{R}^{d \times d} \times \mathbb{R}^{p} \to \mathbb{R}_{\infty} \mathrel{\vcenter{:}}= \mathbb{R}\cup \{\infty\}\), with \(p\in \mathbb{N}\), is defined over a suitable class of admissible functions. The vector \(\zeta \in \mathbb{R}^p\) represents spatially varying parameters, including material history, damage evolution or phase information. In many established engineering models, the density \(W\) is either inherently nonconvex or develops nonconvexity as its parameters evolve, as in damage or phase transition models. Such nonconvexity not only poses significant mathematical challenges, including the potential non-existence of minimisers, but also leads to serious problems in numerical simulation, such as mesh dependence and reduced robustness.

A common approach to overcome these issues is relaxation via the (semi-)convexification of the energy density, where the term semi-convexity covers the notion of polyconvexity, quasiconvexity and rank-one convexity. By replacing the nonconvex energy density \(W\) with its (semi-)convex envelope, a relaxed functional is obtained, which guarantees the existence of minimisers and allows for the identification of bounds of infimising sequences for the original nonconvex problem. On the computational side, these favourable properties translate into robust and mesh-insensitive numerical results.

As shown in [1]–[3], the notion of polyconvexity is well suited for nonlinear elasticity. In general, for a given function \(W\), it is not straightforward to verify whether \(W\) is polyconvex, which has led to extensive investigations of sufficient and necessary conditions in the literature, see e.g.[4]–[16]. However, for the relaxation of nonpolyconvex functions, explicit analytical representations of the polyconvex envelope \(W^{{\operatorname{pc}}}\) are typically not available, necessitating computational relaxation methods. Conventional approaches for relaxation include optimisation-based and computational geometry methods, such as algorithms approximating the rank-one convex envelope [17]–[23] or dedicated algorithms for computing the polyconvex envelope [24]–[27]. However, these methods are often subject to the curse of high-dimensionality, since they require the mesh-based discretisation of the (\(d \times d\))-dimensional space of deformation gradients. Polyconvexification algorithms further increase the dimensionality issue because they perform the convexification in the minors space of the deformation gradient. For isotropic functions, this dimension can be reduced by characterising the polyconvexity via signed singular values, as in [15] (see [28] for a corresponding polyconvexification algorithm). However, despite their efficiency, these algorithms remain resource and time consuming, especially in spatial dimension \(d=3\). As a result, conventional relaxation algorithms are often impractical for engineering applications that require iterative computations and multi-query evaluations of the polyconvex envelope for varying parameter values.

For example, frameworks such as the isotropic pseudo-time incremental damage model introduced in [29] involve a (\(d\times d + 1\))-parameter-dependent family of convexification problems, consisting of the deformation gradient and the scalar internal variable from the previous time step that captures the material history (see 2). Such parameter dependence introduces additional dimensionality as the material parameter \(\zeta \in \mathbb{R}^{p}\) (for the damage example it holds \(p = d^2 + 1\)) can vary at any material point in \(\Omega\). This requires the polyconvex envelope to be computed from scratch for each parameter configuration, resulting in a (\(d\times d+ p\))-dimensional problem. Overall, the accurate approximation of the polyconvex envelope for parameter-dependent energy densities is of great importance for the practical application of computational relaxation techniques.

To overcome the high computational cost of conventional relaxation schemes and make them practically feasible, we propose a neural network-based approach that compresses the parameter-dependent polyconvex envelope into the parameters of a small-scale artificial neural network that incorporates intrinsic properties of the polyconvexification problem into its architecture. Our method exploits a design that encodes polyconvexity by combining a feature extraction layer which computes the minors function based on a signed singular value characterisation of isotropic functions, with Input Convex Neural Networks (ICNNs) [30]. In addition, symmetry properties and envelope inequality are systematically enforced during training through penalty terms in the loss functional, ensuring that the network satisfies the necessary physical and mathematical constraints. This tailored approach not only speeds up the approximation of the polyconvex envelope, but also increases the reliability of the approximations in practical engineering applications.

Recent advances in neural networks, particularly through the development of ICNNs [30], have enabled architectures that enforce convexity in energy density functions, a key requirement for accurate material modelling. For example, [31], [32] use ANNs to learn constitutive laws from stress–strain data, while preserving fundamental principles of solid mechanics through ICNNs that ensure material stability. It should be emphasised that in these two works the ANNs are convex with respect to the strain tensor \(E=\frac{1}{2}(F^TF - \mathbb{I})\) instead of the deformation gradient \(F\). These models are thus not guaranteed to be polyconvex and are therefore not elliptic for all states. In addition, [33] and [34] propose machine-learning-based constitutive models for finite deformation and electro-mechanically coupled behaviour respectively, using ICNNs to enforce hyperelasticity, anisotropy and polyconvexity. In [35], a neural network-based hyperelastic constitutive model is introduced that inherently satisfies standard constitutive conditions using ICNNs, while [36] designs multi-scale heterogeneous structures with spatially varying microstructures by incorporating physical principles such as polyconvexity, objectivity and thermodynamic consistency.

Without using ICNNs, some works have also explored ANNs for modelling the behaviour of mechanically sound materials [37]–[42]. More closely to our approach, other works focus on (Partially) Input Convex Neural Network formulations in the invariants of the right Cauchy–Green tensor \(C\), sufficient but not necessary for polyconvexity, see e.g. [43]–[45]. Recently, [46] introduced convex signed singular value neural networks (CSSV-NN) for isotropic hyperelastic energy formulations and in [47] attention was drawn to the representation of incompressible materials. Unlike [46], which aims to learn polyconvex hulls directly from values of potentially nonpolyconvex functions, our approach learns and compresses polyconvex envelopes from reference solutions obtained either analytically or via established numerical algorithms [28]. This ensures polyconvex data as input and provides stability in the training process, as well as accuracy and control over the learning target within the prescribed domain. The approach aims to construct a compressed representation of the polyconvex envelope on the entire signed singular value space, facilitating its efficient evaluation in complex engineering loading scenarios. In addition, the method extends to the more general case of parameter-dependent classes of reference envelopes: the trained network can thus be applied to entire families of polyconvexification tasks, for example within damage simulations, rather than being limited to a single envelope.

The remaining parts of this paper are structured as follows. In 2, we introduce a guiding engineering application by presenting a parameter-dependent family of energy densities arising from an isotropic damage problem. In 3, we recall the signed singular value formulation of the polyconvex envelope for isotropic functions, while 4 introduces the properties-preserving neural networks that enforce polyconvexity, symmetries and the envelope inequality. Finally, [sec:numericalexperimentsmath] [sec:numericalexperimentsisodamage] provide numerical examples in two and three spatial dimensions including both mathematical benchmarks and the engineering isotropic damage formulation, rewritten in the signed singular value framework and parameter-reduced via a splitting approach.

2 Isotropic Damage Problem: A Guiding Example↩︎

The numerical simulation of the isotropic damage model of [29] is a prime representative of a problem where parameter-dependent polyconvexification is encountered in an engineering application and serves as a guiding example for our work.

Let \(F\) denote the deformation gradient, i.e. \(F = \nabla u\) for a sufficiently smooth deformation \(u\colon \Omega \to \mathbb{R}^d\) of a body \(\Omega\) in spatial dimension \(d \in \{2,3\}\) as in 1 . In what follows, we consider two different effective hyperelastic strain energy density functions \(\psi^{0}\) in the undamaged state: the Saint Venant–Kirchhoff (STVK) model \[\label{eq:STVK} \psi^{0}_{\rm STVK}(F) = \frac{\mu}{4} \, \lvert F^T F - \mathbb{I}\rvert^2 + \frac{\lambda}{8} \, \left(\lvert F \rvert^2 - d\right)^2,\tag{2}\] representing the special case of a linear stress–strain relation in a Lagrangian formulation, and a compressible neo-Hookean (NH) model \[\label{eq:NH} \psi^0_{\rm NH}(F) = \frac{\mu}{2} (\operatorname{tr}(F^T F) - d) - \mu \ln(\det F) + \frac{\lambda}{2} \ln(\det F)^2 ,\tag{3}\] showing a nonlinear response, which also incorporates the determinant constraint for the penalisation of negative determinants of the deformation gradient. Here, \(\lambda\) and \(\mu\) denote the Lamé constants.

The pseudo-time incremental energy density function capturing isotropic damage for the time step \(k\to k+1\) is formulated in [29] as \[\label{eq:Wdamage} \begin{align} W(F_{k+1}; F_{k}, \alpha_{k}) & = \psi(F_{k+1},p(F_{k + 1}; \alpha_{k})) - \psi(F_k,\alpha_k) \\ &\qquad + p(F_{k + 1}; \alpha_{k})\, D(p(F_{k + 1}; \alpha_{k})) - \alpha_k \, D(\alpha_k) - \overline{D}(p(F_{k + 1}; \alpha_{k})) + \overline{D}(\alpha_k). \end{align}\tag{4}\] In 4 , the scalar \(\alpha_k \in \mathbb{R}\) and the matrix \(F_k \in \mathbb{R}^{d\times d}\) denote the internal variable and the deformation gradient from the previous pseudo-time step \(k\), respectively, and serve as parameters for the current step \(k+1\); they are considered fixed and known. The strain energy density function \(\psi\) is given by \[\label{eq:psi1minusD} \psi(F, \alpha) = (1 - D(\alpha)) \, \psi^{0}(F),\tag{5}\] where \(\psi^{0}\) is an isotropic energy density such as the STVK 2 or NH 3 examples above. Through this formulation, the characteristics of the undamaged energy density, such as isotropy or determinant constraints, enter the model. The function \(D\colon \mathbb{R}\to [0,1)\) is a non-decreasing exponential damage function, as considered in [48], which maps the internal variable \(\alpha\) to the interval \([0,1)\), with \(0\) corresponding to the undamaged state and values close to \(1\) indicating full damage. Specifically, \[\label{eq:damagefunction} D(\alpha) = d_\infty \left(1 - \exp\left(- \frac{\alpha}{d_0}\right)\right),\tag{6}\] where \(d_\infty \in (0,1)\) is the asymptotic limit (the maximum possible damage) and \(d_0 \in \mathbb{R}^+\) is the damage saturation parameter. The reciprocal of \(d_0\) reflects the rate at which the asymptotic limit \(d_\infty\) is reached. The function \[\overline{D}(\alpha) = d_\infty \left(\alpha + d_0 \, \exp\left(-\frac{\alpha}{d_0}\right)\right)\] is the antiderivative of \(D\). The evolution of the internal variable \(\alpha_{k+1}\), and hence the damage \(D(\alpha_{k+1})\), is determined by a path function \(p\), defined as \[\label{eq:pathfunctionW} \alpha_{k+1} = p(F_{k+1}; \alpha_{k}) = \begin{cases} \psi^{0}(F_{k+1}) & \text{if } \, \psi^{0}(F_{k+1}) > \alpha_{k}, \\ \alpha_{k} & \text{else}. \end{cases}\tag{7}\] The function \(W\) in 4 is nonconvex because the damage evolution process is governed by the path function \(p\). Theoretically, nonconvexity prevents the existence of minimisers. Numerically, it poses challenges such as stability problems, nonconvergence and oscillatory behaviour in the simulation of boundary value problems; see [20]–[22], [29], where these issues have been addressed by rank-one convexification, for further discussion of these issues.

3 Polyconvexification of Isotropic Functions↩︎

In this section, we present a polyconvexification approach based on a reformulation in terms of the signed singular values of the deformation gradient, which underpins our neural network design. For ease of notation, we omit the explicit parameter dependence and consider functions \(W\colon\mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\). The parameter-dependent formulation is obtained by direct analogy, by including the parameter vector \(\zeta\) and considering the function \(W(\,\cdot\,; \zeta)\).

Assuming that \(W\) is isotropic, we use its characterisation by signed singular values to reduce the dimension of the problem. In this formulation, polyconvexification of an isotropic function reduces to convexification of a function defined on the manifold determined by the minors of the signed singular values. This reformulation is advantageous because it reduces the domain from \(\mathbb{R}^{d \times d}\) to \(\mathbb{R}^d\), thereby reducing the computational cost of a grid-based representation of the polyconvex envelope. This reduced approach has been successfully applied in an efficient conventional polyconvexification algorithm [28] (see also 4.6). In the context of this work, it serves as the basis for learning and predicting the polyconvex envelope via properties-preserving artificial neural networks.

Let \(d \in \{2,3\}\) and let \(W \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\) be a function which maps \((d \times d)\)-matrices to real scalars or infinity. The notion of polyconvexity relies on the minors of the matrices \(F \in \mathbb{R}^{d \times d}\). Given the determinant \(\det(F)\) and the adjugate \(\operatorname{adj}(F)\) of \(F\), let \[\label{eq:M40F41} \mathcal{M}(F)= \begin{cases} (F, \det(F)) & \text{if } \; d=2,\\ (F, \operatorname{adj}(F), \det(F)) & \text{if } \;d=3 \end{cases}\tag{8}\] denote the minors of \(F\), i.e. \(\mathcal{M}(F)\) is a vector of dimension \(K_d = 5\) if \(d =2\) and \(K_d = 19\) if \(d=3\). A function \(V \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\) is said to be polyconvex if there exists a convex function \(G\colon \mathbb{R}^{K_d} \to \mathbb{R}_{\infty}\) such that for all \(F \in \mathbb{R}^{d \times d}\), \[V(F) = G(\mathcal{M}(F)).\] The polyconvex envelope \(W^{{\operatorname{pc}}}(F) \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\) of \(W\), defined by the pointwise supremum \[\label{eq:defWpc} W^{{\operatorname{pc}}}(F) = \sup \{V(F) \;\lvert \;V \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty} \;\text{polyconvex}, V \leq W \},\tag{9}\] is the largest polyconvex function below \(W\).

We will restrict ourselves to the class of isotropic densities and call \(W\) isotropic if and only if \[W(F) = W(R_1 F R_2)\] for all \(F\in \mathbb{R}^{d \times d}\) and all \(R_1, R_2 \in \mathcal{SO}(d)\) and we say that \(W\) is \(\mathcal{SO}(d)\times \mathcal{SO}(d)\)-invariant in this case, where \(\mathcal{SO}(d)\) denotes the special orthogonal group of \((d \times d)\)-matrices. Note that we define isotropy in the sense as introduced in [1], the left \(\mathcal{SO}(d)\)-invariance includes objectivity, while the right invariance reflects full material symmetry. Following [15], isotropic functions can be characterised by the signed singular values of their arguments. Let \(0\leq \sigma_1,\ldots,\sigma_d\in\mathbb{R}\) denote the singular values of a matrix \(F\in\mathbb{R}^{d\times d}\). Then \(\nu_1,\ldots, \nu_d \in \mathbb{R}\) are called signed singular values of \(F\), if they have the same absolute values as the singular values of \(F\), i.e. up to permutation it holds \(|\nu_i| = \sigma_i\), and their signs satisfy \(\operatorname{sign}(\nu_1 \cdot \ldots \cdot \nu_d) = \operatorname{sign}(\det(F))\). In this definition, the signed singular values are only unique up to permutations in \[\Pi_d= \left\{P \operatorname{diag}(\varepsilon) \in \mathcal{O}(d) \mid P \in \operatorname{Perm}(d), \varepsilon \in \{-1,1\}^d, \varepsilon_1\cdot\ldots\cdot\varepsilon_d = 1 \right\},\] where \(\operatorname{diag}\) refers to the diagonal matrix with entries given by the vector of its argument and \(\operatorname{Perm}(d) \subset \{0,1\}^{d \times d}\) denotes the set of permutation matrices. In what follows, we denote by \(\nu \colon \mathbb{R}^{d \times d} \to \mathbb{R}^{d}\) the signed singular value mapping assumed to be well-defined, e.g. by fixing the order of the entries according to the magnitude of the absolute values and assuming positivity of the entries up to only one entry, i.e. assigning the sign of the determinant to a single fixed entry.

It is possible to identify the set of isotropic functions \(W\colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\) with the set of \(\Pi_d\)-invariant functions \(\Phi\colon\mathbb{R}^d \to \mathbb{R}_{\infty}\), i.e. \({\Phi(\hat{\nu}) = \Phi(S \hat{\nu})}\) for all \(\hat{\nu}\in \mathbb{R}^{d}\) and all \(S \in \Pi_d\). The identifications are given by \[\label{eq:W61PhiPhi61W} W(F) = \Phi(\nu(F)) \qquad \text{and} \qquad \Phi(\hat{\nu}) = W(\operatorname{diag}(\hat{\nu}))\tag{10}\] for all \(F\in \mathbb{R}^{d\times d}\) and for all vectors \(\hat{\nu} \in \mathbb{R}^{d}\).

As shown in [15], the polyconvexity of isotropic functions can be characterised directly using a lower dimensional mapping \(\Phi\). For this purpose, in analogy to the minors \(\mathcal{M}\) in 8 , we define \(k_d \mathrel{\vcenter{:}}= 2^d - 1\), and introduce the mapping \(m\colon \mathbb{R}^{d} \to \mathbb{R}^{k_d}\) by \[m(\hat{\nu}) = \begin{cases} ({\nu}_1, {\nu}_2, {\nu}_1 \, {\nu}_2) & \text{if } \;d = 2, \\ ({\nu}_1, {\nu}_2, {\nu}_3, {\nu}_2 \, {\nu}_3, {\nu}_3 \, {\nu}_1, {\nu}_1 \, {\nu}_2, {\nu}_1 \, {\nu}_2 \, {\nu}_3) & \text{if } \;d = 3. \end{cases}\] We call \(m(\hat{\nu})\) the vector of minors of \(\hat{\nu} = (\nu_1,\ldots, \nu_d)\in\mathbb{R}^d\); in the literature, this vector is also called the elementary polynomials. A \(\Pi_d\)-invariant function \(\Psi\) is called (signed singular value) polyconvex if there exists a convex function \(g\colon \mathbb{R}^{k_d} \to \mathbb{R}_{\infty}\) such that \(\Psi = g\circ m\), see [15]. The (signed singular value) polyconvex envelope of \(\Phi\) is then defined by \[\label{eq:defPhipc} \Phi^{{\operatorname{pc}}}(\hat{\nu}) = \sup\{\Psi(\hat{\nu}) \mid \Psi \colon \mathbb{R}^{d} \to \mathbb{R}_{\infty} \text{ polyconvex }, \Psi \leq \Phi\}.\tag{11}\] According to [28], the polyconvex envelope \(W^{{\operatorname{pc}}}\) of the original density can be identified with \(\Phi^{{\operatorname{pc}}}\), similar to 10 , in the following sense: \[\label{eq:Wpc61PhipcPhipc61Wpc} W^{{\operatorname{pc}}}(F) = \Phi^{{\operatorname{pc}}}(\nu(F)) \qquad \text{and} \qquad \Phi^{{\operatorname{pc}}}(\hat{\nu}) = W^{{\operatorname{pc}}}(\operatorname{diag}(\hat{\nu})).\tag{12}\] So we can limit ourselves to the approximation of \(\Phi^{{\operatorname{pc}}}\). Thanks to [28], it can be obtained by \[\label{eq:Phipc61hc} \Phi^{{\operatorname{pc}}}(\hat{\nu}) = h^{{\operatorname{c}}} (m(\hat{\nu})),\tag{13}\] where \(h^{{\operatorname{c}}}\) is the convex envelope of the function \[h\colon \mathbb{R}^{k_d} \to \mathbb{R}_\infty, \qquad x \mapsto \begin{cases} \Phi(\hat{\nu}) &\text{if } \;x = m(\hat{\nu}), \\ \infty &\text{else}. \end{cases}\] The relation 13 transforms the polyconvexification problem of \(\Phi\) into a convexification problem of the function \(h\) in the lifted signed singular value space. Compared to the original polyconvexification problem defined on \(d \times d\) matrices, this formulation reduces the dimension of the manifold to be convexified from \(d \times d\) to \(d\). In particular, as shown in [28], the dimension of the ambient space is reduced from \(19\) and \(5\) to \(7\) and \(3\) for \(d=3\) and \(d=2\), respectively, allowing an efficient numerical treatment.

4 Properties-Preserving Neural Networks↩︎

Given the notion of polyconvexity for isotropic functions introduced in the previous section, we want to approximate the parameter-dependent polyconvex envelope \(\Phi^{{\operatorname{pc}}}\colon \mathbb{R}^{d}\times \mathbb{R}^{p} \to \mathbb{R}\) using a neural network denoted by \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\colon \mathbb{R}^{d}\times \mathbb{R}^{p} \to \mathbb{R}\). This is achieved by approximating the convex envelope \(h^{{\operatorname{c}}}\) of the function \(h\colon \mathbb{R}^{k_d}\times \mathbb{R}^{p} \to \mathbb{R}_{\infty}\) defined in 13 . For clarity, we restrict ourselves to the case where both \(\Phi^{{\operatorname{pc}}}\) and \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) are finite and do not attain the value infinity. In some expressions, we simplify the notation by writing \(\hat{m} \mathrel{\vcenter{:}}= m(\hat{\nu})\) for the minors of the signed singular values vector \(\hat{\nu}\).

4.1 Neural Network Concepts↩︎

An artificial neural network is a mathematical model defined as a composition of functions, each of which is typically formed by further compositions of functions [49]. This model can be conveniently represented as a network structure. In what follows, we present the basic notions of layers, weights, inputs and outputs of a general neural network.

Consider a neural network consisting of \(k\) layers. We denote the index of layers by the superscript \(i \in \{0,\ldots,k\}\), where \(i = 0\) denotes the input layer, \(i = 1\) the first hidden layer and \(i = k\) the output layer. The total number of hidden layers is thus \(k-1\). The number of neurons in the \(i\)-th layer is denoted by \(d_{i}\); in particular, \(d_{0}\) is the number of inputs and \(d_{k}\) is the number of outputs. Let \(W_{i} \in \mathbb{R}^{d_{i+1} \times d_{i}}\) denote the weight matrix associated with the connection from the \(i\)-th to the \((i+1)\)-th layer, \(b_{i} \in \mathbb{R}^{d_{i + 1}}\) the corresponding bias vector and \(g_i\) the activation function, acting component-wisely. For \(i = 0, \dots, k-1\), the network recurrence on the layer outputs \(z_i \in \mathbb{R}^{d_{i}}\) can be written in matrix–vector form as \[z_{i + 1} = g_{i}\!\left(W_{i} \, z_{i} + b_{i} \right),\] where \(z_{0}\) denotes the input and \(z_{k}\) the overall output of the network.

Building on the characterisation of the (signed singular value) polyconvex envelope \(\Phi^{{\operatorname{pc}}}\) in 3, we aim to construct a neural network approximation \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) that reliably preserves its physical properties. The architecture must ensure polyconvexity, i.e. convexity with respect to the minors of the signed singular values, and satisfy the envelope inequality \(\Phi^{{\operatorname{pc}}} \leq \Phi\) (pointwisely) as in 11 . In addition, the network must respect the \(\Pi_{d}\)- invariance of \(\Phi\) representing symmetries. The following sections outline strategies for enforcing these three key properties into the design of neural networks.

4.2 Enforcing Polyconvexity↩︎

To enforce the convexity in the minors \(m(\hat{\nu})\) of the input vector \(\hat{\nu}\), a particular class of neural networks is employed: Input Convex Neural Networks (ICNN) introduced in [30]. In what follows, two variants of ICNN are considered: the Fully Input Convex Neural Networks (FICNN) and the Partially Input Convex Neural Networks (PICNN).

4.2.1 Fully Input Convex Neural Networks (FICNN)↩︎

The model illustrated in 1 defines a neural network over the input \(\hat{m}\) using the recurrence for \(i=0, \dotsc, k-1\), \[z_{i+1} = g_i(W_i^{(z)} z_i + W_i^{(\hat{m})} \hat{m} + b_i),\] where \(z_i\) denotes the layer activations (with \(z_0 = 0\) , \(W_0^{(z)} \equiv 0\)) and \(g_i\) are nonlinear activation functions. The overall network evaluation reads \[\Phi^{{\operatorname{pc}}}_{{\rm pred}}(\hat{\nu}) \mathrel{\vcenter{:}}= h^{\rm c}_{\rm pred}(\hat{m};\theta) = z_k,\] where \(h^{\rm c}_{\rm pred}\) is the neural network approximation of the function \(h^{\rm c}\) defined in 13 and depends on \(\theta = \{ W_{0:k-1}^{(\hat{m})}, W_{1:k-1}^{(z)}, b_{0:k-1} \}\), the parameter vector collecting the weights and biases.

The convexity of \(h^{{\operatorname{c}}}_{\rm pred}\) is addressed first. In order for \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) to be polyconvex as defined in 3, the \(\Pi_d\)-invariance is also required, which will be addressed in 4.4.

The proof follows from the fact that non-negative linear combinations of convex functions are also convex, and that the composition of a convex and convex non-decreasing function is also convex. The constraint that the \(g_i\) are convex and non-decreasing is not particularly restrictive, since common nonlinear activation functions such as the ReLU or SoftPlus activation functions satisfy this constraint. A review of convex activation functions is presented in [50].

4.2.2 Partially Input Convex Neural Networks (PICNN)↩︎

The previous model, i.e. FICNN, provides convexity over the entire input of the neural network, which may in fact be a restriction on the allowable class of models. Furthermore, this full joint convexity is unnecessary in the setting where the convexity is required only over some inputs. This is why we present the Partially Input Convex Neural Networks (PICNN), as introduced in [30], which are convex over only some inputs of the network. The \(k\)-layer PICNN architecture (2) is defined by the recurrence for \(i=0, \dotsc, k-1\), \[\begin{align} &u_{i+1} = \tilde{g}_i \left (\tilde{W}_i u_i + \tilde{b}_i \right), \\ &z_{i+1} = g_i\left( W_i^{(z)}\Big(z_i \circ \max([W_i^{(zu)}u_i+b_i^{(z)}],0)\Big) + W_i^{(\hat{m})} \Big(\hat{m} \circ [W_i^{(\hat{m}u)}u_i+b_i^{(\hat{m})}]\Big) + W_i^{(u)}u_i + b_i \right), \end{align}\] with \(u_0 =\zeta\), \(z_0 =0\) and \(W_0^{(z)}=0\) and \[\Phi^{\operatorname{pc}}_{\rm pred}(\hat{\nu}; \zeta) \mathrel{\vcenter{:}}= h_{\rm pred}^{\rm c}(\hat{m};\zeta;\theta) = z_k\] denotes the overall evaluation of the network. Moreover, \(h_{\rm pred}^{{\operatorname{c}}}\) is the neural network approximation of the function \(h^{{\operatorname{c}}}\) defined in 13 and depends on \(\theta\), the parameter vector collecting the weights and biases. The \(u_i \in \mathbb{R}^{n_i}\) and \(z_i \in \mathbb{R}^{m_i}\) denote the hidden units for the \(\zeta\)-path and \(\hat{m}\)-path, and \(\circ\) denotes the Hadamard product, i.e. the element-wise product between two vectors. The functions \(g_i\) and \(\tilde{g}_i\) are activation functions.

Note that the schemes shown in 1 2 already reflect the convexity with respect to the minors \(\hat{m}\), i.e. the formulation to ensure polyconvexity. However, the change from the original input \((\hat{\nu}, \zeta)\) of the function \(\Phi^{\operatorname{pc}}\) to the input \((m(\hat{\nu}), \zeta)\) of the function \(h^{{\operatorname{c}}}\), i.e. the transitioning from the signed singular values to the minors lifted signed singular values, can be interpreted as a hard-coded feature extracting layer in the overall network architecture.

In the present work, FICNN is employed when the function \(\Phi^{{\operatorname{pc}}}\) to be predicted does not depend on parameters but only on the signed singular values vector \(\hat{\nu}\), and PICNN is employed when the function to be predicted depends on the signed singular values vector \(\hat{\nu}\) along with additional parameters \(\zeta\). Indeed, the use of PICNN allows ensuring convexity with respect to the minors \(m(\hat{\nu})\) of the signed singular values vector \(\hat{\nu}\) while relaxing the constraints on the \(\zeta\)-path, allowing a better representation.

4.3 Enforcing the Envelope Inequality↩︎

By definition 11 , the polyconvex envelope \(\Phi^{\operatorname{pc}}\) of the function \(\Phi\) is the largest polyconvex function below \(\Phi\). Consequently, the goal is to ensure that the polyconvex envelopes predicted by the neural network are below the function \(\Phi\), i.e. \(\Phi^{{\operatorname{pc}}}_{{\rm pred}} \leq \Phi\) which we refer to as envelope inequality. In the present work, the envelope inequality is imposed in a weak sense by penalisation, avoiding the addition of constraints directly in the architecture. To achieve this, a custom loss function is designed, which penalises predictions of the neural network located above the function \(\Phi\) during the training process. By minimising this loss function, the network ultimately achieves predictions below \(\Phi\). The loss function is thus written as \[\label{eq:loss1} \mathcal{L}= \mathcal{L}_{{\rm mse}} + \lambda_{{\rm ineq}} \, \mathcal{L}_{{\rm ineq}}\tag{14}\] with \[\mathcal{L}_{{\rm mse}} = \frac{1}{N} \sum_{i=1}^N (\Phi^{{\operatorname{pc}}}_{i} - \Phi^{{\operatorname{pc}}}_{{\rm pred},i})^2 , \label{eq:Lmse}\tag{15}\] the classical mean square error, and \[\label{eq:lossineq} \mathcal{L}_{{\rm ineq}} = \frac{1}{N} \sum_{i=1}^N \max \{\Phi^{{\operatorname{pc}}}_{{\rm pred},i} - \Phi_i, 0 \}^2 ,\tag{16}\] where \(N\) is the size of the learning data, \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) is the neural network prediction of the polyconvex envelope, and \(\Phi^{{\operatorname{pc}}}\) is the target value to be approximated. The subscript \(i\) denotes the evaluation of these functions at the \(i\)th learning input data point given as a tuple of the form \((m(\hat{\nu}), \zeta)\). The parameter \(\lambda_{{\rm ineq}} \geq 0\) is the so-called penalty parameter, which determines how the two terms are weighted against each other and determines how strongly the network should be penalised if the prediction lies above \(\Phi\).

4.4 Enforcing the \(\Pi_{d}\)-Symmetry↩︎

The function \(\Phi\) from 10 is subject to the invariance of \(\Pi_{d}\), and the same holds for \(\Phi^{{\operatorname{pc}}}\) from 12 , that is it holds \[\Phi^{{\operatorname{pc}}}(\pi(\hat{\nu})) = \Phi^{{\operatorname{pc}}}(\hat{\nu})\] for all \(\pi \in \Pi_{d}\). Consequently, this symmetry should be preserved by a neural network approximation to ensure the \(\Pi_d\)-invariance as required for the notion of (signed singular value) polyconvexity, as defined in 3.

As for the envelope inequality in 4.3, the symmetry of the prediction can be imposed in a weak sense by penalisation. To achieve this, a custom loss function which penalises more severely if the prediction of the neural network is not symmetric during the training process is designed. By minimising this loss function, the network should ultimately achieve symmetric predictions. Building upon the previously defined loss function 14 (keeping the same notations), already facilitating the envelope inequality, we add a term for the penalisation of the non-symmetry, leading to the loss function defined by \[\label{eq:loss} \mathcal{L}= \mathcal{L}_{{\rm mse}} + \lambda_{{\rm ineq}} \, \mathcal{L}_{{\rm ineq}} + \lambda_{{\rm sym}} \, \mathcal{L}_{{\rm sym}}\tag{17}\] with \[\label{eq:losssym} \mathcal{L}_{{\rm sym}} = \frac{1}{\lvert \Pi_d \rvert} \sum_{\pi \in \Pi_d} \frac{1}{N} \sum_{i=1}^{N} \left(\Phi^{{\operatorname{pc}}}_{{\rm pred},i}(\hat{\nu}) -\Phi^{{\operatorname{pc}}}_{{\rm pred},i}(\pi(\hat{\nu}))\right)^2,\tag{18}\] where the normalisation constant is motivated by the fact that the term corresponding to \(\pi = \operatorname{id}_{{}}\) does not contribute to the overall loss. The parameter \(\lambda_{{\rm sym}} \geq 0\) is the penalty parameter, which determines how strongly the network should be penalised if the network is not symmetric with respect to the inputs. In addition, symmetry is further fostered by data augmentation. We incorporate symmetric data in the dataset, i.e. for a given point \(\hat{\nu}\), the points \(\pi(\hat{\nu})\) for \(\pi \in \Pi_d\) are also included in the dataset.

4.5 Implementation and Hyperparameters↩︎

The neural network architectures are implemented using PyTorch, and all the trainings are performed on an Intel^® Core^TM i9-11900K machine, using only one core at a base frequency of 3.50 GHz.

The input data of the neural networks consist of the minors \(m(\hat{\nu}) \in \mathbb{R}^{k_d}\) of the signed singular values \(\hat{\nu}\) along with additional parameters \(\zeta \in \mathbb{R}^{p}\) depending on the definition of the function \(\Phi\), i.e. the input data consists of vectors \((m(\hat{\nu}), \zeta) \in \mathbb{R}^{k_d + p}\). The target values correspond to the evaluation of the polyconvex envelope \(\Phi^{{\operatorname{pc}}}\) at the points \(\hat{\nu}\) and the parameters \(\zeta\), expressed as \(\Phi^{{\operatorname{pc}}}(\hat{\nu}; \zeta)\). In summary, the learning data (gathering both training and validation data) consists of tuples of the form \((m(\hat{\nu}), \zeta, \Phi^{{\operatorname{pc}}}(\hat{\nu}; \zeta), \Phi(\hat{\nu}; \zeta))\) containing the input features, the target value \(\Phi^{{\operatorname{pc}}}(\hat{\nu}; \zeta)\), and the value of the function \(\Phi(\hat{\nu}; \zeta)\), which is necessary for the evaluation of the loss function 17 .

In all the numerical examples, we design the neural network architectures to be as small as possible while maintaining high predictive accuracy. Indeed, preliminary numerical experiments have shown that increasing the number of layers or neurons of the neural networks presented below does not necessarily improve prediction performance; on the contrary, it may even degrade them. Also, we aim to use the smallest dataset size that ensures accurate predictions, while making training computationally efficient and feasible on a standard laptop, which has been achieved thanks to several investigations. We employ the ReLU activation function in all hidden layers, as it satisfies the convexity condition required by ICNN. The output layer consists of a single neuron with a linear activation function. Numerical experiments indicate that selecting ReLU as the convex activation function enhances learning efficiency and improves predictive accuracy compared to the SoftPlus function, another commonly employed activation function in ICNNs.

The weights \(W^{(z)}\) which belong to the convex part of the network are initialised from a normal distribution \(\mathcal{N}(0.1,0.1)\) with \(0.1\) mean and \(0.1\) standard deviation, and then projected onto \(\mathbb{R}_+\). Also, after each training step, the weights \(W^{(z)}\) are projected onto the \(\mathbb{R}_+\) halfspace to ensure their positiveness. This projection is essential to ensure the positivity of the weights. The naive application of ReLU for weight clipping can result in zero entries in the weight matrices and would waste approximation potential. For this projection, a shifted version of ReLu, i.e. \(x \to \max(x, 0) + \varepsilon\) with the small offset \(\varepsilon = 10^{-6}\) is employed. This initialisation and projection strategy has shown better performance compared to a uniform distribution and projections using exponential or SoftPlus functions. The other weights are initialised from a uniform distribution \(\mathcal{U}(- 1/\sqrt{n}, 1/\sqrt{n})\) with \(n\) the input size and all the biases are initialised according to \(\mathcal{N}(0, 0.1)\).

We employ the Adam optimiser with a learning rate of \(\eta = 0.001\), a batch size of \(128\), with the data shuffled at each epoch to improve model generalisation, and a patience of \(5\), i.e. the training process is stopped if no improvement of the validation loss is achieved for \(5\) successive epochs. Although more complex training strategies, e.g. with the use of learning rate decay, could have been implemented, numerical experiments demonstrate that such an approach seems unnecessary in this context. The training and validation losses are computed with the loss function 17 . The hyperparameters \(\lambda_{{\rm ineq}}\) and \(\lambda_{{\rm sym}}\) involved in \(\mathcal{L}\) have been determined by performing several preliminary tests with different parameter values for \(\lambda_{{\rm ineq}}\) and \(\lambda_{{\rm sym}}\). In each numerical example presented below, the given values of the \((\lambda_{{\rm ineq}}, \lambda_{{\rm sym}})\) achieved favourable training performance and predictions.

To obtain more reliable numerical results, for each example, we conduct ten or twenty independent runs of training, using the same training and validation data. The variations arise solely from the initialisation of model weights and the stochastic nature of data loading during training, since a shuffle is used. The numerical results presented are obtained by averaging the outputs of the ten network realisations. To assess the accuracy of the neural network predictions \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) in comparison to the ground truth \(\Phi^{{\operatorname{pc}}}\), obtained either analytically or by numerical algorithms, we compute the following error metrics, namely the mean error \[{\rm Mean \;err} = \frac{1}{N}\sum_{i=1}^N \lvert\Phi^{{\operatorname{pc}}}_{{\rm pred},i}-\Phi^{{\operatorname{pc}}}_{i}\rvert ,\] the relative quadratic error \[{\rm Rel\;quad \;err} = \sqrt{\frac{\sum_{i=1}^N \lvert \Phi^{{\operatorname{pc}}}_{{\rm pred},i}-\Phi^{{\operatorname{pc}}}_{i}\rvert^2}{\sum_{i=1}^N \lvert\Phi^{{\operatorname{pc}}}_{i}\rvert^2}} ,\] and the relative maximum error \[{\rm Rel\;max\;err} = \frac{\max_{i=1,\ldots, N} \lvert \Phi^{{\operatorname{pc}}}_{{\rm pred},i}-\Phi^{{\operatorname{pc}}}_{i} \rvert}{\max_{i=1,\ldots, N} \lvert\Phi^{{\operatorname{pc}}}_{i}\rvert},\] where \(N\) denotes the number of evaluation points.

4.6 Data Generation↩︎

For the physically relevant benchmark examples considered in the next sections, no closed-form analytical representations of the polyconvex envelopes are available. This is why, for the computation of the target values the polyconvex envelopes are approximated numerically utilising the algorithm for isotropic functions based on a linear programming approach presented in [28]. For this approximation procedure, the pointwise characterisation of the polyconvex envelope \(\Phi^{{\operatorname{pc}}}\) at \(\hat{\nu} \in \mathbb{R}^{d}\) by the optimisation problem \[\label{eq:poly-opt-prob-iso} \Phi^{{\operatorname{pc}}}(\hat{\nu}) = \inf \left\{\sum_{i = 1}^{k_d + 1} \xi_i \, \Phi(\nu_i)\,\biggl\vert\, \xi_{i} \in [0, 1],\, \nu_i \in \mathbb{R}^{d},\, \sum_{i = 1}^{k_d + 1} \xi_{i} = 1,\, \sum_{i = 1}^{k_d + 1} \xi_i \, m(\nu_i) = m(\hat{\nu})\right\}\tag{19}\] is employed. For an algorithmic approximation of this problem, consider a discretisation of the signed singular value space by a point cloud, denoted by \(\Sigma_{\delta} = \{\nu_1, \ldots, \nu_{N_\delta}\} \subset \mathbb{R}^{d}\), a possible choice is a structured lattice on \([-r, r]^{d}\) with lattice size \(\delta\) and discretisation radius \(r\). The lifting of this point cloud \(m(\Sigma_\delta)\) can be employed to turn the nonlinear optimisation problem 19 into the following linear program \[\label{eq:linearprogram} \Phi^{{\operatorname{pc}}}_{\delta} (\hat{\nu}) = \min \left\{\sum_{i = 1}^{N_\delta} \xi_{i} \, \Phi(\nu_i) \;\biggl\vert\; \xi_{i} \geq 0,\, \sum_{i = 1}^{N_\delta} \xi_{i} = 1,\, \sum_{i = 1}^{N_\delta} \xi_{i}\,m(\nu_{i}) = m(\hat{\nu}) \right\},\tag{20}\] which can be efficiently solved numerically using standard algorithms for linear programming. The overall procedure is referred to as signed singular value polyconvexification by linear programming (SVPC LP). In exact arithmetic, there exists a minimiser \(\xi \in \mathbb{R}^{N_\delta}\) of 20 with at most \(k_d+1\) non-zero entries reflecting the convex coefficients of the supporting points of the polyconvex envelope at \(\hat{\nu}\) in the set \(\Sigma_{\delta}\) which are sometimes also referred to as volume fractions. For a detailed description of the algorithm see [28], a MATLAB and Python implementation can be found under https://github.com/TmNmr/SVPC.

In this paper, we use a python implementation of this algorithm which uses the linprog function from scipy.optimize for solving the linear program. As discretisation of the signed singular value space, the shifted lattice generated by \(\Sigma_{\delta} = (\delta\, \mathbb{Z}^{d} + \tfrac{\delta}{2} \, \mathbb{1}_{d}) \cap [-r, r]^d\) is employed. The lattice parameters \(r\) and \(\delta\) were chosen such that an absolute error of order \(10^{-5}\) can be predicted, see e.g. [28]. Moreover the symmetry induced by the \(\Pi_d\)-invariance is exploited, reducing the computation to grid points \(\hat{\nu} \in \Sigma_{\delta}\) which are located in the signed singular value cone \(0 \leq \lvert \nu_1 \rvert \leq \nu_2 \leq \ldots \leq \nu_d \leq r\) and exploitation of determinant constraints reduces the evaluation domain further to \(0 \leq \nu_1 \leq \nu_2 \leq \ldots \leq \nu_d \leq r\). It should be noted that for points close to the coordinate axes, a sufficiently fine resolution of the signed singular value space is required due to the determinant constraints, that is why finer lattice widths \(\delta\) are employed for these regions. Nevertheless, in what follows, we denote by \(\Phi^{{\operatorname{pc}}}_{\delta}\) simply the approximation of the polyconvex envelope by the SVPC LP algorithm. In addition, since SVPC LP delivers independent point approximations, the computation of the learning data is performed in parallel.

5 Mathematical Benchmark Examples↩︎

To show the potential of the approximation approaches of the property-preserving neural networks as introduced 4, we consider well-known examples in the context of relaxation algorithms in both two and three spatial dimensions. We consider experiments of increasing complexity: parameter-independent cases in two and three spatial dimensions, followed by parameter-dependent problems. Throughout all of the numerical experiments, we denote the analytical polyconvex envelope by \(\Phi^{{\operatorname{pc}}}\) and by \(\Phi^{{\operatorname{pc}}}_{{\rm pred}}\) the neural network prediction.

5.1 The Kohn–Strang–Dolzmann Function↩︎

This function was studied in [54], [55], subsequently modified to achieve continuity in [17], [18] and further studied in [24] and serves as a first benchmark example to illustrate the capability of the neural network approach to approximate the polyconvex envelope. We consider the function \(W \colon \mathbb{R}^{2\times 2} \to \mathbb{R}\) , defined as \[W(F)= \begin{cases} 1+\lvert F \rvert^2 & \text{if } \, \lvert F \rvert \geq \sqrt{2}-1,\\ 2\sqrt{2} \, \lvert F \rvert & \text{else}, \end{cases}\] where \(\lvert F \rvert \mathrel{\vcenter{:}}= (\sum_{i,j=1}^{d} F_{ij}^2)^{1/2}\) denotes the Frobenius norm of the matrix \(F\). The polyconvex envelope of \(W\) is explicitly known, see e.g. [17], and reads \[W^{{\operatorname{pc}}}(F)= \begin{cases} 1+\lvert F \rvert^2 & \text{if } \, \rho(F) \geq 1, \\ 2 \,(\rho(F) - \lvert {\det}(F) \rvert) & \text{else}, \end{cases}\] where \(\rho(F) \mathrel{\vcenter{:}}=\sqrt{\lvert F \rvert^2 + 2 \lvert {\det}(F) \rvert}\). The functions \(W\) and \(W^{{\operatorname{pc}}}\) are isotropic; consequently, they can be rewritten in terms of the signed singular values and reduce to \(\Phi\), \(\Phi^{{\operatorname{pc}}} \colon \mathbb{R}^2 \to \mathbb{R}\) with \[\Phi(\hat{\nu}) = \begin{cases} 1+ \nu_1^2 +\nu_2^2 & \text{if } \, \sqrt{\nu_1^2 +\nu_2^2} \geq \sqrt{2}-1,\\ 2\sqrt{2} \, \sqrt{\nu_1^2 +\nu_2^2} & \text{else}, \end{cases}\] and \[\Phi^{{\operatorname{pc}}}(\hat{\nu})= \begin{cases} 1+ \nu_1^2 +\nu_2^2 & \text{if } \, \rho(\hat{\nu}) \geq 1, \\ 2 \, (\lvert \nu_1 \rvert + \lvert \nu_2 \rvert - \lvert \nu_1 \nu_2 \rvert) & \text{else}, \end{cases} \label{eq:PhipcKSD}\tag{21}\] where \(\rho(\hat{\nu}) = \lvert \nu_1 \rvert + \lvert \nu_2 \rvert\).

5.1.1 Network Architecture and Learning Data↩︎

In this example, we implement a FICNN which consists of two hidden layers, with 10 and 20 neurons, respectively, as presented in 3. The learning domain for \(\hat{\nu}\) is defined as the box \([-\overline{\nu},\overline{\nu}]^2\) with \(\overline{\nu}=1.05\). For the training, the set of signed singular values \(\hat{\nu}\) used as input is generated by discretising each axis with \(751\) points, so that the point 0 is included, making \(751^2\) training data points in total. Instead of using a uniform distribution of the points, we use a local refinement towards the origin, utilising a quadratic transformation. This refinement approach increases the data density near the point of non-differentiability of both functions \(\Phi\) and \(\Phi^{{\operatorname{pc}}}\). For the validation dataset, we randomly sample 169200 points within the domain \([-\overline{\nu},\overline{\nu}]^2\), accounting for \(30 \%\) of the training dataset size. The target values, i.e. the polyconvex envelope, are obtained by evaluation of the analytical function \(\Phi^{\operatorname{pc}}\) as in 21 . For the loss function, we choose the penalty parameters \(\lambda_{{\rm ineq}} = 50\) and \(\lambda_{{\rm sym}} = 10\).

5.1.2 Numerical Results↩︎

On average, for a single realisation, the training process requires \(6 \pm 1\) minutes to complete \(32 \pm 7\) epochs. The final training loss is \(1.3 \times 10^{-3}\pm 2.1 \times 10^{-4}\), while the validation loss reaches \(1.3 \times 10^{-3} \pm 2.3 \times 10^{-4}\). The learning curves for both training and validation losses on one exemplary realisation are presented in 3. In addition to the loss function 17 , the individual contributions to the overall validation loss \(\mathcal{L}\), i.e. \(\mathcal{L}_{{\rm mse}}\) 15 , \(\mathcal{L}_{{\rm ineq}}\) 16 and \(\mathcal{L}_{{\rm sym}}\) 18 , including the scaling by the penalisation parameters, are plotted, showing that the symmetry as well as the inequality property is taken into account and improved during the training process. In particular, the contributions of \(\mathcal{L}_{{\rm ineq}}\) and \(\mathcal{L}_{{\rm sym}}\) are significantly smaller than that of \(\mathcal{L}_{{\rm mse}}\), substantiating the weak enforcement of the corresponding properties.

1 presents the approximation errors of the predictions, where the analytical polyconvex envelope is used as a ground truth. These errors are computed on a uniform \(100 \times 100\) discretisation of the domain \([-\overline{\nu},\overline{\nu}]^2\). The results indicate that the predicted envelopes deviate by about 2% to 3% from the analytical polyconvex envelope \(\Phi^{{\operatorname{pc}}}\), demonstrating an accuracy sufficient for the intended engineering applications. The approximation quality is further illustrated in 4, where two one-dimensional cross-sections of the predicted polyconvex envelope are depicted. It is important to note that the neural network successfully captures the kink of the function \(\Phi^{{\operatorname{pc}}}\) at the point \((0,0)\) with high accuracy. Further, it is notable that the neural network based approximation captures the nonconvexity along the diagonal cross-section of the envelope, rendering it a consistent polyconvex function. It is observable that the standard deviation \(\sigma\) of the different network realisations is quite negligible when it comes to the approximation accuracy. Specifically, the neural network implemented in this example contains 344 parameters, comprising both weights and biases. With the same number of parameters using the conventional approach, it would only be possible to store \(19 \times 19\) grid values, and then approximate the value at any point by interpolation. While this standard method remains feasible for simple cases, it can quickly become intractable for more complex scenarios, such as parameter-dependent functions or real engineering problems.

5.2 Saint Venant–Kirchhoff Energy Density↩︎

Building on the previous two-dimensional example, we now move to a function well-known in elasticity simulations and illustrate the feasibility of our approach in three spatial dimensions. Specifically, we consider the Saint Venant–Kirchhoff model in three spatial dimensions with determinant constraints given by \[\label{eq:WSTVKdetconstr} W(F) = \begin{cases} \frac{\mu}{4} \, \lvert F^T F - \mathbb{I}\rvert^2 + \frac{\lambda}{8} \left(\lvert F\rvert^2 - 3 \right)^2 & \text{ if } \det F > 0, \\ \infty & \text{ else.} \end{cases}\tag{22}\] In contrast to the unconstrained case 2 , where the polyconvex envelope is known analytically [56], no analytical polyconvex envelope is available in the determinant-constrained setting, making numerical approximations necessary. Accordingly, its evaluation within an offline–online framework necessitates the storage of a large set of values on a three-dimensional grid, followed by interpolation. This example thus allows to illustrate a practical application of our approach, namely the efficient compression of a precomputed polyconvex envelope by a properties-preserving neural network, which in turn enables accurate and computationally efficient online prediction.

In order to learn \(\Phi^{{\operatorname{pc}}}\), we focus on the signed singular value reformulation of 22 , which is given by \[\label{eq:PhiSTVKdetconstr} \Phi(\hat{\nu}) = \begin{cases} \frac{\mu}{4} \, \sum_{i = 1}^{3} (\nu_i^2 - 1)^2 + \frac{\lambda}{8} \, (\lvert\hat{\nu}\rvert^2 - 3)^2 & \text{ if } \nu_1 \,\nu_2\,\nu_3 > 0, \\ \infty & \text{ else.} \end{cases}\tag{23}\] Within this formulation, we consider the material parameters \(\mu = 0.4\) and \(\lambda = 0.4\) for the numerical experiments.

5.2.1 Network Architecture and Learning Data↩︎

We implement a FICNN as described in 1, noting that the input layer is now of dimension \(k_3 = 7\), with three hidden layers consisting of 20, 25 and 30 neurons, respectively, leading to a total of 1888 trainable parameters. To generate the learning data, we approximate the polyconvex envelope of the three-dimensional function \(\Phi\) 23 by the SVPC LP algorithm described in 4.6, with the parameters \(\delta = 0.047\) and \(r = 1.5\). The learning domain for \(\hat{\nu}\) is defined as the box \([\nu_{\min}, \nu_{\max}]^3\) with \(\nu_{\min} = 0.4\) and \(\nu_{\max} = 1.4\). The learning data is generated by uniformly discretising each coordinate axis within the box by \(50\) points and considering all possible \(\Pi_d\) permutations of this discretisation, resulting in 1000000 data points covering the entire signed singular value space, half of which yield finite function values of \(\Phi\). The penalty parameters in the loss function 17 are chosen as \(\lambda_{{\rm sym}} = 1.0\) and \(\lambda_{{\rm ineq}} = 2.0\).

5.2.2 Numerical Results↩︎

On average, for a single realisation, the training process requires \(17 \pm 6\) minutes to complete \(32 \pm 10\) epochs. The final training loss is \(1.15 \times 10^{-5} \pm 2.4 \times 10^{-6}\), while the final validation loss reaches \(1.73 \times 10^{-5} \pm 1.2 \times 10^{-5}\). The mean squared error 15 between the neural network approximation and the learning data, i.e. reference polyconvex envelope \(\Phi^{{\operatorname{pc}}}_{\delta}\), accounts to \(5.03 \times 10^{-6} \pm 1.3 \times 10^{-6}\). The mean error on the learning data is \(0.002 \pm 2.3 \times 10^{-4}\), the relative quadratic error is \(0.020 \pm 2.5 \times 10^{-3}\) and the relative maximum error is \(0.027 \pm 5.9 \times 10^{-3}\).

5 shows three cross-sections, \((\nu_1, \nu_1, \nu_1)\)-axis (left), \((\nu_1, \nu_1, 1)\)-axis (middle) and \((\nu_1, 1, 1)\)-axis (right), for the predictions of the neural network model, corresponding to triaxial, biaxial and uniaxial deformations in the signed singular value space, respectively. The reference polyconvex envelope \(\Phi^{{\operatorname{pc}}}_{\delta}\) in 5 is computed by the SVPC LP algorithm, see [28] and illustrated in 4.6 using the discretisation parameters \(r = 1.5\) and \(\delta = 0.025\), and evaluated for \(\nu_1 \in [\nu_{\min}, \nu_{\max}] = [0.4, 1.4]\) uniformly discretised in 200 points. The resulting points on the cross-sections are not part of the learning data set, thus demonstrating the interpolation quality of the network. On all of the three illustrated cross-sections, the maximum absolute error is found to be equal to \(0.0092\), which corresponds to a relative maximum error of \(2\%\). This error is of the same order as the error observed in the learning data for the entire box \([\nu_{\min}, \nu_{\max}]^3\). These numerical results underline the compression potential of the neural-network-based representation approach. In particular, the proposed neural network compresses the large learning data set, consisting of \(\sim 100^3\) point values of \(\Phi^{{\operatorname{pc}}}_{\delta}\), into only 1888 trainable parameters, corresponding to merely a \(12^3\) grid.

5.3 A Parametric Example: The Generalised Kohn–Strang–Dolzmann Function↩︎

After the previous initial benchmark examples, we consider a more complex case, i.e. a two-parameter-dependent family of functions. The following example was studied in [57], [58] and modified to achieve continuity. We consider the function \(W \colon \mathbb{R}^{2\times 2} \times \mathbb{R}_{+} \times \mathbb{R}_{+}\setminus\{0\} \rightarrow \mathbb{R}\), defined as \[W(F; \lambda, \alpha)= \begin{cases} \lambda + \alpha \, \lvert F \rvert^2 & \text{if } \, \lvert F \rvert \geq \sqrt{\frac{\lambda}{\alpha}} \, (\sqrt{2}-1), \\ 2\sqrt{2\, \lambda \,\alpha} \, \lvert F \rvert & \text{else}, \end{cases}\] where \(\lvert F \rvert \mathrel{\vcenter{:}}= (\sum_{i,j=1}^{d} F_{ij}^2)^{1/2}\). The polyconvex envelope of \(W\) is known analytically in closed form and reads \[W^{{\operatorname{pc}}}(F; \lambda, \alpha)= \begin{cases} \lambda + \alpha \, \lvert F \rvert^2 & \text{if } \, \rho(F) \geq \sqrt{\frac{\lambda}{\alpha}}, \\ 2 \sqrt{\lambda \, \alpha} \, \rho(F) - 2 \alpha \, \lvert {\det}(F) \rvert & \text{else}, \end{cases}\] where \(\rho(F) \mathrel{\vcenter{:}}= \sqrt{\lvert F \rvert^2 + 2\, \lvert {\det}(F) \rvert}\). The functions \(W\) and \(W^{{\operatorname{pc}}}\) are isotropic, and, rewritten in terms of the signed singular values, they reduce to \(\Phi\), \(\Phi^{{\operatorname{pc}}}\colon \mathbb{R}^2 \times \mathbb{R}_{+} \times \mathbb{R}_{+} \setminus\{0\}\to \mathbb{R}\) with \[\Phi(\hat{\nu}; \lambda, \alpha) = \begin{cases} \lambda + \alpha \, (\nu_1^2 +\nu_2^2) & \text{if } \, \sqrt{\nu_1^2 +\nu_2^2} \geq \sqrt{\frac{\lambda}{\alpha}} \, (\sqrt{2}-1),\\ 2\sqrt{2\, \lambda \, \alpha} \, \sqrt{\nu_1^2 +\nu_2^2} & \text{else}, \end{cases}\] and \[\Phi^{{\operatorname{pc}}}(\hat{\nu}; \lambda, \alpha)= \begin{cases} \lambda + \alpha(\nu_1^2 +\nu_2^2) & \text{if } \;\rho(\hat{\nu}) \geq \sqrt{\frac{\lambda}{\alpha}},\\ 2 \sqrt{\lambda \, \alpha} \, (\lvert \nu_1 \rvert + \lvert \nu_2 \rvert) - 2 \alpha \, \lvert \nu_1 \nu_2 \rvert & \text{else}, \end{cases} \label{eq:PhipcGKSD}\tag{24}\] where \(\rho(\hat{\nu}) = \lvert \nu_1 \rvert + \lvert \nu_2 \rvert\).

5.3.1 Network Architecture and Learning Data↩︎

For this parameter-dependent example, we implement a PICNN where each path consists of three hidden layers with 10, 20 and 20 neurons, respectively. The convex input, denoted by \(\hat{m}\), represents the minors of the signed singular values, i.e. \(m(\hat{\nu})\), while the nonconvex inputs, denoted as vector \(\zeta\), correspond to the two parameters \(\lambda\) and \(\alpha\). The learning domain for the parameters \(\lambda\) and \(\alpha\) is \([1,2]\), while for the signed singular value input \(\hat{\nu}\) the learning domain is defined as \([-\overline{\nu}, \overline{\nu}]^2\), with \(\overline{\nu}=1.5\), ensuring in particular that \(\overline{\nu} \geq \sqrt{\lambda/\alpha}\). For the training dataset, the set of points in the signed singular value space is obtained by discretising each axis into 251 points, so that the point \(0\) is included, with a local refinement towards the origin using a quadratic transformation, as in the example in 5.1. The parameter \(\lambda\) takes values from the discrete set \(\{1,1.2, \ldots, 1.8,2\}\) and \(\alpha\) follows the same parameter discretisation. These discretisations lead to a total of 2268036 points in the learning dataset. For the validation dataset, \(30 \%\) of these points are randomly selected, leaving \(70 \%\) for the training data set. The target values are computed using the analytical function \(\Phi^{{\operatorname{pc}}}\), cf. 24 .

6 illustrates a one-dimensional representation of the target values for the training dataset along the diagonal \((\nu_1, \nu_1)\)-axis. For the loss function \(\mathcal{L}\) from 17 , we choose the penalty parameters \(\lambda_{{\rm ineq}}=50\) and \(\lambda_{{\rm sym}}=20\).

5.3.2 Numerical Results↩︎

On average, for a single realisation, the training process requires \(27 \pm 5\) minutes to complete \(22 \pm 4\) epochs. The final training loss is \(8.07 \times 10^{-3} \pm 1.6 \times 10^{-3}\), while the validation loss reaches \(8.57 \times 10^{-3} \pm 1.4 \times 10^{-3}\). The learning curves for one selected realisation are plotted in 6.

7 illustrates the relative errors between the predictions and the analytical polyconvex envelopes for different sets of parameters \((\lambda, \alpha)\), including both values present in the training set and those outside of it. These errors are computed on a uniform \((100 \times 100)\)-discretisation of the domain \([-\overline{\nu},\overline{\nu}]^2\). Across all considered parameter sets, the errors remain consistently low, ranging between 2% and 4%. Additionally, 8 provides one-dimensional cross-sections of these predictions. Notably, the numerical experiments indicate that errors are primarily localised at the domain boundaries, while the neural network successfully captures the kink at the point \((0,0)\) with an accuracy sufficient for typical engineering applications. Furthermore, 8 demonstrates that the neural network is capable of making accurate predictions even for values of \(\hat{\nu}\) outside the training domain, defined as \([-1.5, 1.5]^2\). In particular, the predictions remain accurate for \(\hat{\nu}\) in \([-1.8, 1.8]^2\). This highlights the network’s ability to extrapolate reliably in the signed singular values argument. The behaviour in the extrapolation regime is particularly relevant in finite element simulations, where unreasonably large strains may arise during Newton iterations.

A further extrapolation analysis shows that the implemented neural network performs well even for parameter pairs \(\{\lambda, \alpha\}\) in the range \([1, 2]^2 \cup [1.5, 2.5]^2\), i.e. also for pairs outside the training domain, thereby indicating a degree of extrapolation capability in terms of parameters as well (but inevitably with a loss in the precision). An example of such behaviour is illustrated in 9. Beyond this range, however, the prediction quality deteriorates noticeably. Nevertheless, extrapolation in these parameters is less critical than extrapolation in the signed singular values since the parameters are either known (and possibly constant) or bounded (see e.g. 6) for a given model. Apart from the extrapolation capacity, the crucial aspect here is its compression capability within the training range, which will become even more apparent in higher-dimensional parameter spaces. In particular, the neural network implemented in this example contains 3291 parameters, comprising both weights and biases. With the same number of parameters using the conventional approach, it would only be possible to store a \(19 \times 19\) grid, with 3 values for each \(\lambda\) and \(\alpha\). With this limited grid size, it is evident that the full range of data cannot be recovered with the same precision as is achievable with the neural networks presented here.

6 Engineering Example: Isotropic Damage Problem↩︎

Having validated our approach with (mathematical) benchmark problems, we return to the isotropic damage problem introduced in 2. In order to fit this model into our setting, we rephrase it in signed singular value formulation.

6.1 Reformulation of the Isotropic Damage Problem↩︎

The function \(W\) from 4 is dependent on \(d\times d + 1\) parameters and on the \(d \times d\)-deformation gradient. Due to the isotropy of \(W\), it can be recast into a signed singular value formulation utilising \(\hat{\nu} \in \mathbb{R}^d\). To stress the dependence on the signed singular values, the function \(\varphi\) is now employed to rewrite 5 as \[\varphi(\hat{\nu}, \alpha) = (1 - D(\alpha)) \, \varphi^{0}(\hat{\nu}).\] Within this formulation, the function \(\varphi^{0}\) denotes the signed singular value formulation of the isotropic undamaged energy density \(\psi^{0}\) from 5 , i.e. \(\varphi^{0}\) and \(\psi^{0}\) are related by \[\varphi^0(\hat{\nu}) = \psi^0(\operatorname{diag}(\hat{\nu})) \qquad\quad \text{and} \qquad\quad \psi^0(F) = \varphi^0(\nu(F)).\] Consequently, the pseudo-time incremental energy density \(W\) in the signed singular value formulation, denoted by \(\Phi\), for the time step \(k + 1\) reads \[\begin{align} \label{eq:Phidamage} \begin{aligned} \Phi(\hat{\nu}_{k+1}; \hat{\nu}_{k}, \alpha_{k}) & = \varphi(\hat{\nu}_{k+1},p(\hat{\nu}_{k + 1}; \alpha_{k})) - \varphi(\hat{\nu}_k,\alpha_k) \\ &\qquad + p(\hat{\nu}_{k + 1}; \alpha_{k})\, D(p(\hat{\nu}_{k + 1}; \alpha_{k})) - \alpha_k \, D(\alpha_k) - \overline{D}(p(\hat{\nu}_{k + 1}; \alpha_{k})) + \overline{D}(\alpha_k). \end{aligned} \end{align}\tag{25}\] The parameter dependence in this function can be given the following interpretation. The scalar parameter \(\alpha_{k}\) plays still the same role as the internal variable, while the vector \(\hat{\nu}_{k} \in \mathbb{R}^{d}\) belongs to the signed singular values of the deformation gradient \(F_{k}\) from the previous pseudo-time step. Within 25 , the internal variable evolution is written explicitly using the path function in the signed singular value formulation as \[\label{eq:pathfunctionIsotropic} \alpha_{k+1} = p(\hat{\nu}_{k + 1}; \alpha_{k}) = \begin{cases} \varphi^{0}(\hat{\nu}_{k + 1}) & \text{if } \, \varphi^{0}(\hat{\nu}_{k + 1}) > \alpha_{k}, \\ \alpha_{k} & \text{else.} \end{cases}\tag{26}\] The formulation as stated in 25 significantly reduces dimensionality in both the \(\hat{\nu}_{k+1}\) argument as well as the \(\hat{\nu}_{k}\) parameter dependence, opening the possibility for efficient parameter-dependent polyconvexification of the function \(\Phi(\hat{\nu}_{k+1}; \hat{\nu}_{k}, \alpha_{k})\) in the argument \(\hat{\nu}_{k+1}\). For the neural network, this leads to a reduction in parameter space from \(d \times d + 1\) to \(d + 1\) dimensions, and a reduction in the minors input argument from dimension \(K_d\) to \(k_d\), corresponding to the dimension of the vector \(m(\hat{\nu}_{k+1})\).

The damage parameters \(d_0\) and \(d_\infty\) are set to \(d_{0}=0.5, d_{\infty}=0.99\) in our numerical experiments. The Lamé constants \(\lambda\) and \(\mu\) of the materials in 2 and 3 are set to \(\lambda=0\), \(\mu=0.5\), as in [22].

6.2 Normalisation by a Splitting Approach↩︎

At this stage, the incremental energy density \(\Phi\) depends on \(d + 1\) parameters and shows significant variation over the parameter domain \(\hat{\nu}_{k}\) and \(\alpha_{k}\). This pronounced separation between function curves is a challenge for neural networks as it hinders efficient learning. Large gaps between function values can prevent smooth interpolation and generalisation, making it difficult to capture underlying patterns during training—unless a large amount of data is used, which becomes intractable even in two spatial dimensions. Although \(\Phi\) depends only on \(d + 1\) parameters, the parameter \(\hat{\nu}_{k}\) requires a discretisation as fine as the discretisation for the argument \(\hat{\nu}_{k+1}\) since they play a similar role, making the learning computationally infeasible.

To overcome these difficulties, we take advantage of the structure of the pseudo-time incremental energy density function. The function \(\Phi\) from 25 can be expressed as \[\Phi(\hat{\nu}_{k+1}; \hat{\nu}_{k}, \alpha_{k}) = \tilde{\Phi}(\hat{\nu}_{k+1}; \alpha_{k}) + \Phi_{\rm shift}(\hat{\nu}_{k}, \alpha_{k}),\] where \(\tilde{\Phi} \colon \mathbb{R}^{d} \times \mathbb{R}\to \mathbb{R}_{\infty}\) and \(\Phi_{\rm shift}\colon \mathbb{R}^{d} \times \mathbb{R}\to \mathbb{R}_{\infty}\) are defined as \[\label{eq:Phinormalised} \begin{align} \tilde{\Phi} (\hat{\nu}_{k+1}; \alpha_{k}) & \mathrel{\vcenter{:}}= \varphi(\hat{\nu}_{k+1}, p(\hat{\nu}_{k + 1}; \alpha_{k})) \\ & \qquad + p(\hat{\nu}_{k + 1}; \alpha_{k})\, D(p(\hat{\nu}_{k + 1}; \alpha_{k})) - \alpha_k \, D(\alpha_k) - \overline{D}(p(\hat{\nu}_{k + 1}; \alpha_{k})) + \overline{D}(\alpha_k) \end{align}\tag{27}\] and \[\Phi_{\rm shift}(\hat{\nu}_{k}, \alpha_{k}) \mathrel{\vcenter{:}}= - \varphi(\hat{\nu}_k,\alpha_k),\] respectively. It should be stressed that the function \(\Phi_{\rm shift}\) is independent of \(\hat{\nu}_{k+1}\), hence only dependent on the parameters and constant in the convexification argument \(\hat{\nu}_{k + 1}\). Assuming the function \(\varphi^{0}\) is normalised in the sense that \(\inf_{\hat{\nu}} \varphi^{0}(\hat{\nu}) = \varphi^{0}(\mathbb{1}_d) = 0\), the function \(\tilde{\Phi}\) is also normalised, i.e. \(\inf_{\hat{\nu}_{k+1}} \tilde{\Phi}(\hat{\nu}_{k+1}; \alpha_{k}) = \tilde{\Phi}(\mathbb{1}_d) = 0\), where \(\mathbb{1}_d \in \mathbb{R}^{d}\) denotes the vector containing only ones. Consequently, the polyconvex envelope of \(\Phi\) can be obtained from the polyconvexification of the function \(\tilde{\Phi}\) by \[\label{eq:splitpc} \Phi^{{\operatorname{pc}}}(\hat{\nu}_{k+1}; \hat{\nu}_{k}, \alpha_{k}) = \tilde{\Phi}^{{\operatorname{pc}}}(\hat{\nu}_{k+1}; \alpha_{k}) + \Phi_{{\rm shift}}(\hat{\nu}_{k}, \alpha_{k}).\tag{28}\] Therefore, \(\tilde{\Phi}^{{\operatorname{pc}}}\) should be the focus of an approximation by a neural network or a standard algorithm. Note that this normalisation is domain independent, and just relies on the split 28 , considering the contribution \(\Phi_{{\rm shift}}\) as a shift. Removing the dependence on the previous time step \(\hat{\nu}_{k}\) reduces the polyconvexifaction problem to a one-parameter-dependent family, making the learning feasible. In what follows, the neural networks are trained to predict the function \(\tilde{\Phi}^{{\operatorname{pc}}}\) and the function \(\Phi^{{\operatorname{pc}}}\) is recovered a posteriori by applying the shift \(\Phi_{\rm shift}\) as stated in 28 . Since the polyconvex envelopes for both \(\Phi\) and \(\tilde{\Phi}\) are not known analytically, the ground truth is computed according to 4.6.

6.3 Numerical Results for the Two-Dimensional Saint Venant–Kirchhoff Model↩︎

We consider the function \(\Phi\) from 25 in the Saint Venant–Kirchhoff-based formulation, i.e. \(\psi^{0}\) from 2 is chosen with the material parameters as before. We aim for the representation of the function \(\tilde{\Phi}^{{\operatorname{pc}}}\colon\mathbb{R}^d \times \mathbb{R}_{+} \to \mathbb{R}\), i.e. the polyconvex envelope of the normalised version 27 , by a neural network.

6.3.1 Network Architecture and Learning Data↩︎

In this example, we implement a PICNN consisting of three hidden layers, where the \(\hat{m}\)-path consists of layers with \(30\), \(60\) and \(60\) neurons and the \(\zeta\)-path, i.e. the parameter path, consists of \(15\), \(30\) and \(30\) neurons, respectively. The convex input \(\hat{m}\) represents the minors of the signed singular values \(\hat{\nu}_{k+1}\), i.e. \(m(\hat{\nu}_{k+1})\), while the nonconvex inputs, denoted as \(\zeta\), correspond to the parameter \(\alpha_k\). The learning domain for the parameter \(\alpha_k\) is set to \([0, \alpha_\infty]\) while the learning domain for \(\hat{\nu}_{k+1}\) is defined as \([\nu_{\min},\nu_{\max}]^2\) with \(\nu_{\min}=0.1\) and \(\nu_{\max}=5\) and all permutations included in \(\Pi_d\). The data points in the singed singular value space are obtained by discretising each axis on the interval \([\nu_{\min}, \nu_{\max}]\) by \([\nu_{\min} : 0.005 : 1.2] \cup [1.2 : 0.02 : \nu_{\max}]\), leading to a \(411 \times 411\) grid in \(\mathbb{R}_{+}^d\). Additionally all possible permutations of these points induced by transformations in \(\Pi_d\) are included to extend the data to full \(\mathbb{R}^d\). Since \(\tilde{\Phi}^{{\operatorname{pc}}}(\hat{\nu}; \alpha_{k}) = + \infty\), for \(\nu_1 \, \nu_2 \leq 0\), only the quadrants with positive product \(\nu_1 \, \nu_2\) are considered in the learning data set, resulting in 337842 grid points in the signed singular value space. Notably, the learning domain for \(\hat{\nu}_{k+1}\) covers the part of the signed singular value space associated to positive determinant. For \(\alpha_k\), we choose the set of learning values as \(\mathcal{I}_{\alpha_k} = \{0, 0.1, \dotsc, 1.5, 1.75, 2, 2.25, 2.5, 3, 3.5, 4\}\), i.e. 23 values. These discretisations lead to a learning dataset consisting of 7770366 points. The target values \(\tilde{\Phi}^{{\operatorname{pc}}}_{\delta}\) are computed with the SVPC LP algorithm described in 4.6, with discretisation radius \(r = 5.1\). The lattice width is varied, since points close to the origin require a finer resolution of the computational grid to resolve the growth of the function towards the determinant constrained regime. That is why for \(\hat{\nu}\) with \(\min(\nu_1, \nu_2) \leq 0.125\), a lattice size of \(\delta \approx 0.01\), for \(\hat{\nu}\) with \(\min(\nu_1, \nu_2) \leq 1.25\), a lattice size \(\delta \approx 0.02\) and for all other evaluation points \(\delta = 0.04\) is chosen. For this discretisation strategy the computation of the target values takes approximately 400 CPU hours. 10 illustrates the target values of the learning data along the positive part of the diagonal axis \((\nu_1, \nu_1)\). From the learning dataset, \(70 \%\) of the samples are randomly assigned to the training set, while the remaining \(30 \%\) constitute the validation set. For this experiment the patience parameter is set to \(10\) and the penalty parameters included in the loss function are set to \(\lambda_{{\rm ineq}}=4\) and \(\lambda_{{\rm sym}}=2\).

6.3.2 Numerical Results↩︎

On average, for a single realisation, the training process requires \(1.5\) h \(\pm\) \(18\) minutes to complete \(38 \pm 8\) epochs. The final training loss is \(3.11 \times 10^{-5}\pm 4.9 \times 10^{-6}\), while the validation loss reaches \(3.22 \times 10^{-5} \pm 6.6 \times 10^{-6}\). The learning curves for one of these realisations are presented in 10, including the individual contributions to the validation loss evaluation.

11 illustrates the relative errors between the predictions and the polyconvex envelopes computed by the algorithm from 4.6 for different values of \(\alpha_k\), including values that are outside the training domain. These errors are computed on a uniform \(100 \times 100\) lattice of the domain \([0.1, 5]^2 \subset [\nu_{\min},\nu_{\max}]^2\). Across all considered values \(\alpha_k\), the errors remain consistently low, ranging between 1% and 2%. In particular, it is important to note that the hypothesis stating that the energies become independent of \(\alpha_k\) for \(\alpha_k\geq \alpha_{\infty}\) is verified and validates the choice of \(\alpha_\infty=4\).

Additionally, 12 provides examples of re-shifted predictions, following 28 , which are in good agreement with the polyconvex envelopes computed by the SVPC LP algorithm. In particular, these predictions are obtained for values of \(\alpha_k\) that are not included in the learning set \(\mathcal{I}_{\alpha_k}\) and lie outside the parameter’s training domain, demonstrating good predictability and generalisation capabilities of the network.

For comparison, the computation of one polyconvex envelope on a \(100 \times 100\)-grid on the box \([0.1, 5]^2\) using the SVPC LP algorithm of 4.6, with the same discretisation parameters as above, takes \(49\) minutes on a single CPU (exploiting the symmetry due to \(\Pi_d\)-invariance) while its prediction via the neural network takes only \(0.05\) seconds, emphasising once again the benefits of using a neural-network compression approach for the polyconvexification in engineering applications. In addition, the neural network implemented in this example contains 14800 parameters, comprising both weights and biases. Once again, this number of parameters remains smaller than the storage required for only two polyconvex envelopes on a \(100 \times 100\) grid.

6.4 Numerical Results for the Two-Dimensional Neo-Hookean Model↩︎

We consider the function \(\Phi\) from 25 in the neo-Hookean based formulation, i.e. \(\psi^{0}\) from 3 is chosen and the material parameters are set as before. We aim for the representation of the function \(\tilde{\Phi}^{{\operatorname{pc}}}\colon\mathbb{R}^d \times \mathbb{R}_{+} \to \mathbb{R}\), i.e. the polyconvex envelope of the normalised version 27 , by a neural network.

6.4.1 Network Architecture and Learning Data↩︎

As before, we implement a PICNN consisting of three hidden layers, where the \(\hat{m}\)-path consists of layers with \(30\), \(60\) and \(60\) neurons and the \(\zeta\)-path, i.e. the parameter path, consists of \(15\), \(30\) and \(30\) neurons, respectively. The convex input \(\hat{m}\) represents the minors of the signed singular values \(\hat{\nu}_{k+1}\), i.e. \(m(\hat{\nu}_{k+1})\), while the nonconvex inputs, denoted as \(\zeta\), correspond to the parameter \(\alpha_k\). The learning domain for the parameter \(\alpha_k\) is set to \([0, \alpha_\infty]\) while the learning domain for \(\hat{\nu}_{k+1}\) is defined as \([\nu_{\min},\nu_{\max}]^2\) with \(\nu_{\min}=0.55\) and \(\nu_{\max}=18\) and all permutations included in \(\Pi_d\). The data points in the singed singular value space are obtained by discretising each axis on the interval \([\nu_{\min}, \nu_{\max}]\) by \([\nu_{\min} : 0.015 : 1.75] \cup [1.75 : 0.05 : 18]\), leading to a \(406 \times 406\) grid in \(\mathbb{R}_{+}^d\). Additionally, all possible permutations of these points induced by transformations in the symmetry group \(\Pi_d\) are included to extend the data to full \(\mathbb{R}^d\). Since \(\tilde{\Phi}^{{\operatorname{pc}}}(\hat{\nu}; \alpha_{k}) = + \infty\), for \(\nu_1 \, \nu_2 \leq 0\), only the quadrants with positive product \(\nu_1 \, \nu_2\) are considered in the learning data set, resulting in 329672 grid points in the signed singular value space. For \(\alpha_k\), we choose the set of learning values as \(\mathcal{I}_{\alpha_k} = \{0, 0.1, \dotsc, 1.5, 1.75, 2, 2.25, 2.5, 3, 3.5, 4\}\), i.e. 23 values. These discretisations lead to a learning dataset consisting of 7582456 points.

The target values \(\tilde{\Phi}^{{\operatorname{pc}}}_{\delta}\) are computed with the SVPC LP algorithm described in 4.6, with discretisation radius \(r = 20\). The lattice width is varied, since points close to the origin require a finer resolution of the computational grid to resolve the growth of the function towards the determinant-constrained regime. That is why for \(\hat{\nu}\) with \(\min(\nu_1, \nu_2) \leq 1.25\), a lattice size \(\delta \approx 0.039\) and for all other evaluation points \(\delta = 0.078\) is chosen. For this discretisation strategy the computation of the target values takes approximately 600 CPU hours. 13 illustrates the target values of the learning data along the positive part of the diagonal axis \((\nu_1, \nu_1)\). From the learning dataset, \(70 \%\) of the samples are randomly assigned to the training set, while the remaining \(30 \%\) constitute the validation set. For this experiment, the patience parameter is set to \(10\) and the penalty parameters included in the loss function are set to \(\lambda_{{\rm ineq}}=4\) and \(\lambda_{{\rm sym}}=2\).

6.4.2 Numerical Results↩︎

On average, for a single realisation, the training process requires \(3.2 \pm 1.2\) h to complete \(63 \pm 22\) epochs. The final training loss is \(3.04 \times 10^{-5}\pm 8.2 \times 10^{-6}\), while the validation loss reaches \(3.57 \times 10^{-5} \pm 1.6 \times 10^{-5}\). The learning curves for one of these realisations are presented in 13.

14 illustrates the relative errors between the predictions and the polyconvex envelopes computed by the SVPC LP algorithm for different values of \(\alpha_k\), including values that are outside the training domain. These errors are computed on a uniform \(100 \times 100\) discretisation of the domain \([\nu_{\min},\nu_{\max}]^2\). Across all considered parameter sets, the errors remain consistently low, ranging between 1% and 2%. In particular, it is important to note that the hypothesis stating that the energies become independent of \(\alpha_k\) for \(\alpha_k\geq \alpha_{\infty}\) is verified and validates our choice of \(\alpha_\infty=4\). Additionally, 15 provides an example of re-shifted predictions, i.e. application of 28 . It has to be noted that the predictions are in good agreement with the polyconvex envelopes computed by the SVPC LP algorithm from 4.6.

For comparison, the computation of one polyconvex envelope on a \(100 \times 100\)-grid on the box \([0.55, 18]^2\) using the SVPC LP algorithm of 4.6, with the same discretisation parameters as employed for the learning data generation, takes \(1\) hour on a single CPU (exploiting the symmetry due to \(\Pi_d\)-invariance) while its prediction via the neural network takes only \(0.05\) seconds, highlighting the benefits of using a neural-network compression approach for the polyconvexification in engineering applications. In addition, the neural network implemented in this example contains 14800 parameters, comprising both weights and biases. As before, this number of parameters remains smaller than the storage required for only two polyconvex envelopes on a \(100 \times 100\) grid. These aspects emphasise once again the benefits of using a neural-network-based representation of the polyconvex envelope for parameter-dependent families of functions for applications in engineering problems.

7 Conclusion↩︎

We have demonstrated the effectiveness of a neural network design in predicting polyconvex envelopes with high accuracy and computational efficiency. Our results show that such neural networks can generalise well beyond the learning dataset, enabling fast, multi-query evaluations and real-time computations. Moreover, we have introduced a splitting strategy which decouples the isotropic damage problem from the previous time step state, thereby improving the feasibility and robustness of the training process. Future research will focus on extending this framework to predict not only the polyconvex envelopes but also their derivatives as well as the incorporation of determinant constraints into the neural networks, as relevant for engineering applications in computational mechanics. The results presented in this paper pave the way for complex material simulations in real engineering contexts.

Acknowledgement↩︎

Many of the ideas in this paper were initially formulated and tested during Helena Althoff’s Master’s thesis [59]. Fruitful discussions with Daniel Balzani and Maximilian Köhler are gratefully acknowledged, and some ideas are the result of joint discussions with David Wiedemann.

References↩︎

---

1. The authors gratefully acknowledge funding from the German Research Foundation (DFG) within the Priority Programme 2256 Variational Methods for Predicting Complex Phenomena in Engineering Structures and Materials (project number 441154176, reference IDs PE1464/7-2 and PE2143/5-2). Furthermore, we would like to thank the Bavarian State Ministry of Science and the Arts for funding the Augsburg AI Production Network as part of the High-Tech Agenda Plus.↩︎