- 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.
For 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