Source code for pydmd.spdmd

"""Derived module from dmdbase.py for sparsity-promoting DMD."""

import numpy as np
from numpy.linalg import solve
from scipy.sparse import csc_matrix as sparse
from scipy.sparse import hstack as sphstack
from scipy.sparse import vstack as spvstack
from scipy.sparse.linalg import spsolve

from .dmd import DMD


def soft_thresholding_operator(v, k):
    """
    Soft-thresholding operator as defined in 10.1063/1.4863670.

    :param np.ndarray v: The vector on which we apply the operator.
    :param float k: The threshold.
    :return np.ndarray: The result of the application of the soft-tresholding
        operator on ´v´.
    """
    return np.multiply(
        np.multiply(np.divide(1 - k, np.abs(v)), v), np.abs(v) > k
    )


[docs]class SpDMD(DMD): """ Sparsity-Promoting Dynamic Mode Decomposition. Promotes solutions having an high number of amplitudes set to zero (i.e. *sparse solutions*). Reference: 10.1063/1.4863670 :param svd_rank: the rank for the truncation; If 0, the method computes the optimal rank and uses it for truncation; if positive interger, 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. :type svd_rank: int or float :param int tlsq_rank: rank truncation computing Total Least Square. Default is 0, that means TLSQ is not applied. :param bool exact: flag to compute either exact DMD or projected DMD. Default is True. :param opt: argument to control the computation of DMD modes amplitudes. See :class:`DMDBase`. Default is False. :type opt: bool or int :param rescale_mode: Scale Atilde as shown in 10.1016/j.jneumeth.2015.10.010 (section 2.4) before computing its eigendecomposition. None means no rescaling, 'auto' means automatic rescaling using singular values, otherwise the scaling factors. :type rescale_mode: {'auto'} or None or numpy.ndarray :param bool forward_backward: If True, the low-rank operator is computed like in fbDMD (reference: https://arxiv.org/abs/1507.02264). Default is False. :param sorted_eigs: Sort eigenvalues (and modes/dynamics accordingly) by magnitude if `sorted_eigs='abs'`, by real part (and then by imaginary part to break ties) if `sorted_eigs='real'`. Default: False. :type sorted_eigs: {'real', 'abs'} or False :param float abs_tolerance: Controls the convergence of ADMM. See :func:`_loop_condition` for more details. :param float rel_tolerance: Controls the convergence of ADMM. See :func:`_loop_condition` for more details. :param int max_iterations: The maximum number of iterations performed by ADMM, after that the algorithm is stopped. :param float rho: Controls the convergence of ADMM. For a reference on the optimal value for `rho` see 10.1109/TAC.2014.2354892 or 10.3182/20120914-2-US-4030.00038. :param float gamma: Controls the level of "promotion" assigned to sparse solution. Increasing `gamma` will result in an higher number of zero-amplitudes. :param bool verbose: If `False`, the information provided by SpDMD (like the number of iterations performed by ADMM) are not shown. :param bool enforce_zero: If `True` the DMD amplitudes which should be set to zero according to the solution of ADMM are manually set to 0 (since we solve a sparse linear system to find the optimal vector of DMD amplitudes very small terms may survive in some cases). :param release_memory: If `True` the intermediate matrices computed by the algorithm are deleted after the termination of a call to :func:`fit`. """ def __init__( self, svd_rank=0, tlsq_rank=0, exact=True, opt=False, rescale_mode=None, forward_backward=False, sorted_eigs=False, abs_tolerance=1.0e-6, rel_tolerance=1.0e-4, max_iterations=10000, rho=1, gamma=10, verbose=True, enforce_zero=True, release_memory=True, zero_absolute_tolerance=1.0e-12, ): super().__init__( svd_rank=svd_rank, tlsq_rank=tlsq_rank, exact=exact, opt=opt, rescale_mode=rescale_mode, forward_backward=forward_backward, sorted_eigs=sorted_eigs, ) self.rho = rho self.gamma = gamma self._max_iterations = max_iterations self._abs_tol = abs_tolerance self._rel_tol = rel_tolerance self._verbose = verbose self._enforce_zero = enforce_zero self._release_memory = release_memory self._zero_absolute_tolerance = zero_absolute_tolerance self._P = None self._q = None self._Plow = None self._modes_activation_bitmask_proxy = None
[docs] def fit(self, X): """ Compute the Dynamic Modes Decomposition of the input data. :param X: the input snapshots. :type X: numpy.ndarray or iterable """ super().fit(X) P, q = self._optimal_dmd_matrices() self._P = sparse(P) self._q = q # Cholesky factorization of matrix P + (rho/2)*I Prho = P + np.identity(len(self.amplitudes)) * self.rho / 2 self._Plow = np.linalg.cholesky(Prho) # find which amplitudes are to be set to 0 zero_amplitudes = self._find_zero_amplitudes() # compute the (sparse) vector of optimal DMD amplitudes self._b = self._optimal_amplitudes(zero_amplitudes) # re-allocate the Proxy to avoid problems due to the fact that we # re-computed the amplitudes self._allocate_modes_bitmask_proxy() # release memory if self._release_memory: self._P = None self._q = None self._Plow = None return self
[docs] def _update_alpha(self, beta, lmbd): """ Update the vector :math:`\\alpha_k` of DMD amplitudes. :param np.ndarray beta: Current value of :math:`\\beta_k` (vector of non-zero amplitudes). :param np.ndarray lmbd: Current value of :math:`\\lambda_k` (vector of Lagrande multipliers). :return: The updated value :math:`\\alpha_{k+1}`. :rtype: np.ndarray """ uk = beta - lmbd / self.rho return solve( self._Plow.conj().T, solve(self._Plow, self._q + uk * self.rho / 2) )
[docs] def _update_beta(self, alpha, lmbd): """ Update the vector :math:`\\beta` of non-zero amplitudes. :param np.ndarray alpha: Updated value of :math:`\\alpha_{k+1}` (vector of DMD amplitudes). :param np.ndarray lmbd: Current value of :math:`\\lambda_k` (vector of Lagrange multipliers). :return: The updated value :math:`\\beta_{k+1}`. :rtype: np.ndarray """ return soft_thresholding_operator( alpha + lmbd / self.rho, self.gamma / self.rho )
[docs] def _update_lagrangian(self, alpha, beta, lmbd): """ Update the vector :math:`\\lambda` of Lagrange multipliers. :param np.ndarray alpha: Updated value of :math:`\\alpha_{k+1}` (vector of DMD amplitudes). :param np.ndarray beta: Updated value of :math:`\\beta_{k+1}` (vector of non-zero amplitudes). :param np.ndarray lmbd: Current value of :math:`\\lambda_k` (vector of Lagrange multipliers). :return: The updated value :math:`\\lambda_{k+1}`. :rtype: np.ndarray """ return lmbd + (alpha - beta) * self.rho
[docs] def _update(self, beta, lmbd): """ Operate an entire step of ADMM. :param np.ndarray beta: Current value of :math:`\\beta_k` (vector of non-zero amplitudes). :param np.ndarray lmbd: Current value of :math:`\\lambda_k` (vector of Lagrande multipliers). :return: A tuple containing the updated values :math:`\\alpha_{k+1},\\beta_{k+1},\\lambda_{k+1}` (in this order). :rtype: tuple """ a_new = self._update_alpha(beta, lmbd) b_new = self._update_beta(a_new, lmbd) l_new = self._update_lagrangian(a_new, b_new, lmbd) return a_new, b_new, l_new
[docs] def _loop_condition(self, alpha, beta, lmbd, old_beta): """ Check whether ADMM can stop now, or should perform another iteration. :param np.ndarray alpha: Current value of :math:`\\alpha_k` (vector of DMD amplitudes). :param np.ndarray beta: Current value of :math:`\\beta_k` (vector of non-zero amplitudes). :param np.ndarray lmbd: Current value of :math:`\\lambda_k` (vector of Lagrange multipliers). :param np.ndarray old_beta: Old value of :math:`\\beta_{k-1}` (vector of non-zero amplitudes). :return bool: `True` if ADMM can stop now, `False` otherwise. """ primal_residual = np.linalg.norm(alpha - beta) dual_residual = self.rho * np.linalg.norm(beta - old_beta) eps_primal = np.sqrt(len(alpha)) * self._abs_tol + self._rel_tol * max( np.linalg.norm(alpha), np.linalg.norm(beta) ) eps_dual = np.sqrt( len(alpha) ) * self._abs_tol + self._rel_tol * np.linalg.norm(lmbd) return primal_residual < eps_primal and dual_residual < eps_dual
[docs] def _find_zero_amplitudes(self): """ Use ADMM to find which amplitudes (i.e. their position in the DMD amplitudes array) which can be set to zero according to the given parameters. Note that this method does not compute amplitudes, but only which amplitudes are to be set to 0. Optimal amplitudes should be computed separately afterwards (see :func:`_find_sparsity_promoting_amplitudes`). :return np.ndarray: A boolean vector whose `True` items correspond to amplitudes which should be set to 0. """ n_amplitudes = len(self.amplitudes) # initial values of lmbd and beta are all 0 beta0 = np.zeros(n_amplitudes, dtype="complex") lmbd0 = np.zeros(n_amplitudes, dtype="complex") # perform a first step of ADMM alpha, beta, lmbd = self._update(beta0, lmbd0) old_beta = beta0 # count the number of iterations of ADMM i = 0 # at the beginning of each iteration check if ADMM can stop (because of # loop_condition or number of iterations) while ( not self._loop_condition(alpha, beta, lmbd, old_beta) and i < self._max_iterations ): i += 1 old_beta = beta alpha, beta, lmbd = self._update(beta, lmbd) if self._verbose: print("ADMM: {} iterations".format(i)) # zero values in beta are associated with DMD amplitudes which can be # set to 0 return np.abs(old_beta) < self._zero_absolute_tolerance
[docs] def _optimal_amplitudes(self, zero_amplitudes): """ Find the optimal DMD amplitudes with the constraint that the given indexes should be set to 0. :param np.ndarray zero_amplitudes: Boolean vector. :return np.ndarray: Vector of optimal DMD amplitudes. Amplitudes at indexes corresponding to `True` indexes in `zero_amplitudes` are set to 0. """ n_amplitudes = len(self.amplitudes) n_of_zero = np.count_nonzero(zero_amplitudes) # vectors of the canonical base of R^n_amplitudes, from which we # extract only those corresponding to items set to 0 in zero_amplitudes E = np.identity(n_amplitudes)[:, zero_amplitudes] # left hand side of the linear system whose solution is the vector of # optimal DMD amplitudes. KKT = spvstack( [ sphstack([self._P, E], format="csc"), sphstack( [ E.conj().T, sparse((n_of_zero, n_of_zero), dtype="complex"), ], format="csc", ), ], format="csc", ) # right hand side of the linear system rhs = np.concatenate( ( self._q, np.zeros((n_of_zero,)), ) ) opt_amps = spsolve(KKT, rhs)[:n_amplitudes] if self._enforce_zero: opt_amps[zero_amplitudes] = 0 return opt_amps