Source-linked AI summary
Variable projection methods for an optimized dynamic mode decomposition
Travis Askham, J. Nathan Kutz
TL;DR
Standard DMD can produce noise-biased eigenvalues, while optimized DMD requires a potentially large nonlinear optimization problem. The paper uses variable projection and exponential data fitting to compute optimized DMD efficiently, including for unevenly spaced sample times. The resulting decomposition overcomes exact-DMD bias issues and can substantially improve robustness to noise and accuracy of identified yearly patterns.
Problem
Standard DMD eigenvalues are biased by sensor noise, while optimized DMD avoids much of this bias but requires a potentially large nonlinear optimization problem.
Method
The paper reformulates DMD as exponential data fitting and applies variable projection to eliminate many unknowns, enabling optimized DMD for arbitrary sample times.
Results
The optimized DMD overcomes exact-DMD bias issues, with sometimes an order of magnitude improvement in error and more accurate yearly and half-yearly wavelengths in sea surface temperature data.
Takeaways & Limitations
Optimized DMD provides a flexible DMD-based analysis that can handle unevenly spaced samples and improve robustness to noise.
Takeaways & Limitations
When eigenvalues are confluent or nearly confluent, exponential fitting becomes rank-deficient or ill-conditioned, making non-diagonalizable dynamics difficult to approximate.
Abstract
from arXiv · showhide
The dynamic mode decomposition (DMD) has become a leading tool for data-driven modeling of dynamical systems, providing a regression framework for fitting linear dynamical models to time-series measurement data. We present a simple algorithm for computing an optimized version of the DMD for data which may be collected at unevenly spaced sample times. By making use of the variable projection method for nonlinear least squares problems, the algorithm is capable of solving the underlying nonlinear optimization problem efficiently. We explore the performance of the algorithm with some numerical examples for synthetic and real data from dynamical systems and find that the resulting decomposition displays less bias in the presence of noise than standard DMD algorithms. Because of the flexibility of the algorithm, we also present some interesting new options for DMD-based analysis.
Key word.
Dynamic mode decomposition is presented as a data-driven framework for modeling dynamical systems through inverse linear systems and variable projection algorithms.
- The paper concerns inverse linear systems and a variable projection algorithm for DMD.
- Dynamic mode decomposition decomposes dynamical-system data into modes associated with oscillation frequencies and growth or decay rates.
- DMD provides dynamical information but does not produce orthogonal modes, unlike proper orthogonal decomposition.
1. Introduction.
The introduction motivates optimized DMD as a response to noise bias in standard DMD and develops efficient variable-projection algorithms that also support uneven sample times.
- Motivation: Sensor noise biases eigenvalues computed by standard DMD algorithms.
- Motivation: Optimized DMD fits all snapshots simultaneously, reducing bias but requiring a potentially large nonlinear optimization problem.
- Contributions: Variable projection eliminates many unknowns from the exponential-fitting optimization underlying optimized DMD.
- Contributions: The proposed algorithms allow sample times that are not equispaced and include an initialization scheme based on standard DMD.
- Evaluation: The paper evaluates optimized DMD on synthetic and real data, including examples with unevenly spaced snapshots.
2. Preliminaries.
The preliminaries establish the notation used throughout the paper, including conventions for matrices and vectors.
- Notation: Matrices are generally denoted by bold capital letters, while vectors are denoted by bold lowercase letters.
- Notation: The Moore-Penrose pseudoinverse is denoted by A†.
2.1. Notation.
The notation defines common matrix and vector operations before reviewing classical variable projection algorithms and their modern extensions.
- Indexing and slicing: Vector entries and matrix entries are indexed with subscripts, while submatrices and column-based vectorizations use MATLAB-style slicing notation.
- Matrix operations: A⊺ denotes transpose, A∗ denotes conjugate transpose, and A† denotes the Moore-Penrose pseudoinverse.
- Matrix operations: The notation A = diag(v) represents a diagonal matrix whose diagonal entries are the entries of v.
- Matrix operations: The Kronecker product is denoted by A ⊗ B.
- Variable projection: The paper reviews classical variable projection algorithms relevant to optimized DMD and notes modern advances in the framework.
2.2. Variable projection.
Variable projection converts separable nonlinear least-squares problems into optimization over nonlinear parameters alone, recovering linear coefficients afterward. Applied to inverse differential equations, this framework supports optimized DMD with arbitrary sample times, while its applicability is constrained by rank and eigenvalue structure.
- Variable projection: For fixed α, variable projection computes the minimizing coefficients as β = Φ(α)†η, reducing the optimization to α alone.After finding α̂, the corresponding coefficients are recovered by β̂ = Φ(α̂)†η.
- Variable projection: Levenberg–Marquardt iteratively updates α using the residual Jacobian and trust-region scaling to seek a nearby local minimizer.The parameter ν_i is adjusted so accepted steps reduce the residual and can approach Gauss–Newton updates.
- Multiple right-hand sides: The multiple-right-hand-side formulation couples several functions through the same nonlinear parameter vector α while solving for coefficient matrix B in blocked form.For a given α, B = Φ†H, and the residual and Jacobian computations can also be blocked.
- Computational considerations: Computational cost depends strongly on efficient SVD, sparse derivative operations, and the least-squares solve for the Jacobian update.Under the stated sparsity and full-rank assumptions, solving for δ costs O(k^2mn), while Jacobian formation is organized into four steps.
- Inverse differential equations: Recasting DMD as inverse differential equations extends the method naturally to arbitrary sample times and allows fitting a lower-dimensional linear system in a higher-dimensional state space.The inverse problem finds A from z(t_i) at m ≥ n sample times, and the framework also applies when n > l.
- Limitations: Nearly confluent eigenvalues make Φ(α) rank-deficient or ill-conditioned, while dynamics from non-diagonalizable systems are difficult to approximate using exponentials alone.A block upper-triangular decomposition is proposed as an alternative direction, but extending the method to rank selection remains future work.
2.3. The DMD.
The DMD fits a linear map between snapshot data and represents that map through eigenvalue–eigenvector pairs. Exact DMD uses pseudoinverse-based least squares and SVD truncation, but noise, rank selection, numerical conditioning, and computational bias motivate optimized and debiased variants.
- Exact DMD: Exact DMD defines the best-fit linear operator as A = YX† for data matrices X and Y.This operator minimizes the Frobenius residual in overdetermined systems and gives the minimum-norm solution in underdetermined systems.
- Exact DMD: DMD is the set of eigenvectors and eigenvalues of the fitted operator, providing modes with associated oscillation frequencies and growth or decay rates.The modes supply dynamical information but are not generally orthogonal, unlike POD modes.
- Rank selection: SVD truncation retains the largest r singular values to avoid fitting dynamics to low-energy modes that may be corrupted by noise.The truncation choice can substantially affect the resulting modes and eigenvalues, and hard-thresholding is only one possible rule.
- Exact DMD: The standard exact-DMD workflow forms X and Y, computes a reduced SVD of X, constructs a reduced operator, and eigendecomposes it.Each reduced eigenpair yields a DMD eigenvalue and mode.
- Bias and debiasing: Sensor noise biases computed DMD eigenvalues, partly because standard algorithms treat snapshot pairs separately and favor the forward-time direction.Debiased alternatives include methods designed to reduce this bias within the standard DMD framework.
- Numerical issues: Square-root nonuniqueness can make reduced-operator construction require an apparently O(2^r) search, and eigenvector ill-conditioning can arise when eigendecomposition is unavailable.The paper advises avoiding numerically difficult Jordan decompositions and notes that alternative decompositions have trade-offs.
- Optimized DMD: Optimized DMD fits all snapshots simultaneously, reducing much of the original DMD bias while requiring a potentially large nonlinear optimization problem.The paper combines this formulation with variable projection to develop efficient algorithms for arbitrary sample times.
3. The optimized DMD.
The optimized DMD fits snapshot data to exponential dynamics by solving a nonlinear least-squares problem, with variable projection eliminating linear unknowns. The section also develops an SVD-based approximation, initialization strategies, and computational improvements for unevenly sampled data.
- Definition: Optimized DMD minimizes the Frobenius residual between the transposed snapshot matrix and an exponential model Φ(α)B over α and B.The optimized eigenvalues are the fitted α values, while the eigenmodes are obtained from normalized columns of the fitted coefficient matrix.
- Limitations: Levenberg–Marquardt is not guaranteed to find the global minimizer, so algorithms 2 and 3 may be viewed as post-processors that improve an initial α through a nearby local minimizer.The paper nevertheless reports significant improvement over other DMD methods in its numerical section.
- Approximate optimized DMD: The approximate optimized DMD replaces X with its optimal rank-r approximation Xr, and algorithm 3 computes the resulting reduced problem.This formulation is useful for denoising or restricting data to low-dimensional structure, while often producing reconstructions comparable to the full problem.
- Approximate optimized DMD: After one rank-r SVD, each variable-projection iteration costs O(r^3m), compared with O(r^2mn) for the original formulation, enabling significant speedups.The initial rank-r SVD costs O(mn min(m, n)) with a standard algorithm or O(r^2(m+n)+rmn) with a randomized algorithm.
- Initialization: Variable projection requires a good initial guess for α, so the paper proposes finite-difference initialization and can use exact DMD eigenvalues after logarithmic scaling for equispaced samples.The solution of the reduced problem can also provide a good initial condition for the full optimized DMD problem.
4. Examples.
Across synthetic noisy-data examples, optimized DMD generally recovers eigenvalues and reconstructions more accurately than standard alternatives, with advantages especially visible at higher noise and limited snapshots. Its performance remains strong for both periodic and growing/decaying systems, although some comparisons favor fbDMD on particular measures.
- Example 1: The reconstructed system-matrix advantage of optimized DMD is strongest at the highest noise levels and fewest snapshots.Exact DMD error eventually levels off at higher noise, whereas the other methods generally improve with more snapshots.
- Example 1: Optimized DMD has lower eigenvalue error than fbDMD and tlsDMD across all noise levels, decreasing faster as snapshots increase.This advantage is more pronounced for eigenvalue recovery than for reconstructed system-matrix error.
- Example 1: The optimized DMD eigenvalue distribution has smaller spread than fbDMD and tlsDMD in the periodic-system confidence-ellipse comparison.Exact DMD remains biased into the left half-plane, while fbDMD and tlsDMD correct that bias.
- Example 1: Optimized DMD improves snapshot reconstruction across noise levels and snapshot counts, while other methods’ errors increase with more snapshots, especially under higher noise.Its reconstruction error is relatively flat as the number of snapshots increases.
- Example 2: For growing-system eigenvalues, optimized DMD is up to an order of magnitude more accurate in some settings and retains an advantage for hidden eigenvalues.For hidden eigenvalues, all methods eventually flatten as the signal decays, but optimized DMD remains advantageous across tested noise levels and snapshot counts.
- Example 2: In the growing-system example, reconstruction errors roughly coincide across methods even though optimized DMD recovers eigenvalues significantly better.The growing modes increasingly dominate later snapshots, limiting the distinction visible in reconstruction error.
- Computational cost: Optimized DMD requires more computation than the other methods, but the increase is modest in the reported runtime comparison.For larger rank r, its computational cost appears to increase like r^2.
- Example 3: For uncertain sample times, fbDMD and optimized DMD perform best in reconstruction, while optimized DMD has lower eigenvalue error and distinct confidence-ellipse spreads.In the confidence ellipses, fbDMD has smaller real-axis spread and optimized DMD smaller imaginary-axis spread.
5. Conclusions and future directions.
The optimized DMD reduces bias and improves recovery of dynamical patterns, including yearly and half-yearly sea-surface-temperature modes. Its flexibility supports uneven sampling and new analysis directions, but it requires nonlinear nonconvex optimization with higher computational cost in larger examples.
- The optimized DMD overcomes bias issues of exact DMD and can improve robustness to noise substantially.In one example, it captured hidden dynamics with sometimes an order-of-magnitude lower error than other DMD methods.
- For sea surface temperature data, yearly and half-yearly wavelength accuracy was comparable with the data, while exact DMD and tlsDMD were an order of magnitude less accurate.
- Computing the optimized DMD requires solving a nonlinear, nonconvex optimization problem, and the experiments do not establish that the global optimum was found.
- Runtime was about six times that of exact DMD for the larger climate example, although for the largest rank it was roughly equal to the omitted SVD projection cost.
- The optimized DMD does not require evenly spaced snapshots, enabling analyses with arbitrary sample times.The paper highlights cross-validation with arbitrary times and possible incoherent sampling for detecting high-frequency signals.
- Future directions include sparsity priors, block Schur decompositions for transient dynamics, and incoherent time sampling for high-frequency signals.