Nonlinear Data Fitting by Total Least Squares Regression

This example illustrates a general approach to nonlinear data fitting by finding the best ellipse through a set of two-dimensional data points.  The desired ellipse is formulated as the solution of a nonlinear constrained optimization problem, and solved with the KNITRO for Mathematica application package.  The KNITRO solver, developed by Ziena Optimization, Inc. (www.ziena.com), is the premier algorithm for large-scale nonlinear optimization, and readily handles problems with thousands of unknowns.
This note shows how to formulate the regression problem in Mathematica, and solve it numerically using the KNITRO package. The regression technique easily generalizes to multi-dimensional data points and arbitrary fitting functions.
Suppose we have the 300 sample data points shown in Figure 1 (they were generated randomly for this example).  We wish to find the "best" ellipse that fits this data.

The ellipse can be specified in many ways, but we choose the general quadratic form  a (x^2) + b (y^2) + c x y  +  d x  +  
          e y  =  1, where {a,b,c,d,e} are unknown parameters to be determined by optimization.
Collect the N data points as a list of pairs:  dataPts = {{dataX[1],dataY[1]}, {dataX[2],dataY[2]}, ... }.  For each data point, define a "shadow" point with coordinates {shX[i],shY[i]}.  The data fitting problem is to find shadow points that are closest to their corresponding data point, such that all shadow points lie on a single quadratic curve.  Both the shadow points and the quadratic parameters are unknowns.  Formally, we write

\(TraditionalForm\`\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(min\ \ \ \ \ \ \ \(∑\+\(i = 1\)\%N\((shX\_i - \ dataX\_i)\)\^2\)\  + \((shY\_i - \ dataY\_i)\)\^2

    	subject\ \ to\ \ \ \ \ \ a\ shX\_i\^2 + b\ shY\_i\^2 + 
        c\ shX\_i\ shY\_i + \ d\ shX\_i + e\ shY\_i = 1\ \ , \ \ i = 
      1.  . N\)\)\)

This optimization problem has 2N+5 unknowns, a quadratic objective function, and N quadratic equality constraints.  We initialize  the shadow point unknowns to equal the data points and call KNITRO to compute a solution.

Here is Mathematica code that constructs the objective function  f and a list of N equality constraints c (the quadratic curve unknowns {a,b,c,d,e} are called {aElps,bElps,...}):
    Table[shX[i], {i,n}];  Table[shY[i], {i,n}];
    f = Sum[(shX[i] - dataX[[i]])^2 + (shY[i] - dataY[[i]])^2, {i,1,n}];
    c = Reap[ Do[ Sow[ aElps shX[i]^2 + bElps shY[i]^2 + cElps shX[i] shY[i]
                       + dElps shX[i] + eElps shY[i] == 1 ], {i, 1, n} ] ];
    c = Last[Last[c]];
Here is Mathematica code to assign initial values and call KNITRO:
    xWithInit = Last[Last[ Reap[ Do[ Sow[ {shX[i], dataX[[i]]} ], {i, 1, n} ] ] ]];
    yWithInit = Last[Last[ Reap[ Do[ Sow[ {shY[i], dataY[[i]]} ], {i, 1, n} ] ] ]];
    varsWithInit = Join[ xWithInit, yWithInit,
                         {{aElps, 1},{bElps, 1},{cElps, 1},{dElps, 1},{eElps, 1}}];
    result = KnitroMinimize[{f,c},{},varsWithInit];

KNITRO solves the problem and returns the optimal ellipse parameters.  The plot below shows where the best ellipse lies.

The form of the unknown ellipse,  a (x^2) + b (y^2) + c x y  +  d x  +  
          e y  =  1, is general enough to admit other quadratic curves.  For example, a different set of data points might be fit best by values of {a,b,c,d,e} that describe a hyperbola (see http://mathworld.wolfram.com/Ellipse.html).  KNITRO will have no difficulty finding the best hyperbola with this same problem formulation.
In general, the optimization problem formulated in this example may have many local minima; that is, points which satisfy the constraints and minimize the objective function in a local neighborhood of (R^(2  N + 5)).  Of course we want the best local minimum (the global minimum), but global optimization of nonconvex problems can be extremely time consuming when the number of unknowns is large.  Instead, we make a good initial guess and let KNITRO solve for a nearby local minimum.  The initial guess is to set the shadow point unknowns equal to the data points.  Initial values for the quadratic curve parameters are less significant, and they are simply set to one.

The KNITRO solver has many tuning parameters, described in the user's manual, but in this example we used the default parameter settings.  To illustrate the computational power of KNITRO on large-scale problems, solve the same regression problem with increasing size N, again using default settings.  Measure the CPU time taken by KnitroMinimize[], which does not count the time taken to generate or display data points.   (Times are for Mathematica 5.2 on Linux, running an AMD Athlon 32-bit processor at 2.1 GHz, with 2G memory.)

N (number data points)   Total CPU time   KNITRO CPU time   Evaluation CPU time
           300                4.5 sec           0.1 sec             1.2 sec
           600               12   sec           0.3 sec             2.6 sec
          1200               44   sec           0.6 sec             4.9 sec
          2400              160   sec           1.5 sec            12   sec

Total CPU time is measured with the Mathematica function Timing[].  Computation time is divided between (1) creation of Mathematica NumericalFunction objects, (2) KNITRO solver calculations, and (3) evaluation of numerical partial derivatives in Mathematica.  KNITRO reaches a local minimum by iteration, and requires the numerical value of first and second partial derivatives at each iterate.  Mathematica supplies this information from NumericalFunction objects.  Both KNITRO and Mathematica use sparse arrays for the partial derivatives; hence, computational effort increases relatively slowly with problem size N.  The CPU time for activities (2) and (3) are measured by including KnitroOptionList->{{"log_knml",1}} as an option to the KnitroMinimize[] call.  In this example, most of the computation time is spent creating NumericalObjects in Mathematica; the fraction of time spent in KNITRO is small.

A full listing of the Mathematica code used for this example is available on request from info@ziena.com.  For more information about KNITRO, go to http://www.ziena.com/

Copyright (c) 2006 by Ziena Optimization, Inc.  All rights reserved.