Elsevier

Computers & Graphics

Volume 37, Issue 6, October 2013, Pages 687-696
Computers & Graphics

SMI 2013
Steepest descent paths on simplicial meshes of arbitrary dimensions

https://doi.org/10.1016/j.cag.2013.05.003Get rights and content

Highlights

  • Calculation of steepest descent paths on simplicial meshes.

  • Suitable for PL functions on arbitrary-dimensional domains.

  • The elementary triangular patch is modeled using a modified side-vertex algorithm.

  • Generalization of existing algorithms for two and three dimensions.

Abstract

This paper introduces an algorithm to compute steepest descent paths on multivariate piecewise-linear functions on Euclidean domains of arbitrary dimensions and topology. The domain of the function is required to be a finite PL-manifold modeled by a simplicial complex. Given a starting point in such a domain, the resulting steepest descent path is represented by a sequence of segments terminating at a local minimum. Existing approaches for two and three dimensions define few ad hoc procedures to calculate these segments within simplices of dimensions one, two and three. Unfortunately, in a dimension-independent setting this case-by-case approach is no longer applicable, and a generalized theory and a corresponding algorithm must be designed. In this paper, the calculation is based on the derivation of the analytical form of the hyperplane containing the simplex, independent of its dimension. Our prototype implementation demonstrates that the algorithm is efficient even for significantly complex domains.

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 f:RnR is a real function and f(x)=(f/x1,f/x2,,f/xn) is its gradient, then the steepest descent direction for x0Rn is given by the vector d=f(x0). When starting from an initial point x0, the steepest descent method computes a new point x1x0αf(x0). Then, in a second step it computes another point x2x1αf(x1), and so on. In general, the algorithm computes a sequence of points x0xk 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 pi is the position of the ith vertex of a closed and manifold triangle mesh M, then the gradient of the function f:MR can be estimated as [12]f(pi)jN(i)[f(pj)f(pi)]wj,where the wjs are proper weights that do not depend on f, and N(i) is the set of vertices connected to i by an edge.

In contrast, the gradient f can be computed exactly for points within triangles. Furthermore, f does not change across the points within a given triangle. So, on a triangle t with vertices (pi,pj,pk) and unit normal N, the gradient is the solution f|t of the 3×3 linear system [13](pjpipkpjN)f|t=(f(pj)f(pi)f(pk)f(pj)0)

Therefore, alternatively to Eq. (1), the steepest direction at the vertex position pi 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 p of the path: (1) p is in the interior of a triangle—the steepest direction is unique and orthogonal to the level lines (e.g. Eq. (2); (2) p is in the interior of an edge—one or two locally steepest directions are possible; and (3) p 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 p0 within a (maximal) n-dimensional simplex σ0, computes the direction of steepest descent 0 and shoots a ray from p0 in the direction of 0. The ray intersects an n1-dimensional simplex which is a face of both σ0 and an opposite maximal simplex σ1. Let p1 be the point of intersection and 1 be the direction of steepest descent calculated within σ1. Once again, the algorithm shoots a ray from p1 in the direction of 1,

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 |K|r. A point p|K|r can either be the image of a vertex or belong to the image of a higher-dimensional simplex. We first consider the case where p 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 (n+1)-uples 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)

  • Pascucci V. Topology diagrams of scalar fields in scientific visualization. In: Topological data structures for...
  • Cited by (6)

    • Protonic effects on the melting behavior of magmas: The Forsterite-Diopside-Silica system at high pressure as an example

      2021, Chemical Geology
      Citation 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 Structures
      Citation 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 Geology
      Citation 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
    View full text