collapsed and enhanced the 2a and 2b notebooks into a single one and linked readme appropriately, sorted out optimization for robustdiff, reran notebooks, re-added tests for robustdiff

Author: pavelkomarov

README.md

@@ -83,8 +83,7 @@ The following heuristic works well for choosing `tvgamma`, where `cutoff_frequen
 
 Much more extensive usage is demonstrated in Jupyter notebooks:
 * Differentiation with different methods: [1_basic_tutorial.ipynb](https://github.com/florisvb/PyNumDiff/blob/master/examples/1_basic_tutorial.ipynb)
-* Parameter Optimization with known ground truth (only for demonstration purpose):  [2a_optimizing_parameters_with_dxdt_known.ipynb](https://github.com/florisvb/PyNumDiff/blob/master/examples/2a_optimizing_parameters_with_dxdt_known.ipynb)
-* Parameter Optimization with unknown ground truth: [2b_optimizing_parameters_with_dxdt_unknown.ipynb](https://github.com/florisvb/PyNumDiff/blob/master/examples/2b_optimizing_parameters_with_dxdt_unknown.ipynb)
+* Parameter Optimization:  [2_optimizing_hyperparameters.ipynb](https://github.com/florisvb/PyNumDiff/blob/master/examples/2_optimizing_hyperparameters.ipynb)
 * Automatic method suggestion: [3_automatic_method_suggestion.ipynb](https://github.com/florisvb/PyNumDiff/blob/master/examples/3_automatic_method_suggestion.ipynb)
 
 ## Repo Structure

examples/1_basic_tutorial.ipynb

@@ -74,7 +74,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "id": "1eba51f1-bf31-4f84-b82b-176461319de6",
    "metadata": {},
    "outputs": [],
@@ -90,9 +90,7 @@
     "\t\t\t(tvrdiff, 'TVRDiff'),\n",
     "\t\t\t(smooth_acceleration, 'SmoothAccelTVR'),\n",
     "\t\t\t(rtsdiff, 'RTSDiff'),\n",
-    "\t\t\t(robustdiff, 'RobustDiff'),\n",
-    "\t\t    (robustdiff, 'RobustDiff2'),\n",
-    "\t\t    (robustdiff, 'RobustDiff3')]\n",
+    "\t\t\t(robustdiff, 'RobustDiff')]\n",
     "sims = [(pi_cruise_control, 'Cruise Control'),\n",
     "\t\t(sine, 'Sum of Sines'),\n",
     "\t\t(triangle, 'Triangles'),\n",

examples/2_optimizing_hyperparameters.ipynb

@@ -276,6 +276,10 @@ def robustdiff(x, dt, order, log_q, log_r, proc_huberM=6, meas_huberM=0):
     norm case, because :math:`c_2` approaches :math:`\\frac{\\sqrt{2}}{M}`, cancelling the :math:`M` multiplying :math:`|\\cdot|`
     and leaving behind :math:`\\sqrt{2}`, the proper normalization.
 
+    Note that :code:`log_q` and :code:`proc_huberM` are coupled, as are :code:`log_r` and :code:`meas_huberM`, via the relation
+    :math:`\\text{Huber}(q^{-1/2}v, M) = q^{-1}\\text{Huber}(v, Mq^{-1/2})`, but they are still independent enough that for
+    the purposes of optimization we cannot collapse them.
+
     :param np.array[float] x: data series to differentiate
     :param float dt: step size
     :param int order: which derivative to stabilize in the constant-derivative model (1=velocity, 2=acceleration, 3=jerk)
@@ -337,7 +341,7 @@ def huber_const(M): # from https://jmlr.org/papers/volume14/aravkin13a/aravkin13
                 else huber_const(meas_huberM)*cvxpy.sum(cvxpy.huber(meas_resids, meas_huberM)) # CVXPY quirk: norm(, 1) != sum(abs()) for matrices
 
     problem = cvxpy.Problem(cvxpy.Minimize(objective))
-    try: problem.solve(solver=cvxpy.CLARABEL); print("CLARABEL succeeded")
+    try: problem.solve(solver=cvxpy.CLARABEL)
     except cvxpy.error.SolverError: pass # Could try another solver here, like SCS, but slows things down
 
     if x_states.value is None: return np.full((A.shape[0], N), np.nan) # There can be solver failure, even without error

examples/2a_optimizing_parameters_with_dxdt_known.ipynb

@@ -60,10 +60,10 @@
                   'pad_to_zero_dxdt': {True, False},
                   'high_freq_cutoff': [1e-3, 5e-2, 1e-2, 5e-2, 1e-1]}, # give numerical params in a list to scipy.optimize over them
                  {'high_freq_cutoff': (1e-5, 1-1e-5)}),
-    rbfdiff: ({'sigma': [1e-3, 1e-2, 1e-1, 1],
+    rbfdiff: ({'sigma': [1e-2, 1e-1, 1],
                 'lmbd': [1e-3, 1e-2, 1e-1]},
-              {'sigma': (1e-3, 1e3),
-                'lmbd': (1e-4, 0.5)}),
+              {'sigma': (1e-2, 1e3),
+                'lmbd': (1e-3, 0.5)}),
     tvrdiff: ({'gamma': [1e-2, 1e-1, 1, 10, 100, 1000],
                'order': {1, 2, 3}}, # warning: order 1 hacks the loss function when tvgamma is used, tends to win but is usually suboptimal choice in terms of true RMSE
               {'gamma': (1e-4, 1e7)}),
@@ -87,23 +87,15 @@
            'forwardbackward': {True, False}},
                         {'q': (1e-10, 1e10),
                          'r': (1e-10, 1e10)}),
-    # robustdiff: ({'order': {1, 2, 3}, # warning: order 1 hacks the loss function when tvgamma is used, tends to win but is usually suboptimal choice in terms of true RMSE
-    #               'log_q': [1., 4, 8, 12], # decimal after first entry ensure this is treated as float type
-    #               'log_r': [-1., 1, 4, 8],
-    #         #'proc_huberM': [0., 2, 6], # 0 is l1 norm, 1.345 is Huber 95% "efficiency", 2 assumes about 5% outliers,
-    #         'meas_huberM': [0., 2, 6]}, # and 6 assumes basically no outliers -> l2 norm. Try (1 - norm.cdf(M))*2 to see outlier portion
-    #              {'log_q': (-1, 18),
-    #               'log_r': (-5, 18),
-    #         'proc_huberM': (0, 6),
-    #         'meas_huberM': (0, 6)}),
     robustdiff: ({'order': {1, 2, 3}, # warning: order 1 hacks the loss function when tvgamma is used, tends to win but is usually suboptimal choice in terms of true RMSE
-                  'log_q': [1., 4, 8, 12], # decimal after first entry ensure this is treated as float type
-                  'log_r': [-1., 1, 4, 8],
-               #'qr_ratio': [10**k for k in range(-1, 16, 4)],
-                 'huberM': [0., 5, 10]}, # 0. so type is float. Good choices here really depend on the data scale
-                 {'log_q': (0, 18),
-                  'log_r': (-5, 10),
-                 'huberM': (0, 20)}), # really only want to use quadratic when nearby; 20sigma is a huge distance
+                  'log_q': [1., 4, 7, 10, 13, 16], # decimal after first entry ensure this is treated as float type
+                  'log_r': [-1., 2, 5, 8, 11],
+            'proc_huberM': [0., 2, 6], # 0 is l1 norm, 1.345 is Huber 95% "efficiency", 2 assumes about 5% outliers,
+            'meas_huberM': [0., 2, 6]}, # 6 assumes basically no outliers per outlier_portion = (1 - norm.cdf(M))*2
+                 {'log_q': (-5, 18),
+                  'log_r': (-5, 18),
+            'proc_huberM': (0, 6),
+            'meas_huberM': (0, 6)}),
     lineardiff: ({'kernel': 'gaussian',
                    'order': 3,
                    'gamma': [1e-1, 1, 10, 100],

examples/2b_optimizing_parameters_with_dxdt_unknown.ipynb

@@ -51,7 +51,7 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs)
     (constant_acceleration, {'r':1e-3, 'q':1e4}), (constant_acceleration, [1e-3, 1e4]),
     (constant_jerk, {'r':1e-4, 'q':1e5}), (constant_jerk, [1e-4, 1e5]),
     (rtsdiff, {'order':2, 'log_qr_ratio':7, 'forwardbackward':True}),
-    #(robustdiff, {'order':3, 'log_q':8, 'log_r':0}), # Add back when design stabilizes
+    (robustdiff, {'order':3, 'log_q':7, 'log_r':2}),
     (velocity, {'gamma':0.5}), (velocity, [0.5]),
     (acceleration, {'gamma':1}), (acceleration, [1]),
     (jerk, {'gamma':10}), (jerk, [10]),
@@ -223,10 +223,10 @@ def spline_irreg_step(*args, **kwargs): return splinediff(*args, **kwargs)
               [(-2, -3), (0, 0), (0, -1), (1, 1)],
               [(-1, -2), (1, 1), (0, -1), (1, 1)],
               [(0, 0), (3, 3), (0, 0), (3, 3)]],
-    robustdiff: [[(-14, -15), (-17, -17), (0, -1), (1, 1)],
-                 [(-14, -14), (-13, -13), (0, -1), (1, 1)],
-                 [(-14, -14), (-13, -13), (0, -1), (1, 1)],
-                 [(-12, -12), (-2, -2), (0, -1), (1, 1)],
+    robustdiff: [[(-15, -15), (-14, -15), (0, -1), (1, 1)],
+                 [(-14, -15), (-13, -13), (0, -1), (1, 1)],
+                 [(-14, -15), (-13, -13), (0, -1), (1, 1)],
+                 [(-9, -9), (-2, -2), (0, -1), (1, 1)],
                  [(0, 0), (2, 2), (0, 0), (2, 2)],
                  [(1, 1), (3, 3), (1, 1), (3, 3)]],
     lineardiff: [[(-7, -8), (-14, -14), (0, -1), (0, 0)],

examples/4_performance_analysis.ipynb

@@ -48,9 +48,9 @@ def plot(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth, xlim=None, show_error=True
     axes[1].legend(loc='lower right', fontsize=12)
 
     if show_error:
-        rms_dxdt = rmse(dxdt_truth, dxdt_hat)
+        rmse_dxdt = rmse(dxdt_truth, dxdt_hat)
         R_sqr = error_correlation(dxdt_truth, dxdt_hat)
-        axes[1].text(0.05, 0.95, f"RMSE = {rms_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",
+        axes[1].text(0.05, 0.95, f"RMSE = {rmse_dxdt:.2f}\n$R^2$ = {R_sqr:.2g}",
                      transform=axes[1].transAxes, fontsize=15, verticalalignment='top')
 
     return fig, axes
@@ -80,7 +80,7 @@ def plot_comparison(dt, dxdt_truth, dxdt_hat1, title1, dxdt_hat2, title2, dxdt_h
 
 def robust_rme(x, x_hat, padding=0, M=6):
     """Robustified/Huberized Root Mean Error metric, used to determine fit between smooth estimate and data.
-    Equals np.linalg.norm(x[s] - x_hat[s]) / np.sqrt(N) if M=float('inf'), and dang close for even M=6 or even 2.
+    Equals np.linalg.norm(x[s] - x_hat[s]) / np.sqrt(N) if M=float('inf'), and dang close for M=6 or even 2.
 
     :param np.array[float] x: noisy data
     :param np.array[float] x_hat: estimated smoothed signal, returned by differentiation algorithms in addition