Numerical explorations of feasibility algorithms for finding points in the intersection of finite sets

TL;DR

This paper investigates the effectiveness of various projection algorithms in solving feasibility problems involving finite sets, highlighting how algorithm choice and tuning parameters influence convergence in nonconvex scenarios.

Contribution

It provides a numerical analysis of four projection methods applied to finite set feasibility problems, exploring their behavior across different configurations.

Findings

Algorithm choice significantly affects convergence behavior.

Tuning parameters play a crucial role in the success of the methods.

Different problem configurations exhibit varied local and global convergence patterns.

Abstract

Projection methods are popular algorithms for iteratively solving feasibility problems in Euclidean or even Hilbert spaces. They employ (selections of) nearest point mappings to generate sequences that are designed to approximate a point in the intersection of a collection of constraint sets. Theoretical properties of projection methods are fairly well understood when the underlying constraint sets are convex; however, convergence results for the nonconvex case are more complicated and typically only local. In this paper, we explore the perhaps simplest instance of a feasibility algorithm, namely when each constraint set consists of only finitely many points. We numerically investigate four constellations: either few or many constraint sets, with either few or many points. Each constellation is tackled by four popular projection methods each of which features a tuning parameter. We examine the behaviour for a single and for a multitude of orbits, and we also consider local and global behaviour. Our findings demonstrate the importance of the choice of the algorithm and that of the tuning parameter.

Equations

$$x\in C:=C_1\cap \cdots \cap C_m$$

$$C_i$$

$$X={\mathbb{R}}^2$$

$$0\in C_i\subset [-10,10]\times [-10,10]$$

$$C=C_1\cap \cdots \cap C_m=\{0\}.$$

$$x_{k+1}\coloneqq\big{(}(1-\lambda)\mathrm{Id}+\lambda P_{C_{m}}\big{)}\circ\cdots\circ\big{(}(1-\lambda)\mathrm{Id}+\lambda P_{C_{1}}\big{)}x_{k}.$$

$$x_{k+1}:=x_k+\lambda \cdot \frac{{\sum }_{i=1}^m\Vert x_k-P_{C_i}{x}_k{\Vert }^2}{\Vert {\sum }_{i=1}^m(x_k-P_{C_i}{x}_k){\Vert }^2}\sum _{i=1}^m(P_{C_i}{x}_k-x_k)$$

$$(\forall i\in\{1,\ldots,m\})\quad x_{k+1,i}\coloneqq x_{k,i}+\lambda\big{(}P_{C_{i}}(2\bar{x}_{k}-x_{k,i})-\bar{x}_{k}\big{)}.$$

$$x_{k+1}\coloneqq\big{(}(1-\tfrac{\lambda}{2})P_{C_{m}}+\tfrac{\lambda}{4}\big{(}\mathrm{Id}+R_{C_{1}}R_{C_{m}}\big{)}\big{)}\circ\cdots\circ\\ \big{(}(1-\tfrac{\lambda}{2})P_{C_{2}}+\tfrac{\lambda}{4}\big{(}\mathrm{Id}+R_{C_{3}}R_{C_{2}}\big{)}\big{)}\circ\big{(}(1-\tfrac{\lambda}{2})P_{C_{1}}+\tfrac{\lambda}{4}\big{(}\mathrm{Id}+R_{C_{2}}R_{C_{1}}\big{)}\big{)}x_{k}.$$

$$d:X\to {\mathbb{R}}_{+}:x\mapsto \sqrt{\frac{{\sum }_{i=1}^m\Vert x-P_{C_i}x{\Vert }^2}{{\sum }_{i=1}^m\Vert x_0-P_{C_i}{x}_0{\Vert }^2}},$$

$$d(y_k)<ϵ:=10^{-6}$$

Full text

yclic Douglas–Rachford algorithm, Douglas–Rachford algorithm, extrapolated parallel projection method, method of cyclic projections, nonconvex feasibility problem, optimization algorithm, projection

AMS 2010 Subject Classification: 49M20, 49M27, 49M37, 65K05, 65K10, 90C25, 90C26, 90C30

.1 Introduction

\runinhead

Background

Let $X$ be a Euclidean space (i.e., a finite-dimensional Hilbert space), with inner product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. The feasibility problem is a common problem in science and engineering: given finitely many closed subsets $C_{1},\ldots,C_{m}$ of $X$, it asks to

[TABLE]

We henceforth assume that the intersection $C$ is nonempty. Algorithms for solving (FP) exist when the constraint sets $C_{i}$ allow for simple projectors $P_{C_{i}}$ (i.e., nearest-point mappings). When $C_{i}$ is convex, then the projector $P_{C_{i}}$ is a nice (firmly nonexpansive and single-valued) operator defined on the entire space $X$; when $C_{i}$ is not convex, then $P_{C_{i}}$ is nonempty and set-valued. For notational simplicity, we will use $P_{C_{i}}$ to denote an arbitrary but fixed selection of the set-valued projector. (If $S$ is a subset of $X$, then $P_{S}(x)$ is a minimizer of the function $s\mapsto\|x-s\|$, where $s\in S$. For other notions not explicitly defined in this paper, we refer the reader to BGMbib-BBor .)

Assuming that the operators $P_{C_{1}},\ldots,P_{C_{m}}$ are readily available and implementable, one may try to solve (FP) iteratively by generating a sequence $(x_{k})_{k\in\mathbb{N}}$ of vectors in $X$ that employs the projection operators $P_{C_{i}}$ in some fashion to produce the next update. There are hundreds of papers dealing with algorithms for solving convex or nonconvex feasibility problems. Thus, we refrain from providing a comprehensive list of references and rather point to the following recent books and “meta” papers as starting points: BGMbib-BBor ; BGMbib-BC ; BGMbib-BK ; BGMbib-Ceg ; BGMbib-CZak ; BGMbib-CZ ; BGMbib-Com97a ; BGMbib-Com97b . (We note that the recent manuscript BGMbib-BDL deals with a feasibility problem where one set is a doubleton.) The convergence theory in the nonconvex case is much more challenging and usually of local character.

\runinhead

Goal of this paper

The goal of this paper is to showcase the surprising numerical complexity of the most simple instance of (FP); namely, when each constraint set {svgraybox}

[TABLE]

In this case, the projection operator is very easy to implement — this is achieved by simply measuring the distance of the point to each point in $C_{i}$ and returning the closest one. Furthermore, we will restrict ourselves to the simple case when the underlying space {svgraybox}

[TABLE]

is simply the Euclidean plane. Even in this setting, the difficulty and richness of the dynamic behaviour is impressively illustrated.

It is our hope that the complexity revealed will spark further analytical research in feasibility algorithms with the goal to explain the observed complexity and ultimately to aid in the design of new algorithms for solving difficult feasibility problems.

\runinhead

Organization of the paper

The remainder of the paper is organized as follows. In Section .2, we present the four constellations we will use for our numerical exploration throughout the remainder of the paper. These constellations correspond to feasibility problems that we will attempt to solve using the algorithms listed in Section .3. Section .4 provides details on the implementation and execution of the numerical experiments. The “best” tuning parameter $\lambda_{\text{\,best}}$ is determined in Section .5. We then track typical orbits of the algorithms in Section .6. Local and global behaviour is investigated in Section .7. Some interesting (and beautiful) behaviour outside of the main numerical experiments are collected in Section .8. The final Section .9 contains some concluding remarks.

.2 The four constellations

Even though we restrict ourselves already to finitely many constraint sets with finitely many points in the Euclidean plane, the infinitely many possibilities to experiment make it a daunting task to explore this space. We opted to probe this universe as follows.

The points in each constraint set $C_{i}$ are chosen randomly. We will ensure that the origin belongs to each set $C_{i}$ {svgraybox}

[TABLE]

to have a consistent feasibly problem with {svgraybox}

[TABLE]

We will focus on two alternatives for the number of constraint sets, either “few” or “many”. We will also consider constraint sets with a maximum number of points in the constraint sets, either “few” or “many”. From now on, we will use the following language: {svgraybox}

• •

The number of few sets is $3$.

• •

The number of many sets is $10$.

• •

The number of few points is $20$.

• •

The number of many points is $100$.

This will give rise to four constellations: few sets with few points, few sets with many points, many sets with few points, and many sets with many points. The four constellations used in our numerical experiments are shown in Figure 1.

.3 The four feasibility algorithms

We will numerically solve instances of (FP) using four algorithms which we briefly review in this section. While there is a myriad of competing algorithms available, our selection consists of trustworthy “work horses” that have been employed elsewhere and for which the convergence theory in the convex case is fairly well understood. Each of these algorithms has a “tuning” parameter $\lambda$ in the range $\left]0,2\right[$. The default value $\lambda_{\text{\,default}}$ is $1$. Guided by experiments, we will also (numerically) look for the “best” value $\lambda_{\text{\,best}}$. We now turn to these four algorithms. Each algorithm will have a governing sequence driving the iteration, and a (possibly different) monitored sequence which is meant to find a solution of (FP).

\runinhead

Cyclic Projections (CycP)

Given $x_{0}\in X$, the governing sequence is defined by {svgraybox}

[TABLE]

The default parameter is $\lambda_{\text{\,default}}=1$, from the range $\left]0,2\right[$. The sequence monitored is \big{(}\tfrac{1}{m}\sum_{i=1}^{m}P_{C_{i}}x_{k}\big{)}_{k\in\mathbb{N}}. Selected references: BGMbib-BBor ; BGMbib-BC ; BGMbib-BK ; BGMbib-Ceg ; BGMbib-CCCDH ; BGMbib-CZ ; BGMbib-Com97b .

\runinhead

Extrapolated Parallel Projections (ExParP)

Given $x_{0}\in X$, the governing and monitored sequence is defined by {svgraybox}

[TABLE]

if $x_{k}\notin C$; $x_{k+1}=x_{k}$ otherwise. The default parameter is $\lambda_{\text{\,default}}=1$, from the range $\left]0,2\right[$. Selected references: BGMbib-BC ; BGMbib-BCK ; BGMbib-Com97a .

\runinhead

Douglas–Rachford (DR)

Given $x_{0}\in X$, $\mathbf{x}_{0}\coloneqq(x_{0,1},\ldots,x_{0,m})=(x_{0},\ldots,x_{0})\in\mathbf{X}\coloneqq X^{m}$, $\mathbf{x}_{k}=(x_{k,1},\ldots,x_{k,m})\in\mathbf{X}$, and $\bar{x}_{k}\coloneqq\tfrac{1}{m}\sum_{i=1}^{m}x_{k,i}$, the next iterate is $\mathbf{x}_{k+1}=(x_{k+1,1},\ldots,x_{k+1,m})$, where {svgraybox}

[TABLE]

The default parameter is $\lambda_{\text{\,default}}=1$, from the range $\left]0,2\right[$. The sequence monitored is $(\bar{x}_{k})_{k\in\mathbb{N}}$. Selected references: BGMbib-BC ; BGMbib-BM ; BGMbib-EckBer ; BGMbib-ERT ; BGMbib-LM .

\runinhead

Cyclic Douglas–Rachford (CycDR)

Given $x_{0}\in X$, the governing sequence is defined by {svgraybox}

[TABLE]

The default parameter is $\lambda_{\text{\,default}}=1$, from the range $\left]0,2\right[$. The sequence monitored is \big{(}\tfrac{1}{m}\sum_{i=1}^{m}P_{C_{i}}x_{k}\big{)}_{k\in\mathbb{N}}. Selected references: BGMbib-Luke ; BGMbib-DP18a ; BGMbib-LST . (For other cyclic version of DR, see BGMbib-BNP ; BGMbib-BT . Also, if $m=2$ and $C_{1}=C_{2}=\{0\}$, then $(x_{k})_{k\in\mathbb{N}}=((\lambda/2)^{k}x_{0})_{k\in\mathbb{N}}$ is actually unbounded when $x_{0}\neq 0$ and $\lambda>2$.)

.4 Setting up the numerical explorations

\runinhead

Stopping criteria The feasibility measure

[TABLE]

where $x_{0}\in X\smallsetminus C$, vanishes exactly when $x\in C$. We stop each algorithm with monitored sequence $(y_{k})_{k\in\mathbb{N}}$ either when

[TABLE]

or when the maximum number of iterations, which we set to 1000, is reached. These values were chosen to allow a reasonable exploration of the feasibility problem while maintaining computational efficiency.

\runinhead

Details on program A program was developed in C++ to run the different experiments, see Figure 2 for two screenshots which we describe next. In the main tab of the user interface one can select the algorithm to be used and set up the problem to be solved by choosing the number of sets and the maximum number of elements per set. By clicking on the diagram showing the current constellation of points, the user can select a starting point and immediately observe the resulting orbit being rendered over the constellation. The graph of the feasibility measure $d$, corresponding to the current orbit, is also displayed.

The Cartographer tab allows the exploration of a very large number of starting points to construct a picture of the performance of a given algorithm. This two-dimensional plot shows for each starting point the number of iterations required to solve the problem, ranging from zero (black) to the maximum number of iterations allowed (white). The plot is generated progressively and uses Quasi-Monte Carlo sampling for the selection of the starting points. This is the most computationally intensive part of the software, and it is fully multi-threaded to take advantage of modern processor architectures.

.5 Determining the “best” parameter $\lambda_{\text{\,best}}$

In this section, we consider our four constellations (see Section .2) and run on each of them the four algorithms (see Section .3) with the parameter $\lambda$ ranging over $\left]0,2\right[$. The curves shown in Fig. 3, 4, 5, and 6 give an estimate of the success rate of each algorithm, evaluated for 200 evenly-spaced values of $\lambda$. For each value of $\lambda$, 5000 starting points are drawn from $[-10,10]\times[-10,10]$ using Quasi-Monte Carlo sampling, and the success rate is estimated by dividing the number of times the algorithm is successful by this number of starting points. Thus, a “best” parameter $\lambda_{\text{\,best}}$ is determined. It is this parameter that we will use to compare with the default parameter $\lambda_{\text{\,default}}$, which is $1$ in all cases.

\runinhead

Discussion

For each of the four constellations considered above, we visually inspected the $\lambda$-curves indicating success rates. We then picked for each algorithm a parameter called $\lambda_{\text{\,best}}$ which improved performance over the default parameter $\lambda_{\text{\,default}}=1$. The results are recorded in the following table.

We will use the parameters $\lambda_{\text{\,best}}$ for the experiments in subsequent sections.

.6 Tracking orbits

In this section, we consider our four given constellations (see Section .2). For each constellation, which is organized in a separate subsection, the same starting point is used. We then consider each of our four fixed algorithms (see Section .3) and show orbits for $\lambda_{\text{\,default}}=1$ and for $\lambda_{\text{\,best}}$ (see Section .5), and the corresponding feasibility measure $d$ (see Section .4).

.6.1 Few sets with few points

.6.2 Few sets with many points

.6.3 Many sets with few points

.6.4 Many sets with many points

.6.5 Discussion

The numerical results in this subsection suggest the following: The most challenging constellation is the one with few sets and many points. The least challenging constellation is the one with many sets and few points for which all algorithms are successful.

The worst algorithm is CycP. ExParP with $\lambda_{\text{\,best}}$ solves all four constellations. DR solves all constellations but $\lambda$ has to be chosen appropriately. CycDR works well with $\lambda_{\text{\,default}}$ in terms of number of iterations required; however, it was not able to solve the constellation with few sets and many points.

The experiments in this section suggest that (i) ExParP, DR, and CycDR are algorithms worthwhile exploring and that (ii) experimenting with $\lambda$ may lead to improved convergence.

Because the results in this section feature a fixed starting point, we will explore in the next section the four constellations for a multitude of starting points.

.7 Local and global behaviour

In this section, we continue to consider our four constellations (see Section .2) which our four algorithms attempt to solve (see Section .3).

In contrast to Section .6 where we tracked a single orbit, we here illustrate local and global behaviour of the algorithms for a multitude of starting points, sampled from $[-10,10]\times[-10,10]$ and $[-100,100]\times[-100,100]$, respectively. We do this for $\lambda_{\text{\,default}}=1$ and for $\lambda_{\text{\,best}}$ (see the table in Figure 7 in Section .5);

For each starting point in the given range, these plots display as gray levels the number of iterations the algorithm needed in its attempt to solve the problem represented by the given constellation. Black corresponds to the minimum number of iterations (zero), and white to the maximum number of iterations (1000). The latter is obtained when the algorithm is unsuccessful. Therefore, the darker the image, the better the performance.

To quantitatively assess the performance of each algorithm, success rates are also provided. These are obtained by dividing the number of times the algorithm is successful by the number of starting points used.

Each of these images was generated using at least 15 million starting points. Depending on the constellation, the time required to generate these pictures ranged between a few minutes to about 3 hours using a quad-core computer.

.7.1 Few sets with few points

.7.2 Few sets with many points

.7.3 Many sets with few points

.7.4 Many sets with many points

.7.5 Discussion

Comparing the success rates reported in the figures above, it appears that ExParP, DR, and CycDR are good choices; we recommend that CycP be not used. The effect of the tuning parameter $\lambda$ is very striking for most algorithms when comparing performance of $\lambda_{\text{\,default}}$ with $\lambda_{\text{\,best}}$.

.8 Divertissements

We experimented also with other constellations and encountered some interesting behaviour of ExParP. This algorithm seems to exhibit fractal-like behaviour for some constellations — whether they are created randomly or not. In the following, we present three images that we found particularly delightful in Figure 16 and Figure 17.

.9 Concluding remarks

We encountered a somewhat surprising complexity in the behaviour of four algorithms for solving feasibility problems in a simple nonconvex case. The importance of the tuning parameter $\lambda$ is apparent as is the proximity to solutions (local vs global behaviour). Further studies are needed to find effective guidelines for users in terms of choice of algorithms and the choice of the parameter $\lambda$. Finally, and similarly to BGMbib-Boretal , we encountered beauty in our numerical explorations. It is our hope that others will join us and explore theoretically and numerically this fascinating universe of constellations.

Acknowledgements

We thank the referee for constructive comments and suggestions. The research of HHB was partially supported by NSERC.