Computation of the Hausdorff Distance between Two Compact Convex Sets

Abstract

The Hausdorff distance between two closed sets has important theoretical and practical applications. Yet apart from finite point clouds, there appear to be no generic algorithms for computing this quantity. Because many infinite sets are defined by algebraic equalities and inequalities, this a huge gap. The current paper constructs Frank–Wolfe and projected gradient ascent algorithms for computing the Hausdorff distance between two compact convex sets. Although these algorithms are guaranteed to go uphill, they can become trapped by local maxima. To avoid this defect, we investigate a homotopy method that gradually deforms two balls into the two target sets. The Frank–Wolfe and projected gradient algorithms are tested on two pairs $(A,B)$ of compact convex sets, where: (1) A is the box $[-\mathbf{1},\mathbf{1}]$ translated by $\mathbf{1}$ and B is the intersection of the unit ball and the non-negative orthant; and (2) A is the probability simplex and B is the ${\ell }_1$ unit ball translated by $\mathbf{1}$. For problem (2), we find the Hausdorff distance analytically. Projected gradient ascent is more reliable than the Frank–Wolfe algorithm and finds the exact solution of problem (2). Homotopy improves the performance of both algorithms when the exact solution is unknown or unattained.

1. Introduction

The Hausdorff distance [1,2] between two compact sets A and B in a Euclidean space ${\mathbb{R}}^p$ is defined as

$$\begin{array}{ccc}\hfill d_H(A,B)& =& max\left\{d\right(A,B),d(B,A\left)\right\},\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill d(A,B)& =& \underset{\mathit{x}\in A}{max}\mathrm{dist}(\mathit{x},B),d(B,A)=\underset{\mathit{x}\in B}{max}\mathrm{dist}(\mathit{x},A),\hfill \end{array}$$

$\mathrm{dist}(\mathit{x},A)={min}_{\mathit{y}\in A}\parallel \mathit{x}-\mathit{y}\parallel$, and $\mathrm{dist}(\mathit{x},B)={min}_{\mathit{y}\in B}\parallel \mathit{x}-\mathit{y}\parallel$. Here, $\parallel ·\parallel$ denotes the standard Euclidean norm in ${\mathbb{R}}^p$. The Blaschke formula

$$\begin{array}{ccc}\hfill d_H(A,B)& =& \underset{\mathit{x}}{max}|\mathrm{dist}(\mathit{x},A)-\mathrm{dist}(\mathit{x},B)|\hfill \end{array}$$

serves as an alternative definition of Hausdorff distance [3]. Wikipedia has a helpful entry on Hausdorff distance with a two-dimensional illustration. The theoretical value of Hausdorff distance stems from the fact that it turns the collection of compact sets into a complete separable metric space. In general, Hausdorff distance is challenging to compute.

Hausdorff distance has many applications. For instance, it is instrumental in defining continuity, compactness, and completeness for integral operators, differential operators, and Fourier transforms in functional analysis [4,5]. These concepts are in turn relevant to the analysis of the existence, uniqueness, and stability of solutions to various equations in mathematics and physics [6]. In computer vision, Hausdorff distance enables object recognition [7] and allows one to quantify the difference between two different representations of the same object [8]. Edge detection and pixelization are usually necessary preprocessing steps. Other applications include robotics [9], the fractal modeling of biological structures [10], and the numerical computation of attractors in dynamical systems [11].

The current paper derives and tests two new algorithms for computing the distance $d_H(A,B)$ between two compact convex sets. The previous work on this intrinsically interesting problem is mostly limited to finite point clouds, usually in two and three dimensions [12,13,14,15,16]. The formulas

$$\begin{array}{ccc}\hfill d(A,B)& =& \underset{\mathit{x}\in A}{max}\underset{\mathit{y}\in B}{min}\parallel \mathit{x}-\mathit{y}\parallel \mathrm{and}d(B,A)=\underset{\mathit{x}\in B}{max}\underset{\mathit{y}\in A}{min}\parallel \mathit{x}-\mathit{y}\parallel \hfill \end{array}$$

can be naively implemented for finite sets A and B. The naive implementation benefits from fast software such as the Julia Distances package, which exploits matrix multiplication to find all Euclidean distances between the column vectors of two matrices. The ImageDistances.jl package (github.com/JuliaImages/Images.jl) appears to rely on the naive method [12] for computing Hausdorff distances. Once the distances $d_{ij}=\parallel {\mathit{x}}_i-{\mathit{y}}_j\parallel$ are computed, the computational complexity of the naive method is $O\left(mn\right)$, where m and n equal the number of points ${\mathit{x}}_i\in A$ and ${\mathit{y}}_j\in B$, respectively. This complexity can be reduced by various devices, as suggested in the cited references.

The algorithms for computing the Hausdorff distance between two polygons [17], two curves [18] in the plane, and a curve and a surface [19] represent exceptions to the discrete point-cloud method. More complicated sets defined by algebraic formulas can be attacked by pixelating the sets or sampling them at a dense set of random points. The methods of continuous optimization offer an attractive alternative to the various current methods. Although the calculation of $d_H(A,B)$ takes us outside the comfortable realm of convex optimization, the tools of convex calculus are highly pertinent. To their credit, these tools perform well in higher dimensions. It remains to be seen whether the Hausdorff distance will have practical value in shape recognition in this regime. It would be prudent to keep the possibility in mind.

To calculate $d_H(A,B)$, it clearly suffices to calculate $d(A,B)$ and $d(B,A)$ separately and take the maximum. For many sets B, $\mathrm{dist}(\mathit{x},B)$ is explicitly known or can be computed by an efficient algorithm [20,21]. The Euclidean distance $\mathrm{dist}(\mathit{x},B)$ can be expressed as

$$\begin{array}{ccc}\hfill \mathrm{dist}(\mathit{x},B)& =& \parallel \mathit{x}-P_B\left(\mathit{x}\right)\parallel ,\hfill \end{array}$$

where $P_B\left(\mathit{x}\right)$ is the projection of $\mathit{x}$ onto B. When B is convex, the projection operator $P_B\left(\mathit{x}\right)$ is single-valued. For closed nonconvex sets, the projection is multiple-valued for some $\mathit{x}$, but these points are very rare; indeed, they are of Lebesgue measure 0 [22]. The scaled squared distance $\frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2$ function is smoother than $\mathrm{dist}(\mathit{x},B)$. One can show that the former function is differentiable with gradient

$$\begin{array}{ccc}\hfill \nabla \frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2& =& \mathit{x}-P_B\left(\mathit{x}\right)\hfill \end{array}$$

(1)

at all points $\mathit{x}$ where $P_B\left(\mathit{x}\right)$ is single-valued [23].

The support function ${\sigma }_A\left(\mathit{v}\right)={max}_{\mathit{x}\in A}{\mathit{v}}^{\top }\mathit{x}$ and the corresponding support set ${\mathrm{supp}}_A\left(\mathit{v}\right)={\mathrm{argmax}}_{\mathit{x}\in A}{\mathit{v}}^{\top }\mathit{x}$ also play key roles in our algorithm development. The maximum of $d(A,B)$ exists and necessarily occurs on the boundary of A. In fact, Bauer’s maximum principle [24,25] implies that the maximum is achieved at an extreme point $\mathit{x}$ of A. The point of B corresponding to $\mathit{x}$ occurs on the boundary of B, but not necessarily at an extreme point of B. The supporting set ${\mathrm{supp}}_A\left(\mathit{v}\right)$ is a singleton if and only if ${\sigma }_A\left(\mathit{v}\right)$ is differentiable at $\mathit{v}$.

Our first algorithm for computing $d(A,B)$,

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& \in & {\mathrm{supp}}_A\left({\mathit{v}}_n\right)={\mathrm{supp}}_A[{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]\hfill \end{array}$$

(2)

is a Frank–Wolfe algorithm [26,27]. Our second algorithm,

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& =& P_A({\mathit{x}}_n+{\mathit{v}}_n)=P_A[2{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]\hfill \end{array}$$

(3)

is a projected gradient algorithm [21,28]. Both algorithms force the objective function $\frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2$ uphill and are iterated until convergence. Because the algorithms can become trapped by local maxima, they are not infallible in finding a global maximum. To overcome this tendency, we introduce a homotopy method that gradually transitions the calculation of the Hausdorff distance from the simple case of two balls to the actual problem of finding $d(A,B)$. Homotopy is one of several heuristics for maximizing multi-modal functions [29]. The crucial difference between projected gradient ascent and the Frank–Wolfe algorithm is that one depends on projection while the other depends on both projection and supporting sets. This difference obviously favors projected gradient ascent.

As a roadmap to the rest of this paper, Section 2 presents (a) the basic notation, (b) a brief overview of the minorization–maximization principle that stands behind the new iterative algorithms (2) and (3), (c) a summary of the support functions and supporting sets, (d) the derivation of both algorithms, (e) a description of the homotopy method, and (f) an explanation of relevant convergence theory. Section 3 tests our two iterative algorithms and the point-cloud method on two representative problems. Both iterative algorithms are orders of magnitude faster than the point-cloud method and benefit from the homotopy heuristic. Projected gradient ascent is also more accurate than the point-cloud method. Section 4 summarizes our conclusions, mentions limitations, and suggests new avenues for research. Appendix A, proves some of the mathematical assertions made in the text and provides the full Julia code for our numerical examples. Note that the code is organized from bottom to top, with the main program occurring at the bottom.

2. Materials and Methods

As a prelude to the derivation of the two algorithms, it would be helpful to clarify our notation and make a few remarks about MM algorithms, support functions, and supporting sets. For projected gradient ascent and its homotopy modification, we provide algorithm flowcharts.

2.1. Notation

Here are the notational conventions used throughout this article. All vectors appear in boldface. All entries of the vectors $\mathit{0}$ and $\mathbf{1}$ equal 0 or 1, respectively. The ${}^{\top }$ superscript indicates a vector transpose. The Euclidean norm of a vector $\mathit{x}$ is denoted by $\parallel \mathit{x}\parallel$. For a smooth real-valued function $f\left(\mathit{x}\right)$, we write its gradient (column vector of partial derivatives) as $\nabla f\left(\mathit{x}\right)$ and its first differential (row vector of partial derivatives) as $df\left(\mathit{x}\right)=\nabla f{\left(\mathit{x}\right)}^{\top }$. Finally, we denote the directional derivative of $f\left(\mathit{x}\right)$ in the direction $\mathit{v}$ by $d_{\mathit{v}}f\left(\mathit{x}\right)$. When $f\left(\mathit{x}\right)$ is differentiable, $d_{\mathit{v}}f\left(\mathit{x}\right)=df\left(\mathit{x}\right)\mathit{v}$.

2.2. MM Algorithms

The algorithms explored here are minorization–maximization (MM) algorithms [23,30]. They depend on surrogate functions $g(\mathit{x}\mid {\mathit{x}}_n)$ that minorize the original objective $f\left(\mathit{x}\right)$ around the current iterate ${\mathit{x}}_n$ in the sense of satisfying the tangency condition $g({\mathit{x}}_n\mid {\mathit{x}}_n)=f\left({\mathit{x}}_n\right)$ and the domination condition $g(\mathit{x}\mid {\mathit{x}}_n)\le f\left(\mathit{x}\right)$ for all $\mathit{x}$. The surrogate balances the two goals of hugging the objective tightly and simplifying maximization. Maximizing the surrogate produces the next iterate ${\mathit{x}}_{n+1}$ and drives the objective uphill because

$$\begin{array}{ccc}\hfill f\left({\mathit{x}}_{n+1}\right)& \ge & g({\mathit{x}}_{n+1}\mid {\mathit{x}}_n)\ge g({\mathit{x}}_n\mid {\mathit{x}}_n)=f\left({\mathit{x}}_n\right).\hfill \end{array}$$

In minimization, the surrogate majorizes the objective and is instead minimized. The tangency condition remains the same, but the domination condition $g(\mathit{x}\mid {\mathit{x}}_n)\ge f\left(\mathit{x}\right)$ is now reversed. The celebrated EM (expectation–maximization) principle for maximum likelihood estimation with missing data [31] is a special case of minorization–maximization. In the EM setting, Jensen’s inequality supplies the surrogate as the expectation of the complete data log-likelihood conditional on the observed data.

2.3. Support Functions and Supporting Sets

The set of supporting points ${\mathrm{supp}}_S\left(\mathit{v}\right)={\mathrm{argmax}}_{\mathit{x}\in S}{\mathit{v}}^{\top }\mathit{x}$ determines the support function ${\sigma }_S\left(\mathit{v}\right)$. For instance, the ${\ell }_1$ unit ball has ${\mathrm{supp}}_S\left(\mathit{v}\right)$ equal to the convex hull of the vertices $\pm {\mathit{e}}_i$ where $|v_i|$ is largest. For the unit simplex, ${\mathrm{supp}}_S\left(\mathit{v}\right)$ equals the convex hull of the vertices ${\mathit{e}}_i$ where $v_i$ is largest. For a Minkowski sum $A+B$, ${\mathrm{supp}}_{A+B}\left(\mathit{v}\right)={\mathrm{supp}}_A\left(\mathit{v}\right)+{\mathrm{supp}}_B\left(\mathit{v}\right)$. If S is either a convex cone or a compact convex set that is symmetric about the origin with a non-empty interior, then its support function ${\sigma }_S\left(\mathit{y}\right)$ has a special form. In the former case, ${\sigma }_S\left(\mathit{y}\right)$ is the indicator of the dual cone, and in the latter case, ${\sigma }_S\left(\mathit{y}\right)$ is a norm. The support function of a Cartesian product is the Cartesian product of the separate support functions. For instance, the support function of rectangle $[\mathit{a},\mathit{b}]$ reduces to the one-dimensional case, where ${\mathrm{supp}}_{[a,b]}\left(\mathit{v}\right)$ is a when $v_i<0$, b when $v_i>0$, and all of $[a,b]$ when $v_i=0$. There are many other known support functions. For instance, one-sided penalties such as $cmax\{y,0\}$ and asymmetric penalties such as ${\sigma }_{[-a,b]}\left(y\right)$ are covered by the current theory. Indeed, the former is the support function generated by the interval $[0,c]$. The latter is the tilted absolute value equal to $by$ for $y\ge 0$ and to $-ay$ for $y<0$. The support function of a singleton $\left\{\mathit{a}\right\}$ is the linear function ${\mathit{a}}^{\top }\mathit{y}$. More generally, the support function of the convex hull of the set $\{{\mathit{a}}_1,\dots ,{\mathit{a}}_d\}$ is ${max}_{1\le i\le d}{\mathit{a}}_i^{\top }\mathit{y}$. The support function of the line segment from $-\mathit{a}$ to $\mathit{a}$ is the absolute value $|{\mathit{a}}^{\top }\mathit{y}|$. Adding a constant vector $\mathit{a}$ to a set S produces the support function ${\sigma }_S\left(\mathit{y}\right)+{\mathit{a}}^{\top }\mathit{y}$. It is trivial to project onto $S+\mathit{a}$ if one can project onto S. For any non-negative scalar c, the set $cS$ has support function $c{\sigma }_S\left(\mathit{y}\right)$. Again, it is trivial to project onto $cS$ if one can project onto S.

2.4. Derivation of the Algorithms

When B is convex, the supporting hyperplane inequality

$$\begin{array}{ccc}\hfill \frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2& \ge & \frac{1}{2}\mathrm{dist}{({\mathit{x}}_n,B)}^2+{\mathit{v}}_n^{\top }(\mathit{x}-{\mathit{x}}_n)\hfill \end{array}$$

for ${\mathit{v}}_n={\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)$ generates our first algorithm. Maximizing this minorization over $\mathit{x}\in A$ is equivalent to calculating the support function ${\sigma }_A\left({\mathit{v}}_n\right)={sup}_{\mathit{x}\in A}{\mathit{v}}_n^{\top }\mathit{x}$. If ${\mathrm{supp}}_A\left(\mathit{v}\right)$ denotes the set of points in A where ${\sigma }_A\left(\mathit{v}\right)$ is attained, then the Frank–Wolfe algorithm just described can be phrased as

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& \in & {\mathrm{supp}}_A\left({\mathit{v}}_n\right)={\mathrm{supp}}_A[{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)].\hfill \end{array}$$

The MM principle guarantees that the next iterate ${\mathit{x}}_{n+1}$ will tend to increase the objective $\frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2$ unless ${\mathit{v}}_n=\mathit{0}$. This exception occurs when ${\mathit{x}}_n\in B$. Fortunately, when the iterates begin in $A\setminus B$, they remain in $A\setminus B$. Indeed, if ${\mathit{x}}_n\in A\setminus B$ but ${\mathit{x}}_{n+1}\in B$, then the obtuse angle condition [23] requires

$$\begin{array}{ccc}\hfill {[{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]}^{\top }{\mathit{x}}_{n+1}& \le & {[{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]}^{\top }{P}_B\left({\mathit{x}}_n\right)\hfill \\ & <& {[{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]}^{\top }{\mathit{x}}_n,\hfill \end{array}$$

contradicting the optimality of ${\mathit{x}}_{n+1}$. To achieve the requirement ${\mathit{x}}_0\in A\setminus B$ of the Frank–Wolfe method, we put ${\mathit{x}}_0=P_A\left(\mathit{r}\right)$ for a random vector $\mathit{r}$ and then check that $\parallel P_B\left({\mathit{x}}_0\right)-\mathit{x}\parallel >0$.

If B is closed and convex, then the gradient of the function $\frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2$ is Lipschitz with constant 1 [23]. This fact plus the outcome of completing the square entails the minorization

$$\begin{array}{ccc}\hfill \frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2& \ge & \frac{1}{2}\mathrm{dist}{({\mathit{x}}_n,B)}^2+{\mathit{v}}_n^{\top }(\mathit{x}-{\mathit{x}}_n)-\frac{1}{2}{\parallel \mathit{x}-{\mathit{x}}_n\parallel }^2\hfill \\ & =& \frac{1}{2}\mathrm{dist}{({\mathit{x}}_n,B)}^2-\frac{1}{2}\parallel \mathit{x}-{\mathit{x}}_n-{\mathit{v}}_n{\parallel }^2+\frac{1}{2}{\parallel {\mathit{v}}_n\parallel }^2.\hfill \end{array}$$

Hence, the MM principle implies that defining

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& =& P_A({\mathit{x}}_n+{\mathit{v}}_n)=P_A[2{\mathit{x}}_n-P_B\left({\mathit{x}}_n\right)]\hfill \end{array}$$

also increases $\frac{1}{2}\mathrm{dist}{(\mathit{x},B)}^2$. This second of our two algorithms is a special case of projected gradient ascent. Its flowchart (Algorithm 1) summarizes this straightforward strategy started from many random points ${\mathit{x}}_0$.

Algorithm 1 Computation of $d(A,B)$ by Projected Gradient Ascent

Require: Projection operators $P_A\left(\mathit{x}\right)$ and $P_B\left(\mathit{y}\right)$, initial point ${\mathit{x}}_0$, maximum iterations n, and convergence criterion $ϵ>0$. 1: $\mathit{x}={\mathit{x}}_0$ 2: for $iter=1:n$  do 3:     ${\mathit{x}}_{\mathrm{new}}=P_A[2\mathit{x}-P_B\left(\mathit{x}\right)]$ 4:     if $\parallel {\mathit{x}}_{\mathrm{new}}-\mathit{x}\parallel <ϵ$ then 5:         break 6:     else 7:         $\mathit{x}={\mathit{x}}_{\mathrm{new}}$ 8:     end if 9: end for 10: Return ${\mathit{x}}_{\mathrm{new}}$

2.5. A Homotopy Method

Although both algorithms are guaranteed to increase the objective, they both suffer from the danger of being trapped by local maxima. One obvious remedy is to launch the algorithms from different random points. A more systematic alternative is to exploit homotopy. The idea is to gradually deform both sets A and B from the unit ball U at the origin, where $d_H(U,U)=0$ is known, into the target sets A and B. In practice, we follow the solution path along the family of set pairs $[tA+(1-t)U,tB+(1-t\left)U\right]$ from $t=0$ to $t=1$. This strategy is viable for projected gradient ascent because we can project points onto the Minkowski convex combination $tC+(1-t)D$ by three devices. First, it is well known that when A and B are balls with radii $r_A$ and $r_B$ and centers ${\mathit{c}}_A$ and ${\mathit{c}}_B$, respectively, the distance $\mathrm{dist}(\mathit{x},B)$ is maximized by taking $\mathit{x}={\mathit{c}}_A-r_A\frac{({\mathit{c}}_B-{\mathit{c}}_A)}{\parallel {\mathit{c}}_B-{\mathit{c}}_A\parallel }$, unless $A\subset B$, in which case the maximum is 0 [32]. For the convenience of the reader, Proposition A1 of Appendix A proves this assertion. Second, one can exploit the projection identity $P_{tS}\left(\mathit{z}\right)=t{P}_S\left(t^{-1}\mathit{z}\right)$ for any $t>0$. Third, there is an effective algorithm for projecting onto a Minkowski sum $C+D$ [33]. The idea is to alternate the minimization of $\parallel \mathit{z}-\mathit{c}-\mathit{d}\parallel$ with respect to $\mathit{c}\in C$ and $\mathit{d}\in D$. The iteration scheme ${\mathit{c}}_{n+1}=P_C(\mathit{z}-{\mathit{d}}_n)$ and ${\mathit{d}}_{n+1}=P_D(\mathit{z}-{\mathit{c}}_{n+1})$ is guaranteed to converge at a linear rate when either set is strongly convex. Recall that a convex K is strongly convex if there exists an $r>0$ such that

$$\begin{array}{ccc}\hfill \alpha \mathit{x}+(1-\alpha )\mathit{y}+\frac{r}{2}\alpha (1-\alpha ){\parallel \mathit{x}-\mathit{y}\parallel }^2\mathit{z}& \in & K\hfill \end{array}$$

for all $\mathit{x}$ and $\mathit{y}$ in K, $\alpha \in [0,1]$, and unit vectors $\mathit{z}$ [34]. In particular, $(1-t)U$ is strongly convex when $t\in [0,1)$. Furthermore, the Cartesian product of two strongly convex sets is strongly convex [35]. The homotopy method succeeds because the early sets are more rounded and the objective generates fewer local maxima. The price for better performance is iterations within iterations and an overall slower algorithm. Algorithms 2 and 3 summarize our homotopy strategy for projected gradient ascent.

Algorithm 2 Minkowski Set Projection

Require: Projection operators $P_A\left(\mathit{x}\right)$ and $P_B\left(\mathit{y}\right)$, external point $\mathit{y}$, convexity constant c, maximum iterations n, and convergence criterion $ϵ>0$. 1: $\mathit{a}=\mathit{b}=\mathit{0}$ 2: for  $iter=1:n$  do 3:     ${\mathit{a}}_{\mathrm{new}}=c{P}_A[(\mathit{y}-\mathit{b})/c]$ 4:     ${\mathit{b}}_{\mathrm{new}}=(1-c)P_B[(\mathit{y}-{\mathit{a}}_{\mathrm{new}})/(1-c)]$ 5:     if $\parallel {\mathit{a}}_{\mathrm{new}}-\mathit{a}\parallel +\parallel {\mathit{b}}_{\mathrm{new}}-\mathit{b}\parallel <ϵ$ then 6:         break 7:     else 8:         $\mathit{a}={\mathit{a}}_{\mathrm{new}}$ 9:         $\mathit{b}={\mathit{b}}_{\mathrm{new}}$ 10:     end if 11: end for 12: Return ${\mathit{a}}_{\mathrm{new}}$ and ${\mathit{b}}_{\mathrm{new}}$

Algorithm 3 Homotopy Modification of Projected Gradient Ascent

Require: Projection operators $P_A\left(\mathit{x}\right)$ and $P_B\left(\mathit{x}\right)$, centers $c_A$ and $c_B$, and homotopy points h. 1: Set $d=\parallel c_A-c_B\parallel$ and $\mathit{x}={\mathit{c}}_A-\frac{({\mathit{c}}_B-{\mathit{c}}_A)}{d}$ ▹ distance between two balls of radius 1 2: Let $P_U$ be projection onto the unit ball 3: for $i=1:h-2$ do ▹ intermediate homotopy phases 4:     $c=\frac{i}{h-1}$ 5:     Put $P_{MA}=$ Minkowski sum projection for $P_A$ and $P_U$ and constant c 6:     Put $P_{MB}=$ Minkowski sum projection for $P_B$ and $P_U$ and constant c 7:     Perform projected gradient ascent with $P_{MA}$, $P_{MB}$, and initial point $\mathit{x}$ 8:     Let $\mathit{x}$ be the outcome 9: end for 10: Perform projected gradient ascent with $P_A$ and $P_B$ and initial point $\mathit{x}$ 11: Return converged value of $\mathit{x}$ ▹ output final phase of homotopy

For the Frank–Wolfe algorithm, similar homotopy tactics apply. For a Minkowski sum $C+D$, ${\mathrm{supp}}_{C+D}\left(\mathit{v}\right)={\mathrm{supp}}_C\left(\mathit{v}\right)+{\mathrm{supp}}_D\left(\mathit{v}\right)$. This fact plus the identity ${\mathrm{supp}}_{tC}\left(\mathit{v}\right)={\mathrm{supp}}_C\left(t\mathit{v}\right)$ for $t\ge 0$ makes it possible to carry out the homotopy method.

2.6. Convergence

Because this topic has been covered in previous studies [21,36,37,38,39], we give an abbreviated treatment here. Each algorithm is summarized by a closed algorithm map ${\mathit{x}}_{n+1}\in M\left({\mathit{x}}_n\right)$ that increases the objective $f\left(\mathit{x}\right)$. The limit points of the map occur among the stationary points of $f\left(\mathit{x}\right)$. By definition, a stationary point $\mathit{x}$ satisfies $d_{\mathit{v}}f\left(\mathit{x}\right)=df\left(\mathit{x}\right)\mathit{v}\le 0$ for all tangent vectors $\mathit{v}$ at $\mathit{x}$. The set of tangent vectors $\mathit{v}$ is the closure of the set of points $c(\mathit{y}-\mathit{x})$ with $\mathit{y}\in C$ and $c>0$. This is a place where the convexity of C comes into play. Hence, $\mathit{x}$ is a stationary point if and only if $df\left(\mathit{x}\right)\mathit{x}\ge df\left(\mathit{x}\right)\mathit{y}$ for all $\mathit{y}\in C$. With this distinction in mind, we state our basic theoretical findings for the Frank–Wolfe method. Homotopy is omitted in these considerations.

Proposition 1. 

The limit points of the Frank–Wolfe iterates (2) are stationary points of the objective $f\left(\mathit{x}\right)={min}_{\mathit{y}\in B}\parallel \mathit{x}-\mathit{y}\parallel$ on A. Furthermore, the bound

$$\begin{array}{ccc}\hfill \underset{0\le k\le n}{min}\underset{\mathit{y}\in A}{max}df\left({\mathit{x}}_k\right)(\mathit{y}-{\mathit{x}}_k)& \le & \frac{1}{n+1}[\underset{\mathit{x}\in A}{max}f\left(\mathit{x}\right)-f\left({\mathit{x}}_0\right)]\hfill \end{array}$$

holds. Thus, the stationary condition ${max}_{\mathit{y}\in A}df\left(\mathit{x}\right)(\mathit{y}-\mathit{x})\le 0$ is reasonable to expect at a limit point $\mathit{x}$ of the Frank–Wolfe algorithm.

Here is the corresponding finding for projected gradient ascent.

Proposition 2. 

The limit points of the projected gradient ascent iterates (3) are also stationary points of the objective $f\left(\mathit{x}\right)={min}_{\mathit{y}\in B}\parallel \mathit{x}-\mathit{y}\parallel$ on A. Furthermore, the bound

$$\begin{array}{ccc}\hfill \underset{0\le k\le n}{min}\parallel {\mathit{x}}_{k+1}-{\mathit{x}}_k\parallel & \le & \sqrt{\frac{2}{(n+1)}[f\left({\mathit{x}}_0\right)-{min}_{\mathit{x}\in A}f\left(\mathit{x}\right)]}\hfill \end{array}$$

holds.

Although the convergence rate $O\left(\frac{1}{\sqrt{n}}\right)$ of projected gradient ascent is slower than the corresponding slow convergence rate $O\left(\frac{1}{n}\right)$ of the Frank–Wolfe method, in practice, both algorithms usually converge in fewer than 100 iterations. In the case of the Frank–Wolfe algorithm, each iterate is an extreme point. Many convex sets possess only a finite number of extreme points, and convergence to one of them is guaranteed. Unfortunately, the converged point often provides just a local maximum.

3. Results

We tested the Frank–Wolfe and projected gradient ascent algorithms on two pairs $(A,B)$ of compact convex sets: (1) where A is the box $[-\mathbf{1},\mathbf{1}]$ translated by $\mathbf{1}$ and B is the intersection of the unit ball and the non-negative orthant; and (2) where A is the probability simplex and B is the ${\ell }_1$ unit ball translated by $\mathbf{1}$. These examples are representative, and for the second pair one can show that $d_H(A,B)=\sqrt{p}$, where p is the dimension of the ambient space. See Proposition A2 of Appendix A.

Table 1. Computation of $d_H(A,B)$ by various methods.

Set Pair	p	Method	Homotopy	Maximum	Mean	Std	Secs
(box, ball ∩ orthant)	2	proj grad	false	1.8284	1.3314	0.40789	0.00163
(box, ball ∩ orthant)	2	proj grad	true	1.8284	1.8284	0.0	0.00349
(box, ball ∩ orthant)	2	Frank-Wolfe	false	0.41421	0.16569	0.20394	0.000564
(box, ball ∩ orthant)	2	Frank-Wolfe	true	0.41421	0.41421	0.0	0.00101
(box, ball ∩ orthant)	2	point cloud	false	1.8159	1.8159	0.0	0.625
(simplex, L1 ball)	2	proj grad	false	1.4142	1.4142	0.0	0.001
(simplex, L1 ball)	2	proj grad	true	1.4142	1.4142	0.0	0.00272
(simplex, L1 ball)	2	Frank-Wolfe	false	1.4142	1.4142	0.0	0.000358
(simplex, L1 ball)	2	Frank-Wolfe	true	1.4142	1.4142	0.0	0.0574
(simplex, L1 ball)	2	point cloud	false	1.4142	1.4142	0.0	0.639
(box, ball ∩ orthant)	3	proj grad	false	2.4641	1.6358	0.50531	$7.34\times {10}^{-5}$
(box, ball ∩ orthant)	3	proj grad	true	2.4641	2.4641	0.0	0.000605
(box, ball ∩ orthant)	3	Frank-Wolfe	false	0.73205	0.31788	0.25266	$7.08\times {10}^{-5}$
(box, ball ∩ orthant)	3	Frank-Wolfe	true	0.73205	0.73205	0.0	0.000145
(box, ball ∩ orthant)	3	point cloud	false	2.4221	2.4221	0.0	0.63
(simplex, L1 ball)	3	proj grad	false	1.7321	1.7321	0.0	$5.86\times {10}^{-6}$
(simplex, L1 ball)	3	proj grad	true	1.7321	1.7321	0.0	0.000158
(simplex, L1 ball)	3	Frank-Wolfe	false	1.2247	1.2247	0.0	$4.44\times {10}^{-6}$
(simplex, L1 ball)	3	Frank-Wolfe	true	1.2247	1.2247	0.0	0.00186
(simplex, L1 ball)	3	point cloud	false	1.732	1.732	0.0	0.615
(box, ball ∩ orthant)	10	proj grad	false	5.0	3.4576	0.73165	0.0001
(box, ball ∩ orthant)	10	proj grad	true	5.3246	5.3246	0.0	0.000666
(box, ball ∩ orthant)	10	Frank-Wolfe	false	2.0	1.2288	0.36583	$9.71\times {10}^{-5}$
(box, ball ∩ orthant)	10	Frank-Wolfe	true	2.1623	2.1623	0.0	0.000183
(box, ball ∩ orthant)	10	point cloud	false	4.8158	4.8158	0.0	0.629
(simplex, L1 ball)	10	proj grad	false	3.1623	3.1623	0.0	$7.66\times {10}^{-6}$
(simplex, L1 ball)	10	proj grad	true	3.1623	3.1623	0.0	0.000749
(simplex, L1 ball)	10	Frank-Wolfe	false	2.6667	2.6667	0.0	$7.83\times {10}^{-6}$
(simplex, L1 ball)	10	Frank-Wolfe	true	2.6667	2.6667	0.0	0.00196
(simplex, L1 ball)	10	point cloud	false	3.1626	3.1626	0.0	0.613
(box, ball ∩ orthant)	100	proj grad	false	15.248	13.016	0.69754	0.000621
(box, ball ∩ orthant)	100	proj grad	true	19.0	19.0	0.0	0.00126
(box, ball ∩ orthant)	100	Frank-Wolfe	false	7.124	6.0078	0.34877	0.00168
(box, ball ∩ orthant)	100	Frank-Wolfe	true	9.0	9.0	0.0	0.000304
(box, ball ∩ orthant)	100	point cloud	false	15.598	15.598	0.0	0.629
(simplex, L1 ball)	100	proj grad	false	10.0	10.0	0.0	$2.77\times {10}^{-5}$
(simplex, L1 ball)	100	proj grad	true	10.0	10.0	0.0	0.000852
(simplex, L1 ball)	100	Frank-Wolfe	false	9.8494	9.8494	0.0	$2.21\times {10}^{-5}$
(simplex, L1 ball)	100	Frank-Wolfe	true	9.8494	9.8494	0.0	0.00299
(simplex, L1 ball)	100	point cloud	false	9.9493	9.9493	0.0	0.64
(box, ball ∩ orthant)	1000	proj grad	false	45.174	43.698	0.6652	0.00591
(box, ball ∩ orthant)	1000	proj grad	true	62.246	62.246	0.0	0.00875
(box, ball ∩ orthant)	1000	Frank-Wolfe	false	22.087	21.349	0.3326	0.00256
(box, ball ∩ orthant)	1000	Frank-Wolfe	true	30.623	30.623	0.0	0.00531
(box, ball ∩ orthant)	1000	point cloud	false	48.627	48.627	0.0	1.13
(simplex, L1 ball)	1000	proj grad	false	31.623	31.623	0.0	0.000377
(simplex, L1 ball)	1000	proj grad	true	31.623	31.623	0.0	0.0113
(simplex, L1 ball)	1000	Frank-Wolfe	false	31.575	31.575	0.0	0.000303
(simplex, L1 ball)	1000	Frank-Wolfe	true	31.575	31.575	0.0	0.013
(simplex, L1 ball)	1000	point cloud	false	31.597	31.597	0.0	1.11

presents our findings. The computation times are in seconds per trial across 100 random initializations and appear to scale affinely (a constant plus linear) in p. The columns Maximum, Mean, and Std convey summary statistics of the converged values of the Hausdorff distance. The point-cloud method generated ${10}^4$ random points to mimic each continuous set. Because the point-cloud method is non-iterative, a single run captured its performance. For the record, all computations were carried out on a MacBook Pro with a 2.3 GHz 8-core i9 chip and 16 GB of memory. Although the algorithms were embarrassingly parallel across trials, our Julia code is completely serial.

The random point-cloud method was not remotely competitive with projected gradient ascent in either accuracy or speed on these sample problems. It did produce approximate distances that confirmed the best results of the iterative methods. As measured by the quality of its solution, projected gradient ascent also outperformed the Frank–Wolfe algorithm. The Frank–Wolfe method was probably too aggressive, perhaps because it moved directly to an extreme point of A in computing $d(A,B)$. On the second problem, projected gradient ascent attained the global maximum across all trials. When the standard deviation of the converged values is positive, it follows that some trials were trapped by inferior local maxima. Both iterative algorithms benefit from the homotopy heuristic, which is fully deterministic. Accordingly, the standard deviations under homotopy equaled 0. Homotopy increased the computation times by less than an order of magnitude for 11 evenly spaced homotopy points. Finally, as p increased, both problems appeared easier to solve by the iterative methods. This behavior was particularly evident when $p=1000$. In contrast, the random point-cloud solutions deteriorated as p increased.

4. Discussion

The Hausdorff distance problem is intrinsically interesting, with theoretical applications throughout mathematics and practical applications in image processing. Given the non-convexity of the problem, it has not received nearly the attention in the mathematical literature as the closest point problem. Exact values of $d_H(A,B)$ are available in a few special cases such as the two highlighted in our appendix. Research on fast algorithms tends to be limited to random point clouds. Infinite sets defined by mathematical formulas have been largely ignored. The current paper partially rectifies this omission and demonstrates the value of continuous optimization techniques. The Frank–Wolfe and projected gradient ascent algorithms are relatively easy to code and extremely fast, even in high dimensions. Our preliminary experiments tilt toward projected gradient ascent as the more reliable of the two options. A naive implementation of the point-cloud method is not remotely competitive with projected gradient ascent. More exhaustive studies are warranted beyond the proof of principle examples presented here.

The standard convergence arguments covered in Section 2.6 guarantee that all limit points of the two algorithm classes are stationary points. In practice, convergence appears much faster than the slow rates mentioned in Propositions 1 and 2. We suspect, but have not proved, that full convergence to a stationary point always occurs. This exercise would require a foray into the difficult terrain of real algebraic geometry [40]. In any event, convergence to a global maximum is not guaranteed. Fortunately, safeguards can be put in place to improve the chances of successful convergence. The homotopy method capitalizes on the exact distance between two balls. Minkowski set rounding smooths the boundary of the target sets and steers iterates in a productive direction.

The computation of the Hausdorff distance $d_H(A,B)$ is apt to be much more challenging when either A or B is non-convex. Many sets can be represented as finite unions of compact convex sets. If $A={\cup }_i{A}_i$ and $B={\cup }_j{B}_j$, then the computation of $d_H(A,B)$ reduces to the computation of $d(A_i,B)$ for each index i and $d(B_j,A)$ for each index j. The identities $d(A,B)={max}_{i}d(A_i,B)$ and $d(B,A)={max}_{i}d(B_j,A)$ make this claim obvious. The further identity $\mathrm{dist}(\mathit{x},B)={min}_j\mathrm{dist}(\mathit{x},B_j)$ implies that

$$\begin{array}{ccc}\hfill d(A,B)& =& \underset{i}{max}d(A_i,B)=\underset{i}{max}\underset{\mathit{x}\in A_i}{max}\underset{j}{min}\mathrm{dist}(\mathit{x},B_j).\hfill \end{array}$$

In general, saddlepoint problems of this sort are hard to solve. One possible line of attack is to minorize ${min}_j\mathrm{dist}(\mathit{x},B_j)$ and then maximize the minorization.

The validation and implementation of this strategy will require substantial effort beyond the introduction to the problem pursued in the current paper. Let us merely add that the fast implementation of a more general Hausdorff distance algorithm will depend on the nature of the candidate sets $A_i$ and $B_j$. In the plane, triangles are appealing [41,42]. It is straightforward to project onto a triangle, and a triangle by definition possesses exactly three extreme points. Furthermore, a great deal of research under the heading of finite elements has identified good algorithms for triangulating complicated regions of the plane. The software Triangle for generating two-dimensional meshes and Delaunay triangulations is surely pertinent [43]. The triangularization of surfaces forms part of the MESH software for Hausdorff distance estimation [44].

We hope this paper will provoke a greater focus on the Hausdorff distance problem. As a prototype non-convex problem, it is worthy of far more attention. Continuous optimization tools can be brought to bear on the problem and may ultimately generate more efficient algorithms than the discrete algorithms designed for finite point clouds. The fact that our algorithms scale well in higher dimensions is a plus. We would also like to highlight the illumination that the MM principle brings to the construction of new high-dimensional optimization algorithms, including those considered here. Although neglected in the past, MM may well be the single most unifying principle of algorithm construction in continuous optimization.