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


Learning Algebraic Representation for Systematic Generalization in Abstract Reasoning

Is intelligence realized by connectionist or classicist? While connectionist approaches have achieved superhuman performance, there has been growing evidence that such task-specific superiority is particularly fragile in systematic generalization. This observation lies in the central debate between connectionist and classicist, wherein the latter continually advocates an algebraic treatment in cognitive architectures. In this work, we follow the classicist's call and propose a hybrid approach to improve systematic generalization in reasoning. Specifically, we showcase a prototype with algebraic representation for the abstract spatial-temporal reasoning task of Raven's Progressive Matrices (RPM) and present the ALgebra-Aware Neuro-Semi-Symbolic (ALANS) learner. The ALANS learner is motivated by abstract algebra and the representation theory. It consists of a neural visual perception frontend and an algebraic abstract reasoning backend: the frontend summarizes the visual information from object-based representation, while the backend transforms it into an algebraic structure and induces the hidden operator on the fly. The induced operator is later executed to predict the answer's representation, and the choice most similar to the prediction is selected as the solution. Extensive experiments show that by incorporating an algebraic treatment, the ALANS learner outperforms various pure connectionist models in domains requiring systematic generalization. We further show the generative nature of the learned algebraic representation; it can be decoded by isomorphism to generate an answer.

Abstract. Is intelligence realized by connectionist or classicist? While connectionist approaches have achieved superhuman performance, there has been growing evidence that such task-specific superiority is particularly fragile in systematic generalization. This observation lies in the central debate between connectionist and classicist, wherein the latter continually advocates an algebraic treatment in cognitive architectures. In this work, we follow the classicist's call and propose a hybrid approach to improve systematic generalization in reasoning. Specifically, we showcase a prototype with algebraic representation for the abstract spatial-temporal reasoning task of Raven's Progressive Matrices (RPM) and present the ALgebra-Aware Neuro-Semi-Symbolic (ALANS) learner. The ALANS learner is motivated by abstract algebra and the representation theory. It consists of a neural visual perception frontend and an algebraic abstract reasoning backend: the frontend summarizes the visual information from object-based representation, while the backend transforms it into an algebraic structure and induces the hidden operator on the fly. The induced operator is later executed to predict the answer's representation, and the choice most similar to the prediction is selected as the solution. Extensive experiments show that by incorporating an algebraic treatment, the ALANS learner outperforms various pure connectionist models in domains requiring systematic generalization. We further show the generative nature of the learned algebraic representation; it can be decoded by isomorphism to generate an answer.

Introduction
"Thought is in fact a kind of Algebra." -William James [25] Imagine you are given two alphabetical sequences of "c, b, a" and "d, c, b", and asked to fill in the missing element in "e, d, ?". In nearly no time will one realize the answer to be c. However, more surprising for human learning is that, effortlessly and instantaneously, we can "freely generalize" [37] the solution to any partial consecutive ordered sequences. While believed to be innate in early development for human infants [39], such systematic generalizability has constantly been missing and proven to be particularly challenging in existing connectionist models [2,29]. In fact, such an ability to entertain a given thought and semantically related contents strongly implies an abstract algebra-like treatment [12]; in literature, it is referred to as the "language of thought" [11], "physical symbol system" [44], and "algebraic mind" [37]. However, in stark contrast, existing connectionist models tend only to capture statistical correlation [7,26,29], rather than providing any account for a structural inductive bias where systematic algebra can be carried out to facilitate generalization.
This contrast instinctively raises a question-what constitutes such an algebraic inductive bias? We argue that the foundation of modeling counterpart to the algebraic treatment in early human development [37,39] lies in algebraic computations set up on mathematical axioms, a form of formalized human intuition and the beginning of modern mathematical reasoning [17,34]. Of particular importance to the algebra's basic building blocks is the Peano Axiom [46]. In the Peano Axiom, the essential components of algebra, the algebraic set, and corresponding operators over it, are governed by three statements: (1) the existence of at least one element in the field to study ("zero" element), (2) a successor function that is recursively applied to all elements and can, therefore, span the entire field, and (3) the principle of mathematical induction. Building on such a fundamental axiom, we begin to form the notion of an algebraic set and induce the operator to construct an algebraic structure. We hypothesize that such an algebraic treatment set up on fundamental axioms is essential for a model's systematic generalizability, the lack of which will only make it sub-optimal.
To demonstrate the benefits of adopting such an algebraic treatment in systematic generalization, we showcase a prototype for Raven's Progressive Matrices (RPM) [48,49], an exemplar task for abstract spatial-temporal reasoning [51,68]. In this task, an agent is given an incomplete 3ˆ3 matrix consisting of eight context panels with the last one missing, and asked to pick one answer from a set of eight choices that best completes the matrix. Human's reasoning capability of solving this abstract reasoning task has been commonly regarded as an indicator of "general intelligence" [4] and "fluid intelligence" [19,24,55,56]. In spite of the task being one that ideally requires abstraction, algebraization, induction, and generalization [4,48,49], recent endeavors unanimously propose pure connectionist models that attempt to circumvent such intrinsic cognitive requirements [21,51,58,64,68,70,74]. However, these methods' inefficiency is also evident in systematic generalization; they struggle to extrapolate to domains beyond training [51,70], shown also in this paper.
To address the issue, we introduce an ALgebra-Aware Neuro-Semi-Symbolic (ALANS) learner. At a high-level, the ALANS learner is embedded in a general neuro-symbolic architecture [14,36,66,67] but has the on-the-fly operator learnability, hence semi-symbolic. Specifically, it consists of a neural visual per-ception frontend and an algebraic abstract reasoning backend. For each RPM instance, the neural visual perception frontend first slides a window over each panel to obtain the object-based representation [26,63] for every object. A belief inference engine latter aggregates all object-based representation in each panel to produce the probabilistic belief state. The algebraic abstract reasoning backend then takes the belief states of the eight context panels, treats them as snapshots on an algebraic structure, lifts them into a matrix-based algebraic representation built on the Peano Axiom and the representation theory [23], and induces the hidden operator in the algebraic structure by solving an inner optimization problem [3,8,72]. The answer's algebraic representation is predicted by executing the induced operator: its corresponding set element is decoded by isomorphism, and the final answer is selected as the one most similar to the prediction.
The ALANS learner enjoys several benefits in abstract reasoning with an algebraic treatment: 1. Unlike previous monolithic models, the ALANS learner offers a more interpretable account of the entire abstract reasoning process: the neural visual perception frontend extracts object-based representation and produces belief states of panels by explicit probability inference, whereas the algebraic abstract reasoning backend induces the hidden operator in the algebraic structure. The final answer's representation is obtained by executing the induced operator, and the choice panel with minimum distance is selected. This process much resembles the top-down bottom-up strategy in human reasoning missed in recent literature [21,51,58,64,68,70,74]: humans reason by inducing the hidden relation, executing it to generate a feasible solution in mind, and choosing the most similar answer available [4]. 2. While keeping the semantic interpretability and end-to-end trainability in existing neuro-symbolic frameworks [14,36,66,67], ALANS is semi-symbolic in the sense that the symbolic operator can be learned and concluded on the fly without manual definition for every one of them. Such an inductive ability also enables a greater extent of the desired generalizability. 3. By decoding the predicted representation in the algebraic structure, we can also generate an answer that satisfies the hidden relation in the context. This work makes three major contributions. (1) We propose the ALANS learner, a neuro-semi-symbolic design, in contrast to existing monolithic models.
(2) To demonstrate the efficacy of incorporating an algebraic treatment in reasoning, we show the superior systematic generalization ability of the proposed ALANS learner in various extrapolatory RPM domains. (3) We present analyses into both neural visual perception and algebraic abstract reasoning.

Quest for Symbolized Manipulation
The idea to treat thinking as a mental language can be dated back to Augustine [1,62]. Since the 1970s, this school of thought has undergone a dramatic revival as the quest for symbolized manipulation in cognitive modeling, such as "language of thought" [11], "physical symbol system" [44], and "algebraic mind" [37]. In their study, connectionist's task-specific superiority and inability to generalize beyond training [7,26,51,68] have been hypothetically linked to a lack of such symbolized algebraic manipulation [7,29,38]. With evidence that an algebraic treatment adopted in early human development [39] can potentially address the issue [2,36,38], classicist [12] approaches for generalizable reasoning used in programs [40] and blocks world [61] have resurrected. As a hybrid approach to bridge connectionist and classicist, recent developments lead to neuro-symbolic architectures. In particular, the community of theorem proving has been one of the earliest to endorse the technique [13,50,53]: BILP [10] and NLM [9] make inductive programming end-to-end, and DeepProbLog [35] connects learning and reasoning. Recently, Hudson and Manning [22] propose NSM for visual question answering where a probabilistic graph is used for reasoning. Yi et al . [67] demonstrate a neuro-symbolic prototype for the same task where a perception module and a language parsing module are separately trained, with the predefined logic operators associated with language tokens chained to process the visual information. Mao et al . [36] soften the predefined operators to afford end-to-end training with only question answers. Han et al . [14] use the hybrid architecture for metaconcept learning. Yi et al . [66] and Chen et al . [6] show how neuro-symbolic models can handle explanatory, predictive, and counterfactual questions in temporal and causal reasoning. Lately, NeSS [5] exemplifies an algorithmic stack machine that can be used to improve generalization in language learning. ALANS follows the classicist's call but adopts a neuro-semi -symbolic architecture: it is end-to-end trainable as opposed to Yi et al . [66,67] and the operator can be learned and concluded on the fly without manual specification.

Abstract Visual Reasoning
Recent works by Santoro et al . [51] and Zhang et al . [68] arouse the community's interest in abstract visual reasoning; the task of Raven's Progressive Matrices (RPM) is introduced as such a measure for intelligent agents. As an intelligence quotient test for humans [48,49], RPM is believed to be strongly correlated with human's general intelligence [4] and fluid intelligence [19,24,55,56]. Early RPM-solving systems employ symbolic representation based on hand-designed features and assume access to the underlying logics [4,31,32,33]. Another stream of research on RPM recruits similarity-based metrics to select the most similar answer from the choices [20,30,41,42,43,54]. However, these visual or semantic features are unable to handle uncertainty from imperfect perception, and directly assuming access to the logic operations simplifies the problem. Recently proposed data-driven approaches arise from the availability of large datasets: Santoro et al . [51] extend a pedagogical RPM generation method [59], whereas Zhang et al . [68] use a stochastic image grammar [75] and introduce structural annotations in it, which Hu et al . [21] further refine to avoid shortcut solutions by statistics in candidate panels. Despite the fact that RPM intrinsically requires one to perform abstraction, algebraization, induction, and generalization, Belief Inference Fig. 1. An overview of the ALANS learner. For an RPM instance, the neural visual perception module produces the belief states for all panels: an object CNN extracts object attribute distributions for each image region, and a belief inference engine marginalizes them out to obtain panel attribute distributions. For each panel attribute, the algebraic abstract reasoning module transforms the belief states into matrix-based algebraic representation and induces hidden operators by solving inner optimizations. The answer representation is obtained by executing the induced operators, and the choice most similar to the prediction is selected as the solution. An example of the underlying discrete algebra and its correspondence is also shown on the right.
existing methods bypass such cognitive requirements using a single feedforward pass in connectionist models: Santoro et al . [51] use a relational module [52], Steenbrugge et al . [57] augment it with a VAE [28], Zhang et al . [68] assemble a dynamic tree, Hill et al . [18] arrange the data in a contrastive manner, Zhang et al . [70] propose a contrast module, Zhang et al . [74] formulate it in a studentteacher setting, Wang et al . [58] build a multiplex graph network, Hu et al . [21] aggregate features from a hierarchical decomposition, and Wu et al . [64] apply a scattering transformation to learn objects, attributes, and relations. Recently, Zhang et al . [71] employ a neuro-symbolic design but requires full knowledge over the hidden relations to perform abduction. While our work adopts the visual perception module and employs a similar training strategy from Zhang et al . [71], the ALANS learner manages to induce the hidden relations, enabling on-the-fly relation induction and systematic generalization on relational learning. The recent work of Neural Interpreter (NI) [47] is a complementary neural approach to our method: Although both NI and ALANS decompose the reasoning process into sub-components and aggregate them, NI focuses more on compositionality, routing new input via different paths of learned modules to generalize, whereas ALANS more on induction, enabling a learned module to adapt on the fly.

The ALANS Learner
In this section, we introduce the ALANS learner for the RPM problem. In each RPM instance, an agent is given an incomplete 3ˆ3 panel matrix with the last entry missing and asked to induce the operator hidden in the matrix and choose from eight choice panels one that follows it. Formally, let the answer variable be denoted as y, the context panels as tI o,i u 8 i"1 , and choice panels as tI c,i u 8 i"1 . Then the problem can be formulated as estimating P py | tI o,i u 8 i"1 , tI c,i u 8 i"1 q. According to the common design [4,51,68], there is one operator that governs each panel attribute. Hence, by assuming independence among attributes, we propose to factorize the probability of P py " n | tI o,i u 8 i"1 , tI c,i u 8 i"1 q as where y a denotes the answer selection based only on attribute a and T a the operator on a.
Overview As shown in Fig. 1, the ALANS learner decomposes the process into perception and reasoning: the neural visual perception frontend is adopted from Zhang et al . [71] and extracts the belief states from each of the sixteen panels, whereas the algebraic abstract reasoning backend views an instance as an example in an abstract algebra structure, transforms belief states into algebraic representation by the representation theory, induces the hidden operators, and executes the operators to predict the representation of the answer. Therefore, in Eq. (1), the operator distribution is modeled by the fitness of an operator and the answer distribution by the distance between the predicted representation and that of a candidate.

Neural Visual Perception
We follow the design in Zhang et al . [71] and decompose visual perception into an object CNN and a belief state inference engine. Specifically, for each panel, we use a sliding window to traverse the spatial domain of the image and feed each image region into an object CNN. The CNN has four branches, producing for each region its object attribute distributions, including objectiveness (if the region contains an object), type, size, and color. The belief inference engine summarizes the panel attribute distributions (over position, number, type, size, and color) by marginalizing out all object attribute distributions (over objectiveness, type, size, and color). As an example, the distribution of the panel attribute of Number can be computed as such: for N image regions and their predicted objectiveness where P pr o j q denotes the jth region's estimated objectiveness distribution, and R o is a binary sequence of length N that sums to k. All panel attribute distributions compose the belief state of a panel. In the following, we denote the belief state as b and the distribution of an attribute a as P pb a q. For more details, please refer to Zhang et al . [71] and the Appendix.

Algebraic Abstract Reasoning
Given the belief states of both context and choice panels, the algebraic abstract reasoning backend concerns the induction of hidden operators and the prediction of answer representation for each attribute. The fitness of induced operators is used for estimating the operator distribution and the difference between the prediction and the choice panel for estimating the answer distribution.
Algebraic Underpinning Without loss of generality, here we assume row-wise operators. For each attribute, under perfect perception, the first two rows in an RPM instance provide snapshots into an example of group [15] constrained to an integer-indexed set, a simple algebra structure that is closed under a binary operator. To see this, note that an accurate perception module would see each panel attribute as a deterministic set element. Therefore, RPM instances with unary operators, such as progression, are group examples with special binary operators where one operand is constant. Instances with binary operators, such as arithmetics, directly follow the group properties. Those with ternary operators are ones defined on a three-tuple set from rows.
< l a t e x i t s h a 1 _ b a s e 6 4 = " g y R l S F L N g D X F R B Y x f 4 F K P v 6 j g a w = " > A A A B 7 3 i c b Z D L S s N A F I Y n 9 V b r r e r S z d A i V I S S W P C y C 7 p  Algebraic Representation A systematic algebraic view allows us to felicitously recruit ideas in the representation theory [23] to glean the hidden properties in the abstract structures: it makes abstract algebra amenable by reducing it onto linear algebra. Following the same spirit, we propose to lift both the set elements and the hidden operators to a learnable matrix space. To encode the set element, we employ the Peano Axiom [46]. According to the Peano Axiom, an integer-indexed set can be constructed by (1) a zero element (0), (2) a successor function (Sp¨q), and (3) the principle of mathematical induction, such that the kth element is encoded as S k p0q. Specifically, we instantiate the zero element as a learnable matrix M 0 and the successor function as the matrix-matrix product parameterized by M .
In an attribute-specific manner, the representation of an attribute taking the kth value is pM a q k M a 0 . For operators, we consider them to live in a learnable matrix group of a corresponding dimension, such that the action of an operator on a set can be represented as matrix multiplication. Such algebraic representation establishes an isomorphism between the matrix space and the abstract algebraic structure: abstract elements on the algebraic structure have a bijective mapping to/from the matrix space, and inducing the abstract relation can be reduced to solving for a matrix operator. See Fig. 2 for a graphical illustration of the isomorphism.
Operator Induction Operator induction concerns about finding a concrete operator in the abstract algebraic structure. By the property of closure, we formulate it as an inner-level regularized linear regression problem: a binary operator T a b for attribute a in a group minimizes ℓ a b pT q defined as where under visual uncertainty, we take the expectation w.r.t. the distributions in the belief states of context panels P pb a o,i q in the first two rows, and denote its algebraic representation as M pb a o,i q. For unary operators, one operand can be treated as constant and absorbed into T . Note that Eq. (3) admits a closed-form solution (please refer to the Appendix for details). Therefore, the operator can be learned and adapted for different instances of binary relations and concluded on the fly. Such a design also simplifies the recent neuro-symbolic approaches, where every single symbol operator needs to be hand-defined [14,36,66,67]. Instead, we only specify an inner-level optimization framework and allow symbolic operators to be quickly induced based on the neural observations, while keeping the semantic interpretability in the neuro-symbolic methods. Therefore, we term such a design semi-symbolic.
The operator probability in Eq. (1) is then modeled by each operator type's fitness, e.g., for binary, Operator Execution To predict the algebraic representation of the answer, we solve another inner-level optimization similar to Eq. (3), but now treating the representation of the answer as a variable: where the expectation is taken w.r.t. context panels in the last row. The optimization also admits a closed-form solution (please refer to the Appendix for details), which corresponds to the execution of the induced operator in Eq. (3). The predicted representation is decoded probabilistically as the predicted belief state of the solution, Answer Selection Based on Eqs. (1) and (4), estimating the answer distribution is now boiled down to estimating the conditional answer distributions for each attribute. Here, we propose to model it based on the Jensen-Shannon Divergence (JSD) of the predicted belief state and that of a choice, P py a " n | T a , tI o,i u 8 i"1 , tI c,i u 8 i"1 q 9 expp´d a n q, where we define d a n as d a n " D JSD pP p p b a | T a q}P pb a c,n qqq.
Discussion Comparing with the possible problem-solving process by humans [4], we argue that the proposed algebraic abstract reasoning module offers a computational and interpretable counterpart to human-like reasoning in RPM. Specifically, the induction component resembles fluid intelligence, where one quickly induces the hidden operator by observing the context panels. The execution component synthesizes an image by executing the induced operator, and the choice most similar to the image is selected as the answer. We also note that by decoding the predicted representation in Eq. (6), a solution can be generated : by sequentially selecting the most probable operator and the most probable attribute value, a rendering engine can directly render the solution. The reasoning backend also enables end-to-end training: by integrating the belief states from neural perception, the module conducts both induction and execution in a soft manner, such that the gradients can be back-propagated and both the visual frontend and the reasoning backend jointly trained.

Training Strategy
We train the entire ALANS learner by minimizing the cross-entropy loss between the estimated answer distribution and the ground-truth selection and an auxiliary loss [51,58,68,71] that shapes the operator distribution from the reasoning engine, i.e., where ℓp¨q denotes the cross-entropy loss, y ‹ the correct choice in candidates, and y a ‹ the ground-truth operator selection for attribute a. The first part of the loss encourages the model to select the right choice for evaluation, while the second part motivates meaningful internal representation to emerge. Compared to Zhang et al . [71], the system requires joint operation from not only a trained perception module θ, but also the algebraic encodings from the zero elements tM a 0 u and the successor functions tM a u, and correspondingly, induced operators T . We notice the three-stage curriculum in Zhang et al . [71] is crucial for such a neuro-semi-symbolic system. In particular, we use λ a to balance the trade-off in the curriculum: in the first stage, we only train parameters regarding objectiveness; in the second stage, we freeze objectiveness parameters and cyclically train parameters involving type, size, and color; in the last stage, we fine-tune all parameters.

Experiments
A cognitive architecture with systematic generalization is believed to demonstrate the following three principles [12,37,38]: (1) systematicity, (2) productivity, and (3) localism. Systematicity requires an architecture to be able to entertain "semantically related" contents after understanding a given thought. Productivity states the awareness of a constituent implies that of a recursive application of the constituent; vice versa for localism.
To verify the effectiveness of an algebraic treatment in systematic generalization, we showcase the superiority of the proposed ALANS learner on the three principles in the abstract spatial-temporal reasoning task of RPM. Specifically, we use the generation methods proposed in Zhang et al . [68] and Hu et al . [21] to generate RPM problems and carefully split training and testing to construct the three regimes. The former generates candidates by perturbing only one attribute of the correct answer while the later modifies attribute values in a hierarchical manner to avoid shortcut solutions by pure statistics. Both methods categorize relations in RPM into three types, according to Carpenter et al . [4]: unary (Constant and Progression), binary (Arithmetic), and ternary (Distribution of Three), each of which comes with several instances. Grounding the principles into learning abstract relations in RPM, we fix the configuration to be 3ˆ3Grid and generate the following data splits for evaluation (please refer to the Appendix for details): -Systematicity: the training set contains only a subset of instances for each type of relation, while the test set all other relation instances. -Productivity: as the binary relation results from a recursive application of the unary relation, the training set contains only unary relations, whereas the test set only binary relations. -Localism: the training and testing sets in the productivity split are swapped to study localism.
We follow Zhang et al . [68] to generate 10, 000 instances for each split and assign 6 folds for training, 2 folds for validation, and 2 folds for testing.

Experimental Setup
We evaluate the systematic generalizability of the proposed ALANS learner on the above three splits, and compare the ALANS learner with other baselines, including ResNet [16], ResNet+DRT [68], WReN [51], CoPINet [70], MXGNet [58], LEN [74], HriNet [21], and SCL [64]. We use either official or public implementations that reproduce the original results. All models are implemented in Py-Torch [45] and optimized using ADAM [27] on an Nvidia Titan Xp GPU. We validate trained models on validation sets and report performance on test sets.

Systematic Generalization
Tab. 1 shows the performance of various models on systematic generalization, i.e., systematicity, productivity, and localism. Compared to results reported in existing works mentioned above, all pure connectionist models experience a devastating performance drop when it comes to the critical cognitive requirements on systematic generalization, indicating that pure connectionist models fail to perform abstraction, algebraization, induction, or generalization needed in solving the abstract reasoning task; instead, they seem to only take a shortcut to bypass them. In particular, MXGNet's [58] superiority is diminishing in systematic generalization. In spite of learning with structural annotations, ResNet+DRT [68] does not fare better than its base model. The recently proposed HriNet [21] slightly improves on ResNet [16] in this aspect, with LEN [74] being only marginally better. WReN [51], on the other hand, shows oscillating performance across three regimes. Evaluated under systematic generation, SCL [64] and CoPINet [70] also far deviate from "superior performance." These observations suggest that pure connectionist models highly likely learn from variation in visual appearance rather than the algebra underlying the problem. Embedded in a neural-semi-symbolic framework, the proposed ALANS learner improves on systematic generalization by a large margin. With an algebra-aware design, the model is considerably stable across different principles of systematic generalization. The algebraic representation learned in relations of either a constituent or a recursive composition naturally supports productivity and localism, while semi-symbolic inner optimization further allows various instances of an operator type to be induced from the algebraic representation and boosts systematicity. The importance of the algebraic representation is made more significant in the ablation study: ALANS-Ind, with algebraic representation replaced by independent encodings and the algebraic isomorphism broken, shows inferior performance. We also examine the performance of the learner with perfect visual annotations (denoted as ALANS-V) to see how the proposed algebraic reasoning module works: the gap despite of accurate perception indicates space for improvement for the inductive reasoning part of the model. In the next section, we further show that the neuro-semi-symbolic decomposition in ALANS's design enables diagnostic tests into its jointly learned perception module and reasoning module. This design is in stark contrast to black-box models.

Analysis into Perception and Reasoning
The neural-semi-symbolic design affords analyses into both perception and reasoning. To evaluate the neural perception and the algebraic reasoning modules, we extract region-based object attribute annotations from the datasets [21,68] and categorize all relations into three types, i.e., unary, binary, and ternary.
Tab. 2 shows the perception module's performance on the test sets in the three regimes of systematic generalization. We note that in order for the ALANS learner to achieve the desired results shown in Tab. 1, ALANS learns to construct the concept of objectiveness perfectly. The model also shows fairly accurate prediction on the attributes of type and size. However, on the texture-related concept of color, ALANS fails to develop a reliable notion on it. Despite that, the general prediction accuracy of the perception module is still surprising, considering that the perception module is jointly learned with ground-truth annotations on answer selections. The relatively lower accuracy on color could be attributed to its larger space compared to other attributes. Tab. 3 lists the reasoning module's performance during testing for the three aspects. Note that on position, the unary operator (shifting) and binary operator (set arithmetics) do not systematically imply each other. Hence, we do not count them as probes into productivity and localism. In general, we notice that the better the perception accuracy on one attribute, the better the performance on reasoning. However, we also note that despite the relatively accurate perception of objectiveness, type, and size, near perfect reasoning is never guaranteed. This deficiency is due to the perception uncertainty handled by expectation in Eq. (3): in spite of correctness when we take arg max, marginalizing by expectation will unavoidably introduce noise into the reasoning process. Therefore, an ideal reasoning module requires the perception frontend to be not only correct but also certain. Computationally, one can sample from the perception module and optimize Eq. (9) using REINFORCE [60]. However, the credit assignment problem and variance in gradient estimation will further complicate training.

In-Distribution Performance
To further evaluate how models perform under the regular Independent and Identically Distributed (I.I.D.) setup, we train the models on the original datasets generated by Zhang et al . [68] and Hu et al . [21] and measure the model accuracy in the test splits. We compare ALANS with published baselines in Tab. 1.
Tab. 4 (left) shows the results on the RAVEN dataset [68]. With the jointly trained vision component, the ALANS learner does not fare better than the best connectionist approaches, making it on par with SCL only. As the dataset is  known to have shortcut solutions, neural approaches like MXGNet and CoPINet could potentially find it easier to solve and hence achieve much superior results in this setup. However, ALANS-V, the variant with a perfect perception component reaches a level of much robustness and accuracy, attaining the best results in grid-like layouts, empirically believed to be the hardest in human evaluation [68]. Tab. 4 (right) shows the results on the I-RAVEN dataset [21]. Apart from ALANS-V's realizing the best performance across all models, we also notice the consistency of performance of the proposed method across datasets, with or without the shortcut issues. All other methods show drastically varying performance, particularly for CoPINet and MXGNet, arguably because of the choice generation strategy that effectively prunes easy solution paths via statistics.
In summary, by analyzing the results from Tabs. 1 and 4 together, we notice that ALANS not only attains reasonable performance on the I.I.D. setup but also generalizes systematically.

Generative Potential
Compared to existing discriminative-only RPM-solving methods, the proposed ALANS learner is unique in its generative potential. As mentioned above, the final panel attribute can be decoded by sequentially selecting the most probable hidden operator and the attribute value. A solution can be generated when equipped with a rendering engine. In Fig. 3, we use the rendering program from Zhang et al . [68] to showcase the generative potential in the ALANS learner. 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " h P e q u C A H 3 r I R w 8 J g W V F U m 0 6 h J a w = " > A A A C Z H i c b V B b S 8 M w G M 3 q b c 7 b d P g k S H E I P o 1 2 4 u 1 t 6 I u P C k 6 F d Y w 0 / T b D k r Q k q W 6 E / h Z f 9 S f 5 B / w d p t 0 Q q x 4 I H M 7 5 b j l h w q j S n v d R c R Y W l 5 Z X q q u 1 t f W N z a 3 6 9 s 6 9 i l N J o E t i F s v H E C t g V E B X U 8 3 g M Z G A e c j g I R x f 5 f 7 D M 0 h F Y 3 G n p w n 0 O R 4 J O q Q E a y s N 6 o 2 g m G E k R J k J J h z L c T a o N 7 2 W V 8 D 9 S / w 5 a a I 5 b g b b l Y s g i k n K Q W j C s F I 9 3 0 t 0 3 2 C p K W G Q 1 Y J U Q Y L J G I + g Z 6 n A H F T f F J s z 9 9 A q k T u M p X 1 C u 4 X 6 s 8 N g r t S U h 7 a S Y / 2 k f n u 5 + J / X S / X w v G + o S F I N g s w W D V P m 6 t j N o 3 A j K o F o N r U E E 0 n t r S 5 5 w h I T b Q M r b U l o f l r p H 2 Y y O 7 8 W C H g h M e d Y R C Y g R Y A m i K g Y m R M / y w 7 L / q T k n 1 i / C P s i x + l 3 t H / J f b v l H 7 e O b 9 v N z u U 8 9 i r a Q w f o C P n o D H X Q N b p B X U T Q F L 2 i N / R e + X T W n Y a z O y t 1 K v O e B i r B 2 f 8 C p U a 8 M A = = < / l a t e x i t > 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " h P e q u C A H 3 r I R w 8 J g W V F U m 0 6 h J a w = " > A A A C Z H i c b V B b S 8 M w G M 3 q b c 7 b d P g k S H E I P o 1 2 4 u 1 t 6 I u P C k 6 F d Y w 0 / T b D k r Q k q W 6 E / h Z f 9 S f 5 B / w d p t 0 Q q x 4 I H M 7 5 b j l h w q j S n v d R c R Y W l 5 Z X q q u 1 t f W N z a 3 6 9 s 6 9 i l N J o E t i F s v H E C t g V E B X U 8 3 g M Z G A e c j g I R x f 5 f 7 D M 0 h F Y 3 G n p w n 0 O R 4 J O q Q E a y s N 6 o 2 g m G E k R J k J J h z L c T a o N 7 2 W V 8 D 9 S / w 5 a a I 5 b g b b l Y s g i k n K Q W j C s F I 9 3 0 t 0 3 2 C p K W G Q 1 Y J U Q Y L J G I + g Z 6 n A H F T f F J s z 9 9 A q k T u M p X 1 C u 4 X 6 s 8 N g r t S U h 7 a S Y / 2 k f n u 5 + J / X S / X w v G + o S F I N g s w W D V P m 6 t j N o 3 A j K o F o N r U E E 0 n t r S 5 5 w h I T b Q M r b U l o f l r p H 2 Y y O 7 8 W C H g h M e d Y R C Y g R Y A m i K g Y m R M / y w 7 L / q T k n 1 i / C P s i x + l 3 t H / J f b v l H 7 e O b 9 v N z u U 8 9 i r a Q w f o C P n o D H X Q N b p B X U T Q F L 2 i N / R e + X T W n Y a z O y t 1 K v O e B i r B 2 f 8 C p U a 8 M A = = < / l a t e x i t > 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " h P e q u C A H 3 r I R w 8 J g W V F U m 0 6 h J a w = " > A A A C Z H i c b V B b S 8 M w G M 3 q b c 7 b d P g k S H E I P o 1 2 4 u 1 t 6 I u P C k 6 F d Y w 0 / T b D k r Q k q W 6 E / h Z f 9 S f 5 B / w d p t 0 Q q x 4 I H M 7 5 b j l h w q j S n v d R c R Y W l 5 Z X q q u 1 t f W N z a 3 6 9 s 6 9 i l N J o E t i F s v H E C t g V E B X U 8 3 g M Z G A e c j g I R x f 5 f 7 D M 0 h F Y 3 G n p w n 0 O R 4 J O q Q E a y s N 6 o 2 g m G E k R J k J J h z L c T a o N 7 2 W V 8 D 9 S / w 5 a a I 5 b g b b l Y s g i k n K Q W j C s F I 9 3 0 t 0 3 2 C p K W G Q 1 Y J U Q Y L J G I + g Z 6 n A H F T f F J s z 9 9 A q k T u M p X 1 C u 4 X 6 s 8 N g r t S U h 7 a S Y / 2 k f n u 5 + J / X S / X w v G + o S F I N g s w W D V P m 6 t j N o 3 A j K o F o N r U E E 0 n t r S 5 5 w h I T b Q M r b U l o f l r p H 2 Y y O 7 8 W C H g h M e d Y R C Y g R Y A m i K g Y m R M / y w 7 L / q T k n 1 i / C P s i x + l 3 t H / J f b v l H 7 e O b 9 v N z u U 8 9 i r a Q w f o C P n o D H X Q N b p B X U T Q F L 2 i N / R e + X T W n Y a z O y t 1 K v O e B i r B 2 f 8 C p U a 8 M A = = < / l a t e x i t > 7 < l a t e x i t s h a 1 _ b a s e 6 4 = " h P e q u C A H 3 r I R w 8 J g W V F U m 0 6 h J a w = " > A A A C Z H i c b V B b S 8 M w G M 3 q b c 7 b d P g k S H E I P o 1 2 4 u 1 t 6 I u P C k 6 F d Y w 0 / T b D k r Q k q W 6 E / h Z f 9 S f 5 B / w d p t 0 Q q x 4 I H M 7 5 b j l h w q j S n v d R c R Y W l 5 Z X q q u 1 t f W N z a 3 6 9 s 6 9 i l N J o E t i F s v H E C t g V E B X U 8 3 g M Z G A e c j g I R x f 5 f 7 D M 0 h F Y 3 G n p w n 0 O R 4 J O q Q E a y s N 6 o 2 g m G E k R J k J J h z L c T a o N 7 2 W V 8 D 9 S / w 5 a a I 5 b g b b l Y s g i k n K Q W j C s F I 9 3 0 t 0 3 2 C p K W G Q 1 Y J U Q Y L J G I + g Z 6 n A H F T f F J s z 9 9 A q k T u M p X 1 C u 4 X 6 s 8 N g r t S U h 7 a S Y / 2 k f n u 5 + J / X S / X w v G + o S F I N g s w W D V P m 6 t j N o Ground-truth relations are also listed. Note the generated results do not look exactly the same as the correct candidate choices due to random rotations during rendering, but they are semantically correct.

Conclusion and Limitation
In this work, we propose the ALgebra-Aware Neuro-Semi-Symbolic (ALANS) learner, echoing a normative theory in the connectionist-classicist debate that an algebraic treatment in a cognitive architecture should improve a model's systematic generalization ability. In particular, the ALANS learner employs a neural-semi-symbolic architecture, where the neural visual perception module is responsible for summarizing visual information and the algebraic abstract reasoning module transforms it into algebraic representation with isomorphism established by the Peano Axiom and the representation theory, conducts operator induction, and executes it to arrive at an answer. In three RPM domains reflective of systematic generalization, the proposed ALANS learner shows superior performance compared to other pure connectionist baselines. The proposed ALANS learner also bears some limitations. For one thing, we make the assumption in our formulation that relations on different attributes are independent and can be factorized. This assumption is not universally correct and could potentially lead to failure in more complex reasoning scenarios when attributes are correlated. For another, we assume a fixed and known space for each attribute in the perception module, while in the real world the space for one attribute could be dynamically changing. In addition, the reasoning module is sensitive to perception uncertainty as has already been discussed in the experimental results. Besides, the gap between perfection and the status quo in reasoning remains to be filled. In this work, we only show how the hidden operator can be induced with regularized linear regression via the representation theory. However, more elaborate differentiable optimization problems can certainly be incorporated for other problems.
With the limitation in mind, we hope that this preliminary study could inspire more research on incorporating algebraic structures into current connectionist models and help address challenging modeling problems [65,69,73,76].