Zeroth-order Nonconvex Stochastic Optimization: Handling Constraints, High-Dimensionality and Saddle-Points

Introduction

In this work, we propose and analyze algorithms for solving the following stochastic optimization problem

Algorithms available for solving problem (1.1) depend crucially on the constraint set $\mathcal{X}$, along with the structure imposed on the objective function, $f$. Despite decades of work in zeroth-order optimization literature, there still exists several challenges, primarily motivated by contemporary statistical machine learning problems. A majority of the existing zeroth-order algorithms are predominantly analyzed in the low-dimensional unconstrained setting. Furthermore, when $f$ is non-convex, apart from the first-order stationarity result for gradient descent (GD) algorithm in [GL13], other meaningful theoretical results are lacking in the zeroth-order optimization literature. In this work, we provide theoretically sound algorithms to address the following three main drawbacks of existing zeroth-order optimization methods.

The first issue we address is that of constrainted zeroth-order stochastic optimization. For the problem in (1.1), depending on the geometry of the constraint set $\mathcal{X}$, the cost of computing the projection to the set might be prohibitive. In the first-order oracle setting, this lead to the re-emergence of Conditional Gradient (CG) algorithms recently [HK12, Jag13]. But the performance of the CG algorithm under the zeroth-order oracle is unexplored in the literature to the best of our knowledge, both under convex and nonconvex settings. Hence it is natural to ask if CG algorithms, with access to zeroth-order oracle has similar convergence rates compared to zeroth-order GD algorithms for the unconstrained case. To address this question, we propose and analyze in Section 2 a classical version of CG algorithm with zeroth-order information and provide convergence results. We then propose a modification in Section 2.2 that has improved rates, when $f$ is convex. Notably, we demonstrate that with zeroth-order information, the complexity of CG algorithms also depend linearly on the dimensionality, similar to the GD algorithms, thereby facilitating constrained zeroth-order optimization.

Finally, we address the issue of avoiding saddle-points in zeroth-order stochastic optimization. When the function $f$ is non-convex, designing algorithms that avoid saddle-points and converge to local minimizers is challenging, as exemplified by worst-case computational hardness results [MK87, CGT18]. Hence, it is necessary to impose further structure on the problem to obtain meaningful results. A particularly interesting structure on $f$ is the so-called strict saddle property, which necessitates that all local minima are global minima. This structure has regained popularity as several useful stochastic optimization problems in statistical machine learning are shown to posses this property; for example, phase retrieval [SQW18], tensor decomposition [GHJY15], matrix completion and sensing [BNS16, GLM16] and training deep neural networks [KK19]. See also the survey article [SQW15]. Motivated by this, algorithms that avoid saddle-points and converge to second-order stationary points have re-gained popularity as well. Indeed, methods based on exact or in-exact second-order oracle naturally converge to second-order stationary points [NP06, CGT11a, CGT11b, XRKM17, TSJ+17, CDHS18, AZ18]. Furthermore, first-order methods escape saddle points by leveraging an additional noise term in each iteration; for example [GHJY15, JGN+17, RZS+18] and the references therein. But to the best of our knowledge, there is no algorithm for efficiently avoiding saddle-points under zeroth-order oracle information. In this work, we propose a zeroth-order cubic regularized Newton method, that converges efficiently to second-order stationary points with just noisy function evaluations. In order to do so, we interpret the Gaussian smoothing for zeroth-order gradient estimation [NS17], as an instantiation of Stein’s identity [Ste72, Ste81]. Based on this interpretation, we develop provable techniques for estimating the Hessian of a function at a point with just function queries, leveraging higher-order Stein’s identity. Notably, our Hessian estimator is based only on inner-product evaluations thereby having a linear-in-dimension time runtime. We also provide a comprehensive complexity analysis of the proposed algorithm in terms of achieving second-order stationary points.

Our contributions: To summarize the above discussion, in this paper we make the following contributions to the literature on zeroth-order stochastic optimization.

We first analyze a classical version of CG algorithm in the nonconvex (and convex) setting, under access to zeroth-order information and provide results on the convergence rates in the low-dimensional setting. We then propose and analyze a modified CG algorithm in the convex setting with zeroth-order information and show that it attains improved rates in the low-dimensional setting.

Next, we consider a zeroth-order stochastic gradient algorithm in the high-dimensional nonconvex setting and illustrate an implicit regularization phenomenon –the algorithm converges to first-order stationary points with rates that depend only poly-logarithmically on dimensionality. We also propose a truncated zeroth-order stochastic gradient algorithm in the convex setting which also depends only poly-logarithmically on the dimensionally but has improved dependence on the error-tolerance.

Finally, we propose a zeroth-order Stochastic cubic regularized Newton method that avoids saddle points and converges to second-order stationary points efficiently. Our algorithm is based on a novel technique for estimating the Hessian of a function from function queries based on Stein’s identities.

Our contributions extend the applicability of zeroth-order stochastic optimization to the constrained, high-dimensional and non-convex settings and also provide theoretical insights in the form of rates of convergence. A summary of the results is provided in Table 1.

We now list the main assumptions we make in this work. Additional assumptions will be introduced in the appropriate sections as needed. We start with the assumption on the zeroth-order oracle.

It should be noted that in the above assumption, we do not observe $\nabla F(x,\xi)$ and we just assume that it is an unbiased estimator of gradient of $f$ and its variance is bounded. Furthermore, we make the following smoothing assumption about the noisy estimation of $f$.

Function $F$ has Lipschitz continuous gradient with constant $L$, almost surely for any $\xi$, i.e., $\|\nabla F(y,\xi)-\nabla F(x,\xi)\|_{*}\leq L\|y-x\|,$ which consequently implies that $|F(y,\xi)-F(x,\xi)-\langle\nabla F(x,\xi),y-x\rangle|\leq\frac{L}{2}\|y-x\|^{2}.$

It is easy to see that the above two assumptions imply that $f$ also has Lipschitz continuous gradient with constant $L$ since

due the Jensen’s inequality for the dual norm. We now collect some facts about a gradient estimator based on the above zeroth-order information. Let $u\sim N(0,I_{d})$ be a standard Gaussian random vector. For some $\nu\in(0,\infty)$ consider the smoothed function $f_{\nu}(x)={\bf E}_{u}\left[f(x+\nu u)\right]$. Nesterov [NS17] has shown that $\nabla f_{\nu}(x)=$

This relation implies that we can estimate gradient of $f_{\nu}$ by only using evaluations of $f$. In particular, one can define stochastic gradient of $f_{\nu}(x)$ as

which is an unbiased estimator of $\nabla f_{\nu}(x)$ under Assumption 1 since

We leverage the following properties of $f_{\nu}$ due to Nesterov [NS17] in our proofs later, which we replicate below for completeness.

For a Gaussian random vector $u\sim N(0,I_{d})$ we have that

for any $k\geq 2$. Moreover, the following statements hold for any function $f$ whose gradient is Lipschitz continuous with constant $L$.

The gradient of $f_{\nu}$ is Lipschitz continuous with constant $L_{\nu}$ such that $L_{\nu}\leq L$.

We next introduce the Stein’s identity, popular in the statistics and probability theory literature.

Furthermore, when the function $g$ has a twice continuously differentiable Hessian, $\nabla^{2}g(\cdot)$, we have the following (where the Expectation is assumed to exist):

Based on the above theorem, the Gaussian smoothing approach of estimating gradients from function queries proposed by [NS17], is indeed based on Stein’s identity. Indeed, if we let $g(u)=f(x+\nu u)$ in Equation 1.9, it is easy to see that the identity in Equation 1.3 holds by simply evaluating the Gaussian Stein’s identity in Equation 1.9. We elaborate more on this connection and extensions in Section 4.1. We conclude the section, by defining the following criterion which are used to analyze the complexity of our proposed algorithms.

Assume that a solution $\bar{x}\in\mathcal{X}$ as output of an algorithm and a target accuracy $\epsilon>0$ are given. Then:

If $f$ is convex, $\bar{x}$ is called an $\epsilon$-optimal point of problem (1.1) if ${\bf E}[f(\bar{x})]-f(x_{*})\leq\epsilon$, where $x_{*}$ denotes an optimal solution of the problem.

If $f$ is nonconvex, $\bar{x}$ is called an $\epsilon$-stationary point of the unconstrained variant of problem (1.1) if ${\bf E}[\|\nabla f(\bar{x})\|_{*}]\leq\epsilon$. For the constrained case, $\bar{x}$ should satisfies ${\bf E}[\langle\nabla f(\bar{x}),\bar{x}-u\rangle]\leq\epsilon$ for all $u\in\mathcal{X}$.

If $f$ is nonconvex, $\bar{x}$ is called an $\epsilon$-local optima of the unconstrained variant of problem (1.1) if

where for a symmetric matrix $A$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denotes the minimum and maximum eigenvalue.

It should be pointed out that while the above performance measures are presented in expectation form, one can also use their high probability counterparts. Since, convergence results in this case can be obtained by making sub-Gaussian tail assumptions on the output of the zeroth-order oracle and using the standard two-stage process presented in [GL13, LZ16], we do not elaborate more on this approach. Furthermore, note that the aforementioned measures for evaluating the algorithms are from the derivative-free optimization point of view. In the literature on optimization with bandit feedback, the preferred performance measure is the so-called regret of the algorithm [BCB12, Sha13] which may have a different behavior than our performance measures.

Handling Constraints: Zeroth-order Stochastic Conditional Gradient Type Method

The feasible set $\mathcal{X}$ is bounded such that $\max_{x,y\in\mathcal{X}}\|y-x\|\leq D_{\mathcal{X}}$ for some $D_{\mathcal{X}}>0$. Moreover, for all $x\in\mathcal{X}$, there exists a constant $B>0$ such that $\|\nabla f(x)\|\leq B$.

We should point out that under Assumptions 1 and 2, the second statement in Assumption 3 follows immediately by the first one and choosing $B:=LD_{\mathcal{X}}+\|\nabla f(x_{*})\|$. However, we just use $B$ in our analysis for simplicity.

The vanilla ZSCG method is formally presented in Algorithm 1 and a few remarks about it follows.

First, note that this algorithm differs from the classical CG method in estimating the gradient using zeroth-order information and in outputting a random solution from the generated trajectory. This randomization scheme is the current practice in the literature to provide convergence results for nonconvex stochastic optimization (see e.g., [GL13, RSPS16]). Second, $\bar{G}_{\nu}^{k}$ is the averaged variant of the gradient estimator presented in Subsection 1.1 and is still an unbiased estimator of $\nabla f_{\nu}(z_{k-1})$. Moreover, it can be easily seen that it has a reduced variance with respect to the individual estimators i.e.,

We emphasize that the use of the above variance reduction technique in stochastic CG methods is standard and has been previously proposed and leveraged in several works (see e.g., [LZ16, HL16, RSPS16, MHK18a, MHK18b, Gha18]). Indeed, when exact gradient is not available, an error term appears in the convergence analysis which should converge to at a certain rate as the algorithm moves forward. Hence, the choice of $m_{k}$ plays a key role in the convergence analysis of Algorithm 1. $\bar{G}_{\nu}^{k}$ can be also viewed as a biased estimator for $\nabla f(z_{k-1})$. Finally, since $f$ is possibly nonconvex, we need a different criteria than the optimality gap to provide convergence analysis of Algorithm 1. The well-known Frank-Wolfe Gap given by

has been widely use in the literature to show rate of convergence of the CG methods when $f$ is convex (see e.g., [FW56, DR70, Hea82]). In this case, it is easy to see that

When $f$ is nonconvex, this criteria is still useful since $\langle\nabla f(z_{k-1}),z_{k-1}-u\rangle\leq g_{{}_{\mathcal{X}}}(z_{k-1}),~{}\forall u\in\mathcal{X}$, which implies that one can obtain an approximate stationary point of problem (1.1) by minimizing $g_{{}_{\mathcal{X}}}^{k}$, in the view of Definition 1.1. Note that in our setting, this quantity is not exactly computable and it is only used to provide convergence analysis of Algorithm 1 as shown in the next result.

Let $\{z_{k}\}_{k\geq 0}$ be generated by Algorithm 1 and Assumptions 1, 2, and 3 hold.

Let $f$ be nonconvex, bounded from below by $f^{*}$, and let the parameters of the algorithm be set as

for some constant $B_{L\sigma}\geq\max\{\sqrt{B^{2}+\sigma^{2}}/L,1\}$ and a given iteration bound $N\geq 1$. Then we have

where $R$ is uniformly distributed over $\{1,\ldots,N\}$ and $g_{k}$ is defined in (2.5). Hence, the total number of calls to the zeroth-order stochastic oracle and linear subproblems required to be solved to find an $\epsilon$-stationary point of problem (1.1) are, respectively, bounded by

Let $f$ be convex and let the parameters be set to

where $R$ is random variable from $\{1,\ldots,N\}$ whose probability distribution is given by

Hence, the total number of calls to the zeroth-order stochastic oracle and linear subproblems required to be solved to find and $\epsilon$-optimal solution of problem (1.1) are, respectively, bounded by

In order to prove Theorem 2.1, we need the following result that provides upper bounds for the variance of our gradient estimator.

Let $\bar{G}_{\nu}^{k}$ be computed by (2.1). Then under Assumptions 1, 2 and 3, we have

Proof. [Proof of Lemma 2.1] First note that using (1.8) for function $F$ instead of $f$, under Assumptions 1 and 2, we obtain

Also noting (1.4), (2.4), and the fact that $\|\nabla f_{\nu}\|\leq B$ under Assumption 3, we have

which together with the above relation clearly imply (2.14). We can then obtain (2.15) by noting (1.7) and the fact that

Proof. [Proof of Theorem 2.1] Denoting $\Delta_{k}=\bar{G}_{\nu}^{k}-\nabla f(z_{k-1})$, noting (1.2), (2.3), and (2.5), we have

where the last inequality follows from boundedness of the feasible set, (2.5), and the fact that

due to the optimality condition of (2.2). Taking expectation from both sides of the above inequality, summing them up, rearranging the terms, and noting Lemma 2.1, we obtain

Hence, choosing $\alpha_{k}=\alpha_{1}$ and $m_{k}=m_{1}$ for all $k\geq 1$, and noting that $R$ is a uniform random variable, we have

which together with (2.7) imply (2.8). Hence, (2.9) follows by noting that the total number of calls to the stochastic oracle is bounded by $\sum_{k=1}^{N}m_{k}$.

Now assume that $f$ is convex. Hence, by (2.6) and (2.16), we have

Taking expectation from both sides of the above inequality, dividing them by $T_{k}$, and summing them up, and noting (2.12), we obtain

Combining the above relations, we get (2.11) and (2.13).

Observe that the complexity bounds in (2.9), in terms of $\epsilon$, match the ones obtained in [Gha18, RSPS16, MHK18b] for stochastic CG method with first-order oracle applied to nonconvex problems. For convex problems, similar observation can be made for terms in (2.13) which match the ones in [HL16, Gha18]. Note that the linear dependence of our complexity bounds on $d$ is unimprovable due to the lower bounds for zeorth-order algorithms applied to convex optimization problems [DJWW15]. We conjecture that this is also the case for nonconvex problems.

2 Improved Rates for Convex Problems

Our goal in this subsection is to improve the complexity bounds of the ZCSG method when $f$ is convex. Recall that the ZSCG method presented in Section 2.1 involves two main steps: the gradient evaluation step and the linear optimization step. Motivated by [LZ16], we now propose a modified algorithm that allows one to skip the gradient evaluation from time to time. Notice that, as our gradients are estimated by calling the zeroth-order oracle, this directly reduces the number of calls to the zeroth-order oracle. We first state a subroutine in Algorithm 2 used in our modified algorithm.

Note that Algorithm 2 is indeed the zeroth-order conditional gradient method for inexactly solving the following quadratic program

which is the standard subproblem of stochastic first-order methods applied to a minimization problem when $g$ is an unbiased stochastic gradient of the objective function at $x$. We now present Algorithm 3 which applies the CG method to inexactly solve subproblems of the stochastic accelerated gradient method. This way of using CG methods can significantly improve the total number of calls to the stochastic oracle. Our next result provides convergence analysis of this algorithm.

Let $\{z_{k}\}_{k\geq 1}$ be generated by Algorithm 3, the function $f$ be convex, and

and for some constants $D_{X}^{0}\geq\|x_{0}-x_{*}\|^{2}$ and $B_{L\sigma}\geq\max\{\sqrt{B^{2}+\sigma^{2}}/L,1\}$. Then under Assumptions 1, 2, and 3, we have

Hence, the total number of calls to the stochastic oracle and linear subproblems solved to find and $\epsilon$-stationary point of problem (1.1) are, respectively, bounded by

Proof. First, note that by (1.2), we have

where the second inequality follows from convexity of $f_{\nu}$, (2.20), and (2.22). Also note that by (2.18) and (2.21), we have

Letting $u=x_{*}$ in the above inequality and multiplying it by $\alpha_{k}$, summing it up with (2.26), and denoting $\bar{\Delta}_{k}=\bar{G}_{\nu}^{k}-\nabla f_{\nu}(w_{k})$, we obtain

subtracting $f_{\nu}(x_{*})$ from both sides of the above inequality, diving them by $\hat{\Gamma}_{k}$, taking expectation, summing them up, noting (1.6) assuming that $\alpha_{1}=1$, $\gamma_{k}\geq 2L\alpha_{k}$, and $\gamma_{k}\alpha_{k}/\hat{\Gamma}_{k}$ is constant for any $k\geq 1$, we obtain

due to (2.2) and (2.29), we obtain (2.24).

Furthermore, note that the function $h_{\gamma}$ defined in Algorithm 2 is indeed negative the FW-gap of the CG method applied to problem (2.19). From classical analysis of the CG method and similar to our result in Theorem 2.1, one can show that the FW-gap is bounded by $LD_{\mathcal{X}}^{2}/T$ if the CG method runs for $T$ iteration. Since the gradient of the objective function in (2.19) is Lipschitz continuous with constant $\gamma$, we have

which together with the choice of $\mu_{k}$ and $\gamma_{k}$ in (2.2), imply that at iteration $k$ of Algorithm 1, we need to run Algorithm 2 for at most $T_{k}=4D_{\mathcal{X}}^{2}N/D_{0}^{2}$ iterations. Therefore, the total number of iterations of Algorithm 2 to find an $\epsilon$-stationary point of problem (1.1) is bounded by $\sum_{k=1}^{N}T_{k}\leq 48LD_{\mathcal{X}}^{2}/\epsilon^{2}$ due to (2.25).

Observe that while the number of linear subproblems required to find an $\epsilon$-optimal solution of problem (1.1) is the same for both Algorithms 1 and 3, the number of calls to the stochastic zeroth-order oracle in Algorithm 3 is significantly smaller than that of Algorithm 1. It is also natural to ask if such an improvement is achievable when $f$ is nonconvex. This situation is more subtle and the answer depends on the performance measure used to measure the rate of convergence. Indeed, we can obtain improved complexity bounds for a different performance measure than the Frank-Wolfe gap with a modified algorithm. However, the complexity bounds are of the same order as (2.9) in terms of the Frank-Wolfe gap for the modified algorithm. For the sake of completeness, we add this algorithm and its convergence analysis in in Section 2.3.

3 Zeroth-order Stochastic Gradient Method with Inexact Updates-Nonconvex case

In this section, we present a zeroth-order stochastic gradient method which applies the CG method to solve the subproblems. This algorithm shares the main idea of Algorithm 3, but for nonconvex problems. We show while this algorithm enjoys better complexity bound than Algorithm 3, it possess the same one when the same performance measure is employed.

Since we are now using the CG method for inexactly solving (2.19), we can provide an alternative termination criterion than the FW-gap given in (2.5) to provide our convergence analysis. In particular, we use the gradient mapping defined as

where $P_{\mathcal{X}}$ is the solution to (2.19). This quantity which has been widely used in the literature as a convergence criteria for solving nonconvex problems (see, e.g., [NY83, Nes04]), plays an analogues role of the gradient in constrained problems. Next result provides some properties for this criteria.

Let $P_{\mathcal{X}}(\cdot)$ be defined in (2.19), $\gamma>0$, and $x\in\mathcal{X}$ are given.

Let $P^{\mu}_{\mathcal{X}}$ be the inexact solution of (2.19) such that

Let $g_{{}_{\mathcal{X}}}(\cdot)$ be the Frank-Wolfe gap defined in (2.5). Then we have

where the last inequality follows from Lipschitz continuity of the Euclidian projection over the feasible set $\Pi_{\mathcal{X}}$. Second, by optimality condition of (2.19), we have

where the last inequality follows from (2.5). Furthermore, (2.32) also implies that

where the last inequality follows from Assumption 3.

Now we are ready to state the main result for the nonconvex case.

Let $\{x_{k}\}$ be generated by Algorithm 4, the function $f$ be nonconvex, and

Then under Assumptions 1, 2, and 3, we have

where $R$ is uniformly distributed over $\{0,\ldots,N-1\}$ and $g_{\mathcal{X}}$ is defined in (2.30). Hence, the total number of calls to the stochastic oracle and linear subproblems solved to find and $\epsilon$-stationary point of problem (1.1) are, respectively, bounded by

Letting $u=x_{k-1}$ in (2.27), summing it up with the above inequality, and denoting $\Delta_{k}=\bar{G}_{\nu}^{k}-\nabla f(x_{k-1})$, we obtain

Taking expectation from the above inequalities, summing them up, re-arranging the terms, and in the view of Lemma 2.1, we have

which together with the facts that $x_{k}=P^{\mu_{k}}_{\mathcal{X}}(x_{k-1},\bar{G}_{\nu}^{k},\gamma_{k})$ and

which implies (2.34). Rest of the proof is similar to that of Theorem 2.2 and hence we skip the details.

We point out that while the complexity bounds in (2.35) are better than those in (2.9) in terms of dependence on the target accuracy $\epsilon$, they have been obtained for a different performance measure. Indeed, if only the Frank-Wolfe gap is considered then it is easy to see that both bounds are of the same order of magnitude due to part c of Lemma 2.2.

Handling High-Dimensionality: Zeroth-order Stochastic Gradient Methods

In this subsection, we consider the zeroth-order stochastic gradient method presented in [GL13] (provided in Algorithm 5 for convenience) and provide a refined convergence analysis for it under the sparsity assumption 1, when $f$ is nonconvex. Our main convergence result for Algorithm 5 under the gradient sparsity assumption is stated below.

Let $\{x_{k}\}_{k\geq 0}$ be generated by Algorithm 5 and stepsizes are chosen such that $\forall k\geq 1$,

for some $\hat{s}\geq s$, $\hat{C}\geq C$ (the universal constant defined in Lemma 3.1), and $D_{0}\geq f(x_{0})-f^{*}$. Assume that $f$ is nonconvex. Then under Assumptions 1, 2, and 4, we have

where $\zeta=\{\xi,u,R\}$ and $R$ is uniformly distributed over $\{0,\ldots,N-1\}$. Hence, the total number of calls to the stochastic oracle (number of iterations) required to find an $\epsilon$-stationary point of problem (1.1), in the view of Definition 1.1, is bounded by

Before proving the theorem, we first present two technical results which play key roles in our convergence analysis.

Let $u\sim N(0,I_{d})$ be a $d$-dimensional standard Gaussian vector. Then for all integer $k\geq 1$ and for some universal constant $C$, we have ${\bf E}\left[\|u\|^{k}_{\infty}\right]\leq C(2\log d)^{k/2}$.

Proof. Let $Z=\|u\|_{\infty}$ and denote by $p(x)$ the standard normal pdf. Note that we have

where we define $x_{d}=\sqrt{2\ln d}$. Now we have

and by l’Hospital’s rule, for large $d$ we have

Hence we have for some universal constant $C$,

The following statements hold for function $f$ and its smooth approximation $f_{\nu}$.

Under Assumptions 1 and 2, gradient of $f$ is Lipschitz continuous with constant $L$ and

where the last inequality follows from Lemma 3.1. Second, noting this lemma again, Assumption 4, and part a), we have

Furthermore, by (1.4), Holder inequality, Lemma 3.1, and under Assumption 4 we have

Proof. [Proof of Theorem 3.1] Noting (1.4), Lemma 3.2.a), and with the notion of $G_{\nu,k}\equiv G_{\nu}(x_{k},\xi_{k},u_{k})$, we have

which after taking expectation imply that

where the last inequality follow from Holder inequality and Lemma 3.2.b). Summing both sides of the above inequality over the iterations and rearranging terms, we get

where $R$ is uniformly distributed over $\{0,\ldots,N-1\}$ since

due to the constant choice of $\gamma_{k}$ in (3.1). Therefore, we have

which together with the choice of smoothing parameter in (3.1) imply (3.2).

Note that the above theorem establishes rate of convergence of Algorithm 5 which only poly-logarithmically depends on the problem dimension $d$, by just selecting the step-size appropriately, under additional assumption that the gradient is sparse. This significantly improves the linear dimensionality dependence of the rate of convergence of this algorithm as presented in [GL13] for general nonconvex smooth problems.

Remarkably, Algorithm 5 does not require any special operation to adapt to the sparsity assumption. This demonstrates an implicit regularization phenomenon exhibited by the zeroth-order stochastic gradient method in the high-dimensional setting when the performance is measured by the size of the gradient in the dual norm. We emphasize that the choice of the performance measure is motivated by the fact that we allow $f$ to be nonconvex. Trivially, the result also applies to the case when $f$ is convex, for the same performance measure.

2 Zeroth-order Stochastic Gradient Method for Convex Problems

We now consider the case when the function $f$ is convex. In this setting, a more natural performance measure is the convergence of optimality gap in terms of the function values. For this situation, we propose and analyze a truncate variant of Algorithm 5 that demonstrates similar poly-logarithmic dependence on the dimensionality. To proceed, in addition to Assumption 4, we also make the following sparsity assumption on the optimal solution of problem (1.1).

Problem (1.1) has a sparse optimal solution $x_{*}$ such that $\|x_{*}\|_{0}\leq s^{*}$, where $s^{*}\approx s$.

Our algorithm for the convex setting is presented in Algorithm 6. Note that this algorithm could be considered as a truncated variant of Algorithm 5 and a zeroth-order stochastic variant of the truncated gradient descent algorithm [JTK14]. In the next result, we present convergence analysis of this algorithm.

Let $\{x_{k}\}_{k\geq 1}$ be generated by Algorithm 5, $f$ is convex, Assumptions 1, 2, 4, and 5 hold. Also assume the stepsizes are chosen such that, $\forall k\geq 1$,

for some $\hat{C}\geq C$, $\hat{s}\geq\max\{s,s^{*}\}$, and $D_{X}^{0}\geq\|x_{0}-x_{*}\|^{2}$.

where $\bar{x}_{N}=\frac{\sum_{k=0}^{N-1}x_{k}}{N}$. Hence, the total number of calls to the stochastic oracle (number of iterations) required to find an $\epsilon$-optimal point of problem (1.1) is bounded by

where the inequality follows from the facts that $|J^{k}|\leq 2\hat{s}+s^{*}$ and $\|G_{\nu,k}^{J^{k}}\|\leq\|G_{\nu,k}\|$. Taking expectation from both sides of the above inequality, summing them up, noting Lemma 3.2, convexity of $f_{\nu}$ (due to convexity of $f$), we have

where the last inequality follows from the fact that $f(x_{k})-f(x_{*})\geq 1/(2Ls)\|\nabla f(x_{k})\|_{2}^{2}$ due to the convexity of $f$ and sparsity of its gradient. Rearranging the terms in the above inequality and noting that $\bar{x}_{N}=\frac{\sum_{k=0}^{N-1}x_{k}}{N}$, we obtain

due to the constant choice of $\gamma_{k}$ in (3.5). Hence, (3.6) follows by using the choice of parameters in (3.5) into the above relation.

While for convex case, similar to the nonconvex case, the complexity of Algorithm 6 depends poly-logarithmically on $d$, it only linearly depends on the choice of $\hat{s}$, facilitating zeroth-order stochastic optimization in high-dimensions under sparsity assumptions.

As discussed in detail in [WDBS18], both Assumption 4 and 5 are implied when we assume the function $f$ depends on only $s$ of the $d$ coordinates. But, both Assumption 4 and 5 are comparatively weaker than that assumption. Furthermore, unlike [WDBS18], we do not make any assumption on the sparsity or smoothness of the second-order derivative of the objective function $f$ for our results.

As mentioned before, [WDBS18] considers only the convex case. Furthermore, their gradient estimator with zeroth-order oracle requires $\text{poly}(s,s^{*},\log d)$ function queries in each iteration whereas our estimator is based on only one function query per iteration. Moreover, [WDBS18] requires computationally expensive debiased Lasso estimators whereas our method requires only simple thresholding operations (for convex case) to handle sparsity.

Handling Saddle-Points: Zeroth-Order Cubic Regularization Method

In this section, we study zeroth-order stochastic cubic regularized Newton method for unconstrained version of Problem 1.1. Throughout this section, we equip our space with the self- dual Euclidean norm, i.e., $\|\cdot\|=\|\|_{2}$. Furthermore, for a matrix $A$, we denote by $\|A\|_{F}$, its Frobenious norm and by $\|A\|$, its operator norm. We also make the following smoothness assumption on the Hessian of the objective function $f$, which is a generalization of the assumption in Equation 1.2.

The function $f$ is twice differentiable and has Lipschitz continuous Hessian i.e., there exists $L_{H}>0$ such that

It can be easily seen that the above assumption is equivalent to

Note that such an assumption in standard in the analysis of second-order optimization techniques [NP06]. We next describe a general technique for estimating the Hessian of a function based on Stein’s identity in Section 4.1 and use it to provide a zeroth-order cubic regularization method and its analysis in Section 4.2.

Charles Stein, in his seminal paper [Ste72], proposed a method for characterizing Gaussian random variables. Specifically, a random vector, $u\sim N(0,I_{d})$, is standard Gaussian if and only if, ${\bf E}\left[u~{}g(u)\right]={\bf E}\left[\nabla g(u)\right]$, for all absolutely continuous function $g$. Note that Stein’s identity, naturally relate function queries (left hand side of Equation 1.9) to gradients (right hand side of Equation 1.9) and thus is naturally suited for zeroth-order optimization. Indeed the Gaussian smoothing technique proposed by [NS17], is based on the Stein’s identity. Indeed, if we let $g(u)=f(x+\nu u)$ in Equation 1.9, it is easy to see that the identity in Equation 1.3 holds by simply evaluating the Gaussian Stein’s identity in Equation 1.9. Recall that the results in Sections 2 and 3 are essentially based on approximately estimating the gradient information based on the Gaussian smoothing technique [NS17]. In this section, we develop techniques for approximately estimating the Hessian using zeroth-order oracle, based on second-order Stein’s identities. It is worth noting that [Erd16] also use Stein’s identities to estimate the Hessian but they only work in the restricted framework of generalized linear models with Gaussian data. Our use of Stein’s identity to estimate Hessians, is completely different and we provide Hessian estimators for a general class of non-covnex, smooth functions, even for deterministic functions.

The second-order Gaussian Stein’s identity, that we provide here informally for convenience, states thats ${\bf E}[(uu^{\top}-I_{d})~{}g(u)]={\bf E}[\nabla^{2}g(u)]$, for all functions $g$ with well-defined Hessians. Similar to first-order Stein’s identity, this naturally relates function queries to Hessians. In order to leverage this, similar to the previous case, we let $g(u)=f(x+\nu u)$ and note that we have

This provides a way of approximately estimating the Hessian of the function $f_{\nu}$ by approximating the expectation on the left hand side using Gaussian samples. Hence, we can leverage this estimate of Hessian of the smoothed function to get an approximate estimate of Hessian of $f$. Similar to the gradient-free setting, we now have the following estimates of the Hessian.

Let the Hessian estimator be defined in (4.4) and Assumption 6 hold for $F(x,\xi)$. Then, we have

Proof. Noting (4.4) and Holder’s inequality, we have

which together with Assumptions 2, assumption (4.2) for $F(x,\xi)$, and the fact that

Moreover, by Holder’s inequality, we have

Under Assumption 6, denoting the Hessian of $f$ by $H_{f}$ for simplicity, we have

Proof. Taking $y=x+\nu u$ in Equation 4.2, note that we have

which together with (4.9) and (1.5), imply that

Note that (4.8) is obtained only under Assumption 6. However one could obtain an improved bound on the approximation error, by making the more restrictive assumption of interchangeability of differentiation and expectation as follows: $\|H_{f_{\nu}}-H_{f}\|=\|{\bf E}[\nabla^{2}f(x+\nu u)]-\nabla^{2}f(x)\|\leq{\bf E}\|\nabla^{2}f(x+\nu u)-\nabla^{2}f(x)\|\leq L_{H}\nu{\bf E}\|u\|\leq L_{H}\nu\sqrt{d}$. While this provides an improved dependency on $d$, we remark that this improvement does not translate to the improvement in the number of zeroth-order oracle calls, at least for the cubic regularized method as discussed in Section 4.2.

Recall the definition of a spectral function below.

Then, we define the robust Hessian estimator as

where $\kappa>0$ is a tuning parameter. This provides us with a robust Hessian estimator that allows for the function $F$ to have heavy tails. Furthermore, the more standard median-of-means estimator [NY83] provides a robust gradient estimator as well. A thorough treatment of the estimation error of the robust Hessian and gradient follows from an analysis similar to that of [Min18] and [NY83] respectively, although we do not outline the details in the current paper. We also remark that while the spectral truncation argument makes the estimator robust, the computational advantage of the vanilla estimator in Equation 4.4 is lost.

2 Zeroth-Order Stochastic Cubic Regularized Newton Method

Our goal in this subsection is to provide a second-order algorithmic framework using the estimated gradient and Hessian based on Stein’s identities. In particular, we present a zeroth-order stochastic cubic regularized Newton method in Algorithm 7. Note that the output of this algorithm, similar to the other algorithms presented in this paper for nonconvex problems, is a random index from the generated trajectory. In order to analyze its complexity, we first state a result due to [NP06] that provides optimality conditions of the cubic regularized subproblem in step 2 of Algorithm 7.

Our next result is the analogous result of Lemma 2.1 for the averaged Hessian matrices.

Let $\bar{H}_{\nu}^{k}$ be computed by (4.11), $b_{k}\geq 4(1+2\log 2d)$. Then under Assumptions 1 and 2, we have

Proof. First, note that by Theorem 1 in [Tro16], we have

where $\Delta_{k,i}=H_{\nu}(x_{k-1},\xi^{H}_{k,i},u^{H}_{k,i})-\nabla^{2}f_{\nu}(x_{k-1})$ and $C(d)=4(1+2\log 2d)$. Now, noting (4.7), we have

which together with the above inequality and the fact that

Combining this inequality with (4.8), we obtain (4.14). Moreover, by Holder’s inequality we have

Now, by vector-valued Rosenthal’s inequality (see, for example, Theorem 5.2 in [Pin94]) and (4.6), we obtain

which together with the above inequality and (4.16) imply (4.15).

We now proceed to provide the complexity results for Algorithm 7. We first require two intermediate results.

Let $\{x_{k}\}$ be computed by Algorithm 7. Then under Assumptions 1 and 2, we have

where $\delta_{k}^{g},\delta_{k}^{H}>0$ are chosen such that

Proof. By the equality condition in Lemma 4.3 and (4.1), we have

Taking expectation from both sides of the above inequality and noting that $\delta_{k}^{g},\delta_{k}^{H}$ given in (4.18) are well-defined by properly choosing $m_{k}$ and $b_{k}$ in Lemmas 2.1 and 4.4, we obtain

Also, by smoothness assumption of the Hessian and the inequality relation in Lemma 4.3

Taking expectation from both sides of the above inequality and noting definition of $\delta_{k}^{H}$ in (4.18), we obtain

Combining the above inequality with (4.19), we obtain (4.17).

Let $\{x_{k}\}$ be computed by Algorithm 7 for a given iteration limit $N\geq 1$. Then under Assumptions 1 and 2, we have

where $R$ is an integer random variable whose probability distribution $P_{R}(\cdot)$ is supported on $\{1,\ldots,N\}$ and given by

and $\delta_{k}^{g},\delta_{k}^{H}>0$ are defined in (4.18).

Proof. First, note that by (4.2), (4.12), and the fact that $\alpha_{k}\geq L_{H}$, we have

Combining the above two relations, we obtain

where the last inequality follows from the Young’s inequality. Taking expectation from both sides, re-arranging the terms, and noting (4.18), we obtain

Summing up the above inequalities, dividing both sides by $\sum_{k=1}^{N}\alpha_{k}$, and noting (4.21), we obtain (4.20).

Let $\{x_{k}\}$ be computed by Algorithm 7 for a given iteration limit $N\geq 1$. Moreover, assume that the parameters are set to

where $R$ is uniformly distributed over $\{1,\ldots,N\}$. As a consequence, to obtain an $\epsilon$ second-order stationary point of the problem, the total number of samples required to compute the gradient and Hessian are, respectively, bounded by

Proof. First, note that by (4.22), Lemmas 2.1, and 4.4, we can ensure that (4.18) is satisfied by $\delta_{k}^{g}=2\epsilon/5$ and $\delta_{k}^{H}=\epsilon/138$. Moreover, by Lemma 4.6, we have

Hence, by choosing $N$ according (4.22), and noting Lemma 4.6, we obtain (4.23). Therefore, $x_{R}$ is an $4\epsilon$ second-order stationary point of the problem. Finally, note that the total number of required samples to obtain such a solution is bounded by

Discussion

In this work, we propose and analyze zeroth-order stochastic approximation algorithms for convex and nonconvex problems motivated by modern machine learning challenges. Specifically, we provide zeroth-order algorithms to deal with constraints, dimensionality and saddle-points in nonconvex stochastic optimization problems. While our focus was on general stochastic optimization problems, one could naturally obtain better rates in the case of finite-sum optimization problems with various variance reduction techniques. Several concrete extensions are possible for future work. The performance of conditional gradient algorithm in the high-dimensional constrained optimization setting is not well-explored; the interaction between the geometry of the constraint set, sparsity structure and zeroth-order information is extremely interesting to explore. Obtaining regret bounds for the non-convex problems considered in this work is more challenging. Furthermore, lower bounds can be explored for the cases considered in this paper when $f$ is nonconvex. Finally, obtaining second-order stationarity results in the constrained setting is more challenging. We plan to extend our results for these setting in the future.

References