SMI 2013Steepest descent paths on simplicial meshes of arbitrary dimensions
Graphical abstract
Introduction
Being able to efficiently calculate the local minima of a real-valued function is of clear importance for numerous application contexts. In some domains, such as physics or chemistry, it is not rare to deal with problems where the function F to be minimized represents the status of a system, and its local minima represent particular configurations for which the system is stable or in equilibrium. It is thus interesting to know how such a system evolves toward a stable state when left unconstrained from an initial configuration. Such an evolution can be represented within the domain of F by a sequence of configurations that form a path connecting the initial state with the final stable state. Normally the system cannot jump from a configuration to another without continuity (though there are exceptions, e.g. in quantum mechanics) and thus the path is continuous as well.
For several natural phenomena, a given non-stable state tends to change toward equilibrium with maximum speed. Such a phenomenon can be modeled by a real-valued function F that associates an energy with each state of the system and, if F is differentiable, the direction of the maximum speed is the direction of the negative gradient of F. When the function is expressed through a closed analytical form, computing the gradient becomes a matter of symbolic calculus and thus is relatively easy. In contrast, when the function is described through a finite set of samples, one may need to employ more sophisticated techniques.
In this paper, we review the methods proposed so far to compute steepest descent paths on piecewise linear surface meshes, and generalize them within a unified approach that works on simplicial complexes of any dimension. Furthermore, we provide a novel solution to resolve the cases of indifferent equilibrium that may arise when tracking the steepest descent path. Finally, we show how our solution can be exploited in the field of Material Sciences to study the crystallization of molten substances.
Note that generalizing existing algorithms is not as easy as it could appear at a first sight. On triangle meshes, for example, there are three possibilities: (1) the path passes through triangles, or (2) it may run along edges or (3) it may cross vertices. Thus, algorithms treating these models can enumerate a few possibilities and handle each of them as a particular case. Furthermore, the treatment of special cases of indifferent equilibrium can also be handled by enumerating the possibilities. Conversely, in a generalized approach such as the one introduced in this paper, a dimension-independent solution must be provided which can treat all the cases based on a common approach. To the best of our knowledge no solution to this problem has been proposed so far, while there are application contexts (e.g. Material Sciences) that call for a solution.
The steepest descent method, also known as gradient descent method, is an algorithm to compute local minima of real-valued functions of n variables. While running, it constructs a piecewise-linear approximation of the steepest descent path connecting an initial point to a local minimum.
Formally, if is a real function and is its gradient, then the steepest descent direction for is given by the vector . When starting from an initial point , the steepest descent method computes a new point . Then, in a second step it computes another point , and so on. In general, the algorithm computes a sequence of points that is expected to converge to a local minimum. The positive parameter might be kept constant or can change from step to step. Small values of prevent the algorithm to oscillate around minima, but too small values might significantly slow down the entire process.
Near minima, the steepest descent method tends to go slowly and in some cases to have an erratic path. For these reasons, other minimization algorithms have been proposed: notable examples are the conjugate gradient that uses conjugate directions instead of the local gradient to take into account the previously chosen directions, and the Newton–Raphson that exploits information of the second derivative to locate the minimum of the function.
With small modifications, the steepest descent method can be adapted to the case of piecewise-linear meshes [1] approximating differentiable functions. In this case, there is no need to choose the parameter because each segment of the path is bounded by the size of mesh elements. As a consequence, the aforementioned drawbacks are no longer an issue. Unfortunately, though there are some methods for calculation on two- [1] and three- [2], [3] dimensional meshes, a complete adaptation is missing for the cases where differentiable functions are defined on n-dimensional simplicial domains.
When the domain is a two-dimensional simplicial mesh, several methods exist to calculate steepest descent paths. A big family of algorithms uses descent (and ascent) paths to bound the cells of the so-called Morse–Smale complexes [1], [2], [3]. Within this family, the methods described in [4], [5], [6] force the steepest descent paths to be sequences of mesh edges. Specifically, these approaches start at a given vertex, and at each step simply choose the neighboring vertex associated with the smallest function value. Though these algorithms are extremely efficient, they produce only approximate solutions. Conversely, in [7], [8] the paths are computed based on [1] and are free to cross the triangles, therefore produce much more precise results. These algorithms have applications in both Computer Graphics [9], [10] and Visualization [11].
Clearly, if the function is piecewise-linearly defined on a simplicial mesh, the gradient is generally undefined for points on edges and for the vertices, so it may be necessary to provide an estimate for these particular cases. If is the position of the ith vertex of a closed and manifold triangle mesh , then the gradient of the function can be estimated as [12]where the are proper weights that do not depend on f, and is the set of vertices connected to i by an edge.
In contrast, the gradient can be computed exactly for points within triangles. Furthermore, does not change across the points within a given triangle. So, on a triangle t with vertices and unit normal , the gradient is the solution of the 3×3 linear system [13]
Therefore, alternatively to Eq. (1), the steepest direction at the vertex position can be estimated by considering its incident triangle's gradient having largest magnitude.
The most detailed description of an algorithm that computes the steepest paths is given for the two-dimensional case in [1]. Therein, the path can go along edges or pass through faces of the triangulation. Three cases occur for a starting point of the path: (1) is in the interior of a triangle—the steepest direction is unique and orthogonal to the level lines (e.g. Eq. (2); (2) is in the interior of an edge—one or two locally steepest directions are possible; and (3) is a vertex—there may be as many locally steepest directions as the number of its incident triangles.
The case of three-dimensional manifold domains is treated in [2], [3], where smooth manifolds are approximated with piecewise linear simplicial complexes that linearly interpolate samples of the function f. Both in [3] and in [2] ascending/descending arcs are defined as piecewise linear curves (1-manifold) between saddles and maxima/minima lying along edges of the input mesh.
Section snippets
Definitions and problem statement
In this section basic definitions are introduced to support a formal description of the algorithm. The following definitions are adapted from [14], [15].
Algorithm description
In the easiest case, the algorithm starts with a given point within a (maximal) n-dimensional simplex , computes the direction of steepest descent and shoots a ray from in the direction of . The ray intersects an simplex which is a face of both and an opposite maximal simplex . Let be the point of intersection and be the direction of steepest descent calculated within . Once again, the algorithm shoots a ray from in the direction of ,
Gradient evaluation
In order to implement the proposed algorithm, we need a procedure to derive the gradient vector of f at each point in its domain . A point can either be the image of a vertex or belong to the image of a higher-dimensional simplex. We first consider the case where belongs to the image of a maximal simplex (i.e. if K is n-dimensional, a simplex is maximal if its dimension is n). Afterwards, we generalize the result to lower-dimensional simplices through dimensionality reduction.
Implementation and results
Since the algorithm is dimension-independent, in our implementation we have employed an extended indexed data structure with adjacencies [16] enriched with function values at vertices. Namely, if the mesh is n-dimensional and made of NV vertices and NS maximal simplices, our data structure explicitly encodes:
- •
an array G of NV n-uples of coordinates representing the vertex positions along with an array F of NV corresponding function values;
- •
an array S of NS of indexes in G representing
Conclusions and future work
This research shows that it is possible to efficiently compute steepest descent paths on piecewise-linear functions defined on Euclidean domains of any dimension. In practice, such functions are defined through a finite sampling and the resulting discretization may lead to flat simplices for which we have proposed an appropriate treatment. In some cases, a vector field may be available instead of the function, and often such a field determines the path's direction at all the points; in these
Acknowledgments
This work was partly supported by the Italian MIUR-PRIN Project N. 2009B3SAFK (Topology of phase diagrams and lines of descent) and by the Petromaks program of the Norwegian Research Council through the Geoillustrator project (N. 200512). Thanks are due to the SMG members at IMATI for helpful discussions.
References (23)
- et al.
Topology preserving data simplification with error bounds
Comput Graph
(1998) - et al.
Shape approximation by differential properties of scalar functions
Comput Graph
(2010) - et al.
On converting sets of tetrahedra to combinatorial and PL manifolds
Comput-Aided Geom. Des
(2009) - et al.
Computing smooth approximations of scalar functions with constraints
Comput Graph
(2009) - Edelsbrunner H, Harer J, Zomorodian A. Hierarchical morse complexes for piecewise linear 2-manifolds. In: SCG '01:...
- Edelsbrunner H, Harer J, Natarajan V, Pascucci V. Morse–Smale complexes for piecewise linear 3-manifolds. In: SCG '03:...
- Natarajan V, Pascucci V. Volumetric data analysis using Morse–Smale complexes. In: SMI '05: proceedings of the...
- et al.
Algorithms for extracting correct critical points and constructing topological graphs from discrete geographic elevation data
Comput Graph Forum
(1995) - Bremer P, Pascucci V. A practical approach to two-dimensional scalar topology. In: Topology-based methods in...
- Bremer P-T, Edelsbrunner H, Hamann B, Pascucci V. A multi-resolution data structure for two-dimensional morse...
Cited by (6)
Protonic effects on the melting behavior of magmas: The Forsterite-Diopside-Silica system at high pressure as an example
2021, Chemical GeologyCitation Excerpt :The thermodynamic and thermophysical properties of molten silica and molten alumina were obtained ab-initio (Ottonello et al., 2010a; Belmonte et al., 2013) as well as those of some high-P polymorphs of crucial interest (Ottonello et al., 2009a, 2009b, 2010b; Belmonte et al., 2014). Because manifold simplicial complexes are the most effective models to represent the topological complexities of multiphase systems (Natali et al., 2011), the loci of chemical potential equality were obtained by a simplicial procedure and the following convex-hull analysis of the segmented de Launay space (SIMPLICIAL code; Attene and Ottonello, 2011; Natali et al., 2013). We may see that the actual topological complexity of the simplex exceeds by far the simplified realizations which may be found in literature.
Improved genetic algorithm for optimization design of a three-dimensional braided composite joint
2018, Composite StructuresCitation Excerpt :Meanwhile, the effect of GAs is significantly influenced by the population size [32], but the population scale in the present study is restricted by the finite computational resources. Thus, an SDM [33–35] (though a more appropriate name here might be the steepest ascend method) is employed to improve the precision of the solution. In the GA, an excellent individual might be buried during the evolution.
The system MgO-Al<inf>2</inf>O<inf>3</inf>-SiO<inf>2</inf> under pressure: A computational study of melting relations and phase diagrams
2017, Chemical GeologyCitation Excerpt :Phase diagram topology is computed in a wide range of P-T conditions by a Gibbs free energy minimisation algorithm which employs the convex-hull method to analyze equipotential surfaces in multi-component systems. The convex-hull technique is described elsewhere (e.g. Attene and Ottonello, 2011; Natali et al., 2013; Ottonello et al., 2013), while some examples on the use of this method to calculate ternary phase diagrams can be found in the literature (e.g. Connolly and Kerrick, 1987; Voskov and Voronin, 2010; Voskov et al., 2015). Heat capacity (CP), standard-state entropy (S0298.15) and enthalpy of formation from the elements (H0f,298.15) of the pure liquid end-member compositions (i.e. MgO, Al2O3 and SiO2) have been obtained from ab initio structure-energy and vibrational calculations performed with the Polarized Continuum Model (PCM) (Tomasi and Persico, 1994; Ottonello et al., 2010a; Belmonte et al., 2013; see also Supplementary material for details).
Ab initio reactivity of Earth's materials
2018, Rivista del Nuovo Cimento