Charles Hermite (1822 - 1901)

Sunday morning in Paris was opportunity to walk to the Montparnasse cemetery and pay respect to some of my personal heroes buried there. While the graves of Beckett, Cortázar and Poincaré attract quite many attention and visitors, it is less known that the great French mathematician Charles Hermite is buried there as well. Not only that he was a mentor to Henri Poincaré, Henri Padé, Thomas Stieltjes and Mihajlo Petrovic Alas, but his work on interpolation and function approximation is at the very core of the modern numerical methods used in computational fluid dynamics and radiative transfer (even when this is not so obvious or properly acknowledged). His work from 1878 ("Sur la formule d'interpolation de Lagrange") should be read by anyone interested in function approximation. The name on his grave stone are barely readable nowadays.

Quadratic splines: Part I

Some basic properties of parabola in the context of the second-order Hermite interpolation equation. 

There is a one-dimensional dataset with $n$ points: $x = (x_0, \dots, x_{n-1})$ and $y = (y_0, \dots, y_{n-1})$, such that $x_{i} > x_{i-1}$ for any $i$. We introduce the following labels: $$h_{i-1} = x_i - x_{i-1}$$ $$\delta_{i-1} = \frac{y_i - y_{i-1}}{x_i - x_{i-1}}$$ $h$ is always positive, while $\delta$ can have any value. On each segment, we approximate the unknown function with a parabola. On the $[i-1, i]$ segment: \begin{equation}f(x) = a_0 + a_1 (x - x_{i-1}) + a_2 (x - x_{i-1})^2\end{equation} The first and the second derivative of the parabola are: \begin{equation}f'(x) = a_1 + 2 a_2 (x - x_{i-1})\end{equation} \begin{equation}f''(x) = 2 a_2 \end{equation} The Hermite interpolation theorem provides us with a tool to determine the $a$ coefficients on the $[i-1, i]$ segment if we know $y_{i-1}$, $y_i$ and one of the first derivative in the end points, either $y'_{i-1}$ or $y'_i$. Let's assume that we know ("know" in the sense "we can compute") $y'_{i-1}$. Then there are three equations that the parabola must obey: \begin{equation}f(x_{i-1}) = y_{i-1} = a_0\end{equation} \begin{equation}f'(x_{i-1}) = y'_{i-1} = a_1\end{equation} \begin{equation}f(x_{i}) = y_i = a_0 + a_1 h_{i-1} + a_2 h_{i-1}^2\end{equation} The solution is then trivial: \begin{equation}a_0 = y_{i-1}\end{equation} \begin{equation}a_1 = y'_{i-1}\end{equation} \begin{equation}a_2 = \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2}\end{equation} Once the $a$ coefficients are known, the parabola is fully determined. That also means that the derivative in the point $i$ is fixed: $$f'(x_{i}) = y'_i = a_1 + 2 a_2 (x_i - x_{i-1}) = y'_{i-1} + 2 \frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} h_{i-1}$$ $$y'_i = \frac{2(y_{i} - y_{i-1}) - y'_{i-1} h_{i-1}}{h_{i-1}} $$ \begin{equation}y'_i = 2\delta_{i-1} - y'_{i-1}\end{equation} Local extremum

The extremum of this parabola is located at $x_m$: $$f'(x_m) = 0 = a_1 + 2 a_2 (x_m - x_{i-1})$$ \begin{equation}x_m = x_{i-1} - \frac{a_1}{2a_2}\end{equation} Monotonicity

The function $f(x)$ is monotone on the $[i-1, i]$ segment if it does not have an extremum on it, i.e. if $$x_m < x_{i-1}\;\;\;\;\;\;\;\;\;\;\mathrm{or}\;\;\;\;\;\;\;\;\;\; x_m > x_i$$ The first condition readily translates into: $$x_{i-1} - \frac{a_1}{2a_2} < x_{i-1}$$ $$ \frac{a_1}{2a_2} > 0$$ or, when we substitute $a_1$ and $a_2$: $$ \frac{y'_{i-1} h_{i-1}^2}{2(y_{i} - y_{i-1} - y'_{i-1} h_{i-1})} > 0$$ $$ \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$ As by definition $h_{i-1} > 0$ $$ \frac{y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} > 0$$ If $y'_{i-1} > 0$ then $\delta_{i-1} - y'_{i-1} > 0$ and $\delta_{i-1} > y'_{i-1}$. If $y'_{i-1} < 0$, then $\delta_{i-1} < y'_{i-1}$. Therefore, the first condition ($x_m < x_{i-1}$) reduces to: \begin{equation}\boxed{|\delta_{i-1}| > |y'_{i-1}|}\end{equation} The second condition leads to: $$x_{i-1} - \frac{a_1}{2a_2} > x_{i}$$ $$-x_{i} + x_{i-1} - \frac{a_1}{2a_2} > 0$$ $$h_{i-1} + \frac{a_1}{2a_2} < 0$$ $$h_{i-1} + \frac{y'_{i-1} h_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$ $$\frac{2\delta_{i-1} - y'_{i-1}}{2(\delta_{i-1} - y'_{i-1})} < 0$$ There are two cases. First, if $2\delta_{i-1} - y'_{i-1} > 0$, i.e. $2\delta_{i-1} > y'_{i-1}$. then it must be $\delta_{i-1} - y'_{i-1}< 0$, i.e.$\delta_{i-1} < y'_{i-1}$. Therefore, the first case leads to $2\delta_{i-1} > y'_{i-1} > \delta_{i-1}$. Secondly, if $2\delta_{i-1} - y'_{i-1} < 0$, i.e. $2\delta_{i-1} < y'_{i-1}$, it must be $\delta_{i-1} - y'_{i-1}> 0$, i.e. $\delta_{i-1} > y'_{i-1}$. So the second case leades to $\delta_{i-1} > y'_{i-1}> 2\delta_{i-1}$. The first case is possible only when both $\delta_{i-1}$ and $y'_{i-1}$ are positive and the second case is possible only when both are negative. Therefore, together, the two cases reduce the second condition ($x_m > x_{i}$) to: \begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}| > |\delta_{i-1}|}\end{equation} Together, the two conditions combine into \begin{equation}\boxed{|2\delta_{i-1}| > |y'_{i-1}|}\end{equation} A much faster way to reach the same result is to realize that the the derivatives $y'_{i-1}$ and $y'_{i}$ must be of the same sign if the function is monotone on the $[i-1, i]$ segment, i.e. $$y'_{i-1} \cdot (2\delta_{i-1} - y'_{i-1}) > 0$$ which obviously leads to $|2\delta_{i-1}| > |y'_{i-1}|$.

Therefore, for a Hermite parabola defined by $y_{i-1}$, $y_i$ and $y'_{i-1}$ to be monotone on the interval $(i-1, i)$, the derivative $y'_{i-1}$ must be of the same sign as the linear slope $\delta_{i-1}$ and smaller in magnitude than $2\delta_{i-1}$.

Convexity

For the function $f$ to be convex on $(i-1, i)$ (if parabola is convex on a segment, it's convex everywhere) the condition is $f''(x) > 0$. In the case of our parabola: $$f''(x) = 2a_2 > 0$$. $$\frac{y_{i} - y_{i-1} - y'_{i-1} h_{i-1}}{h_{i-1}^2} > 0$$ or $$y_{i} - y_{i-1} - y'_{i-1} h_{i-1} > 0$$ \begin{equation}\delta_{i-1} > y'_{i-1}\end{equation} This result is obvious $\delta_{i-1} = y'_{i-1}$ is the straight line going through the end-points of the segment. For the parabola to be concave, it must hold $f''(x) < 0$, i.e. \begin{equation}\delta_{i-1} < y'_{i-1}\end{equation} 

 

Cubic spline interpolation: Periodic splines

If the unknown function is periodic, then the periodicity may be used to specify the boundary conditions. Before (in the case of the natural splines) we had $n+1$ value of x and y and as a result we obtained values of the cubic polynomial at $n$ intervals.

Now our function is periodic with period $p$, thus $S(x) = S(x+p)$. In practice, we have an additional data point $x_{n+1}$, such that
$$x_{n+1} = x_0 + p,$$
and
$$y_{n+1} = y_0.$$

The total number of data points is $n+2$ and the total number of segments (and cubic splines) between them is $n+1$.

Cubic spline interpolation: Natural splines

Last year I published a post on the data interpolation and smoothing using the cubic splines using IDL. (Soon I'll do an update of that post as well.) Here I go a step back and show how the classical spline interpolation works. As before, the derivation follows the book of D.S.G.Pollock (highly recommended to any student interested in the analysis of time-series). I'll first go through some theory general for all the cubic splines, than through the derivation of the natural cubic splines and at the end I'll show some examples. In the next post I'll describe the solution of the clamped cubic splines and do the same examples using that variation of the method.

Cubic splines: Interpolation and smoothing

Cubic splines are widely used technique for piecewise smooth interpolation. Between any two points of a discretely given function, the function is approximated by a third order polynomial. However, splines can also be used for smoothing or fitting of measurements. The theory behind the splines is nicely explained by Pollock (1999). I implemented that algorithm in IDL procedure splinecoeff.pro.

To try it yourself, here is a small example. Download the text file with the values of the sunspot area data, daily_area.txt from solarscience.msfc.nasa.gov, and run the following:

; Load the data
n = FILE_LINES('daily_area.txt')-1
OPENR, 1, 'daily_area.txt'
a = ''
READF, 1, a
area = FLTARR(6, n)
READF, 1, area
CLOSE, 1

; Extract data between 1977-01-01 and 2011-12-31
initial = WHERE(area[0, *] EQ 1977 AND area[1, *] EQ  1 AND area[2, *] EQ  1)
final   = WHERE(area[0, *] EQ 2011 AND area[1, *] EQ 12 AND area[2, *] EQ 31)
time = JULDAY(area[1,*], area[2,*], area[0,*])
x = time[initial:final]
y = REFORM(area[3, initial:final])
nx = N_ELEMENTS(x)

; Compute coefficients of the cubic splines
coeff = SPLINECOEFF(x, y, lambda = 1.d3)

; Plot the original data and the spline function for given lambda
y1 = FLTARR(nx-1)
x1 = x[0:nx-2]
FOR i = 0, N_ELEMENTS(y)-2 DO y1[i] = coeff.d[i] + $
                                      coeff.c[i] * (x[i+1]-x[i]) + $
                                      coeff.b[i] * (x[i+1]-x[i])^2 + $
                                      coeff.a[i] * (x[i+1]-x[i])^3
PLOT, x, y, psym = 3
OPLOT, x1, y1

The only difference between the interpolation, smoothing and fitting of the measurement is in the free parameter $\lambda$. If $\lambda$ is very small, than the result is the cubic spline interpolation. If $\lambda$ is very large, the measurement is smoothed by cubic splines. In the extreme case, the resulting splines become a linear function. The following animation show how the result change with $\lambda$:



Note: See paper of Djurovic and Paquet (1999) for an example where the cubic splines are used for filtering and smoothing of data.