Initial commit

Author: sichinaga

pydmd/preprocessing/__init__.py

@@ -4,5 +4,6 @@
 
 from .hankel import hankel_preprocessing
 from .pre_post_processing import PrePostProcessingDMD
+from .randomized import randomized_preprocessing
 from .svd_projection import svd_projection_preprocessing
 from .zero_mean import zero_mean_preprocessing

pydmd/preprocessing/randomized.py

@@ -0,0 +1,62 @@
+"""
+Randomized pre-processing.
+"""
+
+from functools import partial
+from typing import Dict, Union
+
+import numpy as np
+
+from pydmd.dmdbase import DMDBase
+from pydmd.preprocessing import PrePostProcessingDMD
+from pydmd.utils import compute_rqb
+
+svd_rank_type = Union[int, float]
+seed_type = Union[None, int]
+
+def randomized_preprocessing(
+    dmd: DMDBase,
+    svd_rank: svd_rank_type,
+    oversampling: int,
+    power_iters: int,
+    test_matrix: np.ndarray = None,
+    seed: seed_type = None,
+):
+    """
+    Randomized QB pre-processing.
+
+    :param dmd: DMD instance to be wrapped.
+    """
+    return PrePostProcessingDMD(
+        dmd,
+        partial(
+            _rand_preprocessing,
+            svd_rank=svd_rank,
+            oversampling=oversampling,
+            power_iters=power_iters,
+            test_matrix=test_matrix,
+            seed=seed,
+        ),
+        _rand_postprocessing,
+    )
+
+def _rand_preprocessing(
+    state: Dict,
+    X: np.ndarray,
+    svd_rank: svd_rank_type,
+    oversampling: int,
+    power_iters: int,
+    test_matrix: np.ndarray,
+    seed: seed_type,
+    **kwargs
+):
+    Q = compute_rqb(
+        X, svd_rank, oversampling, power_iters, test_matrix, seed
+    )[0]
+    state["compression_matrix"] = Q.conj().T
+
+    return (state["compression_matrix"].dot(X),) + tuple(kwargs.values())
+
+
+def _rand_postprocessing(state: Dict, X: np.ndarray) -> np.ndarray:
+    return state["compression_matrix"].conj().T.dot(X)

pydmd/rdmd.py

@@ -10,7 +10,7 @@
 import numpy as np
 
 from .cdmd import CDMD
-from .utils import compute_rank
+from .utils import compute_rank, compute_rqb
 
 
 class RDMD(CDMD):
@@ -21,6 +21,9 @@ class RDMD(CDMD):
         the Randomized QB Decomposition. If not provided, the `svd_rank` and
         `oversampling` parameters will be used to compute the random matrix.
     :type test_matrix: numpy.ndarray
+    :param seed: Seed used to initialize the random generator when computing
+        random test matrices.
+    :type seed: int
     :param oversampling: Number of additional samples (beyond the desired rank)
         to use when computing the random test matrix. Note that values {5,10}
         tend to be sufficient.
@@ -34,6 +37,7 @@ class RDMD(CDMD):
     def __init__(
         self,
         test_matrix=None,
+        seed=None,
         oversampling=10,
         power_iters=2,
         svd_rank=0,
@@ -58,6 +62,7 @@ def __init__(
         self._oversampling = oversampling
         self._power_iters = power_iters
         self._test_matrix = test_matrix
+        self._seed = seed
 
     def _compress_snapshots(self):
         """
@@ -68,28 +73,39 @@ def _compress_snapshots(self):
         :return: the compressed snapshots
         :rtype: numpy.ndarray
         """
-        # Define the random test matrix if not provided.
-        if self._test_matrix is None:
-            m = self.snapshots.shape[-1]
-            r = compute_rank(self.snapshots, self._svd_rank)
-            self._test_matrix = np.random.randn(m, r + self._oversampling)
+        Q, B, Omega = compute_rqb(
+            self.snapshots,
+            self._svd_rank,
+            self._oversampling,
+            self._power_iters,
+            self._test_matrix,
+            self._seed,
+        )
+        self._compression_matrix = Q.conj().T
+        self._test_matrix = Omega
 
-        # Compute sampling matrix.
-        Y = self.snapshots.dot(self._test_matrix)
+        # # Define the random test matrix if not provided.
+        # if self._test_matrix is None:
+        #     m = self.snapshots.shape[-1]
+        #     r = compute_rank(self.snapshots, self._svd_rank)
+        #     self._test_matrix = np.random.randn(m, r + self._oversampling)
 
-        # Perform power iterations.
-        for _ in range(self._power_iters):
-            Q = np.linalg.qr(Y)[0]
-            Z = np.linalg.qr(self.snapshots.conj().T.dot(Q))[0]
-            Y = self.snapshots.dot(Z)
+        # # Compute sampling matrix.
+        # Y = self.snapshots.dot(self._test_matrix)
 
-        # Orthonormalize the sampling matrix.
-        Q = np.linalg.qr(Y)[0]
+        # # Perform power iterations.
+        # for _ in range(self._power_iters):
+        #     Q = np.linalg.qr(Y)[0]
+        #     Z = np.linalg.qr(self.snapshots.conj().T.dot(Q))[0]
+        #     Y = self.snapshots.dot(Z)
 
-        # Project the snapshot matrix onto the smaller space.
-        B = Q.conj().T.dot(self.snapshots)
+        # # Orthonormalize the sampling matrix.
+        # Q = np.linalg.qr(Y)[0]
 
-        # Save the compression matrix.
-        self._compression_matrix = Q.conj().T
+        # # Project the snapshot matrix onto the smaller space.
+        # B = Q.conj().T.dot(self.snapshots)
+
+        # # Save the compression matrix.
+        # self._compression_matrix = Q.conj().T
 
         return B

pydmd/utils.py

@@ -2,16 +2,18 @@
 
 import warnings
 from numbers import Number
-from typing import NamedTuple
+from typing import NamedTuple, Union
 from collections import namedtuple
 import numpy as np
 from numpy.lib.stride_tricks import sliding_window_view
 
 #  Named tuples used in functions.
 #  compute_svd uses "SVD",
-#  compute_tlsq uses "TLSQ".
+#  compute_tlsq uses "TLSQ",
+#  compute_rqb uses "RQB".
 SVD = namedtuple("SVD", ["U", "s", "V"])
 TLSQ = namedtuple("TLSQ", ["X_denoised", "Y_denoised"])
+RQB = namedtuple("RQB", ["Q", "B", "Omega"])
 
 
 def _svht(sigma_svd: np.ndarray, rows: int, cols: int) -> int:
@@ -189,6 +191,44 @@ def compute_svd(
     return SVD(U, s, V)
 
 
+def compute_rqb(
+    X: np.ndarray,
+    svd_rank: Number,
+    oversampling: int = 10,
+    power_iters: int = 2,
+    Omega: np.ndarray = None,
+    seed: Union[None, int] = None,
+) -> NamedTuple(
+    "RQB", [("Q", np.ndarray), ("B", np.ndarray), ("Omega", np.ndarray)]
+):
+    if X.ndim != 2:
+        raise ValueError("Please ensure that input data is a 2D array.")
+
+    # Define the random test matrix if not provided.
+    if Omega is None:
+        m = X.shape[-1]
+        r = compute_rank(X, svd_rank)
+        rng = np.random.default_rng(seed)
+        Omega = rng.standard_normal((m, r + oversampling))
+
+    # Compute sampling matrix.
+    Y = X.dot(Omega)
+
+    # Perform power iterations.
+    for _ in range(power_iters):
+        Q = np.linalg.qr(Y)[0]
+        Z = np.linalg.qr(X.conj().T.dot(Q))[0]
+        Y = X.dot(Z)
+
+    # Orthonormalize the sampling matrix.
+    Q = np.linalg.qr(Y)[0]
+
+    # Project the snapshot matrix onto the smaller space.
+    B = Q.conj().T.dot(X)
+
+    return RQB(Q, B, Omega)
+
+
 def pseudo_hankel_matrix(X: np.ndarray, d: int) -> np.ndarray:
     """
     Arrange the snapshot in the matrix `X` into the (pseudo) Hankel