On the Benefits of Memory for Modeling Time-Dependent PDEs
Ricardo Buitrago Ruiz Tanya Marwah Albert Gu Andrej Risteski



Topological characterization of modified Kane-Mele-Rashba models via local spin Chern marker
    Tarik P. Cysne
    September 9, 2024 – Version 1.0
============================================================================================



§ ABSTRACT 
    Data-driven techniques 
have emerged as a promising alternative to traditional numerical methods 
for solving partial differential equations (PDEs).
These techniques frequently offer a better trade-off between computational cost and accuracy for many PDE families of interest. For time-dependent PDEs, existing methodologies typically treat PDEs as Markovian systems, i.e., 
the evolution of the system only depends on the “current state”, and not the past states. However, distortion of the input signals — e.g., due to discretization or low-pass filtering — can render the evolution of the 
distorted
signals non-Markovian. 
In this work, motivated by the Mori-Zwanzig theory of model reduction,
we investigate the impact of architectures with memory 
for modeling PDEs: that is, when past states are explicitly used to predict the future. 
We introduce Memory Neural Operator (MemNO), 
a network based on the recent SSM architectures and Fourier Neural Operator (FNO). 
We empirically demonstrate on a variety of PDE families of interest that when the input is given on a low-resolution grid, MemNO significantly outperforms the baselines without memory, 
achieving more than 6 × less error on unseen PDEs. Via a combination of theory and experiments, we show that the effect of memory is particularly significant when the solution of the PDE has high frequency Fourier components (e.g., low-viscosity fluid dynamics), and it also increases robustness to observation noise.













































§ INTRODUCTION

Time-dependent partial differential equations (PDEs) 
are central to modeling
various scientific and physical phenomena,
necessitating the design of accurate and computationally efficient solvers.
Recently, data-driven neural network based approaches <cit.>  have emerged as an attractive alternative to classic numerical solvers, such as finite element and finite difference methods <cit.>.  


Classical approaches are computationally expensive in high dimension, and struggle with PDEs which are very sensitive to initial conditions. Learned approaches can often negotiate these difficulties better, at least for the family of PDEs they are trained on. 

One example of a data-driven approach is learning a neural solution operator, which for a time-dependent PDE learns a time 
evolution map that predicts the solution of the PDE for future time steps.



The operators are frequently autoregressively parametrized, such that the network predicts future states based on previous ones <cit.>, and the number of past states the model is conditioned on serves as “memory” and is treated as
a tunable hyperparameter.

Recent works <cit.> 
suggest that optimal performance across various PDE families 
can be achieved by conditioning the models only on the immediate past
state—i.e., treating the system as Markovian—though this is in settings in which the training data is very high-resolution. 

In many practical settings, we expect to only observe a part of the system. This could be due to limited resolution
of the measurement devices collecting the data, inherent observational errors in the system, or prohibitive computational difficulty generating high-quality synthetic data. 
This can lead to significant information loss, particularly in systems like turbulent flows <cit.> or shock formation in fluid dynamics <cit.>, 
where PDEs change abruptly in space and time.
In such situations, classical results from dynamical systems (Mori-Zwanzig theory), suggest that the system becomes strongly non-Markovian. 

More precisely, Mori-Zwanzig theory <cit.> is an ansatz to understand the evolution of a subspace of system (e.g., the top k Fourier components). Under certain conditions, this evolution can be divided into a Markovian term (the evolution of the chosen subspace under the PDE), a memory term (which is a weighted sum of the values of all previous iterates in the chosen subspace), and an “unobservable” term, which depends on the values of the initial conditions orthogonal to the chosen subspace.  

we first show that if the input data contains
observes only a part of the PDE—for example, 
PDE defined on a low-resolution grid, if the data contains only 
the top k frequency components of a state—then the effects 
of the unobservable part of the system 
can be captured by 
as a memory of the observed variables in the previous state.
We then show that the dynamics of a PDE
can be divided into 
observed and unobserved terms, thus,
the resulting dynamics of the observed variables 
contain a memory term which is a weighted sum of all the previous observed states
plus an unobserved noise term. this is a strange description of MZ— and “thus” is a strange word. You want to say that the projection of the dynamics to the observable subspace has a Markovian term + a “correction” factor that also lies in the observed subspace but isn't MarkovianWe also prove that
under certain conditions, the memory term can be arbitrary large.
Therefore, 
the resulting dynamics of the system not Markovian, necessitating 
the need to maintain some form of a memory of the past states. 


In this paper, we study the effects of explicit
memory when deploying neural operators for time-dependent PDEs. By memory
we loosely mean a representation of the previous states
of a PDE. 
We focus on PDE families that take the form  ∂_t u(x,t) =  u(x, t)
where u(x, t): Ω× [0, T] → is a time-dependent function 
defined over the domain Ω, and  is a (possibly non-linear) operator.
This is a generic form of a time-dependent PDE system and contains many PDE families of interest in different application domains (e.g., heat diffusion, Navier-Stokes, Kuramoto-Sivashinsky, Black-Scholes, Schrödinger equation to name a few). 





We introduce Memory Neural Operator (MemNO), 
an architecture which combines Fourier neural operator (FNO) <cit.>
and 
the S4 architecture <cit.>. 
The MemNO architecture can be seen as an adaptation of the FNO <cit.>
architecture where the FNO layers model the spatial dynamics of the PDE 
while the S4 layers <cit.>
maintain a compressed memory of the past states. 
We choose S4 models over
recurrent architectures like LSTM <cit.>
due to superior performance in modeling long range dependencies <cit.>, ease of training, 
and favorable memory and computational scaling with both state dimension 
and sequence length.
Through our experiments we 
show that for PDEs observed on low resolution grids and/or with observation noise, 
MemNO outperforms their Markovian (memoryless) baselines—achieving 6 × less
loss on unseen PDEs.  Our contributions are as follows:

    
  * We introduce MemNO, a Memory Neural Operator architecture which uses a combination of FNO layers and S4 layers to model the spatial and temporal dynamics of a PDE. The S4 layers explicitly model memory, by introducing a weighted average over the past (compressed) states. 
    
  * Even in relatively simple linear PDEs, we theoretically show the memory term can result in a solution that is (arbitrarily) closer to the correct solution, compared to the Markovian approximation — in particular when the operator describing the PDE “mixes” the observed and unobserved subspace.  
    
  * Across several families of one-dimensional and two-dimensional PDEs, we show 
    that when the input is supplied on a 
    low-resolution grid, or has additional observation noise, memory based architectures outperform the best performing FNO based baselines by a significant margin. 
    
  * Finally, we empirically show that this effect is more pronounced for PDEs which result in solutions with high order
    frequency modes, and introduce a metric which indicates 
    when memory based models will have most impact.
    
    
    













































 









§ RELATED WORK

Data-driven neural solution operators <cit.>
have emerged as the dominant approach for 
approximating 
PDEs, 
given their ability to model multiple families of PDEs at once, 
and relatively fast computation at inference time.

Recently, many architectures have been proposed to improve the performance of neural 
operators across multiple families of PDEs,


 <cit.> designed
the Fourier Neural Operator (FNO), a resolution invariant 
architecture
that uses a convolution based integral kernels 
evaluated in the Fourier space.  <cit.> later
introduced Factorized FNO (FFNO) architecture, that builds upon and improves
the FNO architecture by adding separable spectral layers and residual connections.
Additionally, they perform extensive ablations that point out
that strategies like the Markov 
assumption—i.e., predicting a state only from its immediate prior state—is optimal 
and outperforms models that use the history of past timesteps as input.  <cit.> performed a similar study for long rollouts of the PDE solution and concluded the the optimal performance is indeed achieved under the Markovian assumption.





However, we show that 
when there is a loss of information in the observation a PDE,


a model that the uses the history 
of past states

outperforms its
Markovian counterpart, often achieving 6× less error on unseen PDEs from the same family.

Our work is motivated by the Mori-Zwanzig formalism <cit.>
which shows that
a partial observation of the current state of the system 
can be compensated using memory of past states. 








Our work is also inspired by a previous study on the effects of memory 
in modeling PDE dynamics by  <cit.>. Here the authors 
draw parallels to the Mori-Zwanzig equations and LSTM <cit.>  
to model the dynamics of the top Fourier components of a time-dependent 1D Kuramoto-Sivashinsky
and 2D shear flow equations, one single PDE at a time.
However, in our work, we study the benefits of memory in neural operator settings, i.e, 
we have a single model that learns the dynamics of an entire family of PDE at once.
Furthermore, we use
the S4 state space model architecture <cit.> to model
the temporal dependencies, which in our experiments has better performance and is more stable than LSTMs.









§ PRELIMINARIES


In this section, we introduce several definitions, as well as background on the Mori-Zwanzig formalism as applied to our setting.






 §.§ Partial Differential Equations (PDEs)


[Space of square integrable functions]
    For integers d, V and an open set Ω⊂ℝ^d, we define L^2(Ω; ℝ^V) as the space of square integrable functions u: Ω→ℝ^V such that ‖ u ‖_L^2≤∞, where ‖ u ‖_L^2 = (∫_Ω‖ u(x) ‖^2_2 dx )^1/2.






[Restriction]
    Given a function u: Ω→ℝ^V and a subset A⊂Ω, we denote u _A as the restriction of u to the domain A, i.e. u_A: A →ℝ^V. 
    


The general form the PDEs we consider in this paper will be as follows: 

[Time-Dependent PDE]
    
    For an open set Ω⊂ℝ^d and an interval [0,T]⊂ℝ, a Time-Dependent PDE is the following expression:
    
    ∂ u/∂ t (t,x)    = ℒ[u](t,x),        ∀ t ∈[0,T], x ∈Ω, 
    
            u(0,x)     = u_0(x),           ∀ x ∈Ω, 
    ℬ[u_∂Ω](t)    = 0,               ∀ t ∈ [0,T]

    where ℒ: L^2(Ω; ℝ^V)→ L^2(Ω; ℝ^V) is a differential operator in x which is independent of time, u_0(x)∈ L^2(Ω; ℝ^V) and ℬ is an operator defined on the boundary of ∂Ω, commonly referred as the boundary condition.
    
    


Unless otherwise stated, both in the experiments  and in the theory we will largely work with periodic boundary conditions: 

[Periodic Boundary Conditions]
    For Ω = [0,L]^d, we define the periodic boundary conditions as the condition: 
    
    u(x_1, ⋯, x_k-1, 0, x_k+1, ⋯ x_d) = u(x_1, ⋯, x_k-1, L, x_k+1, ⋯ x_d)

    for all (x_1, ⋯, x_k-1, x_k+1, ⋯, x_L) ∈ [0,L]^d-1 and all k=1, ⋯, d.


Finally, we will frequently talk about a grid of a given resolution:
 [Equispaced grid with resolution f]
Let  Ω=[0,L]^d. An equispaced grid with resolution f in Ω is the following set 𝒮⊂ℝ^d:

    𝒮 = {.(i_1 L/f, ⋯, i_k L/f) | 0 ≤ i_k ≤ f-1  for 1≤ k ≤ d }.

We will also denote by |𝒮| the number of points in 𝒮.












 §.§ Mori-Zwanzig



The Mori-Zwanzig formalism <cit.> deals with cases where an equation is known for a full system, yet only a part of it is observed. It leverages the knowledge of past observed states of a system to compensates for the loss of information that arises from the partial observation. In our paper, partial observation can refer to observing the solution at a discretized grid in space or only observing the Fourier modes up to a critical frequency. In particular, the Mori-Zwanzig formalism in the context of time-dependent PDEs is well-known in the Physics literature as the Nakajima–Zwanzig equation (<cit.>)





Now, we will apply the Nakajima–Zwanzig equation to our setting. Assume we have a PDE as in Definition <ref>. Let 𝒫: L ^2(Ω; ℝ^V) → L^2(Ω; ℝ^V) be a linear projection operator. We define 𝒬=I - 𝒫, where I is the identity operator. In our setting, for the PDE solution at timestep t u_t ∈  L ^2(Ω; ℝ^V), 𝒫[u_t] is the part of the solution that we observe and 𝒬[u_t] is the unobserved part. Thus, the initial information we receive for the system is 𝒫[u_0]. Applying 𝒫 and 𝒬 to Equation <ref> and using u = 𝒫[u] + 𝒬[u], we get:

    ∂/∂ t𝒫[u] (t,x)    = 𝒫ℒ[u] (t,x) = 𝒫ℒ𝒫 [u] (t,x) + 𝒫ℒ𝒬 [u] (t,x)  
    ∂/∂ t𝒬[u] (t,x)    = 𝒬ℒ[u] (t,x) = 𝒬ℒ𝒫 [u] (t,x) + 𝒬ℒ𝒬 [u]  (t,x)

Solving for <ref> yields 𝒬[u](t,x) = ∫_0^t exp𝒬ℒ(t-s)𝒬ℒ𝒫[u](s,x) ds + e^𝒬ℒ t𝒬[u_0](t,x).

Plugging into <ref>, we obtain a Generalized Langevin Equation <cit.> for 𝒫[u]:


    ∂/∂ t𝒫[u] (t,x)     = 𝒫ℒ𝒫 [u] (t,x) + 𝒫ℒ∫_0^t exp𝒬ℒ(t-s)𝒬ℒ𝒫[u](s,x) ds + 𝒫ℒ e^𝒬ℒ t𝒬[u_0](t,x)


We will refer to the first summand on the right hand side of <ref> as the Markovian term because it only depends on 𝒫[u](t,x), the second summand as the memory term because it depends on 𝒫[u](s,x) for 0≤ s ≤ t, and the third summand as the unobserved residual as it depends on 𝒬[u_0] which is never observed.

We note that (<ref>) is exact, not an approximation, so it is equivalent to solving the full system. Typically, the term that is most difficult to compute is the exponential of the memory term, and thus several methods to approximate it have been proposed. 
In the physics literature, the memory term has been approximated through a perturbation expansion of the exponential <cit.>, or by approximating the operator exp𝒬ℒ(t-s): L ^2(Ω; ℝ^V) → L ^2(Ω; ℝ^V) through operators defined in  𝒫[L ^2(Ω; ℝ^V)] <cit.>. In the machine learning literature, <cit.> develop the equations for the case when the operator 𝒫 kept only the top-k modes, and designed a hybrid approach where the memory term was approximated with an LSTM <cit.>, and then used as an additional input of a numerical solver. In this work, we treat the whole memory term as an operator ℳ: 𝒞([0,T],𝒫[L ^2(Ω; ℝ^V)]) →𝒫[L ^2(Ω; ℝ^V)][
Here 𝒞( A, B ) denotes the space of continuous functions u: A → B.
] to be learnt by a parametrized sequential layer of a Neural Operator.










§ OUR APPROACH












































 §.§ Training procedure




First, we describe the high-level training scaffolding for our method, namely the way the data is generated, and the loss we use.  

Training data: Let u∈𝒞([0,T];L^2(Ω; ℝ^V)) be the solution of the PDE given by Definition <ref>. Let 𝒮 be an equispaced grid in Ω with resolution f, and let 𝒯 be another equispaced grid in [0,T] with N_t+1 points. Given u_0(x)_𝒮, our goal is to predict u(t,x)_𝒮 for t∈𝒯 using a Neural Operator. 

Training loss: As it is standard, we proceed through empirical risk minimization on a dataset of trajectories. More specifically, given a loss function ℓ : (ℝ^|𝒮|,ℝ^|𝒮|)→ℝ, a dataset of training trajectories (u(t,x)^(i))_i=0^N, and parametrized maps 𝒢^Θ_t:ℝ^|𝒮|→ℝ^|𝒮| for t∈𝒯, we define:

    Θ^* = min_Θ1/N∑_i=0^N-11/N_t∑_t=1^N_tℓ( u(t,x) _𝒮, 𝒢^Θ_t[ u_0(x) _𝒮] )

We then aim to find an adequate architecture choice such that 𝒢^Θ^* has low test error on unseen trajectories of the same PDE.



 §.§ The Architecture: Memory Neural Operator


In this section we describe Memory Neural Operator (MemNO), a Deep Learning framework to incorporate memory into Neural Operators.

Let NO^Θ_t be a Neural Operator with L layers,
and denote NO^Θ_t[u_0] the prediction of the solution of the PDE at time t.
We will assume that this Neural Operator follows the Markovian assumption, i.e. we can write:


    NO^Θ_t_i+1[u_0] = r_out∘ℓ_L ∘ℓ_L-1∘ ... ∘ℓ_0 ∘ r_in[NO^Θ_t_i[u_0]]

Where r_in: ℝ^|𝒮|→ℝ^|𝒮| × h_0 and r_out: ℝ^|𝒮|× h_L+1→ℝ^|𝒮| are projector operators; ℓ_j: ℝ^|𝒮| × h_j→ℝ^|𝒮|× h_j+1 are parametrized layers; and h_j is the dimension of the j-th hidden layer. 

Our goal is to define a network 𝒢^Θ_t 
that builds upon NO^Θ_t and can incorporate memory. For this
we
take inspiration from the Mori-Zwanzig theory exposed in Section <ref>.
Comparing (<ref>) with (<ref>), 
we identify ℓ_L ∘ℓ_L-1∘ ... ∘ℓ_0 with the Markov term
which models the spatial dynamics. 
To introduce the memory term, we interleave an additional residual sequential layer ℳ that acts on hidden representations of the solution at previous timesteps. Concretely, the MemNO architecture can be written as:


    𝒢_t_i+1^Θ[u_0] = ℛ_out∘ℒ_L ∘ ... ∘ℒ_k+1∘ℳ∘ℒ_k∘ ... ∘ℒ_0 ∘ℛ_in[𝒢_t_i^Θ[u_0] , 𝒢_t_i-1^Θ[u_0] ,..., u_0 ]

Where -1≤ k ≤ L is a chosen hyperparameter.[
k=L referst to inserting M after all the S layers, and k=-1 refers to inserting M as the first layer.
As we show in Appendix <ref>, our experiments are not very sensitive to the choice of k.
]
Now, the spatial ℒ_j layers are understood to be applied timestep-wise. 
That is, if v^(j)(t') is the hidden representation at the j layer for a timestep t'≤ t_i, then ℒ_j+1[v^(j)(t_i), ..., v^(j)(t_0)] = [ℓ_j[v_i^(j)(t_i)], ..., ℓ_j[v_0^(j)(t_0)]],
and analogously for ℛ_in and ℛ_out. 

Thus, the ℒ_j layers still follow the Markovian assumption. The memory is introduced through ℳ, which consists of a sequential layer m that is applied to the time history of the hidden representation of the k-th layer, that is ℳ:ℝ^i× |𝒮| × h_k ⟶ℝ^|𝒮| × h_k with (ℳ[v^(k)(t_i), ..., v^(k)(t_0)])_sh = m[v_sh^(k)(t_i), ..., v_sh^(k)(t_0)]. Note that m is the same for each of the |𝒮| × h_k  elements of the hidden layer.


The main motivation of our MemNO framework is that it can be utilized with any
existing neural operator layer ℓ, and with any 
(causal) sequential model ℳ.

Thus it provides a modular architecture design 

which we hope  can serve as a useful tool for practitioners.



























§ THEORETICAL MOTIVATION FOR MEMORY: A SIMPLE EXAMPLE



In this section, we provide a simple, but natural example of a (linear) PDE, along with (in the nomenclature of Section <ref>) a natural projection operator given by a Fourier truncation measurement operator, such that the memory term in the generalized Langevin equation (GLE) can have an arbitrarily large impact on the quality of the calculated solution. 
We will work with periodic functions over [0,2π] which have a convenient basis: 
[Basis for 2π-periodic functions] A function f: ℝ→ℝ is 2π-periodic if f(x + 2π) = f(x). We can identify 2π-periodic functions with functions over the torus T:={e^iθ: θ∈ℝ}⊆ℂ by the map f̃(e^ix) = f(x). Note that {e^i x n}_n ∈ℤ is a basis for the set of 2π-periodic functions.


We will define the following measurement operator: 
[Fourier truncation measurement] 
  The operator 𝒫_k: L^2(T;ℝ) → L^2(T;ℝ) acts on f ∈ L^2(T;ℝ), f(x) = ∑_n=-∞^∞ a_n e^i n x as  𝒫_k(f) = ∑_n=-k^k a_n e^i n x.

We will also define for notational convenience the functions {𝐞_n}_n ∈ℤ, where  𝐞_n(x) := e^-inx + e^inx. 
Now we consider the following operator to define a linear time-dependent PDE: 
 Let ℒ: L^2(T;ℝ) → L^2(T;ℝ) be defined as 
ℒ u(x) = -Δ u(x) + B · (e^-ix + e^ix) u(x)
for B > 0. Then, we have: 


    ∀ 1 ≤ n ∈ℕ, ℒ (𝐞_n)    = n^2 𝐞_n + B (𝐞_n-1 + 𝐞_n+1)   &  ℒ (𝐞_0) = 2 B 𝐞_1



The crucial property of this operator is that it acts by “mixing” the n-th Fourier basis  with the (n-1)-th and (n+1)-th: thus information is propagated to both the higher and lower-order part of the spectrum. 
Given the above proposition, we can easily write down the evolution of a PDE with operator ℒ in the basis {𝐞_n}_n ∈ℤ: 
  Let ℒ be defined as in Proposition <ref>. Consider the PDE 

    ∂/∂ t u(t,x)    = ℒ u(t,x) 
     
        u(0,x)    = ∑_n ∈ℕ_0 a_n(0) 𝐞_n

Let u(t,x) = ∑_n ∈ℕ_0 a_n^(t)𝐞_n.
Then, the coefficients a_n^(t) satisfy: 

    ∀ 1 ≤ n ∈ℕ, ∂/∂ t a_n^(t)   =n^2 a_n^(t) + B ( a_n-1^(t) + a_n+1^(t))
    ∂/∂ t a_0^(t)   = 2B a_1^(t)



With this setup in mind, we will show that as B grows, the memory term in (<ref>) can have an arbitrarily large effect on the calculated solution:  

    
   Consider the Fourier truncation operator 𝒫_1 and let 𝒬 = I - 𝒫_1. Let u(0, x) have the form in Proposition <ref> for B > 0 sufficiently large, and let a^(0)_n > 0, ∀ n > 0. Consider the memoryless and memory-augmented PDEs: 

    ∂ u_1/∂ t   = 𝒫_1 ℒ u_1 
    ∂ u_2/∂ t   = 𝒫_1 ℒ u_2 + 𝒫_1ℒ∫_0^t exp𝒬ℒ(t-s)𝒬ℒ u_2(s) ds

with u_1(0,x) = u_2(0,x) = 𝒫_1 u(0,x). Then, u_1 and u_2 satisfy: 
   
    ∀ t > 0, u_1(t) - u_2(t)_L_2   ≳  B t u_1(t)_L_2
    ∀ t > 0, u_1(t) - u_2(t)_L_2   ≳ B texp(√(2) Bt)




    Note that the two conclusions of the theorem mean that both the absolute difference, and the relative difference between the PDE including the memory term (<ref>) and not including the memory term (<ref>) can be arbitrarily large as B, t →∞.  




    The choice of ℒ is made for ease of calculation of the Markov and memory term. Conceptually, we expect the solution to (<ref>) will differ a lot from the solution to (<ref>) if the action of the operator ℒ tends to “mix” components in the span of  𝒫 and the span of 𝒬.
    




If we solve the equation ∂/∂ t u(t,x) = ℒ u(t,x) exactly, we can calculate that u(t)_L_2 will be on the order of exp(2Bt). This can be seen by writing the evolution of the coefficients of u(t) in the basis {𝐞_n}, which looks like: 
∂/∂ t[ a_0; a_1;   … ] = 𝒪[ a_0; a_1;   … ] 
where 𝒪 is roughly a tridiagonal Toeplitz operator 𝒪 =  [               ⋮       ⋮       ⋮       ⋮        ;       …       B     n^2       B       0       …;       …       0       B (n+1)^2       B       …;               ⋮       ⋮       ⋮       ⋮         ].
The largest eigenvalue of this operator can be shown to be on the order of at least 2B (equation (4) in <cit.>). The Markov term results in a solution of order exp(√(2)Bt) ( (<ref>),(<ref>)), which is multiplicatively smaller by a factor of exp((2 - √(2))Bt). The result in this Theorem shows the memory-based PDE (<ref>) results in a multiplicative “first order” correction which can be seen by Taylor expanding exp(√(2)Bt) ≈ 1 + √(2)Bt + 1/2(√(2)B)^2t^2 + ….     




§ MEMORY HELPS WITH LOW-RESOLUTION DATA AND INPUT NOISE: A CASE STUDY


In this section we present a case study for several common PDEs of practical interest, showing that MemNO brings accuracy benefits when the data is supplied in low resolution.


Through our experiments we show 
the difference in the performance between

a baseline “memoryless” architecture, which we choose to be Factorized Fourier Neural Operator (FFNO) <cit.> 
and a memory-augmented architecture using S4 <cit.>, which we denote as the S4-Factorized Fourier Neural Operator (s4FFNO). 
The architectural details for both the 
architectures are elaborated upon in 
Appendix <ref>.



 §.§ Setup: Training and evaluation procedure


To construct our datasets, we first produce discretized trajectories of a PDE for N_t timesteps, i.e. (u(t))_t=0^N_t in a high resolution discretized spatial grid 𝒮^HR⊂ℝ^d, i.e. u(t)∈ℝ^|𝒮^HR|. 
We then produce datasets that consist of lower resolution versions of the above trajectories, i.e. on a grid 𝒮^LR of lower resolution f. The low resolution grid would correspond to the observation grid obtained through a physical measurement device, whereas the high resolution grid is used in the numerical solver to obtain trajectories that follow real-world PDE dynamics.
For 1-dimensional datasets, the discretized trajectory on 𝒮^LR is obtained by cubic interpolation of the trajectory in the high resolution grid. In 2D, the discretized trajectory is obtained by downsampling. 
We will show results in different resolutions, in which case both train and test trajectories are at such resolution, and the loss function is also computed at the chosen resolution. 
Our training loss and evaluation metric is normalized Root Mean Squared Error (nRMSE): 

    nRMSE(u(t,x) _𝒮^LR, û(t) ) = ‖ u(t,x) _𝒮^LR - û(t) ‖_2/‖ u(t,x) _𝒮^LR‖_2 ,

where ∥·∥_2 is the euclidean norm in ℝ^|𝒮^HR|. More details on training are given in appendix <ref>.














 §.§ Kuramoto–Sivashinsky equation (1D): a study in low-resolution




The Kuramoto-Sivashinsky equation (KS) is a nonlinear PDE that is used as a modeling tool in fluid dynamics, chemical reaction dynamics, and ion interactions. Due to its chaotic behavior it can model instabilities in various physical systems. For a viscosity ν, it is written as u_t + uu_x + u_xx + ν u_xxxx = 0.
We generated datasets for KS at different viscosities and resolutions. The results are shown in Table <ref>. We can see s4FFNO outperforms FFNO across these viscosities and resolutions, having an nRMSE that can be more than six times smaller. 

We also note that, since the memory layer is applied element-wise in the time dimension, it has very few parameters compared to the spatial layers, as seen in the difference of parameters between s4FFNO and FFNO in column 3.






The key factor for the improved performance of MemNO over memoryless Neural Operators in not the absolute resolution, but rather the resolution relative to the frequency spectrum of the solution. The lower the viscosity, the higher the frequencies that appear in the spectrum. This can be clearly seen in Figure <ref>: in the top row, as viscosities increase, the resolution at which there is a significant difference between s4FFNO and FFNO decreases. In the second row of the figure, we show a visualization of the frequency spectrum of a solution at those viscosities. 



We note that even if the initial condition does not contain high frequencies, in the KS equation high frequencies will appear as the system evolves—indeed, this dataset was generated with initial conditions whose maximum Fourier mode was 8.  This is in qualitative agreement with the theoretical motivation in Section <ref>—albeit the KS equation is substantially more complicated since it is nonlinear, so it is hard to fully theoretically analyze. 
We provide a similar study on 1D Burgers equation 
in the Appendix <ref>.









































    











 §.§ Navier Stokes equation (2D): study in observation noise


The Navier Stokes equation describes the motion of a viscous fluid. 
Like in <cit.>, we consider the incompressible form in the 2D unit torus, which is given by: 

    ∂ w(x,t)/∂ t + u(x,t) ·∇ w(x,t)    = νΔ w(x,t) + f(x),     x ∈ (0, 1)^2, t ∈ (0,T]
    
    ∇· u(x,t)    = 0,     x ∈ (0, 1)^2, t ∈ [0,T]
    
    
        w(x, 0)    = w_0(x),     x ∈ (0, 1)^2

Where w = ∇× u is the vorticity, w_0 ∈ L^2((0, 1)^2; ℝ) is the initial vorticity, ν∈ℝ_+ is the viscosity coefficient, and f ∈ L^2((0, 1)^2; ℝ) is the forcing function. In general, the lower the viscosity, the more rapid the changes in the solution and the harder it is to solve it numerically and with a Neural Operator.
We investigate the effect of memory when adding IID Gaussian noise to the inputs of our neural networks. Noise would represent the observation noise arising from the intrisic error of the measurement device. The noise ϵ∈ℝ^|𝒯|× |𝒮^LR| is sampled IID from a Gaussian distribution ϵ_ts∼𝒩(0, σ), and then added to training and test inputs. During training, for each trajectory a different noise with the same σ is sampled at each iteration of the optimization algorithm. The targets in training and testing represent our ground truth and are not added noise.

In Figure <ref>, we show the results for ν=10^-3 when adding noise levels from σ=0.0 (no noise) to σ=2.048. s4FFNO-2D outperforms FFNO-2D across most noise levels, and the difference between the two is especially significant for noise levels beyond 0.128, where FFNO-2D is around 50% higher than s4FFNO-2D (note the logarithmic scale). For this viscosity, adding small levels of noise actually helps training, which was also observed in other settings in <cit.>. Figure <ref> shows the same experiment performed with ν=10^-5. Again, s4FFNO-2D outperforms FFNO-2D across most noise levels. FFNO-2D losses are similarly around 50% higher for noise levels above 0.032. In this viscosity, adding these levels of noise does not help performance.









 §.§ Relationship with fraction of unobserved information













In this section, we provide a simple experiment to quantify the 
effect of the fraction of unobserved information on the performance of memory based models.
Given a grid of resolution f, we define the Fourier truncation measurement 𝒫_f/2 as in Section <ref>, 
which simply keeps the top f/2+1 modes and discards the other high frequency modes.
Assume u(t)∈ L^2(Ω; ℝ^V) is the solution of a 1-dimensional PDE at time t,  and a_n^(t) for n∈ℤ is its Fourier Transform. We define the quantity:

    ω_f = 1/N_t∑_i=1^N_t∑_|n|≥f/2 |a_n^(t_i)|^2/∑_n∈ℤ |a_n^(t_i)|^2

ω_f is approximate indicator of the
amount of information that is lost when the solution of the PDE is observed at resolution f across time.

We show that there is a positive correlation between ω_f and the difference in nRMSE between FFNO and s4FFNO for the KS experiment in Figure <ref>, and also the for Burgers' experiments of  Appendix<ref> in Figure <ref>.
This demonstrates the benefits
of memory as a way to compensate missing information in the observation.






















[18]ht0.4
    
    < g r a p h i c s >

    Values of ω_f and the difference in nRMSE between FFNO and s4FFNO for different resolutions in the KS experiment of Section <ref> with ν=0.1.
    ω_f is averaged over all trajectories in the dataset. The value is computed approximating the continuous Fourier modes with Discrete Fourier modes of the solution in the highest resolution available (512 for KS).
    



§ CONCLUSION AND FUTURE WORK

We study the benefits of maintaining memeory while modeling 
time dependent PDE systems. Taking inspiration from the Mori-Zwanzig formulation, we show that when 
we only observe part of the PDE initial condition (for example, PDEs observed on low-resolution or with
input noise), the system is no longer Markovian, and the dynamics depend on a 
memory term.


To this end, we introduce 
MemNO, an architecture that combines Fourier Neural Operator (FNO) and 
the S4 architecture. 
Through our experiments on different 1D and 2D PDEs, we show that the MemNO
architecture outperforms the memoryless baselines. 

We present several avenues for future work. First, our experiments on 
observation noise are limited to the setting where the input noise is IID. 
Further, extending the experiments and observing the effects of memory in
more real-world settings (for example, with non-IID noise or in the presence of aliasing) is a fertile ground for future work, and also necessary to ensure that the application of this method does not have unintended negative consequences when broadly applied in society. 

Lastly, while we limit our 
study of the effects of memory to FNO based architectures, performing 
similar studies for different architectures like Transformer based neural operators <cit.>
and diffusion based operators <cit.> 
is an interesting direction for future work.









§ ACKNOWLEDGEMENTS

RBR is supported by the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/EU22/11930090.
TM is supported in part by CMU Software Engineering Institute via Department of Defense under contract FA8702-15-D-0002.
AR is supported in part by 
NSF awards IIS-2211907, CCF-2238523, and Amazon Research.
The authors also thank Cartesia AI for their generous provison of computational resources.

plainnat




§ EXTENDED RELATED WORK

Neural Operators. The Fourier Neural Operator (FNO) is a Neural Operator that performs a transformation in the frequency space of the input <cit.>. Other models have proposed different inductive biases for Neural Operators, including physics based losses and constraints <cit.>,
using Deep Equilibrium Model (DEQ) <cit.>
 to design specialized
architectures for steady-state (time-independent) 
PDEs <cit.>, and using local message passing 
Graph Neural Networks (GNNs) <cit.>
based encoders to model irregular geometries <cit.>.
Other methodologies to solve for 
PDEs include methods like  <cit.>
that use the U-Net <cit.>
type architectures 
and works like <cit.>
that introduce different Transformer <cit.> based 
neural solution operators for modeling both time-dependent and time-independent PDEs.
While most of these methodology 
are designed for time-dependent PDEs, 
there is no clear consensus of how to model 
the past-states to predict future states, and most
of these methods predict the PDE states over time
in an auto-regressive way by conditioning 
the model on varying lengths of the past states <cit.>.

Foundation models. Lately, there have been community efforts 
towards
creating large scale
foundational models for modeling multiple PDE families <cit.>, and weather prediction <cit.>.
We hope that our study is useful in informing the architectural 
design of future models.


§ NETWORK ARCHITECTURES


Factorized Fourier Neural Operator (FFNO) (<cit.>):  This model is a refinement over the original Fourier Neural Operator (<cit.>). Given a hidden dimension h and a spatial grid 𝒮, its layers ℓ: ℝ^|𝒮|× h→ℝ^|𝒮|× h are defined as:

    ℓ(v) = v + Linear_h,h'∘σ∘Linear_h',h∘𝒦 [v]


where σ is the GeLU activation function <cit.> and h' is an expanded hidden dimension. 𝒦 is a kernel integral operator that performs a linear transformation in the frequency space. Denoting by 

FFT_α, IFFT_α are to the Discrete Fast Fourier Transform and the Discrete Inverse Fast Fourier Transform along dimension α <cit.>, it can be written as:

    𝒦[v] = ∑_α∈{1,...,d}IFFT[R_α·FFT_α[v]]

for learnable matrices of weights R_α∈ℂ^h^2 × k_max. k_max is the maximum number of Fourier modes which are used in 𝒦. We use all Fourier modes by setting k_max = f/2. 

In our experiments, The FFNO model consists of 4 FFNO layers. For experiments in 1D, the hidden dimensions are all 128 (h_j=128 for j=0,1,2,3) and the expanded hidden dimension of FFNO's MLP h' is 4 · 128. For experiments in 2D, the hidden dimensions are all 64 and the expanded hidden dimension is 4 · 64.  



S4 - Factorized Fourier Neural Operator (s4FFNO): This model uses our MemNO framework. To discern the effect of memory, all layers except the memory layer will be the same as FFNO. For the memory layer, we choose an S4 layer <cit.> with a state dimension of 64 and a diagonal S4 (S4D) kernel.[The S4 repository has two available kernels, the diagonal S4 (S4D) and the Normal Plus Low Rank S4 (S4NPLR). In our experiments, we didn't find a significant difference between the two, and chose S4D for simplicity.]



For all our models, we use a simple spatial positional encoding E. In 1D, if the grid has f equispaced points in [0,L], then E ∈ℝ^f and the positional encoding is defined as E_i= i/L for 0≤ i ≤ f-1. In 2D, if we have a f× f points 2D equispaced grid in [0,L_x]×[0,L_y], the positional encoding is defined as E_ij = (i/L_x, j/L_y). In the input lifting operator ℛ_in, the input and the grid are stacked and a Linear layer from ℝ^2 →ℝ^h_0 is applied element-wise. For the decoder ℛ_out, we use another Linear layer (without positional encoding). 





§ BURGERS' EQUATION (1D): A STUDY ON LOW-RESOLUTION


The Burgers' equation with viscosity ν is a nonlinear PDE used as a modeling tool in fluid mechanics, traffic flow, and shock waves analysis. It encapsulates both diffusion and advection processes, making it essential for studying wave propagation and other dynamic phenomena. It is known for exhibiting a rich variety of behaviors, including the formation of shock waves and the transition from laminar to turbulent flow. The viscous Burgers' equation is written as:

    u_t + uu_x = ν u_xx

We used the publicly available dataset of the Burgers' equation in the PDEBench repository (<cit.>) with viscosity e^-3, which is available at resolution 1024. We compare our models at resolutions 64, 128, 256, 512 and 1024 and show the results at figure <ref>. As in the case of KS, s4FFNO outperforms FFNO, especially at low resolutions. Furthermore, we show the differnece in nRMSE at each timestep in figure <ref>. We observe that at the first timestep there is no difference between the two models - this makes sense because s4FFNO has the exact same architecture as FFNO for the first timestep. Yet as the initial condition is rolled out, there is more history of the trajectory and the difference between FFNO and s4FFNO increases.






 §.§ Relationship with fraction of unobserved information

As mentioned in Section <ref>
we measure the correlation of 
ω_f defined in (<ref>) with the difference in the nRMSE 
between FFNO and s4FFNO. The results can be seen in Figure <ref>





§ DATA GENERATION




 §.§ Kuramoto–Sivashinsky equation

For Kuramoto-Sivashinsky (KS) equation is given by: 

    u_t + uu_x + u_xx + ν u_xxxx   = 0    (t,x) ∈ [0,T]×[0,L] 
    
        u(0,x)     = u_0(x)    x ∈ [0,L]



We use periodic boundary conditions. Our data generation method is very similar to the one of <cit.>. We employ the method of lines <cit.>, where the spatial dimension is discretized, and the PDE is transformed to a systems of Ordinary Differential Equations (ODEs), one per point in the grid. In order to compute the spatial derivative of the solution at each point in the grid, a pseudospectral method is used, where derivatives are computed in frequency space and then converted to the original space through a Fast Fourier Transform. This method is implemented in the  method of the  package <cit.>. Similarly, the system of ODEs is solved numerically with a implicit Runge-Kutta method of the Radau IIA family of order 5 <cit.>, which is implemented in the  method of . We refer to the code provided in <cit.> to reproduce this data generation, however certain small modifications have to be made, like using a fixed Δ t per trajectory and increasing the number of modes in the initial condition.

As for the PDE parameters, we use L=64 and T=5. For the initial condition, we use a superposition of sinusodial waves: 

    u_0(x) = ∑_i=0^20 A_i sin(2π k_i/L  x + ϕ_i)

where for each trajectory, the A_i are sampled from a continuous uniform in [-0.5,0.5], the k_i are sampled from a discrete uniform in {1,2,...,8}, and the ϕ_i are sampled from a uniform uniform in [0,2π].  We discretize [0,T] into 26 equispaced points separated by Δ t=0.2[
We reiterate that, as opposed to <cit.>, we don't have a random Δ t per trajectory], keep the 25 temporal points generated by the numerical solver and discard our initial condition.
In the experiments in section <ref>, for each of the four values of the viscosity (0.15, 0.125, 0.1, 0.075), we generated a dataset with spatial resolution 512 with 2048 training samples and 256 test samples. For the experiment in the sequential model ablation in section <ref>, we generated one dataset with viscosity 0.15 in resolution 256, 4096 training samples and 256 test samples.



 §.§ Burgers' 1D equation

The 1D Burgers' equation can be written as:

    u_t + uu_x = ν u_xx   (t,x)∈ [0,T]×[0,L]

For the Burgers' equation, we take the publicly available Burgers' dataset of PDEBench <cit.> with viscosity 0.001. Out of the 10000 samples of the dataset, we use a 10% for testing. For training, we found it sufficient to use 2048 samples. Additionaly, for training and testing we only used the 20 first timesteps, since we observed that after the 20th timestep the diffusion term of the equation u_xx attenuates all high frequencies and the solution changes very slowly.



 §.§ Navier Stokes 2D equation

The incompressible Navier Stokes equation in the 2D unit torus is given by: 

    ∂ w(x,t)/∂ t + u(x,t) ·∇ w(x,t)    = νΔ w(x,t) + f(x),     x ∈ (0, 1)^2, t ∈ (0,T]
    
    ∇· u(x,t)    = 0,     x ∈ (0, 1)^2, t ∈ [0,T]
    
    
        w(x, 0)    = w_0(x),     x ∈ (0, 1)^2


For the data generation, we follow the method of <cit.>, yet with different temporal and spatial grids. The initial conditions w_0 are sampled from a Gaussian Random field 𝒩(0, 7^3/2(-Δ + 49I)^-2.5) with periodic boundary conditions. The forcing term is f(x_1,x_2) = 0.1 (sin2π(x_1+x_2) + cos2π(x_1+x_2)). At each timestep, the velocity is obtained from the vorticity by solving a Poisson equation. Then, spatial derivatives are obtained, and the non-linear term is computed in the physical space and then dealiased. A Crank-Nicholson scheme is used to move forward in time, with a timestep of 10^-4. We use a 512x512 spatial grid which is then downsampled to 64x64 for our experiments. For the viscosity ν=10^-3, we use a final time of 16 seconds and sample every 0.5 seconds. For the viscosity ν=10^-5, we use a final time of 3.2 seconds and sample every 0.1 seconds. For more details on the data generation algorithm, we refer to <cit.>.








§ TRAINING DETAILS


In this section, we will provide a detailed description of the training hyperparameters used in the KS experiments of Section <ref>, in the Burgers experimente of section <ref> and the Navier Stokes experiments of section <ref>. We start with the training hyperparameters. All our experiments used a learning rate of 0.001. For the number of epochs, in KS and Burgers, the training was done over 200 epochs with cosine annealing learning scheduling <cit.>; whereas in Navier Stokes we trained for 300 epochs and halved the learning rate every 90. As for the number of samples, KS and Burgers were trained with 2048 samples and Navier Stokes with 1024 samples. Lastly, we observed that the batch size was a sensitive hyperparameter for both the memory and memoryless models (it seemed to affect both equally) so we run a sweep at each experiment to select the best performing one. In the results shown in the paper, KS and Navier Stokes use a batch size of 32, and Burgers a batch size of 64. 

Another relevant detail is the memory length in training, that is, the number of past states that were fed to the memory layer in the MemNO model. In the KS and Burgers experiments, the maximum memory length was 25 (which is the same as the number of timesteps of the dataset). That means  that for the last timestep, the previous 24 states were fed into the memory layer. However, for GPU memory limitations in Navier Stokes the memory length was 16, half the number of timesteps of each trajectory in the dataset. In this case, the memory was reset after the 16th timestep, i.e. for the 16th timestep the 15 past states were fed to the memory model, yet for the 17th timestep only the 16th timestep was fed. Then, for the 18th timestep, the 17th and 16th were fed, and so on.

As in <cit.>, experiments were trained using teacher forcing. This means that for the prediction of the i-th timestep during training, the ground truth of the i-1 previous steps was fed to the model (as opposed to the prediction of the model for such steps).

We run our experiments on A6000/A6000-Ada GPUs. The Navier Stokes 2D experiments required around 34GB of GPU memory for the batch size of 32 and took around 5 hours to finish, whereas the rest of experiments in 1D required a lower GPU memory (less than 10GB) and each run normally took less than an hour.





§ ABLATIONS ON THE MEMORY LAYER


In this section we present two ablations regarding the memory layer of MemNO.


 §.§ Ablation: Choice of sequential model


In section <ref> we introduced MemNO as an architecture framework which allowed the introduction of memory through any choice of a sequential layer, which we chose as S4 in the previous experiments. In this section, we explore two other candidates for the sequential layers: a transformer and an LSTM. We introduce Transformer-FFNO (T-FFNO) and LSTM-FFNO as two models that are identical to s4FFNO except in the sequential layer, where a transformer and an LSTM are used respectively. The LSTM model only has one layer and the transformer layer includes causal masking and a positional encoding. The positional encoding for pos across the time dimension and i across the hidden dimension is given by:

    PE(pos, 2i)    = sin(pos/10000^2i/dim_model) 
    
        PE(pos, 2i + 1)    = cos(pos/10000^2i/dim_model)




We show results for the KS dataset with viscosity ν=0.15 and different resolutions. This dataset was generated using a resolution of 256 and contains 4096 samples, twice as many compared to the KS datasets of <ref>, given that transformers are known to perform better in high-data regimes. The results are shown in Figure <ref>. TFFNO performs significantly worse than s4FFNO across almost all resolutions, and even performs worse than FFNO. In constrast, LSTM-FFNO outperforms FFNO, which shows that MemNO can work with other sequentials model apart from S4. The memory term in Equation <ref> is a convolution in time, which is equivalent to the S4 layer and very similar to a Recurrent Neural Network (RNN) style layer, as showed in <cit.>. We believe that this inductive bias in the memory layer is the reason why both s4FFNO and LSTM-FFNO outperform FFNO. However, S4 was designed with a bias for continuous signals and has empirically has proven better performance in these kind of tasks <cit.>, which is in agreement with its increased performance over LSTMs in this experiment. Additionally, we observed that LSTMs were instable to train in Navier Stokes 2D datasets.

Lastly, we make two remarks. Firstly, we believe that transformers performed worse due to overfitting, given that the train losses were normally comparable or even smaller than the train losses of the rest of the models at each resolution. Modifications of the transformer model or to the training hyperparameters as in other works <cit.> might solve this issue. Secondly, recently there has been a surge of new sequential models such as Mamba <cit.>, RWQK <cit.>, xLSTM <cit.> or LRU <cit.>. We leave it as future work to study which of these sequential model has better overall performance, and hope that our study on the settings where the memory effect is relevant can help make accurate comparisons between them.







 §.§ Ablation: memory layer configuration


In Section <ref> we introduced the memory layer in MemNO as a single layer to be interleaved with Neural Operator layers. In our experiments, we inserted it after the second layer of a four layer Neural Operator. In this section, we explore the impact of having different layer configurations, including the possibility of having several memory layers. We will denote the configurations with a sequence of S and T letters. S means a Neural Operator layer (some sort of Spatial convolution), and T a memory layer (some sort of Time convolution). For example, SSTSS denotes the architecture of our experiments, were we have 2 Neural Operators layers, followed by a memory layer, followed by other 2 Neural Operator layers. Similarly, SSSST denotes 4 Neural Operators layers followed by a memory layer. In Table <ref>, we present the results for the KS dataset with ν=0.1 and final time of 4 seconds for several models. We include the s4FFNO model we used in previous experiments in the first row (with configuration SSTSS), and the FFNO model in the last row. In the middle rows, we show different configurations of memory and Neural Operator layers. It can be observed that all models with at least a memory layer outperform FFNO. There are slight differences between configurations, yet we focused mainly on the comparison to the memoryless model. For that reason, we fixed SSTSS configuration in our previous experiment, which was the most efficient (only one memory layer) and symmetric. We leave as further work determining if there are settings where a given configuration pattern can be substantially better than the rest.





§ APPENDIX: QUANTIFYING THE EFFECT OF MEMORY
 






We proceed to the (<ref>) first. Note that u_1(t), ∀ t ≥ 0 can be written as u_1(t) = a_0^(t)𝐞_0 + a^(t)_1 𝐞_1. Moreover, by Proposition <ref>, we have 

    ∂ a^(t)_0/∂ t   = 2B a^(t)_1
    ∂ a^(t)_1/∂ t   = a_1^(t) + B a_0^(t)

In matrix form, these equations form a linear matrix ODE: 

    ∂/∂ t[ a_0^(t); a_1^(t) ] = [  0 2B;  B  1 ][ a_0^(t); a_1^(t) ]

The solution of this ODE is given by [ a_0^(t); a_1^(t) ] = exp(t [  0 2B;  B  1 ])  [ a_0^(0); a_1^(0) ]. By the first statement of Lemma <ref> and the non-negativity of a_0^(0), a_1^(0), we get:

    a_0^(t)   ≤ 10 e^√(2)Bt(a_0^(0) + a_1^(0)),
     
    a_1^(t)   ≤ 10 e^√(2)Bt(a_0^(0) + a_1^(0))




We proceed to (<ref>). Note that for any s ≥ 0, we can write u_2(s) = â_0^(s)𝐞_0 + â_1^(s)𝐞_1 with â_0^(0) = a_0^(0) and â_1^(0) = a_1^(0). By Proposition <ref>, we have

    𝒬ℒ u_2(x) = B â_1^(s)𝐞_2(x)

Moreover, given a function v(x), 
the action of the operator exp𝒬ℒ(t̃) on v is given by the solution w(t̃, x) to the PDE 

    ∂/∂ t w(t,x)    = 𝒬ℒ w(t,x) 
     
        w(0,x)    = v(x)



If w(t,x) = ∑_n ∈ℕ_0 b_n^(t)𝐞_n and ∀ n ∈ℕ_0, b_n^(0)≥ 0, we are interested in solving the previous PDE with initial conditions b^(0)_2 = B â_1^(s) and b^(0)_n=0  ∀ n≠ 2.


 We claim that the coefficients â_n^(t)≥ 0 ∀ t>0 and ∀ n∈{0,1}. For t=0 this is by definition, and we will prove it for all t by way of contradiction. Suppose the claim is not true, then there exists a t^*>0, and some n^* ∈{0,1} such that â_n^*^(t^*)=0, and â_n^(s)>0 ∀ n∈{0,1} and ∀ s < t^*. 
But from continuity this implies that there exists 0<t'<t^* such that ∂/∂ tâ_n^*^(t')<0. However, it can be easy to see that if â_n^(s)>0  ∀ s≤ t', then 𝒫_1 ℒ u_2(t')>0 and 𝒫_1ℒ∫_0^t'exp𝒬ℒ(t-s)u_2(s) ds > 0. Therefore, from (<ref>), ∂/∂ tâ_n^*^(t')>0, which is a contradiction.

This claim implies that b_n^(0)≥0  ∀ n ∈ℕ, and in turn it implies that b_n^(t)≥ 0  ∀ n∈ℕ, t>0.
Applying 𝒬ℒ results in the following inequalities for the coefficients b_1^(t), b_2^(t), b_3^(t): 


    ∂/∂ t b^(t)_1    ≥   b^(t)_1 + B b^(t)_2 ≥ B b_2^(t)
    ∂/∂ t b^(t)_2    ≥  B b^(t)_1 + 4 b^(t)_2 + B b^(t)_3 ≥ B b^(t)_1 + B b^(t)_3
    ∂/∂ t b^(t)_3    ≥  B b^(t)_2 + 9 b^(t)_3 ≥  B b^(t)_2

Thus, we can write a linear matrix ODE for the vector (b_1^(t), b_2^(t), b_3^(t)):

    ∂/∂ t[ b_1^(t); b_2^(t); b_3^(t) ]≥[ 0 B 0; B 0 B; 0 B 0 ][ b_1^(t); b_2^(t); b_3^(t) ]


Therefore, using Lemma <ref>, for sufficiently large B we have b_2^(t-s)≥ B e^√(2)B(t-s)/10â_1^(s). 

Hence, if we write ∫_0^t exp𝒬ℒ(t-s)𝒬ℒ u_2(s) ds in the basis {𝐞_n}_n ∈ℕ_0, the coefficient for 𝐞_2 will be lower bounded by 

    ∫_0^t 1/10 B e^B (t-s) a^(s)_1 ds




Applying the second statement of Lemma <ref> and using the non-negativity of a^(0)_0 and a^(0)_1, we have â^(s)_1 ≥1/10  e^√(2)B s(a_0^(0) + a_1^(0)). Hence, the coefficient for 𝐞_2 is lower bounded by 

    ∫_0^t1/10 B e^√(2) B (t-s)1/10  e^√(2)B s(a_0^(0) + a_1^(0)) ds ≥B t/100  e^√(2)Bt (a_0^(0) + a_1^(0))


We finally need to consider what happens after applying the outermost operator 𝒫_1 ℒ. Because of Proposition <ref> again, applying ℒ makes the coefficient in front of 𝐞_1 at least B^2 t/100  e^√(2)Bt(a_0^(0) + a_1^(0)). Finally, applying 𝒫_1 preserves the coefficient in front of 𝐞_1. 

Hence, equation (<ref>) results in the following evolution inequalities:   

    ∂â^(t)_0/∂ t   ≥ 2B â^(t)_1
    ∂â^(t)_1/∂ t   ≥â_1^(t) + B â_0^(t) +  B^2 t/100  e^√(2)Bt(a_0^(0) + a_1^(0))

Using the second statement of Lemma <ref> again we have that  â_0(t) ≥1/10  e^√(2)B s(a_0^(0) + a_1^(0)). Thus, dropping the (positive) term â_1^(t) in equation <ref>, we have:

    ∂â^(t)_1/∂ t   ≥(1/10 + B t/100)B e^√(2)Bt(a_0^(0) + a_1^(0))


Integrating this equations yields:

    â^(t)_1 ≥ a^(0)_1 + 1/200 e^√(2) B t( √(2) B t + 10 √(2) - 1 )  (a_0^(0) + a_1^(0))

Thus, we have a_1^(t)≳ B t  e^√(2)Bt(a_0^(0) + a_1^(0)). Together with equation <ref>, the claim of the Theorem follows. 



 There exists B>0 sufficiently large such that for all t>0 the matrix [   0 2Bt;  Bt   t ] satisfies:

    ∀ i, j ∈{1,2}, exp([   0 2Bt;  Bt   t ])_i,j   ≤ 10 exp(√(2) Bt) 
    ∀ i, j ∈{1,2}, exp([   0 2Bt;  Bt   t ])_i,j   ≥1/10exp(√(2) Bt)




By direct calculation, we have: 


    exp([   0 2Bt;  Bt   t ])
       = 1/2 √(8B^2 + 1)[ √(8 B^2 + 1) g(B,t) - h(B,t)                    4B h(B,t);                    2B h(B,t) √(8 B^2 + 1) g(B,t) + h(B,t) ]

where:

    g(B,t) = e^1/2(√(8 B^2 + 1)+1)t + e^-1/2(√(8 B^2 + 1) - 1)t


    h(B,t) = e^1/2(√(8 B^2 + 1)+1)t - e^-1/2(√(8 B^2 + 1) - 1)t











Thus, the statement follows. 



 For all B > 0, the matrix [ 0 B 0; B 0 B; 0 B 0 ] satisfies:

    ∀ i, j ∈{1,2,3}, exp([ 0 B 0; B 0 B; 0 B 0 ])_i,j   ≥1/10exp(√(2) B)




By direct calculation:

    exp([ 0 B 0; B 0 B; 0 B 0 ])_i,j   = 
       1/4 e^-√(2) B[  2 e^√(2) B + e^2 √(2) B + 1       √(2) e^2 √(2) B - √(2) -2 e^√(2) B + e^2 √(2) B + 1;       √(2) e^2 √(2) B - √(2)           2 (e^2 √(2) B + 1)       √(2) e^2 √(2) B - √(2); -2 e^√(2) B + e^2 √(2) B + 1       √(2) e^2 √(2) B - √(2)  2 e^√(2) B + e^2 √(2) B + 1 ]

Thus, the statement follows.

























































































