- Generated by
1.13.2
|
fitpack
Modern Fortran library for curve and surface fitting with splines
|
This page describes the mathematical foundations for bivariate spline fitting in FITPACK: extending the one-dimensional B-spline framework to approximate surfaces \( z = s(x, y) \) defined over a rectangular domain.
A bivariate spline on the rectangular domain \( [a, b] \times [c, d] \) is defined as a tensor product of univariate B-splines:
\[ s(x, y) = \sum_{i=1}^{n_x - k_x - 1} \sum_{j=1}^{n_y - k_y - 1} c_{ij} \, N_{i, k_x}(x) \, M_{j, k_y}(y) \]
where:
The two knot vectors are independent: the number of knots and their placement can differ between the \( x \)- and \( y \)-directions. In FITPACK, both degrees default to \( k_x = k_y = 3 \) (bicubic).
Given \( m \) scattered observations \( (x_i, y_i, z_i) \) with weights \( w_i > 0 \), the goal is to find a tensor product spline \( s(x, y) \) that approximates the data while remaining smooth.
Each data point \( (x_i, y_i) \) lies in a unique cell of the knot grid. For a given point, let \( \mu \) and \( \nu \) index the non-zero basis functions in \( x \) and \( y \), respectively. The observation matrix \( \mathbf{A} \) has entries:
\[ A_{i, \ell} = N_{\mu, k_x}(x_i) \, M_{\nu, k_y}(y_i), \qquad \ell = (\mu - 1)(n_y - k_y - 1) + \nu \]
where the two-dimensional index pair \( (\mu, \nu) \) is mapped to a single column index \( \ell \). Each row of \( \mathbf{A} \) has at most \( (k_x + 1)(k_y + 1) \) non-zero entries, so the matrix is sparse.
The weighted least-squares problem is:
\[ \min_{\mathbf{c}} \sum_{i=1}^{m} w_i^2 \left( z_i - s(x_i, y_i) \right)^2 \]
which in matrix form becomes \( \min \| \mathbf{W}(\mathbf{z} - \mathbf{A} \mathbf{c}) \|_2^2 \) with \( \mathbf{W} = \mathrm{diag}(w_1, \ldots, w_m) \).
The smoothing spline minimizes a roughness functional subject to a data fidelity constraint:
\[ \min \iint \left[ \left( \frac{\partial^2 s}{\partial x^2} \right)^2 + 2 \left( \frac{\partial^2 s}{\partial x \, \partial y} \right)^2 + \left( \frac{\partial^2 s}{\partial y^2} \right)^2 \right] dx \, dy \quad \text{subject to} \quad \sum_{i=1}^{m} w_i^2 (z_i - s(x_i, y_i))^2 \leq S \]
where \( S \) is the smoothing factor. This is the bivariate analogue of the univariate smoothing spline problem.
When few data points fall within certain cells of the knot grid, the observation matrix \( \mathbf{A} \) can become rank deficient. This is more common in the bivariate case than in one dimension, because adding knots in both directions creates cells that may contain no data at all.
FITPACK handles rank deficiency using Householder QR factorization with column pivoting. When the numerical rank \( r \) of the system is less than the number of unknowns, the algorithm computes a minimum-norm solution among all least-squares solutions.
Starting from a minimal knot set, surfit iteratively refines the approximation:
The alternating strategy ensures balanced refinement across both coordinate directions.
When the data lie on a rectangular grid \( z_{ji} = z(y_j, x_i) \) with \( i = 1, \ldots, m_x \) and \( j = 1, \ldots, m_y \), the problem has a Kronecker product structure that can be exploited for dramatic computational savings.
On a grid, the observation matrix factors as a Kronecker product:
\[ \mathbf{A}_{\text{grid}} = \mathbf{B}_x \otimes \mathbf{B}_y \]
where \( \mathbf{B}_x \) is the \( m_x \times (n_x - k_x - 1) \) matrix of univariate B-spline values at the \( x \)-grid points, and \( \mathbf{B}_y \) is the \( m_y \times (n_y - k_y - 1) \) matrix at the \( y \)-grid points. The coefficient matrix conceptually satisfies:
\[ \mathbf{C} = \mathbf{B}_x^{-1} \, \mathbf{Z} \, (\mathbf{B}_y^T)^{-1} \]
where \( \mathbf{Z} \) is the \( m_x \times m_y \) matrix of observed values and \( \mathbf{C} \) holds the spline coefficients.
The factored form means the normal equations can be solved dimension by dimension rather than simultaneously:
| Approach | Complexity |
|---|---|
Scattered (surfit) | \( \mathcal{O}(m \cdot n_x \cdot n_y) \) |
Gridded (regrid) | \( \mathcal{O}(m_x \cdot m_y \cdot (k_x + k_y)) \) |
For a 100 x 100 grid with 20 knots in each direction, the gridded approach is roughly two orders of magnitude faster. Always prefer fitpack_grid_surface when data lie on a rectangular grid.
The smoothing framework (roughness minimization subject to a fidelity constraint) is identical to the scattered case; only the numerical solution method differs.
The Kronecker structure that makes regrid fast is not special to two dimensions. For data sampled on a \( d \)-dimensional rectangular grid \( z_{i_1 i_2 \cdots i_d} = z\big(x^{(1)}_{i_1}, x^{(2)}_{i_2}, \ldots, x^{(d)}_{i_d}\big) \), FITPACK fits a \( d \)-fold tensor-product B-spline
\[ s(x_1, \ldots, x_d) = \sum_{j_1} \sum_{j_2} \cdots \sum_{j_d} c_{j_1 j_2 \cdots j_d} \prod_{r=1}^{d} N^{(r)}_{j_r, k_r}(x_r) \]
where \( N^{(r)}_{j, k_r} \) is the \( j \)-th B-spline of degree \( k_r \) on the knot vector of axis \( r \). The routines that fit and evaluate the 2-D surface — regrid and ndspev — are written for runtime dimension dims, and the generic fitpack_gridded_spline type wraps them for any \( d \le \) MAX_IDIM.
The observation matrix is a \( d \)-fold Kronecker product
\[ \mathbf{A}_{\text{grid}} = \mathbf{B}_1 \otimes \mathbf{B}_2 \otimes \cdots \otimes \mathbf{B}_d \]
so the least-squares system separates: the coefficient tensor is obtained by solving banded systems one axis at a time, sweeping \( r = 1, \ldots, d \) and contracting each axis in turn. The cost is
\[ \mathcal{O}\!\left( \Big(\textstyle\prod_{r=1}^{d} m_r\Big) \cdot \sum_{r=1}^{d} k_r \right), \]
linear in the total number of grid points \( \prod_r m_r \) — the same per-point efficiency as the 2-D case, generalized.
In the smoothing mode, knots are added adaptively until the residual meets the target. With more than one axis, only one axis is refined per iteration: the arbiter new_knot_dimension_nd picks the axis with the largest residual that has not yet reached its knot ceiling, so refinement alternates between directions rather than exhausting one axis first. This is the \( d \)-way generalization of the two-way \( x \leftrightarrow y \) alternation in regrid.
The value grid is stored as a flat, row-major buffer with domain axis 1 the slowest-varying and axis \( d \) the fastest:
\[ \text{offset}(i_1, \ldots, i_d) = \Big( \cdots \big( (i_1 - 1)\, m_2 + (i_2 - 1) \big)\, m_3 + \cdots \Big)\, m_d + (i_d - 1). \]
The concrete 2-D fitpack_grid_surface stores \( z(i_y, i_x) \), whose column-major memory is exactly this row-major buffer with the axes reversed; fitpack_gridded_spline accepts either convention through an optional row_major flag (permuting only on the one-time ingest copy).
The smoothing factor scales with the grid size: a natural starting choice is \( S \approx \prod_{r=1}^{d} m_r \) (roughly unit expected squared residual per point), decreased to follow more detail. Two distinct "dimension" counts appear and must not be confused:
| Symbol | Meaning |
|---|---|
\( d \) (dims, ndim) | number of domain axes \( x_1, \ldots, x_d \) |
idim | number of codomain components (for vector-valued data) |
The gridded routines fit a scalar field ( \( \texttt{idim} = 1 \)); dims is what the generic type varies at runtime.
Memory, not arithmetic, sets the limit. The coefficient tensor has \( \prod_r (n_r - k_r - 1) \) entries and the value grid \( \prod_r m_r \), both growing exponentially in \( d \). In practice \( d \approx 3\text{–}5 \) is the useful range for dense grids; beyond that the full tensor becomes memory-bound, and sparse or reduced-rank representations are preferable.
regrid, ndspev, fitpack_gridded_spline; Dierckx, Ch. 5, §5.4For bivariate splines, the roughness measure used in the smoothing formulation is the thin-plate energy functional:
\[ R[s] = \iint \left[ \left( \frac{\partial^2 s}{\partial x^2} \right)^2 + 2 \left( \frac{\partial^2 s}{\partial x \, \partial y} \right)^2 + \left( \frac{\partial^2 s}{\partial y^2} \right)^2 \right] dx \, dy \]
This is the natural extension of \( \int (s'')^2 \, dx \) to two dimensions. It penalizes curvature equally in all directions and is invariant under rotation of the coordinate axes. The cross-derivative term \( (\partial^2 s / \partial x \, \partial y)^2 \) ensures that diagonal oscillations are penalized as well.
For a tensor product spline, this integral can be expressed as a quadratic form in the coefficients \( c_{ij} \), involving integrals of products of B-spline second derivatives. FITPACK computes these integrals analytically using the knot vectors.
The smoothing factor \( S \) controls the balance between fidelity and smoothness:
For scattered data with unit weights ( \( w_i = 1 \)), a reasonable starting point is:
\[ S \approx m \]
where \( m \) is the number of data points. This follows from the expectation that, for a good fit with normally distributed residuals, the sum of squared residuals should be \( \mathcal{O}(m) \).
For gridded data with unit weights, the analogous rule is:
\[ S \approx m_x \cdot m_y \]
After fitting, inspect the mean squared error via the mse() method to assess fit quality. If the residual is too large, decrease \( S \); if the surface oscillates, increase \( S \).
When data have non-unit weights, interpret \( S \) relative to the weighted residual. Setting \( w_i = 1 / \sigma_i \) (where \( \sigma_i \) is the measurement standard deviation) makes the residual sum an approximate \( \chi^2 \) statistic, and \( S \approx m \) corresponds to a reduced \( \chi^2 \) of unity.
For surfaces embedded in \( \mathbb{R}^d \) and parameterized by \( (u, v) \), each coordinate component is represented as a separate bivariate spline sharing common knot vectors:
\[ s_l(u, v) = \sum_{i=1}^{n_u - k_u - 1} \sum_{j=1}^{n_v - k_v - 1} c_{ij}^{(l)} \, N_{i, k_u}(u) \, M_{j, k_v}(v), \qquad l = 1, \ldots, d \]
All \( d \) components share the same knot vectors \( \mathbf{t}_u \) and \( \mathbf{t}_v \), and the same spline degrees \( k_u \) and \( k_v \). Only the coefficient arrays \( c_{ij}^{(l)} \) differ between components.
Either parameter direction (or both) may be periodic:
Periodicity is useful for closed surfaces such as cylinders and tori. The FITPACK routine parsur handles the general case.
| Problem | FITPACK Type | Core Routine | Key Advantage |
|---|---|---|---|
| Scattered \( (x_i, y_i, z_i) \) | fitpack_surface | surfit | Handles arbitrary point distributions |
| Gridded \( z(y_j, x_i) \) | fitpack_grid_surface | regrid | Kronecker structure, much faster |
| Parametric \( \mathbf{r}(u, v) \) | fitpack_parametric_surface | parsur | Surfaces in \( \mathbb{R}^d \) |
1.13.2