Huberized TVR

[pavelkomarov]

This one is pretty easy to code. We can replace cvxpy.sum_squares(y - x) in TVR with sigma**2 * cvxpy.sum(cvxpy.huber((y - x)/sigma, M)), where sigma is calculated robustly with median_absolute_deviation as in the integration constant estimation code in utility.py and in evaluate.robust_rme. We might also want to add a huber_const like in kalman_smooth.convex_smooth so the Huber can't flatten out too much as $M \to 0$.

This adds a hyperparameter to the method, $M$, which raises its total to 3 (along with order and $\gamma$). Would need some changes in the default optimization dict at the top of _optimize.py too. But M shouldn't be too hard to optimize over if the $\sigma$ normalization and huber_const are used, because M=0 is an $\ell_1$ norm, M=float('inf') is an $\ell_2$ norm, and M = 2 or 6 or whatever has a standard interpretation of becoming linear at that many sigma from the mean error, such that by $6\sigma$ there are really hardly any inliers that would be treated by the linear region of the Huber, so they get the fit benefits of being treated with the parabolic region. I've accounted for this notion in the default optimization dictionary for robustdiff by letting huberM start at values [0., 2, 6]. (Note the period is important so it's treated as a continuous float by the optimizer and not cast back to discrete ints.)

With that understanding, the code changes are small, and the heaviest part of the work comes after, running experiments using notebook 4 (probably modified to not rerun absolutely everybody) to see how this thing improves on the previous design. More than likely it would slay at the outlier case, but we gotta make sure. Note that when there are outliers, optimize has a huberM of its own and uses evaluate.robust_rme in its loss calculations by default, but with huberM=6 (which shows no serious difference with huberM=float('inf'), interestingly). That should maybe be brought lower to something like 2 if outliers are expected. I believe a ternary to make this selection is already in situ within notebook 4.

[pavelkomarov]

I've realized this is actually mildly complicated by the fact you can't just go and calculate the $\sigma_\text{MAD}$ of y - x, because y isn't a vector of values; it's cvxpy.Variables. sigma could always be given as yet another hyperparam, or I can try to scale in a different way and let M be in those units instead, since the net effect of the two would just be to make the second entry of a huber() into M*sigma instead of just M.

[pavelkomarov]

More word on this, because it turns out to be more complicated than it appeared.

Optimizing $\sum \text{Huber}(\mathbf{y} - \mathbf{\hat{x}}, M) +\gamma\text{TV}(\frac{d^\nu}{dt^\nu} \mathbf{\hat{x}})$ as $M \to 0$ is highly problematic, because it means you're seeking a sparse residual as well as a sparse derivative. The thing is, the $\ell_1$ norm actually penalizes more than the $\ell_2$ norm as your estimate gets close, so an $\ell_1$-based solution is much more likely to fit through (or very close to) more of the data points, $\mathbf{y}$, versus find a neat middle path that sits smoothly between data points like when regressing with the $\ell_2$ norm. The derivatives of such a thing can be sparse if it's composed of piecewise flat, linear, or quadratic bits, and $L(\Phi)$ loves that, leaving you with a terrible "optimal" fit.

Interestingly, when optimizing to derivative RMSE directly, some come out with huberM=0 as the optimal choice, but they absolutely crank $\gamma$ to compensate, e.g. Optimal parameters: {'gamma': 1400.0, 'huberM': 0.0, 'order': 2}.

This same thing doesn't seem to happen as badly when proc_huberM is decreased in RobustDiff, and I've been struggling to imagine why. I think what's going on there is your $\mathbf{Q}^{-1/2}(\mathbf{x}_{1:} - \mathbf{A}\mathbf{x}_{:-1})$ residuals term holds residuals for a bunch of different orders of information, residual position, residual velocity, etc., so when you slap a $\sum |\cdot|_1$ over all the columns on there (equivalently sum(abs(resids_matrix)) in code), it doesn't have the power to force sparsity across all those interlinked quantities in quite the same way as the TVR case, so it doesn't hack $L(\Phi)$ quite as badly.

So how can we combat this? The way I found is to force tvrdiff's huberM to have a floor of 2 or so, to ensure there is enough space for the $\ell_2$ norm to operate and allow the solution to smooth out and not be jerked around as much by nearby data points. But it's kind of a hack. If $\mathbf{y}$ truly is sparse, then it makes sense to allow $M$ to decrease, but this is not really our typical situation here with noisy data; that would be more of a sparse control inputs thing.

The wider issue here is that devising a good $L(\Phi)$ is hard. Even up to now I've observed some functions hacking its $\text{TV}$ component by favoring order=1 solutions, for instance, when those are basically wholly inappropriate for all but the Triangles simulation. I use search_space_updates to limit the order in those cases artificially, reasoning a user would be smart enough to realize what's going on and do the same, but it's some pretty advanced usage going on there, and I don't love how finicky this can be. Especially in this $\ell_1$ process term thing, it has been a nightmare to figure out and experiment around.

[pavelkomarov]

Okay, here are some plots. I gotta say that the before and afters of all these look essentially exactly the same. The way to see it is to open them all in Preview and flip back and forth really quickly. The $\gamma$ move in some cases, but hardly at all. I guess that's good; the optimization still finds the optimum it did before in the absence of outliers.
Archive.zip

The one set of plots that looks different and merits actual display here is for outliers. Before:

and after:

Sick, right? Well, kind of. The method is doing better than before, but not super excellently, and not as well as robustdiff in terms of error correlation. Does it have the expressive power to handle outliers even better than shown here? Probably, because optimizing for RMSE gets pretty good-looking plots, but optimization of $L(\Phi)$ is just pain with this thing, as detailed above, so it's unclear whether this method can really be said to check the box.

But in any event, I'm done with it. It's nice to have the extra hyperparameter as a handhold in case someone wants to use it.