Research Program on Splines
Spline Methods for the Continuous/Discrete Pocessing of Images
SNSF Grant: 200020_184646 / 1Principal Investigator: Prof. Michael Unser
Participants: Pakshal Borha, Shayan Aziznejad, Thomas Debarre, Alexis Goujon, Stanislas Ducotterd, Joaquim Campos, Mehrsa Pourya, Rahul Parhi
Period: July 1, 2019–October 31, 2023.
Research Digest
The work supported by this grant is part of a longterm research program (1998-present) that aims at demonstrating the relevance and usefulness of splines to a number of disciplines, including signal and image processing, statistics and stochastic modeling, mathematics (functional analysis and optimization theory) and, more recently, machine learning. The main effort during the past term has been directed toward the specification of splines as the minimizers of some sparsity-promoting functional and the benefit of this property in the derivation of novel algorithms for signal processing, machine learning, and the resolution of inverse problems.
A. Splines as Minimizers of Sparsity-Promoting Functionals
Splines can often be specified as the minimizers of a loss functional that consists of a data term (which measures the goodness of fit to the available data) plus a regularization that involves a defining operator L and a suitable norm. In the traditional setting, the norm is quadratic and the search space (a.k.a. native space) is a reproducing-kernel Hilbert space (RKHS). We have generalized the concept by switching to more general norms—in particular, the ℳ norm (a.k.a. total variation in measure theory), which promotes sparsity and results in a generalized total-variation (gTV) criterion. Our major finding is that, under very mild hypotheses, the use of gTV regularization induces solutions that are nonuniform L-splines with a few adaptive knots. This is true for a broad category of problems such as data fitting (interpolation), supervised learning, and the resolution of linear inverse problems [Unser 2021 FCOM].
The crucial ingredient for the deployment of our theory (representer theorem) is the completeness (Banach property) of the corresponding native space. We have shown that it hinges upon the construction of a generalized inverse of L [Unser 2019 arXiv:1904.10818]. The critical point here is to ascertain the stability of this generalized inverse, which presently still needs to be done on a case-by-case basis. So far, we have treated the "easy" case where the null space of L is trivial—that is, when the inverse operator corresponds to a classical convolution operator [Unser 2020 IEEE TSP]—and a few other specific instances such as L = D2 (the second derivative) [Appendix B, Unser 2019 JMLR], which is especially relevant for machine learning.
We have also investigated a new regularizer: the ℳ-norm of the pointwise Shatten-1 (or nuclear) norm of the Hessian, denoted by HTV(ƒ) (Hessian total variation) [Aziznejad 2021 NFAO], which we claim to be the proper multidimensional generalization of 2nd-order TV. We could show that, similarly to 2nd-order TV, HTV promotes solutions that are continuous piecewise-linear (CPWL)—the multidimensional analogs of linear splines—which also supports its use as a measure of complexity for neural networks [Aziznejad 2023 SIAM JMDS].
B. Multi-Splines
The traditional representation of cardinal splines involves basis functions (B-splines) that are shifted replicates of a single generator. As natural extension, we investigated multi-splines, which are (direct) sums of ordinary splines. While such multi-splines can be expressed in terms of conventional B-splines, we identified alternative representations that rely on shorter basis functions. In particular, we proposed a mechanism to construct the shortest basis functions for direct sums of splines of successive orders [Goujon 2021 JCAM]. We established the connection between such representations and generalized-sampling techniques, with a special emphasis on derivative sampling, which goes hand-in-hand with Hermite splines [Fageot 2020 JCAM].
The prefix "multi" may also refer to the multidimensional domain of a function, which opens up a new range of possibilities. Of special relevance to machine learning are the multidimensional splines that are CPWL as these are precisely the functions encoded by deep ReLU networks or regularized fits with a HTV regularization. As one may expect, CPWL functions also have a more direct B-spline-like representation in terms of "hat" or nodal basis functions specified on some triangulation. We have characterized the stability properties of such explicit representations, while providing sharp estimates of their Riesz bounds [Goujon 2023 ACHA]. We have also shown how the CPWL nodal functions could be expressed in terms of ReLU-like atoms.
C. Splines and Machine Learning
We like to view the well-documented link between splines and RKHS (de Boor 1966) as the precursor of Scholkopf's famous representer theorem, which supports the use of kernel estimators in classic machine learning. We have now managed to go one step farther by replacing the RHKS regularization by suitable forms of gTV. The foundation for most of the results in this section is our extended representer theorem for abstract Banach spaces [Unser 2021 FCOM] and its variant for seminorm regularizations such as gTV [Unser 2022 ACHA].
Our first proposal is a sparse alternative to (multi-)kernel methods. It yields a representation with fewer atoms (kernels) that are positioned adaptively [Aziznejad 2023 SIAM JMDA]. For the univariate setting (d = 1), we investigated the use of 2nd-order TV regularization to find the best spline approximation of a series of data points with the minimal number of linear segments, and designed a very effective algorithm to find the sparsest solution [Debarre 2022 JCAM]. We could also prove that the inclusion of an additional stability constraint on ƒ (max Lipschitz constant) does not affect the linear-spline character of the solution [Aziznejad 2022 OJSP]. To extend the approach to higher dimensions, we investigated an isotropic variant of 2nd-order TV that involves the Laplacian and the Radon transform of computed tomography, the harder part being the proper characterization of the corresponding native space [Unser 2023 JMLR]. This allowed us to prove that the extremal points (solutions) of the corresponding training problem are parameterized by a ReLU network with one hidden layer, which makes an interesting connection with neural networks.
We found 2nd-order TV regularization to be very helpful in the learning of free-form activations in deep neural networks, as formalized in our deep-spline framework [Unser 2019 JMLR]. This can be done very efficiently with the help of linear B-splines as substitutes for ReLUs [Bohra 2020 OJSP]. We have proposed a variant of the scheme that offers stability control via the enforcement of Lipschitz constraints [Aziznejad 2020 IEEE TSP]. We could also prove that our deep-spline framework is universal in the class of stable neural networks with a Lipchitz-1 constraint on each layer [Neumayer 2023 SIAM JMDS].
As alternative to the neural-network parameterization of CPWLs, we investigated a functional approach to supervised learning that solely relies of HTV regularization. An important finding is that HTV lends itself to an exact discretization in a CWPL basis [Aziznejad 2023 SIAM JMDS]. This allowed us to develop some efficient training algorithms: first for d = 2 with the help of hexagonal box splines [Campos 2022 OJSP]; and then for intermediate dimensions up to d = 7 using the hat functions induced by a Delaunay triangulation [Pourya 2023 OJSP].
D. Splines and Inverse Problems
We are the initiators of a global approach to the resolution of linear inverse problems that bypasses (somewhat arbitrary) discretization steps and directly formulates the recovery problem in the continuum [Unser 2021 FCOM]. While our formulation involves an infinite-dimensional optimization (because functions have an infinite number of degrees of freedom), we could prove that the solution (optimum) generally exists and that, under suitable conditions, it has an explicit parameterization that involves a finite number of atoms (basis elements) that are matched to the regularization (semi-)norm. Another way to put it is that the functional optimization also yields the optimal discretization.
We have applied our paradigm to the reconstruction of 1D signals, using various forms of continuous-domain regularization. In particular, we developed exact reconstruction schemes for scenarios where our signal results from the superposition of two components, each of which being constrained by its own regularization. The models of interest are: (i) hybrid-spline dictionaries where each component has his own gTV with distinct regularization operator [Debarre 2019 IEEE TSP], and (ii) the sparse-plus-smooth model that involves a mix of gTV and Hibertian regularization [Debarre 2021 OJSP]. While continuous-domain optimization problems with gTV regularization generally admit multiple solutions, we have identified configurations involving Fourier domain-measurements, for which the solution is guaranteed to be unique [Debarre 2022 IP], [Debarre 2023 JCAM].
Beside gTV, we investigated the use of other regularizers, such as HTV for the interpolation of data in scanning x-ray transmission microscopy [Debarre ISBI 2020]. Another intriguing option is a spline energy with an Lp-norm with p > 0, for which we were able to derive an exact B-spline-based discretization [Bohra 2020 IEEE TSP]. This allowed us to compare the "best" interpolants in the Lp sense, including cases with p in [1, 0) [Bohra ICASSP 2020].
Last but not least, we have relied on spline modeling (deep-spline framework) to learn convex regularizers for the resolution of inverse problems [Goujon 2023 IEEE TCI]. These data-driven regularizers yield image reconstructions that compare favorably with the classic techniques, including TV.
E. Splines and Stochastic Processes
In previous works, we have established a mathematical connection between splines and stochastic processes. The main idea is that both entities can be expressed as solutions of the differential equation s = L w with a driving term w (the so-called innovation). In the case of deterministic splines, w is a weighted sum of Dirac impulses positioned on the knots, while in the case of stochastic processes, w is a continuous-domain white noise (not necessarily Gaussian).
Based on the property that any white noise can be described as a limit of impulsive noises (i.e., some random linear combination of Dirac impulses which are distributed uniformly over the line), we have proposed a general generation mechanism that exploits the equivalence between the stochastic process s and a random L-spline with sufficiently dense knots [Dadi 2020 IEEE TSP].
On the theoretical front, we have characterized the complete family of complex-valued differential operators L that are both shift- and scale-invariant. These operators are the complex-order extension of the fractional derivatives. They have allowed us to specify the complete (non-Gaussian) family of self-similar processes [Amini arXiv:2206.10688]. We have also investigated the wavelet transform of compound-Poisson processes to get a better handle on their Besov/Hölder regularity [Aziznejad 2022 IEEE IT].
We have considered the problem of the optimal recovery of Lévy processes from linear measurements and have developed a Monte-Carlo sampling procedure for the determination of the MMSE solution [Bohra 2023 IEEE TSP]. We have then used the latter—the best possible reconstruction—as gold-standard for the comparison of signal-reconstruction methods, including approaches that are based on deep neural networks.
F. Spline Methods for Image Analysis
As complement to our initial research objectives, we have improved some of our earlier methods for image analysis. In particular, we have developed a general framework for steerable filtering that relies on splines for the effective encoding of the radial frequency of the component filters [Fageot 2021 IEEE TIP]. We have implemented our scheme as the Steer'n'Detect plugin for ImageJ, which provides a fast and accurate method for the detection of a pattern of interest at any orientation in a 2D image [Uhlmann 2022 BI].
We have designed an efficient scheme to trace and eventually simplify/stylize the outline of an object [Jover 2022 IEEE TIP]. Our method performs a regularized spline fit of a series of points in the plane. It uses a variant of gTV that is specifically designed for parametric curves and is invariant to rotations of the coordinate system.