Discovering interpretable physical models using symbolic regression and discrete exterior calculus

Exterior calculus (EC) deals with objects called differential forms and their manipulations (differentiation, integration) over smooth manifolds. Compared to classical tensor calculus using index notation, it provides a coordinate-independent description of physical models and unifies several differential operators, which results in more compact and readable equations that easily generalize to manifolds of different dimension [10]. DEC offers an inherently discrete counterpart to EC, designed to preserve some fundamental differential properties. In particular, in the context of DEC, cell-complexes and cochains play the role of smooth manifolds and differential forms, respectively. Here, we provide a very brief introduction to DEC, geared toward the applications we will consider in section 3. We report in appendix A more formalism of DEC definitions and notations. Without loss of generality, we limit the presentation to simplicial complexes. A more detailed treatment of DEC can be found in [10, 11, 20].

Henceforth, $\mathcal{E}^N$ is the embedding N-dimensional Euclidean point space, $\mathcal{V}$ is the associated space of translations associated with $\mathcal{E}^N$ and Lin is the space of linear applications from $\mathcal{V}$ to $\mathcal{V}$ (tensors).

From the physical modeling viewpoint, a simplicial complex provides the topological and geometrical structure that supports the physical phenomenon. A generic simplex is indicated with the letter σ and we use a pair of subscripts to identify each of them, such that $\sigma_{p,i}$ is the i-th p-simplex. Cochains are the discrete counterpart of differential forms, hence they are also known as discrete forms. Specifically, a cochain can describe a physical quantity (differential form) integrated on the corresponding simplex and is thus the proxy for a field. The discrete analogue of the exterior derivative in EC is represented by the coboundary d, which plays a key role in expressing balance statements [12] of physical models. A dual complex is needed to represent the behavior of a physical quantity in the neighborhood of a simplex (see section 4.2 in [12]); a role similar to that of the dual complex is played by, e.g. stencils in the finite difference method.

In addition to the topological operators introduced so far, metric-dependent ingredients are needed to complete the formulation of physical theories, especially as concerns constitutive equations. Indeed, these equations relate primal to dual variables [12, 14], so that it is convenient to introduce an isomorphism between primal and dual cochains, called the discrete Hodge star and denoted by $\star$. Finally, in continuum theories the L² inner product allows to express power expenditures by pairing power-conjugated quantities and thus allows to represent potential energy terms in variational models for physical systems, which will be employed in our numerical experiments. In particular, it acts on a pair of primal/dual cochains and in the following is indicated by $\langle , \rangle$. This inner product allows to define the discrete codifferential δ as the adjoint of the coboundary with respect to this operator.

We tackle SR through genetic programming [21] (GP), an evolutionary technique that explores the space of symbolic models starting from an initial population and using genetic operations. The candidate models (individuals) are represented by expression trees, whose nodes can be either primitives (i.e. functions and operators from a prescribed pool) or terminals (i.e. variables and constants). Specifically, for the primitives we use the standard operations on scalars and the DEC operations introduced in section 2.1, see table B.1.

The initial population is generated using the ramped half-and-half method (minimum height = 2, maximum height = 5) that allows to obtain trees with a variety of sizes [22]. We constrain the initial individuals to have length—the number of nodes of its tree—less or equal than 100 and to contain the unknowns of the problem.

2.2.1. Model discovery strategy

The goal of the GP evolution is to maximize the following fitness measure

$$\begin{align} F\left(I\right) : = -\left(\alpha\,\mbox{MSE}\left(I\right) + \eta\,R\left(I\right)\right), \end{align} \tag{ 1 }$$

where α is a scaling factor set equal to $\{1, 10, 1000\}$ for Poisson, Elastica and Linear elasticity, respectively, MSE is the mean-square error of the solution corresponding to the current individual I and R is a regularization term, and η is the associated hyperparameter. The regularization term aims at preventing overfitting, i.e. high fitness on the training dataset and low fitness on the validation dataset. To favor individuals with low complexity among those who accurately fit the training dataset, we represent the regularization term as $R(I) : = \text{length}(I)$.

Starting from the initial population, the program evolves such a population through the operations of crossover and mutation, similarly to genetic algorithms. Individuals eligible for crossover and mutation are chosen through a binary stochastic tournament selection, where the competing individuals have given survival probabilities such that the weakest individual too has a chance of winning the tournament. In all the experiments we used one-point crossover and a mixed mutation operation that is chosen among uniform mutation, node replacement and shrink (see here for details on these operations). The two survival probabilities for the stochastic tournament and the three probabilities associated to the mixed mutation are hyperparameters tuned as a part of the model selection.

As anticipated, the discrete mathematical framework introduced in the previous section naturally suggests a type system for the primitives of the GP procedure. For example, we cannot sum cochains with different dimensions. To enforce this typing, we use a variant of GP called Strongly typed GP (STGP), which generates and manipulates trees there are type-consistent. Crucially, in our approach the physical consistency of the generated expressions is guaranteed by the use of dimensionless quantities. To evolve the population, we use a $(\mu + \lambda)$-genetic algorithm with $\mu = \lambda$ and overlapping generations, such that top µ individuals with the highest fitness are selected (truncation selection) among the total population including parents and offspring [23].

2.2.2. Datasets

2.2.3. Variational models of physical systems

We implement our framework in two open-source Python libraries: dctkit (Discrete Calculus Toolkit) [24] for the DEC concepts described in section 2.1 and alpine [25] for the SR. Specifically, dctkit leverages the JAX [26] library as a high-performance numerical backend for the manipulation of arrays and matrices encoding cochains and DEC operators. For the minimization of the discrete energies found during the SR process, we use the limited-memory Broyden–Fletcher–Goldfarb–Shanno (LBFGS) algorithm implemented in the pygmo [27] library, with the gradients evaluated via the automatic differentiation feature of JAX. The alpine library is based on the STGP algorithms provided by a custom version of the Distributed Evolutionary Algorithms in Python (DEAP) [28] library and on the distributed computing features of the ray library [29].

In the continuous setting, the Poisson equation can be obtained as the stationarity conditions of the Dirichlet functional [33]. Such a functional contains the L² norm of the gradient of the unknown field. In the discrete setting, we represent the dimensionless unknown as a primal scalar-valued zero-cochain u, and application of the coboundary operator d leads to a primal scalar-valued one-cochain that is the discrete equivalent of the gradient of the continuous field integrated along the edges. Denoting by f the primal scalar-valued zero-cochain representing the (dimensionless) source, we can write the potential energy function as

$$\begin{align} \mathcal{E}_{\text{p}}\left(u\right) : = \frac{1}{2} \langle du, du \rangle - \langle u, f \rangle = \frac{1}{2} \langle \delta du, u \rangle - \langle u, f \rangle. \end{align} \tag{ 2 }$$

It can be shown [34] that the stationarity conditions of $\mathcal{E}_{\text{p}}$ correspond to the discrete Poisson equation. The benchmark consists in recovering equation (2) defined over a two-simplicial complex. In particular, a Delaunay triangulation of the unit square with 230 nodes was constructed using gmsh [35].

For the SR procedure, we generated 12 arrays representing the coefficients of the zero-cochain u by sampling the following scalar fields at the nodes of the mesh:

$$\begin{align} &u_{1,i}\left(x,y\right) : = \left(i + 1\right)\mathrm{e}^{\sin\left(x\right)} + \left(i+1\right)^2 \mathrm{e}^{\cos\left(y\right)}, && i = 0,1,2,3, \end{align} \tag{ 3 }$$

$$\begin{align} &u_{2,i}\left(x,y\right) : = \left(i + 1\right)\log\left(1 + x\right) + \frac{1}{i+1} \log\left(1+y\right), && i = 0,1,2,3, \end{align} \tag{ 4 }$$

$$\begin{align} &u_{3,i}\left(x,y\right) : = x^{i+3} + y^{i+3}, && i = 0,1,2,3. \end{align} \tag{ 5 }$$

The corresponding source cochains f were computed taking the discrete Laplacian equation (A.9) of these functions. The full dataset is made of the (u, f) pairs. The values of the functions $u_{i,j}$ evaluated at the boundary of the square domain were used to enforce Dirichlet boundary conditions during the minimization of each candidate energy produced in the SR. As initial guesses for the minimization procedures, we took the null zero-cochain.

As reported in table 2, the set of hyperparameter values includes regularization (η = 0.1). With such a set, we obtained a recovery rate of 66% in the model discovery runs. Notice that the model discovery process is generally effective, even when the exact model is not found, as shown by the decrease of the mean of the fitness of the best individual of each generation, as the evolution proceeds (figure 2(A)).

Zoom In Zoom Out Reset image size

To assess the effect of regularization on model discovery, we performed 50 additional runs with the best combination of hyperparameters without regularization (η = 0) found in the grid search (see table C.1). In this case, the recovery rate dropped to 60%, while the length of the tree corresponding to the best individual found at the end of the SR procedure was $39.97\pm 15.53$, which is significantly larger than that (9) of the shortest string corresponding to the energy equation (2). To compensate for such a growth in model complexity, we set η = 0.1 while keeping all the other hyperparameters fixed and ran 50 additional evolutions seeding the initial population with the 50 individuals found at the end of the model discovery runs without regularization. After 20 generations, all runs converged (recovery rate = 100%) to equivalent expressions of the energy equation (2), with length comprised between 9 and 11 (fitness equal to 0.9 and 1.1, respectively), see table 3 and figure 3(A).

Zoom In Zoom Out Reset image size

Table 3. Discrete energies corresponding to the four best individuals obtained in the model discovery runs $(n = 50)$ with η = 0.1 and the initial population consisting of the best 50 individuals found in the model discovery with η = 0.

#	Energy	Training fit.	Test fit.	Test MSE
1	$\langle u, \delta(d(u/2) - f\,) \rangle$	0.9	0.9	$9.8 \cdot 10^{-10}$
2	$\langle u, \delta(du) - 2f \rangle$	0.9	0.9	$9.8 \cdot 10^{-10}$
3	$\langle \star(\star u), \delta(du) - 2f \rangle$	1.1	1.1	$9.8 \cdot 10^{-10}$
4	$\langle du, du \rangle - \langle\, f, 2u \rangle$	1.1	1.1	$9.8 \cdot 10^{-10}$

Unlike the previous problem, here we emulate a physical experiment involving an elastic cantilever rod subject to an end load by generating synthetic experimental data, instead of starting from an exact solution of a given discrete model. Specifically, we compute numerically the solution of a continuous model for different values of the load, we perturb it with noise and sample it at mesh nodes, to mimic experimental measurements. In the following, $\mathcal{I} = [0,L]\subset \mathbb{R}$ is a one-dimensional (1D) interval that represents the coordinates of the points of the longitudinal axis of the rod in its undeformed, straight configuration, L is the length of the rod, $\tilde{u}(\sigma)$ is the angle formed by the tangent to the axis of the rod and the horizontal direction at the dimensionless arc-length $\sigma\colon = s/L$. The axis of the rod is assumed to be inextensible and its curvature can be readily computed as $\tilde{u}^{\prime}$, where a prime denotes differentiation with respect to σ. For a constant cross-section rod subject to a vertical end-load at the right end and clamped at the left end, the boundary-value problem for the (dimensionless) continuous Euler's elastica equation reads [31]

$$\begin{align} \tilde{u}^{^{\prime\prime}} + f \cos \tilde{u} = 0 \quad \text{in } \left[0,1\right], \quad \tilde{u}\left(0\right) = 0, \quad \tilde{u}^{\prime}\left(1\right) = 0, \end{align} \tag{ 6 }$$

and the corresponding variational formulation is

$$\begin{equation} \min_{\tilde{u}}\ \frac{1}{2}\int_{0}^1 \left(\tilde{u}^{\prime}\right)^2\,\textrm{d}\sigma - \int_{0}^1 f \sin \tilde{u} \,\textrm{d}\sigma ,\quad \tilde{u} \left(0\right) = 0. \end{equation} \tag{ 7 }$$

In equations (6) and (7), $f = PL^2/B$ is the dimensionless load parameter, where P is the vertical component of the load applied at the right end and B is the bending stiffness of the rod.

In the discrete setting, we represent the unknown of the problem as a dual scalar-valued zero-cochain u, where the value at each dual node is the angle formed by the corresponding primal edge and the horizontal direction. The discrete curvature k may be evaluated by first taking the difference between the angles of two consecutive edges ($d^\star u$), then transferring this information to the primal nodes ($\star d^\star u$) and finally considering only the internal nodes ($\unicode{x1D7D9}_{\text{int}} \odot \star {d}^\star u$), where k is meaningful. Hence, the discrete equivalent of (7) reads

$$\begin{equation} \mathcal{E}_{\text{el}}\left(u\right) : = \frac{1}{2} \langle k, k \rangle - \langle f\unicode{x1D7D9}, \sin u \rangle, \end{equation} \tag{ 8 }$$

where $k = \unicode{x1D7D9}_{\text{int}} \odot \star {d}^\star u$. The 1D mesh for this problem was generated directly by code. Precisely, we generated a one-simplicial complex in $[0,1]$ consisting of 11 uniformly distributed nodes.

To generate the datasets, we considered a rod with length L = 1 m and bending stiffness $B \approx 7.854\, \textrm{N}\textrm{m}^2$. We computed numerically the solutions of the continuous problem equation (6) for ten different values of the force $P \in \{-5,-10,\ldots, -45, -50\}$ N and we sampled the (x, y) coordinates of the points of the deformed configurations at the nodes of the one-simplicial complex. Finally, we perturbed them with uniformly distributed noise (amplitude equal to ≈5% of the order of magnitude of the clean data) and we recovered the arrays of the coefficients of u. Linear approximations of these data were used as initial guesses for the minimization of the candidate energies found during the SR. The full dataset consists of the $(u,PL^2)$ pairs.

In a realistic experimental scenario, the load magnitude is prescribed, while the bending stiffness is an unknown material-dependent parameter that must be calibrated to reproduce the measurements. Hence, for each individual of the SR procedure the bending stiffness must be fitted on one solution $\bar{u}$ belonging to the training set. For this purpose, for each candidate energy $\mathcal{E}(u,f)$ we solved the constrained optimization problem

$$\begin{align} \min_{f \unicode{x2A7E} 0} \quad & ||u_f - \bar{u}||^2 \qquad \textrm{s.t.}\ u_f \in \arg \min_{\!\!\!\!\!\!\!\!\!\!\!\!u \text{: } u\left(0\right) = 0} \mathcal{E}\left(u,f\right). \end{align} \tag{ 9 }$$

Then, from the relation $B = PL^2/f$ we recovered the value of B and used it to evaluate the fitness of the individual over the whole training, validation and test sets, as described in section 2.2. In particular, taking the ground-truth model, $\mathcal{E} \equiv \mathcal{E}_{\text{el}}$, and solving the optimization problem we can identify a value $B_{\text{noise}} \approx 7.4062\, \textrm{N}\textrm{m}^2$ for the bending stiffness that accounts for the noise.

Both the average and the standard deviation of the fitness of the best individual of each generation over 50 model discovery runs decrease during the SR evolution (figure 4(A)) and attain the values 0.28 and 0.13, respectively, after 100 generations. Also, the selected models fit the data of the test set well (figures 4(B) and (C)). These observations confirm the effectiveness of the model discovery process along with the accuracy and generalization capabilities of the recovered models.

Zoom In Zoom Out Reset image size

We notice that the first and the last models among top four of table 4 are such that $B = B_{\text{noise}}/2h$, where h is the (primal) edge length. This relation exploits the uniform edge length of the simplicial complex representing the axis of the rod. Therefore, these models coincide with that of equation (8) up to addition and/or multiplication with a constant.

Table 4. Discrete energies corresponding to the four best individuals obtained in the model discovery runs $(n = 50)$ of the Elastica problem using the final values of the hyperparameters reported in table 2.

#	Energy	Training fit.	Test fit.	Test MSE	B
1	$\langle \star \unicode{x1D7D9}_{\text{int}}, (d^\star u)^2\rangle - \langle \sin u, f\unicode{x1D7D9}\rangle$	0.196	0.1939	0.0084	37.0312 Nm2
2	$\langle \arccos(-1) \sin u + \delta^\star(\sin(\sin(d^\star u))), u - f\unicode{x1D7D9} \rangle$	0.2123	0.2247	0.0075	7.1779 Nm2
3	$\langle u - f\unicode{x1D7D9}, \delta^\star (\sin(\sin(d^\star u))) + 1/\exp(-1)\sin u \rangle$	0.214	0.2168	0.0067	6.8262 Nm2
4	$\langle \star \unicode{x1D7D9}_{\text{int}}, (d^\star u)^2 \rangle - \langle f\unicode{x1D7D9}, \sin u \rangle + 1/2$	0.216	0.2139	0.0084	37.0312 Nm2

Differently from the previous benchmarks, this case deals with a vector-valued unknown and thus proves the capability of our framework of handling another class of problems. We collect the dimensionless coordinates of the nodes of the two-simplex representing the reference configuration of an elastic body in a primal vector-valued zero-cochain. We assume plane strain conditions. Let $\boldsymbol{e}_i, i = 1,2$ be the covariant basis for $\mathcal{V}$ made of two edge vectors of the reference configuration of a given two-simplex and let $\boldsymbol{e}^{\prime}_i$ be the covariant basis $\mathcal{V}$ made of two edge vectors of the corresponding current configuration of the two-simplex. The contravariant bases are indicated with superscript indices. We represent the two-dimensional (2D) deformation gradient F as a dual (tensor-valued) zero-cochain whose coefficient on a given two-simplex is

$$\begin{equation} \boldsymbol{F}\left(\star \sigma^2\right)\colon = {\boldsymbol{e}^{\prime}}_i \otimes {\boldsymbol{e}}^i. \end{equation} \tag{ 10 }$$

Then, we define the 2D, infinitesimal strain measure as a dual zero-cochain

$$\begin{equation} \boldsymbol{E}\colon = \frac{1}{2}\left(\boldsymbol{F} + \boldsymbol{F}^\mathrm{T}\right) - \boldsymbol{I}, \end{equation} \tag{ 11 }$$

and by applying the constitutive equation for linearly elastic isotropic bodies we obtain the corresponding dimensionless (2D) stress as a dual zero-cochain

$$\begin{equation} \boldsymbol{S}\colon = 2\boldsymbol{E} + \frac{\lambda\_}{\mu\_}\,\text{tr}\left(\boldsymbol{E}\right)\boldsymbol{I}, \end{equation} \tag{ 12 }$$

where I is the identity dual zero-cochain and $(\mu\_, \lambda\_)$ are the Lamé moduli (here, we take $\lambda\_ / \mu\_ = 10$). Notice that the second term in the constitutive equation involves the coefficient-wise multiplication between a scalar-valued and a tensor-valued cochain.

In the case of zero body forces and assuming Dirichlet boundary conditions (prescribed positions of the nodes), the potential energy function reads

$$\begin{equation} \mathcal{E}_\mathrm{LE}\colon = \frac{1}{2} \langle \boldsymbol{S},\boldsymbol{E} \rangle. \end{equation} \tag{ 13 }$$

The goal of this benchmark is recovering equation (13) as a function of F . The deformed configuration is given in terms of the current positions of the nodes that minimized equation (13). In principle, an admissible candidate energy $\mathcal{E}(\boldsymbol{F})$ must be frame-indifferent, i.e. $\mathcal{E}(\boldsymbol{F} + \boldsymbol{W}) = \mathcal{E}(\boldsymbol{F})$ for any skew-symmetric dual zero-cochain W . To enforce such a constraint, we added the penalty term $- (\mathcal{E}(\boldsymbol{F} + \tilde{\boldsymbol{W}}) - \mathcal{E}(\boldsymbol{F}))^2$ to the fitness equation (1), where

$$\begin{equation} \tilde{\boldsymbol{W}}\colon = \sum_{i = 0}^{n_2} \begin{bmatrix} 0 & e\\ -e & 0 \end{bmatrix} \star\!{\sigma^{2,i}}. \end{equation} \tag{ 14 }$$

This addition is clearly necessary but not sufficient. In practice, we observed that all the energies obtained at the end of the model discovery runs are frame-indifferent.

We considered a two-simplicial complex (dimensionless reference configuration) generated as a Delaunay triangulation of the unit square with 142 nodes. The full dataset (training, validation and test sets) consists of 20 samples, each including the positions of the nodes in the deformed configurations associated to the problems of uniaxial tension and pure shear. For these problems, the (uniform) strain and stress cochains are:

• Uniaxial tension:
  $$\begin{equation} \boldsymbol{S}_\mathrm{pt} : = \sum_{i = 0}^{n_2} \begin{bmatrix} s & 0\\ 0 & 0 \end{bmatrix} \star\!\sigma^{2,i}, \quad \boldsymbol{E}_\mathrm{pt} : = \sum_{i = 0}^{n_2} \begin{bmatrix} \epsilon & 0\\ 0 & \ell \end{bmatrix} \star\!\sigma^{2,i}, \end{equation} \tag{ 15 }$$
  where $s = 2\epsilon + (\lambda\_/\mu\_)(\epsilon + \ell)$, $\ell = - (\lambda\_/\mu\_)\,\epsilon/(2 + (\lambda\_/\mu\_))$, and $\epsilon \in \{0.01, 0.02, \ldots, 0.1\}$;
• Pure shear:
  $$\begin{equation} \boldsymbol{S}_\mathrm{ps} : = \sum_{i = 0}^{n_2} \begin{bmatrix} 0 & \gamma\\ \gamma & 0 \end{bmatrix} \star\!\sigma^{2,i}, \quad \boldsymbol{E}_\mathrm{ps} : = \sum_{i = 0}^{n_2} \frac{1}{2}\begin{bmatrix} 0 & \gamma\\ \gamma & 0 \end{bmatrix} \star\!\sigma^{2,i}, \end{equation} \tag{ 16 }$$
  where $\gamma \in \{0.01, 0.02, \ldots, 0.1\}$.

In both cases, the current positions of the nodes can be easily computed from the strain cochain and from equations (10) and (11). As initial guess for the energy minimization procedures associated to the evaluation of the fitnesses of the candidate energies, we took the reference positions of the nodes.

In the model discovery runs performed using the set of hyperparameters of table 2, we obtained a recovery rate of $92\%$. However, the SR process is always very effective, even in the runs where the correct solution is not attained (figure 5(A)). Notice that the second best individual of table 5 equals $1/2 \, \mathcal{E}_\mathrm{LE}$. For the remaining energies reported in table 5, we can still prove their equivalence to $\mathcal{E}_\mathrm{LE}$, under specific hypothesis. In the case of homogeneous deformations ($\boldsymbol{F}(\star \sigma_{2,i}) = \tilde{\boldsymbol{F}}$ for any i)

$$\begin{align} &\langle \boldsymbol{I}, \boldsymbol{F}^\mathrm{T} \rangle^2 = \left(\sum_{i = 0}^{n_2} |\sigma_{2,i}| \text{tr}\left(\boldsymbol{F}\left(\star \sigma_{2,i}\right)\right)\right)^2 = \text{tr}\left(\tilde{\boldsymbol{F}}\right)^2A^2, \end{align} \tag{ 17 }$$

$$\begin{align} &\langle\text{tr}\left(\boldsymbol{F}\right)\boldsymbol{I}, \text{sym}\left(\boldsymbol{F}\right)\rangle = \sum_{i = 0}^{n_2} |\sigma_{2,i}| \text{tr}\left(\boldsymbol{F}\left(\star \sigma_{2,i}\right)\right)^2 = \text{tr}\left(\tilde{\boldsymbol{F}}\right)^2A, \end{align} \tag{ 18 }$$

$$\begin{align} &\langle \boldsymbol{I}, \boldsymbol{I} \rangle = \sum_{i = 0}^{n_2} |\sigma_{2,i}| \text{tr}\left(\boldsymbol{I}\right) = 2A, \end{align} \tag{ 19 }$$

where A is the total area of the two-complex. Moreover, the following identities hold in general:

$$\begin{align} &\langle \boldsymbol{I} - \boldsymbol{F}, \text{sym}\left(\boldsymbol{F}\right)\rangle = \langle \boldsymbol{I} - \boldsymbol{F}^\mathrm{T}, \text{sym}\left(\boldsymbol{F}\right)\rangle = \langle \boldsymbol{I} - \text{sym}\left(\boldsymbol{F}\right), \text{sym}\left(\boldsymbol{F}\right)\rangle, \end{align} \tag{ 20 }$$

$$\begin{align} &\langle \boldsymbol{I}, \boldsymbol{F} \rangle = \langle \boldsymbol{I}, \boldsymbol{F}^\mathrm{T} \rangle = \langle \boldsymbol{I}, \text{sym}\left(\boldsymbol{F}\right) \rangle, \end{align} \tag{ 21 }$$

$$\begin{align} &\langle \text{tr}\left(\boldsymbol{F}\right)\boldsymbol{I}, \boldsymbol{I} \rangle = 2\langle \boldsymbol{I}, \boldsymbol{F} \rangle. \end{align} \tag{ 22 }$$

Additionally, when A = 1 as in our benchmark, equations (17)–(19) provide the identities $\langle \boldsymbol{I}, \boldsymbol{F}^\mathrm{T} \rangle^2 = \langle\text{tr}(\boldsymbol{F})\boldsymbol{I}, \text{sym}(\boldsymbol{F})\rangle$ and $\langle \boldsymbol{I}, \boldsymbol{I} \rangle = 2$. These facts, together with equations (20)–(22), allow to prove the aforementioned equivalences. Of course, some of these equivalences do not hold in general (i.e. when deformations are not homogeneous), so the recovered models correspond to different expressions of the potential energy, while having the same accuracy (in terms of MSE) with respect to the datasets. To resolve this ambiguity, we compute the associated energy densities by dividing the energies by A, since the deformations are homogeneous. We find that energies #1, #3 and #4 of table 5 are not admissible, because their energy densities linearly depend on the area A due to the term in equation (17). However, these expressions may be fixed by interpreting such a term as that in equation (18), so that they all correspond to the exact energy density that was used to generate the dataset. We remark that this a posteriori discussion can be done since our method produces interpretable models that are amenable to analysis.

Zoom In Zoom Out Reset image size

Table 5. Discrete energies corresponding to the four best individuals obtained from the model discovery runs $(n = 50)$ of the Linear elasticity problem using the final values of the hyperparameters reported in table 2.

#	Energy	Training MSE	Test MSE
1	$\langle \boldsymbol{I} - \boldsymbol{F}^\mathrm{T}, \boldsymbol{I} \rangle^2 - 0.1(\langle \boldsymbol{I} - \boldsymbol{F}^\mathrm{T}, \text{sym}(\boldsymbol{F}) - \boldsymbol{I} \rangle + \langle \boldsymbol{I} - \text{sym}(\boldsymbol{F}), \text{sym}(\boldsymbol{F}) - \boldsymbol{I} \rangle)$	0	$1.1672 \cdot 10^{-16}$
2	$\langle -0.5\boldsymbol{I} + 0.5 \,\text{sym}(\boldsymbol{F}), \text{sym}(\boldsymbol{F}) - \boldsymbol{I} + (\text{tr}(\boldsymbol{F}) - \text{tr}(\boldsymbol{I}))(2\text{tr}(\boldsymbol{I})\boldsymbol{I} + \boldsymbol{I}) \rangle$	0	$2.0785 \cdot 10^{-16}$
3	$\langle \boldsymbol{I} - \boldsymbol{F}^\mathrm{T}, \boldsymbol{I} - \text{sym}(\boldsymbol{F}) + 5\langle \boldsymbol{I} - \boldsymbol{F}, \boldsymbol{I} \rangle \boldsymbol{I} \rangle$	0	$3.7309 \cdot 10^{-16}$
4	$\langle \text{sym}(\boldsymbol{F}) - \boldsymbol{I}, \text{sym}(\boldsymbol{F}) - \boldsymbol{I} \rangle + 5 \langle \boldsymbol{I} - \boldsymbol{F}, \boldsymbol{I} \rangle^2$	0	$8.7033 \cdot 10^{-16}$

We notice that traditional SR can be used when the deformation gradient and the stress are homogeneous. However, its performance is much worse than that of our method, because the former lacks tensor-based operations, which are needed to achieve a compact representation of the energy density. The results of the comparison between the two methods are reported in appendix D, along with an example of a genuine field problem in linear elasticity, which clearly cannot be solved via a traditional SR approach. We believe that these simple, yet effective examples demonstrate the advantage of the proposed approach in situations where conventional SR methods struggle to learn an accurate model.

We have developed our framework based on a special class of problems (variational and steady-state). This choice was motivated by the simpler implementation of the numerical routines compared to time-dependent, non-variational problems. To overcome this limitation, we will implement time-integration schemes and, regarding non-variational problems, we will use trees to represent the residual of the system of governing equations, instead of the energy functional.

A p-simplex is the convex hull of p + 1 independent points in $\mathcal{E}^N$. A collection of simplices defines a simplicial complex K if it satisfies the following conditions:

• every lower-order simplex of a given simplex in K is in K;
• the intersection of any two simplices is either empty or a boundary element of both simplices.

A (scalar-valued) p-chain is a collection of values, one for each of the p-simplices of K. Each p-simplex corresponds to a basic chain $\sigma_{p,i}$ (i.e. a chain assigning 1 to a given simplex and 0 to the others) and we can use the same notation to indicate both the simplex and the corresponding basic chain. The p-chains group is denoted by $\mathcal{C}_p(K, \mathbb{R})$. A p-chain c_p can be written in an unique way as a finite sum of basic chains:

$$\begin{align} c_p = \sum_{i = 1}^{n_p} c_p\left(\sigma_{p,i}\right) \sigma_{p,i}, \end{align} \tag{ A.1 }$$

where n_p is the number of p-simplices of K, $c_p(\sigma_{p,i}) \in \mathbb{R}$, and n_p is the number of p-simplices. In computations, a chain can be represented by a vector where the ith entry is the value of the chain on the ith simplex.

A p-cochain is a linear functional that maps chains to scalars. Analogously to chains, we can represent a cochain c^p in terms of basic cochains

$$\begin{align} c^p = \sum_{i = 1}^{n_p} c^p\left(\sigma^{p,i}\right) \sigma^{p,i}, \end{align} \tag{ A.2 }$$

where $c^p(\sigma^{p,i})$ are coefficients to be chosen in an appropriate space and $\sigma^{p,i}$ is the ith basic p-cochain, i.e. the cochain that assigns 1 to the ith basic chain, which we can again identify with the ith simplex, and 0 to the others. The space of scalar-valued p-cochains on a simplicial complex K is denoted by $\mathcal{C}^p(K, \mathbb{R})$, while the spaces of vector-valued and tensor-valued cochains are denoted by $\mathcal{C}^p(K, \mathcal{V})$ and $\mathcal{C}^p(K, \text{Lin})$, respectively. We indicate any of these spaces with $\mathcal{C}^p(K, \cdot)$. We will use bold symbols for vector-valued and tensor-valued cochains and for their coefficients $\boldsymbol{c}^p(\sigma^p_i) \in \mathcal{V}$ and $\boldsymbol{C}^p(\sigma^p_i) \in \text{Lin}$, respectively. In computations, a cochain can be represented by a multi-dimensionl array, where each row contains its value(s) associated to a given simplex.

A p-cochain can be readily evaluated on a p-chain by exploiting the duality between basic chains and basic cochains. We denote such an evaluation by the pairing

$$\begin{align} \langle \langle c^p, c_p \rangle \rangle : = \sum_{i = 1}^{n_p} c^p\left(\sigma^{p,i}\right) c_p\left(\sigma_{p,i}\right). \end{align} \tag{ A.3 }$$

with $\langle \langle \sigma^{p,i}, \sigma_{p,j} \rangle \rangle = \delta^i_j$. Notice that $\langle \langle c^p, \sigma_{p,i} \rangle \rangle = c^p(\sigma^{p,i})$ and $\langle \langle \sigma^{p,i}, c_p \rangle \rangle = c_p(\sigma_{p,i})$. This operation is a discrete preliminary to integration when regarding cochains as densities with respect to the measures imparted to cells by real-valued chains [16].

The boundary of a p-simplex is the union of its $(p-1)$-faces, that is, $(p-1)$-simplices consisting of subsets of the vertices of the p-simplex. Formally, the boundary operator is a linear map $\partial: \mathcal{C}_p(K, \mathbb{R}) \rightarrow \mathcal{C}_{p-1}(K, \mathbb{R})$ defined by

$$\begin{align} \partial \sigma_p = \sum_{i = 0}^p \left(-1\right)^i \left[\sigma_{0,0}, \ldots, \hat{\sigma}_{0,i}, \ldots, \sigma_{0,p}\right], \end{align} \tag{ A.4 }$$

where $\sigma_p = [\sigma_{0,0}, \ldots, \sigma_{0,p}]$ is a generic p-simplex and $[\sigma_{0,0}, \ldots, \hat{\sigma}_{0,i}, \ldots, \sigma_{0,p}]$ is the $(p-1)$-simplex obtained omitting the ith vertex. The boundary operator encodes fundamental topological information about the simplicial complex and allows to define the coboundary operator $d: \mathcal{C}^{p-1}(K, \mathbb{R}) \rightarrow \mathcal{C}^p(K, \mathbb{R})$ as its adjoint with respect to the pairing defined in equation (A.3):

$$\begin{align} \langle \langle d c^{p-1}, c_p \rangle \rangle = \langle \langle c^{p-1}, \partial c_p \rangle \rangle. \end{align} \tag{ A.5 }$$

The boundary operator $\partial$ may be represented by the incidence matrix [12, 20], while the coboundary operator corresponds to its transpose. Notice that there is one (co)boundary matrix per simplex dimension.

Given a simplicial complex K of dimension n, we can define its dual $\star K$ by associating each p-simplex $\sigma_p \in K$ with its dual $(n-p)$-simplex $\star \sigma_{(n-p)}$, following a specific rule (see Remark 2.5.1 in [11]). In this work, we use the circumcentric (or Voronoi) dual. Once the orientation of the primal is chosen, that of the dual is induced, and the dual boundary and coboundary $d^\star$ account for such orientation (see [10, 11, 36] for more details).

The diagonal discrete Hodge star $\star: \mathcal{C}^p(K, \cdot) \rightarrow \mathcal{C}^{n-p}(\star K, \cdot)$ is defined as follows [10, 11]:

$$\begin{align} \frac{1}{|\star\!\sigma_p|}\star\!c^p\left(\star \sigma^p\right)\colon = \frac{1}{|\sigma_p|} c^p\left(\sigma^p\right). \end{align} \tag{ A.6 }$$

Here, $\star \sigma^p$ is the basic dual $(n-p)$-cochain corresponding to σ^p .

We define the L² inner product on cochains as

$$\begin{align} \langle c_1^p, c_2^{p} \rangle\colon = \sum_{i = 1}^{n_p} c_1^p\left(\sigma^{p,i}\right) \cdot \star c_2^{p}\left(\star \sigma^{p,i}\right), \end{align} \tag{ A.7 }$$

with $c_1^p, c^p_2 \in \mathcal{C}^p(K, \cdot)$, where the dot denotes the standard multiplication between scalars or the standard Euclidean inner product between vectors or tensors.

Further, with respect to this inner product, we may introduce the adjoint of the coboundary operator, the discrete codifferential operator δ, which maps primal scalar-valued p-cochains into primal scalar-valued $(p-1)$-cochains, as

$$\begin{align} \delta\colon = \left(-1\right)^{n\left(p-1\right) + 1} \star d^\star \star, \end{align} \tag{ A.8 }$$

by requiring that $\langle d\alpha, \beta \rangle = \langle \alpha, \delta \beta \rangle$, for any $\alpha \in \mathcal{C}^{k-1}(K, \mathbb{R})$, $\beta \in \mathcal{C}^{k}(K, \mathbb{R})$. In the previous statement, we write $\star$ for both the Hodge star and its inverse and $d^\star$ is the dual coboundary operator. The dual inner product and dual codifferential $\delta^\star$ can be analogously defined.

In order to deal with the Poisson equation (see section 3.1), we introduce the discrete Laplace–de Rham operator:

$$\begin{align} \Delta\colon = d\delta + \delta d. \end{align} \tag{ A.9 }$$

This is the discrete analogue of the Laplace–de Rham operator on differential forms, which coincides with minus the standard Laplacian of a scalar field [37].

Finally, we extend functions operating on scalar fields (such as $\sin, \cos, \exp, \log$) to cochains by component-wise application. In particular, the component-wise multiplication between scalar-valued cochains will be denoted by the symbol $\odot$.