#  Finding best approximation pairs for two intersections of closed convex
sets

Heinz H. Bauschke, Shambhavi Singh, and Xianfu Wang  Mathematics, University
of British Columbia, Kelowna, B.C. V1V 1V7, Canada. E-mail:
<EMAIL_ADDRESS>Mathematics, University of British Columbia, Kelowna,
B.C. V1V 1V7, Canada. E-mail<EMAIL_ADDRESS>Mathematics, University
of British Columbia, Kelowna, B.C. V1V 1V7, Canada. E-mail<EMAIL_ADDRESS>
(October 15, 2021)

###### Abstract

The problem of finding a best approximation pair of two sets, which in turn
generalizes the well known convex feasibility problem, has a long history that
dates back to work by Cheney and Goldstein in 1959.

In 2018, Aharoni, Censor, and Jiang revisited this problem and proposed an
algorithm that can be used when the two sets are finite intersections of
halfspaces. Motivated by their work, we present alternative algorithms that
utilize projection and proximity operators. Our modeling framework is able to
accommodate even convex sets. Numerical experiments indicate that these
methods are competitive and sometimes superior to the one proposed by Aharoni
et al.

2020 Mathematics Subject Classification: Primary 65K05; Secondary 47H09,
90C25.

Keywords: Aharoni–Censor–Jiang algorithm, best approximation pair,
Douglas–Rachford algorithm, dual-based proximal method, proximal distance
algorithm, stochastic subgradient descent.

##  1 Introduction

Throughout this paper, we assume that

$Y$ is a finite-dimensional real Hilbert space with inner product
$\left\langle{\cdot},{\cdot}\right\rangle\colon Y\times Y\to\mathbb{R}$,

and induced norm $\|\cdot\|$. Let $m\in\\{1,2,\ldots\\}$, set
$I:=\\{1,\ldots,m\\}$ and suppose that

$(\forall i\in I)\;\;\text{$A_{i}$ and $B_{i}$ are nonempty closed convex
subsets of $Y$}$

such that

$A:=\bigcap_{i\in I}A_{i}\neq\varnothing\;\;\text{and}\;\;B:=\bigcap_{i\in
I}B_{i}\neq\varnothing.$

(We assume here without loss of generality that there are as many sets $A_{i}$
as $B_{j}$; otherwise, we can either “copy” sets or use the full space $Y$
itself.) It will occasionally be convenient to work with the convention
$A_{m+1}=A_{1}$, $A_{m+2}=A_{2}$, etc.; or, more formally,
$A_{n}=A_{1+\mathrm{rem}(n-1,m)}$ and $B_{n}=B_{1+\mathrm{rem}(n-1,m)}$. We
also assume that the projection operators $P_{A_{i}}$ and $P_{B_{i}}$ are
“easy” to compute while the projections $P_{A}$ and $P_{B}$ are “hard” and not
readily available (unless $m=1$). The problem we are interested in is to find
a _best approximation pair_ , i.e., to

$\text{Find $(\bar{a},\bar{b})\in A\times B$ such
that~{}}\|\bar{a}-\bar{b}\|=\inf_{(a,b)\in A\times B}\|a-b\|.$ (1)

(Note that this problem is actually a generalization of the famous _convex
feasibility problem_ which asks to find a point in $A\cap B$ provided that
this intersection is nonempty which we do not assume here!) This problem has a
long history, and the first systematic study was given by Cheney and Goldstein
in 1959 [12]; see also, e.g., [3], [4], and [6]. These works, however, assume
that the projection operators $P_{A}$ and $P_{B}$ are explicitly available,
which essentially means that $m=1$. Recently, Aharoni, Censor, and Jiang (see
[1]) tackled the general case. Indeed, assuming that the sets $A_{i}$ and
$B_{i}$ are _halfspaces_ , they presented a new algorithm — which we call ACJ
for simplicity — for solving 1 where they do not require knowledge of the
projectors $P_{A}$ and $P_{B}$ onto the corresponding polyhedra $A$ and $B$.

_The purpose of this paper is to provide other approaches to solving 1. We
also provide the required proximity operators as well as numerical
comparisons. The algorithms considered will rely only on the operators
$P_{A_{i}}$ and $P_{B_{i}}$ and some other operators that are available in
closed form._

The algorithms presented will work for general closed convex sets, not just
polyhedra as long as the projection operators onto the individual sets making
up the intersections are available. Implementable formulae for the underlying
algorithmic operators are provided. Numerical experiments, similar to one in
Aharoni et al.’s paper, are also performed. Our results show that other
algorithms should be seriously considered for solving 1, especially if $m$ is
small.

The remainder of the paper is organized as follows. In Section 2, we consider
reformulations of 1 that are more amenable to the algorithms discussed in
Section 3. These algorithms rely on computable formulae which we present in
Section 4. We present a small example on which these algorithms are run and
convergence is observed using different metrics in Section 5. In Section 5.1
and Section 5.2, we consider examples where the solutions are known or
unknown, respectively. We also discuss the (positive) effect of pairing up
constraints (see Section 5.3). We conclude the paper in Section 6 with a brief
summary of our findings.

The notation employed in this paper is fairly standard and follows largely [5]
to which we also refer the reader for general background material. For the
reader’s convenience, let us review some notions from convex analysis that are
of fundamental importance to this paper. The _indicator function_ of a set $C$
is written as $\iota_{C}$; we have $\iota_{C}(x)=0$ if $x\in C$ and
$\iota_{C}(x)=+\infty$ otherwise. The corresponding _distance function_ is
$d_{C}(x)=\inf_{c\in C}\|x-c\|$. If $C$ is convex, closed, and nonempty, then
for every $x$, there exists a unique point $P_{C}(x)\in C$ such that
$d_{C}(x)=\|x-P_{C}(x)\|$. The corresponding operator $P_{C}$ is called the
_projection operator_ or simply _projector_ of $C$. For instance, if $A$ is
the _halfspace_ $A=\big{\\{}{x\in
Y}~{}\big{|}~{}{\left\langle{a},{x}\right\rangle\leq\alpha}\big{\\}}$, where
$a\neq 0$ and $\alpha\in\mathbb{R}$, then

$P_{A}(x)=x-\frac{\max\\{0,\left\langle{a},{x}\right\rangle-\alpha\\}}{\|a\|^{2}}a.$
(2)

More generally, if $f$ is a function that is convex, lower semicontinuous, and
proper, then for every $x$, the function $y\mapsto
f(y)+\tfrac{1}{2}\|x-y\|^{2}$ has a unique minimizer which is the celebrated
_proximal mapping_ or _prox operator_ of $f$, written
$\operatorname{Prox}_{f}$. Note that if $f=\iota_{C}$, then we recover the
projection operator $P_{C}$: $\operatorname{Prox}_{\iota_{C}}=P_{C}$. A vector
$y$ is a _subgradient_ of $f$ at $x$ if for every $h$, we have
$f(x)+\left\langle{y},{h}\right\rangle\leq f(x+h)$; the set of all
subgradients at $x$ is the _subdifferential_ of $f$ at $x$, written $\partial
f(x)$.

##  2 Modeling 1

In the product Hilbert space

$X:=Y\times Y,$

we define nonempty closed convex sets by

$(\forall i\in I)\;\;C_{i}:=A_{i}\times B_{i},$

along with their intersection

$(\forall i\in I)\;\;C:=\bigcap_{i\in I}C_{i}=A\times B.$

Note that the projector onto $C_{i}$ is still easy to compute; indeed,
$P_{C_{i}}\colon(x,y)\mapsto(P_{A_{i}}x,P_{B_{i}}y)$. The problem 1 is thus
equivalent to

minimize $h(x,y)$ subject to $(x,y)\in C$, (3)

where

$h(x,y)=\alpha\|x-y\|^{p},$ (4)

$\alpha>0$, and $p\geq 1$. Note that $h$ has full domain and it is convex —
but $h$ is _not_ strictly convex because it has many minimizers:
$\big{\\{}{(x,x)}~{}\big{|}~{}{x\in Y}\big{\\}}$. Also note that if $p>1$,
then $h$ is differentiable.

The problem 1 can thus also be alternatively thought of as

minimize $h(x,y)+\sum_{i\in I}\iota_{C_{i}}(x,y)$, (5)

which features a nonsmooth objective function. In the next section, we survey
various algorithms that could be used to solve the problem 1 or its
reformulations. We also consider the case when 5 is approximated by

minimize $h(x,y)+\sum_{i\in I}Ld_{C_{i}}(x,y)$, (6)

for some “large” constant $L$.

##  3 Algorithms for solving 1

In this section, we discuss various algorithms. The algorithms in Section 3.1
and Section 3.2 are able to solve the original problem exactly while those in
the remaining subsections can be viewed as attempting to solve a perturbed
problem.

###  3.1 The Aharoni–Censor–Jiang Algorithm

This algorithm was recently proposed by Aharoni, Censor, and Jiang in [1]. We
denote their algorithm as ACJ. ACJ builds on the earlier HLWB algorithm. (The
letters in HLWB signify relevant works by Halpern [15], by Lions [20], by
Wittmann [23], and by Bauschke [2]; the name HLWB was coined by Censor in
[11].) ACJ can be viewed as an alternating version of HLWB to find a solution
of 1. Here is the description of ACJ. First, we fix a sequence
$(\lambda_{k})_{k\in{\mathbb{N}}}$ of positive real numbers such that

$\lambda_{k}\to
0,\quad\sum_{k\in{\mathbb{N}}}\lambda_{k}=\infty,\quad\sum_{k\in{\mathbb{N}}}|\lambda_{k}-\lambda_{k+m}|<\infty$
(7)

and also an increasing (not necessarily strictly though) sequence of natural
numbers $(n_{k})_{k\in{\mathbb{N}}}$ such that

$n_{k}\to\infty\quad\text{and}\quad\sup_{k_{0}\in\mathbb{N}}\sum_{k>k_{0}}\prod_{n>n_{k_{0}}}^{n_{k}}(1-\lambda_{n})<\infty.$
(8)

For instance, 7–8 hold when $(\forall{k\in{\mathbb{N}}})$
$\lambda_{k}=\frac{1}{k+1}$ and $n_{k}=\lfloor 1.1^{k}\rfloor$ (see [1, page
512]). Next, given our sequence of sets $(A_{i})_{i\in\mathbb{N}}$ and
${n\in{\mathbb{N}}}$, we define the operator

$Q_{A,n}\colon Y\times Y\to Y\colon(w,w^{\prime})\mapsto w_{n},$ (9)

where $w_{0}=w^{\prime}$ and $w_{n}$ is computed iteratively via

$(\forall i\in\\{0,1,\ldots,n-1\\})\quad
w_{i+1}=\lambda_{i+1}w+(1-\lambda_{i+1})P_{A_{i+1}}(w_{i}).$ (10)

The operator $Q_{B,n}$, for $(B_{i})_{i\in\mathbb{N}}$ and
${n\in{\mathbb{N}}}$, is defined analogously. Finally, we initialize
$(x_{0},y_{0})\in X\times X$, and iteratively update via

$\displaystyle(\forall{k\in{\mathbb{N}}})\quad(x_{k+1},y_{k+1})$
$\displaystyle:=\begin{cases}\big{(}Q_{B,n_{k}}(y_{k},y_{k}^{\prime}),y_{k}\big{)},&\text{if
$k$ is odd;}\\\
\big{(}x_{k},Q_{A,n_{k}}(x_{k},x_{k}^{\prime})\big{)},&\text{if $k$ is
even,}\end{cases}$ (11)

and where $(x^{\prime}_{k},y^{\prime}_{k})_{k\in{\mathbb{N}}}$ in $X\times X$
is a bounded sequence that can either be fixed beforehand, e.g.,
$(x^{\prime}_{k},y^{\prime}_{k})=(y_{0},x_{0})$, or dynamically updated using,
e.g., $(x^{\prime}_{0},y^{\prime}_{0})=(y_{0},x_{0})$ and
$(x^{\prime}_{k},y^{\prime}_{k})=(y_{k-1},x_{k-1})$ for
$k\in\\{1,2,\ldots\\}$.

The main result of [1] yields the convergence of
$(x_{k},y_{k})_{k\in{\mathbb{N}}}$ to a solution of 1 provided that
$C\neq\emptyset$ and each $A_{i}$ and $B_{i}$ is a _halfspace_. We refer the
reader to 2 for the formula for the projection onto a halfspace which is
required by ACJ.

######  Remark 3.1.

Note that ACJ takes into account the order of the sets while the problem 1
does not. It is a nice feature of ACJ that it works throughout in the “small”
space $X\times X$. On the other hand, we are not aware of any extension of ACJ
to the case when the sets underlying the intersections are not halfspaces.
This is an interesting topic for further research.

###  3.2 Douglas–Rachford Algorithm

This algorithm, abbreviated as DR, can be traced back to the paper by Douglas
and Rachford [14] although its relevance to optimization was brought to light
later in the seminal paper by Lions and Mercier [21]. DR can deal with
problems of the form 5, and it implicitly operates in the space $X^{m+1}$.
First, set $f_{0}(x,y):=\alpha\|x-y\|^{p}$ with $\alpha>0$ and $p\geq 1$, as
well as $f_{1}:=\iota_{C_{1}},\ldots,f_{m}:=\iota_{C_{m}}$ and
$I_{0}:=\\{0\\}\cup I$. Second, fix a parameter $0<\lambda<2$ (the default
being $\lambda=1$).

Now initialize
$\mathbf{z}_{0}:=(z_{0,0},z_{0,1},\ldots,z_{0,m})=\big{(}(x_{0,0},y_{0,0}),(x_{0,1},y_{0,1}),\ldots,(x_{0,m},y_{0,m})\big{)}\in
X^{m+1}$. Given $\mathbf{z}_{k}=(z_{k,0},z_{k,1},\ldots,z_{k,m})\in X^{m+1}$,
set

$\displaystyle\bar{z}_{k}$ $\displaystyle:=\frac{1}{m+1}\sum_{i\in
I_{0}}z_{k,i}$ (12a) $\displaystyle(\forall i\in I_{0})\;\;x_{k,i}$
$\displaystyle:=\operatorname{Prox}_{f_{i}}(2\bar{z}_{k}-z_{k,i})$ (12b)
$\displaystyle(\forall i\in I_{0})\;\;z_{k+1,i}$
$\displaystyle:=z_{k,i}+\lambda(x_{k,i}-\bar{z}_{k})$ (12c)

to obtain the update
$\mathbf{z}_{k+1}:=(z_{k+1,0},z_{k+1,1},\ldots,z_{k+1,m})$.

If $i\in I$, then the prox operators corresponding to $f_{i}$ is simply the
projector $P_{C_{i}}$. In particular, if each $C_{i}=A_{i}\times B_{i}$ is the
Cartesian product of two halfspaces, then we may utilize 2 twice to compute
$P_{C_{i}}=\operatorname{Prox}_{f_{i}}$. The prox operator
$\operatorname{Prox}_{f_{0}}$ will be computed in closed form for
$p\in\\{1,2\\}$ in Section 4.1 below. When $p=1$, which produced better
numerical results, then $f_{0}(x,y)=\alpha\|x-y\|$ and

$\displaystyle\operatorname{Prox}_{f_{0}}(x,y)=(x,y)-\frac{1}{\max\big{\\{}2,\|x-y\|/\alpha\big{\\}}}\big{(}x-y,y-x\big{)}.$
(13)

It is well known (see, e.g., [5, Proposition 28.7]) that the sequence
$(\bar{z}_{k})_{k\in\mathbb{N}}$ will converge to a solution of 5, i.e., of 1.

######  Remark 3.2.

The DR approach does not care about the order of the sets presented — unlike,
ACJ! A downside is that it operates in the larger space $X^{m+1}$ which can
become an issue if $m$ is large. On the positive side, if $C_{1}\cap\cdots\cap
C_{m}=\varnothing$, then $(\bar{z}_{k})_{k\in{\mathbb{N}}}$ will converge to a
minimizer of $f_{0}$ over the set of least-squares solutions (see [7,
Corollary 6.8] for further information). Finally, it does not require the
constraint sets $C_{i}$ to be Cartesian product of halfspaces.

###  3.3 Dual-Based Proximal Method

We largely follow Beck’s [8, Section 12.4.2] (see also [9] and [13] for
further background material) but slightly modify the algorithms presented
there to give two additional methods for solving 1. We will work with the form
given in 5 where $h(x,y)$ needs to be $\varepsilon$-strongly convex, for some
$\varepsilon>0$, which precludes using $\alpha\|x-y\|^{p}$ directly. However,
below we will add $\varepsilon\frac{1}{2}(\|x\|^{2}+\|y\|^{2})$ to this last
function to obtain the required $\varepsilon$-strong convexity. We point out
that by adding this energy term and solving the corresponding new perturbed
optimization problem, the solution obtained does not solve the original
problem exactly.

The first method considered is the Dual Proximal Gradient method, which —
following [8] — we abbreviated as DPG. Because the algorithm requires strong
convexity of the objective function, we consider

$f_{0}(x,y):=\alpha\tfrac{1}{2}\|x-y\|^{2}+\varepsilon\tfrac{1}{2}\big{(}\|x\|^{2}+\|y\|^{2}\big{)}\;\;\text{with
$\alpha>0$ and $\varepsilon>0$}.$ (14)

We also set $f_{1}:=\iota_{C_{1}},\ldots,f_{m}:=\iota_{C_{m}}$. Second, fix a
parameter $L\geq m/\varepsilon$.

Now initialize
$\mathbf{z}_{0}:=(z_{0,1},\dots,z_{0,m})=\big{(}(x_{0,1},y_{0,1}),\ldots,(x_{0,m},y_{0,m})\big{)}\in
X^{m}$, and update it using

$\displaystyle s_{k}$ $\displaystyle:=\sum_{i\in I}z_{k,i}$ (15a)
$\displaystyle x_{k}$ $\displaystyle:=\operatorname*{argmax}_{w\in
X}\big{[}\left\langle{w},{s_{k}}\right\rangle-f_{0}(w)\,\big{]}$ (15b)
$\displaystyle(\forall i\in I)\;\;z_{k+1,i}$
$\displaystyle:=z_{k,i}-\frac{1}{L}x_{k}+\frac{1}{L}P_{f_{i}}(x_{k}-Lz_{k,i})$
(15c)

to obtain $\mathbf{z}_{k+1}:=(z_{k+1,1},\dots,z_{k+1,m})$. This is the primal
representation of DPG, see [8, page 356], which is most convenient in our
setting. Once again, the prox operators corresponding to $f_{i}$ for $i\in I$
are just the projectors $P_{C_{i}}$. If the sets $C_{i}$ are Cartesian
products of halfspaces, we may use 2 to compute $P_{C_{i}}$. The closed form
for the argmax operator in 15b is given by

$\displaystyle
x_{k}=\frac{1}{(2\alpha+\varepsilon)\varepsilon}\big{(}(\alpha+\varepsilon)u_{k}+\alpha
v_{k},(\alpha+\varepsilon)v_{k}+\alpha u_{k}\big{)},\quad\text{where
$s_{k}=(u_{k},v_{k})$.}$ (16)

This formula will be proved in Section 4.2 below. For sufficiently small
$\varepsilon>0$, the primal sequence $(x_{k})_{k\in{\mathbb{N}}}$ approximates
a solution of 5 and hence of 1 provided that the relative interiors of the
sets $C_{i}$ form a nonempty intersection (see [8, page 362]). Note that we do
not expect that the primal sequence converges to an exact solution of 1
because the objective function $f_{0}$ in 14 is not identical to the one
required to tackle 1.

An accelerated version of DPG, known as Fast Dual Proximal Gradient or simply
FDPG, applies a FISTA-type acceleration (see [8, Section 12.3] and [9].) Here
is how FDPG proceeds: Starting with $\mathbf{z}_{0}$ as before, initialize
$\textbf{w}_{0}:=\mathbf{z}_{0}$, $t_{0}:=1$, and update via

$\displaystyle s_{k}^{\prime}$ $\displaystyle:=\sum_{i\in I}w_{k,i}$ (17a)
$\displaystyle u_{k}$ $\displaystyle:=\operatorname*{argmax}_{v\in
X}\big{[}\left\langle{v},{s_{k}^{\prime}}\right\rangle-f_{0}(v)\,\big{]}$
(17b) $\displaystyle(\forall i\in I)\;\;z_{k+1,i}$
$\displaystyle:=w_{k,i}-\frac{1}{L}u_{k}+\frac{1}{L}P_{f_{i}}(u_{k}-Lw_{k,i})$
(17c) $\displaystyle t_{k+1}$ $\displaystyle:=\frac{1+\sqrt{1+4t_{k}^{2}}}{2}$
(17d) $\displaystyle(\forall i\in I)\;\;w_{k+1,i}$
$\displaystyle:=z_{k,i+1}+\left(\frac{t_{k}-1}{t_{k+1}}\right)(z_{k,i+1}-z_{k,i})$
(17e)

to get the primal sequence of interest

$x_{k+1}:=\operatorname*{argmax}_{v\in
X}\big{[}\left\langle{v},{s_{k+1}}\right\rangle-
f_{0}(v)\,\big{]},\quad\text{where}\quad s_{k+1}:=\sum_{i\in I}z_{k+1,i}.$
(18)

Again, for sufficiently small $\varepsilon>0$, the sequence
$(x_{k})_{k\in{\mathbb{N}}}$ approximates a solution “close” to that of 5 —
but not exactly — and (as a consequence of 1) provided that the relative
interior of $C$ is not empty.

######  Remark 3.3.

Note that although a smaller $\varepsilon$ ensures a solution that is closer
to that of the original problem 1, it also increases the lower bound for $L$,
which in turn reduces the step size for each iteration as seen in 15c and 17c,
and so the speed of convergence reduces as well. These algorithms are not
affected by the order of the sets presented, but like DR, they operate in a
“large” space (here $X^{m}$). This may become a problem when $m$ is large.

###  3.4 Proximal Distance Algorithm

The Proximal Distance Algorithm, or PDA for short, was first introduced by
Lange and Keys [19]. It is motivated by the framework of MM algorithms, where
MM stands for majorize/minimize or for minorize/maximize depending on the
underlying problem. This framework was pioneered by Lange; see, e.g., his book
[18] on this topic. It can be interpreted as a prox-gradient method applied to
the function $\tfrac{1}{\rho}h+\tfrac{1}{m}\sum_{i\in
I}\tfrac{1}{2}d^{2}_{C_{i}}$, where the penalty parameter is in theory driven
to $+\infty$ (see [18, Section 5.5] for a gentle introduction). The parameter
$\rho$ has to be carefully driven to infinity. We will apply PDA to the
problem formulation given by 5. Set $h(x,y):=\alpha\|x-y\|$, which is
$\sqrt{2}\alpha$-Lipschitz by Proposition 4.2, for $\alpha>0$. Also, write
$z=(x,y)$. The PDA with starting point $z_{0}\in X$ generates a sequence
$(z_{k})_{k\in{\mathbb{N}}}$ via

$z_{k+1}:=\operatorname{Prox}_{\rho_{k}^{-1}h}\Big{(}\sum_{i=1}^{m}\frac{1}{m}P_{C_{i}}z_{k}\Big{)},$
(19)

where

$\displaystyle\operatorname{Prox}_{\rho_{k}^{-1}h}(x,y)=(x,y)-\frac{1}{\max\big{\\{}2,\rho_{k}\|x-y\|/\alpha\big{\\}}}\big{(}x-y,y-x\big{)}$
(20)

by 36 and where $(\rho_{k})_{k\in\mathbb{N}}$ is a sequence of positive (and
“sufficiently large”) parameters. If the sets $C_{i}$ are Cartesian products
of halfspaces, we may use 2 to compute $P_{C_{i}}$. Under suitable choices of
the parameter sequences, the sequence $(z_{k})_{k\in{\mathbb{N}}}$
approximates a solution of 5. Lange and Keys recommend
$\rho_{k}=\min\\{(1.2)^{k}\rho_{0},\rho_{\text{max}}\\}$ but other choices may
yield better performance (see [19, Sections 4 and 5] and [16] for details).
Keys, Zhou, and Lange also point out a Nesterov-style accelerated version of
PDA (accPDA for short), which proceeds as follows:

$\displaystyle w_{k}$ $\displaystyle:=z_{k}+\frac{k-1}{k+2}(z_{k}-z_{k-1}),$
(21a) $\displaystyle z_{k+1}$
$\displaystyle:=\operatorname{Prox}_{\rho_{k}^{-1}h}\Big{(}\sum_{i=1}^{m}\frac{1}{m}P_{C_{i}}w_{k}\Big{)}.$
(21b)

See [16, Algorithm 1 and Section 3] for further information. Note that because
$\rho_{k}\leq\rho_{\text{max}}<+\infty$, both PDA and accPDA find a solution
of the penalized but not of the original problem.

###  3.5 Stochastic Subgradient Descent

The roots of stochastic gradient descent can be traced back to two key papers
from the early 1950s co-authored by Robbins and Monro [22] and by Kiefer and
Wolfowitz [17]; see also [10] for a recent survey. The method has since been
generalized to many different settings. We follow largely the presentation in
[8]. Set $f_{0}(x,y):=\alpha\|x-y\|$ and $(\forall i\in I)$
$f_{i}:=Ld_{C_{i}}$, where $L>0$. Then $f_{0}$ is $\sqrt{2}\alpha$-Lipschitz
and (by 35 below)

$\displaystyle
f_{0}^{\prime}(z)=f_{0}^{\prime}(x,y)=\alpha\big{(}\operatorname{sign}(x-y),-\operatorname{sign}(x-y)\big{)},$
(22)

where “$\operatorname{sign}$” is defined in 25. The other functions $f_{i}$
are $L$-Lipschitz. Moreover, for $i\in I$, we have

$f_{i}^{\prime}(z)=L\operatorname{sign}(z-P_{C_{i}}z)\in\partial f_{i}(z)$
(23)

by, e.g., [5, Example 16.62]. If the sets $C_{i}$ are Cartesian products of
halfspaces, we may use 2 to compute $P_{C_{i}}$. Now Stochastic Subgradient
Descent, which we abbreviate as SSD (see [8, Section 8.3] for further
information), applied to 6 generates a sequence via

$z_{k+1}:=z_{k}-\eta_{k}f^{\prime}_{i_{k}}(z_{k}),$ (24)

where $(\eta_{k})_{{k\in{\mathbb{N}}}}$ is a sequence of positive parameters
(typically constant or a constant divided by $\sqrt{k+1}$) and where
$f^{\prime}_{i_{k}}(z_{k})\in\partial f_{i_{k}}(z_{k})$ where $i_{k}$ is
uniformly sampled from $I_{0}:=\\{0\\}\cup I$.

Under appropriate conditions, the sequence generated by SSD approximates a
minimizer of the function $\alpha\|x-y\|+L\sum_{i\in I}d_{C_{i}}$. Note that
for large $L$, the distance functions converge pointwise to the corresponding
indicator functions, but they are different for fixed $L$. In this sense, SSD
finds a perturbed but not exact solution of the original problem.

##  4 Useful operators

In this section, we collect formulae for operators that are used later in our
numerical experiments.

###  4.1 Prox and (sub)differential operators

Denote the standard unit ball by $\mathbb{B}$: $\mathbb{B}:=\big{\\{}{y\in
Y}~{}\big{|}~{}{\|y\|\leq 1}\big{\\}}$. It will be convenient to define the
generalized signum functions on $Y$ via

$\operatorname{sign}(x):=\begin{cases}x/\|x\|,&\text{if $x\neq 0$;}\\\
0,&\text{if
$x=0$}\end{cases}\;\;\text{and}\;\;\operatorname{Sign}(x):=\begin{cases}\\{x/\|x\|\\},&\text{if
$x\neq 0$;}\\\ \mathbb{B},&\text{if $x=0$.}\end{cases}$ (25)

By [5, Example 16.32], we have

$(\forall x\in
Y)\;\;\operatorname{sign}(x)\in\operatorname{Sign}(x)=\partial\|\cdot\|(x).$
(26)

######  Proposition 4.1.

Let $\alpha>0$. Then the prox operator of the function

$h\colon Y\times
Y\to\mathbb{R}\colon(x,y)\mapsto\alpha\tfrac{1}{2}\|x-y\|^{2}.$ (27)

is given by

$\operatorname{Prox}_{h}\colon(x,y)\mapsto\frac{1}{2\alpha+1}\big{(}(1+\alpha)x+\alpha
y,\alpha x+(1+\alpha)y\big{)}.$ (28)

Moreover, $\nabla h\colon(x,y)\mapsto\big{(}\alpha(x-y),\alpha(y-x)\big{)}$ is
$(1+2\alpha)$-Lipschitz continuous.

Proof. Set $z=(x,y)\in Y\times Y$ and $B\colon Y\times Y\to
Y\colon(x,y)\mapsto\sqrt{\alpha}(x-y)$. Then

$h(z)=\tfrac{1}{2}\|Bz\|^{2}=\tfrac{1}{2}\left\langle{Bz},{Bz}\right\rangle=\tfrac{1}{2}\left\langle{B^{*}Bz},{z}\right\rangle$
(29)

and thus $\nabla h=B^{*}B$. It follows that

$\operatorname{Prox}_{h}=(\operatorname{Id}+B^{*}B)^{-1}.$ (30)

Write $B$ in block matrix form,
$B=\sqrt{\alpha}\begin{bmatrix}\operatorname{Id},-\operatorname{Id}\end{bmatrix}$.
Then
$B^{*}=\sqrt{\alpha}\begin{bmatrix}\operatorname{Id},-\operatorname{Id}\end{bmatrix}^{\intercal}$,

$\nabla h=B^{*}B=\alpha\begin{bmatrix}\operatorname{Id}&-\operatorname{Id}\\\
-\operatorname{Id}&\operatorname{Id}\end{bmatrix},$ (31)

and

$\operatorname{Id}+\nabla
h=\operatorname{Id}+B^{*}B=\begin{bmatrix}(1+\alpha)\operatorname{Id}&-\alpha\operatorname{Id}\\\
-\alpha\operatorname{Id}&(1+\alpha)\operatorname{Id}\end{bmatrix}.$ (32)

The largest eigenvalue of the very last matrix is $1+2\alpha$ which implies
that $1+2\alpha$ is the sharp Lipschitz constant of $\nabla h$. Finally,

$\operatorname{Prox}_{h}=(\operatorname{Id}+\nabla
h)^{-1}=\big{(}\operatorname{Id}+B^{*}B\big{)}^{-1}=\frac{1}{2\alpha+1}\begin{bmatrix}(1+\alpha)\operatorname{Id}&\alpha\operatorname{Id}\\\
\alpha\operatorname{Id}&(1+\alpha)\operatorname{Id}\end{bmatrix}$ (33)

and the result follows. $\hfill\quad\blacksquare$

######  Proposition 4.2.

Let $\alpha>0$. The function

$h\colon Y\times Y\to\mathbb{R}\colon(x,y)\mapsto\alpha\|x-y\|,$ (34)

is $\alpha\sqrt{2}$-Lipschitz and a convenient selection of $\partial h$ is
given by

$\displaystyle(x,y)$
$\displaystyle\mapsto\alpha\big{(}\operatorname{sign}(x-y),-\operatorname{sign}(x-y)\big{)}$
(35a)
$\displaystyle=\begin{cases}\big{(}\alpha(x-y)/\|x-y\|,\alpha(y-x)/\|y-x\|\big{)},&\text{if
$x\neq y$};\\\ \big{(}0,0\big{)},&\text{if $x=y$}\end{cases}$ (35b)
$\displaystyle\in\big{\\{}{\alpha(s,-s)}~{}\big{|}~{}{s\in\operatorname{Sign}(x-y)}\big{\\}}$
(35c) $\displaystyle=\partial h(x,y).$ (35d)

The prox operator of $h$ is given by

$\operatorname{Prox}_{h}\colon(x,y)\mapsto(x,y)-\frac{1}{\max\big{\\{}2,\|x-y\|/\alpha\big{\\}}}\big{(}x-y,y-x\big{)}.$
(36)

Proof. Set
$A=\begin{bmatrix}\operatorname{Id},-\operatorname{Id}\end{bmatrix}$,
$z=(x,y)$, and $f(z)=\|Az\|$. Then $Az=x-y$ and $h=\alpha f$. Furthermore,
$\|Az\|^{2}=\|x-y\|^{2}=\|x\|^{2}+\|y\|^{2}-2\left\langle{x},{y}\right\rangle\leq
2\|x\|^{2}+2\|y\|^{2}=2(\|x\|^{2}+\|y\|^{2})=2\|z\|^{2}$ which shows that $f$
is $\sqrt{2}$-Lipschitz and therefore $h$ is $\alpha\sqrt{2}$-Lipschitz. The
subdifferential formula follows using [5, Proposition 16.6(i) and Example
16.32]. Note that, for $\beta\geq 0$,

$A^{\intercal}=\begin{bmatrix}\operatorname{Id}\\\
-\operatorname{Id}\end{bmatrix},\;\;AA^{\intercal}=2\operatorname{Id},\;\;(AA^{\intercal}+\beta\operatorname{Id})^{-1}=\frac{1}{2+\beta}\operatorname{Id},$
(37)

and

$\displaystyle\operatorname{Id}-A^{\intercal}(AA^{\intercal}+\beta\operatorname{Id})^{-1}A$
$\displaystyle=\begin{bmatrix}\operatorname{Id}&0\\\
0&\operatorname{Id}\end{bmatrix}-\begin{bmatrix}\operatorname{Id}\\\
-\operatorname{Id}\end{bmatrix}\frac{1}{2+\beta}\begin{bmatrix}\operatorname{Id},-\operatorname{Id}\end{bmatrix}$
(38a) $\displaystyle=\begin{bmatrix}\operatorname{Id}&0\\\
0&\operatorname{Id}\end{bmatrix}-\frac{1}{2+\beta}\begin{bmatrix}\operatorname{Id}&-\operatorname{Id}\\\
-\operatorname{Id}&\operatorname{Id}\end{bmatrix}$ (38b)
$\displaystyle=\frac{1}{2+\beta}\begin{bmatrix}1+\beta&1\\\
1&1+\beta\end{bmatrix}.$ (38c)

We now discuss cases.

_Case 1:_ $\|x-y\|\leq 2\alpha$.
In view of 37 (with $\beta=0$), this is equivalent to
$\|(AA^{\intercal})^{-1}Az\|\leq\alpha$. Using [8, Lemma 6.68] and 38 (with
$\beta=0$), we obtain (switching back to row vector notation for convenience)

$\operatorname{Prox}_{h}(x,y)=\frac{1}{2}(x+y,x+y)=(x,y)-\frac{1}{2}(x-y,y-x).$
(39)

_Case 2:_ $\|x-y\|>2\alpha$.
In view of 37 (with $\beta=0$), this is equivalent to
$\|(AA^{\intercal})^{-1}Az\|>\alpha$. Using again 37 (with general $\beta\geq
0$), we set and obtain

$\displaystyle g(\beta)$
$\displaystyle=\|(AA^{\intercal}+\beta\operatorname{Id})^{-1}Az\|^{2}-\alpha^{2}$
(40a) $\displaystyle=\frac{1}{(2+\beta)^{2}}\|x-y\|^{2}-\alpha^{2}.$ (40b)

Because $\alpha>0$ and $\beta\geq 0$, we obtain the equivalences: $g(\beta)=0$
$\Leftrightarrow$ $\|x-y\|=\alpha(2+\beta)$ $\Leftrightarrow$
$\beta=\|x-y\|/\alpha-2$. Set

$\beta^{*}:=\frac{\|x-y\|}{\alpha}-2>0$ (41)

so that $g(\beta^{*})=0$. Moreover, $2+\beta^{*}=\|x-y\|/{\alpha}$ and
$1+\beta^{*}=(\|x-y\|-\alpha)/\alpha$. Hence 38 (with $\beta=\beta^{*}$) turns
into

$\displaystyle\operatorname{Id}-A^{\intercal}(AA^{\intercal}+\beta^{*}\operatorname{Id})^{-1}A$
$\displaystyle=\frac{1}{2+\beta^{*}}\begin{bmatrix}1+\beta^{*}&1\\\
1&1+\beta^{*}\end{bmatrix}$ (42a)
$\displaystyle=\frac{\alpha}{\|x-y\|}\begin{bmatrix}\frac{\|x-y\|-\alpha}{\alpha}&1\\\
1&\frac{\|x-y\|-\alpha}{\alpha}\end{bmatrix}$ (42b)
$\displaystyle=\begin{bmatrix}1&0\\\
0&1\end{bmatrix}-\frac{\alpha}{\|x-y\|}\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}.$ (42c)

Using [8, Lemma 6.68] and 42, we obtain (switching back to row vector notation
for convenience)

$\operatorname{Prox}_{h}(x,y)=(x,y)-\frac{\alpha}{\|x-y\|}(x-y,y-x).$ (43)

Finally, the formula given in 36 follows by combining 39 and 43.
$\hfill\quad\blacksquare$

###  4.2 Argmax operator

Consider, for $\alpha>0$ and $\varepsilon>0$,

$f_{0}(x,y):=\alpha\tfrac{1}{2}\|x-y\|^{2}+\varepsilon\tfrac{1}{2}\big{(}\|x\|^{2}+\|y\|^{2}),$
(44)

which is a perturbation of $\alpha\tfrac{1}{2}\|x-y\|^{2}$ that is
$\varepsilon$-strongly convex.

Given $(u,v)\in X$, the dual-based proximal methods of Section 3.3 require
from us to find the unique maximizer of

$\displaystyle(x,y)$
$\displaystyle\mapsto\left\langle{(u,v)},{(x,y)}\right\rangle-f_{0}(x,y)$
(45a)
$\displaystyle=\left\langle{u},{x}\right\rangle+\left\langle{v},{y}\right\rangle-\alpha\tfrac{1}{2}\|x-y\|^{2}-\varepsilon\tfrac{1}{2}\big{(}\|x\|^{2}+\|y\|^{2}).$
(45b)

We now derive an explicit formula for this maximizer:

######  Proposition 4.3.

Given $\alpha>0$, $\varepsilon>0$, and $(u,v)\in X$, the unique maximizer of
45 is

$(x,y)=\frac{1}{(2\alpha+\varepsilon)\varepsilon}\big{(}(\alpha+\varepsilon)u+\alpha
v,(\alpha+\varepsilon)v+\alpha u\big{)}.$ (46)

Proof. Because $f_{0}$ is strongly convex, we employ standard convex calculus
to find the maximizer by finding the zero of the gradient of the function in
45. That is, we need to solve

$\displaystyle(0,0)$ $\displaystyle=\big{(}u-\alpha x+\alpha y-\varepsilon
x,v-\alpha y+\alpha x-\varepsilon y\big{)};$ (47)

or equivalently (switching to more formal column vector notation),

$\begin{bmatrix}\alpha+\varepsilon&-\alpha\\\
-\alpha&\alpha+\varepsilon\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}=\begin{bmatrix}u\\\ v\end{bmatrix}.$ (48)

Because

$\begin{bmatrix}\alpha+\varepsilon&-\alpha\\\
-\alpha&\alpha+\varepsilon\end{bmatrix}^{-1}=\frac{1}{(2\alpha+\varepsilon)\varepsilon}\begin{bmatrix}\alpha+\varepsilon&\alpha\\\
\alpha&\alpha+\varepsilon\end{bmatrix}$ (49)

we obtain the announced formula. $\hfill\quad\blacksquare$

##  5 Numerical experiments

We start by discussing metrics — in the sense of “standard of measurement” not
in the sense of topology — to evaluate the quality of the the iterates for the
different algorithms proposed in Section 3. To make the measurement uniform
over the different algorithms, we use this metric once for each prox or proj
evaluation. For the same purpose, we consider $A_{i}$ or $B_{i}$ as unit
inputs in these operators. So, for example, as DR projects to all $A_{i}$ and
$B_{i}$ along with a prox evaluation, we get a total of $2m+1$ operations.
Since it is not possible to obtain the final output of a given iteration
before computing all the proxes, we repeat the final output given by the
iteration $2m+1$ times. Again for the sake of uniformity, we only repeat the
final output of each iteration, regardless of whether intermediate updates can
be calculated. We shall consider two cases: in the first, true solutions are
known (and assumed to be unique) while in the second, they aren’t. The former
case allows us to inspect the progress of the iterates towards the solution,
while in the latter case a metric is needed to gauge the performance of the
algorithms. The latter scenario is the one most realistic for applications.

###  5.1 Two examples

Figure 1: Visualization of the sets A (red) and B (blue) for Example 5.1 and
the solution pair (black dots) for 1.

The convergence plots for the algorithms are straightforward when the solution
is known. Assuming the solution is unique and denoted by $(\bar{x},\bar{y})$,
we use the metric

$(x,y)\mapsto\|(x,y)-(\bar{x},\bar{y})\|$ (50)

applied to the appropriate iterates of the algorithms.

The distance to the solution 50 is evaluated once for each projection or prox
operator evaluation. For example, because each iteration of DR (see Section
3.2) uses $2m+1$ prox evaluation (see the first paragraph of Section 5), one
DR step invokes $2m+1$ “updates”. This approach ensures that the evaluation
has some uniformity/fairness over the different algorithms.

We begin with an example where we know the solution. This (small-scale)
example is motivated by the one provided by Aharoni et al. in [1, Section 5].

######  Example 5.1.

Consider the subsets $A$ and $B$ of $Y=\mathbb{R}^{2}$, defined by the two
systems of $m=4$ linear inequalities

Figure 2: Convergence plot for Example 5.1 using the known-solution metric

$\left[\begin{array}[]{rrr}4&3&-17\\\ 1&0&4\\\ 1&1&11\\\
0&1&5\end{array}\right]\begin{bmatrix}x_{1}\\\ x_{2}\\\
1\end{bmatrix}\leq\left[\begin{array}[]{r}0\\\ 0\\\ 0\\\
0\end{array}\right]\quad\text{and}\quad\left[\begin{array}[]{rrr}5&-4&-30\\\
1&-2&0\\\ -1&-4&24\\\ -2&-1&13\end{array}\right]\begin{bmatrix}x_{1}\\\
x_{2}\\\ 1\end{bmatrix}\leq\left[\begin{array}[]{r}0\\\ 0\\\ 0\\\
0\end{array}\right],$ (51)

respectively. The corresponding problem 1 possesses the unique solution

$\big{(}(-6,-5),(4,5)\big{)},$ (52)

which is also visualized in Fig. 1.

We ran the algorithms from Section 3, and also accelerated versions when
available. The algorithms were run for a total 1800 Prox or Projection
operations per algorithm. The accelerated versions performed clearly better
than the original versions in this case. Therefore, for the clarity of the
exposition, we do not report the DPG and PDA results. We used the following
parameters and also report to how many iterations in each algorithm this
corresponded:

  1. ACJ:

55 iterations, with $\lambda_{k}={1}/{(k+1)}$, $n_{k}=\lfloor 1.1^{k}\rfloor$,
and

$(x^{\prime}_{k},y^{\prime}_{k})=\begin{cases}(y_{0},x_{0}),&\text{if }k=0;\\\
(y_{k-1},x_{k-1}),&\text{otherwise}\end{cases}$ (see Section 3.1).

  2. DR:

200 iterations, with $p=1$, $\alpha=5$, and $\lambda=1$ (see Section 3.2).

  3. FDPG:

200 iterations, with $\alpha=1$, $L=16$, and $\varepsilon=1/4$ (see Section
3.3).

  4. accPDA:

200 iterations, with $\alpha=1$, $\rho_{0}=1$ and $\rho_{\text{max}}=100000$
(see Section 3.4).

  5. SSD:

1800 iterations, with $\alpha=1$, $L=1$, and $\eta_{k}=1/\sqrt{k+1}$ (see
Section 3.5).

In each case, we use the starting point $x_{0}=y_{0}=(8,-13)$. The distance of
the iterates from the solution, calculated using 50, results in the
convergence plot shown in Fig. 2, where the grey-dotted lines marks intervals
of 50 iterations of DR. From the plot, we see that accPDA appears to perform
better than DR for the first 50 iterations, but then DR slowly but steadily
begins to produce the most accurate solution. We note that since FDPG is
solving for a strong convex version of the objective function $\|x-y\|$, it
converges to a solution that is not the same as our original problem.

Note that ACJ and SSD do not perform as well as the other algorithms in
Example 5.1; however, when $m$ becomes larger, they become much more
competitive as the following example illustrates:

Figure 3: Example 5.2 for which $n=m=1000$

######  Example 5.2.

Let $n\in\\{1,2,\ldots\\}$, and consider the subsets $A$ and $B$ of
$Y=\mathbb{R}^{n}$, defined by the two system of $m=n$ linear inequalities

$\begin{bmatrix}x_{1}\\\ \vdots\\\ x_{n}\end{bmatrix}\geq\begin{bmatrix}5\\\
\vdots\\\ 5\end{bmatrix}\quad\text{and}\quad\begin{bmatrix}x_{1}\\\ \vdots\\\
x_{n}\end{bmatrix}\leq\begin{bmatrix}-5\\\ \vdots\\\ -5\end{bmatrix},$ (53)

respectively. Clearly, the unique solution to the corresponding problem 1 is
$(\bar{x},\bar{y})\in Y\times Y$, where $\bar{x}=(5,\ldots,5)$ and
$\bar{y}=(-5,\ldots,-5)$. This time we run the algorithms for around $10^{8}$
prox evaluations with the starting point $x_{0}=y_{0}=(0,\dots,0)$. It turns
out that in this case, PDA outperforms accPDA, so we omit the results of the
latter. We only use the error metric once after every 50000 evaluations. We
used the following parameters and also report to how many iterations in each
algorithm this corresponded:

  1. ACJ:

169 iterations, with $\lambda_{k}={1}/{(k+1)}$, $n_{k}=\lfloor
1.1^{k}\rfloor$, and

$(x^{\prime}_{k},y^{\prime}_{k})=\begin{cases}(y_{0},x_{0}),&\text{if }k=0;\\\
(y_{k-1},x_{k-1}),&\text{otherwise}\end{cases}$ (see Section 3.1).

  2. DR:

50000 iterations, with $p=1$, $\alpha=5$, and $\lambda=1$ (see Section 3.2).

  3. FDPG:

50000 iterations, with $\alpha=1$, $L=10000$, and $\varepsilon=1/10$ (see
Section 3.3).

  4. PDA:

50000 iterations, with $\alpha=1$, $\rho_{0}=1$ and $\rho_{\text{max}}=10^{5}$
(see Section 3.4).

  5. SSD:

$\approx 10^{8}$ iterations, with $\alpha=1$, $L=10$, and
$\eta_{k}=1/\sqrt{k+1}$ (see Section 3.5).

For $n=1000=m$ and the error metric given by 50, we obtain the plot shown in
Fig. 3. Note that ACJ and SSD operate in $X=Y\times Y=\mathbb{R}^{2000}$
while, for instance, DR operates in the much bigger space
$X^{m+1}=\mathbb{R}^{2001000}$! The plot makes it clear that in this situation
ACJ and SSD fare much better than in Example 5.1.

We note that given enough iterations, DR again trumps the other algorithms,
but the initial descent is very slow. In fact,the other algorithms perform
much better in the beginning than DR. This suggest an interesting topic for
further research: one could consider a hybrid approach, where one uses an
algorithm such as SSD, ACJ, PDA or FDPG, and then switches over to DR. Note
that only ACJ and DR are known to converge to a solution of the original — the
algorithms SSD, PDA, and FDPG solve perturbed versions and thus can play a
role in quickly getting “close” to nearby points in a preprocessing capacity.

###  5.2 What to do in the absence of known solutions

Figure 4: Using 54 to measure performance of the algorithms

In general, one has no access to true solutions, so it becomes necessary to
measure performance by a metric different from 50. We propose the measure

$D_{\delta}(x,y):=\|x-y\|+\frac{1}{\delta}\sum_{i\in I}d_{C_{i}}(x,y),$ (54)

where $\delta>0$. Because the problem asks to find a point in $A\times B$,
feasibility is of greater importance than minimizing $\|x-y\|$; thus, a
smaller value of $\delta$ is desirable to stress feasibility. (One could also
envision a “dynamic” metric, where $\delta\to 0$ as the number of iterations
increases, but we have not tested this.) Revisiting the problem considered in
Example 5.1, we show in Fig. 4 the convergence plot using the parameter
$\delta=1/10$. This time, the horizontal axis was taken with a log scaling to
increase readability of the resulting graph. The behaviour of the algorithms
in the plot on the right, which includes the feasibility conditions along with
the distance between $x$ and $y$, resembles the one seen in Fig. 2.

###  5.3 Combining constraints

Figure 5: Comparing variants of DR on Example 5.1 with paired projections

In some cases, it is possible to combine constraints and still be able to
compute the projection onto the intersection. For instance, if all sets
$A_{i}$ are halfspaces, then any two sets $A_{i}$ and $A_{j}$ may be combined
and the projection onto the intersection is explicitly available using, e.g.,
[5, Proposition 2.22–2.24].

Revising Example 5.1 in this light, we ran DR with the paired projection (PP),
with choosing mostly non-adjacent halfspaces (labelled as DR+PP in Fig. 5),
and also with choosing explicitly adjacent halfspaces (labelled as DR+PP
adjacent in Fig. 5) along with the original version of DR (labelled as DR in
Fig. 5) To compare this fairly, we count one “paired” projection (onto the
intersection of two halfspaces) as being equivalent to two regular
projections. The convergence plot shown in Fig. 5 illustrates that adding the
paired projections significantly improves the performance significantly, even
more so when the halfspaces are adjacent. The well known and characteristic
“rippling” seen in typical DR curves is heavily damped in the last case.

Figure 6: Comparing ACJ on Example 5.1 with paired projections Figure 7:
Appearance of the iterations of ACJ for the paired projections

Using this technique on ACJ, we note that the paired-projection variants do
not improve the performance significantly; see Fig. 6. On the other hand, the
approach of the iterates to the true solution looks far less scattered as can
be seen in Fig. 7.

In higher dimensions, further investigations are needed to determine the
“best” way to pair up halfspaces as “adjacent halfspaces”.

##  6 Conclusion

We revisited the recent study by Aharoni, Censor, and Jiang on finding best
approximation pairs of two polyhedra. The framework we proposed works for two
sets that are themselves finite intersections of closed convex sets with
“simple” projections. Several algorithms were proposed and the required prox
operators were computed. Our numerical experiments suggested that other
algorithms deserve serious consideration.

## Acknowledgments

The authors thank the editor and the reviewers for helpful suggestions and
constructive feedback which helped us to improve the presentation of the
results, and Patrick Combettes for pointing out the relevant reference [13].
HHB and XW are supported by the Natural Sciences and Engineering Research
Council of Canada.

## References

  * [1] R. Aharoni, Y. Censor, Z. Jiang, Finding a best approximation pair of points for two polyhedra, _Computational Optimization and Applications_ 71 (2018), 509–523. https://doi.org/10.1007/s10589-018-0021-3
  * [2] H.H. Bauschke, The approximation of fixed points of compositions of nonexpansive mappings in Hilbert space, _Journal of Mathematical Analysis and Applications_ 202 (1996), 150–159. https://doi.org/10.1006/jmaa.1996.0308
  * [3] H.H. Bauschke and J.M. Borwein, On the convergence of von Neumann’s alternating projection algorithm for two sets, _Set-Valued Analysis_ 1 (1993), 185–212. https://doi.org/10.1007/BF01027691
  * [4] H.H. Bauschke and J.M. Borwein, Dykstra’s alternating projection algorithm for two sets, _Journal of Approximation Theory_ 79 (1994), 418–443. https://doi.org/10.1006/jath.1994.1136
  * [5] H.H. Bauschke and P.L. Combettes, _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_ , second edition, Springer, 2017. https://doi.org/10.1007/978-3-319-48311-5
  * [6] H.H. Bauschke, P.L. Combettes, and D.R. Luke, Finding best approximation pairs relative to two closed convex sets in Hilbert spaces, _Journal of Approximation Theory 127_ (2004), 178–192. https://doi.org/10.1016/j.jat.2004.02.006
  * [7] H.H. Bauschke and W.M. Moursi, On the behavior of the Douglas-Rachford algorithm for minimizing a convex function subject to a linear constraint, _SIAM Journal on Optimization_ 30 (2020), 2559–2576. https://doi.org/10.1137/19M1281538
  * [8] A. Beck, _First-Order Methods in Optimization_ , SIAM 2017. https://doi.org/10.1137/1.9781611974997
  * [9] A. Beck and M. Teboulle, A fast dual proximal-gradient algorithm for convex minimization and applications, _Operations Research Letters_ 42 (2014), 1–6. https://doi.org/10.1016/j.orl.2013.10.007
  * [10] L. Bottou, F.E. Curtis, and J. Nocedal, _Optimization methods for large-scale machine learning_ , _SIAM Review_ 60 (2018), 223–311. https://doi.org/10.1137/16M1080173
  * [11] Y. Censor, Computational acceleration of projection algorithms for the linear best approximation problem, _Linear Algebra and its Applications_ 416 (2006), 111–123. https://doi.org/10.1016/j.laa.2005.10.006
  * [12] W. Cheney and A.A. Goldstein, Proximity maps for convex sets, _Proceedings of the AMS_ 10 (1959), 448–450. https://doi.org/10.2307/2032864
  * [13] P.L. Combettes, Ð. Dũng, and B.C. Vũ, Dualization of signal recovery problems, _Set-Valued and Variational Analysis_ 18 (2010), 373–404. https://doi.org/10.1007/s11228-010-0147-7
  * [14] J. Douglas and H.H. Rachford, On the numerical solution of heat conduction problems in two and three space variables, _Transactions of the AMS_ 82 (1956), 421–439. https://doi.org/10.1090/S0002-9947-1956-0084194-4
  * [15] B. Halpern, Fixed points of nonexpanding maps, _Bulletin of the AMS_ 73 (1967), 957–961. https://doi.org/10.1090/S0002-9904-1967-11864-0
  * [16] K.L. Keys, H. Zhou, and K. Lange, Proximal distance algorithms: theory and practice, _Journal of Machine Learning Research_ 20 (2019), 1–38. https://jmlr.csail.mit.edu/papers/v20/17-687.html
  * [17] J. Kiefer and J. Wolfowitz, Stochastic estimation of the maximum of a regression function, _Annals of Mathematical Statistics_ 23(3) (1952), 462–466. https://doi.org/10.1214/aoms/1177729392
  * [18] K. Lange, _MM Optimization Algorithms_ , SIAM 2016. https://doi.org/10.1137/1.9781611974409
  * [19] K. Lange and K.L. Keys, The proximal distance algorithm, in _Proceedings of the 2014 International Congress of Mathematicians_ , pages 96–116, Seoul: Kyung Moon, 4. See also https://arxiv.org/abs/1507.07598
  * [20] P.-L. Lions, Approximation de points fixes de contractions, _Comptes Rendus Hebdomadaires des Séances de l’Académie des Sciences Séries A (Sciences Mathématiques)) et B (Sciences Physiques)_ 284(21), (June 6, 1977), 1357–1359. https://gallica.bnf.fr/ark:/12148/bpt6k5731057m
  * [21] P.-L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators, _SIAM Journal on Numerical Analysis_ 16 (1979), 964–979. https://doi.org/10.1137/0716071
  * [22] H. Robbins and S. Monro, A stochastic approximation method, _Annals of Mathematical Statistics_ 22(3) (1951), 400–407. https://doi.org/10.1214/aoms/1177729586
  * [23] R. Wittmann, Approximation of fixed points of nonexpansive mappings, _Archiv der Mathematik_ 58 (1992), 486–491. https://doi.org/10.1007/BF01190119