Multistart Algorithm for Identifying All Optima of Nonconvex Stochastic Functions

We propose a multistart algorithm to identify all local minima of a constrained, nonconvex stochastic optimization problem. The algorithm uniformly samples points in the domain and then starts a local stochastic optimization run from any point that is the"probabilistically best"point in its neighborhood. Under certain conditions, our algorithm is shown to asymptotically identify all local optima with high probability; this holds even though our algorithm is shown to almost surely start only finitely many local stochastic optimization runs. We demonstrate the performance of an implementation of our algorithm on nonconvex stochastic optimization problems, including identifying optimal variational parameters for the quantum approximate optimization algorithm.


Introduction
We consider the problem of identifying all local minima of the constrained stochastic optimization problem where F (x, ξ) is the observable random function, ξ is a random variable defined on the probability space (Ω, F, P ξ ), and D ⊂ R d is compact. Identifying all local minima is relevant in various engineering and scientific applications. For instance, Gheribi et al. [13] describe a multicomponent chemical thermodynamics system where all low-melting compositions-each corresponding to a local minimum of their objective-are desired. While the system in [13] assumes that the evaluation of temperatures is deterministic, in practice, such temperatures measurements would be are accompanied by stochastic noise. Furthermore, in biophysics and biochemistry, an important problem is to understand the transitions in the shape of a protein molecule as it folds itself from a disordered (high-energy) state to a native (least-energy) state. The energy function has multiple local minima, and each local minimum corresponds to an intermediate stable state of a protein molecule. Identifying each of them can be crucial for understanding the protein folding pathways [1,8]. Li et al. [22] highlight the stochastic nature of the protein folding process. Moreover, as a direct consequence of identifying all local minima, we obtain global minimizers of f (x) as well. Identifying such minima is important in various modern applications, such as parameter tuning of stochastic event simulators [4,17,21,24] and identifying optimal parameters within a quantum approximate optimization algorithm (QAOA) [6,7].
In this paper, we propose a multistart algorithm for nonconvex stochastic optimization (MANSO) to identify all the local minima of f (x) in D using state-of-the-art local stochastic optimization methods. (Definition 2 formally defines the meaning of "identifying" a local minimum.) It is designed to extend the popular multistart algorithm Multi-Level Single Linkage (MLSL) [31,32], which guarantees to find all the local minima of a nonconvex deterministic function. MANSO uniformly samples points in D and starts a local optimization method, such as stochastic gradient descent, from points that are considered the best points in their neighborhood. By design, MANSO seeks to use as few local optimization runs as possible by carefully shortlisting sampled points based on the noisy observations of f (x). In particular, like MLSL, we construct a specific sequence of radii {r k }, where k is the number of iterations of MANSO, which enables us to certify that a uniformly sampled point is indeed probabilistically best to start a local stochastic optimization run.
We also provide theoretical guarantees on the performance of MANSO. Primarily, under assumptions that the stationary points of the true objective function f in D are sufficiently separated and the local stochastic optimization method is guaranteed to converge to the first-order stationary points with high probability, we prove that MANSO identifies all local minima of f on D with high probability. Furthermore, we show that MANSO does so while starting only finitely many local stochastic optimization runs.
In addition, we demonstrate the performance of MANSO on two benchmark nonconvex stochastic optimization problems (up to 10 dimensions) using a derivative-free adaptive sampling trust-region stochastic optimization algorithm (ASTRO-DF) [33] as a local stochastic optimization method.
Related Work In general, a nonconvex optimization problem is NP-hard even in the deterministic case, where the true function f (x) can be evaluated at any point x ∈ D. In the past few decades, extensive research has been conducted to develop stochastic algorithms for solving nonconvex deterministic optimization problems, with the aim of finding global optima. Such stochastic methods can be broadly categorized into multistart [31,32,20] and Bayesian optimization(BO) [19,10] methods. Also, some recent works leverage BO techniques for exploring the domain in the multistart framework [25]. In contrast, the parallel work for stochastic nonconvex optimization still requires much attention.
In the deterministic case, one of the popular approaches is to use a stochastic search technique where randomly sampled points in a compact domain are explored with the help of a local search method and a specific set of rules. These rules are designed to qualify any new randomly sampled point to start a local search, and they depend on the already observed values of the points in the neighborhood of the point being tested. With the help of such rules, these types of algorithms try to avoid finding the same local minima multiple times. Such types of algorithms are known as multistart algorithms [31,32,23]. Among various multistart algorithms, one of the most popular is MLSL [31,32], which provides asymptotic guarantees in identifying all the local minima assuming all stationary points of the objective function are separated by a positive distance. Locatelli [23] improves upon this assumption of separated stationary points and introduces N-MLSL (non-monotonic MLSL). Parallel implementation of MLSL is also discussed in [20] with asymptotic guarantees to find all local minima when the true function f (x) can be concurrently evaluated at any x ∈ D, unlike MLSL. Mathesen et al. [25] propose a novel approach, where the restart points are decided based on ideas developed in BO literature [10,27,35]. Instead of randomly selecting the restart points from the search domain, their method chooses points based on a surrogate Gaussian process model of the objective function. Naturally, their surrogate model is updated sequentially using the observed function values at points that are already evaluated. A detailed review of such algorithms, including significant recent developments, can be found in [10]. Similarly, the methods in [18,29,30] consider other approaches for utilizing a surrogate model within a multistart framework. Multistart methods that seek improved efficiency have considered early termination of local searches [38] and tunneling and evolutionary strategies [36]. We note that these methods are for deterministic objectives. Extending such approaches to the stochastic setting may yield similar improvements.
Our algorithm extends the MLSL algorithm to the case when we have access only to noisy function values. Like MLSL, we also rely on state-of-the-art local stochastic optimization techniques with first-and second-order convergence guarantees [16,12,11,33]. However, to the best of our knowledge, no work has developed a multistart algorithm for stochastic function evaluations.
Here is a brief roadmap of the paper. In the next section, we define notations and definitions used in the paper, followed in Section 3 by the assumptions required to prove the theoretical properties of MANSO. In Section 4 we describe the details of our proposed method MANSO, and in Section 5 we establish its theoretical properties. In Section 6.5 we describe the parameters of the quantum approximate optimization algorithm, which we will use for testing the MANSO algorithm. In Section 6.1 we present our numerical results. We conclude in Section 7 with a summary and brief description of further work.

Notations and Definitions
We use capital calligraphic letters to denote Lebesgue measurable subsets of R d . The volume (Lebesgue measure) of a set A is denoted as m(A). We use |A| to denote the cardinality of a discrete set A. Unless otherwise stated, · denotes the Euclidean norm. We denote the set of natural numbers as N.
We use Var ξ [·] and Cov ξ [·, ·] to represent the variance and covariance, respectively, with respect to the random variable ξ. If {A n (ξ)} n∈N is a sequence of events and P ξ ({A n (ξ)}) is the probability of event n occurring, then this sequence of events occurs with high probability (w.h.p.) if for any δ > 0 there exists an n 0 ∈ N such that P ξ (A n (ξ)) > 1 − δ for any n ≥ n 0 . We index the probability measure of the random variable ξ by ξ itself just to differentiate it from the uniform sampling measure used subsequently in the paper. Moreover, the samples space of the random variable ξ, Ω is arbitrary and the range of ξ is an arbitrary Borel measurable set. When a sequence {a t } is O(b t ), it implies that there exist a K > 0 and t 0 ≥ 1 such that ∀t ≥ t 0 , a t ≤ Kb t . We letf n (·) denote an estimate of f (x) constructed using n i.i.d. measurements of F (x, ξ).
We now establish notation to be used to describe the MANSO algorithm. Let k ∈ N denote the number of iterations of MANSO. Let S k ⊆ D be the collection of uniformly sampled points up to and including iteration k. Let A k ⊆ S k be the collection of sampled points from which a local stochastic optimization method has been started and is still active at iteration k, and let L k denote the set of points generated from all the local stochastic optimization runs up to and including iteration k. Let X * ⊂ D denote the set of local minima of f in D andX * k denote the set of local minima identified before iteration k of MANSO. Let Y * denote the set of stationary points of f on D that are not local minima. We let {X a i , i ∈ N} denote the random sequence of iterates generated by a given local stochastic optimization method started from a ∈ D. {X a i , i ∈ N} is a stochastic process defined on the probability space (Ω, F a , P ξ ). Let {x i , i ∈ N} be its realization. We represent the filtration (information) available at iteration i of the local stochastic optimization method started from a ∈ D as {F a i }, which is an increasing family of σ-algebras of {F a } on which the stochastic process {X a i } is defined. We also call {F a i } a filtration at iteration i of the local stochastic optimization (LSO) run started from a ∈ D.
We now list other important notation: for some fixed σ > 4, where Γ is the gamma function. The radius is used in Lemma 2 to define the neighborhood of a candidate "probabilistically best" point.
• Let ∂D denote the boundary of the set D.
• Let B (y; r) := {x ∈ D : x − y ≤ r} represent a ball of radius r centered at any y ∈ D.
• For τ > 0, let be the points in D within τ of the boundary.
• We denote η as the minimum distance between any two distinct stationary points of f on D, that is, Next, we define the domain of attraction for method M for any local minima.
Definition 1 (Domain of attraction). For any x * ∈ X * , the domain of attraction L x * is defined as a subset of D such that if a local stochastic method M is started from any point in it, then it will converge to the local minimum x * w.h.p. Formally, for any > 0, Note that the domain of attraction is defined only for those points that the local stochastic method M will converge to with high probability. Moreover, it is not required that the respective domain of attraction for each x * ∈ X * partition D. We consider Definition 1 to be a reasonable stochastic extension of the domain of attraction for the deterministic case considered in the original MLSL paper and its extensions (e.g., [31,Theorem 4], [23], [20,Assumption 2]).
We also define the event ω−identifying a local minimum, which we use throughout. In our algorithm, ω is a tuning parameter and is given as input by the user.
Definition 2 (ω-identifying local minima). Let ω ≤ η 2 be fixed, and let {X a i } be the sequence of iterates generated by the local method M started from the point a ∈ L x * . A local method M has identified the local minimum x * when the event { X a i − x * < ω} occurs for all i ≥ i 0 .

Assumptions
We now state the assumptions needed for our theoretical analysis. We group them into assumptions about the objective and domain of the problem (SO), assumptions about the true function f and its estimate computed by using F (x, ξ) at x ∈ D, and assumptions about the local stochastic optimization method used within the multistart framework.
Assumption 1. We first impose the following conditions on the problem (SO).
1. D is a compact set, and f is twice continuously differentiable on D.
2. The minimum distance between any two distinct stationary points is positive, that is, η > 0.

There exists
Because the set D is compact by Assumption 1.1, Assumption 1.2 implies that f has finitely many stationary points in D and that f is not flat in D. Assumption 1.3 ensures that there are no stationary points near the boundary of D. (The parts of Assumption 1 are the same as those considered in the deterministic case [31].) It is useful to have notation for the set of points in D within ω of a stationary point of f : For any ω ∈ (0, η), let Assumption 2. For any x ∈ D, we assume that the estimatef n (x) of f (x), which is constructed by using n i.i.d. measurements of the measurable function F (x, ξ), satisfies Sincef n (·) is constructed by using n i.i.d. measurements of the measurable function F (·, ξ), this assumption requires F (·, ξ) to satisfy some regularity conditions. For instance,f n (·) = 1 Assumption 3. We make the following assumptions about the local stochastic optimization method M . If {X a i } is a random sequence of iterates produced by M when started from a ∈ L x * , then it satisfies the following.
where n is the number of random samples of f (·) at respective points.
, for a given ω < η 2 , whereX * k−1 is the collection of approximate local minima up to iteration (k − 1). S3 a / ∈ Q τ , that is, near the boundary of set D. S4 a has not started any LSO.
4. For x, y ∈ X * and x = y, we also assume that the sequence of iterates {X a i } and {Y b j } generated by LSOs started at a ∈ L x and b ∈ L y , respectively. Then Some of the conditions on the local method M in Assumption 3 are strong conditions that may be difficult to satisfy in practice by a local stochastic optimization method on a nonconvex problem. Yet, we find these assumptions to be a natural stochastic version of their deterministic counterparts: for Assumption 3.3 above is similar to the strictly decent property assumed for the (deterministic) local optimization method [31]. Assumption 3.1 ensures that the domain of attraction for any local minima has positive measure and that L x * does not lie entirely in the boundary set Q τ . Assumption 3.2 ensures that no two local minima have overlapping domains of attraction. In general, local stochastic methods [12,11] guarantee convergence to a stationary point only w.h.p.; that is, they satisfy Assumption 3.3(i). In addition, we need Assumption 3.3.1(ii) to ensure that x * is within an ω-ball of all iterates after a large number of iterations with probability 1, given that the local method is started in that L x * and the norm of the gradient at the last iterate (of an LSO run ) is small enough (less than ν). Furthermore, in Assumption 3.4, we assume that the iterates (realizations of the sequence of iterates) generated by any two LSO runs started in different domains of attractions are at least ω apart. These conditions are necessary for developing the algorithm and for showing that the algorithm MANSO identifies all the local minima w.h.p. using the LSO method M . Table 1 lists the conditions that our algorithm checks when deciding where to start an LSO. Algorithm 1 states our algorithm for ω-identifying all local minima of f (x).

Statement of the Algorithm
Uniformly sample a point a in D, evaluatef n (a) and add it to S k . 5 Start LSO method M from all points in S k satisfying the conditions listed in Table 1 for a given r k , β, σ, ω, and sampling effort n, ∀x ∈ S k . Add those points to A k . 6 Update L k by adding the next iterate generated by each LSO started from all the points in A k . 7 Terminate any LSO run if its current iterate is within 2ω distance of any point in L k−1 (from other LSOs) and remove that LSO run from A k . 8 UpdateX * k by adding any new local minima identified during the local searches, and remove that LSO run from A k .
Notice that in condition S1 of Table 1, we have used n ∈ N as the number of random samples of f (·) at respective points. We show in Theorem 1 later that for any n ∈ N the total number of LSO runs started by Algorithm 1 is finite. Since the results on the diffusion approximation of nonconvex SGD [14] show that the use of smaller batch sizes in batch SGD methods help escape sharp local minima and nondegenerate saddle points, we anticipate that the use of a larger n would be "better and faster" in ω-identifying all the local minima.
Note that MANSO can be viewed as a stochastic analogue of MLSL. That is, if we further assume that (1) (2) Var ξ f n (x) → 0 as n → ∞ , and for any two points {x, y} ∈ D, Cov ξ f n (x),f n (y) → 0 as n → ∞, then the above notion (see Table 1, S1) of not finding a "probabilistically best" point z in its r k −neighborhood converges to the original MLSL [31,32] condition of finding a "better" point z, as n → ∞. Also, while implementing MANSO we estimate Var ξ f n (x) −f n (y) using samples of f (·) for any {x, y} ∈ D. In step 7 of Algorithm 1, since we do not have prior knowledge about the minimum separation between stationary points, η, we must choose a small enough positive value for each ω and τ .

MANSO Asymptotic Analysis
MANSO seeks to use as few LSO runs as possible to ω-identify all local minima of f in D. We will show that the total number of LSO runs started by MANSO is finite, even if the algorithm is run forever. We state our main theoretical results in this section, but to ease presentation, the proofs of lemmas are deferred to the appendix.
We first show a limiting result for the measure of balls around any point a ∈ D \ (Q τ ∪ T ω ), where Q τ and T ω are defined in (1) and (3), respectively. Lemma 1. Under Assumption 1.1 and 2, for any a ∈ D \ (Q τ ∪ T ω ), β ∈ (0, 1/2), and for all n ≥ 1 Specifically, Lemma 1 derives a bound on the measure of the set of points that are not probabilistically best in a ball of radius r around a point a drawn uniformly from the set of points in D not within ω of a stationary point or τ of the boundary of D. (Not being probabilistically best is measured with respect to β, an integer n, and the tolerance n (x; a).) Lemma 1 shows in the limit as r converges to zero that the set of not probabilistically best points has a measure of at least half of the ball around a. Note that it is true for any n ∈ N.
Next, we use Lemma 1 and construct a specific sequence of radius {r k } to show that the probability of starting an LSO run from any previously sampled point is bounded by a term that converges to zero as the number of iterations increases. The proof of the following result uses arguments similar to those used in [31,Theorem 8].
Lemma 2. Let t k be the number of LSO runs started from any point in S k \ (Q τ ∪ T ω ) during iteration k of Algorithm 1. Then under Assumption 1.1 and 2 and for r k as defined in Section 2 and used in Table 1 with Subsequently, we use the summability of the bound obtained in Lemma 2 on the probability of starting any LSO run from the set of sampled points not within ω of any stationary point or τ of the boundary of D to show in Theorem 1 that the total number of LSO runs started by Algorithm 1 is finite. Furthermore, condition S3 in Table 1 ensures that LSO does not start from any point in Q τ , and step 7 in Algorithm 1 ensures that the total number of LSO runs started from points sampled within the ω-ball of any stationary point is finite as k → ∞. (The total number of stationary points is finite because of Assumption 1.

2.)
That is, the number of LSO runs started by MANSO is finite even if MANSO runs forever.
where t k is the number of LSO runs started from points in a ∈ S k \ (Q τ ∪ T ω ) during iteration k of Algorithm 1. Note that using the first Borel-Cantelli lemma [5, Theorem 2.3.1], we have from equation (4) that This is equivalent to Since t k ≥ 0, ∀k ≥ 1, the result in (5) implies that lim k→∞ t k → 0 P − almost surely, that is, with P-probability 1. Since {t k } is a sequence of natural numbers, it implies that Furthermore, if a point in S k belongs to Q τ , then MANSO (see S3 in Table 1) does not start a run at that point, since we assumed in Assumption 1.3 that no minimum lies in Q τ . Therefore, the probability of starting an LSO run from any point a ∈ S k ∩Q τ is zero, and the result in (6) holds true for all a ∈ (S k \(Q τ ∪T ω ))∪Q τ , that is, for all a ∈ S k \ T ω . Now consider the last case when a ∈ S k ∩ T ω . Let t k be the number of LSO runs started by MANSO in the iteration k from any point a ∈ S k ∩ T ω . Since in step 7 of Algorithm 1 we kill the LSO run if its current iterate is within 2ω distance of an already generated iterate from any of the previous LSO, a run can be started from any a ∈ S k ∩ T ω at most once. Since Assumption 1.2 implies that there are only a finite number of local minima, there exists a k 0 ≥ 1 such that for all k ≥ k 0 the number of LSO runs started at points a ∈ S k ∩ T ω will be zero. Therefore the number of LSO runs t k → 0 as the number of samples increases to infinity. Since {t k } is a sequence of natural numbers, it implies that Since t k = t k + t k , the result follows immediately by adding (6) and (7).
Our next goal is to show that under certain assumptions MANSO will ω-identify all the local minima w.h.p. Recall that X * is the collection of local minima of f in D. The Assumption 3.3 guarantees that if MANSO starts an LSO run from any point in the domain of attraction L x * of a local minimum x * ∈ X * , then the LSO method M identifies it w.h.p. In particular, we prove that if a point is sampled in the domain of attraction of a local minimum, then the probability of that local minimum not being identified is sufficiently small for large enough iterations of both MANSO and LSO. Combined with the fact that these domains of attraction are of positive measure and the probability of getting a uniformly sampled point in any domain of attraction approaches 1 as the number of iteration increases (since the number of sampled points increase with each iteration), we show that all the local minima are identified w.h.p.
Theorem 2. Under Assumption 3, and given the sequence of radii {r k } as constructed in Lemma 2 and used in condition S1 of Table 1, MANSO identifies all the local minima of (SO) w.h.p.
Proof. Let a be a point sampled in iteration k of Algorithm 1 in the domain of attraction L x * of the local minimum x * (see Definition 1) for the first time and none of the points sampled before belong to L x * . Since m(L x * ) > 0 because of Assumption 3.1, the probability of obtaining at least a uniformly sampled point in L x * approaches 1 as the number of samples increases to infinity [2]. We also assume that the local minimum x * has not been identified yet.
RecallX * k is the collection of approximate local minima identified up to iteration k and that set A k is the collection of sampled points from which the LSO run has started and is still active up to iteration k. Recall the definition (see Definition 2) of the event that the local minimum is identified at any arbitrary iteration l as I l =: {∀ l > l ∃p ∈ A l : X p l − x * < ω}, where {x p l , l ∈ N} is the sequence of iterates generated by LSO at p. Now, let us compute the probability that the local minimum has been identified up to the iterationk ≥ k. Observe that = 1 − P ξ (LSO has started from a in iteration k, I C k |Sk ∩ L x * = ∅) + P ξ ( LSO has not started from a in iteration k , I C k |Sk ∩ L x * = ∅) .
We first analyze the first probabilistic term in (8). Since LSO starts at a, then by definition of A k , a ∈ A k . Even if we sample another point in L x * at any iterationk > k, we will have a ∈ Ak because we terminate only the latest LSO in step 7 of MANSO. Also, if at some iterationk > k any of the iterates generated by some other sampled point, b (at iteration k < k) in L y * , for x * = y * , jumps into L x * , then because of Assumption 3.2 and Assumption 3.4 we will still have a ∈ Ak on the event I C k , since LSO at a will not be terminated at step 7 of MANSO because Consequently, since a ∈ L x * and it remains in Ak for anyk > k, then because of Assumption 3.3 there exists a k 0 ∈ N such that for any iterationk − k > k 0 of LSO at a, we have P ξ (LSO has started from a in iteration k, I C k |Sk ∩ L x * = ∅) < Note that we use the same iteration counter for the LSO run and MANSO, since we progress one step of each active LSO in Ak in each iteration of MANSO (see Step 6). Choosingk large enough such thatk − k > k 0 , we obtain the last inequality in (9) by using Assumption 3.3, since LSO at a never gets terminated given that the local minimum x * has not been ω-identified.
Next, we analyze the second probabilistic term in (8). In this case LSO may not start from a in iteration k, since it gets rejected because of any of the following conditions from Table 1 not being satisfied by a in the iteration k. Using these conditions in the third term of (8), we have P ξ ( LSO has not started from a in iteration k, where S(1) denotes the event that a does not satisfy condition S1 from Table 1 and similarly for S2, S3, and S4. The last three cases are straightforward to analyze. First consider event S2. Since we assumed that the local minimum x * has not been identified yet, then the event {a ∈ ∪ x∈X * k−1 B (x; ω)} is of probability measure zero, and hence For event S3, a ∈ Q τ , and we assumed at the beginning that a ∈ L x * as well. But because of Assumption 1.3 L x * ∩ Q τ = ∅. Therefore it follows that Since we assumed at the beginning that for anyk < k none of the uniformly sampled points belong to L x * , the final event S4 is an impossible event, and therefore For the first event S1 we assumed that there existsā ∈ B (a; r k ) ∩ (S k ) that does not satisfies S1 for a point a in L x * . Now notice that First consider the case whenā / ∈ L x * . Since r k → 0 as k → ∞ and a ∈ L x * , there must exist a k 0 ≥ k such thatā / ∈ B (a; rk) ∩ (Sk) for allk ≥ k 0 . Therefore, LSO will start from a at iteration k 0 . Using arguments similar to those used in (9), there exists a k ≥ k such that ∀k ≥ k On the other hand, because of our assumption that for anyk < k none of the uniformly sampled points belong to L x * , the caseā ∈ L x * is an impossible event, and thus Since m(L x * ) > 0 due to Assumption 3.1, the probability of obtaining at least a uniformly sampled point in L x * approaches 1 as the number of samples increases to infinity [2]. Therefore limk →∞ P ξ (Sk ∩ L x * = ∅) = 1. Now, substituting equation (9)-(15) into (8), we obtain, for large enoughk, whereΛk = min{Λk −k0 , Λk −k }. Therefore, we have shown that any local minimum x * ∈ X * will be identified w.h.p. Also, note that because of Assumption 3.3,Λk decreases to 0 as the number of iterationsk of LSO increases with the number of iterations of MANSO.

Numerical Experiments
We compare implementations of MANSO in their ability to solve difficult synthetic benchmark problems and to identify optimal variational parameters within the quantum approximate optimization algorithm (QAOA) [6]. Our MANSO implementation and scripts to perform our numerical experiments are available: https://github.com/prat212/MANSO.git

Synthetic Benchmark Experiments
We benchmark our implementation of MANSO on nonconvex optimization problems with large variance in their observations. In particular we fix two non-convex benchmark functions, Branin-Hoo (d = 2) and Shekel (d = 4, 6, 8 and 10), and make each non-convex objective evaluation stochastic by adding a Gaussian noise with variance 1. We add this significant noise only to make the testing of MANSO rigorous and robust. We generate 10 sample paths of Gaussian noise to create a set of 10 problems each for the Branin and Shekel functions. We seek to find all of the local minima for each problem within a fixed budget of function evaluations, B. We use ASTRO-DF [33] as the local method; it is a derivative-free trust-region stochastic optimizer selected because of its theoretical guarantee to converge to first-order critical points of the objective function. Other optimizers with convergence guarantees to first-order critical points (e.g., [11,12]) could naturally be used as the local method in MANSO. To improve performance in our numerical experiments, we ensure that there are at most 10 active LSO runs at any given iteration. That is at each iteration k of MANSO, we sample a point uniformly in D if the total number of active runs is no larger than a fixed threshold, heuristically set to 10 in the experiments. Naturally, MANSO ensures that the conditions listed in Table 1 are satisfied before starting an LSO run from all the sampled points. Note that for condition (S1) in Table 1, for any two points {x, a} ∈ D, our implementation estimates Var ξ f n (x) −f n (a)) using the n samples of f (x) and f (a).
We compare the performance of MANSO with drawing points uniformly in the search domain. We compare only with such a random search method because we are unaware of other methods that aim to find all local optima of stochastic nonconvex functions. Our empirical results demonstrate that MANSO outperforms random search in identifying points that are within a small ball of all local minima. As we would like random search to perform equally well across problems independent of problem dimension, we consider each ball around a local minimum to always have a volume that is a small fraction (e.g., 1 /1000) of the volume of the domain. This number is arbitrarily chosen but gives a sense that how fast MANSO and random search can evaluate points near the local minima. We measure the performance of methods using data profiles [26].

Data profiles
Data profiles present the fraction of problems "solved" from a set of problems P after a certain number of function evaluations by an implementation of method h in a set of methods H. The set of implementations H is created by adjusting β, ω, and n of MANSO. For a given objective function, we create different problems by changing the initial random seed. For our comparisons, the set H contains a uniform random sampling method and different versions of MANSO, obtained by varying its tuning parameters ω, τ, β, and n. Data profiles mark a problem instance p ∈ P as solved based on a user-defined test criterion. We use the criterion proposed in [20] to classify that a problem p ∈ P is solved: a problem is solved when a point is evaluated near each local minima for the problem. Let there be j local minima for a given problem. Mathematically, we define a test that ensures that a local minima x * ∈ {x * 1 , x * 2 , . . . , x * j } is identified at level ρ d (ζ) after e evaluations, as where ρ d (ζ) = 1 √ π d Γ(1 + d/2)m(D)ζ and the set E e is constructed by sequentially adding points the number of times they are being evaluated by a method h ∈ H. That is, all the points in S k ∪ L k till e out of B budget is used. Notice that the volume of a ball of radius ρ d (ζ) is ζ times the volume of the search domain D. Hereafter, we use (17) for a problem p ∈ P, method h ∈ H and for all and define the data profile metric for e > 0 function evaluations as In the next section, we present the details of two benchmark nonconvex functions on which we evaluated the performance of MANSO.

Benchmark problems
We consider the Branin-Hoo and Shekel (4, 6, 8, and 10 dimensions) functions. Results for d = 6 and d = 8 Shekel problems appear in Appendix B.

Shekel function:
The Shekel function is the d-dimensional nonconvex problem with m local minima:

Experimental analysis
We plot the data profiles for the Branin function (d = 2) in Figure 1 and the Shekel d = 4 and d = 10 functions in Figures 2 and 3, respectively. We observe from the data profiles presented in each plot that MANSO outperforms the uniform random search method in finding a point in a ball of volume 10 −4 (10 −3 for the Branin function) times the volume of the domain D centered at the respective true local minima. Next, we discuss the effect of MANSO hyperparameters on its performance. Effect of n: Recall that as n increases the variance in estimates of the function values decrease. However, this confidence is attained at the expense of shedding more budget. Consequently, the number of points evaluated by MANSO decrease. Hence varying n controls the trade-off between exploration and variance in function evaluation. We can observe this effect by comparing the MANSO performance on Shekel-4D function in Figure 2. Nonetheless, it is evident from the plots in Figure 2 that the increasing n does not guarantee that finding the local minima will be faster. Effect of ω: Recall that ω is a hyperparameter used in Step 7 of MANSO, to terminate an LSO run if any of its iterates are within 2ω distance of any iterate generated by some other LSO. Intuitively, a larger ω would result in more termination. In Figure 1, the effect of ω can be observed by comparing the green (ω = 0.05) and blue (ω = 0.5) lines. Note that MANSO under the green experiment was able to explore more points and thus identified 2 local minima of the Branin function faster than the blue experiment. In the blue experiments, due to larger ω, the number of points evaluated by MANSO is less than the green as a large portion of the budget is used for evaluating new sampled points.

Variational Parameter Optimization
Quantum approximate optimization algorithm (QAOA) is a hybrid algorithm that uses a parameterized trial quantum state ψ(x) as defined by the parameters x. (The values in x ∈ R 2p are rotations or angles that parameterize 2p unitary operators.) What is desired is parameters x such that when the trial state is measured, the measurement outcome corresponds to the solution of the optimization problem. This is achieved by finding parameters x that give a large expected value for ψ(x) T Hψ(x), where H is the problem Hamiltonian encoding some classical objective h. Under certain conditions on the Hamiltonian H, f (x) must be evaluated by using a quantum computer. The search for optimal parameters x can therefore be considered as a (classical) numerical optimization problem of the form (SO) with f (x) = ψ(x) T Hψ(x). The  stochasticity in the objective arises from not being able to compute the value of observable ψ(x) T Hψ(x) by using a quantum circuit but rather having to compute the objective from a sample: f (x) = yi∈sample h(y i ).
QAOA has nontrivial performance guarantees [6,7] and requires the execution of only moderately sized quantum circuits, with the depth controlled by the number of steps p. For these reasons, QAOA is an especially promising candidate algorithm for demonstrating quantum advantage on near-term quantum computers. Yet, the quality of the solution produced by QAOA depends critically on the quality of the parameters x used by the algorithm. Identifying such parameters is difficult because the objective landscape is highly nonconvex with many local minima with poor objective values [37,34]. Figure 4 shows an example contour plot with p = 1. While nonglobal optima are not necessarily of interest in the QAOA problem setting, we consider the difficulty of finding a global optimum to be a considerable test of our MANSO implementation.
We consider the problem of using QAOA to find the maximum cut on the Petersen graph with a depth of p = 5, that is, d = 10. The global optimal value for this problem is −12. The performance of MANSO to identify maximum cut with three local solvers (ASTRO-DF [33], BOBYQA [3], and Snobfit [15]) are summarized in Figure 4. Moreover, we also considered non-MANSO global optimizers such as Bayesian optimization [28] to solve a deterministic version of the QAOA problem. However, the method was significantly slow due to large matrix computation and produced the best candidate global minima with a value -11.74 only after 5000 evaluations. Consequently, we are not comparing MANSO with other approaches such as Bayesian optimization. For each local solver, we run MANSO on the MAXCUT problem with a budget of 150,000 function evaluations and we check for the termination condition in step 7 of MANSO after 500 function evaluations have been performed by the local search method. We also fix n = 5, ω = 0.01, β = 0.1 and τ = 0.01. We repeat each experiment with a given local solver 20 times and plot the range of best function value identified in Figure 4. Although BOBYQA is designed to solve deterministic problems, we use it for this stochastic problem as it has been reported that it empirically performs well on problems with stochastic noise [3]. In particular, it is evident from our experimental result in Figure 4 (right) too that BOBYQA performance is competitive with other stochastic solvers. However, we note that BOBYQA failed on 10 out of 20 experiments as it produced singular Hessian matrices of the noisy QAOA objective.

Conclusion
We propose the MANSO algorithm to identify all the local minima of a stochastic nonconvex function. We construct an efficient scheme to judiciously determine when to start a local stochastic optimization run from a sampled point in a compact search domain. We show that under that MANSO starts only finitely many local stochastic optimization runs. We also show that MANSO identifies all the local minima asymptotically with high probability, given that the local stochastic optimization method is guaranteed to converge to a local minimum with high probability. MANSO's theoretical guarantees also require that the sequence of iterates generated from the local stochastic search started in a domain of attraction and cannot leave that domain with high probability. (Certainly, this is a restrictive assumption for a stochastic optimization method, but it is analogous to the assumptions in the foundational MLSL work [31,32].) Our experimental results show that MANSO can display strong performance even when coupled with a local optimization that does not satisfy such a restrictive assumption. The assumption that there are no flat regions in the true objective function may be removed by using the techniques developed in [23] for a multistart algorithm MLSL for deterministic nonconvex objectives.
Furthermore, we demonstrate the efficacy of our algorithm on two benchmark problems with dimensions ranging from 2 to 10, and we compare the performance with that of a uniform random search method. We also use MANSO to find the global minima of a highly nonconvex 10-dimensional Peterson graph. We aim to apply MANSO to more complex and higher-dimension benchmark functions and application problems as part of our future work. Similar to [20], an asynchronously parallel version of MANSO can be developed to improve its computational performance.

A Proofs
Below are the proofs of Lemma 1 and Lemma 2.
Recall that T ω is the union of ω-radius balls centered at stationary points of f in D. Now consider a ∈ D \ (Q τ ∪ T ω )}. Define the set E(a; r; ρ) := x ∈ D : x − a ≤ r and ∇f (a) T (x − a) + 1 2 where ρ is the largest eigenvalue of ∇ 2 f (x) for x ∈ D. Using the Taylor expansion of f around a, we know that for all x ∈ D, with x − a ≤ r, there exists a θ ∈ [0, 1] such that f (x) − f (a) = ∇f (a) T (x − a) + 1 2 (x − a) T ∇ 2 f (a + θ(x − a))(x − a).
For ease of reference, let H = ∇ 2 f (a + θ(x − a)) = H and v = x − a.
Proof of Lemma 2. For any a ∈ S k \ (Q τ ∪ T ω ) and using the fact that S k is sampled uniformly, for a vanishing sequence {r k } the Lemma 1 implies that there exists a k 0 ≥ 1 such that for all k ≥ k 0 , where the first inequality bounds the probability that for any sampled point in S k \ (Q τ ∪ T ω ) none of the remaining sampled points are in A (a; r k ; n; β) (see (S1) of Table 1). The second inequality follows from Boole's inequality. The last inequality in (31) is due to Lemma 1 and the fact that |S k \ (Q τ ∪ T ω )| < |S k |.
Recall that for any a ∈ D and r k sufficiently small such that B (a; r k ) ⊆ D, m(B (a; r k )) = r d k π d/2 Γ(1+d/2) . Combined with this fact, for any σ > 0 and choosing a sequence r k = 1 √ π d Γ(1 + d/2)m(D)σ log |S k | |S k | , we have Since |S k | → ∞ as k → ∞, the fact that e B Shekel with d = 6 and d = 8