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, mode_prox=None, remove_bad_bags=False, bag_warning=100, bag_maxfail=200, varpro_opts_dict=None)[source]

Bases: 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.6.

  • 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.

  • mode_prox (function) – Optional proximal operator function to apply to the DMD modes. If use_proj is False, this function is applied at every iteration of the variable projection routine. If use_proj is True, this function is instead applied at the end of the variable projection routine after the modes have been projected back to the space of the full input data.

  • remove_bad_bags (bool) – Whether or not to exclude results from bagging trials that didn’t converge according to the tolerance used for variable projection. Default is False, all trial results are kept regardless of convergence.

  • 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. This parameter becomes active only when remove_bad_bags=True. Use negative arguments for no warning condition.

  • bag_maxfail (int) – Number of consecutive non-converged trials of BOP-DMD at which to terminate the fit. Default is 200. This parameter becomes active only when remove_bad_bags=True. Use negative arguments for 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:
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

plot_eig_uq(eigs_true=None, xlim=None, ylim=None, figsize=None, dpi=None, flip_axes=False, draw_axes=False)[source]

Plot BOP-DMD eigenvalues against 1 and 2 standard deviations.

Parameters:

  • eigs_true (np.ndarray or iterable) – True continuous-time eigenvalues, if known.

  • xlim (iterable) – Desired limits for the x-axis.

  • ylim (iterable) – Desired limits for the y-axis.

  • figsize (iterable) – Width, height in inches.

  • dpi (int) – Figure resolution.

  • flip_axes (bool) – Whether or not to swap the real and imaginary axes on the eigenvalue plot. If True, the real axis will be vertical and the imaginary axis will be horizontal.

  • draw_axes (bool) – Whether or not to draw the real and imaginary axes.

plot_mode_uq(*, x=None, y=None, d=1, modes_shape=None, order='C', cols=4, figsize=None, dpi=None, plot_modes=None, plot_conjugate_pairs=True)[source]

Plot BOP-DMD modes alongside their standard deviations.

Parameters:

  • x (np.ndarray or iterable) – Points along the 1st spatial dimension where data has been collected.

  • y (np.ndarray or iterable) – Points along the 2nd spatial dimension where data has been collected. This parameter is only applicable when the data snapshots are 2-D, which must be indicated with modes_shape.

  • d (int) – Number of delays applied to the data. If d is greater than 1, then each plotted mode will be the average mode taken across all d delays.

  • modes_shape (iterable) – Shape of the modes. If not provided, the shape is assumed to be the flattened space dim of the snapshot data. Provide as width, height dimension.

  • order ({"C", "F", "A"}) – Read the elements of snapshots using this index order, and place the elements into the reshaped array using this index order. It has to be the same used to store the snapshots.

  • cols (int) – Number of columns to use for the subplot grid.

  • figsize (iterable) – Width, height in inches.

  • dpi (int) – Figure resolution.

  • plot_modes (int or iterable) – Number of leading modes to plot, or the indices of the modes to plot. If None, then all available modes are plotted. Note that if this parameter is given as a list of indices, it will override the plot_complex_pair parameter.

  • plot_conjugate_pairs (bool) – Whether or not to omit one of the modes that correspond with a complex conjugate pair of eigenvalues.

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, mode_prox, remove_bad_bags, 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: 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.

  • mode_prox (function) – Optional proximal operator function to apply to the DMD modes. If use_proj is False, this function is applied at every iteration of the variable projection routine. If use_proj is True, this function is instead applied at the end of the variable projection routine after the modes have been projected back to the space of the full input data.

  • remove_bad_bags (bool) – Whether or not to exclude results from bagging trials that didn’t converge according to the tolerance used for variable projection. Default is False, all trial results are kept regardless of convergence.

  • 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. This parameter becomes active only when remove_bad_bags=True. Use negative arguments for no warning condition.

  • bag_maxfail (int) – Number of consecutive non-converged trials of BOP-DMD at which to terminate the fit. Default is 200. This parameter becomes active only when remove_bad_bags=True. Use negative arguments for 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:
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:
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