¹

Key words. Polyconvexity, input convex neural network, monotonicity, relaxation
AMS subject classifications. 49J45, 49J10, 74G65, 74B20, 68T07

1 Introduction↩︎

Many problems in non-linear elasticity can be formulated as variational minimisation problems of the form \[I(u) = \int_{\Omega} W(\nabla u)\,\mathrm{d}x,\] where \(\Omega \subset \mathbb{R}^d\) in spatial dimension \(d \in \{2,3\}\), denotes the reference domain, \(u\colon \Omega \to \mathbb{R}^d\) an admissible deformation, and \(W \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty} \mathrel{\vcenter{:}}= \mathbb{R}\cup \{\infty\}\) the energy density function. In practical applications, the function value \(W(F)=\infty\) is used to model physically inadmissible deformation states, for example through orientation-preserving determinant constraints.

Many constitutive models arising in non-linear elasticity, phase transformations, damage mechanics and fracture are inherently non-convex, see, for example, [1]–[9]. Consequently, the associated variational problems may fail to admit minimisers, while numerical approximations often exhibit mesh dependence, reduced robustness and pronounced sensitivity with respect to discretisation and material parameters.

A classical remedy is provided by relaxation theory. That is, instead of the original non-convex energy density \(W\), one considers a suitable semiconvex envelope, thereby obtaining a relaxed problem that admits minimisers and captures the effective behaviour of oscillatory minimising sequences. Among the various notions of semiconvexity, polyconvexity, introduced by Ball in the works [1], [10], [11], constitutes a physically meaningful sufficient condition for weak lower semicontinuity of the functional and therefore for the existence of minimisers [12]. Consequently, the polyconvex envelope plays a central role in the relaxation of variational problems.

Explicit analytical representations of semiconvex envelopes are available only in a limited number of special cases, see, for example, [13]–[24]. For general non-convex energy densities, however, closed-form expressions for the corresponding semiconvex envelopes are typically unavailable, necessitating computational relaxation procedures.

This has motivated the development of a broad range of computational semiconvexification methods, see, for example, [25]–[36]. In the context of polyconvexity, dedicated polyconvexification algorithms have been proposed in [37]–[40]. For isotropic functions, dimension-reduced formulations have recently been developed in [41] based on the characterisation of polyconvexity through signed singular values established in [42].

The computation of semiconvex envelopes is computationally demanding, since the convex envelope construction is inherently non-local. Moreover, the associated numerical procedures typically require the discretisation of the \(d \times d\)-dimensional deformation gradient space, leading to a substantial computational burden. In the polyconvex setting, this difficulty is further aggravated by the fact that convexification is performed in a lifted space associated with the minors of the deformation gradient. As a result, the dimensionality of the convexification problem increases substantially, particularly in three spatial dimensions. Although the signed singular value characterisation considerably alleviates this issue for isotropic functions and has enabled efficient dimension-reduced polyconvexification algorithms, repeated evaluations of polyconvex envelopes remain prohibitively expensive. The computational burden becomes particularly pronounced in applications involving concurrent relaxation within parameter-dependent boundary value problem simulations, where polyconvex envelopes must be evaluated repeatedly for varying material parameters and loading states.

Recent developments in machine learning have provided powerful tools for constructing structure-preserving surrogate models. In particular, Input Convex Neural Networks (ICNNs) [43] enable the incorporation of convexity directly into the network architecture and have therefore attracted significant attention in constitutive modelling and computational mechanics. Various neural-network-based formulations enforcing polyconvexity through sufficient criteria have been proposed, see, for example, [44]–[50]. Tailored architectures for isotropic hyperelasticity based on the sufficient and necessary characterisation of polyconvexity established in [42] were introduced in [51] and subsequently extended to incompressible formulations in [52].

While these approaches focus on the direct learning of constitutive models, the compression of precomputed polyconvex envelopes was considered in [53]. There, a property-preserving neural-network architecture was proposed for the compression of polyconvex envelopes of isotropic functions based on the sufficient and necessary characterisation from [42]. The resulting surrogate representation acts on the signed singular value space and preserves the structural properties of the polyconvex envelope through a combination of architectural constraints and suitably designed loss terms.

The present contribution builds upon this framework and introduces an alternative compression strategy based on the classical sufficient polyconvexity criterion by Ball [1] for determinant-constrained isotropic functions. In contrast to [53], the proposed approach operates on the singular value space and therefore benefits from a further reduction of the characterising domain. The corresponding neural-network architecture combines convexity-preserving structures (ICNNs) with additional monotonicity constraints inherent in Ball’s criterion, similar to the monotonicity-preserving architectures considered in [50], [52]. As a consequence, the resulting surrogate automatically satisfies the structural requirements of the underlying sufficient criterion for polyconvexity. However, the additional monotonicity assumptions render the proposed representation more restrictive than the equivalent characterisation from [42]. Consequently, the compressed representation may, in general, provide only a lower bound on the polyconvex envelope. Nevertheless, the reduction from the signed singular value space to the singular value space substantially decreases the computational effort.

The proposed methodology is investigated for determinant-constrained isotropic energy densities whose polyconvex envelopes are not available in analytical form. The numerical experiments demonstrate that the reduced representation yields substantial computational savings while maintaining the accuracy of state-of-the-art approaches based on the full signed singular value characterisation. This makes the proposed surrogate particularly attractive for applications requiring repeated evaluations of polyconvex envelopes.

The remainder of this paper is organised as follows. 2 reviews the relevant criteria for polyconvexity of isotropic functions and the associated envelope characterisations. 3 introduces the corresponding property-preserving neural-network architectures. Finally, 4 presents numerical experiments and compares the proposed compression approach with existing methods.

2 Sufficient and Necessary Criteria for Polyconvexity of Isotropic Functions↩︎

This section introduces the theoretical foundations of polyconvexity as well as sufficient and necessary criteria for the polyconvexity of isotopic functions. Let \(d \in \{2,3\}\) and let \(W \colon \mathbb{R}^{d \times d} \to \mathbb{R}_{\infty}\) denote a energy density function that maps \(d \times d\)-matrices to real scalars or infinity. We denote the determinant of a matrix \(F\) by \(\det(F)\) and its adjugate by \(\operatorname{adj}(F)\). Polyconvexity is characterised through the minors of the deformation gradient, represented by the mapping \(\mathcal{M}\colon \mathbb{R}^{d \times d} \to \mathbb{R}^{K_d}\), where \(K_{2} = 5\) and \(k_{3} = 19\), defined by \[\mathcal{M}(F)= \begin{cases} (F, \, \det(F)) & \text{ if } d = 2,\\ (F, \, \operatorname{adj}(F), \, \det(F)) & \text{ if } d = 3. \end{cases}\] Specifically, 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}\) such that \[V(F) = G(\mathcal{M}(F)) \qquad \text{ for all } F\in\mathbb{R}^{d\times d}.\]

Let \(\mathcal{O}(d)\) denote the orthogonal group, \(\mathcal{SO}(d)\) the special orthogonal group and \(\mathcal{S}(d)\subseteq \{0,1\}^{d\times d}\) the group of permutation matrices. Throughout this work, we focus on isotropic energy densities \(W\colon\mathbb{R}^{d\times d}\to \mathbb{R}_\infty\), i.e. for all \(F\in\mathbb{R}^{d\times d}\) and for all \(R_1, R_2\in \mathcal{SO}(d)\) it holds \[W(F) = W(R_1\, F\, R_2).\] Note that this definition of isotropy incorporates objectivity and full material symmetry, following the convention outlined in [1], and is also referred to as \(\mathcal{SO}(d)\times\mathcal{SO}(d)\)-invariance.

A classical result, see e.g. [12], states that isotropic functions admit a dimension-reduced representation in terms of signed singular values. Based on the signed singular value decomposition, i.e. for all \(F\in\mathbb{R}^{d\times d}\) there exist \(R_1, R_2\in \mathcal{SO}(d)\) and \(\hat{\nu}\in \mathbb{R}^d\) such that \[F = R_1 \, \mathrm{diag}(\hat{\nu}) \, R_2,\] it is possible to derive the dimension reduced representation. In general, the signed singular values \(\hat{\nu}\in \mathbb{R}^d\) are only unique up to transformations included in the symmetry group consisting of orientation preserving signed permutations, formally described by \[\Pi(d)= \left\{\operatorname{diag}(\varepsilon )\, S \in \mathcal{O}(d)\;\middle|\; S \in \mathcal{S}(d),\, \varepsilon \in \{ \pm 1\}^d,\, \varepsilon_{1} \cdots \varepsilon_{d} = 1\right\}.\] A canonical representative might be obtained through the signed singular value mapping. We introduce the singular value mapping \(\sigma\colon\mathbb{R}^{d\times d}\to\mathbb{R}^d_+\), where \(R_{+} = \{x \geq 0\}\), defined by \[\sigma(F) = \left[\sigma_1(F), \, \ldots, \, \sigma_d(F)\right]^\top,\] with \(0\leq \sigma_1(F) \leq \ldots \leq \sigma_d(F)\). The signed singular value mapping \(\nu\colon\mathbb{R}^{d\times d}\to\mathbb{R}^d\) is then defined by \[\nu(F) = \left[\operatorname{sign}(\det(F)) \, \sigma_1(F), \, \sigma_2(F), \, \ldots, \, \sigma_d(F)\right]^\top.\] Every isotropic function \(W\colon\mathbb{R}^{d\times d}\to\mathbb{R}_\infty\) is in one-to-one correspondence with a \(\Pi_d\)-invariant function \(\Phi\colon\mathbb{R}^d\to\mathbb{R}_\infty\), i.e. \(\Phi(\hat{\nu}) = \Phi(S \, \hat{\nu})\) for all \(S\in\Pi(d)\), through the identifications \[W(F) = \Phi(\nu(F)) \quad \text{ for all } F\in\mathbb{R}^{d\times d} \quad \text{ and vice versa } \quad \Phi(\hat{\nu}) = W(\operatorname{diag}(\hat{\nu})) \quad \text{ for all } \hat{\nu} \in \mathbb{R}^{d}.\]

Following [42], polyconvexity can equivalently be transferred to functions \(\Phi\colon\mathbb{R}^d\to\mathbb{R}_\infty\) acting on the signed singular values. To this end, we introduce the vector-valued minor mapping \(m\colon\mathbb{R}^d\to\mathbb{R}^{k_d}\), where \(k_2=3\) and \(k_3=7\), defined by \[m(\hat{\gamma}) = \begin{cases} (\hat{\gamma}_1, \, \hat{\gamma}_2, \, \hat{\gamma}_1 \, \hat{\gamma}_2) & \text{ if } d = 2, \\ (\hat{\gamma}_1, \, \hat{\gamma}_2, \, \hat{\gamma}_3,\, \hat{\gamma}_2 \, \hat{\gamma}_3, \, \hat{\gamma}_3 \, \hat{\gamma}_1, \, \hat{\gamma}_1 \, \hat{\gamma}_2, \, \hat{\gamma}_1 \, \hat{\gamma}_2 \, \hat{\gamma}_3) & \text{ if } d = 3. \end{cases}\] A function \(\Phi\colon\mathbb{R}^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 \[\Phi(\hat{\nu}) = g(m(\hat{\nu})) \qquad \text{ for all } \hat{\nu}\in \mathbb{R}^d.\]

The following theorem provides a necessary and sufficient criterion for the polyconvexity of isotropic functions and establishes the equivalence of polyconvexity in the deformation gradient and signed singular value representations. The theorem below was stated by Wiedemann and Peter in [42], where lower semicontinuity is incorporated explicitly into the characterisation. Since lower semicontinuity is automatically satisfied within our neural-network framework considered here, we omit this aspect and refer to [42] for the details of the statement.

While 1 provides an equivalent characterisation of polyconvex isotropic functions, many explicit constructions in non-linear elasticity are based on the classical sufficient criterion due to Ball [1], which is formulated directly in terms of the singular values.

Moreover, under the assumptions of 2, the result extends naturally to extended-real-valued functions incorporating orientation-preserving determinant constraints. In particular, let \(\Upsilon \colon \mathbb{R}_{+}^{d} \to \mathbb{R}_{\infty}\) be of the form \(\Upsilon(\hat{\sigma}) = \bar{g}(m(\hat{\sigma}))\) for all \(\hat{\sigma} \in \mathbb{R}^d_{+}\) where \(\bar{g}\) satisfies the assumptions from 2. Then the function \[\label{eq:ballcritextended} W(F) = \begin{cases} \Upsilon(\sigma(F)) & \text{ if } \det(F) > 0, \\ \infty & \text{ otherwise} \end{cases}\tag{1}\] is polyconvex on \(\mathbb{R}^{d\times d}\).

While 2 provides a convenient sufficient condition of polyconvex functions, a dimension reduced computation of polyconvex envelopes relies on the characterisation via 1. In particular, the envelope construction can be reduced to the convexification of a lifted function in the signed singular value space. To this end, we briefly recall the definition of the polyconvex envelope.

The polyconvex envelope of the energy density \(W\) is the largest polyconvex function below \(W\), i.e. the mapping \(W^{{\operatorname{pc}}} \colon \mathbb{R}^{d\times d} \to \mathbb{R}_\infty\) defined by \[W^{{\operatorname{pc}}}(F) = \sup \left\{ V(F) \;\middle|\; V \colon \mathbb{R}^{d \times d} \to \mathbb{R}_\infty,\;V \text{ polyconvex, } V \leq W \right\}.\] Correspondingly, for a function \(\Phi \colon \mathbb{R}^d \to \mathbb{R}_\infty\), the polyconvex envelope is defined as the largest polyconvex function below \(\Phi\), namely \[\label{eq:defPhipcenvelope} \Phi^{{\operatorname{pc}}}(\hat{\nu}) = \sup \left\{ \Psi(\hat{\nu}) \;\middle|\; \Psi \colon \mathbb{R}^{d} \to \mathbb{R}_\infty,\;\Psi \text{ polyconvex, } \Psi \leq \Phi \right\}.\tag{2}\]

For isotropic functions \(W\) and the associated \(\Pi(d)\)-invariant functions \(\Phi\), [41] shows that the polyconvex envelopes can be determined through the convexification of the function \(h \colon \mathbb{R}^{k_d} \to \mathbb{R}_{\infty}\) defined on the lifted signed singular values by \[\label{eq:h} h(x) = \begin{cases} W(F) & \text{ if } x = m(\nu(F)), \\ \infty & \text{ otherwise}. \end{cases}\tag{3}\] In particular, for all \(F\in\mathbb{R}^{d\times d}\) the natural identification \[\label{eq:WpcPhipchc} W^{\operatorname{pc}}(F)=\Phi^{\operatorname{pc}}(\nu(F))=h^\mathrm{c}(m(\nu(F)))\tag{4}\] holds.

The characterisation 4 provides a viable strategy for determining the polyconvex envelope via a standard convexification in the lifted signed singular value space. In the following, we consider an alternative representation motivated by the sufficient criterion of 2. Restricting attention to determinant-constrained isotropic functions \(W\colon\mathbb{R}^{d\times d}\to \mathbb{R}_\infty\), we assume that the polyconvex envelope of \(W\) admits a representation of the form \[\label{eq:Wpc61Upsilonpc} W^{{\operatorname{pc}}}(F) = \begin{cases} \Upsilon^{{\operatorname{pc}}}(\sigma(F)) & \text{ if } \det(F) > 0, \\ \infty & \text{ otherwise}, \end{cases}\tag{5}\] where \(\Upsilon^{{\operatorname{pc}}}\) is \(\mathcal{S}(d)\)-invariant and satisfies \[\label{eq:Upsilonpc61barhc} \Upsilon^{{\operatorname{pc}}}(\hat{\sigma}) = \bar{h}^{{\operatorname{c}}}(m(\hat{\sigma})),\tag{6}\] with \(\bar{h}^{{\operatorname{c}}}\colon \mathbb{R}^d_{+} \to \mathbb{R}_{\infty}\) convex and non-decreasing in the first \(k_d-1\) arguments. It should be stressed, that technically \(\Upsilon^{{\operatorname{pc}}}\) is not a polyconvexification. Rather, it serves as the singular value representation of the polyconvex envelope and is a notation used to distinguish it from \(\Phi^{{\operatorname{pc}}}\). Indeed, under the assumption 5 , the singed singular value polyconvex envelope is related by \[\Phi^{{\operatorname{pc}}} (\hat{\nu}) = \begin{cases} \Upsilon^{{\operatorname{pc}}}(\sigma(\operatorname{diag}(\hat{\nu}))) & \text{ if } \nu_{1} \cdot \ldots \cdot \nu_{d} > 0, \\ \infty & \text{ otherwise}. \end{cases}\] The composition with \(\sigma \circ \operatorname{diag}\) ensures that \(\Upsilon^{{\operatorname{pc}}}\) gets the singular values as argument, i.e. the input vector \(\hat{\nu}\) is mapped to the positive octant as it was the convention of the singular value mapping \(\sigma\). Consequently, \(\Upsilon^{{\operatorname{pc}}}\) may be viewed as the restriction of \(\Phi^{{\operatorname{pc}}}\) to the singular value cone.

It is important to note that the representation of the form 5 is an assumption. Indeed, the monotonicity requirements imposed on \(\bar{h}^{{\operatorname{c}}}\) are not implied by the characterisation 4 . If these conditions are not satisfied, the resulting function obtained from the monotone convexification, i.e. a monotonous convex envelope of the function \(h\) from 3 over \(\mathbb{R}^{k_d}_{+}\), the formulation \(\Upsilon^{{\operatorname{pc}}}\) in 5 still characterises a polyconvex function below \(W\), but in general constitutes only a lower bound of the polyconvex envelope \(W^{{\operatorname{pc}}}\).

The representation of the polyconvex envelope of \(W\) via \(\Upsilon^{{\operatorname{pc}}}\) as in 5 provides the foundation for the alternative compression approach considered in this work. While the additional monotonicity requirements may, in general, render the resulting approximation a lower bound on the polyconvex envelope, the associated reduction from the signed singular value space to the singular value space yields a considerable computational advantage and motivates the neural-network architecture introduced in the following section.

3 Property-Preserving Neural Network Architectures↩︎

This section introduces the neural-network-based representations of polyconvex envelopes of isotropic functions based on the characterisations presented in the previous section. In addition to the compression of the signed singular value polyconvex envelope \(\Phi^{{\operatorname{pc}}}\) from 4 , as proposed in [53], we consider the singular value based representation \(\Upsilon^{{\operatorname{pc}}}\) from 5 . The resulting neural network approximations are denoted by \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) and \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), respectively.

The approximation \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) is based on the necessary and sufficient characterisation of 1, whereas \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) follows the representation associated to 2, i.e. the Ball-criterion-based characterisation. Although the latter may, in general, provide only a lower bound on the polyconvex envelope, it benefits from a reduced characterising domain, acting solely on the positive octant in the singed singular value space, the singular value cone. As a consequence, substantially lower computational costs can be achieved while preserving the structural properties required by the corresponding polyconvexity criterion.

Both approximations are realised by the same underlying architecture, so called Input Convex Neural Networks, as proposed in [43]. Following this contribution, convexity is enforced through non-negativity constraints on selected weights together with convex and monotone activation functions. In the present setting for polyconvexity, convexity is required with respect to the minors, i.e. the lifted (signed) singular values.

To describe both variants simultaneously, we denote by \(\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}} \in \{\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}, \Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\}\) the considered approximations in an abstract sense. Furthermore, let \(\hat{\gamma}\in\mathbb{R}^{d}\) denote the corresponding input variable, i.e. \(\hat{\nu} \in \mathbb{R}^{d}\) in the case of \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) and \(\hat{\sigma}\in\mathbb{R}^d_+\) for \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), and define the lifted input of the network by \(\hat{m}= m(\hat{\gamma})\). The application of the minors mapping \(m\) is interpreted as a hard coded feature-extraction layer. For an \(L\)-layer network, the hidden activations \(z_{\ell}\) are defined recursively by \[z_{\ell+1} = \rho_\ell \left(W_\ell^{(z)} z_{\ell} + W_\ell^{(\hat{m})} \hat{m} + b_\ell \right)\] for \(\ell = 1, \ldots, L - 1\) and the output \(z_{\ell} \in \mathbb{R}\). The weights \(W_\ell^{(z)}\), the passthrough layer weights \(W_\ell^{(\hat{m})}\) and biases \(b_\ell\) are collected in the trainable parameter vector \(\theta\), with the typical convention \(z_0 = 0\) and \(W_0^{(z)} \equiv 0\). The overall network output is denoted by \[\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}(\hat{\sigma}; \theta) = \bar{h}^{{\operatorname{c}}}_{\mathcal{NN}}(m(\hat{\sigma});\theta) = z_{L}\] and \[\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}(\hat{\nu}; \theta) = h^{{\operatorname{c}}}_{\mathcal{NN}}(m(\hat{\nu});\theta) = z_{L},\] respectively.

Note that neither \(\Phi_{\mathcal{NN}}^{{\operatorname{pc}}}\) nor \(\Upsilon_{\mathcal{NN}}^{{\operatorname{pc}}}\) explicitly encode the determinant constraint included in the function \(W^{{\operatorname{pc}}}\) as in 5 . Both networks are finite-valued representations of the corresponding characterising functions and the extension by \(\infty\) for \(\det(F)\leq 0\) should be enforced separately during the evaluation of \(W_{\mathcal{NN}}^{{\operatorname{pc}}}\), in agreement with 5 .

The convexity and monotonicity properties of \(h^{{\operatorname{c}}}_{\mathcal{NN}}\) and \(\bar{h}^{{\operatorname{c}}}_{\mathcal{NN}}\) required by 1 and 2 are enforced through suitable restrictions on the network parameters and activations. According to [53] and [43], convexity of the network output with respect to the lifted (signed) singular values input arguments is guaranteed by employing convex and non-decreasing activation functions \(\rho_\ell\) and by constraining the weights in the \(z\)-path to be non-negative, i.e. \(W_\ell^{(z)} \geq 0\) for \(\ell = 1,\ldots,L-1\). Under these assumptions, the functions \(h^{{\operatorname{c}}}_{\mathcal{NN}}\) and \(\bar{h}^{{\operatorname{c}}}_{\mathcal{NN}}\) are convex with respect to their inputs.

For the representation \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), the monotonicity assumptions of 2 must additionally be satisfied. Since the activation functions are already chosen to be non-decreasing, the additional monotonicity of \(\bar{h}^{{\operatorname{c}}}_{\mathcal{NN}}\) in the first \(k_d-1\) arguments is obtained by enforcing the corresponding columns of the weights in the \(\hat{m}\)-path to be non-negative, i.e. \([W_\ell^{(\hat{m})}]_{:,1:k_d-1} \geq 0,\) for \(\ell = 0,\ldots,L-1\).

Consequently, \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) satisfies the convexity requirements of 1, whereas \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) additionally satisfies the monotonicity assumptions required by 2. Therefore, both architectures preserve the structural properties underlying the corresponding polyconvexity criteria by construction.

Following [53], the remaining characteristic properties, namely the symmetry conditions and the envelope inequality, are incorporated into the neural-network models in a weak sense through the loss function. Let \(\mathcal{S}\) denote the associated symmetry group to the considered neural-network representation \(\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}\), i.e. \(\Pi(d)\) in the case of \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) and \(\mathcal{S}(d)\) for \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\). Furthermore, let \((\hat{\gamma}_{i}, \Phi^{{\operatorname{pc}}}(\hat{\gamma}_{i}), \Phi(\hat{\gamma}_{i}))\), for \(i = 1, \ldots, N\) denote the training data points. The symmetry property is enforced through the penalty term \[\mathcal{L}_{{\rm sym}}(\theta; \mathcal{S}) = \frac{1}{N} \sum_{i = 1}^{N} \frac{1}{\lvert \mathcal{S}\rvert} \sum_{S \in \mathcal{S}} \left( \hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}({\hat{\gamma}_{i}};\theta) - \hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}(S \, \hat{\gamma}_i;\theta)\right)^2.\] Moreover, the inequality property in the definition of the polyconvex envelope \(\Phi^{{\operatorname{pc}}} \leq \Phi\), is incorporated through the penalty \[\mathcal{L}_{{\rm ineq}}(\theta) = \frac{1}{N} \sum_{i = 1}^{N}\max \bigl\{\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}(\hat{\gamma}_{i};\theta)- \Phi(\hat{\gamma}_{i}),0\bigr\}^{2}.\] Combining these contributions with the mean squared approximation error for the target envelope function values yields the overall loss function \[\mathcal{L}(\theta) = \mathcal{L}_{{\rm mse}}(\theta) + \lambda_{{\rm sym}} \, \mathcal{L}_{{\rm sym}}(\theta) + \lambda_{{\rm ineq}} \, \mathcal{L}_{{\rm ineq}}(\theta),\] where \(\lambda_{{\rm sym}}\) and \(\lambda_{{\rm ineq}}\) denote scalar penalty parameters associated with the symmetry and inequality contributions, respectively.

4 Numerical Experiments↩︎

In the following section, we focus on a well-established benchmark problem to assess the performance of the proposed compression approach and compare it with the state-of-the-art ansatz introduced in [53]. We consider the Saint Venant–Kirchhoff energy density in three spatial dimensions, incorporating determinant constraints that penalise self-intersection and self-interpenetration. This energy density \(W\colon \mathbb{R}^{3 \times 3} \to \mathbb{R}_{\infty}\) is defined by \[W(F) = \begin{cases} \frac{\mu}{4} \lvert F^\top F - \mathbb{I}_{3} \rvert^2 + \frac{\lambda}{8} \left(\lvert F\rvert^2 - 3 \right)^2 & \text{ if } \det(F) > 0, \\ \infty & \text{ otherwise}, \end{cases}\] where \(\lvert A \rvert = \sqrt{\operatorname{tr}(A^\top A)}\) denotes the Frobenius norm and \(\mathbb{I}_{3} \in \mathbb{R}^{3 \times 3}\) the identity matrix. Exploiting the \(\mathcal{SO}(d)\times\mathcal{SO}(d)\)-invariance of \(W\), the function can be rephrased in the dimension-reduced representation \(\Phi \colon \mathbb{R}^{3} \to \mathbb{R}_{\infty}\), acting on the signed singular values \(\hat{\nu} = [\nu_1, \ldots, \nu_d]^\top \in \mathbb{R}^3\), given by \[\Phi(\hat{\nu}) = \begin{cases} \frac{\mu}{4} \, \sum_{i = 1}^{3} \left(\nu_i^2 - 1\right)^2 + \frac{\lambda}{8} \left(\lvert\hat{\nu}\rvert^2 - 3\right)^2 & \text{ if } \nu_{1} \, \nu_{2} \, \nu_{3} > 0, \\ \infty & \text{ otherwise}. \end{cases}\] For the unconstrained model, the polyconvex envelope is known in closed form, see [22]. However, no analytical representation is known for the determinant-constrained setting considered here, making a numerical approximation of the envelope indispensable.

Owing to the isotropic structure of the energy density, the polyconvex envelope can be approximated efficiently by the signed singular value polyconvexification algorithm based on linear programming (SVPC LP) proposed in [41]. The method performs a convexification of a discrete representation of the function \(h\) from 3 on the minors manifold using a lifted signed singular value discretisation, i.e.  \[W^{{\operatorname{pc}}}_{\delta}(F) = \Phi^{{\operatorname{pc}}}_{\delta}(\nu(F)) = \begin{cases} h^{\operatorname{c}}_{\delta}(m(\nu(F))) & \text{ if } \det(F) > 0, \\ \infty & \text{ otherwise}. \end{cases}\] The resulting numerical approximation of the singed singular value polyconvex envelope is denoted by \(\Phi^{{\operatorname{pc}}}_{\delta}\).

The original energy density \(\Phi\) and its numerical polyconvex envelope approximation \(\Phi^{{\operatorname{pc}}}_{\delta}\) are depicted in 1. Since no analytical envelope is available, this benchmark provides a particularly suitable test case for the proposed compression framework: the numerically computed envelope point values obtained by SVPC LP serve as target values that are subsequently compressed using the structure-preserving neural network approaches introduced before.

Studying the numerical envelope approximation \(\Phi^{{\operatorname{pc}}}_{\delta}\) suggests that the polyconvex envelope admits a representation of the form 1 . This motivates the application of the Ball-criterion-based compression approach, which amounts to learning the associated convex function \(\bar{h}^{{\operatorname{c}}}\) from 6 on \(\mathbb{R}_{+}^{3}\).

Since Ball’s criterion, cf. 2, provides only a sufficient condition for polyconvexity, the additional monotonicity constraints imposed on \(\bar{h}^{{\operatorname{c}}}\), and in the neural network may, in principle, be overly restrictive. Consequently, the corresponding neural network ansatz may represent only a lower bound of the true polyconvex envelope. For the present benchmark, however, the numerical reference approximation \(\Phi^{{\operatorname{pc}}}_{\delta}\) already satisfies the required monotonicity properties and the Ball-criterion-based neural network can be employed without introducing any observable additional approximation error in the experiments presented below.

The numerical setup is as follows. All experiments are conducted in the three-dimensional setting with Lamé parameters \(\mu = 0.4\) and \(\lambda = 0.4\). We consider the proposed compression model \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) from 6 and the state–of–the–art model \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) from 4 for comparison reasons. Both models are implemented in Python using the PyTorch framework. The architectures for both approaches share the same principal structure. After the feature extraction via the hard coded minors layer, the neural network consists of a input layer of dimension \(k_{3} = 7\), corresponding to the lifted signed singular value dimension, followed by three hidden layers with \(10\), \(10\) and \(20\) neurons, respectively, and a scalar (finite) valued output layer.

As activation function, we employ the Softplus function with parameter \(\beta=20\), which is convex and non-decreasing and therefore compatible with the theoretical requirements to ensure convexity. Furthermore, the weights satisfy the structure-preserving architecture constraints introduced in the previous section by projection of the relevant weights onto the positive half space via the Softplus function. With this configuration, each network comprises 648 trainable parameters and apart from the additional monotonicity constraints required for the Ball-criterion-based model \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), both architectures employ essentially the same implementation, ensuring a fair comparison of the respective compression strategies.

Training is performed using the Adam optimiser with learning rate \(\eta=\num{5e-3}\) and batch size \(\num{256}\). Early stopping with a patience of \(10\) epochs is employed, while the maximal number of training epochs is set to 250. The penalty parameters in the loss function \(\mathcal{L}\) are chosen as \(\lambda_{{\rm sym}}=2.0\) and \(\lambda_{{\rm ineq}}=5.0\).

We aim to compress the polyconvex envelope on a bounded subset of the signed singular value space rather than along prescribed deformation paths. Consequently, the learning domain contains the signed singular values corresponding to arbitrary deformation states and is not restricted to specific deformation patterns. The discretisation is based on the box \([\nu_{\min}, \nu_{\max}]^{3} \subset \mathbb{R}_{+}^{3}\) with \(\nu_{\min}=0.4\) and \(\nu_{\max}=1.4\). Each coordinate direction is discretised by \(75\) equidistantly distributed points.

For the Ball-criterion-based architecture \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), the learning domain is given by \([\nu_{\min}, \nu_{\max}]^{3} \subset \mathbb{R}_{+}^{3}\), resulting in a total of 421875 lattice points. In contrast, the architecture \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) operates on the entire signed singular value space, to be more precise, the part of signed singular value space which corresponds to \(\det F > 0\), i.e. \(\{\hat{\nu} \in \mathbb{R}^{d} \mid \nu_1 \, \nu_2 \, \nu_3 > 0\}\). Therefore, all sign-preserving permutations must additionally be considered, i.e. \(\Pi(d)\, [\nu_{\min}, \nu_{\max}]^{3}\), yielding a total of 1687500 data points in \([\nu_{\min}, \nu_{\max}]^{3} \, \subset \mathbb{R}^{3}\). This extension is necessary in order to satisfy the appropriate notion of convexity on the full signed singular value space.

It should be noted that throughout this work, we restrict ourselves to the finite-valued regime of the energy density and do not explicitly consider the extended-real-valued setting. In a practical simulation framework, the determinant constraint could instead be enforced by a dedicated preprocessing layer that detects inadmissible deformation states and returns the penalty value \(\infty\).

The restriction of \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) to the positive quadrant reduces the size of the learning domain by a factor of \(4\), clearly illustrating the domain reduction achieved by the Ball-criterion-based formulation. This reduction substantially decreases the number of data points involved in the training process. The trade-off is the additional monotonicity requirement, which may in general be overly restrictive and consequently yield only a lower-bound approximation of the polyconvex envelope.

The learning data is given by tuples of the form \((\hat{\nu}, \Phi^{{\operatorname{pc}}}_{\delta}(\hat{\nu}), \Phi(\hat{\nu}))\), where the envelope values are obtained by the SVPC LP algorithm. The original function values \(\Phi\) are included in order to enforce the inequality constraint appearing in the loss function. The resulting data set is split randomly into a training set (70%) and a validation set (30%). To reduce the influence of the random data partitions and random weight initialisations, all experiments are repeated for ten independent network realisations.

1 depicts the resulting averaged neural network compressions \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) and \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) on a representative cross-sectional slice in the positive octant. It should be noted that this domain constitutes the natural domain for \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), whereas for \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) it merely represents a restriction of the full signed singular value space used during training. For illustration and comparison reasons, this restriction is used throughout the visualisations below.

Visually, both neural network approximations are nearly indistinguishable and closely match the reference envelope \(\Phi^{{\operatorname{pc}}}_{\delta}\). This observation suggests that, for the determinant-constrained Saint Venant–Kirchhoff model considered here, the Ball-criterion-based ansatz is indeed able to recover the polyconvex envelope. In particular, although \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) is, in general, expected to provide only a lower-bound, no discrepancy between the corresponding compression and the reference envelope is observed.

To complement this qualitative assessment, 1 reports quantitative performance measures of the training process for both approaches. The training and validation losses of both approaches are around the same order of magnitude. However, a direct comparison of the loss values is only partially meaningful due to the differing symmetry notions and the fact that, for instance, the mean squared error is evaluated on different computational domains. For this reason, additional illustrations provided below aim to shed further light on the approximation capabilities of the respective models.

In terms of training duration, the \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) approach typically requires a larger number of epochs, it nevertheless achieves a significantly lower overall training time. In the present configuration, the average training time is approximately 22 minutes for \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) and 35 minutes for \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\), corresponding to a speedup factor of about 1.6. We emphasise that this speedup should be regarded as conservative. Across a range of additional architectures and penalty parameter choices, we observed speedup factors between approximately \(1.4\) and \(6\). The precise value depends on the chosen hyperparameters and training configuration.

2 shows the individual approximation and property errors of the neural network compressions. The top row corresponds to the compression via \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\), while the bottom row shows the results for \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\). The left column depicts the pointwise approximation error, the middle column the symmetry error, and the right column the inequality error, respectively.

As a local measure of non-symmetry of the neural network compression, we consider the pointwise deviation from the symmetrised network output. For a network \(\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}} \in \{\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}, \Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\}\) with associated symmetry group \(\mathcal{S}\in \{\mathcal{S}(d), \Pi(d)\}\), the symmetry error \(e_{{\rm sym}}(\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}; \mathcal{S}) \colon \mathbb{R}^{3} \to \mathbb{R}\) describes the pointwise distance to the symmetrised network output and is defined by \[e_{{\rm sym}}(\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}; \mathcal{S})(\hat{\gamma}) = \hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}(\hat{\gamma}) - \frac{1}{\lvert \mathcal{S}\rvert} \sum_{S \in \mathcal{S}}\hat{\Gamma}^{{\operatorname{pc}}}_{\mathcal{NN}}(S \, \hat{\gamma}).\]

The pointwise approximation errors indicate comparable accuracy for both compression approaches. Moreover, both methods exhibit symmetry and inequality errors of similar magnitude, significantly lower than the pure pointwise approximation error, confirming effective incorporation of the structural properties through the penalty formulation.

3 illustrates three one-dimensional cross sections of the three-dimensional Saint Venant–Kirchhoff example for the neural network compressions \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\) and \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\). Both approaches exhibit a reliable approximation of the reference envelope, with no discernible difference in accuracy along the selected slices. The variability across independent training runs is negligible, as indicated by the virtually vanishing standard deviations \(\sigma_{\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}}\) and \(\sigma_{\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}}\), justifying the presentation of averaged network outputs in the preceding figures. Both compression approaches demonstrate stable and accurate interpolation within the training domain, while maintaining acceptable accuracy in regions outside the training domain.

Overall, both approaches yield accurate and structure-preserving approximations of the polyconvex envelope, with no significant differences in accuracy and in the enforcement of the required physical properties. At the same time, the neural network reduces the 1687500 (resp. 421875) data points for \(\Phi^{{\operatorname{pc}}}_{\mathcal{NN}}\) (resp. \(\Upsilon^{{\operatorname{pc}}}_{\mathcal{NN}}\)) to a representation with only 648 trainable parameters, which corresponds, in terms of storage capacity, to representing the envelope on a \(9 \times 9 \times 9\) lattice and highlights the huge compression potential for precomputed polyconvex envelopes. The Ball-criterion-based formulation is characterised by a consistently reduced training time and, consequently, a significant decrease in computational cost.

5 Conclusion↩︎

The proposed structure-preserving neural network framework provides an efficient and reliable approach for the compression of polyconvex envelopes of isotropic functions based on the classical sufficient criterion for polyconvexity of [1]. Although this compression approach generally provide only a lower bound on the polyconvex envelope, the proposed architecture proved successful for the physically relevant determinant-constrained energy density considered in this work, whose polyconvex envelope is not available in analytical form. Substantial reductions in computational cost can be achieved through the reduction from the signed singular value space to the singular value space while preserving the required structural properties and maintaining the accuracy of state-of-the-art approaches. These computational savings appear particularly attractive in parameter-dependent settings, where they may facilitate concurrent relaxation in boundary value problem simulations, such as damage models, a direction that could be pursued in future work.

Acknowledgement↩︎

The authors would like to express their gratitude for the fruitful scientific discussions with Daniel Peterseim and Malte A. Peter.

References↩︎