\section{Introduction}

\section{Finding the nearest point in a polytope}

In this section, we present a general acceleration technique for an \textit{arbitrary} algorithm for finding the nearest point in a
polytope, defined as the convex hull of a finite number of points in $\mathbb{R}^d$, in the case when the number of points
in the polytope is much greater than the dimension of the space. In other words, we discuss an acceleration
technique for methods of solving the following optimization problem
\[
  \min_x \: \| x - z \| \enspace \text{ subject to } \enspace
  x \in P := \co\big\{ x_1, \ldots, x_{\ell} \big\} \subset \mathbb{R}^d.
  \qquad \eqno{(\mathcal{P})}
\]
in the case when $\ell \gg d$. Here $z, x_1, \ldots, x_{\ell} \in \mathbb{R}^d$ are given points and $\| \cdot \|$ is
the Euclidean norm.

\subsection{General acceleration technique}

Before we proceed to the description of the acceleration technique, let us first recall the
following well-known optimality condition for the problem $(\mathcal{P})$ (see, e.g. \cite{Wolfe,MDM}),
for the sake of completeness.

\begin{proposition} \label{prp:NearPoint_OptCond}
A point $x_* \in P$ is a globally optimal solution of the problem $(\mathcal{P})$ if and only if
\begin{equation} \label{eq:NearPoint_OptCond}
  \langle x_* - z, x_i - x_* \rangle \ge 0 \quad \forall i \in I = \{ 1, \ldots, \ell \}.
\end{equation}
\end{proposition}

Suppose that an algorithm $\mathcal{A}$ for solving the nearest point problem is given. For any point
$w \in \mathbb{R}^d$ and any polytope $Q = \co\{ u_1, \ldots, u_m \}$ the algorithm returns a unique solution
$x_* = \mathcal{A}(w, Q)$ of the problem
\[
  \min_{x \in Q} \: \| x - w \|.
\]
No other assumptions on the algorithm $\mathcal{A}$ are imposed.

Our aim is to design an acceleration technique for the algorithm $\mathcal{A}$ that would improve its performance. This
acceleration technique will be presented as a meta-algorithm that utilises the algorithm $\mathcal{A}$ as a subroutine
on each step.

Recall that we are mainly interested in the case when the number of points $\ell$ is significantly greater than the
dimension of the space $d$. Furthermore, nothing is known about the structure of the algorithm $\mathcal{A}$.
Therefore, perhaps, the only straightforward way to potentially accelerate this algorithm is by applying $\mathcal{A}$
to a polytope generated by a relatively small number of points $\{ x_{i_1}, \ldots, x_{i_s} \}$ from the set
$\{ x_1, \ldots, x_{\ell} \}$ and then replacing some of the extreme points of the polytope
$P_0 = \co\{ x_{i_1}, \ldots, x_{i_s} \}$ by different points from $P$ and repeating the same procedure till
the projection of $z$ onto $P_0$ coincides with the projection of $z$ onto $P$. The key ingredient of this strategy is
\textit{an exchange rule} that replaces points from $P_0$ by different points from $P$.

To define an accelerating meta-algorithm for solving the problem $(\mathcal{P})$, we need to choose parameter
$s \in \mathbb{N}$ that defines the size of the \textit{subpolytope} $P_0$ and an exchange rule $\mathcal{E}$ whose
output consists of two index sets $J_1$ and $J_2$. The first set $J_1 \subseteq \{ i_1, \ldots, i_s \}$ defines
the points in $P_0$ that must be removed from $P_0$, while the second index set
$J_2 \subseteq I \setminus \{ i_1, \ldots, i_s \}$ defines the points that are included into $P_0$ before the next
iteration. In the general case, the index sets $J_1$ and $J_2$ might have arbitrary sizes and $J_1$ can be empty, but
for the sake of simplicity we restrict our consideration to the case when the cardinalities $|J_1|$ and $|J_2|$ of the
sets $J_1$ and $J_2$ coincide and are equal to some number $q \in \mathbb{N}$. Any exchange rule $\mathcal{E}$ that for
a given input consisting of $s$ indices $\{ i_1, \ldots, i_s \}$ returns two collections of indices $(J_1, J_2)$ with
\[
  J_1 \subseteq \{ i_1, \ldots, i_s \}, \quad
  J_2 \subseteq I \setminus \{ i_1, \ldots, i_s \}, \quad
  |J_1| = |J_2| = q
\]
is called an $(s, q)$-\textit{exchange rule}.

\begin{remark}
It should be noted that an actual exchange rule obviously requires more input parameters than just a set of indices
$\{ i_1, \ldots, i_s \}$. In particular, it might need some information about the point $z$, the set
$\{ x_1, \ldots, x_{\ell} \}$, the nearest point $\mathcal{A}(z, P_0)$ in the polytope $P_0$, etc. as its input.
However, for the sake of shortness, below we explicitly indicate only a set of indices as an input of an exchange rule
$\mathcal{E}$.
\end{remark}

Thus, we arrive at the following scheme of the meta-algorithm for solving the problem $(\mathcal{P})$ given in
Algorithm~\ref{alg:NearestPoint}. In its core, this meta-algorithm consists in choosing a \textit{subpolytope} $P_n$ of
the original polytope $P$ and finding the projection $y_n$ of the given point $z$ onto $P_n$. If this projection happens
to coincide with the projection of $z$ onto $P$ (this fact is verified with the use of optimality conditions from
Prp.~\ref{prp:NearPoint_OptCond}), then the meta-algorithm terminates. Otherwise, one replaces some points in $P_n$,
thus constructing a new subpolytope $P_{n + 1}$, and repeats the same procedure till the projection of $z$ onto $P_n$
coincides with the projection of $z$ onto the original polytope $P$.

From the geometric point of view, the meta-algorithm consists in choosing a ``small'' subpolytope $P_n$ in the original
polytope $P$ and gradually ``shifting'' it with each iteration till $P_n$ contains the projection of $z$ onto $P$. This
procedure is performed with the hope that it is much faster to compute a projection of a point onto a small subpolytope
and gradually shift it, rather than to compute the projection of this point onto the original polytope with a very large
number of vertices (i.e. with $\ell \gg d$). The results of numerical experiments reported below demonstrate that this
hope is fully justified.

\begin{algorithm}[t]	\label{alg:NearestPoint}
\caption{Meta-algorithm for finding the nearest point in a polytope.}

\noindent\textbf{Input:} {a point $z \in \mathbb{R}^d$, a collection of points $\{ x_1, \ldots, x_{\ell} \} \subset
\mathbb{R}^d$, an algorithm $\mathcal{A}$ for solving the nearest point problem, parameters
$s, q \in \{ 1, \ldots, \ell \}$ with $s \ge q$, and an $(s, q)$-exchange rule $\mathcal{E}$.}

\noindent\textbf{Initialization:} {Put $n = 0$, choose an index set $I_n \subseteq I$ with $|I_n| = s$, and define
$P_n = \co\{ x_i \mid i \in I_n \}$.}

\noindent\textbf{Step 1.} {Compute $y_n = \mathcal{A}(z, P_n)$. If $y_n$ satisfies the optimality condition
\begin{equation} \label{eq:TrialPointOptCond}
  \langle y_n - z, x_i - y_n \rangle \ge 0 \quad \forall i \in I,
\end{equation}
\textbf{return} $y_n$.}

\noindent\textbf{Step 2.} {Compute $(J_{1n}, J_{2n}) = \mathcal{E}(I_n)$ and define
\[
  I_{n + 1} = \Big( I_n \setminus J_{1n} \Big) \cup J_{2n}, \quad
  P_{n + 1} = \co\big\{ x_i \bigm| i \in I_{n + 1} \big\}.
\]
Set $n = n + 1$ and go to \textbf{Step 1}.}
\end{algorithm}

The two following lemmas describe simple conditions on the exchange rule $\mathcal{E}$ ensuring that the proposed
meta-algorithm indeed solves the problem $(\mathcal{P})$ and, furthermore, terminates after a finite number of
iterations. These lemmas provide one with two convenient criteria for choosing effective exchange rules. Although
(rather awkward) proofs of these lemmas are obvious, we present them for the sake of completeness.

For any set $\Omega \subset \mathbb{R}^d$ and any $x \in \mathbb{R}^d$ denote by
$\dist(x, \Omega) = \inf_{y \in \Omega} \| x - y \|$ the distance from $x$ to $\Omega$.

\begin{lemma} \label{lem:DistDecay}
Suppose that the exchange rule $\mathcal{E}$ satisfies the distance decay condition: if for some $n \in \mathbb{N}$ the
point $y_n$ does not satisfy optimality condition \eqref{eq:TrialPointOptCond}, then
\[
  \dist(z, P_{n + 1}) < \dist(z, P_n).
\]
Then Algorithm~\ref{alg:NearestPoint} terminates after a finite number of steps and returns an optimal solution of the
problem $(\mathcal{P})$.
\end{lemma}

\begin{proof}
From Proposition~\ref{prp:NearPoint_OptCond} and the termination criterion \eqref{eq:TrialPointOptCond} of
Algorithm~\ref{alg:NearestPoint} it follows that if this algorithm terminates in a finite number of steps, then the last
computed point $y_n$ (and the output of Algorithm~\ref{alg:NearestPoint}) is an optimal solution of the problem
$(\mathcal{P})$. Thus, we only need to prove that the algorithm does terminate in a finite number of steps.

To prove the finite termination, note that there is only a finite number of distinct subsets of the set
$\{ x_1, \ldots, x_{\ell} \}$ with cardinality $s$. Moreover, from the distance decay condition
\eqref{eq:DistDecayCondition} it follows that if the algorithm does not terminate in $n \in \mathbb{N}$ iterations, then
all index sets $I_0, I_1, \ldots, I_n$ are distinct. Therefore, in a finite number of steps the algorithm must find a
point $y_n$ satisfying the termination
criterion \eqref{eq:TrialPointOptCond}.
\end{proof}

\begin{lemma} \label{lem:SubpolytopeSize}
Let an $(s, q)$-exchange rule $\mathcal{E}$ with $s \ge q$ satisfy the distance decay condition for any point
$z \in \mathbb{R}^d$ and any polytope $P = \co\{ x_1, \ldots, x_{\ell} \}$. Then $s \ge d + 1$.
\end{lemma}

\begin{proof}
Arguing by reductio ad absurdum, suppose that $s \le d$. Let us furnish a particular point $z$ and a particular polytope
$P$ for which any $(s, q)$-exchange rule with $s \le d$ and $s \ge q$ fails to satisfy the distance decay condition.

Let $z = 0$. If $d = 1$, define $\ell = 2$, $x_1 = 1$, and $x_2 = -1$. In this case $s = q = 1$ and for any choice of
$I_0$ it is obviously impossible to satisfy the condition $\dist(z, P_1) < dist(z, P_2)$.

Let now $d \ge 2$. Put $z = 0$, $\ell = d + 1$, and define the points $x_i$ as follows:
\begin{gather} \label{eq:CounterExample1}
  x_1 = (1, 0, \ldots, 0, -1), \quad x_2 = (0, 1, 0, \ldots, 0, -1), \ldots,
  x_{d - 2} = (0, \ldots, 0, 1, 0, -1),
  \\
  x_{d - 1} = (-1, \ldots, -1, 1, -1), \quad x_d = (-1, \ldots, -1), \quad x_{d + 1} = (0, \ldots, 0, 1).
  \label{eq:CounterExample2}
\end{gather}
As is easily seen, any $d$ points from the set $\{ x_1, \ldots, x_{d + 1} \}$ are linearly independent and $0 \in P$,
since
\[
  \sum_{i = 1}^{d + 1} \alpha_i x_i = 0, \quad \sum_{i = 1}^{d + 1} \alpha_i = 1
\]
for
\[
  \alpha_1 = \ldots = \alpha_{d - 2} = \frac{1}{2(d - 1)}, \quad \alpha_{d - 1} = \alpha_d = \frac{1}{4(d - 1)},
  \quad \alpha_{d + 1} = \frac{1}{2}.
\]
Hence, in particular, for any index set $K \subset I$ with $|K| = s \le d$ one has $\dist(z, P_K) > 0$, where
$P_K = \co\{ x_i \mid i \in K \}$. Therefore, for any choice of an $(s, q)$-exchange rule with $s \le d$ and $s \ge q$
the stopping criterion \eqref{eq:TrialPointOptCond} cannot be satisfied, which by the previous lemma implies that this
exchange rule does not satisfy the distance decay condition.
\end{proof}

\begin{remark}
One can readily verify that if in the lemma above one imposes the additional assumption that $z \notin P$, then the
statement of the lemma holds true for $s \ge d$. To prove this result, one simply needs to put $z = 0$, $\ell = d$, and
define $P$ as the convex hull of the first $d$ points from \eqref{eq:CounterExample1}, \eqref{eq:CounterExample2}.
\end{remark}

\subsection{The steepest descent exchange rule}
\label{subsect:ExchangeRule}

Let us present a particular exchange rule based on the optimality condition from Proposition~\ref{prp:NearPoint_OptCond}
(or, equivalently, the stopping criterion \eqref{eq:TrialPointOptCond}) and satisfying the distance decay condition for
any point $z$ and any polytope $P$.

Bearing in mind Lemma~\ref{lem:SubpolytopeSize}, we propose to consider the following $(d + 1, 1)$-exchange rule that on
each iteration of Algorithm~\ref{alg:NearestPoint} replaces only one point in the current polytope $P_n$. Suppose that
for some $n \in \mathbb{N}$ the stopping criterion \eqref{eq:TrialPointOptCond} is not satisfied. Then we define the new
point $x_{1n}$ from $P$, that is included into $P_{n + 1}$, as any point from $P$ on which the minimum in
\[
  \min_{i \in I} \langle y_n - z, x_i - y_n \rangle \ge 0
\]
is attained (cf. a similar rule for including new points in the major cycle of the Wolfe method \cite{Wolfe}).

To find a point that is removed from $P_n$, note that $z \notin P_n$ due to the definition of $y_n$ and the fact that
condition \eqref{eq:TrialPointOptCond} does not hold true. Therefore, $y_n$ belongs to the boundary of $P_n$, which by
\cite[Lemma~2.8]{Zeigler} implies that $y_n$ is contained in a face of $P_n$ of dimension at most $d - 1$.
Hence by \cite[Prp.~1.15 and 2.3]{Zeigler} the point $y_n$ can be represented as a convex combination of at
most $d$ points from the set $\{ x_i \mid i \in I_n \}$. In other words, there exists $i_{2n} \in I_n$ such that
$y_n \in \co\{ x_i \mid I_n \setminus \{ i_{2n} \} \}$, and it is natural to remove the point $x_{i_{2n}}$
from the polytope $P_n$.

Thus, we arrive at the following theoretical scheme of a $(d + 1, 1)$-exchange rule given in Algorithm
\ref{alg:SteepestDescentExchRule} that we call \textit{the steepest descent exchange rule}. From the discussion above
it follows that this exchange rule is correctly defined, provided $y_n$ does not satisfy the stopping criterion
\eqref{eq:TrialPointOptCond}. Let us also verify that this exchange rule always satisfies the distance decay condition.

\begin{algorithm}[t] \label{alg:SteepestDescentExchRule}
\caption{The steepest descent exchange rule.}

\noindent\textbf{Input:} an index set $I_n \subset I$ with $|I_n| = d + 1$, the point $z$,
the set $\{ x_1, \ldots, x_{\ell} \}$, and the projection $y_n$ of $z$ onto $P_n = \co\{ x_i \mid i \in I_n \}$.

\noindent\textbf{Step 1:} Find $i_{1n} \in I_n$ such that $y_n \in \co\{ x_i \mid I_n \setminus \{ i_{1n} \} \}$.

\noindent\textbf{Step 2:} Find $i_{2n} \in I$ such that
\[
  \langle y_n - z, x_{i_{2n}} - y_n \rangle = \min_{i \in I} \langle y_n - z, x_i - y_n \rangle.
\]
\textbf{Return} $(\{ i_{1n} \}, \{ i_{2n} \})$.
\end{algorithm}

\begin{theorem} \label{thm:StDescExchRule}
For any point $z \in \mathbb{R}^d$ and any polytope $P = \co\{ x_1, \ldots, x_{\ell} \}$ with $\ell > d + 1$
the steepest descent exchange rule satisfies the distance decay condition.
\end{theorem}

\begin{proof}
Suppose that for some $n \in \mathbb{N}$ the stopping criterion \eqref{eq:TrialPointOptCond} does not hold true. Denote
by
\[
  P_n^0 = \co\{ x_i \mid i \in I_n \setminus \{ i_{1n} \} \}
\]
the polytope obtained from $P_n$ after removing the point selected by the steepest descent exchange rule. By
construction $y_n \in P_n^0$. Moreover, due to the definition of the exchange rule and the fact that the stopping
criterion is not satisfied one has
\begin{equation} \label{eq:GradAtNewPoint}
  \langle y_n - z, x_{i_{2n}} - y_n \rangle < 0.
\end{equation}
Define $x_n(\alpha) = (1 - \alpha) y_n + \alpha x_{i_{2n}}$. Clearly, $x_n(\alpha) \in P_{n + 1}$ for any
$\alpha \in [0, 1]$, since $P_{n + 1} = \co\{ P_n^0, x_{i_{2n}} \}$. Moreover, for any $\alpha \in \mathbb{R}$ one has
\begin{align*}
  f(\alpha) &:= \| x_n(\alpha) - z \|^2
  \\
  &= (1 - \alpha)^2 \| y_n - z \|^2 +
  2 \alpha (1 - \alpha) \langle y_n - z, x_{i_{2n}} - z \rangle + \alpha^2 \| x_{i_{2n}} - z \|^2.
\end{align*}
Note that $f(0) = \| y_n - z \|^2$ and
\[
  f'(0) = - 2 \| y_n - z \|^2 + 2 \langle y_n - z, x_{i_{2n}} - z \rangle
  = \langle y_n - z, x_{i_{2n}} - y_n \rangle < 0
\]
due to \eqref{eq:GradAtNewPoint}. Therefore, for any sufficiently small $\alpha \in (0, 1)$ one has
$f(\alpha) < f(0)$, which implies that
\begin{align*}
  \dist(z, P_{n + 1}) &= \min_{x \in P_{n + 1}} \| z - x \| \le \| z - x_n(\alpha) \|
  \\
  &= \sqrt{f(\alpha)} < \sqrt{f(0)} = \| z - y_n \| = \dist(z, P_n).
\end{align*}
Thus, the steepest descent exchange rule satisfies the distance decay condition.
\end{proof}

\begin{corollary}
Let $\ell > d + 1$. Then Meta-algorithm~\ref{alg:NearestPoint} with $s = d + 1$, $q = 1$, and the steepest descent exchange rule
terminates after a finite number of iterations and returns an optimal solution of the problem $(\mathcal{P})$.
\end{corollary}

Let us discuss a possible implementation of the steepest descent exchange rule. Clearly, the most difficult part of this rule
consists in finding an index $i_{1n} \in I_n$ such that $y_n \in \co\{ x_i \mid I_n \setminus \{ i_{1n} \} \}$. This difficulty
can be overcome in the following way.

Namely, suppose that instead of returning a projection of a point $w$ onto a polytope $Q = \co\{ u_1, \ldots, u_m \}$, the output
$\mathcal{A}(w, Q)$ of the algorithm $\mathcal{A}$ is actually a vector $\alpha \in \mathbb{R}^{\ell}$ such that
\[
	u_* = \sum_{i = 1}^m \alpha_i u_i, \quad \sum_{i = 1}^m \alpha_i = 1, \quad \alpha_i \ge 0
	\quad \forall i \in \{ 1, \ldots, m \},
\]
where $u_*$ is an optimal solution of the problem $\min_{u \in Q} \| u - w \|$. Let us note that for the vast majority of existing
methods for finding the nearest point in a polytope this assumption either holds true by default or can be satisfied by slightly
modifying the corresponding method.

Let the points $x_i$, $i \in I_n$, be affinely independent for some $n \in \mathbb{N}$, and suppose that the stopping criterion
\eqref{eq:TrialPointOptCond} is not satisfied. Then, as was noted above $y_n$ does not belong to the interior of $P_n$, which by
the characterization of interior points of a polytope \cite[Lemma~2.8]{Zeigler} implies that $y_n$ \textit{cannot} be represented
in the form
\[
	y_n = \sum_{i \in I_n} \alpha_i x_i, \quad \sum_{i \in I_n} \alpha_i = 1, \quad \alpha_i > 0 \quad \forall i \in I_n.
\]
Consequently, at least of of the coordinates, say $\alpha_{ni_{1n}}$ of the vector $\alpha_n = \mathcal{A}(z, P_n)$ is equal to zero,
and one can choose this $i_{1n}$ of Step~1 of the steepest descent exchange rule.

In the case when the points $x_i$, $i \in I_n$, are affinely dependent and all $\alpha_{ni} > 0$, $i \in I_n$, one can apply
the following simple procedure, inspired the Wolfe method \cite{Wolfe}, to find the required index $i_{1n}$. This procedure
is described in Algorithm~\ref{alg:IndexReplSearch}. The proposition below proves the correctness of the proposed method.

\begin{algorithm}[t]	\label{alg:IndexReplSearch}
\caption{Index removal method.}

\noindent\textbf{Step 1.} {Compute an optimal solution $\beta_n \in \mathbb{R}^{d + 1}$ of the problem
\[
	\min_{\beta} \Big\| \sum_{i \in I_n} \beta_i x_i - z \Big \|^2 \quad \text{subject to} \quad \sum_{i \in I_n} \beta_i = 1,
\]
that is, find a projection of $z$ onto the affine subspace generated by the points $x_i$, $i \in I_n$.}

\noindent\textbf{Step 2.} {If there exists $i \in I_n$ such that $\beta_{ni} < 0$, then find an index $k \in I_n$ on which
the minimum in
\begin{equation} \label{eq:WolfeCoefficient}
	\min\Big\{ \frac{\alpha_{ni}}{\alpha_{ni} - \beta_{ni}} \Bigm| i \in I_n \colon \alpha_{ni} - \beta_{ni} > 0 \Big\}
\end{equation}
is attained and \textbf{return} $j$.}

\noindent\textbf{Step 3.} {For each $j \in I_n$ check whether $y_n \in P_n^j := \co\{ x_i \mid i \in I_n \setminus \{ j \} \}$
(e.g. by finding the distance between $y_n$ and $P_n^j$ via the algorithm $\mathcal{A}$). If $y_n \in P_n^j$, \textbf{return} $j$.
}
\end{algorithm}

\begin{proposition}
Let the points $x_i$, $i \in I_n$, be affinely dependent for some $n \in \mathbb{N}$, $\alpha_{ni} > 0$ for all $i \in I_n$,
where $\alpha_n = \mathcal{A}(z, P_n)$, and $j \in I_n$ be the output of Algorithm~\ref{alg:IndexReplSearch}.
Then $y_n \in P_n^j = \co\{ x_i \mid i \in I_n \setminus \{ j \} \}$.
\end{proposition}

\begin{proof}
Suppose at first that the algorithm executes step 3. By our assumptions on the point $x_i$ the dimension of the polytope
$P_n = \co\{ x_i \mid i \in I_n \}$ is smaller than $d$. Therefore by \cite[Prp.~1.15]{Zeigler} the point $y_n$ can be
represent as the convex combination with positive coefficients of at most $d$ points from the set $\{ x_i \mid i \in I_n \}$.
Therefore, there exists $j \in I_n$ such that $y_n^j \in P_n^j$, which implies the required result.

Suppose now that the algorithm does not execute Step 2, i.e. that there exists $i \in I_n$ such that $\beta_{ni} < 0$.
Denote by $\theta$ the minimum in \eqref{eq:WolfeCoefficient}, and put
\[
	u_n = \sum_{i \in I_n} \beta_{ni} x_i, \quad w_n = (1 - \theta) y_n + \theta u_n.
\]
Note that by definition $\theta > 0$ and $\| u_n \| \le \| y_n \|$, which, in particular, implies that $\| w_n \| \le \| y_n \|$.
Moreover, one has
\[
	w_n = \sum_{i \in I_n} \gamma_{ni} x_i, \quad \gamma_{ni} = (1 - \theta) \alpha_{ni} + \theta \beta_{ni},
	\quad \sum_{i \in I_n} \gamma_{ni} = 1.
\]
In addition, if $\beta_{ni} \ge 0$, then clearly $\gamma_{ni} \ge 0$, while in the case $\beta_{ni} < 0$ one has
\[
	(1 - \theta) \alpha_{ni} + \theta \beta_{ni} = \alpha_{ni} - \theta(\alpha_{ni} - \beta_{ni})
	= (\alpha_{ni} - \beta_{ni}) \Big( \frac{\alpha_{ni}}{\alpha_{ni} - \beta_{ni}} - \theta \Big) \ge 0,
\]
and by definition of $\theta$ one has $\gamma_{nj} = 0$. Thus, $w_n \in P_n^j$. Moreover, from the inequality
$\| w_n \| \le \| y_n \|$ and the uniqueness of the projection onto a convex set it follows that $w_n = y_n$, which implies
the desired result.
\end{proof}

\begin{remark}
In the worst case scenario Algorithm~\ref{alg:IndexReplSearch} has to solve a projection problem and check the inclusion
$y_n \in P_n^j$ $d$ times to find the required index $j$. Such situation occurs only when the algorithm executes step 3, which happens
only in the case when the points $x_i$, $i \in I_n$, are affinely dependent, $\alpha_{ni} > 0$ for all $i \in I_n$, where
$\alpha_n = \mathcal{A}(z, P_n)$, and, furthermore, the projection of $z$ onto the affine subspace generated by the points $x_i$,
$i \in I_n$, coincides with $y_n$. It should be pointed out that in our numerical experiments on various test problems such situation
has never occurred, although one can easily provide an artificial example of such behaviour of Algorithm~\ref{alg:IndexReplSearch}.
\end{remark}

\subsection{Robustness with respect to round-off errors}

It is clear that any algorithm $\mathcal{A}$ for solving the problem $(\mathcal{P})$ can find an optimal solution of this problem
only in theory, while in practice it always returns an approximate solution of this problem due to the fact that the computations
are performed with finite precision and only a finite number of iterations can be performed. Therefore it is an important practical
issue to analyse the performance of Meta-algorithm~\ref{alg:NearestPoint} in the case when only computations with finite precision
are possible.

Assume that instead of an ``ideal'' algorithm $\mathcal{A}$ one uses its ``approximate'' version $\mathcal{A}_{\varepsilon}$,
$\varepsilon > 0$, that returns an approximate, in some sense (that will be specified below), solution of the nearest point problem.
Then it is obvious that the stopping criterion
\[
  \langle y_n - z, x_i - y_n \rangle \ge 0 \quad \forall i \in I
\]
of Meta-algorithm~\ref{alg:NearestPoint} cannot be satisfied, and must be replaced by the inequality
\[
  \langle y_n - z, x_i - y_n \rangle \ge - \eta \quad \forall i \in I
\]
for some small $\eta > 0$. Bearing in mind these changes we arrive at the following robust version of Meta-algorithm~\ref{alg:NearestPoint}
given in Algorithm~\ref{alg:NearestPointRobust}.

\begin{algorithm}[t]	\label{alg:NearestPointRobust}
\caption{Robust meta-algorithm for finding the nearest point in a polytope.}

\noindent\textbf{Input:} {a point $z \in \mathbb{R}^d$, a collection of points $\{ x_1, \ldots, x_{\ell} \} \subset
\mathbb{R}^d$, an algorithm $\mathcal{A}_{\varepsilon}$, $\varepsilon > 0$, for approximately solving the nearest point problem,
parameters $s, q \in \{ 1, \ldots, \ell \}$ with $s \ge q$, an $(s, q)$-exchange rule $\mathcal{E}$, and tolerance $eta > 0$.}

\noindent\textbf{Initialization:} {Put $n = 0$, choose an index set $I_n \subseteq I$ with $|I_n| = s$, and define
$P_n = \co\{ x_i \mid i \in I_n \}$.}

\noindent\textbf{Step 1.} {Compute $y_n^{\varepsilon} = \mathcal{A}_{\varepsilon}(z, P_n)$. If $y_n$ satisfies the optimality condition
\begin{equation} \label{eq:TrialPointApproxOptCond}
  \langle y_n^{\varepsilon} - z, x_i - y_n^{\varepsilon} \rangle \ge -\eta \quad \forall i \in I,
\end{equation}
\textbf{return} $y_n$.}

\noindent\textbf{Step 2.} {Compute $(J_{1n}, J_{2n}) = \mathcal{E}(I_n)$ and define
\[
  I_{n + 1} = \Big( I_n \setminus J_{1n} \Big) \cup J_{2n}, \quad
  P_{n + 1} = \co\big\{ x_i \bigm| i \in I_{n + 1} \big\}.
\]
Set $n = n + 1$ and go to \textbf{Step 1}.}
\end{algorithm}

Let us first show that under some natural assumptions Meta-algorithm~\ref{alg:NearestPointRobust} terminates after
a finite number of iterations and returns an approximate, in some sense, solution of the problem $(\mathcal{P})$.

\begin{lemma} \label{lem:ApproxDistDecay}
Suppose that Meta-algorithm~\ref{alg:NearestPointRobust} is correctly defined and the exchange rule $\mathcal{E}$ satisfies
the distance decay condition: if for some $n \in \mathbb{N}$ the point $y_n^{\varepsilon}$ does not satisfy
the approximate optimality condition \eqref{eq:TrialPointApproxOptCond}, then
\[
  \dist(z, P_{n + 1}) < \dist(z, P_n).
\]
Then Meta-algorithm~\ref{alg:NearestPointRobust} terminates after a finite number of steps and returns a point
$y_n^{\varepsilon}$ such that $\| y_n^{\varepsilon} - x_* \| < \sqrt{\eta}$, where $x_*$ is an optimal solution of the problem
$(\mathcal{P})$.
\end{lemma}

\begin{proof}
The fact that the algorithm terminates in a finite number of steps can be proved in exactly the same way as in the proof
of Lemma~\ref{lem:DistDecay}. Therefore, one only needs to check that the validity of the inequality \eqref{eq:TrialPointApproxOptCond}
implies that $\| y_n^{\varepsilon} - x_* \| < \sqrt{\eta}$.

Indeed, suppose that a point $u \in P$ satisfies the inequality
\begin{equation} \label{eq:ApproxOptPoint}
	\langle u - z, x_i - u \rangle \ge - \eta \quad \forall i \in I.
\end{equation}
Adding and subtracting $z$ one gets
\begin{align*}
	\| u - x_* \|^2 &= \| u - z \|^2 - 2 \langle u - z, x_* - z \rangle + \| x_* - z \|^2
	\\
	&= \Big( \| u - z \|^2 - \langle u - z, x_* - z \rangle \Big)
	+ \Big( \| x_* - z \|^2 - \langle u - z, x_* - z \rangle \Big)
	\\
	&= \langle u - z, u - x_* \rangle - \langle x_* - z, u - x_* \rangle.
\end{align*}
Hence applying Prp.~\ref{prp:NearPoint_OptCond} one obtains
\[
	\| u - x_* \|^2 \le \langle u - z, u - x_* \rangle.
\]
Since $x_* \in P$, there exists $\alpha_i \ge 0$, $i \in I$ such that
\[
	x_* = \sum_{i = 1}^{\ell} \alpha_i x_i, \quad \sum_{i = 1}^{\ell} \alpha_i = 1.
\]
Hence with the use of \eqref{eq:ApproxOptPoint} one finally gets that
\[
	\| u - x_* \|^2 \le \sum_{i = 1}^{\ell} \alpha_i \langle u - z, u - x_i \rangle
	\le \sum_{i = 1}^{\ell} \alpha_i \eta = \eta,
\]
which implies that $\| u - x_* \| \le \sqrt{\eta}$.
\end{proof}

Let us now check that under some natural assumptions on the algorithm $\mathcal{A}_{\varepsilon}$ the steepest
descent exchange rule with $y_n$ replaced by $y_n^{\varepsilon}$ satisfies the distance decay condition from the previous lemma.

\begin{theorem}
Let $\ell > d + 1$ and the algorithm $\mathcal{A}_{\varepsilon}$ satisfy the two following conditions:
\begin{enumerate}
\item{\textbf{Exchange correctness}: for any $n$ one has $y_n \in P_{n + 1}$ and there exists $j \in I_n$ such that
$\mathcal{A}_{\varepsilon}(z, P_n) \in \co\{ x_i \mid i \in I_n \setminus \{ j \} \}$;
}

\item{\textbf{Approximate optimality}: for any point $w$ and any polytope $Q \subset \mathbb{R}^d$ one has
$\mathcal{A}_{\varepsilon}(w, Q) \in Q$ and
\[
	\big\| \mathcal{A}_{\varepsilon}(w, Q) - u_* \big\| < \varepsilon,
\]
where $u_*$ is an optimal solution of the problem $\min_{u \in Q} \| u - w \|$.
}
\end{enumerate}
Suppose also that
\begin{equation} \label{eq:ParametersConsistency}
	\eta > \varepsilon D(P) + \varepsilon\big( \dist(z, P_0) + \varepsilon \big), \quad
	D(P) = \max_{i, j \in I} \| x_i - x_j \|.
\end{equation}
Then Meta-algorithm~\ref{alg:NearestPointRobust} with $s = d + 1$, $q = 1$, and the steepest descent exchange rule
with $y_n$ replaced by $y_n^{\varepsilon}$ is correctly defined, terminates after a finite number of iterations, and
returns a point $y_n^{\varepsilon}$ such that $\| y_n^{\varepsilon} - x_* \| < \sqrt{\eta}$, where $x_*$ is an optimal
solution of the problem $(\mathcal{P})$.
\end{theorem}

\begin{proof}
Clearly, the exchange correctness condition ensures that the steepest descent exchange rule and, therefore,
Meta-algorithm~\ref{alg:NearestPointRobust} itself, are correctly defined. Therefore, by Lemma~\ref{lem:ApproxDistDecay}
it is sufficient to prove that the distance decay condition holds true.

Suppose that for some $n \in \mathbb{N}$ the point $y_n^{\varepsilon}$ does not satisfy
the stopping criterion \eqref{eq:TrialPointApproxOptCond}. Let also $(i_{1n}, i_{2n}) \in I_n \times I$ be
the output of the steepest descent exchange rule. By definition of $i_{2n}$ and the fact that condition
\eqref{eq:TrialPointApproxOptCond} does not hold true one has
\begin{equation} \label{eq:NotApproxOptimal}
	\langle y_n^{\varepsilon} - z, x_{i_{2n}} - y_n^{\varepsilon} \rangle
	= \min_{i \in I} \langle y_n^{\varepsilon} - z, x_i - y_n^{\varepsilon} \rangle \le - \eta.
\end{equation}
Define $x_n(\alpha) = (1 - \alpha) y_n + \alpha x_{i_{2n}}$. Clearly, $x_n(\alpha) \in P_{n + 1}$ for all
$\alpha \in [0, 1]$, since $y_n \in P_{n + 1}$ by the exchange correctness condition and $x_{i_{2n}} \in P_{n + 1}$
by the definition of the exchange rule.

For any $\alpha \in \mathbb{R}$ one has
\[
  f(\alpha) := \| x_n(\alpha) - z \|^2 = (1 - \alpha)^2 \| y_n^{\varepsilon} - z \|^2 +
  2 \alpha (1 - \alpha) \langle y_n^{\varepsilon} - z, x_{i_{2n}} - z \rangle + \alpha^2 \| x_{i_{2n}} - z \|^2.
\]
Note that $f(0) = \| y_n - z \|^2 = \dist(z, P_n)^2$ and
\[
  f'(0) = - 2 \| y_n - z \|^2 + 2 \langle y_n - z, x_{i_{2n}} - z \rangle
  = \langle y_n - z, x_{i_{2n}} - y_n \rangle.
\]
Let us show that from inequality \eqref{eq:ParametersConsistency} it follows that $f'(0) < 0$. Then for any sufficiently small
$\alpha \in (0, 1)$ one has $f(\alpha) < f(0)$ and
\begin{align*}
  \dist(z, P_{n + 1}) &= \min_{x \in P_{n + 1}} \| z - x \| \le \| z - x_n(\alpha) \|
  \\
  &= \sqrt{f(\alpha)} < \sqrt{f(0)} = \| z - y_n \| = \dist(z, P_n),
\end{align*}
that is, the distance decay condition holds true.

Let us check the inequality $f'(0) < 0$. Indeed, adding and subtracting $y_n^{\varepsilon}$ twice one gets that
\begin{align*}
	f'(0) &= \langle y_n - z, x_{i_{2n}} - y_n \rangle
	= \langle y_n - y_n^{\varepsilon}, x_{i_{2n}} - y_n \rangle + \langle y_n^{\varepsilon} - z, x_{i_{2n}} - y_n \rangle
	\\
	&= \langle y_n - y_n^{\varepsilon}, x_{i_{2n}} - y_n \rangle + \langle y_n^{\varepsilon} - z, x_{i_{2n}} - y_n^{\varepsilon} \rangle
	+ \langle y_n^{\varepsilon} - z, y_n^{\varepsilon} - y_n \rangle.
\end{align*}
Consequently, by approximate optimality condition on $\mathcal{A}_{\varepsilon}$ and inequality \eqref{eq:NotApproxOptimal} one has
\[
	f'(0) \le \varepsilon \| x_{i_{2n}} - y_n \| - \eta + \| y_n^{\varepsilon} - z \| \varepsilon
	\le \varepsilon D(P) - \eta + \varepsilon\big( \dist(z, P_n) + \varepsilon \big).
\]
Hence with the use of \eqref{eq:ParametersConsistency} one gets that $f'(0) < 0$ for $n = 0$. Therefore $\dist(z, P_1) < \dist(z, P_0)$,
which implies that $f'(0) < 0$ for $n = 1$ and $\dist(z, P_2) < \dist(z, P_1) < \dist(z, P_0)$.
Now, arguing by induction one can check that $f'(0) < 0$ for any $n$ such that the stopping criterion \eqref{eq:TrialPointApproxOptCond}
does not hold true.
\end{proof}

Let us comment on the assumptions of the theorem above.
\begin{enumerate}
\item{One can replace the approximate optimality condition on the algorithm $\mathcal{A}_{\varepsilon}$ from the theorem
above by the following condition: for any point $w$ and any polytope $Q = \co\{ u_1, \ldots, u_m \} \subset \mathbb{R}^d$ one has
\[
	\big\langle \mathcal{A}_{\varepsilon}(w, Q) - z, u_i - \mathcal{A}_{\varepsilon}(w, Q) \big\rangle
	\ge - \varepsilon \quad \forall i \in \{ 1, \ldots, m \},
\]
that is, the algorithm $\mathcal{A}_{\varepsilon}$ returns a point approximately satisfying the optimality condition for
the problem $\min_{u \in Q} \| u - w \|$. If $u_*$ is an optimal solution of this problem, then as was shown in the proof of
Lemma~\ref{lem:ApproxDistDecay} one has $\| \mathcal{A}_{\varepsilon}(w, Q) - u_* \| \le \sqrt{\varepsilon}$. Consequently,
the theorem above remains to hold true in this case provided inequality \eqref{eq:ParametersConsistency} is replaced by the inequality
$\eta > \sqrt{\varepsilon} D(P) + \sqrt{\varepsilon}\big( \dist(z, P_0) + \sqrt{\varepsilon} \big)$.
}

\item{Roughly speaking, condition \eqref{eq:ParametersConsistency} means that to solve the problem $(\mathcal{P})$ with
a pre-specified accuracy with the use of Meta-algorithm~\ref{alg:NearestPointRobust} one must assume that the algorithm $\mathcal{A}$
solves the corresponding reduced nearest point subproblems $\min_{x \in P_n} \| x - z \|$ with a significantly higher accuracy.
In other words, condition \eqref{eq:ParametersConsistency} can be viewed as a rigorous formulation of an intuitively obvious fact that
Meta-algorithm~\ref{alg:NearestPointRobust} has lower accuracy than the algorithm $\mathcal{A}$ that is used as
a subroutine on each iterations.
}

\item{Finally, the exchange correctness condition on the algorithm $\mathcal{A}_{\varepsilon}$ can be interpreted as
the implicit assumption that the steepest descent exchange rule is able to remove (and actually removes) such a point $x_{i_{1n}}$
from $P_n$ that is not needed for computing the actual projection $y_n$ of $z$ onto $P_n$. In the case when the points $x_i$,
$i \in I_n$, are affinely independent this condition, in essence, means that the point $y_n^{\varepsilon}$ belongs to the same face
of $P_n$ as $y_n$. Once can also interpret the exchange correctness condition as the assumption that the algorithm $\mathcal{A}_{\varepsilon}$
returns a point having the same ``structure'', i.e. the same nonzero coefficients in the convex combination
\[
	y_n^{\varepsilon} = \sum_{i \in I_n} \alpha_i x_i, \quad \sum_{i \in I_n} \alpha_i = 1, \quad \alpha_i \ge 0,
\]
as the projection $y_n$ of $z$ onto $P_n$.
}
\end{enumerate}

\subsection{Numerical experiments}

To demonstrate the effectiveness of the proposed acceleration technique, it was verified numerically via multiple experiments with various
values of $d$ and $\ell$. For each choice of $d$ and $\ell$ we randomly generated 10 problems and an average computation time for these
problems was used as a performance measure for each method.

The problem data was chosen as follows. First, we randomly generated $\ell$ points $\{ u_1, \ldots, u_{\ell}$ in $\mathbb{R}^d$,
uniformly distributed over the $d$-dimensional cube $[-1, 1]^d$. Then these points were compressed and shifted as follows:
\[
	x_i = (1 + 0.01 u_{i1}, u_{i2}, \ldots, u_{id})^T \quad \forall i \in \{ 1, \ldots, \ell \}.
\]
The point $z$ was chosen as $z = 0$. As was noted in \cite{Wolfe} and demonstrated by our numerical experiments, this particular problem
is very challenging for methods for finding the nearest point in a polytope (especially in the case when $\ell$ is much greater than $d$)
both in terms of computation time and finding an optimal solution with high accuracy. For our numerical experiments we used the values
$d = 3$, $d = 10$, and $d = 50$, while values of $\ell$ was chosen from the set
$\{ 100, 200, 300, 500, 1000, 2000, 3000, 5000, 10000, 20000 \}$ and depended on $d$.

Without trying to conduct exhaustive numerical experiments, we tested our acceleration technique on only 3 classic methods: namely,
the MDM method \cite{MDM}, the Wolfe method \cite{Wolfe}, and a method based on quadratic programming. All methods were implemented in
\textsc{Matlab}. The last method was based on solving the problem
\[
	\min_{\alpha} \frac{1}{2} \Big\| \sum_{i = 1}^{\ell} \alpha_i x_i \Big\|^2 \quad \text{subject to} \quad
	\sum_{i = 1}^{\ell} \alpha_i = 1, \quad \alpha_i \ge 0, \quad i \in I.
\]
with the use of \texttt{quadprog}, the standard \textsc{Matlab} routine for solving quadratic programming problems. We used this routine
with default settings. The number of iterations of the MDM and Wolfe method was limited to $10^6$, and we used the inequalities
\[
	\delta(v_k) < \varepsilon, \quad \langle X, P_J \rangle > \langle X, X \rangle - \varepsilon
\]
with $\varepsilon = 10^{-4}$ as termination criteria for the MDM method and the Wolfe method respectively (see the descriptions of these
methods in \cite{MDM,Wolfe}). The value $10^{-4}$ was used, since occasionally both methods failed to terminate with smaller value of
$\varepsilon$ for large $\ell$ (this statement was especially true for the MDM method).

Finally, we implemented each method ``on its own'' and also incorporated each method within the robust acceleration technique, i.e.
Meta-algorithm~\ref{alg:NearestPointRobust}, with the steepest descent exchange rule given in Algorithm~\ref{alg:SteepestDescentExchRule},
and Algorithm~\ref{alg:IndexReplSearch} for searching the index of a point to be removed from the polytope $P_n$ in the case when
all coefficients in the expression of $y_n^{\varepsilon}$ as the convex combination of the points $x_i$, $i \in I_n$, were positive (note,
however, that such situation never occurred in our numerical experiments).

The number of iterations of the acceleration technique was limited to $\ell$. We used the value $\eta = 10^{-4}$ for $d = 3$ and $d = 10$
and $\eta = 5 \cdot 10^{-4}$ for $d = 50$, since the acceleration technique occasionally failed to terminate for $\eta = 10^{-4}$ in this
case. In addition, we terminated the computations if the computation time exceeded 1 minute.

\begin{figure}[t] \label{fig:NearestPoint3}
\caption{The results of numerical experiments in the case $d = 3$ for \texttt{quadprog} routine (left figure) and
the MDM method (right figure).}
\includegraphics[width=0.42\textwidth]{quadprog_3}
\includegraphics[width=0.42\textwidth]{MDM_3}
\end{figure}

\begin{figure}[t] \label{fig:NearestPoint10}
\includegraphics[width=0.42\textwidth]{quadprog_10}
\includegraphics[width=0.42\textwidth]{MDM_10}
\caption{The results of numerical experiments in the case $d = 10$ for \texttt{quadprog} routine (left figure) and
the MDM method (right figure).}
\end{figure}

\begin{figure}[t] \label{fig:NearestPoint50}
\includegraphics[width=0.42\textwidth]{quadprog_50}
\includegraphics[width=0.42\textwidth]{Wolfe_50}
\caption{The results of numerical experiments in the case $d = 50$ for \texttt{quadprog} routine (left figure) and
the Wolfe method (right figure).}
\end{figure}

The results of numerical experiments are given in Fig.~\ref{fig:NearestPoint3}--\ref{fig:NearestPoint50}. Let us note that both the MDM
method and its accelerated version were extremely slow in comparison with other methods in the case $d = 50$. That is why we do not present
the corresponding results of the numerical experiments. Furthermore, the results of numerical experiments demonstrated that the difference
in performance between the Wolfe method and its accelerated version significantly increases with the growth of $d$. Since the results were
quantitatively the same for all $d$, here we present them only in the case $d = 50$, which allows one to more clearly assess the efficiency
of the acceleration technique.

The numerical experiments clearly demonstrate that the proposed acceleration technique with the steepest descent exchange rule allows one
to significantly reduce the computation time for method of finding the nearest point in a polytope. In the case of \texttt{quadprog} routine
and the MDM method the reduction in time is proportional to the number of points $\ell$. Furthermore, numerical experiments also show that
for the accelerated versions of these methods the computation time increases \textit{linearly} with dimension for the problem under consideration
(but we do not claim the linear in $\ell$ complexity of the acceleration technique for all types of problem data).

In the case of the Wolfe method the situation is more complicated. For relatively small values of $\ell$ the ``pure'' Wolfe method outperforms
its accelerated version, but for larger $\ell$ the accelerated version is faster than the pure method and the reduction of computation time
increases with the increase of $\ell$.

\begin{remark}
In our implementation of the acceleration technique (Meta-algo\-rithm~\ref{alg:NearestPointRobust}) the algorithm $\mathcal{A}$
was applied afresh on each iteration (i.e. without using any information from the previous iteration). It should be noted that, in particular,
the performance of the accelerated version of the Wolfe method can be significantly improved, if one uses the \textit{corral}
(see \cite{Wolfe}) obtained on the previous iteration as the initial guess for the next iteration (a similar remark is true for accelerated
versions of other methods). However, since our main goal was to demonstrate the performance of the acceleration technique on its own,
here we do not discuss potential ways this technique can be efficiently integrated with a particular method for finding the
nearest point in a polytope and do not present any results of numerical experiments for such fully integrated accelerated methods.
\end{remark}

\section{Computing the distance between two polytopes}

The acceleration technique for method for finding the nearest point in a polytope studied in the previous section admits a straightforward
extension to the case of problem of computing the distance between two convex polytopes defined as the convex hulls of finite sets
of points. This section is devoted to a detailed discussion of this extensions, that is, to a detailed discussion of an acceleration
technique for method of solving the following optimization problem
\[
  \min_{(x, y)} \: \| x - y \| \enspace \text{ subject to } \enspace
  x \in P := \co\big\{ x_1, \ldots, x_{\ell} \big\}, \quad
  y \in Q = \co\big\{ y_1, \ldots, y_m \}.
  \qquad \eqno{(\mathcal{D})}
\]
in the case when $\ell \gg d$ and $m \gg d$. Here $x_i \in \mathbb{R}^d$, $i \in I = \{ 1, \ldots, \ell \}$, and
$y_j \in \mathbb{R}^d$, $j \in J = \{ 1, \ldots, m \}$, are given points.

\subsection{An extension of the acceleration technique}

Before we proceed to the description of an acceleration technique, let us present convenient optimality conditions
for the problem $(\mathcal{D})$. These conditions are well-known, but we include their short proofs for the sake of completeness.
Denote by $Pr_A(x)$ the Euclidean projection of a point $x \in \mathbb{R}^d$ on the closed convex set $A$, i.e. an optimal
solution of the problem $\min_{y \in A} \| x - y \|$.

\begin{proposition} \label{prp:PolytopesDistOptCond}
A pair $(x_*, y_*) \in P \times Q$ is an optimal solution of the problem $(\mathcal{D})$ if and only if
$Pr_P(y_*) = x_*$ and $Pr_Q(x_*) = y_*$ or, equivalently,
\begin{equation} \label{eq:SeparateOptimality}
	\langle x_* - y_*, x_i - x_* \rangle \ge 0 \quad \forall i \in I, \quad
	\langle y_* - x_*, y_j - y_* \rangle \ge 0 \quad \forall j \in J.
\end{equation}
\end{proposition}

\begin{proof}
The fact that conditions \eqref{eq:SeparateOptimality} are equivalent to the equalities $Pr_P(y_*) = x_*$ and $Pr_Q(x_*) = y_*$
follows directly from Prp.~\ref{prp:NearPoint_OptCond}. Let us check that these conditions are equivalent to the optimality of
$(x_*, y_*)$.

Indeed, as is easily seen, inequalities \eqref{eq:SeparateOptimality} are satisfied if and only if
\[
	\langle x_* - y_*, x - x_* \rangle \ge 0 \quad \forall x \in P, \quad
	\langle y_* - x_*, y - y_* \rangle \ge 0 \quad \forall y \in Q.
\]
In turn, these conditions are satisfied if and only if
\[
	\langle x_* - y_*, x - x_* \rangle + \langle y_* - x_*, y - y_* \rangle \ge 0
	\quad \forall x \in P, \: y \in Q.
\]
It remains to note that by standard optimality conditions for convex programming problems the above inequality
is equivalent to the optimality of $(x_*, y_*)$.
\end{proof}

Suppose that an algorithm $\mathcal{A}$ for solving the problem $(\mathcal{D})$ is given. For any two polytopes
$U, W \subset \mathbb{R}^d$, defined as the convex hulls of finite collections of points, it returns optimal solution
$(x_*, y_*) = \mathcal{A}(U, W)$ of the problem
\[
	\min_{(x, y)} \: \| x - y \| \quad \text{subject to} \quad x \in P, \quad y \in Q.
\]
We propose to accelerate this method in precisely the same way as the method for solving the nearest point problem.
Namely, one chooses ``small'' subpolytopes $P_0 \subset P$ and $Q_0 \subset Q$ and applies the algorithm $\mathcal{A}$
to find the distance between these subpolytopes. If an optimal solution of this problem coincides with an optimal solution
of the problem $(\mathcal{D})$ (this fact is verified via the optimality conditions from Prp.~\ref{prp:PolytopesDistOptCond}), then
the computations are terminated. Otherwise, one shifts these polytopes and repeats the same procedure till an optimal solution
of the distance problem for subpolytopes $P_0$ and $Q_0$ coincides with an optimal solution of the problem $(\mathcal{D})$. A general
structure of this acceleration technique (meta-algorithm) is essentially the same as the structure of Meta-algorithm~\ref{alg:NearestPoint}
and is given in Algorithm~\ref{alg:PolytopesDistance}.

\begin{algorithm}[t]	\label{alg:PolytopesDistance}
\caption{Meta-algorithm for finding the distance between two polytopes.}

\noindent\textbf{Input:} {two collection of points $\{ x_1, \ldots, x_{\ell} \} \subset
\mathbb{R}^d$ and $\{ y_1, \ldots, y_m \} \subset \mathbb{R}^d$, an algorithm $\mathcal{A}$ for solving the problem of computing
the distance between two polytopes, parameters $s, q \in \{ 1, \ldots, \ell \}$ with $s \ge q$, and
an $(s, q)$-exchange rule $\mathcal{E}$.}

\noindent\textbf{Initialization:} {Put $n = 0$, choose index sets $I_n \subseteq I$ and $J_n \subseteq J$ with
$|I_n| = |J_n| = s$. Define $P_n = \co\{ x_i \mid i \in I_n \}$ and $Q_n = \co\{ y_j \mid j \in J_n \}$.}

\noindent\textbf{Step 1.} {Compute $(x_n, y_n) = \mathcal{A}(P_n, Q_n)$. Compute
\[
	\theta_{xn} = \min_{i \in I} \langle x_* - y_*, x_i - x_* \rangle, \quad
	\theta_{yn} = \min_{j \in J} \langle y_* - x_*, y_j - y_* \rangle.
\]
If $\theta_{xn} \ge 0$ and $\theta_{yn} \ge 0$, \textbf{return} $(x_n, y_n)$.}

\noindent\textbf{Step 2.} {If $\theta_{xn} < 0$, compute $(I_{1n}, I_{2n}) = \mathcal{E}(I_n)$ and define
\[
  I_{n + 1} = \Big( I_n \setminus J_{1n} \Big) \cup J_{2n}, \quad
  P_{n + 1} = \co\big\{ x_i \bigm| i \in I_{n + 1} \big\}.
\]
If $\theta_{yn} < 0$, compute $(J_{1n}, J_{2n}) = \mathcal{E}(J_n)$ and define
\[
  J_{n + 1} = \Big( J_n \setminus J_{1n} \Big) \cup J_{2n}, \quad
  Q_{n + 1} = \co\big\{ y_j \bigm| j \in J_{n + 1} \big\}.
\]
Set $n = n + 1$ and go to \textbf{Step 1}.}
\end{algorithm}

Arguing in essentially the same way as in the proofs of Lemmas~\ref{lem:DistDecay} and \ref{lem:SubpolytopeSize} one
can verify that the following results hold true. These results, just like the aforementioned lemmas, can be viewed
as convenient criteria for choosing an efficient exchange rule for Algorithm~\ref{alg:PolytopesDistance}.

For any sets $A, B \subset \mathbb{R}^d$ denote by
\[
	\dist(A, B) = \inf\big\{ \| x - y \| \bigm| x \in A, \: y \in B \big\}
\]
the Euclidean distance between these sets.

\begin{lemma}
Suppose that the exchange rule $\mathcal{E}$ satisfies the distance decay condition for the problem $(\mathcal{D}$:
if for some $n \in \mathbb{N}$ the pair $(x_n, y_n)$ does not satisfy optimality condition $\theta_{xn} \ge 0$ and
$\theta_{yn} \ge 0$, then
\begin{equation} \label{eq:DistDecayCondition}
  \dist(P_{n + 1}, Q_{n + 1}) < \dist(P_n, Q_n).
\end{equation}
Then Algorithm~\ref{alg:PolytopesDistance} terminates after a finite number of steps and returns an optimal solution of the
problem $(\mathcal{D})$.
\end{lemma}

\begin{lemma}
Let an $(s, q)$-exchange rule $\mathcal{E}$ with $s \ge q$ satisfy the distance decay condition for any polytopes
$U, W \subset \mathbb{R}^d$. Then $s \ge d + 1$.
\end{lemma}

The following theorem, whose prove is similar to the proof of Theorem~\ref{thm:StDescExchRule}, shows that the steepest
descent exchange rule, proposed in Subsection~\ref{subsect:ExchangeRule}, satisfies the distance decay condition for
the problem $(\mathcal{D})$.

\begin{theorem}
For any polytopes $P = \co\{ x_1, \ldots, x_{\ell} \}$ and $Q = co\{ y_1, \ldots, y_m \}$ with $\ell \ge d + 1$ and
$m \ge d + 1$ the steepest descent exchange rule satisfies the distance decay condition for the problem $(\mathcal{D})$.
\end{theorem}

Let us finally note that the results of numerical experiments, that will be reported in later versions of this article,
demonstrate that Meta-algorithm~\ref{alg:PolytopesDistance} with the steepest descent exchange rule allows one to significantly
accelerate existing methods for computing the distance between two polytopes in the case when $\ell \gg d$ and $m \gg d$.

\bibliographystyle{abbrv}