Tchebycheff splines

[Matt Piercy]

Has anyone tried to program these in mata? I would be interested in the code if they are willing to share it.

[Choonjoo Lee]

I have a Stata/Mata implementation that I would be happy to share.

The current ado and help files are still research/development code. Although they have passed a set of internal numerical and reproducibility tests, they have not yet undergone extensive independent validation. Feedback, test cases, and comparisons with other implementations would therefore be very welcome.

The appropriate routines depend on the problem you wish to solve. If you could briefly describe your application—including the spline space, knot structure, continuity conditions, and any shape restrictions—I will try to provide the most relevant code.

The current implementation supports:

Piecewise-linear splines, including monotonicity, convexity, and concavity constraints, as well as optional free-knot estimation.
Polynomial multi-degree B-splines (MDB-splines), with interval-specific polynomial degrees and knot-specific continuity orders.
Multi-degree Tchebycheffian B-splines (MDTB-splines), with local spaces specified through the characteristic roots of constant-coefficient linear differential operators. This provides a common framework for polynomial, exponential, and trigonometric spline spaces.

Conditional on a fixed basis and fixed knots, the coefficient problem is formulated and solved as a finite linear program. Free-knot estimation is currently available for the piecewise-linear basis and uses deterministic profile searches with multiple starting values. The resulting knot locations should therefore be interpreted as local solutions rather than guaranteed global optima.

Shape constraints are imposed through linear constraints and adaptive numerical separation, depending on the basis. In particular, shape verification for the MDTB basis is numerical rather than based on a formal interval-arithmetic certificate.

If you let me know what kind of spline problem you are considering, I will be happy to share the relevant ado and help files and assist as much as I can.

Choonjoo.

[Matt Piercy]

Hi Choonjoo,
I tried the splines out on a piecewise linear spline for a dataset to visualise c-section times. These worked OK. Then I read about Tchbycheff splines and got them to work in python but the result was completely different. I would appreciate if you could take a look. Thanks, Matt

[Matt Piercy]

sorry I forgot to put in the graph code.

[Matt Piercy]

This is a do file of the spline and a Tchebycheff file in python for same data.

[Choonjoo Lee]

Hi Matt,

I traced the discrepancy to the fact that the Python code and Mata's
spline3() solve different problems.

Your Python code creates three independent cubic Chebyshev polynomial blocks.
There are 12 coefficients but only 4 observations, and no continuity
constraints are imposed at the internal knots. The shared endpoints are also
included in both adjacent blocks. np.linalg.lstsq() therefore selects one
solution to an underdetermined problem.

Mata spline3(), in contrast, computes a natural cubic interpolating spline:
it passes through all observations, is C2 at the interior knots, and has zero
second derivative at the two boundary knots. The curves therefore need not
agree.

I have attached a standalone Python script that reproduces the original
Chebyshev-block calculation, and a second small script that reproduces the
natural cubic spline used by spline3(). Both use the same Emma data. The PNG
figures show the resulting curves.

The original Python code can run outside Stata as a .py file. If it is run
inside Stata 16, the Python executable must be configured first, and the
Python block must be closed with end. The attached standalone scripts avoid
that Stata/Python setup issue.

Best,
Choonjoo

1.T-spline emma's data.py

import numpy as np
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt


def cheb_poly(n, x):
"""Compute Chebyshev polynomial T_n(x)."""
if n == 0:
return np.ones_like(x)
if n == 1:
return x
t0, t1 = np.ones_like(x), x
for _ in range(2, n + 1):
t0, t1 = t1, 2 * x * t1 - t0
return t1


def cheb_basis(m, x):
"""Return [T0(x), T1(x), ..., T_m(x)]."""
return np.column_stack([cheb_poly(k, x) for k in range(m + 1)])


def cheb_spline_basis(x, knots, m):
"""Build the independent piecewise Chebyshev blocks in the original code."""
blocks = []
for a, b in zip(knots[:-1], knots[1:]):
mask = (x >= a) & (x <= b)
block = np.zeros((x.size, m + 1), dtype=float)
if np.any(mask):
xi = 2 * (x[mask] - a) / (b - a) - 1
block[mask, :] = cheb_basis(m, xi)
blocks.append(block)
return np.hstack(blocks)


def cheb_spline_interpolate(x_data, y_data, m=3):
knots = np.linspace(min(x_data), max(x_data), len(x_data))
design = cheb_spline_basis(x_data, knots, m)
coeffs, *_ = np.linalg.lstsq(design, y_data, rcond=None)
return coeffs, knots


def cheb_spline_eval(x, coeffs, knots, m=3):
with np.errstate(divide="ignore", over="ignore", invalid="ignore"):
return cheb_spline_basis(x, knots, m) @ coeffs


x_data = np.array([1, 2, 3, 4], dtype=float)
y_data = np.array([0.392, 0.323, 0.201, 0.083], dtype=float)

coeffs, knots = cheb_spline_interpolate(x_data, y_data, m=3)
x_grid = np.linspace(x_data.min(), x_data.max(), 400)
y_grid = cheb_spline_eval(x_grid, coeffs, knots, m=3)

plt.figure(figsize=(8, 5))
plt.plot(x_grid, y_grid, label="Original Python Chebyshev blocks", color="blue")
plt.scatter(x_data, y_data, color="red", s=60, label="Data points")
plt.title("Original Python code: independent Chebyshev blocks")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig("emma_python_original.png", dpi=160)
print("wrote emma_python_original.png")
print("knots:", knots)
print("number of coefficients:", len(coeffs))
print("interpolation residuals:", y_data - cheb_spline_eval(x_data, coeffs, knots, m=3))

2. reproduce_stata_spline3.py

#!/usr/bin/env python3
"""Reproduce Mata spline3() natural cubic interpolation in Python.

The local representation matches Mata's spline3() output:

s_i(x) = y_i + b_i (x-x_i) + c_i (x-x_i)^2 + d_i (x-x_i)^3

for x_i <= x <= x_{i+1}. The spline interpolates every input point, is C2
at interior knots, and has zero second derivative at both endpoints.

Input CSV format: one x column and one or more response columns, e.g.

x,y,z
1,.392,.512
2,.323,.269
3,.201,.138
4,.083,.081

Example:
python3 reproduce_stata_spline3.py --input emma_data.csv \
--target y z --output emma_spline3.png
"""

from __future__ import annotations

import argparse
import csv
from pathlib import Path
from typing import Iterable

def _solve_tridiagonal(lower: list[float], diagonal: list[float],
upper: list[float], rhs: list[float]) -> list[float]:
"""Solve a tridiagonal system using the Thomas algorithm."""
n = len(diagonal)
if n == 0:
return []
if len(lower) != n - 1 or len(upper) != n - 1 or len(rhs) != n:
raise ValueError("invalid tridiagonal system")
c = upper[:]
d = rhs[:]
b = diagonal[:]
for i in range(1, n):
if b[i - 1] == 0:
raise ValueError("singular tridiagonal system")
factor = lower[i - 1] / b[i - 1]
b[i] -= factor * c[i - 1]
d[i] -= factor * d[i - 1]
if b[-1] == 0:
raise ValueError("singular tridiagonal system")
solution = [0.0] * n
solution[-1] = d[-1] / b[-1]
for i in range(n - 2, -1, -1):
solution[i] = (d[i] - c[i] * solution[i + 1]) / b[i]
return solution


def natural_cubic_coefficients(x: Iterable[float], y: Iterable[float]) -> list[list[float]]:
"""Return [b, c, d] rows in the same local form as spline3()."""
x = [float(value) for value in x]
y = [float(value) for value in y]
if len(x) != len(y):
raise ValueError("x and y must be one-dimensional arrays of equal length")
if len(x) < 2 or any(right <= left for left, right in zip(x, x[1:])):
raise ValueError("x must contain at least two strictly increasing values")

h = [right - left for left, right in zip(x, x[1:])]
n = len(x)

# M_i is the second derivative at knot i. Natural boundaries: M_0=M_n-1=0.
M = [0.0] * n
if n > 2:
lower = h[1:-1]
diagonal = [2.0 * (h[i - 1] + h[i]) for i in range(1, n - 1)]
upper = h[1:-1]
rhs = [6.0 * ((y[i + 1] - y[i]) / h[i] -
(y[i] - y[i - 1]) / h[i - 1])
for i in range(1, n - 1)]
interior = _solve_tridiagonal(lower, diagonal, upper, rhs)
M[1:-1] = interior

coefficients = []
for i in range(n - 1):
b = (y[i + 1] - y[i]) / h[i] - h[i] * (2.0 * M[i] + M[i + 1]) / 6.0
c = M[i] / 2.0
d = (M[i + 1] - M[i]) / (6.0 * h[i])
coefficients.append([b, c, d])
return coefficients


def spline3_eval(x: Iterable[float], y: Iterable[float],
coefficients: list[list[float]], query: Iterable[float]) -> list[float]:
"""Evaluate a spline represented by natural_cubic_coefficients()."""
x = [float(value) for value in x]
y = [float(value) for value in y]
q = [float(value) for value in query]
if any(value < x[0] or value > x[-1] for value in q):
raise ValueError("query points must lie within the input x range")
values = []
for value in q:
interval = 0
while interval < len(x) - 2 and value >= x[interval + 1]:
interval += 1
dx = value - x[interval]
b, c, d = coefficients[interval]
values.append(y[interval] + b * dx + c * dx**2 + d * dx**3)
return values


def read_csv(path: Path) -> dict[str, list[float]]:
with path.open(newline="") as stream:
rows = list(csv.DictReader(stream))
if not rows or "x" not in rows[0]:
raise ValueError("CSV must contain an x column and at least one response column")
columns = {name: [float(row[name]) for row in rows]
for name in rows[0]}
return columns


def main() -> None:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--input", required=True, type=Path)
parser.add_argument("--target", nargs="+", required=True)
parser.add_argument("--grid", type=int, default=400)
parser.add_argument("--output", type=Path, default=None,
help="optional PNG plot; requires matplotlib")
args = parser.parse_args()

data = read_csv(args.input)
x = data["x"]
if args.grid < 2:
raise ValueError("--grid must be at least 2")
query = [x[0] + (x[-1] - x[0]) * i / (args.grid - 1)
for i in range(args.grid)]
fitted = {}

for target in args.target:
if target not in data:
raise ValueError(f"target column not found: {target}")
y = data[target]
coefficients = natural_cubic_coefficients(x, y)
fitted[target] = spline3_eval(x, y, coefficients, query)

print(f"[{target}] local coefficients [b, c, d]")
for i, row in enumerate(coefficients, start=1):
print(f" interval {i}: {row[0]: .12g} {row[1]: .12g} {row[2]: .12g}")
at_data = spline3_eval(x, y, coefficients, x)
max_error = max(abs(fitted_value - observed_value)
for fitted_value, observed_value in zip(at_data, y))
print(f" max interpolation error: {max_error:.3e}")
print()

if args.output is not None:
try:
import matplotlib.pyplot as plt
except ImportError as exc:
raise SystemExit("--output requires matplotlib") from exc
for target in args.target:
plt.plot(query, fitted[target], label=f"spline3: {target}")
plt.scatter(x, data[target], s=35)
plt.xlabel("x")
plt.ylabel("value")
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig(args.output, dpi=160)
print(f"wrote {args.output}")


if __name__ == "__main__":
main()


//. end.

[Choonjoo Lee]

Hi Matt,

To clarify, Mata's spline3() is a natural cubic interpolating spline, not a
general Tchebycheffian spline implementation. It is useful as a
cubic-polynomial baseline, but it does not implement the exponential or
trigonometric local spaces that can be used in a Tchebycheffian spline.

I also have a separate Stata/Mata implementation of multi-degree
Tchebycheffian B-splines (MDTB-splines). The local spaces are specified
through characteristic roots, with fixed knots, continuity conditions, and
minimax fitting. The implementation can represent polynomial, exponential,
trigonometric, and mixed local spaces.

If your goal is to fit a genuinely Tchebycheffian spline, I can provide the
relevant user-written commands and examples(let me know your email address). It would be helpful to know:

1. the intended local spaces or characteristic roots;
2. the knot locations;
3. the desired continuity conditions; and
4. whether you want interpolation, least-squares fitting, or minimax fitting.

Please note that the MDTB implementation is still research/development code.
It has passed internal numerical and reproducibility tests, but it has not yet
undergone extensive independent validation. The shape-restriction checks for
the MDTB basis are numerical rather than formal interval-arithmetic
certificates. Therefore, any results should be treated as experimental and
independently checked for the specific application.

Best,
Choonjoo

[Matt Piercy]

Hi Choonjoo,
Thanks. I gather your python code is reproducing in python a cubic spline similar to what I have in stata. I wanted to see what the TB splines looked like. I couldn't understand why they looked so different from the cubic splines. The explanation copilot gave me was that they are better for solving multidimensional problems. I would be grateful if you could send me you stata/mata code or do file for the T-spline to see whether they looked different from the python ones. I would be interested in interpolation spline with the `emma's data' data point, so the knots should be the data points and the ends. I understand it is not for publication. I am trying to extend my understanding to see if this is something useful. I did fit an L-spline to the same data in python, and it looked pretty similar to the cubic spline. I have attached the do file for this.
Thanks,
Matt
E. [email protected]

[Choonjoo Lee]

Hi Matt,
I’ve sent the details and files by email.
Cheers!