BOPDMD: Optimized DMD and Bagging, Optimized DMD¶

Derived module from dmdbase.py for Optimized DMD and Bagging, Optimized DMD (BOP-DMD).

References: - Travis Askham and J. Nathan Kutz. Variable projection methods for an optimized dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems, 17, 2018. - Diya Sashidhar and J. Nathan Kutz. Bagging, optimized dynamic mode decomposition (bop-dmd) for robust, stable forecasting with spatial and temporal uncertainty-quantification. 2021. arXiv:2107.10878.

class BOPDMD(svd_rank=0, compute_A=False, use_proj=True, init_alpha=None, proj_basis=None, num_trials=0, trial_size=0.6, eig_sort='auto', eig_constraints=None, bag_warning=100, bag_maxfail=-1, varpro_opts_dict=None)[source]

Bases: pydmd.dmdbase.DMDBase

Bagging, Optimized Dynamic Mode Decomposition.

Parameters

svd_rank (int or float) – The rank for the truncation; If 0, the method computes the optimal rank and uses it for truncation; if positive integer, the method uses the argument for the truncation; if float between 0 and 1, the rank is the number of the biggest singular values that are needed to reach the ‘energy’ specified by svd_rank; if -1, the method does not compute truncation.
compute_A (bool) – Flag that determines whether or not to compute the full Koopman operator A. Default is False, do not compute the full operator. Note that the full operator is potentially prohibitively expensive to compute.
use_proj (bool) – Flag that determines the type of computation to perform. If True, fit input data projected onto the first svd_rank POD modes or columns of proj_basis if provided. If False, fit the full input data. Default is True, fit projected data.
init_alpha (numpy.ndarray) – Initial guess for the continuous-time DMD eigenvalues. If not provided, one is computed via a trapezoidal rule approximation. Default is None (alpha not provided).
proj_basis (numpy.ndarray) – Orthogonal basis for projection, where each column of proj_basis contains a basis mode. If not provided, POD modes are used. Default is None (basis not provided).
num_trials (int) – Number of BOP-DMD trials to perform. If num_trials is a positive integer, num_trials BOP-DMD trials are performed. Otherwise, standard optimized dmd is performed. Default is 0.
trial_size (int or float) – Size of the randomly selected subset of observations to use for each trial of bagged optimized dmd (BOP-DMD). If trial_size is a positive integer, trial_size many observations will be used per trial. If trial_size is a float between 0 and 1, int(trial_size * m) many observations will be used per trial, where m denotes the total number of data points observed. Note that any other type of input for trial_size will yield an error. Default is 0.2.
eig_sort ({"real", "imag", "abs", "auto"}) – Method used to sort eigenvalues (and modes accordingly) when performing BOP-DMD. Eigenvalues will be sorted by real part and then by imaginary part to break ties if eig_sort=”real”, by imaginary part and then by real part to break ties if eig_sort=”imag”, or by magnitude if eig_sort=”abs”. If eig_sort=”auto”, one of the previously-mentioned sorting methods is chosen depending on eigenvalue variance. Default is “auto”.
eig_constraints (set(str) or function) – Set containing desired DMD operator eigenvalue constraints. Currently available constraints are “stable”, which constrains eigenvalues to the left half of the complex plane, “imag”, which constrains eigenvalues to the imaginary axis, and “conjugate_pairs”, which enforces that eigenvalues are always present with their complex conjugate. Note that constraints may be combined if valid. May alternatively be a custom eigenvalue constraint function that will be applied to the computed eigenvalues at each step of the variable projection routine.
bag_warning (int) – Number of consecutive non-converged trials of BOP-DMD at which to produce a warning message for the user. Default is 100. Use arguments less than or equal to zero for no warning condition.
bag_maxfail (int) – Number of consecutive non-converged trials of BOP-DMD at which to terminate the fit. Default is -1, no stopping condition.
varpro_opts_dict (dict) – Dictionary containing the desired parameter values for variable projection. The following parameters may be specified: init_lambda, maxlam, lamup, use_levmarq, maxiter, tol, eps_stall, use_fulljac, verbose. Default values will be used for any parameters not specified in varpro_opts_dict. See BOPDMDOperator documentation for default values and descriptions for each parameter.

property A

Get the full Koopman operator A.

Returns: the full Koopman operator A.
Return type: numpy.ndarray

_check_eig_constraints(eig_constraints)[source]

Function that verifies that…

if the given eig_constraints is a function, the function is able to

take a (n,) numpy.ndarray and return a (n,) numpy.ndarray - if the given eig_constraints is a set, it does not contain an unsupported constraint class, and does not contain an invalid combination of eigenvalue constraints

_initialize_alpha()[source]

Uses projected trapezoidal rule to approximate the eigenvalues of A in: z’ = Az.

The computed eigenvalues will serve as our initial guess for alpha.

Returns: Approximated eigenvalues of the matrix A.
Return type: numpy.ndarray

property amplitudes_std

Get the amplitudes standard deviation.

Returns: amplitudes standard deviation.
Return type: numpy.ndarray

property atilde

Get the reduced Koopman operator A, called Atilde.

Returns: the reduced Koopman operator A.
Return type: numpy.ndarray

property compute_A

Returns: flag that determines whether to compute the full operator A.
Return type: bool

property dynamics

Get the time evolution of each mode.

Returns: matrix that contains all the time evolution, stored by row.
Return type: numpy.ndarray

property eig_constraints

Get the eigenvalue constraints.

Returns: eigenvalue constraints.
Return type: set(str)

property eigenvalues_std

Get the eigenvalues standard deviation.

Returns: eigenvalues standard deviation.
Return type: numpy.ndarray

fit(X, t)[source]

Compute the Optimized Dynamic Mode Decomposition.

Parameters

X (numpy.ndarray or iterable) – the input snapshots.
t (numpy.ndarray or iterable) – the input time vector.

forecast(t)[source]

Predict the output X given the input time t using the fitted DMD model. If model has been fitted using multiple ensembles, an average prediction and its variance is returned.

Parameters: t (numpy.ndarray or iterable) – the input time vector.
Returns: system prediction at times given by vector t.
Return type: numpy.ndarray or numpy.ndarray, numpy.ndarray

property init_alpha

Returns: initial guess used for the continuous-time DMD eigenvalues.
Return type: numpy.ndarray

property modes_std

Get the modes standard deviation.

Returns: modes standard deviation.
Return type: numpy.ndarray

property num_trials

Returns: the number of BOP-DMD trials to perform.
Return type: int

print_varpro_opts()[source]: Prints a formatted information string that displays all chosen variable projection parameter values.

property proj_basis

Returns: the projection basis used, with modes stored by column.
Return type: numpy.ndarray

property svd_rank

Returns: the rank used for the svd truncation.
Return type: int or float

property time

Get the vector that contains the time points of the fitted snapshots.

Returns: the vector that contains the original time points.
Return type: numpy.ndarray

property trial_size

Returns: size of the data subsets used during each BOP-DMD trial.
Return type: int or float

property use_proj

Returns: flag that determines whether to fit projected or full data.
Return type: bool

class BOPDMDOperator(compute_A, use_proj, init_alpha, proj_basis, num_trials, trial_size, eig_sort, eig_constraints, bag_warning, bag_maxfail, init_lambda=1.0, maxlam=52, lamup=2.0, use_levmarq=True, maxiter=30, tol=1e-06, eps_stall=1e-12, use_fulljac=True, verbose=False)[source]

Bases: pydmd.dmdoperator.DMDOperator

BOP-DMD operator.

Parameters

compute_A (bool) – Flag that determines whether or not to compute the full Koopman operator A.
use_proj (bool) – Flag that determines the type of computation to perform. If True, fit input data projected onto the first svd_rank POD modes or columns of proj_basis if provided. If False, fit the full input data.
init_alpha (numpy.ndarray) – Initial guess for the continuous-time DMD eigenvalues.
proj_basis (numpy.ndarray) – Orthogonal basis for projection, where each column of proj_basis contains a basis mode.
num_trials (int) – Number of BOP-DMD trials to perform. If num_trials is a positive integer, num_trials BOP-DMD trials are performed. Otherwise, standard optimized dmd is performed.
trial_size (int or float) – Size of the randomly selected subset of observations to use for each trial of bagged optimized dmd (BOP-DMD). If trial_size is a positive integer, trial_size many observations will be used per trial. If trial_size is a float between 0 and 1, int(trial_size * m) many observations will be used per trial, where m denotes the total number of data points observed. Note that any other type of input for trial_size will yield an error.
eig_sort ({"real", "imag", "abs", "auto"}) – Method used to sort eigenvalues (and modes accordingly) when performing BOP-DMD. Eigenvalues will be sorted by real part and then by imaginary part to break ties if eig_sort=”real”, by imaginary part and then by real part to break ties if eig_sort=”imag”, or by magnitude if eig_sort=”abs”. If eig_sort=”auto”, one of the previously-mentioned sorting methods is chosen depending on eigenvalue variance.
eig_constraints (set(str) or function) – Set containing desired DMD operator eigenvalue constraints. Currently available constraints are “stable”, which constrains eigenvalues to the left half of the complex plane, “imag”, which constrains eigenvalues to the imaginary axis, and “conjugate_pairs”, which enforces that eigenvalues are always present with their complex conjugate. Note that constraints may be combined if valid. May alternatively be a custom eigenvalue constraint function that will be applied to the computed eigenvalues at each step of the variable projection routine.
bag_warning (int) – Number of consecutive non-converged trials of BOP-DMD at which to produce a warning message for the user. Default is 100. Use arguments less than or equal to zero for no warning condition.
bag_maxfail (int) – Number of consecutive non-converged trials of BOP-DMD at which to terminate the fit. Default is -1, no stopping condition.
init_lambda (float) – Initial value used for the regularization parameter in the Levenberg method. Default is 1.0. Note: Larger lambda values make the method more like gradient descent.
maxlam (int) – Maximum number of of steps used in the inner Levenberg loop, i.e. the number of times you increase lambda before quitting. Default is 52.
lamup (float) – The factor by which you increase lambda when searching for an appropriate step. Default is 2.0.
use_levmarq (bool) – Flag that determines whether you use the Levenberg algorithm or the Levenberg-Marquardt algorithm. Default is True, use Levenberg-Marquardt.
maxiter (int) – The maximum number of outer loop iterations to use before quitting. Default is 30.
tol (float) –
The tolerance for the relative error in the residual. i.e. the program will terminate if

norm(y-Phi(alpha)*b,’fro’)/norm(y,’fro’) < tol

is achieved. Default is 1e-6.
eps_stall (float) –
The tolerance for detecting a stall. i.e. if

error(iter-1)-error(iter) < eps_stall*err(iter-1)

the program halts. Default is 1e-12.
use_fulljac (bool) – Flag that determines whether or not to use the full expression for the Jacobian or Kaufman’s approximation. Default is True, use full expression.
verbose (bool) – Flag that determines whether or not to print warning messages that arise during the variable projection routine, and whether or not to print information regarding the method’s iterative progress. Default is False, don’t print information.

property A

Get the full Koopman operator A.

Returns: the full Koopman operator A.
Return type: numpy.ndarray

_argsort_eigenvalues(eigs)[source]

Helper function that computes and returns the indices that sort the given array of eigenvalues according to the operator’s eig_sort attribute. Sets eig_sort according to eigs if not already done so.

Parameters: eigs (numpy.ndarray) – array of eigenvalues to sort
Returns: array of indices that sort the given eigenvalues
Return type: numpy.ndarray

_bag(H, trial_size)[source]

Given a 2D array of data X, where each row contains a data snapshot, randomly sub-selects and returns data snapshots while preserving the original snapshot order. Note that if trial_size is a positive integer, trial_size many observations will be used per trial. If trial_size is a float between 0 and 1, int(trial_size * m) many observations will be used per trial, where m denotes the total number of snapshots in X. The indices of the sub-selected snapshots are also returned.

Parameters

H (numpy.ndarray) – Full data matrix to be sub-selected from.
trial_size (int or float) – Size of the sub-selection from H.

Returns

Matrix of sub-selected data snapshots, stored in each row, and a vector of each snapshots’s row index location in H.

Return type

numpy.ndarray, numpy.ndarray

static _compute_irank_svd(X, tolrank)[source]

Helper function that computes and returns the SVD of X with a rank truncation of irank, which denotes the number of singular values of X greater than tolrank * s1, where s1 is the largest singular value of the matrix X.

Parameters

X (numpy.ndarray) – Matrix to decompose.
tolrank (float) – Determines the rank of the returned SVD.

Returns

irank truncated SVD of X.

Return type

numpy.ndarray, numpy.ndarray, numpy.ndarray

_exp_function(alpha, t)[source]

Matrix of exponentials.

Parameters

alpha (numpy.ndarray) – Vector of time scalings in the exponent.
t (numpy.ndarray) – Vector of time values.

Returns

Matrix A such that A[i, j] = exp(t_i * alpha_j).

Return type

numpy.ndarray

_exp_function_deriv(alpha, t, i)[source]

Derivatives of the matrix of exponentials.

Parameters

alpha (numpy.ndarray) – Vector of time scalings in the exponent.
t (numpy.ndarray) – Vector of time values.
i (int) – Index in alpha of the derivative variable.

Returns

Derivatives of Phi(alpha, t) with respect to alpha[i].

Return type

scipy.sparse.csr_matrix

_push_eigenvalues(eigenvalues)[source]

Helper function that constrains the given eigenvalues according to the arguments found in self._eig_constraints. If no constraints were given, this function simply returns the given eigenvalues. Applies the provided eig_constraints function if a function was provided instead of a set of constraints.

Parameters: eigenvalues (numpy.ndarray) – Vector of original eigenvalues.
Returns: Vector of constrained eigenvalues.
Return type: numpy.ndarray

_single_trial_compute_operator(H, t, init_alpha)[source]: Helper function that computes the standard optimized dmd operator. Returns the resulting DMD modes, eigenvalues, amplitudes, reduced system matrix, full system matrix, and whether or not convergence of the variable projection routine was reached.

_variable_projection(H, t, init_alpha, Phi, dPhi)[source]

Variable projection routine for multivariate data. Attempts to fit the columns of H as linear combinations of the columns of Phi(alpha,t) such that H = Phi(alpha,t)B. Note that M denotes the number of data samples, N denotes the number of columns of Phi, IS denotes the number of functions to fit, and IA denotes the length of the alpha vector.

Parameters

H (numpy.ndarray) – (M, IS) matrix of data.
t (numpy.ndarray) – (M,) vector of sample times.
init_alpha (numpy.ndarray) – initial guess for alpha.
Phi (function) – (M, N) matrix-valued function Phi(alpha,t).
dPhi (function) – (M, N) matrix-valued function dPhi(alpha,t,i) that contains the derivatives of Phi wrt the ith component of alpha.

Returns

Tuple of two numpy arrays and a boolean representing: 1. (N, IS) best-fit matrix B. 2. (N,) best-fit vector alpha. 3. Flag indicating whether or not convergence was reached.

Return type

Tuple[numpy.ndarray, numpy.ndarray, bool]

References: - Extensions and Uses of the Variable Projection Algorith for Solving Nonlinear Least Squares Problems by G. H. Golub and R. J. LeVeque ARO Report 79-3, Proceedings of the 1979 Army Numerical Analsysis and Computers Conference. - Variable projection for nonlinear least squares problems. Computational Optimization and Applications 54.3 (2013): 579-593 by Dianne P. O’Leary and Bert W. Rust.

_varpro_opts_warn()[source]: Checks the validity of the parameter values in _varpro_opts. Throws an error if any parameter value has an invalid type and generates a warning if any value lies outside of the recommended range.

property amplitudes_std

Get the amplitudes standard deviation.

Returns: amplitudes standard deviation.
Return type: numpy.ndarray

compute_operator(H, t)[source]

Compute the low-rank and the full BOP-DMD operators.

Parameters

H (numpy.ndarray) – Matrix of data to fit.
t (numpy.ndarray) – Vector of sample times.

Returns

The BOP-DMD amplitudes.

Return type

numpy.ndarray

property eigenvalues_std

Get the eigenvalues standard deviation.

Returns: eigenvalues standard deviation.
Return type: numpy.ndarray

property modes_std

Get the modes standard deviation.

Returns: modes standard deviation.
Return type: numpy.ndarray

property varpro_opts

Get the variable projection options.

Returns: the variable projection options.
Return type: tuple