add regularization because of sparse case (@Abdu-Hekal #88)

EthanJamesLew · Apr 13, 2024 · d216b75 · d216b75
1 parent b08476e
commit d216b75
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 10 deletions.
diff --git a/autokoopman/estimator/koopman.py b/autokoopman/estimator/koopman.py
@@ -55,6 +55,8 @@ def wdmdc(X, Xp, U, r, W):
 def swdmdc(X, Xp, U, r, Js, W):
     """State Weighted Dynamic Mode Decomposition with Control (wDMDC)"""
     import cvxpy as cp
+    import warnings
+
     assert len(W.shape) == 2, "weights must be 2D for snapshot x state"
 
     if U is not None:
@@ -76,7 +78,8 @@ def swdmdc(X, Xp, U, r, Js, W):
     # SW-eDMD objective
     weights_obj = np.vstack([(np.abs(J) @ w)[:n_obs] for J, w in zip(Js, W)]).T 
     P = sf * cp.multiply(weights_obj, Xp[:, :n_obs].T - K @ X[:, :n_obs].T)
-    objective = cp.Minimize(cp.sum_squares(P))
+    # add regularization 
+    objective = cp.Minimize(cp.sum_squares(P) + 1E-4 * 1.0 / (n_obs**2) * cp.norm(K, "fro"))
 
     # unconstrained problem
     constraints = None
@@ -85,8 +88,14 @@ def swdmdc(X, Xp, U, r, Js, W):
     prob = cp.Problem(objective, constraints)
 
     # solve for the SW-eDMD Koopman operator
-    result = prob.solve()
-    # TODO: check the result
+    _ = prob.solve(solver=cp.CLARABEL)
+    #_ = prob.solve(solver=cp.ECOS)
+
+    # backup case
+    if K.value is None:
+        # give a warning about the optimization failure
+        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 

diff --git a/notebooks/weighted-cost-func.ipynb b/notebooks/weighted-cost-func.ipynb
@@ -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=(10, 2)),\n",
+    "    initial_states=np.random.uniform(low=-2.0, high=2.0, size=(1, 2)),\n",
     "    tspan=[0.0, 6.0],\n",
     "    sampling_period=0.1\n",
     ")\n",
@@ -90,7 +90,7 @@
     "    #w = np.sum(traj.abs().states, axis=1)\n",
     "    #w = 1/(traj.abs().states+1.0)\n",
     "    w = np.ones(traj.states.shape)\n",
-    "    w[:, -3:] = 1.0\n",
+    "    w[:, -3:] = 0.0\n",
     "    w[:, :2] = 1.0\n",
     "    weights.append(w)\n",
     "\n",
@@ -136,8 +136,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=(20, 21, 1)       # rank range (start, stop, step) DMD hyperparameter\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",
     ")"
    ]
   },
@@ -168,9 +168,9 @@
     "    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=5,             # k-folds validation for tuning, helps stabilize the scoring\n",
+    "    n_splits=None,             # k-folds validation for tuning, helps stabilize the scoring\n",
     "    lengthscale=(0.1, 1.0),\n",
-    "    rank=(20, 21, 1)       # rank range (start, stop, step) DMD hyperparameter\n",
+    "    rank=(1, 41, 5)       # rank range (start, stop, step) DMD hyperparameter\n",
     ")"
    ]
   },
@@ -290,7 +290,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.9.0"
   }
  },
  "nbformat": 4,