Documents Live, a web authoring and publishing system
If you see this, something is wrong
Table of contents
First published on Tuesday, Oct 27, 2026 and last modified on Tuesday, Oct 27, 2026 by François Chaplais.
Flatiron Institute, Simons Foundation
Abstract
Coherent sets—regions of phase space that resist mixing—organize transport in stochastic dynamical systems and are encoded in the leading singular functions of the stochastic Koopman operator. These are typically computed numerically and globally; a local, analytic theory near a fixed point has been missing, although fixed points and the transitions between them organize the flow. We supply the local theory for a saddle, the transition state through which probability flows between basins. Naively the construction fails—probability escapes along unstable directions, the invariant measure is non-normalizable, the singular functions leave \( L^2\) , and singular values exceed one. We show that conditioning on non-escape (a Doob’s \( h\) -transform) repairs this; that the required \( h\) is the unique stabilizing solution of an algebraic Riccati equation; and that the resulting slow modes decay at the escape rate along the unstable manifold. We obtain these modes in closed form—products of Hermite polynomials in canonical coordinates times the \( h\) -function—and use their unitary structure to derive a data-driven estimator: canonical correlation analysis on conditioned trajectories. When a single unstable direction dominates, the leading mode is the reaction coordinate, locally recovering the stable manifold and committor of transition-path theory. The construction connects coherent sets to transition-path theory, to the stabilizing Riccati solution of linear control, and to the time-reversed dynamics of generative diffusion models.
1 Introduction
A nonlinear dynamical system is organized by its fixed points and by the transitions between them: trajectories dwell near stable equilibria and cross between them through the neighborhoods of unstable equilibria [1, 2]. Local behavior around fixed points therefore controls the global dynamics, and characterizing it is a natural route to understanding the system. For deterministic dynamics, this local analysis is standard—the Jacobian’s eigenvectors give normal modes and its eigenvalues give growth rates. For stochastic dynamics one would like the analogous local theory of the slow, collective modes that govern how probability is transported through the neighborhood.
Those modes are the singular functions of the stochastic Koopman operator whose singular values lie close to one, the spectral objects underlying coherent sets [3, 4, 5]: subsets of phase space that remain dynamically coherent rather than dispersing. At present, coherent sets are studied almost exclusively numerically and globally: one discretizes a transfer or Koopman operator over the entire phase space of a nonlinear system—using Ulam-type discretizations [6, 7], finite elements [8], extended dynamic mode decomposition [9], or diffusion maps [10]—and reads off a partition [3, 5]. This yields neither a closed form nor a local law: no statement of what the coherent modes look like near a given fixed point. The present paper supplies that missing local, analytic description.
A stable fixed point (a sink) linearizes to a contracting Ornstein–Uhlenbeck (OU) process: an isolated metastable well with no transport through it. The dynamically interesting fixed points are saddles—hyperbolic points with at least one unstable direction. A saddle is a transition state: the bottleneck through which trajectories pass between basins, where reaction coordinates and committor functions are defined [11, 12, 13]. Its linearization is a continuous OU process whose drift matrix has eigenvalues of both signs, and this is the regime we solve.
The breakdown is easy to see. When the drift matrix has unstable eigenvalues, there is no normalizable stationary measure: probability escapes along the unstable manifold. As a consequence, the naive extension of singular functions from stable fixed points to saddles leads to non-integrable functions and singular values greater than one [14] (Fig. 6). This breaks the singular-value problem at its foundation. The singular value represents the correlation coefficient between the left singular function evaluated at the initial time and the right singular function evaluated at time \( \tau\) . In the coherent-set picture, where these functions approximate indicators of regions, this correlation is essentially the fraction of trajectories that begin in the coherent set and remain there after time \( \tau\) —a survival probability, which cannot exceed one. The breakdown is not merely analytic: the data-driven route fails in tandem, because the empirical covariances that any CCA-type estimator [15, 16, 5] relies on diverge along the unstable directions as the trajectories escape, so the correlation matrix has no finite limit and the canonical problem is ill-posed on the raw data (Fig. 6). Any correct theory of the saddle case must first restore a proper Hilbert space.
Our resolution, and the conceptual core of the paper, is to condition the dynamics: restricting to trajectories that have not yet escaped is a Doob’s \( h\) -transform [17] that stabilizes the drift and makes the slow modes well defined. We show (i) that the required conditioning is the unique stabilizing solution of an algebraic Riccati equation [18], selected automatically as the attractor of a Riccati flow (Sec. 3); (ii) that the conditioned slow modes admit closed form whose decay rate is the escape rate along the unstable manifold, and that multiplication by \( h\) is a unitary map relating them to ordinary Hermite modes (Sec. 4); (iii) that this unitary structure yields a data-driven estimator—canonical correlation analysis on conditioned trajectories—which recovers the analytic modes exactly (Sec. 5); (iv) that in the special case with a single slow unstable mode, the leading predictive direction converges to the unstable left eigenvector, justifying the reaction coordinate in classical theories (Sec. 6); and (v) that the construction establishes bridges to transition-path theory, stochastic thermodynamics, linear control, and diffusion-model time reversal (Sec. 7).
2 Setup and main result
2.1 Local linearization of stochastic nonlinear dynamics
Consider a stochastic nonlinear system on \( \mathbb{R}^d\) ,
(1)
where \( \mathbf b:\mathbb{R}^d\to\mathbb{R}^d\) is the drift and \( D\succ0\) is the constant diffusion matrix. Let \( \mathbf x_*\) be a fixed point, \( \mathbf b(\mathbf x_*)=0\) . Writing \( A=\nabla\mathbf b(\mathbf x_*)\) for the Jacobian and shifting the origin to \( \mathbf x_*\) , the local dynamics is the OU process
(2)
We assume for clarity that the eigenvalues and eigenvectors of \( A\) are real (the complex case is discussed in Appendix S4 D). Denote by \( \lambda_i\) the eigenvalues of \( A\) , by \( r_i\) the right eigenvectors (\( A r_i=\lambda_i r_i\) ), and by \( \ell_i\) the left eigenvectors (\( A^{\top}\ell_i=\lambda_i\ell_i\) ), normalized biorthogonally, \( \ell_i^{\top}r_j=\delta_{ij}\) . We split the spectrum into a stable part (eigenvalues \( \lambda_i^-<0\) , eigenvectors \( r_i^-,\ell_i^-\) ) and an unstable part (eigenvalues \( \lambda_i^+>0\) , eigenvectors \( r_i^+,\ell_i^+\) ), collect the unstable left eigenvectors into the columns of \( L_{+}\) and the unstable eigenvalues into \( \Lambda_{+}=\mathrm{diag}(\lambda_i^+)\) , so that \( A^{\top}L_{+}=L_{+}\Lambda_{+}\) . A saddle is the case where both parts are nonempty.
2.2 Fokker–Planck operator and the Koopman generator
The density \( p(\mathbf x,t)\) obeys the Fokker–Planck equation \( \partial_t p=\mathcal{L} p\) with the forward generator
(3)
Its adjoint with respect to the Euclidean inner product \( \langle f,g\rangle=\int f g\,d\mathbf x\) is the backward Kolmogorov operator, which is exactly the generator of the finite-time stochastic Koopman operator:
(4)
where the Koopman operator \( \mathcal{K}_\tau\) propagates observables \( g:\mathbb{R}^d\to\mathbb{R}\) by the conditional expectation [19]
(5)
A stationary measure satisfying \( \mathcal{L}\rho_0=0\) is given by \( \rho_0(\mathbf x)\propto\exp(-\tfrac12\mathbf x^{\top}\Sigma^{-1}\mathbf x)\) , with \( \Sigma\) solving the Lyapunov equation
(6)
Eq. (6) has a unique solution iff no two eigenvalues of \( A\) sum to zero [20], and \( \Sigma\succ0\) is the covariance iff \( A\) is Hurwitz (all \( \lambda_i<0\) ). If \( A\) has unstable eigenvalues, \( \Sigma\) is indefinite (negative definite if all modes are unstable); then \( \rho_0\) grows in the unstable directions and is not normalizable. Although such a non-normalizable \( \rho_0\) does not define a probability distribution on \( \mathbb{R}^d\) , it remains a well-defined local stationary profile for the linearized dynamics near the fixed point. This is the root of the obstruction described in Sec. I.
2.3 Singular functions of the Koopman operator
Coherent modes are defined variationally. Given a reference measure \( \mu\) with \( \langle f,g\rangle_\mu\coloneq\int fg\,d\mu\) and \( \|f\|_\mu\coloneq\sqrt{\langle f,f\rangle_\mu}\) , the leading singular value and functions solve
(7)
where the maximizers are the leading left singular function \( f_1\) and right singular function \( g_1\) . Further modes maximize the same ratio subject to \( \mu\) -orthogonality to the previous ones. Introducing the \( \mu\) -adjoint \( \mathcal{K}_\tau^{\dagger}\) via \( \langle f,\mathcal{K}_\tau g\rangle_\mu=\langle\mathcal{K}_\tau^{\dagger}f,g\rangle_\mu\) , and the forward–backward and backward–forward operators \( \mathcal F_\tau=\mathcal{K}_\tau\mathcal{K}_\tau^{\dagger}\) and \( \mathcal B_\tau=\mathcal{K}_\tau^{\dagger}\mathcal{K}_\tau\) , the singular functions are equivalently given by the following eigenfunctions [19]:
(8)
Typically, \( \mu\) is taken as the invariant measure. However, when \( A\) has unstable modes, the invariant measure is non-normalizable, the naive stable-OU polynomial singular functions are not integrable with respect to this measure, and the ratio in Eq. (7) can exceed one—the \( \sigma>1\) pathology of Sec. I.
2.4 Summary of the main result
We show that conditioning replaces the ill-defined problem by a well-posed one on a stabilized process (Sec. 3), whose singular functions are then obtained exactly. The result, derived in Sec. 4 and Appendix S8, is the following. Conditioning produces a stabilized drift \( \widetilde{A}=A+DP_+\) (Hurwitz) with \( P_+=L_{+}(L_{+}^{\top}\Sigma L_{+})^{-1}L_{+}^{\top}\) , a stabilized covariance \( \widetilde{\Sigma }\) satisfying \( \widetilde{A}\widetilde{\Sigma }+\widetilde{\Sigma }\widetilde{A}^{\top}=-D\) , a Gaussian factor \( h(\mathbf x)=\exp(\tfrac12\mathbf x^{\top}P_+\mathbf x)\) , and canonical directions \( v_i,u_i\) with canonical correlations \( \gamma_i\in(0,1)\) (the singular values of the stabilized process, defined in Sec. 5.1). We introduce a reference weight \( \widetilde{\rho}_0\) defined in Sec. 3.3, the natural extension of \( \rho_0\) to saddles, on which \( \mathcal{K}_\tau\) is a contraction. The saddle Koopman singular functions and singular values, with respect to \( \widetilde{\rho}_0\) , are indexed by a multi-index \( \boldsymbol n\in\mathbb N_0^{\,d}\) and given in closed form by
(9)
where \( \mathrm{He}_n\) is the probabilists’ Hermite polynomial and \( \boldsymbol n!=\prod_i n_i!\) . The Gaussian factor \( h\) decays along the unstable manifold, with the decaying prefactor of the singular value \( e^{-\rm{Tr}(\Lambda_{+})\tau}\) set by the escape rate \( \rm{Tr}(\Lambda_{+})=\sum_i\lambda_i^+\) . The remainder of the paper constructs each ingredient and proves Eq. (9).
3 Stabilization by conditioning
3.1 Doob’s \( h\) -transform
Conditioning on a single endpoint. Retain only trajectories arriving at a prescribed point \( \mathbf x(T)=\mathbf y\) . The conditioned process is again a diffusion—Doob’s \( h\) -transform—with the same \( D\) and a drift modified by the gradient of the log transition density [17]. Writing \( h(\mathbf x,t)=p(\mathbf y,T\mid\mathbf x,t)\) for the probability of reaching \( \mathbf y\) from \( (\mathbf x,t)\) , the drift is modified as
(10)
Conditioning on a final distribution. More generally, one can condition on a prescribed terminal distribution \( \pi(\mathbf x_T)\) . The \( h\) -function in this case is given by the average of the single-point \( h\) -function weighted by the conditional-to-original distribution ratio (Appendix S1),
(11)
where \( w\coloneq \frac{\pi}{p}\) is the reweighting coefficient. Thus, the conditioning \( h\) -function is itself a Koopman image of the reweighting \( w\) . For a Gaussian \( w(\mathbf x)\propto\exp(-\tfrac12\mathbf x^{\top}\Sigma'^{-1}\mathbf x)\) , the long-horizon limit, \( T\gg t\) , yields a Gaussian \( h\) ,
(12)
with \( L_{+}\) the unstable left eigenvectors of Sec. 2. The proof and explanation are given in the next subsection. We stress that this long-horizon limit is independent of the reweighting window \( \Sigma'\) . Following Eq. (10), the conditioned drift matrix is modified as
(13)
The two extremes are transparent:
(14)
(15)
In the all-stable case, conditioning does nothing. In the all-unstable case, the conditioned drift is exactly the score-based reverse-time SDE of generative diffusion modeling [21, 22]: \( \widetilde{A} \mathbf x=A\mathbf x-D\nabla\log\rho_0(\mathbf{x})\) . The saddle interpolates between them.
Proposition 1 (Selective time reversal)
With \( \widetilde{A}\) from Eq. (13),
(16)
the stable eigenvalues are unchanged and the unstable ones are reflected across zero. In particular, \( \widetilde{A}\) is Hurwitz.
Proof in Appendix S3. Conditioning is thus a selective time reversal, flipping the arrow of time only along the escaping directions. The numerical verification is given in Fig. 2, and Fig. 5 shows the original and stabilized drift of an example linear saddle system.
3.2 The Riccati flow and its stabilizing attractor
We now show why the conditioning selects precisely \( P_+\) in Eq. (12). From Eq. (11), \( h_s= \mathcal{K}_s w\) with \( s=T-t\) ; the OU semigroup preserves the Gaussian class, so \( h_s\propto\exp(\tfrac12\mathbf x^{\top}P_s\mathbf x)\) with \( P_s\) symmetric, and we have (Appendix S2)
(17)
Substituting into \( \partial_s h_s=\mathcal{L}^{\dagger} h_s\) gives a matrix Riccati flow for the quadratic form (Appendix S5),
(18)
As \( s\to\infty\) this flow converges to a fixed point \( P\) solving the algebraic Riccati equation [18]
(19)
and thus \( \exp\big(\tfrac12\mathbf x^{\top}P\mathbf x\big)\) is an eigenfunction of \( \mathcal{L}^{\dagger}\) with eigenvalue \( \frac{1}{2}\rm{Tr}(DP)\) .
To see which fixed point is reached, we first describe the structure of all solutions, then identify the attractor.
Structure of the solutions. Since \( P\) enters through the quadratic form, we restrict throughout to symmetric solutions \( P=P^\top\) . Each symmetric solution \( P\) of Eq. (19) produces a conditioned drift \( \widetilde{A}=A+DP\) whose eigenvalues have a simple relation to those of \( A\) : every eigenvalue \( \lambda_i\) of \( A\) appears in \( \widetilde{A}\) either unchanged (“kept”) or negated to \( -\lambda_i\) (“flipped”). The one-dimensional case makes this concrete: for \( \dot x=\lambda x+\xi\) the Riccati equation \( 2\lambda P+DP^2=0\) has the two roots \( P=0\) (giving \( \widetilde{A}=\lambda\) , kept) and \( P=-2\lambda/D=\Sigma^{-1}\) (giving \( \widetilde{A}=-\lambda\) , flipped). In \( d\) dimensions, a solution is specified by choosing, independently for each of the \( d\) modes, whether to keep or flip it. There are therefore \( 2^d\) solutions, given explicitly by
(20)
where \( \mathcal F\subseteq\{1,\dots,d\}\) is the set of flipped modes and \( L_{\mathcal F}\) collects their left eigenvectors. Three landmarks: \( \mathcal F=\varnothing\) gives \( P=0\) ; the full set gives \( P=\Sigma^{-1}\) (so \( \widetilde{A}=-\Sigma A^{\top}\Sigma^{-1}\) ); and \( \mathcal F=\{\text{unstable modes}\}\) gives \( P_+\) in Eq. (12). The systematic derivation via the Hamiltonian matrix \( \bigl[\begin{smallmatrix}A&D\\ 0&-A^{\top}\end{smallmatrix}\bigr]\) and its invariant subspaces is in Appendix S4 [18].
Which one does the flow reach? A conditioned process that does not escape must be contracting, i.e. \( \widetilde{A}\) must be Hurwitz (all eigenvalues negative). Among the \( 2^d\) choices, exactly one achieves this: flip every unstable mode (\( \lambda_i^+>0\mapsto-\lambda_i^+<0\) ) and keep every stable mode. This is the stabilizing solution, Eq. (12). The mathematical reason for the Riccati flow (18) converging to it and to no other is the following. Linearizing the flow about any solution \( P_{\mathcal F}\) gives growth rates \( \widetilde\lambda_i+\widetilde\lambda_j\) , where \( \{\widetilde\lambda_i\}\) are the eigenvalues of \( \widetilde{A}\) , and these are all negative—so \( P_{\mathcal F}\) is attracting—only when every \( \widetilde\lambda_i<0\) , i.e. only for the stabilizing choice; every other solution has a growing direction and is a saddle of the flow (Appendix S5). Hence, long-horizon conditioning lands on Eq. (12) regardless of the conditioning window \( \Sigma'\) . We numerically verify the convergence for different windows \( \Sigma'\) , as shown in Fig. 3.
3.3 The reference measure \( \widetilde{\rho}_0\)
The stabilized process \( \dot{\mathbf x}=\widetilde{A}\mathbf x+\boldsymbol\xi\) has a stationary Gaussian density \( \rho_{\widetilde{A}}\propto\exp(-\tfrac12\mathbf x^{\top}\widetilde{\Sigma }^{-1}\mathbf x)\) , \( \widetilde{A}\widetilde{\Sigma }+\widetilde{\Sigma }\widetilde{A}^{\top}=-D\) , \( \widetilde{\Sigma }\succ0\) . In the all-stable case, \( \widetilde{\Sigma }=\Sigma\) ; in the all-unstable case, \( \widetilde{\Sigma }=-\Sigma\) . As illustrated by the following proposition, \( \widetilde{\Sigma }\) can generically be expressed in closed form, through the pair of complementary Riccati solutions.
Proposition 2 (Closed form for the stabilized covariance)
(21)
where \( P_\pm=L_\pm\big(L_\pm^{\top}\Sigma L_\pm\big)^{-1}L_\pm^{\top}\) , with \( L_\pm\) the collection of unstable/stable left eigenvectors of \( A\) .
Proof in Appendix S6. In contrast to \( P_+\) , \( P_-\) solely flips all the stable modes (Eq. (20) with \( \mathcal F=\{\text{stable}\}\) ); correspondingly \( A+DP_-\) is purely unstable. Eq. (21) is the covariance-level counterpart of the selective time reversal of Proposition 1. Thus, \( \rho_{\widetilde{A}}\propto \exp\Big[-\tfrac12\mathbf x^{\top}\big(P_--P_+\big)\mathbf x\Big]\) .
We now consider the reference measure defined as
(22)
The all-stable case gives \( P_-=\Sigma^{-1},\: P_+=0\) and the all-unstable case gives \( P_-=0,\: P_+=\Sigma^{-1}\) . Thus, in both all-stable and all-unstable extremes, \( \widetilde{\rho}_0=\rho_0\) . In the genuine saddle case, they are generically not identical because the stable and unstable directions demand opposite curvatures, and \( \widetilde{\rho}_0=\rho_{\widetilde{A}}h^{-2}\) is the minimal deformation that honestly keeps the stable manifold’s Gaussian and the unstable manifold’s reciprocal Gaussian. Thus, \( \widetilde{\rho}_0\) is the natural extension of \( \rho_0\) from stable/unstable extremes to saddles, which decouples the stable and unstable covariances. We restrict the variational maximization for singular functions, Eq. (7), to the Hilbert space \( L^2(\widetilde{\rho}_0)\) . The factorization Eq. (22) makes the relevant integrability transparent: for any polynomial \( \mathcal H\) ,
(23)
In particular, the singular functions in Eq. (9) are square-integrable in \( L^2(\widetilde{\rho}_0)\) . As shown next, this specific measure is exactly the one making conditioning a unitary change of representation, and it is what the data-driven method computes.
4 Exact singular modes and their unitary origin
The main result is Eq. (9): the saddle Koopman singular functions with respect to \( \widetilde{\rho}_0\) are products of univariate Hermite polynomials in the canonical coordinates \( v_i^{\top}\mathbf x\) , \( u_i^{\top}\mathbf x\) , multiplied by the Gaussian factor \( h\) , with singular values \( e^{-\rm{Tr}(\Lambda_{+})\tau}\prod_i\gamma_i^{n_i}\) . A direct derivation is given in Appendix S8. Here we explain the result through the structure that makes it true, which also motivates the estimator of Sec. 5.
As shown in Sec. 3.2, \( h\) is a positive eigenfunction of the Koopman generator, for \( P_+\) solving Eq. (19),
(24)
It follows that the Koopman generator of the stabilized process is given by
(25)
which leads to the following proposition.
Proposition 3 (Conditioning is a unitary change of representation)
Let \( M_h f=hf\) . Because \( \widetilde{\rho}_0=\rho_{\widetilde{A}}h^{-2}\) , \( M_h: L^2(\rho_{\widetilde{A}})\to L^2(\widetilde{\rho}_0) \) is unitary:
(26)
and it intertwines the Koopman operators of the original and stabilized processes,
(27)
Proof in Appendix S7. Since the singular value decomposition is invariant under unitary conjugation and a scalar prefactor merely rescales singular values, the singular functions of the saddle Koopman operator (with respect to \( \widetilde{\rho}_0\) ) are \( M_h\) applied to those of the stabilized process (with respect to \( \rho_{\widetilde{A}}\) ), with singular values scaled by \( e^{-\rm{Tr}(\Lambda_{+})\tau}\) . It also follows from the unitarity that \( \mathcal{K}_\tau\) is a contraction on \( L^2(\widetilde{\rho}_0)\) , so the singular value problem is well-posed (Appendix S9).
The stabilized problem is a stable-OU singular-value problem, and its singular functions are products of univariate Hermite polynomials [14]. The first-order left and right singular functions are the canonical scores \( v_i^{\top}\mathbf x\) and \( u_i^{\top}\mathbf x\) , respectively, where \( v_i,u_i\) solve
(28)
and \( \gamma_i\in(0,1)\) are the canonical correlations in descending order \( \gamma_1\ge\gamma_2\ge…\) . Crucially, in the canonical coordinates, the stationary measure factorizes into independent standard Gaussians, so the higher-order singular functions are simply products of univariate Hermite polynomials (proof in Appendix S8).
The leading coherent mode is linear. Because \( \gamma_i\in(0,1)\) , excluding the trivial mode \( \boldsymbol n=\boldsymbol 0\) , the product \( \prod_i\gamma_i^{n_i}\) is largest only if there is a single mode, i.e., \( \sum_i n_i=1\) . That is, the largest nontrivial singular value is \( \sigma_{\boldsymbol e_1}=e^{-\rm{Tr}(\Lambda_{+})\tau}\gamma_1\) , attained by a first-order mode: \( \Phi^L_1=h\,v_1^{\top}\mathbf x\) . Higher-order modes are strictly less coherent. Because coherent sets are partitioned by the sign pattern of top singular functions, and \( h>0\) , this explains why the linear filter—and hence canonical correlation analysis, which computes exactly the first-order modes (Sec. 5)—suffices to identify the dominant coherent partition.
5 Data-driven estimator
5.1 Koopman singular functions and canonical correlation analysis
The mainstream data-driven methods for approximating transfer operators from trajectory data—extended dynamic mode decomposition (EDMD) [9], the variational approach for Markov processes (VAMP) [16], and other variants—are, at their core, correlation methods: they lift the state through a dictionary of observables and solve a generalized eigenproblem of empirical instantaneous and time-lagged covariance matrices. When the object sought are the leading singular functions rather than eigenfunctions, this reduction is exactly a canonical correlation analysis (CCA) between the observables at the initial and final times [16, 5]. Our estimator is of the same family. The point of this section is that, for the saddle case, the correct CCA is performed on the conditioned (stabilized) trajectories rather than the raw ones. We first make the CCA–Koopman correspondence precise.
CCA [15] finds linear projections of two random vectors that are maximally correlated. Applied in time to \( \mathbf x(0)\) and \( \mathbf x(\tau)\) of a stationary process, it seeks \( v,u\) maximizing
(29)
with the canonical directions and the canonical correlations given by the eigenproblems of \( C_0^{-1}C_\tau C_0^{-1}C_{-\tau}\) and \( C_0^{-1}C_{-\tau} C_0^{-1}C_{\tau}\) . For a stable OU process, this is exactly the first-order singular-value problem, Eqs. (8) and (28), with \( \mu=\rho_0\) : substituting \( C_0=\Sigma\) , \( C_\tau=\Sigma e^{A^{\top}\tau}\) , \( C_{-\tau}=e^{A\tau}\Sigma\) ,
(30)
Thus, the CCA canonical scores \( v^{\top}\mathbf x\) , \( u^{\top}\mathbf x\) are the first-order Koopman singular functions and the canonical correlations \( \gamma_i\) are the singular values [23, 16, 14]. In short, CCA is the empirical singular value decomposition of the Koopman operator for stable OU processes, which is why a data-driven procedure can recover the analytic modes. However, one cannot naively extend such CCA on raw trajectories to saddle cases, because the correlation matrices, \( C_0\) and \( C_\tau\) , diverge along unstable directions (Fig. 6).
5.2 Algorithm for coherent-set partition
In this subsection, we present the correct data-driven algorithm for estimating the coherent-set partition in the local neighborhood of a fixed point. Coherent sets are partitioned by the sign pattern of top Koopman singular functions [3]. The estimation of canonical coordinates, \( \xi_i=v_i^{\top}\mathbf x\) and \( \eta_i=u_i^{\top}\mathbf x\) in Eq. (9), is all that the data-driven method needs: every singular function is a product of Hermite polynomials in them, \( \Phi^{L}_{\boldsymbol n}\propto h\prod_i\mathrm{He}_{n_i}(\xi_i)\) and \( \Phi^{R}_{\boldsymbol n}\propto h\prod_i\mathrm{He}_{n_i}(\eta_i)\) . Because the Gaussian factor \( h>0\) never changes sign, coherent sets are the sign sets of the Hermite part alone, i.e., they are determined entirely by the estimated canonical coordinates. Thus, once the \( \{(v_i,u_i)\}\) are estimated, any coherent-set partition is available in closed form without further fitting. According to Eq. (28), the estimation of \( \{(v_i,u_i)\}\) is a linear CCA on the stabilized process, obtained from the original (unstable) trajectories by the conditioning reweighting, \( w\) in Eq. (11). This gives the following estimator.
Algorithm: CCA with conditioning.
- Take the fixed point of interest as the center.
- Collect trajectories \( \{\mathbf x^{(m)}(t)\}_{t\in[0,T]}\) from any normalizable initial distribution.
- Choose a confining covariance \( \Sigma'\) ; assign trajectory \( m\) the weight \( w_m\propto\mathcal N\big(\mathbf x^{(m)}(T);0,\Sigma'\big)\) , implementing conditioning on non-escape.
- Choose lags \( 0\ll t_1<t_2\ll T\) . Under the reweighted empirical measure \( \{w_m\}\) , form the weighted instantaneous and time-lagged covariances of \( (\mathbf x(t_1),\mathbf x(t_2))\) and solve the CCA eigenproblem for the leading canonical pairs \( \{(v_i,u_i)\}\) .
- Build any desired singular function from the estimated canonical coordinates via Eq. (9), and partition phase space by its sign pattern.
Leading linear partition. The most coherent nontrivial modes are first order (Sec. 4): \( \Phi^{L}_1\propto h\,\xi_1\) and \( \Phi^{R}_1\propto h\,\eta_1\) , whose sign sets \( \{\xi_1\gtrless0\}\) and \( \{\eta_1\gtrless0\}\) are hyperplane partitions. The partitions are thus read off from the signs of the projections onto \( v_1\) and \( u_1\) , the direct analogue of the Laplacian’s leading eigenvector in spectral clustering [24]. We tested the algorithm on both stochastic linear saddle systems and nonlinear bistable dynamics with saddle fixed points. The numerically estimated leading predictive vector aligns with the theoretical expectation as the sample size increases (Fig. 7 and the right panels of Fig. 8).
5.3 Practical regime
Conditioning an escaping process has a cost: the weights concentrate on rare surviving trajectories, so the effective sample size decays roughly as \( e^{-c\,\rm{Tr}(\Lambda_{+})\,T}\) . The estimator is thus reliable for moderate horizons and modest escape rates. The window-independence guaranteed by the Riccati attractor (Sec. 3.2) provides a plateau in \( T\) over which the estimate is unbiased before the variance penalty dominates; one selects \( T\) at the short end of the plateau.
6 Special case: single, dominant unstable mode
A case of particular physical importance is a saddle with a single unstable direction whose escape rate is much lower than the relaxation rate of every stable direction, \( \lambda^+\ll\min_i|\lambda_i^-|\) (equivalently, \( \lambda^+<\min_i|\lambda_i^-|\) and the long-lag limit \( \tau\to\infty\) ). This is precisely the separation of timescales underlying Kramers’ and Eyring’s classical reaction-rate theory [11] and transition-path theory [12]: a rare, slow escape along a single reaction coordinate against a fast-relaxing bath of transverse modes. Conventionally, the reaction coordinate is identified with the unstable direction of the linearized dynamics at the transition state [25, 26, 11]; we show here that in this regime, the full saddle construction of Secs. 3–5 justifies that identification with exact closed-form solutions.
The left and right singular functions of the Koopman operator correspond to predictive and retrospective coherent sets, respectively. Given the time-reversal symmetry, we focus on prediction in the rest of this section. Write \( \ell^+,r^+\) for the unique unstable left and right eigenvectors of \( A\) , \( \lambda^+\) for the unstable eigenvalue, and \( \kappa\coloneq \ell^{+\top}\Sigma\ell^+\) . The stabilizing solution (Eq. (12)) is rank one, \( P_+=\ell^+\ell^{+\top}/\kappa\) .
Proposition 4 (Left eigenvectors of \( \widetilde{A}\) for single unstable mode)
For a saddle \( A\) with a single unstable mode,
(31)
In particular, the unstable left eigenvector of \( A\) is a left eigenvector of \( \widetilde{A}\) with flipped eigenvalue.
Proof in Appendix S10. Proposition 1 already gives the matching statement for right eigenvectors (\( \widetilde{A} r_i^-=\lambda_i^-r_i^-\) , \( \widetilde{A}\Sigma\ell^+=-\lambda^+\Sigma\ell^+\) ); Eq. (31) is its left-eigenvector counterpart. The flip preserves \( \ell^+\) exactly whenever there is only one unstable mode to flip. With the additional assumption of timescale separation such that the unstable mode dominates the long-lag limit of the stabilized transfer operator, the leading local predictive direction simply converges to \( \ell^+\) , which is orthogonal to the stable subspace (Fig. 5).
Proposition 5 (Long-lag convergence of the leading predictive direction)
If \( \lambda^+<\min_i|\lambda_i^-|\) , then \( -\lambda^+\) is the eigenvalue of \( \widetilde{A}\) closest to zero, so \( \ell^+\) is the left eigenvector of the slowest mode of \( \widetilde{A}\) and
(32)
Proof in Appendix S10. We now interpret the physical meaning of Proposition 5 by connecting it to the stable manifold and the committor. The stable manifold of a fixed point \( \mathbf x_*\) is the set of points whose deterministic orbit converges to \( \mathbf x_*\) rather than escaping [27]. In the single-unstable regime of this section, the stable manifold partitions the phase space into two disjoint regions whose states possess distinct asymptotic future behaviors. For a linear saddle, the hyperplane orthogonal to the left unstable eigenvector, \( \{\mathbf{x}:(\mathbf{x}-\mathbf{x}_*)^\top \ell^+=0\}\) , is exactly the stable manifold. For a saddle fixed point in a genuinely nonlinear system, this hyperplane is tangent to the stable manifold, serving as a local linear approximation. Thus, the projection of states near the saddle onto the local predictive direction \( v_1\) forecasts the asymptotic future evolution. We numerically verified that the estimated \( \hat{v}_1\) is perpendicular to the stable manifold at the saddle (middle panels of Fig. 8).
For a stochastic system, the mathematical object that more directly quantifies the probabilistic future evolution is the committor \( q(\mathbf{x})\) , which is defined as the probability of reaching basin A before basin B conditioned on the initial state \( \mathbf x\) [12]. The committor depends on the global definitions of the two basins and therefore need not, in general, be determined by the local linearization. Nevertheless, for a linear saddle with basin boundaries compatible with the unstable coordinate, the committor is exactly an affinely normalized error function of \( \ell^{+\top}(\mathbf{x}-\mathbf{x}^*)\) , and its committor contours are hyperplanes normal to \( \ell^+\) (Appendix S11 and Ref. [28]). We empirically observed that the local predictive direction is orthogonal to the contour of the committor at the saddle (middle panels of Fig. 8). These results verify that the singular functions we construct, along with the numerical estimator, can yield local predictive representations consistent with classical theories.
The single-unstable case also closes the loop on the stabilized covariance with an inverse-free formula:
Proposition 6 (Stabilized covariance for single unstable mode)
For a saddle \( A\) with a single unstable mode,
(33)
Proof in Appendix S10. Eq. (33) is the covariance-level flip of the unstable mode. As a check, in the all-unstable corner (\( d=1\) , \( r^+=\ell^+=1\) ) this reduces to \( \widetilde{\Sigma }=\Sigma-2\Sigma=-\Sigma\) , matching Eq. (15). We note that Eqs. (31) and (33) do not extend to multiple unstable modes by the naive matrix analogue.
Equations (32) and (33) together give \( v_1,\widetilde{\Sigma }\) in closed form. In chemical reaction rates, this provides the stabilized transverse covariance directly, wherever transition-state theory’s separation of timescales holds [25, 26, 11]. In enhanced-sampling methods, a central practical problem is identifying the reaction coordinate to bias simulations along [29, 30], and our algorithm provides a data-driven estimator for that coordinate with a proven guarantee of what it converges to, rather than a heuristic collective variable chosen by hand [31]. The same closed form applies to any metastable system with this structure, independent of the specific physical mechanism.
7 Discussion and connections
Transition-path theory. A saddle is a transition state; its unstable manifold is the linearized reaction coordinate, and the conditioned process is the local linear theory of the reactive ensemble. When the saddle has a single unstable direction whose escape is the slowest process in the system—the classical rare-event regime of Kramers’ and Eyring’s rate theory [11, 12]—the leading singular function is aligned with the linear committor direction (Sec. 6).
Linear control and stochastic realization. The stabilization solution is the attractor of an algebraic Riccati flow [18], and the eigenvalue reflection of Proposition 1 is familiar from acausal stochastic realization [32, 23]. The agreement is not only formal: \( \widetilde{A}\) is the optimal closed-loop generator of the linear-quadratic regulator, and the same solution is selected there because the regulated dynamics must not diverge—the control-theoretic counterpart of requiring the conditioned process not to escape. Canonical correlations between past and future of a stable stochastic system are likewise classical [23, 32]; what is new here is the escaping case, where no invariant probability measure exists, and the correlation problem must first be repaired by conditioning.
Stochastic thermodynamics. The selective time reversal of Proposition 1 reverses the time-irreversible part of the drift only on the unstable manifold, flipping the entropy-producing arrow on the escaping directions while leaving the recurrent ones untouched—a subspace-resolved statement about dissipation.
Generative diffusion models. The all-unstable case, Eq. (15), is exactly the reverse-time score-based SDE of generative sampling [21, 22]. It is the extreme point of a family whose interior is the saddle, with conditioning acting only on the escaping subspace: reverse-time generative dynamics and transport through a transition state are two ends of one construction.
Escape rate and survival. The decay \( e^{-\rm{Tr}(\Lambda_{+})\tau}\) sets the rate at which coherence is lost, and \( (\mathcal{L}^{\dagger},h,-\rm{Tr}(\Lambda_{+}))\) is the principal eigen-triple of a killed (conditioned) semigroup: \( h\) is the survival function, small where escape is imminent, and \( \rho_{\widetilde{A}}=h\nu\) is the stationary law of the associated \( Q\) -process, with \( \nu\) its dual eigenfunction living on the unstable side [33]. The same mechanism underlies the metastable spectra of open systems more broadly: a decay rate becomes visible only once an invariant probability measure is given up in favor of a survival-weighted or otherwise non-normalizable reference, as in the Pollicott–Ruelle resonances of chaotic flows and the Gamow states of quantum scattering [34, 35]. Here that role is played by the reference weight \( \widetilde{\rho}_0\) itself.
Limitations. Real eigenvalues are assumed for clarity only; the complex case extends by flipping conjugate pairs together (Appendix S4 D). Timescale separation is needed for the committor correspondence, not for the closed forms. Linearity is a genuine restriction: the estimator applies to nonlinear systems only within a neighborhood whose radius must be established empirically, and its variance grows with the escape rate, limiting the usable horizons. With several unstable directions, the flip preserves the unstable left subspace but distinguishes no direction within it, leaving open what plays the role of the reaction coordinate. The framework is local by design, complementing global numerical coherent-set computations.
Acknowledgment
We thank Joshua L. Pughe-Sanford, Charles L. Epstein, and Ruoyu Shang for helpful discussions and comments, which helped improve this paper. The Flatiron Institute is a division of the Simons Foundation.
References
[1] Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering Chapman and Hall/CRC 2024
[2] Random Perturbations of Dynamical Systems Springer 2012 3
[3] An analytic framework for identifying finite-time coherent sets in time-dependent dynamical systems Physica D: Nonlinear Phenomena 2013 250 1–19
[4] Dynamic isoperimetry and the geometry of Lagrangian coherent structures Nonlinearity 2015 28 10 3587–3622
[5] Kernel methods for detecting coherent structures in dynamical data Chaos: An Interdisciplinary Journal of Nonlinear Science 2019 29 12
[6] A Collection of Mathematical Problems Interscience Publishers 1960 New York
[7] Finite approximation for the Frobenius-Perron operator. A solution to Ulam's conjecture Journal of Approximation theory 1976 17 2 177–186
[8] Robust FEM-based extraction of finite-time coherent sets using scattered, sparse, and incomplete trajectories SIAM Journal on Applied Dynamical Systems 2018 17 2 1891–1924
[9] A data–driven approximation of the koopman operator: Extending dynamic mode decomposition Journal of Nonlinear Science 2015 25 6 1307–1346
[10] Diffusion maps Applied and computational harmonic analysis 2006 21 1 5–30
[11] Reaction-rate theory: fifty years after Kramers Reviews of modern physics 1990 62 2 251
[12] Towards a theory of transition paths Journal of statistical physics 2006 123 3 503–523
[13] Transition path sampling: Throwing ropes over rough mountain passes, in the dark Annual review of physical chemistry 2002 53 1 291–318
[14] Neurons as Detectors of Coherent Sets in Sensory Dynamics Advances in Neural Information Processing Systems 2026 38 158904–158924
[15] Relations between two sets of variates Breakthroughs in statistics: methodology and distribution Springer 1992 162–190
[16] Variational approach for learning Markov processes from time series data Journal of Nonlinear Science 2020 30 1 23–66
[17] Conditional Brownian motion and the boundary limits of harmonic functions Bulletin de la Société mathématique de France 1957 85 431–458
[18] Algebraic riccati equations Clarendon press 1995
[19] Dynamical systems and complex networks: A Koopman operator perspective Journal of Physics: Complexity 2024 5 4 041001
[20] Matrix analysis Cambridge university press 2012
[21] Reverse-time diffusion equation models Stochastic Processes and their Applications 1982 12 3 313–326
[22] Score-Based Generative Modeling through Stochastic Differential Equations International Conference on Learning Representations 2021
[23] Markovian representation of stochastic processes by canonical variables SIAM journal on control 1975 13 1 162–173
[24] Algebraic connectivity of graphs Czechoslovak mathematical journal 1973 23 2 298–305
[25] Brownian motion in a field of force and the diffusion model of chemical reactions physica 1940 7 4 284–304
[26] Statistical theory of the decay of metastable states Annals of Physics 1969 54 2 258–275
[27] Differential equations and dynamical systems Springer Science & Business Media 2013 7
[28] One-dimensional reaction coordinates for diffusive activated rate processes in many dimensions The Journal of chemical physics 2005 122 1 014503
[29] Escaping free-energy minima Proceedings of the national academy of sciences 2002 99 20 12562–12566
[30] Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint Annual review of physical chemistry 2016 67 1 159–184
[31] Spectral gap optimization of order parameters for sampling complex molecular systems Proceedings of the National Academy of Sciences 2016 113 11 2839–2844
[32] Canonical correlations of past inputs and future outputs for linear stochastic systems Systems & control letters 1984 4 4 209–215
[33] Quasi-stationary distributions Markov chains, diffusions and dynamical systems 2013
[34] Resonances of chaotic dynamical systems Physical review letters 1986 56 5 405
[35] Dirac Kets, Gamow Vectors and Gel'fand Triplets: The Rigged Hilbert Space Formulation of Quantum Mechanics Lectures in Mathematical Physics at the University of Texas at Austin Springer 1989