remove clipping function

change to rank adaptation problem for sw-edmd
modify weighted cost func for new rank values

Author: EthanJamesLew

autokoopman/estimator/koopman.py

@@ -81,14 +81,21 @@ def swdmdc(X, Xp, U, r, Js, W):
     # so the objective isn't numerically unstable
     sf = (1.0 / n_snap)
 
+    # check that rank is less than the number of observations
+    if r > n_obs:
+        warnings.warn("Rank must be less than the number of observations. Reducing rank to n_obs.")
+        r = n_obs
+    
     # koopman operator
-    K = cp.Variable((n_obs, n_obs + n_inps))
+    # Variables for the low-rank approximation
+    K_U = cp.Variable((n_obs, r))
+    K_V = cp.Variable((r, n_obs + n_inps))
 
     # SW-eDMD objective
-    weights_obj = np.vstack([(np.clip(np.abs(J), 0.0, 1.0) @ w) for J, w in zip(Js, W)]).T 
-    P = sf * cp.multiply(weights_obj, Yp.T - K @ Y.T)
+    weights_obj = np.vstack([(np.abs(J) @ w) for J, w in zip(Js, W)]).T 
+    P = sf * cp.multiply(weights_obj, Yp.T - (K_U @ K_V) @ Y.T)
     # add regularization 
-    objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K, "fro"))
+    objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K_U @ K_V, "fro"))
 
     # unconstrained problem
     constraints = None
@@ -100,14 +107,14 @@ def swdmdc(X, Xp, U, r, Js, W):
     try:
         _ = prob.solve(solver=cp.CLARABEL)
         #_ = prob.solve(solver=cp.ECOS)
-        if K.value is None:
+        if K_U.value is None or K_V.value is None:
             raise Exception("SW-eDMD (cvxpy) Optimization failed to converge.")
     except:
         warnings.warn("SW-eDMD (cvxpy) Optimization failed to converge. Switching to unweighted DMDc.")
         return dmdc(X, Xp, U, r)
 
     # get the transformation
-    Atilde = K.value 
+    Atilde = K_U.value @ K_V.value 
     return Atilde[:, :state_size], Atilde[:, state_size:]
 
 

notebooks/weighted-cost-func.ipynb

@@ -16,7 +16,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "dfdb47e2-0ea0-4bf8-8279-8500ff3cf21f",
    "metadata": {},
    "outputs": [],
@@ -36,7 +36,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "id": "291d3409-1c8c-44cb-8380-44f08019b57d",
    "metadata": {},
    "outputs": [],
@@ -45,7 +45,7 @@
     "import autokoopman.benchmark.fhn as fhn\n",
     "fhn = fhn.FitzHughNagumo()\n",
     "training_data = fhn.solve_ivps(\n",
-    "    initial_states=np.random.uniform(low=-2.0, high=2.0, size=(1, 2)),\n",
+    "    initial_states=np.random.uniform(low=-2.0, high=2.0, size=(10, 2)),\n",
     "    tspan=[0.0, 6.0],\n",
     "    sampling_period=0.1\n",
     ")\n",
@@ -68,7 +68,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "id": "e2d42e41-46c2-467c-9ce7-9bd6a7c509a1",
    "metadata": {},
    "outputs": [],
@@ -92,6 +92,7 @@
     "    w = np.ones(traj.states.shape)\n",
     "    w[:, -3:] = 0.0\n",
     "    w[:, :2] = 1.0\n",
+    "    w[:, 0] = 1.0\n",
     "    weights.append(w)\n",
     "\n",
     "    # weight garbage trajectory to zero\n",
@@ -103,6 +104,25 @@
     "#weights = {idx: w for idx, w in enumerate(weights)}"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "ddd76415-6d19-4a38-a2b0-84eb48d0fdfc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from autokoopman.observable import ReweightedRFFObservable\n",
+    "import autokoopman.observable as akobs\n",
+    "\n",
+    "X, WX = list(zip(*list((trajectories[i], w) for i, w in enumerate(weights))))\n",
+    "X, WX = np.vstack(X), np.vstack(WX)\n",
+    "X, WX = np.tile(X, (3, 1)), np.tile(WX, (3, 1))\n",
+    "idxs = np.random.permutation(np.arange(len(X)))\n",
+    "Y, WY = X[idxs], WX[idxs]\n",
+    "\n",
+    "reweight_obs =  akobs.IdentityObservable() | akobs.ReweightedRFFObservable(5, 40, 1.0, X, Y, WX, WY, 'rff')"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "a706f212-36cd-4203-b209-cb7c5ce4ad94",
@@ -122,7 +142,19 @@
    "execution_count": null,
    "id": "98510aa7-3416-4181-a493-00500be53f61",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "Tuning GridSearchTuner:   0%|              | 0/40 [00:00<?, ?it/s]/home/elew/AutoKoopman/notebooks/../autokoopman/estimator/koopman.py:113: UserWarning: SW-eDMD (cvxpy) Optimization failed to converge. Switching to unweighted DMDc.\n",
+      "  warnings.warn(\"SW-eDMD (cvxpy) Optimization failed to converge. Switching to unweighted DMDc.\")\n",
+      "Tuning GridSearchTuner:  12%|▊     | 5/40 [00:23<03:07,  5.37s/it]/home/elew/anaconda3/envs/autokoopman/lib/python3.12/site-packages/numpy/linalg/linalg.py:2582: RuntimeWarning: overflow encountered in multiply\n",
+      "  s = (x.conj() * x).real\n",
+      "Tuning GridSearchTuner:  30%|█▌   | 12/40 [00:52<02:04,  4.44s/it]"
+     ]
+    }
+   ],
    "source": [
     "# learn model from weighted data\n",
     "experiment_results = auto_koopman(\n",
@@ -136,8 +168,8 @@
     "    n_obs=40,              # maximum number of observables to try\n",
     "    max_opt_iter=200,       # maximum number of optimization iterations\n",
     "    grid_param_slices=5,   # for grid search, number of slices for each parameter\n",
-    "    n_splits=None,         # k-folds validation for tuning, helps stabilize the scoring\n",
-    "    rank=(40, 41, 1)       # NOTE: don't rank tune (for now)--SW-eDMD already optimizes for this\n",
+    "    n_splits=5,            # k-folds validation for tuning, helps stabilize the scoring\n",
+    "    rank=(1, 41, 5)       # rank (SW-eDMD now uses rank adaptation)\n",
     ")"
    ]
   },
@@ -290,7 +322,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.0"
+   "version": "3.12.3"
   }
  },
  "nbformat": 4,