Sampling-based model order reduction for stochastic differential equations driven by fractional Brownian motion

In this paper, we study large-scale linear fractional stochastic systems representing, e.g., spatially discretized stochastic partial differential equations (SPDEs) driven by fractional Brownian motion (fBm) with Hurst parameter H > 1 / 2 . Such equations are more realistic in modeling real-world phenomena in comparison to frameworks not capturing memory effects. To the best of our knowledge, dimension reduction schemes for such settings have not been studied so far. In this work, we investigate empirical reduced order methods that are either based on snapshots (e.g. proper orthogonal decomposition) or on approximated Gramians. In each case, dominant subspaces are learned from data. Such model reduction techniques are introduced and analyzed for stochastic systems with fractional noise and later applied to spatially discretized SPDEs driven by fBm in order to reduce the computational cost arising from both the high dimension of the considered stochastic system and a large number of required Monte Carlo runs. We validate our proposed techniques with numerical experiments for some large-scale stochastic differential equations driven by fBm.


Introduction
The proper orthogonal decomposition (POD) is a data-driven strategy for the reduction of large-scale models that is based on the SVD of snapshot matrices and utilizes an offline-online splitting.During the offline stage, the available empirical data regarding the solution of the system is used to determine a dominant subspace of the full configuration space.After that, the high-dimensional system is projected onto that subspace with the aim of capturing the original dynamics, and in the online stage, such a reduced model is solved numerically at a lower cost of computation than the original model.The POD method has been successfully utilized in a wide variety of scientific and engineering problems, such as fluid dynamics [7,12,19], electric circuit analysis [16], or structural dynamics [2]; see also [4,5], for additional references.
There have been other applications in which POD has been used to investigate stochastic systems; for example, see [11], and [20].However, POD techniques for stochastic differential equations driven by fractional Brownian motion (fBm) have not been studied yet.
The fBm is an excellent candidate to simulate numerous phenomena thanks to its self-similar and long-range dependency qualities.Some examples of these phenomena include financial markets (see, for example, [8,13]), and traffic networks (see, e.g., [10,21]).It is vital to point out that if H ̸ = 1  2 the process W H is neither a semi-martingale nor a Markov process, and standard stochastic calculus approaches are no longer applicable.Therefore, in order to construct a stochastic calculus framework for fBm, it will be necessary to use methods that are distinct from the standard Ito calculus.In recent years, numerous developments of stochastic calculus and stochastic differential equations with respect to the fBm have taken place, particularly for H ∈ ( 1 2 , 1) (see, for example, [1,9,14,15]).Throughout this entire work, we focus on the Hurst parameter H ∈ ( 1 2 , 1), and we also employ Young's integral [22] to define an integral with respect to fBm.This concept is discussed in detail in Section 2.
In this research, we investigate large-scale linear fractional stochastic systems, resulting from the spatial discretization of SPDEs driven by fBm.We employ three approaches for this goal.The first two approaches are based on the POD, and the third method is based on an empirical Gramian.The difference between the first and second POD techniques is that the second method divides the system into two subsystems.Subsequently, POD applies to reduce each subsystem individually.The sum of these two reduced systems is then the reduced order candidate for the original system.To employ the third technique, a Gramian is introduced which characterizes dominant subspaces in terms of its eigenspaces.Such approaches have the advantage of finding a good quality reduced system uniformly in the control.We approximate the Gramian using samples of the solution since it cannot be determined via a Lyapunov equation as it can be done having Wiener or Lévy noise [3].This is one of the main obstacles to studying model reduction for fractional stochastic systems.
The paper is structured as follows.In Section 2, we provide a quick introduction to fBm and describe some of its key characteristics before defining Young's integral.In Section 3, we briefly discuss the setting and the general reduced system structure.In Section 4, we present the POD and empirical Gramian approaches.We end the paper in Section 5 with numerical experiments illustrating our results.
A mean zero Gaussian process W H = {W H (t), t ≥ 0} is called fractional Brownian motion (fBm) with Hurst parameter H ∈ (0, 1) if it has the covariance function Mandelbrot and Van Ness investigated this process after it was proposed by Kolmogorov, and they constructed a stochastic integral representation of it in terms of a standard Brownian motion.The climatologist Hurst created a statistical analysis of annual Nile river runoffs, and the resulting parameter is known as the Hurst index.The fBm has the following properties: • Self-similarity: The probability distributions of the processes {a −H W H (at), t ≥ 0} and {W H (t), t ≥ 0} are the same for any constant a > 0. This is a direct consequence of (1).
• Stationary increments: By definition, the increments of a fBm are normally distributed with zero mean.From (1), it further follows that , so that the distribution only depends on t − s.
• Hölder continuity: [6] allows us to calculate the following modulus of continuity for fBm trajectories: There exists a non negative random variable G ϵ,T for all ϵ > 0 and T > 0 such that, almost surely, The Hurst parameter H accounts not only for the sign of the correlation of the increments, but also for the regularity of the sample paths, i.e., for any ϵ > 0, the trajectories are Hölder continuous of order H − ϵ.
We will now consider the integral defined by Young [22] in 1936 in order to define a stochastic integral w.r.t.W H , H > 1/2.This scenario covers integrands and integrators with a certain Hölder regularity.
• As pointed out above, the paths of W H are a.s.Hölder continuous for α = H − ϵ.We can define T 0 Y (s)dW H (s) for a process Y with (H − ϵ)-Hölder continuous trajectories path-wise in the sense of Young if H > 1/2.

Setting and projection based reduced system
We consider the following stochastic system with control u satisfying ∥u∥ 2 where 2) is defined as an integral equation using Definition 2.1 to make sense of t 0 N i x(s)dW H i (s).We aim to find a space that is spanned by the columns of V ∈ R n×r such that x(t) ≈ V x(t).Inserting this into (2) yields We enforce the error e(t) to be orthogonal to some space spanned by columns of W ∈ R n×r .Multiplying (3) with where Often W ⊤ V = I is chosen.If W = V , we have a Galerkin approximation.In this paper, we focus and empirical Galerkin based methods, where W = V has orthonormal columns.We start with the POD that is based on snapshots and subsequently discuss techniques based on empirical Gramians representing a uniform reduced model approach in u.
4 Empirical reduced order methods

Proper orthogonal decomposition
The first strategy As a first strategy, we apply the POD technique to problem (2), meaning that we create a snapshots matrix and pick V from the SVD of this matrix in order to find (4) with W = V .To be more precise, we sample the solution for a fixed u and x 0 , so we get where N, N s > 0 are the number of time points and samples.We want this method to be computationally inexpensive.Therefore, N s should be small.We consider and calculate the SVD of Z: Therefore, V Z = V * is a weighted basis for the snapshots.If Σ contains all large singular values, then V contains the dominant basis vectors and we obtain a good reduced system (ROM) from that V if the data is representative.
The second strategy We split the main system (2) into subsystems as follows: and observe that y = y 1 + y 2 .In contrast to the first strategy, the POD method is applied to each subsystem individually to receive a ROM.The sum of the subsystem reduced order outputs then is an approximation for y.

Empirical Gramian method
In this subsection, we introduce a Gramian P that allows us to identify the dominant subspaces of (2).In particular, the projection matrix V shall represent the eigenspaces of P corresponding to the large eigenvalues of the Gramian.This approach has the advantage of providing good quality ROMs uniformly in the control u.However, since fBm with H > 1/2 does not have independent increments, no Lyapunov equations associated with P can be derived.Therefore, we consider an empirical version of this Gramian meaning that we discretize the integral representation of P by a Monte-Carlo method.The fundamental solution to ( 6) is defined as the matrix valued stochastic process Φ solving for t ≥ s ≥ 0. We further set Φ(t) := Φ(t, 0).Then, we can show that the state of ( 2) is given by We define a (time-limited) Gramian P := P 1 + P 2 for (2), where where P 1 and P 2 can be interpreted as Gramians for ( 5) and ( 6).Since P can not be computed directly, we approximate it by x 2 (s i , 0, X 0 , w j )x 2 (s i , 0, X 0 , w j ) ⊤ + x 2 (T, s i , B, w j )x 2 (T, s i , B, w j ) ⊤ =: P , (10) where x 2 (t, s, X 0 , •) denotes the solution to (6) with initial time s and initial state X 0 .P is called empirical Gramian.As mentioned above, the ROM (4) with W = V based on P is obtained by conducting an eigenvalue decomposition of the empirical Gramian and subsequently choosing V to be the matrix of eigenvectors associated with the large eigenvalues of P .We begin with a modified version of an example studied in [3] and consider the following controlled stochastic partial differential equation: where a, b > 0, γ ∈ R and W H is an one-dimensional fBm and m = 1.We set Due to the fact that the Dirichlet Laplacian generates a C 0 -semigroup, the H-valued solution X of ( 11) is meant to be in the mild sense.The eigenfunctions (h k ) k∈N of this Laplacian demonstrate an orthonormal basis of H which is ordered based on the corresponding eigenvalues (−λ k ) k∈N .The average temperature in the uncontrolled area is the quantity of interest, i.e., We define the associated output operator as 4 ] 2 x(ζ)dζ for x ∈ H.A numerical solution to this SPDE can be found by first considering a spatial discretization as a possible first step.We exploit a spectral Galerkin technique in this context to build an approximation X n to X by an orthogonal projection of (11) onto the subspace H n = span{h 1 , • • • , h n }.X n converges to the SPDE solution with n → ∞ and is determined by its Fourier coefficients x(t) = (⟨X n (t), h 1 ⟩ H , . . ., ⟨X n (t), h n ⟩ H ), that satisfy (2) with coefficients where we fix a = 0.2 and b = 0.1.We refer to [3] for more details on the discretization.We now fix the dimension of spatial discretization to be n = 225.For the discretization in time, the Runge-Kutta scheme in [18] is used, where the number of time steps is N = 100.The control u(t) = e −bt and the initial condition x 0 = b (⟨sin(•), h k ⟩ H ) k=1•••n were used to calculate snapshots for the POD methods.The quality of the resulting ROMs is calculated for u(t) = 2 π cos(t) and x 0 = 0, respectively.The reduction error is measured by the quantity sup t∈[0,1] E∥y(t) − ȳ(t)∥ 2 which is approximated based on N s = 10 4 samples.We apply POD with its two different strategies, see Section 4.1, and the empirical Gramian approach, see Section 4.2, to obtain the ROM.We present the error for these three approaches for two different Hurst parameters H = 0.75 and H = 0.95 in Figure 1.According to this figure, although the performance of the first and second POD strategies differs somewhat in a few of the resulting errors, these two ways produce almost identical performance, which is superior to the result of the empirical Gramian method.In addition, the quantity of errors caused by the three aforementioned approaches with H = 0.95 for the system (11) can be observed in Table 1.
Figure 2 depicts the first 50 singular values of the snapshot matrices and the eigenvalues of the empirical Gramian P .For simplicity, these values are called Hankel singular values (HSVs) of system (2) for the three aforementioned methods.The second strategy involves plotting σ i,1 + σ i,2 , where σ i,1 and σ i,2 represent the POD singular values of subsystems ( 5) and ( 6), respectively.It is understood that small truncated HSV correspond to a small reduction error in the MOR scheme.Due to the rapid decay of these values, we can conclude that even for small reduced dimensions r, small errors can be achieved.
Example 2 We consider the following controlled stochastic partial differential equation which is a modification of the example studied in [17] is investigated and the output is an approximation of the string's midpoint, which is where a, b > 0 and ϵ > 0. We know from [17] that the Galerkin system looks as follows: where the coefficients are represented by for j = 1 . . ., n and ν = 1, . . ., n 2 , • the output matrix is given by C ⊤ = (c k ) k=1,...,n with c 2l = 0 and c 2l−1 = 1 √ 2πl 2 cos l π 2 − cos l π 2 + ϵ .In this case, we assume b = 0.1, a = 2 and H = 0.75.Further, the sizes of spatial and time discretization are n = 100 and N = 100, respectively.In this example, we consider the same scenario as we did in the first example; which means we calculate snapshots and errors using different control functions and different initial values.We use N s = 100, u(t) = 2 1−e −2π e −t and the initial state defined above in the computation of the snapshots.We calculated the problem's error using the following values: N s = 10 4 , u(t) = 2 π cos(t) and x 0 = 0.As can be seen in Figure 3 and Table 2, the performance of the three techniques is nearly identical.With a more thorough perspective, however, the results of the first and second POD strategies are totally consistent, and these two ways are superior to the empirical Gramian method by a minor difference.

Fig. 1 :
Fig. 1: The error for the three approaches with two different Hurst parameters H = 0.75 (a) and H = 0.95 (b).

Fig. 3 :
Fig.3: The error for the three approaches with Hurst parameters H = 0.75.