On low complexity Acceleration Techniques for Randomized Optimization:\\ Supplementary Online Material
Sebastian U. Stich\footnotemark[1]
Recently it was shown by Nesterov (2011) that techniques form convex optimization can be used to successfully accelerate simple derivative-free randomized optimization methods. The appeal of those schemes lies in their low complexity, which is only $\Theta(n)$ per iteration---compared to $\Theta(n^2)$ for algorithms storing second-order information or covariance matrices. From a high-level point of view, those accelerated schemes employ correlations between successive iterates---a concept looking similar to the evolution path used in Covariance Matrix Adaptation Evolution Strategies (CMA-ES). In this contribution, we (i) implement and empirically test a simple accelerated random search scheme (SARP). Our study is the first to provide numerical evidence that SARP can effectively be implemented with adaptive step size control and does not require access to gradient or advanced line search oracles. We (ii) try to empirically verify the supposed analogy between the evolution path and SARP. We propose an algorithm CMA-EP that uses only the evolution path to bias the search. %By storing $m$ directions, This algorithm can be generalized to a family of low memory schemes, with complexity $\Theta(mn)$ per iteration, following a recent approach by Loshchilov (2014). The study shows that the performance of CMA-EP heavily depends on the spectra of the objective function and thus it cannot accelerate as consistently as SARP. ~\\ \noindentKeywords: gradient-free optimization, accelerated random search, evolution path, adaptive step size, Covariance Matrix Adaptation, spectra
ppsn14_supplement.tex
%
%\documentclass{llncs}
\documentclass[10pt, a4paper]{article}
\usepackage{graphicx}
\usepackage{amssymb, amsmath}
\usepackage{color}
\usepackage{grffile}
\usepackage{epstopdf}
\epstopdfsetup{update}
\usepackage[ruled,vlined]{algorithm2e}
\setlength{\algomargin}{1em}
\usepackage[lofdepth,lotdepth, caption=false]{subfig}
\usepackage{dsfont}
%\usepackage[all=normal,floats]{savetrees}

\newcommand*\red{\color{red}}

\newcommand{\R}{\mathds{R}}
\newcommand{\id}{\mathds{1}}
\newcommand{\0}{\mathbf{0}}
\newcommand{\1}{\mathbf{1}}
\newcommand{\E}{\ensuremath{\mathbb{E}}}
\newcommand{\Tr}{\mathrm{Tr}}

\newcommand{\x}{\mathbf{x}}
\newcommand{\xx}{\mathbf{x}}
\newcommand{\yy}{\mathbf{y}}
\newcommand{\pp}{\mathbf{p}}
\newcommand{\m}{\mathbf{m}}
\newcommand{\uu}{\mathbf{u}}
\newcommand{\vv}{\mathbf{v}}
\newcommand{\Ss}{\mathbb{S}}

\DeclareMathOperator*{\argmin}{arg\,min}

%save space
%\renewcommand\floatpagefraction{.9}
%\renewcommand\topfraction{.9}
%\renewcommand\bottomfraction{.9}
%\renewcommand\textfraction{.1}   
%\addtolength{\intextsep}{-4ex}
%\addtolength{\intextsep}{-3ex}
\addtolength{\abovecaptionskip}{-1ex}
%\addtolength{\floatsep }{-1.5ex}
\addtolength{\belowcaptionskip }{-2ex}
%\addtolength{\textfloatsep}{-2ex}
 \setlength{\abovedisplayskip}{0pt}
 \setlength{\abovedisplayshortskip}{0pt}
 \setlength{\belowdisplayskip}{0pt}
 \setlength{\belowdisplayshortskip}{0pt}
%

%
%\usepackage{makeidx}  % allows for indexgeneration
%
\begin{document}
%
\pagestyle{headings}  % switches on printing of running heads

%
\title{On low complexity Acceleration Techniques for Randomized Optimization:\\ Supplementary Online Material}
%
%\titlerunning{On low complexity Acceleration Techniques for Randomized Optimization}  % abbreviated title (for running head)
%                                     also used for the TOC unless
%                                     \toctitle is used
%
\author{Sebastian U. Stich\footnotemark[1]}
%
%\authorrunning{Stich} % abbreviated author list (for running head)
%
%%%% list of authors for the TOC (use if author list has to be modified)
%\tocauthor{Sebastian Urban Stich}
%
%\institute{
%\affiliation{
%Institute of Theoretical Computer Science\\ 
%ETH Z\"{u}rich, CH 8092 Z\"{u}rich, Switzerland\\
%\texttt{sstich@inf.ethz.ch}
%\\ WWW home page:
%\texttt{http://users/\homedir iekeland/web/welcome.html}
%}

\date{June 8, 2014}

\renewcommand{\thefootnote}{\fnsymbol{footnote}}
\footnotetext[1]{Institute of Theoretical Computer Science, ETH Z\"{u}rich, \texttt{sstich@inf.ethz.ch}}
\renewcommand{\thefootnote}{\arabic{footnote}}

\maketitle              % typeset the title of the contribution
%
%
%
\begin{abstract}
Recently it was shown by Nesterov (2011) that techniques form convex optimization can be used to successfully accelerate simple derivative-free randomized optimization methods. The appeal of those schemes lies in their low complexity, which is only $\Theta(n)$ per iteration---compared to $\Theta(n^2)$ for algorithms storing second-order information or covariance matrices. From a high-level point of view, those accelerated schemes employ correlations between successive iterates---a concept looking similar to the evolution path used in Covariance Matrix Adaptation Evolution Strategies (CMA-ES). 
In this contribution, we (i) implement and empirically test a simple accelerated random search scheme (SARP). Our study is the first to provide numerical evidence that SARP can effectively be implemented with adaptive step size control and does not require access to gradient or advanced line search oracles. We (ii) try to empirically verify the supposed analogy between the evolution path and SARP. We propose an algorithm CMA-EP that uses only the evolution path to bias the search. 
%By storing $m$ directions, 
This algorithm can be generalized to a family of low memory schemes, with complexity $\Theta(mn)$ per iteration, following a recent approach by Loshchilov (2014).  The study shows that the performance of CMA-EP heavily depends on the spectra of the objective function and thus it cannot accelerate as consistently as SARP.
~\\
\noindent\textbf{Keywords:} gradient-free optimization, accelerated random search, evolution path, adaptive step size, Covariance Matrix Adaptation, spectra
\end{abstract}
%
\section{Introduction}
\label{sec:intro}
The Gradient Method~\cite{Polyak:1987,Nesterov:2004}---one of the most fundamental schemes in convex optimization---has iteration complexity $\Theta(n)$, where $n$ is the dimension. On strongly convex functions its convergence rate is linear, depending only on the condition number of the objective function. To overcome the difficulty imposed by ill-conditioned problems, second-order methods like Newton's method % However, obtaining second-order information turns out to be very costly or infeasible on many practical problems. 
or first order Quasi-Newton methods such as the BFGS scheme \cite{Broyden:1970,Fletcher:1970,Goldfarb:1970,Shanno:1970} are a welcome alternative.
% build a quadratic model of the objective function only based on first-order information. 
%Those schemes have become are ubiquitous in all areas of science and engineering. 
Those schemes maintain a quadratic model of the objective function and their complexity is bounded by $\Omega(n^2)$. Limited memory schemes like L-BFGS~\cite{Nocedal:1980,Liu:1989fu} trade-off linear iteration complexity $\Theta(mn)$ (where $m$ is a fixed parameter), versus convergence rate. Accelerated versions of the Gradient Method have linear complexity $\Theta(n)$ per iteration and converge with optimal rate among all first-order methods. On strongly convex problems the convergence rate is  proportional to the square root of the condition number~\cite{Polyak:1987,Nesterov:2004,nesterov:1983,nesterov:2007,tseng:2008}.

Randomized (gradient-free) schemes do not require first-order information, they operate by only querying function values. Such schemes are nowadays a ubiquitous tool for solving many practical problems in science and engineering  where first-order information is difficult to compute or does not exist. Among the first proposed schemes that are still of considerable (theoretical) importance are Adaptive Step Size Random Search (aSS)~\cite{Schumer:1968} and the (almost identical) well-known (1+1)-Evolution Strategy (ES)~\cite{Rechenberg:1973} in Evolutionary Computation. More recent schemes comprise Random Pursuit (RP)~\cite{Mutseniyeks:1964,Stich:2011}, or Random Gradient Descent~\cite{nesterov:11}. Those schemes can be viewed as generalizations of the Gradient Method to zeroth-order, with iteration complexity $\Theta(n)$. Likewise, analogues of the second-order schemes try to estimate an approximation of the Hessian by finite difference computations~\cite{Leventhal:2011,Stich:2012} or by estimating correlations among search directions. A very popular algorithm of this kind is the Covariance Matrix Adaptation Evolution Strategy (CMA-ES)~\cite{Hansen:2001,Hansen:2003}. Limited memory variants have been proposed in~\cite{knight:2007,Loshchilov:2014}, with iteration complexity $\Theta(mn)$. Especially, the later variant due to Loshchilov shows excellent convergence in high dimensions also for small values of $m$.
Only recently, also zeroth-order analogues of the accelerated gradient schemes have been introduced~\cite{nesterov:11,lee:2013}. Those schemes massively outperform the simple random search schemes on convex problems. This performance gain does not come for free, as those schemes require valid bounds on the condition number as input parameters. However, their low iteration complexity of $\Theta(n)$ could make them a promising choice for large scale problems, where the fully-quadratic schemes inherently fail. We focus on a very simple accelerated random search scheme, which we call SARP. 

By inspecting closely the accelerated search schemes, one could conclude that the difference to the classical schemes can be explained by an additional ``drift'' term~\cite{Nesterov:2004,tseng:2008} that takes into account correlations of the last iterates. In Evolutionary Computation, correlations between successive iterates are often expressed by the evolution path~\cite{ostermeier1994}. In the popular CMA-ES, the evolution path accumulates the correlation of successive iterates over a finite horizon of the order of $\Theta(n)$ steps~\cite{Hansen:2001}. In this work, we are interested, if the evolution path can be used for acceleration, competitive to the accelerated zeroth-order schemes form convex optimization. 
To this end, we introduce a variant of CMA-ES, called EP-CMA, that only uses the information stored in the evolution path to bias the direction of the search.
Similar to the approach proposed in~\cite{Loshchilov:2014}, this scheme can be generalized to a family of schemes, which we call EP-CMA-$m$. 
 For $m=1$, the approach is similar to~\cite{sun:2013}. For $m>1$, the schemes are similar in spirit to the LM-BFGS schemes~\cite{Nocedal:1980,Liu:1989fu} and for $m$ large enough, the scheme is intended to approach $(1+1)$-CMA-ES. EM-CMA-$m$ could be implemented with $\Theta(mn)$ complexity per iteration (not presented here), although this was not required, as the dimension $n \leq 100$ in our empirical study.

The remainder of this paper is structured as follows. In Section~\ref{sec:algorithms} we present the accelerated random search scheme SARP and detail the EP-CMA-$m$ schemes. In Section~\ref{sec:emp} we empirically test the performance of all schemes on three quadratic and the non-convex Rosenbrock function, and highlight the key results. We discuss these results and conclude the paper in Section~\ref{sec:disc}.
%
%
%%plots line serach
\def \sfactor{0.83}
%
 \begin{figure}[!ht]
    \centering
    \null\hfill %
     \scalebox{\sfactor}{
     \begin{minipage}[t]{.55\linewidth}
     \vspace{0pt}
     \begin{algorithm}[H]
     %
     % lineS
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl \lIf{exact}{
          \Return $\mathrm{exactLS}(\x, \uu/\|\uu\|)$
        }
        \nl \lElse{
           \Return $\mathrm{aSS}(\x, \uu, \sigma, p)$
        }%
      \caption{{lineSearch$(\x, \uu, [\sigma, p])$}}%
     \end{algorithm}%
     \begin{algorithm}[H]
     %
     % lineS
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl $\sigma_+ \leftarrow \min_{\lambda}f(\x+\lambda \uu);\, \x_+ \leftarrow \x + \sigma_+ \uu$\;
      \nl \Return{$(\x_+, \sigma_+)$}
      \caption{{exactLS$(\x, \uu, [\sigma])$}}%
     \end{algorithm}%
    \end{minipage}%
    }% end scalebox
    %
    \hfill
    %
    \scalebox{\sfactor}{
    \begin{minipage}[t]{.55\linewidth}
    \vspace{0pt}
        \begin{algorithm}[H]
     %
     % aS
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl \eIf{ $f(\x+ \sigma \uu) \leq f(\x)$}{
      \nl $\x_+ \leftarrow \x + \sigma \uu$; $\sigma_+ \leftarrow \sigma \cdot \exp(1/3)$\;
      %\nl $\sigma_+ \leftarrow \sigma \cdot \exp(1/3)$\;
      }{
      \nl $\x_+ \leftarrow \x$; $\sigma_+ \leftarrow \sigma \cdot \exp\left(-\frac{p}{3(1-p)}\right)$\;
      %\nl $\sigma_+ \leftarrow \sigma \cdot \exp\left(-\frac{0.27}{3(1-0.27)}\right)$\;
      }%
      \nl \Return{$(\x_+, \sigma_+)$}\;%
      \caption{{aSS$(\x, \uu, \sigma,p)$ \textit{(adaptive step size)}}}%
     \end{algorithm}%
    \end{minipage}%
    }%end scalebox
     \hfill\null\par
    \caption{Line search oracles for gradient-free optimization.}\label{algo:ls}
   \end{figure} 
%
%
%
% 
   %
%\begin{figure}
%\centering
%\includegraphics[width=.99\textwidth]{matlab/plots/felli4.dim_100.all=0.ls=1-eps-converted-to.pdf}
%
%\includegraphics[width=.99\textwidth]{matlab/plots/felli4lin.dim_100.all=0.ls=1-eps-converted-to.pdf}   
%\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ (top) and $f_{\rm lin}$ (bottom) with $L=1$\textsc{e}4 in $n=100$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
%\label{fig:4withls}
%\end{figure}  
%
%\begin{figure}
%\centering
%\includegraphics[width=.99\textwidth]{matlab/plots/felli4.dim_100.all=1.ls=0-eps-converted-to.pdf}
%
%\includegraphics[width=.99\textwidth]{matlab/plots/felli4lin.dim_100.all=1.ls=0-eps-converted-to.pdf}   
%\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ (top) and $f_{\rm lin}$ (bottom) with $L=1$\textsc{e}4 in $n=100$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
%\label{fig:4all}
%\end{figure}  
%
%
%\begin{figure}
%\centering
%\includegraphics[width=.99\textwidth]{matlab/plots/felli6.dim_100.all=1.ls=0-eps-converted-to.pdf}
%
%\includegraphics[width=.99\textwidth]{matlab/plots/felli6lin.dim_100.all=1.ls=0-eps-converted-to.pdf}   
%\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ (top) and $f_{\rm lin}$ (bottom) with $L=1$\textsc{e}6 in $n=100$ dimensions. For 51 runs (11 for RP) we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9 (1\textsc{e}-2 for RP). The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
%\label{fig:6all}
%\end{figure}  
%
\section{Algorithms}
\label{sec:algorithms}
We here present the optimization schemes considered in this study. We first detail two Random Pursuit algorithms and a simplistic variant of a standard (1+1)-CMA-ES. Then we introduce the new EP-CMA-$m$ schemes.

\textbf{RP.}
Random Pursuit is a basic optimization scheme that iteratively generates a sequence of approximate solutions to the global optimization problem $\min_{\xx \in \R^n} f(\xx)$. In each step a search direction is drawn $\uu_k \sim \mathcal{N}(0, \id_n)$. In Random Pursuit with exact line search (RP-exact), first proposed in~\cite{Mutseniyeks:1964} and analyzed in~\cite{Stich:2011}, the step size $\sigma$ is determined by minimizing the objective function in direction $\uu$, i.e. $\sigma = \argmin_\lambda f(\x + \lambda \uu)$. For quadratic functions $f(\x):= \frac{1}{2} \x^T A \x$ with Hessian $A$, the expected one-step progress can be estimated as:
\begin{align}
\label{eq:generalRP}
\E \left[f(\x_+) \mid \x\right] \leq \left(1- \lambda_{\rm min}(A)/\Tr[A] \right) f(\x) \,,
\end{align}
where $\x$ is the current iterate, $\x_+ :=\x + \sigma \uu$ denotes the next iterate. This statement can also be generalized to arbitrary smooth convex functions \cite{Stich:2011}.
Stich et al. \cite{Stich:2011} show that RP-exact still converges if the line search is not performed exactly, but allowing relative errors. Therefore, we also consider Random Pursuit with adaptive step sizes (RP) instead of exact line search. In RP the step size is dynamically controlled such as to approximately guarantee a certain probability $p$ of finding an improving iterate. Depending on %the specific random direction generator and 
the underlying test function, different optimality conditions can be formulated for the value $p$. 
Schumer and Steiglitz~\cite{Schumer:1968} suggest the setting $p = 0.27$ which is considered throughout this work. We use immediate exponential step size control as explicitly formulated in the aSS sub-routine in Fig.~\ref{algo:ls}. RP is identical to the well known (1+1)-ES. 
%J\"agersk\"upper \cite{Jaegerskuepper:2005} showed that, for quadratic functions, the expected one-step progress of the (1+1)-ES is almost identical to the one shown in Eq.~(\ref{eq:generalRP}).    
%   
%   %
%
\def \sfactor{0.83}
%
 \begin{figure}[!ht]
    \centering
    \null\hfill % 
     \scalebox{\sfactor}{
     \begin{minipage}[t]{.55\linewidth}
     \vspace{0pt} 
     \begin{algorithm}[H]
     %
     % RP
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl \For{$k=1$ to $N$}{
       \nl $\uu_{k} \sim \mathcal{N}(0, I_n)$\;
      \nl $(\x_{k},\sigma_k) \leftarrow \mathrm{lineSearch}(\x_{k-1}, \uu_{k}, [\sigma_{k-1}])$\;
      }%
      \nl \Return{$\x_N$}\;
      \caption{{RP$(\x_0, N, [\sigma_0 , p])$}}%
     \end{algorithm}%  
      \begin{algorithm}[H]
     %
     % CMA Revolver
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl $\hat{\pp}_0 \leftarrow \mathbf{0};\, \hat{\pp}_{1},\dots,\hat{\pp}_{m-1} \leftarrow \mathbf{0};\,q=0$\;
      \nl \For{$k=1$ to $N$}{
       \nl $C_k \leftarrow I_n$\; 
       \nl \For{$i =1$ to $m-1$}{$C_{k} \leftarrow (1-c_{\rm{cov}}) C_{k}+c_{\rm{cov}} \hat{\pp}_{i} \hat{\pp}_{i}^T$}
       \nl $C_{k} \leftarrow (1-c_{\rm{cov}}) C_{k}+c_{\rm{cov}} \pp_{k-1} \pp_{k-1}^T$
       \nl $\uu_{k} \sim \mathcal{N}(0, C_{k})$\;     
       \nl $(\x_k, \sigma_k) \leftarrow \mathrm{aSS}(\x_{k-1}, \uu_k, \sigma_{k-1})$\;
       \nl $\yy_{k} \leftarrow (\xx_k - \xx_{k-1})/\sigma_{k-1}$\;
       \nl \uIf{$\yy_{k} \neq \mathbf{0}$ (\textit{success})}{
       \nl $\pp_{k} \leftarrow (1-c_{\rm{c}}) \pp_{k-1} + \sqrt{c_{\rm{c}}(2-c_{\rm{c}})}\yy_{k}$ \;
       }
       \nl \lElse{$\pp_{k} \leftarrow (1-c_{\rm p})\pp_{k-1}$}
       \nl \If{$k > q + n^2/m$}{
       $\hat{\pp}_{1} \leftarrow \hat{\pp}_{2}, \dots, \hat{\pp}_{m-2} \leftarrow \hat{\pp}_{m-1};\, q=k$
       }
      }%
      \nl \Return{$\x_N$}\;
      \caption{{EP-CMA-$m$ $(\x_0, N, \sigma_0, p, c_{\rm{c}}, c_{\rm{cov}})$}}%
     \end{algorithm}%       
    \end{minipage}%
    }% end scalebox
    %
    \hfill
    %
    \scalebox{\sfactor}{
    \begin{minipage}[t]{.55\linewidth}
    \vspace{0pt} 
      \begin{algorithm}[H]
     %
     % updateH
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
      %
      \nl $\yy_0 \leftarrow \x_0;\, \vv_0 \leftarrow x_0;\, \theta \leftarrow \sqrt{\frac{m}{2 n^2 L}}$\;
      \nl \For{$k=1$ to $N$}{
       \nl $\uu_{k} \sim \mathcal{N}(0, I_n)$\;
      \nl $(\x_{k},\sigma_k) \leftarrow \mathrm{lineSearch}(\yy_{k-1}, \uu_{k}, [\sigma_{k-1}])$\;
       \nl $\yy_k \leftarrow (\theta \vv_{k-1} + \x_{k}) /(1+\theta)$\;
       \nl $\vv_k \leftarrow (1-\theta)\vv_{k-1} + \theta \yy_{k} + \theta n \frac{L}{m} \sigma_k \uu_{k}$\;
      }%
      \nl \Return{$\x_N$}\;
      \caption{{SARP$(\x_0, N, m, L, [\sigma_0])$}}%
     \end{algorithm}%
     \vspace{0.255cm}
     \begin{algorithm}[H]
     %
     % CMA
     %
      \NoCaptionOfAlgo
      \SetAlgoVlined
      \DontPrintSemicolon
      \SetInd{0.1em}{0.5em}
     %
      \nl $C_0 \leftarrow I_n;\, \pp_0 \leftarrow \mathbf{0}$\;
      \nl \For{$k=1$ to $N$}{
       \nl $\uu_{k} \sim \mathcal{N}(0, C_{k-1})$\;     
       \nl $(\x_k, \sigma_k) \leftarrow \mathrm{aSS}(\x_{k-1}, \uu_k, \sigma_{k-1})$\;
       \nl $\yy_{k} \leftarrow (\xx_k - \xx_{k-1})/\sigma_{k-1}$\;
       \nl \eIf{$\yy_{k} \neq \mathbf{0}$ (\textit{success})}{
       \nl $\pp_{k} \leftarrow (1-c_{\rm{c}}) \pp_{k-1} + \sqrt{c_{\rm{c}}(2-c_{\rm{c}})}\yy_{k}$ \;
       \nl $C_{k} \leftarrow (1-c_{\rm{cov}}) C_{k-1} + c_{\rm{cov}} \pp_{k} \pp_{k}^T$
       }{
       \nl $C_{k} \leftarrow C_{k-1};\, \pp_{k} \leftarrow (1-c_{\rm p})\pp_{k-1}$       
       }
      }%
      \nl \Return{$\x_N$}\;
      \caption{{(1+1)-CMA$(\x_0, N, \sigma_0, p, c_{\rm{c}}, c_{\rm{cov}})$}}%
     \end{algorithm}%      
    \end{minipage}%
    }%end scalebox
     \hfill\null\par
    \caption{RP, EP-CMA and CMA-ES schemes.}\label{algo:all}
   \end{figure}

\textbf{SARP.} Accelerated random search schemes are fundamentally different from the simple random search schemes. Instead of generating only  one sequence of iterates, those algorithms typically maintain two or more sequences simultaneously (here essentially $\xx_k$ and $\yy_k$, see Fig.~\ref{algo:all}). Those sequences allow to store gathered knowledge on the objective function which yields better performance. In Fig.~\ref{algo:all} we present a simple version of the accelerated random search scheme proposed in~\cite{Stich:2011} and refer to it as simple accelerated random search (SARP). Like RP, SARP can (in practice) be used with exact line search oracles or with adaptive step size control, although convergence for those oracles has not been proven yet. For Nesterov's accelerated random search scheme~\cite{nesterov:11}, the expected one-step progress can be estimated as
\begin{align}
\label{eq:SARP}
\E \left[f(\x_+) \mid \x\right] \leq \left(1- (n\sqrt{\kappa})^{-1} \right) f(\x) \,,
\end{align}
where condition $\kappa = L/m$ and the two parameters $m \leq \lambda_{\rm min}(A)$ and $L \geq \lambda_{\rm max}(A)$ are required as input to the algorithm (and always provided in our numerical study). This rate is much better than~\eqref{eq:generalRP} and we hope to see that SARP attains comparable performance. SARP is not a monotone scheme, that is, the function values of the iterates are not monotonically decreasing. SARP is closely related to the first-order accelerated search scheme of Nesterov~\cite{Nesterov:2004}. This scheme also simultaneously maintains two sequences $\xx'_k$ and $\yy'_k$ of iterates (but requires access to the gradient in every iteration). For Nesterov's first-order scheme it is known~\cite[p.79]{Nesterov:2004}
%, see also~\cite{drori2013performance}, 
that the sequence $\yy'_{k}$ obeys
\begin{align}
\begin{split}
\label{eq:ARP}
\yy'_{k+1} %&= \xx'_{k+1} + \left(\left(\yy'_k + \sigma_k \nabla f(\yy'_k) - \xx'_k\right) - \theta \left(\xx'_{k+1} - \xx'_{k}\right)\right)/(1+\theta) \\
 &= \xx'_{k+1} + \beta' \left(\xx'_{k+1}- \xx'_k\right)\,,
\end{split}
\end{align}
for $\beta'=1 - 2/\sqrt{\kappa} + O(1/\kappa)$. Thus the additional $(\xx'_{k+1}-\xx'_{k})$ acts like a drift term, cf.~\cite{Polyak:1987}. For SARP with parameter $\theta'= \sqrt{1/(n^2 \kappa)}$ (only slightly different from $\theta$ in Fig.~\ref{algo:all}) the same reformulation of the update reveals
\begin{align}
\label{eq:sarp}
\yy_{k+1} = \xx_k + \beta \left(\xx_{k+1} - \xx_{k}\right)\,,
\end{align}
for $\beta= (1-\theta')/(1+\theta') = 1 - 2/(n \sqrt{\kappa}) + O(1/\kappa)$.
The main term contributing to the drift is approximately only an $1/n$-fraction of the step, accounting for the uncertainty emerging form the randomness.

\textbf{(1+1)-CMA-ES.}
In contrast to the presented Random Pursuit schemes, in CMA-ES new search points are sampled from a multivariate normal distribution $\uu_k \sim \mathcal{N}(0, C_k)$ whose parameter $C_k$ is updated in each iteration based on the evaluation of the samples. The covariance matrix can be adapted using different rank-1~\cite{Hansen:2001,igel:06} or rank-k updates~\cite{Hansen:2003}. 
In addition, the CMA-ES scheme is augmented by an auxiliary variable called evolution path that takes into account the correlation of successive means taken over a finite horizon. In~\cite{Hansen:2001,ostermeier1994}, the evolution path $\pp_k$ is updated as 
\begin{align}
\label{eq:path}
\pp_{k+1} = (1-c_{\rm c}) \pp_k + \sqrt{c_{\rm c}(1-c_{\rm c})} \uu_k\,.
\end{align}
Cumulative information about successive steps is stored in the variable $\pp_{k}$. %, whereas the new information is damped by a factor of roughly $\sqrt{2}/\sqrt{n}$ for the chosen parameters. 
We use a simplistic CMA-ES variant, closely following~\cite{Hansen:2001}, see Fig.~\ref{algo:all}. We use the simple Adaptive Step Size control aSS to determine the step size $\sigma_k$, the covariance matrix update solely uses the information of the evolution path like in~\cite{igel:06} and for simplicity we refrain from implementing any regularization features, in contrast to~\cite{igel:06}. We use the same parameters that were proposed in~\cite{igel:06} for the (1+1)-CMA-ES, namely $c_{\rm c} = 2/(n+2)$, $c_{\rm p} = 1/12$ and $c_{\rm cov}=2/(n^2+6)$.

\textbf{EP-CMA-1.}
The evolution path $\pp_k$ accumulates information over successful steps. This accumulation can be seen as a smoothing of the noisy information obtained in single steps, at the effect that the evolution path points into direction of more promising function values~\cite{Hansen:2001}. In this study, we are interested if the evolution path can be used in a similar way as the drift term in~\eqref{eq:ARP} or~\eqref{eq:sarp}, respectively, to accelerate the search. There are several ways to incorporate the evolution path $\pp_k$ into the update scheme. We suggest to use the path $\pp_k$ in the following way: in the simple random search scheme RP (equivalent to (1+1)-ES), we sample in iteration $k$ a direction from $\uu_k \sim \mathcal{N}(0, (1-c_{\rm{cov}}) I_n+c_{\rm{cov}} \pp_{k} \pp_{k}^T)$, with bias along the direction indicated by $\pp_k$. This has the effect that we only follow successful steps, but the drift imposed by the evolution path might be smaller than it ideally should be. The scheme EP-CMA-1 is detailed in Fig~\ref{algo:all}, we used $c_{\rm c}$ and $c_{\rm p}$ as above, and $c_{\rm cov} = 1/5$. This approach is similar to~\cite{sun:2013}.

\textbf{EP-CMA-$\mathbf{m}$.} The proposed EP-CMA-1 can easily be generalized to a whole family of optimization schemes by an approach presented in~\cite{Loshchilov:2014}. In EP-CMA-1, only the information stored in the current evolution path $\pp_k$ is used to bias the search direction. But we could also afford to temporarily store a small number $m$ of past $\pp_{k'}$ for $k'<k$, and use the information collectively to bias the search. As two successive evolution paths are likely highly correlated, we propose to store the evolution path only every $n^2\!/m$-th generation (and up to at most $(m-1)$ copies simultaneously). The resulting scheme is detailed in Fig.~\ref{algo:all}. We used $c_{\rm c} = 2/(n+2)$ as in CMA-ES, and for $m > 1$, $c_{\rm cov} = 2/(6+m)$ for EP-CMA-$m$. % for $m>1$. 
If implemented carefully, EP-CMA-$m$ has $\Theta(mn)$ complexity per iteration (not shown in Fig~\ref{algo:all}). For $m = n^2$, the updates of EP-CMA-$m$ are identical to the updates of (1+1)-CMA-ES, if \emph{limited} to a finite horizon of $n^2$ steps. In contrast, the low memory method proposed in~\cite{knight:2007} behaves similar to CMA-ES already for $m=n$, but has iteration complexity $\Theta(nm^2)$.

\section{Empirical Study}
\label{sec:emp}
We now present the setup of our empirical study. We focus on the following schemes: (i) the two Random Pursuit schemes with adaptive step size control (denoted as RP and SARP) and with exact line search (denoted as RP-exact and SARP-exact), (ii) the simplified (1+1)-CMA-ES and (iii) the EP-CMA-$m$ schemes as introduced in Sec.~\ref{sec:algorithms}, see Fig.~\ref{algo:all}. We use EP-CMA-$m$ with parameters $m=1,2,4,\sqrt{n},n$. This totals in 10 different schemes, all of which were implemented in MATLAB and will be made available on the authors website.
%
%
%


We tested the performance of all algorithms on three variants of the ellipsoidal benchmark function~\cite{Hansen:2001} and the non-convex Rosenbrock function, detailed in Table~\ref{tab:functions}. The quadratic functions were chosen in such a way that the extremal values of their spectra (1 and $L$) both agree. 
%The eigenvalues of  $f_{\rm lin}$ are equally spaced between $1$ and $L$, the eigenvalues of $f_{\rm exp}$ exponentially, and the eigenvalues of $f_{\rm two}$ take only two different values. 
We considered the quadratic functions with parameters $L=1$\textsc{e}4 and $L=1$\textsc{e}6 each, and repeat the experiments in dimensions $n=20,40,60,80,100$.

For all experiments, initial settings were $\xx_0 = \mathbf{1}$, $\sigma_0 = 1$ and $p=0.27$ (for schemes with the aSS routine). We count the number of iterations (\#ITS) needed to decrease the function value (FVAL) below $1$\textsc{e}-9.
%
A graphical summary of our results can be found in Fig.~\ref{fig:all}-\ref{fig:rosen} and in the appendix.
%Results not depicted here are reported in the supplementary online material~\cite{Stich:2014s}.
We now proceed by discussing some of the key results.
%
%

\begin{figure}[!t]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_100.all=1.ls=1-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_100.all=1.ls=1-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots/felli4lin.dim_100.all=1.ls=1-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots/felli6lin.dim_100.all=1.ls=1-eps-converted-to.pdf}
\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ (top) and $f_{\rm lin}$ (bottom) with $L=1$\textsc{e}4 (left) and $L=1$\textsc{e}6 (right) in $n=100$ dimensions. For 51 (11 for RP on $f_{\rm lin}$ with $L=1$\textsc{e}6) runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers. (RP on $f_{\rm lin}$ with $L=1$\textsc{e}6 reaching FVAL $<$ 1\textsc{e}-2 after 1\textsc{e}6.5$n$ \#ITS.)}
\label{fig:all}
\end{figure} 
%
%
\begin{table}[b]
  \centering
	\caption{List of benchmark functions.}
	\hrule
	\vspace{-.5em}
    \begin{align*}
    &f_{\rm exp} (\xx) = \frac{1}{2} \sum_{i=1}^n L^\frac{i-1}{n-1} x_i^2 &
    f_{\rm rosen} (\xx) = \sum_{i=1}^{n-1} \left(100 \cdot \left(x_i^2 - x_{i+1} \right)^2 + \left(x_i - 1\right)^2 \right)& \\
    &f_{\rm lin} (\xx) = \frac{1}{2} \sum_{i=1}^n \left(1 + i\frac{(L-1)}{(n-1)}\right) x_i^2     &
    f_{\rm two} (\xx) = \frac{1}{2}\sum_{i=1}^{\lfloor n/2 \rfloor} x_i^2 + \frac{L}{2} \sum_{i = \lceil n/2 \rceil} x_i^2&
    \end{align*}
    \vspace{-.5em}
 	\hrule
 \label{tab:functions}
\end{table}   
%
%
%
\textbf{Line search.} 
Both RP and SARP were tested with exact line search oracle and adaptive step size control. In Fig.~\ref{fig:all} we see that the exact schemes outperform their adaptive variants in $n=100$ dimensions by a factor of roughly 2-3. This pattern is observed throughout the whole benchmark in all dimensions. Thus we omit to display the results for exact line search in subsequent Figs.~\ref{fig:dimfexp}-\ref{fig:rosen}. 
%RP showed consistently better performance on $f_{\rm exp}$ than on $f_{\rm lin}$ (Fig.~\ref{fig:all}-\ref{fig:dimfexp}), in agreement to~\eqref{eq:generalRP}. SARP performed slightly better on $f_{\rm lin}$ than on $f_{\rm exp}$ (Fig.~\ref{fig:all}-\ref{fig:dimfexp}).

\textbf{SARP vs. EP-CMA-1.} 
The picture is twofold. In Fig.~\ref{fig:all} we see 
that EP-CMA-1 outperforms SARP by a factor of roughly 5 on $f_{\rm lin}$ with $L=1$\textsc{e}4 (factor 24 for $L=1$\textsc{e}6). The smallest eigenvalue of this function is separated form the second largest by a gap of roughly $n$. Hence, knowledge of one important direction reduces the conditioning of the function by a large factor. This factor becomes smaller in higher dimension. This scaling in the dimensions is indeed observed empirically, and depicted in the appendix. %~\cite{Stich:2014s}.

On the other three functions, SARP performs consistently better than EP-CMA-1. On $f_{\rm exp}$ with $L=1$\textsc{e}4 in $n=100$ dimensions the factor is roughly 3, its roughly 14 for $L=1$\textsc{e}6  (Fig.~\ref{fig:all}), and exceeds 10 on both $f_{\rm two}$ and $f_{\rm rosen}$ (Fig.~\ref{fig:rosen}).
%
% 3 times more iterations than SARP on $f_{\rm exp}$ with $L=1$\textsc{e}4 in $n=100$ dimensions to reach the target accuracy (14 times more for $L=1$\textsc{e}6). But 
Considering the scaling in dimension (Figs.~\ref{fig:dimfexp}-\ref{fig:rosen}; and the appendix), %
%~\cite{Stich:2014s}), 
we observe that the relative performance (\#ITS$/n$) of SARP remains constant on all four benchmark functions, as predicted by theory for a similar method~\cite{nesterov:11,lee:2013}.
% On $f_{\rm exp}$, the advantage of SARP over EP-CMA-1 for parameter $L=1$\textsc{e}6 is roughly a factor of 10-14 over all dimensions (2-3 for $L=1$\textsc{e}4). Our results on  $f_{\rm lin}$ (not depicted here) show that the advantage of EP-CMA-1 over SARP in $n=20$ dimensions is roughly 20 for parameter $L=1$\textsc{e}4 (120 for parameter $L=1$\textsc{e}6) and then steadily drop to roughly 6 (30 for parameter $L=1$\textsc{e}6) in $n=80$ dimensions.

%\textbf{SARP vs. EP-CMA-1.} 
%In Fig.~\ref{fig:all} we see that EP-CMA-1 requires roughly 3 times more iterations that SARP on $f_{\rm exp}$ with $L=1$\textsc{e}4 in $n=100$ dimensions to reach the target accuracy (14 times more for $L=1$\textsc{e}6). But on $f_{\rm lin}$ it outperforms SARP by a factor of roughly 5 (factor 24 for $L=1$\textsc{e}6).
%Considering the scaling in dimension (Fig.~\ref{fig:dimfexp}), we observe that the relative performance (\#ITS$/n$) of SARP remains constant on all four benchmark functions, as predicted by theory for similar method~\cite{nesterov:11,lee:2013}. On $f_{\rm exp}$, the advantage of SARP over EP-CMA-1 for parameter $L=1$\textsc{e}6 is roughly a factor of 10-14 over all dimensions (2-3 for $L=1$\textsc{e}4). Our results on  $f_{\rm lin}$ (not depicted here) show that the advantage of EP-CMA-1 over SARP in $n=20$ dimensions is roughly 20 for parameter $L=1$\textsc{e}4 (120 for parameter $L=1$\textsc{e}6) and then steadily drop to roughly 6 (30 for parameter $L=1$\textsc{e}6) in $n=80$ dimensions.

\begin{figure}[!t]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_20.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_20.all=1.ls=0-eps-converted-to.pdf}

%\vspace{-1ex}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_40.all=1.ls=0-eps-converted-to.pdf}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_40.all=1.ls=0-eps-converted-to.pdf}
%
%\vspace{-1ex}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_60.all=1.ls=0-eps-converted-to.pdf}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_60.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_80.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_80.all=1.ls=0-eps-converted-to.pdf}

\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ with $L=1$\textsc{e}4 (left) and $L=1$\textsc{e}6 (right) in $n=20$ and $n=80$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
\label{fig:dimfexp}
\end{figure}

\textbf{EP-CMA-schemes.}
The EP-CMA-$m$ schemes consistently work better for increasing values of $m$ throughout the whole benchmark (Figs.~\ref{fig:all}-\ref{fig:rosen}). On $f_{\rm exp}$ with $L=1$\textsc{e}4 the difference in \#ITS between EP-CMA-$n$ and EP-CMA-1 is roughly a factor of 10, and 20 for $L=1$\textsc{e}6  (Fig.~\ref{fig:all}). The gap becomes gradually smaller on $f_{\rm lin}$, $f_{\rm two}$ (especially for $L=1$\textsc{e}6, see appendix), %~\cite{Stich:2014s}), 
and is insignificant on $f_{\rm rosen}$ (Fig.~\ref{fig:rosen}).
%
On $f_{\rm lin}$ the EP-CMA-$m$ schemes perform extremely well, already for small $m$. EP-CMA-4 performs approximately as good as CMA-ES, for both parameters $L=1$\textsc{e}4 and $L=1$\textsc{e}6 (Fig.~\ref{fig:all}).
%
On $f_{\rm exp}$ in $n=100$ dimensions and parameter $L=1$\textsc{e}4, both SARP and EP-CMA-4 need about the same \#ITS. For parameter $L=1$\textsc{e}4  the performance of SARP is the same as the performance of EP-CMA\nobreakdash-$\sqrt{n}$ (Fig.~\ref{fig:all}).
%
On both $f_{\rm two}$ and $f_{\rm rosen}$, the EP-CMA-$m$ scheme cannot reach the performance of SARP, though on $f_{\rm rosen}$ the EP-CMA-$m$ schemes perform as good as CMA-ES (Fig.~\ref{fig:rosen}).
%For $f_{\rm exp}$ in $n=100$ dimensions and parameter $L=1$\textsc{e}4, EP-CMA-$n$
%needs only 1.5 times more \#ITS as CMA-ES to reach the target accuracy (Fig.~\ref{fig:all}). EP-CMA-4 need about the same \#ITS as SARP. The same observations hold for dimensions $n=80$ and $n=60$, however, for $n=40$ SARP performs approximately as good as EP-CMA-2, and for $n=20$ EP-CMA-1 outperforms SARP (Fig.~\ref{fig:dimfexp}). 
%For parameter $L=1$\textsc{e}6 the picture is different. EP-CMA-$n$ performs approximately 10 times worse than CMA-ES and 10 times better than SARP, throughout the dimensions $n=40,\dots, 100$ (Fig.~\ref{fig:all}, \ref{fig:dimfexp}). The performance of SARP is the same as the performance of EP-CMA-$\sqrt{n}$, except for dimension $n=20$, where EP-CMA-4 performs already better than SARP. 
%On $f_{\rm lin}$ the EP-CMA-$m$ schemes perform extremely well, already for small $m$. EP-CMA-4 performs approximately as good as CMA-ES, for both parameters $L=1$\textsc{e}4 and $L=1$\textsc{e}6 (Fig.~\ref{fig:all}), and on all dimensions (not depicted).

%\textbf{EP-CMA-schemes.}
%The EP-CMA-$m$ schemes consistently work better for increasing values of $m$ throughout the whole benchmark (Fig.~\ref{fig:all}-\ref{fig:dimfexp}). For $f_{\rm exp}$ in $n=100$ dimensions and parameter $L=1$\textsc{e}4, EP-CMA-$n$
%needs only 1.5 times more \#ITS as CMA-ES to reach the target accuracy (Fig.~\ref{fig:all}). EP-CMA-4 need about the same \#ITS as SARP. The same observations hold for dimensions $n=80$ and $n=60$, however, for $n=40$ SARP performs approximately as good as EP-CMA-2, and for $n=20$ EP-CMA-1 outperforms SARP (Fig.~\ref{fig:dimfexp}). 
%For parameter $L=1$\textsc{e}6 the picture is different. EP-CMA-$n$ performs approximately 10 times worse than CMA-ES and 10 times better than SARP, throughout the dimensions $n=40,\dots, 100$ (Fig.~\ref{fig:all}, \ref{fig:dimfexp}). The performance of SARP is the same as the performance of EP-CMA-$\sqrt{n}$, except for dimension $n=20$, where EP-CMA-4 performs already better than SARP. On $f_{\rm lin}$ the EP-CMA-$m$ schemes perform extremely well, already for small $m$. EP-CMA-4 performs approximately as good as CMA-ES, for both parameters $L=1$\textsc{e}4 and $L=1$\textsc{e}6 (Fig.~\ref{fig:all}), and on all dimensions (not depicted).

\textbf{CMA-ES.}
Fig.~\ref{fig:dimfexp} shows nicely the quadratic dependence of the performance of CMA-ES on the dimension $n$, see appendix %also~\cite{Stich:2014s} 
where we report the data for all considered dimensions.
The \#ITS of the Random Pursuit schemes (RP, SARP) to reach the target accuracy increases only linearly (the relative performance (\#ITS$/n$) is constant over the dimensions).
In the dimensions $n \leq 100$ considered here, CMA-ES is the best performing scheme on $f_{\rm exp}$ (Fig.~\ref{fig:all}) and $f_{\rm two}$ (Fig.~\ref{fig:rosen}); on $f_{\rm lin}$ the EP-CMA-$m$ schemes match its performance for $m \geq 4$ (Fig.~\ref{fig:all}). A notable exception is the behavior on the non-convex $f_{\rm rosen}$, where only SARP can accelerate and the other schemes, including CMA-ES, require over 10 times more \#ITS to reach the same accuracy (Fig.~\ref{fig:rosen}).
% the performance advantage of CMA-ES is slowly reduced (but sill significantly large in $n=100$ dimensions).
%At last but not least, we like to note that  
%

\begin{figure}[!t]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_20.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_20.all=1.ls=0-eps-converted-to.pdf}

%\vspace{-1ex}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_40.all=1.ls=0-eps-converted-to.pdf}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_40.all=1.ls=0-eps-converted-to.pdf}
%
%\vspace{-1ex}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli4.dim_60.all=1.ls=0-eps-converted-to.pdf}
%\includegraphics[width=.49\textwidth]{matlab/plots/felli6.dim_60.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_80.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_80.all=1.ls=0_withlabel-eps-converted-to.pdf}

\caption{Evolution of FVAL vs. \#ITS on $f_{\rm two}$ with $L=1$\textsc{e}4 (left) and $f_{\rm rosen}$ (right) in $n=20$ and $n=80$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
\label{fig:rosen}
\end{figure}
%
%
\section{Discussion and Conclusions}
\label{sec:disc}
In this contribution we emphasize the importance of accelerated random gradient schemes~\cite{nesterov:11,lee:2013}.
%, i.e. zeroth-order counterparts of the fastest possible first-order methods for convex optimization problems. 
Each iteration in SARP has only linear complexity, yet the scheme takes correlations between successive iterates into account.
%These schemes take correlations between successive iterates into account and have only linear complexity per iteration.
In CMA-ES, such correlations are collected in the evolution path~\cite{Hansen:2001,ostermeier1994} and stored in the covariance matrix. This requires $\Theta(n^2)$ simple operations per iteration. 
The proposed EP-CMA-1 uses as well the information of the evolution path to bias the search, but does not store a full-rank covariance matrix. 
%This yields linear iteration complexity, at the expense of slower convergence. 
%By storing more than one vector at a time, this scheme can easily be generalized to limited memory variants of CMA-ES. We tested the performance of the newly proposed EP-CMA variants and one accelerated random pursuit method (SARP) on two quadratic functions.

\textbf{Line search.}
We empirically tested two Random Pursuit algorithms with an exact line search oracle.
%For the two given benchmark functions the minimizer along the search direction could be computed analytically. 
%This was possible on the given benchmark functions, however, 
Such an oracle is in general not available for general black-box optimization problems and the line search must for instance be implemented as bisection search~(cf.~\cite{Stich:2011,Stich:2012a}) at the expense of additional function evaluations per iteration. The empirical data shows that both Random Pursuit algorithms do perform well if a simple adaptive step size scheme is used instead of the line search. This makes both schemes (especially SARP) promising candidates for black-box optimization, also in high dimension, as the runtime scales only linearly with the dimension. 
%Both schemes scale only linearly with the dimension and can therefore also be used in higher dimensions.
% problems. 
%ue to the linear iteration complexity, these schemes can be used in high dimensions.
%For RP-exact, our findings are in agreement with~\eqref{eq:generalRP}. 
Up to our knowledge, no experimental results for SARP with adaptive step size have been published yet (the authors in~\cite{Stich:2011} considered a line search with high accuracy, almost like SARP-exact). 

\textbf{Acceleration in EP schemes.} 
Our empirical results show that the sole use of the evolution path
%in EP-CMA-1 
can lead to astonishingly good performance---depending on the problem and its eigenvalue spectrum.
%The performance gain was substantial on a function whose smallest eigenvalue is well separated form the second smallest in the spectra. There EP-CMA-1 did not only outperform RP, but also SARP by an order of magnitude. Observe, that knowing one essential direction reduces the condition number of this problem by a factor of roughly $1/n$, thus the findings are in agreement to~\eqref{eq:generalRP}. 
%This effect is much smaller on $f_{\rm exp}$ where most eigenvalues are concentrated around the smallest one. There, the performance gain of EP-CMA-1 over RP was rather small (factor 2-3), and EP-CMA-1 performed an order of magnitude worse than SARP. 
%Our results show that exclusively using the evolution path can lead to better performance, but this increase can be explained by the spectra of the objective function and 
The speed-up of EP-CMA-1 on $f_{\rm lin}$ can be explained by the fact that the condition number of the problem drops once 
the algorithm has learned the most insensitive direction. %(most elongated direction of the level sets) 
%and can compensate for it. 
Hence, the acceleration can be explained by formula~\eqref{eq:generalRP} rather than~\eqref{eq:SARP}. 
For SARP the situation is more promising. 
The data indicates that the convergence on $f_{\rm exp}$ and $f_{\rm two}$ is as described in~\eqref{eq:SARP}. The same seems to be true on the non-convex $f_{\rm rosen}$ where SARP needs an order of magnitude less \#ITS than all other schemes, including CMA-ES. Only on
$f_{\rm lin}$ this does not to hold, as SARP is only one order of magnitude faster than RP.  %The \#ITS of SARP to reach the target accuracy scale only linearly in the dimension (see also~\cite{Stich:2014s}. 
Consider the update~\eqref{eq:sarp}. By expansion we obtain
%. By expanding the update we can write
\begin{align}
\yy_{k+1} %&= \xx_k + \beta \left(\xx_{k+1} - \xx_{k}\right) = \xx_k + \beta \left(\yy_k + \sigma_{k+1} \uu_{k+1} - \xx_k \right) \\
 &= \xx_k + \beta \left(\sigma_{k+1} \uu_{k+1} + \beta \left(\xx_{k} - \xx_{k-1}\right) \right) = \xx_k + \sum_{i=1}^{k+1} \beta^{k+2-i} \sigma_i \uu_i\,.
\end{align}
We see that the drift is a weighted average of the previous steps $\sigma_i \uu_i$. The discount factor $\beta$ is the expected convergence rate. Therefore, the influence of a step $\sigma_i \uu_i$ on $\yy_{k+1}$ is roughly the same for all $i=1,\dots, k$. In contrast to this, the evolution path $\pp_k$ stores only information of the directions of the last steps (but no step sizes). The discount factor is approximately $1-2/n$. Although the evolution path $\pp_k$ is a cumulation of all old steps, the weigh of old steps is exponentially small compared to the influence of the newest steps. We might conclude that the mechanism of  accelerated random schemes like SARP is therefore inherently different to the concept of the evolution path, supporting reports in~\cite{Stich:2012a}. However, we cannot rule out the possibility, that with a different choice of internal parameters of EP-CMA-1 the difference to SARP could be reduced.

\textbf{Limited Memory schemes.}
The performance of the  proposed EP-CMA-$m$ schemes 
uniformly increases for larger parameters $m$, as well as the complexity of each single iteration. 
An optimal trade-off for the parameter $m$ has to be found, depending on the dimension $n$ and the cost of individual function evaluations. 
The data shows that the EP-CMA-$m$ schemes can dramatically improve the performance of simple random search already for small values of $m$. The speed-up depends crucially on the eigenvalue spectra of the objective function. 
It seems that these schemes can not reach the performance of the related variants in~\cite{Loshchilov:2014}.% However, they do not reach the performance of CMA-ES, even for $m=n$.

% For exponential spectra, EP-CMA-$n$
%with $\Theta(n^2)$ iteration complexity
%was still an order of magnitude slower than CMA-ES. The EP-CMA-4 scheme %(with only linear complexity per iteration), 
%was more than two order of magnitudes slower than CMA-ES, but one order of magnitude faster than the plain RP.


We generally conclude, that the here
% proposed limited memory variants of CMA-ES could 
proposed algorithmic schemes with linear iteration complexity could
be a promising way to handle high dimensional black-box optimization problems.
However, the empirical data suggest that there is an intrinsic limitation for the EP schemes, as they depend on the eigenvalue spectrum of the objective function. 
%based only on the evolution path. 
%that use only information from the evolution path. 
%
This behavior is not observed for SARP. We like to advocate that features of accelerated schemes (like SARP) should therefore be taken seriously into account when facing high dimensional problems.
%, as those schemes optimal (linear) scaling of the iteration complexity on the dimension. %Further research should determine if true acceleration can be obtained by using the evolution path in a different way.

%Supplementary online material:~\cite{Stich:2014s}



% ---- Bibliography ----
%
%\bibliographystyle{unsrt}
\bibliographystyle{splncs}
\bibliography{ppsn2014}


\newpage
\appendix
\section{Appendix}
~
\vfill
\begin{figure}[h!]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4.dim_20.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6.dim_20.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4.dim_40.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6.dim_40.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4.dim_60.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6.dim_60.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4.dim_80.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6.dim_80.all=1.ls=0-eps-converted-to.pdf}

\vspace{-1ex}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4.dim_100.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6.dim_100.all=1.ls=0-eps-converted-to.pdf}

\caption{Evolution of FVAL vs. \#ITS on $f_{\rm exp}$ with $L=1$\textsc{e}4 (left) and $L=1$\textsc{e}6 (right) in $n=\{20,40,60,80,100\}$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
\label{fig:supplement:exp}
\end{figure}
\vfill


\begin{figure}[h!]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4lin.dim_20.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6lin.dim_20.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4lin.dim_40.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6lin.dim_40.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4lin.dim_60.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6lin.dim_60.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4lin.dim_80.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6lin.dim_80.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli4lin.dim_100.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/felli6lin.dim_100.all=1.ls=0-eps-converted-to.pdf}

\caption{Evolution of FVAL vs. \#ITS on $f_{\rm lin}$ with $L=1$\textsc{e}4 (left) and $L=1$\textsc{e}6 (right) in $n=\{20,40,60,80,100\}$ dimensions. For 51 runs (11 for RP and $L=1$\textsc{e}6) we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
\label{fig:supplement:lin}
\end{figure} 

\begin{figure}[h!]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_20.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo6.dim_20.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_40.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo6.dim_40.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_60.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo6.dim_60.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_80.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo6.dim_80.all=1.ls=0-eps-converted-to.pdf}

\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo4.dim_100.all=1.ls=0-eps-converted-to.pdf}
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/ftwo6.dim_100.all=1.ls=0-eps-converted-to.pdf}

\caption{Evolution of FVAL vs. \#ITS on $f_{\rm two}$ with $L=1$\textsc{e}4 (left) and $L=1$\textsc{e}6 (right) in $n=\{20,40,60,80,100\}$ dimensions. For 51 runs (see below) we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.
(For $L=1$\textsc{e}6 in dimension $n=100$: only 11 runs of all schemes, except EP-CMA-100 which is omitted.)}
\label{fig:supplement:two}
\end{figure} 

\begin{figure}[h!]
\centering
\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_20.all=1.ls=0-eps-converted-to.pdf}


\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_40.all=1.ls=0-eps-converted-to.pdf}


\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_60.all=1.ls=0-eps-converted-to.pdf}


\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_80.all=1.ls=0-eps-converted-to.pdf}


\includegraphics[width=.49\textwidth]{matlab/plots_supplement/frosen.dim_100.all=1.ls=0-eps-converted-to.pdf}


\caption{Evolution of FVAL vs. \#ITS on $f_{\rm rosen}$ in $n=\{20,40,60,80,100\}$ dimensions. For 51 runs we recorded \#ITS needed to reach FVAL of 1\textsc{e}-9. The trajectory realizing the median values is depicted, mean and one standard deviation are indicated by markers.}
\label{fig:supplement:rosen}
\end{figure}





\end{document}
ppsn14_supplement.bbl
\begin{thebibliography}{10}

\bibitem{Polyak:1987}
Polyak, B.:
\newblock Introduction to Optimization.
\newblock Optimization Software - Inc. (1987)

\bibitem{Nesterov:2004}
Nesterov, Y.:
\newblock {I}ntroductory {L}ectures on {C}onvex {O}ptimization.
\newblock Kluwer (2004)

\bibitem{Broyden:1970}
Broyden, C.G.:
\newblock {The Convergence of a Class of Double-rank Minimization Algorithms 1.
  General Considerations}.
\newblock IMA J. of Appl. Math. \textbf{6}(1) (1970)  76--90

\bibitem{Fletcher:1970}
Fletcher, R.:
\newblock {A new approach to variable metric algorithms}.
\newblock The Computer Journal \textbf{13}(3) (1970)  317--322

\bibitem{Goldfarb:1970}
Goldfarb, D.:
\newblock {A Family of Variable-Metric Methods Derived by Variational Means}.
\newblock Mathematics of Computation \textbf{24}(109) (1970)  23--26

\bibitem{Shanno:1970}
Shanno, D.F.:
\newblock {Conditioning of Quasi-Newton Methods for Function Minimization}.
\newblock Mathematics of Computation \textbf{24}(111) (1970)  647--656

\bibitem{Nocedal:1980}
Nocedal, J.:
\newblock {Updating Quasi-Newton Matrices with Limited Storage}.
\newblock Mathematics of Computation \textbf{35}(151) (1980)  773--782

\bibitem{Liu:1989fu}
Liu, D.C., Nocedal, J.:
\newblock On the limited memory {BFGS} method for large scale optimization.
\newblock Mathematical Programming \textbf{45}(1--3) (1989)  503--528

\bibitem{nesterov:1983}
Nesterov, Y.:
\newblock A method of solving a convex programming problem with convergence
  rate {O}$(1/k^2)$.
\newblock Soviet Mathematics Doklady \textbf{27}(2) (1983)  372--376

\bibitem{nesterov:2007}
Nesterov, Y.:
\newblock Smoothing technique and its applications in semidefinite
  optimization.
\newblock Mathematical Programming \textbf{110}(2) (2007)  245--259

\bibitem{tseng:2008}
Tseng, P.:
\newblock On accelerated proximal gradient methods for convex-concave
  optimization.
\newblock submitted to SIAM Journal on Optimization (2008)

\bibitem{Schumer:1968}
Schumer, M., Steiglitz, K.:
\newblock Adaptive step size random search.
\newblock Automatic Control, IEEE Transactions on \textbf{13}(3) (1968)
  270--276

\bibitem{Rechenberg:1973}
Rechenberg, I.:
\newblock Evolutionsstrategie; Optimierung technischer Systeme nach Prinzipien
  der biologischen Evolution.
\newblock Frommann-Holzboog (1973)

\bibitem{Mutseniyeks:1964}
Mutseniyeks, V.A., Rastrigin, L.A.:
\newblock Extremal control of continuous multi-parameter systems by the method
  of random search.
\newblock Eng.Cyb. \textbf{1} (1964)  82--90

\bibitem{Stich:2011}
Stich, S.U., M\"uller, C.L., G\"artner, B.:
\newblock {Optimization of convex functions with Random Pursuit}.
\newblock SIAM Journal on Optimization \textbf{23}(2) (2013)  1284--1309

\bibitem{nesterov:11}
Nesterov, Y.:
\newblock Random {G}radient-{F}ree {M}inimization of {C}onvex {F}unctions.
\newblock Technical report, ECORE (2011)

\bibitem{Leventhal:2011}
Leventhal, D., Lewis, A.S.:
\newblock {Randomized Hessian estimation and directional search}.
\newblock Optimization \textbf{60}(3) (2011)  329--345

\bibitem{Stich:2012}
Stich, S.U., G\"artner, B., M\"uller, C.L.:
\newblock {Variable Metric Random Pursuit}.
\newblock submitted, http://arxiv.org/abs/1210.5114 (2012)

\bibitem{Hansen:2001}
Hansen, N., Ostermeier, A.:
\newblock {Completely Derandomized Self-Adaption in Evolution Strategies}.
\newblock {Evolutionary Computation} \textbf{9}(2) (2001)  159--195

\bibitem{Hansen:2003}
Hansen, N., Muller, S.D., Koumoutsakos, P.:
\newblock {Reducing the time complexity of the derandomized evolution strategy
  with covariance matrix adaptation (CMA-ES)}.
\newblock Evol.\ Comput. \textbf{11}(1) (2003)  1--18

\bibitem{knight:2007}
Knight, J.N., Lunacek, M.:
\newblock Reducing the {S}pace-time {C}omplexity of the {CMA-ES}.
\newblock In: GECCO '07, ACM (2007)  658--665

\bibitem{Loshchilov:2014}
Loshchilov, I.:
\newblock A {Computationally} {Efficient} {Limited} {Memory} {CMA-ES} for
  {Large} {Scale} {Optimization}.
\newblock to appear, GECCO'14, http://arxiv.org/abs/1404.5520 (2014)

\bibitem{lee:2013}
Lee, Y.T., Sidford, A.:
\newblock Efficient accelerated coordinate descent methods and faster
  algorithms for solving linear systems.
\newblock In: FOCS. IEEE (2013)  147--156

\bibitem{ostermeier1994}
Ostermeier, A., Gawelczyk, A., Hansen, N.:
\newblock Step-size adaptation based on non-local use of selection information.
\newblock In: PPSN III. LNCS.
\newblock (1994)  189--198

\bibitem{sun:2013}
Sun, Y., Schaul, T., Gomez, F., Schmidhuber, J.:
\newblock A linear time natural evolution strategy for non-separable functions.
\newblock In: Proc. 15th Genetic and evolutionary computation conference
  companion, ACM (2013)  61--62

\bibitem{igel:06}
Igel, C., Suttorp, T., Hansen, N.:
\newblock {A Computational Efficient Covariance Matrix Update and a (1+1)-CMA
  for Evolution Strategies}.
\newblock In: {GECCO}. (2006)  453--460

\bibitem{Stich:2012a}
Stich, S.U., M\"uller, C.L.:
\newblock {On Spectral Invariance of Randomized Hessian and Covariance Matrix
  Adaptation Schemes}.
\newblock In: PPSN XII.
\newblock (2012)  448--457

\end{thebibliography}
