Closest Farthest Widest

Abstract

The current paper proposes and tests algorithms for finding the diameter of a compact convex set and the farthest point in the set to another point. For these two nonconvex problems, I construct Frank–Wolfe and projected gradient ascent algorithms. Although these algorithms are guaranteed to go uphill, they can become trapped by local maxima. To avoid this defect, I investigate a homotopy method that gradually deforms a ball into the target set. Motivated by the Frank–Wolfe algorithm, I also find the support function of the intersection of a convex cone and a ball centered at the origin and elaborate a known bisection algorithm for calculating the support function of a convex sublevel set. The Frank–Wolfe and projected gradient algorithms are tested on five compact convex sets: (a) the box whose coordinates range between −1 and 1, (b) the intersection of the unit ball and the non-negative orthant, (c) the probability simplex, (d) the Manhattan-norm unit ball, and (e) a sublevel set of the elastic net penalty. Frank–Wolfe and projected gradient ascent are about equally fast on these test problems. Ignoring homotopy, the Frank–Wolfe algorithm is more reliable. However, homotopy allows projected gradient ascent to recover from its failures.

1. Introduction

Let S be a compact convex set and $\mathit{p}$ any external or internal point. This paper investigates algorithms for computing the two functions ${\mathrm{far}}_S\left(\mathit{p}\right)={max}_{\mathit{x}\in S}\parallel \mathit{x}-\mathit{p}\parallel$ and $\mathrm{diam}\left(S\right)={max}_{(\mathit{x},\mathit{y})\in S\times S}\parallel \mathit{x}-\mathit{y}\parallel$. The first is the farthest distance from $\mathit{p}$ to $\mathit{x}\in S$, and the second is the diameter of S. Both objective functions are convex and Lipschitz under the Euclidean norm, the first with constant 1 and the second with constant 2. Furthermore, if S is merely bounded, and $\mathrm{conv}\left(S\right)$ denotes its closed convex hull, and $\mathrm{ext}\left(S\right)$ denotes the extreme points of $\mathrm{conv}\left(S\right)$, then

$$\begin{array}{ccc}\hfill {\mathrm{far}}_S\left(\mathit{p}\right)& =& {\mathrm{far}}_{\mathrm{conv}\left(S\right)}\left(\mathit{p}\right) = {\mathrm{far}}_{\mathrm{ext}\left(S\right)}\left(\mathit{p}\right)\hfill \\ \hfill \mathrm{diam}\left(S\right)& =& \mathrm{diam}\left[\mathrm{conv}\right(S\left)\right] = \mathrm{diam}\left[\mathrm{ext}\right(S)\times \mathrm{ext}(S\left)\right].\hfill \end{array}$$

The fact that $\mathrm{conv}\left(S\right)$ is compact is crucial in the reaching these conclusions. Thus, without loss of generality, we can assume that S is a compact convex set. All boundary points S are extreme points when S is strictly convex [1,2].

If S is the convex hull of a finite cloud of points $\{{\mathit{x}}_1,\dots ,{\mathit{x}}_m\}$, then ${\mathrm{far}}_S\left(\mathit{p}\right)$ can be found by identifying the point ${\mathit{x}}_i$ that maximizes $\parallel {\mathit{x}}_i-\mathit{p}\parallel$. For a ball S and a point $\mathit{y}\notin S$, ${\mathrm{far}}_S\left(\mathit{p}\right)$ is found by projecting $\mathit{y}$ onto S and then extending the line through $\mathit{y}$ and its projection $\mathit{z}$ to the antipodal point of $\mathit{z}$. If the set S is a Cartesian product, then calculation of ${\mathrm{far}}_S\left(\mathit{p}\right)$ and $\mathrm{diam}\left(S\right)$ reduce to easy lower-dimensional problems. Many diameters are known. For a ball, ellipsoid $\{\mathit{x}:\frac{1}{2}{\mathit{x}}^{\top }\mathit{A}\mathit{x}\le r\}$, rectangle $[\mathit{a},\mathit{b}]$, probability simplex, and ${\mathsf{\ell }}_1$ unit ball, the diameter is, respectively, twice its radius, $\sqrt{8r{\lambda }_{min}^{-1}}$, $\parallel \mathit{a}-\mathit{b}\parallel$, $\sqrt{2}$, and 2, where ${\lambda }_{min}$ is the smallest eigenvalue of the positive definite matrix $\mathit{A}$ [3].

The distance function $\mathrm{dist}(\mathit{p},S)={min}_{\mathit{x}\in S}\parallel \mathit{x}-\mathit{p}\parallel$ is extremely well studied [4,5,6,7]. The unique point where the minimum is attained is the projection $P_S\left(\mathit{p}\right)$ of $\mathit{p}$ onto S. Many projection operators are known [4,6,7,8,9]. The web sites [10,11] provide Julia, Python, and Matlab implementations of the most commonly encountered projection operators. Projection onto a convex sublevel set $S=\{\mathit{x}:g(\mathit{x})\le 0\}$ is often more computationally demanding than other projection problems. For the special case when the proximal operator ${\mathrm{prox}}_{\lambda g}\left(\mathit{y}\right)={\mathrm{argmin}}_{\mathit{x}}\left[\lambda g\left(\mathit{x}\right)+\frac{1}{2}{\parallel \mathit{y}-\mathit{x}\parallel }^2\right]$ is known for all $\lambda >0$, one can solve the projection problem by bisection [12]. Bisection also turns out to be an attractive strategy for computing ${\mathrm{far}}_S\left(\mathit{p}\right)$ and $\mathrm{diam}\left(S\right)$ for sublevel sets.

The corresponding operators for ${\mathrm{far}}_S\left(\mathit{p}\right)$ and $\mathrm{diam}\left(S\right)$ can be multi-valued. Finding them is more challenging because the underlying optimization problem is no longer convex. The current paper proposes two algorithms for each problem. These algorithms are not infallible. They are ascent algorithms, so they tend to find local maxima. Convex hull algorithms can be harnessed to solve the farthest and diameter problems for finite point clouds in dimensions greater than 3 [13,14]. To our knowledge, these are the only competitive algorithms in common use in higher dimensions. Hence, any progress in solving these two fundamental geometric problems should be welcome.

Unfortunately, convex hull algorithms scale poorly in high dimensions. For n points in ${\mathbb{R}}^p$, an n-vertex polytope can have as many as $O\left(n^{\lfloor \frac{p}{2}\rfloor }\right)$ facets [15]. This translates into a worst-case computational computational complexity of $O\left(n^{\lfloor \frac{p}{2}\rfloor }\right)$ for traditional convex hull algorithms. It is true that once the m extreme points are extracted from the convex hull, the farthest and diameter problems require computing just m and $m^2$ distances, respectively. In contrast, we circumvent the convex hull problem entirely and attack the farthest and diameter problems directly. Thus, we easily handle problems in dimension $p=1000$.

Given the naturalness of the farthest and diameter problems, a few comments should suffice to motivate each. The minimum enclosing ball problem reduces to find

$$\begin{array}{ccc}\hfill \underset{\mathit{c}\in S}{\mathrm{argmin}}\underset{\mathit{x}\in S}{max}\parallel \mathit{x}-\mathit{c}\parallel & =& \underset{\mathit{c}\in S}{\mathrm{argmin}}\underset{S}{\mathrm{far}}\left(\mathit{c}\right).\hfill \end{array}$$

The solution point $\mathit{c}$ is called the Chebyshev center. Article [16] explains the pertinence to statistics. As for the diameter problem, convergence estimates for the Frank–Wolfe method of optimization depend on the squared diameter of the underlying set [17,18]. In general for an L-Lipschitz function $f\left(\mathit{x}\right)$, the inequality $\left|f\right(\mathit{x}\left)\right|\le \left|f\right(\mathit{y}\left)\right|+L\mathrm{diam}\left(S\right)$ bounds objective values. Bounding objective values is crucial in deriving complexity estimates for optimization algorithms such as projected gradient descent. If instead $\nabla f\left(\mathit{x}\right)$ possesses an L-Lipschitz gradient, then the bound

$$\begin{array}{ccc}\hfill \parallel \nabla f\left(\mathit{x}\right)-\nabla f\left(\mathit{y}\right)\parallel & \le & L\mathrm{diam}\left(S\right)\hfill \end{array}$$

is available for proving convergence to a stationary point. The diameter of a set also plays a crucial role in estimating the concentration of probability measures [19]. Finally, let us draw attention to a theorem of Jung to the effect that the radius r of the minimum enclosing ball problem satisfies $r\le \mathrm{diam}\left(S\right)\sqrt{\frac{p}{2(p+1)}}$ for $S\subset {\mathbb{R}}^p$ [20].

Our algorithms are minorization–maximization (MM) algorithms [7,21]. Such algorithms depend on a surrogate function that minorizes 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 two goals, 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 now the domination condition $g(\mathit{x}\mid {\mathit{x}}_n)\ge f\left(\mathit{x}\right)$ is reversed. The acronym MM also applies to majorization–minimization.

The celebrated EM (expectation–maximization) principle for maximum likelihood estimation with missing data [22] 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. Projected gradient descent, proximal gradient descent, and the convex–concave procedure [23] can also be viewed as MM algorithms.

The convexity of $\frac{1}{2}{\parallel \mathit{x}-\mathit{p}\parallel }^2$ yields the supporting hyperplane minorization

$$\begin{array}{ccc}\hfill \frac{1}{2}{\parallel \mathit{x}-\mathit{p}\parallel }^2& \ge & \frac{1}{2}{\parallel {\mathit{x}}_n-\mathit{p}\parallel }^2+{({\mathit{x}}_n-\mathit{p})}^{\top }(\mathit{x}-{\mathit{x}}_n).\hfill \end{array}$$

To improve ${\parallel \mathit{x}-\mathit{p}\parallel }^2$, we take

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& \in & \underset{\mathit{x}\in S}{\mathrm{argmax}}{\mathit{v}}^{\top }\mathit{x}\hfill \end{array}$$

for $\mathit{v}={\mathit{x}}_n-\mathit{p}$. This is the (naive) Frank–Wolfe algorithm [17,24] with full steps operating on the support function ${\sigma }_S\left(\mathit{v}\right)={sup}_{\mathit{x}\in S}{\mathit{v}}^{\top }\mathit{x}$ defined by S. The collection of points ${\mathrm{supp}}_S\left(\mathit{v}\right)={\mathrm{argmin}}_{\mathit{x}\in S}{\mathit{v}}^{\top }\mathit{x}$ at which the maximum is attained constitute the support map. Let us emphasize that the Frank–Wolfe algorithm is pertinent to the maximization of any differentiable convex function $f\left(\mathit{x}\right)$ over a compact convex set S because the supporting hyperplane minorization

$$\begin{array}{ccc}\hfill f\left(\mathit{x}\right)& \ge & f\left({\mathit{x}}_n\right)+df\left({\mathit{x}}_n\right)(\mathit{x}-{\mathit{x}}_n)\hfill \end{array}$$

generalizes to this context. The ascent property follows immediately from this minorization and the MM principle.

The Frank–Wolfe method applies to a host of problems. For example, suppose $g\left(\mathit{x}\right)$ maps ${\mathbb{R}}^p$ into itself. Consider the problem of minimizing the loss $f\left(\mathit{x}\right)=\frac{1}{2}{\parallel g\left(\mathit{x}\right)\parallel }^2$ over a compact convex set S. The function $f\left(\mathit{x}\right)$ has differential $df\left(\mathit{x}\right)=g{\left(\mathit{x}\right)}^{\top }dg\left(\mathit{x}\right)$ and gradient $\nabla f\left(\mathit{x}\right)=dg{\left(\mathit{x}\right)}^{\top }g\left(\mathit{x}\right)$. At iteration n, Frank–Wolfe decreases the linear approximation

$$\begin{array}{ccc}\hfill f\left(\mathit{y}\right)& \approx & f\left({\mathit{x}}_n\right)+g{\left({\mathit{x}}_n\right)}^{\top }dg\left({\mathit{x}}_n\right)(\mathit{y}-{\mathit{x}}_n)\hfill \end{array}$$

to $f\left(\mathit{y}\right)$. If the support function of S is simple to compute, then Frank–Wolfe additionally needs to compute just a single matrix times vector product $dg{\left(\mathit{x}\right)}^{\top }g\left(\mathit{x}\right)$ per iteration. The projection operator $P_S\left(\mathit{x}\right)$ is never invoked in updating $\mathit{x}$. In this minimization example, the MM principle is no longer operative, so a more sophisticated choice of step length is advisable. The problem of finding a root of $g\left(\mathit{x}\right)=\mathbf{0}$ on S also succumbs to a conjugate gradient method [25].

There are many known examples of support functions. All are closed, convex, positively homogeneous, and satisfy ${\sigma }_S\left(\mathbf{0}\right)=0$. The support function ${\sigma }_S\left(\mathit{y}\right)$ is by definition the Fenchel conjugate ${sup}_{\mathit{x}}[{\mathit{y}}^{\top }\mathit{x}-i_S\left(\mathit{x}\right)]$ of the $0/\infty$ indicator $i_S\left(\mathit{y}\right)$ of S. It is worth mentioning a few known support functions. The support function of a closed convex cone is the convex indicator of its polar cone. The support function of a polyhedral set can be found by linear programming. The support function of a union of two sets is the maximum of the two support functions. The support function of a Minkowski sum is the sum of the two support functions. The support function of the unit ball of a norm is the dual norm. The support function of a Cartesian product is the sum of the support functions of the two product sets.

The Frank–Wolfe version of the diameter problem exploits the two supporting hyperplane minorizations

$$\begin{array}{ccc}\hfill {\parallel \mathit{x}-\mathit{y}\parallel }^2& \ge & \parallel {\mathit{x}}_n-{\mathit{y}}_n{\parallel }^2+2{({\mathit{x}}_n-{\mathit{y}}_n)}^{\top }(\mathit{x}-{\mathit{x}}_n)\hfill \\ \hfill {\parallel \mathit{x}-\mathit{y}\parallel }^2& \ge & \parallel {\mathit{x}}_n-{\mathit{y}}_n{\parallel }^2+2{({\mathit{y}}_n-{\mathit{x}}_n)}^{\top }(\mathit{y}-{\mathit{y}}_n).\hfill \end{array}$$

Adding these and dividing by 2 yields the minorization

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

To improve ${\parallel \mathit{x}-\mathit{y}\parallel }^2$, we take

$$\begin{array}{ccc}\hfill ({\mathit{x}}_{n+1},{\mathit{y}}_{n+1})& \in & \underset{(\mathit{x},\mathit{y})\in S\times S}{\mathrm{argmax}}{\left(\genfrac{}{}{0.0pt}{}{{\mathit{x}}_n-{\mathit{y}}_n}{{\mathit{y}}_n-{\mathit{x}}_n}\right)}^{\top }\left(\genfrac{}{}{0.0pt}{}{\mathit{x}}{\mathit{y}}\right)\hfill \\ & =& \underset{\mathit{x}\in S}{\mathrm{argmax}}{({\mathit{x}}_n-{\mathit{y}}_n)}^{\top }\mathit{x}\times \underset{\mathit{y}\in S}{\mathrm{argmax}}{({\mathit{y}}_n-{\mathit{x}}_n)}^{\top }\mathit{y},\hfill \end{array}$$

which is a Frank–Wolfe update for the function $f(\mathit{x},\mathit{y})=\frac{1}{2}{\parallel \mathit{x}-\mathit{y}\parallel }^2$ on the Cartesian product $S\times S$.

For the sake of completeness, the Julia functions for our two Frank–Wolfe algorithms follow. The maximum number of iterations, 100, and the convergence tolerances, ${10}^{-8}$ and $2\times {10}^{-8}$, can be reset by the user.

function farthest(Supp::Function, x, p)

  tol = 1.0e-8

  for iter = 1:100

    xnew = Supp(x - p)

    conv = norm(xnew - x)

    x .= xnew

    if conv < tol

      break

    end

  end

  return (norm(x - p), x)

end

function widest(Supp::Function, x, y)

  tol = 2.0e-8

  for iter = 1:100

    (xnew, ynew) = (Supp(x - y), Supp(y - x))

    conv = norm(xnew - x) + norm(ynew - y)

    x .= xnew

    y .= ynew

    if conv < tol

      break

    end

  end

  return (norm(x - y), x, y)

end

Projected gradient ascent offers another avenue for solving the farthest and diameter problems. The two algorithms we later derive and apply are

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& =& P_S(2{\mathit{x}}_n-\mathit{p})\mathrm{for};\underset{S}{\mathrm{far}}\left(\mathit{p}\right)\hfill \\ \hfill \left(\begin{array}{c}{\mathit{x}}_{n+1}\\ {\mathit{y}}_{n+1}\end{array}\right)& =& \left[\begin{array}{c}{P}_S(\frac{3}{2}{\mathit{x}}_n-\frac{1}{2}{\mathit{y}}_n)\\ P_S(\frac{3}{2}{\mathit{y}}_n-\frac{1}{2}{\mathit{x}}_n)\end{array}\right]\mathrm{for};\mathrm{diam}\left(S\right).\hfill \end{array}$$

As MM algorithms, these two algorithms are guaranteed to increase the objective function. Fortunately, projected gradient ascent does not require the objective function to be convex. Once again the Julia code is straightforward.

function farthest(Proj::Function, x, p)

  tol = 1.0e-8

  for iter = 1:100

    xnew = Proj(2x - p)

    conv = norm(xnew - x)

    x .= xnew

    if conv < tol

      break

    end

  end

  return (norm(x - p), x)

end

function widest(Proj::Function, x, y)

  tol = 2.0e-8

  for iter = 1:100

    (xnew, ynew) = (Proj(3x / 2 - y / 2), Proj(3y / 2 - x / 2))

    conv = norm(xnew - x) + norm(ynew - y)

    x .= xnew

    y .= ynew

    if conv < tol

      break

    end

  end

  return (norm(x - y), x, y)

end

The contributions of this article include (a) deriving and testing these algorithms, (b) investigating simplifications arising from symmetry, (c) describing a homotopy method for avoiding local maxima, (d) finding the support function of the intersection of a convex cone and a ball centered at the origin, and (e) elaborating a known bisection algorithm for calculating the support function of a convex sublevel set. The next section derives the projected gradient ascent algorithm and lays out our contributions to items (b) through (e). Section 3 briefly tackles convergence of the various algorithms. Section 4 describes a few numerical experiments, and Section 5 discusses overall conclusions, limitations, and new directions for research. The author’s earlier article [26] applies similar techniques to the problem of computing the Hausdorff distance between two compact convex sets

Here are the notational conventions used throughout this article. All vectors appear in boldface. All entries of the vector $\mathbf{0}$ equal 0. 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)$, I 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, I 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. Derivations

2.1. Projected Gradient Ascent and Homotopy

Let us first tackle the closest, farthest, and diameter problems when the projection operator $P_S\left(\mathit{p}\right)$ is available. These three problems are equivalent to minimizing the functions

$$\begin{array}{ccc}\hfill c\left(\mathit{x}\right)& =& \frac{1}{2}{\parallel \mathit{x}-\mathit{p}\parallel }^2\hfill \\ \hfill f\left(\mathit{x}\right)& =& -\frac{1}{2}{\parallel \mathit{x}-\mathit{p}\parallel }^2\hfill \\ \hfill w(\mathit{x},\mathit{y})& =& -\frac{1}{2}{\parallel \mathit{x}-\mathit{y}\parallel }^2\hfill \end{array}$$

with gradients

$$\begin{array}{ccc}\hfill \nabla c\left(\mathit{x}\right)& =& \mathit{x}-\mathit{p}\hfill \\ \hfill \nabla f\left(\mathit{x}\right)& =& -(\mathit{x}-\mathit{p})\hfill \\ \hfill \nabla w(\mathit{x},\mathit{y})& =& \left(\begin{array}{c}\mathit{y}-\mathit{x}\\ \mathit{x}-\mathit{y}\end{array}\right)\hfill \end{array}$$

over the sets S and $S\times S$. One possibility for this task is projected gradient descent.

Because the Lipschitz constants of $\nabla c\left(\mathit{x}\right)$ and $\nabla f\left(\mathit{x}\right)$ are both 1, the projected gradient steps for the closest and farthest problems are

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& =& P_S[{\mathit{x}}_n-L^{-1}\nabla c\left({\mathit{x}}_n\right)]\hfill \\ & =& P_S({\mathit{x}}_n-{\mathit{x}}_n+\mathit{p})\hfill \\ & =& P_S\left(\mathit{p}\right)\hfill \\ \hfill {\mathit{x}}_{n+1}& =& P_S[{\mathit{x}}_n-L^{-1}\nabla f\left({\mathit{x}}_n\right)]\hfill \\ & =& P_S({\mathit{x}}_n-\mathit{p}+{\mathit{x}}_n)\hfill \\ & =& P_S(2{\mathit{x}}_n-\mathit{p}).\hfill \end{array}$$

As expected, the closest algorithm ${\mathit{x}}_{n+1}=P_S\left(\mathit{p}\right)$ reduces to ordinary projection. The Lipschitz constant for $\nabla w(\mathit{x},\mathit{y})$ is determined by

$$\begin{array}{ccc}\hfill ∥\left(\begin{array}{c}\mathit{y}-\mathit{x}\\ \mathit{x}-\mathit{y}\end{array}\right)-\left(\begin{array}{c}\mathit{v}-\mathit{u}\\ \mathit{u}-\mathit{v}\end{array}\right)∥& =& \sqrt{{2\parallel \mathit{y}-\mathit{x}-\mathit{v}+\mathit{u}\parallel }^2}\hfill \\ & \le & \sqrt{2{(\parallel \mathit{x}-\mathit{u}\parallel +\parallel \mathit{y}-\mathit{v}\parallel )}^2}\hfill \\ & \le & \sqrt{4(\parallel \mathit{x}-\mathit{u}{\parallel }^2+\parallel \mathit{y}-\mathit{v}{\parallel }^2)}\hfill \\ & =& 2∥\left(\begin{array}{c}\mathit{x}-\mathit{u}\\ \mathit{y}-\mathit{v}\end{array}\right)∥\hfill \end{array}$$

as 2. Consequently, the projected gradient step for the diameter problem is

$$\begin{array}{ccc}\hfill \left(\begin{array}{c}{\mathit{x}}_{n+1}\\ {\mathit{y}}_{n+1}\end{array}\right)& =& P_{S\times S}\left[\left(\begin{array}{c}{\mathit{x}}_n\\ {\mathit{y}}_n\end{array}\right)-\frac{1}{2}\nabla w({\mathit{x}}_n,{\mathit{y}}_n)\right]\hfill \\ & =& P_{S\times S}\left[\left(\begin{array}{c}{\mathit{x}}_n\\ {\mathit{y}}_n\end{array}\right)-\frac{1}{2}\left(\begin{array}{c}{\mathit{y}}_n-{\mathit{x}}_n\\ {\mathit{x}}_n-{\mathit{y}}_n\end{array}\right)\right]\hfill \\ & =& \left[\begin{array}{c}{P}_S(\frac{3}{2}{\mathit{x}}_n-\frac{1}{2}{\mathit{y}}_n)\\ P_S(\frac{3}{2}{\mathit{y}}_n-\frac{1}{2}{\mathit{x}}_n)\end{array}\right].\hfill \end{array}$$

This diameter update relies on the facts that minimization of a Euclidean distance is equivalent to minimization of a squared Euclidean distance and that the squared objective splits over the $\mathit{x}$ and $\mathit{y}$ parameters.

The chief problem with all of the proposed algorithms is their propensity to veer toward local maxima. One remedy is good initialization. Given a point $\mathit{p}\notin S$, a general tactic for initialization in the farthest problem is to start with the projection $P_S\left(\mathit{p}\right)$. This defines the line segment $[\mathit{p},P_S\left(\mathit{p}\right)]$, which can be extended until it hits the boundary of S at a second point ${\mathit{x}}_0$. Replacing $P_S\left(\mathit{p}\right)$ by ${\mathit{x}}_0$ produces a more distant point of S and a better starting value in the farthest problem. If S is an even set in the sense that $S=-S$, then as discussed in Section 2.3, the second point reduces to $P_S(-\mathit{p})$.

Beyond good initialization, it is worth considering the expensive alternative of homotopy [27] for the diameter problem. The idea is to gradually deform the unit ball B, where the diameter problem is trivial to solve, into the target set S. Thus, we follow the solution path along the family of sets $tS+(1-t)B$ from $t=0$ to $t=1$. For the Frank–Wolfe method, I exploit the fact that the Minkowski convex combination $tS+(1-t)B$ has support map $t{\mathrm{supp}}_S\left(\mathit{z}\right)+(1-t){\mathrm{supp}}_B\left(\mathit{z}\right)$.

For projected gradient ascent, one can project points onto the Minkowski convex combination $tS+(1-t)B$ by two devices. First, it is well known that $P_{tS}\left(\mathit{z}\right)=t{P}_S\left(t^{-1}\mathit{z}\right)$ for any $t>0$. Second, there is an effective algorithm for projecting onto a Minkowski sum $A+B$ [28]. The idea is to alternate minimization of $\parallel \mathit{z}-\mathit{a}-\mathit{b}\parallel$ with respect to $\mathit{a}\in A$ and $\mathit{b}\in B$. The iterative scheme ${\mathit{a}}_{n+1}=P_A(\mathit{z}-{\mathit{b}}_n)$ and ${\mathit{b}}_{n+1}=P_B(\mathit{z}-{\mathit{a}}_{n+1})$ is guaranteed to converge at a linear rate when either set is strongly convex. In particular, $(1-t)B$ is strongly convex. The homotopy method is motivated by the intuition that the early sets are more rounded and that the objective possesses fewer local maxima. The price for better performance is iterations within iterations and an overall slower algorithm.

The following Julia code implements the homotopy method for the diameter problem under Frank–Wolfe. The map ${\mathrm{supp}}_S\left(\mathit{y}\right)$ is passed to each of the functions. The code for projected gradient ascent is similar.

function MinkowskiSupp(Supp::Function, z, t)

  return t × Supp(z) + (1 - t) × (1 / norm(z)) × z

end

function widest_homotopy(Supp::Function, n)

  x = randn(n)

  x = x /norm(x) # random point on unit sphere

  y = -x  # point on opposite side of unit sphere

  (homotopy_points, tol) = (10, 2.0e-8)

  for iter = 0:homotopy_points

    t = iter / homotopy_points

    for i = 1:100

      xnew = MinkowskiSupp(Supp, x - y, t)

      ynew = MinkowskiSupp(Supp, y - x, t)

      conv = norm(x - xnew) + norm(y - ynew)

      x .= xnew

      y .= ynew

      if conv < tol

        break

      end

    end

  end

  return (norm(x - y), x, y)

end

2.2. Supporting Points and Sublevel 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)$. Our numerical tests require knowing ${\mathrm{supp}}_S\left(\mathit{v}\right)$ in some specific examples. For instance, the ${\mathsf{\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. Here, ${\mathit{e}}_i$ is included when $v_i>0$, and $-{\mathit{e}}_i$ is included when $v_i<0$. 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. The rectangle $[\mathit{a},\mathit{b}]$ is a Cartesian product. Hence, ${\mathrm{supp}}_S\left(\mathit{v}\right)$ is also a Cartesian product. In the one-dimensional case, ${\mathrm{supp}}_S\left(\mathit{v}\right)$ is a when $v_i<0$, b when $v_i>0$, and all of $[a,b]$ when $v_i=0$. 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)$. This fact plus the identity ${\mathrm{supp}}_{tA}\left(\mathit{v}\right)=t{\mathrm{supp}}_A\left(\mathit{v}\right)$ for $t\ge 0$ makes it easy to carry out the homotopy method with our Frank–Wolfe algorithms.

If ${\mathrm{supp}}_S\left(\mathit{v}\right)$ is a singleton, then ${\sigma }_S\left(\mathit{v}\right)$ is differentiable at $\mathit{v}$ by Danskin’s theorem [7]. Conversely, if ${\sigma }_S\left(\mathit{v}\right)$ is differentiable, then ${\mathrm{supp}}_S\left(\mathit{v}\right)$ is a singleton by Corollary 25.3.1 of Rockafellar [29]. Because ${\sigma }_S\left(\mathit{v}\right)$ is convex, finite, and locally Lipschits, Rademacher’s theorem [30,31] implies that it is differentiable almost everywhere. Hence, ${\sigma }_S\left(\mathit{v}\right)$ is a singleton almost everywhere.

Later, we will need ${\mathrm{supp}}_A\left(\mathit{v}\right)$ for the intersection of the non-negative orthant and the ball $B_r$ of radius r around the origin. This is a special case of projection onto the intersection of an arbitrary closed convex cone K and the ball. In this general setting, $\mathit{w}={\mathrm{supp}}_{K\cap B_r}\left(\mathit{v}\right)=\frac{r}{\parallel P_K\left(\mathit{v}\right)\parallel }{P}_K\left(\mathit{v}\right)$ when $P_K\left(\mathit{v}\right)\ne \mathbf{0}$, and $\mathit{w}=\mathbf{0}$, otherwise. To prove this assertion, it suffices to show that

$$\begin{array}{ccc}\hfill d_{\mathit{u}-\mathit{w}}{\sigma }_S\left(\mathit{v}\right)& =& {\mathit{v}}^{\top }(\mathit{u}-\mathit{w}) \le 0\hfill \end{array}$$

for any point $\mathit{u}\in K\cap B_r$. When $P_K\left(\mathit{v}\right)\ne \mathbf{0}$, the Moreau decomposition and the Cauchy-Schwarz inequality imply that

$$\begin{array}{ccc}\hfill {\mathit{v}}^{\top }(\mathit{w}-\mathit{u})& =& {[P_K\left(\mathit{v}\right)+P_{K^{\circ }}\left(\mathit{v}\right)]}^{\top }(\mathit{w}-\mathit{u})\hfill \\ & =& \frac{r}{\parallel P_K\left(\mathit{v}\right)\parallel }{\parallel P_K\left(\mathit{v}\right)\parallel }^2-{[P_K\left(\mathit{v}\right)+P_{K^{\circ }}\left(\mathit{v}\right)]}^{\top }\mathit{u}\hfill \\ & \ge & \frac{r}{\parallel P_K\left(\mathit{v}\right)\parallel }{\parallel P_K\left(\mathit{v}\right)\parallel }^2-P_K{\left(\mathit{v}\right)}^{\top }\mathit{u}\hfill \\ & \ge & r\parallel P_K\left(\mathit{v}\right)\parallel -r\parallel P_K\left(\mathit{v}\right)\parallel \hfill \\ & =& 0,\hfill \end{array}$$

where $K^{\circ }$ is the polar cone of K. Otherwise, $P_K\left(\mathit{v}\right)=\mathbf{0}$, and

$$\begin{array}{ccc}\hfill {\mathit{v}}^{\top }(\mathit{w}-\mathit{u})& =& -{[P_K\left(\mathit{v}\right)+P_{K^{\circ }}\left(\mathit{v}\right)]}^{\top }\mathit{u} = -P_{K^{\circ }}{\left(\mathit{v}\right)}^{\top }\mathit{u} \ge 0\hfill \end{array}$$

by definition.

The polar cone of the non-negative orthant consists of those vectors $\mathit{v}$ with $\mathit{v}\le \mathbf{0}$. The current intersection support map can also be deduced less directly from Proposition 2.2 of paper [32]. The corresponding problem of projecting onto the intersection of a ball and cone is treated by Lange [7] and Bauschke et al. [33].

As already mentioned, the support function of a sublevel set is generally challenging to compute. As a supplement to the brief discussion in Section 6.4.2 of the survey [12], consider the Lagrangian $-{\mathit{v}}^{\top }\mathit{x}+\mu [g\left(\mathit{x}\right)-c]$. If $g\left(\mathit{x}\right)$ is differentiable, then the Lagrangian satisfies the KKT stationary conditions

$$\begin{array}{ccc}\hfill \mathbf{0}& =& -\mathit{v}+\mu \nabla g\left(\mathit{x}\right)\hfill \\ \hfill \mathbf{0}& =& g\left(\mathit{x}\right)-c.\hfill \end{array}$$

Slater’s condition postulates the existence of a point $\mathit{y}$ with $g\left(\mathit{y}\right)<c$. Under this condition, there is $\mu \ge 0$ satisfying the KKT conditions. The solution is unique and can be found by bisection when $\mathit{v}\ne \mathbf{0}$ and $g\left(\mathit{x}\right)$ is not just convex but also strongly convex. If $h\left(\mathit{u}\right)$ is the inverse function of $\nabla g\left(\mathit{x}\right)$, then one seeks a root of the equation $\varphi \left(\mu \right)=g\left[h\left({\mu }^{-1}\mathit{v}\right)\right]=c$ by bisection. Recall that the equation $\nabla g\left(\mathit{x}\right)=\mathit{y}$ is uniquely solvable for all $\mathit{y}$ when $g\left(\mathit{x}\right)$ is strongly convex and differentiable. The function $\varphi \left(\mu \right)$ is strictly decreasing in $\mu$ because

$$\begin{array}{ccc}\hfill {\varphi }^{\prime }\left(\mu \right)& =& dg\left[h\left({\mu }^{-1}\mathit{v}\right)\right]dh\left({\mu }^{-1}\mathit{v}\right)(-{\mu }^{-2}\mathit{v})\hfill \\ & =& {\mu }^{-1}{\mathit{v}}^{\top }dh\left({\mu }^{-1}\mathit{v}\right)(-{\mu }^{-2}\mathit{v})\hfill \\ & =& -{\mu }^{-3}{\mathit{v}}^{\top }{d}^{2}g{\left[h\left({\mu }^{-1}\mathit{v}\right)\right]}^{-1}\mathit{v}\hfill \\ & <& 0\hfill \end{array}$$

when $g\left(\mathit{x}\right)$ is also twice differentiable. Finally, under strong convexity, ${\varphi }^{\prime }\left(\mu \right)$ tends to $-\infty$ as $\mu \downarrow 0$, so $\varphi \left(\mu \right)\to \infty$ as $\mu \downarrow 0$.

As an example, consider the sublevel set $\{\mathit{x}:g(\mathit{x})\le c\}$ determined by the elastic net function $g\left(\mathit{x}\right)={\parallel \mathit{x}\parallel }_1+\frac{\rho }{2}{\parallel \mathit{x}\parallel }^2$, which is separable but not fully differentiable. There is at least a partial inverse. Because the partial derivative $\frac{\partial }{\partial x_i}g\left(\mathit{x}\right)=\mathrm{sgn}\left(x_i\right)+\rho x_i$, the ith component of the inverse function $h\left(\mathit{u}\right)$ should satisfy

$$\begin{array}{ccc}\hfill u_i& =& \nabla g{\left[h\left(\mathit{u}\right)\right]}_i=\mathrm{sgn}\left[h{\left(\mathit{u}\right)}_i\right]+\rho h{\left(\mathit{u}\right)}_i.\hfill \end{array}$$

Hence, the function $h\left(\mathit{u}\right)$ with ith component $h{\left(\mathit{u}\right)}_i=\frac{\mathrm{sgn}\left(u_i\right)}{\rho }\left(\right|u_i{|-1)}_{+}$ serves as a partial inverse. For $\mathit{v}\ne \mathbf{0}$ and $c>0$, the function $\varphi \left(\mu \right)=g\left[h\left({\mu }^{-1}\mathit{v}\right)\right]$ is continuous, tends to 0 as $\mu \to \infty$ and to ∞ as $\mu \to 0$. Hence, the equation $\varphi \left(\mu \right)=c$ is solvable by the intermediate value theorem.

2.3. Symmetry

Let us first consider a permutationally invariant set S. In S, we can swap coordinates and remain within S. For a support function $\mathit{x}={\mathrm{supp}}_S\left(\mathit{y}\right)$, the swap criterion dictates that

$$\begin{array}{ccc}\hfill y_i{x}_j+y_j{x}_i& \le & y_i{x}_i+y_j{x}_j\iff \hfill \\ \hfill 0& \le & (y_i-y_j)(x_i-x_j).\hfill \end{array}$$

Thus, the components of the map $\mathit{x}={\mathrm{supp}}_S\left(\mathit{y}\right)$ should satisfy the inequality

$$\begin{array}{ccc}\hfill (y_i-y_j)(x_i-x_j)& \ge & 0\hfill \end{array}$$

for all i and j. In other words, the components of the output $\mathit{x}={\mathrm{supp}}_S\left(\mathit{y}\right)$ should be consistently ordered with the components of the input $\mathit{y}$ of ${\mathrm{supp}}_S\left(\mathit{y}\right)$.

For the distance function, the swap criterion is

$$\begin{array}{ccc}\hfill {(y_i-x_j)}^2+{(y_j-x_i)}^2& \ge & {(y_i-x_i)}^2+{(y_j-x_j)}^2\iff \hfill \\ \hfill 0& \le & (y_i-y_j)(x_i-x_j).\hfill \end{array}$$

Hence, the components of the output $\mathit{x}=P_S\left(\mathit{y}\right)$ should also be consistently ordered with the components of the input $\mathit{y}$ of $P_S\left(\mathit{y}\right)$.

For the farthest function, the swap criterion is

$$\begin{array}{ccc}\hfill {(y_i-x_j)}^2+{(y_j-x_i)}^2& \le & {(y_i-x_i)}^2+{(y_j-x_j)}^2\iff \hfill \\ \hfill 0& \ge & (y_i-y_j)(x_i-x_j).\hfill \end{array}$$

Hence, the components of the map $\mathit{x}={\mathrm{far}}_S\left(\mathit{y}\right)$ should satisfy the inequality

$$\begin{array}{ccc}\hfill (y_i-y_j)(x_i-x_j)& \le & 0\hfill \end{array}$$

for all i and j. In other words, the components should be reverse consistently ordered.

With each of these three functions, permutational invariance should inform our choice of algorithm initial points. Alternatively, if one opts to test the extreme points of a set S to find ${\mathrm{far}}_S\left(\mathit{y}\right)$, then most of these can be eliminated from contention by incompatibility with permutational invariance. A sublevel set $\{\mathit{x}:g(\mathit{x})\le c\}$ is permutationally invariant when $g\left(\mathit{x}\right)$ is permutationally invariant.

Finally, let us turn to the diameter problem for even sets $S=-S$. Choose $\mathit{x}\in S$ farthest from the origin. The Cauchy–Schwarz inequality $-{\mathit{y}}^{\top }\mathit{x}\le {\parallel \mathit{x}\parallel }^2$ for $\mathit{y}\in S$ implies that

$$\begin{array}{ccc}\hfill {\parallel \mathit{y}\parallel }^2-2{\mathit{y}}^{\top }\mathit{x}& \le & {3\parallel \mathit{x}\parallel }^2\iff \hfill \\ \hfill {\parallel \mathit{y}-\mathit{x}\parallel }^2& \le & {4\parallel \mathit{x}\parallel }^2 = {\parallel (-\mathit{x})-\mathit{x}\parallel }^2.\hfill \end{array}$$

Thus, the diameter $2\parallel \mathit{x}\parallel$ is achieved when $\parallel \mathit{x}\parallel$ is maximal. For instance, the diameter of the elastic net sublevel set $\{\mathit{x}:\parallel \mathit{x}{\parallel }_1+\frac{1}{2}\parallel \mathit{x}{\parallel }^2\le 1\}$ is $2(\sqrt{3}-1)$, regardless of the dimension. The maximum value $\sqrt{3}-1=0.73205$ of $\parallel \mathit{x}\parallel$ is achieved when all but one component of $\mathit{x}$ is 0. Sometimes a translate $S+\mathit{b}$ of S is even when S itself is not even. If this is the case, the diameter can be found by maximizing $\parallel \mathit{x}+\mathit{b}\parallel$ over S. A sublevel set $\{\mathit{x}:g(\mathit{x}\le c\}$ is even when $g\left(\mathit{x}\right)$ is an even function.

3. Convergence

Zangwill’s [34] theorem offers the quickest route to proving convergence of both projected gradient ascent and our simplified Frank–Wolfe algorithm. Unfortunately, Zangwill’s theorem says nothing about the rate of convergence. The theorem involves a solution set $\mathrm{\Gamma }\subset S$, a descent function $f\left(\mathit{x}\right)$, and an algorithm map $M\left(\mathit{x}\right)$. In our case, $\mathrm{\Gamma }$ consists of the stationary points $\mathit{x}\in S$ where $d_{\mathit{v}}f\left(\mathit{x}\right)\ge 0$ for all tangent vectors $\mathit{v}$. Note that this condition is necessary but not sufficient for $\mathit{x}$ to be a local minimum of $f\left(\mathit{x}\right)$. The algorithm map $M\left(\mathit{x}\right)$ is said to be closed at $\mathit{x}$ if whenever ${\mathit{x}}_n\in S$ converges to $\mathit{x}$ and ${\mathit{y}}_n\in M\left({\mathit{x}}_n\right)$ converges to $\mathit{y}$, then $\mathit{y}\in M\left(\mathit{x}\right)$. If $M\left(\mathit{z}\right)$ is single-valued and continuous at $\mathit{x}$, then $M\left(\mathit{z}\right)$ is certainly closed at $\mathit{x}$. Here is Zangwill’s theorem.

Proposition 1. 

Suppose that

• All iterates ${\mathit{x}}_{n+1}\in M\left({\mathit{x}}_n\right)$ fall in the compact set S.
• The map M is closed at $\mathit{x}$ when $\mathit{x}\notin \mathrm{\Gamma }$.
• The function $f\left(\mathit{x}\right)$ is continuous on S and satisfies $f\left[M\right(\mathit{x}\left)\right]\le f\left(\mathit{x}\right)$, with strict inequality for $\mathit{x}\notin \mathrm{\Gamma }$.

Then, the limit of any convergent sub-sequence ${\mathit{x}}_{n_m}$ of ${\mathit{x}}_n$ belongs to Γ.

Zangwill’s theorem also applies to maximization provided that the progress condition $f\left[M\right(\mathit{x}\left)\right]\le f\left(\mathit{x}\right)$ is replaced by $f\left[M\right(\mathit{x}\left)\right]\ge f\left(\mathit{x}\right)$, and the stationary condition $d_{\mathit{v}}f\left(\mathit{x}\right)\ge 0$ is replaced by $d_{\mathit{v}}f\left(\mathit{x}\right)\le 0$.

3.1. Convergence of Frank–Wolfe

To prove that the algorithm map is closed for our version of Frank–Wolfe, suppose that $df\left({\mathit{x}}_n\right){\mathit{y}}_n\ge df\left({\mathit{x}}_n\right)\mathit{y}$ for every $\mathit{y}\in S$ and that ${\mathit{x}}_n$ and ${\mathit{y}}_n$ converge to ${\mathit{x}}_{\infty }$ and ${\mathit{y}}_{\infty }$. Then assuming $df\left(\mathit{x}\right)$ is continuous, $df\left({\mathit{x}}_{\infty }\right){\mathit{y}}_{\infty }\ge df\left({\mathit{x}}_{\infty }\right)\mathit{y}$ for every $\mathit{y}\in S$. The ascent property is baked into Frank–Wolfe because it is a minorization–maximization algorithm. Finally, a stationary point $\mathit{x}$ satisfies $df\left(\mathit{x}\right)\mathit{v}\le 0$ for all tangent vectors $\mathit{v}$. The set of tangent vectors $\mathit{v}$ is the closure of the set of points $c(\mathit{y}-\mathit{x})$ with $\mathit{y}\in S$ and $c>0$. This is where the convexity of S 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 S$, which is equivalent to $\mathit{x}\in M\left(\mathit{x}\right)$. If $\mathit{x}\notin M\left(\mathit{x}\right)$, then the objective strictly increases. Thus, Zangwill’s theorem applies to Frank–Wolfe.

These arguments shed light on the rate of convergence as measured by closeness to stationarity [24,35]. Indeed, adding the inequality

$$\begin{array}{ccc}\hfill f\left({\mathit{x}}_{n+1}\right)-f\left({\mathit{x}}_n\right)& \ge & df\left({\mathit{x}}_n\right)({\mathit{x}}_{n+1}-{\mathit{x}}_n) \ge df\left({\mathit{x}}_n\right)(\mathit{y}-{\mathit{x}}_n)\hfill \end{array}$$

for arbitrary $\mathit{y}\in S$ leads by telescoping to

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

This in turn implies

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

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

Mangasarian [36] demonstrates that Frank–Wolfe converges in a finite number of iterations when the objective $f\left(\mathit{x}\right)$ is differentiable, convex, and bounded from above on a polyhedral set S. This result is more or less obvious given that the number of extreme points of S is finite. Yurtsever and Sra [37] show that the well-known convex–concave procedure (CCCP) [23] and its generalization are special cases of the Frank–Wolfe method with step-size selection. This insight quantifies the convergence rate of the CCCP.

3.2. Convergence of Projected Gradient Descent

In this scenario, let $f\left(\mathit{x}\right)$ be the function to be minimized. Furthermore, let L be the Lipschitz constant of $\nabla f\left(\mathit{x}\right)$ on S. Projected gradient descent minimizes the surrogate function $g(\mathit{x}\mid {\mathit{x}}_n)=f\left({\mathit{x}}_n\right)+df\left({\mathit{x}}_n\right)(\mathit{x}-{\mathit{x}}_n)+\frac{L}{2}{\parallel \mathit{x}-{\mathit{x}}_n\parallel }^2$ over S. The algorithm map

$$\begin{array}{ccc}\hfill {\mathit{x}}_{n+1}& =& M\left({\mathit{x}}_n\right) = P_S[{\mathit{x}}_n-L^{-1}\nabla f\left({\mathit{x}}_n\right)]\hfill \end{array}$$

is single-valued and continuous. The descent condition is again automatic. Proposition 7.3.2 of [7] now shows that Zangwill’s theorem applies provided that we identify stationary points as satisfying either of the equivalent conditions $M\left(\mathit{x}\right)=\mathit{x}$ or $f\left[M\right(\mathit{x}\left)\right]=f\left(\mathit{x}\right)$. These conditions are also equivalent to our postulated stationary condition because

$$\begin{array}{ccc}\hfill \nabla g({\mathit{x}}_n\mid {\mathit{x}}_n)& =& \nabla f\left({\mathit{x}}_n\right).\hfill \end{array}$$

In this regard, observe that $g(\mathit{x}\mid {\mathit{x}}_n)$ is strongly convex, so the requirement $dg({\mathit{x}}_n\mid {\mathit{x}}_n)\mathit{v}\ge 0$ for all tangent vectors $\mathit{v}$ is both necessary and sufficient for ${\mathit{x}}_n$ to be a global minimum of $g(\mathit{x}\mid {\mathit{x}}_n)$. Proposition 7.3.4 of my book [7] additionally proves that the collection of limit points of the MM sequence ${\mathit{x}}_{n+1}=M\left({\mathit{x}}_n\right)$ is compact and connected. Thus, when the extreme points of S are isolated, the MM sequence actually converges to one of them. Alternatively, in the setting of semi-algebraic sets and functions, Attouch et al. [38] prove full convergence using the tools of algebraic geometry.

To attack the rate of convergence of non-convex projected gradient descent, I present a simplified version of an argument featured by Beck [5]. The references [6,35,39,40,41] provide further background. Our point of departure is the observation that the stationary condition $df\left(\mathit{x}\right)(\mathit{y}-\mathit{x})\ge 0$ for all $\mathit{y}\in S$ is equivalent to satisfaction of the equation $\mathit{x}=P_S[\mathit{x}-s\nabla f\left(x\right)]$ for any $s>0$. The obtuse angle criterion

$$\begin{array}{ccc}\hfill {[\mathit{x}-s\nabla f\left(\mathit{x}\right)-\mathit{x}]}^{\top }(\mathit{y}-\mathit{x})& \le & 0\hfill \end{array}$$

for all $\mathit{y}\in S$ is both necessary and sufficient for $\mathit{x}=P_S[\mathit{x}-s\nabla f\left(x\right)]$. However, the obtuse angle criterion is just a disguised version of $df\left(\mathit{x}\right)(\mathit{y}-\mathit{x})\ge 0$ for all $\mathit{y}\in S$.

The quantity $\parallel {\mathit{x}}_n-P_S[{\mathit{x}}_n-L^{-1}\nabla f\left({\mathit{x}}_n\right)]\parallel =\parallel {\mathit{x}}_n-{\mathit{x}}_{n+1}\parallel$ serves as a measure of how far ${\mathit{x}}_n$ is from stationarity. The obtuse angle condition implies that

$$\begin{array}{ccc}\hfill {[{\mathit{x}}_n-L^{-1}\nabla f\left({\mathit{x}}_n\right)-{\mathit{x}}_{n+1}]}^{\top }({\mathit{x}}_n-{\mathit{x}}_{n+1})& \le & 0,\hfill \end{array}$$

which is equivalent to

$$\begin{array}{ccc}\hfill df\left({\mathit{x}}_n\right)({\mathit{x}}_{n+1}-{\mathit{x}}_n)& \le & -L\parallel {\mathit{x}}_n-{\mathit{x}}_{n+1}{\parallel }^2.\hfill \end{array}$$

It follows that

$$\begin{array}{ccc}\hfill f\left({\mathit{x}}_{n+1}\right)-f\left({\mathit{x}}_n\right)& \le & g({\mathit{x}}_{n+1}\mid {\mathit{x}}_n)-g({\mathit{x}}_n\mid {\mathit{x}}_n)\hfill \\ & \le & dg({\mathit{x}}_n\mid {\mathit{x}}_n)({\mathit{x}}_{n+1}-{\mathit{x}}_n)+\frac{L}{2}{\parallel {\mathit{x}}_{n+1}-{\mathit{x}}_n\parallel }^2\hfill \\ & =& df\left({\mathit{x}}_n\right)({\mathit{x}}_{n+1}-{\mathit{x}}_n)+\frac{L}{2}{\parallel {\mathit{x}}_n-{\mathit{x}}_{n+1}\parallel }^2\hfill \\ & \le & -L\parallel {\mathit{x}}_{n+1}-{\mathit{x}}_n{\parallel }^2+\frac{L}{2}{\parallel {\mathit{x}}_n-{\mathit{x}}_{n+1}\parallel }^2.\hfill \end{array}$$

Rearrangement of this inequality gives

$$\begin{array}{ccc}\hfill \parallel {\mathit{x}}_{n+1}-{\mathit{x}}_n{\parallel }^2& \le & \frac{2}{L}[f\left({\mathit{x}}_n\right)-f\left({\mathit{x}}_{n+1}\right)].\hfill \end{array}$$

Adding these inequalities and telescoping yield

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

In other words, the convergence rate is $O\left(\frac{1}{\sqrt{n}}\right)$ for the distance from stationarity.

4. Numerical Experiments

We tested the Frank–Wolfe and projected gradient ascent algorithms on five compact convex sets: (a) the box $[-\mathbf{1},\mathbf{1}]=\{\mathit{x}:\parallel \mathit{x}{\parallel }_{\infty }\le 1\}$, (b) the intersection of the unit ball and the non-negative orthant, (c) the probability simplex, (d) the ${\mathsf{\ell }}_1$ unit ball, and (e) the sublevel set $\{\mathit{x}:\parallel \mathit{x}{\parallel }_1+\frac{1}{2}\parallel \mathit{x}{\parallel }^2\le 1\}$ determined by the elastic net penalty. All five sets are permutationally invariant. Sets (a), (d), and (e) are even. These examples are representative, and all five exact diameters are available for comparison.

Table 1. Evaluation of the Frank–Wolfe algorithms.

Set	Dimension	Type	Homotopy	Fraction	Maximum	Seconds
box	2	farthest	no		2.9624	0.13
box	2	widest	no	1.0	2.8284	0.0708
box	2	widest	yes	1.0	2.8284	0.0709
ball ∩ orthant	2	farthest	no		2.5994	0.0817
ball ∩ orthant	2	widest	no	0.52	1.4142	0.0756
ball ∩ orthant	2	widest	yes	0.52	1.4142	0.0612
simplex	2	farthest	no		2.5465	0.0734
simplex	2	widest	no	1.0	1.4142	0.0697
simplex	2	widest	yes	1.0	1.4142	0.0707
L1 ball	2	farthest	no		2.5465	0.0754
L1 ball	2	widest	no	1.0	2.0	0.0764
L1 ball	2	widest	yes	1.0	2.0	0.0761
elastic net	2	farthest	no		2.2883	0.113
elastic net	2	widest	no	1.0	1.4641	0.0818
elastic net	2	widest	yes	1.0	1.4641	0.277
box	3	farthest	no		3.9572	0.0472
box	3	widest	no	1.0	3.4641	0.0458
box	3	widest	yes	1.0	3.4641	0.0651
ball ∩ orthant	3	farthest	no		3.2791	0.0443
ball ∩ orthant	3	widest	no	0.78	1.4142	0.0436
ball ∩ orthant	3	widest	yes	0.78	1.4142	0.0694
simplex	3	farthest	no		3.0727	0.045
simplex	3	widest	no	1.0	1.4142	0.0461
simplex	3	widest	yes	1.0	1.4142	0.0763
L1 ball	3	farthest	no		3.0727	0.0439
L1 ball	3	widest	no	1.0	2.0	0.0465
L1 ball	3	widest	yes	1.0	2.0	0.0671
elastic net	3	farthest	no		2.8509	0.0705
elastic net	3	widest	no	1.0	1.4641	0.0489
elastic net	3	widest	yes	1.0	1.4641	0.249
box	10	farthest	no		5.6758	0.0475
box	10	widest	no	1.0	6.3246	0.0454
box	10	widest	yes	1.0	6.3246	0.0612
ball ∩ orthant	10	farthest	no		3.6022	0.043
ball ∩ orthant	10	widest	no	1.0	1.4142	0.0449
ball ∩ orthant	10	widest	yes	1.0	1.4142	0.0848
simplex	10	farthest	no		3.2937	0.0449
simplex	10	widest	no	1.0	1.4142	0.0442
simplex	10	widest	yes	1.0	1.4142	0.0797
L1 ball	10	farthest	no		3.2937	0.0434
L1 ball	10	widest	no	1.0	2.0	0.0427
L1 ball	10	widest	yes	1.0	2.0	0.0696
elastic net	10	farthest	no		3.1077	0.0496
elastic net	10	widest	no	1.0	1.4641	0.0537
elastic net	10	widest	yes	1.0	1.4641	0.328
box	1000	farthest	no		59.61	0.0444
box	1000	widest	no	1.0	63.246	0.0505
box	1000	widest	yes	1.0	63.246	0.175
ball ∩ orthant	1000	farthest	no		31.962	0.0425
ball ∩ orthant	1000	widest	no	1.0	1.4142	0.0447
ball ∩ orthant	1000	widest	yes	1.0	1.4142	0.469
simplex	1000	farthest	no		31.325	0.0456
simplex	1000	widest	no	1.0	1.4142	0.046
simplex	1000	widest	yes	1.0	1.4142	0.548
L1 ball	1000	farthest	no		31.325	0.0427
L1 ball	1000	widest	no	1.0	2.0	0.0445
L1 ball	1000	widest	yes	1.0	2.0	0.436
elastic net	1000	farthest	no		31.3	0.523
elastic net	1000	widest	no	1.0	1.4641	0.277
elastic net	1000	widest	yes	1.0	1.4641	9.21

and

Table 2. Evaluation of the projected gradient ascent algorithms.

Set	Dimension	Type	Homotopy	Fraction	Maximum	Seconds
box	2	farthest	no		2.9624	0.802
box	2	widest	no	1.0	2.8284	0.0291
box	2	widest	yes	0.99	2.8284	0.129
ball ∩ orthant	2	farthest	no		2.5994	0.128
ball ∩ orthant	2	widest	no	0.78	1.0	0.0278
ball ∩ orthant	2	widest	yes	0.51	1.4142	0.09
simplex	2	farthest	no		2.5465	0.761
simplex	2	widest	no	1.0	1.4142	0.0259
simplex	2	widest	yes	0.7	1.4142	0.0906
L1 ball	2	farthest	no		2.5465	0.154
L1 ball	2	widest	no	1.0	2.0	0.0273
L1 ball	2	widest	yes	0.99	2.0	0.0912
elastic net	2	farthest	no		2.2883	0.159
elastic net	2	widest	no	1.0	1.4641	0.0279
elastic net	2	widest	yes	0.82	1.4641	0.112
box	3	farthest	no		3.9572	0.0269
box	3	widest	no	1.0	3.4641	0.0246
box	3	widest	yes	0.99	3.4641	0.0357
ball ∩ orthant	3	farthest	no		3.2791	0.0257
ball ∩ orthant	3	widest	no	0.9	1.0	0.0259
ball ∩ orthant	3	widest	yes	0.77	1.4142	0.0286
simplex	3	farthest	no		3.0727	0.0248
simplex	3	widest	no	0.64	1.4142	0.0258
simplex	3	widest	yes	0.93	1.4142	0.0291
L1 ball	3	farthest	no		3.0727	0.0255
L1 ball	3	widest	no	1.0	2.0	0.025
L1 ball	3	widest	yes	1.0	2.0	0.0296
elastic net	3	farthest	no		2.8509	0.0574
elastic net	3	widest	no	1.0	1.4641	0.0303
elastic net	3	widest	yes	0.74	1.4641	0.0616
box	10	farthest	no		5.6758	0.0283
box	10	widest	no	1.0	6.3246	0.0314
box	10	widest	yes	0.95	6.3246	0.0401
ball ∩ orthant	10	farthest	no		3.6022	0.0268
ball ∩ orthant	10	widest	no	1.0	1.0	0.0279
ball ∩ orthant	10	widest	yes	1.0	1.4142	0.0396
simplex	10	farthest	no		3.2937	0.0272
simplex	10	widest	no	0.08	1.0801	0.0271
simplex	10	widest	yes	0.88	1.4142	0.0327
L1 ball	10	farthest	no		3.2937	0.0271
L1 ball	10	widest	no	1.0	2.0	0.0272
L1 ball	10	widest	yes	0.97	2.0	0.0363
elastic net	10	farthest	no		3.1077	0.0415
elastic net	10	widest	no	1.0	1.4641	0.0289
elastic net	10	widest	yes	0.44	1.4641	0.0962
box	1000	farthest	no		59.61	0.0515
box	1000	widest	no	1.0	63.246	0.0914
box	1000	widest	yes	0.01	63.019	1.11
ball ∩ orthant	1000	farthest	no		31.962	0.0284
ball ∩ orthant	1000	widest	no	1.0	1.0	0.0377
ball ∩ orthant	1000	widest	yes	1.0	1.4142	0.209
simplex	1000	farthest	no		31.325	0.0466
simplex	1000	widest	no	0.02	1.0005	0.0588
simplex	1000	widest	yes	0.02	1.4142	0.632
L1 ball	1000	farthest	no		31.325	0.0435
L1 ball	1000	widest	no	1.0	2.0	0.0722
L1 ball	1000	widest	yes	0.52	2.0	0.566
elastic net	1000	farthest	no		31.3	0.609
elastic net	1000	widest	no	1.0	1.4641	0.29
elastic net	1000	widest	yes	0.03	1.4641	4.2

present our findings. The column headed “Fraction” records the fraction of the 100 trials that achieve the maximum objective. The column headed “Seconds” records cumulative execution times over 100 trials.

Homotopy offers no advantage in the farthest problem and is limited to the diameter problem. The farthest problem is initialized by $-\mathit{p}$ under Frank–Wolfe and by $P_C(-\mathit{p})$ under projected gradient ascent. Thus, all 100 trials are identical in each scenario. The diameter problem is initialized by two random opposing points on the unit sphere. Given the nature of the support and projection maps, all subsequent iterations fall within the set C. Both the Frank–Wolfe algorithms and the projected gradient descent algorithms perform well. They attain the same solutions, which in the diameter problems are identical to the known solution. Contrary to our expectation, projected gradient ascent is not noticeably slower than Frank–Wolfe. Homotopy over 11 intervening points rescues projected gradient descent on the diameter problem for sets (b) and (c). Homotopy is noticeably slower than the unadorned algorithms. The elastic net problem tends to take the most time owing to the inefficiency of bisection. Computation times scale well when the dimension d of the ambient space reaches 1000. 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 are embarrassingly parallel across trials, our Julia code is completely serial.

On the basis of reliability, these problems favor Frank–Wolfe. Although Frank–Wolfe algorithms possess a theoretically faster rate of convergence than projected gradient ascent, in the empirical trials, the two methods are quite comparable in speed. The “Seconds” columns of Table 1 and Table 2 have respective means of 0.260 and 0.240 and respective medians of 0.065 and 0.038. The corresponding standard deviations 1.192 and 0.697 are relatively high, so for this and other obvious reasons, one should probably not read too much into these crude comparisons.

5. Discussion

The farthest and diameter problems are natural problems of intrinsic interest. Given their non-convexity, they have not received nearly the attention in the mathematical literature as the closest problem. Exact mathematical solutions are available in some special cases. Research on fast algorithms appears 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 in high dimensions. Our preliminary experiments tilt toward the Frank–Wolfe algorithms as the more reliable of the two options. The full Julia code for our numerical examples appears at https://github.com/KennethLange/ClosestFarthestWidest (accessed on 19 January 2024).

Standard convergence arguments covered here guarantee that all limit points of the two algorithm classes are stationary points. I suspect, but have not proven, that full convergence to a stationary point always occurs for the Frank–Wolfe algorithms. This exercise would require a foray into the difficult terrain of real algebraic geometry [42]. 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. I have suggested symmetry tactics for choosing good starting points. The homotopy method capitalizes on exact solutions for the unit ball. Minkowski set rounding smooths the boundary of the target set and steers iterates in a productive direction. For semi-algebraic sets and functions, Attouch et al. [38] prove that the projected gradient algorithms always converge for objectives like ours with Lipschitz gradients. Again, convergence to a global maximum is not guaranteed.

I hope that this paper will provoke greater focus on the farthest and diameter problems. As prototype non-convex problems, they are worthy of further attention. In view of the connections to other Frank–Wolfe algorithms, I also encourage more community effort in finding and cataloging support functions ${\sigma }_S\left(\mathit{y}\right)$ and their corresponding maps ${\mathrm{supp}}_S\left(\mathit{v}\right)$. The effort put into this task so far is weaker than the effort put into devising projection maps. Once in possession of such maps, construction of fast algorithms is immensely easier. Finally, I would like to highlight the illumination that the MM principle brings to the construction of new high-dimensional optimization algorithms, including the ones considered here.