Non-polynomial splines, a non-linear problem

Question

I'm looking for references on how to construct spline-like functions from a basis that does not include piecewise polynomials.

To be specific, given a class of functions such as "decaying exponentials" or "sines and cosines" (which are parameterized by a single parameter, e.g. the decay rate or the frequency), is there an efficient and numerically stable method to construct a function that is piecewise a linear combination of $N$ such functions (whose parameters are to be determined), interpolates given data $(x_k,f_k)$ [which are assumed to be such that they can be interpolated using the given function class, i.e. e.g. monotonically decreasing for decaying exponentials] and has continuous derivatives at $x_k$ up to order $2N-2$ (in order to fix $N$ linear coefficients and $N$ non-linear parameters)?

I can of course write down the explicit equations needed to satisfy these conditions, but directly solving those using a non-linear equation system solver does not look too promising as an approach.

What literature I could find so far on splines with non-polynomial components referred to spaces spanned by polynomials and some given non-polynomial. Here I'm looking for the case where there are no polynomials and the non-polynomial functions are parameterized by a parameter whose values are to be determined by the interpolation and smoothness conditions.

Answer 1

This is not a full answer, but I guess it is more than a comment.

One way to reduce the task is to apply the spirit of VARPRO to separate the linear coefficients and the non-linear parameters, i.e. alternatingly solve $N$ of the conditions as a linear equation system for the coefficients of the solution in a given $N$-dimensional subspace, and solve the remaining $N$ conditions as a non-linear system for $N$ parameters giving a different subspace (in which the coefficients, being kept fixed, are those of the solution to those other $N$ conditions). This reduces the problem size by a factor of two.

Answer 2

The interpolant
$$A\left(e^{\frac{a}{A}(x-x_i)}-1\right)+B\cdot\left(\cosh\left(\frac{a}{A}(x-x_i\right)-1\right)+y_i\ =\ (A+B)\mathbf{e^{\frac{a}{A}(x-x_i)}}+B\cdot \mathbf{e^{-\frac{a}{A}(x-x_i)}}-(A+B)+y_i\\ a=\frac{d}{dx}S(x_i),\\ \Delta x_i=x_{i+1}-x_i,\\ \Delta y_i=y_{i+1}-y_i,\\ B=\frac{\Delta y_i-A\left(e^{\frac{a}{A}\Delta x_i}-1\right)}{\cosh\left(\frac{a}{A}\Delta x_i\right)-1}$$ solves the $C^1$ interpolation problem albeit not the periodic one.

The construction of the spline $S(x)$ happens from left to right and requires knowldege of the slope at $x_0$ to be able to calculate the Interpolant that connects $(x_0,y_0)$ with $(x_1,y_1)$; knowing the interpolant we can determine the slope at $x_1$ and we are in the same situation as before so that eventually $S(x)$ is determined.

The underlying ideas that led to identifying the interpolant for $C^1$ continuous interpolation can most likely be generalized to higher degrees of smoothness but I'm still on my way.

Addendum:
parameter $A$ provides limited control over the shape of the interpolant:

• if we want the "purest" expontial functions, the $A$ should be chosen to minimize $B$
• if we strive for preserving shape we can chose $A$ to control the abscissa of the (unique) local extremum in cases of $y_i\lt y_{i+1}\land y_{i+2}\lt y_{i+1}$
  • align it with the vertex of the parabola through $\left(\,(x_i,y_i),\,(x_{i+1},y_{i+1}),\,(x_{i+2},y_{i+2})\,\right)$
  • align it with the abscissa of the intersection of the line with slope $a$ that contains $(x_i,y_i)$ and the line through $(x_{i+1},y_{i+1})$ and $(x_{i+2},y_{i+2})$

Answer 3

if the functions to be "splined" are invertible on the intervals $[x_i,x_{i+1}]$ and, if $f_{i-1}(x)$ and $f_i(x)$ are two such functions that are invertible on $[x_{i-1},x_i]$, resp. on $[x_i,x_{i+1}]$, additionally $f_{i-1}^{-1}(x_i)=f_i^{-1}(x_i)$
then $$S(x)\ :=\ \mathcal{F}^{-1}\Big(\mathcal{S}\big(x_i,\mathcal{F}\left(y_i+\gamma(x_i)\,\right)\big)\Big)-\gamma(x) $$ yields an entire class of interpolating non-polynomial splines. I used the "strange" notation in order to stress the analogy to signal processing: if the original interpolation task is non-polynomial in the "time"-domain then "transforming" the $y_i$ via the inverse functions makes the interpolation polynomial in the "frequency" domain; interpreting interpolation as filtering completes the analogy to signal processing.
The key concept is Homomorphism and was introduced to signal processing in R.W. Schafer "Echo removal by discrete generalized linear filtering". Res. Lab. Electron. MIT, Tech. Rep. (1969)

The $\gamma(x)$ is a way to allow for mixed interpolation, e.g. exponentials plus polynomials; the case $\gamma(x)\equiv \mathrm{const}$ seem interesting in its own right when investigating the spline's limit behavior for e.g. $\gamma(x)\equiv \mathrm{const}\,\to\,\infty$, especially in the case of exponential splines.

Maybe its worth mentioning albeit trivial that this "homomorphic" splining directly carries over higher dimensions or parametric interpolation and one isn't limited to polynomials in the "frequency" domain; rational interpolants would be the "next bigger thing" to use.

Addendum:
the notation I used is aimed at emphasizing the analogy to transformation from the time domain to the frequency domain e.g. via a Fourier transform $\mathcal{F}$; the analogy to interpolating $(x_i,y_i)$ with exponential splines $S(x)\in C^{k-1}$ of the form $e^{\sum_{j=0}^k a_{ij}x^j}, S(x_i)=y_i$ is that the tranformation $\mathcal{F}$ corresponds to applying the analogy of the Fourier transformation $\mathcal{F}$, i.e. the inverse $\ln(\cdot)$ of $e^{(\cdot)}$, to the ordinates $y_i$ of the data to be interpolated and then calculate the interpolating polynomial spline with these transformed ordinates; that spline is denoted by $\mathcal{S}(x)$ because it is calculated in the analogy of the frequency domain.
the final step is to go back to the analogy of the time domain by applying the analogy $\mathcal{F}^{-1}$, i.e. $e^{(\cdot)}$ of the inverse Fourier transformation to $\mathcal{S}(x)$ to obtain the interpolating spline $S(x)$ in the analogy of the time domain.

The $\gamma(x)$ allow for modeling situations where the data is "almost" polynomial, i.e. if for some order $h$ the magnitudes of the divided differences of order $h$ become small, then it may, depending on the model of the origin of the data, improve the quality of the interpolating spline if the values of a "lower polynomial hull" at the abscissas $x_i$ are subtracted from the ordinates prior to the transition to the analogy of the frequency domain.
Another usecase for introducing $\gamma(x)$ is if one is only interested in the shape of the interpolating curve but the ordinates of the data are outside the value range of the interpolation curve, e.g. negative $y_i$ in the case of exponential splines; then adding $\gamma(x)$ can fix these issues.

An example may make things clearer:
assuming we have sampled an empirical distribution of which we assume that is the sum of Gaussian's that we would like to recover. In that case an Ansatz could be to take the logarithm of the ordinates, calculate a quadratic spline and take the abscissas of apices of parabolas with negative leading coefficient as an initial guess for the central values of the Gaussians.