Adjoint sensitivity analysis on chaotic dynamical systems by Non-Intrusive Least Squares Adjoint Shadowing (NILSAS)
Angxiu Ni Chaitanya Talnikar
physics.comp-ph nlin.CD









Angxiu]Angxiu Nicor
[cor]Corresponding author.
niangxiu@gmail.com
[url]https://math.berkeley.edu/ niangxiu/
[Angxiu]Department of Mathematics, University of California, Berkeley, CA 94720, USA

Chai]Chaitanya Talnikar
chaitukca@gmail.com
[Chai]Nvidia Corporation, Santa Clara, CA 95051, USA


We develop the NILSAS algorithm, 
which performs adjoint sensitivity analysis of chaotic systems via computing the adjoint shadowing direction.
NILSAS constrains its minimization to the adjoint unstable subspace, and can be implemented with little modification to existing adjoint solvers.
The computational cost of NILSAS is independent of the number of parameters.
We demonstrate NILSAS on the Lorenz 63 system and a weakly turbulent three-dimensional flow over a cylinder.




Chaos Adjoint methods Shadowing methods Sensitivity analysis 






§ INTRODUCTION



Sensitivity analysis is a powerful tool in helping scientists and engineers design products <cit.>,
control processes and systems <cit.>, solve inverse problems <cit.>, 
estimate simulation errors <cit.>, 
assimilate measurement data <cit.> 
and quantify uncertainties <cit.>.


Conventional sensitivity analysis works well for systems with a stable fixed point or a periodic orbit.
However, when the system is chaotic and the design objective is a long-time-averaged quantity, 
conventional sensitivity methods, be it tangent or adjoint, fails to provide useful sensitivity information.
Chaotic dynamical systems are typically modeled as hyperbolic systems,
for which many sensitivity analysis methods have been developed, such as by Lea et al. in <cit.>,
by Abramov and Majda in <cit.>,
and by Lucarini et al. in <cit.>.


Sensitivity of hyperbolic dynamical systems can be computed via the shadowing direction.
Bowen <cit.> proved that for a trajectory of a uniform hyperbolic systems, if we introduce a small perturbation in the governing equation,
there exists a shadowing trajectory, which satisfies the perturbed equations, but still lies close to the base trajectory.
Later Pilyugin <cit.> gave a formula for the first order difference 
between the shadowing trajectory and the base trajectory.
In this paper we call such first order difference the shadowing direction.
The Least Squares Shadowing (LSS) approach, developed by Wang et al., <cit.>, 
computes such shadowing direction through an L^2 minimization, and use it for sensitivity analysis.


The Non-Intrusive Least Squares Shadowing method (NILSS) developed by Ni et al. <cit.>
finds a `non-intrusive' formulation of the LSS problem which allows constraining the minimization problem in LSS to the unstable subspace.
For many real-life problems, the dimension of unstable subspace is much lower than the dimension of the dynamical system, 
and NILSS can be thousands times faster than LSS. 
A major variant of NILSS is the Finite Difference NILSS (FD-NILSS) algorithm <cit.>,
whose implementation requires only primal solvers, but not tangent solvers.
FD-NILSS has been applied to several complicated flow problems 
<cit.> which were too expensive for previous sensitivity analysis methods.


The marginal cost for a new parameter in NILSS is only computing one extra inhomogeneous tangent solution.
Yet for cases where there are many parameters and only a few design objectives, an adjoint version of NILSS is desired.
A continuous adjoint version appeared in the first publication of NILSS <cit.>; 
however, this version of lacks the constraint on the neutral subspace, which will be explained in our current paper.
Blonigan <cit.> developed a discrete adjoint version of NILSS, 
which was later implemented by Chandramoorthy et al. <cit.> using automatic differentiation.
In comparison to this discrete adjoint NILSS, NILSAS does not require tangent solvers, and requires less modification to existing adjoint solvers,
and the simplicity of the formula of NILSAS should also give it more robustness and perhaps better convergence.


Recently, we defined the adjoint shadowing direction for both hyperbolic flows and diffeomorphisms,
which is a bounded inhomogeneous adjoint solution with several other properties <cit.>.
We showed that the adjoint shadowing direction exists uniquely on a given trajectory, and can be used for adjoint sensitivity analysis.
Adjoint shadowing direction is defined using only adjoint flows, giving us a chance to get rid of tangent solvers in our algorithm, and arrive at a neat formula.


This paper presents the Non-Intrusive Least Squares Adjoint Shadowing (NILSAS) algorithm, 
where we construct a least squares problem to approximate the adjoint shadowing direction and then compute adjoint sensitivity.
The `non-intrusive' formulation of NILSAS allows it be built using existing adjoint solvers,
and more importantly, it allows the minimization be constrained to the unstable adjoint subspace.
The main body of this paper will be about continuous dynamical systems, 
and we briefly discuss NILSAS for discrete systems in <ref>.


We organize the rest of this paper as follows.
First, we prepare our study by defining our problem and reviewing adjoint flows and adjoint shadowing directions; we also provide some intuitions to help understanding adjoint shadowing directions.
Then we derive the NILSAS algorithm.
Then we present a detailed procedure list for our algorithm and give several remarks on the algorithm.
Finally, we demonstrate NILSAS on the Lorenz 63 system and a weakly turbulent three-dimensional (3D) flow over a cylinder.




§ PREPARATIONS



The governing equation of a hyperbolic flow, which models a time-evolving chaotic process, is: 

    ut = f(u,s),    u(t=0) = u_0 .

This differential equation is called the primal system, and a solution u(t) is the primal solution.
Here f(u,s):^m×→^m is a smooth function, u∈^m is the state of the dynamical system, u_0 is the initial condition, and s∈ is the parameter.
For now we assume there is only one parameter;
and in section <ref> we will explain how we can compute sensitivities with respect to several parameters with almost no additional cost. 


In this paper, we assume there is only one objective, which is a long-time-averaged quantity.
To define it, we first let J(u,s):^m×→ be a continuous function that represents the instantaneous objective.
The objective is obtained by averaging over a semi-infinite trajectory:

    J:= lim_T→∞1/T J(u,s)dt  .

We assume the system has a global attractor <cit.>, hence J only depends on s.
Our goal of this paper is to develop an algorithm computing the sensitivity, dJ / d s, whose marginal cost for a new parameter is negligible.




 §.§ Adjoint flow



 
  A homogeneous adjoint solution w(t):→^m is a function which solves the homogeneous adjoint equation:
  
    wt + f_u ^T w =0,

  where ·^T is the matrix transpose.
  An inhomogeneous adjoint solution is a function v(t):→^m which solves:
  
    vt + f_u ^T v =g(t),

  where g(t):→^m is a vector-valued function of time.



In numerical implementations, we typically solve adjoint equations backward in time.
This is because, as shown in <cit.>, when solving backward in time, 
the dimension of the unstable adjoint subspace is the same as the unstable tangent subspace, which is typically much lower than m.
On the other hand, if we solve the adjoint equation forward in time,
the unstable subspace has much higher dimension, causing strong numerical instability.



  In this paper, an adjoint covariant Lyapunov vector (CLV) with adjoint Lyapunov exponent (LE) λ
  is a homogeneous adjoint solution (t) such that there is a constant C, for any t_1,t_2∈, 
  
    (t_1)≤ C e^λ(t_2-t_1)(t_2) .





Note that the time direction in the above definition is reversed: 
if the adjoint CLV grows exponentially backward in time, its exponent is positive.
Adjoint CLVs with positive exponents are called unstable, those with negative exponents are stable, 
and those with zero exponent is are neutral.
In this paper, we sort adjoint CLVs by descending order of their exponents.
  The earliest mention of adjoint CLVs was by Kuptsov and Parlitz in <cit.>.
  For the purpose of defining adjoint shadowing directions and deriving the NILSAS method, 
  we recently proved the existence of adjoint CLVs under the same assumptions of adjoint shadowing theorem, 
  and found some relation between CLVs and adjoint CLVs.



  Adjoint CLVs are homogeneous adjoint solutions whose norm grows exponentially, and the adjoint LE is measured backward in time.
  This is similar but also different from tangent CLVs, which are tangent solutions measured forward in time.
  The CLV structure for the adjoint flow is the same as the tangent flow.
  That is, the adjoint LE spectrum is the same as the tangent LE spectrum.
  Moreover, the subspace of CLVs with an exponent λ is perpendicular to the subspace of all adjoint CLVs with exponents not λ, and vice versa.
  If we can write the full set of CLVs as a matrix valued function of time, W(t), then W^-T(t), where ·^-T is the inverse of transpose, is a matrix whose columns are adjoint CLVs:
  readers can verify that W^-T(t) satisfies the properties listed above.
  Note that we do not know if CLVs and adjoint CLVs with the same exponent are perpendicular or parallel.



We assume our system is uniform hyperbolic and it has a bounded global attractor.
Definition of hyperbolicity can be found in most textbook on dynamical system such as <cit.>,
and readers may also refer to <cit.> for a definition using the same notation as this paper.
Uniform hyperbolicity requires that the tangent space can be split into stable subspace, unstable subspace, and a neutral subspace of dimension one.
Together with the boundedness of the attractor, 
we can show the angles between two subspaces of different sets of tangent CLVs are always larger than some positive angle.
Since the adjoint equations have the same structure as the tangent ones,
there is only one neutral adjoint CLV, and adjoint CLVs are always bounded away from each other <cit.>. 




 §.§ Adjoint shadowing directions



In <cit.>, the author defined adjoint shadowing directions, 
proved their unique existence on a given trajectory, 
and showed how to use them for adjoint sensitivity analysis.
We briefly restate the main results in this subsection.




  On a trajectory u(t) on the attractor, for t≥ 0,
  the adjoint shadowing direction ^∞: _+→^m is defined as a function with the following properties:
  
    
  * ^∞ solves the inhomogeneous adjoint equation:
      
    ^∞t + f_u^T ^∞ = - J_u  ,

      where subscripts are partial derivatives, that is, f_u = ∂ f /∂ u, J_u = ∂ J /∂ u.
    
  * ^∞(t=0) has zero component in the unstable adjoint subspace. 
    
  * ^∞(t) is bounded by a constant for all t∈_+.
    
  * The averaged inner-product of ^∞ and f is zero:
      
    ^∞, f_avg := lim_T→∞1/T∫_0^T ^∞(t), f(t) = 0  ,
 
      where ·,· is the inner-product on the Euclidean space.
  



  We remind readers to distinguish the three different kinds of adjoint solutions we mentioned: homogeneous adjoint solutions, inhomogeneous adjoint solutions and adjoint shadowing directions.
  Homogeneous adjoint solutions are different from inhomogeneous ones, since homogeneous adjoint equations must have zero right-hand-sides.
  The adjoint shadowing direction is an inhomogeneous adjoint solution, but not any inhomogeneous adjoint solution:
  it must in extra have three more properties listed in the definition.
  In fact, one way to view NILSAS is that we search the space of all inhomogeneous adjoint solutions to find one such that it mimics the other three properties.
  More specifically, we minimize the L^2 norm, and constrain the inner product with f: this derivation will be revealed in later sections.




  For a uniform hyperbolic system with a global compact attractor,
  on a trajectory on the attractor,
  there exists a unique adjoint shadowing direction.
  Further, we have the adjoint sensitivity formula:
  
    Js = lim_T→∞1/T^∞, f_s + J_s   dt  .



  We explain the assumption of the adjoint shadowing theorem.
  First, if a dynamical system has a compact global attractor, it means there is a bounded set of states, or the attractor, 
  such that no matter what initial condition the system starts from, the trajectory will eventually enter the attractor and never leave.
  Second, uniform hyperbolicity here mainly means that there is only one neutral CLV.
  Third, by the compactness, the angles between all CLVs are larger than a positive angle, regardless of where we are on the attractor.

  Why do we make above assumptions in theories for shadowing methods?
  The main reason for assuming only one neutral CLV in shadowing methods is to prevent linear growth in inhomogeneous tangent/adjoint solutions.
  The main reason for global attractability is to ensure that shadowing trajectories are representative of the averaged behavior of the system.
  The main reason for compactness is because we want a bound for the projection operators projecting onto a particular subspace.
  Still, we remind readers that, in practice, shadowing methods may be effective beyond above assumptions, 
  as to be discussed in section <ref>.


Rather than giving an explicit expression of adjoint shadowing directions, which can be found as well in <cit.>,
the definition is stated as a criterion,
where we check several properties to determine if a function is indeed the adjoint shadowing direction.
In fact, we forged this definition for designing the NILSAS algorithm, which will be revealed in the next section.




 §.§ Interpreting adjoint shadowing directions



We give several different perspectives to help readers build intuitions on adjoint shadowing directions.
To start with, we revisit some formal descriptions of shadowing operators.
The shadowing operator, denoted as 𝒮, 
can be viewed roughly as mapping a vector-valued function f_s(t) to another vector-valued function v^∞(t),
Here f_s(t) is the perturbation on f due to parameter perturbations;
v^∞ is the (tangent) shadowing direction, 
which is first order approximation of the difference between the shadowing and the base trajectory.
Note that 𝒮 is a linear operator, and both f_s(t) and v^∞(t) are linear approximations.
If we neglect the subtleties due to the neutral CLV, we have roughly

    Js≈v^∞, J_u_avg  =  𝒮(f_s), J_u_avg .

Where ·, ·_avg is an inner product.

  
First we provide a utility point of view for adjoint shadowing directions, which is also an algebraic point of view.
Riesz's representation theorem tells us that there is an adjoint operator  such that 

    𝒮(f_s), J_u_avg = f_s, (J_u)_avg .

The adjoint shadowing direction, ^∞, can be viewed as (J_u).
Suppose now that we have a computer program which approximately functions as .
If we have two parameters s_1, s_2, then they can perturb the governing equation by f_s_1 and f_s_2, whereas J_u keeps the same.
This means that we only need to run our adjoint program once, and use the result to inner-product with both f_s_1 and f_s_2.
For cases where there are many parameters s and a few objectives J, we only need to run our program a few times.


Then a physical point of view.
Formally dJ_avg/ds = ∂ f /∂ s, ^∞,
and we can formally cancel ds on both side and divide by df.
Thus we get ^∞ = ∂ J_avg/∂ f.
This means that the adjoint shadowing direction can be viewed as how perturbations in the output of f affects the objective.
For example, in our fluid examples, if the ρ E component of ^∞ has value β for all u, 
it means that if we can somehow keep infusing unit energy into the flow, the objective will be changed by β.
Note that in ^∞ we do not prescribe how perturbations in f is generated from changing a particular parameter; the dependence f on s is described in another term, f_s.


Then an implementation point of view.
The adjoint operator of a matrix is simply its transpose.
For many other cases, we can see that matrix transposition often appears as a crucial step in the formulation of adjoint operators.
So once we derive an new adjoint formula of something, we may ask if it can be presented as neatly as transposing a matrix.
Adjoint shadowing directions satisfy inhomogeneous adjoint equations, which is a linear ODE whose matrix is the transpose of the Jacobian matrix.
Hence there is chance that algorithms, such as NILSAS, do not differ too much from existing adjoint solvers.
In fact, giving a neat recipe for adjoint shadowing direction is the main contribution of both the adjoint shadowing theorem and the NILSAS algorithm,
since the existence is already given by Riesz's representation theorem.


Above different views on adjoint shadowing directions are the best intuitions we can come up with now.
Still, we acknowledge that it is not yet as intuitive as shadowing direction,
which can be viewed as vectors `starting from' a base trajectory, and `pointing to' a shadowing trajectory, 
with parameter s+Δ s, that truly exists.
Moreover, our intuition is not quite refined yet, as we still need to draw upon the abstract notion of duality with their tangent counterparts,
in order to describe the subtle difference among all those adjoint solutions: 
homogeneous adjoint solutions, inhomogeneous adjoint solutions and adjoint shadowing directions.
We humbly ask readers to inform us better intuitions.




§ DERIVING NILSAS





 §.§ The non-intrusive formulation




On a finite trajectory of time span [0,T], the NILSAS algorithm computes a  which approximates ^∞.
Since the definition of the adjoint shadowing direction is similar to the tangent shadowing direction,
it is not surprising that NILSAS has similar formulation as NILSS.

In NILSAS, we strictly enforce the first property of adjoint shadowing directions by constraining our solutions to inhomogeneous adjoint solutions.
The second property is changed to a symmetric statement that stable component in (T) should be O(1), 
which can be easily satisfied.
The third property is approximated by minimizing the L^2 norm of the inhomogeneous adjoint solution.
The fourth property is strictly enforced by adding a constraint to our minimization problem.
In this subsection, we explain why the  given by this reverse-engineering is a good approximation of ^∞.



Our algorithm strictly enforces the first property of ^∞.
To do this, we represent the solution set of equation (<ref>) as a particular solution plus the space of homogeneous solutions.
We select the particular solution as the conventional inhomogeneous adjoint solution ^*, which is defined as the solution of:

    ^*t + f_u ^T ^* =-J_u ,  ^*(T)=0 .

Then we select the collection of all adjoint CLVs, = [_1,⋯,_m], 
all of which have terminal condition _j(T)=1, as the basis of the space of homogeneous solutions.
Hence we can enforce the first property by considering candidates only in the following form for some {a_j}_j=1^m:

    = ^* + ∑_j=1^m a_j _j = ^* +  a  ,

where the coefficients a = [a_1,...,a_m]^T is a column vector.
Another interpretation of our way of enforcing the first property is that, 
we want to start  from ^*, and modify by adding adjoint CLVs to approximate ^∞;
in other words, the coefficients a should be such that Za ≈^∞ - ^*.
We then use other properties to determine the coefficients for stable, unstable and neutral CLVs. 



When defining adjoint shadowing directions in <cit.>, the author was considering functions defined starting from time zero,
whereas in our case here, adjoint solutions are solved from T backward in time.
Hence, in order to keep the sensitivity formula, we change the second property to a symmetric statement,
that is, we want the stable component in (T) be in the order O(1), meaning be bounded by a constant independent of T.
Equivalently, we require coefficients for stable CLVs be O(1).
Another way to interpret is that, if this O(1) condition is true, 
then since ^∞(T) is O(1), the stable component in (T)-^∞(T) is O(1); 
now since stable CLVs decay exponentially fast, their contribution in - ^∞ can be neglected.
This O(1) condition is a loose requirement and, as we will see, it can be easily satisfied.



To mimic the boundedness in the third property of adjoint shadowing directions, we minimize _L^2,
which determines the coefficients for the unstable CLVs. 
Indeed, this minimization removes significant unstable CLVs from -^∞, 
since otherwise this difference would grow exponentially, 
and since ^∞ is bounded,  would have large L^2 norm.



We strictly enforce the fourth property of .
This determines the coefficient for the neutral adjoint CLV, 
since as shown in <cit.>,
f is always orthogonal to non-neutral adjoint CLVs.
Note also that the norm of the neutral adjoint CLV is bounded, unlike neutral tangent CLV, which can have linear growth.
Hence we can allow its coefficient be O(1) without jeopardizing the boundedness property we used earlier.
In fact, the adjoint NILSS in <cit.> lacks exactly this constraint on the neutral adjoint CLV.



To summarize, we determine coefficients for unstable adjoint CLVs via a minimization,
the coefficient for the neutral adjoint CLV via equation (<ref>),
and we do not care coefficients for stable adjoint CLVs too much.
Hence there is no need to provide stable CLVs to our algorithm;
it is even unnecessary to provide accurate non-stable CLVs, they can contain some stable components at T.
Further, we care not individual CLVs but only their span.
Hence we can replace  by =[_1,⋯,_M], 
with M≥+1,  being the number of unstable CLVs,
and {_j}_j=1^M are homogeneous adjoint solutions whose non-stable components at T span the entire non-stable subspace.
Such a set of solutions can be obtained by solving homogeneous adjoint equations from almost all terminal conditions of M randomized unit vectors. 



With above discussions, we see that our algorithm should solve the NILSAS problem on one segment:

    min_a∈^M    1/2^* +  a, ^* +  a ,
       s.t. ∫_0^T ^* +  a, f = 0  .

This is simply a least squares problem with arguments a∈^M.
Note that all adjoint solutions can be computed after making little modifications to existing adjoint solvers, 
and our minimization is constrained to essentially only the unstable adjoint subspace:
these are the benefits of the non-intrusive formulation. 
Letting = + a, we can compute sensitivity via equation (<ref>) on a finite trajectory:

    Js≈1/T, f_s + J_s   dt  .







 §.§ Dividing trajectory into segments




An issue in numerical stability is that, as the trajectory gets longer, 
all adjoint solutions become dominated by the fastest growing adjoint CLV; 
as a result, the minimization problem in equation (<ref>) becomes ill-conditioned.
This issue also happened in NILSS <cit.> and in the algorithm for computing CLVs <cit.>,
and here we use a similar technique to resolve it, that is, dividing the whole trajectory into multiple segments, and rescaling at interfaces.



Roughly speaking, at the end of each segment,
we orthogonalize and rescale adjoint solutions  and ^* so that they are no longer dominated by the first CLV.
Note here since adjoint solutions are integrated backward in time, the initial condition is at the end of a segment.
Despite that now  and ^* are discontinuous across segments, we can still construct  as their linear combinations,
and keep  continuous across all segments.
This continuous  computed from multiple segments should be identical to that solved on one large segment containing the entire trajectory.



We first define some notations, as shown in figure <ref>.
Let T be the time length of the entire trajectory, and K the total number of segments.
We denote the time span of the i-th segment  by [t_i, t_i+1], where t_0=0, t_K = T.
For quantities defined on a entire segment such as _i, ^*_i, C_i, d_i and a_i,
their subscripts are the same as the segment they are defined on.
For quantities defined only at the interfaces between segments such as Q_i, R_i, b_i, p_i and λ_i,
their subscripts are the same as the time point they are defined at.
Some of the notations are used immediately below, the others will be used in section <ref> and <ref>.






At time t_i, we perform QR factorization to _i(t) = [_i1(t), ⋯_iM(t)], 
which is a m× M matrix whose column vectors are homogeneous adjoints on segment i.
We use the Q-matrix, the matrix with orthonormal columns, as the terminal condition for _i-1 on segment i-1.
More specifically, 

    _i(t_i) = Q_i R_i,   and  _i-1(t_i) = Q_i  .




At time t_i, we also rescale ^*_i(t), which is the particular inhomogeneous adjoint solution on segment i.
To do this, we subtract from ^*_i its orthogonal projection onto homogeneous adjoint solutions.
More specifically, 

    p_i := ^*_i(t_i) - Q_i b_i  ,
       where  
      b_i = Q_i ^T ^*_i(t_i)  ,
       and ^*_i-1(t_i)=p_i  .

This rescaling maintains the continuity of the affine space _i^*+span{_ij}_j=1^M across different segments.



The continuity of affine space allows us to impose continuity condition for _i, which is the adjoint shadowing direction on segment i.
On each segment,  _i = ^*_i + _i a_i for some a_i∈^M.
The continuity condition can now be expressed via a relation between a_i and a_i-1:

    ^*_i(t_i) + _i(t_i) a_i = ^*_i-1(t_i) + _i-1(t_i) a_i-1 .

Apply equation (<ref>) and (<ref>), cancel ^*_i(t_i) on each side, we get:

    Q_i R_i a_i = -Q_i b_i + Q_i a_i-1

Since Q_i has orthonormal columns, Q_i^TQ_i=I∈^M× M. 
Multiplying Q_i^T to the left of both sides, we have the continuity condition for :

    a_i-1 = R_i a_i + b_i  .






§ THE NILSAS ALGORITHM
 




 §.§ Procedure list of the algorithm




Now we give a procedure list of the NILSAS algorithm.
To start with, we need to have an inhomogeneous adjoint solver and a homogeneous adjoint solver, both can take arbitrary terminal conditions.
The inhomogeneous adjoint equation we solve in NILSAS has right-hand-side -J_u, which is the same as many existing adjoint solvers.
Hence for inhomogeneous adjoint solvers in NILSAS, we only need to change existing solvers to be able to take arbitrary terminal conditions.
For homogeneous adjoint solvers, we only need to further change the right-hand-side of existing solvers to zero.



We provide the following data to NILSAS: 
1) the number of homogeneous adjoint solutions, M≥+1, where  is the number of unstable CLVs; 
2) the total number of segments, K; 
3) for convenience, we assume that the length of all time segments are the same, denoted by Δ T.
The  total time length is determined by T =K Δ T.
Moreover, in the procedure list below, inner products are written in matrix notations, and by default vectors are in column forms.




  
  * Integrate the primal system for sufficiently long time before t=0 so that u(t=0) is on the attractor.

  
  * Compute the trajectory u(t), t∈ [0,T], by integrating the primal system.

  
  * Generate terminal conditions for _i and ^*_i on the last segment i=K-1:
  
    
  * Randomly generate a m× M full rank matrix, Q'.
      Perform QR factorization: Q_K R_K = Q'.
    
  * Set p_K = 0.
  

  
  * Compute _i and v^*_i on all segments.
  For i=K-1 to i=0 do:
  
    
  * To get _i(t), whose columns are homogeneous adjoint solutions on segment i, solve:
      
    _it + f_u ^T _i = 0,   _i(t_i+1) = Q_i+1 .

      To get ^*_i(t), solve the inhomogeneous adjoint equation:
      
    ^*_it + f_u ^T ^*_i = -J_u  ,  ^*_i (t_i+1) = p_i+1 .


    
  * Compute the following integrations.
      
    C_i = ∫_t_i^t_i+1_i^T _i dt  ,   
            d^wv^*_i = ∫_t_i^t_i+1_i^T ^* dt  , 
       d^wf_i = ∫_t_i^t_i+1_i^T f dt  ,   
            d^v^*f_i = ∫_t_i^t_i+1^*T f dt  , 
       d^wf_s_i = ∫_t_i^t_i+1_i^T f_s dt  ,   
            d^v^*f_s_i = ∫_t_i^t_i+1^*T f_s dt  , 
       d^J_s_i = ∫_t_i^t_i+1 J_s dt  ,

      where d^wv^*_i, d^wf_i, d^wf_s_i∈^M; d^v^*f_i, d^v^*f_s_i, d^J_s_i∈;
      C_i∈^M× M is the covariant matrix.

    
  * Orthonormalize homogeneous adjoint solutions via QR factorization:
      
    Q_i R_i = _i(t_i)


    
  * Rescale the inhomogeneous adjoint solution using Q_i: 
      
    p_i = ^*_i(t_i) - Q_i b_i  ,  where  b_i = Q_i^T ^*_i(t_i)  .

  


  * Compute the adjoint shadowing direction {_i}_i=0^K-1.
  
    
  * Solve the NILSAS problem on multiple segments:
      
    min_a_0,⋯,a_K-1∈^M∑_i=0^K-11/2 (a_i)^T C_i a_i +  (d^wv^*_i)^T a_i,   
         a)   a_i-1 = R_i a_i + b_i  ,     i=1,⋯,K-1  ,
         b)  ∑_i=0^K-1 (d^wf_i)^T a_i + ∑_i=0^K-1 d^v^*f_i = 0  .

      This is a least squares problem in {a_i}_i=0^K-1⊂^M.
      In <ref> we suggest one way to solve this problem.

    
  * On each time segment i, _i is given by
      
    _i(t) = ^*_i(t) + _i(t) a_i .

  

  
  * Compute the derivative by:
    
    dJ/ds≈1/T∑_i=0^K-1∫_t_i^t_i+1( _i^T f_s  +J_s ) dt
          = 1/T∑_i=0^K-1( d^v^*f_s_i + a_i^T d^wf_s_i + d^J_s_i )






 §.§ Remarks about NILSAS





  §.§.§ Miscellaneous




We remark that if we are not interested in obtaining (t) for all t, 
there is no need to store adjoint solutions  and ^* in computers, which is typically lots of data.
To compute sensitivity, we only need to store d^wv^*_i, d^wf_i, d^wf_s_i,
d^v^*f_i, d^v^*f_s_i, d^J_s_i, C_i given in equation (<ref>),
R_i given in equation (<ref>), and b_i given in equation (<ref>).



Moreover, when computing quantities in equation (<ref>), 
we can estimate the integration using particular values of the integrands evaluated at several snapshots,
to further reduce the storage management cost.
For example, similar to <cit.>, we can estimate C_i by the terminal value of W_i, which is Q_i,
and thus C_i=I∈^M× M.
We suggest further research be done to determine which estimation is best practice.



NILSAS has the benefit of typical adjoint algorithms,
that is, for a new parameter s,  does not change, 
so we only need to give new f_s, J_s, and recompute equation (<ref>) and (<ref>).
Hence the extra cost for a new parameter is only performing an L^2 inner product, which is negligible in comparison with the total cost of the algorithm.
More specifically, assume there are n parameters s=[s_1,⋯,s_n], we can define 

    J s    = [Js_1, ⋯, Js_n ] ∈ R^1× n ;
    
      J_s    = [ J s_1, ⋯,  J s_n]  ∈ R^1× n ;
    
      f_s    = [ f s_1, ⋯,  f s_n]  ∈ R^m× n .

With these definitions, the NILSAS algorithm, in particular equation (<ref>) and (<ref>),
extend to several parameters with almost no extra cost.
An extreme example is where f_s is unknown a priori and we can now use  to design an optimal control, f_s.


  
  The assumptions in theorem <ref> are made for theoretically proving the unique existence of adjoint shadowing directions and convergence of NILSAS.
  In practice, it is possible that NILSS/NILSAS are still valid on a chaotic system which fails these assumptions.
  For example, the 3D cylinder flow we investigate later in this paper has at least two neutral CLVs, 
  corresponding to translations in time and in the span-wise directions, due to the periodic boundary condition.
  In fact, 
  in <cit.> we also showed that the smallest angle between tangent CLVs depends on meshes and may fall below a threshold value:
  this further violates our assumptions.
  However, we did found the trend that angles between tangent CLVs gets larger when their indices are further apart: 
  this property is related to hyperbolicity, but has not been well investigated yet.
  As we shall see, both NILSS and NILSAS compute correct sensitivities on this 3D flow.
  The generality of shadowing methods is as suggested by the chaotic hypothesis <cit.>,
  that is, theoretical tools may still valid for non-uniform hyperbolic chaotic systems, 
  even though those tools can only be rigorously proved with a stricter assumption.
  We do not expect NILSS and NILSAS be valid for all chaotic systems;
  however, they are valid somewhere beyond our current assumptions.
  
  We call for more research to identify the limit of shadowing methods, especially in real-life problems.



In <ref>, we discuss in detail NILSAS for discrete systems, more specifically, hyperbolic diffeomorphisms.
Adjoint shadowing directions for diffeomorphisms were also defined in <cit.>.
Because of the absence of the neutral subspace, the NILSAS algorithm for hyperbolic diffeomorphisms is easier than flows.
To obtain NILSAS for diffeomorphisms, we no longer compute d^wf_i, d^v^*f_i, 
and no longer impose the second constraint in the NILSAS problem.
Of course, we should change integrations to summations, and adjoint equations to their discrete counterparts.




  §.§.§ Number of homogeneous adjoint solutions



Since the number of homogeneous adjoint solutions should be strictly larger than then number of unstable adjoint CLVs, 
which equals the number of unstable tangent CLVs, 
we first discuss the number of unstable CLVs, about which there are two questions: 
(1)  whether the absolute number can be large; (2) whether the number is significantly lower than the dimension of the system.
We are interested in (1) because we want to estimate the cost of NILSAS.
We are interested in (2) because we want to determine whether computational efficiency can benefit from the non-intrusive formulation, which restricts minimization to unstable subspaces.
Roughly, the efficiency improvement due to the non-intrusive formulation is proportional to the ratio of the system dimension to unstable subspace dimension.


First, the absolute number of unstable CLVs can be large and the cost of NILSAS increases.
We think maybe this is the price to pay for chaos, 
that is, for more chaotic systems, numerical methods should be more expensive, not only for NILSS/NILSAS, 
but also for other common methods such as computing long-time averages, which should take longer time to converge for more chaotic systems.
Second, in a recent paper <cit.>, based on observations on flow past a 3D cylinder, we conjectured that for open flows, CLVs active in the freestream or less turbulent regions are stable.
At least for these open flows, where there are large areas of freestream, a large fraction of CLVs should be stable.
For these cases, we can benefit from the non-intrusive formulation by restricting minimization to unstable subspaces.
On the other hand, there are cases where significant part of all CLVs are unstable, such as Hamiltonian systems.
For such cases the non-intrusive formulation does not help reduce the computational cost.
The other possible case in fluid mechanics might be wall bounded turbulence,
where for higher Reynolds number there are more unstable CLVs, but also more CLVs since finer meshes should be used, 
and we do not know yet if, in the limit, a significant fraction of CLVs can be unstable.


We provide some examples in computational fluid mechanics on the number of unstable CLVs and the dimension of the system.
For a 2D incompressible channel flow over a backward facing step at Reynolds number Re=2.5× 10^4, 
there are 13 unstable CLVs in a system of dimension 4× 10^4 <cit.>.
For a 2D NACA 0012 airfoil at Mach number Ma=0.2, angle of attack 20, Re=2400,
there are less than 5 unstable CLVs for different implementations with system dimension
ranging from 7× 10^3 to 8× 10^5 <cit.>.
For a 3D turbulent channel flow with Ma=0.3 and Re_τ = 180 on a domain of size 4π× 2× 2π, 
there are about 1.5× 10^3 unstable CLVs out of a system of dimension 2.2× 10^6 <cit.>.
For a 3D weakly turbulent flow over a cylinder at Re=5.2× 10^2, 
there are 20 unstable out of 1.9× 10^6 <cit.>.
To conclude, we believe that although the cost of NILSS and NILSAS can be high, 
for many cases, the non-intrusive formulation is beneficial for computational efficiency.


In our procedure list we listed M in the setting of NILSAS and required it strictly larger than m_us, the number of unstable CLVs.
How should we know m_us before running NILSAS? And what if we chose initial M smaller than required?
First, the number of unstable modes is roughly positively related to how chaotic the flow is.
This is not a rigorous criterion but readers can look at some test cases to have a rough sense.
But there is not any precise method that allows us to know the exact number at the first glance.
Second, even if we started with an insufficient M, we can add adjoint solutions inductively in NILSAS,
rather than recomputing everything all over again.
More specifically, in the NILSAS problem in equation (<ref>) 
and the sensitivity formula in equation (<ref>),
say we want to add k more adjoint solutions,
then coefficients arrays d^wv^*_i, d^wf_i, d^wf_s_i and b_i,
should be augmented by k more entries, 
while the old coefficient arrays are not changed inside the new arrays;
similarly, the coefficient matrices C_i, R_i should be augmented by k rows and k columns.


The headache of choosing an initial M is further relieved by the fact that adjoint solutions can be more efficiently computed in batches.
Within each segment, we can accelerate NILSAS by taking advantage of the fact that all adjoint solutions, 
both homogeneous and inhomogeneous, use the same Jacobian f_u.
If the numerical integration is vectorized, 
we can integrate all adjoint solutions simultaneously without repeatedly loading f_u into the computer CPU, 
which is the most time-consuming procedure in the numerical integration.
At each time step, instead of several matrix-vector products, we can perform one matrix-matrix products,
where the second matrix is composed of several adjoint solutions;
then we add the right-hand-side to the inhomogeneous adjoint solution.
For example, for a 4th order IEDG solver, the marginal cost for one more adjoint solution can be only, say 0.037, 
of the first adjoint solution.
In this scenario we should start NILSAS and then add adjoint solutions by batches on the order of 1/0.037≈ 27 adjoint solutions per batch.
This should be further faster than adding adjoint solutions one by one.
[The ideas of taking advantage of vectorized integration 
  and the estimation on IEDG were both given during private discussion by Pablo Fernandez, 
  who co-authored with us on finite difference NILSS (FD-NILSS), see arXiv:1711.06633.]




  §.§.§ Other settings of NILSAS



It is required by the algorithm that we run primal system long enough before the main part of NILSAS, so that our initial condition is on the attractor.
In general, we can not know very well what is `long enough' before we do any computations, 
and this run-up time is determined empirically as the time when the flow field starts to repeat itself.
In a typical scenario, we would run a primal simulation before taking interest in any sensitivities.
When running that primal simulation, there is the same question of when we reach the stage that enough long-time behavior has been captured:
typically this is indicated by that several objective functions began to oscillate around some averaged values.


We should also determine the time length T on which we run NILSAS.
In practice T is determined empirically as the time when the sensitivity computed by NILSAS converges to within the uncertainty bound we desire.
However, there is one caveat that the adjoint solutions are computed backwards in time.
Now if we find T insufficient, we can not add time after T without recomputing all adjoint solutions, since integration adjoint solutions forward in time is very unstable.
Rather, we should add time before our current trajectory.
In practice, we should run our primal simulation till enough long-time behavior has been captured, then start computing adjoint solutions from the end of that primal trajectory.
We have found that typically NILSAS requires a shorter trajectory to compute sensitivity than that required to reflect average behavior.


Then we discuss the choice of segment length Δ T.
Similar to NILSS, Δ T is determined by that within one segment, the leading adjoint CLV does not dominate the M-th adjoint CLV.
This is because otherwise we would have covariant matrices, {C_i}_i=0^K-1, with small condition number, 
which would lead to eventually the poor condition of the NILSAS problem in equation (<ref>).
We recommend Δ T (λ_1-λ_M) to be O(1), in which case within one segment,
the leading CLV would grow to be about e^1=2.7 time larger than the M-th CLV.


A related question is that if the leading LE is large and we select a small Δ T, will the cost of frequent rescaling offset the cost reduction due to non-intrusive formulation?
First, as we discussed in section <ref>, there are a lot of fluid systems whose CLVs are mostly stable, in which case the non-intrusive formulation is beneficial.
Second, our understanding is that the numerical methods should have smaller time steps for more chaotic systems, to capture accurate motions on all scales.
Hence one segment, although is shorter in physical time, may still contain many small time steps.
As a result, the rescaling may not be more frequent for more chaotic systems. 
Again, we call for more research on numerical schemes and LE spectrum, 
especially for systems other than fluid or extremely chaotic systems.




  §.§.§ Comparison with other shadowing algorithms




There are currently several variants of NILSS <cit.>,
such as the Finite-Difference NILSS (FD-NILSS) <cit.> and discrete adjoint NILSS <cit.>.
NILSAS, as well as these NILSS variants, bears part of the merit of `non-intrusive' formulation, that is, in comparison to LSS,
the minimization problems in these algorithms are constrained to the unstable subspaces.
Hence, for many real-life problems, where the unstable subspaces have significantly lower dimension than the dynamical systems,
these algorithms should be significantly faster than LSS.



We compare NILSAS with variants of NILSS in table <ref>.
In particular, we want to compare in more detail the two adjoint algorithms: discrete adjoint NILSS versus NILSAS.



NILSAS should be easier to implement than the discrete adjoint NILSS since it does not require tangent solvers, 
and requires less modification to existing adjoint solvers. 
Furthermore, unlike discrete adjoint NILSS, 
NILSAS does not explicitly depend on the fact that inner-products between adjoint and tangent homogeneous solutions are constants:
since this property holds true only for analytic solutions but is typically false for numerical solutions,
we think NILSAS should be more robust to implementations of tangent and adjoint solvers, and should typically have better convergence.
We suggest more numerical comparison be done to compare the two methods.







§ APPLICATIONS





 §.§ Application on Lorenz 63 system




In this subsection we apply NILSAS to the Lorenz 63 system as an illustration.
[The python code used for this section is at: https://github.com/niangxiu/nilsas.]
Lorenz 63 is an ordinary differential equations system with three states u = [x, y, z]:

    xt = σ(y-x),
        yt=x(ρ-z)-y,
        zt = xy- β z.

We fix β=8/3.
This system models the heat transfer in a fluid layer heated from below and cooled from above.
In particular, x is the convection rate, y the horizontal temperature variation, and z the vertical temperature variation.
The parameters σ and  ρ are proportional to the Prandtl number and Rayleigh number.
We select the instantaneous objective function as J(u) = z, and hence our objective J is the averaged vertical temperature variation.



The primal system and adjoint equations are integrated via the explicit time-stepping scheme:

    u_k+1   = u_k + f(u_k) Δ t 
    _k   = _k+1  + f_u(u_k)^T _k+1Δ t
    ^*_k   = ^*_k+1  + f_u(u_k)^T ^*_k+1Δ t + J_u(u_k) Δ t

where the subscript k denotes the time step number in numerical integration.
The time step size is Δ t = 0.001.
For NILSAS, time segment length is Δ T=0.2, thus there are 200 time steps per segment.



We want to determine the number of unstable CLVs for the Lorenz system.
The Lyapunov exponents, λ_1, λ_2, λ_3, satisfy the following constraints <cit.>:

    λ_1 + λ_2 + λ_3 = -(1+σ+β) < 0  .

Moreover, one of these exponents corresponds to the neutral CLV so it is zero.
Hence there is at most one positive exponent, so in NILSAS we set the number of homogeneous solutions M=2.



We verify that NILSAS gives correct sensitivities by computing J and ∂J/ ∂ρ for different ρ, while fixing σ=10.
The Lorenz system has one quasi-hyperbolic strange attractor when 25≤ρ <31, and one non-hyperbolic attractor when 31 ≤ρ≤ 50:
none of these cases strictly satisfies our uniform hyperbolic assumption.
As shown in figure <ref>, as ρ becomes larger, 
the system becomes non-hyperbolic, and the sensitivity results given by NILSAS begin to oscillate.
Nevertheless, NILSAS gives that ∂J/∂ρ is approximately 1 for all ρ, 
which matches the trend between J and ρ:
this again shows that NILSAS can be effective for systems not satisfying assumptions of theorem <ref>, 
as we discussed in section <ref>.






Then we show that both J and the sensitivities computed by NILSAS converge as the trajectory length T gets larger,
while fixing σ=10 and ρ=28.
Figure <ref> shows that the standard deviation of J reduces at the rate of T^-0.5.
Figure <ref> shows that the sensitivities computed by NILSAS, 
with respect to both ρ and σ, converge faster than the rate of T^-0.5.










NILSAS computes sensitivities with respect to multiple parameters with almost no additional cost, 
since the adjoint shadowing solution  does not depend on the choice of parameters.
Figure <ref> illustrates the contour of J with respect to ρ and σ, 
and the gradient, [∂J /∂ρ, ∂J /∂σ], 
is computed by NILSAS.
Since we use the same length unit for both parameters, gradients should be perpendicular to the level sets of the objective:
this is indeed the case, and it shows NILSAS gives correct gradient information.






Finally, we draw the norm of an adjoint shadowing direction in figure <ref>.
As we can see from the left plot,
the norm of the adjoint shadowing direction does not grow exponentially,
satisfying the third property of adjoint shadowing directions.
Moreover, as shown in the right plot,
the adjoint shadowing direction computed by NILSAS is continuous.
This shows that our dividing trajectory technique indeed allows us to recover a continuous adjoint shadowing direction.







 §.§ Application on a weakly turbulent flow past a three-dimensional cylinder




In this subsection, we apply NILSAS to a 3D subsonic flow over a cylinder at Reynolds number Re=1100 and Mach number Ma=0.093.
[The flow solver, adFVM, used for this section is at: https://github.com/chaitan3/adFVM, the particular file that implements the NILSAS algorithm used in this case is apps/nilsas.py.]
The flow-wise length of the domain is 60d, where d=0.25 mm  is the diameter of the cylinder. 
The Reynolds number is defined using the diameter of the cylinder and the density, velocity and viscosity of inflow. 
The span-wise extent, at z = 2d, is sufficient to capture most of the important flow features, like a turbulent wake and flow separation.
The front view of our fluid problem is shown in figure <ref>.






We use compressible Navier-Stokes equations with the ideal gas law approximating the thermodynamic state equation <cit.>.
The gas is assumed to be air.
More specifically, the governing equations are:

    ∂ρ/∂ t + ∇· (ρ𝐮) = 0  , 
       ∂ (ρ𝐮)/∂ t + ∇· (ρ𝐮𝐮) + ∇ p = ∇·σ +   , 
       ∂ (ρ E)/∂ t + ∇· (ρ E 𝐮 + p 𝐮) = ∇· (𝐮·σ + αργ∇ e)  ,
       σ = μ (∇𝐮 + ∇𝐮^T) - 2μ/3(∇·𝐮)𝐈 ,
           c = √(γ p/ρ) ,
       p = (γ - 1)ρ e  ,
           e = E - 𝐮·𝐮/2 .

Here ρ is the density, 𝐮 is the velocity vector, ρ E is the total energy,
p is pressure, e is internal energy of the fluid, c is the speed of sound,
γ = 1.4 is the isentropic expansion factor and
μ is the viscosity field modeled using Sutherland's law for air 

    μ = C_s T^3/2/T + T_s

where T_s = 110.4 K and C_s = 1.458 × 10^-6kg/m s √(K).
α is the thermal diffusivity modeled using

    α = μ/ρ Pr

where Pr = 0.71 is the Prandtl number.



We use an unstructured hexahedral mesh with approximately 7× 10^5 cells, with 50 cells in the span-wise direction.
The front view of our mesh is shown in figure <ref>.






We use a second order finite volume method (FVM) <cit.> for unstructured hexahedral meshes.
The central differencing scheme is used to interpolate cell averages of the flow solution onto faces of the mesh <cit.>. 
The numerical fluxes for the conservative flow variables are computed using the Roe approximate Riemann solver <cit.>. 
An explicit time integration scheme, the strong stability preserving third order Runge-Kutta method <cit.>, 
is used for time marching the numerical flow solution.
The size of the time step is determined using the acoustic Courant-Friedrichs-Lewy (CFL) condition <cit.>,
and we choose our CFL number to be 1.2.
The flow solver is implemented in Python using the adFVM <cit.> library,
which provides a high-level abstract application programming interface for writing efficient CFD applications.
The flow solver is parallelized using the Message Passing Interface (MPI) library. 



We use implicit Large Eddy Simulations (LES) in our numerical simulation.
In an LES, the large scale eddies of the flow are resolved by the grid, 
while the contribution from the small scale eddies to the filtered Navier-Stokes equations 
are modeled using a sub-grid scale Reynolds stress model <cit.>.
In this paper, the numerical error of the discretization scheme serves as the LES model.
It has been shown that when using a relatively dissipative discretization method,
the numerical viscosity from the grid can be of the same order of magnitude as the sub-grid scale viscosity
<cit.>,
and thus can be regarded as an implicit LES model.



On the inlet boundary, we specify stagnation pressure and temperature, corresponding to a fixed Reynolds number Re=1100 
and a Mach number which we choose to be the system parameter.
For the base case, we choose Mach number Ma=0.093.
Periodic boundary condition is used in the span-wise direction. 
The surface of the cylinder is maintained at a constant temperature of 300 K.
Static pressure of 1 atmosphere unit is prescribed on the outlet boundary. 



A snapshot of the flow field simulated with above settings is shown in figure <ref>.
As we can see, this flow exhibits weak turbulence in the wake.
In particular, the top view shows that this flow is 3D.

 
 



We choose our system parameter as the Mach number of the incoming flow.
The objective function is the time-averaged normalized drag over the cylinder. 
More specifically,

    J_avg = 2/ρ_r u_r^2 zd∫ (pn_x -μn·∇ u_x)  dS dt,.

Here the second integral is over the surface of the cylinder;
u_r = 31.4 m/s and ρ_r = 1.3 kg/m^3 are the reference velocity and density of the base case, where Ma=0.093.
For the base case, the normalized drag and the drag coefficient are the same, whereas for other Mach numbers they are different.



We run the flow simulation for 10^6 time steps, which corresponds to approximately 720 t_r.
Here the time unit t_r is the amount of time that the flow takes to traverse the length of the cylinder, that is, t_r = d/u_r. 
This time interval is sufficient to obtain a statistically converged estimate of the design objective.
The standard deviation of the time-averaged objective is computed using
the autoregressive time series analysis techniques described in <cit.> and <cit.>,
and we use one standard deviation as the confidence interval.
The normalized drag for the base case is 1.2 ± 0.03.
Our results reasonably matches the results from experiments <cit.>,
which is approximately 1.0 ± 0.15.
Figure <ref> shows different objectives for different incoming Mach number.






We first estimate the sensitivity by the linear least-squares regression method using 5 data points with different parameter values.
To use this linear regression method, we need to make the assumption that the relation between parameters and objectives is linear.
We select one standard deviation of the relevant estimator as our confidence interval.
Note that one shortcoming of the linear regression method is that the linear assumption may not be true when parameters are spaced far apart
such that a linear approximation no longer holds on the dataset;
on the other hand, when parameters are too close, the uncertainty in the objectives will lead to large error in the sensitivity.
In the base case, the sensitivity of drag with respect to the inflow Mach number, given by linear regression, is 25.0 ± 2.1. 
This sensitivity is visualized in figure <ref>.




Adjoint LEs are shown in figure <ref>,
where the confidence interval is also selected as one standard deviation given by autoregressive time series analysis.
The subsonic flow over a 3D cylinder has m_us = 9 unstable adjoint CLVs.





In NILSAS, the number of homogeneous adjoint solutions computed is set to M = 20.
The number of segments in NILSAS is K=100 and the number of time steps per segment is 500. 
Each segment roughly corresponds to 0.4t_r. 
Consequently, the time length of trajectory used in NILSAS is 40t_r, 
which is much lower than that required to obtain a reasonably accurate sensitivity using the linear regression method.
In this particular implementation, corresponding to the discussion in section <ref>,
on segment [t_i, t_i+1], we approximate integrations in equation (<ref>) using snapshots at t_i+1.
As a result, we have C_i=I, d_i^wv^* = 0:
this approximation eases the implementation responsible for storing adjoint solutions.








The sensitivity computed by NILSAS is 20.8 ± 3.5, which is visualized in figure <ref>.
Here the confidence interval is also selected as one standard deviation given by autoregressive time series analysis.
Comparing to the sensitivity estimated via linear regression methods, the relative difference is less than 20%. 
As we can see, the sensitivity computed by NILSAS correctly reflects the trend between parameters and objectives.



We remark that NILSAS may work for systems do not strictly satisfy the assumptions in theorem <ref>,
and our fluid problem is such an example.
First, our system has at least two neutral CLVs:
the first one corresponds to the common time translation of continuous dynamical systems,
and the second corresponds to span-wise translations due to the periodic boundary conditions.
Second, due to the similarity of this fluid problem with the one investigated in <cit.>,
whose tangent CLVs appear to have occasional tangencies,
it is reasonable to assume that adjoint CLVs in our current system may also have occasional tangencies.
Still, like NILSS did in <cit.>,
NILSAS computes a correct sensitivity: this encourages us to test NILSS and NILSAS on more general chaotic systems.



We compare computational cost for sensitivity analysis via the linear regression method and NILSAS.
The linear regression method runs the primal solver for a total of 5× 10 ^6 steps.
NILSAS runs the primal solver and 21 adjoint solvers, each for 5× 10 ^4 steps, which leads to
 1.1× 10^6 steps in total.
To build more favor towards NILSAS, note that,
first, adjoint solvers can be further accelerated due to the vectorization we discussed in section <ref>;
second, NILSAS has no additional cost for sensitivities to multiple parameters.
For chaotic problems with a higher number of positive LEs, 
the cost of NILSAS increases;
however, if the percentage of positive LEs is still low,
the non-intrusive formulation can still be a key technique for designing fast sensitivity algorithms.




§ CONCLUSIONS




To compute the gradient of long-time averaged objectives in chaotic systems,
we develop the Non-Intrusive Least Squares Adjoint Shadowing (NILSAS) algorithm, 
which approximates the adjoint shadowing direction by a `non-intrusive' formulation of a least squares problem.
NILSAS is demonstrated on the Lorenz 63 system and a weakly turbulent 3D flow over a cylinder,
where it gives accurate sensitivities for both cases.



Similar to NILSS <cit.>, NILSAS can be implemented with little modification to existing adjoint solvers,
and its minimization is carried out only in the unstable adjoint subspace.
Unlike NILSS, NILSAS has the benefit of adjoint approaches that its cost does not increase with the number of parameters;
thus making NILSAS ideal for applications where there are many parameters, or where f_s is unknown a priori.
NILSAS does not require tangent solvers, and is easy to implement.






§ SOLVING THE NILSAS PROBLEM


We discuss one way to solve the NILSAS problem in equation (<ref>).
The corresponding Lagrange function is:

    ∑_i=0^K-11/2 (a_i)^T C_i a_i + (d^wv^*_i)^T a_i  
    
      +    ∑_i=1^K-1λ_i^T ( a_i-1 - R_i a_i - b_i )
      + λ' ( ∑_i=0^K-1 (d^wf_i)^T a_i + ∑_i=0^K-1 d^v^*f_i )  ,

where λ_i is the Lagrange multiplier for the continuity condition at t_i.
By the Lagrange multiplier method,
the minimizer for the NILSAS problem is at the solution of the following linear equation systems:

    C     B^T
     B     0  a
    λ 
      = -d 
     b  ,

where the block matrices C∈^MK × MK, B∈^(MK-M+1)× MK,
vectors a, d∈^MK and λ, b ∈^MK-M+1.
More specifically, 

    C = C_0 
       C_1 
          ⋱
             C_K-1
      ,   
      B = I    -R_1
       I    -R_2 
          ⋱   ⋱
              I    -R_K-1
    
      (d^wf_0)^T       ⋯      (d^wf_K-1)^T, 
        a = a_0 
    ⋮
    a_K-1
      ,   λ = λ_1 
    ⋮
    λ_K-1
    λ'
      ,   
      d = d^wv^*_0 
    ⋮
    d^wv^*_K-1
      ,   
      b = b_1 
    ⋮
    b_K-1
     -  ∑_i=0^K-1 d^v^*f_i
      ,
 
where {C_i}_i=0^K-1, {R_i}_i=1^K-1⊂^M× M;
{a_i}_i=0^K-1, {d^wf_i}_i=0^K-1, {d^wv^*_i}_i=0^K-1, {λ_i}_i=1^K-1, {b_i}_i=1^K-1⊂^M;
λ', {d^v^*f_i}_i=0^K-1⊂.


We can solve the Schur complement of equation (<ref>) for λ:

    - B C^-1 B^T λ = B C^-1d + b  ,

where C^-1 can be computed via inverting each diagonal block in C.
Then we can compute a by:

    a = - C^-1(B^T λ +d)  .





§ NILSAS ON DISCRETE SYSTEMS





 §.§ Backgrounds and notations



We provide a brief introduction on discrete dynamical systems, in particular, hyperbolic diffeomorphisms.
More details are provided in <cit.>.
First, the governing equation for a discrete dynamical system is: 

    u_l+1 = f(u_l,s),    l≥ 0 .

The objective is:

    J:= lim_N→∞1/N∑_l=0^N-1 J(u_l,s).

Similar to flows, we assume u_l∈^m, and f(u,s) and J(u,s) are smooth.
We also assume f is a diffeomorphism in u, that is, for each fixed s, f has a smooth inverse.
Also, for simplicity of notations, we assume there is only one parameter s∈.


We first look at the tangent equations.
The homogeneous tangent diffeomorphism is: 

    w_l+1 = f_ul w_l  .

where the second subscript of f_ul indicate where the partial derivative is evaluated, that is,
f_ul:= ∂ f/ ∂ u (u_l,s).
A tangent CLV with exponent λ is a homogeneous tangent solution {ζ_l}^∞_i=0 such that there is constant C,
for any integer l_1, l_2, ζ_l_2≤ C e^λ(l_2-l_1)ζ_l_1. 
The uniform hyperbolicity for diffeomorphisms is defined as that all LEs are not 1.



On the adjoint side, the homogeneous adjoint diffeomorphism is defined as:

    _l = f_ul^T _l+1 ,

where ·^T is the matrix transpose.
The particular inhomogeneous adjoint diffeomorphism we will be using is:

    _l = f_ul^T _l+1 + J_ul .

On a trajectory {u_l}_l=0^∞ on the attractor,
the adjoint shadowing direction {_l}_l=0^∞ is a sequence with the following properties:

  
  * {_l}_l=0^∞ solves an inhomogeneous adjoint equation:
    
    _l = f_ul^T _l+1 + J_ul ,

  
  * _0 has zero component in the unstable adjoint subspace. 
  
  * _l is bounded by a constant for all l≥ 0.




It was proved in <cit.> that for a uniform hyperbolic diffeomorphism with a global compact attractor,
on a trajectory on the attractor, there exists a unique adjoint shadowing direction.
Further, we have the adjoint sensitivity formula:

    Js = lim_N→∞1/N∑_l=0^N-1( _l+1, f_sl + J_sl)  .







 §.§ Procedure list of NILSAS




We provide a procedure list for the NILSAS algorithm on discrete chaotic systems, more specifically, hyperbolic diffeomorphisms.
To start with, we need an inhomogeneous adjoint solver and a homogeneous adjoint solver, both can take arbitrary terminal conditions.
The inhomogeneous adjoint equation still has right-hand-side -J_u, same as many existing adjoint solvers for discrete systems.
We provide the following data: 
1) the number of homogeneous adjoint solutions, M≥, where  is the number of unstable CLVs, 
note that because the lack of neutral CLV, we can use one less homogeneous adjoint solution; 
2) the total number of segments, K; 
3) number of steps in one segment, L.



We can have three subscripts, the first, typically being u or s, indicates this term is a partial derivative;
the second, typically being i, 0, or K, indicates the segment number;
the third, typically being l, 0, or L, indicates the step number inside a segment.
Disappearance of the first subscript means that term is not a partial derivative.
Disappearance of the third subscript means either we are considering all steps in a segment, 
or that term is defined only once per segment interface.




  
  * Integrate the primal system for sufficiently many steps so that the initial condition, u_00, is on the attractor.

  
  * Compute the trajectory u_il for 0≤ i≤ K-1 and 0≤ l ≤ L.
    Here we assume the step at end of each segment overlaps with the start of next segment,
    that is, u_iL = u_i+1,0.

  
  * Generate terminal conditions for _i and ^*_i on the last segment i=K-1:
  
    
  * Randomly generate a m× M full rank matrix, Q'.
      Perform QR factorization: Q_K R_K = Q'.
    
  * Set p_K = 0.
  

  
  * Compute _i and v^*_i on all segments.
  For i=K-1 to i=0 do:
  
    
  * To get _il, whose columns are homogeneous adjoint solutions on segment i, solve:
      
    _il = f_uil^T _i,l+1 ,   _iL = Q_i+1 .

      To get ^*_i(t), solve the inhomogeneous adjoint equation:
      
    _il = f_uil^T _i,l+1 + J_uil ,  ^*_iL = p_i+1 .


    
  * Compute the following integrations.
      
    C_i = ∑_l=1^L_il^T _il dt  ,   
            d^wv^*_i = ∑_l=1^L_il^T ^*_il dt  ,   
            d^J_s_i = ∑_l=1^L J_sil dt  ,
       d^wf_s_i =∑_l=1^L_il^T f_si,l-1 dt  ,   
            d^v^*f_s_i = ∑_l=1^L^*T_il f_si,l-1 dt  ,

      where d^wv^*_i, d^wf_s_i∈^M; d^v^*f_s_i, d^J_s_i∈;
      C_i∈^M× M is the covariant matrix.
      Note that when multiplying adjoint solutions with f_s, their time steps are not the same:
      this asymmetry is the same as that in equation (<ref>).
      We are not sure yet if this technical detail can be neglected in practice.
    
    
  * Orthonormalize homogeneous adjoint solutions via QR factorization:
      
    Q_i R_i = _i0


    
  * Rescale the inhomogeneous adjoint solution using Q_i: 
      
    p_i = ^*_i0 - Q_i b_i  ,  where  b_i = Q_i^T ^*_i0 .

  


  * Compute the adjoint shadowing direction {_il} for  0≤ i≤ K-1 and 0≤ l ≤ L.
  
    
  * Solve the NILSAS problem on multiple segments:
      
    min_a_0,⋯,a_K-1∈^M   ∑_i=0^K-11/2 (a_i)^T C_i a_i +  (d^wv^*_i)^T a_i,   
       a_i-1 = R_i a_i + b_i  ,     i=1,⋯,K-1  .

      This is a least squares problem in {a_i}_i=0^K-1⊂^M.
      Note we do not have the other constraint as NILSAS in the continuous case.


    
  * On each time segment i, _il is given by
      
    _il = ^*_il + _il a_i .

  

  
  * Compute the derivative by:
    
    dJ/ds≈1/KL∑_i=0^K-1( d^v^*f_s_i + a_i^T d^wf_s_i + d^J_s_i )




model1-num-names


