Searching the solution landscape by generalized high-index saddle dynamics

Given an autonomous dynamical system [44]

where $\bm{F}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a $\mathcal{C}^{r}(r\geqslant 2)$ function, a point ${\bm{x}}^{*}\in\mathbb{R}^{n}$ is called a stationary point (or equilibrium solution) of Eq. 2 if $\bm{F}(\bm{x}^{*})=\bm{0}$. Let $J(\bm{x})=\nabla\bm{F}(\bm{x})$ denote the Jacobian of $\bm{F}(\bm{x})$, and $\langle\cdot,\cdot\rangle$ denote an inner product on $\mathbb{R}^{n}$, i.e. $\langle\bm{x},\bm{y}\rangle=\bm{y}^{\top}\bm{x}$.

For a stationary point $\bm{x}^{\ast}$, taking $\bm{x}=\bm{x}^{*}+\bm{y}$ in Eq. 2, we have

where $\|\cdot\|$ denotes the norm induced by the inner product. The associated linear system

is used to determine the stability of $\bm{x}^{*}$. Let $\{\bm{w}_{1},\cdots,\bm{w}_{k_{\mathrm{u}}}\}$, $\{\bm{w}_{k_{\mathrm{u}}+1},\cdots,\bm{w}_{k_{\mathrm{u}}+k_{\mathrm{s}}}\}$ and $\{\bm{w}_{k_{\mathrm{u}}+k_{\mathrm{s}}+1},\cdots,\bm{w}_{k_{\mathrm{u}}+k_{\mathrm{s}}+k_{\mathrm{c}}}\}\subset\mathbb{C}^{n}$ denote the (generalized) right eigenvectors of $J(\bm{x})$ corresponding to the eigenvalues of $J(\bm{x})$ with positive, negative and zero real parts, respectively, where $k_{\mathrm{u}}+k_{\mathrm{s}}+k_{\mathrm{c}}=n$. With these eigenvectors, we define unstable, stable and center subspaces of the Jacobian $J(\bm{x})$, respectively, as $\mathcal{W}^{\mathrm{u}}(\bm{x})$, $\mathcal{W}^{\mathrm{s}}(\bm{x})$ and $\mathcal{W}^{\mathrm{c}}(\bm{x})$:

From the primary decomposition theorem [27], $\mathbb{R}^{n}$ can be decomposed as a direct sum:

By treating $J(\bm{x}^{\ast})$ as a linear operator on $\mathbb{R}^{n}$, $\mathcal{W}^{\mathrm{u}}(\bm{x}^{*})$, $\mathcal{W}^{\mathrm{s}}(\bm{x}^{*})$ and $\mathcal{W}^{\mathrm{c}}(\bm{x}^{*})$ are the invariant subspaces of the linear system Eq. 4.

The nonlinear system Eq. 3 possesses a $k_{\mathrm{u}}$-dimensional local invariant unstable manifold $\mathcal{M}^{\mathrm{u}}_{\mathrm{loc}}(\bm{x}^{*})$, a $k_{\mathrm{s}}$-dimensional local invariant stable manifold $\mathcal{M}^{\mathrm{s}}_{\mathrm{loc}}(\bm{x}^{*})$, and a $k_{\mathrm{c}}$-dimensional local invariant center manifold $\mathcal{M}^{\mathrm{c}}_{\mathrm{loc}}(\bm{x}^{*})$ near the origin. They are intersecting and tangent to the respective invariant subspaces $\mathcal{W}^{\mathrm{u}}(\bm{x}^{*})$, $\mathcal{W}^{\mathrm{s}}(\bm{x}^{*})$, $\mathcal{W}^{\mathrm{c}}(\bm{x}^{*})$ [44]. “Local” here refers to that these manifolds with boundaries are defined only in a neighborhood of $\bm{x}^{*}$. The invariant subspaces of the linear system Eq. 4 can be used to study the stationary point of the nonlinear system Eq. 3.

A stationary point $\bm{x}^{*}$ of Eq. 2 is called hyperbolic if no eigenvalues of $J(x^{*})$ have zero real part, which means $\mathcal{W}^{\mathrm{c}}(\bm{x}^{*})=\{\bm{0}\}$ and the decomposition Eq. 6 can be rewritten as

A hyperbolic stationary point is called a saddle if $\mathcal{W}^{\mathrm{u}}(\bm{x}^{*})$ and $\mathcal{W}^{\mathrm{s}}(\bm{x}^{*})$ are nontrivial. The hyperbolic stationary point $x^{*}$ is called a sink (source) if all the eigenvalues of $J(\bm{x}^{*})$ have negative (positive) real parts. The index of a stationary point $\bm{x}^{\ast}$ is defined as the dimension of the unstable subspace $\mathcal{W}^{\mathrm{u}}(\bm{x}^{*})$ [44].

The HiOSD method is designed for finding index-$k$ saddles of an energy function $E(\bm{x})$ [46]. For gradient systems $\bm{F}(\bm{x})=-\nabla E(\bm{x})$, the Jacobian $J(\bm{x})=-\nabla^{2}E(\bm{x})=-G(\bm{x})$ is self-adjoint with all eigenvalues real, where $G(\bm{x})$ denotes the Hessian of $E(\bm{x})$. The hyperbolic $k$-saddle $\bm{x}^{*}$ is a local maximum on the linear manifold $\bm{x}^{*}+\mathcal{W}^{\mathrm{u}}(\bm{x}^{*})$ and a local minimum on $\bm{x}^{*}+\mathcal{W}^{\mathrm{s}}(\bm{x}^{*})$.

Regardless of the dimer approximation temporarily, the high-index saddle dynamics (HiSD) for a $k$-saddle ($k$-HiSD) has the following form:

which involves a position variable $\bm{x}$ and $k$ direction variables $\bm{v}_{i}$. The dynamics Eq. 8 is coupled with an initial condition:

The dynamics Eq. 8a represents a transformed gradient flow:

where $\mathcal{P}_{\mathcal{V}}$ denotes the orthogonal projection operator on a finite-dimensional subspace $\mathcal{V}$. Here, $-\mathcal{P}_{\mathcal{W}^{\mathrm{u}}(\bm{x})}\bm{F}(\bm{x})$ is taken as an ascent direction on the subspace $\mathcal{W}^{\mathrm{u}}(\bm{x})$ and $\bm{F}(\bm{x})-\mathcal{P}_{\mathcal{W}^{\mathrm{u}}(\bm{x})}\bm{F}(\bm{x})$ is a descent direction on the subspace $\mathcal{W}^{\mathrm{s}}(\bm{x})$. The dynamics Eq. 8b finds an orthonormal basis of $\mathcal{W}^{\mathrm{u}}(\bm{x})$. Thanks to the self-adjoint Jacobian $J(\bm{x})$, we can simply take $\bm{v}_{i}$ as a unit eigenvector corresponding to the $i$th smallest eigenvalue of $G(\bm{x})$. Therefore, the $i$th eigenvector $\bm{v}_{i}$ can be obtained by a constrained optimization problem with the knowledge of $\bm{v}_{1},\cdots,\bm{v}_{i-1}$:

Then the $k$ Rayleigh quotients Eq. 11 are minimized simultaneously using the gradient flow of $\bm{v}_{i}$:

as the dynamics Eq. 8b.

The GHiSD for a $k$-saddle ($k$-GHiSD) of the dynamical system Eq. 2 has the following form:

with an initial condition Eq. 9.

The dynamics Eq. 13a is the same as the dynamics Eq. 8a of HiSD. For non-gradient systems, the Jacobian $J(\bm{x})$ is not self-adjoint, and the eigenvectors may not be orthogonal. The loss of orthogonality between $\mathcal{W}^{\mathrm{u}}(\bm{x})$ and $\mathcal{W}^{\mathrm{s}}(\bm{x})$ indicates that it may be hard to decompose $\bm{F}(\bm{x})$ according to the direct sum Eq. 7. For this reason, the gentlest ascent dynamics (GAD) for searching $1$-saddles in non-gradient systems requires two direction variables to approximate the left and right eigenvectors corresponding to the largest eigenvalue of $J(\bm{x})$ [16]. Nevertheless, Gu and Zhou later proposed a simplified GAD in non-gradient systems that uses only one direction variable and proved that its stable critical point corresponds to a $1$-saddle of the original system [20]. Thus, we adopt the similar idea that the dynamics Eq. 10 could be used to find high-index saddles of non-gradient systems and the HiSD Eq. 8a is preserved in GHiSD Eq. 13a, as long as $\bm{v}_{1},\cdots,\bm{v}_{k}$ approximate an orthonormal basis of the subspace $\mathcal{W}^{\mathrm{u}}(\bm{x})$.

The dynamics Eq. 13b aims to find an orthonormal basis $\bm{v}_{1},\cdots,\bm{v}_{k}$ of the subspace $\mathcal{W}^{\mathrm{u}}(\bm{x})$. The dynamical flow

should preserve the orthonormal condition

which indicates

Thus, the dynamics Eq. 13b is derived from Eq. 14 and Eq. 16, which also accords with Eq. 8b under the condition $J(\bm{x})=J^{\top}(\bm{x})$.

Next, we show the following theorem for linear stability for the $k$-GHiSD Eq. 13 under some assumptions.

Remark: In the above theorem, all eigenvalues of $J(\bm{x}^{*})$ are assumed to be real to simplify the derivations. These assumptions can be relaxed as the eigenvalues $\lambda_{k+1},\cdots,\lambda_{n}$ have nonzero real parts. Specifically, assumed that $\lambda_{1}>\cdots>\lambda_{k}>\mathrm{Re}\lambda_{k+1}\geqslant\cdots\geqslant\mathrm{Re}\lambda_{n}$, are nonzero, the proof can be accomplished similarly.

Assumed that the eigenvalues $\lambda_{1},\cdots,\lambda_{k}$ are real, the only difference between GHiSD Eq. 13 and the HiSD Eq. 8 results from $J(\bm{x})\neq J(\bm{x})^{\top}$. However, if $J^{*}$ has a pair of complex eigenvalues with positive real parts, $(\bm{x}^{*},\bm{v}_{1},\cdots,\bm{v}_{k})$ cannot be a stable stationary point of GHiSD Eq. 13 for any $(\bm{v}_{1},\cdots,\bm{v}_{k})$. In spite of the loss of stable convergence of each $\bm{v}_{i}$, the subspace $\mathrm{span}\{\bm{v}_{1},\cdots,\bm{v}_{k}\}$ can still converge stably to $\mathcal{W}^{\mathrm{u}}(\bm{x})$, which also achieves our goal. Therefore, GHiSD Eq. 13 can also be applied to search hyperbolic high-index saddle points with complex eigenvalues.

Now we consider the numerical implementation of GHiSD. The dynamics Eq. 13b of GHiSD iterates $\bm{v}_{i}$ along $J(\bm{x})\bm{v}_{i}$ while maintaining the orthonormal condition Eq. 15. Applying an explicit Euler scheme with a sufficiently small step size $\beta>0$, we can establish $\mathcal{W}^{\mathrm{u}}(\bm{x})$ via

which can be regarded as a power method of $I+\beta J(\bm{x})$ to calculate the eigenvectors of the eigenvalues of $J(\bm{x})$ with positive real parts at an exponential rate. Here $\mathrm{orth}(\cdot)$ is a normalized orthogonalization function, which can be realized by a Gram–Schmidt procedure.

Since the Jacobian is usually expensive to compute and store in practice, we adopt the dimer method to avoid evaluating the Jacobian explicitly [25]. A dimer centered at $\bm{x}$ with a direction of $\bm{v}_{i}$ and length $2l$ is applied to approximate $J(\bm{x})\bm{v}_{i}$ in Eq. 26:

which is essentially a central difference for directional derivatives and $l>0$ is a small constant. Eventually, by applying an explicit Euler scheme to the dynamics Eq. 13a of GHiSD, we obtain a numerical algorithm for searching $k$-saddle of non-gradient systems:

where $\alpha,\beta>0$ are step sizes. In practice, the discretization scheme Eq. 28 is implemented with sufficiently small step sizes to ensure numerical stability and convergence.

In this subsection, we introduce a systematic numerical procedure to construct the solution landscape of a dynamical system. This procedure consists of a downward search algorithm and an upward search algorithm, which is a generalization of constructing the pathway map [45].

The downward search is the core procedure to search stationary points starting from a high-index saddle. Given a $k$-saddle $\hat{\bm{x}}$, let $\hat{\bm{v}}_{1},\cdots,\hat{\bm{v}}_{k}$ denote the orthonormal basis of $\mathcal{W}^{\mathrm{u}}(\hat{\bm{x}})$ in Theorem 2.1. We choose a direction $\hat{\bm{v}}_{i}$ from the unstable directions $\{\hat{\bm{v}}_{1},\cdots,\hat{\bm{v}}_{k}\}$, and slightly perturb the parent state $\hat{\bm{x}}$ along this direction. An $m$-GHiSD $(m<k)$ is started from the point $\hat{\bm{x}}\pm\epsilon\hat{\bm{v}}_{i}$ and the $m$ initial directions from the unstable directions need to exclude $\hat{\bm{v}}_{i}$. A typical choice for $m$-GHiSD in a downward search is $(\hat{\bm{x}}\pm\epsilon\hat{\bm{v}}_{m+1},\hat{\bm{v}}_{1},\cdots,\hat{\bm{v}}_{m})$. This procedure is repeated to newly-found saddles until sinks are found. The downward search algorithm is presented in detail as Algorithm 1.

The downward search algorithm drives a systematic search of stationary points starting from a given parent state in a controlled procedure. On the other hand, if the parent state is unknown or multiple parent states exist, an upward search algorithm is adopted to search high-index saddles from a sink (or a saddle). For an upward search, more directions $\hat{\bm{v}}_{1},\cdots,\hat{\bm{v}}_{K}$ at $\hat{\bm{x}}$ need to be computed using Eq. 27, where $K>k$ is the highest index of the saddle to search. Similarly, we choose a direction $\hat{\bm{v}}_{i}$ from the stable directions $\{\hat{\bm{v}}_{k+1},\cdots,\hat{\bm{v}}_{K}\}$, and start $m$-GHiSD from $\hat{\bm{x}}\pm\epsilon\hat{\bm{v}}_{i}$, where the $m$ initial directions should include $\hat{v}_{i}$. A typical choice for $m$-GHiSD $(m>k)$ in an upward search is $(\hat{\bm{x}}\pm\epsilon\hat{\bm{v}}_{m},\hat{\bm{v}}_{1},\cdots,\hat{\bm{v}}_{m})$. Algorithm 2 presents the upward search algorithm.

By using a combination of downward searches and upward searches, the entire search can navigate up and down systematically to find the complete solution landscape, as long as the saddle points are connected somewhere.

We first study a simple three-dimensional (3D) dynamical system to illustrate the solution landscape:

All stationary points of Eq. 29 are listed in Table 1 and the complete solution landscape is shown in Fig. 2. To search the solution landscape, we first find the source $a_{1}$ using 3-GHiSD, which is the inverse dynamics of Eq. 29. Then we apply the downward search algorithm to construct the complete solution landscape, as shown in Fig. 2(B). Fig. 2(A) shows a 3D diagram of the solution landscape, which cannot be visualized in a high-dimensional space.