A Scalable Observation-Driven Time-Dependent Basis for a Reduced Description of Transient Systems

TL;DR

This paper introduces a scalable, observation-driven variational method to extract a time-dependent basis for transient systems, enabling reduced-order modeling without requiring knowledge of the underlying system model.

Contribution

It presents a novel variational principle-based approach for deriving a time-dependent orthonormal basis solely from system observations, applicable to unknown systems.

Findings

Method effectively captures transient dynamics in diverse systems

Outperforms static basis methods like POD and DMD in transient scenarios

Computationally efficient with linear complexity and memory requirements

Abstract

We present a variational principle for the extraction of a time-dependent orthonormal basis from random realizations of transient systems. The optimality condition of the variational principle leads to a closed-form evolution equation for the orthonormal basis and its coefficients. The extracted modes are associated with the most transient subspace of the system, and they provide a reduced description of the transient dynamics that may be used for reduced-order modeling, filtering and prediction. The presented method is matrix-free and relies only on the observables of the system and ignores any information about the underlying system. In that sense, the presented reduction is purely observation-driven and may be applied to systems whose models are not known. The presented method has linear computational complexity and memory storage requirement with respect to the number of observables and the number of random realizations. Therefore, it may be used for a large number of observations and samples. The effectiveness of the proposed method is tested on three examples: (i) stochastic advection equation, (ii) a reduced description of transient instability of Kuramoto-Sivashinsky, and (iii) a transient vertical jet governed by incompressible Navier-Stokes equation. In these examples, we contrast the performance of the time-dependent basis versus static basis such as proper orthogonal decomposition, dynamic mode decomposition and polynomial chaos expansion.

Figures

Tables

Table 1. Table 1: Reconstruction error of Kuramoto-Sivashinsky equation: comparison between dynamic mode decomposition (DMD) and dynamic basis (DB).

r	2	4	6	8
Error (DMD)	9.01$\times {10}^{-2}$	9.027$\times {10}^{-2}$	9.18$\times {10}^{-2}$	9.24$\times {10}^{-2}$
Error (DB)	2.10$\times {10}^{-1}$	6.92$\times {10}^{-3}$	$1.88\times {10}^{-3}$	8.92$\times {10}^{-4}$

Equations

$$\mathbf{x}=\mathbf{Uz}+ϵ$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t;\mathbf{\xi }),$$

$$\mathbf{T}(t)$$

$$\mathbb{E}[\mathbf{T}]=\mathbf{T}{\mathbf{w}}_{\mathbf{\xi }},$$

$$\displaystyle\big{<}\mathbf{u}_{i},\mathbf{u}_{j}\big{>}$$

$$\mathbb{E}[{\mathbf{y}}_i^T{\mathbf{y}}_j]$$

$$\mathbf{T}(t)=\mathbf{U}(t)\mathbf{Y}(t{)}^T+\mathcal{E}(t),$$

$$\mathbf{Y}(t)$$

$$\mathbf{U}(t)$$

$$\left\|\mathbb{E}\big{[}\mathbf{T}\big{]}^{2}\right\|_{F}^{1/2}=\bigg{(}\sum_{i=1}^{n}\sum_{j=1}^{s}\mathbf{w}_{{\mathbf{x}}_{i}}\mathbf{w}_{{\boldsymbol{\xi}}_{j}}\mathbf{T}_{i,j}^{2}\bigg{)}^{1/2},$$

$$\mathcal{F}(\dot{\mathbf{U}},\dot{\mathbf{Y}})=\left\|\mathbb{E}\bigg{[}\frac{d(\mathbf{U}\mathbf{Y}^{T})}{dt}-\dot{\mathbf{T}}\bigg{]}^{2}\right\|_{F},$$

$$\big{<}\dot{\mathbf{u}}_{i},\mathbf{u}_{j}\big{>}+\big{<}\mathbf{u}_{i},\dot{\mathbf{u}}_{j}\big{>}=0,$$

$$\mathcal{G}(\dot{\mathbf{U}},\dot{\mathbf{Y}})=\left\|\mathbb{E}\bigg{[}\sum_{i=1}^{r}\frac{d(\mathbf{u}_{i}\mathbf{y}^{T}_{i})}{dt}-\dot{\mathbf{T}}\bigg{]}^{2}\right\|_{F}+\sum_{i,j=1}^{r}\lambda_{ij}(t)\big{(}\big{<}\mathbf{u}_{i},\dot{\mathbf{u}}_{j}\big{>}-\boldsymbol{\phi}_{ij}\big{)},$$

$$\dot{\mathbf{U}}$$

$$\dot{\mathbf{Y}}$$

$$\dot{\mathbf{U}}$$

$$\dot{\mathbf{Y}}$$

$$\prod _{\perp \mathbf{U}}\mathbb{E}[\dot{\mathbf{T}}\mathbf{Y}]{\mathbf{C}}^{-1}=\left(\mathbf{I}-\mathbf{U}{\mathbf{U}}^T{\mathbf{W}}_{\mathbf{x}}\right)\dot{\mathbf{T}}{\mathbf{W}}_{\mathbf{\xi }}\mathbf{Y}{\mathbf{C}}^{-1}=\left(\dot{\mathbf{T}}{\mathbf{W}}_{\mathbf{\xi }}\mathbf{Y}-\mathbf{U}\dot{{\mathbf{T}}_r}\right){\mathbf{C}}^{-1},$$

$$\dot{\mathbf{U}}$$

$$\dot{\mathbf{Y}}$$

$$\dot{\mathbf{U}}=\prod _{\perp \mathbf{U}}\mathbb{E}[\dot{\mathbf{T}}\mathbf{Y}]{\mathbf{C}}^{-1}=0,$$

$${\mathbf{T}}_0=\sum _{i=1}^r\sqrt{{\lambda }_i}{\mathbf{v}}_i{\mathbf{z}}_i^T+{\mathbf{E}}_0,$$

$${\dot{\mathbf{T}}}_k=\frac{-{\mathbf{T}}_{k+2}+8{\mathbf{T}}_{k+1}-8{\mathbf{T}}_{k-1}+{\mathbf{T}}_{k-2}}{12\Delta t}+\mathcal{O}(\Delta t^4).$$

$$\mathbf{C}\mathbf{\psi }=\mathbf{\psi }\mathbf{\Lambda },$$

$$\hat{\mathbf{Y}}=\mathbf{Y}\mathbf{\psi },\hat{\mathbf{U}}=\mathbf{U}\mathbf{\psi }.$$

$$\frac{\partial u}{\partial t}+V(\xi )\frac{\partial u}{\partial x}=0,x\in [-1,1],t\in [0,T_f],$$

$$V(\xi )=\overline{v}+\sigma \xi ,u_0(x;\xi )=\sin \pi x,$$

$$u(x,t;\xi)=\sin\pi\big{(}x-(\overline{v}+\sigma\xi)t\big{)},$$

$$\mathbb{E}[u(x,t;\xi )]=\frac{1}{\sigma \pi t}\sin \pi (x-\overline{v}t)\sin \sigma \pi t,$$

$$\mathbb{E}[u(x,t;\xi {)}^2]-\mathbb{E}[u(x,t;\xi ){]}^2=\frac{1}{2}-\frac{1}{4\sigma \pi t}\cos 2\pi (x-\overline{v}t)\sin 2\sigma \pi t-\frac{1}{{\sigma }^2{\pi }^2{t}^2}{\sin }^2\pi (x-\overline{v}t){\sin }^2\sigma \pi t,$$

$$\Sigma(t)=\int_{-1}^{1}\big{(}\mathbb{E}[u(x,t;\xi)^{2}]-\mathbb{E}[u(x,t;\xi)]^{2}\big{)}dx=1-\frac{\sin^{2}\sigma\pi t}{(\sigma\pi t)^{2}}.$$

$$u(x,t;\xi )-\mathbb{E}[u(x,t;\xi )]=u_1{y}_1+u_2{y}_2,$$

$$u_1(x,t)$$

Full text

1 Introduction

The ability to discover intrinsic structure from data has been called the the fourth paradigm of scientific discovery [1]. Dimension reduction enables such discoveries and clearly belongs to the fourth paradigm. The goal of dimension reduction is to map high-dimensional data into a meaningful representation of reduced dimensionality. As a result, dimension reduction has been widely used in classification, visualization, and compression of high-dimensional data [2]. In the scientific and engineering applications dimension reduction has been used to build reduced order models to propagate uncertainty, control and optimization [3, 4].

This paper focuses on dimension reduction of vector-valued observables $\mathbf{x}(t;\omega)\in\mathbb{R}^{n}$ generated by time-dependent stochastic system where $\omega$ denotes a random observation. A linear reduction model seeks to approximate the target $\mathbf{x}$ with

[TABLE]

where $\mathbf{U}\in\mathbb{R}^{n\times r}$ is a deterministic matrix whose columns can be regarded as modes and $\mathbf{z}\in\mathbb{R}^{r}$ are the reduced stochastic variable, where $r$ is the reduction size and $\epsilon$ is the reduction error. Both $\mathbf{U}$ and $\mathbf{z}$ are in general time-dependent. In the view of Equation 1, polynomial chaos expansion (PCE) can be seen as a special case where $\mathbf{z}$ represented by polynomial chaos – a static (i.e. time-independent) basis – while $\mathbf{U}=\mathbf{U}(t)$ is time-dependent. PCE has reached widespread popularity especially for steady-state and parabolic stochastic partial differential equations (SPDE) [5]. However, for time-dependent dynamical systems with transient instabilities the number of PCE terms required for a given level of accuracy must rapidly increase [6]. Similarly, it was shown in [7] that the reduction size of PCE must increase in time to maintain a given accuracy. On the other hand, data-driven methods such as proper orthogonal decomposition (POD) and dynamic mode decomposition (DMD) belong to a class of methods where the modes $\mathbf{U}$ are static and they are computed from numerical/experimental observations of the dynamical system and $\mathbf{z}(t;\omega)$ is obtained by often projection of the data to these modes [8]. For highly transient systems, DMD and POD require large number of modes to capture transient time behavior [9, 10]. This is partly because both DMD and POD decompositions seek to obtain an optimal reduction of the data in a time-average sense. In the case of DMD a best-fit linear dynamical system matrix $\mathbf{A}$ is sought so that $\|\mathbf{x}^{k+1}-\mathbf{A}\mathbf{x}^{k}\|_{2}$ is minimized across all snapshots, where $k$ shows the time snapshots. Similarly, the POD is an optimal basis for all snapshots. However, transient instabilities, for example, by definition have a short life time. The effect of these transient instabilities may be lost in a reduction based on time average.

More recently dynamically orthogonal (DO) [11] and bi-orthogonal (BO) [12] expansions have proven promising in tackling time-dependent SPDEs. In these reduction techniques both $\mathbf{U}(t)$ and $\mathbf{z}(t;\omega)$ are dynamic and they evolve with time. In DO/BO framework, exact evolution equations for the evolution of $\mathbf{U}(t)$ and $\mathbf{z}(t;\omega)$ are derived. In the case of SPDEs, the evolution equation for $\mathbf{U}(t)$ is a system of deterministic PDEs that are solved fully coupled with a system of stochastic ordinary differential equations (SODE) and an evolution equation for the mean equation. Although DO and BO algorithms lead to different evolution equations for the stochastic coefficient and the spatial basis, it is shown in [13, 14] that both these methods are equivalent, i.e. the BO and DO reductions span the same subspaces. In deterministic applications, time-dependent basis have been used for dimension reduction of deterministic Schrodinger equations by the Multi Configuration Time Dependent Hartree (MCTDH) method [15, 16], and later by Koch and Lubich [17]. Recently, optimally time-dependent reduction (OTD) was introduced [18], in which an evolution equation for a set of orthonormal modes is found that capture the most instantaneously transient subspace of a dynamical system [19].

The motivation for this paper is to obtain a reduced description of transient dynamics that relies only on random observations — i.e. an equation-free reduction. This is in contrast to model-dependent reduction techniques, for example DO and BO methodologies, that require solving high-dimensional ODE/PDE for the determination of the time-dependent modes. The model-dependence of these techniques prevents the application of these methods to complex systems whose underlying models are not known or suffer from physics deficiency due to multi-scale physics effects, for example, atmospheric sciences and weather forecasting, computational neuroscience and complex biological systems. Moreover, the evolution equation for the DO/BO modes is tied to the equation that governs the dynamical system at hand, which would require a case by case implementation.

The objective of this paper is to find a reduced description of transient dynamics that is observation driven and discards any information about the underlying model. To this end, we present a variational principle for the determination of the set of time-dependent orthonormal modes that capture the low-dimensional structure from observable of the stochastic system. The structure of the paper is as follows: In Section 2, we present the variational principle for the extraction of the time-dependent modes and their coefficients. In Section 3, we present three demonstration cases: (a) the stochastic advection equation, (b) Kuramoto-Sivashinsky equation, and (c) a transient vertical jet in quiescent background flow governed by incompressible Navier-Stokes equations. In Section 4, we present the conclusions of our work.

2 Methodology

2.1 Setup

We consider a stochastic dynamical system in the form of:

[TABLE]

where $\mathbf{x}$ is the state space and the above dynamical system is finite-dimensional, i.e. $\mathbf{x}\in\mathbb{R}^{m}$, or infinite dimensional, expressed, for example, by SPDEs and $\boldsymbol{\xi}=(\xi_{1},\xi_{2},\dots,\xi_{d})\in\mathbb{R}^{d}$ is a $d$-dimensional random space with the joint probability given by $\rho(\boldsymbol{\xi})$. The dynamical system may not be known and the proposed reduction relies only on vector-valued observations of the above system. We base our formulation on continuous-time observations, represented by matrix $\mathbf{T}(t)\in\mathbb{R}^{n\times s}$, whose columns contain vector-valued observations of the stochastic system at time $t$:

[TABLE]

Each column of this matrix is obtained by first drawing a random sample of $\boldsymbol{\xi}$ with joint probability distribution of $\rho(\boldsymbol{\xi})$ and solving the forward model given by equation 2. Before we proceed we introduce the sample mean operator:

[TABLE]

where $\mathbf{w}_{\boldsymbol{\xi}}=[\mathbf{w}_{{\boldsymbol{\xi}}_{1}},\mathbf{w}_{{\boldsymbol{\xi}}_{2}},\dots,\mathbf{w}_{{\boldsymbol{\xi}}_{s}}]^{T}\in\mathbb{R}^{s\times 1}$ represents the sample weights. For example, for standard Monte Carlo samples, $\mathbf{w}_{\boldsymbol{\xi}}=\mathbf{1}/s$, where $\mathbf{1}\in\mathbb{R}^{s\times 1}$ is a column vector with all elements equal to 1. For cases where controlled observations of system 2 is possible, the probabilistic collocation strategy may be used for sampling. In that case $\mathbf{w}_{\boldsymbol{\xi}}$ represents the associated probabilistic collocation weights. Similarly, $\mathbf{w}_{\mathbf{x}}=[\mathbf{w}_{{\mathbf{x}}_{1}},\mathbf{w}_{{\mathbf{x}}_{2}},\dots,\mathbf{w}_{{\mathbf{x}}_{n}}]^{T}\in\mathbb{R}^{n\times 1}$ represents the associated weights in the state space. For applications governed by SPDEs, $\mathbf{w}_{\mathbf{x}}$ may be chosen as space integration weights (e.g. quadrature weights), such that $\mathbf{w}_{\mathbf{x}}^{T}\mathbf{u}$ approximates $\int_{\Omega}ud\Omega$, where $\Omega$ is the spatial domain and $\mathbf{u}$ contains vector-valued observations of spatial function $u$. We assume that $\mathbf{T}(t)$ has zero mean, i.e. $\mathbb{E}[\mathbf{T}(t)]=0$. This is accomplished by subtracting the ensemble mean from all the realizations. Let’s define the inner product in both state and random spaces:

[TABLE]

where $\mathbf{W}_{\mathbf{x}}=\mbox{diag}(\mathbf{w}_{\mathbf{x}})$ and $\mathbf{W}_{\boldsymbol{\xi}}=\mbox{diag}(\mathbf{w}_{\boldsymbol{\xi}})$ are matrices of size $n\times n$ and $s\times s$, respectively. Therefore, \big{<}\mathbf{u}_{i},\mathbf{u}_{j}\big{>} is an inner product of $\mathbf{u}_{i}$ and $\mathbf{u}_{j}$ in the state space, and $\mathbb{E}[\mathbf{y}^{T}_{i}\mathbf{y}_{j}]$ is an inner product of $\mathbf{y}_{i}$ and $\mathbf{y}_{j}$ in the random space. Our goal is to build a low-rank approximation that can capture transient behavior in the stochastic data $\mathbf{T}(t)$. To this end, we seek a reduction in the form of:

[TABLE]

with

[TABLE]

where $\mathbf{Y}(t)\in\mathbb{R}^{s\times r}$ and \mathbf{Y}(t)=\big{[}\mathbf{y}_{1}(t)\ |\ \mathbf{y}_{2}(t)\ |\ \dots\ |\ \mathbf{y}_{r}(t)\big{]} and $\mathbf{U}(t)\in\mathbb{R}^{n\times r}$ is a set of orthonormal modes, i.e. \big{<}\mathbf{u}_{i}(t),\mathbf{u}_{j}(t)\big{>}=\delta_{ij} and $\mathcal{E}$ is the reduction error. In the proposed form $\mathbf{U}(t)$ represents a time-dependent low rank subspace which we will refer to as dynamic basis and $\mathbf{Y}(t)$ are the stochastic coefficients of the dynamic basis and we will refer to $\mathbf{Y}(t)$ as the stochastic coefficients. In the following we present a variational principle that results in closed-form evolution equations for the dynamic basis $\mathbf{U}(t)$ and the stochastic coefficients $\mathbf{Y}(t)$.

2.2 Variational Principle

In this section we present a variational principle from which the evolution equations of the low-rank approximation are derived. Our goal is to build a low-rank approximation that can capture transient data. To this end, we first introduce a weighted Frobenius norm:

[TABLE]

where $\mathbf{T}_{i,j}$ is the entry of matrix $\mathbf{T}$ at the $i$th row and the $j$th column. Now we consider the variational principle in the form of:

[TABLE]

subject to the orthonormality condition of the modes, i.e. \big{<}\mathbf{u}_{i}(t),\mathbf{u}_{j}(t)\big{>}=\delta_{ij}. The variational principle 5 seeks to minimize the difference between the rate of change of all observations and the dynamic basis reduction, and the control parameters are $\dot{\mathbf{U}}$ and $\dot{\mathbf{Y}}$. To incorporate the orthonormality condition , we observe that:

[TABLE]

and we denote \boldsymbol{\phi}_{ij}(t):=\big{<}\mathbf{u}_{i},\dot{\mathbf{u}}_{j}\big{>}. It is clear from the above equation that $\boldsymbol{\phi}_{ij}(t)$ is a skew-symmetric matrix, i.e. $\boldsymbol{\phi}_{ij}(t)=-\boldsymbol{\phi}_{ji}(t)$. We incorporate the orthonormality constraint into the variational principle via Lagrange multipliers. This follows:

[TABLE]

where $\lambda_{ij}(t),i,j=1,\dots,r$ are the Lagrange multipliers. In Appendix A, we show that the first-order optimality condition of the above functional leads to the closed-form evolution equation of $\mathbf{U}$ and $\mathbf{Y}$ as given by:

[TABLE]

where $\mathbf{C}=\mathbb{E}[\mathbf{Y}^{T}\mathbf{Y}]$ is the reduced covariance matrix ($\mathbf{C}\in\mathbb{R}^{r\times r}$), and $\underset{\perp\mathbf{U}}{\prod}$ is the orthogonal projection matrix to the space spanned by the complement of $\mathbf{U}$, i.e. $\underset{\perp\mathbf{U}}{\prod}\mathbf{V}=\mathbf{V}-\mathbf{U}\mathbf{U}^{T}\mathbf{W}_{\mathbf{x}}\mathbf{V}=(\mathbf{I}-\mathbf{U}\mathbf{U}^{T}\mathbf{W}_{\mathbf{x}})\mathbf{V}$.

Solving equations 8 and 9 require determination of $\boldsymbol{\phi}$. We show in Appendix B that for any skew-symmetric matrix $\boldsymbol{\phi}$ the above reduction remains equivalent for all times. Equivalent reductions can be defined as follows: The two reduced subspaces $\{\mathbf{Y},\mathbf{U}\}$ and $\{\widetilde{\mathbf{Y}},\widetilde{\mathbf{U}}\}$, where $\mathbf{U},\widetilde{\mathbf{U}}\in\mathbb{R}^{n\times r}$ and $\mathbf{Y},\widetilde{\mathbf{Y}}\in\mathbb{R}^{s\times r}$, are equivalent if there exists a transformation matrix $\mathbf{R}\in\mathbb{R}^{r\times r}$ such that $\mathbf{U}=\widetilde{\mathbf{U}}\mathbf{R}$ and $\mathbf{Y}=\widetilde{\mathbf{Y}}\mathbf{R}$, where $\mathbf{R}$ is an orthogonal rotation matrix, i.e. $\mathbf{R}^{T}\mathbf{R}=\mathbf{I}$. Therefore, changing $\boldsymbol{\phi}$ does not affect the subspace the both $\mathbf{U}$ and $\mathbf{Y}$ span, and it will only lead to an in-subpsace rotation. In the rest of this paper, we take the simplest case of $\boldsymbol{\phi}\equiv\mathbf{0}$. In that case the evolution of the dynamic basis reduces to:

[TABLE]

For the computation of the right hand side of equation 10, $\mathbf{U}\mathbf{U}^{T}$ may not be computed nor stored explicitly. This can be done by observing that:

[TABLE]

where $\dot{\mathbf{T}_{r}}=\mathbf{U}^{T}\mathbf{W}_{\mathbf{x}}\dot{\mathbf{T}}\mathbf{W}_{\boldsymbol{\xi}}\mathbf{Y}$ and $\dot{\mathbf{T}_{r}}\in\mathbb{R}^{r\times r}$. Therefore, the evolution equations for the dynamic basis and the reduced manifold are given by:

[TABLE]

Besides the simplicity of the above evolution equations, the choice of $\boldsymbol{\phi}=\mathbf{0}$ also leads to an insightful interpretation of the dynamic basis reduction, namely the update of the dynamic basis $\dot{\mathbf{U}}$ is always orthogonal to the space spanned by $\mathbf{U}$, i.e. $\dot{\mathbf{U}}\perp\mathbf{U}$, and this update is driven by the projection of $\mathbb{E}[\dot{\mathbf{T}}\mathbf{Y}]$ onto to the orthogonal complement of $\mathbf{U}$ (see equation 12). This means that if $\mathbb{E}[\dot{\mathbf{T}}\mathbf{Y}]$ is in the span of $\mathbf{U}$, then

[TABLE]

However, when $\mathbb{E}[\dot{\mathbf{T}}\mathbf{Y}]$ is not in the span of $\mathbf{U}$, it results in $\dot{\mathbf{U}}\neq\mathbf{0}$. This built-in adpativity of the dynamic basis allows $\mathbf{U}$ to “follow” the data instantaneously. The $\mathbf{Y}$ coefficients are updated by projection of $\dot{\mathbf{T}}$ onto the space spanned by $\mathbf{U}$ (see equation 13). Therefore, the evolution equation of $\mathbf{Y}$ captures the evolution of $\mathbf{T}$ within the subspace $\mathbf{U}$. We also observe that the above scheme is agnostic with regard to the underlying model that generates the observation and observations may come from a system whose model is not known.

2.3 Workflow and computational complexity

The dynamic basis reduction is governed by time-continuous evolution equations 12 and 13. For time-discrete data these equations must be discretized in time. Here we summarize all the steps required for solving equations 12 and 13:

Step 1: Given the time-discrete realizations in the form of $\widetilde{\mathbf{T}}_{k},k=0,1,2,\dots,K$ we compute zero-mean by subtracting the instantaneous sample mean from the data: $\mathbf{T}_{k}=\widetilde{\mathbf{T}}_{k}-\mathbb{E}[\widetilde{\mathbf{T}}_{k}]\mathbf{1}_{1\times s}$, where $\mathbf{1}_{1\times s}$ is a row vector of size $s$ whose all elements are equal to one.

Step 2: We compute the initial condition for $\mathbf{Y}_{0}$ and $\mathbf{U}_{0}$ by performing Karhunen-Loève decomposition of the $\mathbf{T}_{0}$ at $t=0$. This follows:

[TABLE]

where vectors $\mathbf{v}_{i}\in\mathbb{R}^{n}$ are an orthonormal basis, $\mathbf{z}_{i}\in\mathbb{R}^{s}$ are the stochastic coefficients, and $\lambda_{i}$ are their correspond eigenvalues and $\mathbf{E}_{0}$ is the reduction error. The initial conditions for $\mathbf{U}(t)$ and the stochastic coefficients $\mathbf{Y}(t)$ are then given by: $\mathbf{u}_{i_{0}}=\mathbf{v}_{i}$ and $\mathbf{y}_{i_{0}}=\sqrt{\lambda_{i}}\mathbf{z}_{i}$ for $i=1,2,\dots,r$.

Step 3: The time derivative of the observations $\dot{\mathbf{T}}(t)$ can be estimated via finite difference schemes. For example the fourth-order central finite difference approximation takes the form of:

[TABLE]

For the first two time steps one sided finite-difference schemes may be used.

Step 4: Evolve the dynamic basis and the stochastic coefficients to the next time step $k$ using a time-integration scheme. We use fourth-order Runge-Kutta scheme. A comparison with lower-order time integration is scheme is given in the last demonstration case.

Step 5: Perform a Gram–Schmidt orthonormalization on the dynamic basis $\mathbf{U}$. Although the dynamic basis remains orthonormal by construction, the numerical orthonormality of $\mathbf{U}$ may not be preserved due to the nonlinearity of the evolution equations 12 and 13. An enforcement of orthonormality of $\mathbf{U}$ resolves this issue. Now steps 3-5 can be repeated in a time-loop till the completion of all time steps.

A quick inspection of the evolution equations 12-13 reveals two favorable properties with regard to the computational cost and memory storage requirement of the dynamic basis reduction: (1) The computational complexity of solving the evolution equations 12-13 is $\mathcal{O}(rns)$, i.e. linear computational complexity with respect to the reduction size $r$, the size of the observations in state space $n$ and the sample size $s$. We note that given that $r$ is often a small number, i.e. $r<<n$ and $r<<s$, the cost of inverting the covariance matrix $\mathbf{C}$ is often negligible. (2) The memory storage requirement is also linear with respect to $n$, $r$ and $s$. This is because at each time instant only computation of $\dot{\mathbf{T}}$ is required, which can be computed from few snapshots of $\mathbf{T}$. This allows for loading matrix $\mathbf{T}$ in small time batches. The linear computational and memory complexity of the dynamic basis reduction may be exploited for applications with very large data sets.

2.4 Mode ranking within the subspace

In this section we present a ranking criterion within the dynamic subspace $\mathbf{U}$. This can be accomplished by using the reduced covaraince matrix $\mathbf{C}$, which, in general, is a full matrix and this implies that the $\mathbf{Y}$ coefficients are in general correlated. The uncorrelated $\mathbf{Y}$ coefficients may be computed via the eigen-decomposition of the covariance matrix $\mathbf{C}$. To this end, let

[TABLE]

be the eigenvalue decomposition of matrix $\mathbf{C},$ where $\boldsymbol{\psi}\in\mathbb{R}^{r\times r}$ is matrix of eigenvectors of $\bf{C}$ and $\boldsymbol{\Lambda}=\mbox{diag}(\lambda_{1},\lambda_{2},\dots,\lambda_{r})$ is a diagonal matrix of the eigenvalues of matrix $\mathbf{C}$. The uncorrelated $\mathbf{Y}$ coefficients and the dynamic basis associated with those uncorrelated coefficients, denoted by $\hat{\mathbf{Y}}$ and $\hat{\mathbf{U}}$, are obtained by an in-subspace rotation:

[TABLE]

It is easy to verify that $\hat{\mathbf{Y}}$ coefficients are uncorrelated, i.e. $\mathbb{E}[\hat{\mathbf{Y}}\hat{\mathbf{Y}}^{T}]=\boldsymbol{\Lambda}$. We adopt a variance-based ranking for $\{\hat{\mathbf{Y}},\hat{\mathbf{U}}\}$. Since matrix $\mathbf{C}$ is positive semidefinite, without loss of generality, let $\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{r}\geq 0$, where $\lambda_{i}$ is the variance of $\hat{\mathbf{y}}_{i}$.

Since $\{\hat{\mathbf{Y}},\hat{\mathbf{U}}\}$ are obtained by an in-subspace rotation of $\{\mathbf{Y},\mathbf{U}\}$, they are equivalent representations of the low rank dynamics, i.e. $\hat{\mathbf{U}}\hat{\mathbf{Y}}^{T}=\mathbf{U}\boldsymbol{\psi}\boldsymbol{\psi}^{T}\mathbf{Y}^{T}=\mathbf{U}\mathbf{Y}^{T}$, where we have used $\boldsymbol{\psi}\boldsymbol{\psi}^{T}=\mathbf{I}$, since eigenvectors of a symmetric matrix ($\mathbf{C}$) are orthonormal. In all the demonstration cases in this paper, we use the ranked uncorrelated representation of the dynamic basis for the analysis of the results. For the sake of simplicity in the notation, we use $\{\mathbf{Y},\mathbf{U}\}$ to refer to the ranked uncorrelated representation instead of $\{\hat{\mathbf{Y}},\hat{\mathbf{U}}\}$, noting that both of these representations are equivalent.

3 Demonstration cases

3.1 Stochastic advection equation

As the first example, we apply dynamic basis reduction to a one-dimensional stochastic advection equation. As we will show in this section, the stochastic advection equation chosen here is exactly two-dimensional, when expressed versus the dynamic basis. It also has an exact analytical solution for the dynamic basis $\mathbf{U}$ and the stochastic coefficients $\mathbf{Y}$, which enables an appropriate verification of the dynamic basis algorithm. The stochastic advection equation is given by:

[TABLE]

subject to periodic boundary conditions at the ends of the spatial domain $x\pm 1$ and the initial condition of $u_{0}(x;\xi)$. In the above equation $T_{f}$ is the final time. The solution to the above equation is a stochastic random field in the form of $u(x,t;\xi)$. The stochastic advection speed $V(\xi)$ and the deterministic initial condition are chosen to be

[TABLE]

where $\xi$ is a one-dimensional uniform random variable defined on $[-1,1]$ and $\sigma$ is a constant and $\overline{v}$ is the mean of the transport velocity. The exact solution of equation 15 is

[TABLE]

and its mean and variance are:

[TABLE]

The expectation operator is given by $\mathbb{E}[u(x,t;\xi)]=\displaystyle\int_{-1}^{1}u(x,t;\xi)\rho(\xi)d\xi$, where $\rho(\xi)=1/2$ is the probability density function. The total variance is equal to:

[TABLE]

The mean-subtracted stochastic field can be expressed with the following two dimensional expansion:

[TABLE]

where

[TABLE]

It is easy to verify that the two-dimensional expansion of the above random field satisfies the dynamic basis orthonormality condition in the continuous form, i.e.

[TABLE]

The expansion of the stochastic field as expressed above in Equation 16 is in bi-orthogonal form since $\mathbb{E}[y_{i}y_{j}]=0,i\neq j$.

The details of solving the evolution equations of the dynamic basis $\mathbf{U}(t)$ and the stochastic coefficients $\mathbf{Y}(t)$ are as follows. The matrix-valued random observations $\mathbf{T}(t)$ are generated by solving the stochastic advection equation 15 at $s=64$ collocation points. Theses sampling points are chosen to be the roots of the Legendre polynomial of order 64 defined on the interval of $\xi\in[-1,1]$. The 64 realizations of equation 15 are generated using a deterministic solver that uses Fourier modes of order $n=128$ for space discretization and fourth-order Runge-Kutta with $\Delta t=0.001$ for time-advancement. The final time is set to $T_{f}=10$ for all realizations.

We first performed dimension reduction of size $r=2$. In Figure 1, the total variance $\Sigma(t)$ obtained from the reduction is compared against the exact variance. The total variance of the low-rank approximation is obtained by $\Sigma(t)=\sum_{i=1}^{r}\lambda_{i}(t)$ where $\lambda_{i}(t)$ are the eigenvalues of the covariance matrix $\mathbf{C}$. We also solved the stochastic advection equation 15 using probabilistic collocation method (PCM), in which the stochastic solution is approximated as:

[TABLE]

where $\phi_{i}(\xi)$ are Legendre-chaos polynomials of order $i$, $v_{i}(x,t)$ are space-time coefficients and $e(x,t;\xi)$ is the expansion error. In Figure 1, the total variance obtained from different reduction sizes $r$ of the PCM are also shown where it can be seen that as the size of PCM expansion increases the PCM variance follows the exact solution for a longer time, but eventually the PCM error increases and even a 20-dimensional PCM reduction under-performs the dynamic basis reduction with $r=2$. For the long-term error analysis for this problem using the polynomial chaos expansion, the reader may refer to reference [7]. We note that in the PCM reduction the stochastic basis is fixed * a priori*, i.e. the polynomial chaos $\phi_{i}(\xi)$ while in the dynamic basis reduction the stochastic coefficients $\mathbf{y}_{i}(\xi,t)$ are time-dependent and intuitively this enables the stochastic coefficient adapt instantaneously to the data according to the variational principle 7.

In Figure 1 the covariance eigenvalues of the dynamic basis reduction with $r=3$ are compared against the exact solution. Since the solution is two-dimensional when expressed versus the dynamic basis, the third mode must remain inactive. As shown in Figure 1, this is confirmed numerically as the third eigenvalue does not pick up energy, while the first two eigenvalues are in excellent agreement with the exact solution. In general small eigenvalues of the covariance matrix may cause numerical instability due to the computation of $\mathbf{C}^{-1}$ in equation 12. The issue of singularity of matrix $\mathbf{C}$ can be resolved by computing the pseudo-inverse of the covaraince matrix. See reference [14] for similar treatment of singularity of the covaraince matrix. In Figure 2, the dynamic basis and the stochastic coefficients are compared against the analytical values for two different times. In both cases excellent agreement is observed.

3.2 Transient Instability: Kuramoto-Sivashinsky

One of the main utilities of a low-rank approximation of high-dimensional dynamical systems is to provide interpretable explanation of complex spatio-temporal dynamics. In this section, we demonstrate the application of the dynamic basis reduction to the one-dimensional stochastic Kuramoto-Sivashinsky equation. This equation has been widely studied for chaotic spatio-temporal dynamics in wide-ranging applications including reaction-diffusion systems [20], fluctuations in fluid films on inclines [21] and instabilities in laminar flame fronts [22]. The Kuramoto-Sivashinsky equation exhibits intrinsic instabilities and as a result the stochastic perturbations may undergo a significant growth. The spontaneity of the instability makes it a challenging problem for any reduction technique as the reduced subspace during the transition could be very different from that of the equilibrium states. We apply the dynamic basis reduction to realizations of the Kuramoto-Sivashinsky system to demonstrate how the $\mathbf{U}(t)$ “follows” the dynamics as the observations undergo a transient instability. This transient instability is visualized in Figure 3(a) where a sudden transition occurs near $t=0.75$. The solution shown in Figure 3(a) is one realization of the Kuramoto-Sivashinsky equation with the following setup. We consider Kuramoto-Sivashinsky equation

[TABLE]

with $\epsilon=0.01$ and subject to random initial condition

[TABLE]

In the above equation $u_{b}(x)$ is chosen to be the terminal solution of equation 17 after long-time integration. The random perturbation of the initial condition, denoted by $u^{\prime}(x;\xi)$, is set to have the periodic covariance given by:

[TABLE]

where $\sigma$ is the standard deviation and $l_{c}$ is the correlation length. We perform Karhunen-Loève decomposition

[TABLE]

where $\lambda_{i}$ and $\phi_{i}(x)$ are the $i$th eigenvalue and eigenfunction of the covariance kernel. The random perturbation is then approximated as:

[TABLE]

where $\xi_{i}\in[0,1]$ are independent uniform random numbers with unit variance. The dimension $d$ is chosen such that 99% of the energy is captured by the Karhunen-Loève decomposition. Decreasing $l_{c}$ increases the number of dimensions required to capture the energy threshold. We consider two cases of low-dimensional and high-dimensional random perturbations with $l_{c}=2.5$ and $l_{c}=0.05$, respectively. The covariance kernels for these two correlation lengths are shown in Figures 3(b)-3(c).

Three-dimensional random perturbation: We consider the correlation length of $l_{c}=2.5$ and $\sigma=0.1$. The three-dimensional ($d=3$) expansion of the KL decomposition captures more than $99\%$ of the energy of the covariance kernel. We sample the Kuramoto-Sivashinsky equation at probabilistic $Q=5$ collocation points in each $\xi_{i}$ direction, and therefore, the total number of samples is $s=5^{3}$. Correspondingly, the sample weights $\mathbf{w_{\xi}}$ are the probabilistic collocation weights. We use chebfun [23] as a blackbox solver to solve the Kuramoto-Sivashinsky equation 17 for $s=5^{3}$ realizations of the random initial condition. The solution is post-processed with fixed Chebyshev polynomial order $n=256$. Matrix $\mathbf{T}$ is sampled at discrete times with constant time interval of $dt=0.001$ from $t=0$ to $T_{f}=1.2$. This time interval encompasses the transition period that occurs $0.6<t<1.0$.

We perform dynamic basis reduction with size $r=4$. In Figure 4(a), the eigenvalues of the reduced covariance matrix versus time are shown. The eigenvalues of the covariance matrix represent variance or the energy associated with the dynamic basis. The first three dynamic basis associated with three largest eigenvalues of the covaraince matrix are shown in Figures 4(b)-4(d). We first observe that the energy of the first mode remains roughly constant for $0<t<0.5$. However, the second and third modes undergo four order of magnitude growth during the same period. Correspondingly, the first dynamic basis remain nearly invariant for $0<t<0.5$ (see Figure 4(b)). However, $\mathbf{u}_{2}$ and $\mathbf{u}_{3}$ are evolving during the same time interval (see Figures 4(c) and 4(d)). The second and third modes have near eigenvalue crossings periodically in this interval which causes very rapid change of the modes when their eigenvalues become close. The significant growth of the energy of the second and third mode reveals the transient unstable directions of the dynamical system. However, these two modes are “buried” underneath the most energetic mode, e.g. the first mode. We note that detecting the unstable directions is crucial for control of these systems [24]. Detecting these directions can be utilized in building precursors of sudden transitions or rare events in transient systems.

We also computed DMD modes to compare against the dynamic basis reduction. DMD has been very successful as a diagnostic tool to characterize complex dynamical systems. However, it is difficult for DMD to handle transient phenomena. To this end, we computed the DMD modes for a deterministic realization of Kuramoto-Sivashinsky equation. The realization is generated by solving the Kuramoto-Sivashinsky equation for the collocation point $\boldsymbol{\xi}_{0}=(-0.91,-0.54,0.91)$. This collocation point is chosen arbitrarily and any other collocation point would result in the same qualitative observations that we make in this section. To compute the DMD modes we first performed singular value decomposition (SVD) of the mean-subtracted snapshots of the Kuramoto-Sivashinsky equation in the span of $t\in[0,T_{f}]$ for $\boldsymbol{\xi}_{0}=(-0.91,-0.54,0.91)$. The SVD truncation threshold $e=\sum_{i=1}^{n}\sigma_{i}/\sum_{i=1}^{N}\sigma_{i}$ was set to $99\%$, where $\sigma_{i}$ are the singular values, and $N$ is the total number of snapshots. The $99\%$ threshold resulted in $n=11$ low-rank approximation. We then computed the DMD modes $\boldsymbol{\Phi}$ according to procedure explained in reference [8]. The DMD-reconstructed data was then computed using: $\mathbf{u}(t)\simeq\sum_{i=1}^{n}\boldsymbol{\Phi}_{k}\exp(\omega_{k}t)b_{k}$, where $b_{k}$ is the initial amplitude of each mode, $\boldsymbol{\Phi}_{k}$’s are the DMD eigenvectors and $\omega_{k}$’s are the eigenvalues associated to DMD modes. In Figure 5, we compare the reconstruction of the DMD with $n=11$ and DB with $r=4$. The DB reconstruction obtained by: $\mathbf{u}(t)\simeq\sum_{i=1}^{r}\mathbf{Y}_{s_{0},i}(t)\mathbf{u}_{i}(t)$, where $\mathbf{y}_{i}(t,\boldsymbol{\xi}_{0})$ is the stochastic coefficients associated with $\boldsymbol{\xi}_{0}=(-0.91,-0.54,0.91)$, i.e. a fixed row ($s$) of matrix $\mathbf{Y}$. It is clear that DMD fails to capture the transient instability, while DB reconstruction looks qualitatively accurate.

We perform a quantitative comparison between DMD and DB. We rank the $n=11$ DMD modes based on their growth rates from increasing to decreasing, and we perform DMD reconstruction $\mathbf{u}(t)\simeq\sum_{i=1}^{r}\boldsymbol{\Phi}_{k}\exp(\omega_{k}t)b_{k}$ for $r=2,4,6$ and 8. We also perform DB for the same reduction sizes. We then compute the reconstruction error: $\epsilon^{2}=\sum_{k=1}^{K}\|\mathbf{u}_{true}(t_{k})-\mathbf{u}(t_{k})\|^{2}/K$ where $K$ is the total number of snapshots and $\mathbf{u}(t_{k})$ is the reconstructed solution based on DMD or DB. In Table 1, the reconstruction errors $\epsilon$ for DB and DMD are shown. DMD is incapable of converging to the true solution as the size of DMD modes increases, however, DB converges to the true solution as the reductions size increases. The limitation of DMD for capturing transient dynamics has been recognized before; see for example [9], Chapter 1, Section 1.5.2.

One-hundred-dimensional random perturbation: In this section we present the results of dynamic basis reduction of one-hundred dimensional random initial condition with correlation length of $l_{c}=0.05$ (see Figure 3(c)). To capture the small scales of the random initial condition the state space is represented with Chebyshev polynomial of order $n=1024$. We draw Monte Carlo samples of the initial condition with sample sizes of $s=1000$ and $s=5000$. We note these sample size are quite small given the high dimensionality of the random space. However, as we demonstrate here the dynamic basis reduction extracts the low-dimensional structure of the observables despite the severely low number of samples. We choose reduction size of $r=6$. In Figure 6 the eigenvalues of the covariance matrix for the samples sizes are shown. As it is clear, the dynamic basis reduction with sample sizes of 1000 and 5000 are very close before and during the transition. This implies that despite the high dimensionality of the random space, the transient response has a low-dimensional structure — albeit a highly time-dependent one — that is effectively captured with the proposed reduction with a severely sparse samples.

3.3 Transient Vertical Jet

In this section we perform dynamic reduction for a two-dimensional transient vertical jet flushed into quiescent background flow. The jet is driven by random boundary conditions as described in this section. The schematic of the problem is shown in Figure 7.

To generate the samples, we performed direct numerical simulation (DNS) of the two-dimensional incompressible Navier-Stokes equations:

[TABLE]

where $\mathbf{u}(\mathbf{x},t;\boldsymbol{\xi})$ is the random velocity field with two components of $\mathbf{u}=(u,v)$ and the spatial coordinates of $\mathbf{x}=(x,y)$, and $p(\mathbf{x},t;\boldsymbol{\xi})$ is the random pressure field and $Re=500$ is the Reynolds number. We consider a stochastic jet velocity given by:

[TABLE]

The mean of the prescribed jet boundary condition is similar to previous studies [25, 26]. The random component of the boundary condition is represented by a uniform random field with the covariance kernel specified by:

[TABLE]

where $l_{c}=0.5$ and $\sigma=0.01$. The dimension $d=2$ is chosen such that 99% of the energy is captured by the Karhunen-Loève decomposition. We performed DNS using spectral/hp element method with $Ne=1604$ and spectral polynomial order of $p=5$. With these settings, the size of the discrete phase space is roughly $n=115,000$ when accounting for the two components of the velocity vector field. The time integration was performed using third-order stiffly stable scheme [27] with the time intervals of $\Delta t=5\times 10^{-3}$. We sample the DNS at PCM points with $Q=6$ quadrature points in each random dimension. Therefore, the total number of samples are $s=Q^{2}=36$.

As a reference reduction, we first computed the proper orthogonal decomposition (POD) modes. Since POD modes are often used in deterministic flows, we explain how these modes are computed in a stochastic setting. To do so, we first formed the matrix $\mathbf{A}$ that contains all mean-subtracted snapshots of all 36 samples, i.e. the entire DNS data. We then performed singular value decomposition of this matrix. The POD modes are then computed as the left singular eigenvectors of matrix $\mathbf{A}$. Therefore, the computed POD modes are an optimal basis for the entire DNS data in the second norm sense. In the middle three panels of Figure 7, the first three POD modes associated with largest singular values are shown.

We solved the evolution Eqs. 12 and 13 with $r=4$ with fourth order Runge-Kutta integration scheme. The time snapshots of matrix $\mathbf{T}$ were obtained by skipping every three time steps of the DNS, i.e. the time interval of dynamic basis evolution is $\Delta t=2\times 10^{-2}$ .

In Figure 8, we show the results of the state of the dynamics basis for three time snapshots at $t_{1}=2$, $t_{2}=6$ and $t_{3}=10$ time units, with each row showing the results at a fixed time. For the purpose of visualizing the mean, we solved a passive scalar transport equation along with Navier-Stokes equation with Schmidt number $Sc=1$. See [18] for more details for similar approach for visualization. In the leftmost column, the mean data is visualized by the scalar transport intensity. The mean shows the puff growing in size as the jet entrains the quiescent background flow while advecting upward. In the middle three columns, the first three most dominant dynamic bases are visualized by their vorticity: $\omega_{z_{i}}=\partial u_{i}/\partial y-\partial v_{i}/\partial x$. The evolution of the dynamic basis along with the data is clearly seen. The dynamic basis reveals the spatio-temporal coherent structures in the flow. Similar to the POD modes (see Figure 7), the higher modes of the dynamic basis show finer spatial structures. However, the POD modes are static and they are optimal in a time-average sense, as opposed to the dynamic basis that evolves according to $\mathbf{\dot{T}}$. We will compare the POD reduction error with the dynamic basis in the section. In the rightmost column, the low-dimensional attractor is visualized by plotting the iso-surface of joint probability density function $\pi_{y_{1}y_{2}y_{3}}$ of the three most dominant stochastic coefficients: $\mathbf{y}_{1},\mathbf{y}_{2}$ and $\mathbf{y}_{3}$. Since the DNS are sampled at PCM points, the polynomial chaos expansion is used to generate more samples for visualizing $\pi_{y_{1}y_{2}y_{3}}$ in the post-processing stage. The attractor clearly shows a low-dimensional manifold as it evolves in time.

In Figure 9, we make a comparison between POD and dynamic basis reduction errors. This error is obtained by:

[TABLE]

where for POD modes $\mathbf{U}$ represents the first $r$ dominant POD modes and $\mathbf{Y}$ is obtained via the orthogonal projection of the data onto $\mathbf{U}$, i.e.: \mathbf{Y}(t)=\big{<}\mathbf{U},\mathbf{T}(t)\big{>}. As shown in Figure 9, the reduction of dynamic basis $r=2$ is comparable to that of POD with $r=20$ for long time, and in the early evolution (i.e. $t<4$), the error of dynamic basis ($r=2)$ is smaller than POD with $r=50$. Similar behaviour is observed for dynamic basis $r=3,4$ and 5, where a ten-fold increase in POD reduction size is required to have a comparable reduction error with dynamic basis. We also observe that the dynamic basis error increases with time. This is due two effects: (1) as the jet advects upward, the uncertainty grows and more modes are required to better estimate the observations; and (2) the effect of unresolved modes in the dynamic basis reduction. In the case of drastic reduction, for example $r=2$, discarding the unresolved modes will lead to error increase over time. The effect of unresolved modes can be properly studied in the Mori-Zwanzig formalism [28]. As it can be seen in Figure 9 increasing the number of dynamic basis significantly decreases the error and the growth rate of the error.

In Figure 9, we show the results of different time discretization schemes for: (1) computing $\dot{\mathbf{T}}$ with two different schemes of first-order Euler explicit (EE1) and fourth-order central finite difference (FD4) and (2) the time advancement of $\mathbf{U}(t)$ and $\mathbf{Y}(t)$ via Equations 12 and 13 with EE1 and fourth-order Runge-Kutta (RK4). All four combinations of the above schemes are considered with the identifiers in the form of X-Y, where the first part of the identifier (X) shows the discretization scheme used for estimating $\dot{\mathbf{T}}$ and the second part of the identifier (Y) shows the time-advancement scheme for solving Equations 12 and 13. In particular we considered these schemes: EE1-EE1, EE1-RK4, FD4-EE1 and FD4-RK4. We used dynamic basis with the size of $r=5$ for all these schemes. The error of schemes EE1-EE1 and FD4-EE1 are very close – suggesting that the time-advancement error of solving solving Equations 12 and 13 is the primary source of error and increasing the accuracy of estimating $\dot{\mathbf{T}}$ does not improve the overall accuracy. Comparing the errors of EE1-EE1 with EE1-RK4 confirms that the error can be reduced by increasing the accuracy of time integrator of the evolution equations for $\mathbf{U}$ and $\mathbf{Y}$. However, comparing the errors of EE1-RK4 with FD4-RK4 demonstrates that the primary contributor of error for the scheme EE1-RK4 is the estimation of $\dot{\mathbf{T}}$. The totality of these observations show that both the accuracy of estimating $\dot{\mathbf{T}}$ and solving the evolution equations of the dynamic basis affect the reduction error potentially significantly. Obviously, if $\Delta t$ is decreased, the contribution of time-discretization error will become less significant compared to reduction error.

4 Conclusion

The objective of this paper is to develop a reduced description of the transient dynamics that relies on the observables of the dynamics. To this end, we present a variational principle whose optimality conditions lead to closed-form evolution equation for a set of time-dependent orthonormal modes and their associated coefficients, referred to as dynamic basis. Dynamic basis is driven by the time derivative of the random observations, and it does not require the underlying dynamical system that generates the observables. In that sense, the proposed reduction is purely observation-driven and may apply to systems whose models are not known.

The application of the dynamic basis was demonstrated for a variety of examples: (1) advection equation with random wave speed for which exact solution for the dynamic basis was derived and it was used for validation purposes. (2) Kuramoto-Sivashisky equation subject to random initial condition where it was shown that the dynamic basis captures the transient intrinsic instability - even in the presence of high-dimensional random initial condition. The performance of the dynamic basis reduction was compared against the dynamic mode decomposition. (3) Navier-Stokes equation subject to random boundary condition for which we compared the effectiveness of the dynamic basis versus proper orthogonal decomposition, as an example of static-basis reductions.

The dynamic basis evolution equation has a linear computational complexity with respect to the number of observables and the number of random realizations. This is a favorable feature if applications of large dimension are of interest, for example weather forecast, ocean models and turbulent flows.

Acknowledgement

The authors would like to thank Prof. Karniadakis for useful comments. The author gratefully acknowledge the allocated computer time on the Stampede supercomputers, awarded by the XSEDE program, allocation number ECS140006. The author has been supported through a NASA Grant number 80NSSC18M0150, and an award by American Chemical Society, Grant number 59081-DN19.

Appendix A. Optimality condition of the variational principle

In the following we show that the first-order optimality condition of the variational principle 7 leads to closed-form evolution equations for $\mathbf{U}$ and $\mathbf{Y}$. We first observe that:

[TABLE]

First order optimality condition requires $\dfrac{\partial\mathcal{G}}{\partial\dot{\mathbf{u}}_{k}}=0$ and $\dfrac{\partial\mathcal{G}}{\partial\dot{\mathbf{y}}_{k}}=0$ for $k=1,2,\dots,r$. We first observe that the gradient of the $\mathcal{G}$ with respect to $\dot{\mathbf{u}}_{k}$ can be expressed as:

[TABLE]

To eliminate $\lambda_{kj}$, we take the inner product of the above equation by $\mathbf{u}_{l}$:

[TABLE]

where we have made use of the definition of \boldsymbol{\phi}_{lj}(t):=\big{<}\mathbf{u}_{l},\dot{\mathbf{u}}_{j}\big{>}. Rearranging the above equation results in:

[TABLE]

After multiplying equation 20 from left by $\dfrac{\mathbf{W}_{\mathbf{x}}^{-1}}{2}$ and substituting $\lambda_{jk}$ from the above equation we obtain:

[TABLE]

In the vectorized form the above equation can be written as:

[TABLE]

where $\mathbf{C}_{ij}=\mathbf{y}_{i}^{T}\mathbf{W}_{\boldsymbol{\xi}}\mathbf{y}_{j}$ is the covariance matrix. Similarly,

[TABLE]

In the above expression we have made use of $\left(\dot{\mathbf{u}}^{T}_{i}\mathbf{W}_{\mathbf{x}}\mathbf{u}_{j}\right)\left(\mathbf{y}^{T}_{i}\mathbf{W}_{\boldsymbol{\xi}}\dot{\mathbf{y}}_{j}\right)=\left(\mathbf{y}^{T}_{i}\mathbf{W}_{\boldsymbol{\xi}}\dot{\mathbf{y}}_{j}\right)\left(\dot{\mathbf{u}}^{T}_{i}\mathbf{W}_{\mathbf{x}}\mathbf{u}_{j}\right)$ since both $\left(\dot{\mathbf{u}}^{T}_{i}\mathbf{W}_{\mathbf{x}}\mathbf{u}_{j}\right)$ and $\left(\mathbf{y}^{T}_{i}\mathbf{W}_{\boldsymbol{\xi}}\dot{\mathbf{y}}_{j}\right)$ are scalars. By multiplying the above equation 21 by $\dfrac{\mathbf{W}_{\boldsymbol{\xi}}^{-1}}{2}$ from left and replacing $\mathbf{u}^{T}_{k}\mathbf{W}_{\mathbf{x}}\mathbf{u}_{j}=\delta_{kj}$ and $\dot{\mathbf{u}}^{T}_{j}\mathbf{W}_{\mathbf{x}}\mathbf{u}_{k}=\boldsymbol{\phi}_{kj}=-\boldsymbol{\phi}_{jk}$, we obtain:

[TABLE]

In the vectorized form the above equation can be expressed as:

[TABLE]

Appendix B. Proof of equivalence of reductions

Theorem .1

*Suppose that $\{\mathbf{Y},\mathbf{U}\}$ and $\{\widetilde{\mathbf{Y}},\widetilde{\mathbf{U}}\}$ satisfy the evolution equations (8 and 9) with different choices of $\boldsymbol{\phi}$ functions denoted by $\boldsymbol{\phi}\in\mathbb{R}^{r\times r}$ and $\widetilde{\boldsymbol{\phi}}\in\mathbb{R}^{r\times r}$, respectively. We also assume that the two reductions are initially equivalent, *i.e.

• $\mathbf{U}(0)=\widetilde{\mathbf{U}}(0)\mathbf{R}_{0}$ and $\mathbf{Y}(0)=\widetilde{\mathbf{Y}}(0)\mathbf{R}_{0}$, where $R_{0}\in\mathbb{R}^{r\times r}$ is an orthogonal rotation matrix. Then $\{\mathbf{Y},\mathbf{U}\}$ and $\{\widetilde{\mathbf{Y}},\widetilde{\mathbf{U}}\}$ are equivalent for $t>0,$ with a rotation matrix $\mathbf{R}(t)\in\mathbb{R}^{r\times r}$ governed by the matrix differential equation $\dot{\mathbf{R}}=\mathbf{R}\boldsymbol{\boldsymbol{\phi}}-\widetilde{\boldsymbol{\phi}}\mathbf{R}$.*

proof: We replace $\mathbf{U}=\widetilde{\mathbf{U}}\mathbf{R}$ and $\mathbf{Y}=\widetilde{\mathbf{Y}}\mathbf{R}$ into the evolution equations (8 and 9). The evolution of the dynamic basis can be written as:

[TABLE]

Therefore:

[TABLE]

where $\widetilde{\mathbf{C}}=\mathbb{E}[\widetilde{\mathbf{Y}}^{T}\widetilde{\mathbf{Y}}]=\mathbf{R}\mathbf{C}\mathbf{R}^{T}$ and therefore, $\widetilde{\mathbf{C}}^{-1}=\mathbf{R}\mathbf{C}^{-1}\mathbf{R}^{T}$, where $\mathbf{C}$ are $\widetilde{\mathbf{C}}$ are similar matrices and have the same eigenvalues. In the last equation we have used $\dot{\mathbf{R}}=\mathbf{R}\boldsymbol{\boldsymbol{\phi}}-\widetilde{\boldsymbol{\phi}}\mathbf{R}$, and made use of the identity $\mathbf{R}^{-1}=\mathbf{R}^{T}$. The evolution of the stochastic coefficients can be written as:

[TABLE]

Therefore:

[TABLE]

This completes the proof.