A Stochastic Derivative-Free Optimization Method with Importance Sampling: Theory and Learning to Control

TL;DR

This paper introduces a novel derivative-free optimization method that incorporates importance sampling, providing theoretical complexity improvements and demonstrating practical effectiveness in high-dimensional control tasks.

Contribution

It is the first to combine importance sampling with derivative-free optimization and offers new theoretical complexity bounds for non-convex, convex, and strongly convex functions.

Findings

The method achieves improved complexity bounds.

Experiments confirm theoretical predictions on synthetic and real datasets.

Significant sample complexity reduction in high-dimensional control problems.

Abstract

We consider the problem of unconstrained minimization of a smooth objective function in $\R^n$ in a setting where only function evaluations are possible. While importance sampling is one of the most popular techniques used by machine learning practitioners to accelerate the convergence of their models when applicable, there is not much existing theory for this acceleration in the derivative-free setting. In this paper, we propose the first derivative free optimization method with importance sampling and derive new improved complexity results on non-convex, convex and strongly convex functions. We conduct extensive experiments on various synthetic and real LIBSVM datasets confirming our theoretical results. We further test our method on a collection of continuous control tasks on MuJoCo environments with varying difficulty. Experiments suggest that our algorithm is practical for high dimensional continuous control problems where importance sampling results in a significant sample complexity improvement.

Figures

Tables

Table 1. Table 1: Summary of the new derived complexity results as opposed to uniform sampling where r 0 = f ​ ( x 0 ) − f ​ ( x ∗ ) subscript 𝑟 0 𝑓 subscript 𝑥 0 𝑓 subscript 𝑥 r_{0}=f(x_{0})-f(x_{*}) . The assumptions listed are in addition to Assumption 1 . R 0 < ∞ subscript 𝑅 0 R_{0}<\infty indicates a bounded level set where the exact definition is given in Assumption 2 . The key differences in complexity between the uniform and importance sampling are detailed in text.

Assumptions on $f$	Uniform Sampling Complexity	Importance Sampling	Importance Sampling Complexity [NEW]	Theorem
None	$\frac{4\sqrt{2}{r}_0{n}^{2}L}{{ϵ}^2}$	$p_i=\frac{\sqrt{L_i}}{{\sum }_{i=1}^n\sqrt{L_i}}$	$\frac{4\sqrt{2}{r}_0{\left({\sum }_{i=1}^n\sqrt{L_i}\right)}^2}{{ϵ}^2}$	1
None	$\frac{4\sqrt{2}{r}_0{n}^{2}L}{{ϵ}^2}$	$p_i=\frac{L_i}{{\sum }_{i=1}^n{L}_i}$	$\frac{4\sqrt{2}{r}_{0}n\left({\sum }_{i=1}^n{L}_i\right)}{{ϵ}^2}$	2
Convex,	$8R_0^2{n}^{2}L\left(\frac{1}{ϵ}-\frac{1}{r_0}\right)$	$p_i=\frac{L_i}{{\sum }_{i=1}^n{L}_i}$	$8R_0^{2}n{\sum }_{i=1}^n{L}_i\left(\frac{1}{ϵ}-\frac{1}{r_0}\right)$	3
$\lambda$-strongly convex	$\frac{nL}{\lambda }\log \left(\frac{r_0}{ϵ}\right)$	$p_i=\frac{L_i}{{\sum }_{i=1}^n{L}_i}$	$\frac{{\sum }_{i=1}^n{L}_i}{\lambda }\log \left(\frac{r_0}{ϵ}\right)$	4

Table 2. Table 2: For each MuJoCo task, we report the average number of episodes required to achieve a predefined reward threshold. Results for our method is averaged over five random seeds, the rest is copied from ( ? ) (N/A means the method failed to reach the threshold. UNK means the results is unknown since they are not reported in the literature.)

		Fixed Step Size		Adaptive Step Size
	Threshold	STP	STP${}_{\text{IS}}$	STP	STP${}_{\text{IS}}$	ARS(V1-t)	ARS(V2-t)	NG-lin	TRPO-nn
Swimmer-v1	325	320	110	200	90	100	427	1450	N/A
Hopper-v1	3120	3970	2400	3720	1870	51840	1973	13920	10000
HalfCheetah-v1	3430	13760	4420	5040	2710	8106	1707	11250	4250
Ant-v1	3580	107220	43860	96980	26480	58133	20800	39240	73500
Humanoid-v1	6000	N/A	530200	N/A	296800	N/A	142600	130000	UNK

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x),$$

$$f(x_{\ast })\le f(x+t{e}_i)\le f(x)+{\nabla }_{i}f(x)t+\frac{L_i}{2}{t}^2$$

$$x_{k+1}=\arg \min \{f(x_k),f(x_{+}),f(x_{-})\}.$$

$$\underset{y\in {\mathbb{R}}^n}{\min }{f}_B(y)\stackrel{\text{def}}{=}f(By)$$

$$K\ge \frac{2{\left(\frac{\sqrt{2}(f(x_0)-f(x_{\ast }))}{{\alpha }_0}+\frac{{\alpha }_0{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2}\right)}^2}{{\left({\min }_i\frac{p_i}{v_i}\right)}^2{ϵ}^2},$$

$$\frac{4\sqrt{2}\left(f(x_0)-f(x_{\ast })\right){\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{{\left({\min }_i\frac{p_i}{v_i}\right)}^2{ϵ}^2}.$$

$$\frac{4\sqrt{2}(f(x_0)-f(x_{\ast })){\left({\sum }_{i=1}^n\sqrt{L_i}\right)}^2}{{ϵ}^2}.$$

$$\frac{4\sqrt{2}(f(x_0)-f(x_{\ast }))n\left({\sum }_{i=1}^n{L}_i\right)}{{ϵ}^2}.$$

$$\frac{4\sqrt{2}(f(x_0)-f(x_{\ast }))n^{2}L}{{ϵ}^2}.$$

$$K\ge \frac{2n(f(x_0)-f(x_{\ast }))}{\left({\min }_i\frac{p_i}{v_i}\right)\left(1-\frac{{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2n\left({\min }_i\frac{p_i}{v_i}\right)}\right){ϵ}^2},$$

$$\frac{4\left(f(x_0)-f(x_{\ast })\right)n\left({\sum }_{i=1}^n{L}_i\right)}{{ϵ}^2}.$$

$$R_0\stackrel{\text{def}}{=}\underset{x}{\max }\{\Vert x-x_{\ast }{\Vert }_{\infty }:f(x)\le f(x_0)\}<\infty$$

$$f(x)-f(x_{\ast })$$

$$\le \Vert \nabla f(x){\Vert }_1\Vert x-x_{\ast }{\Vert }_{\infty }\le R_0\Vert \nabla f(x){\Vert }_1.$$

$$k\ge \frac{8R_0^{2}n}{{\min }_i\frac{p_i}{v_i}}\left(\frac{1}{ϵ}-\frac{1}{r_0}\right),$$

$$k\ge 8R_0^{2}n\left(\sum _{i=1}^n{L}_i\right)\left(\frac{1}{ϵ}-\frac{1}{r_0}\right).$$

$$f(x)\ge f(y)+⟨\nabla f(y),x-y⟩+\frac{\lambda }{2}\Vert x-y{\Vert }_2^2.$$

$$k\ge \frac{{\max }_i\frac{v_i}{p_i}}{\lambda }\log \left(\frac{2\left(f(x_0)-f(x_{\ast })\right)}{ϵ}\right),$$

$$\mathbb{E}\left[f(x_{k+1})|x_k\right]\le f(x_k)-\mathbb{E}\left[{\alpha }_{i_k}|⟨\nabla f(x_k),e_{i_k}⟩||x_k\right]+\frac{1}{2}\mathbb{E}\left[{\alpha }_{i_k}^2{L}_{i_k}|x_k\right].$$

$$f(x_{k+1})$$

$$K\ge \frac{2{\left(\frac{\sqrt{2}(f(x_0)-f_{\ast })}{{\alpha }_0}+\frac{{\alpha }_0{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2}\right)}^2}{{\left({\min }_i\frac{p_i}{v_i}\right)}^2{ϵ}^2},$$

$$\mathbb{E}\left[f(x_{k+1})|x_k\right]\le f(x_k)-\underset{\text{¬}}{\underbrace{\mathbb{E}\left[{\alpha }_{i_k}|⟨\nabla f(x_k),e_{i_k}⟩||x_k\right]}}+\underset{\text{­}}{\underbrace{\frac{1}{2}\mathbb{E}\left[{\alpha }_{i_k}^2{L}_{i_k}|x_k\right]}},$$

$$\mathbb{E}\left[{\alpha }_{i_k}|⟨\nabla f(x_k),e_{i_k}⟩||x_k\right]$$

$$\frac{1}{2}\mathbb{E}\left[{\alpha }_{i_k}^2{L}_{i_k}|x_k\right]$$

$${\theta }_{k+1}\le {\theta }_k-\frac{{\alpha }_0{\min }_i\frac{p_i}{v_i}}{\sqrt{k+1}}{g}_k+\frac{{\alpha }_0^2{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2(k+1)},$$

$$g_k\le \frac{1}{{\min }_i\frac{p_i}{v_i}}\left(\frac{({\theta }_k-{\theta }_{k+1})\sqrt{k+1}}{{\alpha }_0}+\frac{{\alpha }_0{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2\sqrt{k+1}}\right).$$

$$\sum _{j=l}^{2l}{\theta }_j-{\theta }_{j+1}={\theta }_l-{\theta }_{2l+1}\le f(x_0)-f(x_{\ast })\stackrel{\text{def}}{=}C,$$

$$g_j$$

$$K\ge \frac{2n\left(f(x_0)-f(x_{\ast })\right)}{{\min }_i\frac{p_i}{v_i}\left(1-\frac{{\sum }_{i=1}^n\frac{p_i{L}_i}{v_i^2}}{2n{\min }_i\frac{p_i}{v_i}}\right){ϵ}^2},$$

$$\mathbb{E}\left[f(x_{k+1})|x_k\right]$$

$$\mathbb{E}\left[{\alpha }_{i_k}|⟨\nabla f(x_k),e_{i_k}⟩||x_k\right]$$

Full text

Introduction

In this paper, we consider the optimization problem

[TABLE]

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a “smooth” but not necessarily convex function, bounded from below and it achieves its global minimum at some $x_{*}\in\mathbb{R}^{n}$. In particular, we enforce throughout the paper the following smoothness assumption:

Assumption 1**.**

The objective function $f$ has coordinate-wise Lipschitz gradient, with Lispchitz constants $L_{1},\dots,L_{n}>0$. Moreover, $f$ is bounded from below by $f(x_{*})\in\mathbb{R}$. That is, $f$ satisfies

[TABLE]

for all $x\in\mathbb{R}^{n}$ and $t\in\mathbb{R}$, where $\nabla_{i}f(x)$ is the $i$th partial derivative of $f$ at $x$.

DFO. We consider the Derivative-Free Optimization (DFO) (?; ?) setting. That is, we assume that the derivatives of $f$ are numerically impractical to obtain, unreliable (e.g., noisy function evaluations (?)), or not available at all. In typical DFO applications, evaluations of $f$ are possible through runs/simulations of some black-box software only. Optimization problems of this type appear in many applications, including computational medicine (?), fluid-dynamics (?; ?), localization (?; ?) and continuous control (?; ?).

Literature on DFO methods for solving (1) has a long history. Some of the first approaches were based on deterministic direct search (DDS) (?). Subsequently, additional variants of DDS, including randomized approaches, were proposed in (?; ?; ?; ?; ?; ?). However, complexity bounds for deterministic direct search methods have only been established recently by the works of (?; ?; ?; ?). Recently, complexity bounds have also been derived for randomized methods (?; ?; ?). For instance, the work of (?; ?) imposes a decrease condition on whether to accept or reject a step of a set of random directions. Moreover, (?; ?) derived new complexity bounds for accelerated random search.

More recently, Bergou et. al. proposed a new randomized direct search method called Stochastic Three Points (STP) method. STP, in each iteration $k$, generates a random search direction $s_{k}$ according to a certain probability law, then compares the objective function at three points: current iterate $x_{k}$, a point $x_{+}=x_{k}+\alpha_{k}s_{k}$ in the direction of $s_{k}$ and a point $x_{-}=x_{k}-\alpha_{k}s_{k}$ in the direction of $-s_{k}$. The method then chooses the best of these three points as the new iterate:

[TABLE]

**Notation: ** As for the notations, $\mathbb{E}\left[\cdot\right]$ denotes the expectation operator. The standard inner product is defined as $\langle x,y\rangle=x^{\top}y$. We also denote the $\ell_{1}$-norm and $\ell_{2}$-norm by $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$, respectively. We define $L=\max_{i}L_{i}$ for a given sequence of scalars $L_{1},\ldots,L_{n}$.

Paper Overview and Contributions

While importance sampling, a term that typically refers to the nonuniform sampling of random directions in stochastic algorithms, has been widely investigated in gradient based methods (?; ?; ?; ?), to the best of our knowledge there exists no work on importance sampling in the random direct search setting. To this end, we study STP and analyze its complexity with arbitrary probabilities. In particular, we restrict the random directions to be sampled from discrete distributions, i.e., in each iteration of STP a random direction $s_{k}$ from a finite set of independent directions $\{b_{1},\dots,b_{n}\}\subset\mathbb{R}^{n}$ is sampled. That is, we set $s_{k}=b_{i}$ with probability $p_{i}>0$. We then propose new sampling strategies that are either optimal or at least improve the complexity bounds, i.e., importance sampling.

Coordinate directions

Without loss of generality, we only consider directions in the canonical basis of $\mathbb{R}^{n}$, i.e., $e_{1},\dots,e_{n}$. The general case can be recovered via a linear change of variables: $x=By$, where $B\in\mathbb{R}^{n\times n}$. Indeed, consider the problem

[TABLE]

instead. A coordinate update $y_{k+1}=y_{k}+\alpha_{k}e_{i}$ for the re-parameterized problem (2) corresponds to updates of the form $x_{k+1}=x_{k}+\alpha_{k}b_{i}$, where $b_{i}$ is the $i$th column of $B$, for the original problem (1). In light of the above discussion, the newly proposed algorithm dubbed STP${}_{\text{{IS}}}$ is formally described as Algorithm 1.

Complexity bounds

To the best of our knowledge, ours are the first complexity bounds (bounds on the number of iterations) for a DFO method with importance sampling. We design importance sampling that improves the worst-case iteration complexity bounds compared to state-of-the-art algorithms. These bounds have the same dependence on the precision $\epsilon$ as classical bounds in the literature, i.e. $1/\epsilon^{2}$ for non-convex $f$, $1/\epsilon$ for convex $f$ and $\log(1/\epsilon)$ for strongly convex $f$; see for instance (?; ?). However, the leading constant, which is often the bottleneck in practical performance, especially when low or medium accuracy solutions are acceptable, is improved and often dramatically so. Typically, the improvement is via replacing the maximum Lipschitz constant of the gradient by the average Lipschitz constants of all coordinates (see Theorems 1, 2, 3, and 4). The improvement we obtain is similar to the improvement obtained by importance sampling in stochastic coordinate (gradient) descent methods  (?; ?; ?). Table 1 summarizes complexity results obtained in this paper for STP and for STP${}_{\text{{IS}}}$. The assumptions in Table 1 are in addition to Assumption 1.

Empirical results

In addition to our theoretical analysis, we conduct extensive testing to show the efficiency of the proposed method in practice. We use both synthetic and real111We use several LIBSVM datasets (?). datasets for ridge regression and squared SVM problems. In the non-convex case, we use continuous control tasks from the MuJoCo (?) suite following the recent success of DFO compared to model-free RL (?; ?). Results show that our approach leads to huge speedups compared against uniform sampling, the improvement can reach several orders of magnitude and comparable or better than state-of-art policy gradient methods.

Non-Convex Case

This section describes our complexity results for Algorithm 1 in the case when $f$ is allowed to be non-convex. We show that this method guarantees complexity bounds with the same order in $\epsilon$ as classical bounds in the literature, i.e., $1/\epsilon^{2}$ with an improved dependence on the Lipschitz constant. All proofs are left for the appendix.

Theorem 1**.**

Let Assumption 1 be satisfied. Choose $\alpha_{i_{k}}=\tfrac{\alpha_{0}}{v_{i_{k}}\sqrt{k+1}}$ where $\alpha_{0}>0$. If

[TABLE]

then $\displaystyle\min_{k=0,1,\dots,K}\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]\leq\epsilon.$

Note that the complexity depends on $\alpha_{0}$. The optimal choice of $\alpha_{0}$ minimizing (3) is $\alpha_{0}^{*}=8^{1/4}\sqrt{\frac{{f(x_{0})-f(x_{*})}}{\sum_{i=1}^{n}\frac{p_{i}{L_{i}}}{{v_{i}^{2}}}}},$ in which case the complexity bound (3) takes the form

[TABLE]

Importance sampling. The complexity depends on the choice of the probabilities $\left\{p_{i}\right\}_{i=1}^{n}$ and the quantities $\left\{v_{i}\right\}_{i=1}^{n}$. For instance, if $p_{i}=\frac{\sqrt{L_{i}}}{\sum_{i=1}^{n}\sqrt{L_{i}}}$ and $v_{i}=\sqrt{L_{i}}$, then the complexity becomes

[TABLE]

On the other hand, if $p_{i}=\frac{{L_{i}}}{\sum_{i=1}^{n}{L_{i}}}$ and $v_{i}={L_{i}}$ then the complexity becomes

[TABLE]

Under the choice of uniform sampling, i.e. $p_{i}=\frac{1}{n}$ and the choice $v_{i}=L$, we recover the uniform sampling complexity of (?)

[TABLE]

Note that $\left(\sum_{i=1}^{n}\sqrt{L_{i}}\right)^{2}\leq n^{2}L$ and $n\left(\sum_{i=1}^{n}{L_{i}}\right)\leq n^{2}L$. Therefore, complexity bounds in the number of iterations is improved with the proposed importance sampling strategies (4) and (5). We now state a complexity theorem for Algorithm 1 when using non-decreasing stepsizes.

Theorem 2**.**

Let Assumption 1 be satisfied. Choose $\alpha_{i_{k}}=\frac{\epsilon}{nv_{i_{k}}}$ where $\sum_{i=1}^{n}\frac{p_{i}L_{i}}{v_{i}^{2}}<2n\left(\min_{i}\frac{p_{i}}{v_{i}}\right)$. If

[TABLE]

then $\displaystyle\min_{k=0,1,\dots,K}\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]\leq\epsilon.$

Under the choice of importance sampling $p_{i}=\frac{L_{i}}{\sum_{i=1}^{n}{L_{i}}}$ and $v_{i}={L_{i}}$, the complexity (6) becomes

[TABLE]

Similar to Theorem 1, the uniform sampling complexity of (?) can be recovered with $p_{i}=\frac{1}{n}$ and $v_{i}=L$. Note that the uniform sampling complexity is proportional to $n^{2}L$ and since $n\sum_{i=1}^{n}{L_{i}}\leq n^{2}L$, the worst case complexity of the number of iterations is also improved for this choice of importance sampling.

Convex Case

This section describes our complexity results for Algorithm 1 when the objective function $f$ is convex. We show that this method guarantees complexity bounds with the same order in $\epsilon$ as classical bounds in the literature, i.e., $1/\epsilon$ with an improved dependence on the Lipschitz constant. We will need the following additional assumption in the sequel.

Assumption 2**.**

The function $f$ is convex and has a bounded level set at $x_{0}$. That is, f satisfies:

[TABLE]

Note that if $f$ is convex and has bounded level sets, the following holds:

[TABLE]

Theorem 3**.**

Let Assumptions 1 and 2 be satisfied. Choose $\alpha_{i_{k}}=\frac{|f(x_{k}+te_{i_{k}})-f(x_{k})|}{tv_{i_{k}}}$ and sufficiently small $t$ (see the appendix for the bound on $t$). If

[TABLE]

then $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$.

Here, the complexity bound in (9) depends on the choice of the probabilities $\left\{p_{i}\right\}_{i=1}^{n}$ and the quantities $\left\{v_{i}\right\}_{i=1}^{n}$. Note that if $v_{i}=L_{i}$, it is easy to show that the minimum of the complexity bound (9) in $p_{i}$, i.e. minimizing $1/\min_{i}\frac{p_{i}}{v_{i}}$ in $p_{i}$, over a probability simplex is attained at $p_{i}=\frac{L_{i}}{\sum_{i=1}^{n}L_{i}}$. Thus, the complexity bound of this this importance sampling is

[TABLE]

Since uniform sampling is proportional to $n^{2}L$ and that $n\sum_{i=1}^{n}L_{i}\leq n^{2}L$, importance sampling is clearly better than uniform sampling.

Strongly Convex Case

This section describes our complexity results for Algorithm 1 when objective function $f$ is $\lambda$-strongly convex. We show that this method guarantees complexity bounds with the same order in $\epsilon$ as classical bounds in the literature, i.e., $\log(1/\epsilon)$ with an improved dependence on the Lipschitz constant. First, we define $\lambda$-strongly convexity functions.

Assumption 3**.**

The function $f$ is $\lambda$-strongly convex. That is, for some $\lambda>0$, the following holds

[TABLE]

Theorem 4**.**

Let Assumptions 1 and 3 be satisfied. Choose $\alpha_{i_{k}}=\frac{|f(x_{k}+te_{i_{k}})-f(x_{k})|}{tv_{i_{k}}}$ and a sufficiently small $t$ (see the appendix for the bound on $t$). If

[TABLE]

then $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$.

The complexity bound in (10) is minimized, in $p_{i}$ for $v_{i}=L_{i}$, with $p_{i}=\frac{L_{i}}{\sum_{i=1}^{n}L_{i}}$. Importance sampling improves over uniform sampling, since $\sum_{i=1}^{n}L_{i}\leq nL$.

Experiments

We conduct extensive experiments on synthetic and real datasets comparing the uniform sampling STP against the importance sampling version STP${}_{\text{{IS}}}$. The experiments are conducted on several choices of the function $f$. In particular, we perform experiments on regularized ridge regression on synthetic data and squared SVM loss on real data from the LIBSVM dataset (?). Moreover, for non-convex problems, we compare STP and STP${}_{\text{{IS}}}$ on various continuous control environments on MuJoCo (?). We also compare against state-of-art solvers for the continuous control task.

Ridge regression on synthetic data

We compare STP against STP${}_{\text{{IS}}}$ on synthetic data on the regularized ridge regression problem: $f(x)=\frac{1}{2m}\|Ax-y\|_{2}^{2}+\frac{\lambda}{2}\|x\|_{2}^{2},$ where $A\in\mathbb{R}^{m\times n}$, $y\in\mathbb{R}^{m}$ are the data and $\lambda>0$ is the regularization parameter. The elements of $A$ and $y$ were sampled from the standard Gaussian distribution $\mathcal{N}\left(0,1\right)$. Note that for ridge regression, $L_{i}=\frac{1}{m}\|A(:,i)\|_{2}^{2}+\lambda$ where, following (?), we normalize data such that $\|A(:,1)\|_{2}=1$ and $\|A(:,i)\|_{2}=\frac{1}{m}$, $i=2,\ldots,m$ and set $\lambda=\frac{1}{m}$. We compute a high accuracy solution $x_{*}$ by solving the ridge regression problem exactly with a linear solver. Thereafter, the metric used is the difference between the current objective value and the optimal one, i.e. $f(x)-f(x_{*})$. Since the objective is $\lambda$-strongly convex, we use the stepsize suggested by Theorem 4. In all experiments, we set $t=10^{-4}$. We perform experiments across difference choices of $m$ and $n$. In the first row of Figure 1, we compare both methods with a fixed $n=10$ and a varying $m$, i.e. $m\in\{10^{3},10^{4},10^{5}\}$. The superior performance of $\texttt{STP}_{\texttt{IS}}$ over STP is evident from Figure 1. Moreover, we conduct further experiments where $m$ is fixed such that $m=100$ but with a varying dimension, i.e. $n\in\{10^{1},10^{2},10^{3}\}$. All experiments are conducted 10 times and we report the average, worst and best performances. A similar behaviour is also present as seen in the second row of Figure 1 where $\texttt{STP}_{\texttt{IS}}$ is far more superior to STP. In all experiments, the stopping criterion is set such that both STP and $\texttt{STP}_{\texttt{IS}}$ run for exactly $5\times 10^{2}$ iterations for small problems, i.e. $n=10$, while for problems of size $n=10^{2}$ and $n=10^{3}$, both methods are terminated at $5\times 10^{3}$ and $15\times 10^{3}$ iterations, respectively.

Ridge regression and squared SVM on real data

We also conduct experiments on the regularized ridge regression problem on real datasets where $A$ and $y$ are from LIBSVM data. We follow the same protocol as the experiments on synthetic data. We compare both algorithms on 6 different datasets, namely, australian, mushrooms, a9a, heart, cov1 and ijcnn1. In addition to ridge regression, we conduct experiments on the same real datasets on the regularized squared SVM loss: $f(x)=\frac{1}{2}\sum_{i=1}^{m}\max\left(0,1-y_{i}a_{i}^{\top}x\right)^{2}+\frac{\lambda}{2}\|x\|_{2}^{2},$ where $a_{i}$ is the $i^{\text{th}}$ row of $A$. Note that $L_{i}=\|A(:,i)\|_{2}^{2}+\lambda$. Since the squared SVM problem does not exhibit a closed form solution, we compare both STP and $\texttt{STP}_{\texttt{IS}}$ in terms of the objective value $f(x)$. In Figure 2, we show the comparison between both STP and $\texttt{STP}_{\texttt{IS}}$ on the ridge regression on all 6 datasets. It is clear that using the proposed importance sampling is far more superior to standard uniform sampling. The improvement is also consistently present on the squared SVM problem as seen in Figure 3.

Continuous control experiments

Here, we address the problem of model-free control of a dynamical system. Model-free reinforcement learning algorithms (especially policy gradient methods), provide an off-the-shelf model-free approach to learn how to control a dynamical system. Such models have been typically benchmarked in a simulator. Thus, we adopt the MuJoCo (?) continuous control suite following its wide adaptation. We choose 5 problems with various difficulty Swimmer-v1, Hopper-v1, HalfCheetah-v1, Ant-v1, and Humanoid-v1. In all experiments, we use linear policies similar to (?; ?).

Considering the stochastic nature of the dynamical systems, i.e. $f$ is stochastic, we take multiple ($K$) measurements for $f(x_{k}),f(x_{+})$ and $f(x_{-})$ and use their mean as the function values. Considering the varying dimensionality of the state space, we use different $K$ for each problem, in particular, we set $K=2$ for Swimmer-v1, $K=4$ for Hopper-v1 and HalfCheetah-v1, $K=40$ for Ant-v1 and $K=120$ for Humanoid-v1. These values are decided using grid search over the set of $K\in\{1,2,4,8,16\}$ for low dimensional problems and $K\in\{20,40,120,240\}$ for high dimensional Ant-v1 and Humanoid-v1 problems. Following our remark given by Equation 2, we use a square matrix $B$ sampled from a standard Gaussian distribution $\mathcal{N}(0,1)$. In our experiments, this coordinate transform resulted in a better performance. Since the Lipschitz constants are not available for continuous control, we learn an estimate of the function using a parametric family (specifically multi-layer perceptron) and use its Lipschitz constants as the estimates to decide importance sampling weights. This is similar to actor-critic methods(?) used in the policy gradient literature. Similar to us, actor-critic methods learn an estimate of the value function and use it to decide which point to evaluate. We defer the details of estimation procedure to the appendix. Following the common practice, we perform all experiments with 5 random initialization and measure the mean average reward at each iteration. We give detailed comparison of STP and our proposed importance sampling variant STP${}_{\text{{IS}}}$ in terms of reward vs. sample complexity in Figure 4 for both adaptive and fixed step size cases (see Theorems 1 and 2). Shaded regions in figures show standard deviations.

As seen from Figure 4, our proposed importance sampling version STP${}_{\text{{IS}}}$ significantly improves sample complexity when compared to STP. Moreover, the difference is significant for high dimensional problems like HalfCheetah, Ant and Humanoid. The results suggest that STP fails to scale to very high dimensional problems like Humanoid. Our method tackles this and improves the sample complexity of STP. Such results also suggest that it is feasible to estimate the coordinate-wise Lipschitz gradient constants, detailed in Assumption 1, of a complicated non-convex function using a data-driven approach. An interesting conclusion from Figure 4 is that adaptive step size performs better than fixed step size even after a large hyper-parameter search, particularly, for higher dimensional problems.

In order to compare our method with the existing state-of-art DFO and policy gradient methods, we also tabulate the sample complexity of our method and several existing baselines. Similar to (?), we compute the average number of episodes needed to reach a predefined threshold. Although there are many DFO and policy gradient methods in literature, we report ARS (?) as a representative DFO method since it outperforms other baselines. As for policy gradient approaches, we report TRPO (?) as a representative policy gradient method since it is widely used in the community. Moreover, we use NG (?) as a policy gradient method using linear policies. We list sample complexity results for all methods and tasks in Table 2.

As the results suggest, STP is competitive with existing solutions for low dimensional problems (Swimmer, Hopper and HalfCheetah) whereas it underperforms existing solutions for Ant and fails to solve the Humanoid problem. Our theoretically proposed importance sampling version STP${}_{\text{{IS}}}$ significantly improves STP and results in a performance either competitive with or better than existing baselines in all problems except for Humanoid. Although our method successfully solves the Humanoid problem, it has worse sample complexity than other solutions. Hence, scaling STP to very large dimensional continuous control problems (e.g. Humanoid-v1 state space has more than 1000 dimensions) is still an open problem. Moreover, for lower dimensional problems like (Swimmer and Hopper), our method outperforms all existing methods.

An interesting question is whether we can use first order methods utilizing the estimated function. In order to estimate the Lipschitz smoothness constant of the function, we utilize a parametric family and one can argue that its gradients can be used as a surrogate gradient in first order optimization. We compare our method with a simple first order baseline using the surrogate gradient in SGD; it fails in all environments. Hence, we do not show the results in the paper.

Conclusion

We propose and analyze a DFO algorithm with importance sampling STP${}_{\text{{IS}}}$ enjoying the best known complexity bounds known in DFO literature. Experiments on ridge regression and squared SVM objectives for both synthetic and LIBSVM datasets demonstrate the superiority of STP${}_{\text{{IS}}}$ over its uniform version. We also conduct experiments on a collection of continuous control tasks on several MuJoCO environments. We are orders of magnitudes better than the uniform sampling and comparable or better than the state-of-art methods.

Acknowledgments. This work was supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research.

Appendix A Preliminaries

We establish the key lemma which will be used to prove the theorems stated in the paper.

Lemma 1**.**

If $f$ satisfies Assumption 1 and following the STP${}_{\text{{IS}}}$ update, the following holds:

[TABLE]

Proof.

Since

[TABLE]

Then the result follows by taking conditional expectation on $x_{k}$. ∎

Appendix B Non convex case

Theorem 1**.**

Let Assumption 1 hold. Choose $\alpha_{i_{k}}=\tfrac{\alpha_{0}}{v_{i_{k}}\sqrt{k+1}}$, where $\alpha_{0}>0$. If

[TABLE]

then $\displaystyle\min_{k=0,1,\dots,K}\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]\leq\epsilon.$

Proof.

We have from Lemma 1

[TABLE]

From ① we have

[TABLE]

From ② we have

[TABLE]

By injecting ① and ② in (12) and taking the expectation we get

[TABLE]

where $\theta_{k}=\mathbb{E}\left[f(x_{k})\right]$ and $g_{k}=\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]$. By re-arranging the terms we get:

[TABLE]

We have that the sequence $\left\{f(x_{k})\right\}_{k\geq 0}$ is monotonically decreasing and $f$ is bounded from below by $f(x_{*})$, hence $f(x_{*})\leq\theta_{k+1}\leq\theta_{k}\leq f(x_{0})$ for all $k$. Letting $l=\lfloor K/2\rfloor$, this implies that

[TABLE]

from which we conclude that there must exist $j\in\{l,\dots,2l\}$ such that $\theta_{j}-\theta_{j+1}\leq C/(l+1)$. This implies that

[TABLE]

∎

Theorem 2**.**

Let Assumption 1 be satisfied. Choose $\alpha_{i_{k}}=\frac{\epsilon}{nv_{{i_{k}}}}$ where $\sum_{i=1}^{n}\frac{p_{i}L_{i}}{v_{i}^{2}}<2n\left(\min_{i}\frac{p_{i}}{v_{i}}\right)$. If

[TABLE]

then $\displaystyle\min_{k=0,1,\dots,K}\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]\leq\epsilon$.

Proof.

From equation (12) we have

[TABLE]

From ① we have

[TABLE]

From ② we have

[TABLE]

By injecting ① and ② in (16) and taking the expectation we get

[TABLE]

where $\theta_{k}=\mathbb{E}\left[f(x_{k})\right]$ and $g_{k}=\mathbb{E}\left[\|\nabla f(x_{k})\|_{1}\right]$. By re-arranging the terms we get:

[TABLE]

We have that the sequence $\left\{f(x_{k})\right\}_{k\geq 0}$ is monotonically decreasing and $f(x)$ is bounded from below by $f(x_{*})$, hence $f(x_{*})\leq\theta_{k+1}\leq\theta_{k}\leq f(x_{0})$ for all $k$. Letting $l=\lfloor K/2\rfloor$, this implies that

[TABLE]

from which we conclude that there must exist $j\in\{l,\dots,2l\}$ such that $\theta_{j}-\theta_{j+1}\leq C/(l+1)$. This implies that

[TABLE]

∎

Appendix C Convex case

We state a lemma which will be useful latter on in the analysis

Lemma 2**.**

Let assumption 1 be satisfied. Let $\alpha_{i_{k}}=\frac{|f(x_{k}+te_{i_{k}})-f(x_{k})|}{tv_{i_{k}}}$, then the following inequality holds:

[TABLE]

Proof.

It follows directly by noting that

[TABLE]

which follows in a similar fashion to Theorem [13] in (?). ∎

Theorem 3**.**

Let Assumptions 1 and 2 be satisfied. Choose $\alpha_{i_{k}}=\frac{\alpha_{0}\left(f(x_{k})-f(x_{*})\right)}{v_{i_{k}}}$, where $0<\alpha_{0}<\frac{2\min_{i}\frac{p_{i}}{v_{i}}}{R_{0}\sum_{i=1}^{n}\frac{p_{i}L_{i}}{v_{i}^{2}}}$. If

[TABLE]

then $\mathbb{E}\left[f(x_{k})\right]-f(x_{*})\leq\epsilon$.

Proof.

We have from Lemma 1

[TABLE]

From ① we have

[TABLE]

From ② we have

[TABLE]

Let $r_{k}=\mathbb{E}\left[f(x_{k})\right]-f(x_{*})$. By substituting ① and ② in (21) and then taking the expectation we get

[TABLE]

where $a=\frac{\alpha_{0}}{R_{0}}\min_{i}\frac{p_{i}}{v_{i}}-\frac{\alpha_{0}^{2}}{2}\sum_{i=1}^{n}\frac{p_{i}L_{i}}{v_{i}^{2}}.$ Therefore,

[TABLE]

Therefore, we have $\frac{1}{r_{k}}\geq\frac{1}{r_{0}}+ka$ and hence $r_{k}\leq\frac{1}{\frac{1}{r_{0}}+ka}$. Therefore, if $k\geq\frac{1}{a}\left(\frac{1}{\epsilon}-\frac{1}{r_{0}}\right),$ then $r_{k}\leq\frac{1}{\frac{1}{r_{0}}+ka}\leq\epsilon$. ∎

Note that the stepsizes in the previous theorem depend on the optimal value $f(x_{*})$. In practice, we cannot always use these stepsizes as we usually do not know $f(x_{*})$. Next theorem will suggest stepsizes that are independent from $f(x_{*})$ for which we get an optimized complexity as well.

Theorem 4**.**

Let Assumptions 1 and 2 be satisfied. Choose $\alpha_{i_{k}}=\frac{|f({x}_{k}+te{i_{k}})-f({x}_{k})|}{tv_{i_{k}}}$, then for a small t that satisfies

[TABLE]

where $r_{0}=f(x_{0})-f(x_{*}),$ and if

[TABLE]

then $r_{k}=\mathbb{E}[f({x}_{k})]-f(x_{*})\leq\epsilon$.

Proof.

Taking the expectation of (19) on $x_{k}$ with the choice of $\alpha_{i_{k}}$ we have:

[TABLE]

As for ① , taking the total expectation we have

[TABLE]

where $r_{0}=\mathbb{E}\left[f(x_{0})\right]-f(x_{*})$. Note that the first inequality follows by Jensen’s inequality. The last inequality follows from the fact that $f$ has Lipschitz gradient, with Lispchitz constant $nL$ (this is a direct property from Assumption 1).

As for ② , taking total expectation, we have

[TABLE]

where $r_{k}=\mathbb{E}\left[f(x_{k})-f(x_{*})\right]$ Note that the second inequality follows since for $h\in\mathbb{R}^{n}$ the following holds $\|h\|_{2}\geq\frac{1}{\sqrt{n}}\|h\|_{1}$. Lastly, by subtracting $f(x_{*})$ and taking expectation of (24), we have

[TABLE]

Thus, we have that

[TABLE]

Note that by setting $t$ to be the the first two upper bounds in (22), the numerator of (26) is lower bounded as

[TABLE]

Moreover, by setting $t$ to the last two upper bounds in (22), the denominator of (26) is lower bounded as

[TABLE]

Then $\frac{1}{r_{k+1}}-\frac{1}{r_{k}}\geq\frac{\min_{i}\frac{p_{i}}{v_{i}}}{8R_{0}^{2}n}=a$ then $\frac{1}{r_{k}}\geq\frac{1}{r_{0}}+ka$. Hence if

[TABLE]

we have $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$. ∎

Appendix D Strongly Convex Case

Theorem 5**.**

Let Assumptions 1 and 3 be satisfied. Choose $\alpha_{i_{k}}=\frac{\alpha_{0}}{v_{{i_{k}}}}\sqrt{2\mu\left(f(x_{k})-f(x_{*})\right)}$, where $0<\alpha_{0}<\frac{2}{\sum_{i=1}^{n}\frac{p_{i}L_{i}}{v_{i}^{2}}}$. If

[TABLE]

then $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$.

Proof.

Since $f$ is $\lambda$-strongly convex, then:

[TABLE]

where $\|h\|_{p^{-2}\circ v^{2}}^{2}=\sum_{i=1}^{n}\frac{v_{i}^{2}}{p_{i}^{2}}h_{i}^{2}$. Minimizing the right hand side in $h$ and substituting again, we have the inequality:

[TABLE]

By subtracting $f(x_{*})$ from both sides of Lemma 1 we get

[TABLE]

From ① , we have

[TABLE]

Since the following holds:

[TABLE]

Thus, combining the previous result with (31), we have

[TABLE]

From ② , we have

[TABLE]

Lastly, substituting the results of ① and ② in (30) and taking the expectation with respect to $x_{k}$ where $r_{k}=\mathbb{E}[f(x_{k})-f(x_{*})]$, we get

[TABLE]

[TABLE]

then $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$. ∎

Theorem 6**.**

Let Assumptions 1 and 3 be satisfied. Choose $\alpha_{i_{k}}=\frac{|f(x_{k}+te_{i_{k}})-f(x_{k})|}{tv_{i_{k}}}$ for some small t, $\mu=\frac{\lambda}{\max_{i}\frac{v_{i}}{p_{i}^{2}}}$. If

[TABLE]

then $\mathbb{E}\left[f(x_{k})-f(x_{*})\right]\leq\epsilon$.

Proof.

Taking the expectation of (19) conditioned on $x_{k}$, we have:

[TABLE]

Taking total expectation in ① we have

[TABLE]

where $r_{0}=\mathbb{E}\left[f(x_{0})\right]-f(x_{*})$. Note that the first inequality follows by Jensen’s inequality. The last inequality follows from the fact that $f$ has Lipschitz gradient, with Lispchitz constant $nL$. As for ② , note that since $\|h\|_{v\circ p^{-1}}^{2}\leq\max_{i}\frac{v_{i}}{p_{i}}\|h\|_{2}^{2}$; thus, by strong convexity we have that

[TABLE]

where $\|h\|_{p^{-1}\circ v}^{2}=\sum_{i=1}^{n}\frac{v_{i}}{p_{i}}h_{i}^{2}$ and $\|\nabla f(x_{k})\|_{p\circ v^{-1}}\stackrel{{\scriptstyle\text{def}}}{{=}}\sqrt{\sum_{i=1}^{n}\frac{p_{i}}{v_{i}}\nabla_{i}f(x_{k})^{2}}$. Thus, it follows that $\|\nabla f(x_{k})\|^{2}_{p\circ v^{-1}}\geq 2\mu(f(x_{k})-f(x_{*}))$. Thereafter, taking the expectation of ② and with the use of the total expectation rule, ② can be lower bounded as follows

[TABLE]

By subtracting $f(x_{*})$ in (36) and taking expectation, we get

[TABLE]

Thus, by choosing $t$ such that $\left(t\max_{i}\sqrt{p_{i}}\sqrt{2n{L}r_{0}}+\frac{t^{2}}{8}\sum_{i=1}^{n}\frac{p_{i}L_{i}^{3}}{v_{i}^{2}}\right)\frac{1}{\mu}\leq\frac{\epsilon}{2}$ and if

[TABLE]

we have $\mathbb{E}\left[f(x_{k})\right]-f(x_{*})\leq\epsilon$. ∎

Appendix E Estimating the Lipschitz Smoothness Constant

As we explain in the Section Experiments, we do not have an access to the Lipschitz smoothness constants of the function for the continuous control case. Instead, we estimate them using the points which have been evaluated. Our experiments suggest that, estimating the function directly using a parametric family works better than estimating smoothness constants directly. In other words, we estimate the function to be optimized directly using a parametric family ($\hat{f}(\cdot,\theta)$) and use its Lipschitz smoothness constants as an estimate.

Consider the DFO step $n$; using the queried sampled $\{x_{i},f(x_{i})\}_{i\in[n-1]}$, we can estimate the function of interest by solving empirical risk minimization problem as:

[TABLE]

We use the resulting $\hat{f}(\cdot,\theta_{n})$ to compute Lipschitz smoothness constants and importance sampling weights. The parametric family we consider is a multi-layer perceptron with single hidden layer and tanh non-linearity. Input dimension is the policy size, output dimension is 1 and hidden layer dimension is chosen as 16 for Swimmer-v1 and Hopper-v1, 64 for HalfCheetah-v1 and Ant-v1, and 256 for Humanoid-v1. Learning has been performed using ADAM (?) optimizer at each iteration. Step size (learning rate) has been chosen as $10^{-3}$ for all experiments. At each iteration we choose a $32$ random samples in uniform from $x_{0},\dots,x_{n-1}$ and use it as a batch.