\section{Introduction}
Derivative-free or zeroth-order optimization \cite{rosenbrock1960automatic,brent1973algorithms,spall2003introduction} is one of the oldest areas in optimization, which constantly attracts attention of the learning community, mostly in connection to the online learning in the bandit setup \cite{bubeck2012regret}. 
We study stochastic derivative-free optimization problems in a two-point feedback situation, considered by \cite{agarwal2010optimal,duchi2015optimal,shamir2017optimal} in the learning community and by \cite{nesterov2017random,stich2011linear,ghadimi2013stochastic,ghadimi2016mini-batch,gasnikov2016gradient-free} in the optimization community. 
Two-point setup allows to prove complexity bounds, which typically coincide with the complexity bounds for gradient-based algorithms up to a small-degree polynomial of $n$, where $n$ is the dimension of the decision variable. 
On the contrary, problems with one-point feedback are harder and complexity bounds for such problems either have worse dependence on $n$, or worse dependence on the desired accuracy of the solution, see \cite{nemirovsky1983problem,protasov1996algorithms,flaxman2005online, agarwal2011stochastic,jamieson2012query,shamir2013complexity,liang2014lzeroth-order,bach16highly-smooth,bubeck2017kernel-based} and the references therein.

More precisely, we consider the following optimization problem
\begin{equation}
\label{eq:PrSt}
   \min_{x\in\mathbb R^n} \left\{ f(x) := \mathbb E_{\xi}[F(x,\xi)] = \int_{\mathcal{X}}F(x,\xi)dP(x) \right\}, 
\end{equation}
where $\xi$ is a random vector with probability distribution $P(\xi)$, $\xi \in \mathcal{X}$, and the function $f(x)$ is closed  and  convex. Note that $F(x,\xi)$ can be non-convex at $x$ with positive probability. Moreover, we assume that, for $P$ almost every $\xi$, the function $F(x,\xi)$ has gradient $g(x,\xi)$, which is $L(\xi)$-Lipschitz continuous with respect to the Euclidean norm and $L_2:=\sqrt{\mathbb E_{\xi} L(\xi)^2 } < +\infty$. Under this assumptions, $\mathbb E_{\xi}g(x,\xi) = \nabla f(x)$ and $f$ has $L_2$-Lipschitz continuous gradient with respect to the Euclidean norm. Also we assume that
\begin{equation}
\label{stoch_assumption_on_variance}
    \mathbb E_{\xi}[\|g(x,\xi) - \nabla f(x)\|_2^2] \leqslant \sigma^2,
\end{equation}
where $\|\cdot\|_2$ is the Euclidean norm. We emphasize that, unlike \cite{duchi2015optimal}, we \textit{do not assume that $\mathbb E_\xi\left[\|g(x,\xi)\|_2^2\right]$ is bounded} since it is not the case for many unconstrained optimization problems, e.g. for deterministic quadratic optimization problems.

Finally, we assume that an optimization procedure, given a pair of points $(x,y) \in \mathbb R^{2n}$, can obtain a pair of noisy stochastic realizations $(\widetilde{f}(x,\xi),\widetilde{f}(y,\xi))$ of the objective value $f$, which we refer to as \textit{oracle call}. Here
\begin{equation}
\label{eq:tf_def}
\widetilde{f}(x,\xi) = F(x,\xi) + \eta(x,\xi), \quad |\eta(x,\xi)| \leqslant  \Delta, \;\forall x \in \mathbb R^n, \; \text{a.s.\ in } \xi,
\end{equation}
and $\xi$ is independently drawn from $P$. This makes our problem more complicated than problems studied in the literature. Not only we have  stochastic noise in the problem \eqref{eq:PrSt}, but an additional noise $\eta(x,\xi)$, which can be adversarial (see Section~\ref{sec:role_of_parameters} for the detailed explanation).

We notice that our model of two-point feedback oracle is pretty general and covers deterministic exact oracle and even specific types of one-point feedback oracle. For example, if the function $F(x,\xi)$ is separable, i.e.\ $F(x,\xi) = f(x) + h(\xi)$, where $\mathbb E_\xi\left[h(\xi)\right] = 0$, $|h(\xi)|\le \frac{\Delta}{2}$ for all $\xi$ and the oracle gives us $F(x,\xi)$ in the given point $x$, then for all $\xi_1,\xi_2$ we can define $\widetilde{f}(x,\xi_1) = F(x,\xi_1)$ and $\widetilde{f}(y,\xi_2) = F(y,\xi_2) = F(y,\xi_1) + h(\xi_2) - h(\xi_1)$. Since $|h(\xi_2)-h(\xi_1)| \le |h(\xi_2)| + |h(\xi_1)| \le \Delta$ we can use representation \eqref{eq:tf_def} omitting dependence of $\eta(x,\xi_1)$ on $\xi_2$, because in our analysis we only rely on the fact that $|\eta(x,\xi)| \le \Delta$ for all $x\in\mathbb R^n$ almost surely in $\xi$. Moreover, such oracle can be met in practice, since rounding errors can be put into this context. e.g.\ by artificial adding random bit modulo $2$ to the last bit in machine's number representation format (see \cite{gasnikov2016stochastic} for details).

As it is known \cite{lan2012optimal,devolder2011stochastic,dvurechensky2016stochastic}, if the stochastic approximation $g(x,\xi)$ for the gradient of $f$ is available, an accelerated gradient method has oracle complexity bound (i.e. overall number of first-order oracle calls) $O\left(\max\left\{\sqrt{L_2R_2^2/\varepsilon},\sigma^2R_2^2/\varepsilon^2\right\}\right)$, where $\varepsilon$ is the target optimization error, the goal being to find such $\hat{x}$ that $\mathbb E f(\hat{x}) - f^* \leq \varepsilon$. Here $f^*$ is the global optimal value of $f$.
The question, to which we give a positive answer in this paper, is as follows.

\textit{Is it possible to solve a stochastic optimization problem with the same $\varepsilon$-dependence in the iteration and sample complexity and only noisy observations of the objective value?} 

Problem \eqref{eq:PrSt} is an unconstrained problem and the first choice of geometric setup is usually given by Euclidean norm. Surprisingly, as we show below, using a non-standard proximal setup given by $\|\cdot\|_1$-norm can give some benefits. So, we solve the problem \eqref{eq:PrSt} using two proximal setups \cite{ben-tal2015lectures}, characterized by the value
$p\in \{1,2\}$ and its conjugate $q \in \{2,\infty\}$, given by the identity $\frac{1}{p}+\frac{1}{q} = 1$. The case $p=1$ corresponds to the choice of $\|\cdot\|_1$-norm in $\mathbb R^n$ and corresponding prox-function, which is strongly convex with respect to this norm (we provide the details below). The case $p=2$ corresponds to the choice of the Euclidean $\|\cdot\|_2$-norm in $\mathbb R^n$ and squared Euclidean norm as the prox-function.

\subsection{Related Work}

Non-smooth deterministic and stochastic problems in the two-point derivative-free setting was considered by \cite{nesterov2017random}.\footnote{We list the references in the order of the date of the first appearance, but not in the order of the date of official publication.} 
Non-smooth stochastic problems were considered by \cite{shamir2017optimal} and independently by \cite{bayandina2017gradient-free}, the latter paper considering also problems with additional noise of an unknown nature in the objective value.
\cite{duchi2015optimal} consider the smooth stochastic optimization problems, yet under additional quite restrictive assumption $\mathbb E_\xi\left[\|g(x,\xi)\|_*^2\right] \le +\infty$. Their bound was improved by \cite{gasnikov2016gradient-free,gasnikov2017stochastic} for the problems with non-Euclidean proximal setup and noise in the objective value.  Strongly convex problems with different smoothness assumptions were considered by \cite{gasnikov2017stochastic,bayandina2017gradient-free}. 
Smooth stochastic convex optimization problems, without the assumption that $\mathbb E \|g(x,\xi)\|_2^2 < +\infty$, were studied by \cite{ghadimi2016mini-batch,ghadimi2013stochastic} for the Euclidean case.
Accelerated and non-accelerated derivative-free method, but for deterministic problems was proposed in \cite{nesterov2017random} and extended in \cite{bogolubsky2016learning,dvurechensky2017randomized} for the case of additional bounded noise in the function value. 
Table~\ref{tbl:summ_and_comp} presents a detailed comparison of our results and results in the literature on two-point feedback derivative-free optimization and assumptions, under which they are obtained.  



{
\renewcommand{\arraystretch}{1.7}
\begin{table}
\small
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Method & \begin{tabular}{c}\iffalse Additional \fi Assumptions \end{tabular} & Oracle complexity, $\widetilde{O}\left( \cdot \right)$ & $p=1$ & Stoch. & Noise\\
\hline
\makecell{MD \\ {\tiny[\cite{duchi2015optimal,gasnikov2016gradient-free, gasnikov2017stochastic,gasnikov2016stochastic}]}\\
{\tiny[\cite{shamir2017optimal,bayandina2017gradient-free}]}}
&bound. gr. &$\frac{{\color{blue}n^{\frac{2}{q}}}{\color{red}M_2^2}R_p^2}{\varepsilon^2}$ & $\surd$ & $\surd$ & $\surd$ \\
\hline
\makecell{RSPGF\\ {\tiny[\cite{ghadimi2016mini-batch,ghadimi2013stochastic}]}} & bound. var.
& $\max\left\{\frac{{\color{red}n}L_2R_2^2}{\varepsilon},\frac{{\color{red}n}{\color{blue}\sigma^2}R_2^2}{\varepsilon^2}\right\}$ & $\times$ & $\surd$ & $\times$ \\
\hline
\makecell{RS \\ {\tiny\cite{nesterov2017random,bogolubsky2016learning}} } & 
& $\frac{{\color{red}n}L_2R_2^2}{\varepsilon}$& $\times$ & $\times$ & $\surd$  \\
\hline
\makecell{RDFDS {\color{blue}(New!)}\\ {[This paper]}} & bound. var.
& $\max\left\{\frac{{\color{blue}n^{\frac{2}{q}}}L_2R_p^2}{\varepsilon},\frac{{\color{blue}n^{\frac{2}{q}}}{\color{blue}\sigma^2}R_p^2}{\varepsilon^2}\right\}$ & $\surd$ & $\surd$ & $\surd$  \\
\hline
\makecell{
AccRS\\ {\tiny[\cite{nesterov2017random,dvurechensky2017randomized}]}
}
&
& ${\color{red}n}\sqrt{\frac{L_2R_2^2}{\varepsilon}}$ & $\times$ & $\times$ & $\surd$  \\
\hline
\makecell{ARDFDS {\color{blue}(New!)}\\ {[This paper]}} & bound. var. 
& $\max\left\{{\color{blue}n^{\frac{1}{2}+\frac{1}{q}}}\sqrt{\frac{L_2R_p^2}{\varepsilon}},\frac{{\color{blue}n^{\frac{2}{q}}}{\color{blue}\sigma^2}R_p^2}{\varepsilon^2}\right\}$ & $\surd$ & $\surd$ & $\surd$   \\
\hline
\end{tabular}
\caption{Comparison of oracle complexity (total number of oracle calls) of different methods with two point feedback discussed in the paper for convex optimization problems. In the column "Assumptions" we use "bound. gr." to define $\mathbb E_\xi\left[\|g(x,\xi)\|_2^2\right] \le M_2^2$ and "bound. var." to define $\mathbb E\|g(x,\xi) - \nabla f(x)\|_2^2 \leqslant \sigma^2$. Column $p=1$ corresponds to the support of non-Euclidean setup, column "Stoch." to the support of stochastic optimization methods, "Noise" corresponds to the support of additional noise of an unknown nature.
\label{tbl:summ_and_comp}
}
\end{center}
\end{table}
}

We also mention the works \cite{nemirovsky1983problem,protasov1996algorithms,flaxman2005online,saha2011improved,dekel2015bandit,gasnikov2017stochastic,agarwal2011stochastic,liang2014lzeroth-order,belloni2015escaping,bubeck2017kernel-based,shamir2013complexity,jamieson2012query,hazan2014bandit,bach16highly-smooth,jamieson2012query} who study derivative-free optimization with one-point feedback in different settings, and works \cite{nesterov2005smooth,allen2014linear} on coupling non-accelerated methods to obtain acceleration, which inspired our work. A more detailed description of the related work can be found in the Appendix.


\subsection{Our Contributions}
As our main contribution, we propose an accelerated method for smooth stochastic derivative-free optimization, which we call Accelerated Randomized Derivative-Free Directional Search (ARDFDS). Our method has the complexity bound 
\begin{equation}
\label{eq:ARDFDSComplInformal}
\widetilde{O}\left(\max\left\{n^{\frac12+\frac{1}{q}}\sqrt{\frac{L_2R_p^2}{\varepsilon}},\frac{n^{\frac{2}{q}}\sigma^2R_p^2}{\varepsilon^2}\right\}\right),
\end{equation}
where $R_p$ characterizes the distance in $\|\cdot\|_p$ - norm between the starting point of the algorithm and a solution to \eqref{eq:PrSt}. 
In the Euclidean case $p=q=2$, the first term in the above bound has better dependence on $n$, $\varepsilon$,  $L_2$ and $R_2$ than the bound in \cite{ghadimi2016mini-batch,ghadimi2013stochastic}. Unlike these papers, our bound also covers the non-euclidean case $p=1$, $q=\infty$ and due to that allows to obtain better complexity bounds. To illustrate this, let us start method from a point $x_0$ and define the sparsity $s$ of the vector $x_0 - x^*$, i.e. $\|x_0-x^*\|_1 \leq s \cdot \|x_0-x^*\|_2$ and $1 \leq s \leq \sqrt{n}$. Then the complexity of our method for $p=1,q=\infty$ is $\widetilde{O}\left(\max\left\{\sqrt{\frac{ns^2L_2\|x_0-x^*\|_2^2}{\varepsilon}},\frac{s^2\sigma^2\|x_0-x^*\|_2^2}{\varepsilon^2}\right\}\right)$, which is always no worse than the complexity for $p=q=2$ $\widetilde{O}\left(\max\left\{\sqrt{\frac{n^2L_2\|x_0-x^*\|_2^2}{\varepsilon}},\frac{n\sigma^2\|x_0-x^*\|_2^2}{\varepsilon^2}\right\}\right)$ and allows to gain up to $\sqrt{n}$ if $s$ is close to 1. Notably, this is done automatically, without any prior knowledge of $s$.

Unlike other authors, we consider additional, possibly adversarial noise $\eta(x,\xi)$ in the objective value and analyze how this noise affects the convergence rate estimates. We emphasize that even if noise is \textit{uncontrolled}, e.g. we do not know noise level $\Delta$ or we cannot reduce $\Delta$, we can run our algorithms and still expect the convergence (see Section~\ref{sec:role_of_parameters} for the detailed explanation). This is important when the objective is given as a solution to some auxiliary problem which can't be solved exactly, e.g. in bi-level optimization or reinforcement learning. 
It should also be mentioned that our assumption $\mathbb E_{\xi} [L(\xi)^2] < +\infty$ is weaker than the assumption $L(\xi) \leq L_2$ a.s. in $\xi$, which is used in \cite{ghadimi2016mini-batch,ghadimi2013stochastic}. 

As our second contribution, we propose a non-accelerated Randomized Derivative-Free Directional Search (RDFDS) method with the complexity bound
\begin{equation}
\label{eq:RDFDSComplInformal}
\widetilde{O}\left(\max\left\{\frac{n^{\frac{2}{q}}L_2R_p^2}{\varepsilon},\frac{n^{\frac{2}{q}}\sigma^2R_p^2}{\varepsilon^2}\right\}\right),
\end{equation}
where unlike \cite{ghadimi2016mini-batch,ghadimi2013stochastic} a non-euclidean case $p=1$, $q=\infty$ with the gain in the complexity up to the factor of $n$ is possible. Interestingly, in this case, we obtain a nearly dimension independent complexity bound despite we use only noisy function value observations.
\subsubsection{\textbf{Why it is important to improve the first term of the maximum?}}
\begin{enumerate}
    \item\textbf{Acceleration when $n$ is big.} First of all, we notice that the first term dominates the second term when $\sigma^2 \le \frac{\varepsilon^{\frac{3}{2}}n^{\frac{1}{2}-\frac{1}{q}}\sqrt{L_2}}{R_p}$ in the accelerated case and $\sigma^2 \le \varepsilon L_2$ in the non-accelerated case which could be met in practice if $\varepsilon$, $L_2$ and $n$ are big enough in comparison with $R_p$. For ARDFDS with $p=1,q=\infty$ it means that if we want to get $\varepsilon$-solution with $\varepsilon = 10^{-3}$ and $L_2 = 100$, $R_p = 10$, $n = 10 000$ (or bigger), then the variance should satisfy $\sigma^2 \le 10^{-1}$ in order to have accelerated rate for ARDFDS, which, actually, is not a very restrictive assumption.
    \item\textbf{Better dimension-dependence in the deterministic case.} We underline also that even in the deterministic case $\sigma = 0$ without additive noise both our non-accelerated and accelerated complexity bounds for $p<2$ and, in particular, for $p=1$ are new. Moreover, disregarding $\ln n$ factors, the existing bounds \cite{nesterov2017random} are $n/s^2$ and $n/s$ times worse than our new bounds respectively in non-accelerated and accelerated cases, where $s\in [1,\sqrt{n}]$. Importantly, in the non-accelerated case our bound is dimension-independent up to a $\ln n$ factor.
    \item\textbf{Parallel computation of mini-batches makes acceleration reasonable when $\sigma^2$ is not small.} Imagine the situation when one can have an access to $k \ge m = n^{\frac{1}{q}- \frac{1}{2}}\sigma^2\sqrt{\frac{R_p^2}{\varepsilon^3 L_2}}$ processors (or workers). For example, if $\sigma^2=1$, which is not small, $n=10000$, $R_p = 10$, $\varepsilon = 10^{-3}$ and $L_2 = 100$, then we need $k=10^{2.5}\approx 316$ processors which is just nothing for modern supercomputers and clusters that often have $\sim 10^5-10^6$ processors. It means that one can call zeroth-order oracle on each processor, compute own stochastic approximation of the gradient which is described below and after that send it to one chosen processor (directly or via others processors) which just take an average of received data. In other words, one can compute mini-batches in parallel. Such technique gives us an opportunity to reduce real time of computing stochastic approximation of the gradient on each iteration and to make running time equal to $O(N)$ (i.e.\ iteration complexity). But $N$ is equal exactly to the first term in our bounds (see Table~\ref{tbl:ARDFDS_params} for the details) which means that in such situation it is possible to make real running time of the method proportional to the first term in bounds we obtained. More rigorously, it means that if we have $k$ processors than due to parallel computing of mini-batches we can obtain the real working time equal to $\widetilde{O}\left(\max\left\{n^{\frac12+\frac{1}{q}}\sqrt{\frac{L_2R_p^2}{\varepsilon}},\frac{1}{k}\frac{n^{\frac{2}{q}}\sigma^2R_p^2}{\varepsilon^2}\right\}\right)$ and $\widetilde{O}\left(\max\left\{\frac{n^{\frac{2}{q}}L_2R_p^2}{\varepsilon},\frac{1}{k}\frac{n^{\frac{2}{q}}\sigma^2R_p^2}{\varepsilon^2}\right\}\right)$ for ARDFDS and RDFDS respectively. We notice that $k$ could be even smaller that $m = n^{\frac{1}{q}- \frac{1}{2}}\sigma^2\sqrt{\frac{R_p^2}{\varepsilon^3 L_2}}$ to have that the first term of maximum is bigger than the second one (see the first reason in this list).
\end{enumerate}

\section{Algorithms for Stochastic Convex Optimization}
\subsection{Preliminaries}\label{sec:preliminaries}

\noindent \textbf{Proximal setup.}
Let $p\in[1,2]$ and $\|x\|_p$ be the $p$-norm in $\mathbb R^n$ defined as $\|x\|_p^p = \sum\limits_{i=1}^n|x_i|^p$
and $\|\cdot\|_{q}$ be its dual, defined by $\|g\|_{q} = \max\limits_{x} \big\{ \langle g, x \rangle, \| x \|_p \leq 1 \big\}$, where $q \in [2,\infty]$ is the conjugate number to $p$, given by $\frac{1}{p} + \frac{1}{q} = 1$, and, for $q = \infty$, we define $\|x\|_\infty = \max\limits_{i=1,\ldots,n}|x_i|$. We choose a \textit{prox-function} $d(x)$, which is continuous, convex on $\mathbb R^n$ and is $1$-strongly convex on $\mathbb R^n$ with respect to $\|\cdot\|_p$, i.e., for any $x, y \in \mathbb R^n$ $d(y)-d(x) -\langle \nabla d(x) ,y-x \rangle \geq \frac12\|y-x\|_p^2$.
Without loss of generality, we assume that $\min\limits_{x\in \mathbb R^n} d(x) = 0$.
We define also the corresponding \textit{Bregman divergence} $V[z] (x) = d(x) - d(z) - \langle \nabla d(z), x - z \rangle$, $x, z \in \mathbb R^n$. Note that, by the strong convexity of $d$,
\begin{equation}
\label{eq:VStrConv}
V[z] (x) \geq \frac{1}{2}\|x-z\|_p^2, \quad x, z \in \mathbb R^n.
\end{equation}
For the case $p=1$, we choose the following prox function \cite{ben-tal2015lectures} $d(x) = \frac{{\rm e}n^{(\kappa-1)(2-\kappa)/\kappa}\ln n}{2} \|x\|_\kappa^2$, where $\kappa=1 + \frac{1}{\ln n}$ and, for the case $p=2$, we choose the prox-function to be proportional to the Euclidean norm: $d(x) = \frac{1}{2}\|x\|_2^2.$


\noindent \textbf{Main technical lemma.}
In our proofs of complexity bounds, we rely on the following lemma. The proof is rather technical and is provided in the appendix. 
\begin{lemma}
		\label{Lm:MainTechLM}
    Let $e \in RS_2(1)$, i.e be a random vector uniformly distributed on the surface of the unit Euclidean sphere in $\mathbb R^n$, $p\in[1,2]$ and $q$ be given by $\frac{1}{p}+\frac{1}{q} = 1$. Then, for $n \geqslant8$,		
		\begin{equation}\label{assumption_e_norm}
        \mathbb E_e\|e\|_q^2 \leq \rho_n,
    \end{equation}
    \begin{equation}\label{assumption_e_product}
        \mathbb E_e\left(\langle s,e \rangle^2 \|e\|_q^2 \right) \leq \frac{6\rho_n}{n} \|s\|_2^2,     \quad \forall s \in \mathbb R^n.
    \end{equation}
		with $\rho_n \overset{\text{eqdef}}{=} \min\{q-1,\,16\ln n - 8\}n^{\frac{2}{q}-1}$.
\end{lemma}

\noindent \textbf{Stochastic approximation of the gradient.} 
Based on the noisy observations \eqref{eq:tf_def} of the objective value, we form the following stochastic approximation of $\nabla f(x)$
\begin{equation}
\label{eq:FinDiffStocGrad}
\widetilde{\nabla}^m f^t(x) = \frac{1}{m}\sum\limits_{i=1}^m\frac{\widetilde{f}(x+te,\xi_i)-\widetilde{f}(x,\xi_i)}{t}e,
\end{equation}
where $e \in RS_2(1)$, $\xi_i$, $i=1,...,m$ are independent realizations of $\xi$, $m$ is the \textit{batch size}, $t$ is some small positive parameter, which we call \textit{smoothing parameter}.

\subsection{Algorithms and Main Theorems}
Our Accelerated Randomized Derivative-Free Directional Search (ARDFDS) method is listed as Algorithm~\ref{Alg:ARDFDS}.
\begin{algorithm}
	\caption{Accelerated Randomized Derivative-Free Directional Search (ARDFDS)}
	\label{Alg:ARDFDS
	\begin{algorithmic}[1]
		\REQUIRE $x_0$~---starting point; $N$~--- number of iterations; $m \geq 1$~--- batch size; $t > 0$~--- smoothing parameter, $\{\alpha_k\}_{k=1}^N$~--- stepsizes.
		\ENSURE point $y_N$.
		\STATE $y_0 \leftarrow x_0, \, z_0 \leftarrow x_0$.
		\FOR{$k=0,\, \dots, \, N-1$.}
	
		\STATE $\tau_k \leftarrow \frac{2}{k+2}$
		\STATE Generate $e_{k+1} \in RS_2\left( 1 \right)$ independently from previous iterations and $\xi_i$, $i=1,...,m$ -- independent realizations of $\xi$. 
		\STATE $x_{k+1} \leftarrow \tau_kz_k + (1-\tau_k)y_k $.
		\STATE Calculate $\widetilde{\nabla}^m f^t(x_{k+1})$ given in \eqref{eq:FinDiffStocGrad}.
		\STATE $y_{k+1} \leftarrow x_{k+1}-\frac{1}{2L_2}\widetilde{\nabla}^m f^t(x_{k+1})$.
		\STATE $z_{k+1} \leftarrow \mathop{\arg\!\min}\limits_{z\in\mathbb R^n} \Big\{ \alpha_{k+1} n \left\langle \widetilde{\nabla}^m f^t(x_{k+1}), \, z-z_k \right\rangle+V[z_{k}] \left( z \right)\Big\}$
		\ENDFOR
		\RETURN $y_N$
	\end{algorithmic}
\end{algorithm}
\begin{theorem}
\label{Th:ARDFDSConv}
    Let ARDFDS be applied to solve problem \eqref{eq:PrSt} with $\alpha_{k+1} = \frac{k+2}{96n^2\rho_n L_2}$ which implies $\tau_k = \frac{2}{k+2} = \frac{1}{48\alpha_{k+1}n^2\rho_n L_2}$. If we set $\Theta_p = V[z_0](x^*)$ which is defined by the chosen proximal setup, then
    \begin{equation}\label{eq:ARDFDSConv}
        \begin{array}{rcll}
            \mathbb E[f(y_N)] - f(x^*) 
        &\leqslant& \frac{384 n^2\rho_nL_2\Theta_p}{N^2} +\frac{384N}{nL_2}\frac{\sigma^2}{m} + \frac{12\sqrt{2n\Theta_p}}{N^2}\left( \frac{L_2t}{2} + \frac{2\Delta }{t} \right)&\\ &&\quad+ 
				 \frac{6N}{L_2} \left( L_2^2t^2+\frac{16\Delta^2}{t^2} \right)+ \frac{N^2}{24n\rho_nL_2} \left( L_2^2t^2 + \frac{16\Delta^2 }{t^2} \right),& \forall n\ge 8.
        \end{array}
    \end{equation}
\end{theorem}

\textbf{Sketch of the proof of the Theorem~\ref{Th:ARDFDSConv} for the deterministic case with exact oracle.} To introduce the main idea and intuition behind our proof of the Theorem~\ref{Th:ARDFDSConv} we consider here the most simple situation when $\sigma^2 = 0$ (i.e.\ deterministic problem) and $\Delta = 0$ (i.e.\ exact oracle). In this case technical details corresponding to the control of variance and noise do not disturb us to show what is really important to get better dimension-dependence in the oracle complexity.

First of all, we notice that in this special case we don't need to do mini-batches of size $m>1$, so, we can set $m=1$ and rewrite $\widetilde{\nabla}^m f(x)$ as
\begin{equation}\label{eq:stochastic_approx_deterministic}
    \widetilde{\nabla}^m f(x) \equiv \widetilde{\nabla}f(x) = \frac{f(x+te) - f(x)}{t}e = \left(\langle\nabla f(x), e\rangle + \theta(x,t,e)\right)e,
\end{equation}
where $\theta(x,t,e) = \frac{f(x+te) - f(x)}{t} - \langle\nabla f(x), e\rangle = \frac{1}{t}\left(f(x+te) - f(x) - \langle\nabla f(x),te\rangle\right)$. Note that due to $L_2$-smoothness we have
\begin{equation}\label{eq:theta_upper_bound_deterministic}
    |\theta(x,t,e)| \le \frac{1}{t} \cdot\frac{L_2t^2}{2} = \frac{L_2 t}{2}.
\end{equation}

The technique of linear coupling from \cite{allen2014linear} inspired us to study how this method behaves when we use gradient approximation by finite differences like in \eqref{eq:stochastic_approx_deterministic} instead of true gradient in the method. In the linear coupling method one can choose any norm to define $L$-smoothness of the objective. However, we are also motivated to have clear comparison in the deterministic case with other methods like RSPGF from \cite{ghadimi2013stochastic,ghadimi2016mini-batch}, RS from \cite{nesterov2017random,bogolubsky2016learning} and AccRS from \cite{nesterov2017random,dvurechensky2017randomized}, which use $L_2$-smoothness with respect to Euclidean norm. Therefore, we decide to use $L_2$-smoothness of the objective in the $\ell_2$-norm and, as the consequence, it is natural for our purposes to consider gradient step in the linear coupling method with respect to Euclidean norm, since it is a well-known fact the gradient descent step with constant stepsize $\frac{1}{L_2}$ could be obtained in view of minimization of the quadratic upper bound coming from Lipschitz-continuity of the gradient: $\mathop{\arg\!\min}_{x\in\mathbb R^n}\left\{f(x_k) + \langle\nabla f(x_k), x - x_k\rangle + \frac{L_2}{2}\|x_k - x_*\|_2^2\right\} = x_k - \frac{1}{L_2}\nabla f(x_k).$ Thus, in the line $7$ of the ARDFDS we set $y_{k+1} = x_{k+1} - \frac{1}{2L_2}\widetilde{\nabla}f(x)$. We substitute classic coefficient $-\frac{1}{L_2}$ by $-\frac{1}{2L_2}$ just because of some technical issues that we analyse accurately in the full prove which one can find in the appendix. But if we also chose mirror descent step with Bregman divergence corresponding to the $\ell_2$ proximal setup, we will obtain the same bound as for AccRS from \cite{nesterov2017random,dvurechensky2017randomized}. Then we switch gears and come up with the following counter-intuitive idea. Let consider \textit{mirror descent step with respect to $\ell_p$ proximal setup, where $p\in[1,2)$, but remain the gradient descent step the same as for choosing Lipschitz-continuity of the gradient with respect to $\ell_2$ norm.} Next we emphasize the place in our proof where this idea naturally appears.

The first step of the proof is the standard analysis of the one iteration of mirror descent. We have
\begin{eqnarray*}
    \alpha_{k+1} n\langle \widetilde{\nabla}f^t(x_{k+1}), \, z_k-x_*\rangle
		=  \alpha_{k+1} n \langle \widetilde{\nabla}f^t(x_{k+1}), \,  z_k-z_{k+1}\rangle + \alpha_{k+1} n \langle \widetilde{\nabla}f^t(x_{k+1}), \, z_{k+1}-x_*\rangle\\
		\overset{\text{\ding{172}}}{\leqslant}  \alpha_{k+1} n \langle \widetilde{\nabla}f^t(x_{k+1}), \, z_k-z_{k+1}\rangle + \langle -\nabla V[z_k](z_{k+1}), \,  z_{k+1}-x_*\rangle
		\overset{\text{\ding{173}}}{=}  \alpha_{k+1} n \langle \widetilde{\nabla}f^t(x_{k+1}), \, z_k-z_{k+1}\rangle\\
		+ V[z_k](x_*) - V[z_{k+1}](x_*)
		- V[z_k](z_{k+1})
		\overset{\text{\ding{174}}}{\leqslant} \left(\alpha_{k+1} n \langle \widetilde{\nabla}f^t(x_{k+1}), \,  z_k-z_{k+1}\rangle - \frac{1}{2}\|z_k-z_{k+1}\|_p^2\right)\\ + V[z_k](x_*)
		- V[z_{k+1}](x_*)
		\overset{\text{\ding{175}}}{\leqslant} \frac{\alpha_{k+1}^2n^2}{2}\| \widetilde{\nabla}f^t(x_{k+1})\|_q^2 + V[z_k](x_*) - V[z_{k+1}](x_*),
\end{eqnarray*}
where $\text{\ding{172}}$ follows from the definition of $z_{k+1}$, whence $\langle \nabla V[z_k](z_{k+1}) +  \alpha_{k+1} n \widetilde{\nabla}f^t(x_{k+1}), \, u - z_{k+1}\rangle \geqslant 0$ for all $u\in \mathbb R^n$ and, as a consequence, for $u=x_*$; $\text{\ding{173}}$ follows from the "magic identity" Fact 5.3.3 in \cite{ben-tal2015lectures} for the Bregman divergence; $\text{\ding{174}}$ follows from \eqref{eq:VStrConv}; and $\text{\ding{175}}$ follows from the Fenchel inequality $\zeta\langle s,z \rangle - \frac{1}{2}\|z\|_p^2 \leq \frac{\zeta^2}{2}\|s\|_q^2$. Letting $\theta_{k+1} = \theta(x_{k+1},t,e_{k+1})$, we continue
\begin{eqnarray*}
    \alpha_{k+1} n\langle \widetilde{\nabla}f^t(x_{k+1}), \, z_k-x_*\rangle =\alpha_{k+1}n\langle\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}, z_k - x_*\rangle
    + \alpha_{k+1}n\theta_{k+1}\langle e_{k+1},z_k - x_*\rangle\\
    \ge \alpha_{k+1}n\langle\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}, z_k - x_*\rangle
    - \alpha_{k+1}n|\theta_{k+1}|\cdot |\langle e_{k+1},z_k - x_*\rangle|\\
    \overset{\eqref{eq:theta_upper_bound_deterministic}}{\ge} \alpha_{k+1}n\langle\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}, z_k - x_*\rangle
    - \frac{\alpha_{k+1}nL_2t}{2}|\langle e_{k+1},z_k - x_*\rangle|.
\end{eqnarray*}
Using triangle inequality and $\|a+b\|_q^2 \le \left(\|a\|_q + \|b\|_q\right)^2 \le 2\|a\|_q^2 + 2\|b\|_q^2$ we obtain $\|\widetilde{\nabla}f(x_{k+1})\|_q^2 \overset{\eqref{eq:stochastic_approx_deterministic}}{\le} 2\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2 + 2\|\theta_{k+1}e_{k+1}\|_q^2 \overset{\eqref{eq:theta_upper_bound_deterministic}}{\le} 2\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2 + \frac{L_2^2 t^2}{2}\|e_{k+1}\|_q^2.$
Putting all together we get 
\begin{eqnarray*}
    \alpha_{k+1}n\langle\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}, z_k - x_*\rangle &\le& \alpha_{k+1}^2n^2\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2 + V[z_k](x_*)\\
    &&- V[z_{k+1}](x_*) + tW_k,
\end{eqnarray*}
where $W_k = \frac{\alpha_{k+1}nL_2}{2}|\langle e_{k+1}, z_k-x_*\rangle| + \frac{\alpha_{k+1}^2n^2L_2^2t}{4}\|e_{k+1}\|_q^2$. If we take the expectation $\mathbb E_{e_{k+1}}[\cdot]$ with respect randomness coming from $e_{k+1}$ from both sides of the previous inequality, we will establish
\begin{eqnarray*}
    \alpha_{k+1}\langle\nabla f(x_{k+1}),z_k - x_*\rangle &\le& \alpha_{k+1}^2n^2\mathbb E_{e_{k+1}}\left[\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2\right] + V[z_k](x_*)\\
    &&- \mathbb E_{e_{k+1}}\left[V[z_{k+1}](x_*)\right] + t\mathbb E_{e_{k+1}}[W_k],
\end{eqnarray*}
since $\mathbb E_{e_{k+1}}[\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}] = \frac{\nabla f(x_{k+1})}{n}$ that can be obtained from $\mathbb E_{e_{k+1}}[(e_{k+1}^i)^2] = \frac{1}{n}$ and\newline $\mathbb E_{e_{k+1}}[e_{k+1}^ie_{k+1}^j] = 0$ for $i\neq j$, where $e_{k+1}^i$ is the $i$-th component of the $e_{k+1}$. \textit{Note that before this moment our proof looks similar for all $p\in [1,2]$.} \textit{Now it is the right time to reveal the main trick in our proof.}

\textit{The second step of the proof consists of upper bounding $\mathbb E_{e_{k+1}}\left[\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2\right]$ by $f(x_{k+1}) - \mathbb E_{e_{k+1}}[f(y_{k+1})]$ with some multiplicative and additive terms to connect progresses of mirror descent and gradient descent steps. But since we have $L_2$-smoothness defined with respect to Euclidean norm it is more convenient first to upper bound $\mathbb E_{e_{k+1}}\left[\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2\right]$ by $\|\nabla f(x_{k+1})\|_2^2$ with some multiplicative factor \textbf{which is the most important step in the whole proof since it influences the dimension-dependence in the final bound} and after that to upper bound $\|\nabla f(x_{k+1})\|_2^2$ by $f(x_{k+1}) - \mathbb E_{e_{k+1}}[f(y_{k+1})]$ with some multiplicative and additive terms using standard gradient descent analysis.} Our main technical lemma (see Lemma~\ref{Lm:MainTechLM}) helps us here:
\begin{equation*}
    \mathbb E_{e_{k+1}}\left[\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2\right] \overset{\eqref{assumption_e_product}}{\le} \frac{6\rho_n}{n}\|\nabla f(x_{k+1})\|_2^2, \quad \rho_n = \min\left\{q-1,16\ln n - 8\right\}n^{\frac{2}{q}-1},\; \forall n\ge 8.
\end{equation*}
It remains to estimate $\|\nabla f(x_{k+1})\|_2^2$. Using $L_2$-smoothness and invoking gradient step we have
\begin{eqnarray*}
    f(y_{k+1}) \le f(x_{k+1}) - \frac{1}{2L_2}\langle\nabla f(x_{k+1}), \widetilde{\nabla}f(x_{k+1})\rangle + \frac{L_2}{2}\left\|\frac{1}{2L_2}\widetilde{\nabla}f(x_{k+1})\right\|_2^2
    \overset{\eqref{eq:stochastic_approx_deterministic}}{\le} f(x_{k+1})\\ - \frac{1}{2L_2}\langle\nabla f(x_{k+1}),e_{k+1}\rangle^2 - \frac{1}{2L_2}\theta(x_{k+1},t,e_{k+1})\langle\nabla f(x_{k+1}),e_{k+1}\rangle\\
    +L_2\left(\frac{1}{2L_2}\langle\nabla f(x_{k+1}),e_{k+1}\rangle +\frac{1}{2L_2}\theta(x_{k+1},t,e_{k+1})\right)^2
    \overset{\eqref{eq:theta_upper_bound_deterministic}}{\le} f(x_{k+1}) - \frac{1}{4L_2}\langle\nabla f(x_{k+1}),e_{k+1}\rangle^2 + \frac{L_2t^2}{16}.
\end{eqnarray*}
Using Lemma~B.10 from \cite{bogolubsky2016learning} and rearranging the terms we get $\|\nabla f(x_{k+1})\|_2^2 \le 4nL_2\left(f(x_{k+1})-\mathbb E_{e_{k+1}}[f(y_{k+1})]\right) + \frac{L_2nt^2}{4}.$
If we multiply the right hand side of the previous inequality by $2$ we will still obtain right inequality $\|\nabla f(x_{k+1})\|_2^2 \le 8nL_2\left(f(x_{k+1})-\mathbb E_{e_{k+1}}[f(y_{k+1})]\right) + \frac{L_2nt^2}{2}.$ For the deterministic case it is no reason to do so, but we do it to have clear matching with the general case. Putting all together we get
\begin{eqnarray}
    \alpha_{k+1}\langle\nabla f(x_{k+1}),z_k - x_*\rangle &\le& 48\alpha_{k+1}^2n^2\rho_nL_2\left(f(x_{k+1})-\mathbb E_{e_{k+1}}[f(y_{k+1})]\right) + V[z_k](x_*) \notag \\
    &&- \mathbb E_{e_{k+1}}\left[V[z_{k+1}](x_*)\right] + t\widetilde{W}_{k}, \label{eq:one_step_pr1_informal}
\end{eqnarray}
where $\widetilde{W}_k = W_k + 3\alpha_{k+1}^2n^2\rho_nL_2t$.

The remaining part of the proof is very similar to the original proof of linear coupling method, but still hides non-trivial technique how to estimate sum of $\mathbb E\left[\widetilde{W}_k\right]$. However, for the simplicity in this sketch of the proof we omit the details about estimating terms connecting with $\widetilde{W}_k$ (see appendix for the detailed proof). Further,
\begin{eqnarray*}
    \alpha_{k+1}(f(x_{k+1}) - f(x_*)) \leqslant \alpha_{k+1}\langle \nabla f(x_{k+1}) ,\, x_{k+1} - x_* \rangle
    = \alpha_{k+1}\langle\nabla f(x_{k+1}), \, x_{k+1}-z_k\rangle\\ + \alpha_{k+1}\langle\nabla f(x_{k+1}), \, z_{k}-x_*\rangle
	\overset{\text{\ding{172}}}{=} \frac{(1-\tau_k)\alpha_{k+1}}{\tau_k}\langle\nabla f(x_{k+1}), \, y_{k}-x_{k+1}\rangle+\alpha_{k+1}\langle\nabla f(x_{k+1}), \, z_{k}-x_*\rangle\\
	\overset{\text{\ding{173}}}{\leqslant} \frac{(1-\tau_k)\alpha_{k+1}}{\tau_k} (f(y_k)-f(x_{k+1})) + \alpha_{k+1}\langle\nabla f(x_{k+1}), \, z_{k}-x_*\rangle
	\overset{\eqref{eq:one_step_pr1_informal}}{\leqslant} \frac{(1-\tau_k)\alpha_{k+1}}{\tau_k} (f(y_k)-f(x_{k+1}))\\ +  48 \alpha_{k+1}^2 n^2 \rho_n L_2 \left(f(x_{k+1}) - \mathbb E_{e_{k+1}}[f(y_{k+1})]\right)
	+ V[z_k](x_*) - \mathbb E_{e_{k+1}}[V[z_{k+1}](x_*)] + t\widetilde{W}_k \\
	\overset{\text{\ding{174}}}{=} (48 \alpha_{k+1}^2 n^2 \rho_n L_2-\alpha_{k+1})f(y_k) - 48 \alpha_{k+1}^2 n^2 \rho_n L_2\mathbb E_{e_{k+1}}[f(y_{k+1})]\\
	+\alpha_{k+1}f(x_{k+1}) + V[z_k](x_*) - \mathbb E_{e_{k+1}}[V[z_{k+1}](x_*)] + t\widetilde{W}_k.
\end{eqnarray*}
That is, $\text{\ding{172}}$ is since $x_{k+1} := \tau_kz_k+(1-\tau_k)y_k\; \Leftrightarrow\; \tau_k(x_{k+1}-z_k) = (1-\tau_k)(y_k - x_{k+1})$, $\text{\ding{173}}$ follows from the convexity of $f$ and inequality $1-\tau_k\geqslant0$ and $\text{\ding{174}}$ is since $\tau_k = \frac{1}{48\alpha_{k+1}n^2\rho_n L_2}$. Rearranging the terms we get $48 n^2 \rho_n L_2 \alpha_{k+1}^2\mathbb E_{e_{k+1}}[f(y_{k+1})] - (48 n^2 \rho_n L_2 \alpha_{k+1}^2-\alpha_{k+1})f(y_k) - V[z_k](x_*) + \mathbb E_{e_{k+1}}[V[z_{k+1}](x_*)] - t\widetilde{W}_k \leqslant \alpha_{k+1}f(x_*).$
Taking the full expectation from the both sides of the previous inequality and summing the results for $k=0,1,\ldots,N$ we get $48n^2\rho_nL_2\alpha_{l}^2\mathbb E[f(y_N)] + \sum\limits_{k=1}^{N-1}\frac{1}{192n^2\rho_nL_2}\mathbb E[f(y_k)] - V[z_0](x_*) + \mathbb E[V[z_N](x_*)]-t\sum\limits_{k=0}^{N-1}\mathbb E\left[\widetilde{W}_k\right]  \leqslant \sum\limits_{k=0}^{N-1}\alpha_{k+1}f(x_*),$ where we use the simple recursion $48n^2\rho_nL_2\alpha_{k+1}^2 - \alpha_{k+1} + \frac{1}{192n^2\rho_nL_2} = 48n^2\rho_nL_2\alpha_{k}^2$ which trivially follows from our choice of $\alpha_{k+1}$. We define $\Theta_p := V[z_0](x^*)$. Since $\sum\limits_{k=0}^{N-1}\alpha_{k+1} = \frac{N(N+3)}{192n^2\rho_nL_2}$, $V[z_N](x_*) \ge 0$ and for all $k=1,\ldots,N$, $f(y_k) \leqslant f(x^*)$, we obtain
\begin{eqnarray*}
    \frac{(N+1)^2}{192n^2\rho_nL_2}\mathbb E[f(y_N)] &\leqslant& f(x_*)\left(\frac{(N+3)N}{192n^2\rho_nL_2} - \frac{N-1}{192n^2\rho_nL_2}\right) + \Theta_p + t\sum\limits_{k=0}^{N-1}\mathbb E\left[\widetilde{W}_k\right],
\end{eqnarray*}
and after rearranging the terms we get
\begin{equation}\label{eq:main_result_deterministic_informal}
    \mathbb E[f(y_N)] - f(x_*) \le \frac{192n^2\rho_nL_2\Theta_p}{(N+1)^2} + t\cdot\frac{192n^2\rho_nL_2}{(N+1)^2}\sum\limits_{k=0}^{N-1}\mathbb E\left[\widetilde{W}_k\right].
\end{equation}
In the appendix we have a detailed analysis and estimate the second term in the right-hand side of \eqref{eq:main_result_deterministic_informal}. In general, the second term could be made arbitrary small by choosing proper $t$. Therefore, the leading term in the complexity bound \eqref{eq:main_result_deterministic_informal} for deterministic case with exact oracle and small enough $t$ is always $\frac{192n^2\rho_nL_2\Theta_p}{(N+1)^2}$, which means that in this case in order to achieve $\mathbb E\left[f(y_N)\right] - f(x_*) \le \varepsilon$ we need to perform $O\left(\sqrt{n^2\rho_nL_2\Theta_p/\varepsilon}\right)$ iterations of ARDFDS. Furthermore, in this case at each iteration we use only one call of our oracle which means that the oracle complexity is also $O\left(\sqrt{n^2\rho_nL_2\Theta_p/\varepsilon}\right)$. From this perspective it is clear why \textit{choosing $p=1$ is always better than $p=2$}. That is, for our choice of prox-functions described in the Section~\ref{sec:preliminaries} we have $\Theta_1 = O\left(\|x_0 - x_*\|_1^2\right)$ and $\Theta_2 = \frac{1}{2}\|x_0 - x_*\|_2^2 = O\left(\|x_0 - x_*\|_2^2\right)$. Moreover, if the dimension $n$ is big enough then $\rho_n \overset{\text{eqdef}}{=} \min\{q-1,\,16\ln n - 8\}n^{\frac{2}{q}-1}$ is $\widetilde{O}\left(\frac{1}{n^2}\right)$ for the choice of $p=1$ and $\widetilde{O}\left(\frac{1}{n}\right)$ for the choice of $p=2$. What is more, due to Cauchy-Schwartz inequality $\|x_0-x_*\|_1 \le \sqrt{n}\|x_0-x_*\|_2$. Taking it all into account we obtain that in general situation oracle complexity of ARDFDS with $p=1$ is not worse then for ARDFDS with $p=2$. But if we are lucky and set $x_0$ in such a way that $x_0 - x_*$ is sparse (i.e.\ we can choose $x_0 = 0$ if we minimize the loss for over-parameterized models, since the solution $x_*$ is sparse), we will have $\|x_0-x_*\|_1 = O(\|x_0-x_*\|_2)$ and $\widetilde{O}\left(\sqrt{nL_2\Theta_2/\varepsilon}\right)$ oracle complexity for the choice of $p=1$ which is $\sqrt{n}$ better then $\widetilde{O}\left(n\sqrt{L_2\Theta_2/\varepsilon}\right)$. 

\textit{Now it is clearly why it was so important to get a tight upper bound of $\mathbb E_{e_{k+1}}\left[\|\langle\nabla f(x_{k+1}),e_{k+1}\rangle e_{k+1}\|_q^2\right]$ by $\|\nabla f(x_{k+1})\|_2^2$, since it is the only place in our proof where the constant $\rho_n$ appears and $\rho_n$ is, actually, the only thing that improves the dimension-dependence in the oracle complexity.}

In the Table~\ref{tbl:ARDFDS_params} we give the appropriate choice of the ARDFDS parameters $N$, $m$, $t$ and accuracy of the function values evaluation $\Delta$ where $\alpha_{k+1} = \frac{k+2}{96n^2\rho_n L_2}$ (see Corollary~\ref{cor:ardfds_conv} in the appendix). The last row represents the total number $Nm$ of function evaluations, which was advertised in \eqref{eq:ARDFDSComplInformal}.

{
{\small
\begin{table*}[t]
\centering
    \centering
    \begin{tabular}{|c|c|c|}
        \hline
         & $p=1$ & $p=2$ \\
        \hline
        $N$ & $O\left(\sqrt{\frac{ n\ln nL_2 \Theta_1}{\varepsilon}}\right)$ & $O\left(\sqrt{\frac{ n^{2}L_2\Theta_2}{\varepsilon}}\right)$ \\
			
        \hline
        $m$ & $O\left(\max\left\{1, \frac{\sigma^2}{\varepsilon^{3/2}}\cdot\sqrt{\frac{\Theta_1\ln n}{nL_2}}\right\}\right)$ & $O\left(\max\left\{1,\frac{\sigma^2}{\varepsilon^{3/2}}\cdot\sqrt{\frac{\Theta_2}{L_2}}\right\}\right)$ \\ 
		\hline	
        $\Delta$ & $O\left(\min\left\{\frac{\varepsilon^{3/2}}{\sqrt{ L_2\Theta_1n\ln n}},\, \frac{\varepsilon^2}{nL_2\Theta_1}\right\}\right)$ & $O\left(\min\left\{\frac{\varepsilon^{3/2}}{n\sqrt{L_2\Theta_2}},\,  \frac{\varepsilon^2}{nL_2\Theta_2}\right\}\right)$\\
        \hline	
		$t$ & $O\left(\min\left\{\frac{\varepsilon^{3/4}}{\sqrt[4]{L_2^3\Theta_1 n\ln n}},\, \frac{\varepsilon}{L_2\sqrt{n\Theta_1}}\right\}\right)$ & $O\left(\min\left\{\frac{\varepsilon^{3/4}}{\sqrt[4]{n^2L_2^3\Theta_2}},\,  \frac{\varepsilon}{L_2\sqrt{n\Theta_2}}\right\}\right)$ \\
        \hline
        Func. eval-s & $O\left(\max\left\{\sqrt{\frac{ n\ln nL_2\Theta_1}{\varepsilon}}, \frac{\sigma^2\Theta_1 \ln n}{\varepsilon^2}\right\}\right)$ & $O\left(\max\left\{\sqrt{\frac{ n^{2}L_2\Theta_2}{\varepsilon}}, \frac{\sigma^2\Theta_2 n}{\varepsilon^2}\right\}\right)$\\
			
        \hline
    \end{tabular}
    \caption{Algorithm~\ref{Alg:ARDFDS} parameters for the cases $p=1$ and $p=2$.\label{tbl:ARDFDS_params}}
\end{table*}}
}
Our Randomized Derivative-Free Directional Search (RDFDS) method is listed as Algorithm~\ref{Alg:RDFDS}.
\begin{algorithm}
	\caption{Randomized Derivative-Free Directional Search (RDFDS)}
	\label{Alg:RDFDS
	\begin{algorithmic}[1]
		\REQUIRE $x_0$~---starting point; $N$~--- number of iterations; $m \geq 1$~--- batch size; $t > 0$~--- smoothing parameter, $\alpha$~--- stepsize.
		\ENSURE point $\bar{x}_N$.
		\FOR{$k=0,\, \dots, \, N-1$.}
	
		\STATE Generate $e_{k+1} \in RS_2\left( 1 \right)$ independently from previous iterations and $\xi_i$, $i=1,...,m$ -- independent realizations of $\xi$. 
		\STATE $x_{k+1} \leftarrow \mathop{\arg\!\min}\limits_{x\in\mathbb R^n} \Big\{ \alpha n \left\langle \widetilde{\nabla}^m f^t(x_{k}), \, x-x_k \right\rangle+V[x_{k}] \left( x \right)\Big\}.$
		\STATE Calculate $\widetilde{\nabla}^m f^t(x_{k+1})$ given in \eqref{eq:FinDiffStocGrad}.
		\ENDFOR
		\RETURN$\bar{x}_N \leftarrow \frac{1}{N}\sum\limits_{k=0}^{N-1}x_k$
	\end{algorithmic}
\end{algorithm}

\begin{theorem}\label{theorem_convergence_mini_gr_free_non_acc}
    Let RDFDS with $\alpha = \frac{1}{48n\rho_n L_2}$ be applied to solve problem \eqref{eq:PrSt} and $\Theta_p = V[x_0](x^*)$. Then
    \begin{equation}\label{theo_main_result_mini_gr_free}
        \begin{array}{rcll}
        \mathbb E[f(\bar{x}_N)] - f(x_*) &\leqslant& \frac{384n\rho_nL_2\Theta_p}{N} + \frac{2\sigma^2}{L_2m} + \left(\frac{n}{6L_2} + \frac{N}{3L_2\rho_n} \right)\left(\frac{L_2^2t^2}{2} + \frac{8\Delta^2}{t^2}\right)& \\
        &&\quad+ \frac{8\sqrt{2n\Theta_p}}{N}\left(\frac{L_2t}{2} + \frac{2\Delta}{t}\right),& \forall n\ge 8.
    \end{array}
    \end{equation}
   
\end{theorem}
The proof of the theorem and appropriate choice of parameters are given in the Section~\ref{sec:rdfds} in the Appendix. The proof of this result is pretty similar to the proof of the Theorem~\ref{Th:ARDFDSConv}.

\subsection{Role of the parameters $\Delta, t$ and $\sigma^2$ in our algorithms}\label{sec:role_of_parameters}
\textbf{Role of $\Delta$ and $t$.} We want to mention that \textit{there is no need to know the noise level $\Delta$ to run our algorithms}. As it can be seen from \eqref{eq:ARDFDSConv}, the ARDFDS method is robust (in the sense of \cite{nemirovski2009robust}) to the choice of the smoothing parameter $t$. Namely, if we under/overestimate $t$ by a constant factor, the corresponding terms in the convergence rate will increase only by a constant factor. 

Our Theorems~\ref{Th:ARDFDSConv} and \ref{theorem_convergence_mini_gr_free_non_acc} are applicable in two situations, the noise being a) controlled and b) uncontrolled.
\begin{enumerate}
    \item[a)] Our assumptions on the noise level in Tables~\ref{tbl:ARDFDS_params} and \ref{tbl:RDFDS_params} can be met in practice. For example, in \cite{bogolubsky2016learning}, the objective function is defined by some auxiliary problem and its value can be calculated with accuracy $\Delta$ at the cost proportional to $\ln\frac{1}{\Delta}$, which would result in only a $\ln\frac{1}{\varepsilon}$ factor in the total complexity of our methods in this paper. 
    \item[b)] One can minimize the r.h.s. of \eqref{eq:ARDFDSConv} and \eqref{theo_main_result_mini_gr_free} in $N$ and obtain the minimum possible accuracy of the solution. This minimum error can not be arbitrarily small. But, this is reasonable: one can not solve the problem with better accuracy than the accuracy of the available information. And here we come to the curious phenomenon: it can happen that accelerated algorithm works worse than a non-accelerated one. This is not an issue of our approach, but a general drawback of accelerated methods, as for the full gradient methods (see \cite{devolder2014first}) it was shown that accelerated gradient method accumulates the error and it was proved that it is impossible to have acceleration without error accumulation. We conjecture that, in our derivative-free setting, it is impossible to have slower rate of error accumulation than we have.
\end{enumerate}

\textbf{Role of $\sigma^2$.} Although, all the related works, which we are aware of, assume $\sigma^2$ to be known, adaptivity to the variance $\sigma^2$ is a very important direction of future work.

\section{Experiments}\label{sec:experiments_main}

We considered the following artificial problem in order to emphasize main theoretical contributions of our work. That is, we construct the objective function such that it reveals advantages of using $\|\cdot\|_1$-norm setup instead of classical Euclidean setup.

The objective function is Nesterov's function $f(x) = \frac{L_2}{4}\left(\frac{1}{2}\left[(x^1)^2 + \sum\limits_{i=1}^{n-1}(x^i - x^{i+1})^2 + (x^n)^2\right] - x^1\right)$ from \cite{nesterov2004introduction}, where $x^i$ is $i$-th component of vector $x\in\mathbb R^n$. It is well-known and one can easily check that $f$ is convex, $L_2$-smooth w.r.t.\ $\|\cdot\|_2$-norm and the objective attains its minimal value $f_* = \frac{L_2}{8}\left(-1+\frac{1}{n+1}\right)$ at the point $x_* = (x_*^1,\ldots,x_*^n)^\top$ such that $x_*^i = 1 - \frac{i}{n+1}$. We introduce stochastic noise in our function and consider $F(x,\xi) = f(x) + \xi\langle a, x\rangle$, where $\xi$ is Gaussian random variable with mean $\mu = 0$ and variance $\sigma^2$ and $a\in\mathbb R^n$ is some vector from the unit Euclidean sphere ($\|a\|_2^2 = 1$). It implies that $f(x) = \mathbb E_\xi\left[F(x,\xi)\right]$ and $F(x,\xi)$ is $L_2$-smooth in $x$ w.r.t.\ $\|\cdot\|_2$-norm since $g(x,\xi) - g(y,\xi) = \nabla f(x) - \nabla f(y)$. Moreover, $\mathbb E_\xi g(x,\xi) = \nabla f(x)$ and $\mathbb E_\xi\left[\|g(x,\xi) - \nabla f(x)\|_2^2\right] = \mathbb E_\xi\left[\|\xi a\|_2^2\right] = \|a\|_2^2\mathbb E_\xi\left[\xi^2\right] = \sigma^2$ for all $x\in\mathbb R^n$. Also we introduce the additive noise $\eta(x) = \Delta\sin\|x\|_2$. It is clear that $|\eta(x)| \le \Delta$ for all $x\in\mathbb R^n$. Thus, we consider the situation when oracle takes as an input a pair of points $(x,y)$ and returns noisy stochastic realizations $(\widetilde{f}(x,\xi),\widetilde{f}(y,\xi))$, where $\widetilde{f}(x,\xi) = F(x,\xi) + \eta(x) = f(x) + \xi\langle a,x\rangle + \Delta\sin\|x\|_2.$ We compare our methods in Euclidean and $\|\cdot\|_1$-norm proximal setup with the RSPGF method from \cite{ghadimi2013stochastic} on the problem defined above. As the theory predicts our methods with $p=1$ work better and better when $n$ is growing (see Figure~\ref{fig:nesterov}). See other details of implementation in Section~\ref{sec:experiments} of the Appendix. One can find our code here: \url{https://github.com/eduardgorbunov/ardfds}.

\begin{figure}
\centering 
\includegraphics[scale=0.20]{nesterov_100dim_50000it_8_2000_32_1000_stepsz_compare_methods_iterations.pdf}
\includegraphics[scale=0.20]{nesterov_1000dim_500000it_8_2000_32_1000_stepsz_compare_methods_iterations.pdf}
\includegraphics[scale=0.20]{nesterov_5000dim_1000000it_8_2000_32_1000_stepsz_compare_methods_iterations.pdf}
\includegraphics[scale=0.20]{nesterov_100dim_50000it_8_2000_32_1000_stepsz_compare_methods_runtimes.pdf}
\includegraphics[scale=0.20]{nesterov_1000dim_500000it_8_2000_32_1000_stepsz_compare_methods_runtimes.pdf}
\includegraphics[scale=0.20]{nesterov_5000dim_1000000it_8_2000_32_1000_stepsz_compare_methods_runtimes.pdf}
\caption{The $y$-axis corresponds to the relative accuracy $\frac{f(x_k) - f(x_*)}{f(x_0)-f(x_*)}$ in both rows. In the first row the $x$-axis shows number of iterations and in the second row~--- time in seconds. We use \_E and \_NE in the graphs to define $\|\cdot\|_2$-norm and $\|\cdot\|_1$-norm proximal setups relatively and RSPGF to define the corresponding method from \cite{ghadimi2013stochastic}.}
\label{fig:nesterov}
\end{figure}


\section{Conclusion}
In this paper, we propose two new algorithms for smooth stochastic derivative-free optimization with two-point feedback inexact oracle. Our first algorithm is an accelerated one and the second one is a non-accelerated one. Interestingly, despite the traditional choice of $\ell_2$-norm proximal setup for unconstrained optimization problems, our analysis shows that the method with $\ell_1$-norm proximal setup has better complexity.

\acks{The work of of E.~Gorbunov was supported by the grant of Russian's President MD-1320.2018.1, the work of P.~Dvurechensky was supported by RFBR 18-29-03071 mk and the work of A.~Gasnikov was supported by RFBR 18-31-20005 mol-a-ved and was prepared within the framework of the HSE University Basic Research Program and funded by the Russian Academic Excellence Project '5-100'.}

