High-Fidelity Digital Twin Data Models by Randomized Dynamic Mode Decomposition and Deep Learning with Applications in Fluid Dynamics

In plain words

A detailed simulation of a fluid can produce a huge table of numbers over time. This paper builds a much smaller stand-in for that simulation, which the author calls a digital twin data model. The stand-in is meant to behave like the original, so it can answer questions for times and settings where running the full simulation is slow.

The stand-in has two parts. The first finds a short list of spatial patterns that capture most of the motion in the simulation output, and for each pattern it records how it grows or shrinks and oscillates. The second part is a small neural network that learns how the strength of each pattern changes over time, so the stand-in can be run without the full simulation.

The method is tested on three cases of a simple one-dimensional flow equation with sharp fronts, at three settings of a viscosity-related parameter. The author reports small errors against the simulation and short run times, though no full-order run time is given for comparison. The evidence is in-sample and comes from these three cases only.

TL;DR

The paper uses a randomized singular value decomposition inside dynamic mode decomposition (DMD) to get a low-rank model, and it picks the rank by minimising a time-averaged error §4, alg. 1, eq. 29. A cascade-forward neural network, called an NLARX estimator, then models the time series of the DMD amplitudes, and the pair forms the digital twin data model §5, eq. 31–32. The method is tested on the 1D viscous Burgers equation at Re = 10², 10³ and 10⁴ §2. The reported errors are in-sample, measured on the same data the models were built from eq. 33 §6 tab. 1–4.

Problem

Full-order simulations of shock-like flows are costly, and the paper wants a reduced model that keeps the dynamics but is cheap to evaluate §1. Proper orthogonal decomposition has truncation errors that limit its use on complex systems §1. DMD links each spatial mode to one frequency and growth rate, which is useful for flow analysis §1, §3, but choosing which modes to keep is a known difficulty §4. The Moore–Penrose approach to the DMD operator is said not to be feasible for high-dimensional data §3. The paper does not require Galerkin projection of the governing equations [abstract].

Method

Data. The snapshots come from a finite-difference solver for the viscous Burgers equation with a discontinuous initial state §2, eq. 1–2. The snapshot matrix V has N_t + 1 columns eq. 10. The shifted matrices V0 and V1 hold the first and last N_t columns eq. 12.

DMD operator. The paper writes the DMD operator as the minimiser of ‖V1 − V0 S‖ eq. 18. The full-size solution S = V0⁺ V1 eq. 19 is set aside as infeasible in high dimension §3. The SVD route instead works from V0 = U Σ Wᴴ eq. 20. Each mode φ_j has a Ritz value λ_j, a growth rate σ_j and a frequency ω_j eq. 24. The reduced model keeps N_DMD leading modes, with amplitudes a_j(t_i) = ã_j λ_j^(i−1) eq. 26.

Randomized SVD (k-RSVD). A random test matrix with r = min(N_t, 2k) columns samples the range of V0. The sample is orthonormalised, V0 is projected onto the resulting basis, and an economy SVD of the small matrix is computed. The left vectors are lifted back to the full space alg. 2. fig. 2 gives a schematic of this routine.

Randomized DMD (Algorithm 1). For each rank k from 2 to N_t − 1, the algorithm runs k-RSVD, solves the reduced problem of eq. 18, takes the eigen-decomposition of the reduced operator, forms the dynamic modes Φ = U X, projects them onto the first snapshot to get amplitudes, and rebuilds the snapshots through a Vandermonde factor alg. 1 steps 1–6, eq. 30.

# Randomized SVD of rank k of the snapshot matrix V0  [alg. 2]
# Input:  V0 (Nx, Nt); integer k with 2 <= k < Nt
# Output: U (Nx, k); Sigma (k, k); W (Nt, k), as the paper writes them
def k_rsvd(V0, k):
    # Random test matrix with r = min(Nt, 2k) columns  [alg. 2 step 1]
    r = min(Nt, 2 * k)
    test = random_matrix(Nt, r)   # distribution and seed not stated
    # Sample the range of V0  [alg. 2 step 2]
    sample = V0 @ test   # (Nx, r)
    # Orthonormalise the sample by Gram-Schmidt  [alg. 2 step 3]
    basis = gram_schmidt(sample)   # (Nx, r)
    # Project the snapshots onto the small basis  [alg. 2 step 4]
    small = basis.conj().T @ V0   # (r, Nt)
    # Economy SVD of the small matrix  [alg. 2 step 5]
    # Truncation to k is not shown in the steps
    T, Sigma, W = svd_economy(small)
    # Lift the left vectors back to the full space  [alg. 2 step 6]
    U = basis @ T
    return U, Sigma, W

# Randomized DMD for one candidate rank k  [alg. 1 steps 2-6]
# Input:  V0, V1 (Nx, Nt) shifted snapshot matrices  [eq. 12]
# Output: eigenvalues, modes Phi (Nx, k), amplitudes, rebuilt snapshots
def rdmd_model(V0, V1, k):
    # Low-rank factors of V0  [alg. 1 step 2]
    U, Sigma, W = k_rsvd(V0, k)
    # Reduced operator S minimising ||V1 - V0 S||  [alg. 1 step 3] [eq. 18]
    # The paper gives no formula for S in terms of U, Sigma, W
    S = reduced_operator(U, Sigma, W, V1)
    # Eigen-decomposition S X = X Lambda  [alg. 1 step 4]
    eigenvalues, X = eig(S)
    # Dynamic modes Phi = U X  [alg. 1 step 4]
    Phi = U @ X
    # Amplitudes from projecting the modes onto the first snapshot
    # [alg. 1 step 5]; the projection rule is not stated
    amplitudes = project(Phi, V0[:, 0])
    # Rebuild the snapshots with a Vandermonde factor  [alg. 1 step 6] [eq. 30]
    # The paper's Vandermonde runs to power N-2, but N is not defined
    vandermonde = powers_of(eigenvalues)
    V_dmd = Phi @ diag(amplitudes) @ vandermonde
    return eigenvalues, Phi, amplitudes, V_dmd

# Pick the rank: smallest time-averaged error over the snapshots  [eq. 27-29]
# Input:  V0, V1, snapshots (Nx, Nt + 1), Nt
def select_rank(V0, V1, snapshots, Nt):
    # Try every rank from 2 to Nt - 1  [alg. 1 step 1]
    best_error = infinity
    best_model = None
    for k in range(2, Nt):
        model = rdmd_model(V0, V1, k)
        # L2 error averaged over the Nt sample times  [eq. 27] [eq. 28]
        error = mean_l2_error(snapshots, model.V_dmd)
        if error < best_error:
            best_error = error
            best_model = model
    # Rank and model with the smallest error  [eq. 29]
    return best_model

Online stage (NLARX). Each DMD amplitude series a(t) is the input to a nonlinear autoregressive estimator, and its output â(t) is the temporal coefficient of the digital twin data model §5, eq. 31–32. The estimator f is a cascade-forward network with 10 hidden units §5, fig. 3. The regressors are n_a past outputs and n_b past inputs with a pure delay n_k eq. 32. The three orders are chosen per coefficient group (see Hyperparameters) tab. 2–4.

# Online stage: NLARX estimate of one temporal coefficient  [§5] [eq. 32]
# a_past: DMD amplitude series (input); a_hat_past: estimator outputs
def nlarx_step(a_hat_past, a_past, n_a, n_b, n_k, f):
    # Past outputs a_hat(t-1) .. a_hat(t-n_a)  [eq. 32]
    # Delayed inputs a(t-n_k) .. a(t-n_k-n_b+1)  [eq. 32]
    # Lag ordering inside the input vector is not stated
    regressors = concatenate(
        a_hat_past[:n_a],
        a_past[n_k:n_k + n_b],
    )
    # Nonlinear network f plus its linear block  [§5]
    return f(regressors)

# Training cost as printed  [eq. 33]
# Reading the norm of a scalar series as its absolute value (our reading)
# Input:  a (Nt,) amplitudes; a_hat (Nt,) estimator output
def nlarx_cost(a, a_hat, Nt):
    total = 0.0
    for i in range(Nt):
        total = total + abs(a[i] - a_hat[i])
    return total / (2 * Nt)

Evaluation metrics. The paper defines an error and a correlation for the DMD model eq. 27, eq. 34 and the same two for the digital twin data model eq. 35, eq. 36. Both correlation formulas are discussed in Pitfalls.

Hyperparameters

Results

Author-reported, in-sample. Setup: the snapshot data of §2 for each Reynolds number (301 snapshots, N_x = 101). The text names no external baseline. The only comparison is between two rank-search methods, GA-SQP and SA, which the author says took similar effort §6. The text reports no full-order solver run time.

Reduced DMD model tab. 1

Digital twin data model, NLARX per coefficient group. Each row of tab. 2–tab. 4 lists the coefficient indices j sharing one (n_a, n_b, n_k) and carries one metric. A blank cell means the paper gives no value in that row.

Timing §6, fig. 13: offline stage at most about 2 s; online stage 8 to 17 s, longer for the highest Reynolds number. Hardware is not stated.

Qualitative fig. 7–9, fig. 10–12, fig. 14–16: the leading modes and their growth rates and frequencies are plotted; the digital twin solutions are plotted for each case; the first temporal coefficient is shown against its estimator output for each case.

Paper vs code

No official code.

The paper has no code link, no data availability statement and no project page. The arXiv metadata has no comment field with a link. The DOI in the front matter is the journal page, not a code repository. GitHub searches for "randomized dynamic mode decomposition", "digital twin data model dmd" and "digital twin dmd burgers" returned repositories owned by other people. A user search for "bistrian" and for "diana bistrian" found no account that could be verified as the author's. Nothing was cloned, and no code claims are made.

Pitfalls

• Snapshot count. Section 2 gives N_t + 1 = 301 snapshots §2, but the Table 1 caption says 300 data snapshots were processed tab. 1. The DMD matrices use N_t = 300 columns eq. 12, and the error average runs over N_t samples eq. 28. The text does not say which count the experiments used.
• Data source is not clear. Section 2 says the data come from a finite-difference solver §2. The paper also gives an exact Cole–Hopf solution eq. 9, fig. 1. The text does not say the exact solution was used for anything. The grid spacing and time stepping of the solver are not stated.
• The "shock" initial state. The paper calls these cases shock waves §2, but the initial state has u_L = 0.1 below u_R = 0.5 eq. 2. Our own note, not the paper's: for the inviscid equation, a left state below the right state gives a rarefaction, not a shock. Reproduce the solver output rather than assume the wave type.
• Frequency formula. Equation 24 writes ω_j = arg(|λ_j|)/Δt. The argument of a modulus is not a phase. The intended formula is presumably arg(λ_j) divided by Δt. Growth rate is log|λ_j|/Δt. Use the phase of λ_j.
• The reduced operator is missing. Algorithm 1 step 3 says to solve the minimisation of eq. 18 alg. 1 step 3. The full solution eq. 19 is called infeasible for high-dimensional data §3. The paper never writes S in terms of U, Σ and W, so an implementer has to choose a reduction. This is the most important gap for reproduction.
• Amplitude projection. Algorithm 1 step 5 projects the modes onto the first snapshot, but the projection rule (least squares or another) is not stated.
• Vandermonde size. The Vandermonde matrix in step 6 runs to power N − 2, but N is not defined anywhere in the text alg. 1 step 6.
• Random projection details. The test matrix is drawn with rand() alg. 2 step 1. The distribution and random seed are not stated, so the randomized results may differ between runs. The steps do not show truncation of the r-column SVD back to k, although the output is stated as rank k alg. 2 step 5, output. The step that computes U = QT is called "right singular vectors", but the matrix is used as the left (mode) matrix U (Nx, k) alg. 2 step 6; alg. 1 step 2.
• Rank search. Algorithm 1 loops over every k from 2 to N_t − 1 alg. 1 step 1 and then calls for solving the optimisation of eq. 29. Section 6 says GA-SQP and SA gave similar effort §6, and the conclusions say hybrid simulated annealing was used §7. Which method produced Table 1 is not stated. The claim that the optimal rank is "the unique solution" §6 is not tested in the text.
• The error is not relative. The figure captions call the error "relative" fig. 4–6 captions, but equation 27 divides nothing by ‖u‖. Use the equation as printed, and check the figure scale.
• Correlation formulas. Equations 34 and 36 read as a norm of an elementwise product, but the text says the inner product is Hermitian §6. For real data the intended definition is unclear. Table 1 gives C_DMD = 1.0000 for all three cases, so the metric does not separate the cases, although the text says the method gives "the best correlation coefficient" §6.
• NLARX cost. Equation 33 sums the norms of scalar differences, which is an absolute error, scaled by 1/(2 N_t). A squared-error loss, which many implementations use, would differ. Training details are not stated, so the choice has to be made by the implementer.
• In-sample evaluation. The NLARX input (the DMD amplitudes) and its training target are the same series eq. 32–33. The validation figures show the estimator's output over the trained window §6, fig. 14–16. No held-out time window, initial condition or Reynolds number is tested. The text does not say whether the estimator runs one step ahead or in free simulation.
• Network structure. Figure 3 shows an input of size 3 fig. 3. The tables use n_a + n_b values of 2, 3 and 4 across groups tab. 2–4. The paper does not say whether one network is trained per coefficient or per group.
• Table layout. Tables 2–4 put E_DTM and C_DTM in different rows, and the cell is blank where the row carries the other metric tab. 2–4. Read each value against its own row's coefficient indices.
• Cost claims are unquantified. The paper says the DTM gives reduced cost in CPU time and hardware §1, §7. It gives offline and online times only, on unstated hardware §6. It does not report the time of the full-order solver.
• Journal version. The LaTeX front matter says the paper was originally published in Modelling 3 (2022), no. 3, pages 314–332 [front matter]. The arXiv record is version 1 dated 2026-09-15. The card is based on the arXiv source.
• Naming. The conclusions call the method "adaptive randomized DMD" §7, but Section 4 describes randomized DMD with a rank search only. Nothing in the text defines an adaptive step beyond the search.

Lineage

• Builds on: the author's randomized DMD for non-intrusive reduced-order modelling, Bistrian and Navon (2017) §1, §4, and their earlier DMD vs POD comparison, Bistrian and Navon (2014) §1, §3. These are not cards in this wiki.
• Compared against: no external baseline. The paper compares GA-SQP and simulated annealing for the rank search only §6.
• Superseded by: not known from the content of this wiki.