Numerical issues in scipy's Savitzky Golay filter coefficients for large polynomial order

Question

Consider the design of a Savitzky-Golay filter of window length 101 and (high) polynomial order 20. Using scipy version 1.10.1, the filter coefficients can be obtained in python as:

from scipy import signal
h = signal.savgol_coeffs(window_length=101, polyorder=20)

leading to

For reference, in Matlab r2022b

order=20; framelen=101; b = sgolay(order,framelen); h = b((framelen+1)/2,:);

gives

which looks more reasonable. For lower order polynomials, the results coincide much more closely.

It looks like scipy's implementation solves a least squares problem where the design matrix is computed as:

 order = np.arange(polyorder + 1).reshape(-1, 1)
 A = x ** order

leading to a condition number > 1e29!

Is there some other way to get these coefficients in python?

Answer 1

Some languages offers higher accuracy data types.
For instance, on Python, you may use the BigFloat package.
There is always enough bits you can throw at it and get something working.

Yet a better approach is to use a different approach.
Since you're after a polynomial like model, you could use Kalman Filter to get similar results with better numerical properties.
The trick here will be tweaking the ratio between $ \boldsymbol{Q} $ and $ \boldsymbol{R} $ to have similar results.

With some tweaking and hyper parameters grid search you will be able to easily do so.

Remark: If you share some signals we can try to implement the approach.