From 98e9ee3fde6276d89d25891a99bdbcc960d2ed82 Mon Sep 17 00:00:00 2001 From: MothNik Date: Thu, 20 Jun 2024 22:58:13 +0200 Subject: [PATCH 1/9] python/scipy/savgol-precision: Added a numerically stabilised version of the Savitzky Golay filter computation as described in Scipy Issue #20825 https://github.com/scipy/scipy/issues/20825#issuecomment-2181539118 --- python/scipy/savgol-precision/savgol_work.py | 288 +++++++++++++++++-- 1 file changed, 268 insertions(+), 20 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index 6026da9..62494f4 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -1,10 +1,216 @@ +from math import comb, fsum +from typing import Literal, Optional, Union -from mpmath import mp -import mpsig import numpy as np from matplotlib import pyplot as plt +from mpmath import mp from scipy._lib._util import float_factorial -from scipy.linalg import lstsq + +import mpsig + + +def _pinv_correction_factors_normalized_polyvander( + pinv_row_idx: int, + polyorder: int, + polyvander_column_scales: np.ndarray, + x_center: float, + x_scale: float, +) -> np.ndarray: + """ + Computes the correction factors for the pseudo-inverse of a normalized polynomial + Vandermonde matrix where the x-variable was normalized as ``x_scaled = (x - x_center) / x_scale``. + After this normalization, the columns of the Vandermonde matrix ``A`` given by + ``A = polyvander(x_scaled, polyorder)`` are scaled to have unit Euclidean norm. + These steps ensure the numerical stability of the pseudo-inverse computation, but + distorts the pseudo-inverse of the polynomial fit. It will then be matched to the + scaled columns of ``A`` who were themselves normalized before, rather than the + original columns of ``A`` that are based on ``x``. + + This function computes the correction factors to be applied to the rows of the + pseudo-inverse of ``A`` to recover the pseudo-inverse of the polynomial fit based + on the original columns of ``A``. + Multiplying out all scaling and centering factors and rearranging terms gives the + following equation for the correction factors: + + ``` + phi_ik = ((-1)**(i - k)) * comb(i, k) * ((1.0 / x_scale)**i) * 1.0 / col_scale[i] * (x_center ** (i - k)) + ``` + + where + + - ``k`` is the row index of the pseudo-inverse matrix to be corrected (the ``k``-th + row corresponds to the ``k``-th coefficient of the polynomial fit, e.g., ``k=0`` + corresponds to the constant term ``a0``), + - ``i`` is the iterator for the polynomial order + - ``comb(i, k)`` is the binomial coefficient, and + - ``col_scale[i]`` is the Euclidean norm of the ``i``-th column of the normalized + polynomial Vandermonde matrix + + These can be applied to correct the element in the ``k``-th row and the ``j``-th + column of the distorted pseudo-inverse ``JD+`` as follows: + + ``` + J+[k, j] = sum(phi_ik * JD+[i, j] for i in range(k, polyorder + 1)) + ``` + + to obtain the corrected pseudo-inverse ``J+``. + + Parameters + ---------- + pinv_row_idx : int + The row index of the pseudo-inverse matrix to be corrected. This corresponds to + ``k`` in the above equation. + polyorder : int + The polynomial order of the polynomial fit. + polyvander_column_scales : np.ndarray of shape (polyorder + 1,) + The Euclidean norms of the columns of the normalized polynomial Vandermonde matrix. + x_center, x_scale : float + The centering and scaling factors used to normalize the x-variable. + + Returns + ------- + pinv_correction_factors : np.ndarray of shape (polyorder + 1 - pinv_row_idx,) + The correction factors for the pseudo-inverse matrix row corresponding to + ``pinv_row_idx``. + + """ # noqa: E501 + + # first, the signs are determined by the (-1)**(i - k) factor ... + i_vect = np.arange(start=pinv_row_idx, stop=polyorder + 1, dtype=np.int64) + signs = (-1) ** (i_vect - pinv_row_idx) + # ... followed by the binomial coefficients ... + binomials = np.array([comb(i, pinv_row_idx) for i in i_vect], dtype=np.int64) + # ... then the factors for the x-variable + x_factors = ( + (x_center ** (i_vect - pinv_row_idx)) + / (x_scale**i_vect) + / polyvander_column_scales[pinv_row_idx::] + ) + + return signs * binomials * x_factors + + +def _correct_pinv_row_normalized_polyvander( + pinv_row_idx: int, + pinv_scaled_polyvander: np.ndarray, + polyorder: int, + polyvander_column_scales: np.ndarray, + x_center: float, + x_scale: float, +) -> np.ndarray: + """ + Corrects the ``k``-th row of pseudo-inverse of a normalized polynomial Vandermonde + matrix to obtain the ``k``-th row of the pseudo-inverse of the polynomial fit based + on the original columns of the Vandermonde matrix. + For the correction it relies on floating point accurate summation to avoid loss of + precision due to cancellation errors with high polynomial orders. + Please refer to the documentation of the function :func:`_pinv_correction_factors_normalized_polyvander` + for a detailed explanation of the correction factors. + + Parameters + ---------- + pinv_row_idx : int + The row index of the pseudo-inverse matrix to be corrected. This corresponds to + ``k`` in the above equation. + pinv_scaled_polyvander : np.ndarray of shape (polyorder + 1, m) + The pseudo-inverse of the normalized polynomial Vandermonde matrix. + polyorder : int + The polynomial order of the polynomial fit. + polyvander_column_scales : np.ndarray of shape (polyorder + 1,) + The Euclidean norms of the columns of the normalized polynomial Vandermonde matrix. + x_center, x_scale : float + The centering and scaling factors used to normalize the x-variable. + + Returns + ------- + pinv_corrected_polyvander : np.ndarray of shape (polyorder + 1, m) + The corrected pseudo-inverse of the polynomial fit that corresponds to the + original columns of the polynomial Vandermonde matrix. + + """ # noqa: E501 + + correction_factors = _pinv_correction_factors_normalized_polyvander( + pinv_row_idx=pinv_row_idx, + polyorder=polyorder, + polyvander_column_scales=polyvander_column_scales, + x_center=x_center, + x_scale=x_scale, + ) + return np.array( + [ + fsum(pinv_scaled_polyvander[pinv_row_idx::, j] * correction_factors) + for j in range(pinv_scaled_polyvander.shape[1]) + ] + ) + + +def super_stabilised_savgol_coeffs( + window_length: int, + polyorder: int, + deriv: int = 0, + delta: Union[float, int] = 1.0, + pos: Optional[int] = None, + use: Literal["conv", "dot"] = "conv", +) -> np.ndarray: + + if polyorder >= window_length: + raise ValueError("polyorder must be less than window_length.") + + halflen, rem = divmod(window_length, 2) + + if pos is None: + if rem == 0: + pos = halflen - 0.5 + else: + pos = halflen + + if not (0 <= pos < window_length): + raise ValueError("pos must be nonnegative and less than " "window_length.") + + if use not in ["conv", "dot"]: + raise ValueError("`use` must be 'conv' or 'dot'") + + if deriv > polyorder: + coeffs = np.zeros(window_length) + return coeffs + + # Form the design matrix A. The columns of A are powers of the integers + # from -pos to window_length - pos - 1. The powers (i.e., rows) range + # from 0 to polyorder. (That is, A is a vandermonde matrix, but not + # necessarily square.) + x = np.arange(-pos, window_length - pos, dtype=float) + + if use == "conv": + # Reverse so that result can be used in a convolution. + x = x[::-1] + + # the stable version of the polynomial Vandermonde matrix is computed + x_min, x_max = x.min(), x.max() + x_center = 0.5 * (x_min + x_max) + x_scale = 0.5 * (x_max - x_min) + x_normalized = (x - x_center) / x_scale + A_scaled_polyvander = np.polynomial.polynomial.polyvander(x_normalized, polyorder) + + # to stabilize the pseudo-inverse computation, the columns of the polynomial + # Vandermonde matrix are normalized to have unit Euclidean norm + a_column_scales = np.linalg.norm(A_scaled_polyvander, axis=0) + A_scaled_polyvander /= a_column_scales[np.newaxis, ::] + + # then, the pseudo-inverse is computed and corrected to obtain the pseudo-inverse + # of the polynomial fit based on the columns the original polynomial Vandermonde + # matrix based on the original x-values + A_pinv_distorted = np.linalg.pinv(A_scaled_polyvander) + + return ( + float_factorial(deriv) / delta**deriv + ) * _correct_pinv_row_normalized_polyvander( + pinv_row_idx=deriv, + pinv_scaled_polyvander=A_pinv_distorted, + polyorder=polyorder, + polyvander_column_scales=a_column_scales, + x_center=x_center, + x_scale=x_scale, + ) def stabilised_savgol_coeffs( @@ -33,8 +239,7 @@ def stabilised_savgol_coeffs( pos = halflen if not (0 <= pos < window_length): - raise ValueError("pos must be nonnegative and less than " - "window_length.") + raise ValueError("pos must be nonnegative and less than " "window_length.") if use not in ["conv", "dot"]: raise ValueError("`use` must be 'conv' or 'dot'") @@ -76,7 +281,7 @@ def stabilised_savgol_coeffs( def simple_savgol_coeffs(window_length, polyorder, pos=None): if window_length % 2 != 1: - raise ValueError('window_length must be odd when pos is None') + raise ValueError("window_length must be odd when pos is None") if pos is None: pos = window_length // 2 t = np.arange(window_length) @@ -87,8 +292,12 @@ def simple_savgol_coeffs(window_length, polyorder, pos=None): def check_savgol_coeffs(c, cref): - relerr = np.array([float(abs((c0 - cref0)/cref0)) if cref != 0 else np.inf - for c0, cref0 in zip(c, cref)]) + relerr = np.array( + [ + float(abs((c0 - cref0) / cref0)) if cref != 0 else np.inf + for c0, cref0 in zip(c, cref) + ] + ) return relerr @@ -98,21 +307,60 @@ def check_savgol_coeffs(c, cref): window_len = 51 order = 10 pos = 35 + deriv = 5 + delta = 0.5 + + coeffs = mpsig.savgol_coeffs(window_len, order, pos=pos, deriv=deriv, delta=delta) + if deriv == 0: + cnp = simple_savgol_coeffs(window_len, order, pos=pos) - coeffs = mpsig.savgol_coeffs(window_len, order, pos=pos) - cnp = simple_savgol_coeffs(window_len, order, pos=pos) - cs = stabilised_savgol_coeffs(window_len, order, pos=pos) + c_sup_s = super_stabilised_savgol_coeffs( + window_len, order, pos=pos, deriv=deriv, delta=delta + ) + # c_s = stabilised_savgol_coeffs(window_len, order, pos=pos, deriv=deriv, delta=delta) - e_cnp = check_savgol_coeffs(cnp, coeffs) - e_cs = check_savgol_coeffs(cs, coeffs) + if deriv == 0: + e_cnp = check_savgol_coeffs(cnp, coeffs) + e_c_sup_s = check_savgol_coeffs(c_sup_s, coeffs) + # e_c_s = check_savgol_coeffs(c_s, coeffs) - plt.plot(e_cnp, 'o', alpha=0.65, - label=f'numpy Polynomial.fit\n(rel err: max {e_cnp.max():8.2e}, mean {np.mean(e_cnp):8.2e})') - plt.plot(e_cs, 'd', alpha=0.65, - label=f'stabilised_savgol_coeffs\n(rel err: max {e_cs.max():8.2e}, mean {np.mean(e_cs):8.2e}') + if deriv == 0: + plt.plot( + e_cnp, + "8", + alpha=0.65, + label=f"numpy Polynomial.fit\n(rel err: max {e_cnp.max():8.2e}, median {np.median(e_cnp):8.2e})", + ) + + plt.plot( + e_c_sup_s, + "d", + alpha=0.65, + label=f"super stabilised_savgol_coeffs\n(rel err: max {e_c_sup_s.max():8.2e}, median {np.median(e_c_sup_s):8.2e}", + ) + # plt.plot(e_c_s, 'x', alpha=0.65, + # label=f'stabilised_savgol_coeffs\n(rel err: max {e_c_s.max():8.2e}, median {np.median(e_c_s):8.2e}') plt.legend(shadow=True, framealpha=1) plt.grid() - plt.title(f'Relative Error of Savitzky-Golay Coefficients\n{window_len=} {order=} {pos=}') - plt.xlabel('Coefficient index') - plt.ylabel('Relative error of coefficient') + plt.title( + f"Relative Error of Savitzky-Golay Coefficients\n{window_len=} {order=} {pos=}, {deriv=}" + ) + plt.xlabel("Coefficient index") + plt.ylabel("Relative error of coefficient") + + fig, ax = plt.subplots() + + ax.plot(coeffs, "o", alpha=0.65, label="mpsig.savgol_coeffs") + if deriv == 0: + ax.plot(cnp, "8", alpha=0.65, label="numpy Polynomial.fit") + + ax.plot(c_sup_s, "d", alpha=0.65, label="super stabilised_savgol_coeffs") + # ax.plot(c_s, 'x', alpha=0.65, label='stabilised_savgol_coeffs') + + ax.legend(shadow=True, framealpha=1) + ax.grid() + ax.set_title(f"Savitzky-Golay Coefficients\n{window_len=} {order=} {pos=}") + ax.set_xlabel("Coefficient index") + ax.set_ylabel("Coefficient value") + plt.show() From 3555e9bf3ce394b02dc1f323a6f25ce04535ab03 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 22 Jun 2024 20:23:00 +0200 Subject: [PATCH 2/9] enhances speed of stabilized Savitzky-Golay coefficient computation by avoiding the explicit computation of the pseudoinverse --- python/scipy/savgol-precision/savgol_work.py | 101 +++++++++++-------- 1 file changed, 59 insertions(+), 42 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index 62494f4..fe48311 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -1,4 +1,4 @@ -from math import comb, fsum +from math import comb from typing import Literal, Optional, Union import numpy as np @@ -9,7 +9,7 @@ import mpsig -def _pinv_correction_factors_normalized_polyvander( +def _pinv_correction_factors_stabilized_polyvander( pinv_row_idx: int, polyorder: int, polyvander_column_scales: np.ndarray, @@ -17,7 +17,7 @@ def _pinv_correction_factors_normalized_polyvander( x_scale: float, ) -> np.ndarray: """ - Computes the correction factors for the pseudo-inverse of a normalized polynomial + Computes the correction factors for the pseudo-inverse of a stabilized polynomial Vandermonde matrix where the x-variable was normalized as ``x_scaled = (x - x_center) / x_scale``. After this normalization, the columns of the Vandermonde matrix ``A`` given by ``A = polyvander(x_scaled, polyorder)`` are scaled to have unit Euclidean norm. @@ -75,35 +75,43 @@ def _pinv_correction_factors_normalized_polyvander( """ # noqa: E501 - # first, the signs are determined by the (-1)**(i - k) factor ... + # first, the signs are computed + # incorporating the signs as -1 and 1 is not efficient, especially not when they + # are computed by raising -1 to the power of the iterator + # therefore, the signs are already pre-multiplied with the only constant factor in + # the formula, namely ``x_center ** (-pinv_row_idx)`` + x_center_to_pinv_row_idx = x_center ** (-pinv_row_idx) + # the following is equivalent to a multiplication of the factor with the signs, + # however, since the signs are alternating, they are simply repeated as often as + # necessary + prefactors = [x_center_to_pinv_row_idx, -x_center_to_pinv_row_idx] + n_full_repetitions, n_rest_repetitions = divmod(polyorder + 1 - pinv_row_idx, 2) + prefactors = np.array( + prefactors * n_full_repetitions + prefactors[:n_rest_repetitions] + ) + + # then, the binomial coefficients are computed ... i_vect = np.arange(start=pinv_row_idx, stop=polyorder + 1, dtype=np.int64) - signs = (-1) ** (i_vect - pinv_row_idx) - # ... followed by the binomial coefficients ... binomials = np.array([comb(i, pinv_row_idx) for i in i_vect], dtype=np.int64) - # ... then the factors for the x-variable - x_factors = ( - (x_center ** (i_vect - pinv_row_idx)) - / (x_scale**i_vect) - / polyvander_column_scales[pinv_row_idx::] - ) + # ... followed by the x-factors + x_factors = ((x_center / x_scale) ** i_vect) / polyvander_column_scales[ + pinv_row_idx:: + ] - return signs * binomials * x_factors + return prefactors * binomials * x_factors -def _correct_pinv_row_normalized_polyvander( +def _get_corrected_pinv_row_from_stabilized_polyvander( pinv_row_idx: int, - pinv_scaled_polyvander: np.ndarray, + normalized_polyvander: np.ndarray, polyorder: int, polyvander_column_scales: np.ndarray, x_center: float, x_scale: float, ) -> np.ndarray: """ - Corrects the ``k``-th row of pseudo-inverse of a normalized polynomial Vandermonde - matrix to obtain the ``k``-th row of the pseudo-inverse of the polynomial fit based - on the original columns of the Vandermonde matrix. - For the correction it relies on floating point accurate summation to avoid loss of - precision due to cancellation errors with high polynomial orders. + Obtains the corrected ``k``-th row of the pseudo-inverse of a polynomial Vandermonde + matrix from a stabilized version of the latter. Please refer to the documentation of the function :func:`_pinv_correction_factors_normalized_polyvander` for a detailed explanation of the correction factors. @@ -112,8 +120,9 @@ def _correct_pinv_row_normalized_polyvander( pinv_row_idx : int The row index of the pseudo-inverse matrix to be corrected. This corresponds to ``k`` in the above equation. - pinv_scaled_polyvander : np.ndarray of shape (polyorder + 1, m) - The pseudo-inverse of the normalized polynomial Vandermonde matrix. + normalized_polyvander : np.ndarray of shape (m, polyorder + 1) + The normalized polynomial Vandermonde matrix. The columns of this matrix are + scaled to have unit Euclidean norm after the x-variable was normalized. polyorder : int The polynomial order of the polynomial fit. polyvander_column_scales : np.ndarray of shape (polyorder + 1,) @@ -123,25 +132,36 @@ def _correct_pinv_row_normalized_polyvander( Returns ------- - pinv_corrected_polyvander : np.ndarray of shape (polyorder + 1, m) - The corrected pseudo-inverse of the polynomial fit that corresponds to the - original columns of the polynomial Vandermonde matrix. + pinv_corrected_row : np.ndarray of shape (m,) + The corrected row of the pseudo-inverse of the polynomial fit that corresponds + to the original columns of the polynomial Vandermonde matrix. """ # noqa: E501 - correction_factors = _pinv_correction_factors_normalized_polyvander( + # first, the correction factors are computed + correction_factors = _pinv_correction_factors_stabilized_polyvander( pinv_row_idx=pinv_row_idx, polyorder=polyorder, polyvander_column_scales=polyvander_column_scales, x_center=x_center, x_scale=x_scale, ) - return np.array( - [ - fsum(pinv_scaled_polyvander[pinv_row_idx::, j] * correction_factors) - for j in range(pinv_scaled_polyvander.shape[1]) - ] - ) + + # then, the Singular Value Decomposition is obtained and prepared for the + # pseudo-inverse computation + u, s, vt = np.linalg.svd(normalized_polyvander, full_matrices=False) + rcond = np.finfo(normalized_polyvander.dtype).eps * max(normalized_polyvander.shape) + with np.errstate(divide="ignore", invalid="ignore"): + s = np.where(s > rcond * s.max(), np.reciprocal(s), 0.0) + + # finally, the pseudo-inversion is applied without explicitly forming the + # pseudo-inverse matrix + # NOTE: the pseudo-inverse is the one of the transposed Vandermonde matrix, hence + # its not V @ inv(S) @ U.T, but U @ inv(S) @ V.T + # NOTE: the order here is very important to avoid unnecessary operations by ensuring + # that only matrix-vector products are computed rather than matrix-matrix + # products + return u @ (s * (vt[::, pinv_row_idx::] @ correction_factors)) def super_stabilised_savgol_coeffs( @@ -189,23 +209,20 @@ def super_stabilised_savgol_coeffs( x_center = 0.5 * (x_min + x_max) x_scale = 0.5 * (x_max - x_min) x_normalized = (x - x_center) / x_scale - A_scaled_polyvander = np.polynomial.polynomial.polyvander(x_normalized, polyorder) + A_normalized_polyvander = np.polynomial.polynomial.polyvander( + x_normalized, polyorder + ) # to stabilize the pseudo-inverse computation, the columns of the polynomial # Vandermonde matrix are normalized to have unit Euclidean norm - a_column_scales = np.linalg.norm(A_scaled_polyvander, axis=0) - A_scaled_polyvander /= a_column_scales[np.newaxis, ::] - - # then, the pseudo-inverse is computed and corrected to obtain the pseudo-inverse - # of the polynomial fit based on the columns the original polynomial Vandermonde - # matrix based on the original x-values - A_pinv_distorted = np.linalg.pinv(A_scaled_polyvander) + a_column_scales = np.linalg.norm(A_normalized_polyvander, axis=0) + A_normalized_polyvander /= a_column_scales[np.newaxis, ::] return ( float_factorial(deriv) / delta**deriv - ) * _correct_pinv_row_normalized_polyvander( + ) * _get_corrected_pinv_row_from_stabilized_polyvander( pinv_row_idx=deriv, - pinv_scaled_polyvander=A_pinv_distorted, + normalized_polyvander=A_normalized_polyvander, polyorder=polyorder, polyvander_column_scales=a_column_scales, x_center=x_center, From de949bac0e3b41863b1413eb875b674114dfbb42 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sun, 23 Jun 2024 11:35:50 +0200 Subject: [PATCH 3/9] replaced custom SVD-logic for row correction of the pseudoinverse with `np.linalg.lstsq` --- python/scipy/savgol-precision/savgol_work.py | 31 +++++++++----------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index fe48311..64b29aa 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -69,9 +69,10 @@ def _pinv_correction_factors_stabilized_polyvander( Returns ------- - pinv_correction_factors : np.ndarray of shape (polyorder + 1 - pinv_row_idx,) + pinv_correction_factors : np.ndarray of shape (polyorder + 1) The correction factors for the pseudo-inverse matrix row corresponding to ``pinv_row_idx``. + All the elements that are not required are set to zero. """ # noqa: E501 @@ -98,7 +99,9 @@ def _pinv_correction_factors_stabilized_polyvander( pinv_row_idx:: ] - return prefactors * binomials * x_factors + pinv_correction_factors = np.zeros(shape=(polyorder + 1,), dtype=np.float64) + pinv_correction_factors[pinv_row_idx::] = prefactors * binomials * x_factors + return pinv_correction_factors def _get_corrected_pinv_row_from_stabilized_polyvander( @@ -147,21 +150,15 @@ def _get_corrected_pinv_row_from_stabilized_polyvander( x_scale=x_scale, ) - # then, the Singular Value Decomposition is obtained and prepared for the - # pseudo-inverse computation - u, s, vt = np.linalg.svd(normalized_polyvander, full_matrices=False) - rcond = np.finfo(normalized_polyvander.dtype).eps * max(normalized_polyvander.shape) - with np.errstate(divide="ignore", invalid="ignore"): - s = np.where(s > rcond * s.max(), np.reciprocal(s), 0.0) - - # finally, the pseudo-inversion is applied without explicitly forming the - # pseudo-inverse matrix - # NOTE: the pseudo-inverse is the one of the transposed Vandermonde matrix, hence - # its not V @ inv(S) @ U.T, but U @ inv(S) @ V.T - # NOTE: the order here is very important to avoid unnecessary operations by ensuring - # that only matrix-vector products are computed rather than matrix-matrix - # products - return u @ (s * (vt[::, pinv_row_idx::] @ correction_factors)) + # then, the corrected row is computed by solving the Least Squares problem + # normalized_polyvander.T @ pinv_corrected_row = correction_factors + pinv_corrected_row, _, _, _ = np.linalg.lstsq( + normalized_polyvander.transpose(), + correction_factors, + rcond=None, + ) + + return pinv_corrected_row def super_stabilised_savgol_coeffs( From 75bfd579c02a960d28c43c11a8921c7d4d33201b Mon Sep 17 00:00:00 2001 From: MothNik Date: Sun, 23 Jun 2024 12:10:28 +0200 Subject: [PATCH 4/9] - removed the intermediate function to apply the pseudoinverse correction for the super stabilised Savitzky Golay coefficients - in this context, the derivative pre-factor was incorporated into the correction factors as a single rather than multiple multiplications --- python/scipy/savgol-precision/savgol_work.py | 199 ++++++++----------- 1 file changed, 79 insertions(+), 120 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index 64b29aa..89dc257 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -9,63 +9,31 @@ import mpsig -def _pinv_correction_factors_stabilized_polyvander( - pinv_row_idx: int, +def _savgol_pinv_correction_factors( + deriv: int, polyorder: int, - polyvander_column_scales: np.ndarray, x_center: float, x_scale: float, + polyvander_column_scales: np.ndarray, + delta: float, ) -> np.ndarray: """ - Computes the correction factors for the pseudo-inverse of a stabilized polynomial - Vandermonde matrix where the x-variable was normalized as ``x_scaled = (x - x_center) / x_scale``. - After this normalization, the columns of the Vandermonde matrix ``A`` given by - ``A = polyvander(x_scaled, polyorder)`` are scaled to have unit Euclidean norm. - These steps ensure the numerical stability of the pseudo-inverse computation, but - distorts the pseudo-inverse of the polynomial fit. It will then be matched to the - scaled columns of ``A`` who were themselves normalized before, rather than the - original columns of ``A`` that are based on ``x``. - - This function computes the correction factors to be applied to the rows of the - pseudo-inverse of ``A`` to recover the pseudo-inverse of the polynomial fit based - on the original columns of ``A``. - Multiplying out all scaling and centering factors and rearranging terms gives the - following equation for the correction factors: - - ``` - phi_ik = ((-1)**(i - k)) * comb(i, k) * ((1.0 / x_scale)**i) * 1.0 / col_scale[i] * (x_center ** (i - k)) - ``` - - where - - - ``k`` is the row index of the pseudo-inverse matrix to be corrected (the ``k``-th - row corresponds to the ``k``-th coefficient of the polynomial fit, e.g., ``k=0`` - corresponds to the constant term ``a0``), - - ``i`` is the iterator for the polynomial order - - ``comb(i, k)`` is the binomial coefficient, and - - ``col_scale[i]`` is the Euclidean norm of the ``i``-th column of the normalized - polynomial Vandermonde matrix - - These can be applied to correct the element in the ``k``-th row and the ``j``-th - column of the distorted pseudo-inverse ``JD+`` as follows: - - ``` - J+[k, j] = sum(phi_ik * JD+[i, j] for i in range(k, polyorder + 1)) - ``` - - to obtain the corrected pseudo-inverse ``J+``. + Computes the correction factors for the pseudoinversion for the Savitzky Golay + coefficient computations. Please refer to the Notes section for more details. Parameters ---------- - pinv_row_idx : int - The row index of the pseudo-inverse matrix to be corrected. This corresponds to - ``k`` in the above equation. + deriv : int + The derivative order and also row index of the pseudo-inverse matrix to be + corrected. This corresponds to ``d`` in the Notes section. polyorder : int The polynomial order of the polynomial fit. polyvander_column_scales : np.ndarray of shape (polyorder + 1,) The Euclidean norms of the columns of the normalized polynomial Vandermonde matrix. x_center, x_scale : float The centering and scaling factors used to normalize the x-variable. + delta : float + The spacing between the grid points of the signal to be filtered. Returns ------- @@ -74,93 +42,77 @@ def _pinv_correction_factors_stabilized_polyvander( ``pinv_row_idx``. All the elements that are not required are set to zero. + Notes + ----- + The correction takes into account that + + - the x-variable for the Vandermonde matrix was normalized as + ``x_normalized = (x - x_center) / x_scale``, + - the columns of the Vandermonde matrix ``J`` given by + ``J = polyvander(x_normalized, polyorder)`` are scaled to have unit Euclidean norm + by dividing by the column scales ``polyvander_column_scales``, + - that the ``d``-th derivative order introduces a pre-factor of ``d! / delta**d`` + to the coefficients. + + While the first two steps ensure the numerical stability of the pseudoinverse + computation, they distort the pseudoinverse of the polynomial fit. The correction + factors are applied to the rows of the pseudo-inverse of ``J`` to recover the + pseudoinverse of the polynomial fit based on the original columns of ``J``. + + The correction factors are computed as follows: + + ``` + phi_id = (d! / delta**d) * ((-1)**(i - d)) * comb(i, d) * ((1.0 / x_scale)**i) * 1.0 / col_scale[i] * (x_center ** (i - d)) + ``` + + where + + - ``d`` is the derivative order and also the row index of the pseudo-inverse matrix + to be corrected (the ``d``-th row corresponds to the ``d``-th derivative of the + polynomial fit), + - ``i`` is the iterator for the polynomial order, + - ``comb(i, d)`` is the binomial coefficient, + - ``col_scale[i]`` is the Euclidean norm of the ``i``-th column of the normalized + polynomial Vandermonde matrix, and + - ``delta`` is the spacing between the grid points of the signal to be filtered. + + These can be applied to correct the element in the ``d``-th row and the ``j``-th + column of the distorted pseudo-inverse ``JD+`` as follows: + + ``` + J+[d, j] = sum(phi_id * JD+[i, j] for i in range(d, polyorder + 1)) + ``` + + to obtain the corrected pseudo-inverse ``J+``. + """ # noqa: E501 # first, the signs are computed # incorporating the signs as -1 and 1 is not efficient, especially not when they # are computed by raising -1 to the power of the iterator # therefore, the signs are already pre-multiplied with the only constant factor in - # the formula, namely ``x_center ** (-pinv_row_idx)`` - x_center_to_pinv_row_idx = x_center ** (-pinv_row_idx) + # the formula, namely ``deriv! * ((x_center * delta) ** (-deriv))`` + x_center_modified_for_deriv = float_factorial(deriv) * ((x_center * delta) ** (-deriv)) # the following is equivalent to a multiplication of the factor with the signs, # however, since the signs are alternating, they are simply repeated as often as # necessary - prefactors = [x_center_to_pinv_row_idx, -x_center_to_pinv_row_idx] - n_full_repetitions, n_rest_repetitions = divmod(polyorder + 1 - pinv_row_idx, 2) + prefactors = [x_center_modified_for_deriv, -x_center_modified_for_deriv] + n_full_repetitions, n_rest_repetitions = divmod(polyorder + 1 - deriv, 2) prefactors = np.array( prefactors * n_full_repetitions + prefactors[:n_rest_repetitions] ) # then, the binomial coefficients are computed ... - i_vect = np.arange(start=pinv_row_idx, stop=polyorder + 1, dtype=np.int64) - binomials = np.array([comb(i, pinv_row_idx) for i in i_vect], dtype=np.int64) + i_vect = np.arange(start=deriv, stop=polyorder + 1, dtype=np.int64) + binomials = np.array([comb(i, deriv) for i in i_vect], dtype=np.int64) # ... followed by the x-factors - x_factors = ((x_center / x_scale) ** i_vect) / polyvander_column_scales[ - pinv_row_idx:: - ] + x_factors = ((x_center / x_scale) ** i_vect) / polyvander_column_scales[deriv::] pinv_correction_factors = np.zeros(shape=(polyorder + 1,), dtype=np.float64) - pinv_correction_factors[pinv_row_idx::] = prefactors * binomials * x_factors + pinv_correction_factors[deriv::] = prefactors * binomials * x_factors return pinv_correction_factors -def _get_corrected_pinv_row_from_stabilized_polyvander( - pinv_row_idx: int, - normalized_polyvander: np.ndarray, - polyorder: int, - polyvander_column_scales: np.ndarray, - x_center: float, - x_scale: float, -) -> np.ndarray: - """ - Obtains the corrected ``k``-th row of the pseudo-inverse of a polynomial Vandermonde - matrix from a stabilized version of the latter. - Please refer to the documentation of the function :func:`_pinv_correction_factors_normalized_polyvander` - for a detailed explanation of the correction factors. - - Parameters - ---------- - pinv_row_idx : int - The row index of the pseudo-inverse matrix to be corrected. This corresponds to - ``k`` in the above equation. - normalized_polyvander : np.ndarray of shape (m, polyorder + 1) - The normalized polynomial Vandermonde matrix. The columns of this matrix are - scaled to have unit Euclidean norm after the x-variable was normalized. - polyorder : int - The polynomial order of the polynomial fit. - polyvander_column_scales : np.ndarray of shape (polyorder + 1,) - The Euclidean norms of the columns of the normalized polynomial Vandermonde matrix. - x_center, x_scale : float - The centering and scaling factors used to normalize the x-variable. - - Returns - ------- - pinv_corrected_row : np.ndarray of shape (m,) - The corrected row of the pseudo-inverse of the polynomial fit that corresponds - to the original columns of the polynomial Vandermonde matrix. - - """ # noqa: E501 - - # first, the correction factors are computed - correction_factors = _pinv_correction_factors_stabilized_polyvander( - pinv_row_idx=pinv_row_idx, - polyorder=polyorder, - polyvander_column_scales=polyvander_column_scales, - x_center=x_center, - x_scale=x_scale, - ) - - # then, the corrected row is computed by solving the Least Squares problem - # normalized_polyvander.T @ pinv_corrected_row = correction_factors - pinv_corrected_row, _, _, _ = np.linalg.lstsq( - normalized_polyvander.transpose(), - correction_factors, - rcond=None, - ) - - return pinv_corrected_row - - def super_stabilised_savgol_coeffs( window_length: int, polyorder: int, @@ -206,26 +158,33 @@ def super_stabilised_savgol_coeffs( x_center = 0.5 * (x_min + x_max) x_scale = 0.5 * (x_max - x_min) x_normalized = (x - x_center) / x_scale - A_normalized_polyvander = np.polynomial.polynomial.polyvander( + J_normalized_polyvander = np.polynomial.polynomial.polyvander( x_normalized, polyorder ) # to stabilize the pseudo-inverse computation, the columns of the polynomial # Vandermonde matrix are normalized to have unit Euclidean norm - a_column_scales = np.linalg.norm(A_normalized_polyvander, axis=0) - A_normalized_polyvander /= a_column_scales[np.newaxis, ::] - - return ( - float_factorial(deriv) / delta**deriv - ) * _get_corrected_pinv_row_from_stabilized_polyvander( - pinv_row_idx=deriv, - normalized_polyvander=A_normalized_polyvander, + j_column_scales = np.linalg.norm(J_normalized_polyvander, axis=0) + J_normalized_polyvander /= j_column_scales[np.newaxis, ::] + + # then, the correction factors for the subsequent least squares problem are computed + correction_factors = _savgol_pinv_correction_factors( + deriv=deriv, polyorder=polyorder, - polyvander_column_scales=a_column_scales, + polyvander_column_scales=j_column_scales, x_center=x_center, x_scale=x_scale, + delta=delta, ) + # finally, the coefficients are obtained from solving the least squares problem + # J.T @ coeffs = correction_factors + return np.linalg.lstsq( + J_normalized_polyvander.transpose(), + correction_factors, + rcond=None, + )[0] + def stabilised_savgol_coeffs( window_length, From 49a040ce50262ed0b79be0bcfbd8d4ba14486579 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sun, 23 Jun 2024 12:32:05 +0200 Subject: [PATCH 5/9] replaced SciPy `float_factorial` by `math.factorial` for the correction coefficients of the super stabilised Savitzky Golay coefficients. It does not provide any benefit and `math.factorial` is more direct. Yet, it's a micro-optimization. --- python/scipy/savgol-precision/savgol_work.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index 89dc257..d3136ca 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -1,10 +1,9 @@ -from math import comb +from math import comb, factorial from typing import Literal, Optional, Union import numpy as np from matplotlib import pyplot as plt from mpmath import mp -from scipy._lib._util import float_factorial import mpsig @@ -92,7 +91,7 @@ def _savgol_pinv_correction_factors( # are computed by raising -1 to the power of the iterator # therefore, the signs are already pre-multiplied with the only constant factor in # the formula, namely ``deriv! * ((x_center * delta) ** (-deriv))`` - x_center_modified_for_deriv = float_factorial(deriv) * ((x_center * delta) ** (-deriv)) + x_center_modified_for_deriv = factorial(deriv) * ((x_center * delta) ** (-deriv)) # the following is equivalent to a multiplication of the factor with the signs, # however, since the signs are alternating, they are simply repeated as often as # necessary From 1655e7886490757da7115a2ca6f268d8ee160083 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sun, 30 Jun 2024 17:10:15 +0200 Subject: [PATCH 6/9] made `pos` more type checker friendly --- python/scipy/savgol-precision/savgol_work.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index d3136ca..fa1daf4 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -126,13 +126,14 @@ def super_stabilised_savgol_coeffs( halflen, rem = divmod(window_length, 2) - if pos is None: + pos_internal = pos + if pos_internal is None: if rem == 0: - pos = halflen - 0.5 + pos_internal = halflen - 0.5 else: - pos = halflen + pos_internal = halflen - if not (0 <= pos < window_length): + if not (0 <= pos_internal < window_length): raise ValueError("pos must be nonnegative and less than " "window_length.") if use not in ["conv", "dot"]: @@ -146,7 +147,12 @@ def super_stabilised_savgol_coeffs( # from -pos to window_length - pos - 1. The powers (i.e., rows) range # from 0 to polyorder. (That is, A is a vandermonde matrix, but not # necessarily square.) - x = np.arange(-pos, window_length - pos, dtype=float) + x = np.arange( + start=-pos_internal, + stop=window_length - pos_internal, + step=1.0, + dtype=float, + ) if use == "conv": # Reverse so that result can be used in a convolution. From b795f4b028fa1dd764cb7c792f531c06e896f38d Mon Sep 17 00:00:00 2001 From: MothNik Date: Sun, 30 Jun 2024 18:49:07 +0200 Subject: [PATCH 7/9] fixed potential zero-division bug even though it's not clear whether this could ever be encountered in this context --- python/scipy/savgol-precision/savgol_work.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index fa1daf4..b05646a 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -170,6 +170,7 @@ def super_stabilised_savgol_coeffs( # to stabilize the pseudo-inverse computation, the columns of the polynomial # Vandermonde matrix are normalized to have unit Euclidean norm j_column_scales = np.linalg.norm(J_normalized_polyvander, axis=0) + j_column_scales[j_column_scales == 0.0] = 1.0 J_normalized_polyvander /= j_column_scales[np.newaxis, ::] # then, the correction factors for the subsequent least squares problem are computed From 3c175a04a78d5672b4fc2b073b7630ca3209f4ea Mon Sep 17 00:00:00 2001 From: MothNik Date: Mon, 1 Jul 2024 23:05:34 +0200 Subject: [PATCH 8/9] - fixed accidentally untested failing case of ``pos = window_length // 2`` in super stabilised Savitzky Golay coefficient computation by using a simplified solution in this case - made relative error test in main truly insensitive to ZeroDivision (`inf` was never used in case of zero division because Python does not evaluate lazily) - added a more extensive test for the coefficients against the correct solution to main to avoid thing like the first bullet point happening --- python/scipy/savgol-precision/savgol_work.py | 71 ++++++++++++++------ 1 file changed, 52 insertions(+), 19 deletions(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index b05646a..3f8c42d 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -8,7 +8,7 @@ import mpsig -def _savgol_pinv_correction_factors( +def _savgol_pinv_correction_factors_with_x_center( deriv: int, polyorder: int, x_center: float, @@ -18,7 +18,8 @@ def _savgol_pinv_correction_factors( ) -> np.ndarray: """ Computes the correction factors for the pseudoinversion for the Savitzky Golay - coefficient computations. Please refer to the Notes section for more details. + coefficient computations when the filter is not centered within its window, i.e., + ``x_center != 0``. Please refer to the Notes section for more details. Parameters ---------- @@ -173,15 +174,27 @@ def super_stabilised_savgol_coeffs( j_column_scales[j_column_scales == 0.0] = 1.0 J_normalized_polyvander /= j_column_scales[np.newaxis, ::] - # then, the correction factors for the subsequent least squares problem are computed - correction_factors = _savgol_pinv_correction_factors( - deriv=deriv, - polyorder=polyorder, - polyvander_column_scales=j_column_scales, - x_center=x_center, - x_scale=x_scale, - delta=delta, - ) + # of ``x_center != 0``, the correction factors for the subsequent least squares + # problem are a bit more involved to compute + if x_center != 0.0: + correction_factors = _savgol_pinv_correction_factors_with_x_center( + deriv=deriv, + polyorder=polyorder, + polyvander_column_scales=j_column_scales, + x_center=x_center, + x_scale=x_scale, + delta=delta, + ) + + # for ``x_center == 0``, the correction factors are highly simplified + else: + # in this case, the correction factors just a column vector whose ``d``-th + # entry is given by + # ``(d! / delta**d) * 1.0 / col_scale[d] * 1.0 / (x_scale ** d)`` + correction_factors = np.zeros(shape=(polyorder + 1,), dtype=np.float64) + correction_factors[deriv] = factorial(deriv) / ( + j_column_scales[deriv] * ((delta * x_scale) ** deriv) + ) # finally, the coefficients are obtained from solving the least squares problem # J.T @ coeffs = correction_factors @@ -271,21 +284,22 @@ def simple_savgol_coeffs(window_length, polyorder, pos=None): def check_savgol_coeffs(c, cref): - relerr = np.array( - [ - float(abs((c0 - cref0) / cref0)) if cref != 0 else np.inf - for c0, cref0 in zip(c, cref) - ] - ) + relerr = np.empty_like(c) + for i, (c0, cref0) in enumerate(zip(c, cref)): + if cref0 != 0: + relerr[i] = abs((c0 - cref0) / cref0) + else: + relerr[i] = np.nan + return relerr if __name__ == "__main__": mp.dps = 150 - window_len = 51 + window_len = 55 order = 10 - pos = 35 + pos = 27 deriv = 5 delta = 0.5 @@ -343,3 +357,22 @@ def check_savgol_coeffs(c, cref): ax.set_ylabel("Coefficient value") plt.show() + + +# Long Test section + +if __name__ == "__main__" and True: + + for order in [0, 1, 9, 10]: # odd and even + for deriv in [0, 1, 2, 3, 4, 5, 6]: + for window_len in [11, 31]: # to keep computation time reasonable + for pos in range(0, window_len): # all possible positions + print(f"Checking {window_len=} {order=} {deriv=} {pos=}") + coeffs = mpsig.savgol_coeffs( + window_len, order, pos=pos, deriv=deriv + ) + c_sup_s = super_stabilised_savgol_coeffs( + window_len, order, pos=pos, deriv=deriv + ) + + assert np.allclose(coeffs, c_sup_s) \ No newline at end of file From 20042c628b2f852ab2c6b4202e4a6ae3fde8af7e Mon Sep 17 00:00:00 2001 From: MothNik Date: Mon, 1 Jul 2024 23:20:02 +0200 Subject: [PATCH 9/9] made `pos == window_length // 2` check for super stabilised Savitzky-Golay coefficients safe against floating point errors in the zero-check --- python/scipy/savgol-precision/savgol_work.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/python/scipy/savgol-precision/savgol_work.py b/python/scipy/savgol-precision/savgol_work.py index 3f8c42d..4057463 100644 --- a/python/scipy/savgol-precision/savgol_work.py +++ b/python/scipy/savgol-precision/savgol_work.py @@ -176,7 +176,11 @@ def super_stabilised_savgol_coeffs( # of ``x_center != 0``, the correction factors for the subsequent least squares # problem are a bit more involved to compute - if x_center != 0.0: + # NOTE: this is not checked via ``x_center`` which can suffer from floating point + # errors, but via the internal variable ``pos_internal`` which is ensured to + # be an integer for the cases where ``x_center == 0`` (so this is a safe + # check) + if pos_internal != halflen: correction_factors = _savgol_pinv_correction_factors_with_x_center( deriv=deriv, polyorder=polyorder,