Merge pull request PyDMD#568 from ClimeTrend/fit-econ

Economy fit method for the BOPDMD class

Author: mtezzele

pydmd/bopdmd.py

@@ -1351,17 +1351,31 @@ def _check_eig_constraints(self, eig_constraints):
                     msg = "Invalid eigenvalue constraint combination provided."
                     raise ValueError(msg)
 
-    def _initialize_alpha(self):
+    def _initialize_alpha(self, s=None, V=None):
         """
         Uses projected trapezoidal rule to approximate the eigenvalues of A in
             z' = Az.
         The computed eigenvalues will serve as our initial guess for alpha.
-
+        If the singular values and right singular vectors are provided, they
+        will be used to compute the projection. Otherwise, the snapshot data
+        will be projected onto the projection basis.
+
+        :param s: the singular values of the snapshot matrix. If provided, they
+            will be used along with V to compute the projection.
+        :type s: numpy.ndarray
+        :param V: the right singular vectors of the snapshot matrix. If
+            provided, they will be used along with s to compute the projection.
+        :type V: numpy.ndarray
         :return: Approximated eigenvalues of the matrix A.
         :rtype: numpy.ndarray
         """
-        # Project the snapshot data onto the projection basis.
-        ux = self._proj_basis.conj().T.dot(self.snapshots)
+        if s is not None and V is not None:
+            # If the singular values and right singular vectors are provided,
+            # use them to compute the projection.
+            ux = np.diag(s).dot(V)
+        else:
+            # Project the snapshot data onto the projection basis.
+            ux = self._proj_basis.conj().T.dot(self.snapshots)
         ux1 = ux[:, :-1]
         ux2 = ux[:, 1:]
 
@@ -1461,6 +1475,105 @@ def fit(self, X, t):
 
         return self
 
+    def fit_econ(self, s, V, t):
+        """
+        Compute the Optimized Dynamic Mode Decomposition using the
+        SVD results of the snapshot matrix.
+        Use this method if you have pre-computed the SVD and the
+        original snapshot matrix is not avaialble or too large to
+        store in memory.
+        The left singular vectors (U) must be specified as the projection
+        basis when initializing the BOPDMD object.
+
+        :param s: the singular values.
+        :type s: array_like
+        :param V: matrix of right singular vectors, where each row is a
+            different singular vector.
+        :type V: numpy.ndarray
+        :param t: the input time vector.
+        :type t: numpy.ndarray
+        """
+        self._reset()
+        self._time = np.asarray(t).squeeze()
+        s = np.asarray(s).squeeze()
+
+        if self._proj_basis is None or not self._use_proj:
+            msg = """
+            proj_basis must be provided when using fit_econ,
+            and use_proj must be set to True.
+            """
+            raise ValueError(msg)
+        if (
+            not isinstance(self._proj_basis, np.ndarray)
+            or self._proj_basis.ndim != 2
+            or self._proj_basis.shape[1] != self._svd_rank
+        ):
+            msg = (
+                "proj_basis must be a 2D np.ndarray with "
+                f"{self._svd_rank} columns."
+            )
+            raise ValueError(msg)
+
+        if self._time.ndim > 1:
+            msg = "Input time vector t must be one-dimensional."
+            raise ValueError(msg)
+
+        if s.ndim != 1 or len(s) != self._svd_rank:
+            msg = f"s must be one-dimensional and of length {self._svd_rank}."
+            raise ValueError(msg)
+
+        if (
+            not isinstance(V, np.ndarray)
+            or V.ndim != 2
+            or V.shape[0] != self._svd_rank
+            or V.shape[1] != len(self._time)
+        ):
+            msg = (
+                "V must be a 2D numpy.ndarray with shape "
+                f"({self._svd_rank}, {len(self._time)})."
+            )
+            raise ValueError(msg)
+
+        # Set/check the initial guess for the continuous-time DMD eigenvalues.
+        if self._init_alpha is None:
+            self._init_alpha = self._initialize_alpha(s=s, V=V)
+        elif (
+            not isinstance(self._init_alpha, np.ndarray)
+            or self._init_alpha.ndim > 1
+            or len(self._init_alpha) != self._svd_rank
+        ):
+            msg = "init_alpha must be a 1D np.ndarray with {} entries."
+            raise ValueError(msg.format(self._svd_rank))
+
+        # Build the BOP-DMD operator now that the initial alpha and
+        # the projection basis have been defined.
+        self._Atilde = BOPDMDOperator(
+            self._compute_A,
+            self._use_proj,
+            self._init_alpha,
+            self._proj_basis,
+            self._num_trials,
+            self._trial_size,
+            self._eig_sort,
+            self._eig_constraints,
+            self._mode_prox,
+            self._remove_bad_bags,
+            self._bag_warning,
+            self._bag_maxfail,
+            **self._varpro_opts_dict,
+        )
+
+        # Define the snapshots that will be used for fitting.
+        snp = np.diag(s).dot(V)
+
+        # Fit the data.
+        self._b = self.operator.compute_operator(snp.T, self._time)
+
+        # Store the singular values
+        self._singular_values = s
+
+        return self
+
     def forecast(self, t):
         """
         Predict the output X given the input time t using the fitted DMD model.

pydmd/plotter.py

@@ -650,7 +650,10 @@ def plot_summary(
         sizes of all other eigenvalues are then scaled according to eigenvalue
         significance.
     :type max_eig_ms: int
-    :param max_sval_plot: Maximum number of singular values to plot.
+    :param max_sval_plot: Maximum number of singular values to plot. If the DMD
+        instance was fitted via the `fit_econ` method, the maximum number of
+        singular values that can be plotted is determined by the length of the
+        array of singular values passed to `fit_econ`.
     :type max_sval_plot: int
     :param title_fontsize: Fontsize used for subplot titles.
     :type title_fontsize: int
@@ -781,7 +784,14 @@ def plot_summary(
     else:
         # Use input data matrix to compute singular values.
         snp = dmd.snapshots
-    s = np.linalg.svd(snp, full_matrices=False, compute_uv=False)
+    if snp is not None:
+        s = np.linalg.svd(snp, full_matrices=False, compute_uv=False)
+    else:
+        try:
+            s = dmd._singular_values
+        except AttributeError as e:
+            msg = "Snapshots and singular values are not available."
+            raise ValueError(msg) from e
     # Compute the percent of data variance captured by each singular value.
     s_var = s * (100 / np.sum(s))
     s_var = s_var[:max_sval_plot]

tests/test_bopdmd.py

@@ -641,3 +641,22 @@ def test_getter_errors():
 
     with raises(ValueError):
         _ = bopdmd.amplitudes_std
+
+
+def test_fit_econ():
+    """
+    Test that the fit_econ method gives the same results as the fit method.
+    """
+    svd_rank = 2
+    bopdmd = BOPDMD(svd_rank=svd_rank, use_proj=True)
+    bopdmd.fit(Z, t)
+
+    U, s, V = np.linalg.svd(Z)
+    U = U[:, :svd_rank]
+    s = s[:svd_rank]
+    V = V[:svd_rank, :]
+    bopdmd_econ = BOPDMD(svd_rank=svd_rank, use_proj=True, proj_basis=U)
+    bopdmd_econ.fit_econ(s, V, t)
+
+    np.testing.assert_allclose(bopdmd_econ.eigs, bopdmd.eigs)
+    np.testing.assert_allclose(bopdmd_econ.modes, bopdmd.modes)

tutorials/user-manual1/user-manual-bopdmd.ipynb

@@ -95,6 +95,25 @@
     flip_continuous_axes=True,  # Rotate the continuous-time eig plot.
 )
 
+# Alternatively, pre-compute the SVD and use the `fit_econ` function.
+# This is useful when the SVD is already computed and the snapshots matrix
+# X is not available or is too large to store in memory.
+U, s, V = np.linalg.svd(X, full_matrices=False)
+U = U[:, :11]
+s = s[:11]
+V = V[:11, :]
+
+bopdmd_econ = BOPDMD(svd_rank=11, num_trials=0, use_proj=True, proj_basis=U)
+bopdmd_econ.fit_econ(s, V, t)
+
+plot_summary(
+    bopdmd_econ,
+    figsize=(12, 6),  # Figure size.
+    index_modes=(0, 1, 3),  # Indices of the modes to plot.
+    snapshots_shape=(ny, nx),  # Shape of the modes.
+    order="F",  # Order to use when reshaping the modes.
+    flip_continuous_axes=True,  # Rotate the continuous-time eig plot.
+)
 
 # In[4]: