Outlier-Robust Differentiation Method

[pavelkomarov]

Right now we have a standard Kalman filter, which can work with given dynamics or assumed constant-velocity or constant acceleration dynamics, and we have forward-backward smoothing. We do not have anything resembling a nonlinear or robust filter.

[pavelkomarov]

When addressing this, we'll have to face potential irregularity across methods: The EKF, UKF, RKF, and $\mathcal{H}_\infty$ are similar cases where the user provides system information via A, and the code doesn't natively know which index of the state holds the original signal and which holds the derivative. In fact, there is no guarantee that the state even keeps track of the derivative in this general formulation (although the way we set up most problems does include derivatives in the state). As a consequence, this function can only return the RTS-smoothed state, can't pick out the signal and its derivative, so the return value of this function doesn't match the common x_hat, dxdt_hat standard.

With a clever enough design, we may be able to enforce the user give us information with enough regularity that we can return x_hat, dxdt_hat like we're used to. But it's also somewhat challenging to get the interface right when we may want to support the user passing observation and state evolution functions instead of matrices and when we may want to have unfixed dt per #1 and #80. I'm open to interface suggestions.

[pavelkomarov]

I've been doing a little more reading, and it seems like EKF is considered pretty rudimentary. UKF generally does better at handling greater nonlinearity and is cheaper for systems that don't have really high dimensional states. Nonlinear models are a pretty rare case in simple differentiation, so it's probably not worth implementing both, might only barely be worth implementing one, and then largely the benefit is the edifying education of the the collaborator who takes it on. Here are some papers on the topic to draw from https://groups.seas.harvard.edu/courses/cs281/papers/unscented.pdf, https://ieeexplore.ieee.org/document/1271397, https://www.cs.unc.edu/~welch/kalman/media/pdf/Julier1997_SPIE_KF.pdf

Regarding the robust case, $\mathcal{H}_\infty$ filtering is often much more robust than necessary, useful when you truly know nothing about the noise. It really isn't made to handle outliers and tends to have less optimal results in terms of norms we intuit than the vanilla Kalman filtering result. Indeed, the Kalman filter is $\mathcal{H}_2$ optimal regardless of the disturbance, so long as you know its covariances; it's just that when the disturbance is Gaussian noise $\mathcal{H}_2$ happens to find you the Maximum Likelihood Estimator, which will be different with different noise assumptions. So honestly in our situation the kind of robustness we want is probably more along the lines of https://sites.math.washington.edu/~burke/papers/reprints/89-Survey-GKS.pdf, https://arxiv.org/pdf/2405.05646, https://www.sciencedirect.com/science/article/pii/0167691194901066/pdf?md5=9928110b2442a4b7dbb5a216a172d89b&pid=1-s2.0-0167691194901066-main.pdf, https://www.kybernetika.cz/content/1991/6/481/paper.pdf, or https://pmc.ncbi.nlm.nih.gov/articles/PMC4435213/. That's a lot to sort through, and I'll do what I can, but we certainly don't have the bandwidth or the need to implement all approaches. We can pick one, or we can even just have some docs that point users to this issue or to those papers to show we've been thinking about the possibility of generalizing in this direction, so there is a place to get started if and when the need ever eventually arises.

As far as interface concerns, I believe we can simply demand that users model their system so the smoothed variable they care about is in the first dimension, and its derivative is in the second. Rows can be flip flopped around without a problem in a linear model, and nonlinear functions can work on flip flopped indices, so this is an easy condition to meet. Variable dt I don't think will get in our way at the interface level, because it's very easy to make the second parameter optionally a vector of data locations, and all the complicated business can be handled internally by the module. Similarly, taking state evolution and measurement functions rather than matrices is a very easy and intuitive swap-out in a method signature, since functions can be passed around like first-class objects in Python, so I foresee very little road bumps there.

[pavelkomarov]

Here is an implementation of $\mathcal{H}_\infty$ filtering compared with Kalman filtering, although I'm not 100% sure this implementation is correct, because the old papers and textbooks on $\mathcal{H}_\infty$ filtering are really abstruse, mostly deal with continuous time formulations, and seem to have pretty varied notation that doesn't seamlessly agree with the common Kalman filter formulations. The example system I'm using is the constant-acceleration one, which turns out not to have a bounded $\mathcal{H}_\infty$ norm, nor $\mathcal{H}_2$ norm, despite being still workable with Kalman Filter. The $\mathcal{H}_p$ paradigm is over-general and weaker than desired in this case.

filt_comparison.ipynb

[pavelkomarov]

The more I think about it, the more I realize this situation could be a whole git repo of its own. We should maybe pull in @saravkin to talk about what repos of his already exist to do some of this. What is the right level of detail if we want to address the problem of differentiation in the presence of outliers? I imagine we don't truly need all the Piecewise Linear-Quadratic Modeling machinery, but we could use some of it and use an assumed linear constant-derivative model to give a nice robustdiff(x, t, **kwargs)-type interface here, backed by CVXPY to solve $\ell_1$ norms, which is already a dependency. And whichever supporting functions make sense to expose can be exposed to enable more sophisticated modeling should users decide they want to go full bananas, like how we presently make kalman_filter and rts_smooth public.

[pavelkomarov]

Lying in shavasana after yoga on Saturday, noticing my mind racing like normal, I decided that I'd be willing to close this issue if we made an $\ell_1$-based robustdiff that worked okay. I'll look in to coding it up, because it shouldn't be so hard from where we're at. Everything beyond handling outliers can be considered out of scope for this project.