# Computation of the Hausdorff Distance between Two Compact Convex Sets
## 1 Introduction
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\), \(\operatorname{dist}(\mathbf{x},B)\) is explicitly known or can be computed by an efficient algorithm [20; 21]. The Euclidean distance \(\operatorname{dist}(\mathbf{x},B)\) can be expressed as  
\[\operatorname{dist}(\mathbf{x},B) = \|\mathbf{x}-P_{B}(\mathbf{x})\|,\]  
where \(P_{B}(\mathbf{x})\) is the projection of \(\mathbf{x}\) onto \(B\). When \(B\) is convex, the projection operator \(P_{B}(\mathbf{x})\) is single-valued. For closed nonconvex sets, the projection is multiple-valued for some \(\mathbf{x}\), but these points are very rare; indeed, they are of Lebesgue measure 0 [22]. The scaled squared distance \(\frac{1}{2}\operatorname{dist}(\mathbf{x},B)^{2}\) function is smoother than \(\operatorname{dist}(\mathbf{x},B)\). One can show that the former function is differentiable with gradient  
\[\
abla\frac{1}{2}\operatorname{dist}(\mathbf{x},B)^{2} = \mathbf{x}-P_{B}(\mathbf{x}) \tag{1}\]  
at all points \(\mathbf{x}\) where \(P_{B}(\mathbf{x})\) is single-valued [23]. The support function \(\sigma_{A}(\mathbf{v})=\max_{\mathbf{x}\in A}\mathbf{v}^{\top}\mathbf{x}\) and the corresponding support set \(\operatorname{supp}_{A}(\mathbf{v})=\operatorname{argmax}_{\mathbf{x}\in A}\mathbf{v}^{ \top}\mathbf{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, [PERSON]'s maximum principle [24; 25] implies that the maximum is achieved at an extreme point \(\mathbf{x}\) of \(A\). The point of \(B\) corresponding to \(\mathbf{x}\) occurs on the boundary of \(B\), but not necessarily at an extreme point of \(B\). The supporting set \(\operatorname{supp}_{A}(\mathbf{v})\) is a singleton if and only if \(\sigma_{A}(\mathbf{v})\) is differentiable at \(\mathbf{v}\). Our first algorithm for computing \(d(A,B)\),  
\[\mathbf{x}_{n+1} \in \text{supp}_{A}(\mathbf{v}_{n}) = \text{supp}_{A}[\mathbf{x}_{n}-P_{B}(\mathbf{x}_{n})] \tag{2}\]  
is a Frank-Wolfe algorithm [26; 27]. Our second algorithm,  
\[\mathbf{x}_{n+1} = P_{A}(\mathbf{x}_{n}+\mathbf{v}_{n}) = P_{A}[2\mathbf{x}_{n}-P_{B}(\mathbf{x}_{n})] \tag{3}\]  
is a projected gradient algorithm [21; 28]. Both algorithms force the objective function \(\frac{1}{2}\operatorname{dist}(\mathbf{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.