Optimization Loss isn't Outlier Robust

[pavelkomarov]

I've been improving RobustDiff for #165, but as I've done so I've found that the optimization results are actually getting worse. The degradation starts when I decouple qr_ratio into log_q and log_r, because the scale of the two matters relative to huberM. Suddenly the method has more expressive power and can find lower-cost solutions, but with worse RMSE and R^2 against the true underlying derivative. I've realized $L(\Phi)$ isn't itself robust, so the hyperparameter optimization is favoring solutions that bend in the direction of outliers.

The solution? Pop a Huber loss in the RMSE evaluation, to get something we might call "Root Mean Huber Error" or "Robust Root Mean Error". We now have another parameter to pick, where that Huber switches from quadratic to linear. If we normalize the inputs by their standard deviation, then we can choose the parameter as some number of sigma, like 2 or 3, which would then count ~95% or ~99.73% of inliers as inliers, assuming a Gaussian distribution. This will affect the scale of the first term in the loss function, so we'll have to account for this against the total variation smoothing term.

The TV term is potentially problematic itself too, because if we are allowed to independently drive down RobustDiff's process term's Huber parameter, optimization might favor approximating the 1-norm to make it artificially sparse. I've run in to this in other cases too, and disallow order 1 for TVR for most datasets, because it can "hack" the loss function.

[pavelkomarov]

This issue goes kind of deep. Actually even estimating integration constants needs to be made robust, otherwise the loss estimate will be off.

def estimate_integration_constant(x, x_hat, M=float('inf')):
    """Integration leaves an unknown integration constant. This function finds a best fit integration
    constant given x and x_hat (the integral of dxdt_hat) by optimizing :math:`\\min_c ||x - \\hat{x} + c||_2`.

    :param np.array[float] x: timeseries of measurements
    :param np.array[float] x_hat: smoothed estimate of x
    :param float M: If given, robustifies constant estimation using Huber loss

    :return: **integration constant** (float) -- initial condition that best aligns x_hat with x
    """
    if M == float('inf'): # calculates the constant to be mean(diff(x, x_hat)), equivalent to argmin_{x0} ||x_hat + x0 - x||_2^2
        return np.mean(x - x_hat) # Solves the L2 distance minimization
    elif M == 0:
        return np.median(x - x_hat) # Solves the L1 distance minimization
    else:
        sigma = np.std(x - x_hat) # for normalizing the input to huber
        return minimize(lambda x0: np.mean(huber((x - (x_hat+x0))/sigma, M)), # fn to minimize in 1st argument
            0, method='SLSQP').x[0] # result is a vector, even if initial guess is just a scalar

Even this isn't quite right, because the standard deviation is thrown off a lot by outliers, so in order to normalize huber inputs, we really have to use a winsorized stddev or something. What a pain.

[pavelkomarov]

Note that when I replaced $\text{trapz}(\mathbf{\hat{\dot{x}}}(\Phi)) + c$ with $\mathbf{\hat{x}}$, which seems simpler and equally valid, the optimization results were not as good. I believe this is because $\mathbf{\hat{\dot{x}}}$ and $\mathbf{\hat{x}}$ are related through a complicated differentiation method, and optimizing with the former is more direct. This makes the full loss function:

$$L(\Phi) = \sigma_\text{MAD} \sqrt{\frac{2}{N}\sum\text{Huber}\bigg( \frac{\text{trapz}(\mathbf{\hat{\dot{x}}}(\Phi)) + c - \mathbf{y}}{\sigma_\text{MAD}}, M \bigg)} + \gamma \bigg(\text{TV}\big(\mathbf{\hat{\dot{x}}}(\Phi)\big)\bigg)$$

where $\Phi$ are the hyperparameter settings of a differentiation method, $\mathbf{\hat{\dot{x}}}(\Phi)$ is an estimate of the derivative, the left term is a generalization of $\text{RMSE}$ that uses the Huber loss, $\sigma_\text{MAD}$ is a robust scatter measure (median absolute deviation), $M$ is a robustness parameter in units of $\sigma_\text{MAD}$ (often set to 2 or 3 so the vast majority of inliers are measured with RMSE), $\text{trapz}$ performs trapezoidal integration, $c$ corrects for a lost constant of integration (outlier-robustly), $\mathbf{y}$ is noisy data, $\gamma$ is a hyper-parameter, and $TV$ is the total variation:

$$\text{TV}(\mathbf{\hat{\dot{x}}}) = \frac{1}{m}\left\lVert\mathbf{\hat{\dot{x}}}_{0:m-1}-\mathbf{\hat{\dot{x}}}_{1:m}\right\rVert_{1}.$$