Intentionally picked Huber M in the optimization loss to be 2, which appears not to hurt performance when there are no outliers but helps the optimization better balance fidelity with smoothness when there are outliers, reran notebook and improved some of the explanatory math

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)}),
@@ -81,29 +81,21 @@
     rtsdiff: ({'forwardbackward': {True, False},
                          'order': {1, 2, 3}, # for this few options, the optimization works better if this is categorical
                   'log_qr_ratio': [float(k) for k in range(-9, 10, 2)] + [12, 16]},
-                 {'log_qr_ratio': [-10, 20]}), # qr_ratio is usually >>1
+                 {'log_qr_ratio': (-10, 20)}), # qr_ratio is usually >>1
     constant_velocity: ({'q': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8], # Deprecated method
                          'r': [1e-8, 1e-4, 1e-1, 1e1, 1e4, 1e8],
            '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], # 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, 16),
+                  'log_r': (-5, 16),
+            'proc_huberM': (0, 6),
+            'meas_huberM': (0, 6)}),
     lineardiff: ({'kernel': 'gaussian',
                    'order': 3,
                    'gamma': [1e-1, 1, 10, 100],
@@ -161,7 +153,8 @@ def _objective_function(point, func, x, dt, singleton_params, categorical_params
             ec = evaluate.error_correlation(dxdt_truth, dxdt_hat, padding=padding)
             cache[key] = ec; return ec
     else: # then minimize sqrt{2*Mean(Huber((x_hat- x)/sigma))}*sigma + gamma*TV(dxdt_hat)
-        cost = evaluate.robust_rme(x, x_hat, padding=padding) + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding)
+        # Huber M=2 means more than 95% of inliers (assuming Gaussianity) are treated with RMSE, while 
+        cost = evaluate.robust_rme(x, x_hat, padding=padding, M=2) + tvgamma*evaluate.total_variation(dxdt_hat, padding=padding)
         cache[key] = cost; return cost
 
 

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

@@ -101,8 +101,7 @@ def test_simulations(request):
         from matplotlib import pyplot
         fig, axes = pyplot.subplots(2, 3, figsize=(18,7), constrained_layout=True)
 
-    for i,(sim,title) in enumerate(zip(
-        [pi_cruise_control, sine, triangle, pop_dyn, linear_autonomous, lorenz_x],
+    for i,(sim,title) in enumerate(zip([pi_cruise_control, sine, triangle, pop_dyn, linear_autonomous, lorenz_x],
         ["Cruise Control", "Sum of Sines", "Triangles", "Logistic Growth", "Linear Autonomous", "Lorenz First Dimension"])):
 
         y, x, dxdt = sim(duration=4, dt=0.01, noise_type='normal', noise_parameters=[0,0.1])