# Spectral Projected Subgradient Method for Nonsmooth Convex Optimization
Problems

Nataša Krejić111Department of Mathematics and Informatics, Faculty of
Sciences, University of Novi Sad, Trg Dositeja Obradovića 4, 21000 Novi Sad,
Serbia. e-mail<EMAIL_ADDRESS>Nataša Krklec Jerinkić,222Department of
Mathematics and Informatics, Faculty of Sciences, University of Novi Sad, Trg
Dositeja Obradovića 4, 21000 Novi Sad, Serbia. e-mail:
<EMAIL_ADDRESS>Tijana Ostojić 333Department of Fundamental
Sciences, Faculty of Technical Sciences, University of Novi Sad, Trg Dositeja
Obradovića 6, 21000 Novi Sad, Serbia. e-mail<EMAIL_ADDRESS>
(August 8, 2022)

###### Abstract

We consider constrained optimization problems with a nonsmooth objective
function in the form of mathematical expectation. The Sample Average
Approximation (SAA) is used to estimate the objective function and variable
sample size strategy is employed. The proposed algorithm combines an SAA
subgradient with the spectral coefficient in order to provide a suitable
direction which improves the performance of the first order method as shown by
numerical results. The step sizes are chosen from the predefined interval and
the almost sure convergence of the method is proved under the standard
assumptions in stochastic environment. To enhance the performance of the
proposed algorithm, we further specify the choice of the step size by
introducing an Armijo-like procedure adapted to this framework. Considering
the computational cost on machine learning problems, we conclude that the line
search improves the performance significantly. Numerical experiments conducted
on finite sum problems also reveal that the variable sample strategy
outperforms the full sample approach.

Key words: Nonsmooth Optimization, Subgradient, Spectral Projected Gradient
Methods, Sample Average Approximation, Variable Sample Size Methods, Line
Search.

##  1 Introduction

We consider the following constrained optimization problem

$\min_{x\in\Omega}f(x)=E(F(x,\xi)),$ (1)

where $\Omega\subset\mathbb{R}^{n}$ is a convex, closed set,
$F:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ is continuous and
convex function with respect to $x$, bounded from below,
$\xi:\mathcal{A}\rightarrow\mathbb{R}^{m}$ is random vector and
$(\mathcal{A},\mathcal{F},P)$ is a probability space. Notice that F is locally
Lipschitz as a consequence of convexity, [2], but possibly nonsmooth. The
importance of the stochastic optimization problem arising from various
scientific fields generated a large amount of literature in the recent years.
Due to difficulty in computing the mathematical expectation in general, the
common approach is to approximate the original objective function $f(x)$ by
applying the Sample Average Approximation (SAA) function

$f_{\mathcal{N}}(x)=\frac{1}{N}\sum_{i\in\mathcal{N}}f_{i}(x),$ (2)

where $f_{i}(x)=F(x,\xi_{i})$ and $N=|\mathcal{N}|$ determines the size of a
sample used for approximation. The sample vectors $\xi_{i},i\in\mathcal{N}$
are assumed to be independent and identically distributed (i.i.d.).

The sample size $N$ determines the precision of the approximation, but it also
influences the computational cost. In order to achieve the a.s. convergence,
one needs to push the sample size to infinity in general. Even if the original
problem is already in the SAA form, i.e., if we are dealing with finite sum
problems, the costs of employing the full sample at each iteration can be
large and thus the variable sample size (VSS) strategy is often applied.
Finding a good way of varying the sample size as well as choosing a sample is
a problem itself and it has been the subject of many research efforts (e.g.
[3], [4], [13], [16]). Although this may affect the algorithm a lot, a
suitable sample size strategy will not be the main concern of this paper. In
order to prove a.s. convergence we assume here that the sample size tends to
infinity and leave the problem of determining the optimal VSS strategy in this
framework for future work.

In this paper we propose a framework for solving nonsmooth constrained
optimization problem (1) assuming that the feasible set $\Omega$ is easy to
project on (for example a box or ball in $\mathbb{R}^{n}$). This allows us to
apply a method of the Spectral Projected Gradient type. The Spectral Projected
Gradient (SPG) method, originally proposed in [6], is well known for its
efficiency and simplicity and it has been widely used and developed as a
solver of constrained optimization problems [5], [11], [16], [27]. The step
length selection strategy in SPG method is crucial for faster convergence with
respect to classical gradient projection methods because it involves second-
order information related to the spectrum of the Hessian matrix. SPG methods
for finite sums problem have been investigated in [5], [27]. In [27] they are
used in combination with the stochastic gradient method and the convergence is
proved assuming that the full gradient is calculated in every $m$ iterations.
In [5] the subsampled spectral gradient methods are analyzed and the effect of
the choice of the spectral coefficient is investigated. In [11] the SPG
direction is employed within Inexact Restoration framework to address
nonlinear optimization problems with nonconvex constraints. The SPG methods
for problems with continuously differentiable objective function given in form
of mathematical expectation have been analyzed in [16].

Another important class of nonsmooth optimization methods are bundle methods,
especially if accuracy in the solution and reliability are a concern (see
[19], [23] and the references therein). The main idea of bundle methods is to
make an approximation of the whole subdifferential of the objective function
instead of using only one arbitrary subgradient at each point. The advantage
of these methods compared to the classical subgradient methods is that they
use more information about the local behavior of the function by approximating
the subdifferential of the objective function with subgradients from previous
iterations that form a bundle. A great number of bundle methods in combination
with Newton-type methods [9] and the trust-region method [1] have been
developed. The modification of bundle methods has been used for nonconvex
problems [22], constrained problems [25] and multi-criteria problems [21], as
well. The drawback of these methods is requirement of solving at least one
quadratic programming subproblem in each iteration, which can be time-
consuming, especially for large scale problems. Proximal bundle methods (see
for example [8], [12] [20]) are based on the bundle methodology and have the
ability to provide exact solutions even if most of the time the available
information is inaccurate, unlike their forerunner variants. The proximal
bundle method takes ideas from subgradient method and proximal method [15],
[23]. First one can be seen as extension of gradient methods in smooth
optimization and the second one is a variant of proximal point method, which
minimizes the original function plus a quadratic part.

The method proposed in this paper is a subgradient method but differs from the
existing ones in the literature in several ways. We propose a way to plug VSS-
SPG ideas into the nonsmooth framework. Since the objective function may be
nonsmooth, we have to use subgradients instead of gradients and thus we refer
to the core algorithm as SPS - Spectral Projected Subgradient method. The
spectral coefficient is calculated by employing consecutive subgradients of
possibly different SAA functions and the safeguard which provides positive,
bounded spectral coefficients is used. We prove a.s. convergence under the
standard assumptions for stochastic environment. Moreover, in order to improve
the performance of the algorithm, we also propose a line search variant of SPS
named LS-SPS. By specifying the line search technique we ensure that LS-SPS
falls into the SPS framework and thus the same convergence results hold.
Although the descent property of the search direction is desirable, it is not
necessary in each iteration to ensure the convergence result. The proposed
line search is well defined and the a.s. convergence is achieved even if the
search direction is not a descent one for the SAA function.

Although the proposed algorithms are constructed to cope with unbounded sample
sizes, they can also be applied on finite sum problems and we devote a part of
considerations to this important special class as well. Preliminary numerical
results on machine learning problems modeled by Hinge Loss functions reveal
several advantages of the proposed method: a) spectral method outperforms the
plain subgradient method; b) VSS method reduces the costs of the full sample
method when applicable; c) line search improves the behavior of the SPS method
with predefined step size sequence.

Therefore, the main contributions of this work are the following:

  * i)

The stochastic SPG method is adapted to the nonsmooth framework;

  * ii)

The a.s. convergence of the proposed SPS method is proved under the standard
assumptions;

  * iii)

The SPS is further upgraded by introducing a specific line search technique
resulting in LS-SPS;

  * iv)

Numerical results on machine learning problems show the efficiency of the
proposed method, especially LS-SPS.

This paper is structured as follows. Details of the proposed algorithm are
presented in the next section. Section 3 deals with the convergence of the
resulting algorithm and also introduces a line search algorithm. In Section 4
we provide some implementation details and discuss results of the relevant
numerical experiments. We conclude in Section 5.

##  2 The Method

In this section, we describe the subsampled spectral projected subgradient
framework algorithm for nonsmooth probems - SPS. For any given
$z\in\mathbb{R}^{n}$, we denote by $P_{\Omega}(z)$ the orthogonal projection
of $z$ onto $\Omega.$ Throughout the paper $\mathcal{N}_{k}$ denotes the
sample used to approximate the objective function and $N_{k}$ denotes its
cardinality.

Algorithm 1: SPS
(Spectral Projected Subgradient Method for Nonsmooth Optimization)

  * S0

Initialization. Given $N_{0}\in\mathbb{N},$ $x_{0}\in\Omega,$
$0<C_{1}<1<C_{2}<\infty,$ $0<\underline{\zeta}\leq\overline{\zeta}<\infty,$
$\zeta_{0}\in\left[\underline{\zeta},\overline{\zeta}\right].$ Set k = 0.

  * S1

Direction. Choose $\bar{g}_{k}\in\partial f_{\mathcal{N}_{k}}(x_{k})$ and set
$p_{k}=-\zeta_{k}\bar{g}_{k}$.

  * S2

Step size. Choose $\alpha_{k}\in\left[C_{1}/k,C_{2}/k\right]$.

  * S3

Main update. Set $x_{k+1}=P_{\Omega}(x_{k}+\alpha_{k}p_{k})$ and
$s_{k}=x_{k+1}-x_{k}$.

  * S4

Sample size update. Chose $N_{k+1}\in\mathbb{N}$.

  * S5

Spectral coefficient update. Calculate
$y_{k}=g_{\mathcal{N}_{k}}(x_{k+1})-\bar{g}_{k}$ where
$g_{\mathcal{N}_{k}}(x_{k+1})\in\partial f_{\mathcal{N}_{k}}(x_{k+1})$. Set
$\zeta_{k+1}=\min\\{\overline{\zeta},\max\\{\underline{\zeta},\frac{s_{k}^{T}s_{k}}{s_{k}^{T}y_{k}}\\}\\}.$

  * S6

Set $k:=k+1$ and go to S1.

Let us now comment the algorithm above. In Step S1 we calculate the direction
by choosing a subgradient of the SAA function $f_{\mathcal{N}_{k}}$ and taking
the opposite direction multiplied by the spectral coefficient. Notice that the
safeguard in Step S5 ensures that the negative subgradient direction is
retained. However, this direction does not have to be descent for the function
$f_{\mathcal{N}_{k}}$ since we take an arbitrary subgradient.

Within Step S2 we choose the step size $\alpha_{k}$ from the given interval.
The constants $C_{1}$ and $C_{2}$ can be arbitrary small and large,
respectively, allowing a wide range of feasible step sizes. This choice was
motivated by the common assumption on the step size sequence in stochastic
algorithms:

$\sum_{k=1}^{\infty}\alpha_{k}^{2}<\infty,\quad\sum_{k=1}^{\infty}\alpha_{k}=\infty.$
(3)

Notice that the choice in S2 ensures that the sequence of step sizes of SPS
algorithm satisfies (3). After finding the direction and the step size, we
project the point $x_{k}+\alpha_{k}p_{k}$ onto the feasible set $\Omega$ and
thus we retain feasibility in all the iterations of the algorithm.

In Step S4 we chose the sample size to be used in the subsequent iteration. In
order to prove convergence result, we will assume that $N_{k}$ tends to
infinity or achieves and retains at the maximal sample size in the case of
finite sums. So, the simplest way to ensure this is to increase the sample
size at each iteration. However, we state this step in the most generic way to
emphasize that other choices are feasible as well, including some adaptive
strategies.

In Step S5 $y_{k}$ is calculated as a difference of two subgradients of the
same approximate function $f_{\mathcal{N}_{k}}$, but different approaches are
feasible as well. For instance, one can use subgradients of different
functions $y_{k}=g_{\mathcal{N}_{k+1}}(x_{k+1})-g_{\mathcal{N}_{k}}(x_{k})$.
This can reduce the costs, especially if the sample is not cumulative, but
also brings additional noise into the spectral coefficient since the
subgradients are calculated for two different functions in general. However,
if we have a finite sum problem and the full sample is reached, the cost of
calculating the subgradient may be reduced to one subgradient per iteration
since one can obviously take $\bar{g}_{k+1}=g_{\mathcal{N}_{k}}(x_{k+1})$. In
general, another choice could be
$y_{k}=g_{\mathcal{N}_{k+1}}(x_{k+1})-g_{\mathcal{N}_{k+1}}(x_{k}).$ This
reduces the influence of a noise and usually provides better approximation of
the spectral coefficient of the true objective function, but it requires
additional evaluations. Although the choice of $y_{k}$ was addressed in the
literature (see [5] for example), in general it remains an open question which
requires thorough analysis before drawing the final conclusions. It is
important to point out that the choice of $y_{k}$ does not affect the
convergence analysis and the theoretical results obtained in next section, but
it may affect the algorithm’s performance significantly.

##  3 Convergence

In this section we analyze conditions needed for a.s. convergence of the SPS
algorithm. A standard assumptions for stochastic environment is stated below.
Recall that samples are assumed to be i.i.d.

######  Assumption A 1.

Assume that $f_{i}(x)=F(x,\xi_{i}),$ $i=1,2,\ldots,$ are continuous, convex
and bounded from below with a constant $C.$ Moreover, assume that the function
$F$ is dominated by a P-integrable function on any compact subset of
$\mathbb{R}^{n}$.

Notice that Assumption A1 implies that $f$ is convex and continuous function
as well as $f_{\mathcal{N}}$ for any given $N$. Moreover, the Uniform Law of
Large Numbers (ULLN) implies that $f_{\mathcal{N}}(x)$ a.s converges uniformly
to $f(x)$ on any compact subset $S\subseteq\mathbb{R}^{n}$ (see Theorem 7.48
in [26] for instance), i.e.,

$\lim_{N\rightarrow\infty}\sup_{x\in S}|f_{N}(x)-f(x)|=0\quad\mbox{a.s.}$ (4)

Also notice that (4) holds trivially if the sample is finite and the full
sample is eventually achieved and retained.

The main result, a.s. convergence of Algorithm SPS, is stated in the following
theorem. We assume that the feasible set is compact, although this assumption
may be relaxed as we will show in the sequel. Moreover, recall that the
convexity implies that the functions $f_{i}$ are locally Lipschitz continuous
and thus for every $x$ and $i$ there exists $L_{i}(x)$ such that for all
$g\in\partial f_{i}(x)$ there holds $\|g\|\leq L_{i}(x)$. However, since there
can be infinitely many functions $f_{i}$ in general we assume that the chosen
subgradients are uniformly bounded. This may be accomplished by scaling the
subgradient by its norm for instance. Let $X^{*}$ and $f^{*}$ be the set of
solutions and the optimal value of problem (1), respectively. The convergence
result is as follows.

######  Theorem 3.1.

Suppose that Assumption A1 holds and $\\{x_{k}\\}$ is a sequence generated by
Algorithm SPS where $N_{k}\to\infty$. Assume also that $\Omega$ is compact and
convex and there exists $G$ such that $||\bar{g}_{k}||\leq G$ for all $k$.
Then

$\liminf_{k\to\infty}f(x_{k})=f^{*}\quad\text{a.s.}.$ (5)

Moreover,

$\lim_{k\rightarrow\infty}x_{k}=x^{*}\quad\text{a.s.}$ (6)

for some $x^{*}\in X^{*},$ provided that $\sum_{k=0}^{\infty}e_{k}/k<\infty$,
where $e_{k}:=\max_{x\in\Omega}|f_{\mathcal{N}_{k}}(x)-f(x)|.$

###### Proof.

Denote by $\mathcal{W}$ the set of all possible sample paths of SPS algorithm.
Suppose that (5) does not hold, i.e., $\liminf_{k\to\infty}f(x_{k})=f^{*}$
does not happen with probability 1. In that case there exists a subset of
sample paths $\tilde{\mathcal{W}}\subseteq\mathcal{W}$ such that
$P(\tilde{\mathcal{W}})>0$ and for every $w\in\tilde{\mathcal{W}}$ there holds

$\liminf_{k\rightarrow\infty}f(x_{k}(w))>f^{*},$

i.e., there exists $\varepsilon(w)>0$ small enough such that
$f(x_{k}(w))-f^{*}\geq 2\varepsilon(w)$ for all $k$. Since $f$ is continuous
on the feasible set $\Omega$, there exists $\tilde{y}(w)\in\Omega$ such that
$f(\tilde{y}(w))=f^{*}+\varepsilon(w).$ This further implies

$f(x_{k}(w))-f(\tilde{y}(w))=f(x_{k}(w))-f^{*}-\varepsilon(w)\geq
2\varepsilon(w)-\varepsilon(w)=\varepsilon(w).$

Let us take an arbitrary $w\in\tilde{\mathcal{W}}.$ Denote
$z_{k+1}(w):=x_{k}(w)+\alpha_{k}(w)p_{k}(w)$. Notice that nonexpansivity of
orthogonal projection and the fact that $\tilde{y}\in\Omega$ together imply

$||x_{k+1}(w)-\tilde{y}(w)||=||P_{\Omega}(z_{k+1}(w))-P_{\Omega}(\tilde{y}(w))||\leq||z_{k+1}(w)-\tilde{y}(w)||.$
(7)

Furthermore, using the fact that $\bar{g}_{k}$ is subgradient of the convex
function $f_{\mathcal{N}_{k}}$, $\bar{g}_{k}\in\partial
f_{\mathcal{N}_{k}}(x_{k})$, we have
$f_{\mathcal{N}_{k}}(x_{k})-f_{\mathcal{N}_{k}}(\tilde{y})\leq\overline{g}^{T}_{k}(x_{k}-\tilde{y})$
and dropping the $w$ in order to facilitate the reading we obtain

$\displaystyle||z_{k+1}-\tilde{y}||^{2}$ $\displaystyle=$
$\displaystyle||x_{k}+\alpha_{k}p_{k}-\tilde{y}||^{2}=||x_{k}-\alpha_{k}\zeta_{k}\overline{g}_{k}-\tilde{y}||^{2}$
(8) $\displaystyle=$
$\displaystyle||x_{k}-\tilde{y}||^{2}-2\alpha_{k}\zeta_{k}\overline{g}^{T}_{k}\left(x_{k}-\tilde{y}\right)+\alpha_{k}^{2}\zeta^{2}_{k}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-\tilde{y}||^{2}+2\alpha_{k}\zeta_{k}(f_{\mathcal{N}_{k}}(\tilde{y})-f_{\mathcal{N}_{k}}(x_{k}))+\alpha_{k}^{2}\zeta_{k}^{2}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-\tilde{y}||^{2}+2\alpha_{k}\zeta_{k}(f(\tilde{y})-f(x_{k})+2e_{k})+\alpha_{k}^{2}\zeta_{k}^{2}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-\tilde{y}||^{2}-2\alpha_{k}\zeta_{k}(f(x_{k})-f(\tilde{y}))+4e_{k}\alpha_{k}\overline{\zeta}+\alpha_{k}^{2}\overline{\zeta}^{2}G^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-\tilde{y}||^{2}-2\alpha_{k}\underline{\zeta}\varepsilon+4e_{k}\alpha_{k}\overline{\zeta}+\alpha_{k}^{2}\overline{\zeta}^{2}G^{2}$
$\displaystyle=$
$\displaystyle||x_{k}-\tilde{y}||^{2}-\alpha_{k}\left(2\underline{\zeta}\varepsilon-4e_{k}\overline{\zeta}-\alpha_{k}\overline{\zeta}^{2}G^{2}\right).$

By ULLN we have $\lim_{k\rightarrow\infty}e_{k}=0\text{ a.s.}$, or more
precisely, $\lim_{k\rightarrow\infty}e_{k}(w)=0$ for almost every
$w\in\mathcal{W}$. Since $P(\tilde{\mathcal{W}})>0$, there must exist a sample
path $\tilde{w}\in\tilde{\mathcal{W}}$ such that

$\lim_{k\rightarrow\infty}e_{k}(\tilde{w})=0.$

This further implies the existence of $\tilde{k}(\tilde{w})\in\mathbb{N}$ such
that for all $k\geq\tilde{k}(\tilde{w})$ we have

$\alpha_{k}(\tilde{w})\overline{\zeta}^{2}G^{2}+4e_{k}(\tilde{w})\overline{\zeta}\leq\varepsilon(\tilde{w})\underline{\zeta}$
(9)

because step S2 of SPS algorithm implies that $\lim_{k\to\infty}\alpha_{k}=0$
for any sample path. Furthermore, since (8) holds for all
$w\in\tilde{\mathcal{W}}$ and thus for $\tilde{w}$ as well, from (7)-(9) we
obtain

$\displaystyle||x_{k+1}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}$ $\displaystyle-$
$\displaystyle\alpha_{k}(\tilde{w})\left(2\underline{\zeta}\varepsilon(\tilde{w})-4e_{k}(\tilde{w})\overline{\zeta}-\alpha_{k}(\tilde{w})\overline{\zeta}^{2}G^{2}\right)$
$\displaystyle=$ $\displaystyle||x_{k}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}$
$\displaystyle+$
$\displaystyle\alpha_{k}(\tilde{w})\left(4e_{k}(\tilde{w})\overline{\zeta}+\alpha_{k}(\tilde{w})\overline{\zeta}^{2}G^{2}-2\underline{\zeta}\varepsilon(\tilde{w})\right)$
$\displaystyle\leq$
$\displaystyle||x_{k}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}+\alpha_{k}(\tilde{w})\left(\underline{\zeta}\varepsilon(\tilde{w})-2\underline{\zeta}\varepsilon(\tilde{w})\right)$
$\displaystyle=$
$\displaystyle||x_{k}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}-\alpha_{k}(\tilde{w})\underline{\zeta}\varepsilon(\tilde{w}).$

and

$||x_{k+s}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}\leq||x_{k}(\tilde{w})-\tilde{y}(\tilde{w})||^{2}-\varepsilon(\tilde{w})\underline{\zeta}\sum_{j=0}^{s-1}\alpha_{j}(\tilde{w}).$

Letting $s\to\infty$ yields a contradiction since
$\sum_{k=0}^{\infty}\alpha_{k}\geq\sum_{k=0}^{\infty}C_{1}/k=\infty$ for any
sample path and we conclude that (5) holds.

Now, let us prove (6) under the additional assumption
$\sum_{k=0}^{\infty}e_{k}/k<\infty$. Notice that this assumption implies that
$\sum_{k=0}^{\infty}\alpha_{k}e_{k}<\infty$ since $\alpha_{k}\leq C_{2}/k$.
Since (5) holds, we know that

$\liminf_{k\rightarrow\infty}f(x_{k}(w))=f^{*},$ (10)

for almost every $w\in\mathcal{W}$. In other words, there exists
$\overline{\mathcal{W}}\subseteq\mathcal{W}$ such that
$P(\overline{\mathcal{W}})=1$ and (10) holds for all
$w\in\overline{\mathcal{W}}$. Let us consider arbitrary
$w\in\overline{\mathcal{W}}$. We will show that
$\lim_{k\rightarrow\infty}x_{k}(w)=x^{*}(w)\in X^{*}$ which will imply the
result (6). Once again let us drop $w$ to facilitate the notation.

Let $K_{1}\subseteq\mathbb{N}$ be a subsequence of iterations such that

$\lim_{k\in K_{1}}f(x_{k})=f^{*}.$

Since $\\{{x_{k}\\}}_{k\in K_{1}}\subseteq\\{{x_{k}\\}}_{k\in\mathbb{N}}$ and
$\\{{x_{k}\\}}_{k\in\mathbb{N}}$ is bounded because of feasibility and the
compactness of the feasible set $\Omega$, there follows that
$\\{{x_{k}\\}}_{k\in K_{1}}$ is also bounded and there exist $K_{2}\subseteq
K_{1}$ and $\tilde{x}$ such that

$\lim_{k\in K_{2}}x_{k}=\tilde{x}.$

Then, we have

$f^{*}=\lim_{k\in K_{1}}f(x_{k})=\lim_{k\in K_{2}}f(x_{k})=f(\lim_{k\in
K_{2}}x_{k})=f(\tilde{x}).$

Therefore, $f(\tilde{x})=f^{*}$ and we have $\tilde{x}\in X^{*}$. Now, we show
that the whole sequence of iterates converges. Let $\\{{x_{k}\\}}_{k\in
K_{2}}:=\\{{x_{k_{i}}\\}}_{i\in\mathbb{N}}.$ Following the steps of (8) and
using the fact that $f(x_{k})\geq f(\tilde{x})$ for all $k$, we obtain that
the following holds for any $s\in\mathbb{N}$

$\displaystyle||x_{k_{i}+s}-\tilde{x}||^{2}$ $\displaystyle\leq$
$\displaystyle||x_{k_{i}}-\tilde{x}||^{2}+4\overline{\zeta}\sum_{j=0}^{s-1}e_{k_{i}+j}\alpha_{k_{i}+j}+\overline{\zeta}^{2}G^{2}\sum_{j=0}^{s-1}\alpha_{k_{i}+j}^{2}$
(11) $\displaystyle\leq$
$\displaystyle||x_{k_{i}}-\tilde{x}||^{2}+4\overline{\zeta}\sum_{j=0}^{\infty}e_{k_{i}+j}\alpha_{k_{i}+j}+\overline{\zeta}^{2}G^{2}\sum_{j=0}^{\infty}\alpha_{k_{i}+j}^{2}$
(12) $\displaystyle=$
$\displaystyle||x_{k_{i}}-\tilde{x}||^{2}+4\overline{\zeta}\sum_{j=k_{i}}^{\infty}e_{j}\alpha_{j}+\overline{\zeta}^{2}G^{2}\sum_{j=k_{i}}^{\infty}\alpha_{j}^{2}.$
(13)

Due to the fact that $\sum_{j=k_{i}}^{\infty}e_{j}\alpha_{j}$ and
$\sum_{j=k_{i}}^{\infty}\alpha_{j}^{2}$ are the residuals of convergent sums,
for any $s\in\mathbb{N}$ there holds

$||x_{k_{i}+s}-\tilde{x}||^{2}\leq a_{i},\text{ where
}\lim_{i\to\infty}a_{i}=0.$

Thus, for all $s\in\mathbb{N}$ we have

$\limsup_{i\to\infty}||x_{k_{i}+s}-\tilde{x}||^{2}\leq\limsup_{i\to\infty}a_{i}=\lim_{i\to\infty}a_{i}=0.$

Since $s\in\mathbb{N}$ is arbitrary, there follows
$\limsup_{k\to\infty}||x_{k}-\tilde{x}||^{2}=0,$ i.e.
$\lim_{k\to\infty}x_{k}=\tilde{x}$ which completes the proof. ∎

Let us comment on $e_{k}$ first. We obtain (5) provided that the sample size
tends to infinity in an arbitrary manner. The stronger result is achieved
under assumption of fast enough increase of the sample size. Having in mind
the interval for step size $\alpha_{k}$, we conclude that the assumption of
summability needed for (6) is satisfied if $e_{k}\leq C_{3}k^{-\nu}$ holds for
all $k$ large enough and arbitrary $C_{3}>0$, where $\nu>0$ can be arbitrary
small. While in general this can be hard to guarantee, for some classes of
functions $F$ (e.g. function plus noise with finite variance), the error bound
for cummulative samples derived in [13] yields $e_{k}\leq C\sqrt{\frac{\ln\ln
N_{k}}{N_{k}}},$ for all $k$ large enough and some positive constant $C$
directly dependent on the noise variance. In that case, it can be shown that
the simple choice of $N_{k}=k$ provides the sufficient growth needed for (6).

An important class of the problems that we consider is a finite sum problem

$\min_{x\in\Omega}f_{\mathcal{N}}(x)=\frac{1}{N}\sum_{i=1}^{N}f_{i}(x),$ (14)

where the functions $f_{i}$ are continuous, convex and possibly nonsmooth. In
that case, we do not need a dominance assumption since (4) is trivially
satisfied if the full sample size is eventually achieved and retained. Thus,
$e_{k}=0$ for all $k$ large enough and $\sum_{k=0}^{\infty}e_{k}/k<\infty$
trivially holds. Moreover, the compactness of the feasible set and convexity
of $f$ imply the Lipschitz continuouity of each $f_{i}$ on $\Omega$ and thus
the subgradients $\bar{g}_{k}$ are uniformly bounded. Furthermore, we also
know that the functions are uniformly bounded from below on $\Omega$. Although
the sample paths may differ, the convergence result is deterministic since the
original objective function $f_{N}$ is eventually used. We summarize the
result in the following corollary of the previous theorem.

######  Corollary 3.2.

Suppose that the functions $f_{i},i=1,...,N$ are continuous and convex and
$\\{x_{k}\\}$ is a sequence generated by Algorithm SPS applied to (14).
Suppose that $N_{k}=N$ for all $k$ large enough and $\Omega$ is compact and
convex. Then $\lim_{k\rightarrow\infty}x_{k}=x^{*}\in X^{*}$.

Since the compactness of $\Omega$ excludes the important class of constrained
problems such as $x\geq 0$ which appear as subproblems in many cases, for
instance in penalty methods, it is important to comment on the alternatives.
The assumption of bounded $\Omega$ may be replaced with the assumption of
bounded iterate sequence $\\{x_{k}\\}_{k}$ which is common in stochastic
analysis. We state the result for completeness.

######  Theorem 3.3.

Suppose that Assumption A1 holds and $\\{x_{k}\\}\subseteq\bar{\Omega}$ is a
sequence generated by Algorithm SPS where $N_{k}\rightarrow\infty$ and
$\bar{\Omega}\subseteq\Omega$ is bounded. Assume that $\Omega$ is closed and
convex and there exists $G$ such that $||\bar{g}_{k}||\leq G$ for all $k$.
Then $\liminf_{k\to\infty}f(x_{k})=f^{*}\text{a.s.}.$ Moreover,
$\lim_{k\rightarrow\infty}x_{k}=x^{*}\in X^{*}$ a.s. provided that
$\sum_{k=0}^{\infty}e_{k}/k<\infty$, where
$e_{k}:=\max_{x\in\bar{\Omega}}|f_{\mathcal{N}_{k}}(x)-f(x)|.$

Now, let us see under which conditions we obtain the boundedness of iterates.
Define an SAA error sequence as

$\bar{e}_{k}=|f_{\mathcal{N}_{k}}(x_{k})-f(x_{k})|+|f_{\mathcal{N}_{k}}(x^{*})-f(x^{*})|,$
(15)

where $x^{*}\in X^{*}$ is an arbitrary solution point. We have the following
result.

######  Proposition 3.4.

Suppose that Assumption A1 holds and $\\{x_{k}\\}$ is a sequence generated by
Algorithm SPS where $N_{k}\to\infty$. Assume that $\Omega$ is closed and
convex and there exists $G$ such that $||\bar{g}_{k}||\leq G$ for all $k$.
Then there exists a compact set $\bar{\Omega}\subseteq\Omega$ such that
$\\{x_{k}\\}\subseteq\bar{\Omega}$ provided that
$\sum_{k=0}^{\infty}\bar{e}_{k}/k\leq C_{4}<\infty.$

###### Proof.

Let $x^{*}$ be an arbitrary solution of the problem (1). Then, by following
similar steps as in (8) and using the nonexpansivity of the projection (7) we
obtain that the following holds for an arbitrary $k$

$\displaystyle||x_{k+1}-x^{*}||^{2}$ $\displaystyle\leq$
$\displaystyle||z_{k+1}-x^{*}||^{2}=||x_{k}-\alpha_{k}\zeta_{k}\overline{g}_{k}-x^{*}||^{2}$
$\displaystyle=$
$\displaystyle||x_{k}-x^{*}||^{2}-2\alpha_{k}\zeta_{k}\overline{g}^{T}_{k}\left(x_{k}-x^{*}\right)+\alpha_{k}^{2}\zeta^{2}_{k}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-x^{*}||^{2}+2\alpha_{k}\zeta_{k}(f_{\mathcal{N}_{k}}(x^{*})-f_{\mathcal{N}_{k}}(x_{k}))+\alpha_{k}^{2}\zeta_{k}^{2}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-x^{*}||^{2}+2\alpha_{k}\zeta_{k}(f(x^{*})-f(x_{k})+\bar{e}_{k})+\alpha_{k}^{2}\zeta_{k}^{2}||\overline{g}_{k}||^{2}$
$\displaystyle\leq$
$\displaystyle||x_{k}-x^{*}||^{2}+2\bar{e}_{k}\overline{\zeta}C_{2}/k+\alpha_{k}^{2}\overline{\zeta}^{2}G^{2}$
$\displaystyle\leq$
$\displaystyle||x_{0}-x^{*}||^{2}+2C_{2}\overline{\zeta}\sum_{k=0}^{\infty}\frac{\bar{e}_{k}}{k}+\overline{\zeta}^{2}G^{2}\sum_{k=0}^{\infty}\frac{C^{2}_{2}}{k^{2}}.$

Thus, there exists a constant $C_{5}$ such that $||x_{k}-x^{*}||\leq C_{5}$
for all $k$, which completes the proof. ∎

We summarize the convergence result for unbounded feasible set in the
following theorem.

######  Theorem 3.5.

Suppose that Assumption A1 holds, the feasible set $\Omega$ is convex and
closed and $\\{x_{k}\\}$ is a sequence generated by Algorithm SPS where
$N_{k}$ tends to infinity fast enough to provide
$\sum_{k=0}^{\infty}\bar{e}_{k}/k\leq C_{4}<\infty$ with $\bar{e}_{k}$ given
by (15). Then

$\liminf_{k\to\infty}f(x_{k})=f^{*}\text{a.s.}.$

Moreover,

$\lim_{k\rightarrow\infty}x_{k}=x^{*}\text{a.s.}$

for some $x^{*}\in X^{*}$ provided that $\sum_{k=0}^{\infty}e_{k}/k<\infty$,
where $e_{k}:=\max_{x\in\bar{\Omega}}|f_{\mathcal{N}_{k}}(x)-f(x)|$ and
$\bar{\Omega}$ is a compact set containing $\\{x_{k}\\}$.

###  3.1 Improving the efficiency - Line Search SPS

Notice that SPS algorithm works with an arbitrary subgradient direction
related to the current SAA function. However, in some applications such as
Hinge Loss binary clasification, it is possible to provide a descent direction
with respect to the SAA function [17] by applying the procedure proposed in
[31] or gradient subsampling technique [7] for instance. On the other hand, it
is well known that applying the line search may improve the performance of the
algorithm significantly, even in the stochastic environment. In order to make
the SPS algorithm more efficient, we propose a line search technique adapted
to the nonsmooth variable sample size framework to fit the SPS convergence
analysis. The proposed line search does not require a descent search direction
in order to be well defined, nor the convergence analysis depends on the
descent property. So, the following property (16) of $\bar{g_{k}}\in\partial
f_{\mathcal{N}_{k}}(x_{k})$ is desirable, but not necessary in order to prove
the convergence of the Line Search SPS (LS-SPS) algorithm presented in the
sequel,

$\sup_{g\in\partial
f_{\mathcal{N}_{k}}(x_{k})}g^{T}p_{k}\leq-\frac{m}{2}\|p_{k}\|^{2}\quad\mbox{for
some}\;m>0.$ (16)

The LS procedure. Since we employ the spectral subgradient method, we use
nonmonotone Armijo-type line search condition

$f_{\mathcal{N}_{k}}(x_{k}+\alpha_{k}p_{k})\leq\max_{i\in[\max\\{1,k-c\\},k]}{f_{\mathcal{N}_{i}}(x_{i})}-\eta\alpha_{k}||p_{k}||^{2},$
(17)

where $p_{k}$ is the search direction as in Step S1 of Algorithm 1. The
candidates for $\alpha_{k}$ that we consider are: $d_{k}$ and $(d_{k}+1/k)/2$,
where $d_{k}=\min\\{1,C_{2}/k\\}$. The reasoning behind this is the following.
The choice of $\alpha_{k}=1/k$ is a typical choice that is suitable for
obtaining a.s. convergence. The line search is employed to estimate if the
larger value of $\alpha_{k}$ may be used. Since the backtracking techniques
usually start with 1, we take the minimum of 1 and $C_{2}/k$ as the initial
choice. Although $C_{2}/k$ must be included to ensure the theoretical
requirements of Step S2, one can take $C_{2}$ arbitrary large such that
$d_{k}=1$ even for the large values of $k$. Thus, practically, 1 would be the
initial choice in all practical applications. We set the middle of the
interval $[1/k,d_{k}]$ as the second possible choice for step size in line
search. Although other strategies are feasible as well, we reduce to these two
guesses to avoid the computational costs of unsuccessful line search attempts.
Thus, if $\alpha_{k}=d_{k}$ satisfies (17), we take this as a step size. If
not, we check (17) with the medium value $\alpha_{k}=(d_{k}+1/k)/2$. If the
condition is satisfied, we retain this choice, otherwise we set
$\alpha_{k}=1/k$.

Algorithm 1: LS-SPS
(Line Search Spectral Projected Subgradient Method for Nonsmooth Optimization)

  * S0

Initialization. Given $N_{0}\in\mathbb{N},$ $x_{0}\in\Omega,$
$0<C_{1}<1<C_{2}<\infty,$ $0<\underline{\zeta}\leq\overline{\zeta}<\infty,$
$\zeta_{0}\in\left[\underline{\zeta},\overline{\zeta}\right],\eta\in(0,1),c\in\mathbb{N}.$
Set k = 0.

  * S1

Direction. Choose $\bar{g}_{k}\in\partial f_{\mathcal{N}_{k}}(x_{k})$
satisfying (16) if possible and set $p_{k}=-\zeta_{k}\bar{g}_{k}$.

  * S2

Step size. Choose $\alpha_{k}\in\\{d_{k},(d_{k}+1/k)/2\\}$ such that (17)
holds if possible. Otherwise set $\alpha_{k}=\frac{1}{k}.$

  * S3

Main update. Set $x_{k+1}=P_{\Omega}(x_{k}+\alpha_{k}p_{k})$ and
$s_{k}=x_{k+1}-x_{k}$.

  * S4

Sample size update. Chose $N_{k+1}\in\mathbb{N}$.

  * S5

Spectral coefficient update. Calculate
$y_{k}=g_{\mathcal{N}_{k}}(x_{k+1})-\bar{g}_{k}$ where
$g_{\mathcal{N}_{k}}(x_{k+1})\in\partial f_{\mathcal{N}_{k}}(x_{k+1})$. Set
$\zeta_{k+1}=\min\\{\overline{\zeta},\max\\{\underline{\zeta},\frac{s_{k}^{T}s_{k}}{s_{k}^{T}y_{k}}\\}\\}.$

  * S6

Set $k:=k+1$ and go to S1.

Remark. LS-SPS algorithm falls into the framework of SPS algorithm as
$\alpha_{k}$ satisfies the condition (17). Thus, the whole convergence
analysis presented for the SPS algorithm also holds for LS-SPS.

##  4 Numerical results

We performed preliminary numerical experiments on the set of binary
classification problems listed in Table 1. The problems are modeled by the
$L_{2}$-regularized Hinge Loss. More precisely, we consider the following
optimization problem for learning with a Support Vector Machine introduced in
[24]

$\min_{x\in\Omega}f(x):=10||x||^{2}+\frac{1}{N}\sum_{i=1}^{N}\max\\{{0,1-z_{i}x^{T}w_{i}\\}},$
$\Omega:=\\{{x\in\mathbb{R}^{n}\;:||x||^{2}\leq 0.1\\}},$

where $w_{i}\in\mathbb{R}^{n}$ are the input features and
$z_{i}\in\\{{1,-1\\}}$ are the corresponding labels. Thus, we have a convex
problem with the compact feasible set easy to project on. Moreover, for these
kind of problems it is possible to calculate the descent direction and we use
the procedure proposed in [31, Algorithm 2, p. 1155] as a subroutine that
provides the descent property (16). We employ this subroutine in all the
tested algorithms to ensure the fair comparison.

| Data set | $N$ | $n$ | $N_{train}$ | $N_{test}$
---|---|---|---|---|---
1 | SPLICE [28] | 3175 | 60 | 2540 | 635
2 | MUSHROOMS [14] | 8124 | 112 | 6500 | 1624
3 | ADULT9 [28] | 32561 | 123 | 26049 | 6512
4 | MNIST(binary) [29] | 70000 | 784 | 60000 | 10000
Table 1: Properties of the datasets used in the experiments.

Our numerical study has several goals. It is designed to investigate:
a) whether the variable sample size approach remains beneficial in the
nonsmooth environment with bounded full sample;
b) whether introducing the line search pays off;
c) whether the spectral coefficient improves the efficiency of the projected
subgradient method.

We set the experiments as follows. The main criterion for comparison of the
methods will be the computational cost modeled by $FEV$ (number of function
evaluations). More precisely, $FEV^{m}_{k}$ represents the number of scalar
products needed for method $m$ to calculate $x_{k}$ (starting from $x_{0}$).
We also track the value of the true objective function across the iterations
to observe the progress of the considered method.

To answer the question a), we compare the VSS methods to their full sample
counterparts. The extension $-F$ (e.g. SPS-F) will indicate that the full
sample is used at every iteration, i.e., $N_{k}=N$ for all $k$. On the other
hand, we assume the following sample size increase for the VSS methods:
$N_{k+1}=\lceil\min\\{1.1N_{k},N\\}\rceil$ with $N_{0}=0.1N$. Obviously there
are many other choices which can be more efficient, but we choose this simple
increase to be tested in the initial phase of the method evaluation. To
address the question b) we compare LS-SPS algorithm to SPS algorithm with the
standard choice of the step size $\alpha_{k}=1/k$. Finally, to address c), we
compare the proposed methods to the first order subgradient method denoted by
LS-PS (Line Search Projected Subgradient), which can be viewed as a special
case of LS-SPS with $\underline{\zeta}=\overline{\zeta}=1.$ We also test the
VSS and the full sample alternatives of the projected subgradient method: LS-
PS and LS-PS-F, respectively. The results of the subgradient method with the
choice of $\alpha_{k}=1/k$ were poor and thus not reported here.

The relevant parameters are as follows. The initial points are chosen randomly
from $(0,1)$ interval and we use the same initial point for all the tested
methods within one run. The step size parameters are
$C_{1}=10^{-2},C_{2}=10^{2}$ and the line search is performed with
$\eta=10^{-4}$ and $c=5$. For the proposed spectral methods, the safeguard
parameters are $\underline{\zeta}=10^{-4}$ and $\overline{\zeta}=10^{4}$.

We perform 5 independent runs for each of the data sets which yields 20 runs
of each method in total. A demonstrative run is presented in Figure 1 where
the objective function $f(x_{k})$ is plotted against the $FEV_{k}$. It reveals
that LS-SPS methods outperform other methods in a sense that reach tighter
vicinity of the solution. Even if we use the spectral coefficient,
predetermined step size was not enough to bring the sequence to the same
vicinity as obtained by the LS counterparts. On the other hand, observing the
LS-PS method which uses line search but without a spectral coefficient, we can
see that the line search itself (without second order information) was not
enough to push the subgradient method towards the solution. Finally, notice
that the computational cost is reduced significantly by employing the VSS
scheme in LS-SPS method.

Figure 1: MNIST data set

In order to compare the tested methods taking into account all the runs and
data sets, we employ the metric based on the ideas of the performance profile
[10] adapted to the stochastic environment in [18]. Roughly speaking, we
estimate the probability of winning for each of the tested methods. For the
considered run, the method $m$ wins if it reaches the vicinity $\tau$ of the
solution with the smallest costs. Since the theoretical stopping criterion is
non-existing, we stop the methods if the maximal number of scalar products is
achieved. At each iteration $k$ we measure the distance from the solution by
observing the relative error of method $m$ with respect to the optimal value,
i.e., $r^{m}_{k}:=(f(x^{m}_{k})-f^{*})/f^{*}$. For each method $m$ and each
run $l$ we register the first iteration $k(m,l)$ at which we have
$r^{m}_{k(m,l)}\leq\tau$ and read the corresponding $FEV^{m}_{k(m,l)}$. Then,
the method $m$ earns a point in run $l$ if
$FEV^{m}_{k(m,l)}=\min_{j}FEV^{j}_{k(j,l)}$. Finally, we estimate the
probability of winning, denoted by $\pi$ by

$\pi=\frac{t}{T},$

where $t$ is the number of earned points and $T$ is the total number of runs.
Notice that in the described situation we can have more winners, in other
words more methods can share the first place if they reached the goal with the
same costs.

The results are presented in Figure 2 for different relative errors
$\tau\in[0.01,3.5]$. They reveal that the VSS methods clearly outperform their
full sample counterparts and that LS-SPS method turns out to be the best
possible choice according to the conducted experiments. The algorithms LS-PS
and SPS reach 1 for very large values of the relative error $\tau$ which is
not relevant, so we do not show this part of the graph.

The Figure 3 represents classical performance profile (PP) graph for fixed
relative error $\tau=1$. The FEV is kept as the criterion for the classical PP
as well. On the $y$-axes we plot the probability that the method is close
enough to the best one, where ”close enough” is determined by the value on
$x$-axes denoted by $q$. More precisely, retaining the same notation as above,
the method $m$ earns a PP(q) point in run $l$ if $FEV^{m}_{k(m,l)}\leq
q\min_{j}FEV^{j}_{k(j,l)}$ and the plotted values correspond to

$\pi_{PP}(q)=\frac{t_{PP}(q)}{T},$

where $t_{PP}(q)$ is the number of earned PP(q) points for the considered
method. Again, from this figure it is clear that LS-SPS is the most robust,
i.e., it has the highest probability of being the optimal solver.

Figure 2: Empirical probabilities of winning ($\pi$) for different relative
errors ($\tau$). Figure 3: Performance profile for level of accuracy $\tau=1$

###  4.1 Additional comparison

In this subsection we show additional numerical results in order to compare
proposed algorithm with the proximal bundle method (PBM). It is known that PBM
gives the best result under a fixed number of iterations. The reason behind
this is the fact that the number of constraints in the quadratic program
solved by PBM may grow linearly with the number of iterations [15].
Accordingly, PBM may become significantly slower when the number of iterations
become larger. For that reasons, we compare fixed number of iterations of LS-
SPS-F with PBM that use full sample size in all iterations, and after that LS-
SPS with VSS PBM (where the sample is changed in the same way as in LS-SPS
through iterations).

In order to ensure the fair comparison, we choose several different
combinations of initial parameters for PBM. Table 2 summarizes the properties
of the observed methods, where $\gamma$ is proximity control parameter, $m$ is
descent coefficient, $\epsilon$ is tolerance parameter and $\omega$ is decay
coefficient. Detailed information about these parameters and implementations
in Matlab of PBM are available at [30].

The Figure 4 shows the results for Full and VSS version of LS-SPS and PBM
algorithms on MNIST data set, while the results on other three data sets are
similar. The objective function $f(x_{k})$ is plotted against the iteration
$k$ and the y-axes are in logarithmic scale. The results show that the
proposed algorithms (LS-SPS-F and LS-SPS) outperform the observed PBM
counterparts.

PBM | $1$ | $2$ | $3$ | $4$ | $5$
---|---|---|---|---|---
$\gamma$ | 1 | 1 | 1 | 1 | 1
$m$ | 0.01 | 0.1 | 0.01 | 0.01 | 0.01
$\epsilon$ | 0.1 | 0.1 | 0.01 | 0.1 | 0.1
$\omega$ | 0.5 | 0.5 | 0.5 | 0.9 | 0.1
Table 2: The initial parameters for PBM.

Figure 4: LS-SPS-F against Full PBM (left) and LS-SPS against VSS PBM (right).
MNIST data set.

##  5 Conclusions

The method SPS for solving constrained optimization problems with a nonsmooth
objective function in the form of mathematical expectation is proposed. We
assume that the feasible set is easy to project on and employ orthogonal
projections in order to ensure feasibility of the iterates. The SAA is used to
estimate the objective function and Variable Sample Size strategy is applied.
SPS combines an SAA subgradient with the spectral coefficient in order to
provide a suitable direction. The step sizes are chosen from the predefined
interval, while the choice is further specified in the line search version -
LS-SPS. The almost sure convergence of the method is proved under the standard
assumptions in stochastic environment. We provide convergence analysis both
for bounded and unbounded sample sizes. In the later case we show that the
sequence of iterates is bounded if the sample size growth is fast enough. We
also analyze the important special case - finite sum problems, for which we
prove deterministic convergence under the reduced set of assumptions.
Numerical experiments are conducted on the set of machine learning problems
modeled by the Hinge Loss binary classification. The methods are compared
through performance profile type of graphs where the main criterion is the
computational cost modeled by the number of scalar products. The initial
results yield several conclusions, including the one that VSS outperforms the
full sample as expected. We also conclude that the spectral coefficient is
beneficial, but it achieves the best performance when combined with the
Armijo-like line search procedure. Additionally, numerical results show clear
advantages of the proposed LS-SPS method with respect to Proximal Bundle
Method. Since the current results are promising, the future work should
include adaptive VSS schemes, inexact projections and further investigation of
the step size selection and the spectral coefficient in stochastic
environment.

Acknowledgement. We are grateful to the associate editor and the anonymous
referee for their comments that helped us to improve the paper.
Funding This work is supported by the Ministry of Education, Science and
Technological Development, Republic of Serbia.
Availability statement The datasets analysed during the current study are
available in the the MNIST database of handwritten digits [29], LIBSVM Data:
Classification (Binary Class) [28] and UCI Machine Learning Repository [14].

## Declarations

Conflict of interest The authors declare no competing interests.

## References

  * [1] P. Apkarian, D. Noll, L. Ravanbod, Nonsmooth bundle trust-region algorithm with applications to robust stability, Set-Valued and Variational Analysis 24(1) (2016), pp. 115 – 148. https://doi.org/10.1007/s11228-015-0352-5
  * [2] A. Bagirov, N.Karmitsa, M. Mäkelä, Introduction to Nonsmooth Optimization, Springer, (2014). https://doi.org/10.1007/978-3-319-08114-4
  * [3] F. Bastin, C. Cirillo, P.L. Toint, An adaptive Monte Carlo algorithm for computing mixed logit estimators, Comput. Manag. Sci. 3(1) (2006), pp. 55-79. https://doi.org/10.1007/s10287-005-0044-y
  * [4] S. Bellavia, N. Krejić, N. Krklec Jerinkić, Subsampled Inexact Newton methods for minimizing large sums of convex function, IMA Journal of Numerical Analysis 40(4) (2018), pp. 2309-2341. https://doi.org/10.1093/imanum/drz027
  * [5] S. Bellavia, N.K. Jerinkić, G. Malaspina, Subsampled nonmonotone spectral gradient methods, Communications in Applied and Industrial Mathematics, 11(1) (2020), pp. 19-34. https://doi.org/10.2478/caim-2020-0002
  * [6] E.G. Birgin, J.M. Martínez, M. Raydan, Nonmonotone Spectral Projected Gradients on Convex Sets, SIAM Journal on Optimization. 10 (2000) pp. 1196-1211. https://doi.org/10.1137/S1052623497330963
  * [7] J. Burke, A. Lewis, M. Overton, Approximating subdifferentials by random sampling of gradients, Mathematics of Operations Research 27(3) (2002), pp. 567-584. https://doi.org/10.1287/moor.27.3.567.317
  * [8] C. Lemarechal, C. Sagastizabal, Variable metric bundle methods: From conceptual to implementable forms, Mathematical Programming 76(3) (1997), pp. 393 – 410. https://doi.org/10.1007/BF02614390
  * [9] L. Lukšan, J. Vlček, A bundle-Newton method for nonsmooth unconstrained minimization, Mathematical Programming 83(1) (1998), pp. 373 – 391. https://doi.org/10.1007/BF02680566
  * [10] E. D. Dolan, J. J. Moré, Benchmarking optimization software with performance profiles, Math. Program., Ser. A 91 (2002), pp. 201-213. https://doi.org/10.1007/s101070100263
  * [11] M. A. Gomes-Ruggiero, J. M. Martínez, S. A. Santos, Spectral projected gradient method with inexact restoration for minimization with nonconvex constraints, SIAM Journal on Scientific Computing, 31(3) (2009), pp. 1628-1652. https://doi.org/10.1137/070707828
  * [12] W. Hare, C. Sagastizábal, M. Solodov, A proximal bundle method for nonsmooth nonconvex functions with inexact information, Computational Optimization and Applications, 63(1) (2016), pp. 1-28. https://doi.org/10.1007/s10589-015-9762-4
  * [13] T. Homem-de-Mello, Variable-Sample Methods for Stochastic Optimization, ACM Trans. Model. Comput. Simul. 13(2) (2003), pp. 108–133. https://doi.org/10.1145/858481.858483
  * [14] M. Lichman, UCI machine learning repository, https://archive.ics. uci.edu/ml/index.php, (2013).
  * [15] Q. Jiang, Proximal Bundle Methods and Nonlinear Acceleration: An Exploration, Senior Thesis, (2021).
  * [16] N. Krejić, N. Krklec Jerinkić, Spectral projected gradient method for stochastic optimization, Journal of Global Optimization, 73 (2018), pp. 59–81. https://doi.org/10.1007/s10898-018-0682-6
  * [17] N. Krejić, N. Krklec Jerinkić, T. Ostojić, Minimizing Nonsmooth Convex Functions with Variable Accuracy, arXiv:2103.13651v3, (2021).
  * [18] N. Krklec Jerinkić, A. Rožnjik, Penalty variable sample size method for solving optimization problems with equality constraints in a form of mathematical expectation, Numerical Algorithms, 83 (2020), pp. 701-718. https://doi.org/10.1007/s11075-019-00699-6
  * [19] M. Mäkelä, Survey of bundle methods for nonsmooth optimization, Optimization methods and software, 17(1) (2002), pp. 1-29. https://doi.org/10.1080/10556780290027828
  * [20] M. Mäkelä, P. Neittaanmaki, Nonsmooth Optimization: Analysis and Algorithms with Applications to Optimal Control, World Scientific Publishing Co. (1992). https://doi.org/10.1142/1493
  * [21] K. Miettinen, Nonlinear Multiobjective Optimization, Springer, (1998). https://doi.org/10.1007/978-1-4615-5563-6
  * [22] R. Mifflin, A modification and an extension of Lemarechal’s algorithm for nonsmooth optimization, Mathematical Programming Studies 17 (1982), pp. 77 – 90. https://doi.org/10.1007/BFb0120960
  * [23] W. D. Oliveira, C. Sagastizábal, Bundle methods in the XXIst century: A bird’s-eye view, Pesquisa Operacional, 34 (2014), pp. 647-670. https://doi.org/10.1590/0101-7438.2014.034.03.0647
  * [24] S. Shalev-Shwartz, Y. Singer, N. Srebro, A. Cotter,  Pegasos: primal estimated sub-gradient solver for SVM, Mathematical programming, 127(1) (2011), pp. 3-30. https://doi.org/10.1007/s10107-010-0420-4
  * [25] C. Sagastizabal, M. Solodov, An infeasible bundle method for nonsmooth convex constrained optimization without penalty function or a filter, SIAM Journal on Optimization, 140(1) (2005), pp. 146 – 169. https://doi.org/10.1137/040603875
  * [26] A. Shapiro, D. Dentcheva, A. Ruszczynski, Lectures on Stochastic Programming: Modeling and Theory, MPS-SIAM Series on Optimization (2009). https://doi.org/10.1137/1.9780898718751
  * [27] C. Tan, S. Ma, Y. H. Dai, Y.Qian, Barzilai-borwein step size for stochastic gradient descent, Advances in neural information processing systems, 29, (2016). https://doi.org/10.48550/arXiv.1605.04131
  * [28] https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/binary.html
  * [29] http://yann.lecun.com/exdb/mnist/
  * [30] https://github.com/ritchie-xl/Bundle-Method-Matlab/blob/master/bundle_method.m
  * [31] J. Yu, S. Vishwanathan, S. Guenter, N. Schraudolph, A Quasi-Newton Approach to Nonsmooth Convex Optimization Problems in Machine Learning, Journal of Machine Learning Research 11 (2010), pp. 1145-1200.


# A Quantum Computing Approach for Multi-robot Coverage Path Planning

Poojith U Rao, Florian Speelman, Balwinder Sodhi, Sachin Kinge

###### Abstract

This paper tackles the multi-vehicle Coverage Path Planning (CPP) problem,
crucial for applications like search and rescue or environmental monitoring.
Due to its NP-hard nature, finding optimal solutions becomes infeasible with
larger problem sizes. This motivates the development of heuristic approaches
that enhance efficiency even marginally.

We propose a novel approach for exploring paths in a 2D grid, specifically
designed for easy integration with the Quantum Alternating Operator Ansatz
(QAOA), a powerful quantum heuristic. Our contribution includes:

  * •

An objective function tailored to solve the multi-vehicle CPP using QAOA.

  * •

Theoretical proofs guaranteeing the validity of the proposed approach.

  * •

Efficient construction of QAOA operators for practical implementation.

  * •

Resource estimation to assess the feasibility of QAOA execution on current
hardware.

  * •

Performance comparison against established algorithms like the Depth First
Search, demonstrating significant runtime improvements.

This work paves the way for leveraging quantum computing in optimizing multi-
vehicle path planning, potentially leading to real-world advancements in
various applications.

###### Index Terms:

Quantum alternating Operator Ansatz, Multi-vehicle Coverage Path Planning, NP-
Hard problem

##  I Introduction

Advancements in robotics, drones, sensors, and navigation technology have led
to the widespread adoption of multi-vehicle systems for numerous applications.
The common objectives include surveillance, monitoring, surveying, and
modeling in landmine detection, lawn mowing, search and rescue operations, and
natural disaster detection. Such objectives can be achieved more efficiently
when multiple vehicles communicate and cooperate. The common aim of such
applications is to scan an area of interest while optimizing various objective
functions. Finding an optimal and collision-free route to cover an area of
interest is the Coverage Path Planning (CPP) problem.

The two main components of CPP include i) viewpoint generation and ii) path
generation [1]. Viewpoint generation refers to identifying important positions
in the area to be covered. The viewpoints help in collecting data that aid in
covering the entire area. Most of the algorithms in the literature take a
uniform distribution of viewpoints across the Area of Interest (AoI). A simple
example is the 2D grid-based approach, where viewpoints are modeled as a grid
with different shapes, such as triangles, squares, hexagons, etc. A more
detailed study on different viewpoint generation approaches is discussed in
[1].

Path generation refers to finding an optimal route that covers these
viewpoints. This step starts with dividing the AoI into sub-areas and
initializing the robot’s location. Then, the covering direction is determined
for each sub-area. When this is known, each sub-area is covered using simple
movements. The efficacy of the path identified can be measured in terms of
numerous parameters such as the percentage of covered area, path overlap rate,
and energy consumption of robots [2].

CPP can be viewed as a combinatorial optimization problem with multiple
objectives and constraints depending on its application. Commonly used
objective functions include the number of rotations or turns in the path,
number of viewpoints covered, computing time, path length, path overlap,
energy consumption, smooth trajectories/maneuvers, etc., and constraints
include avoiding obstacles, no backtracking, a continuous path, usage of
simple trajectories such as straight lines and circles, etc.

Developing new algorithms for solving combinatorial optimization problems with
hard constraints, such as CPP, that can improve the runtime performance even
by a small fraction is of great importance. It is because, at scale, even a
slight improvement in resource optimization accrues significantly in a large-
scale real-world application scenario.

###  I-A Related work

Multiple works aim to solve the Coverage path planning using a single robot.
In the case of a large area of interest (AoI), single robot solutions are not
preferred due to the high possibility of malfunction, mechanical or electronic
breakdown, sensor and actuator faults, and battery drainage. Thus, recent
research uses multiple robots to overcome most of these shortcomings [3].
Numerous challenges associated with the CPP problem along with an exhaustive,
comprehensive review of existing algorithms used for solving it exist in
literature [2]. Some of the relevant works where a multi-vehicle version of
the CPP problem are discussed here.

Spanning trees have been used in numerous works for solving the multi-vehicle
version of CPP [4, 5, 6]. They can be used as a base for creating coverage
paths and a polynomial-time tree construction algorithm is used to
dramatically improves the coverage time [5]. However, the traveled path in
this algorithm depends on the initial positions of each robot, which might
lead to issues such as backtracking among other robots, a high overlap rate,
significantly smaller energy efficiency.

Artificial Potential Fiend (APF) method has also been used to cover the area
and find collision-free paths for the multi-robot system [7]. Despite all the
research efforts, there is still a lack of planning for collision avoidance
between multiple robots when simultaneously accessing the goal under the
potential field.

Dijkstra’s algorithm is also proposed for task planning in a multi-robot
system. Nevertheless, the search path is not optimal in terms of travel
distance [8, 9]. Theta $A^{*}$ algorithm was used in [10] to reduce the
coverage time. However, the algorithm failed to generate a global optimization
solution in terms of path length. The Voronoi partition approach is another
common modeling technique applied in the distributed coordination for a multi-
robot system [11].

Genetic Algorithm [3], Bee Colony Optimization (BCO) [12], Firefly algorithm
(FA) [13, 14], Coordinated Multi-robot Exploration (CME), Gray Wolf
Optimization (GWO) algorithm [15] and numerous other evolutionary algorithm
based solutions have been proposed for the CPP. The performance of FA, PSO,
and BCO in coordinating the swarm robotics system in terms of energy
consumption was compared by [16]. The prevalence of heuristic approaches for
the CPP problem underscores their real-world practicality and value. Hence, in
this work, heuristic based approach capable of taking advantage of quantum
mechanics is proposed.

The multi-vehicle version of the Coverage Path Planning problem is NP-hard.
Thus, mostly heuristic-based [17] solutions are proposed that can neither
guarantee optimality nor shorten running times. There is a recent trend [17]
in using quantum algorithms on near-term quantum devices that promise
potential quantum advantage for solving combinatorial problems. While more
attention is given to quantum annealing due to the near-term availability of
hardware, recent works also show the applicability of gate-based quantum
computers for solving such problems. In particular, for solving problems with
hard constraints on gate-based computers, the quantum alternating operator
ansatz (QAOA) algorithm was proposed by [18].

In this work, an approach to explore paths between two points in a 2D grid is
proposed. The approach involves starting from a trivial initial path (random
state) and exploring nearby solutions until a satisfactory solution is found
within the feasible solution space. The methodology proposed to explore new
solutions can be adapted to other well-known heuristic algorithms. For
example, it can be implemented as mutations/crossover in Genetic Algorithms
(GA), random neighbors in Simulated Annealing (SA), and mixers in Quantum
Alternating Operator Ansatz (QAOA).

In this work, the focus is on the Quantum Alternating Operator Ansatz (QAOA),
which is a meta-heuristic framework proposed by [18] used to perform
approximate optimization on gate-based quantum computers. This is an extended
version of the Quantum Approximation Optimization Algorithm proposed by [19]
that produces approximate solutions for combinatorial optimization problems.
The version proposed in [18] is well suited for problems such as CPP, where
the feasible space is smaller than the complete solution space. The QAOA
algorithm can be summarized as follows:

  1. 1.

Set the quantum system to an initial state, which is trivial to implement and
lies inside the feasible subspace of the problem.

  2. 2.

Apply suitable parameterized quantum operators $p$ times on the initial state.

  3. 3.

Measure the system in the computational basis to find a candidate solution to
the problem.

  4. 4.

Update the parameters of the operators to optimize the cost function.

  5. 5.

Repeat 2-3-4 until desired convergence is obtained.

  6. 6.

Measure in the computational basis to find a solution to the problem.

In the following sections, the methodology to explore paths in a 2D grid is
described. The methodology is described in terms of a Genetic Algorithm.
Later, the design of mixer operators and phase separation operators for the
QAOA algorithm, along with approximate resource estimates, are discussed. The
proposed methodology is evaluated by implementing the Simulated Annealing and
QAOA algorithm and comparing it with the DFS algorithm.

##  II Coverage Path Planning problem

In this section, an approach based on Genetic Algorithms (GAs) is proposed for
exploring all paths between two nodes in a 2D grid. As with any GA, the
approach starts with an initial population and explores new feasible paths by
selecting individuals from the current population to be parents and produce
children for the next generation. An objective function evaluates the quality
of the population. Over successive generations, the population evolves towards
an optimal solution. The important steps of the algorithm are described in the
following subsections.

###  II-A Problem definition

The Area of Interest (AOI) is defined by a rectangular grid of dimensions
$n\times m$ and characterized by a set $E$ of edges. There are $r$ robots,
denoted as $r_{i}$, each initially stationed at $s_{r_{i}}$ (source station).
This can be visualized as a cube-shaped graph of size $r\times n\times m$. The
robots navigate the AOI along the grid edges from the source node to the
destination node. Every node within the grid serves as a viewpoint and is to
be covered by the robots. The path of a robot is represented by a sequence of
edges from the source node to the destination node.

The objective of the multi-vehicle CPP is to allocate paths to the robots with
the following criteria: i) maximize coverage of viewpoints, ii) minimize the
total time required by all robots for coverage, and iii) minimize redundant
visits to a node by different robots.

The paths assigned to each robot must also adhere to the following conditions:
i) avoid backtracking while exploring viewpoints, ii) be continuous and
sequential, and iii) navigate around obstacles on the map.

####  II-A1 Decision variable

Robots move along the edges of the graph to cover all viewpoints. Each edge in
the grid is considered as a decision variable (Equation 1). The variable
represents if an edge is covered by a robot. In a grid of size $n\times m$,
there are $m(n-1)+n(m-1)$ edges. Hence, the total number of variables needed
is $r(m(n-1)+n(m-1))$.

$x_{i,e}=\begin{cases}1&\text{if the robot i covers edge e}\\\
0&\text{otherwise}\\\ &i=1,\cdot\cdot\cdot,r;e\in E;\\\ \end{cases}$ (1)

####  II-A2 Objective function

The objective function has to satisfy four requirements: R1) avoids obstacles
along the path, R2) maximizes the number of viewpoints covered, R3) minimizes
the overall cover time needed and R4) minimize overlapping of paths assigned
to robots.

Obstacles along the path can be avoided by avoiding all the edges leading
towards it. To satisfy R1) and R2), a weight is assigned to each edge in the
graph. A positive weight is assigned to all edges leading to obstacles, and a
negative weight is assigned to the remaining edges. When the objective
function is minimized, the positively weighted edges are preferred over
negatively weighted edges. The cost function for R1) and R2) is given in
Equation 2 where $w_{i}$ is the weight assigned to $i^{th}$ edge.

$c_{1}=\sum_{r,i}w_{i}x_{r,i}$ (2)

The total time the robots take to travel the assigned paths is the cover time.
An optimal time coverage algorithm for a system with $r$ robots and $m\times
n$ viewpoints will result in a total coverage time of $\frac{m\times n}{r}$
[5].

When robots have equal speed, they cover equal distances in equal time. Hence,
under the assumption of equal speed, it is sufficient to minimize the
difference in distance covered between all pairs of robots. The cost function
to compute the sum of all differences is given in Equation 3 where $d_{j}$ is
the length of edge $j$.

$c_{2}=\sum_{r}\left[\sum_{j}d_{j}x_{r,j}-\sum_{j}d_{j}x_{r+1,j}\right]^{2}$
(3)

Two edges are required to cover a node: one incoming edge and one outgoing
edge. Hence, the sum of all decision variables representing edges connected to
a node (which includes all robots) must be two (apart from source and
destination nodes). This is formulated as a cost function in Equation 4 where
$j$ denotes all edges connected to node $i$.

$c_{3}=\sum_{i\in\text{node}}\left[\sum_{j\in\text{i}}\sum_{r}x_{r,j}-2\right]^{2}$
(4)

The final objective function is the weighted sum of all the cost functions.
This is given in Equation 5.

$c=\alpha_{0}c_{1}+\alpha_{1}c_{2}+\alpha_{2}c_{3}$ (5)

###  II-B Initial population

The algorithm starts with an initial population consisting of paths for each
robot from the source node to the destination node. A heuristic to generate
the initial path is to start from the source node and move horizontally to the
column of the destination node and then vertically to the row of the
destination node.

Figure 1: On the left - Scenarios where SBF operation searches for paths
between two points. Top-left is the $S_{1}$ operation, and bottom-left is the
$S_{2}$ operation. On the right - Scenarios where SBF operation must be
avoided Figure 2: Scenario where SBF operation on a four-node sub-grid results
in a loop

###  II-C Mutation

The mutation operators are designed to explore all feasible solutions for the
problem. A feasible solution is one where the path assigned to each robot
meets the following conditions: i) It connects the source and destination
nodes. ii) There is no backtracking involved. iii) An edge/node is covered at
most once. iv) No loops are present. Mutations are performed only on a four-
node sub-grid at a time. The following subsections show how mutations explore
paths in a four-node grid and how repeated mutations on different four-node
sub-grids explore all paths in the grid.

####  II-C1 Four-node graph

A simple four-node graph is shown in Figure 1. It has four edges $ab$, $bc$,
$cd$, and $ad$ denoted by decision variables $x_{ab},{x_{bc}},{x_{cd}}\text{
and }{x_{da}}$. There are exactly two paths between any two nodes in this
graph. The two paths between two nodes are complements of each other. For
example, if $x_{ab}{x_{bc}}{x_{cd}}{x_{da}}$ is a valid path from $a$ to $b$
then $\overline{x_{ab}}\ \overline{x_{ab}}\ \overline{{x_{bc}}}\
\overline{{x_{cd}}}\ \overline{{x_{da}}}$ is also a valid path from $a$ to
$b$. Hence, if a bit string $x_{ab}{x_{bc}}{x_{cd}}{x_{da}}$ denotes a
feasible path between two nodes in a four-node graph, the other feasible path
is obtained by performing a NOT operation on the bit string. This operation of
simultaneously flipping all the edges in a four-node graph is termed as
Simultaneous Bit Flip (SBF) operation.

####  II-C2 Scaling to larger grids

All feasible paths between two nodes in an $m\times n$ grid can be explored by
applying SBF operation on smaller $(m-1)(n-1)$ four-node sub-grids. On
applying an SBF operation on a four-node grid, while the path inside the four-
node grid remains valid, the path outside the four-node grid can become
invalid (a discontinuous or a loopy path). Hence, specific conditions must be
satisfied before an SBF operation can be applied to a four-node sub-grid.
Examples of valid and invalid scenarios are shown in Figure 1. In the
subsequent sections, the logic function to identify valid scenarios is
discussed.

####  II-C3 Logic function representing valid/invalid paths inside four node
sub-grid

The logic function to identify valid scenarios inside a four-node sub-grid
depends on the following scenarios:

1\. There is no path in the sub-grid: Applying an SBF operation on a four-node
grid with no active edges (Figure 1 \- right bottom) forms a loop. Hence, if
$x_{abcd}={0000}$, the operation must not be applied. The logic expression for
the same is $f_{1}(a,b,c,d)=x_{ab}\lor x_{bc}\lor x_{cd}\lor x_{ad}$.

2\. There are two paths in the sub-grid: Applying an SBF operation on a four-
node grid with two parallel paths (Figure 1 \- right top) alters the source
and destination of the paths. Hence if
${x_{abcd}}={0101}\lor{x_{abcd}}={1010}$ the operation must not be applied.
The logic expression for the same is
$f_{2}(a,b,c,d)=\overline{x_{ab}\overline{x_{bc}}x_{cd}\overline{x_{ad}}\lor\overline{x_{ab}}x_{bc}\overline{x_{cd}}x_{ad}}$.

3\. There is a single path in the sub-grid: Applying the SBF operation on a
four-node sub-grid in this scenario may create loops in its adjacent sub-
grids. Hence, additional cases need to be verified. If the path passes through
any of the nodes in the four-node grid but does not enter it, then applying
the SBF operation on the four-node grid crisscrosses two paths. An example
covering this scenario is shown in Figure 2. Let $oe$ denote the edges of the
nodes in the sub-grid but not in the sub-grid. If any of the $oe$ edges are
active, the operation must not be applied on the sub-grid. The logical
expression for the same is
$f_{3}(a,b,c,d)=\overline{x_{a,oe_{1}}x_{a,oe_{1}}\lor
x_{b,oe_{3}}x_{b,oe_{4}}\lor x_{c,oe_{5}}x_{c,oe_{6}}\lor
x_{d,oe_{7}}x_{d,oe_{8}}}$.

The final logic function that decides if the operation can be applied on the
four-node sub-grid is the binary and operation of all the logical expressions.
This is given as

$f(a,b,c,d)=f_{1}(a,b,c,d)\land f_{2}(a,b,c,d)\land f_{3}(a,b,c,d)$ (6)

Figure 3: Exploring all feasible paths between two points $a$ and $b$ by
applying SBF operation on 4-node sub-grids.

The order of mutations controls the exploration of paths. As the population
evolves with more mutations, more feasible solutions are explored. The proof
that SBF operations explore all feasible solutions on four-node subgrids is
discussed in the next section. An example exploration in a 3x3 grid is shown
in Figure 3.

The pseudo-code of the GA explained above is as follows:

  1. 1.

Initial population: Start with bit strings that assign a valid path to each
robot.

  2. 2.

Compute the quality of the population using the objective function.

  3. 3.

Until a desired quality is achieved, repeat the following:

    1. (a)

Mutation: Apply the SBF operation on all four-node sub-grids and add the
offspring to the population.

    2. (b)

Evaluate the quality of the population using the objective function.

The above methodology can be adopted into the Simulated Annealing (SA)
algorithm by using the SBF operation to search for nearby neighbors.
Similarly, for QAOA, the SBF operation can be incorporated into the mixer
design, and the objective function can be incorporated into the phase
operators. The design and construction of these operators are discussed in the
subsequent sections.

##  III Proof of explorability

(a) Corner where SBF operation can be applied in a closed region.

(b) Corner where SBF operation cannot be applied in a closed region

Figure 4: Corner in a closed region

A path between two nodes is trivial if it only moves (either horizontally or
vertically via the edges) in the direction of the destination node. Trivial
paths are also the shortest paths from the source node to the destination
node.

Non-trivial paths move towards and away from the destination (i.e., reverse
their direction) multiple times before reaching the destination node. A path
reverses its direction by turning right or left twice consecutively. In such
cases, a straight line along the edges can be drawn to connect two nodes on
the path such that a closed region is formed. A non-trivial path can be
transformed into a trivial path by removing all such closed regions. The
procedure to reduce a closed region to a straight line is as follows:

While a closed region exists, perform the following steps.

  1. 1.

Apply SBF operation $S_{1}$ to all four-node grids inside the closed region.

  2. 2.

Apply SBF operation $S_{2}$ to all corners in the closed region wherever it is
valid.

In the above procedure, both steps reduce the area inside the closed region.
The $S_{1}$ and $S_{2}$ operations are shown in the Figure 1. The above
procedure fails to produce an output only when no valid operations can be
performed. The following argument shows that there always exists valid
operations $S_{1}$ or $S_{2}$ till a closed region is present.

In both operations $S_{1}$ and $S_{2}$, corners are considered for reducing
the area. Let $G$ (refer to Figure 4 ) be a corner inside the closed region
where the SBF operation is performed. The four-node grid under consideration
is $DEGH$. The border of the region through $G$ passes through nodes $B$ and
$H$ (shown as a solid black line in Figure 4), after which it can continue in
any direction (dotted black lines in Figure 4(a)). In either direction, if the
border passes through $E$, the operation $S_{1}$ can be applied.

If the border does not pass through $E$ and if there is no other segment of
the border passing through $E$, then the $S_{2}$ operation operation can be
applied to reduce the area. The SBF operation $S_{2}$ becomes invalid on the
four-node grid $DEGH$ if there is another segment of the border passing
through $E$.

If there is another border segment through $E$, it must pass through nodes $B$
and $F$. As the border is continuous and encloses a closed region, $D$ must
meet $B$ and $H$ must meet $F$. In both cases, there has to be a corner where
they can meet. The $S_{2}$ operation cannot be invalid because the borders
have to meet. Hence, there has to be a corner where the SBF operation is
valid. The above procedure ends when no area is left to reduce, i.e., only a
straight connecting two nodes is left.

With the above procedure, any nontrivial path can be transformed into a
trivial path. Also, any trivial path can be transformed into any other one by
applying multiple $S_{2}$ operations. As the SBF operations are reversible,
any trivial path can be transformed into a non-trivial path. Hence, any valid
path can be transformed into any other one, allowing all paths to be
explorable by applying SBF operations.

##  IV QAOA for Coverage Path Planning

In this section, the constructions of QAOA for the multi-vehicle CPP problem
are discussed. The techniques of [18] are used to design suitable problem
encodings and mixing operators such that the probability amplitude is
restricted to quantum states encoding feasible solutions. The circuits yield
parameterized quantum states that represent possible solutions for multiple
robots (via repeated state preparations and measurements). Measuring such a
state in the computational basis returns paths for each robot. The methodology
for designing initial state, mixing, and phase separation operators is
presented in the following sections.

###  IV-A Initial state

The initial state is created by applying basis encoding on all bit strings in
the initial population. Thus, the initial state is a superposition of all bit
strings in the initial population (refer to Section II-B).

###  IV-B Mixers

Mixers explore feasible solutions for the problem at hand. Hence, mixers that
can imitate the mutation process discussed earlier are designed. Specifically,
a partial mixer that performs SBF operation and a full mixer for applying SBF
operation on all four-node graphs are needed. As the operation must be applied
in valid scenarios, controlled-unitary partial mixers[17] are preferable. The
design of mixers is discussed in the following sub-section:

@C=1em @R=1em &H 1 1 H
H 1 1 H
H 1 1 H
H R_Z(β) H

Figure 5: Simultaneous Bit Flip mixer that flips all the qubits simultaneously
and is controlled by a parameter $\beta$.

####  IV-B1 Simultaneous Bit Flip mixer

A simultaneous Bit Flip (SBF) mixer performs an SBF operation defined earlier.
It has a trainable parameter $\beta$ that controls the probability of the
operation being applied. The circuit of the SBF mixer is shown in Figure 5.
The SBF mixer is defined as follows:

$U_{ab,bc,cd,da}^{SBF}(\beta)=exp(-i\beta(X_{ab}X_{bc}X_{cd}X_{da})/2)$ (7)

####  IV-B2 Partial controlled four-node mixer $U_{a,b,c,d}^{4N}$

A partial mixer $U_{a,b,c,d}^{4N}$ is designed which is the controlled SBF
mixer with the logic function $f(a,b,c,d)$ as the controlling logic. This is
defined as follows:

$U_{a,b,c,d}^{4N}(\beta)=\wedge_{f(a,b,c,b)}U_{ab,bc,cd,da}^{SBF}(\beta)$ (8)

Here, $U_{ab,bc,cd,da}^{SBF}(\beta)$ denotes the SBF mixer, and $f(a,b,c,d)$
is the controlling logic function.

####  IV-B3 Full mixer design

The full mixer for a grid of size $m\times n$ is defined as the product of all
partially controlled $U_{a,b,c,d}^{4N}$ Unitaries over all four-node grids
present. This is defined as follows:

$U_{CPP}(\beta)=\prod_{i}^{(m-1)(n-1)}U_{a,b,c,d}^{4N}(\beta)$ (9)

As it is possible to transform a random feasible path into any other feasible
path, the full mixer can cover the entire feasible search subspace. The
Boolean condition for each controlled mixer ensures that the new states
explored are always inside the feasible subspace. Hence, the full mixer
satisfies the conditions required by a mixing operator mentioned in [18].

###  IV-C Phase separation operator

The phase separation operator is encoded as a diagonal cost Hamiltonian $C$
that acts on state $x$ as $C\ket{x}=c(x)\ket{x}$. The problem cost function
(see Equation-5) is mapped to the Pauli Z operator using the transformations
described in [20]. The relevant equations after applying the transformations
are shown below.

$C=\sum_{r,i}w_{i}(1-Z_{r,i})+\\\
\sum_{r}\left[\sum_{j}d_{j}(1-Z_{r,j})-\sum_{j}d_{j}(1-Z_{r+1,j})\right]^{2}+\\\
\sum_{r}\sum_{j}(1-Z_{r,j})(1-Z_{r+l,j})$ (10)

This can be written as Equation 11 and implemented using $R_{Z}$ and $R_{ZZ}$
gates.

$C=c_{0}I-\sum_{r,i}c_{1}Z_{ri}+\sum_{r,i}c_{2}Z_{r,i}Z_{r+1,j}$ (11)

###  IV-D Resource estimates

Figure 6: A schematic diagram of a partial QAOA mixer circuit

As optimal circuit compilation is hardware-dependent, approximate resource
estimates of the proposed QAOA solution are shown here. CNOT and single-qubit
gates are considered universal gate sets[21]. We show an approximate resource
estimate in terms of the number of single-qubit gates $\mathbb{N}_{S}$ and the
number of CNOT gates $\mathbb{N}_{C}$. With simple circuit optimizations, the
resource estimates can be further reduced.

The problem formulation requires $r\times e$ qubits where $e=(m(n-1)+n(m-1))$
is the number of edges in the grid, and $r$ is the number of robots. We also
use eight ancillary qubits per robot to compute the logic function controlling
the mixer. As we reuse the ancillary qubits for all controlled partial mixers,
the number of qubits required is significantly reduced but at the cost of
increased circuit depth. Hence, the number of qubits needed is
$\mathbb{N}_{Q}\approx re$.

The resource estimates for the phase separation operator are as follows. Each
linear term in Equation 11 is implemented using an RZ gate requiring a total
of $re$ single-qubit gates. Similarly, each quadratic term is implemented
using 2 CNOT gates and 1 RZ gate, requiring a total of $re(re+1)$ CNOT gates
and $re(re+1)/2$ single-qubit gates.

A Full mixer consists of $r(m-1)(n-1)$ partial mixers. The schematic circuit
diagram of a partial mixer is shown in Figure 6. The controlled Z operation
used in the partial mixer can be decomposed into 2 Z rotations and 2 CNOT
gates [22]. $f_{1}(a,b,c,d)$ needs one 5-qubit Tofolli gate, $f_{2}(a,b,c,d)$
needs two 5-qubit Tofolli gates, and $f_{3}(a,b,c,d)$ needs four 3-qubit
Tofolli gates. A 7-qubit Toffoli gate is used for computing $f(a,b,c,d)$ on
the ancilla qubit $anc_{6}$, which controls the partial mixer (SBF mixer).

An $(l+1)$-qubit Tofolli gate can be decomposed into $4(l-2)$ 3-qubit Tofolli
gates using $l-2$ work qubits, when $l\geq 3$ [22]. Finally, each 3-qubit
Tofolli gate can be decomposed into 6 CNOT and nine single-qubit gates. The
ancillary reuse reduces the number of qubits required but requires recomputing
all the functions. Hence, the total gate requirement is twice of the
controlled partial mixer.

##  V Experiments

The proposed approach was verified by implementing the QAOA and Simulated
Annealing algorithm. The results were compared with the Depth First Search
(DFS) algorithm. DFS algorithm finds the best solution by performing an
exhaustive search on the entire search space. Due to computational
constraints, the QAOA algorithm could be tested for a single robot on a 3x3
grid. Hence, only the SA algorithm was tested for multiple robots on larger
grids.

###  V-A Collision free path for a single robot on 3x3 grid

(a) Initial state

(b) Optimal path

Figure 7: Initial state and optimal path as computed by QAOA, DFS, and SA

In this experiment, it was verified if the QAOA algorithm can find an optimal
path for a single robot on a 3x3 grid. The QAOA algorithm was simulated using
the Pennylane library. The mixer circuit in the QAOA is independent for each
robot. Hence, verifying if the mixer circuit can explore paths for a single
robot is sufficient to ensure that it can explore paths for multiple robots.

Twenty qubits (12 decision variable qubits + 8 ancillary qubits) were needed
to simulate the QAOA algorithm for this case. The parameters were optimized
using gradient descent with the Nesterov momentum algorithm (keeping step size
= 0.1). The initial state and the optimal path returned by the QAOA algorithm
are shown in Figure 7(a) and Figure 7(b). The square nodes represent the
endpoints, while the triangular nodes represent the obstacles.

(a) QAOA cost plot.

(b) SA cost plot.

Figure 8: Optimization plots

The quality of the solution returned by QAOA can be improved by increasing the
number of layers $p$. But this also increases the circuit depth and the number
of training parameters. The loss curve for different numbers of layers $p$ is
plotted in Figure 8(a). Due to the simplicity of the scenario, the optimal
solution is found with just a single layer. The result is verified to be the
same as the result returned by DFS and SA algorithm.

###  V-B Collision-free path for two robots

(a) Optimal path for robot 1

(b) Optimal path for robot 2

Figure 9: Optimal paths for two robots on 4x4 grid

The number of qubits needed to run QAOA algorithm for two robots on a 3x3 grid
($\sim\,$40) is already high enough to be simulated classically. Hence, only
the SA algorithm was tested and compared with the DFS algorithm for multiple
robots. The solutions evaluated by SA for an example scenario of two robots in
a 4x4 grid and 5x5 grids are shown in Figure 9 and Figure 10. The SA loss
curve for the 4x4 scenario is shown in Figure 8(b). These are the same results
returned by the DFS algorithm, thus verifying the results.

(a) Optimal path for robot 1

(b) Optimal path for robot 2

Figure 10: Optimal paths for two robots on 5x5 grid

##  VI Discussion and conclusion

A novel method for exploring paths in a 2D grid was presented in this work.
The methodology was designed such that it could be adapted to quantum
heuristic algorithms such as QAOA and classical algorithms such as
Evolutionary Algorithms. The circuit constructions for the QAOA algorithm were
presented along with approximate resource estimates. The proposed methodology
was evaluated by implementing QAOA and SA algorithms and compared with the DFS
algorithm.

This work may be extended by formulating this problem as a Quadratic
Unconstrained Binary Optimization (QUBO) problem and comparing the results
with annealer-based solutions.

## References

  * [1] R. Almadhoun, T. Taha, L. Seneviratne, and Y. Zweiri, “A survey on multi-robot coverage path planning for model reconstruction and mapping,” SN Applied Sciences, vol. 1, no. 8, pp. 1–24, 2019.
  * [2] C. S. Tan, R. Mohd-Mokhtar, and M. R. Arshad, “A comprehensive review of coverage path planning in robotics using classical and heuristic algorithms,” IEEE Access, 2021.
  * [3] R. Sun, C. Tang, J. Zheng, Y. Zhou, and S. Yu, “Multi-robot path planning for complete coverage with genetic algorithms,” in International Conference on Intelligent Robotics and Applications, pp. 349–361, Springer, 2019.
  * [4] K. Senthilkumar and K. Bharadwaj, “Spanning tree based terrain coverage by multi robots in unknown environments,” in 2008 Annual IEEE India Conference, vol. 1, pp. 120–125, IEEE, 2008.
  * [5] N. Agmon, N. Hazon, and G. A. Kaminka, “Constructing spanning trees for efficient multi-robot coverage,” in Proceedings 2006 IEEE International Conference on Robotics and Automation, 2006. ICRA 2006., pp. 1698–1703, IEEE, 2006.
  * [6] G.-Q. Gao and B. Xin, “A-stc: auction-based spanning tree coverage algorithm formotion planning of cooperative robots,” Frontiers of Information Technology & Electronic Engineering, vol. 20, no. 1, pp. 18–31, 2019.
  * [7] C. Huang, W. Li, C. Xiao, B. Liang, and S. Han, “Potential field method for persistent surveillance of multiple unmanned aerial vehicle sensors,” International journal of distributed sensor networks, vol. 14, no. 1, p. 1550147718755069, 2018.
  * [8] R. Rosa, T. Brito, A. I. Pereira, J. Lima, and M. A. Wehrmeister, “Using multi-uav for rescue environment mapping: Task planning optimization approach,” in Portuguese Conference on Automatic Control, pp. 507–517, Springer, 2020.
  * [9] X. Zhang, S. Liu, and Z. Xiang, “Optimal inspection path planning of substation robot in the complex substation environment,” in 2019 Chinese Automation Congress (CAC), pp. 5064–5068, IEEE, 2019.
  * [10] Y.-C. Huang and H.-Y. Lin, “Development and implementation of a multi-robot system for collaborative exploration and complete coverage,” in 2018 14th International Conference on Signal-Image Technology & Internet-Based Systems (SITIS), pp. 472–479, IEEE, 2018.
  * [11] V. G. Nair and K. Guruprasad, “Gm-vpc: An algorithm for multi-robot coverage of known spaces using generalized voronoi partition,” Robotica, vol. 38, no. 5, pp. 845–860, 2020.
  * [12] I. Caliskanelli, B. Broecker, and K. Tuyls, “Multi-robot coverage: A bee pheromone signalling approach,” in Artificial Life and Intelligent Agents Symposium, pp. 124–140, Springer, 2014.
  * [13] F. De Rango, N. Palmieri, X. S. Yang, and S. Marano, “Bio-inspired exploring and recruiting tasks in a team of distributed robots over mined regions,” in 2015 International Symposium on Performance Evaluation of Computer and Telecommunication Systems (SPECTS), pp. 1–8, IEEE, 2015.
  * [14] N. Palmieri, F. De Rango, X. S. Yang, and S. Marano, “Multi-robot cooperative tasks using combined nature-inspired techniques,” in 2015 7th International Joint Conference on Computational Intelligence (IJCCI), vol. 1, pp. 74–82, IEEE, 2015.
  * [15] K. Albina and S. G. Lee, “Hybrid stochastic exploration using grey wolf optimizer and coordinated multi-robot exploration algorithms,” IEEE Access, vol. 7, pp. 14246–14255, 2019.
  * [16] N. Palmieri, X.-S. Yang, F. De Rango, and S. Marano, “Comparison of bio-inspired algorithms applied to the coordination of mobile robots considering the energy consumption,” Neural Computing and Applications, vol. 31, no. 1, pp. 263–286, 2019.
  * [17] T. Stollenwerk, S. Hadfield, and Z. Wang, “Toward quantum gate-model heuristics for real-world planning problems,” IEEE Transactions on Quantum Engineering, vol. 1, pp. 1–16, 2020.
  * [18] S. Hadfield, Z. Wang, B. O’Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas, “From the quantum approximate optimization algorithm to a quantum alternating operator ansatz,” Algorithms, vol. 12, no. 2, p. 34, 2019.
  * [19] E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization algorithm,” arXiv preprint arXiv:1411.4028, 2014.
  * [20] S. Hadfield, “On the representation of boolean and real functions as hamiltonians for quantum computing,” ACM Transactions on Quantum Computing, vol. 2, no. 4, pp. 1–21, 2021.
  * [21] Z. Wang, N. C. Rubin, J. M. Dominy, and E. G. Rieffel, “X y mixers: Analytical and numerical results for the quantum alternating operator ansatz,” Physical Review A, vol. 101, no. 1, p. 012320, 2020.
  * [22] A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor, T. Sleator, J. A. Smolin, and H. Weinfurter, “Elementary gates for quantum computation,” Physical review A, vol. 52, no. 5, p. 3457, 1995.