Automatic computation of hierarchical biquadratic smoothing splines with minimum GCV

A computationally efficient numerical strategy for fitting approximate minimum GCV bivariate thin plate smoothing splines to large noisy data sets was developed. The procedure discretises the bivariate thin plate smoothing spline equations using biquadratic B-splines and uses a nested grid SOR iterative strategy to solve the discretised system. For efficient optimisation, the process incorporates a double iteration that simultaneously updates both the discretised solution and the estimate of the minimum GCV smoothing parameter. The GCV was estimated using a minimum variance stochastic estimator of the trace of the influence matrix associated with the fitted spline surface. A Taylor series expansion was used to estimate the smoothing parameter that minimises the GCV estimate. The computational cost of the procedure is optimal in the sense that it is proportional to the number of grid points supporting the fitted biquadratic spline. Convergence was improved by adding a first order correction to the solution estimate after each smoothing parameter update. The algorithm was tested on several simulated data sets with varying spatial complexity and noise level. An accurate approximation to the analytic minimum GCV thin plate smoothing spline was obtained in all cases.