Space-time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems

1 Introduction

Many computational models for physics simulations are formulated as linear dynamical systems. Examples of linear dynamical systems include the computational model for the signal propagation and interference in electric circuits, storm surge prediction models before an advancing hurricane, vibration analysis in large structures, thermal analysis in various media, neuro-transmission models in the nervous system, various computational models for micro-electro-mechanical systems, and various particle transport simulations. Depending on the complexity of geometries and desirable fidelity level, these problems become easily large-scale problems. For example, the Boltzmann Transport Equation (BTE) has seven independent variables, i.e., three spatial variables, two directional variables, one energy variable, and one time variable. It is not hard to see that the BTE can easily lead to a high dimensional discretized problem. Additionally, the complex geometry (e.g., reactors with thousands of pins and shields) can lead to a large-scale problem. As an example, a problem with 20 angular directions, a cubit spatial domain of 100 x 100 x 100 elements, 16 energy groups, and 100 time steps leads to 32 billion unknowns. The large-scale hinders a fast forward solve and prevents the multi-query setting problems, such as uncertainty quantification, design optimization, and parameter study, from being tractable. Therefore, developing a Reduced Order Model (ROM) that accelerates the solution process without losing much accuracy is essential.

There are several model order reduction approaches available for linear dynamical systems: (i) Balanced truncation [29, 28]

in control theory community is the most famous one. It has explicit error bounds and guarantees stability. However, it requires the solution of two Lyapunov equations to construct bases, which is a formidable task in large-scale problems. (ii) The moment-matching methods

[4, 19]

provide a computationally efficient framework using Krylov subspace techniques in an iterative fashion where only matrix-vector multiplications are required. The optimal

H_{2}

tangential interpolation for nonparametric systems

[19] is also available. The most crucial part of the moment-matching methods is location of samples where moments are matched. Also, it is not a data-driven approach, meaning that no data is used to construct ROM. (iii) Proper Generalized Decomposition (PGD) [2] was first developed as a numerical method of solving boundary value problems, later extended to dynamical problems [3]. The main assumption of the method is a separated solution representation in space and time, which gives a way for an efficient solution procedure. Therefore, it is considered as a model reduction technique. However, PGD is not a data-driven approach.

Often there are many data available either from experiments or high-fidelity simulations. Those data contain valuable information about the system of interest. Therefore, data-driven methods can maximize the usage of the existing data and enable the construction of an optimal ROM. Dynamic Mode Decomposition (DMD) is a data-driven approach that generates reduced modes that embed an intrinsic temporal behavior. It was first developed by Peter Schmid in [34] and populated by many other scientists. For more detailed information about DMD, we refer to this preprint [37]. As another data-driven approach, Proper Orthogonal decomposition (POD) [6] gives the data-driven optimal basis through the method of snapshots. However, most of POD-based ROMs for linear dynamical systems apply spatial projection only. Temporal complexity is still proportional to temporal discretization of its high-fidelity model. In order to achieve an optimal reduction, a space–time ROM needs to be built where both spatial and temporal projections are applied. In literature, some space–time ROMs are available [14, 38, 41, 40], but they are only applied to small-scale problems.

Recently, several ROM techniques have been applied to various types of transport equations. Starting with Wols’s work that used the POD to simulate the dynamics of an accelerator driven system (ADS) in 2010 [39]

, the publications on ROMs in transport problems have increased in number. For example, a POD-based ROM for eigenvalue problems to calculate dominant eigenvalues in reactor physics applications is developed by Buchan, et al. in

[10]. The corresponding Boltzmann transport equation was re-casted into its diffusion form, where some of dimensions were eliminated. Sartori, et al. in [33] also applied the POD-based ROM to the diffusion form and compared it with a modal method to show that the POD was superior to the modal method. Reed and Robert in [32] also used POD to expand energy that replaced the traditional Discrete Legendre Polynomials (DLP) or modified DLP. They showed that a small number of POD energy modes could capture many-group fidelity. Buchan, et al. in [11] also developed a POD-based reduced order model to efficiently resolve the angular dimension of the steady-state, mono-energetic Boltzmann transport equation. Behne, Ragusa, and Morel [5] applied POD-based reduced order model to accelerate steady state

S_{n}

radiation multi-group energy transport problem. A Petrov-Galerkin projection was used to close the reduced system. Coale and Anistratov in [15] replaces a high-order (HO) system with POD-based ROM in high-order low-order (HOLO) approach for Thermal Radiative Transfer (TRT) problems.

There have been some interesting DMD works for transport problems. McClarren and Haut in [27]

used DMD to estimate the slowly decaying modes from Richardson iteration and remove them from the solution. Hardy, Morel, and Cory

[20] also explored DMD to accelerate the kinetics of subcritical metal systems without losing much accuracy. It was applied to a three-group diffusion model in a bare homogeneous fissioning sphere. An interesting work by Star, et al. in [36] exists, using the DMD method to identify non-intrusive POD-based ROM for the unsteady convection-diffusion scalar transport equation. Their approach applied the DMD method to reduced coordinate systems.

Several papers are found to use the PGD for transport problems. For example, Prince and Ragusa [30] applied the Proper Generalized Decomposition (PGD) to steady-state mono-energetic neutron transport equations where $S_{n}$ angular flux was sought as a finite sum of separable one-dimensional functions. However, the PGD was found to be ineffective for pure absorption problems because a large number of terms were required in the separated representation. Prince and Ragusa [31] also applied the PGD for uncertainty quantification process in the neutron diffusion–reaction problems. Dominesey and Ji used the PGD method to separate space and angle in [16] and to separate space and energy in [17].

However, all these model order reduction techniques for transport equations apply only spatial projection, ignoring the potential reduction in temporal dimension. In this paper, a Space–Time Reduced Order Model (ST-ROM) is developed and applied to large-scale linear dynamical problems. Our ST-ROM achieves complexity reduction in both space and time dimension, which enables a great maximal speed-up and accuracy. It is amenable to any time integrators. It follows the framework initially published in [14], but makes a new contribution by discovering a block structure in space–time basis that enables efficient implementation of the ST-ROM. The block structure in space–time basis allows us not to build a space–time basis explicitly. It enables the construction of space–time reduced operators with small additional costs to the space-only reduced operators. In turn, this allows us to apply the space–time ROM to a large-scale linear dynamical problem.

The paper is organized in the following way: Section 1.1 introduces useful mathematical notations that will be used throughout the paper. Section 2 describes a parametric linear dynamical system and how to solve high-fidelity model in a classical time marching fashion. The full-order space–time formulation is also presented in Section 2 to be reduced to form our space–time ROM in Section 3.2. Section 3 introduces both spatial and spatiotemporal ROMs. The basis generation is described in Section 4 where the traditional POD and incremental POD are explained in Sections 4.1 and 4.2, respectively. Section 4.3 reveals a block structure of the space–time reduced basis and derive each space–time reduced operators in terms of the blocks. We apply our space–time ROM to a large-scale linear dynamical problem, i.e., a neutron transport simulation of solving BTE. Section 6 explains a discretization derivation of the Boltzmann transport equation, using multigroup energy discretization, associated Legendre polynomials for surface harmonic, simple corner balance discretization for space and direction, and the discrete ordinates method. Finally, we present our numerical results in Section 7 and conclude the paper with summary and future works in Section 8.

1.1 Notations

We review some of the notation used throughout the paper. An $ℓ_{2}$ norm is denoted as $∥ \cdot ∥$ . For matrices $A \in R^{m \times n}$ and $B \in R^{k \times l}$ , the Kronecker (or tensor) product of $A$ and $B$ is the $m k \times n l$ matrix denoted by

A \otimes B \equiv ⎛ ⎜ ⎜ ⎝ \begin{matrix} a_{11} B & \dots & a_{1 n} B ⋮ & ⋱ & ⋮ a_{m 1} B & \dots & a_{m n} B \end{matrix} ⎞ ⎟ ⎟ ⎠,

where $A = (a_{i j})$ . Kronecker products have many interesting properties. We list here the ones relevant to our discussion:

If $A$ and $B$ are nonsingular, then $A \otimes B$ is nonsingular with $(A \otimes B)^{- 1} = A^{- 1} \otimes B^{- 1}$ ,
$(A \otimes B)^{T} = A^{T} \otimes B^{T}$ ,
Given matrices $A$ , $B$ , $C$ , and $D$ , $(A \otimes B) \cdot (C \otimes D) = A C \otimes B D$ , as long as both sides of the equation make sense,
$(A + B) \otimes C = A \otimes C + B \otimes C$ , and
$A \otimes (B + C) = A \otimes B + A \otimes C$ .

2 Linear dynamical systems

Parametric continuous dynamical systems that are linear in state are considered:

	$˙ u (t; μ)$	$= A (μ) u (t; μ) + B (μ) f (t; μ), u (0; μ) = u_{0} (μ),$		(1)
	$y (t; μ)$	$= C (μ)^{T} u (t; μ),$		(2)

where $μ \in Ω_{μ} \subset R^{n_{μ}}$ denotes a parameter vector, $u : [0, T] \times R^{n_{μ}} \to R^{N_{s}}$ denotes a time dependent state variable function, $u_{0} : R^{n_{μ}} \to R^{N_{s}}$ denotes an initial state, $f : [0, T] \times R^{n_{μ}} \to R^{N_{i}}$ denotes a time dependent input variable function, and $y : [0, T] \times R^{n_{μ}} \to R^{N_{o}}$ denotes a time dependent output variable function. The system operations, i.e., $A : R^{n_{μ}} \to R^{N_{s} \times N_{s}}$ , $B : R^{n_{μ}} \to R^{N_{s} \times N_{i}}$ , and $C : R^{n_{μ}} \to R^{N_{s} \times N_{o}}$ , are real valued matrices, independent of state variables. We assume that the dynamical system above is stable, i.e., the eigenvalues of $A$ have strictly negative real parts.

Our methodology works for any time integrators, but for the illustration purpose, we apply a backward Euler time integrator to Eq. (1). At $k$ th time step, the following system of equations is solved:

(I_{N_{s}} - Δ t^{(k)} A (μ)) u^{(k)} = u^{(k - 1)} + Δ t^{(k)} B (μ) f^{(k)} (μ),

(3)

where $I_{N_{s}} \in R^{N_{s} \times N_{s}}$

denotes an identity matrix,

Δ t^{(k)}

denotes

k

th time step size with

T = \sum_{k = 1}^{N_{t}} Δ t^{(k)}

, and

u^{(k)} (μ)

and

f^{(k)} (μ)

denote state and input vectors at

k

th time step,

t_{(k)} = \sum_{j = 1}^{k} Δ t^{(j)}

, respectively. A Full Order Model (FOM) solves Eq. (3) every time step. The spatial dimension,

N_{s}

, and the temporal dimension,

N_{t}

can be very large, which leads to a large-scale problem. We introduce how to reduce the high dimensionality in Section 3.

The single time step formulation in Eq. (3) can be equivalently re-written in the following discretized space-time formulation:

A^{st} (μ) u^{st} (μ) = f^{st} (μ) + u_{0}^{st} (μ),

(4)

where the space–-time system matrix, $A^{st} : R^{n_{μ}} \to R^{N_{s} N_{t} \times N_{s} N_{t}}$ , the space–time state vector, $u^{st} : R^{n_{μ}} \to R^{N_{s} N_{t}}$ , the space–time input vector, $f^{st} : R^{n_{μ}} \to R^{N_{s} N_{t}}$ , and the space–time initial vector, $u_{0}^{st} : R^{n_{μ}} \to R^{N_{s} N_{t}}$ , are defined respectively as

A^{st} (μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} I_{N_{s}} - Δ t^{(1)} A (μ) - I_{N_{s}} & I_{N_{s}} - Δ t^{(2)} A (μ) ⋱ & ⋱ - I_{N_{s}} & I_{N_{s}} - Δ t^{(N_{t})} A (μ) \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦,

(5)

u^{st} (μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} u^{(1)} (μ) u^{(2)} (μ) ⋮ u^{(N_{t})} (μ) \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦, f^{st} (μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} Δ t^{(1)} B (μ) f^{(1)} (μ) Δ t^{(2)} B (μ) f^{(2)} (μ) ⋮ Δ t^{(N_{t})} B (μ) f N_{t} (μ) \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦, u_{0}^{st} (μ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} u_{0} (μ) 0 ⋮ 0 \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ .

(6)

A lower block-triangular matrix structure of $A^{st}$ comes from the backward Euler time integration scheme. Other time integrators will give other sparse block structures. No one will solve this space–time system directly because the specific block structure of $A^{st}$ lets one to solve the system in time-marching fashion. However, if the space–time formulation in Eq. (4) can be reduced and solved efficiently, then one might be interested in solving the reduced space–time system in its whole. Section 3.2 shows such a reduction is possible.

3 Reduced order models

We consider a projection-based reduced order model for linear dynamical systems. Section 3.1 shows a typical spatial reduced order model. A space–time reduced order model is described in Section 3.2.

3.1 Spatial reduced order models

A projection-based spatial reduced order model approximates the state variables as a linear combination of a small number of spatial basis vectors, ${ϕ_{1}^{s}, \dots, ϕ_{n_{s}}^{s}}$ , where $ϕ_{k}^{s} \in R^{N_{s}}$ , $k \in N (n_{s})$ with $N (n_{s}) \equiv {1, \dots, n_{s}}$ , $n_{s} ≪ N_{s}$ , i.e.,

u (t; μ) \approx ~ u (t) \equiv u_{ref} (μ) + Φ_{s}^u (t; μ),

(7)

where the spatial basis, $Φ_{s} \in R^{N_{s} * n_{s}}$ with $Φ_{s}^{T} Φ_{s} = I_{n_{s}}$ , is defined as

Φ_{s} \equiv [\begin{matrix} ϕ_{1}^{s} & \dots & ϕ_{n_{s}}^{s} \end{matrix}],

(8)

a reference state is denoted as $u_{ref} (μ) \in R^{N_{s}}$ , and a time-dependent reduced coordinate vector function is defined as $^u : R^{n_{μ}} \to R^{n_{s}}$ . Substituting (7) to (1), gives an over-determined system of equations:

Φ_{s} ˙^u (t; μ) = A (μ) u_{ref} (μ) + A (μ) Φ_{s}^u (t; μ) + B (μ) f (t, μ),

(9)

which can be closed, for example, by Galerkin projection, i.e., left-multiplying both sides of (9) and initial condition in (1) by $Φ_{s}^{T}$ , giving the reduced system of equations and initial conditions:

˙^u (t; μ) =^A (μ)^u (t; μ) +^B (μ) f (t; μ) + {^u}_{ref} (μ),^u (0; μ) = {^u}_{0} (μ)

(10)

where $^A (μ) \equiv Φ_{s}^{T} A (μ) Φ_{s} \in R^{n_{s} \times n_{s}}$ denotes a reduced system matrix, $^B (μ) \equiv Φ_{s}^{T} B (μ) \in R^{n_{s} \times N_{i}}$ denotes a reduced input matrix, ${^u}_{ref} (μ) \equiv Φ_{s}^{T} A (μ) u_{ref} (μ) \in R^{n_{s}}$ denotes a reduced reference state vector, and ${^u}_{0} (μ) \equiv Φ_{s}^{T} (u_{0} (μ) - u_{ref} (μ)) \in R^{n_{s}}$ denotes a reduced initial condition. Once $Φ_{s}$ and $u_{ref}$ are known, the reduced operators, $^A$ , $^B$ , ${^u}_{ref}$ , and ${^u}_{0}$ can be pre-computed. With the pre-computed operators, the system (10) can be solved fairly quickly, for example, by applying the backward Euler time integrator:

(I_{n_{s}} - Δ t^{(k)}^A (μ)) {^u}^{(k)} (μ) = {^u}^{(k - 1)} (μ) + Δ t^{(k)}^B (μ) f^{(k)} (μ) + Δ t^{(k)} {^u}_{ref} (μ),

(11)

where ${^u}^{(k)} (μ) \equiv^u (t_{(k)}; μ) \in R^{n_{s}}$ . Then the output vector, $y^{(k)} (μ) \equiv y (t_{(k)}; μ)$ , can be computed as

y^{(k)} (μ) \equiv^C (μ)^{T} {^u}^{(k)} (μ) + C (μ)^{T} u_{ref} (μ),

(12)

where $^C (μ) \equiv Φ_{s}^{T} C (μ) \in R^{n_{s} \times N_{o}}$ denotes a reduced output matrix. The usual choices for $u_{ref}$ include $0$ , $u_{0}$ , and some kind of average quantities. Note that if $u_{0}$ is used as $u_{ref}$ , then ${^u}_{0} = 0$ , independent of $μ$ , which is convenient. The POD for generating the spatial basis is described in Section 4.

3.2 Space–time reduced order models

The space–time formulation, (4), can be reduced by approximating the space–time state variables as a linear combination of a small number of space–time basis vectors, ${ϕ_{1}^{st}, \dots, ϕ_{n_{s} n_{t}}^{st}}$ , where $ϕ_{k}^{st} \in R^{N_{s} N_{t}}$ , $k \in N (n_{s} n_{t})$ with $n_{s} n_{t} ≪ N_{s} N_{t}$ , i.e.,

u^{st} (μ) \approx {~ u}^{st} (μ) \equiv Φ_{st} {^u}^{st} (μ),

(13)

where the space–time basis, $Φ_{st} \in R^{N_{s} N_{t} * n_{s} n_{t}}$ is defined as

Φ_{st} \equiv [\begin{matrix} ϕ_{1}^{st} & \dots & ϕ_{i + n_{s} (j - 1)}^{st} & \dots ϕ_{n_{s} n_{t}}^{st} \end{matrix}],

(14)

where $i \in N (n_{s})$ , $j \in N (n_{t})$ . The space–time reduced coordinate vector function is denoted as ${^u}^{st} : R^{n_{μ}} \to R^{n_{s} n_{t}}$ . Substituting (13) to (4), gives an over-determined system of equations:

A^{st} (μ) Φ_{st} {^u}^{st} (μ) = f^{st} (μ) + u_{0}^{st} (μ)

(15)

which can be closed, for example, by Galerkin projection, i.e., left-multiplying both sides of (15) by $Φ_{st}^{T}$ , giving the reduced system of equations:

{^A}^{st} (μ) {^u}^{% st} (μ) = {^f}^{st} (μ) +^u_{0}^{st} (μ)

(16)

where ${^A}^{st} (μ) \equiv Φ_{st}^{T} A^{st} (μ) Φ_{st} \in R^{n_{s} n_{t} \times n_{s} n_{t}}$ denotes a reduced space–time system matrix, ${^f}^{st} (μ) \equiv Φ_{st}^{T} f^{st} (μ) \in R^{n_{s} n_{t}}$ denotes a reduced space–time input vector, and ${^u}_{0}^{st} (μ) \equiv Φ_{st}^{T} u_{0}^{st} (μ) \in R^{n_{s} n_{t}}$ denotes a reduced space–time initial state vector. Once $Φ_{st}$ is known, the space–time reduced operators, ${^A}^{st}$ , ${^f}^{st}$ , and ${^u}_{0}^{st}$ , can be pre-computed, but its computation involves the number of operations in $O (N_{s} N_{t})$ , which can be large. Sec. 4.3 explores a block structure of $Φ_{st}$ that shows an efficient way of constructing reduced space–time operators without explicitly forming the full-size space–time operators, such as $Φ_{st}$ , $A^{st}$ , $f^{st}$ , and $u_{0}^{st}$ .

4 Basis generation

4.1 Proper orthogonal decomposition

Figure 1: Illustration of spatial and temporal bases construction, using SVD with $n_{μ} = 3$ . The right singular vector, $v_{i}$ , describes three different temporal behaviors of a left singular basis vector $w_{i}$ , i.e., three different temporal behaviors of a spatial mode. Each temporal behavior is denoted as $v_{i}^{1}$ , $v_{i}^{2}$ , and $v_{i}^{3}$ .

We follow the method of snapshots first introduced by Sirovich [35]. Let $P = {μ_{1}, \dots, μ_{n_{μ}}}$ be a set of parameter samples where we run full order model simulations. Let $U_{p} \equiv [\begin{matrix} u^{(1)} (μ_{p}) & \dots & u^{(N_{t})} (μ_{p}) \end{matrix}] \in R^{N_{s} \times N_{t}}$ , $p \in N (n_{μ})$ , be a full order model state solution matrix for a sample parameter value, $μ_{p} \in Ω_{μ}$ . Then a snapshot matrix, $U \in R^{N_{s} \times n_{μ} N_{t}}$ , is defined by concatenating all the state solution matrices, i.e.,

U \equiv [\begin{matrix} U_{1} & \dots & U_{n_{μ}} \end{matrix}] .

(17)

The spatial basis from POD is an optimally compressed representation of $range (U)$ in a sense that it minimizes the difference between the original snapshot matrix and the projected one onto the subspace spanned by the basis, $Φ_{s}$ :

minimize Φ_{s} \in R^{N_{s} \times n_{s}}, Φ_{s}^{T} Φ_{s} = I_{n_{s}}

{∥ ∥ U - Φ_{s} Φ_{s}^{T} U ∥ ∥}_{s}^{2},

(18)

where $∥ \cdot ∥_{F}$ denotes the Frobenius norm. The solution of POD can be obtained by setting $Φ_{s} = W (:, 1 : n_{s})$ , $n_{s} < n_{μ} N_{t}$ , in MATLAB notation, where $W$ is the left singular matrix of the following thin Singular Value Decomposition (SVD) with $ℓ \equiv min (N_{s}, n_{μ} N_{t})$ :

	$U$	$= W Σ V^{T}$		(19)
		$= ℓ \sum i = 1 σ_{i} w_{i} v_{i}^{T}$		(20)

where $W \in R^{N_{s} \times ℓ}$ and $V \in R^{n_{μ} N_{t} \times ℓ}$ are orthogonal matrices and $Σ \in R^{ℓ \times ℓ}$ is a diagonal matrix with singular values on its diagonal. The equivalent summation form is written in (20), where $σ_{i} \in R$ is $i$ th singular value, $w_{i}$ and $v_{i}$ are $i$ th left and right singular vectors, respectively. Note that $v_{i}$ describes $n_{μ}$ different temporal behavior of $w_{i}$ . For example, Figure 1 illustrates the case of $n_{μ} = 3$ , where $v_{i}^{1}$ , $v_{i}^{2}$ , and $v_{i}^{3}$ describe three different temporal behavior of a specific spatial basis vector, i.e., $w_{i}$ . For general $n_{μ}$ , we note that $v_{i}$ describes $n_{μ}$ different temporal behavior of $i$ th spatial basis vector, i.e., $ϕ_{i}^{s} = w_{i}$ . We set $Υ_{i} = [\begin{matrix} v_{i}^{1} & \dots & v_{i}^{n_{μ}} \end{matrix}]$ to be $i$ th temporal snapshot matrix, where $v_{i}^{k} \equiv v_{i} (1 + (k - 1) N_{t} : k N_{t})$ for $k \in N (n_{μ})$ . We apply SVD on $Υ_{i}$ :

Υ_{i} = Λ_{i} Σ_{i} Ψ_{i}^{T} .

(21)

Then, the temporal basis for $i$ th spatial basis vector can be set $Φ_{t}^{i} = Λ_{i} (:, 1 : n_{t})$ in MATLAB notation. Finally, a space–time basis vector, $ϕ_{i + n_{s} (j - 1)}^{st} \in R^{N_{s} N_{t}}$ , in (14) can be constructed as

ϕ_{i + n_{s} (j - 1)}^{st} = ϕ_{i j}^{t} \otimes ϕ_{i}^{s},

(22)

where $ϕ_{i}^{s} \equiv Φ_{s} (:, i) \in R^{N_{s}}$ denotes $i$ th spatial basis vector and $ϕ_{i j}^{t} \equiv Φ_{t}^{i} (:, j) \in R^{N_{t}}$ denotes $j$ th temporal basis vector that describes a temporal behavior of $ϕ_{i}^{s}$ . The computational cost of SVD for the snapshot matrix, $U \in R^{N_{s} \times n_{μ} N_{t}}$ , assuming $N_{s} ≫ n_{μ} N_{t}$ , is $O (N_{s}^{2} n_{μ} N_{t})$ and the computational cost of SVD for $n_{s}$ temporal snapshot matrices, $Υ_{i} \in R^{N_{t} \times n_{μ}}$ , $i \in N (n_{s})$ , $N_{t} ≫ n_{μ}$ is $O (N_{t}^{2} n_{s} n_{μ})$ . For a large-scale problem, this may be a formidable task. Thus, we use an incremental SVD where a rank one update of existing SVD is achieved with much more memory-efficient way than the thin SVD in Eq. (19). The incremental SVD procedure is explained in Section 4.2.

POD is related to the principal component analysis in statistics

[22] and Karhunen–Loève expansion [26] in stochastic analysis. Since the objective function in (18) does not change even though

Φ_{s}

is post-multiplied by an arbitrary

n_{s} \times n_{s}

orthogonal matrix, the POD procedure seeks the optimal

n_{s}

-dimensional subspace that captures the snapshots in the least-squares sense. For more details on POD, we refer to [6, 21, 23].

4.2 Incremental space–time reduced basis

An incremental SVD is an efficient way of updating the existing singular value decomposition when a new snapshot vector, i.e., a column vector, is added. For a time dependent problem, we start with a first time step solution with a first parameter vector, i.e., $u^{(1)} (μ_{1})$ . If its norm is big enough (i.e., $∥ u^{(1)} (μ_{1}) ∥ > ϵ_{SVD}$ ), then we set the first singular value $σ_{1} = ∥ u^{(1)} (μ_{1}) ∥$ , the first left singular vector be the normalized first snapshot vector, i.e., $w_{1} = u^{(1)} (μ_{1}) / σ_{1}$ , and the right singular vector be $v_{1} = 1$ . Otherwise, we set them empty, i.e., $σ_{1} = []$ , $w_{1} = []$ , and $v_{1} = []$ . This initializing process is described in Algorithm 1. We pass $k$ to initializingIncrementalSVD function as an input argument to indicate $k$ th snapshot vector is being handled. Also, the rank of $W_{k}$ is denoted as $r_{k}$ . In general, $r_{k} \neq k$ because a snapshot vector will not be included if it is too small (i.e., Line 1 in Algorithm 1) or it is linearly dependent on the existing basis (i.e., Line 9 and 13 in Algorithm 2) or it generates a small eigenvalue (i.e., Line 18 in Algorithm 2).

Let’s assume that we have $(k - 1)$ th SVD from previous $k - 1$ snapshot vectors, i.e., $W_{k - 1} Σ_{k - 1} V_{k - 1}^{T}$ , whose rank is $r_{k - 1}$ . If a new snapshot vector, $u$ (e.g., $k$ th time step solution with the first sample parameter value, $u^{(k)} (μ_{1})$ ) needs to be added to the existing SVD, the following factorization can be used [9]:

	$[\begin{matrix} W_{k - 1} Σ_{k - 1} V_{k - 1}^{T} & u \end{matrix}]$			(23)
		$= [\begin{matrix} W_{k - 1} & j \end{matrix}] [\begin{matrix} Σ_{k - 1} & ℓ 0 & p \end{matrix}] {[\begin{matrix} V_{k - 1} & 0 0 & 1 \end{matrix}]}^{T},$		(24)

where $ℓ = W_{k - 1}^{T} u \in R^{r_{k - 1}}$ denotes a reduced coordinate of $u$ that is projected onto the subspace spanned by $W_{k - 1}$ , $p = ∥ u - W_{k - 1} ℓ ∥$ denotes the norm of the difference between $u$ and the projected one, and $j = (u - W_{k - 1} ℓ) / p \in R^{N_{s}}$ denotes a new orthogonal vector due to the incoming vector, $u$ . Note that the left and right matrices of the factorization, i.e., $[\begin{matrix} W_{k - 1} & j \end{matrix}] \in R^{N_{s} \times (r_{k - 1} + 1)}$ and $[\begin{matrix} V_{k - 1} & 0 0 & 1 \end{matrix}] \in R^{k \times (r_{k - 1} + 1)}$ are orthogonal matrices. Let $Q \in R^{(r_{k - 1} + 1) \times (r_{k - 1} + 1)}$ denote the middle matrix of the factorization, i.e.,

Q = [\begin{matrix} Σ_{k - 1} & ℓ 0 & p \end{matrix}] .

(25)

The matrix, $Q$ , is almost diagonal except for $ℓ$ in the upper right block, i.e., one column bordered diagonal. Its size is not in $O (N_{s})$ . Thus, the SVD of $Q$ is computationally cheap, i.e., $O ((r_{k - 1} + 1)^{3})$ . Let the SVD of $Q$ be

(26)

where $¯ ¯¯¯¯ ¯ W \in R^{(r_{k - 1} + 1) \times (r_{k - 1} + 1)}$ denotes the left singular matrix, denotes the singular value matrix, and $¯ ¯¯ ¯ V \in R^{(r_{k - 1} + 1) \times (r_{k - 1} + 1)}$ denotes the right singular matrix of $Q$ . Replacing $Q$ in Eq. (24) with (26) gives

	$[\begin{matrix} W_{k - 1} Σ_{k - 1} V_{k - 1}^{T} & u \end{matrix}]$			(27)
		$= W_{k} Σ_{k} V_{k}^{T},$		(28)

where denotes the updated left singular matrix, $Σ_{k} = ¯ ¯¯ ¯ Σ \in R^{(r_{k - 1} + 1) \times (r_{k - 1} + 1)}$ denotes the updated singular value matrix, and $V_{k} = [\begin{matrix} V_{k - 1} & 0 0 & 1 \end{matrix}] ¯ ¯¯ ¯ V \in R^{k \times (r_{k - 1} + 1)}$ denotes the updated right singular matrix. This updating algorithm is described in Algorithm 2.

Algorithm 2 also checks if $u$ is linearly dependent on the current left singular vectors numerically. If $p < ϵ_{SVD}$ , then we consider that it is linearly dependent. Thus, we set $p = 0$ in $Q$ , i.e., Line 10 of Algorithm 2. Then we only update the first $r_{k - 1}$ components of the singular matrices in Line 14 of Algorithm 2. Otherwise, we follow the update form in Eq. (28) as in Line 16 in Algorithm 2.

Line 18-20 in Algorithm 2 checks if the updated singular value has a small value. If it does, we neglect that particular singular value and corresponding component in left and right singular matrices. It is because a small singular value causes a large error in left and right singular matrices [18].

Although the orthogonality of the updated left singular matrix, $W_{k}$ , must be guaranteed in infinite precision by the product of two orthogonal matrices in Line 14 or 16 of Algorithm 2

, it is not guaranteed in finite precision. Thus, we heuristically check the orthogonality in Lines 21-24 of Algorithm

2 by checking the inner product of the first and last columns of

Φ_{k}

. If the orthogonality is not shown, then we orthogonalize them by the QR factorization. Here

ϵ

denotes unit roundoff (e.g., eps in MATLAB).

The spatial basis can be set after $n_{μ} N_{t}$ incremental steps:

Φ_{s} = W_{n_{μ} N_{t}} (:, 1 : n_{s}) .

(29)

If all the time step solutions are taken incrementally and sequentially from $n_{μ}$ different high-fidelity time dependent simulations, then the right singular matrix, $V_{N_{t} n_{μ}} \in R^{N_{t} n_{μ} \times r_{N_{t} n_{μ}}}$ , holds $n_{μ}$ different temporal behavior for each spatial basis vector. For example, $v_{i}$ describes $n_{μ}$ different temporal behavior of $w_{i}$ . As in Section 4.1, $i$ th temporal snapshot matrix $Υ_{i} \in R^{N_{t} \times n_{μ}}$ can be defined as

Υ_{i} \equiv [\begin{matrix} v_{i}^{1} & \dots & v_{i}^{n_{μ}} \end{matrix}],

(30)

where $v_{i}^{k} \equiv V_{N_{t} n_{μ}} (1 + (k - 1) N_{t} : k N_{t}, i)$ for $k \in N (n_{μ})$ . If we take the SVD of $Υ_{i} = Λ_{i} Σ_{i} Ψ_{i}^{T}$ , then the temporal basis for $i$ th spatial basis vector can be set

Φ_{t}^{i} = Λ_{i} (:, 1 : n_{t}) .

(31)

[ $W_{k}$ , $σ_{k}$ , $V_{k}$ ] = initializingIncrementalSVD( $u$ , $ϵ_{SVD}$ , $k$ )
Input: $u$ , $ϵ_{SVD}$ , $k$
Output: $W_{k}$ , $σ_{k}$ , $V_{k}$

1: if

∥ u ∥ > ϵ_{SVD}

then

σ_{k} \leftarrow [\begin{matrix} ∥ u ∥ \end{matrix}]

W_{k} \leftarrow u / σ_{1}

, and

V_{k} \leftarrow [\begin{matrix} 1 \end{matrix}]

3: else

σ_{k} \leftarrow []

W_{k} \leftarrow []

, and

V_{k} \leftarrow []

5: end if

Algorithm 1 Initializing incremental SVD

[ $W_{k}$ , $σ_{k}$ , $V_{k}$ ] = incrementalSVD( $u$ , $ϵ_{SVD}$ , $ϵ_{SV}$ , $W_{k - 1}$ , $σ_{k}$ , $V_{k - 1}$ , $k$ )
Input: $u$ , $ϵ_{SVD}$ , $ϵ_{SV}$ , $W_{k - 1}$ , $σ_{k - 1}$ , $V_{k - 1}$ , $k$
Output: $W_{k}$ , $σ_{k}$ , $V_{k}$

1: if

r_{k - 1} = 0

r_{k - 1} = r_{max}

then

2: [

W_{k}

σ_{k}

V_{k}

] = initializingIncrementalSVD(

u

ϵ_{SVD}

k

), i.e., apply Algorithm 1

3: return

4: end if

ℓ \leftarrow W_{k - 1}^{T} u

p \leftarrow \sqrt{u^{T} u - ℓ^{T} ℓ}

j \leftarrow (u - W_{k - 1} ℓ) / p

Q \leftarrow [\begin{matrix} diag (s) & ℓ 0 & p \end{matrix}]

9: if

p < ϵ_{SVD}

then

10:

Q_{end, end} \leftarrow 0

11: end if

12:

#

SVD update

13: if

p < ϵ_{SVD}

then

14: , , and

V_{k} \leftarrow [\begin{matrix} V_{k} & 0 0 & 1 \end{matrix}] {¯ ¯¯ ¯ V}_{(:, 1 : r_{k - 1})}

15: else

16:

W_{k} \leftarrow [\begin{matrix} W_{k - 1} & j \end{matrix}] ¯ ¯¯¯¯ ¯ W

, , and

V_{k} \leftarrow [\begin{matrix} V_{k} & 0 0 & 1 \end{matrix}] ¯ ¯¯ ¯ V

17: end if

#

Neglect small singular values: truncation

18: if

{σ_{k}}_{(r_{k})} < ϵ_{SV}

then

19:

σ_{k} \leftarrow {σ_{k}}_{(1 : r_{k} - 1)}

W_{k} \leftarrow {W_{k}}_{(:, 1 : r_{k} - 1)}

V_{k} \leftarrow {V_{k}}_{(:, 1 : r_{k} - 1)}

20: end if

#

Orthogonalize if necessary

21: if

{W_{k}}_{(:, 1)}^{T} {W_{k}}_{(:, end)} > min {ϵ_{SVD}, ϵ \cdot N_{s}}

then

22:

[¯ ¯¯ ¯ Q, R] \leftarrow Q R {W_{k}}

23:

W_{k} \leftarrow ¯ ¯¯ ¯ Q

24: end if

Algorithm 2 Incremental SVD,

W_{- 1} = []

4.3 Space–time reduced basis in block structure

Forming the space–time basis in Eq. (14) through the the Kronecker product in Eq. (22) requires $N_{s} N_{t} n_{s} n_{t}$ multiplications. This is troublesome, not only because it is computationally costly, but also it requires too much memory. Fortunately, a block structure of the space–time basis in (14) is available:

Φ_{s t} = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} Φ_{s} D_{1}^{1} & \dots & \dots & \dots & Φ_{s} D_{1}^{n_{t}} ⋮ & ⋱ & ⋮ & \cdot^{\cdot^{\cdot}} & ⋮ ⋮ & \dots & Φ_{s} D_{k}^{j} & \dots & ⋮ ⋮ & \cdot^{\cdot^{\cdot}} & ⋮ & ⋱ & ⋮ Φ_{s} D_{N_{t}}^{1} & \dots & \dots & \dots & Φ_{s} D_{N_{t}}^{n_{t}} \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ \in R^{N_{s} N_{t} \times n_{s} n_{t}},

(32)

where $k$ th time step temporal basis matrix, $D_{k}^{j} \in R^{n_{s} \times n_{s}}$ , is defined as

D_{k}^{j} \equiv diag ([\begin{matrix} ϕ_{1 j, k}^{t} & \dots & ϕ_{n_{s} j, k}^{t} \end{matrix}]),

(33)

where $ϕ_{i j, k}^{t} \in R$ denotes a $k$ th element of $ϕ_{i j}^{t}$ . Thanks to this block structure, the space–time reduced order operators, such as ${^A}^{st}$ , ${^f}^{st}$ , and ${^u}_{0}^{st}$ can be formed without explicitly forming $Φ_{st}$ . For example, the reduced space–time system matrix, ${^A}^{st}$ can be computed, using the block structures, as

(34)

where ( $j^{'}$ , $j$ )th block matrix, ${^A}_{(j^{'}, j)}^{st} (μ) \in R^{n_{s} \times n_{s}}$ , $j^{'}$ , $j \in N (n_{t})$ is defined as

{^A}_{(j^{'}, j)}^{st} (μ) = N_{t} \sum k = 1 (D_{k}^{j^{'}} D_{k}^{j} - Δ t^{(k)} D_{k}^{j^{'}}^A (μ) D_{k}^{j}) - N_{t} - 1 \sum k = 1 D_{k + 1}^{j^{'}} D_{k}^{j} .

(35)

Note that the computations of $D_{k}^{j^{'}} D_{k}^{j}$ and $D_{k + 1}^{j^{'}} D_{k}^{j}$ are trivial because they are diagonal matrix-products whose individual product requires $n_{s}$ scalar products. Additionally, $^A : R^{n_{μ}} \to R^{n_{s} \times n_{s}}$ , is a reduced order system operator that is used for the spatial ROMs, e.g., see Eq. (10). This can be pre-computed. It implies that the construction of ${^A}^{st} (μ)$ requires computational cost that is a bit larger than the one for the spatial ROM system matrix. The additional cost $O (n_{s}^{2} n_{t}^{2} N_{t})$ to the spatial ROM system matrix construction is required.

Similarly, the reduced space–time input vector, ${^f}^{st} \in R^{n_{s} n_{t}}$ , can be computed as

{^f}^{st} = Φ_{st}^{T} f^{st} = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ \begin{matrix} ⋮ {^f}_{(j)}^{st} ⋮ \end{matrix} ⎤ ⎥ ⎥ ⎥ ⎥ ⎦,

(36)

where the $j$ th block vector, ${^f}_{(j)}^{st} \in R^{n_{s}}$ , $j \in N (n_{t})$ , is given as

{^f}_{(j)}^{st} (μ) = N_{t} \sum k = 1 Δ t^{(k)} D_{k}^{j}^B (μ) f^{(k)} .

(37)

Note that $^B : R^{n_{μ}} \to R^{n_{s} \times N_{i}}$ is used for the spatial ROMs, e.g., see Eq. (10). Also, $^B (μ) f^{(k)}$ needs to be computed in the spatial ROMs. These can be pre-computed. Other operations are related to the row-wise scaling with diagonal term, $D_{k}^{j}$ , whose computational cost is $O (n_{s} n_{t} N_{t})$ . If $f^{(k)}$ is constant throughout the whole time steps, i.e., $f^{(k)} = f$ , then Eq. (37) can be further reduced to

{^f}_{(j)}^{st} (μ) = (N_{t} \sum k = 1 Δ t^{(k)} D_{k}^{j})^B (μ) f,

(38)

where you can compute the summation term first, then multiply the diagonal term with the precomputed term,

Space-time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems

Authors

Efficient space-time reduced order model for linear dynamical systems in Python using less than 120 lines of code

Mixed interpolatory and inference non-intrusive reduced order modeling with application to pollutants dispersion

A reduction methodology using free-free component eigenmodes and Arnoldi enrichment

Space-time reduced-order modeling for uncertainty quantification

Pressio: Enabling projection-based model reduction for large-scale nonlinear dynamical systems

Spectral Methods - Part 2: A comparative study of reduced order models for moisture transfer diffusive problems

Maximum Entropy Snapshot Sampling for Reduced Basis Generation

1 Introduction

1.1 Notations

2 Linear dynamical systems

3 Reduced order models

3.1 Spatial reduced order models

3.2 Space–time reduced order models

4 Basis generation

4.1 Proper orthogonal decomposition

4.2 Incremental space–time reduced basis

4.3 Space–time reduced basis in block structure

Space-time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems

Authors

Related Research

Efficient space-time reduced order model for linear dynamical systems in Python using less than 120 lines of code

Mixed interpolatory and inference non-intrusive reduced order modeling with application to pollutants dispersion

A reduction methodology using free-free component eigenmodes and Arnoldi enrichment

Space-time reduced-order modeling for uncertainty quantification

Pressio: Enabling projection-based model reduction for large-scale nonlinear dynamical systems

Spectral Methods - Part 2: A comparative study of reduced order models for moisture transfer diffusive problems

Maximum Entropy Snapshot Sampling for Reduced Basis Generation

This week in AI

1 Introduction

1.1 Notations

2 Linear dynamical systems

3 Reduced order models

3.1 Spatial reduced order models

3.2 Space–time reduced order models

4 Basis generation

4.1 Proper orthogonal decomposition

4.2 Incremental space–time reduced basis

4.3 Space–time reduced basis in block structure