Computer Science Mathematics
Computer Science-2021.gz:143246
Computer Science
713acb000509d7823cf915fd9099d68cf6dfed65
Curvature-Aware Derivative-Free Optimization

We propose a new line-search method, coined Curvature-Aware Random Search (CARS), for derivative-free optimization. CARS exploits approximate curvature information to estimate the optimal step-size given a search direction. We prove that for strongly convex objective functions, CARS converges linearly if the search direction is drawn from a distribution satisfying very mild conditions. We also explore a variant, CARS-NQ, which uses Numerical Quadrature instead of a Monte Carlo method when approximating curvature along the search direction. We show CARS-NQ is effective on highly non-convex problems of the form $f = f_{\mathrm{cvx}} + f_{\mathrm{osc}}$ where $f_{\mathrm{cvx}}$ is strongly convex and $f_{\mathrm{osc}}$ is rapidly oscillating. Experimental results show that CARS and CARS-NQ match or exceed the state-of-the-arts on benchmark problem sets.


Introduction
We consider minimizing a function f : R d → R, using only function queries f (x), without access to gradients or directional derivatives. This setting is known as derivative-free optimization (DFO). DFO has a long tradition and recently revived for reinforcement learning [1,2,3], hyperparameter tuning [4], and adversarial attacks on neural network based classifiers [5,6]. Common to all such applications is the assumption that evaluating f (x) is expensive, time consuming or inconvenient. Thus, it is desirable for DFO algorithms to use as few function queries as possible.
Classical approaches to DFO include the Nelder-Mead simplex method [7], direct search methods [8] and model-based methods [9]. However, such approaches scale poorly with the problem dimension d. Due in part to the demands of large-scale applications in machine learning, zeroth-order (ZO) approaches to DFO have been gaining increasing attention [10]. Such approaches mimic first-order methods such as gradient descent, but with all derivative information approximated using function queries. ZO methods are iterative, and at each iteration the algorithm selects a direction u k and then takes a step x k+1 = x k + α k u k . While selection of u k has been well studied (see [11] and citations therein), this paper focuses on the selection of α k . Choosing α k intelligently yields faster convergence and may improve the overall query efficiency. When u k ≈ ∇f (x k ) techniques from first-order optimization such as backtracking line search may be employed [12]. However, obtaining a sufficiently accurate approximation to the gradient requires Ω(d) queries per iteration [11], which is impractical when d is large. On the other hand, when u k is a random vector, α k may be taken to be Methods employing (c) typically converge in fewer iterations but have a higher per-iteration query complexity. On the other hand methods using (a) or (b) have fixed per-iteration query complexity but may require many more iterations. We propose Curvature-Aware Random Search (CARS), which adaptively selects α k by approximately applying a single step of Newton's method to (1). Thus CARS enjoys the rapid convergence of (c) yet has the fixed per-iteration query complexity of (a) and (b). In fact, each iteration of CARS requires only three function queries and is guaranteed to decrease the objective function value. CARS is scale-invariant (i.e. minimizing λf (x) is identical to minimizing f (x) for any λ > 0) and works with a wide variety of rules for selecting u k .
The step-size α k selected by CARS may be viewed as a ratio of two stochastic estimators, both of which are high-variance. Inspired by [15], we propose to use numerical quadrature to construct more accurate estimators. The resulting algorithm, dubbed CARS-NQ, performs significantly better on problems of the form f = f cvx + f osc where f cvx is strongly convex and f osc is rapidly oscillating (see Section 5). Heuristically, this is because numerical quadrature "smoothes out" f osc , so CARS-NQ avoids spurious local minima.
To summarize, our main contributions are: 1. A lightweight yet effective and scale-invariant ZO algorithm, CARS, which guarantees descent on every iteration.
2. Guarantees that, in expectation, CARS requires only O(log(1/ε)) function queries to reach an ε-optimal solution (see Section 1.1) when f is strongly convex. In addition, we show CARS is scale-invariant and affine-invariant with a suitable transformation of the distribution (see Section 4).
3. The CARS-NQ algorithm, a heuristic enhancement to CARS designed for non-convex functions of the form f = f cvx + f osc where f cvx is strongly convex and f osc is rapidly oscillating. 4. Empirical evidence that CARS outperforms other ZO line-search algorithms on challenging synthetic benchmark problems [16] and on the real-world problem of generating black-box adversarial attacks.

Assumptions and Notation
Throughout this paper, we assume f is twice continuously differentiable. We denote g(x) = ∇f (x) and H(x) = ∇ 2 f (x) but emphasize we do not assume access to either of these, except briefly as a step in the theoretical analysis of Section 4. For a fixed initial The Euclidean norm is denoted by · . For a vector u ∈ R d , we use the hat notation for its normalization:û = u/ u . We denote f := min x∈R d f (x) and say x k is an ε-optimal solution if f (x k ) − f ≤ ε. We use D to denote a probability distribution on R d . For any set S ⊂ R d by Unif(S) we denote the uniform distribution over S. By S d−1 we mean the unit sphere, S d−1 := {u : u 2 = 1} ⊂ R d . By e 1 , · · · e d we mean the canonical basis vectors in R d . For two matrices A and B we say Assumption 1 (Strong Convexity and Smoothness). f is µ-strongly convex and L-smooth, i.e. there are µ > 0 and L > 0 such that Under Assumption 1, H(z) is positive definite for all z ∈ Q; hence, the following inner product and induced norm are well-defined for all z ∈ Q: x, y H(z) := H(z)x, y , Assumption 1 implies the following relative smoothness and convexity properties [17,Prop. 2]. Lemma 1.1 (L-Relative Smoothness andμ-Relative Convexity). If f satisfies Assumption 1, then f isμ-relatively convex andL-relatively smooth for someL ≥μ > 0, i.e. for all x, y ∈ Q We also make the following regularity assumption on H.
H is a-Hölder continuous for some a > 0, i.e.
for any unit vectorû.
Hölder continuity reduces to Lipschitz continuity when a = 1. Assumption 2 can be used to refine the relative smoothness and relative convexity constants in a smaller region.

Prior Art
For a comprehensive introduction to DFO, we refer the reader to [9] or the more recent survey articles [18,10]. Here we discuss DFO algorithms based on line search: All of the complexity results discussed below assume noise-free access to f (x). The noisy case is more complicated, see [19]. To the best of the authors' knowledge, the first papers to use this scheme were [20,21], where convergence is discussed under the assumptions that u k is a descent direction 1 direction for all k and (3b) is solved sufficiently accurately. Assuming (3b) is solved exactly, [22] proves this scheme finds an ε-optimal solution in O(log(1/ε)) iterations when f is a quadratic function (see also [23] for a discussion of these results in English).
In [24], O(log(1/ε)) iteration complexity was proved assuming access to an approximate line search oracle that solves (3b) sufficiently accurately, for any strongly convex f , as long as u k are cyclically sampled coordinate vectors. This was extended to more general distributions by [25], albeit at the cost of a more opaque condition the approximate line search oracle must satisfy. Similar ideas can be found in [26,27,28]. More recently, [29] studied (3) under the name Random Pursuit which assumes access to an approximate line search oracle satisfying either additive (α − δ ≤α ≤ α + δ) or multiplicative ((1 − δ)α ≤α ≤ α and sign(α) = sign(α )) error bounds. They show Random Pursuit finds an ε-optimal solution in O(log(1/ε) (resp. O(1/ε)) iterations if f is strongly convex (resp. convex). The use of O(·) above suppresses the dependence of the query complexity on the dimension d. In all results stated, the query complexity scales at least linearly with d. This is unavoidable in DFO for generic f , although see [30,31,6,32,33,34] for recent progress in overcoming this.
We highlight a shortcoming of the aforementioned works: Although they provide essentially optimal bounds on the iteration complexity, they do not bound the query complexity. Indeed, the true query complexity will depend on the inner workings of the solver employed to solve (3b). For example, [29] reports each call to the line search oracle requires an average of 4 function queries when d ≤ 128 which increases to 7 when d = 1024. In contrast, CARS requires only three queries per iteration, independent of d. The recently introduced Stochastic Three Point (STP) method [13,35] also uses only three queries per iteration. However, STP is not scale invariant, and in practice we find its performance compares poorly against CARS (see Section 6).
We note [5] proposes a method, ZOO-Newton, which is essentially CARS with D = Unif({e 1 , · · · , e d }). They demonstrate empirically that ZOO-Newton performs well but do not provide any theoretical guarantees. Our convergence guarantees for CARS imply convergence of ZOO-Newton as a special case.
We are certainly not the first to consider adapting Newton's method to the derivativefree setting. However, obtaining an estimate of the d × d Hessian ∇ 2 f (x k ) for general (i.e. unstructured) f (x) is difficult. Thus, one needs to either use Ω(n 2 ) queries [36] in order to obtain an accurate estimate of ∇ 2 f (x k )-far too many for most applications-or use a high-variance approximation to ∇ 2 f (x k ) [37,38,39,40,41]. CARS sidesteps this dichotomy, as it applies Newton's method to a one-dimensional function. Thus the "Hessian" to be estimated is 1 × 1, i.e. a scalar second derivative.

Curvature-Aware Random Search
Given u k sampled from D, consider the one-dimensional Taylor expansion: CARS selects α k ≈ arg min α T 2 (α; x k , u k ). More specifically, observing that the exact minimizer is u k g k /u k H k u k , CARS uses α k = d r /h r , where d r and h r are finite difference approximations: (We write d r and h r , in place of d r (x k ; u k ) and h r (x k ; u k ) when x k and u k are clear from context.) Thus each iteration of CARS is a zeroth-order analogue of Newton's method applied to f restricted to the line spanned by u k . As is well-known [42], pure Newton's method may not converge globally. Thus, following [17] we augment our method with a fixed step-size 1/L and define: We allow the distribution to be iteration dependent, i.e. we allow for u k to be sampled from D k . In computing d r (x k ) and h r (x k ), CARS queries f at the symmetric points x k + ru k and x k − ru k . There is thus a natural synergy between CARS and the Stochastic Three Point Method [13] where x k+1 is chosen as: We exploit this by using STP as a safeguarding mechanism within CARS. Specifically, after computing the CARS point x CARS , we choose the next iterate as This safeguarding ensures the sequence of iterates selected by CARS is non-increasing: We note that the Nesterov-Spokoiny Random Search (NSRS) algorithm chooses the next iterate as x k+1 = x k − dr L u k . Because CARS uses (approximate) second-order information it may be thought of as a "curvature-aware" version of random search.

6:
Evaluate and store f (x k ± ru k ).

Convergence Guarantees
Before proceeding we list several necessary assumptions on D k . All proofs are deferred to Appendix A. To describe the convergence result concisely, we introduce a function η: By Cauchy-Schwarz inequality, η(g, H; D) ≤ 1 for all g, D, and positive definite H.

(Approximate gradient direction) If D k satisfies
We remark that when D k approximates the gradient or Newton direction η D k does not depend on d, leading to dimension independent convergence rates. When D k is modified via a linear transformation η D k changes in predictable manner: This function will be used to show that Assumption 4 guarantees an enough descent with certain probability. In Appendix A.1, we provide lower bounds of p D for various D. Assumption 4 connects the choice of sampling radius r to properties of D k . From (5) and (6) it is apparent that scaling the sampling radius r and the distribution D k by the same factor β > 0 (i.e. replacing (r, D k ) by (β −1 r, βD k )) will not affect CARS. Hence, we use the following notion of scale-free radius.
Assumption 4. Let ε > 0 be the target accuracy. The scale-free sampling radius r u satisfies for some γ > 0 with p D (γ) > 0, where Assumption 4 states that one of the upper bounds of r is a constant depending only on the function's convexity and regularity, and the other bound additionally depends on the target accuracy and the sample distribution. For instance, when H is Lipschitz continuous, we have r = O(ε 1/4 ). Also notice that both bounds in Assumption 4 are scale-invariant in the sense that replacing the pair (f, ε) by (λf, λε) for any λ > 0 does not change the bounds because the λ factors cancel out.
We provide the convergence analysis for CARS for strongly convex f . We begin by temporarily assuming access to Exact Derivatives (ED). In particular, at the k-th iteration, ED takes the step without safeguarding: whereL is the relative convexity parameter.
We use Proposition 4.3 to prove the convergence of CARS, i.e. Algorithm 1, by bounding the distance between the update assuming exact derivatives (ED) and the CARS point: Lemma 4.4 (Expected descent of CARS). Let f satisfy Assumptions 1 and 2, and D satisfy Assumption 3. If for any k ≥ 0, then the update scheme (10) gives We cannot choose r k such that (11) holds deterministically, as this would require knowledge of g k . Instead, we show (11) holds for a positive fraction of CARS iterations. When (11) fails the safeguarding mechanism guarantees we still have descent: Carefully quantifying this intuition we show: Theorem 4.5 (Convergence of CARS). Let f satisfy Assumptions 1 and 2, and D satisfy Assumption 3. Choose r k upon u k such that Assumption 4 holds. If Limitations. Currently our theoretical results for CARS are limited to the case where f is strongly convex. The difficulty in extending our analysis to convex or even non-convex functions is analogous to that in extending convergence results for Newton's methods to non-strongly-convex functions. To do so, Newton's method must be modified to include cubic regularization [42] or a trust region [43]. We leave such extensions for future work while noting that empirically CARS performs excellently on non-convex functions (see Section 6).

Incorporating Numerical Quadrature
Recall the Gaussian smoothing of f is defined as As discussed in [14] and elsewhere . This fact is key in proving the convergence of NS random search [14]. The CARS iteration can be interpreted as (approximately) applying a single step of Newton's method to a smoothed, one-dimensional function. To see this, assume D = N (0, I) and decompose u k ∈ D as a product of two independent random variables: . Suppose u k is drawn first and then fixed. Define the one-dimensional functionf (t; x k ,û k ) = f (x k + tû k ) and its Gaussian smoothing: For ease of notation, we omit x k andû k when they are clear from context. Because d r is an In this section we propose a variant of CARS which uses better estimators of G r [f ] (0) and G r [f ] (0). By integration by parts: As these integrals are one-dimensional they may be accurately approximated using as few as three queries using Gauss-Hermite (GH) quadrature. The following is a simple consequence of [44].
be the GH quadrature points and weights. Define Here the constant for the O-notation depends onf and q. See [44] for numerical expressions for t i and w i . Notice that when q is odd, we only need q − 1 new function queries, since t = 0 is one of the quadrature points. Using d NQ r,q , h NQ r,q allows one to obtain extremely accurate approximations toF r (s) andF r (s) while using a much larger r. Although large r is not desirable when f is strongly convex (c.f. Lemma 4.4) it can be very useful when f is non-convex of the form where f cvx is strongly convex while f osc is rapidly oscillating, i.e. f osc (x) = ψ(x) cos(λφ(x)) where ψ and φ are smooth and ψ ∞ < ∞.
Theorem 5.2. Suppose f is non-convex of the form (17). For any fixedû consider Thus a judicious choice of r "smoothes out" the oscillatory part, while still accurately estimating the first and second derivatives of the strongly convex part. Each iterate of CARS with Numerical Quadrature (CARS-NQ, see Algorithm 2) effectively applies (one step of) Newton's method to f cvx , while ignoring f osc . This suggests, but does not prove, CARS-NQ is robust towards local minima induced by f osc and will converge towards the global minimum of f (assuming it is close to the global minimum of f cvx ). Our empirical results (see Section 6) support this intuition. [45,15] also suggest using GH quadrature in DFO. However the underlying principle of their proposed algorithms, DGS and AdaDGS, is quite different from CARS-NQ. They apply GH quadrature in each coordinate direction and stack these estimates into a vector which they call the Directional Gaussian Smoothed (DGS) gradient. The DGS gradient is not easily interpretable as the gradient of a function, hence analyzing the convergence of AdaDGS and DGS is tricky. Nevertheless in practice both AdaDGS and DGS converge in relatively few iterations, although we note their per-iteration query complexity is Ω(d), making them unsuitable for high-dimensional problems.
Bounding the relative smoothness. Both CARS and CARS-NQ require the knowledge of the relative smoothness constant,L, so as to set the step-size appropriately. Using NQ a proper value ofL may be approximated.
Although Proposition 5.3 does not upper bound the relative smoothness parameter over the whole domain, it guarantees the same inequality at the point where we need it. As in (13) and (14), also the higher order derivatives of G r [f ] can be easily estimated using NQ. Hence, we suggest the following approximation to the relative smoothness parameter for each direction u:L u ≈ where m NQ r,q denote the NQ estimator of G r [f ] (0). The derivation can be found in Appendix A.7.
Algorithm 2 CARS with Numerical Quadrature (CARS-NQ) 1: Input: x 0 : initial point; r: sampling radius; q: number of quadrature points. 2: Get the oracle f (x 0 ). 3: for k = 0 to K do 4: Sample u k from D.

Experimental Results
For a detailed description of all experimental settings and hyperparameters see Appendix B. The code for all the experiments can be found in https://github.com/bumsu-kim/CARS.

Performance Profiles
Performance profiles [46] are a commonly used tool for comparing the performance of multiple algorithms over a suite of test problems. Performance profiles tend to be more informative than single-dimensional summaries (e.g. average number of iterations required to solve a problem). Formally, consider fixed sets of problems P and algorithms S. For each p ∈ P and s ∈ S the performance ratio r p,s is defined by where t p,s is the number of function queries required for s to solve p. This is the relative performance of s on p compared to the best algorithm in S for p. The performance profile of s, ρ s : Therefore, ρ s (1) is the fraction of problems for which s performs the best, while ρ s (τ ) for large τ measures the robustness of s. For all τ , a higher value of ρ s (τ ) is better. We use a log-scale on the horizontal access when plotting ρ s (τ ).

Problems with Highly Oscillatory Noise
We verified that CARS-NQ works especially well when f = f cvx + f osc , as described in Section 5. We added highly oscillatory noise of the form f osc (x) = ψ d i=1 (1 − cos(φ i x i )) to (i) the same Moré-Garbow-Hillstrom problems with additional noise, and (ii) a convex quartic function where α i > 0 and A 0. Details on the noise level ψ and the frequency φ are presented in Appendix B. Performance profiles for (i) and the function values versus number of evaluations for (ii) in Figure 2 illustrates that CARS-NQ works well when faced with oscillatory noise.

Black-box Adversarial Attacks
Suppose N is a convolutional neural network (CNN) trained to classify images. The problem of generating small perturbations x that, when added to a natural image x nat , fool the classifier (i.e. N (x nat + x) = N (x nat )) is known as finding an adversarial attack [49]. As described in [5], when no access to the internal workings of the CNN is available, this problem becomes a black-box, or derivative free, optimization problem. In order to ensure the attacked image x nat + x appears natural, a pixel-wise bound x ∞ ≤ ε atk is usually enforced. CARS showed state-of-the-art performance in generating black-box adversarial attacks for N trained on the MNIST digit classification dataset [50]. In our experiments, N is a two-layer CNN achieving 99% test accuracy on unperturbed images. We use ε atk = 0.2 and consider all 10,000 images from the test set of MNIST. We consider an attack a success if it fools N before a budget of 10,000 queries is met. The success rates, median and average queries for successful attacks are shown in Table 1. The results from ZOO [5], PGD-NES [51], and ZOHA-type algorithms [38] are cited from [38]. As pointed out in Section 4, the choice of sampling directions for CARS is not restrictive. Hence we used a similar initialization and D as the Square Attack [52], which is known to be particularly well-suited for attacking CNN models. Detailed settings and visualization of attacked images can be found in Appendix B.

Conclusion and Future Work
We present a novel DFO algorithm, namely CARS, which enjoys a linear rate of convergence for strongly convex f . We incorporate numerical quadrature with CARS to obtain CARS-NQ, which is a superior solver when f is perturbed by non-convex, oscillatory noise. As both algorithms make a predetermined number of queries per iteration, their convergence behaviors translate into state-of-the-art query efficiency. Our experimental results verify this query efficiency, even when f is not strongly convex. Future work will extend our analysis of CARS to such non-convex and merely convex functions. [

A.1 Estimation and Lower Bounds of p D (γ) for Various Distributions
Recall that for a distribution D in R n , p D : [0, ∞) → [0, 1] is defined as Lemma A.1 (Estimation and Lower Bounds of p D (γ) for Various Distributions). We have the following estimation and lower bounds of p D .

(Uniform on sphere/ Gaussian) When
where I is the incomplete Beta function, and B is the Beta function and y) is the Beta function.
Remark A.2. For the Gaussian distribution and uniform distribution on the sphere, the estimate suggested in Lemma A.1 seems hard to use. In fact, however, simple numerical calculation reveals that it has even better bound than the random coordinate direction scenario when γ = 1/ √ d. In Figure 4, we provide the values of , 000 turns out to be around 0.31733.
Remark A.3. p D plays an important role in the proof of Theorem 4.5. However, we don't need to take an infimum over all v ∈ S d−1 for the proof. In fact, what we need for the proof of Theorem 4.5 is that γ ≤ |û kĝ k | with a positive probability, so we only need to analyze And we can replace p D (γ) in the Theorem by p D k (γ, g k ). Then we also get the following results: then p D (γ,ĝ k ) ≥ p(γ). In particular, when u k = cg k for some scalar c (i.e. gradient direction), we have p D (1, g k ) = 1.

A.2 Lower Bounds of η for Various Distributions
Recall that the function η was defined in Section 1.1 as

(Approximate gradient direction) If D k satisfies
then η(g k , H k ; D k ) ≥ βµ/L.

(Newton direction)
When u k is parallel to H(x k ) −1 g(x k ) with probability 1, we have η(g k , H k ; D k ) = 1.
Proof. First notice that u Hu ≤ L u 2 and g H −1 g ≤ µ −1 g 2 . Thus we can lower bound the expectation in terms of a more familiar quantity, the cosine of the angle between u and g: 1. When D = N (0, I d ) or Unif(S d−1 ), we can replace g by the standard basis vector e 1 by symmetry. Denoting the first entry of u by u 1 and using (21), 2. when D = Unif({e 1 , · · · , e d }), Hence, again from (21), we have the same lower bound µ/(dL) as in the previous case.
3. Similarly, when u is an approximate gradient direction, in the sense that holds, then it immediately follows that η ≥ βµ/L. In particular, when u = cg for some scalar c (i.e. gradient direction), then η ≥ µ/L.

4.
Finally, when u is the Newton direction, that is, u is parallel to H −1 g with probability 1, then u g = u Hu = g H −1 g, and η(g, H; D) = 1. This is the best possible η, since η ≤ 1 from Cauchy-Schwarz inequality.
This finishes the proof. Proof.
This finishes the proof.
Remark A.4. From Lemma 4.2, when B ≈ H −1/2 , and D 0 is either Unif({e 1 , · · · , e d }), N (0, I d ), or Unif(S d−1 ), we have η ≈ 1/d. Note that this means we can (almost) remove the curvature factor µ/L from η and thus the convergence speed are less affected by the curvature. Hence, as a future work, we can consider an algorithm that adaptively changes the sampling distribution so that η increases up to 1/d.

A.3 Proof of Proposition 4.3
Recall that the iteration of CARS assuming exact derivatives (ED) is Proof. Fromμ-relative strong convexity we have where the second inequality follows by taking x = x − x k and c =μ in the following general inequality: arg min Lemma 9 of [17]. Rearranging (22) yields FromL-relative smoothness and Lemma 5 of [17] f where M k := uu u H k u . We can simplify the last term: Now let E k [·] := E[·|x k ] and take the conditional expectation of both sides of (24): (23) and the assumed lower bound on η(g k , H k ; D k )) Subtracting f from both sides yields the desired result.
Remark A.5. Proposition 4.3 is superficially similar to Corollary 1 of [17]. However in [17] the conditions on D are quite restrictive. Specifically D must be a distribution over a set of vectors {d 1 , · · · , d n } where the probability of selecting d i depends on the Hessian of f . In contrast, Proposition 4.3 allows for a wide array of D.
A.4 Proof of Lemma 4.4 for some k ≥ 0, then the update scheme (10) gives Proof. In this proof, for notational convenience let d 0 = g k u k andd 0 = g kû k = d 0 / u k for the first-order directional derivative, and h 0 = u k H k u k andĥ 0 =û k H kûk = h 0 / u k 2 for the second-order. Also define a quadratic function q(t) = −d 0 t −Lh 0 t 2 /2. Since we assume 0 < r k ≤ C 1 |d 0 | 1/(1+a) , d 0 is not 0. FromL-relative smoothness, As in the exact derivatives case, the maximizer of q is t = −d 0 /(Lh 0 ), with corresponding maximum q(t ) = d 2 0 /(2Lh 0 ) = g k M k /(2L), for M k as defined in A.3. We first show t r := −d r /(Lh r ) approximates t well enough so that q(t r ) ≥ q(t )/2. Observe that if holds for some 0 < c < 1, then by completing the square in q(t): Hence, it suffices to show Let e d and e h be the absolute errors in estimating d 0 and h 0 , respectively, i.e. e d = |d 0 − d r | and e h = |h 0 − h r |. Then and thus, for (27) we only need to prove Now we bound e d and e h using Taylor's theorem and Assumption 2. First, we have and thus By Assumption 2, and therefore, Similarly, for the second-order directional derivative, Again using Assumption 2, One sufficient condition for (28) is and simultaneously Using (29), we see that (31) is met when Similarly, from (30) and h 0 =ĥ 0 u k 2 ≥ µ u k 2 , we get (32) when Therefore, (25) implies (28).
If A k does not occur then, as CARS is non-increasing, Thus completes the proof.

B.1 Experimental Settings
In this Section, we list the hyperparameters we used for each experiment. The code for all experiments can be found in https://github.com/bumsu-kim/CARS. We ran experiments on several machines to distribute the load. We used Intel i5-9400F with Nvidia RTX 2060 and i9-9940X with two RTX 2080.
Moré-Garbow-Hillstrom Problems. The Moré-Garbow-Hillstrom Problem set consists of 34 non-convex smooth functions, where the problem dimension lies between 2 and 100. This experiment is conducted in Matlab.
We consider a problem solved when f ( The target accuracies used here are ε = 10 −1 , 10 −3 and 10 −5 . We used the recommended x 0 for each problem. For CARS and CARS-NQ, we used the sampling radius r k = 0.01/ √ k + 1,L = 2. For CARS-NQ, the number of quadrature points q = 5 is used. Having larger q increases the accuracy of the approximation of smoothed derivatives, but it also increases the cost (i.e. function queries required) at each step. Recall from Section 5 that using odd q is always more efficient than even q. Since q = 3 is essentially less accurate than CARS, we used several values of q ≥ 5 and found that q = 5 is the best.
For AdaDGS [15], we used the code provided by the authors, by implementing the original Python code in Matlab. Some modifications on hyperparameters are made due to the difference in the scale of problem dimension, and the lack of domain width. First, the original AdaDGS code performs experiments on high dimensional problems (e.g. d = 1000), whereas 2 ≤ d ≤ 100 in this experiment. Also, the problems are unconstrained, and x 0 −x varies from order of 10 0 to 10 6 . Thus we used the following modified hyperparameters (following the notation of [15]): 1. The number of points used for line search S = 100, since the suggested value 0.05d(M − 1) is too small for our experiments.
For plotting the performance profile, we set the performance ratio r p,s = r M when p is not solved by s. Having r M = ∞ is ideal, but setting it by a sufficiently large number does not make any difference. We used r M = 10 20 .
Problems with Highly Oscillatory Noise. For the presence of oscillatory noise scenario, we performed two experiments, by adding the noise to Moré-Garbow-Hillstrom Problems and a convex quartic function. For both experiments, we used the same noise For Moré-Garbow-Hillstrom Problems with oscillatory noise, we used the same hyperparameters as in the previous experiment, except for CARS-NQ, where we used 50 times larger r.
For the convex quartic function, the problem dimension d = 20 and α i x 4 i + 1 2 x Ax, α i = 0.01 is used and A is set by dBB /Tr(BB ) + 0.01I d , where B ij ∼ N (0, 1).
Black-box Adversarial Attacks. In this section, we explain the experiment setting for black-box adversarial attacks and also provide the hyperparameters that we used. The CNN model we attack has two 5 × 5 convolutional layers with 6 and 16 output channels, followed by a 4 × 4 convolutional layer with 120 output channels. Then two fully connected layers with 84 and 10 units follows. Between layers we use ReLU, and between convolutional layers we use 2 × 2 max-pooling as well. Finally we apply log softmax to the output layer. The test accuracy of the trained model is 98.99%. For more readable description, let the images have width and height both equal to s. Then the problem dimension d is s 2 . For this particular experiment, we make three modifications to CARS. First, since the problem is highly non-convex (h r < 0 at around 50% of the iteration), we don't compute x CARS when h r < 0 at k-th iteration. The second modification is due to the constraint of the problem. Let F be the feasible set: Inspired by [52], we also compute x bdry = x k − t max d r u k , where t max = max{t > 0 : x k − td r u k ∈ F}. This is especially useful when h r < 0 so we don't compute x CARS . Therefore, CARS now requires either 3 (f (x k ± r k u k ) and f (x bdry )) function queries or 4 (f (x CARS ) in addition) per iteration. To sum up, x k+1 = arg min{f (x k ± r k u k ), f (x CARS ), (x bdry )} if h r > 0 arg min{f (x k ± r k u k ), (x bdry )} otherwise.
Lastly, we perturbed x 0 by adding horizontal stripes. That is, where v i,: = ±(1, · · · , 1) for i = 1, · · · , s and x orig is the image to be attacked. This choice of initialization is found to be very effective in [52]. We use the same sampling distribution as Square Attack [52], which is known to be particularly well-suited for attacking CNN models. Identifying s × s images and d-vectors, this distribution generates a square block of all 1 or −1, i.e. u r+1:r+w,c+1:c+w = ±1, where 0 ≤ r, c ≤ s − w and w is the window size. See [52] for more details.
The window size w is determined by p, the fraction of the pixels changed by attack. Namely, w is the closest positive integer to √ ps, and thus this p is an important hyperparameter in this distribution. In [52], the authors used p = 0.8 in the initial stage and halve the value of p when the number of iteration k reaches 10, 50, 200, 1000, 2000, 4000, 6000 and 8000, respectively. For CARS, we used the initial p = 0.2 and halved it at k = 2, 10, 40, 250, 500, 800, 1200 and 1600. Note that the Square Attack only uses one function query at each iteration and CARS for black-box attack uses 3 or 4 queries, and thus p is being halved at similar stage of the algorithms, in terms of the function queries.
Finally, the sampling radius r k is fixed by 1 for every iteration.

B.2 Visualization of Attacked Images
In this section we present the visualization of the attacked images, comparing them with the original images. To be more clear, the images in Figure 5 have Testset ID (TID)'s from 100 to 109. When an image with label n has TID t, this means it is the t-th n appearing in the test set. See the code for more details on TID.