Coarse-graining Langevin dynamics using reduced-order techniques
Lina Ma Xiantao Li Chun Liu

Stillman's conjecture via generic initial ideals
    Jan Draisma and Michał Lasoń and Anton Leykin
    December 30, 2023
=================================================


This paper considers the reduction of the Langevin equation arising from bio-molecular models. To facilitate the construction and implementation of the reduced models,  the problem is formulated as a reduced-order modeling problem. The reduced models can then be directly obtained from a Galerkin projection to appropriately defined Krylov subspaces.  The equivalence to  a moment-matching procedure, previously implemented in <cit.>,   is proved.  A particular emphasis is placed on 
the reduction of the stochastic noise, which is absent in many order-reduction problems. In particular, for order less than six we can show the reduced model obtained from the subspace projection automatically satisfies the fluctuation-dissipation theorem.
Details for the implementations, including a bi-orthogonalization procedure and the minimization of the number of matrix multiplications, will be discussed as well.

Stillman's conjecture via generic initial ideals
    Jan Draisma and Michał Lasoń and Anton Leykin
    December 30, 2023
=================================================





§ INTRODUCTION


Langevin dynamics models arise from a wide variety of problems, especially where a mechanical system is subject to random forces that can be modeled by white noise, e.g.,  as in <cit.>. 
A practical issue arises when the dimension of system is large, in which the computational cost can be overwhelming. For example, in bio-molecular models, the 
degrees of freedom are associated with the position and momentum of the constituting atoms, and the large dimensionality makes it difficult to probe large-scale biological processes over an extended period of time.   In this case, it is of great  interest to develop reduced models, which in bio-molecular modeling, is known as coarse-graining <cit.>. 

There are multiple benefits from such an approach. For example, reduced models can capture directly the dynamics of certain quantities of interest. Secondly, with the reduction of the dimension, the computational cost can be reduced dramatically.  In addition, the quantities of interest often correspond to slow variables. By eliminating fast variables, the time step can also be increased considerably. This allows one to access longer time scales <cit.>. 

There has been tremendous  recent progress in the development of coarse-grained models <cit.>. Most effort, however, is thermodynamics based. Namely, one aims to construct the free energy associated with the reduced variables, which then yields the driving force for the reduced dynamics, known as the potential of mean forces (PMF) <cit.>. As pointed out in <cit.>, the damping mechanics, which also plays an important role in the reduced dynamics, is not  part of the construction. 

In this work, we are interested in an equation-based derivation, where the reduced model can be derived directly from the Langevin dynamics.  Deriving reduced models from a  stochastic dynamical system has been a subject of extensive studies, the most well known of which is the homogenization approach <cit.>. Another important approach is to employ a coordinate transformation using normal forms to separate out the degrees of freedom that are less relevant <cit.>.  More recently, Legoll and Lelievre proposed to use conditional expectations to derived reduced models <cit.>.  Overall, these methods  require either significant scale separation assumption, or simple functions forms in the stochastic differential equations, which for bi-molecular models, does not apply.  For example, the force-field for biomolecular models typically involves complicated function forms.  

Meanwhile,  in the field  of molecular modeling there are also many methods that were proposed to coarse-grain a molecular dynamics model. Most of these methods are derived from a Hamiltonian system of ODEs <cit.>, either motivated by or directly  obtained, from the Mori-Zwanzig projection formalism <cit.>.  Strictly speaking, such a procedure will break down for stochastic models, due to the absence of the semi-group evolution operator.
 For Langevin dynamics, one empirical coarse-graining approach is the partition method  <cit.>, in which the variables are projected into appropriate subspaces.  However, the approach proposed in  <cit.> does not reduced the number of variables. Rather, it is a numerical integration algorithm.   The main reduction comes from filtering out high frequency modes in the numerical algorithm. In our previous work <cit.>, we have furthered this approach, by eliminating the fast-variables. This gives rise to a generalized Langevin equation (GLE) for the reduced variables. In principle, the GLE, under proper assumptions, is an exact model. 
After this reduction of the spatial dimensions, a temporal reduction was introduced to represent the memory term with a small number of auxiliary variable. Known as Markovian embedding, this procedure approximates the GLEs by using an extended system of stochastic differential equations (SDEs) with white noise.  The main idea is using a rational approximation for the Laplace transform, and the coefficients are determined based on a Hermite interpolation. The important advantage is that the approximation can be written as an extended system of SDE with no memory. 

A  well known issue in Padè type of approximations

is that when more conditions are incorporated, the resulting models tend to be ill-posed. In particular, the coefficient matrices are usually ill-conditioned, making it impractical.  Therefore, an important focus of this paper is on  re-formulating the coarse-graining procedure into a reduced-order problem, which has been widely studied <cit.>. 
In particular, we observe a feedback loop between the coarse-grain variables and the additional degrees of freedom, i.e., the fast variables. More specifically, the slow variables impose a mechanical force on the fast dynamics, and in turn, such influences will be propagated back as a force on the slow variables. As a result, the elimination of fast variables can be viewed as an order reduction problem, in that it is a large-dimensional dynamical system with low-dimensional input and low-dimensional output.  We will show that with an appropriate reformulation of the fast dynamics,  the transfer function from the order-reduction problem corresponds precisely to the memory kernel in the GLE. For such problems, one robust numerical method is the Krylov subspace projection <cit.>, which uses a Galerkin projection onto Krylov subspaces. The subspaces can be orthogonalized using the Lanczos algorithm <cit.>.  As a result, instead of mannually constructing the auxiliary system on a case-by-case basis as in the moment matching approach <cit.>, we can automate the procedure numerically.

More importantly, the bi-orthogonalization alleviate the problem of having ill-conditioned matrices.
 
For the current problem, the presence of the noise presents another critical issue. Namely, the random noise in the GLE must satisfy the second fluctuation-dissipation theorem (FDT) <cit.>, a necessary condition for the solution of the GLE to be stationary and to have the correct variance. In the Galerkin projection method, both the noise and the kernel function are being approximated. In general,  they do not satisfy the second FDT, unless the subspaces are properly selected.  We will provide two conditions that ensure such consistency, and we will show Krylov subspaces that fullfill these conditions. 













This paper is organized as follows. Section 2  describes the derivation the GLE system.   The classical approach of  approximating the  Laplace transform of the  memory kernel function θ(t) with a rational function will be presented. 
 Section 3 presents a formulation using the Galerkin  projection to general subspaces. Criteria will be provided in order to maintain the FDT in the reduced system. In Section 4, we introduce appropriate  Krylov subspaces to fulfill the criteria.  The resulting system will also be compared to a moment-matching procedure and the equivalence is proved in this section. Section 5 addresses two important issues in the numerical implementation.  Numerical examples  are shown in Section 6.  






§ MATHEMATICAL DERIVATION



 §.§ The Reduction of the Full Langevin Dynamics Model

We start with the full Langevin dynamics model with N atoms. After proper mass scaling <cit.>,  the system can be expressed as follows,

    {ẋ(t)=     v(t),
    v̇(t)=    F(x) - Γ v(t) + f(t),   
     .








where x=(x_1,x_2,…,x_N) denotes the displacement of all the atoms,  F(x) is the force derived from an empirical potentials V(x) with F=-∇ V, Γ denotes the damping coefficient for the friction term with dimension ℝ^3N× 3N, and f(t) is a stochastic force, usually modeled by a Gaussian white noise, which satisfies the fluctuation-dissipation theorem (FDT),

    f(t), f(t')^⊺=2k_B TΓδ(t-t').




For example, the random force can be written in the conventional form: df(t)= σ dW(t) with W(t) being the standard Brownian motion, and σσ^⊺= 2k_B TΓ.
Here, k_B is the Boltzmann constant, and T is the temperature of the system. This FDT is crucial to ensure that the system reaches the correct equilibrium state  <cit.> . 



















Implementing the full Langevin dynamics model can be very expensive, due to the large number of atoms involved in the entire system. 
Here we briefly go over a reduction procedure. More details can be found in   <cit.>.




The first step in the reduction procedure is to identify slow variables, which at the same time, are sufficient to describe the overall dynamics. In principle, these variables can be selected 
by transforming the system into normal forms  <cit.>.  For bio-molecules, a more intuitive and more efficient approach is based on the residues,  the building blocks of proteins, by choosing the center of mass of each amino acid.  Mathematically, this can be expressed as a small number of basis functions <cit.>, which span a subspace, denoted here by Y, with its orthogonal complement denoted by  Y^⊥.      Y has dimension m and Y^⊥ has dimension 3N-m: m≪ 3N.  We denote the basis vectors by {ϕ_i} and  {ψ_i}, respectively,  as follows,

    Y= span{ϕ_1,ϕ_2,…,ϕ_m},   
    Y^⊥ =  span{ψ_1,ψ_2,…,ψ_3N-m}.

Taking these basis vectors as columns and forming matrices Φ and Ψ,  one can  decompose the solution x in the following form,

    x(t)=Φ q(t) + Ψξ(t),

where q ∈ℝ^m and ξ∈ℝ^3N-m are nodal values associated with the basis vectors. Similarly, 

    v(t)=Φ p(t) + Ψη(t).
 



Meanwhile, a linearization of the force F≈ -Ax is considered, e.g., by principal component analysis (PCA)  <cit.>:

    x(t), x(t)^⊺ = k_BT A^-1,

which ensures that the covariance of the displacement is correct.

Now define the following projected matrices and vectors, 

    A_11 = Φ^⊺  A Φ,    A_12 = Φ^⊺  A Ψ,   Γ_11 = Φ^⊺ΓΦ,   Γ_12 = Φ^⊺ΓΨ,      f_1=Φ^⊺  f(t),


    A_21=Ψ^⊺  A Φ,    A_22=Ψ^⊺ A Ψ,   Γ_21 = Ψ^⊺ΓΦ,   Γ_22=Ψ^⊺ΓΨ,    f_2=Ψ^⊺  f(t).


By using this partition of variables, the original Langevin dynamics can be written in terms of the following first order stochastic differential equations (SDEs),

    {q̇(t)=      p(t), 
    ṗ(t)=   Φ^⊺  F(Φ q) - A_12ξ(t)  - Γ_11 p(t) - Γ_12η (t)+f_1(t),
       .
    {ξ̇(t)=    η(t), 
    η̇(t)=     -A_21 q(t)  - A_22ξ(t) - Γ_21 p(t) - Γ_22η(t) +  f_2(t).
       .

The linearization of the high-frequency modes has been based on numerous observations, e.g., <cit.>.
Essentially, we assume that the low frequency can be well captured by the basis functions in Φ, and the high frequency is nearly Gaussian. 
For example, in the rotation-translation block (RTB) approach, each residue is allowed to move as a rigid body. There is overwhelming evidence that the low-frequency normal modes are well represented by the subspace spanned by such basis functions <cit.>.  




Here (q, p) are the reduced/coarse-grained variables. Notice that the interactions involving the fast variables ξ have been 
linearized. By eliminating  (ξ, η), we have derived a low-dimensional reduced model <cit.>,

    {q̇(t)=      p(t), 
    ṗ(t)=    F_eff(q) - Γ_11 p (t)- ∫_0^t  θ(t-τ) p(τ)dτ + f(t).
    .

The effective force for the reduced system is

    F_eff(q)= Φ^⊺  F(Φ q) -A_12A_22^-1A_21q.


Compared to system (<ref>), the force F_ eff has an extra term -A_12A_22^-1A_21q from the derivation.

θ(t) is the memory kernel function, which is expressed in terms of a matrix exponential, 


    θ(t)= 
     [ A_12,  Γ_12] e^Dt[
    [ A_22^-1       0;       0      -I ]]
     [
    [ A_21; Γ_21 ]],

 where the matrix D ∈ℝ^(6N-2m)× (6N-2m)  is defined as,

    D=
    [
    [     0     I; -A_22 -Γ_22 ]] .


It has also been shown in  <cit.> 

    f =  f_1(t)
      - [ A_12,  Γ_12]  
     ∫_0^⊺  e^D(t-s)[
    [      0; f_2(s) ]] ds
     - [ A_12,  Γ_12]
    e^D t[
    [ ξ(0)+A_22^-1A_21q(0);                 η(0) ]]
    .

This random force  is a stationary Gaussian random process with mean zero, satisfying the second fluctuation-dissipation theorem:

    f(t) f(t')^⊺ = 2k_B T Γ_11δ(t-t') + k_B T θ(t-t').



Equation (<ref>) is known as the generalized Langevin equation (GLE). Currently there are primarily  three existing methods to solve the GLE numerically.
The first approach is to directly approximate the memory term, either by using quadrature formula, or by approximating the kernel function with a sum of exponentials. Known as the Prony sum, the later approach replaces the memory integral by additional variables that can be updated using certain recurrence formulas or by solving an ODEs system <cit.>.  The random noise can be approximated by introducing noises in those ODEs <cit.>. However, the approximation of the sum of exponentials requires the values of the kernel function (<ref>), which is difficult to compute due to the large dimensionality of the matrix D in the matrix exponential. The second approach is to eliminate the memory effect by approximating the kernel function with a delta function in time <cit.>. This approximation can be quite effective when the memory effect is not strong. But in general, the accuracy is quite limited.  The third approach is to approximate the memory effect by introducing auxiliary variables. This has been motivated by the Mori's continued-fraction approach <cit.>,  and has been pursued by many groups <cit.>.




 For example, in <cit.>, the  first order approximation leads to an extended dynamics with auxiliary variable z,

    {q̇(t)=     p(t), 
    ṗ(t)=   F_eff(q)- Γ_11 p(t) - z(t) +  f_1(t),
    ż(t)=     Bz(t) +  Cp(t) + ζ(t).
    .

The coefficients B and C can be found by using a  `moment matching' procedure, and we will elaborate on such procedures in section <ref>. At the same time, methods have been established to sample the additive noise ζ(t)  to ensure the FDT (<ref>). 

In theory, it is possible to advance to high order approximations using the above methods, e.g., a third order method  <cit.>. However, in practice, the matrices generated from the moment matching procedure tend to become ill-conditioned as the order of approximation increases. Moreover, the covariance of the noise and the covariance of the auxiliary variable z  need to be constructed specifically  for  each  order of approximation to ensure the FDT (<ref>), which is nontrivial. Therefore, it is  important to develop an alternative method to improve the robustness and automate the procedure. Inspired by order reduction methods for large-scale dynamical system, we will formulate the current problem as an order reduction problem with stochastic noise. The key  is to identify the low-dimension input and low-dimension output.



§ MODEL REDUCTION FOR THE  STOCHASTIC MODEL





 §.§ A Reformulation of the Orthogonal Dynamics

We will first introduce vector and matrix notations to rewrite the system (<ref>) in a more compact form.
Let y=(ξ, η)^⊺  represents   the partitioned variables, and u(t) = (q,p)^⊺  represents the coarse grained variables. System (<ref>)  can be rewritten as    

the following   SDEs:

    ẏ(t)= D  y (t)+  R  u(t) +g(t),     y(0) ∼N(0,k_B T Q).

with 

    R=
     [
     [     0     0; -A_21 -Γ_21 ]],
         g=[
     [      0; f_2(t) ]].



The matrix Q determines the initial covariance of y, given by,

    Q=
    [
    [ A_22^-1       0;       0       I ]].

Further we let Σ be the variance of the Gaussian noise g(t). It follows the Lyapunov equation, to ensure the stationarity of the solution,
  
    Σ = -k_BT ( DQ + Q D^⊺).

  It can be directly verified that,
 
    Σ=[
    [            0            0;            0 2 k_B TΓ_2,2 ]].

At the same time, we define

    L=[ A_12,  Γ_12],    R=[A_22^-1A_21, -Γ_21]^⊺.

Now the equation (<ref>) can be written as

    {q̇(t)=      p(t), 
    ṗ(t)=   Φ^⊺  F(Φ q)   - Γ_11 p(t) - L y+f_1(t).
       .



The corresponding memory kernel in (<ref>) is given by,

    θ(t)=L e^Dt R.


It is at this point that we recognize the similarity to an order reduction problem: The large-dimensional dynamics (<ref>) contains an input variable u(t), which is low-dimensional. Moreover, of direct importance to the coarse-grained dynamics (<ref>) is Ly, which again is low-dimensional. Also observed, however, is that the dimensions of L and R are different. Fortunately, we can reformulate the problem into the following equivalent dynamics (<ref>), where the input and output dimensions are the same.


    {q̇(t)    =p(t)
    ṗ(t)    =F_ eff(q)-Γ_11 p(t)-L  y + f_1(t),
    ẏ(t)    = D  y(t) + Rp(t) + g(t),      y(0) ∼N(0,k_B T Q).
    .


 





































Consider the following dynamics:

    ẏ_1(t)= D y_1(t) + Rp(t) +  g(t),     y_1(0) ∼ N(0, k_B TQ).

With a substitution into the first two equations in (<ref>) (in which y is replaced by y_1), one obtains a GLE that is equivalent to (<ref>).



































































Using a variation of constant formula, we find,

    y_1(t) = e^D t  y_1(0) + ∫_0^t  e^D(t-τ)Rp(τ) dτ
     +  ∫_0^t  e^D(t-τ)g(τ)dτ.

Next we define the out quantity w_1(t) from (<ref>),

    w_1(t) = L  y_1 = ∫_0^t θ(t-τ) p(τ) dτ+ Le^Dt y_1(0)+∫_0^t  Le^D(t-τ)g(τ) dτ
    =:  ∫_0^t  θ(t-τ) p(τ) dτ+  ζ(t).

Here ζ is the sum of the last two terms. For t>t', we have,

    ⟨ζ(t) ζ^⊺(t')⟩ = 
    k_BT Le^DtQe^D^⊺ t' L^⊺+ ⟨∫_0^t∫_0^t'Le^D(t-τ) g(τ) g^⊺(τ') e^D^⊺ (t-τ') L^⊺ dτ' dτ⟩ 
    =k_BTLe^D(t-t') QL^⊺.

The second step can be carried out by using the Itô's isometry.

Now we replace the term -A_12ξ(t) - Γ_12η (t) by w_1(t)  in system (<ref>). We  have,

    {q̇(t)=      p(t), 
    ṗ(t)=    F_eff(q)  - Γ_11 p(t) - ∫_0^t  θ(t-τ) p(τ) dτ- ζ(t) + f_1(t).
    .

Let f_1(t)=f_1(t)-ζ(t). With the assumption that the initial data of y_1 is uncorrelated with the noise term, we get, 

    ⟨  f_1(t)   f_1^⊺(t')⟩ =   
    2k_BT Γ_11δ(t-t')
    -2k_BTLe^D(t-t')[
    [    0; Γ_21 ]]
    +k_BTLe^D(t-t')QL^⊺
    
    =   2k_BT Γ_11δ(t-t')
    +k_BT θ(t-t').

The last step requires that

    QL^⊺ -2[
    [    0; Γ_21 ]] = R,

which can be easily verified. Now, according to theory of Gaussian processes <cit.>, the processes f(t) and f_1(t) are equivalent. 

Finally,  the memory terms in  (<ref>) with (<ref>) are the same,  the proof of equivalence is thus completed. 


It is clear that the dynamics (<ref>) is very similar to dynamics (<ref>), with subtle modification: p(t) instead of u(t) is involved in the system. More importantly,  in  (<ref>)   the input and the output of the dynamics  have the same dimension.  Our following discussion will be based on y_1, and instead, we will denote this term as y due to the equivalence.



We now have formulated the problem as a reduce-order problem: The dynamics of y involves a large-dimensional dynamical system, in which the variable p(t) is acting as a control variable. Meanwhile, what is of interest to the coarse-grained dynamics is the quantity L  y. As a result, we have at hand a large dynamical system with low-dimensional input and a low-dimensional output. 





 §.§ Properties of General Galerkin Projections

A remarkable success in order reduction problems is the Galerkin projection method to appropriately defined subspaces  <cit.>.
 Motivated by such success, 

we first consider a general Galerkin projection of the SDEs (<ref>),

    ẏ(t) = D  y(t) + Rp(t)+g(t).

More specifically, we  seek y(t) in the subspace X_n=span{ V_1, V_2, … V_n},  with each basis having m columns. We denote the space of test functions by X_n=span{W_1, …, W_n}. Now the projection can be stated as follows: find y(t)∈ X_n, such that for any  χ(t) ∈  X_n,

    (ẏ(t) - D   y(t) - Rp(t)-g(t), χ(t) )=0.


To put it in  a matrix-vector form, let V=[ V_1, V_2, … V_n] and W=[W_1, W_2, …, W_n], and we choose the columns as the basis for the two subspaces. The approximate solution is written as, 

    y(t) = V  z(t),

with z(t) being the nodal values. 
Then the Galerkin projection yields, 

    Mż(t) = D   z (t)+ W^⊺ R  p(t) + W^⊺ g(t),

where we have defined, 
    M=W^⊺ V,   D=W^⊺D V.

With the assumption that M is nonsingular, we can write

    ż(t) = M^-1D  z(t) + M^-1 W^⊺ R p(t) + f(t),

where 

    f(t)= M^-1 W^⊺ f(t),

 and its covariance matrix is given by,

    ⟨ f(t)  f(t')^⊺⟩ = Σδ(t-t'),   Σ=  M^-1 W^⊺Σ W M^-⊺.


With this reduction, we can now write down the reduced model  involving  the variables (p, q, z),


    q̇(t)=      p(t), 
    ṗ(t)=    F_eff(q) - Γ_11 p(t) - LV z(t) + f_1(t),
       
    
    ż(t) =    M^-1D  z(t) + M^-1 W^⊺ R p(t) + f(t),





In contrast to the conventional order  reduction problems  <cit.>, the current approach yields a noise term. Its presence  brings up an important issue: appropriate conditions are needed to ensure that the solution reaches  correct equilibrium, which will be addressed here. 

Due to ergodicity, the solution of the original SDE,  y(t), will evolve into a stationary process, and we expect the approximate solution to become a  stationary process as well. Assuming that the initial variance of z is k_BTQ,  that is, 
    ⟨ z(0) z(0)^⊺⟩ = k_BTQ,

then the stationarity implies that Q must satisfy the Lyapunov equation   <cit.>,

    k_BT(M^-1DQ + QD^⊺M^-⊺) = - Σ. Condition  A

This condition, as one of the necessary conditions to ensure the second FDT, will be referred to as Condition A.


Meanwhile, the projected dynamics (<ref>) corresponds to an approximation of the GLEs (<ref>). 
This can be verified by directly solving (<ref>), and then substitute L y into the equation for p.
With direct calculations, we find that  the approximated kernel can be expressed as,

    θ(t) ≈θ(t) :=L V e^M^-1D tM^-1 W^⊺ R.

Moreover, the low dimensional output is approximated by,

    w(t) ≈ w(t) = L y = ∫_0^T  θ(t-τ) p(τ) dτ + ζ(t),

where

    ζ(t) = LVe^M^-1Dt z(0) + ∫_0^T  LVe^M^-1D(t-τ) f(τ)dτ.

As a result, we obtain an approximate GLE model, 

    {q̇(t)=      p(t), 
    ṗ(t)=    F_eff(q) - Γ_11 p (t)- ∫_0^T θ(t-τ) p(τ)dτ + ζ(t) + f_1(t).
    .



The term ζ(t)  introduces an added Gaussian noise to the coarse-grained dynamics. 
 Together with the Lyapunov equation (<ref>), we can express its time  correlation as follows,

    ⟨ζ(t) ζ(t')^⊺⟩ = k_B T  L V e^M^-1D (t-t')Q V^⊺ L^⊺=k_BT θ(t-t'),  θ(t):=  L V e^M^-1D tQ V^⊺ L^⊺.



Clearly, in general the correlation of the noise θ(t) in (<ref>) might not be consistent with the memory kernel θ(t)  in (<ref>) and (<ref>). Namely, the second FDT, a necessary condition for the reduced model to have the correct statistics, may not be fulfilled.   The following theorem identifies the condition under which such consistency can be guaranteed. 



The coarse grained dynamics (<ref>) and (<ref>) derived from the 
Petrov-Galerkin projection will obey the second FDT (<ref>), if the following condition is satisfied:

    MQ V^⊺ L^⊺=W^⊺ Q L^⊺.




Recall that w(t)=Ly from (35b) needs to be injected into the dynamics of the reduced variables (35a).
The resulting random noise  is  ζ=-ζ(t) + f_1(t), with time correlation,

    ⟨ζ(t) ζ^⊺ (t') ⟩=2k_BT Γ_11δ(t-t') + 
     k_B T  L V e^M^-1D (t-t')Q V^⊺ L^⊺
     -2k_BT LV e^M^-1D (t-t')M^-1W^⊺[
    [    0; Γ_21 ]].

It is clear that if,

    Q V^⊺ L^⊺ - 2 M^-1W^⊺[
    [    0; Γ_21 ]] = M^-1W^⊺ R,

 this will result in the second FDT:

    ⟨ζ (t) ζ^⊺ (t') ⟩ = 2 k_BTΓ_11δ (t-t')+ k_BT θ(t-t').

In light of Equation (<ref>), Equation (<ref>) is equivalent to

    MQ V^⊺ L^⊺=W^⊺ Q L^⊺.
    Condition  B


This equation will be referred to as condition B.








Conditions A and B constitute the basis for constructing consistent stochastic reduced models. While condition A can be enforced by solving the Lyapunov equation,    condition B may not be satisfied by an arbitrary Galerkin projection. Therefore, we need to choose appropriate subspaces  for this to hold automatically. 




§ THE PROJECTION TO KRYLOV SUBSPACES





In this section, we will construct Krylov subspaces for the Galerkin projection procedure, which subsequently leads to approximations of the memory kernel function and 
random force.  


We will also discuss several issues related to the practical implementations. 

It turns out that the Krylov subspace approach has a close connection to a two-point Padè approximation, previously studied in   <cit.> to incorporate both long time and short time statistics. We will review this approach briefly, which will be referred to as moment-matching, and then make connections to the Krylov subspace projection approach.

 
We will consider the case where the damping coefficient is constant, i.e., Γ=γ I.



 §.§ The Moment Matching Approach

Define the moments,

    M_0= θ(0),   M_1= θ'(0),   ⋯,  M_ℓ=  θ^(ℓ)(0), ⋯,   M_∞=∫_0^∞θ(t)dt.

Notice that the moment M_∞ corresponds to the correlation time. 
With the moments, the memory function at t=0 can be expanded as:

    θ(t)=M_0+M_1t+M_2/2t^2+…+M_ℓ/ℓ!t^ℓ+….

Since the exact memory kernel is Le^DtR, it is clear that the moments are given by,
    M_0=LR,    M_1=LDR,  …   M_ℓ=LD^ℓ R,   M_∞ = - LD^-1R.


Meanwhile, the Laplace transform can be expanded near zero, 

    Θ(s)=M_0/s+M_1/s^2 + …+M_ℓ/s^ℓ+1+…,

which can be obtained by repeated integration by parts  <cit.>.




The moment matching procedure is essentially a rational approximation  for the Laplace transform of the memory kernel,

    Θ_n(s)=(s^nI - s^n-1 B_0 - s^n-2B_1-…-B_n-1)^-1(s^n-1 C_0 + s^n-2C_1+…+C_n-1),

such that,

    Θ_n(0)=Θ(0) (= M_∞),   θ_n^(ℓ)(0)=θ^(ℓ)(0) (=M_ℓ)   for   i=0,…, 2n-2,

To solve for the coefficients B_i, one needs to solve a linear system,

    [
    [   -M_∞    M_0      …  M_n-2;    M_0    M_1      …  M_n-1;      …;  M_n-2  M_n-1      … M_2n-3 ]]
    [
    [ B_n-1; B_n-2;     …;   B_0 ]]=
    [
    [  M_n-1;    M_n;      …; M_2n-2 ]]


We will use the second order approximation as an  example (n=2). In this case, the approximation would proceed as follows,


  * Set the Laplace transform of the approximated kernel to,

    Θ_2(s)=(s^2 - s B_0 - B_1)^-1(s C_0 + C_1).


  * Solve for the coefficients using the moments:

    [
    [ -M_∞  M_0;  M_0  M_1 ]]
    [ [ B_1; B_0 ]]=
    [[ M_1; M_2 ]],
       C_0=M_0,   C_1 = -B_1 M_∞.


  * The approximate kernel function in the real time domain can be expressed as:

    θ_2(t) ≈ [0    I]e^ Bt
     C,
       where  
     B = 
    [[   0 B_1;   I B_0 ]],     C=[ [ C_1; C_0 ]]
    .





Remark1 : Once B and C are computed, the variance of the random noise ζ in the stochastic equation ż= B z + C p + ζ(t), as well as the variance of the z(0), will be chosen based on these two matrices to satisfy the FDT. Such computation is quite involved in general.  Fortunately, as we will show, the subspace projection approach simplifies this effort considerably. 


Remark 2: Although one can increase the order of the approximation by simply introducing more moments, there remains an important practical problem, that is,  the condition number of the matrix in equation (<ref>) increases rapidly as the order increases. We hereby list the condition numbers in the following Table <ref> for a test problem.





We now turn to the Krylov subspace projections. 



 §.§ First Order Subspace Projection n=1
 

As the first approximation, 
we choose the subspaces  

    V=R, and   W=D^-⊺L^⊺.


We  show that the resulting approximate kernel function is the same as that from the moment matching approach.




By taking  V=R and  W=D^-⊺L^⊺ in the Galerkin projection, 
the memory kernel θ_1(t) in the projected dynamics (<ref>) 

is equivalent to that from the first order moment matching method.
 In particular,  two moments are matched exactly by the approximate kernel functions,

 


    θ(0) = M_0,   and∫_0^+∞θ_1(t) dt =M_∞.
















With direct computation, we get from (<ref>) that,

    θ_1 (0) = LVM^-1 W^⊺ R, 
      ∫_0^+∞θ_1(t) dt= -LVD^-1W^⊺ R.


By the particular choice of V and W (<ref>), we have M=W^⊺ R. Therefore,

    θ_1(0)=LR=M_0,  ∫_0^+∞θ_1(t) dt=W^⊺ R =LD^-1R=M_∞.





By taking  V=R and  W=D^-⊺L^⊺ in the Galerkin projection, 
the projected dynamics  (<ref>) will automatically satisfy the second FDT (<ref>).


 We need to show that <ref> is satisfied by this choice of W and V in this case.
Given  (<ref>) and <ref>, we have

    -W^⊺Σ W=-MΣM^-⊺= k_BT (DQM^⊺+ MQD^⊺).

Notice that since W=D^-⊺L^⊺,  one has D=LV. In addition, from Equation (<ref>),  we have

    LQW + W^⊺ Q L^⊺ = LVQM^⊺ + MQV^⊺ L^⊺.

It is clear that on both sides, it is a summation of a matrix and its transpose. By moving terms we find,

    LQW -LVQM^⊺ =   MQV^⊺ L^⊺ - W^⊺ Q L^⊺

and  <ref> would hold if either  side equals to zero. We will examine the two terms on the right hand side.  

Since Γ≡γ I, we have Γ_12=0. Further 
notice that,

    D^-1=[
    [ -A_22^-1Γ_22     -A_22^-1;            I            0 ]].

By direct calculations, 
the second term on the right hand side can be simplified to,

    W^⊺ Q L^⊺  = LD^-1QL^⊺=-A_12A_22^-1Γ_22A_22^-1A_21=-M_∞.
 

Regarding the first term on the right hand side of (<ref>), it can be directly verified that,

    M = LD^-1R = - M_∞,    D = LR=M_0,

which are both symmetric matrices. Further, by using Equation (<ref>), we get that,

    Σ = 2k_B TM_∞^-1.

Then the Lyapunov Equation (<ref>) becomes

    k_B T (M_∞^-1 M_0  Q +  Q M_0^⊺  M_∞^-1)=2k_BTM_∞^-1,
from which we obtain the solution Q=M_0^-1. 

Therefore the first term on the right hand side of (<ref>) becomes

    M  Q V^⊺ L^⊺=-M_∞ M_0^-1 R^⊺ L^⊺ = -M_∞,


which would cancel the second term and complete the proof.

























 §.§ Second Order Subspace Projection n=2








We now extend the subspace by choosing,

    V=[R, DR],  and  W=[D^-⊺ L^⊺, L^⊺].

 
 As a result, 
the two matrices M and D in the Galerkin formulation are given by,

    M =[[ - M_∞   M_0;   M_0   M_1 ]],   D =[[ M_0 M_1; M_1 M_2 ]].

It's easy to check that M^-1D= B, as in equation (<ref>). Within this extended approximation, the approximate memory function is given by,

    θ_2(t)= [M_0    M_1] e^ BtM^-1[[ -M_∞;  M_0 ]]
    =[M_0    M_1] e^ Bt[[ I; 0 ]].




We first show that this approximation is equivalent to the moment matching procedure.  It is straightforward to verify that the approximate kernel, denoted by θ_2, from the moment matching procedure, should satisfy the following second order differential equation:

    θ̈_2(t) = B_0θ̇_2(t) + B_1θ_2(t),   θ_2(0)=M_0,   θ̇_2(0)=M_1.


We now show that the kernel function θ_2(t) follows the same equation. Thanks to the uniqueness, we can then conclude the equivalence.  The key observation is that,

    [M_0    M_1] B = [M_1    M_2].

As a result, it can be quickly verified that

    B_1[M_0   M_1]+B_0[M_0   M_1] B=B_1[M_0    M_1]+B_0[M_1   M_2] = [M_1   M_2] B = [M_0   M_1] B^2,

which combined with (<ref>) would lead to 

    θ̈_2(t) = B_0θ̇_2(t) + B_1θ_2(t),   θ_2(0)=M_0,   θ̇_2(0)=M_1.












Therefore, we have this following theorem.

The reduced model    (<ref>) from the Galerkin projection with the choice of V=[R, DR], and W=[D^-⊺L^⊺, L^⊺] produces an approximate memory kernel function, which  is equivalent to  that from the second order moment matching procedure.

Furthermore, we have,

The projected system (<ref>) with the choice of V=[R, DR], and W=[D^-⊺L^⊺, L^⊺] will automatically satisfy the second FDT (<ref>).




We only need to justify <ref>. It is straightforward to show that,
 
    W^⊺ Q L^⊺ = 
      [
      [ LD^-1;     L ]]
      [
      [ A_22^-1       0;       0       I ]]
      [
      [ A_21;    0 ]] = W^⊺ R =  [
      [ -M_∞;  M_0 ]].

 With the choice of V,  we have
 
    V^⊺ L^⊺ =[[   M_0; M_1^⊺ ]]
     =[[ M_0;   0 ]].

 Notice M_1=0, which can be verified  by direct calculation.
Then by some direct calculations with  the representation of the covariance matrix, we have

    MΣ M^⊺ = W^⊺Σ W = [
    [ 2k_BT M_∞         0;         0         0 ]]
    =- k_B T ( D  Q  M^⊺ +  M  Q  D^⊺).

Meanwhile, we have, 

    M  =[ [ -M_∞  M_0;  M_0    0 ]],

which gives,

    M^-1  =[ [               0          M_0^-1;          M_0^-1 M_0^-1M_∞M_0^-1 ]],  Σ=2k_BT[[               0               0;               0 M_0^-1M_∞M_0^-1 ]].

Now we  solve the Lyapunov equation and we find that,

    Q = [[   M_0^-1        0;        0 - M_2^-1 ]].



With Q available,  it can be verified that 

    M  Q V^⊺ L^⊺ = [[ M_0;   0 ]] = W^⊺ Q L^⊺,

which is our condition B, thus it completes the proof. 












 §.§ Generalization to High Order Approximation (n ≥ 2)




Inspired by the previous choices, we consider

    V=[R, DR, …, D^n-1R],   W=[D^-⊺L^⊺, L^⊺,D^⊺ L^⊺,…,(D^⊺)^n-2L^⊺],
 
and apply Galerkin projection to the two subspaces generated by the columns of these two matrices.

 The  corresponding matrices M, D, B and W^⊺ R are given by, respectively,

    M=   [
    [   -M_∞    M_0      …  M_n-2;    M_0    M_1      …  M_n-1;      ⋮;  M_n-2      … M_2n-3 ]],
       D=[
    [    M_0    M_1      …  M_n-1;    M_1    M_2      …    M_n;      ⋮;  M_n-1      … M_2n-2 ]],
    
    
    
     B = M^-1 D=   [
    [     0     0     …     0 B_n-1;     I     0     …     0 B_n-2;     0     I     …     0;     ⋮;     0     0     …     I   B_0 ]],  
    W^⊺ R = [[  -M_∞;   M_0;     ⋮; M_n-2 ]].

Therefore, the approximate kernel under the Galerkin projection can be expressed as,

    θ_n(t)= [M_0    M_1  …   M_n-1 ] e^ Bt[[ I; 0; ⋮; 0 ]].


Meanwhile, the high order approximate memory kernel from the moment matching procedure satisfies the n^th order differential equation:

    θ^(n)_n(t)=B_0 θ^(n-1)_n(t)+B_1θ^(n-2)_n(t)+… B_n-1θ_n(t),   θ_n(0)=M_0,   …,  θ^(n-1)_n(0)=M_n-1.



We first show that these approximate kernel functions are the same.

The function θ_n(t) in equation (<ref>) is equivalent to the function θ_n(t) generated from moment matching procedure as described in section 4.1. In particular, it also satisfies the initial-value problem,

    θ^(n)_n(t)=B_0 θ^(n-1)_n(t)+B_1θ^(n-2)_n(t)+… B_n-1θ_n(t),   θ_n(0)=M_0,   …,  θ^(n-1)_n(0)=M_n-1.







Each M_i is a m by m matrix, and the dimension of   θ_n(t) is also m× m.
For simpler notations, we will denote [M_i   M_i+1  …   M_i+n-1]=G_i.
If  we can show that

    G_0 B^n = B_0G_0 B^n-1 + B_1G_0 B^n-2 + … + B_n-1 G_0,

this will prove θ_n(t) satisfies the same differential equation.
Notice that the recursive relation 
    G_i  B=G_i+1,   for   i=0,… n-2,
 comes straightforward since M  B =  D.


Then it remains to check that

    G_n-1 B = B_0 G_n-1+B_1 G_n-2 + … +B_n-1G_0.

We will take a closer look at each block elements. The first block on the left hand side is M_n, and on the right hand side, we have ∑_i=0^n-1M_iB_n-1-i. They are equal  due to the equation M  B =  D. In fact, all other blocks except the last one can be shown from the same equation.

 The last block automatically equal to each other since they have exactly the same representation. 
 
 For the initial conditions, they can be easily verified using  M  B =  D.


What we will study next is whether this formulation also obeys the second FDT. However, we are not able to prove the general case due to the lengthy calculations involved. We are able to prove the consistency for  n ≤ 5. The following few results are useful for the verification. Numerical tests suggest that the consistency holds also for higher order cases.

 
The moments of the memory function are all symmetric matrices. As a result, M, and D as defined in equation (<ref>) are also symmetric matrices. 


 We only need to show all moments M_i are symmetric. 
Recall that 

    D=[[     0     I; -A_22 -Γ_22 ] = 
    ]
    =[[    0   -I;   -I Γ_22 ]][[ A_22    0;    0   -I ]],

and 

    R=[[ A_22^-1 A_21;        -Γ_21 ]]=[[ A_22^-1       0;       0      -I ]][[ A_21; Γ_21 ]].

Therefore for i>0,

    M_i=LD^iR =   [ [A_12 Γ_12] ][[    0   -I;   -I Γ_22 ]]
    [[ A_22    0;    0   -I ]]⋯[[    0   -I;   -I Γ_22 ]]
    [[ A_22    0;    0   -I ]]
    [[ A_22^-1       0;       0      -I ]]
    [[ A_21; Γ_21 ]] 
    
    =    [ [A_12 Γ_12] ][[    0   -I;   -I Γ_22 ]]
    [[ A_22    0;    0   -I ]]⋯[[    0   -I;   -I Γ_22 ]]
    [
    [ A_21; Γ_21 ]],

which is clearly symmetric.
At the same time, it is straightforward to see that M_0 is symmetric by direct calculation. Finally, 

    M_∞   =[ [A_12 Γ_12] ][
    [ -A_22^-1Γ_22     -A_22^-1;            I            0 ]][[ A_22^-1       0;       0      -I ]][[ A_21; Γ_21 ]]
       =[ [A_12 Γ_12] ][
    [ -A_22^-1Γ_22A_22^-1             A_22^-1;             A_22^-1                   0 ]]
    [[ A_21; Γ_21 ]],

 is symmetric as well.



Assume that M is invertible. 


Condition B is equivalent to, 

    Q [[   M_0;   M_1;     ⋮; M_n-1 ]]=[[ I; 0; ⋮; 0 ]].





When Γ=γ I, Γ_12=0, we have the identity QL^⊺ =R. It is also easy to see that due to symmetry from Lemma 1, one has,  

    V^⊺ L^⊺ = [[   M_0^⊺;   M_1^⊺;       …; M_n-1^⊺ ]]= [[   M_0;   M_1;     …; M_n-1 ]].

Therefore, <ref> becomes

    M Q [[   M_0;   M_1;     …; M_n-1 ]]= W^⊺ R=
    [[  -M_∞;   M_0;     …; M_n-2 ]].

Multiplying both sides by M^-1 (with the assumption that M is invertible), we arrive at equation (<ref>). 







Let Σ = W^⊺Σ W, which has dimension nm × nm. If it is partitioned into a block matrix with each block having dimension m× m , then the  block elements have the following recurrence relations:

    Σ_1,1 = 2M_∞,  Σ_i,2=0,   Σ_2,i=0,
       Σ_ij = -1/γΣ_i+1,j- 1/γΣ_i,j+1-2k_BT M_i+j-3,    i,j >2.



As a result, the elements of Σ can be constructed column by column using the recurrence relation. They can be expressed in terms of the moments M_is. The next lemma shows that the moments also exhibit a recurrence relation, which can be exploited to make the calculation a bit easier.


The moments M_i=LD^iR can be written as a linear combination of matrices A_12A_22^kA_21,

    M_i = ∑_k=0^⌊i/2⌋-1 c_i,k A_12 A_22^k A_21.








The proof of these lemmas can be found in the appendix.



The reduced system  (<ref>) from the  Petrov-Galerkin projection obeys the second FDT for orders   n≤ 5.


It now becomes clear that in order to check wether  the second FDT holds for high order approximation system, one only needs to show equation (<ref>). On the other hand, we know Q is the solution to Lyapunov equation (<ref>), which uniquely determined. Therefore under the assumption that D is nonsingular, Q= Q  D is also uniquely determined. This also leads to the following equation based on the fact that M is symmetric.

    k_B T ( Q M +  M  Q^⊺) = - M Σ M = -W^⊺Σ W=-Σ.

Now the goal is to compute the exact form of Q. We will present the expression of Q and W^⊺Σ W for n=3, 4, 5, and readers can substitute those forms into the equation above to verify. There are some identities needed in order to complete the verification, which we will present in the Appendix. 


For n=3,

    Σ= k_B T[
    [    2M_∞       0 -2γ M_0;       0       0       0; -2γ M_0       0 -2γ M_2 ]],
       Q=[[    I    0    0;    0   -I    0;    0 2γ I    I ]].



For n=4,

    Σ = k_B T[
    [     2M_∞        0  -2γ M_0 2γ^2 M_0;        0        0        0        0;  -2γ M_0        0  -2γ M_2  -2γ M_3;  2γ^2M_0        0  -2γ M_3  2γ^2M_3;          ]],
       Q=[[       I       0       0       0;       0      -I       0       0;       0    2γ I       I       0;       0 -2γ^2 I   -2γ I      -I ]].


For n=5,

    Σ = k_B T[
    [                          2M_∞                             0                       -2γ M_0                      2γ^2 M_0            -2γ^3 M_0 - 2γ M_2;                             0                             0                             0                             0                             0;                       -2γ M_0                             0                       -2γ M_2                       -2γ M_3                       2γ^2M_3;                       2γ^2M_0                             0                       -2γ M_3                       2γ^2M_3                     2 γ^2 M_4;           -2γ^3 M_0 - 2 γ M_2                             0                       -2γ M_4                       2γ^2M_4 -2γ^3 M_4 - 2γ^2 M_5 - 2γ M_6 ]],
    
        Q=[[       I       0       0       0       0;       0      -I       0       0       0;       0    2γ I       I       0       0;       0 -2γ^2 I   -2γ I      -I       0;       0  2γ^3 I  4γ^2 I    4γ I       I ]].






§ NUMERICAL IMPLEMENTATION


In this section, we will describe  the numerical implementation of the Krylov subspace projection method.


In the previous section, we have studied properties of the projected dynamics with particular choices of   V and W. 
However,  as is well known <cit.>, a direct implementation using those matrices often leads to ill-conditioned matrices. This has clearly been shown in Table 1. 
A much more robust approach is to obtain orthogonal basis by using appropriate orthogonalization algorithms.



Let us first introduce the notations for these two Krylov subspaces for an nth order approximation. 

    𝒦_n(D, R) = span{R, DR, …, D^n-1 R},   𝒦_n(D^⊺,D^-⊺ L^⊺)=span{D^-⊺L^⊺, L^⊺, …, (D^n-2)^⊺ L^⊺}.




 §.§ Block Lanczos Algorithms (BLBIO)

 We will adopt  the non-symmetric block Lanczos algorithms from  <cit.>
 to generate orthogonal basis V=[ V_1, …, V_n] and W=[ W_1, …, W_n] for 𝒦_n(D, R) and 𝒦_n(D^⊺,D^-⊺ L^⊺), respectively. 

The Lanczos algorithm proceeds as follows. 
Choose V_1=R, W_1= D^-⊺ L^⊺, and let δ_1 = W_1^⊺ V_1, and for k=1, 2, … compute 


    δ_k^A = W_k ^⊺ D V_k,
       α_k = δ_k^-1δ_k^A,  α_k=δ_k^-⊺(δ_k^A)^⊺,
       β_k-1=δ_k-1^-1γ_k-1^⊺δ_k,  β_k-1=δ_k-1^-⊺γ_k-1^⊺δ_k^⊺,   ( if  n>0)
       V_tmp=DV_k - V_kα_k-V_k-1β_k-1,   
     W_tmp=D^⊺ W_k-W_k α_k - W_k-1β_k-1,
       δ_tmp=W_tmp^⊺ V_tmp
        choose γ_k, γ_k  and δ_k+1,    s.t.  γ_k^⊺δ_k+1γ_k = δ_tmp


Several possible choices have been recommended in  <cit.> for γ_k, γ_kand δ_k+1.  We found that the QR factorization with column pivoting for V_tmp and W_tmp is quite robust. Namely,  


    V_tmp P=UR,   W_tmp P =  U  R.

Then we choose 

    V_k+1=U,    W_k+1= U,  γ_k=RP^⊺,  γ_k= R  P^⊺.


By following this algorithm, we obtain the orthogonality properties among the basis vectors of the Krylov subspaces. In particular, the matrix M is diagonal, and the matrix D is block-tridiagonal. 









As a result, the SDEs for the auxiliary variable z (<ref>) involves sparse matrices.







 §.§ Implementation without Ψ

The implementation of the algorithm requires the matrices L, V, W, and R,  all involving the Ψ matrix as part of the construction. Constructing Ψ  is usually not feasible for large systems. Here we present an algorithm that does not involve   Ψ.



Let's first derive a few useful identities involving the Ψ matrix.
We start with,

    [[  Φ^⊺; Ψ ^⊺ ]] A [Φ  Ψ]=
    [[ A_11 A_12; A_21 A_22 ]],
        ([[  Φ^⊺; Ψ ^⊺ ]] A [Φ  Ψ])^-1
    =[[  Φ^⊺; Ψ ^⊺ ]] A^-1 [Φ  Ψ].


Using a block inversion formula, we get

    (Φ^⊺ A^-1Φ)^-1 = A_11 - A_12A_22^-1A_21.

By left multiplying the equation by A^-1Φ, together with the identity ΦΦ^⊺ = I - ΨΨ^⊺,  we find that,

    A^-1Φ(Φ^⊺ A^-1Φ)^-1 = Φ - Ψ A_22^-1 A_21.

Next, right multiplying the above equation by Φ^⊺ A^-1, we arrive at,

    A^-1Φ(Φ^⊺ A^-1Φ)^-1Φ^⊺ A^-1= ΦΦ^⊺ A^-1 - Ψ A_22^-1Ψ^⊺ A ΦΦ^⊺ A^-1 = ΦΦ^⊺ A^-1 - Ψ A_22^-1Ψ^⊺+ (I-ΦΦ^⊺)A^-1.



Now we define, 
    D = [[ Ψ 0; 0 Ψ ]] D [[ Ψ^⊺   0;   0 Ψ^⊺ ]],      R = [[ Ψ 0; 0 Ψ ]] R,    L =Φ^⊺[ A,  Γ] [[ Ψ 0; 0 Ψ ]]  D^-1[[ Ψ^⊺   0;   0 Ψ^⊺ ]]


We start with the following observation,
 
The following relation holds between the two Krylov subspaces, 𝒦_n(D,R) and 𝒦_n( D,  R):

    [[ Ψ 0; 0 Ψ ]] 
    𝒦_n(D,R) = 
    
    𝒦_n( D,  R).

Similarly,

    [[ Ψ 0; 0 Ψ ]] 𝒦_n (D^⊺ , D^-⊺L^⊺)
    =𝒦_n(  D^⊺, L^⊺).



With these observations, we show that: 







The Lanczos algorithm, the Galerkin projection, and  the sampling of the noise, can be done without Ψ.





First it can be directly shown that, 
 
    D = [ [            0         ΨΨ^⊺; -ΨΨ^⊺ A ΨΨ^⊺   -ΨΨ^⊺ΓΨΨ^⊺ ]].

Thanks again to the identity 

    ΨΨ^⊺ =  I - ΦΦ^⊺,

 we  can evaluate ΨΨ^⊺ through the matrix Φ. Therefore the calculation of D can be done without Ψ.
 
Secondly, to compute L, we notice that the terms involving Ψ are Ψ A_22^-1Γ_22Ψ^⊺, ΨΨ^⊺, and Ψ A_22^-1Ψ^⊺, and these terms can be represented without Ψ from Equation (<ref>) and (<ref>).  The calculation of R is similar.
 
 Thirdly, we see that the solution of the projected dynamics (<ref>) enters the coarse-grained dynamics (<ref>) via a matrix multiplication by LV.
It is straightforward to write L as 
    L=Φ^⊺[ A,  Γ][[ Ψ 0; 0 Ψ ]],
 which means that for the term LV, we can actually compute  Φ^⊺[ A,  Γ] V, where V is constructed using  using the block Lanczos from space  𝒦_n( D,  R).
 
 Now let V and 
 W be the basis generated from the orthogonalization of the new Krylov subspaces 𝒦_n( D,  R) and 𝒦_n( D^⊺,  L^⊺ ), respectively.
 Therefore, the matrices  M = W^⊺ V =  W^⊺ V, D = W^⊺ D V= W^⊺ D  V  and W^⊺R= W^⊺ R can all be generated without introducing Ψ.

 Finally, it remains to show that the sampling of the noise does not have to involve Ψ, which is clearly true since the noise is represented as M^-1 W^⊺[[      0; Ψ^⊺  f ]], and W^⊺Ψ^⊺ =  W^⊺. 



 It is a trivial, but important point in practice that 
 in the numerical implementation, it is not necessary to store the full matrix ΦΦ^⊺. For a given vector u, the multiplication ΦΦ^⊺ u can be done through Φ (Φ^⊺ u).



 §.§ A Summary of the Galerkin Projection

The Galerkin projection method can be summarized as follows,


  * Choose appropriate basis matrix Φ. 

  * Pick the order of approximation n ≥ 1. Use the block-Lanczos algorithm to determine the orthogonal basis V and W,
 for the Krylov subspaces 𝒦_n(D, R) and 𝒦_n(D^⊺,D^-⊺ L^⊺), respectively.

  * Solve the stochastic differential equations (<ref>), where M,  D,  f are defined from equations (<ref>) and (<ref>).
The initial variance of z(t) is determine from <ref>.

Clearly, this procedure avoided manual constructions of the reduced model. This choice of the Krylov subspaces guarantees that the FDT is satisfied (through <ref>), at least till the fifth order of approximation (numerical tests indicate that this is true for higher order cases). 




§ NUMERICAL TEST

We test our algorithm on the example considered in  <cit.>. We simulate the dynamics of the protein Chignolin (PDB id 1uao) at temperature T=298 for .4 nano seconds. 
 The system is set up in solvation, modeled by the generalized Born (GB) model and simulations have been conducted  in TINKER <cit.> using force field CHARMM22.
For the surrounding bath, we considered the case γ=91 ps^-1 which corresponds to water solvant <cit.> and a  low friction  case γ=5 ps^-1. 
In the latter case,  the kernel function exhibits nontrivial behavior <cit.>: it tends to be more oscillatory  compared to the former case. 
By calculating the eigenvalues of A we have identified the under-damped regime to be γ < 13.4 and the over-damped regime to γ > 997.7.   

  Data are collected to compute the PCA matrix A=k_BT⟨ x, x^⊺⟩^-1. The projection matrix are composed of RTB basis <cit.>, since there are 10 residues in Chignolin, the dimension of the coarse-grained variables is 60. The explicit forms of the basis functions in Φ each each translational and rotational mode can be found in <cit.>. 
 



 We first  present the numerical result for γ = 91 in Figure <ref>. 
 On the left panel, we showed the comparison of approximating memory function, from order two to order seven. The right panel of the figure provides the comparison of time correlation of the momentum. Both exact plots are obtained  by running the full model. The order of approximations, n,  is defined as the order of Krylov subspaces, which is equivalent to the order of the rational functions in the moment matching approach. Since the kernel function θ(t) is matrix-valued, we chose the sixth diagonal, θ_6,6(t) for the comparison, this index corresponds to the third rotational component of the first residue. We can observe the approximation is satisfactory for n≤ 5. 


In Figure <ref>, we present a comparison for γ = 5. 

The small damping constant leads to a underdamped system, making  the approximation difficult due to the rapid and non-trivial oscillation. However we can observe substantial improvement of  the accuracy on the memory kernel. The memory effect on auto correlation is evident compared to system with high damping constant. 

Though improvement is significant for the memory kernel,  the velocity time correlation exhibits noticeable error. 








In Figure <ref>, we provide a close-up view over the time interval [0,0.2] ps, and show results from secon order to seventh order approximations.  We observe  increased accuracy as the order of the approximation is increased within this time period. 









Finally, we present  the relative L_2 error for both memory kernel and time correlation, comparing the results of second order and seventh order, for both γ=91 and γ=5 in Figure <ref>. This relative L_2 error is  computed for the time period [0,1]. We showed error for each coarse grained variables, and improvement of accuracy is significant. 




In addition to the Krylov subspaces that were presented in the previous section, we also implemented inverse Krylov subspaces and shifted-inverse Krylov subspaces in the Galerkin projection. These variations can often offer better approximations to the transfer function in order-reduction problems <cit.>, which in our case, corresponds to the memory function. However, through our numerical computations, we found that none of these choices satisfies (<ref>). This implies that the second FDT is not fulfilled, and the reduced dynamics (<ref>) does not produce stationary processes <cit.>. In fact, the variance of the solution will follow the dynamic Lyapunov equation (<cit.> Eqn 3.103), and it will not converge to the steady-state.  























§ CONCLUSION

We adopted reduced-order modeling techniques to reduce the Langevin dynamics model.  We consider reduced models obtained from Petrov-Galerkin projections. By selecting appropriate Krylov spaces,  we  show the mathematical equivalence of the proposed model to the reduced models derived from moment matching procedure. Another emphasis is placed on the statistical consistency, i.e., the fluctuation-dissipation theorem. We are able to identify two conditions that ensure such consistency. We also showed that the Galerkin projections to the selected subspaces automatically satisfy the FDT, at least for n≤ 5. With the block Lanczos algorithm, the models derived this way are more robust. 

One open issue is the case when the damping coefficient Γ is not proportional to an identity matrix. In this case, both condition A and condition B are still sufficient to ensure the FDT. But the Krylov subspaces construction in section 4 may not satisfy condition B. Another open question is whether one can bypass the linear approximations used in (<ref>) and (<ref>). It seems that a different methodology is needed to derive the generalized Langevin equation (<ref>). These issues will be addressed in future works.

 

§ ACKNOWLEDGEMENT

 The research of Li is supported by  NSF Grant DMS-1522617 and DMS-1619661.  The research of Liu is supported by NSF grant NSF Grant DMS-1759535 and DMS-1759536.




§ RECURRENCE FORMULA FOR Σ=W^⊺Σ W

This is the proof of lemma 3. 
For the nth approximation, 

    Σ = W^⊺Σ W = [
    [        LD^-1Σ D^-⊺L^⊺            LD^-1Σ L^⊺        LD^-1Σ D^⊺ L^⊺                     …  LD^-1Σ (D^n-2)^⊺ L^⊺;           LΣ D^-⊺ L^⊺                LΣ L^⊺                     …       LΣ(D^n-2)^⊺ L^⊺;                     …;      LD^n-2Σ D^-T L^⊺                     … LD^n-2Σ (D^n-2)^⊺ L^⊺ ]]


The block element of Σ on the ith row and jth column is given by  

    Σ_ij = LD^i-2Σ (D^⊺)^j-2 L^⊺.

By using equation (<ref>), we arrive at,

    Σ_ij = -k_BT [LD^i-1 Q (D^⊺)^j-2 L^⊺+ LD^i-2Q(D^⊺)^j-1 L^⊺].


Next we define matrix S=[[ A_22^-1       0;       0      -I ]]. It can be easily seen that since Γ=γ I, we have 
    SD^⊺ = DS,    SL^⊺=R,    Q-S=1/2γ k_BTΣ.


With these identities, we are able to manipulate terms, and get,

    Σ_ij   = -k_BT LD^i-1 S (D^⊺)^j-2L^⊺ - k_BT LD^i-2 S (D^⊺)^j-1L^⊺ -  1/γ LD^i-1Σ (D^⊺)^j-2L^⊺- 1/γLD^i-1Σ (D^⊺)^j-2L^⊺
       = -2 k_BT LD^i+j-3 R - 1/γΣ_i+1,j - 1/γΣ_i,j+1 =  - 1/γΣ_i+1,j - 1/γΣ_i,j+1 - 2 k_BTM_i+j-3,

where we have used the notation M_-1=-M_∞.

Meanwhile, the block elements of Σ in the second column and the second row are all zeros, since by direct calculation, Σ L^⊺=0.  Furthermore, we have Σ_11=2k_BTM_∞. For example when j=1, we have

    Σ_i1 = -γΣ_i-1,j-2k_BT M_i-3  , i≥ 3.

Therefore we are able to write out entries of Σ using M_is. 



§ REPRESENTATION OF M_IS

This is the proof of lemma 4.
This calculation is based on the formulas in 
(<ref>) and (<ref>) for the memory function θ(t). The Laplace transform will be given by,

    Θ = L [ λ^-1 I - D] R = 
      λ L [
      [       I    -λ I;  λ A_22 (1+λγ)I ]]^-1  R = λ L 
        [
      [ ( I + λ^2 (1+λγ)^-1A_22)^-1                           *;                                                       * ]] R.

Here we have used a block inversion formula, and the fact that the second block in both L and R are zero. 

At this point, we can invoke the Neumann series of the first diagonal block and we have,


    Θ = λ A_12A^-1_22A_21 -  λ^3 A_12  (1 -λγ + λ^2γ^2 + ⋯) A_21 + 
      -  λ^5 A_12  (1 -λγ + λ^2γ^2 + ⋯)^2 A_22A_21 +
     ⋯.


Therefore, the patterns in the representation of M_i's can be observed.
 



































As an example, the first few moments are listed below

 

    M_∞=γ A_12(A_22^-1)^2A_21, M_0=A_12A^-1_22A_21,    M_1=0, 
       M_2=-A_12A_21,   M_3=γ A_12A_21,   M_4=A_12A_22A_21-γ^2 A_12A_21,
       M_5=-2γ A_12A_22A_21+γ^3A_12A_21,  
    M_6 = -A_12A_22A_21+3γ^2 A_12A_22A_21-γ^4 A_12A_21,
       M_7 = 3γ A_12A_22^2A_21 - 4γ^3 A_12A_22A_21 + γ^5 A_12A_21.


In addition, here are a few identities that is used to prove second FDT for order n=3, 4, 5. 

    γ M_2+M_3=0,
      γ^2 M_3 + 2 γ M_4 +M_5 = 0,  γ^3 M_4 +3 γ^2 M_5 +3γ M_6 + M_7 =0.




unsrt


 