Monotone Piecewise Cubic Interpolation

Statistics
Published

March 13, 2025

Modified

March 30, 2025

Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) is a cubic spline-based interpolation method designed to preserve monotonicity. See MATLAB or SciPy for the implementation details.

1 Interpolation Function

Given n data points (x_1,y_1),\dots,(x_n,y_n) with x_1<\cdots<x_n, where y is monotonic (either y_i\le y_{i+1} or y_{i+1}\ge y_i), define:

h_i=x_{i+1}-x_i,\ s_i=\frac{y_{i+1}-y_i}{h_i} \tag{1}

For x_i<x<x_{i+1}, the interpolation function f(x) is a cubic polynomial:

f(x)=c_0+c_1(x-x_i)+c_2(x-x_i)^2+c_3(x-x_i)^3 \tag{2}

satisfying:

\begin{align*} f(x_i)&=y_i, &f(x_{i+1})&=y_{i+1}\\ f'(x_i)&=d_i, &f'(x_{i+1})&=d_{i+1} \end{align*} \tag{3}

Solving for the coefficients:

\begin{align*} c_0&=y_i, &c_2&=\frac{3s_i-2d_i-d_{i+1}}{h_i}\\ c_1&=d_i, &c_3&=\frac{-2s_i+d_i+d_{i+1}}{h_i^2}\\ \end{align*} \tag{4}

Thus, computing derivatives d_1,\dots,d_n determines c_0,c_1,c_2,c_3 for each interval.

2 Derivative Computation

The derivative at x_i is computed using local information from three neighboring points (Fritsch and Butland 1984):

d_i=G(s_{i-1},s_i,h_{i-1},h_i)= \begin{cases} \frac{s_{i-1}s_i}{rs_i+(1-r)s_{i-1}} & \mathrm{if~}s_{i-1}s_i>0, \\ 0 & \text{otherwise} & \end{cases} \tag{5}

where the ratio r (1/3<r<2/3) is given by:

r=\frac{h_{i-1}+2h_i}{3h_{i-1}+3h_i}=\frac{1}{3}\Big(1+\frac{h_i}{h_{i+1}+h_i}\Big) \tag{6}

If s_{i-1} and s_i have opposite signs (indicating non-monotonicity) or one is zero, then d_i=0. Otherwise, d_i is a weighted harmonic mean:

\frac{1}{d_i}=\frac{r}{s_{i-1}}+\frac{1-r}{s_i} \tag{7}

For endpoints x_1 and x_n, derivatives are computed separately (Moler 2004, chap. 3):

d_1= \begin{cases} 0 & \mathrm{if~}\text{sgn}(\hat{d}_1)\neq \text{sgn}(s_1), \\ 3s_1 & \mathrm{if~}\text{sgn}(s_1)\neq \text{sgn}(s_2) \land|\hat{d}_1|>3|s_1|,\\ \hat{d}_1 & \text{otherwise} & \end{cases} \tag{8}

where:

\hat{d}_1=\frac{(2h_1+h_2)s_1 - h_1s_2}{{h_1+h_2}} \tag{9}

A quadratic polynomial \hat{f}(x)=\hat{c}_0+\hat{c}_1x+\hat{c}_2x^2 is fit through the first three points, and its derivative at x_1 is computed to obtain \hat{d}_1. Additional rules are then applied to preserve monotonicity. Similar rules apply for d_n.

3 Monotonicity Conditions

To ensure monotonicity, define:

\alpha_i=\frac{d_i}{s_i},\ \beta_i=\frac{d_{i+1}}{s_i} \tag{10}

Lemma 1 A sufficient condition for monotonicity is (Fritsch and Carlson 1980):

\alpha_i,\beta_i\ge0 \land \Big( \alpha_i,\beta_i\le3 \lor \phi(\alpha_i,\beta_i)\ge0 \Big) \tag{11}

where:

\phi(\alpha,\beta)=\alpha-\frac{(2\alpha+\beta-3)^2}{3(\alpha+\beta-2)} \tag{12}

If s_{i-1}s_i>0, then \alpha_i>0; otherwise, \alpha_i=0. Since the ratio r satisfies 1/3<r<2/3, \alpha_i is upper-bounded by 3:

\alpha_i =\frac{1}{rs_i/s_{i-1}+(1-r)} <\frac{1}{1-r}<3 \tag{13}

For endpoint \alpha_1, we only need to show the condition of \text{sgn}(\hat{d}_1)=\text{sgn}(s_1)=\text{sgn}(s_2), since other conditions already lie within the region [0,3].

\alpha_1 =1+\frac{1-s_2/s_1}{1+h_2/h_1}<2 \tag{14}

Similarly, \beta_i and endpoint \beta_{n-1} all satisfy the monotonicity condition.

4 Proof of Monotonicity Conditions

Proof. To preserve monotonicity, the derivatives d_i and d_{i+1} must align with the direction of the slope of the interval s_i. This is a necessary condition (Fritsch and Carlson 1980):

\text{sgn}(d_i)=\text{sgn}(d_{i+1})=\text{sgn}(s_i) \Leftrightarrow \alpha_i,\beta_i\ge0 \tag{15}

The derivative of f(x) is a quadratic polynomial:

f'(x)=c_1+ 2c_2(x-x_i)+ 3c_3(x-x_i)^2 \tag{16}

It has a unique extremum at:

x^*=x_i+\frac{h_i}{3}\cdot\frac{2\alpha_i+\beta_i-3}{\alpha_i+\beta_i-2} \tag{17}

and

f'(x^*)=\phi(\alpha_i,\beta_i)s_i \tag{18}

There are three conditions to check: 1) x^*<x_i; 2) x^*>x_{i+1}; 3) x_i\le x^*\le x_{i+1}.

Condition 1) is equivalent to \alpha_i,\beta_i\ge0 \land \frac{2\alpha_i+\beta_i-3}{\alpha_i+\beta_i-2}<0 \tag{19}

The analysis of Condition 2) and 3) are similar, leading to:

\alpha_i,\beta_i\ge0 \land \frac{\alpha_i+2\beta_i-3}{\alpha_i+\beta_i-2}<0 \tag{20}

and

\alpha_i,\beta_i\ge0 \land \frac{2\alpha_i+\beta_i-3}{\alpha_i+\beta_i-2}\ge0 \land \frac{\alpha_i+2\beta_i-3}{\alpha_i+\beta_i-2}\ge0 \land \phi(\alpha_i,\beta_i)\ge0 \tag{21}

Figure 1 illustrates these conditions separately: Condition 1) — blue, Condition 2) — green, Condition 3) — red.

Code 
import numpy as np
import matplotlib.pyplot as plt

x_vals = np.linspace(0, 4, 1000)
y_vals = np.linspace(0, 4, 1000)
x, y = np.meshgrid(x_vals, y_vals)

cond1 = (x > 0) & (y > 0) & ((2 * x + y - 3) / (x + y - 2) < 0)
cond2 = (x > 0) & (y > 0) & ((x + 2 * y - 3) / (x + y - 2) < 0)
cond3 = (
    (x > 0)
    & (y > 0)
    & ((2 * x + y - 3) / (x + y - 2) > 0)
    & ((x + 2 * y - 3) / (x + y - 2) > 0)
    & (x - ((2 * x + y - 3) ** 2) / (3 * (x + y - 2)) > 0)
)

fig, ax = plt.subplots(figsize=(4, 4))

ax.contourf(x, y, cond1, levels=1, colors=["white", "blue"], alpha=1)
ax.contourf(x, y, cond2, levels=1, colors=["white", "green"], alpha=0.5)
ax.contourf(x, y, cond3, levels=1, colors=["white", "red"], alpha=0.25)

ax.set_xlabel(r"$\alpha$")
ax.set_ylabel(r"$\beta$")
ax.set_xlim([0, 4])
ax.set_ylim([0, 4])
ax.xaxis.set_ticks(list(np.arange(0, 4.5, 0.5)))
ax.yaxis.set_ticks(list(np.arange(0, 4.5, 0.5)))

ax.grid(True)
ax.tick_params(
    bottom=True,
    top=True,
    left=True,
    right=True,
    labelbottom=True,
    labeltop=True,
    labelleft=True,
    labelright=True,
)
plt.show()
Figure 1: The monotonicity region

Finally, simplifying yields the final monotonicity condition (Lemma 1).

5 Cubic Hermite Interpolation

Cubic Hermite Interpolation constructs f(x) using both function values and derivatives:

f(x)=y_iH_1(x)+y_{i+1}H_2(x)+d_iH_3(x)+d_{i+1}H_4(x) \tag{22}

where basis functions are:

\begin{align*} H_{1}(x)&=\phi\Big(\frac{x_{i+1}-x}{h_{i}}\Big), &H_{2}(x)&=\phi\Big(\frac{x-x_{i}}{h_{i}}\Big)\\ H_{3}(x)&=-h_{i}\psi\Big(\frac{x_{i+1}-x}{h_{i}}\Big), &H_{4}(x)&=h_{i}\psi\Big(\frac{x-x_{i}}{h_{i}}\Big) \end{align*} \tag{23}

with:

\phi(t)=3t^{2}-2t^{3},\ \psi(t)=t^{3}-t^{2} \tag{24}

This formulation aligns with the cubic polynomial definition above.

References

Fritsch, F. N., and J. Butland. 1984. “A Method for Constructing Local Monotone Piecewise Cubic Interpolants.” SIAM Journal on Scientific and Statistical Computing 5 (2): 300–304. https://doi.org/10.1137/0905021.
Fritsch, F. N., and R. E. Carlson. 1980. “Monotone Piecewise Cubic Interpolation.” SIAM Journal on Numerical Analysis 17 (2): 238–46. https://doi.org/10.1137/0717021.
Moler, Cleve B. 2004. Numerical Computing with Matlab. Society for Industrial; Applied Mathematics. https://doi.org/10.1137/1.9780898717952.ch3.