Reinclude Chebyshev, or actually don't and remove the traces and tell users why

[pavelkomarov]

Hi Floris, I've finally finished all work on my spectral-derivatives repo and am turning full attention to pynumdiff. I'm using issues as a task tracker for myself, as I'm wont to do, and decided the upstream repo is better for this rather than my own fork, both so you can feed back as I go if you want to and so the record of these things is long-lived after I finish and nuke my own fork. Hope the emails don't get too annoying.

First up, I spent a long while working on trying to make my spectral methods better in the presence of noise, but the end result of that is mostly the realization that polynomial bases are just not very good for this, because, as you discovered in the past, there is systematic sensitivity and blowup at the domain edges for theoretical reasons.

In the past you removed Chebyshev-based methods because (1) performance was bad in the presence of noise and (2) your code (now remaining as comments) was based on a long-unsupported package. I can offer to solve the second problem with spectral-derivatives, and it might be nice to offer this method in the case someone wants to numerically differentiate clean data. I propose we do that to clean up the ghostly code, or we make the firm decision not to support this method and then put an argument for why somewhere in the documentation or papers.

[pavelkomarov]

Another potential solution is to remove all support for spectral methods in this repo and instead direct users to spectral-derivatives, where I already support multidimensional input and arbitrary sampling points.

[pavelkomarov]

As of #115 the linear module is looking a touch better, but we need to figure out what to do about this one to really complete its makeover. Now that the unit tests can --plot stuff for us nicely, it will be pretty easy to see what this thing is doing and toy with it. I'll move over to local dependency on spectral-derivatives and post some results plots here. If I can get it to work under any circumstances, great. If not, we'll have a record of why it's a bad idea. I didn't try any sliding window stuff in my repo, so I'm curious whether that helps.

[pavelkomarov]

I sat down to do this and realized that a sliding window with Chebyshev makes no sense, because the whole advantage of Chebyshev is being able to fit a polynomial basis in $O(N \log N)$ via FFT when points are cosine-spaced, and if you take a portion of cosine-spaced points, they're no longer cosine-spaced. You also can't add padding to cosine-spaced points without again breaking cosine-spacing, so if we wanted to do cosine-spacing, we'd end up having to pad and then interpolate. Gross.

If points are equispaced, as the assumption in the code, then you get no savings fitting Chebyshev polynomials vs any other kind of polynomials, so a sliding-window implementation would essentially be redundant with polydiff.

Here is such an implementation, making use of the numpy.polynomial.chebyshev module:

def chebydiff(x, dt, poly_order, window_size=None, step_size=1, kernel='friedrichs'):
    """Fit Chebyshev polynomials to the data, and differentiate those

    :param np.array[float] x: data to differentiate
    :param float dt: step size
    :param int poly_order: keep polynomials up to this order
    :param int window_size: size of the sliding window, if not given no sliding
    :param int step_size: step size for sliding
    :param str kernel: name of kernel to use for weighting and smoothing windows ('gaussian' or 'friedrichs')

    :return: tuple[np.array, np.array] of\n
             - **x_hat** -- estimated (smoothed) x
             - **dxdt_hat** -- estimated derivative of x
    """
    if window_size % 2 == 0:
        window_size += 1
        warn("Kernel window size should be odd. Added 1 to length.")

    def _chebdiff(x, dt, poly_order, weights=None):
        t = np.arange(len(x))*dt

        r = np.polynomial.chebyshev.chebfit(t, x, poly_order, w=weights) # chebfit returns lowest order first
        dr = np.polynomial.chebyshev.chebder(r) # series derivative rule already implemented for us

        dxdt_hat = np.polynomial.chebyshev.chebval(t, dr) # evaluate the derivative and original polynomials at points t
        x_hat = np.polynomial.chebyshev.chebval(t, r) # smoothed x

        return x_hat, dxdt_hat

    if not window_size:
        return _chebdiff(x, dt, poly_order)

    kernel = {'gaussian':utility.gaussian_kernel, 'friedrichs':utility.friedrichs_kernel}[kernel](window_size)
    return utility.slide_function(_chebdiff, x, dt, kernel, poly_order, stride=step_size, pass_weights=True)

A not-sliding-window approach I tried ad nauseam in a notebook on my spectral-derivatives repo, and it doesn't really work, although I could try to ape spectaldiff from this file as a last attempt.

[pavelkomarov]

Okay, here are plots against the test functions for the version above with {'poly_order':3, 'window_size':11, 'step_size':1, 'kernel':'friedrichs'}.

So it works fine. Interestingly, if I run the optimization code in notebook 2a, using the same default search_space for both (see polydiff's entry in the dictionary at the top of _optimize.py), since they have the same parameters

params, val = optimize(chebydiff, x, dt, dxdt_truth=dxdt_truth)
print('Optimal parameters: ', params)
x_hat, dxdt_hat = chebydiff(x, dt, **params)
evaluate.plot(x, dt, x_hat, dxdt_hat, x_truth, dxdt_truth)

then I get practically exactly the same answers for polydiff:

Optimal parameters: {'step_size': 2, 'poly_order': 2, 'window_size': 29, 'kernel': 'friedrichs'} RMS error in velocity: 0.7277649206991202

and this version of chebydiff:

Optimal parameters: {'step_size': 2, 'poly_order': 2, 'window_size': 29, 'kernel': 'friedrichs'} RMS error in velocity: 0.7277649206991037

So that's indeed redundant. We'll move on to a version that's shaped more like spectraldiff with its edge-padding logic. I already tried edge padding in may own repo and didn't have much success, but I also didn't think it should work with Fourier-spectral, yet it seems to do okay, so maybe I missed something. Worth a quick look.

[pavelkomarov]

Okay, here is an alternative implementation inspired by spectraldiff:

def chebydiff(x, dt, relative_cutoff=None, pad_to_zero_dxdt=True):
    """Fit Chebyshev polynomials to the data, and differentiate those

    :param np.array[float] x: data to differentiate
    :param float dt: step size
    :param float relative_cutoff: in 0 to 1
    :param bool pad_to_zero_dxdt: if True, extend the data to smoothly force the derivative to zero.

    :return: tuple[np.array, np.array] of\n
             - **x_hat** -- estimated (smoothed) x
             - **dxdt_hat** -- estimated derivative of x
    """
    L = len(x)

    # make derivative go to zero at ends (optional)
    if pad_to_zero_dxdt:
        padding = min(100, L)
        pre = x[0]*np.ones(padding) # extend the edges
        post = x[-1]*np.ones(padding)
        x = np.hstack((pre, x, post))
        kernel = utility.mean_kernel(padding//2)
        x_hat = utility.convolutional_smoother(x, kernel) # smooth the edges in
        x_hat[padding:-padding] = x[padding:-padding] # replace middle with original signal
        x = x_hat
    else:
        padding = 0

    cutoff_wavenumber = len(x)-1 # need degree N-1 to go through N points
    if relative_cutoff != None: cutoff_wavenumber = int(cutoff_wavenumber*relative_cutoff)

    t = np.arange(len(x))*dt
    r = np.polynomial.chebyshev.chebfit(t, x, cutoff_wavenumber)
    dr = np.polynomial.chebyshev.chebder(r) # series derivative rule already implemented for us

    dxdt_hat = np.polynomial.chebyshev.chebval(t, dr) # evaluate the derivative and original polynomials at points t
    dxdt_hat = dxdt_hat[padding:L+padding]

    x_hat = np.polynomial.chebyshev.chebval(t, r) # smoothed x
    x_hat = x_hat[padding:L+padding]

    # Integrate to get x_hat, alternative way
    #x_hat = utility.integrate_dxdt_hat(dxdt_hat, dt)
    #x0 = utility.estimate_integration_constant(x[padding:L+padding], x_hat)
    #x_hat = x_hat + x0

    return x_hat, dxdt_hat

If I test it with pad_to_zero_dxdt=False, I see it fits the function well in the absence of noise and blows up at the edges in the presence of noise, just as in my earlier experimentation on my own repo:

If I pad and use arbitrarily-chosen freq_cutof=0.7, then I get something that looks kinda okay-ish, but not spectacular:

If I run it though the optimizer as in notebook 2a using the following addition to _optimize.py

chebydiff: ({'pad_to_zero_dxdt': [True, False],
             'freq_cutoff': [0.1, 0.25, 0.5, 0.7, 1]},
            {'freq_cutoff': (0, 1)})

then I see:

** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
[/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626](https://file+.vscode-resource.vscode-cdn.net/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626): RankWarning: The fit may be poorly conditioned
  return pu._fit(chebvander, x, y, deg, rcond, full, w)
[/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626](https://file+.vscode-resource.vscode-cdn.net/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626): RankWarning: The fit may be poorly conditioned
  return pu._fit(chebvander, x, y, deg, rcond, full, w)
[/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626](https://file+.vscode-resource.vscode-cdn.net/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626): RankWarning: The fit may be poorly conditioned
  return pu._fit(chebvander, x, y, deg, rcond, full, w)
[/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626](https://file+.vscode-resource.vscode-cdn.net/usr/local/lib/python3.13/site-packages/numpy/polynomial/chebyshev.py:1626): RankWarning: The fit may be poorly conditioned
  return pu._fit(chebvander, x, y, deg, rcond, full, w)
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value
** On entry to DLASCL, parameter number  4 had an illegal value

Optimal parameters: {'pad_to_zero_dxdt': np.False_, 'freq_cutoff': np.float64(0.25625)} RMS error in velocity: 1.69281054090109

Gross. I think that result speaks for itself.

I'm going to rip out all Chebyshev stuff from this repo, because it's either:

• completely redundant with polydiff if we limit it to the local case or
• less numerically stable and less efficient than spectraldiff for a result that looks worse in the global case.

The non-uniform, global Chebyshev basis functions turn out to be great for expressive power, but expressive power is what you want to regularize away when dealing with noise, so Chebyshev polynomials are exactly the wrong tool. Fundamentally, even if you could pad as much as you wish and sample at the optimal Chebyshev-Lobatto nodes, the sensitivity of the method in the center is different from that as you go toward the domain edges, so high-frequency noise filtration isn't consistent across the domain, resulting in fits that grow more erroneous away from the center. It's just overkill when simple iterated FD works so well.