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.
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). The easiest interpolation method which preserves monotonicity is linear interpolation. However, it is not smooth. Polynomial interpolation is smooth but may introduce overshoots, violating monotonicity.
def plot_interp(interp_func, x, y):
grid_x = np.linspace(min(x), max(x), 100)
grid_y = interp_func(grid_x)
fig, ax = plt.subplots(figsize=(3, 3))
ax.plot(grid_x, grid_y)
ax.scatter(x, y, color='C1', zorder=2)
ax.set_xlabel(r"$x$")
ax.set_ylabel(r"$y$")import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
x = [-3, -2, -1, 0, 1, 2, 3]
y = [-2, -2, -2, 0, 2, 2, 2]
linear_interp = interp1d(x, y, kind="linear")
plot_interp(linear_interp, x, y)
cubic_interp = interp1d(x, y, kind="cubic")
plot_interp(cubic_interp, x, y)
One of the most effective interpolation methods that preserves monotonicity while maintaining smoothness is PCHIP. Some people also refer to it as shape-preserving cubic interpolation. The following figures illustrate PCHIP interpolation, where it also appliable to non-monotonic data points.
from scipy.interpolate import PchipInterpolator
x = [-3, -2, -1, 0, 1, 2, 3]
y = [-2, -2, -2, 0, 2, 2, 2]
pchip_interp = PchipInterpolator(x, y)
plot_interp(pchip_interp, x, y)
x = [-3, -2, -1, 0, 1, 2, 3]
y = [-2, -1, -1.5, 0, 1, 0.5, 2]
pchip_interp = PchipInterpolator(x, y)
plot_interp(pchip_interp, x, y)
How to ensure cubic interpolation preserves monotonicity? Let’s define some notations first. For each interval [x_i,x_{i+1}] (i=1,\dots,n-1), define the step size h_i and slope s_i as:
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.
The derivative d_i is computed using local information from three neighboring points \{(x_{i-1},y_{i-1}),(x_i,y_i),(x_{i+1},y_{i+1})\}. Replaceing with the slopes s_{i-1} and s_i, and step sizes h_{i-1} and h_i, the formula is (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}\left(1+\frac{h_i}{h_{i+1}+h_i}\right) \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 derivatives d_1 and d_n, they are computed separately (Moler 2004, chap. 3). For d_1:
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 \hat{d}_1=\hat{f}'(x_1). Additional rules are then added to preserve monotonicity. Similar procedures apply for d_n.
To ensure monotonicity, the derivatives d_i must satisfy certain conditions. 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 \left( \alpha_i,\beta_i\le3 \lor \phi(\alpha_i,\beta_i)\ge0 \right) \tag{11}
where:
\phi(\alpha,\beta)=\alpha-\frac{(2\alpha+\beta-3)^2}{3(\alpha+\beta-2)} \tag{12}
To verify Equation 5 and Equation 8 satisfy the monotonicity condition, we analyze them separately. 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.
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: (i) x^*<x_i; (ii) x^*>x_{i+1}; (iii) x_i\le x^*\le x_{i+1}.
Condition (i) 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 (ii) and (iii) 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 3 illustrates these conditions: Condition (i) — blue, Condition (ii) — green, Condition (iii) — red.
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=(3, 3))
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()
Finally, simplifying yields the final monotonicity condition (Lemma 1).
@misc{pchip,
author = {Xu, Xuefeng},
title = {Monotone {Piecewise} {Cubic} {Interpolation}},
date = {2025-03-13},
url = {https://xuefeng-xu.github.io/blog/pchip.html},
langid = {en}
}