Curve and surface fitting

The ALGLIB numerical library offers a comprehensive suite of curve and surface fitting algorithms, capable of efficiently handling various small- and large-scale fitting problems.

The library supports 1D curve fitting algorithms, such as polynomial, rational, penalized spline, and 4PL/5PL fitting. It also features N-dimensional fitting methods, such as penalized large-scale 2D splines, large-scale smoothing RBF splines, large-scale TPS and IDW, as well as general N-dimensional linear and nonlinear least squares solvers.

ALGLIB is available in C++, C#, Java and other languages under free and commercial licenses.

ALGLIB curve/surface fitting functionality

Introduction

Curve fitting is intricately related to interpolation and least squares problems, yet it exhibits distinct differences. Interpolation involves functions that precisely fit data points, representing a very specific subset of curve fitting. ALGLIB offers numerous interpolation methods, many of which partially overlap with its curve fitting functionality.

Linear or nonlinear least squares can be viewed as a generalization of curve fitting. The difference is that curve fitting usually fits the same curve (parametric function) computed at various points. In contrast, least squares problems can involve a diverse set of distinct squared functions to be minimized. ALGLIB has a comprehensive suite of linear and nonlinear least squares solvers. Often, the same algorithms are applicable both in curve fitting and generic least squares, depending on the context.

Programming languages supported

ALGLIB supports many programming languages, including C++, C#, Java, Python, and others:

• C++ version. Highly optimized, portable, self-contained library (x86/x64/ARM; works on Windows, Linux, and POSIX systems)
• C# version. Unified C# API for two backends: a fully managed C# implementation and a high-performance native computational core written in C
• Java, Python, and other languages. A wrapper for a high-performance native computational core written in C. Seamless integration with the target programming language

A distinctive feature of ALGLIB is that it provides the same API in all programming languages. This is achieved through our exclusive automatic code translation and wrapper generation technology.

1D curve fitting

Penalized smoothing cubic splines

Penalized (smoothing) cubic splines are a widely used method of 1D curve fitting. In this approach, a cubic spline with M equidistant nodes is fitted to N scattered points. The fitting process involves minimizing a sum of squared deviations from target values, with an additional nonlinearity penalty applied. This functionality is provided by the spline1dfit function from the spline1d subpackage.

Straightforward approach to the problem involves solution of a big dense system, with O(M²) memory and O((N+M)·M²) running time requirements. ALGLIB includes sophisticated implementation of the smoothing spline solver which has O(max(M,N)) memory and time requirements, making it scalable to big datasets with millions of points. The function can easily handle both M≪N and M≫N cases, e.g. it is possible to process 100 points with a 1000-node spline or vice versa.

Polynomial fitting

Fitting by polynomials is another popular method of 1D curve fitting. ALGLIB incorporates an advanced solver for constrained polynomial fitting, which employs numerous mathematical improvements to do robust fits even with high-degree polynomials. Internally, this solver utilizes a stable Chebyshev basis for performing the fit, before converting the polynomial into a stable and efficient barycentric form.

The unconstrained and constrained polynomial fitting functionality is provided by the polynomialfit and polynomialfitwc functions of the lsfit subpackage. Additionally, the resulting polynomial can be converted into a power basis with the polynomialbar2pow function.

Rational fitting

Rational functions are often used for high-accuracy fits of noise-free data. Among these, the Floater-Hormann class of rational functions is particularly noted for its stability characteristics. This model ensures the absence of poles (infinities) within the interval of fitting.

The barycentricfitfloaterhormann and barycentricfitfloaterhormannwc functions perform a Floater-Hormann least squares fit, producing a rational function in the barycentric form.

4-parameter and 5-parameter logistic (4PL and 5PL)

The Four Parameter Logistic (4PL) and Five Parameter Logistic (5PL) curves are widely used in biochemistry. ALGLIB provides a comprehensive set of functions for both unconstrained and constrained 4PL/5PL fits:

• logisticfit4 and logisticfit5 perform an unconstrained 4PL/5PL fit
• logisticfit4ec and logisticfit5ec perform an equality constrained 4PL/5PL fit
• logisticcalc4 and logisticcalc5 convenience functions can be used to evaluate the curve given its coefficients

Despite its apparent simplicity, the 4PL/5PL fitting problem is challenging due to issues like ill conditioning and multiple bad local extrema. ALGLIB successfully addresses these issues by performing several rounds of regularized optimization with gradually relaxed constraints, and by making several restarts from random initial solutions.

N-dimensional curve fitting

Penalized bicubic splines

Smoothing bicubic splines are widely used for 2D surface fitting. Similarly to the 1D case, an equidistant bicubic spline with P·Q nodes is fitted to N scattered points, subject to an additional nonlinearity penalty. This fitting process is notoriously challenging. When solved with a dense solver, the computational demands become prohibitive for grids larger than 50x50 nodes, due to the O((PQ)³) running time.

An unique feature of ALGLIB are two high-performance bicubic fitting solvers capable of handling large-scale problems:

• The BlockLLS solver, which leverages the sparsity of the problem. It generates C2-continuous splines. It is the most robust and controllable solver, capable of handling problems with up to 512x512 nodes and up to a 1.000.000 points.
• The FastDDM solver, which divides the fitting problem into smaller segments and then combines the results. This method, known as domain decomposition, results in fast parallelized computation. However, it produces merely C1 continuous models, as the decomposition process disrupts C2 continuity. This solver is optimized for large-scale problems, starting from 256x256 grids, and up to 10000x10000 grids.

The bicubic spline fitting functionality is provided by the spline2d subpackage.

RBF and thin plate spline fitting

Radial Basis Functions (RBFs), including Gaussian RBFs, biharmonic splines, thin plate splines and multiquadrics, are among the most powerful methods for N-dimensional surface fitting. They support both interpolation and smoothing, can easily generalize to 3 or more dimensions, and naturally support scattered interpolation nodes.

ALGLIB includes two state-of-the-art RBF construction algorithms, both of them suited for smoothing and interpolation of large-scale datasets:

• The HRBF (Hierarchical RBF) algorithm fits models with a Gaussian basis function, leveraging the sparsity of the problem. When used with small radii (up to ten times the mean interpoint distance), it scales well to hundreds of thousands of points.
• The DDM-RBF algorithm fits biharmonic splines, thin plate splines and multiquadrics. It employs a domain decomposition method to split a problem into smaller subproblems. The algorithm naturally supports parallelism and scales well to hundreds of thousands of points.

Both algorithms are provided by the rbf subpackage of the Interpolation package.

IDW models

Inverse Distance Weighting (IDW) models, including the modified Shepard's method and the modified-stabilized algorithm (MSTAB) are frequently used for large-scale surface fitting. This functionality is provided by the idw subpackage of the Interpolation package.

Because these models were primarily designed for interpolation, they do not support controlled smoothing. Additionally, they lack the extrapolation properties of thin plate splines, which makes them less suitable for tasks such as shape or volume reconstruction. On the other hand, they excel on large-scale uniform datasets due to fast model construction and evaluation.

Fitting by general curves/surfaces

General models can be fit to known data points using either a curve/surface fitting paradigm, or a nonlinear least squares paradigm.

General curve/surface fitting involves a single function, f(X|C), with two separate sets of arguments. Here X represents the coordinates of dataset points, while C denotes the parameters being fitted. Typically, X and C differ in dimensionality: X's dimensionality matches the space we're working in (usually 1, 2, or 3 dimensions), and C's dimensionality equals the number of parameters in the model we're fitting.

The lsfit subpackage, part of the Interpolation package, provides this functionality. It can fit differentiable models using the Levenberg-Marquardt algorithm, which employs either analytic derivatives or numerical differentiation. The solver, in its commercial version, is capable of performing parallel model evaluations to speed up computations. The same subpackage also provides linear least squares solvers.

Curve/surface fitting can be also seen as a nonlinear least squares problem. The former focuses on minimizing the sum of squared terms which are essentially the same function computed at different points. The latter minimizes the sum of squared functions, which may be entirely unrelated. Obviously, this includes curve/surface fitting as a special case.

Nonlinear least squares fits are conducted using the nls and minlm subpackages. nls offers a derivative-free solver, whereas minlm provides a Levenberg-Marquardt algorithm for differentiable models (with both analytic gradients and numerical differentiation). The commercial versions of both solvers are capable of parallel model evaluations, greatly enhancing optimization efficiency.

However, the preferable way of fitting curves/surfaces is to use lsfit subpackage, because it provides more parallelism opportunities due to its awareness about model structure, and because it provides fitting-specific features like errors in parameters.

Special cases

Fitting circle/sphere to data (LS, MC, MI, MZ)

Metrology and geospatial applications often require one to fit a circle or sphere (generally, N-dimensional sphere) to data. Depending on the situation, one may want to perform:

• a least squares (LS) fit, i.e. to minimize sum of squared distances from the sphere surface
• a minimum circumscribed (MC) fit, i.e. to minimize the sphere radius subject to the constraint that all points are located inside the sphere
• a maximum inscribed (MI) fit, i.e. to maximize the sphere radius subject to the constraint that all points are located outside the sphere
• a minimum zone (MZ) fit, i.e. to minimize the difference between circumscribed and inscribed radii

The fitsphere subpackage of the Interpolation package solves all these tasks.

Ramer-Douglas-Peucker piecewise linear fit

The Ramer–Douglas–Peucker algorithm simplifies a curve composed of line segments into a simpler curve with fewer points. ALGLIB offers several functions to perform an Ramer–Douglas–Peucker fit:

• the lstfitpiecewiselinearrdp and lstfitpiecewiselinearrdpfixed functions perform a non-parametric Ramer–Douglas–Peucker fit. 'Non-parametric' in this context means that these functions work with curves which can be represented as Y(X). Figures like circles are impossible to model with these functions.
• the parametricrdpfixed function performs a parametric Ramer–Douglas–Peucker fit

This article is licensed for personal use only.