Chapter 4 CURVE FITTING and LEAST SQUARES Summary The curve fitting functions support a full spectrum of techniques for fitting curves to data. The areas covered are: o Rational Approximations o Spline Curve Fits o Least Square Fits Rational approximation functions use a Tchebycheff Pade' method to obtain nearly optimal approximations. Simple cubic and tensioned splines are supported. The least square functions estimate polynomial, general linear, and nonlinear least square fits. Notes on Contents This segment of the library covers methods for approximating functions efficiently, computation of smooth spline approximation curves from a discrete set of values, and a comprehensive set of least squares estimation techniques. o Rational Approximations: The evaluation of elementary functions is typically implemented by using an optimized rational approximation. chcof -------- compute Tchebycheff expansion coefficients. chpade ------- compute a rational approximation. ftch -------- evaluate a rational Tchebycheff approximation. The 'chpade' function fits a rational Tchebycheff Pade' approximation to a specified function. This form is known to yield an approximation very close to the optimal (ie. minimum error) rational approximation. Thus, it provides a computationally efficient technique for generating function values. o Spline Curve Fits: cspl -------- compute a cubic spline interpolation. csplp ------- compute a cubic spline interpolation to a periodic curve. csfit ------- evaluate a cubic spline fit. tnsfit ------ evaluate a "tensioned" spline fit. dcspl ------- evaluate the derivatives of a cubic spline fit. Cubic splines are an excellent class of smooth interpolation functions. These splines are continuous as are their first and second derivatives. The library functions fit both open and closed or periodic curves. Splines in tension can be substituted for cubic splines in cases where the cubic spline interpolation introduces spurious oscillations in the fitted curve. This technique is frequently useful in graphics applications such as contour plotting. o Least Square Fits: Linear Least Squares qrlsq ------ linear least squares using QR reduction. qrvar ------ compute the parameter variance matrix resulting from a QR least squares solution. lsqsv ------ perform a rank analysis on a singular value decomposition of a linear least squares problem. svdlsq ----- compute a SVD analysis of a linear least squares system. sv2lsq ----- compute SVD efficiently for measurements >> number of parameters. The QR reduction is a numerically stable algorithm for computing a least squares solution when the design matrix of the system is known to have full rank. Singular value decomposition is the technique that should be applied when this assumption of full rank breaks down. It provides the maximum amount of information on the character of the system, and ensures a reasonable solution in all cases. Polynomial Least Squares plsq ------ compute a least squares solution using orthogonal polynomials. pplsq ----- compute a least squares fit to a polynomial of specified degree. evpsq ----- evaluate the fit function of an orthogonal polynomial fit. evpsqv ---- evaluate the fit function and rms sigma of an orthogonal polynomial fit. psqcf ----- extract the polynomial coefficient vector of specified degree from an orthogonal fit. psqvar ---- compute the variance matrix of a polynomial coefficient vector. The polynomial least square estimators are based on a numerically stable approach. It returns a set of parameters that support fits with polynomials of any degree up to the maximum specified. It also permits computation of the parameter variance matrix and rms sigmas of fit functions, in cases where this statistical data is desired. Nonlinear Least Squares seqlsq ------- perform a sequential iteration of least square parameter estimation. gnlsq -------- compute a batch Gauss-Newton least square fit from specified initial parameter values. fitval ------- evaluate the value and rms sigma of a nonlinear least square fit function. The library supports a novel mix of sequential and batch mode (Gauss- Newton) estimation of nonlinear least-square fits. The sequential fit is highly effective in the initial parameter search phase of nonlinear least squares estimation. A procedure that starts with a few sequential search passes and refines the resulting estimate using the Gauss-Newton estimator, is highly effective in nonlinear least square fits. ------------------------------------------------------------------------------- General Technical Comments The availability of a nearly optimal rational approximation for function evaluation is an important feature of the library. Efficient computation can be based on the use of this approximation so that the code implemented for initial evaluation can be biased towards accuracy and range of validity. Use of this technique is recommended for any computationally intensive application that requires frequent evaluations of a function. The main novelty in this segment is the implementation of a sequential estimator for nonlinear least squares problems. Experience indicates that this represents a powerful approach to searching parameter space in these very difficult problems. The use of adaptive least squares techniques goes back to Gauss, and it is also featured in the Kalman-Bucy filters employed in modern control theory. Nonlinear estimation is yet another domain where this approach is recommended. ------------------------------------------------------------------------------ FUNCTION SYNOPSES ------------------------------------------------------------------------------ Tchebycheff Polynomial Approximations: ------------------------------------------------------------------------------ chcof Compute the coefficients of a Tchebycheff expansion of the function func(x). (These coefficients can be used as input to chpade.) chcof(double c[],int m,double (*func)()) c = array containing m+1 values of computed coefficients m = maximum degree of fit (polynomials 0 --- m used) func = pointer to user supplied function ( Range of (*func)(x) evaluation -1 <= x <= 1. ) The Tchebycheff expansion is given by f = c[0]/2 + Sum(i=1 to m) c[i]*Ti(x) , where Ti(x) is the ith Tchebycheff polynomial. ----------------------------------------------------------- chpade Compute the coefficients of a rational Tchebycheff approximation. chpade(double c[],double a[],int m,double b[],int n) c = array of dimension m+2*n+2 containing the coefficients of a Tchebycheff series expansion of the function f a = m+1 dimensional array of numerator coefficients m = degree of numerator polynomial b = n+1 dimensional array of denominator coefficients ( b[0]=1.0 , by convention ) n = degree of denominator polynomial Here the input expansion gives f = c[0]/2 + Sum(i=1 to N) c[i]*Ti(x) , with N=m+2n+1. The output approximation to is: f = {Sum(i=0 to m) a[i]*Ti(x)}/{1+Sum(j=1 to m) b[i]Ti(x))} . ------------------------------------------------------------------- ftch Evaluate a rational Tchebycheff Pade approximation. double ftch(double x,double a[],int m,double b[],int n) x = value of independent variable ( -1 <=x <= 1 ) a = m+1 dimensional array of numerator coefficients m = degree of numerator polynomial b = n+1 dimensional array of denominator coefficients n = degree of denominator polynomial return value: f = rational approximation of f(x), with f = {Sum(i=0 to m) a[i]*Ti(x)}/{1+Sum(j=1 to m) b[i]Ti(x))} . ----------------------------------------------------------------------------- Spline Curve Fits: ----------------------------------------------------------------------------- cspl Compute a cubic or tensioned spline fit to an open curve. cspl(double x[],double y[],double z[],int m,double tn) x = array containing m+1 x-coordinates of the curve y = array containing m+1 y-coordinates of the curve z = output array of m+1 second derivatives ( the boundary condition assumes z[0]=z[m]=0. ) m = number of intervals (m+1 input points) tn = tension parameter, with: tn = 0. -> cubic spline tn > 0. -> spline in tension ----------------------------------------------------------------- csplp Compute a cubic or tensioned spline approximation to a closed curve. csplp(double x[],double y[],double z[],int m,double tn) x = array containing m+1 x-coordinates of the curve y = array containing m+1 y-coordinates of the curve z = output array of m+1 second derivatives ( The boundary condition is periodic, with equal first and second derivatives at x[0] and x[m]. ) m = number of intervals (m+1 input points) tn = tension parameter, with: tn = 0. -> cubic spline tn > 0. -> spline in tension --------------------------------------------------------------- csfit Evaluate a cubic spline curve at any interior point. double csfit(double w,double x[],double y[],double z[],int m) w = value of the independent variable x,y,z = m+1 dimensional arrays define a cubic spline m = number of intervals return value: cs(w) = cubic spline interpolation if x[0]<=w<=x[m] ( 0.0 is returned for wx[m] ) This function is used to evaluate splines computed with tn=0.0. ------------------------------------------------------------------ dcspl Evaluate the derivative of a cubic spline fit at interior points. double dcspl(double x,double u[],double v[],double z[],int m) x = value of the independent variable u,v,z = m+1 dimensional arrays define a cubic spline m = number of intervals return value: dcs/dx = derivative of cubic spline interpolation if u[0] <= x <=u[m] ( 0.0 is returned for x < u[0] or x > u[m] ) -------------------------------------------------------------------- tnsfit Evaluate a tensioned spline curve at any interior point. double tnsfit(double w,double x[],double y[],double z[],int m,double tn) w = value of independent variable x,y,z = m+1 dimensional arrays of tensioned spline coefficients m = number of intervals tn = spline tension parameter (should match value tn > 0.0 used in computing the spline parameters) return value: tns(w) = value of tensioned spline interpolation when x[0] <= w <= x[m] ( 0.0 is returned if w < x[0] or w > x[m] ) Increasing the tension parameter from tn=0 produces a sequence of curve fits that varies continuously from the cubic spline (tn=0) to a linear interpolation between anchor points (tn >> 1). ------------------------------------------------------------------------------ Linear Least Squares: ------------------------------------------------------------------------------ qrlsq Compute a linear least squares solution for A*x = b using a QR reduction of the matrix A. double qrlsq(double *a,double *b,int m,int n,int *f) a = pointer to m by n design matrix array A This is altered to upper right triangular form by the computation. b = pointer to array of measurement values b The first n components of b are overloaded by the solution vector x. m = number of measurements (dim(b)=m) n = number of least squares parameters (dim(x)=n) (dim(a)=m*n) f = pointer to store of status flag, with *f = 0 -> solution valid *f = -1 -> rank of A < n (no solution) return value: ssq = sum of squared fit residuals The QR algorithm employs an orthogonal transformation to reduce the matrix A to upper right triangular form, with A = Q*R and R = Q~*A . The matrix R has non-zero components confined to the range 0 <= i <= j < n. The system vector b is transformed to b' = U~*b and Sum(k=j to n-1) R[j,k]*x[k] = b'*[j] is solved for x, using 'solvru' (see Chapter 1). The sum of squared residuals is ssq = Sum(i=n to m-1){b~[i]^2} . ---------------------------------------------------------------- qrvar Compute the parameter variance matrix V for QR least squares. double qrvar(double *a,int m,int n,double ssq) a = pointer to a-matrix output of qrlsq Only the n by n upper right triangular portion R is used. The first n rows are replaced by the parameter variance matrix V. m = number of measurements (m > n assumed) n = number of parameters (dim(a)>=n*n) ssq = sum of squared fit residuals return value: s = ssq/(m-n), the estimated measurement variance The parameter variance matrix is computed by assuming that the measurements are independent, with measurement variance sigsq = ssq/(m-n) and matrix V given by V = sigsq*Inv(R~*R) , with R an upper right triangular matrix whose non-zero elements are equal to corresponding elements of the matrix A output by qrlsq. -------------------------------------------------------------------- The singular value decomposition of the design matrix A can be used to derive a least squares solution in cases where the QR method breaks down because A does not have full rank. In this sense, SVD is a sure-fire approach to linear least squares. The method is based on the decomposition A = U*D*V~ , with U and V orthogonal matrices and D an m by n matrix (m>=n) whose non-zero elements are D[i,i]=d[i] for i=0 to n-1. The resulting least squares solution for x is given by x = V*Inv(D)*U~*b for x[i] with i=0 to n-1. The function 'lsqsv' employs a threshold test to identify small singular values that correspond to rank deficiency in the design matrix. This threshold is specified as a user input. When D[i,i] < th , the corresponding inverse element Inv(D)[i,i] is replaced by zero, so that the ith measurement component is not used in the solution. ------------------------------------------------------------------ lsqsv Generate a SVD based solution for the linear least squares system A*x = b. double lsqsv(double *x,int *pr,double *var,double *d,double *b, \ double *v,int m,int n,double th) x = pointer to array to be loaded with least squares parameters x pr = pointer to store for rank of system solution r (r<=n, and if r normal exit -1 -> input error m < n This form of least squares SVD is designed for use when the number of measurements m is comparable to the number of parameters n. The quantities returned are the orthogonal matrix V, the modified system vector b~, and the singular values d[i]. ------------------------------------------------------------------- sv2lsq Perform a singular value decomposition ot a linear least squares system efficiently when m >> n. int sv2lsq(d,a,b,m,v,n) double *d,*a,*b,*v; int m,n; d = pointer to array of computed singular values a = pointer to array of design matrix A This matrix is altered by the computation. b = pointer to array of measurement values b The computation replaces b by the vector b'=U~*b. m = number of measurements (dim(b)=m) n = number of fit parameters x (dim(d)=n,dim(a)=m*n) return value: status flag 0 -> normal exit -1 -> input error m < n This algorithm is more efficient for SVD when m > (3/2)n, and it should certainly be used in place of 'svdlsq' when m > 2n. -------------------------------------------------------------------- Auxiliary function for linear least squares: Perform the QR part of the SVD least squares reduction. (Called by both svdlsq and sv2lsq.) int qrbdbv(double *d,double *e,double *b,double *v,int n) d = pointer to array of diagonal elements The computation loads this array with singular values. e = pointer to array with superdiagonal elements in the first n-1 positions. These values are altered. b = pointer to array of system vector b This vector is transformed by the QR rotations. v = pointer to orthogonal transformation matrix V This matrix is transformed by the QR rotations. n = dimension (dim(d)=dim(e)=n,dim(v)=n*n,dim(b)>=n) return value: N = number of QR iterations required ------------------------------------------------------------------------------- Polynomial Least Squares: ------------------------------------------------------------------------------- Polynomial fits are one of the most common forms of least squares analysis. A comprehensive set of functions, based on developing a set of polynomials orthogonal on the set of measurements, is provided to support this form of analysis. The orthogonal polynomials op(x) satisfy the orthogonality relation Sum(i=0 to n-1){opj(x[i])*opk(x[i])} = h[k] if j=k and 0 otherwise. Orthogonal polynomials are computed using the recurrence op{k+1}(x) = (x - f[k]*)opk(x)-(h[k]/h[k-1])*op{k-1}(x) , where op{-1} = 0 and op0 =1. The coefficients are determined by h[k] = Sum(i=0 to n-1){(opk(x[i]))^2} and f[k] = Sum(i=0 to n-1){x[i]*opk(x[i])^2} . The fit residuals e[i] are given by e[i] = y[i] - Sum(i=0 to m-1) c[k]*opk(x[i]) for a fit of degree n in x. The sum of squared residuals ssq = Sum(i=0 to n-1) e[i]^2 provides a measure of the fit, and can be used to estimate measurement variance in cases where this statistic makes sense. The variance matrix for the coefficients c[k] of the least squares polynomials takes the form var(c[j]) = ssq/(h[j]*(n-m-1)) for an orthogonal polynomial fit of degree m. ------------------------------------------------------------------------ plsq Compute a polynomial least squares fits to data, using polynomials orthogonal on the x-values. plsq(double *x,double *y,int n,Opol *cf,double *ssq,int m) x = pointer to array of x-values y = pointer to array of measurements converted to array of fit residuals by the computation n = number of measurements (dim(x)=dim(y)= n) cf = pointer to array of structures specifying the orthogonal polynomials ssq = pointer to array of squared residuals for fits of order 1 to m m = maximum order of fit desired (dim(cf)=dim(ssq)=m) (Fits with polynomials of degree 0 to m-1 are computed by this function.) The structure Opol is defined in the header file ccmath.h by the following line. typedef struct orpol {double cf,hs,df;} Opol; This structure holds the coefficient of the polynomial in the least squares fit, and the recurrence coefficients needed to generate the polynomial. cf[k].cf = c[k], cf[k].hs = h[k], and cf[k].df = f[k] . The least squares fit for a polynomial of degree m is given by the expansion fit_n(x) = Sum(k=0 to n){c[k]*opk(x)} for n = 0,1, - - ,m-1 . -------------------------------------------------------------------- pplsq Compute a least squares fit to a polynomial of degree m-1. double pplsq(double *x,double *y,int n,double *bc,int m) x = pointer to array of x-values of measurements y = pointer to array of measurements converted to array of fit residuals by the computation n = number of measurements (dim(x)=dim(y)= n) bc = pointer to vector of polynomial coefficients m = order of fit (dim(bc)= m) (The fit polynomial is of degree m-1 in x.) This function is used in cases where the fit is to be made with a polynomial of known degree, and the output is to be expressed in terms of coefficients of powers of x. fit(x) = Sum(i=0 to m-1) bc[i]*x^i . --------------------------------------------------------------- evpsq Evaluate a computed least squares polynomial fit function. double evpsq(double x,Opol *c,int m) x = argument where fit function is to be evaluated c = array of orthogonal polynomial coefficient structures m = order of fit (dim(c)=m) return value: y = fit(x) The fit functions generated by a call of 'plsq' are evaluated by this function, with fit(x) = Sum(i=0 to m-1) c[k]*opk(x) . --------------------------------------------------------------- evpsqv Evaluate a least squares polynomial fit function and compute the standard deviation of the fit. double evpsqv(double x,Opol *c,int m,double *v,double sig) x = argument where fit function is to be evaluated c = array of orthogonal polynomial coefficient structures sig = estimate of sigma squared for the measurements v = pointer to store for rms sigma of returned value ( y = fit(x) , v ) m = order of fit (dim(c)=m) return value: y = fit(x) The input sig is an estimate of the variance of individual measurements, based on the sum of squared fit residuals. sig = ssq(m)/(n-m-1) . The accuracy of the fit at a point is expressed as y(x) = fit(x) +- *v , with *v the one sigma rms variance. ---------------------------------------------------------------- psqcf Compute the polynomial coefficient vector for a fit by a polynomial of degree m-1. psqcf(double *b,Opol *c,int m) b = pointer to array for polynomial coefficient output c = pointer to array of structures specifying the polynomials orthogonal on x-values m = order of fit (dim(cf)=dim(b)=m) This function expresses a fit obtained with 'plsq' in terms of the coefficients of powers of x, with fit[m-1](x) = Sum(i=0 to m-1) b[i]*x^i . ---------------------------------------------------------------- psqvar Compute the variance matrix of the polynomial coefficients. psqvar(double *v,double sig,Opol *c,int m) v = pointer to array for variance matrix output sig = estimate of sigma squared for the measurements c = pointer to array of structures specifying the orthogonal polynomials m = order of fit (dim(c)= m, dim(v)=m*m) This function extracts the variance matrix of the coefficients of powers of x from a least squares fit computed by 'plsq'. The input sig is an estimated measurement variance sig = ssq(m)/(n-m-1) and the variance matrix V = b*b~ , where the components of the the vector b are computed by 'psqcf'. ------------------------------------------------------------------------------- Nonlinear Least Squares: ------------------------------------------------------------------------------- The basic problem encountered in nonlinear least squares is to minimize ssq = Sum(i=0 to n-1){y[i] - f(x[i],b[1] - b[m])}^2 , where the fit function f is not linear in the parameters b[k]. The range of possible nonlinear functions is so large that a completely general solution to this problem is unlikely. Nevertheless, experience with a combination of sequential and batch estimation steps indicates that this approach is suitable and efficient for a wide range of nonlinear estimation problems. The sequential method is one that was originally applied by Gauss to the reduction of astronomical orbit data. Each measurement in a sequence is used to update the current parameter estimate. These updates employ a linearization of the nonlinear system about the current parameter estimate and an error term of the form e[i] = y[i] - f(x[i],B'[i-1]) , where B'[i-1] is the parameter vector estimated at the previous update. This algorithm is very effective in searching parameter space for a parameter vector near the minimum of ssq. ------------------------------------------------------------------------------- seqlsq Perform a sequential iteration of a least square fit to the function y=(*func)(x,par). double seqlsq(double x[],double y[],int n,double par[],double var[],\ int m,double de,double (*func)(),int kf) x = n dimensional array of independent variables y = n dimensional array of dependent variables n = number of points fitted m = dimension of parameter vector par = m-vector of fit function parameters (The input values represent an initial estimate.) var = array containing the m by m matrix proportional to the parameter variance matrix (constant of proportionality = ssq/(n-m) ) de = interval for derivative computation (dfunc/dpar) func = pointer to the fit function (y=(*func)(x,par) ) kf = control flag, with: kf=0 -> initialize var as a unit matrix kf>0 -> accept input values of var matrix return value: ssq = sum of squares at fit points This function can be used, in conjunction with the Gauss-Newton estimator gnlsq, to estimate nonlinear least square fits. Sequential iterations are effective in the initial search stage. The parameter vector estimate is updated for each measurement y[i] in sequence, ssq = Sum(i=0 to n-1) (y[i] - f(x[i],par))^2 should be interpreted with this sequential variation of par in mind. Iterations of this function are used to search parameter space. Normally the output of an iteration is employed as the input to the next estimation stage. Several other options are outlined below. Termination of a sequence of sequential estimation passes should be based on the statistical significance of the decrease in ssq in successive passes. Input options: 1) The use of a unit matrix for initializing var is usually acceptable. 2) Parameter scale differences in components of par can be addressed by using an initial pass through the Gauss_Newton estimator to obtain an input var. 3) The output matrix var can be multiplied by a scale factor, with s >1 between iterations. This is used to avoid a search with variance too 'stiff' to support significant changes in the parameter vector. Numerical derivatives of the fit function with respect to the parameters are used. The interval used to compute these derivatives should be selected to ensure that (f(x,par'')-f(x,par'))/de , with de = par''[i]-par'[i] yields an effective approximation to the partial derivative. This also applies to the function 'gnlsq'. ---------------------------------------------------------------------- The Gauss-Newton least squares estimator employs a batch process to estimate a parameter vector, with the linearization performed at the input point for each measurement. This estimator converges quite well when the starting point in parameter space is fairly close to a minimum. Thus, it can be used effectively in the final iterations of a nonlinear fit. The output of the previous sequential or batch pass is used as input to 'gnlsq'. gnlsq Perform a Gauss-Newton iteration of a least square fit to the function y=(*func)(x,par). double gnlsq(double x[],double y[],int n,double par[],double var[],\ int m,double de,double (*func)()) x = n dimensional array of independent variables y = n dimensional array of dependent variables n = number of points fitted m = dimension of parameter vector par = m-vector of fit function parameters (The input values represent an initial estimate.) var = array containing the m by m matrix proportional to the parameter variance matrix (constant of proportionality = ssq/(n-m) ) de = interval for derivative computation (dfunc/dpar) func = pointer to the fit function (y=(*func)(x,par) ) return value: ssq = sum of squares at fit points ssq > 0. -> normal exit ssq = -1.0 -> singular matrix var detected (ie. an input error) A single iteration of this function can be used for well-posed linear least square estimates, where it is equivalent to the use of the normal equation method. In nonlinear problems it is most effective in the final iterations, however, an initial iteration is sometimes used to determine the relative scaling of the variance matrix. The sum of squared residuals is given by ssq = Sum(i=0 to n-1) (y[i] - f(x[i],par))^2 , where the initial input parameter vector par is used at each step. The estimated measurement variance, based on this sum, is sigsq = ssq/(n-m). This value should be used to scale the final variance matrix output var when the statistical uncertainty in parameter estimates is of interest. Termination of a sequence of Gauss-Newton iterations is best based on the magnitude of changes in the estimated parameter vector from pass to pass. ---------------------------------------------------------------- fitval Evaluate the value and rms sigma of a general least square fit at x. double fitval(double x,double *s,double *par,double (*fun)(),\ double *v,int n) x = value of independent variable s = pointer to store for computed fit's rms sigma at x par = pointer to estimated fit parameter vector fun = pointer to fit function (y=(*fun)(x,par)) v = pointer to estimated fit variance matrix n = dimension of parameter vector return value: f(x,par) = value of the estimated fit function at x The variance matrix input to fitval must be scaled by the estimated s2 of the fit so that a valid rms sigma can returned in *s. y = f(x,par) +- *s . setfval Allocate/free the storage required by fitval. void setfval(int i,int n) int i,n; i = control flag, with i=0 -> allocate internal storage i=1 -> free storage n = dimension of fit parameter vector The function 'setfval' must be called with i=0 before any calls to 'fitval' are made.