High-index Saddle Dynamics (HiSD)

• High-index Saddle Dynamics (HiSD) is a method for identifying index‑k saddle points by reversing the gradient flow on the unstable eigenspace while descending along stable directions.
• The approach evolves both the state and an orthonormal frame that spans the unstable subspace, ensuring accurate tracking of negative curvature directions and Hessian eigenstructure.
• HiSD extensions include manifold-constrained, Hessian-free, accelerated, preconditioned, and neural-surrogate variants, all aimed at constructing comprehensive solution landscapes.

High-index Saddle Dynamics (HiSD) is a dynamical-systems approach for computing saddle points of a prescribed Morse index and constructing solution landscapes. In its standard finite-dimensional form, HiSD evolves a state variable together with an orthonormal frame spanning the unstable eigenspace, and modifies gradient flow by reversing the sign of the force along that frame while retaining descent on its complement. The literature represented here includes analyses of explicit Euler discretizations, manifold-constrained systems, Hessian-free dimer realizations, accelerated and preconditioned variants, degenerate critical manifolds, neural-surrogate formulations, and PDE-constrained or continuous-in-space extensions (Luo et al., 2022, E et al., 2010).

1. Precursor dynamics and conceptual origin

The clearest precursor to HiSD is the gentlest ascent dynamics (GAD), introduced as a modified flow in which saddle points of the original system, especially index-1 saddles, become stable attractors of an augmented dynamical system (E et al., 2010). For a gradient system

$\dot{x}=-\nabla V(x),$

the basic GAD is

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$

under the normalization $(v,v)=1$. The variable $v$ evolves toward the eigenvector of the Hessian corresponding to the smallest eigenvalue, and the $x$-equation reverses the force component in that direction. In geometric terms, GAD performs descent in all directions except the selected unstable one, where it ascends.

The same paper also formulates GAD for non-gradient systems $\dot{x}=F(x)$, where one must track both right and left eigendirections because the Jacobian need not be symmetric. In that setting the augmented dynamics preserves the normalizations $(v,v)=1$ and $(w,v)=1$, and the $x$-equation reverses the component of $F(x)$ associated with the unstable eigendirection. Its linearized spectrum at a fixed point is given explicitly by

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$0

which yields the criterion that index-1 saddles of the original system become linearly stable fixed points of GAD.

For high-index saddles, the same paper already states the essential HiSD idea: identify the unstable subspace and reverse the flow on that entire subspace while leaving the stable complement unchanged (E et al., 2010). It gives explicit index-2 constructions for both a complex-conjugate unstable pair and two positive real unstable eigenvalues, including a deflation strategy for successive unstable modes. However, those formulations remain a precursor rather than a full general $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$1-index theory, because the paper does not write out a complete general $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$2-index system or prove a general high-index stability theorem.

2. Core mathematical formulation

In the later HiSD literature, the target is typically a smooth energy $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$3. A point $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$4 is an index-$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$5 saddle if

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$6

and the Hessian $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$7 has exactly $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$8 negative eigenvalues and no zero eigenvalue; equivalently, the index is the dimension of the unstable subspace of the Hessian (Luo et al., 2022). In a sign convention frequently used in analysis, the continuous HiSD system is

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$9

with $(v,v)=1$0 (Luo et al., 2022). Equivalent formulations use $(v,v)=1$1 and $(v,v)=1$2, so that

$(v,v)=1$3

(Zhang et al., 2021).

The operator

$(v,v)=1$4

is the defining mechanism of HiSD. It reflects the gradient flow along the unstable eigenspace: along stable directions the method descends, while along unstable directions it reverses sign and effectively ascends (Luo et al., 2022). The direction variables $(v,v)=1$5 are intended to form an orthonormal basis of the unstable eigenspace associated with the $(v,v)=1$6 smallest Hessian eigenvalues. Within the local theory used in later convergence analyses, linearly stable steady states of this dynamics are exactly index-$(v,v)=1$7 saddles (Luo et al., 2022).

A generalized non-gradient version also appears in the numerical-analysis literature. For autonomous systems $(v,v)=1$8 with Jacobian $(v,v)=1$9, the generalized HiSD evolves

$v$0

together with

$v$1

so that the role of the symmetric Hessian is replaced by a symmetrized Jacobian contribution in the cross-direction terms (Zhang et al., 2021). This formulation extends the ascent–descent splitting beyond gradient systems while preserving the central unstable-subspace interpretation.

3. Discretization and local convergence theory

The practically analyzed discrete HiSD is an explicit Euler update in the state variable, coupled with an eigensolver step at each iteration (Luo et al., 2022). In the notation of the convergence theory,

$v$2

followed by

$v$3

The theory is local near a target saddle $v$4, under $v$5 regularity, local Lipschitz continuity of the Hessian, and a fixed sign pattern

$v$6

in a neighborhood of $v$7 (Luo et al., 2022).

A basic technical device is the recursion

$v$8

where $v$9 is the linearized contraction matrix and $x$0 is a higher-order Taylor remainder (Luo et al., 2022). For exact eigenspaces and the step size

$x$1

the one-step estimate is

$x$2

and this extends unchanged from index $x$3 to general index $x$4 because, after reflection of the $x$5 unstable directions, the effective linearized operator has eigenvalues in $x$6 independent of $x$7 (Luo et al., 2022). Consequently, in the exact-eigenspace theory, the saddle index itself does not worsen the local rate; the asymptotic contraction is controlled by curvature through $x$8.

For inexact eigendirections, the same paper measures error either by an angle condition in the index-1 case or, for general $x$9, by projector distance

$\dot{x}=F(x)$0

The resulting contraction factor degrades explicitly with both $\dot{x}=F(x)$1 and $\dot{x}=F(x)$2, which makes precise the paper’s main message: the local behavior of discrete HiSD is governed primarily by local Hessian curvature and eigenspace accuracy (Luo et al., 2022).

A complementary line of analysis treats the discrete scheme as a time discretization of the continuous HiSD trajectory rather than only as a local solver near a saddle. For explicit Euler with normalization or Gram–Schmidt orthonormalization, the global-in-time finite-interval error estimates are

$\dot{x}=F(x)$3

and

$\dot{x}=F(x)$4

showing first-order trajectory accuracy despite the strong nonlinearity and the nonstandard orthonormalization step (Zhang et al., 2021). A later analysis explains why the commonly used implementation that drops the continuous Lagrange multiplier terms and instead applies Gram–Schmidt still works: to first order in the step size, the combined effect of the provisional update and Gram–Schmidt recovers the same directional dynamics, so the continuous Stiefel-manifold preservation mechanism is “nearly a Gram–Schmidt process” (Zhang et al., 2024).

4. Constrained and manifold HiSD

A major extension concerns saddle search on manifolds, particularly the unit sphere $\dot{x}=F(x)$5. For sphere-constrained HiSD, the state variable must remain on the sphere and the direction variables must remain tangent and orthonormal. Writing $\dot{x}=F(x)$6 and $\dot{x}=F(x)$7, the continuous system is

$\dot{x}=F(x)$8

and

$\dot{x}=F(x)$9

with constraints

$(v,v)=1$0

preserved by the flow (Zhang et al., 2023).

The corresponding discrete geometry uses three operations. First, after computing an off-sphere predictor $(v,v)=1$1, one retracts by normalization: $(v,v)=1$2 Second, one transports each direction to the new tangent space by tangent projection: $(v,v)=1$3 Third, one applies Gram–Schmidt orthonormalization: $(v,v)=1$4 A semi-implicit scheme based on these steps is proved to be first-order accurate in time,

$(v,v)=1$5

despite the coupling introduced by retraction, vector transport, and current-step orthonormalization (Zhang et al., 2023).

An earlier explicit-scheme analysis on the sphere established the same $(v,v)=1$6 trajectory accuracy and then refined it to an index-robust estimate by changing both the norm and the parameter scaling (Zhang et al., 2022). In that work the averaged direction-error norm

$(v,v)=1$7

is paired with the scaling

$(v,v)=1$8

which yields

$(v,v)=1$9

with $(w,v)=1$0 independent of the saddle index $(w,v)=1$1 (Zhang et al., 2022). The same sphere paper is also presented as a special case of a broader equality-constrained HiSD framework in which the constraint manifold is defined by $(w,v)=1$2 (Zhang et al., 2023).

5. Hessian-free, accelerated, and preconditioned variants

One major practical direction replaces exact Hessian actions by finite-difference dimer approximations. In the shrinking-dimer saddle dynamics (SSD), the any-index HiSD system is modified by replacing $(w,v)=1$3 with

$(w,v)=1$4

and by evolving the dimer half-length according to

$(w,v)=1$5

The resulting method is Hessian-free, but its consistency defect is $(w,v)=1$6, so the numerical analysis shows that to retain first-order temporal accuracy one must match the dimer length and time step through

$(w,v)=1$7

Under this condition, the explicit Euler plus dimer plus Gram–Schmidt scheme has global first-order accuracy, and Richardson extrapolation upgrades the approximation from $(w,v)=1$8 to $(w,v)=1$9 (Zhang et al., 2022).

Another line of work adds momentum. The accelerated HiSD (A-HiSD) modifies only the $x$0-update by a heavy-ball term,

$x$1

while keeping the eigensolver step for the direction vectors (Luo et al., 2023). The local theory shows that the reflected Hessian at the saddle is positive definite, so the method behaves like heavy-ball iteration on a positive-definite operator after unstable-direction reflection. The resulting local factor improves from the conventional HiSD dependence $x$2 to essentially

$x$3

with $x$4, and the paper proves that A-HiSD is especially advantageous for ill-conditioned problems (Luo et al., 2023).

A more geometric conditioning remedy is preconditioned HiSD (p-HiSD), which equips $x$5 with the $x$6-inner product

$x$7

for an SPD preconditioner $x$8 (Huang et al., 26 Mar 2026). In that metric, unstable directions are generalized eigenvectors of

$x$9

and the reflection becomes

$F(x)$0

The continuous and discrete theory shows that critical points and Morse indices are invariant under this metric change, while the local convergence rate is governed by the preconditioned condition number

$F(x)$1

for the generalized spectrum (Huang et al., 26 Mar 2026). In exact eigenspaces the discrete contraction is

$F(x)$2

and the iteration complexity improves from

$F(x)$3

(Huang et al., 26 Mar 2026).

6. Nonlocal search and degenerate saddles

Classical HiSD is fundamentally local, which motivates methods that either weaken initialization dependence or extend the theory to degenerate saddles. The improved HiSD (iHiSD) introduces a crossover from gradient flow to traditional HiSD,

$F(x)$4

with $F(x)$5 evolving from $F(x)$6 to $F(x)$7 (Su et al., 6 Feb 2025). At $F(x)$8 the dynamics is pure gradient ascent or descent; at $F(x)$9 it is HiSD. Under Morse–Smale assumptions, the paper proves a nonlocal convergence statement: if there exists a path of $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$00 from a stationary point $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$01 to an index-$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$02 stationary point $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$03, then there exists an iHiSD trajectory with

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$04

It further proves that any two stationary points can be connected by a sequence of iHiSD trajectories, thereby partly answering the completeness issue for solution-landscape construction (Su et al., 6 Feb 2025).

Degeneracy caused by nullspaces leads to a different modification. The nullspace-preserving HiSD (NPHiSD) is designed for degenerate multiple-solution problems in which the Hessian has a nontrivial nullspace in addition to unstable and stable subspaces (Jiang et al., 28 Oct 2025). Its central idea is to select ascent directions from the orthogonal complement of a fixed nullspace basis, so that the upward search excludes symmetry-induced zero modes. To reduce the cost of repeated nullspace updates, the search is divided into segments, and the paper proves that a segment admits $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$05 efficient ascent directions if and only if

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$06

where $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$07 measures principal-angle drift between current and frozen nullspaces, and $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$08 accumulates directional mismatch (Jiang et al., 28 Oct 2025).

A complementary convergence theory for degenerate critical manifolds uses Morse–Bott structure. For a connected critical manifold $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$09 of degenerate index-$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$10 saddle points with nullity $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$11, HiSD is proved to converge locally provided

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$12

and the selected frame contains the true negative modes together with arbitrary orthonormal directions in the near-zero eigenspace (Luo et al., 2 Feb 2026). The continuous dynamics is asymptotically stable toward the critical manifold, and the discrete Euler scheme satisfies a linear recursion for the distance $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$13 to the manifold,

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$14

which yields local linear convergence (Luo et al., 2 Feb 2026). The same paper also proves a gradient alignment tendency: asymptotically, the gradient direction aligns with the Hessian eigenvector associated with the smallest nonzero curvature magnitude.

7. PDE, optimal-control, and neural-surrogate extensions

HiSD has also been transplanted beyond finite-dimensional explicit energy landscapes. In one direction, continuous-in-space PDE formulations replace the usual “discretize first” viewpoint. For semilinear elliptic PDEs, the index-1 case was formulated as a coupled parabolic system equivalent to gentlest ascent dynamics, and its semi-discrete and fully discrete finite element schemes were shown to be well posed, $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$15-stable, and first-order accurate in time with index-preservation consequences (Zhang et al., 3 Aug 2025). A later general-index spatiotemporal HiSD for semilinear elliptic problems introduced a fully discrete retraction-free, orthonormality-preserving scheme, established gradient stability and full space–time error estimates, and extended the framework to semilinear advection-reaction-diffusion equations (Zhang et al., 13 Jan 2026).

A different extension treats non-convex optimal control of elliptic equations. In PDE-constrained HiSD (PCHiSD), the reduced control functional is differentiated through a state equation and an adjoint equation, and the saddle search is run in the control variable using reduced gradients and dimer-type approximations of reduced Hessian-vector products (Du et al., 30 Nov 2025). The paper uses these index-$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$16 saddle searches to build a control landscape, with higher-index saddles serving as parent states for downward searches to lower-index saddles and minima.

When the energy is unavailable explicitly or too expensive to evaluate, HiSD can be run on a neural surrogate. In neural network-based HiSD (NN-HiSD), one trains

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$17

from sampled energies using

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$18

optionally augmented by a gradient-matching term

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$19

and then applies HiSD to the surrogate using automatic differentiation for $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$20 and $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$21 (Liu et al., 2024). The local theory proves that if the surrogate error is sufficiently small in gradient and Hessian norms, then there exists a nearby surrogate saddle $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$22 and NN-HiSD converges locally to it with a rate controlled by the perturbed condition number

$$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$23

(Liu et al., 2024). This formulation is used not only on analytic benchmarks but also on data-driven landscapes such as alanine dipeptide and bacterial ribosomal assembly intermediates.

Across these extensions, the defining structure of HiSD remains unchanged: an orthonormal $$\dot{x}=-\nabla V(x)+2(\nabla V(x),v)\,v,\qquad \dot{v}=-\nabla^2V(x)v+(v,\nabla^2V(x)v)\,v,$$24-frame tracks unstable directions, the state dynamics performs ascent on that frame and descent on its complement, and the computational objective is not merely one saddle point but the organization of stationary states into a solution landscape. The modern literature shows that this structure can be analyzed on Euclidean spaces, manifolds, degenerate critical manifolds, surrogate models, and PDE-governed systems, while retaining explicit control of curvature, eigenspace accuracy, and Morse index (Zhang et al., 13 Jan 2026, Liu et al., 2024).