Does the Best-Fitting Curve Always Exist?

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 ℱ(𝑆)=^𝑛_𝑖=1dist𝑃_𝑖,𝑆²,(1.1) 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).

Table 1 

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 𝑛=4 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 𝑑≥2, but for the sake of notational simplicity and ease of illustrations we mostly restrict the exposition to the planar case 𝑑=2.

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 (0,0) and their slopes 1/𝑛 decrease as 𝑛 grows (Figure 1). Naturally, we would consider this sequence as convergent; it converges to the 𝑥-axis, that is, to the line 𝐿={𝑦=0}, as 𝑛→∞.

Figure 1 

Sequence of lines 𝐿𝑛 and the limit line 𝐿.

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 𝑆₁and𝑆₂⊂ℝ² are. It is defined by dist_H𝑆₁,𝑆₂=maxsup_{𝑃1∈𝑆1}dist𝑃₁,𝑆₂,sup_{𝑃2∈𝑆2}dist𝑃₂,𝑆₁.(2.1) If the two sets 𝑆₁and𝑆₂ are closed and the Hausdorff distance between them is zero, that is, dist_H(𝑆₁,𝑆₂)=0, then they coincide 𝑆₁=𝑆₂. If the Hausdorff distance is small, the two sets 𝑆₁and𝑆₂ 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, dist_H(𝐿_𝑛,𝐿)=∞ 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 𝑅={−𝐴≤𝑥≤𝐴,−𝐵≤𝑦≤𝐵},(2.2) 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: dist_H(𝐿_𝑛∩𝑅,𝐿∩𝑅)→0 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: dist_H𝑆₁,𝑆₂;𝑅=maxsup_{𝑃∈𝑆1∩𝑅}dist𝑃,𝑆₂,sup_{𝑄∈𝑆2∩𝑅}dist𝑄,𝑆₁.(2.3) The formula (2.3) applies whenever both sets, 𝑆₁ and 𝑆₂, intersect the window 𝑅. If only one set, say 𝑆₁, intersects 𝑅, we modify (2.3) as follows: dist_H𝑆₁,𝑆₂;𝑅=sup_{𝑃∈𝑆1∩𝑅}dist𝑃,𝑆₂.(2.4) A similar modification is used if only 𝑆₂ intersects 𝑅. If neither set intersects the window 𝑅, we simply set dist_H(𝑆₁,𝑆₂;𝑅)=0 (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 dist_H𝑆_𝑛,𝑆;𝑅→0as𝑛→∞.(2.5)

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 𝑆₁and𝑆₂⊂ℝ² is defined as follows: dist_W𝑆₁,𝑆₂=^∞_𝑘=12^−𝑘dist_H𝑆₁,𝑆₂;𝑅_𝑘,(2.6) where 𝑅_𝑘 is a square window of size 2𝑘×2𝑘, 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 2^−𝑘. 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 2√2𝑘) the above sum is bounded by 2√2∑^∞𝑘=1𝑘2^−𝑘<6.

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 dist(𝑃,𝑆)≤dist(𝑃,𝑆^); 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 dist(𝑃_𝑖,𝑆) to individual points 𝑃_𝑖, see (1.1), it is enough to verify that the distance dist(𝑃,𝑆) 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 dist(𝑃,𝑄)=dist(𝑃,𝑆); such a point 𝑄 exists because 𝑆 is closed. Denote by 𝐷 the disk centered on 𝑃 of radius 1+dist(𝑃,𝑄); it contains 𝑄. Let 𝑅 be a window containing the disk 𝐷.
Since 𝑅 contains 𝑄, it intersects with 𝑆, that is, 𝑅∩𝑆≠∅. This guarantees that dist_H(𝑆_𝑚,𝑆;𝑅)→0. Thus, there are points 𝑄_𝑚∈𝑆_𝑚 such that 𝑄_𝑚→𝑄. Since dist(𝑃,𝑆_𝑚)≤dist(𝑃,𝑄_𝑚), we conclude that the upper limit of the sequence dist(𝑃,𝑆_𝑚) does not exceed dist(𝑃,𝑆), that is, limsupdist𝑃,𝑆_𝑚≤dist(𝑃,𝑆).(3.1) On the other hand, we will show that the lower limit of the sequence dist(𝑃,𝑆_𝑚) cannot be smaller than dist(𝑃,𝑆), that is, liminfdist𝑃,𝑆_𝑚≥dist(𝑃,𝑆).(3.2) The estimates (3.1) and (3.2) together imply that dist(𝑃,𝑆_𝑚)→dist(𝑃,𝑆), as desired, and hence the distance function dist(𝑃,𝑆) 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 lim_𝑘→∞dist𝑃,𝑆_𝑚𝑘=liminfdist𝑃,𝑆_𝑚<dist(𝑃,𝑆).(3.3) Denote by 𝑄_𝑚∈𝑆_𝑚 the point in 𝑆_𝑚 closest to 𝑃, that is, such that dist(𝑃,𝑄_𝑚)=dist(𝑃,𝑆_𝑚). Then we have lim_𝑘→∞dist𝑃,𝑄_𝑚𝑘=lim_𝑘→∞dist𝑃,𝑆_𝑚𝑘<dist(𝑃,𝑆)=dist(𝑃,𝑄).(3.4) Since the points 𝑄_𝑚𝑘 are closer to 𝑃 than the point 𝑄 is, we have 𝑄_𝑚𝑘∈𝐷⊂𝑅. Recall that dist_H(𝑆_𝑚,𝑆;𝑅)→0, hence dist𝑄_𝑚𝑘,𝑆→0as𝑘→∞.(3.5) Denote by 𝐻_𝑚𝑘∈𝑆 the point in 𝑆 closest to 𝑄_𝑚𝑘, that is, such that dist(𝑄_𝑚𝑘,𝐻_𝑚𝑘)=dist(𝑄_𝑚𝑘,𝑆). Now we have, by triangle inequality, dist(𝑃,𝑆)≤dist𝑃,𝐻_𝑚𝑘≤=ddiisstt𝑃𝑃,,𝑄𝑄_{𝑚𝑘𝑚𝑘}+dist𝑄_𝑚𝑘,𝐻_𝑚𝑘+dist𝑄_𝑚𝑘,𝑆.(3.6) The limit of the first term on the right hand side of (3.6) is <dist(𝑃,𝑆) by (3.4), and the limit of the second term is zero by (3.5). This implies dist(𝑃,𝑆)<dist(𝑃,𝑆), 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 𝑆_best∈𝕄 corresponds to the (global) minimum of the objective function ℱ, that is, 𝑆_best=argmin_𝑆∈𝕄ℱ(𝑆). The function ℱ defined by (1.1) cannot be negative; thus, it always has an infimum ℱ₀=inf_𝑆∈𝕄ℱ(𝑆).(4.1) 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 ℱ(𝑆_𝑚)<ℱ(𝑆_𝑚−1); 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 𝑆₀=𝑆_best.

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 𝑃₁=(−1,0), 𝑃₂=(0,0) and 𝑃₃=(1,0). Then the sequence of circles 𝑆_𝑚 defined by 𝑥²+(𝑦−𝑚)²=𝑚² will fit the points progressively better (tighter) as 𝑚 grows, so that ℱ(𝑆_𝑚)→0 as 𝑚→∞ (Figure 2). On the other hand, no circle can pass through three collinear points; hence no circle 𝑆 satisfies ℱ(𝑆)=0. Thus the circle-fitting problem has no solution.

Figure 2 

Sequence of circles and the limit line.

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 𝑦=0.

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, 𝕄_𝐶=𝕄_𝐶∪𝕄_𝐿.(4.2) (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 𝑆_best∈𝕄, 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 (0,0). 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 2𝑟, hence ℱ(𝑆₀)≤(2𝑟)²𝑛.
Now we define our subcollection 𝕄₀⊂𝕄; it consists of all model objects 𝑆∈𝕄 that intersect the larger disk 𝐷_3𝑟 of radius 3𝑟. Objects that lie entirely outside 𝐷_3𝑟 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 𝐷_3𝑟 by a “no man’s land’’—the ring 𝐷_3𝑟⧵𝐷_𝑟, which is 2𝑟 wide. Thus the distances from the given points to any object 𝑆 which was not included in 𝕄₀ are greater than 2𝑟; hence for such objects we have ℱ(𝑆)>(2𝑟)²𝑛. 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 𝑆_best within the subcollection 𝕄₀, then for any other object 𝑆∈𝕄⧵𝕄₀ we will have ℱ(𝑆)>(2𝑟)²𝑛≥ℱ𝑆₀≥ℱ𝑆_best,(4.3) which shows that 𝑆_best 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 𝑗=[3𝑟]+1 be the smallest integer greater than 3𝑟. Recall that all the objects in 𝕄₀ are required to intersect 𝐷_3𝑟; 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 dist_H(𝑆^𝑘(𝑗)∩𝐷_𝑗,𝑆^∗𝑗)→0 as 𝑘→∞. Next, from the subsequence 𝑆^𝑘(𝑗) we extract a subsubsequence, call it 𝑆^𝑘(𝑗+1) that converges in the larger disk 𝐷_𝑗+1; that is, such that dist_H(𝑆^𝑘(𝑗+1)∩𝐷_𝑗+1,𝑆^∗𝑗+1)→0 as 𝑘→∞ for some compact subset 𝑆^∗𝑗+1⊂𝐷_𝑗+1. Since dist_H(𝑆^𝑘(𝑗+1)∩𝐷_𝑗,𝑆^∗𝑗)→0, we see that 𝑆^∗𝑗+1∩𝐷_𝑗=𝑆^∗𝑗; that is, the limit sets 𝑆^∗𝑗 and 𝑆^∗𝑗+1 “agree’’ within 𝐷_𝑗.
Then we continue this procedure inductively for the progressively larger disks 𝐷_𝑗+2,𝐷_𝑗+3,…. In the end we use standard Cantor’s diagonal argument to construct a single subsequence 𝑆_𝑚𝑘 such that for every 𝑖≥𝑗 we have dist_H(𝑆_𝑚𝑘∩𝐷_𝑖,𝑆^∗𝑖)→0 as 𝑘→∞, and the limiting subsets 𝑆^∗𝑖⊂𝐷_𝑖 “agree’’ in the sense 𝑆^∗𝑖+1∩𝐷_𝑖=𝑆^∗𝑖 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 𝐷_3𝑟, 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 1≤𝑘<𝑑.

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, 𝕄_𝐶=𝕄_𝐶∪𝕄_𝐿.(5.1) Now the (extended) circle fitting problem always has a solution; for any data points 𝑃₁,…,𝑃_𝑛 there is a best-fitting object 𝑆_best∈𝕄_𝐶 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.

Figure 3 

Sequences of ellipses converging to various limit objects.

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 𝕄_𝐸=𝕄_𝐸∪𝕄_𝐿∪𝕄_‖∪𝕄_∪,(5.2) 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 𝑆_best∈𝕄_𝐸 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.

Figure 4 

Sequences of hyperbolas converging to various limit objects.

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) 𝕄_𝑄=𝕄_𝑄∪𝕄_𝐿𝐿,(5.3) 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 𝑆_best∈𝕄_𝑄 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 𝜌(𝑥₁,𝑦₁,…,𝑥_𝑛,𝑦_𝑛)>0.

Let us first see how frequently the best-fitting circle fails to exist under this assumption. We begin with the simplest case of 𝑛=3 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 (𝑥_𝑖,𝑦_𝑖), 1≤𝑖≤3 are collinear if and only if 𝑥₁−𝑥₂𝑦₁−𝑦₃=𝑥₁−𝑥₃𝑦₁−𝑦₂.(6.1) 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 𝑛>3 data points, the best-fitting circle fails to exist when the points are collinear, but there are other instances too. For example, let 𝑛=4 points be at (0,1), (0,−1), (𝐴,0), and (−𝐵,0) for some large 𝐴,𝐵≫1. 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 𝑦=0.

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 𝑛>3, the best-fitting circle exists unless the (2𝑛)-vector (𝑥₁,𝑦₁,…,𝑥_𝑛,𝑦_𝑛) belongs to a certain algebraic submanifold in ℝ^2𝑛. 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 𝑆_best 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 𝐴⊂ℝ^2𝑛 of Lebesgue measure zero, that is, Leb(𝐴)=0. 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 𝑛=5 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 𝑆_best∈𝕄_𝐸 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 𝑛=5). 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 𝑛=5 points are interpolated by a hyperbola or a pair of intersecting lines. If there were a best-fitting ellipse 𝐸_best∈𝕄_𝐸, then no other ellipse could provide a better fit; that is, for any other ellipse 𝐸∈𝕄_𝐸 we would have ℱ(𝐸)≥ℱ(𝐸_best), which contradicts the above theorem. Thus the best-fitting object 𝑆_best∈𝕄_𝐸 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 𝑛=5 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 𝑖=1,…,5) was generated independently according to a standard normal law 𝒩(0,1). 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 [0,1]×[0,1] with a uniform distribution; that is, we selected each coordinate 𝑥_𝑖 and 𝑦_𝑖 independently from the unit interval [0,1]. 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 [0,1]×[0,1] by a linear map of the plane, and then the same argument as before applies.

We also conducted numerical tests for 𝑛>5 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 [0,1]×[0,1], uniform in the rectangle [0,2]×[0,1]. Table 2 shows the percentage of randomly generated samples of 𝑛 points for which the best-fitting ellipse exists.

Table 2 

Percentage (%) of randomly generated samples of n 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 [0,1]×[0,1] 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 𝑛≥6 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 ℱ(𝐸_𝑚)<ℱ(𝐸_𝑚−1); 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 𝑆_best in the entire collection of conics, then 𝑆_best 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 ℱ_𝑤(𝑆)=^𝑛_𝑖=1𝑤_𝑖dist𝑃_𝑖,𝑆²,(7.1) for some fixed 𝑤_𝑖>0. 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 dP(𝑥,𝑦) on ℝ², one can “fit’’ an object 𝑆 to it (or approximate it) by minimizing ℱ_𝑃(𝑆)=𝑑²(𝑥,𝑦)dP(𝑥,𝑦),(7.2) where 𝑑(𝑥,𝑦)=dist[(𝑥,𝑦),𝑆]. Note that (7.1) is a particular version of (7.2) in which dP 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 dP 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 dP(𝑥,𝑦) on ℝ². As 𝑛 grows, due to the law of large numbers of probability theory that the normalized objective function (1/𝑛)ℱ 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 dP, which would be the best “asymptotic’’ fit to a random sample of size 𝑛 selected from the distribution dP, as 𝑛→∞. This relation helps us explain some experimental results reported in Section 6.

Suppose that dP is a uniform distribution in a rectangle 𝑅=[0,𝐿]×[0,1]. Then the integral in (7.2) becomes ℱ_𝑃(𝑆)=𝐿⁻¹^𝐿0¹⁰𝑑²(𝑥,𝑦)𝑑𝑦𝑑𝑥.(7.3) We have computed it numerically in order to find the best-fitting conic 𝑆. The results are given below.

Perfect Square (𝑅=[0,1]×[0,1])
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 2𝑎=0.9489 and the minor axis is 2𝑏=0.6445 (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).

Figure 5 

The best fit to a uniform distribution in a square.

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 𝐿>1 increases the best-fitting conic changes. We found that for 𝐿<𝐿₁≈1.2, the best-fitting conic is still an ellipse (though it gets elongated compared to the one found for the square, when 𝐿=1). 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 𝑦=0.25 and 𝑦=0.75 (See Figure 6).

Figure 6 

The best fit to a uniform distribution in a rectangle.

This fact might explain why for longer rectangles, such as 𝑅=[0,2]×[0,1] 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 [0,𝐿]×[0,1] 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 dist(𝑃_𝑖,𝑆) we can write dist𝑃_𝑖,𝑆²=min_{(𝑥′𝑖,𝑦′𝑖)∈𝑆}𝑥^𝑖−𝑥^𝑖²+𝑦_𝑖−𝑦^𝑖².(8.1) Thus we can express the objective function ℱ(𝑆) as follows: ℱ(𝑆)=min^𝑛_𝑖=1𝑥^𝑖−𝑥^𝑖²+𝑦_𝑖−𝑦^𝑖²;𝑥^𝑖,𝑦^𝑖∈𝑆∀𝑖.(8.2)

Now let us represent the 𝑛 data points 𝑃₁=(𝑥₁,𝑦₁),…,𝑃_𝑛=(𝑥_𝑛,𝑦_𝑛) by one point (“megapoint’’) 𝒫 in the 2𝑛-dimensional space ℝ^2𝑛 with coordinates 𝑥₁,𝑦₁,…,𝑥_𝑛,𝑦_𝑛. For the given object 𝑆, let us also define a multidimensional set “megaset’’ 𝔐_𝑆⊂ℝ^2𝑛 as follows: 𝒫′=𝑥^1,𝑦^1,…,𝑥^𝑛,𝑦^𝑛∈𝔐_𝑆⟺𝑥^𝑖,𝑦^𝑖∈𝑆∀𝑖.(8.3) Note that ∑^𝑛𝑖=1[(𝑥_𝑖−𝑥^𝑖)²+(𝑦_𝑖−𝑦^𝑖)²] in (8.2) is the square of the distance from 𝒫 to 𝒫′∈𝔐_𝑆, in the megaspace ℝ^2𝑛. Therefore consider ℱ(𝑆)=min_{𝒫∈𝔐𝑆}dist𝒫,𝒫^²=dist𝒫,𝔐_𝑆².(8.4) Next given a collection 𝕄 of model objects we define a large megaset 𝔐(𝕄)⊂ℝ^2𝑛 as follows: 𝔐(𝕄)=∪_𝑆∈𝕄𝔐_𝑆. Alternatively, it can be defined as 𝑥^1,𝑦^1,…,𝑥^𝑛,𝑦^𝑛∈𝔐(𝕄)⟺∃𝑆∈𝕄∶𝑥^𝑖,𝑦^𝑖∈𝑆∀𝑖.(8.5) The best-fitting object 𝑆_best minimizes the function ℱ(𝑆). Thus, due to (8.4), 𝑆_best 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 𝑆_best 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 𝑆_best.

We conclude that the best-fitting object 𝑆_best∈𝕄 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 𝑑≥1, 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 𝑆_best 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 ℝ^2𝑛.

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 𝒫^(𝑘)=𝑥^1(𝑘),𝑦^1(𝑘),…,𝑥^𝑛(𝑘),𝑦^𝑛(𝑘),(8.6) all belonging to the megaset 𝔐(𝕄), converges, as 𝑘→∞, to a megapoint 𝒫^(∞)=𝑥^1(∞),𝑦^1(∞),…,𝑥^𝑛(∞),𝑦^𝑛(∞)(8.7) in the usual topology of ℝ^2𝑛. Equivalently, for every 𝑖=1,…,𝑛 the point (𝑥^𝑖(𝑘),𝑦^𝑖(𝑘)) converges, as 𝑘→∞, to the point (𝑥^𝑖(∞),𝑦^𝑖(∞)). We need to show that 𝒫^(∞)∈𝔐(𝕄).
Now for every 𝑘 there exists a model object 𝑆^(𝑘)∈𝕄 that contains the 𝑛 points (𝑥^1(𝑘),𝑦^1(𝑘)),…,(𝑥^𝑛(𝑘),𝑦^𝑛(𝑘)). 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 (𝑥^1(∞),𝑦^1(∞)),…,(𝑥^𝑛(∞),𝑦^𝑛(∞)). 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 𝔐(𝕄_𝐿)⊂ℝ^2𝑛 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 𝐶_𝑛,3 algebraic relations as follow: det⎡⎢⎢⎣𝑥^𝑖−𝑥_𝑗𝑦_𝑖−𝑦_𝑗𝑥_𝑖−𝑥_𝑘𝑦_𝑖−𝑦_𝑘⎤⎥⎥⎦=0(8.8) for all 1≤𝑖<𝑗<𝑘≤𝑛. 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 𝑛−2-independent relations; hence it is an (𝑛+2)-dimensional manifold (algebraic variety) in ℝ^2𝑛. 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 𝔐(𝕄_𝐶)⊂ℝ^2𝑛 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: 𝐴𝑥²+𝑦²+𝐵𝑥+𝐶𝑦+𝐷=0.(8.9) This condition can be expressed by 𝐶_𝑛,4 algebraic relations as follow: det⎡⎢⎢⎢⎢⎣𝑥^𝑖−𝑥_𝑗𝑦_𝑖−𝑦_𝑗𝑥^2𝑖−𝑥^2𝑗+𝑦^2𝑖−𝑦^2𝑗𝑥_𝑖−𝑥_𝑘𝑦_𝑖−𝑦_𝑘𝑥^2𝑖−𝑥^2𝑘+𝑦^2𝑖−𝑦^2𝑘𝑥_𝑖−𝑥_𝑚𝑦_𝑖−𝑦_𝑚𝑥^2𝑖−𝑥^2𝑚+𝑦^2𝑖−𝑦^2𝑚⎤⎥⎥⎥⎥⎦=0(8.10) for 1≤𝑖<𝑗<𝑘<𝑚≤𝑛. 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 (𝑛+3)-dimensional algebraic variety (manifold) in ℝ^2𝑛 defined by quartic polynomial equations. Note that the dimension of 𝔐(𝕄_𝐶)∪𝔐(𝕄_𝐿) is one higher than that of 𝔐(𝕄_𝐿), that is, dim𝔐𝕄_𝐶∪𝔐𝕄_𝐿=dim𝔐𝕄_𝐿+1.(8.11) 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 𝔐(𝕄_𝐸)⊂ℝ^2𝑛 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: 𝐴𝑥²+2𝐵𝑥𝑦+𝐶𝑦²+2𝐷𝑥+2𝐸𝑦+𝐹=0.(8.12) This equation actually means that the points belong to one conic (either regular or degenerate). This condition can be expressed by 𝐶_𝑛,6 algebraic relations as follow: det⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑥^𝑖−𝑥_𝑗𝑦_𝑖−𝑦_𝑗𝑥^2𝑖−𝑥^2𝑗𝑦^2𝑖−𝑦^2𝑗𝑥_𝑖𝑦_𝑖−𝑥_𝑗𝑦_𝑗𝑥_𝑖−𝑥_𝑘𝑦_𝑖−𝑦_𝑘𝑥^2𝑖−𝑥^2𝑘𝑦^2𝑖−𝑦^2𝑘𝑥_𝑖𝑦_𝑖−𝑥_𝑘𝑦_𝑘𝑥_𝑖−𝑥_𝑚𝑦_𝑖−𝑦_𝑚𝑥^2𝑖−𝑥^2𝑚𝑦^2𝑖−𝑦^2𝑚𝑥_𝑖𝑦_𝑖−𝑥_𝑚𝑦_𝑚𝑥_𝑖−𝑥_𝑙𝑦_𝑖−𝑦_𝑙𝑥^2𝑖−𝑥^2𝑙𝑦^2𝑖−𝑦^2𝑙𝑥_𝑖𝑦_𝑖−𝑥_𝑙𝑦_𝑙𝑥_𝑖−𝑥_𝑟𝑦_𝑖−𝑦_𝑟𝑥^2𝑖−𝑥^2𝑟𝑦^2𝑖−𝑦^2𝑟𝑥_𝑖𝑦_𝑖−𝑥_𝑟𝑦_𝑟⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦=0(8.13) for 1≤𝑖<𝑗<𝑘<𝑚<𝑙<𝑟≤𝑛. 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 𝔐(𝕄_Conics) corresponding to the collection of all conics, regular and degenerate, that is, 𝕄_Conics=𝕄_𝐸∪𝕄_𝐻∪𝕄_∪∪𝕄_𝐿𝐿(8.14) using our previous notation.
The determinant in (8.13) is a polynomial of the eighth degree, and 𝔐(𝕄_Conics) is a closed (𝑛+5)-dimensional algebraic manifold in ℝ^2𝑛. It is mostly made of two big megasets: 𝔐(𝕄_𝐸) and 𝔐(𝕄_𝐻) and they both are (𝑛+5)-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 𝔐(𝕄_Conics) is schematically depicted in Figure 7, where it is shown as the 𝑥𝑦 plane {𝑧=0} in the 3D space (the latter plays the role of the megaspace ℝ^2𝑛). The positive half-plane 𝐻₊={𝑦>0,𝑧=0} represents the elliptic megaset 𝔐(𝕄_𝐸), and the negative half-plane 𝐻₋={𝑦<0,𝑧=0} represents the hyperbolic megaset 𝔐(𝕄_𝐻). The 𝑥-axis {𝑦=𝑧=0} separating these two half-planes represents the lower-dimensional megasets 𝔐(𝕄_∪∪𝕄_𝐿𝐿) in the decomposition (8.14). True, the real structure of 𝔐(𝕄_Conics) is much more complicated, but our simplified picture still shows its most basic features.

Figure 7 

Schematic illustration of the megaset 𝔐(𝕄Conics).

Now recall that finding the best-fitting conic corresponds to an orthogonal projection of the given megapoint 𝒫 in the megaspace ℝ^2𝑛 (in our illustration, it would be a point (𝑥,𝑦,𝑧)∈ℝ³) onto the megaset 𝔐(𝕄_Conics) (in our illustration—onto the 𝑥𝑦 plane). In Figure 7 the point (𝑥,𝑦,𝑧) is simply projected onto (𝑥,𝑦,0). 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 {𝑦=0} are projected onto the line {𝑦=𝑧=0}. 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 𝐻₊={𝑦>0,𝑧=0}. Finding the best-fitting ellipse would correspond to an orthogonal projection of the given point (𝑥,𝑦,𝑧)∈ℝ³ onto the positive half-plane 𝐻₊={𝑦>0,𝑧=0}. Now if the given point (𝑥,𝑦,𝑧) has a positive 𝑦-coordinate, then it is projected onto (𝑥,𝑦,0), as before, and we get the desired best-fitting ellipse. But if it has a negative 𝑦-coordinate, then it is projected onto (𝑥,0,0), 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 𝑦<0 (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).

Figure 8 

Schematic illustration of the megaset 𝔐(𝕄𝐶)⊂𝔐(𝕄𝐸).

Suppose that the megaset 𝔐(𝕄_𝐶) for the collection of circles is represented by the set 𝑈={𝑦=𝑥²,𝑥≠0,𝑧=0} in our illustration. Note that 𝑈 consists of two curves (branches of a parabola on the 𝑥𝑦 plane), both lie in the half-plane 𝐻₊={𝑦>0,𝑧=0} that corresponds to the collection of ellipses. So the required inclusion 𝑈⊂𝐻₊ does take place. The two curves making 𝑈 terminate at the point (0,0,0), 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 (0,0,0)? It is not hard to see (and prove by elementary geometry) that only points on the 𝑦𝑧 plane may be projected onto the origin (0,0,0) (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 𝑛=5 (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 𝑄₂,𝑄₃,𝑄₄,and𝑄₅. Then let 𝑄^1 be a point near 𝑄₁ such that dist(𝑄^1,𝑃₁)<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 𝑄^1,𝑄₂,𝑄₃,𝑄₄,and𝑄₅. It is easy to see that disdtist𝑃_𝑖𝑃₁,𝑆^≤dist𝑃₁,𝑄^1<dist𝑃₁,𝑄₁=dist𝑃₁,𝑆,,𝑆^≤dist𝑃_𝑖,𝑄_𝑖=dist𝑃_𝑖,𝑆forevery𝑖=2,…,5.(A.1) As a result, ℱ𝑆^=⁵_𝑖=1dist𝑃_𝑖,𝑆^²<⁵_𝑖=1dist𝑃_𝑖,𝑆²=ℱ(𝑆),(A.2) as desired. Since 𝑄^1 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 dist𝑃₁,𝑆^<dist𝑃₁,𝑄₁=dist𝑃₁,𝑆.(A.3) At the same time we have dist𝑃_𝑖,𝑆^≤dist𝑃_𝑖,𝑄_𝑖=dist𝑃_𝑖,𝑆forevery𝑖=2,…,5.(A.4) As a result, consider ℱ𝑆^=⁵_𝑖=1dist𝑃_𝑖,𝑆^²<⁵_𝑖=1dist𝑃_𝑖,𝑆²=ℱ(𝑆),(A.5) 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.