L2CXCV: A Fortran 77 package for least squares convex/concave data smoothing

Fortran 77 software is given for least squares smoothing to data values contaminated by random errors subject to one sign change in the second divided differences of the smoothed values, where the location of the sign change is also unknown of the optimization problem. A highly useful description of the constraints is that they follow from the assumption of initially increasing and subsequently decreasing rates of change, or vice versa, of the process considered. The underlying algorithm partitions the data into two disjoint sets of adjacent data and calculates the required fit by solving a strictly convex quadratic programming problem for each set. The piecewise linear interpolant to the fit is convex on the first set and concave on the other one. The partition into suitable sets is achieved by a finite iterative algorithm, which is made quite efficient because of the interactions of the quadratic programming problems on consecutive data. The algorithm obtains the solution by employing no more quadratic programming calculations over subranges of data than twice the number of the divided differences constraints. The quadratic programming technique makes use of active sets and takes advantage of a B-spline representation of the smoothed values that allows some efficient updating procedures. The entire code required to implement the method is 2920 Fortran lines. The package has been tested on a variety of data sets and it has performed very efficiently, terminating in an overall number of active set changes over subranges of data that is only proportional to the number of data. The results suggest that the package can be used for very large numbers of data values. Some examples with output are provided to help new users and exhibit certain features of the software. Important applications of the smoothing technique may be found in calculating a sigmoid approximation, which is a common topic in various contexts in applications in disciplines like physics, economics, biology and engineering. Distribution material that includes single and double precision versions of the code, driver programs, technical details of the implementation of the software package and test examples that demonstrate the use of the software is available in an accompanying ASCII file.Program summaryTitle of program:L2CXCVCatalogue identifier:ADXM_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXM_v1_0Program obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandComputer:PC Intel Pentium, Sun Sparc Ultra 5, Hewlett-Packard HP UX 11.0Operating system:WINDOWS 98, 2000, Unix/Solaris 7, Unix/HP UX 11.0Programming language used:FORTRAN 77Memory required to execute with typical data:  O(n)O(n), where n is the number of dataNo. of bits in a byte:8No. of lines in distributed program, including test data, etc.:29 349No. of bytes in distributed program, including test data, etc.:1 276 663No. of processors used:1Has the code been vectorized or parallelized?:noDistribution format:default tar.gzSeparate documentation available:YesNature of physical problem:
• Analysis of processes that show initially increasing and then decreasing rates of change (sigmoid shape), as, for example, in heat curves, reactor stability conditions, evolution curves, photoemission yields, growth models, utility functions, etc.
• Identifying an unknown convex/concave (sigmoid) function from some measurements of its values that contain random errors. Also, identifying the inflection point of this sigmoid function.Method of solution:Univariate data smoothing by minimizing the sum of the squares of the residuals (least squares approximation) subject to the condition that the second order divided differences of the smoothed values change sign at most once. Ideally, this is the number of sign changes in the second derivative of the underlying function. The remarkable property of the smoothed values is that they consist of one separate section of optimal components that give nonnegative second divided differences (convexity) and one separate section of optimal components that give nonpositive second divided differences (concavity). The solution process finds the joint (that is the inflection point estimate of the underlying function) of the sections automatically. The underlying method is iterative, each iteration solving a structured strictly convex quadratic programming problem in order to obtain a convex or a concave section over a subrange of data.Restrictions on the complexity of the problem:Number of data, n  , is not limited in the software package, but is limited to 2000 in the main driver. The total work of the method requires 2n−22n−2 structured quadratic programming calculations over subranges of data, which in practice does not exceed the amount of O(n2)O(n2) computer operations.Typical running times:  CPU time on a PC with an Intel 733 MHz processor operating in Windows 98: About 2 s to smooth n=1000n=1000 noisy measurements that follow the shape of the sine function over one period.Summary:L2CXCV is a package of Fortran 77 subroutines for least squares smoothing to n univariate data values contaminated by random errors subject to one sign change in the second divided differences of the smoothed values, where the location of the sign change is unknown. The piecewise linear interpolant to the smoothed values gives a convex/concave fit to the data. The underlying algorithm is based on the property that in this best convex/concave fit, the convex and the concave section are both optimal and separate. The algorithm is iterative, each iteration solving a strictly convex quadratic programming problem for the best convex fit to the first k   data, starting from the best convex fit to the first k−1k−1 data. By reversing the order and sign of the data, the algorithm obtains the best concave fit to the last n−kn−k data. Then it chooses that k as the optimal position of the required sign change (which defines the inflection point of the fit), if the convex and the concave components to the first k   and the last n−kn−k data, respectively, form a convex/concave vector that gives the least sum of squares of residuals. In effect the algorithm requires at most 2n−22n−2 quadratic programming calculations over subranges of data. The package employs a technique for quadratic programming, which takes advantage of a B  -spline representation of the smoothed values and makes use of some efficient O(k)O(k) updating procedures, where k is the number of data of a subrange. The package has been tested on a variety of data sets and it has performed very efficiently, terminating in an overall number of active set changes that is about n, thus exhibiting quadratic performance in n. The Fortran codes have been designed to minimize the use of computing resources. Attention has been given to computer rounding errors details, which are essential to the robustness of the software package. Numerical examples with output are provided to help the use of the software and exhibit certain features of the method. Distribution material that includes driver programs, technical details of the installation of the package and test examples that demonstrate the use of the software is available in an ASCII file that accompanies this work.