1 Introduction

In metrology applications (e.g., quality and wearing assessment of industrial rollers) an object is considered round enough, if all measurement points on its surface lie between two spheres (cirlces) around a common center, such that the difference of their radii is smaller than some prescribed tolerance. This complies with ANSI Y14.5 and ISO 1101 standards [1, 2].

The roundness of any set \(S \subset {\mathbb {R}}^d\) may be defined as the smallest possible difference of radii for a suitable center. A precise definition of the corresponding optimization problem will be given in the next section.

The problem of computing the roundness is also known as the problem of computing the minimum radial separation, the minimum zone sphericity (circularity) error, or the minimum-width shell (annulus).

Computational geometry algorithms that compute exact solutions for finite sets S with n points are mostly based on Voronoi diagrams [3,4,5,6,7]. In the plane, they have worst case \({\mathcal {O}}(n^2)\)-time behavior, but there are also subquadratic \({\mathcal {O}}(n^{3/2+\beta })\)-time algorithms [8] for arbitrarily small \(\beta >0\). Under special assumptions on the set S even (expected) linear time can be achieved with linear programming techniques [9], see also [10]. In [11] it was noted that the general complexity to compute an exact solution in any dimension d is \({\mathcal {O}}(n^{{{\,\mathrm{floor}\,}}(d/2)+1})\). No estimates on the constants involved are given, but they generally depend exponentially on d. Indeed, as pointed out in [12], even in the plane, algorithms, that compute exact global solutions, are rather inefficient for large n. Therefore, it is of practical interest to develop algorithms, that efficiently compute approximate global minimizers of the roundness problem for any prescribed accuracy \(\varepsilon >0\), especially, if the roundness of many manufactured objects has to be evaluated. Since the roundness problem may have several local minima, standard local optimization methods cannot guarantee to find global minimizers. To our knowledge, the algorithm with the best complexity estimate to compute an approximate global minimizer is the one in [11] with \({\mathcal {O}}\big (n \cdot \ln (1/\varepsilon )/\varepsilon ^d\big )\), which for \(d=2\) improves the \({\mathcal {O}}\big (n \cdot \ln (n)+ \frac{n}{\varepsilon ^2}\big )\) result of [12]. Again, no estimates on the constants involved are given. Furthermore, the case \(d=2\) is often treated in a special manner, and the algorithms cannot directly be used to deal with sets S, that contain infinitely many points such as line segments. But, for the case that S is the boundary of a convex polygon with n vertices in the plane, an exact solution can also be computed with the help of Voronoi diagrams in \({\mathcal {O}}(n)\)-time [13].

Here, we analyze a global search algorithm to solve the roundness problem, which adaptively divides the initial search space into smaller parts until an approximate global minimizer with the desired accuracy is found. The main idea is in the spirit of the general branch-and-bound methods for global optimization discussed in [14], but here we make use of the special structure of the roundness problem to reduce the search space. The algorithm is simple to implement, because it is based only on the comparison of objective function values, and neither derivatives nor Voronoi diagrams have to be computed, which also makes it less sensitive to perturbations in the data. The formulation is valid in any dimension d, and for both finite sets S of data points as well as infinite sets like polygonal chains or triangulations of surfaces. An approximate global minimizer is found in \({\mathcal {O}}\big (n \cdot \ln (1/\varepsilon ) \big )\)-time (with n beeing the number of points or vertices), which also improves the \({\mathcal {O}}\big (n \cdot \ln (1/\varepsilon )/\varepsilon ^d\big )\) result of [11]. The constants involved depend on the number of objective function evaluations (the most costly operation), which in the worst case also grows exponentially with the dimension d. But, under mild assumptions on the roundness of S and the initial search space, we compute explicit upper bounds on the number of objective function evaluations that are needed to find an approximate global minimizer, whose objective value is not larger than twice the minimal value. Similar ideas might be used to formulate and analyze global search algorithms for related problems like cylindricity or flatness problems under appropriate assumptions on the data.

In the next section we give a precise definition of the roundness problem, and recall a result proven in [9], which shows that global minimizers of the roundness problem exist for almost round sets. We will need this result for the complexity estimate of our algorithm, which is analyzed in Sect. 3. The performance of the algorithm is tested in Sect. 4 with various numerical experiments.

2 Global Solutions of the Roundness Problem

For vectors \(\varvec{x}=(x_1,\ldots ,x_d), \varvec{y}=(y_1,\ldots ,y_d) \in {\mathbb {R}}^d\) we denote their scalar product by \(\langle \varvec{x} \,,\, \varvec{y}\rangle :=\sum _{i=1}^d x_i \cdot y_i\), and the Euclidean norm by \(\Vert \varvec{x}\Vert _{}:=\sqrt{\sum _{i=1}^d x_i^2}\). By introducing the roundness function for a given set \(S \subset {\mathbb {R}}^d\),

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{x}):=\max _{\varvec{p} \in S}\Vert \varvec{x}-\varvec{p}\Vert _{} - \min _{\varvec{q} \in S}\Vert \varvec{x}-\varvec{q}\Vert _{} \,, \end{aligned}$$

the problem of computing the roundness \({{\,\mathrm{rd}\,}}(S)\) of S may be formulated as the unconstrained optimization problem

$$\begin{aligned} {{\,\mathrm{rd}\,}}(S):=\min _{\varvec{x} \in {\mathbb {R}}^d} {{\,\mathrm{rd}\,}}_S(\varvec{x}) \,. \end{aligned}$$
(RP)

We assume throughout the paper that \(S \subset {\mathbb {R}}^d\) is compact. This implies that the maximum distance \(r_{\max }(\varvec{x}):=\max _{\varvec{p} \in S}\Vert \varvec{x}-\varvec{p}\Vert _{}\) and the minimum distance \(r_{\min }(\varvec{x}):=\min _{\varvec{q} \in S}\Vert \varvec{x}-\varvec{q}\Vert _{}\) of a given center \(\varvec{x} \in {\mathbb {R}}^d\) to the set S are attained, and that the roundness function \({{\,\mathrm{rd}\,}}_S(\varvec{x})=r_{\max }(\varvec{x})-r_{\min }(\varvec{x})\) is continuous. Even so, it was shown in [15, 16] that the minimum in (RP) need not exist in general. As a simple counterexample, consider the case that S consists of 3 distinct points on a straight line: For any center \(\varvec{x} \in {\mathbb {R}}^d\) we have \({{\,\mathrm{rd}\,}}_S(\varvec{x}) >0\), but, by moving the center perpendicularly to and further away from the line, in the limit we get \(\inf _{\varvec{x} \in {\mathbb {R}}^d} {{\,\mathrm{rd}\,}}_S(\varvec{x}) =0\). But a positive result was obtained in [9] for sets which there are called almost round. Here we need this result for slightly more general parameter values. Similar to [9], in addition to compactness, we make the following two assumptions about the set S with respect to some given center \(\varvec{x_0}\in {\mathbb {R}}^d\) and angle \(\alpha \in [0,\pi [\).

  1. (A1)

    The set S contains sufficiently many points sucht that there is at least one point of S inside any cone with apex at \(\varvec{x_0}\) and angle \(\alpha \); see Fig. 1.

  2. (A2)

    The ratio \(\eta \) of the value of the roundness function

    $$\begin{aligned} {{{\,\mathrm{rd}\,}}_S(\varvec{x_0})=r_{\max }(\varvec{x_0})-r_{\min }(\varvec{x_0})} \end{aligned}$$

    to the mean radius \(r(\varvec{x_0}):=\frac{r_{\max }(\varvec{x_0})+r_{\min }(\varvec{x_0})}{2}\) is sufficiently small; more precisely it holds that

    $$\begin{aligned} \eta :=\frac{{{\,\mathrm{rd}\,}}_S(\varvec{x_0})}{r(\varvec{x_0})} < 2 \cdot \cos \left( \tfrac{\alpha }{2}\right) \,. \end{aligned}$$
Fig. 1
figure 1

Left. Blue: Data points of an almost round set S. Red: Each cone with apex at \(\varvec{x_0}\) and angle \(\alpha \) contains at least one point of S. Right. Typical graph of the roundness function \({{\,\mathrm{rd}\,}}_S(\varvec{x})\) of an almost round set S

Full size image

Small values of \(\alpha \) in  Assumption (A1) prevent “big holes” in the data. The ideal case \(\alpha =0\) in particular corresponds to closed surfaces, for instance if S is the boundary of a polygon in the plane. Larger values of \(\alpha \) force the ratio \(\eta \) in Assumption (A2) to be smaller. The next theorem shows that under these assumptions the global minimum exists, and moreover, the search for a minimizer can be restricted to a ball around \(\varvec{x_0}\).

Theorem 2.1

(cf. Lemma 1 in [9]) Under Assumptions (A1) and (A2) the roundness problem (RP) has a global minimum, which is attained at some center in a ball around \(\varvec{x_0}\) with radius

$$\begin{aligned} \rho :=\frac{\cos \left( \tfrac{\alpha }{2}\right) \cdot (4-\eta ^2)}{4 \cdot \cos ^2\left( \tfrac{\alpha }{2}\right) -\eta ^2} \cdot {{\,\mathrm{rd}\,}}_S(\varvec{x_0}) \,. \end{aligned}$$

If Assumption (A1) holds with \(\alpha =0\), then Assumption (A2) is not necessary, and the assertion is valid with \(\rho :={{\,\mathrm{rd}\,}}_S(\varvec{x_0})\).

Proof

At first, we note that Assumption (A2) implies \(0 \le \rho < \infty \) in case \(\alpha >0\). Let \(\varvec{x}\in {\mathbb {R}}^d\) be any center with \(\Vert \varvec{x}-\varvec{x_0}\Vert _{} \ge \rho \). We will show that the inequality \({{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge {{\,\mathrm{rd}\,}}_S(\varvec{x_0})\) holds. Since the roundness function is continuous, its global minimum must then be attained on the ball around \(\varvec{x_0}\) with radius \(\rho \). By  Assumption (A1) there exist \(\varvec{p},\varvec{q} \in S\) with \(\frac{\langle \varvec{q}-\varvec{x_0} \,,\, \varvec{x}-\varvec{x_0}\rangle }{\Vert \varvec{q}-\varvec{x_0}\Vert _{}\cdot \Vert \varvec{x}-\varvec{x_0}\Vert _{}} \ge \cos \left( \tfrac{\alpha }{2}\right) \) and \(\frac{\langle \varvec{p}-\varvec{x_0} \,,\, -(\varvec{x}-\varvec{x_0})\rangle }{\Vert \varvec{p}-\varvec{x_0}\Vert _{}\cdot \Vert \varvec{x}-\varvec{x_0}\Vert _{}} \ge \cos \left( \tfrac{\alpha }{2}\right) \). We set \(t:=\Vert \varvec{x}-\varvec{x_0}\Vert _{}\) and estimate

$$\begin{aligned} \Vert \varvec{p}-\varvec{x}\Vert _{}= \Vert (\varvec{p}-\varvec{x_0})-(\varvec{x}-\varvec{x_0})\Vert _{} \ge \sqrt{r_{\min }^2 + 2\cos \left( \tfrac{\alpha }{2}\right) \cdot r_{\min } \cdot t +t^2} \end{aligned}$$

and

$$\begin{aligned} \Vert \varvec{q}-\varvec{x}\Vert _{} = \Vert (\varvec{q}-\varvec{x_0})-(\varvec{x}-\varvec{x_0})\Vert _{} \le \sqrt{r_{\max }^2 - 2\cos \left( \tfrac{\alpha }{2}\right) \cdot r_{\max } \cdot t +t^2 }\,. \end{aligned}$$

For \(t\ge 0\) the function

$$\begin{aligned} f(t):=\sqrt{r_{\min }^2 + 2\cos \left( \tfrac{\alpha }{2}\right) \cdot r_{\min } \cdot t +t^2 }-\sqrt{r_{\max }^2 - 2\cos \left( \tfrac{\alpha }{2}\right) \cdot r_{\max } \cdot t +t^2 } \end{aligned}$$
(1)

is increasing with \(f(\rho )=r_{\max }-r_{\min }={{\,\mathrm{rd}\,}}_S(\varvec{x_0})\). For \(\Vert \varvec{x}-\varvec{x_0}\Vert _{}=t \ge \rho \) it follows that we get

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge \Vert \varvec{p}-\varvec{x}\Vert _{}- \Vert \varvec{q}-\varvec{x}\Vert _{} \ge f(t) \ge f(\rho )={{\,\mathrm{rd}\,}}_S(\varvec{x_0}) \,. \end{aligned}$$

\(\square \)

We will make use of this theorem in two ways. First, an appropriate initial center \(\varvec{x_0}\) allows us to restrict the initial search space for global minimizers. For example, if \(\alpha \le \frac{\pi }{2}\) and \(\eta \le 0.4\), then it suffices to search in a ball around \(\varvec{x_0}\) with radius \(\rho < \frac{3}{2}\cdot {{\,\mathrm{rd}\,}}_S(\varvec{x_0})\). This may be useful even if the roundness value \({{\,\mathrm{rd}\,}}_S(\varvec{x_0})\) at the initial center is not very small. Second, by applying the theorem with \(\varvec{x_0}\) being a global minimizer, we will derive complexity estimates for our search algorithm in the next section.

To ensure existence of the global minimum in general, it suffices to restrict the choice of possible centers to a compact set \(Q \subset {\mathbb {R}}^d\). Therefore, in the following, we consider the constrained roundness problem

$$\begin{aligned} {{\,\mathrm{rd}\,}}(S|Q):=\min _{\varvec{x} \in Q} {{\,\mathrm{rd}\,}}_S(\varvec{x}) \,, \end{aligned}$$
(CRP)

According to Theorem 2.1, for almost round sets we can choose Q large enough so that solutions of (CRP) are also solutions of (RP). Note that the roundness function in general is not convex. Hence it may have several local minima, and the global minimum may be attained at different centers. Uniqueness of the optimal center cannot so easily be guaranteed, even for almost round sets. As a pathological example, consider the case that S consists of the origin together with a perfect sphere with radius r around the origin: Then the roundness of S is \({{\,\mathrm{rd}\,}}(S)=r\), and the global minimum is attained at every point in the ball with radius \(\frac{r}{2}\) around the origin. Since we are mainly interested in the minimal objective value, we do not pursue any further the issue of uniqueness, but instead refer to [9, 17] for a positive result in dimension \(d=2\) and a discussion why for almost round sets it is unlikely in practice to encounter local minima which are not also global minima.

3 Roundness Evaluation Algorithm

We propose to solve the constrained roundness problem (CRP) with a global search algorithm, which adaptively divides the initial search space into smaller parts until an approximate global minimizer with the desired accuracy is found. It consists of the basic global search Algorithm 1 together with Algorithm 2 to speed up convergence. For reasons of simplicity and efficiency, we initially use a cube (square) Q in (CRP), which is then adaptively divided into smaller cubes in each iteration. The centers of these cubes are candidates for potential minimizers.

3.1 Basic Global Search

We denote a d-dimensional axis-parallel cube with center \(\varvec{x} \in {\mathbb {R}}^d\) and edge length \(h>0\) by

$$\begin{aligned} Q_h(\varvec{x}):=\{\varvec{y} \in {\mathbb {R}}^d\, :\, \max _{i=1,\ldots ,d}|y_i-x_i| \le \tfrac{h}{2}\} \,, \end{aligned}$$

and its set of \(2^d\) vertices by \(V_h(\varvec{x})\). The vertices in \(V_{h/2}(\varvec{x})\) then coincide with the centers of the subcubes that we obtain by dividing \(Q_h(\varvec{x})\) into \(2^d\) subcubes with edge length \(\frac{h}{2}\), see Fig. 2.

Fig. 2
figure 2

Black: Cube \(Q_{h}(\varvec{x})\) with center \(\varvec{x}\) and edge length h. Red: Vertices \(\varvec{v_i} \in V_{h/2}(\varvec{x})\) are centers of the subcubes with edge length \(\frac{h}{2}\)

Full size image

We will prove the convergence of the algorithm with the help of the following lemma, which gives a sharp bound on the growth of the roundness function on a cube with respect to its center.

Lemma 3.1

For all \(\varvec{y} \in Q_h(\varvec{x})\) we have \({{\,\mathrm{rd}\,}}_S(\varvec{x}) \le {{\,\mathrm{rd}\,}}_S(\varvec{y}) +\sqrt{d} \cdot h\).

Proof

To \(\varvec{x}\) we find \(\varvec{p},\varvec{q} \in S\) such that \(\Vert \varvec{x}-\varvec{p}\Vert _{}-\Vert \varvec{x}-\varvec{q}\Vert _{} = {{\,\mathrm{rd}\,}}_S(\varvec{x})\). Using the triangle inequality, we get \({{\,\mathrm{rd}\,}}_S(\varvec{y}) \ge \Vert \varvec{y}-\varvec{p}\Vert _{}-\Vert \varvec{y}-\varvec{q}\Vert _{} \ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) - 2\cdot \Vert \varvec{y}-\varvec{x}\Vert _{}\). Now the assertion follows from the estimate \(\displaystyle \Vert \varvec{y}-\varvec{x}\Vert _{} \le \sqrt{d} \cdot \max _{i=1,\ldots ,d}|y_i-x_i| \le \sqrt{d} \cdot \tfrac{h}{2}\) for \(\varvec{y} \in Q_h(\varvec{x})\). \(\square \)

In each iteration of Algorithm 1, we use a finite set C of centers of cubes \(Q_h(\varvec{x})\) which have the same current edge length h. These centers are candidates for potential minimizers. The center \(\varvec{{\hat{x}}}\) with so far minimal roundness value \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\) is kept as current approximate solution. In line 1 of Algorithm 1 we initialize with \(\varvec{{\hat{x}}}:=\varvec{x_0}\), \(h:=h_0\), and \(C:=\{\varvec{x_0}\}\). New candidate centers are obtained as follows: We start with \(C_\mathrm{new}:=\emptyset \) (line 3) and run through all current centers \(\varvec{x} \in C\) (line 4). Based on Lemma 3.1, we make in line 5 what we call a simple test to decide which of the current cubes can be discarded in the next iterations, and which cubes have to be searched further. If the current center \(\varvec{x} \in C\) passes the test, then we add all centers \(\varvec{x_\mathrm{new}} \in V_{h/2}(\varvec{x})\) of the subcubes to \(C_\mathrm{new}\) (lines 6, 7). Furthermore, if one of the new centers has a smaller roundness value than \(\varvec{{\hat{x}}}\), then we take it as the new approximate solution (lines 8, 9). In the next iteration, we replace C with the new centers (line 14), and halve the edge length h (line 15). The iterations are stopped as soon as the desired accuracy \({{\,\mathrm{rd}\,}}(\varvec{{\hat{x}}})\le {{\,\mathrm{rd}\,}}(S|Q)+\varepsilon \) is reached, which by the next theorem will be guaranteed to be the case if either the current edge length is small enough, or no cubes remain to be searched (line 16).

figure a

Theorem 3.1

Given an initial center \(\varvec{x_0}\in {\mathbb {R}}^d\), edge length \(h_0>0\), and desired accuracy \(\varepsilon >0\), the basic global search Algorithm 1 finds an approximate global minimizer \(\hat{\varvec{x}}\) of the constrained roundness problem (CRP) on the cube \(Q:=Q_{h_0}(\varvec{x_0})\) with \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\le {{\,\mathrm{rd}\,}}(S|Q)+\varepsilon \) after at most \({{\,\mathrm{ceil}\,}}\left( \log _2\left( \frac{\sqrt{d} \cdot h_0}{\varepsilon }\right) \right) \) iterations.

Proof

We denote the union of all cubes with current centers C and edge length h at the beginning of an iteration by

$$\begin{aligned} R_{C,h}:=\bigcup _{\varvec{x} \in C} Q_{h}(\varvec{x}) \,. \end{aligned}$$

By Lemma 3.1, for all \(\varvec{y} \in R_{C,h}\) and corresponding center \(\varvec{x} \in C\) such that \(\varvec{y} \in Q_{h}(\varvec{x})\), we have \({{\,\mathrm{rd}\,}}_S(\varvec{x}) \le {{\,\mathrm{rd}\,}}_S(\varvec{y})+\sqrt{d} \cdot h\). Since the current approximate solution fulfills \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{x})\), we infer that

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{y})+\sqrt{d} \cdot h \quad \text{ for } \text{ all } \varvec{y} \in R_{C,h}. \end{aligned}$$
(2)

We prove inductively that on the complement of the region \(R_{C,h}\) we have

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{y}) + \varepsilon \quad \text{ for } \text{ all } \varvec{y} \in Q\setminus R_{C,h}. \end{aligned}$$
(3)

(This also shows that outside of the region \(R_{C,h}\) the current approximate solution \(\varvec{{\hat{x}}}\) has already reached the desired accuracy, so that the search for new, potentially better, approximate minimizers can be restricted to the region \(R_{C,h}\), see also Fig. 3.)

Initially we have \(C=\{\varvec{x_0}\}\) and \(h=h_0\) (cf. line 1 of Algorithm 1), so that \(R_{C,h}=Q_{h_0}(\varvec{{\hat{x}}})=Q\). Hence, (3) holds trivially since \(Q\setminus R_{C,h}=\emptyset \). Now we assume that (3) holds for the current centers C and edge length h at the beginning of the for-loop in line 4, and show that it also holds for \(R_{C_\mathrm{new},h/2}\) at the end (line 13) with the new centers \(C_\mathrm{new}\) and edge length \(\frac{h}{2}\). Let \(\varvec{y} \in Q\setminus R_{C_\mathrm{new},h/2}\) be arbitrary. If \(\varvec{y}\) is also contained in \(Q\setminus R_{C,h}\), then the induction hypothesis (3) immediately implies \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{y}) + \varepsilon \). Otherwise, we have \(\varvec{y} \in R_{C,h}\). According to the construction of the set of new centers \(C_\mathrm{new}\) (lines 6, 7) this means that \(\varvec{y} \in Q_{h}(\varvec{x})\) for some center \(\varvec{x} \in C\) which must have failed the test in line 5, i.e., \({{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) + \sqrt{d} \cdot h-\varepsilon \). Hence, together with Lemma 3.1, in this case we also get

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{x}) - \sqrt{d} \cdot h + \varepsilon \le {{\,\mathrm{rd}\,}}_S(\varvec{y}) + \varepsilon \,, \end{aligned}$$

so that (3) indeed holds for all \(\varvec{y} \in Q\setminus R_{C_\mathrm{new},h/2}\).

Note that if \(C_\mathrm{new}=\emptyset \) then the new region \(R_{C_\mathrm{new},h/2}\) is empty, too, so that by (3) we can stop iterating (line 16), because then \(\varvec{{\hat{x}}}\) already is an approximate minimizer with the desired accuracy.

By (2) and (3), we conclude that the current approximate solution fulfills

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\le {{\,\mathrm{rd}\,}}(S|Q)+ \max \{\varepsilon ,\sqrt{d} \cdot h\} \,, \end{aligned}$$
(4)

so that we can also stop iterating as soon as the edge length h of the cubes is small enough, i.e., \(\sqrt{d} \cdot h \le \varepsilon \) (line 16). Since at the end of the K-th iteration the edge length of the cubes is \(h=\frac{h_0}{2^K}\), the desired accuracy \(\sqrt{d} \cdot h \le \varepsilon \) is reached after at most \({{\,\mathrm{ceil}\,}}\left( \log _2\left( \frac{\sqrt{d} \cdot h_0}{\varepsilon }\right) \right) \) iterations. \(\square \)

Fig. 3
figure 3

Left. Blue: Data points S, approximately on a circular arc. Red: Currently best center \(\varvec{{\hat{x}}}\) after some iterations with Algorithm 1, and region \(R_{C,h}\) where centers \(\varvec{x}\) with \({{\,\mathrm{rd}\,}}_S(\varvec{x})<{{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) - \varepsilon \) may be found. Right. Typical graph of the roundness function \({{\,\mathrm{rd}\,}}_S(\varvec{x})\) of an approximate circular arc S

Full size image

The most costly operation is the evaluation of the roundness function for the encountered centers. In the worst case, no cubes can be discarded, i.e., after K iterations there are at most \(\sum _{k=0}^K (2^d)^k =\frac{2^{(K+1)\cdot d}-1}{2^d-1}\) function evaluations. On the one hand, this implies that the complexity of the algorithm to compute an approximate solution for finite sets S with n data points is \({\mathcal {O}}(n)\), i.e., it increases only linearly with the number n of data points. On the other hand, there is a potentially huge constant hidden in the \({\mathcal {O}}\)-notation, and for higher dimensions even storing the current centers in computer memory will be a problem. But Theorem 2.1 indicates that many cubes will be discarded up to a small neighborhood around a minimizer. The next lemma shows how to derive much smaller upper bounds on the number of roundness function evaluations for almost round sets S by counting the number of certain integer lattice points.

Lemma 3.2

Let \(\varvec{x_{\min }} \in Q\) be a globlal minimizer of the constrained roundness problem (CRP), and let Assumptions (A1) and (A2) be fullfilled for \(\varvec{x_{\min }}\) in place of \(\varvec{x_0}\) with corresponding mean radius \(r=\frac{r_{\max }(\varvec{x_{\min }})+r_{\min }(\varvec{x_{\min }})}{2}\) and ratio \(\eta =\frac{{{\,\mathrm{rd}\,}}_S(\varvec{x_{\min }})}{r}\). For simplicity we assume that Algorithm 1 is run with initial edge length \(h_0=r\) and desired accuracy \(\varepsilon \le {{\,\mathrm{rd}\,}}(S|Q)\). Then after at most \(K:={{\,\mathrm{ceil}\,}}\left( \log _2\left( \frac{\sqrt{d}}{\eta }\right) \right) \) iterations the current approximate solution fulfills

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\le 2 \cdot {{\,\mathrm{rd}\,}}(S|Q) \,. \end{aligned}$$

We set \(c:=\cos \left( \tfrac{\alpha }{2}\right) \), and for \(k=0,\ldots ,K\) and integer lattice points \(\varvec{z} \in {\mathbb {Z}}^d\) we define

$$\begin{aligned} g_{k}(\varvec{z})&:= \sqrt{\big (1-\tfrac{\eta }{2}\big )^2 + 2 c\cdot \big (1-\tfrac{\eta }{2}\big ) \cdot \tfrac{\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2} }{2^{k-1}} +\tfrac{\big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big )^2 }{2^{2k-2}}}\nonumber \\&\quad -\sqrt{\big (1+\tfrac{\eta }{2}\big )^2 - 2 c\cdot \big (1+\tfrac{\eta }{2}\big ) \cdot \tfrac{\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2} }{2^{k-1}} +\tfrac{\big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big )^2 }{2^{2k-2}}}\,. \end{aligned}$$
(5)

Let \(a_k\) be the number of lattice points \(\varvec{z} \in {\mathbb {Z}}^d\) which fulfill either \(\Vert \varvec{z}\Vert _{} < \tfrac{\sqrt{d}}{2}\), or \(g_{k}(\varvec{z}) < \eta + \tfrac{3\cdot \sqrt{d}}{2^{k-1}}\) in case \(\Vert \varvec{z}\Vert _{} \ge \tfrac{\sqrt{d}}{2}\) (it may be \(a_k=\infty \)). Set \(m_0:=1\) and \(m_k:=\min \{2^d \cdot m_{k-1},a_k\}\) for \(k\ge 1\). Then the number of roundness function evaluations in iteration k is at most \(m_k\), i.e., after K iterations there are at most \(\sum _{k=0}^K m_k\) roundness function evaluations.

Proof

Let \(\varvec{{\hat{x}}}\) be the current approximate solution and \(\varvec{x} \in C\) be any of the centers in line 4 of Algorithm 1 at the beginning of iteration \(k \ge 1\). The cube \(Q_h(\varvec{x})\) is discarded if the test in line 5 fails, i.e., if

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) + \sqrt{d} \cdot h-\varepsilon \,. \end{aligned}$$
(6)

We will derive a sufficient condition for the test to fail.

Let us divide Q into equally sized subcubes, which have current edge length \(h=\frac{h_0}{2^{k-1}}=\frac{r}{2^{k-1}}\) (at the beginning of iteration k). The minimizer must lie in one of those cubes, say \(\varvec{x_{\min }} \in Q_h(\varvec{{\tilde{x}}})\) for some suitable center \(\varvec{{\tilde{x}}} \in Q\) (not necessarily \(\varvec{{\tilde{x}}} \in C\)). By Lemma 3.1, we get

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\tilde{x}}}) \le {{\,\mathrm{rd}\,}}_S(\varvec{x_{\min }}) + \sqrt{d} \cdot h ={{\,\mathrm{rd}\,}}(S|Q) + \sqrt{d} \cdot h \,. \end{aligned}$$
(7)

Furthermore, as was shown in (4) in the proof of Theorem 3.1, we have

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le {{\,\mathrm{rd}\,}}(S|Q) + \max \{\varepsilon ,\sqrt{d} \cdot h\} \le rd_S(\varvec{{\tilde{x}}}) + \max \{\varepsilon ,\sqrt{d} \cdot h\} \,. \end{aligned}$$
(8)

By definition of K it is the smallest integer with \(2^K \ge \frac{\sqrt{d}}{\eta }\). Hence, for \(k \le K\) we get \(\varepsilon \le {{\,\mathrm{rd}\,}}(S|Q)=\eta \cdot r < \sqrt{d} \cdot \tfrac{r}{2^{k-1}} =\sqrt{d} \cdot h\), so that (8) and (7) yield the estimate

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le rd_S(\varvec{{\tilde{x}}}) +\sqrt{d} \cdot h \le {{\,\mathrm{rd}\,}}(S|Q) + 2 \cdot \sqrt{d} \cdot h \,. \end{aligned}$$
(9)

For the current center \(\varvec{x} \in C\), we can write \(\varvec{x}-\varvec{{\tilde{x}}} = h \cdot \varvec{z}\) with some suitable \(\varvec{z} \in {\mathbb {Z}}^d\). Since \(\varvec{x_{\min }} \in Q_h(\varvec{{\tilde{x}}})\) we conclude that

$$\begin{aligned} t:=\Vert \varvec{x} -\varvec{x_{\min }}\Vert _{} \ge \Vert \varvec{x} -\varvec{{\tilde{x}}}\Vert _{} - \Vert \varvec{{\tilde{x}}} -\varvec{x_{\min }}\Vert _{} \ge h \cdot \big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big ) =:t_z\,. \end{aligned}$$

For \(\Vert \varvec{z}\Vert _{} \ge \tfrac{\sqrt{d}}{2}\), we have \(t_z \ge 0\), so that we can argue as in the proof of Theorem 2.1 for \(\varvec{x_{\min }}\) instead of \(\varvec{x_0}\), and with the increasing function f(t) defined there in (1), we get

$$\begin{aligned} {{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge f\big (\Vert \varvec{x} -\varvec{x_{\min }}\Vert _{}\big ) = f(t)\ge f(t_z)=f\big (h \cdot \big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big )\big ) \,. \end{aligned}$$
(10)

Hence, by (9) and (10), the inequality (6) is surely fulfilled in case

$$\begin{aligned} f\big (h \cdot \big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big )\big ) \ge {{\,\mathrm{rd}\,}}(S|Q) + 3 \cdot \sqrt{d} \cdot h \,. \end{aligned}$$

By recalling that \({{\,\mathrm{rd}\,}}(S|Q)=\eta \cdot r\), \(h=\frac{r}{2^{k-1}}\), \(r_{\min }=\big (1-\frac{\eta }{2}\big ) \cdot r\), \(r_{\max }=\big (1+\frac{\eta }{2}\big ) \cdot r\), and the definition of \(g_k\) in (5), we conclude that the cube \(Q_h(\varvec{x})\) is discarded if

$$\begin{aligned} g_{k}(\varvec{z}) = \tfrac{f\big (h \cdot \big (\Vert \varvec{z}\Vert _{}-\tfrac{\sqrt{d}}{2}\big )\big )}{r} \ge \eta + \tfrac{3\cdot \sqrt{d}}{2^{k-1}} \,. \end{aligned}$$

Since the number of new centers is at most \(2^d\)-times the number of old centers, we infer that the number of roundness function evaluations in iteration k is at most \(m_k=\min \{2^d \cdot m_{k-1},a_k\}\) for \(k\ge 1\). Finally, because \(\varepsilon \le {{\,\mathrm{rd}\,}}(S|Q)\) and at the end of iteration K, we also have \(\sqrt{d} \cdot h=\sqrt{d} \cdot \frac{r}{2^{K}} \le {{\,\mathrm{rd}\,}}(S|Q)\), we infer from the first inequality in (8) that the approximate solution at the end of iteration K fulfills \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le 2 \cdot {{\,\mathrm{rd}\,}}(S|Q)\). \(\square \)

Some of the results that we computed with the help of Lemma 3.2 are given in Table 1, where the values in brackets correspond to the order of the general worst case estimates if no cubes would be discarded (note that the worst case estimate depends on K and thus on \(\eta \), but not on \(\alpha \)).

Table 1 Maximum number of roundness function evaluations after K iterations, computed with the help of Lemma 3.2
Full size table

Lemma 3.2 explains well the good initial behavior of Algorithm 1, which only uses the simple test in line 5 to decide which cubes can be discarded. Moreover, it shows that it is not really important that the initial search space is small, because the algorithm itself rapidly decreases the size of the initial search space.

3.2 Acceleration of the Basic Global Search

To speed up the convergence by increasing the number of discarded cubes, we use Algorithm 1 together with Algorithm 2.

figure b

Algorithm 2 augments the simple test by a second one, which exploits the special structure of the roundness function. If the simple test succeeded for the current center \(\varvec{x}\), then, instead of definitely keeping all \(2^d\) subcubes for the next iteration, the second test in line 5 of Algorithm 2 decides whether some of the subcubes can be discarded as well.

Lemma 3.3

If the additional test in line 5 of Algorithm 2 fails, then the subcube \(Q_{h/2}(\varvec{x_\mathrm{new}})\) can be discarded in subsequent iterations.

Proof

If the test in line 5 fails, then we have \(m \ge 0\) and \(m \ge {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) - \varepsilon \). By definition of m (line 4) all vertices \(\varvec{v} \in V_{h/2}(\varvec{x_\mathrm{new}})\) of the subcube \(Q_{h/2}(\varvec{x_\mathrm{new}})\) lie in the set \(H:=\{\varvec{y} \in {\mathbb {R}}^d\, :\, \Vert \varvec{y}-\varvec{p}\Vert _{}-\Vert \varvec{y}-\varvec{q}\Vert _{} \ge m\}\), which for \(m \ge 0\) is either a hyperboloid or a halfspace, and thus convex. This implies that the whole subcube \(Q_{h/2}(\varvec{x_\mathrm{new}})\) is contained in H. Hence, for all \(\varvec{y} \in Q_{h/2}(\varvec{x_\mathrm{new}})\) we get \({{\,\mathrm{rd}\,}}_S(\varvec{y}) \ge \Vert \varvec{y}-\varvec{p}\Vert _{}-\Vert \varvec{y}-\varvec{q}\Vert _{} \ge m \ge {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) - \varepsilon \), i.e., the subcube \(Q_{h/2}(\varvec{x_\mathrm{new}})\) can be discarded in subsequent iterations. \(\square \)

We point out that we do not evaluate the whole costly roundness function at the vertices, but only the values \(\Vert \varvec{v}-\varvec{p}\Vert _{}-\Vert \varvec{v}-\varvec{q}\Vert _{}\) for the single pair \(\varvec{p}, \varvec{q} \in S\) chosen in line 2 of Algorithm 2. Furthermore, since the subcubes have several vertices in common, and the value at the center \(\varvec{x}\) is already known, it actually suffices to compute the values \(\Vert \varvec{v}-\varvec{p}\Vert _{}-\Vert \varvec{v}-\varvec{q}\Vert _{}\) for only \(3^d-1\) vectors \(\varvec{v} \in {\mathbb {R}}^d\) (instead of \(2^d \cdot 2^d=4^d\)), cf. Fig. 4.

Fig. 4
figure 4

Red: Subcubes \(Q_{h/2}(\varvec{x_\mathrm{new}})\) of \(Q_{h}(\varvec{x})\) and their vertices. Blue: Data points \(\varvec{p},\varvec{q} \in S\) and region with \({{\,\mathrm{rd}\,}}_S(\varvec{y}) \ge \Vert \varvec{y}-\varvec{p}\Vert _{}-\Vert \varvec{y}-\varvec{q}\Vert _{}\ge {{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) - \varepsilon \) for all \(\varvec{y} \in Q_{h/2}(\varvec{x_\mathrm{new,2}})\), so that \(Q_{h/2}(\varvec{x_\mathrm{new,2}})\) can be discarded

Full size image

Obviously, the assertions of Theorem 3.1 and Lemma 3.2 remain valid, if we use Algorithm 1 together with Algorithm 2. But since the second test uses the smaller value \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) - \varepsilon \) instead of \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) + \sqrt{d} \cdot h - \varepsilon \), potentially more cubes can be discarded, especially in the final iterations, where the simple functions \(\varvec{y} \mapsto \Vert \varvec{y}-\varvec{p}\Vert _{}-\Vert \varvec{y}-\varvec{q}\Vert _{}\) with fixed \(\varvec{p},\varvec{q} \in S\) are locally good approximations of the roundness function. This claim is supported by the next lemma and our numerical experiments.

Lemma 3.4

Let r be the mean radius at a minimizer \(\varvec{x_{\min }}\) of (CRP). We make the following assumptions about the ratio \(\eta :=\frac{{{\,\mathrm{rd}\,}}(S|Q)}{r}\) and the desired accuracy \(\varepsilon \):

$$\begin{aligned} \eta \le \frac{1}{10} \quad {and} \quad 21 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) \le \varepsilon \le {{\,\mathrm{rd}\,}}(S|Q) \,. \end{aligned}$$
(11)

If the current edge length h in line 1 of Algorithm 2 is small enough, and the distance of the current center \(\varvec{x}\) to its corresponding data point \(\varvec{q} \in S\) chosen in line 2 is large enough such that

$$\begin{aligned} \sqrt{d} \cdot h \le {{\,\mathrm{rd}\,}}(S|Q) \quad {and} \quad \Vert \varvec{x}-\varvec{q}\Vert _{} \ge \tfrac{r}{2} \,, \end{aligned}$$
(12)

then at least one subcube of \(Q_{h}(\varvec{x})\) will be discarded by the test in line 5. In particular, the condition \(\Vert \varvec{x}-\varvec{q}\Vert _{} \ge \tfrac{r}{2}\) is automatically fulfilled, if \(\varvec{x}\) passed the test in line 1, and if Assumption (A1) holds with \(\alpha \le \frac{\pi }{2}\) for \(\varvec{x_{\min }}\) in place of \(\varvec{x_0}\).

Proof

As shown in (4) in the proof of Theorem 3.1, the current approximate solution fulfills \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\le {{\,\mathrm{rd}\,}}(S|Q)+ \max \{\varepsilon ,\sqrt{d} \cdot h\}\). Hence, (11) and (12) imply \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}}) \le 2 \cdot {{\,\mathrm{rd}\,}}(S|Q)\). If the current center \(\varvec{x}\) passed the first test in line 1 of Algorithm 2, then it follows that

$$\begin{aligned} {{\,\mathrm{rd}\,}}(\varvec{x}) < 3 \cdot {{\,\mathrm{rd}\,}}(S|Q)\,. \end{aligned}$$
(13)

We show that the second test in line 5 fails for the center \(\varvec{x_\mathrm{new}}\) of the subcube which lies in the halfspace \(H:=\{\varvec{y} \in {\mathbb {R}}^d\, :\, \langle \frac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}}-\frac{\varvec{x}-\varvec{q}}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \,,\, \varvec{y}-\varvec{x}\rangle \ge 0\}\) (note that at least one of the subcubes must lie in H because the center \(\varvec{x}\) lies in H). Each vertex \(\varvec{v} \in V_{h/2}(\varvec{x_\mathrm{new}})\) can be written as \(\varvec{v}=\varvec{x}+t \cdot \varvec{u}\) for some \(0 \le t \le \sqrt{d} \cdot h\) and a unit vector \(\varvec{u} \in {\mathbb {R}}^d\) with

$$\begin{aligned} \langle \tfrac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}}-\tfrac{\varvec{x}-\varvec{q}}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \,,\, \varvec{u}\rangle \ge 0 \,. \end{aligned}$$
(14)

We set \(h(t):=\Vert \varvec{x}+t \cdot \varvec{u}-\varvec{p}\Vert _{}-\Vert \varvec{x}+t \cdot \varvec{u}-\varvec{q}\Vert _{}\). To see that the test in line 5 fails, it is then sufficient to show that \(h(t) \ge {{\,\mathrm{rd}\,}}_S(\varvec{x})-\varepsilon \) (recall that \({{\,\mathrm{rd}\,}}_S(\varvec{x})-\varepsilon \ge {{\,\mathrm{rd}\,}}(S|Q) -\varepsilon \ge 0\)).

By Taylor expansion around \(t=0\) there exists some \({\tilde{t}} \in (0,t)\) and a function \(g({\tilde{t}})\) such that

$$\begin{aligned} h(t)&= {{\,\mathrm{rd}\,}}_S(\varvec{x}) + t \cdot \langle \tfrac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}}-\tfrac{\varvec{x}-\varvec{q}}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \,,\, \varvec{u}\rangle \cdot \left( 1- \tfrac{t}{2 \cdot \Vert \varvec{x}-\varvec{q}\Vert _{}} \cdot \langle \tfrac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}}+\tfrac{\varvec{x}-\varvec{q}}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \,,\, \varvec{u}\rangle \right) \nonumber \\&\quad -t^2 \cdot \tfrac{{{\,\mathrm{rd}\,}}_S(\varvec{x})}{2 \cdot \Vert \varvec{x}-\varvec{p}\Vert _{}\cdot \Vert \varvec{x}-\varvec{q}\Vert _{}} \cdot \left( 1- \tfrac{t}{2 \cdot \Vert \varvec{x}-\varvec{q}\Vert _{}} \cdot \left( \langle \tfrac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}} \,,\, \varvec{u}\rangle \right) ^2\right) - g({\tilde{t}}) \cdot t^3 \,. \end{aligned}$$
(15)

We will derive lower bounds for the summands in the Taylor expansion.

From (11) and (12) we infer that

$$\begin{aligned} {\tilde{t}} \le t \le \sqrt{d} \cdot h \le {{\,\mathrm{rd}\,}}(S|Q)=\eta \cdot r \le \frac{r}{10} \,, \end{aligned}$$
(16)

which implies \(\tfrac{t}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \le \frac{r/10}{r/2} = \frac{1}{5}\). It follows that the second summand in (15) is non-negative, because \(t\ge 0\) and by (14) we have \(0 \le \langle \frac{\varvec{x}-\varvec{p}}{\Vert \varvec{x}-\varvec{p}\Vert _{}}-\frac{\varvec{x}-\varvec{q}}{\Vert \varvec{x}-\varvec{q}\Vert _{}} \,,\, \varvec{u}\rangle \le 2\). Similarly we see that in the third summand the whole expression in the big brackets is less than or equal to 1, and by (13) and (16) we then can estimate the third summand from below by

$$\begin{aligned} -t^2 \cdot \tfrac{{{\,\mathrm{rd}\,}}_S(\varvec{x})}{2 \cdot \Vert \varvec{x}-\varvec{p}\Vert _{}\cdot \Vert \varvec{x}-\varvec{q}\Vert _{}} \ge -{{\,\mathrm{rd}\,}}(S|Q)^2 \cdot \tfrac{3 \cdot {{\,\mathrm{rd}\,}}(S|Q)}{2 \cdot r^2/4} = -6 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) \,. \end{aligned}$$

Hence we have

$$\begin{aligned} h(t) \ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) -6 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) - g({\tilde{t}}) \cdot t^3 \,. \end{aligned}$$
(17)

An elementary but lengthy calculation leads to the estimate

$$\begin{aligned} |g({\tilde{t}})| \le \frac{2}{\sqrt{3}} \cdot \left( \tfrac{1}{\Vert \varvec{x}-\varvec{p}+{\tilde{t}} \cdot \varvec{u}\Vert _{}^2} + \tfrac{1}{\Vert \varvec{x}-\varvec{q}+{\tilde{t}} \cdot \varvec{u}\Vert _{}^2}\right) \,. \end{aligned}$$

From (16) and (12) we infer that

$$\begin{aligned} \Vert \varvec{x}-\varvec{q}+{\tilde{t}} \cdot \varvec{u}\Vert _{} \ge \Vert \varvec{x}-\varvec{q}\Vert _{} - {\tilde{t}} \ge \tfrac{r}{2} - {\tilde{t}} \ge \tfrac{2r}{5}\,. \end{aligned}$$

Since \(\Vert \varvec{x}-\varvec{p}\Vert _{} \ge \Vert \varvec{x}-\varvec{q}\Vert _{}\) we likewise get \(\Vert \varvec{x}-\varvec{p}+{\tilde{t}} \cdot \varvec{u}\Vert _{} \ge \tfrac{2r}{5}\) and thus

$$\begin{aligned} |g({\tilde{t}})| \le \frac{2}{\sqrt{3}} \cdot \big (\tfrac{25}{4r^2} + \tfrac{25}{4r^2} \big ) =\tfrac{25}{\sqrt{3}} \cdot \tfrac{1}{r^2}\,. \end{aligned}$$

Together with (11) and (16) we finally get the desired lower bound for h(t) in (17), i.e.,

$$\begin{aligned} h(t)&\ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) -6 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) - \tfrac{25}{\sqrt{3}} \cdot \tfrac{t^3}{r^2}\\&\ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) -6 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) - \tfrac{25}{\sqrt{3}} \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q)\\&\ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) - 21 \cdot \eta ^2 \cdot {{\,\mathrm{rd}\,}}(S|Q) \ge {{\,\mathrm{rd}\,}}_S(\varvec{x}) -\varepsilon \,. \end{aligned}$$

It remains to show that the condition \(\Vert \varvec{x}-\varvec{q}\Vert _{} \ge \frac{r}{2}\) is fulfilled if Assumption (A1) holds with \(\alpha \le \frac{\pi }{2}\) for \(\varvec{x_{\min }}\) in place of \(\varvec{x_0}\). At first we note that by (11) we have

$$\begin{aligned} \Vert \varvec{x_{\min }}-\varvec{q}\Vert _{} \ge r_{\min }=r-\frac{{{\,\mathrm{rd}\,}}(S|Q)}{2} \ge \frac{19r}{20} \,. \end{aligned}$$
(18)

As we have seen above, (13) holds if \(\varvec{x}\) passed the first test in line 1 of Algorithm 2. We will show that this implies

$$\begin{aligned} \Vert \varvec{x}-\varvec{x_{\min }}\Vert _{} \le \frac{9r}{20} \,, \end{aligned}$$
(19)

because together with (18) we then can indeed conclude that

$$\begin{aligned} \Vert \varvec{x}-\varvec{q}\Vert _{} \ge \Vert \varvec{x_{\min }}-\varvec{q}\Vert _{}-\Vert \varvec{x}-\varvec{x_{\min }}\Vert _{} \ge \frac{19r}{20}-\frac{9r}{20} = \frac{r}{2}\,. \end{aligned}$$

Let us define

$$\begin{aligned} {\tilde{r}}_{\max } :=r+\frac{3}{2} \cdot {{\,\mathrm{rd}\,}}(S|Q) \ge r_{\max } \quad \text{ and } \quad {\tilde{r}}_{\min } :=r-\frac{3}{2} \cdot {{\,\mathrm{rd}\,}}(S|Q) \le r_{\min } \,,\\ {\tilde{{{\,\mathrm{rd}\,}}}}(S|Q) :={\tilde{{{\,\mathrm{rd}\,}}}}_S(\varvec{x_{\min }}):={\tilde{r}}_{\max }-{\tilde{r}}_{\min }=3\cdot {{\,\mathrm{rd}\,}}(S|Q) \,,\\ {\tilde{\eta }} :=\frac{{\tilde{{{\,\mathrm{rd}\,}}}}(S|Q)}{r} =3 \cdot \eta \quad \text{ and } \quad {\tilde{\rho }} :=\frac{\cos \left( \tfrac{\alpha }{2}\right) \cdot (4-{\tilde{\eta }}^2)}{4 \cdot \cos ^2\left( \tfrac{\alpha }{2}\right) -{\tilde{\eta }}^2} \cdot {\tilde{{{\,\mathrm{rd}\,}}}}(S|Q) \,. \end{aligned}$$

Condition (11) implies \({\tilde{\eta }} \le 0.3\), so that for \(\alpha \le \frac{\pi }{2}\) we have

$$\begin{aligned} {\tilde{\rho }}\le \frac{3}{2}\cdot {\tilde{{{\,\mathrm{rd}\,}}}}(S|Q)=\frac{9}{2} \cdot \eta \cdot r \le \frac{9r}{20}\,. \end{aligned}$$

Now assume that (19) does not hold, which implies \(\Vert \varvec{x}-\varvec{x_{\min }}\Vert _{} > \frac{9r}{20} \ge {\tilde{\rho }}\). But then a similar argument as in the proof of Theorem 2.1 leads to the estimate \({{\,\mathrm{rd}\,}}_S(\varvec{x}) \ge {\tilde{{{\,\mathrm{rd}\,}}}}_S(\varvec{x_{\min }})=3 \cdot {{\,\mathrm{rd}\,}}(S|Q)\), which contradicts (13). \(\square \)

3.3 Initial Center

In many applications, an adequate initial center is at hand. But if none is given, we compute an initial center for finite datasets \(S=\{\varvec{p_1},\ldots ,\varvec{p_n}\}\) by solving

$$\begin{aligned} \min _{\varvec{x} \in {\mathbb {R}}^d, r \ge 0} \sum _{i=1}^n \big (\Vert \varvec{p_i}-\varvec{x}\Vert _{}^2-r^2\big )^2 \,. \end{aligned}$$

In [18] it was shown that this problem is equivalent to the following linear least squares problem by setting \(t:=r^2-\Vert \varvec{x}\Vert _{}^2\),

$$\begin{aligned} \min _{\varvec{x} \in {\mathbb {R}}^d, t \in {\mathbb {R}}} \sum _{i=1}^n \big (\Vert \varvec{p_i}\Vert _{}^2-2 \cdot \langle \varvec{p_i} \,,\, \varvec{x}\rangle -t\big )^2 \,. \end{aligned}$$
(LSQ)

4 Numerical Experiments

We have implemented Algorithm 1 together with Algorithm 2 in Matlab (R2017b), and test its performance with different kind of sets S.

4.1 Random Sets with Predefined Roundness and Unique Optimal Center

At first we generate random sets of test points, but with predefined roundness and unique optimal center. The following lemma explains the construction of such sets.

Lemma 4.1

For given \(\delta \in ]0,1[\) we define the positive numbers \(u:=\sqrt{\frac{4+\delta ^2\cdot d}{4+4d}}\), \(v:=\sqrt{\frac{4d-\delta ^2 \cdot d}{4+4d}}\) and \(w:=\sqrt{d+1}\cdot \frac{\delta }{2}\). Let \(\varvec{e_i} \in {\mathbb {R}}^d\) for \(i=1,\ldots ,d\) be the standard basis vectors. For all \(i,j \in \{1,\ldots ,d\}\), with \(i\not =j\), and \(\sigma _u, \sigma _v, \sigma _w \in \{-1,+1\}\) we define the \(8\cdot d \cdot (d-1)\) different points \(\varvec{p_{i,j,\sigma _u, \sigma _v, \sigma _w}} \in {\mathbb {R}}^d\) by

$$\begin{aligned} \varvec{p_{i,j,\sigma _u, \sigma _v, \sigma _w}}:=\sigma _u \cdot (u+\sigma _w \cdot w) \cdot \varvec{e_i} + \sigma _v \cdot v \cdot \varvec{e_j} \,. \end{aligned}$$

Then all \(\varvec{p_{i,j,\sigma _u, \sigma _v, -1}}\) lie on the sphere with radius \(r_{\min }:=\sqrt{1+\frac{d}{4} \cdot \delta ^2}-\frac{\delta }{2}\) around the origin, and all points \(\varvec{p_{i,j,\sigma _u, \sigma _v, +1}}\) lie on the sphere with radius \(r_{\max }:=\sqrt{1+\frac{d}{4} \cdot \delta ^2}+\frac{\delta }{2}\) around the origin, cf. Fig. 5. The set S of all these points has roundness \({{\,\mathrm{rd}\,}}(S)=\delta \) with the origin as unique optimal center. Obviously, any superset of S obtained by adding arbitrarily many points \(\varvec{p} \in {\mathbb {R}}^d\) with \(r_{\min } \le \Vert \varvec{p}\Vert _{} \le r_{\max }\) to S then also has roundness \(\delta \) with the origin as unique optimal center.

Fig. 5
figure 5

Set S of test points \(\varvec{p_{i,j,\sigma _u, \sigma _v, \sigma _w}}\) for \(i,j \in \{1,2\}\) with \(i\not =j\) and \(\sigma _u, \sigma _v, \sigma _w \in \{-1,+1\}\), with predefined roundness \({{\,\mathrm{rd}\,}}(S)= 0.2\)

Full size image

Proof

Since \({{\,\mathrm{rd}\,}}_S(\varvec{0})=r_{\max }-r_{\min }=\delta \), we have \({{\,\mathrm{rd}\,}}(S)\le \delta \), and the assertion follows by proving that there can be no center \(\varvec{x} \in {\mathbb {R}}^d\) with \({{\,\mathrm{rd}\,}}_S(\varvec{x})< \delta \).

We define the functions \(h, h_{i,j,\sigma _u, \sigma _v}: {\mathbb {R}}^d \rightarrow {\mathbb {R}}\) by

$$\begin{aligned} h_{i,j,\sigma _u, \sigma _v}(\varvec{x}) := d \cdot (x_i-\sigma _u \cdot u)^2 - (x_j-\sigma _v \cdot v)^2 - \sum _{\begin{array}{c} k=1 \\ k\not =i,j \end{array}}^d x_k^2 \end{aligned}$$

and

$$\begin{aligned} h(\varvec{x}):=\sum _{\sigma _u, \sigma _v \in \{\pm 1\}} \sum _{\begin{array}{c} i,j=1 \\ i\not = j \end{array}}^d h_{i,j,\sigma _u, \sigma _v}(\varvec{x}) \,. \end{aligned}$$

Then the two points \(\varvec{p_{i,j,\sigma _u, \sigma _v, -1}}\) and \(\varvec{p_{i,j,\sigma _u, \sigma _v, +1}}\) are foci of the hyperboloid defined by the equation

$$\begin{aligned} h_{i,j,\sigma _u, \sigma _v}(\varvec{x}) = d \cdot \tfrac{\delta ^2}{4} \,, \end{aligned}$$

and which contains the origin (this can easily be confirmed by elementary calculations, using the relations \(u^2+v^2=1\) and \(d \cdot u^2-v^2=d \cdot \tfrac{\delta ^2}{4}\)). The function h is quadratic and strictly convex with gradient \(h'(\varvec{x})=8 \cdot \varvec{x}\) and Hessian \(h''(\varvec{x})=8\cdot \varvec{I}\) with \(d \times d\)-identity matrix \(\varvec{I}\), and hence h has the origin as unique minimizer with \(h(\varvec{0})=d^2 \cdot (d-1)\cdot \delta ^2\).

Now assume that there exists a center \(\tilde{\varvec{x}} \in {\mathbb {R}}^d\) with \({{\,\mathrm{rd}\,}}_S(\tilde{\varvec{x}})< \delta \). This implies that \(|\Vert \varvec{p_{i,j,\sigma _u, \sigma _v, -1}}-\tilde{\varvec{x}}\Vert _{}-\Vert \varvec{p_{i,j,\sigma _u, \sigma _v, +1}}-\tilde{\varvec{x}}\Vert _{}|< \delta \), which is equivalent to \(h_{i,j,\sigma _u, \sigma _v}(\tilde{\varvec{x}}) < d \cdot \tfrac{\delta ^2}{4}\). Summing up all \(4\cdot d \cdot (d-1)\) inequalities yields the inequality \(h(\tilde{\varvec{x}}) < d^2 \cdot (d-1)\cdot \delta ^2=h(\varvec{0})\), which is a contradition to the origin beeing the unique minimizer of h. \(\square \)

We create examples according to Lemma 4.1 for different dimensions d and numbers n of data points (where random points are added to the points \(\varvec{p_{i,j,\sigma _u, \sigma _v, \sigma _w}}\)). In all cases, our algorithm is run with accuracy \(\frac{{{\,\mathrm{rd}\,}}(S)}{100}\), initial center \(\varvec{x_0}=(0.3,\ldots ,0.3)\), and initial edge length \(h_0=1\). The initial roundness is \({{\,\mathrm{rd}\,}}_S(\varvec{x_0}) \approx 1\). For each value of d, n and \({{\,\mathrm{rd}\,}}(S)\) we compute the average number of roundness function evaluations, taken over 100 runs with different random points.

At first, we examine the behavior for increasing dimension with a fixed number \(n=1000\) of data points. The results are displayed in Table 2. The values in brackets are obtained with the simple test of Algorithm 1 only (which only worked well for \(d\le 5\)), and the values without brackets are obtained together with Algorithm 2. This confirms that Algorithm 2 greatly accelerates the basic global search. On our PC (Intel(R) Core(TM) i7 CPU with 2.8 GHz and 12 GB RAM) the average computing times with Algorithm 2 are negligible for \(d=2,3\), less than 0.1 s for \(d=4,5\), about 3 s for \(d=7\), about 20 s for \(d=8\), and about 3 minutes for \(d=9\). In the following we always use Algorithm 1 together with Algorithm 2.

Table 2 Average number of roundness function evaluations for examples created according to Lemma 4.1
Full size table

Next we examine the behavior for an increasing number of data points with fixed dimension and roundness \({{\,\mathrm{rd}\,}}(S)= 0.01\). As can be seen in Fig. 6, the average number of roundness function evaluations is almost constant.

Fig. 6
figure 6

Average number of roundness function evaluations for \({{\,\mathrm{rd}\,}}(S)= 0.01\) and an increasing number of data points

Full size image

Figure 7 confirms that the average computing time increases at most linearly with the number of data points, and that the slope increases with the dimension.

Fig. 7
figure 7

Average computing time in seconds for \({{\,\mathrm{rd}\,}}(S)= 0.01\) and an increasing number of data points

Full size image

4.2 Publicly Available Test Data

We check the performance of our algorithm for the publicly available 2D datasets with CMM data (cartesian coordinates obtained with a Coordinate Measuring Machine) that were used in Table 5 of [19], where a comparative study of existing roundness evaluation algorithms was done. The number of data points is rather small, ranging from \(n=16\) to \(n=32\). The roundness values obtained in [19] range from 0.03 [mm] to 36 [mm], and are given with 4 decimal digits. Therefore we run our algorithm with accuracy \(\varepsilon =10^{-6}\). We compute an initial center \(\varvec{x_0}\) by solving the linear least squares problem (LSQ), and choose the initial edge length \(h_0\) in two ways: Once as \(h_0=r_0\) with the mean radius \(r_0\) at the initial center, and once as \(h_0=2 \rho \) with the radius \(\rho \) according to Theorem 2.1. In all cases, we obtain the same minimal roundness values as given in [19] and computing times are negligible. Table 3 lists the initial roundness value \({{\,\mathrm{rd}\,}}_S(\varvec{x_0})\), the computed roundness value \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})\), the ratio \({\hat{\eta }}\) at the approximate minimizer, the initial edge length \(h_0\), and the number of roundness function evaluations, where the values without brackets correspond to the choice \(h_0=2 \rho \), and the values in brackets correspond to the choice \(h_0=r_0\). With the smaller initial edge length less function evaluations are needed, but, as expected, the difference is not that great.

Table 3 Results for the publicly available 2D-datasets used in Table 5 of [19]
Full size table

4.3 Polygonal Chains and Triangulations

Here we consider examples for the case that S is not a set of finitely many points, but a union of simple sets like line segments or triangles. The roundness function for such sets can be evaluated easily: At first, we recall that for any compact polyhedral set P the maximum distance \(\max _{\varvec{p} \in P}\Vert \varvec{x}-\varvec{p}\Vert _{}\) of a given center \(\varvec{x}\) to any point of P is attained at one of the vertices of P. Consequently this also holds for lines and triangles.

Let \(L:=\{\varvec{v}+ t \cdot \varvec{u}\, :\, t \in [0,1]\}\) be a line segment with \(\varvec{v}, \varvec{u} \in {\mathbb {R}}^d\). Then the minimum distance of \(\varvec{x}\) to any point of L can be computed as

$$\begin{aligned} \min _{\varvec{p} \in L}\Vert \varvec{x}-\varvec{p}\Vert _{}=\min _{t \in [0,1]}\Vert \varvec{x}-(\varvec{v}+ t \cdot \varvec{u})\Vert _{} \,, \end{aligned}$$

and is attained at \(\varvec{v}+ {\hat{t}} \cdot \varvec{u}\) with \({\hat{t}}=\min \big \{1,\max \{0,\tfrac{\langle \varvec{x}-\varvec{v} \,,\, \varvec{u}\rangle }{\Vert \varvec{u}\Vert _{}^2}\}\big \}\).

Now let \(D:=\{\varvec{v}+ t_1 \cdot \varvec{u_1}+ t_2 \cdot \varvec{u_2}\, :\, t_1,t_2 \in [0,1]\,,\, t_1+t_2 \le 1\}\) be a triangle with \(\varvec{v}, \varvec{u_1}, \varvec{u_2} \in {\mathbb {R}}^d\). Then the minimum distance of \(\varvec{x}\) to any point of D can be computed as

$$\begin{aligned} m:=\min _{\varvec{p} \in D}\Vert \varvec{x}-\varvec{p}\Vert _{}=\min _{t_1,t_2 \in [0,1]\,,\, t_1+t_2 \le 1}\Vert \varvec{x}-(\varvec{v}+t_1 \cdot \varvec{u_1}+ t_2 \cdot \varvec{u_2})\Vert _{}\,. \end{aligned}$$

Let \(({\hat{t}}_1,{\hat{t}}_2)\) be the solution of the unconstrained linear least squares problem \(\min _{t_1,t_2 \in {\mathbb {R}}}\Vert \varvec{x}-(\varvec{v}+t_1 \cdot \varvec{u_1}+ t_2 \cdot \varvec{u_2})\Vert _{}^2\). If it fulfills \({\hat{t}}_1,{\hat{t}}_2 \in [0,1]\) and \({\hat{t}}_1+{\hat{t}}_2 \le 1\), then the minimum distance m is attained at \(\varvec{{\hat{p}}}:=\varvec{v}+{\hat{t}}_1 \cdot \varvec{u_1}+ {\hat{t}}_2 \cdot \varvec{u_2}\); otherwise, it is attained at one of the edges of D (which are line segments as above).

This also shows that the roundness function can be evaluated in time \({\mathcal {O}}(n)\), where n is the number of vertices of S.

In the first example, the vertices \(\varvec{v_1},\ldots ,\varvec{v_6}\) of the polygonal chain S in Fig. 8 lie exactly on the circle with radius 1 around the origin, and have coordinates \(\varvec{v_i}=\big ( \cos (\phi _i), \sin (\phi _i)\big )\) with \(\phi _i \in \{0,\frac{\pi }{8},\frac{\pi }{4},\frac{\pi }{2},\frac{5\pi }{8},\frac{13\pi }{16}\}\) for \(i=1,\ldots ,6\). We run our algorithm with initial center \(\varvec{x_0}=(0.3,0.3)\) and initial edge length \(h_0=r_0=0.8459\). The initial roundness is \({{\,\mathrm{rd}\,}}_S(\varvec{x_0})=0.6281\). With accuracy \(\varepsilon =10^{-6}\) we obtain the approximate minimizer \(\varvec{{\hat{x}}}=(-0.0239,-0.0578)\) with roundness \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})=0.0716\) after 168 roundness function evaluations in about 0.02 s.

Fig. 8
figure 8

Blue: The vertices of the polygonal chain S lie exactly on a circle. Red: Computed optimal center \(\varvec{{\hat{x}}}\)

Full size image

In the second example, the 63 vertices \(\varvec{v_{i,j}}\) of the triangulation S shown in Fig. 9 lie exactly on the sphere with radius 1 around the origin, and have coordinates \(\varvec{v_{i,j}}=\big ( \cos (\phi _i) \cdot \sin (\theta _j), \sin (\phi _i)\cdot \sin (\theta _j),\cos (\theta _j)\big )\) with \(\phi _i =i\cdot \frac{\pi }{24}\) for \(i=0,\ldots ,8\), and \(\theta _j =\frac{\pi }{4}+ j \cdot \frac{\pi }{24}\) for \(j=0,\ldots ,6\). We run our algorithm with initial center \(\varvec{x_0}=(0.3,0.3,0.3)\) and initial edge length \(h_0=r_0=0.6482\). The initial roundness is \({{\,\mathrm{rd}\,}}_S(\varvec{x_0})=0.3406\). With accuracy \(\varepsilon =10^{-6}\) we obtain the origin as approximate minimizer \(\varvec{{\hat{x}}}=\varvec{0}\) (rounded to 4 decimal digits) with roundness \({{\,\mathrm{rd}\,}}_S(\varvec{{\hat{x}}})=0.0043\) after 698 roundness function evaluations in about 0.4 s.

Fig. 9
figure 9

The vertices of the triangulation S lie exactly on a sphere

Full size image

These examples further demonstrate that the algorithm also copes well in cases where data are given only on a part of a sphere (circle), so that Assumption (A1) is not fulfilled.

5 Conclusions

We have analyzed a global search algorithm for the solution of the roundness problem. In each iteration, it uses two tests to reduce the search region. For almost round sets the first test leads to relatively small upper bounds on the number of roundness function evaluations that are needed to find an approximate global minimizer, whose roundness function value is not larger than twice the minimal value. The second test is shown to improve the performance in the final iterations to reach a better accuracy. Since our numerical experiments indicate that the improvement in practice is even much better then predicted by the theoretical analysis, it would be interesting to know whether the theoretical result can be strengthened. Furthermore, similar strategies might be used to formulate and analyze global search algorithms for related problems like cylindricity or flatness problems under appropriate assumptions on the data.