CSSD — Cubic smoothing splines for discontinuous signals

Cubic smoothing splines for signals with a-priori unknown discontinuities — joint estimation of the spline and the discontinuity set via dynamic programming, with K-fold cross-validation for the regularisation parameters.

This is a reference implementation for the algorithms described in the paper

M. Storath, A. Weinmann, "Smoothing splines for discontinuous signals", Journal of Computational and Graphical Statistics, 2023, [Preprint]

Overview of main functionalities

1. cssd.m computes a cubic smoothing spline with discontinuities (CSSD) for data (x,y). It is a solution of the following model of a smoothing spline $f$ with a-priori unknown discontinuities $J$

$$\min_{f, J} p \sum_{i=1}^N \left(\frac{y_i - f(x_i)}{\delta_i}\right)^2 + (1-p) \int_{[x_1, x_N] \setminus J} (f''(t))^2 dt + \gamma |J|.$$

where

• $y_i = g(x_i) + \epsilon_i$ are samples of piecewise smooth function $g$ at data sites $x_1, \ldots, x_N$, and an estimate $\delta_i$ of the standard deviation of the errors $\epsilon_i$
• the minimum is taken over all possible sets of discontinuities between two data sites $J \subset [x_1, x_N]\setminus {x_1, \ldots, x_N}$ and all functions $f$ that are twice continuously differentiable away from the discontinuities.
• The model parameter $p \in (0, 1)$ controls the relative weight of the smoothness term (second term) and the data fidelity term.
• The last term is a penalty for the number of discontinuities $|J|$ weighted by a parameter $\gamma > 0.$

2. cssd_cv.m automatically determines values for the model parameters $p$ and $\gamma$ based on K-fold cross validation.

Quickstart

Python (Rust core)

pip install cssd

import numpy as np
from cssd import cssd, cssd_cv

x = np.linspace(0, 1, 100)
y = np.sin(4 * np.pi * x) - np.sign(x - 0.3) - np.sign(0.72 - x)

out = cssd(x, y, p=0.999, gamma=8.0)
out.discont      # detected jump locations
out.pp(x)        # evaluate the piecewise spline

cv = cssd_cv(x, y, cv_type="random", cv_arg=5)
cv.p, cv.gamma, cv.fit.discont

The Python package wraps a Rust extension built with PyO3 and maturin; see crates/cssd-core for the algorithm crate and crates/cssd-py for the bindings.

MATLAB (reference implementation, unchanged)

1. Execute "install_cssd.m" which adds the folder and all subfolders to the Matlab path.
2. Execute any m-file from the demos folder

Examples

Synthetic data

Stock data

Geyser data

Reference

M. Storath, A. Weinmann, "Smoothing splines for discontinuous signals", Journal of Computational and Graphical Statistics, 2023

How to cite

If you use this software, please cite the paper above. GitHub's "Cite this repository" button on the repo page reads the version and date-released fields from CITATION.cff and renders BibTeX/APA.

See also

Sibling projects from the same research program on variational methods for signal and image processing:

• Pottslab — multilabel image segmentation via the Potts / piecewise-constant Mumford-Shah model
• L1TV — exact L1-TV regularisation of real- or circle-valued signals
• MumfordShah2D — edge-preserving image restoration via the Mumford-Shah model
• CircleMedianFilter — fast median filtering for phase or orientation data
• DCEBE — bolus arrival time estimation for DCE-MRI signals

Related external projects:

• Higher order Mumford-Shah models — discrete splines with discontinuities.
• pcw-regrs — degrees-of-freedom penalised piecewise regression.

License

Released under the MIT License. See LICENSE.