Non-Intrusive LS Adjoint Shadowing

Updated 3 October 2025

Non-Intrusive Least Squares Adjoint Shadowing (NILSAS) is a technique that computes sensitivities of long-time averages in chaotic systems by utilizing a constrained least squares approach focused on unstable adjoint modes.
It segments trajectories and applies QR factorizations to ensure numerical stability while integrating with legacy solvers for high-dimensional applications.
Its computational cost scales with the number of unstable adjoint modes rather than full state dimensions, enabling efficient multi-parameter sensitivity analysis in complex turbulent flows.

Non-Intrusive Least Squares Adjoint Shadowing (NILSAS) is an adjoint-based algorithm for computing sensitivities of long-time averaged outputs in chaotic dynamical systems. Unlike conventional tangent/adjoint sensitivity methods, which fail in chaotic regimes due to the exponential divergence of linearized solutions along unstable directions, NILSAS constructs a bounded adjoint shadowing solution via a non-intrusive, constrained least squares approach focused on the unstable adjoint subspace. This allows for efficient, robust sensitivity estimation with computational cost essentially independent of the number of design parameters, making it particularly suitable for large-scale, high-dimensional, and legacy simulation codes.

1. Mathematical Framework of Adjoint Shadowing

In a typical chaotic system governed by

$\frac{du}{dt} = f(u, s), \quad u \in \mathbb{R}^m, \quad s \in \mathbb{R}^p,$

the long-time averaged quantity of interest is

$\overline{J} = \lim_{T \to \infty} \frac{1}{T} \int_0^T J(u(t), s)\, dt,$

where $J(u, s)$ is the instantaneous observable. Naïve application of tangent or adjoint sensitivity methods to $\overline{J}$ fails in chaotic regimes, producing unbounded derivatives due to sensitivity with respect to initial conditions along unstable manifolds.

NILSAS addresses this by computing an adjoint shadowing direction $v^\infty(t)$ satisfying

$\frac{dv}{dt} + f_u(u(t), s)^{\mathsf{T}} v = -J_u(u(t), s),$

with the crucial condition that $v^\infty(t)$ is bounded and its unstable adjoint components vanish. An additional constraint,

$\lim_{T \to \infty} \frac{1}{T} \int_0^T \langle v^\infty(t), f(u(t), s) \rangle \, dt = 0,$

removes the non-uniqueness arising from the neutral (flow-tangent) direction (Ni, 2018). Existence, uniqueness, and boundedness under hyperbolicity have been rigorously established (Ni, 2018, Thakur et al., 13 Feb 2025).

The sensitivity formula in NILSAS is then

$\frac{d\overline{J}}{ds} = \lim_{T \to \infty} \frac{1}{T} \int_0^T \left[ \langle v^\infty(t), f_s(u(t), s) \rangle + J_s(u(t), s) \right] dt,$

where $f_s = \partial f/\partial s$ and $J_s = \partial J/\partial s$ .

2. Non-Intrusive Least Squares Formulation

NILSAS computes the adjoint shadowing direction by representing any inhomogeneous adjoint solution as a linear combination: $v(t) = v^*(t) + W(t) a,$ where $v^*(t)$ is a particular inhomogeneous solution and $W(t)$ are homogeneous adjoint solutions spanning the low-dimensional unstable adjoint subspace (Ni et al., 2018). The coefficient vector $a$ is determined by the constrained least squares problem: $\min_{a \in \mathbb{R}^M} \frac{1}{2} \int_0^T \| v^*(t) + W(t) a \|^2 dt,$ subject to

$\int_0^T (v^*(t) + W(t) a)^{\mathsf{T}} f(u(t), s) dt = 0.$

This restriction to the unstable adjoint subspace is key— $M$ is chosen slightly above the number of positive Lyapunov exponents (typically $M \ll m$ ), making the solution tractable for high-dimensional systems. The subspace is identified by backward-integrating homogeneous adjoint equations and applying periodic QR factorizations to maintain numerical conditioning (Ni et al., 2018, Ni et al., 2016).

Segmenting the trajectory, e.g., $[0, T]$ split into $K$ intervals, further improves stability. On each segment, constrained minimization for $a_i$ (the coefficients for that segment) involves only low-dimensional block matrices, and continuity across segments is enforced via recursion relations involving QR factors (Ni et al., 2018).

3. Implementation and Computational Advantages

Because the NILSAS algorithm operates by wrapping existing adjoint solvers with minimal modifications—primarily to supply arbitrary terminal conditions and record homogeneous adjoint trajectories—it is non-intrusive by design. The main steps are:

Integrate the primal system until the trajectory is on the attractor.
Divide the time domain into segments; on each, integrate the inhomogeneous and $M$ homogeneous adjoint equations backward in time.
At each segment interface, perform QR factorization to identify the local unstable adjoint directions, rescale the basis, and ensure continuity in the representation (Ni et al., 2018).
Assemble and solve the resulting low-dimensional coupled system for the coefficient vectors $\{a_i\}$ .
Once the continuous adjoint shadowing direction $v(t)$ is constructed, compute the sensitivity via quadrature using the formula above for all parameters simultaneously.

The computational cost of NILSAS scales with the number of unstable adjoint modes $M$ rather than the full state dimension $m$ , and is almost independent of the number of design parameters. The cost of additional parameter sensitivities is negligible, as $v(t)$ is reused and only extra inner products need to be evaluated (Ni et al., 2018, Blonigan et al., 2017).

4. Theoretical Foundations and Convergence

NILSAS rests on the mathematical framework of differentiability of SRB measures—ergodic averages in uniformly hyperbolic systems—established in the foundational work of Ruelle and subsequent shadowing theory. The existence and uniqueness of bounded shadowing solutions for the adjoint equation under hyperbolicity, and their ability to deliver the correct linear response, have been rigorously proven (Ni, 2018, Thakur et al., 13 Feb 2025).

The shadowing-based approach allows for robust computation of $d\overline{J}/ds$ with error that decays as $\mathcal{O}(1/\sqrt{T})$ as $T \to \infty$ , insensitive to the choice of adjoint boundary condition as long as it is bounded (Thakur et al., 13 Feb 2025). These properties are maintained even when the underlying system exhibits only approximate or nonuniform hyperbolicity, provided tangencies between CLVs are rare or controlled (Ni, 2017).

5. Comparison with Other Shadowing-Based Sensitivity Methods

NILSAS is closely related to tangent non-intrusive least squares shadowing (NILSS) (Ni et al., 2016) and finite-difference NILSS (FD-NILSS) (Ni et al., 2017), which operate in the tangent rather than adjoint formalism. FD-NILSS uses multiple primal solves to approximate tangent solutions, while NILSAS leverages the adjoint (dual) perspective. The key advantages of NILSAS over these tangent-based approaches include:

Independence of computational cost with respect to the number of parameters: The same constructed adjoint shadowing direction $v(t)$ may be paired with each $f_{s_j}$ .
Suitability for high-dimensional problems: The algorithm's effort scales with the dimension of the unstable adjoint subspace, which is often far smaller than the total system dimension in chaotic but physically dissipative systems (e.g., resolved fluid dynamics).
Integration with legacy codes: Only minor modification of adjoint solvers is required, facilitating coupling with industrial or legacy CFD frameworks (Ni et al., 2018, Blonigan et al., 2014).

Compared to stabilized-march adjoint shadowing (Thakur et al., 1 May 2025), NILSAS solves a global least-squares or KKT system rather than a sequence of triangular systems via backward substitution. The stabilized-march approach may thus have lower cost per segment when $n_u$ (number of unstable adjoint modes) is moderate and the chain of segments is long.

6. Applications and Demonstrations

NILSAS has been demonstrated on:

The Lorenz 63 system: For a long-time averaged observable $J(u) = z$ , NILSAS yields sensitivities $\partial\overline{J}/\partial\rho$ in close agreement with expected theoretical results, with variance decreasing as $T^{-0.5}$ (Ni et al., 2018).
Compressible turbulent cylinder flow: NILSAS delivered accurate sensitivities of mean drag to inflow Mach number, matching linear regression results computed from ensemble parameter variations, despite the underlying system state dimension being $O(10^6)$ and unstable adjoint dimension being $O(10)$ (Ni et al., 2018, Ni, 2017).
Large-scale CFD: The non-intrusive wrapper paradigm enables application to existing, highly optimized codes for large-eddy simulation or DNS (Blonigan et al., 2014, Blonigan et al., 2017).

The computational and memory overhead is significantly lower than that of conventional LSS or adjoint LSS formulations that do not focus exclusively on the unstable subspace.

7. Extensions and Theoretical Developments

The formalism underlying NILSAS extends naturally to spaces of hyperbolic diffeomorphisms as well as continuous-time flows (Ni, 2018). The adjoint shadowing operator $\mathcal{S}$ provides a rigorous means of isolating the shadowing contribution to the linear response. The overall linear response formula can be decomposed as

$\delta \rho(\Phi) = \text{SC} - \text{UC},$

where the shadowing contribution (SC) is captured by adjoint shadowing, and the unstable contribution (UC) accounts for the change in the SRB measure under parameter perturbation (Ni, 2022). In systems with few unstable directions, UC is often negligible and the NILSAS formula provides a complete linear response.

The method is robust to the violation of strict uniform hyperbolicity, as in open or weakly turbulent flows, provided that the adjoint shadowing solution remains bounded and the unstable/neutral subspaces are sufficiently separated (Ni, 2017). Windowed shadowing formulations (Chater et al., 2016), frequency-domain approaches (Kantarakias et al., 2022), and randomized descent algorithms for large-scale least squares (Lorenz et al., 2023) further enhance the scalability of shadowing-based adjoint sensitivity analysis.

In summary, Non-Intrusive Least Squares Adjoint Shadowing (NILSAS) provides a mathematically rigorous, computationally efficient, and implementation-friendly framework for computing sensitivities of time-averaged quantities in chaotic dynamical systems. By focusing on the unstable adjoint subspace and employing a constrained least-squares construction, it overcomes the limitations of classical tangent/adjoint methods and generalizes robustly to large-scale and legacy scientific computing environments (Ni et al., 2018, Ni, 2018, Ni, 2017).