Merge pull request #110 from decargroup/feature/78-implement-random-f…

…ourier-feature-rff-lifting-functions

Feature/78 implement random fourier feature rff lifting functions

.github/workflows/test-package.yml

@@ -10,7 +10,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ['3.7', '3.8', '3.9', '3.10']
+        python-version: ['3.8', '3.9', '3.10']
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{matrix.python-version}}

doc/examples.rst

@@ -60,3 +60,13 @@ initial velocity). Latin hypercube sampling is used to generate 100 centers.
 
 .. plot:: ../examples/6_example_rbf_pendulum.py
    :include-source:
+
+Random Fourier features on a Duffing oscillator
+-----------------------------------------------
+
+This example shows how random Fourier features (and randomly binned features)
+can be used as lifting functions to identify Duffing oscillator dynamics.
+For more details on how these features are generated, see [RR07]_.
+
+.. plot:: ../examples/7_example_rff_duffing.py
+   :include-source:

doc/pykoop.rst

@@ -137,6 +137,7 @@ namespace for convenience.
    pykoop.BilinearInputLiftingFn
    pykoop.ConstantLiftingFn
    pykoop.DelayLiftingFn
+   pykoop.KernelApproxLiftingFn
    pykoop.PolynomialLiftingFn
    pykoop.RbfLiftingFn
    pykoop.SkLearnLiftingFn
@@ -161,6 +162,19 @@ convenience.
    pykoop.EdmdMeta
 
 
+Kernel approximation methods
+============================
+
+The following classes are used to generate random feature maps from kernels
+for kernel approximation lifting functions (i.e., random Fourier feature
+lifting functions).
+
+.. autosummary::
+   :toctree: _autosummary/
+
+   pykoop.RandomBinningKernelApprox
+   pykoop.RandomFourierKernelApprox
+
 Radial basis function centers
 =============================
 
@@ -203,6 +217,7 @@ convenience.
    :toctree: _autosummary/
 
    pykoop.AnglePreprocessor
+   pykoop.example_data_duffing
    pykoop.example_data_msd
    pykoop.example_data_pendulum
    pykoop.example_data_vdp
@@ -269,6 +284,7 @@ ones have been imported into the ``pykoop`` namespace.
    pykoop.Centers
    pykoop.EpisodeDependentLiftingFn
    pykoop.EpisodeIndependentLiftingFn
+   pykoop.KernelApproximation
    pykoop.KoopmanLiftingFn
    pykoop.KoopmanRegressor
    pykoop.dynamic_models.ContinuousDynamicModel

doc/references.rst

@@ -30,3 +30,10 @@ References
    identification of vehicle dynamics using Koopman operator," in 22nd
    International Conference on Process Control, 2019.
    https://doi.org/10.1109%2Fpc.2019.8815104
+.. [RR07] Ali Rahimi and Benjamin Recht. "Random Features for Large-Scale
+   Kernel Machines." in Proc. 20th Int. Conf. Neural. Inf. Process. Syst.,
+   2007. https://proceedings.neurips.cc/paper/2007/hash/013a006f03dbc5392effeb8f18fda755-Abstract.html
+.. [WHDKR16] Matthew O. Williams, Maziar S. Hemati, Scott T. M. Dawson, Ioannis
+   G. Kevrekidis, and Clarence W. Rowley. "Extending data-driven Koopman
+   analysis to actuated systems." IFAC-PapersOnLine, vol. 49, no. 18, pp.
+   704-709, 2016. https://doi.org/10.1016/j.ifacol.2016.10.248

examples/7_example_rff_duffing.py

@@ -0,0 +1,96 @@
+"""Example of random Fourier features on a Duffing oscillator."""
+import numpy as np
+import scipy.stats
+from matplotlib import pyplot as plt
+
+import pykoop
+
+plt.rc('lines', linewidth=2)
+plt.rc('axes', grid=True)
+plt.rc('grid', linestyle='--')
+
+
+def example_random_fourier_features() -> None:
+    """Demonstrate random Fourier featuers on a Duffing oscillator."""
+    # Get example Duffing oscillator data
+    eg = pykoop.example_data_duffing()
+
+    # Create RFF pipeline
+    kp_rff = pykoop.KoopmanPipeline(
+        lifting_functions=[(
+            'rff',
+            pykoop.KernelApproxLiftingFn(
+                kernel_approx=pykoop.RandomFourierKernelApprox(
+                    n_components=100,
+                    random_state=1234,
+                )),
+        )],
+        regressor=pykoop.Edmd(),
+    )
+
+    # Create random binning pipeline for comparison
+    kp_bin = pykoop.KoopmanPipeline(
+        lifting_functions=[(
+            'bin',
+            pykoop.KernelApproxLiftingFn(
+                kernel_approx=pykoop.RandomBinningKernelApprox(
+                    n_components=10,
+                    random_state=1234,
+                )),
+        )],
+        regressor=pykoop.Edmd(),
+    )
+
+    for label, kp in [
+        ('Random Fourier features', kp_rff),
+        ('Randomly binned features', kp_bin),
+    ]:
+        # Fit the pipeline
+        kp.fit(
+            eg['X_train'],
+            n_inputs=eg['n_inputs'],
+            episode_feature=eg['episode_feature'],
+        )
+
+        # Get training and validation episodes
+        ep_train = np.unique(eg['X_train'][:, 0])
+        ep_valid = np.unique(eg['X_valid'][:, 0])
+
+        # Predict new trajectories
+        X_pred = kp.predict_trajectory(
+            eg['x0_valid'],
+            eg['u_valid'],
+        )
+
+        # Plot validation trajectories
+        fig, ax = plt.subplots()
+        tab20 = plt.cm.tab20(np.arange(0, 1, 0.05))
+        color = [
+            (tab20[1], tab20[0]),
+            (tab20[3], tab20[2]),
+            (tab20[5], tab20[4]),
+        ]
+        for (i, ep) in enumerate(ep_valid):
+            idx = eg['X_valid'][:, 0] == ep
+            cax_valid = ax.plot(
+                eg['X_valid'][idx, 1],
+                eg['X_valid'][idx, 2],
+                color=color[i][0],
+                label=f'True traj. {i}',
+            )
+            cax_pred = ax.plot(
+                X_pred[idx, 1],
+                X_pred[idx, 2],
+                label=f'Pred. traj. {i}',
+            )
+        # Set legend
+        ax.legend(loc='lower center', ncol=3)
+        # Set labels
+        ax.set_title(f'{label}: true and predicted trajectories')
+        ax.set_xlabel('$x_1[k]$')
+        ax.set_ylabel('$x_2[k]$')
+
+
+if __name__ == '__main__':
+    example_random_fourier_features()
+    plt.show()

notebooks/6_example_rbf_pendulum.ipynb

@@ -35,7 +35,7 @@
    "id": "61af60a4-e5b5-4d11-93db-38642e4e657f",
    "metadata": {},
    "source": [
-    "Load example data from the library. `eg` is a `dict` containing training data, validation ata, and a few related parameters."
+    "Load example data from the library. `eg` is a `dict` containing training data, validation data, and a few related parameters."
    ]
   },
   {
@@ -312,7 +312,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.7"
+   "version": "3.10.8"
   }
  },
  "nbformat": 4,

pykoop/__init__.py

@@ -1,5 +1,20 @@
 """Koopman operator identification library in Python."""
 
+from .centers import (
+    Centers,
+    ClusterCenters,
+    DataCenters,
+    GaussianMixtureRandomCenters,
+    GaussianRandomCenters,
+    GridCenters,
+    QmcCenters,
+    UniformRandomCenters,
+)
+from .kernel_approximation import (
+    KernelApproximation,
+    RandomBinningKernelApprox,
+    RandomFourierKernelApprox,
+)
 from .koopman_pipeline import (
     EpisodeDependentLiftingFn,
     EpisodeIndependentLiftingFn,
@@ -21,25 +36,17 @@
     DelayLiftingFn,
     PolynomialLiftingFn,
     RbfLiftingFn,
+    KernelApproxLiftingFn,
     SkLearnLiftingFn,
 )
 from .regressors import Dmd, Dmdc, Edmd, EdmdMeta
 from .tsvd import Tsvd
 from .util import (
     AnglePreprocessor,
+    example_data_duffing,
     example_data_msd,
-    example_data_vdp,
     example_data_pendulum,
+    example_data_vdp,
     random_input,
     random_state,
 )
-from .centers import (
-    Centers,
-    GridCenters,
-    UniformRandomCenters,
-    GaussianRandomCenters,
-    GaussianMixtureRandomCenters,
-    QmcCenters,
-    ClusterCenters,
-    DataCenters,
-)

pykoop/dynamic_models.py

@@ -257,6 +257,43 @@ def f(self, t, x, u):
         return x_dot
 
 
+class DuffingOscillator(ContinuousDynamicModel):
+    r"""Duffing oscillator model.
+
+    Equation is ``\ddot{x} + \delta \dot{x} + \beta x + \alpha x^3 = u(t)``
+    where usually ``u(t) = a \cos(\omega t)``.
+    """
+
+    def __init__(
+        self,
+        alpha: float = 1,
+        beta: float = -1,
+        delta: float = 0.1,
+    ) -> None:
+        """Instantiate :class:`DuffingOscillator`.
+
+        Parameters
+        ----------
+        alpha : float
+            Coefficient of cubic term.
+        beta : float
+            Coefficient of linear term.
+        delta : float
+            Coefficient of first derivative.
+        """
+        self.alpha = alpha
+        self.beta = beta
+        self.delta = delta
+
+    def f(self, t: float, x: np.ndarray, u: np.ndarray):
+        # noqa: D102
+        x_dot = np.array([
+            x[1],
+            u - self.delta * x[1] - self.beta * x[0] - self.alpha * x[0]**3
+        ])
+        return x_dot
+
+
 class DiscreteVanDerPol(DiscreteDynamicModel):
     """Van der Pol oscillator.