Abstract
Fitting geometric shapes to observed data points (images) is a popular task in computer vision and modern statistics (errors-in-variables regression). We investigate the problem of existence of the best fit using geometric and probabilistic approaches.
1. Introduction
In many areas of human practice, one needs to approximate a set of points (representing experimental data or observations) by a simple geometric shape, such as a line, plane, circular arc, elliptic arc, and spherical surface. This problem is known as fitting a model object (line, circle, sphere, etc.) to observed data points.
The best fit is achieved when the geometric distances from the given points to the model object are minimized, in the least squares sense. Finding the best fit reduces to the minimization of the objective function where denotes the model object (line, circle, sphere, etc.).
While other fitting criteria are used as well, it is the minimum of (1.1) that is universally recognized as the most desirable solution of the fitting problem. It has been adopted as a standard in industrial applications [1, 2]. In statistics, the minimum of (1.1) corresponds to the maximum likelihood estimate under usual assumptions known as functional model; see [3, 4]. Minimization of (1.1) is also called orthogonal distance regression (ODR).
Geometric fitting in the sense of minimizing (1.1) has a long history. First publications on fitting linear objects (lines, planes, etc.) to given points in 2D and 3D date back to the late nineteenth century [5]. The problem of fitting lines and planes was solved analytically in the 1870s [6], and the statistical properties of the resulting estimates were studies throughout the twenties century [7β9], the most important discoveries being made in the 1970s [3, 10, 11]. These studies gave rise to a new branch of mathematical statistics now known as errors-in-variable (EIV) regression analysis [3, 4].
Since the 1950s, fitting circles and spheres to data points became popular in archaeology, industry, and so forth, [12, 13]. In the 1970s, researchers started fitting ellipses and hyperbolas to data points [14β17].
The interest to the problem of geometric fitting surged in the 1990s when it became an agenda issue for the rapidly growing computer science community. Indeed, fitting simple contours to digitized images is one of the basic tasks in pattern recognition and computer vision. In most cases, those contours are lines, circles, ellipses and other conic sections (called simply conics, for brevity) in 2D, and planes, spheres, ellipsoids in 3D. More complicated curves and surfaces are used occasionally too [1, 18], but it is more common to divide complex images into smaller segments and approximate each of those by a line or by a circular arc. This way one can approximate a complex contour by a polygonal line or a circular spline (see, e.g., [19, 20]).
Most publications on the geometric fitting problem are devoted to analytic solutions (which are only possible for lines and planes), or practical algorithms for finding the best-fitting object, or statistical properties of the resulting estimates. Very rarely one addresses fundamental issues such as the existence and uniqueness of the best fit. If these issues do come up, one either assumes that the best fit exists and is unique, or just points out examples to the contrary without deep investigation.
In this paper we address the issue of existence of the best fit in a general and comprehensive manner. The issue of uniqueness will be treated in a separate paper. These issues turn out to be quite nontrivial and lead to unexpected conclusions. As a glimpse of our results, we provide Table 1 summarizing the state of affairs in the problem of fitting 2D objects (here Yes means the best-fitting object exists or is unique in all respective cases; No means the existence/uniqueness fails in some of the respective cases).
We see that the existence and uniqueness of the best-fitting object cannot be just taken for granted. Actually 2/3 of the answers in Table 1 are negative. The uniqueness can never be guaranteed. The existence is guaranteed only for lines. In typical cases (i.e., for typical sets of data points, in probabilistic sense; see precise definition in Section 6), the best-fitting line and circle do exist, but the best-fitting ellipse does not.
The existence and uniqueness of the best fit are not only of theoretical interest but also practically relevant. For example, knowing under what conditions the problem does not have a solution might help us understand why the computer algorithm keeps diverging, or returns nonsense, or crashes altogether. While the cases where the best fit does not theoretically exist may be exceptional, nearby cases may be practically hard to handle, as the best-fitting object may be extremely difficult to find.
The nonuniqueness has its practical implications too. It means that the best-fitting object may not be stable under slight perturbations of the data points. An example is described by Nievergelt [21]; he presented a set of points that can be fitted by three different circles equally well. Then by arbitrarily small changes in the coordinates of the points one can make any of these three circles fit the points a bit better than the other two circles, thus the best-fitting circle will change abruptly. See also a similar example in [4, Section 2.2].
Here we develop a general approach to the studies of existence of the best fit. Our motivation primarily comes from popular applications where one fits lines, circles, ellipses, and other conics, but our methods and ideas can be applied to much more general models. Our approach works in any dimension; that is, we can treat data points and model objects for any , but for the sake of notational simplicity and ease of illustrations we mostly restrict the exposition to the planar case .
2. W-Convergence and the Induced Topology
A crucial notion in our analysis is that of convergence (of model objects).
Motivating Example
Consider a sequence of lines . They all pass through the origin and their slopes decrease as grows (Figure 1). Naturally, we would consider this sequence as convergent; it converges to the -axis, that is, to the line , as .
However, it is not easy to define convergence so that the above sequence of lines would satisfy it. For example, we may try to use Hausdorff distance to measure how close two objects are. It is defined by If the two sets are closed and the Hausdorff distance between them is zero, that is, , then they coincide . If the Hausdorff distance is small, the two sets nearly coincide with each other. So the Hausdorff distance seems appropriate for our purposes.
But it turns out that the Hausdorff distance between the lines and in our example is infinite, that is, for every . Thus the Hausdorff distance cannot be used to characterize convergence of lines as we would like to see it.
Windows
So why do we think that the line converges to the line , despite an infinite Hausdorff distance between them? It is because we do not really see an infinite line, we only βseeββ objects in a certain finite area, like in Figure 1. Suppose we see objects in some rectangle
which for the moment will play the role of our βwindowββ through which we look at the plane. Then we see segments of our lines within ; that is, we see intersections and . Now clearly the segment gets closer to , as grows, and in the limit they become identical. This is why we see the lines converging to . We see this convergence no matter how large the window is. Note that the Hausdorff distance between and indeed converges to zero: as .
Window-Restricted Hausdorff Distance
Taking cue from the above example, we define the Hausdorff distance between sets and ββwithin a finite window as follows:
The formula (2.3) applies whenever both sets, and , intersect the window . If only one set, say , intersects , we modify (2.3) as follows:
A similar modification is used if only intersects . If neither set intersects the window , we simply set (because we βseeββ two empty sets, which are not distinguishable).
We now define our notion of convergence for a sequence of sets.
Definition 2.1. Let be some sets and another set. We say that the sequence ββconverges to if for any finite window we have
In this definition the use of finite windows is crucial; hence we will call the resulting notion Window convergence, or W convergence, for short. According to this definition, the sequence of lines in the above example indeed converges to the limit line .
For finite-size (bounded) objects, like circles or ellipses, the W-convergence is equivalent to the convergence with respect to the Hausdorff distance. However for unbounded objects, such as lines and hyperbolas, we have to use our window-restricted Hausdorff distance and formula (2.5).
Our definition of convergence (i.e., W-convergence) is intuitively clear, but some complications may arise if it is used too widely. For example, if the limit object is not necessarily closed, then the same sequence may have more than one limit (though the closure of every limit set will be the same). Thus, from now on, to avoid pathological situations, we will assume that all our sets are closed. All the standard model objectsβlines, circles, ellipses, hyperbolasβare closed.
Convergence induces topology in the collection of all closed sets , which we denote by . A set is closed if for any sequence of sets converging to a limit set the limit set also belongs to , that is, . A set is open if its complement is closed.
The above topology on is metrizable. This means we can define a metric on , that is, a distance between closed sets in , in such a way that the W-convergence means exactly that the distance between and goes down to zero. Such a distance can be defined in many ways; here is one of them.
Definition 2.2. The W distance (or Window distance) between two closed sets is defined as follows: where is a square window of size , that is, .
In this formula, we use a growing sequence of nested windows and the Hausdorff distances between and within those windows balanced by the factors . In the formula (2.6), the first nonzero term corresponds to the smallest window that intersects at least one of the two sets, or . The sum in (2.6) is always finite. Indeed, let us suppose, for simplicity, that both and intersect each window . Then (since the distance between any two points in is at most ) the above sum is bounded by .
3. Continuity of Objective Function
Recall finding the best-fitting object for a given set of points consists of minimization of the objective function (1.1).
Let denote the collection of model objects. We put no restrictions on other than all are assumed to be closed sets. The collection is then a subset of the space of all closed subsets . The topology and metric on defined above now induce a topology and metric on .
Redundancy Principle
If an object is a subset of another object , that is, , then for any point we have ; thus cannot fit any set of data points better than does. So for the purpose of minimizing , that is, finding the best-fitting object, we may ignore all model objects that are proper subsets of other model objects (they are redundant). This may reduce the collection somewhat. Such a reduction is not necessary, it is just a matter of convenience, and we will occasionally apply it below.
Conversely, there is no harm in considering any subset of an object as a (smaller) object, too. Indeed, if provides a best fit (i.e., minimizes the objective function ), then so does , because . Hence including into the collection will not really be an extension of , its inclusion will not change the best fit.
Theorem 3.1 (continuity of ). For any given points and any collection of model objects the function defined by (1.1) that is continuous on . This means that if a sequence of objects converges (i.e., W converges) to another object , then .
Proof. Since is the sum of squares of distances to individual points , see (1.1), it is enough to verify that the distance is a continuous function of for any given point .
Suppose and a sequence of closed sets W converging to a closed set . We denote by the point in closest to , that is, such that ; such a point exists because is closed. Denote by the disk centered on of radius ; it contains . Let be a window containing the disk .
Since contains , it intersects with , that is, . This guarantees that . Thus, there are points such that . Since , we conclude that the upper limit of the sequence does not exceed , that is,
On the other hand, we will show that the lower limit of the sequence cannot be smaller than , that is,
The estimates (3.1) and (3.2) together imply that , as desired, and hence the distance function will be continuous on . It remains to prove (3.2).
To prove (3.2) we assume that it is false. Then there is a subsequence in our sequence of sets such that
Denote by the point in closest to , that is, such that . Then we have
Since the points are closer to than the point is, we have . Recall that , hence
Denote by the point in closest to , that is, such that . Now we have, by triangle inequality,
The limit of the first term on the right hand side of (3.6) is by (3.4), and the limit of the second term is zero by (3.5). This implies , which is absurd. The contradiction proves (3.2). And the proof of (3.2) completes the proof of the theorem.
4. Existence of the Best Fit
The best-fitting object corresponds to the (global) minimum of the objective function , that is, . The function defined by (1.1) cannot be negative; thus, it always has an infimum This means that one cannot find objects such that , but one can find objects such that is arbitrarily close to . More precisely, there is a sequence of objects such that for each and as .
In practical terms, one usually runs a computer algorithm that executes a certain iterative procedure. It produces a sequence of objects (here denotes the iteration number) such that ; that is, the quality of approximations improves at every step. If the procedure is successful, the value converges to the minimal possible value, , and the sequence of objects converges (i.e., W-converges) to some limit object . Then the continuity of the objective function (proven earlier) guarantees that ; that is, indeed provides the global minimum of the objective function, so it is the best-fitting object .
A problem arises if the limit object does not belong to the given collection ; hence it is not admissible. Then we end up with a sequence of objects , each of which fits (approximates) the given points better than the previous one, but not as good as the next one. None of them would qualify as the best fit, and in fact the best fit would not exist, so the fitting problem would have no solution.
Illustrative Example
Suppose that our model objects are circles, and our given points are , and . Then the sequence of circles defined by will fit the points progressively better (tighter) as grows, so that as (Figure 2). On the other hand, no circle can pass through three collinear points; hence no circle satisfies . Thus the circle-fitting problem has no solution.
The above sequence of circles converges (in our terms, W-converges) to the axis, which is a line, so it is natural to declare that line βthe best fit.ββ This may not be satisfactory in some practical applications, where one really needs to produce an estimate of the circleβs center and radius. But if we want to present the best-fitting object here, it is clearly and undeniably the line .
Thus, in order to guarantee the existence of the best-fitting object in all cases, we need to include in our collection all objects that can be obtained as limits of sequences of objects from . Such βlimit objectsββ are called limit points of , in the language of topology.
Definition 4.1. A collection which already contains all its βlimit pointsββ is said to be closed. If the collection is not closed, then an extended collection that includes all the limit points of is called the closure of .
For example, the collection of all lines in is closed, as a sequence of lines can only converge to a line. The collection of all circles in is not closed, as a sequence of circles may converge to a line. The closure of the collection of circles includes all circles and all lines, that is, (Strictly speaking, a sequence of circles may also converge to a single point, see Section 5, so singletons need to be included too.)
We see that the collection of model objects must be closed if we want the best-fitting object to exist in all cases. If is not closed, we have to extend it by adding all its limit points and thus make it closed.
Theorem 4.2 (existence of the best fit). Suppose that the given collection of model objects is closed. Then for any given points there exists the best-fitting object , that is, the objective function attains its global minimum on .
Proof. The key ingredients of our proof will be the continuity of the objective function and the compactness of a restricted domain of that function. We recall that a continuous real-valued function on a compact set always takes its maximum value and its minimum value on that set. In metric spaces (like our ), a subset is compact if every sequence of its elements has a subsequence that converges to another element , that is, as .
We recall that our objective function is continuous, and its domain is now assumed to be closed. If it was compact, the above general fact would guarantee that had a global minimum, as desired. But in most practical settings is not compact. Indeed, if it was compact, the function would not only have a minimum, but also a maximum. And this is usually impossible; one can usually find model objects arbitrarily far from the given points, thus making the distances from the points to the object arbitrarily large. Thus we need to find a smaller (restricted) subcollection which will be compact, and then we will apply the above general fact.
Let denote the disk of radius centered on the origin . We can find a large enough disk that satisfies two conditions: (i) covers all the given points and (ii) intersects at least one object , that is, . The distances from the given points to cannot exceed the diameter of , which is , hence .
Now we define our subcollection ; it consists of all model objects that intersect the larger disk of radius . Objects that lie entirely outside are not included in . Note that the subcollection contains at least one object ( mentioned above); hence is not empty.
Recall that all our given points lie in . They are separated from the region outside the larger disk by a βno manβs landβββthe ring , which is wide. Thus the distances from the given points to any object which was not included in are greater than ; hence for such objects we have . So objects not included in cannot fit our points better than does. Thus they can be ignored in the process of minimization of . More precisely, if we find the best-fitting object within the subcollection , then for any other object we will have
which shows that will be also the best-fitting object within the entire collection .
It remains to verify that the subcollection is compact; that is, every sequence of objects has a subsequence converging to another object . We will use the following standard fact. In a compact metric space any sequence of compact subsets has a convergent subsequence with respect to the Hausdorff metric.
Now let be the smallest integer greater than . Recall that all the objects in are required to intersect ; thus, they all intersect as well. The sets are compact. By the above general fact, there is a subsequence in the sequence and a compact subset such that as . Next, from the subsequence we extract a subsubsequence, call it that converges in the larger disk ; that is, such that as for some compact subset . Since , we see that ; that is, the limit sets and βagreeββ within .
Then we continue this procedure inductively for the progressively larger disks . In the end we use standard Cantorβs diagonal argument to construct a single subsequence such that for every we have as , and the limiting subsets βagreeββ in the sense for every . Then it follows that the sequence of objects converges (i.e., W converges) to the closed set . The limit set must belong to our collection because that collection was assumed to be closed. Lastly, intersects the disk , so it also belongs to the subcollection . Our proof is now complete.
5. Some Popular Models
Lines and Planes
The most basic (and oldest) fitting application is fitting lines in 2D. The model collection consists of all lines . It is easy to see that a sequence of lines can only converge to another line; hence, the collection of lines is closed. This guarantees that the fitting problem always has a solution.
Similarly, the collection of lines in 3D is closed. The collection of planes in 3D is closed, too. Thus the problems of fitting lines and planes always have a solution. The same is true for the problem of fitting -dimensional planes in a -dimensional space for any .
Circles and Spheres
A less-trivial task is fitting circles in 2D. Now the model collection consists of all circles . A sequence of circles can converge to an object of three possible types: a circle , a line , or a single point (singleton) . Singletons can be regarded as degenerate circles (whose radius is zero), or they can be ignored based on the redundancy principle introduced in Section 3. But the lines are objects of a different type; thus, the collection of circles is not closed.
As a result, the circle fitting problem does not always have a solution; a simple example was given in Section 4. In order to guarantee the existence of the best-fitting object, one needs to extend the given collection by adding all its βlimit points,ββ in this case-lines. The extended collection will include all circles and all lines, that is,
Now the (extended) circle fitting problem always has a solution; for any data points there is a best-fitting object that minimizes the sum of squares of the distances to . One has to keep in mind, though, that the best-fitting object may be a line, rather than a circle.
The problem of fitting circles to data points has been around since the 1950βs [13]. Occasional nonexistence of the best-fitting circle has been first noticed empirically; see, for example, Berman and Culpin [22]. The first theoretical analysis of the non-existence phenomenon was done by Nievergelt [21] who traced it to the noncompactness of the underlying model space and concluded that lines needed to be included in the model space to guarantee the existence of the best fit. He actually studied this phenomenon in many dimensions, so he dealt with hyperspheres and hyperplanes. He called hyperplanes βgeneralized hyperspheres.ββ Similar analysis and conclusions, for circles and lines, were published later, independently, by Chernov and Lesort [23] and by Zelniker and Clarkson [24]; see also a recent book [4].
Ellipses
A popular task in computer vision is fitting ellipses. Now the model collection consists of all ellipses . A sequence of ellipses may converge to an object of many possible types: an ellipse, a parabola, a line, a pair of parallel lines, a segment of a line, a half-line (a ray); see Figure 3.
Line segments can be regarded as degenerate ellipses whose minor axis is zero. Alternatively, line segments and half-lines can be ignored based on the redundancy principle (Section 3). Thus we have lines, pairs of parallel lines, and parabolas as limit objects of different types.
We see that the collection of ellipses is not closed. As a result, the ellipse-fitting problem does not always have a solution. For example, this happens when given points that are placed on a line or on a parabola. (Actually, there are many more examples; see Section 6.)
In order to guarantee the existence of the best-fitting object, one needs to extend the collection of ellipses by adding all its βlimit points,ββ in this caseβlines, pairs of parallel lines, and parabolas. We denote that extended collection by that where denotes the collection of pairs of parallel lines, and the collection of parabolas (the symbol resembles a parabola).
Now the (extended) ellipse-fitting problem always has a solution, for any data points there will be a best-fitting object that minimizes the sum of squares of the distances to . One has to keep in mind, though, that the best-fitting object may be a line, or a pair of parallel lines, or a parabola, rather than an ellipse.
The ellipse fitting problem has been around since the 1970s. The need to deal with limiting cases was first notice by Bookstein [16] who wrote βThe fitting of a parabola is a limiting case, exactly transitional between ellipse and hyperbola. As the center of ellipse moves off toward infinity while its major axis and the curvature of one end are held constantβ¦.ββ A first theoretical analysis on the non-existence of the best ellipse was done by Nievergelt in [25] who traced it to the noncompactness of the underlying model space and concluded that parabolas needed to be included in the model space to guarantee the existence of the best fit.
Conics
We now turn to a more general task of fitting quadratic curves (also known as conic sections, or simply conics). Now the model collection consists of all nondegenerate quadratic curves, by which we mean ellipses, parabolas, and hyperbolas. A sequence of conics may converge to an object of many possible types: a conic, a line, a line segment, a ray, two opposite rays, and a pair of parallel lines, and a pair of intersecting lines; see Figures 3 and 4.
Every line segment, ray, and two opposite rays are a part (subset) of a full line, so these objects can be ignored based on the redundancy principle (Section 3). But pairs of parallel lines and pairs of intersecting lines constitute limit objects of new types. We see that the collection of quadratic curves (conics) is not closed. As a result, the conic fitting problem does not always have a solution.
In order to guarantee the existence of the best-fitting object, one needs to extend the collection of quadratic curves by adding all its βlimit pointsβββpairs of lines (which may be intersecting, parallel, or coincident) where denotes the collection of pairs of lines.
Now the (extended) conic fitting problem always has a solution; for any data points there is a best-fitting object that minimizes the sum of squares of the distances to . One has to keep in mind, though, that the best-fitting object may be a pair of lines, rather than a genuine conic.
6. Sufficient and Deficient Models
We see that if the collection of adopted model objects is not closed, then in order to guarantee the existence of the best fit we have to close it up, that is, add to all objects that are obtained as limits of existing objects. The new, extended collection will be closed and will always provide the best fit. Here we discuss practical implications of extending .
Definition 6.1. The original objects are called primary objects and the objects that have been added, that is, are called secondary objects.
The above extension of by secondary objects, though necessary, may not be completely harmless. Suppose again that one needs to fit circles to observed points. Since the collection of circles is not closed, one has to add lines to it. That is, one really has to operate with the extended collection in which circles are primary objects and lines are secondary objects. The best-fitting object to a given set of points will then belong to ; hence occasionally it will be a line, rather than a circle.
In practical applications, though, a line may not be totally acceptable as the best-fitting object. For example, if one needs to produce an estimate of the circleβs radius and center, a line will be of little help; it has no radius or center. Thus, when a secondary object happens to be the best fit to the given data points, various unwanted complications may arise.
Next we investigate how frequently such unwanted events occur.
In most applications points are obtained (measured or observed) with some random noise. The noise has a probability distribution, which in statistical studies is usually assumed to be normal [3, 4]. For our purposes it is enough to just assume that given points have an absolutely continuous distribution:
Assumption 6.2. The points are obtained randomly with a probability distribution that has a joint probability density .
Let us first see how frequently the best-fitting circle fails to exist under this assumption. We begin with the simplest case of points. If they are not collinear, then they can be interpolated by a circle. If they are collinear (and distinct), there is no interpolating circle, so the best fit is achieved by a line. The three points , are collinear if and only if This equation defines a hypersurface in . Under the above assumption, the probability of that hypersurface is zero; that is, the collinearity occurs with probability zero. In simple practical terms, it is βimpossible,ββ it βnever happens.ββ If we generate three data points by using a computer random number generator, we will practically never see collinear points. All our practical experience tells us; the best-fitting circle always exists. Perhaps for this reason the circle-fitting problem is usually studied without including lines in the collection of model objects.
In the case of data points, the best-fitting circle fails to exist when the points are collinear, but there are other instances too. For example, let points be at , , , and for some large . Then it is not hard to check, by direct inspection, that the best-fitting circle fails to exist, and the best fit is achieved by the line .
It is not easy to describe all sets of points for which the best-fitting circle fails to exist, but such a description existsβit is given in [4] (see Theorem 8 on page 68 there). We skip technical details, only state the final conclusion; for every , the best-fitting circle exists unless the -vector belongs to a certain algebraic submanifold in . So, under our assumption, the best-fitting circle exists with probability one.
In other words, no matter how large or small the data set is, the best-fitting circle will exist with probability one, so practically we never have to resort to secondary objects (lines); that is, we never have to worry about existence. The model of circles is adequate and sufficient; it does not really require an extension.
Definition 6.3. We say that a collection of model objects is sufficient (for fitting purposes) if the best-fitting object exists with probability one under the above assumption. Otherwise the model collection will be called deficient.
We note that under our assumption an event occurs with probability zero if and only if it occurs on a subset of Lebesgue measure zero, that is, Leb. Therefore, our notion of sufficiency does not depend on the choice of the probability density .
As we have just seen, the collection of circles is sufficient for fitting purposes. Next we examine the model collection of ellipses . As we know, this collection is not closed; its closure must include all lines, pairs of parallel lines, and parabolas.
Ellipses
Let us see how frequently the best-fitting ellipse fails to exist. It is particularly easy to deal with data points. For any set of distinct 5 points in a general linear position (which means that no three points are collinear), there exists a unique quadratic curve (conic) that passes through all of them (i.e., interpolates them). That conic may be an ellipse, or a parabola, or a hyperbola. If 5 points are not in generic linear position (i.e., at least three of them are collinear), then they can be interpolated by a degenerate conic (a pair of lines).
If the interpolating conic is an ellipse, then obviously that ellipse is the best fit (the objective function takes it absolute minimum value, zero). What if the interpolating conic is a parabola or a pair of parallel lines? Then it is a secondary object in the extended model (5.2), and we have an unwanted event; a secondary object provides the best fit. We remark, however, that this event occurs with probability zero, so it is not a real concern. But there are other ways in which unwanted events occur, see below.
Suppose that the interpolating conic is a hyperbola or a pair of intersecting lines. Now the situation is less clear, as such an interpolating conic would not belong to the extended collection ; see (5.2). Is it possible then that the best-fitting object is an ellipse? The answer is βnoββ, and our argument is based on the following theorem that is interesting itself:
Theorem 6.4 (no local minima for ). For any set of five points the objective function has no local minima on conics. In other words, for any non interpolating conic there exist other conics , arbitrarily close to , which fit our data points better than does, that is, .
The proof is given in the appendix, and here we answer the previous question. Suppose points are interpolated by a hyperbola or a pair of intersecting lines. If there were a best-fitting ellipse , then no other ellipse could provide a better fit; that is, for any other ellipse we would have , which contradicts the above theorem. Thus the best-fitting object is not an ellipse but a secondary object (a parabola, or a line, or a pair of parallel lines), that is, an unwanted event occurs.
One can easily estimate the probability of the above unwanted event numerically. For a given probability distribution one can generate random samples of points and for each sample find the interpolating conic (by elementary geometry). Every time that conic happens to be other than an ellipse, an unwanted event occurs; that is, the best fit from within the collection (5.2) is provided by a secondary object.
We have run the above experiment with a standard normal distribution, where each coordinate and (for ) was generated independently according to a standard normal law . We found that the random points were interpolated by an ellipse with probability 22% and by a hyperbola with probability 78%. (Other conics, including parabolas, never turned up; in fact, it is not hard to prove that they occur with probability zero.) This is a striking result; hyperbolas actually dominate over ellipses! Thus the best-fitting ellipse fails to exist with a probability as high as 78%.
We note that this percentage will remain the same if we generate points independently according to any other 2D normal distribution. Indeed, any 2D normal distribution can be transformed to a standard 2D normal distribution by a linear map of the plane . Now under linear transformations conics are transformed to conics, and their types are preserved (i.e., ellipses are transformed to ellipses, hyperbolas to hyperbolas, etc.). Therefore the type of the interpolating conic cannot change.
In another experiment we sampled points from the unit square with a uniform distribution; that is, we selected each coordinate and independently from the unit interval . In this experiment ellipses turned up with probability 28% and hyperbolas with probability 72%, again an overwhelming domination of hyperbolas!
And again, these percentages will remain the same if we generate points according to any other rectangular (uniform) distribution. Indeed, any rectangle can be transformed to the square by a linear map of the plane, and then the same argument as before applies.
We also conducted numerical tests for points. In these tests we generated points , independently according to four predefined distributions: normal with covariance matrix , normal with covariance matrix , uniform in the square , uniform in the rectangle . Table 2 shows the percentage of randomly generated samples of points for which the best-fitting ellipse exists.
We see that in most cases the best-fitting ellipse fails to exist with a high probability. Only for the uniform distribution is the square the failures decrease as , and apparently the probability of existence of the best-fitting ellipse grows to 100%. We will explain this in Section 7.
We note that for points there is no interpolating conic, so one has to find a global minimum of via an extensive search over the model space of conics. This is a very time-consuming procedure and its results are never totally reliable. Thus our estimates, especially for larger βs, are quite approximate, but we believe that they are accurate to within 5%.
Insufficiency of Ellipses
We conclude that the collection of ellipses is not sufficient for fitting purposes. This means that there is a real chance that for a given set of data points no ellipse could be selected as the best fit to the points; that is, the ellipse-fitting problem would have no solution. More precisely, for any ellipse there will be another ellipse providing a better fit, in the sense . If one constructs a sequence of ellipses that fit the given points progressively better and on which the objective function converges to its infimum, then those ellipses will grow in size and converge to something different than an ellipse (most likely, to a parabola).
In practical terms, one usually runs a computer algorithm that executes an iterative procedure such as Gauss-Newton or Levenberg-Marquardt. It produces a sequence of ellipses (here denotes the iteration number) such that ; that is, the quality of approximations improves at every step, but those ellipses would keep growing in size and approach a parabola. This situation commonly occurs in practice [16, 25, 26], and one either accepts a parabola as the best fit or choose an ellipse by an alternative procedure (e.g., by sirect ellipse fit [27]).
We see that whenever the best-fitting ellipse fails to exist, the ellipse-fitting procedure will attempt to move beyond the collection of ellipses, so it will end up on the border of that collection and then return a secondary object (a parabola or a pair of lines). In a sense, the scope of the collection of ellipses is too narrow for fitting purposes. In other words, this collection is deficient, or badly incomplete. Its deficiency indicates that it should be substantially extended for the fitting problem to have a reasonable (adequate) solution with probability one.
Conics
Such an extension is achieved by the collection of all conics (ellipses, hyperbolas, and parabolas) denoted by ; see Section 5. If one searches for the best-fitting object in the entire collection of conics, then will exist with probability one. This was confirmed in numerous computer experiments that we have conducted. In fact, the best-fitting object has always been either an ellipse or hyperbola, so parabolas apparently can be excluded. Thus our numerical tests strongly suggest that the collection of ellipses, and hyperbolas is sufficient (though we do not have a proof).
7. Generalizations
In the basic formula (1.1) all the points make equal contributions to . In some applications, though, one needs to minimize a weighted sum for some fixed . This problem is known as weighted least squares. All our results on the existence of the best fit apply to it as well.
More generally, given a probability distribution on , one can βfitββ an object to it (or approximate it) by minimizing where . Note that (7.1) is a particular version of (7.2) in which is a discrete measure with a finite support. In order to guarantee the finiteness of the integral in (7.2) one needs to assume that the distribution has finite second moments. After that all our results apply.
The generalization (7.2) is relevant to the numerical tests described in Section 6. Suppose again that points are selected independently according to a probability distribution on . As grows, due to the law of large numbers of probability theory that the normalized objective function converges exactly to the integral (7.2).
Now the global minimum of the function would give us the best-fitting object corresponding to the given probability distribution , which would be the best βasymptoticββ fit to a random sample of size selected from the distribution , as . This relation helps us explain some experimental results reported in Section 6.
Suppose that is a uniform distribution in a rectangle . Then the integral in (7.2) becomes We have computed it numerically in order to find the best-fitting conic . The results are given below.
Perfect Square
One would naturally expect that the best fit to a mass uniformly distributed in a perfect square would be either a circle or a pair of lines (say, the diagonals of the square).
What we found was totally unexpected; the best-fitting conic is an ellipse! Its center coincides with the center of the square, and its axes are aligned with the sides of the square, but its axes are different; the major axis is and the minor axis is (assuming that the square has unit side). In fact, there are exactly two such ellipses; one is oriented horizontally and the otherβvertically; see the left panel of Figure 5. These two ellipses beat any circle and any pair of lines (see the table in Figure 5); they provide the global minimum of the objective function (7.3).
The fact that the best-fitting conic to a mass distributed in a square is an ellipse explains why for large samples generated from the uniform distribution in a square ellipses dominate over hyperbolas (Section 6).
As the rectangle gets extended horizontally, that is, as increases the best-fitting conic changes. We found that for , the best-fitting conic is still an ellipse (though it gets elongated compared to the one found for the square, when ). But for , the best-fitting conic abruptly changes from an ellipse to a pair of parallel lines. Those two lines are running through the rectangle horizontally; their equations are and (See Figure 6).
This fact might explain why for longer rectangles, such as used in Section 6, the percentage of samples for which the best-fitting ellipse exists does not grow much as increases. Though it remains unclear which type of conics (ellipses or hyperbolas) dominates for large samples taken from longer rectangles with .
8. Megaspace
Our conclusions can be illustrated by an interesting construction in a multidimensional space βmegaspaceββ. It was first used by Malinvaud [28, Chapter 10] and then by Chernov [4, Sections 1.5 and 3.4].
Recall that given data points , and a model object (a closed subset of ), our objective function is defined by (1.1) and for the distances we can write Thus we can express the objective function as follows:
Now let us represent the data points by one point (βmegapointββ) in the -dimensional space with coordinates . For the given object , let us also define a multidimensional set βmegasetββ as follows: Note that in (8.2) is the square of the distance from to , in the megaspace . Therefore consider Next given a collection of model objects we define a large megaset as follows: . Alternatively, it can be defined as The best-fitting object minimizes the function . Thus, due to (8.4), minimizes the distance from the megapoint representing the data set to the megaset representing the collection .
Thus, the problem of finding the best-fitting object reduces to the problem of finding the megapoint that is closest to the given megapoint (representing the given points). In other words, we need to project the megapoint onto the megaset , and the footpoint of the projection would correspond to the best-fitting object .
We conclude that the best-fitting object exists if and only if there exists a megapoint that is closest to the given megapoint . It is a simple fact that given a set with , the closest point to a given point exists for any if and only if the set is closed. Thus the existence of the best-fitting object requires the megaset to be topologically closed. Again we see that the property of closedness is vital for the fitting problem to have a solution for every set of data points .
Theorem 8.1 (closedness of megasets). If the model collection is closed (in the sense of W convergence), then the megaset is closed in the natural topology of .
Clearly, this theorem provides an alternative proof of the existence of the best-fitting object, provided that the model collection is closed.
Proof. Suppose a sequence of megapoints
all belonging to the megaset , converges, as , to a megapoint
in the usual topology of . Equivalently, for every the point converges, as , to the point . We need to show that .
Now for every there exists a model object that contains the points . The sequence contains a convergent subsequence , that is, such that as for some closed set . (The existence of a convergent subsequence can be verified by the diagonal argument as in the proof of the existence theorem in Section 4.)
As we assumed that the collection is closed, it follows that contains the limit object , that is, . It is intuitively clear (and can be verified by a routine calculus-type argument) that contains all the limit points . Therefore the limit megapoint belongs to the megaset . This proves that the latter is closed.
Megaset for Lines
Let consist of all lines in . The corresponding megaset is described in [28, Chapter 10]. A megapoint belongs to if and only if all the planar points belong to one line (i.e., they are collinear). This condition can be expressed by algebraic relations as follow:
for all . Each of these relations means that the three points , , and are collinear; cf. (6.1). All of these relations together mean that all the points are collinear.
Note that is specified by -independent relations; hence it is an -dimensional manifold (algebraic variety) in . The relations (8.8) are quadratic, so is a quadratic surface. It is closed in topological sense; hence the problem of finding the best-fitting line always has a solution.
Megaset for Circles
Let consist of all circles in . The corresponding megaset is described in [4, Section 3.4], and we briefly repeat it here. A megapoint belongs to if and only if all the planar points belong to one circle (we say that they are βcocircularββ). In that case all these points satisfy one quadratic equation of a special type as follow:
This condition can be expressed by algebraic relations as follow:
for . Note that (8.10) include megapoints satisfying (8.8); hence they actually describe the union .
The determinant in (8.10) is a polynomial of degree four, and is an -dimensional algebraic variety (manifold) in defined by quartic polynomial equations. Note that the dimension of is one higher than that of , that is,
A closer examination shows that plays the role of the boundary of , that is, terminates on . The megaset is not closed, but if we add its boundary to it, it will become closed.
Megaset for Ellipses
Let consist of all ellipses in . The corresponding megaset can be described in a similar manner as above. A point belongs in if and only if all the planar points belong to one ellipse (we say that they are βcoellipticalββ). In that case all these points satisfy one quadratic equation of a general type as follow:
This equation actually means that the points belong to one conic (either regular or degenerate). This condition can be expressed by algebraic relations as follow:
for . Each of these relations means that the six points , , , , , and satisfy one quadratic equation of type (8.12); that is, they belong to one conic (either regular or degenerate). All of these relations together mean that all the points satisfy one quadratic equation of type (8.12); that is, they all belong to one conic (regular or degenerate). Therefore the relations (8.13) describe a much larger megaset corresponding to the collection of all conics, regular and degenerate, that is,
using our previous notation.
The determinant in (8.13) is a polynomial of the eighth degree, and is a closed -dimensional algebraic manifold in . It is mostly made of two big megasets: and and they both are -dimensional. Other megasets (parabolas and pairs of lines) have smaller dimensions and play the role of the boundaries of the bigger megasets and .
Illustrations
The structure of the megaset is schematically depicted in Figure 7, where it is shown as the plane in the 3D space (the latter plays the role of the megaspace ). The positive half-plane represents the elliptic megaset , and the negative half-plane represents the hyperbolic megaset . The -axis separating these two half-planes represents the lower-dimensional megasets in the decomposition (8.14). True, the real structure of is much more complicated, but our simplified picture still shows its most basic features.
Now recall that finding the best-fitting conic corresponds to an orthogonal projection of the given megapoint in the megaspace (in our illustration, it would be a point ) onto the megaset (in our illustrationβonto the plane). In Figure 7 the point is simply projected onto . What are the chances that the footpoint of the projection corresponds to a βboundaryββ object in (i.e., to a secondary object)? Clearly, only the points of the plane are projected onto the line . If the point is selected randomly with an absolutely continuous distribution, then a point on the plane would be chosen with probability zero. This fact illustrates the sufficiency of the model collection of nondegenerate conics (even the sufficiency of ellipses and hyperbolas alone, without parabolas); the best-fitting object will be a primary object with probability one.
But what if our model collection consists of ellipses only, without hyperbolas? Then in our illustration, the corresponding megaset would be the positive half-plane . Finding the best-fitting ellipse would correspond to an orthogonal projection of the given point onto the positive half-plane . Now if the given point has a positive -coordinate, then it is projected onto , as before, and we get the desired best-fitting ellipse. But if it has a negative -coordinate, then it is projected onto , which is on the boundary of the half-plane, so we get a boundary footpoint; that is, a secondary object will be the best fit. We see that all the points with (making a whole half-space!) are projected onto the boundary line, hence for all those points the best-fitting ellipse would not exist! This fact clearly illustrates the deficiency of the model collection of ellipses.
One may wonder: How is it possible that the collection of circles is sufficient (as we proved in Section 5), while the larger collection of ellipses is not? Indeed every circle is an ellipse, hence the collection of ellipses contains all the circles. So why the sufficiency of circles does not guarantee the sufficiency of the bigger, inclusive collection of ellipses? Well, this seemingly counterintuitive fact can be illustrated too (see Figure 8).
Suppose that the megaset for the collection of circles is represented by the set in our illustration. Note that consists of two curves (branches of a parabola on the plane), both lie in the half-plane that corresponds to the collection of ellipses. So the required inclusion does take place. The two curves making terminate at the point , which does not belong to , so it plays the role of the boundary of .
Now suppose a randomly selected point is to be projected onto the set . What are the chances that its projection will end up on the boundary of , that is, at the origin ? It is not hard to see (and prove by elementary geometry) that only points on the plane may be projected onto the origin (and not even all of them; points with large positive -coordinates would be projected onto some interior points of ). So the chance that the footpoint of the projection ends up at the boundary of is zero. This illustrates the sufficiency of the smaller model collection of circles, despite the deficiency of the larger model collection of ellipses (which we have seen in Figure 7).
Appendix
Here we prove the theorem on no local minima for (Section 6).
Proof. Suppose that is a conic on which takes a local minimum. We will assume that is nondegenerate; that is, is an ellipse, or hyperbola, or parabola. Degenerate conics are easier to handle, so we omit them.
Let denote the given points and , , , , and their projections onto . Since does not interpolate all the given points simultaneously, there exists at least one which is not on . Suppose that it is , hence .
Let us first assume that is different from the other projections . Then let be a point near such that dist. In other words, we perturb slightly by moving it in the direction of . Then there exists another non-degenerate conic that interpolates the points . It is easy to see that
As a result,
as desired. Since can be placed arbitrarily close to , the new conic can be made arbitrarily close to .
Now suppose that the projection point coincides with some other projection point(s). In that case there are at most four distinct projection points. We will not move the point , but we will rotate (slightly) the tangent line to the conic at the point . Denote that tangent line by . Note that is perpendicular to the line passing through and .
Now let be another line passing through the point and making an arbitrarily small angle with . It is not hard to show, by elementary geometry, that there exists another nondegenerate conic passing though the same (β€4) projection points and having the tangent line at .
Since the tangent line to the conic at the point is not orthogonal to the line , we easily have
At the same time we have
As a result, consider
as desired. Since the line can be selected arbitrarily close to , the new conic will be arbitrarily close to .
Acknowledgment
N. Chernov was partially supported by National Science Foundation, Grant DMS-0969187.