\section{Introduction}
There is a growing interest in multimodal techniques for solving  Music Information Retrieval (MIR) problems.
Music performances have a highly multimodal content and the different modalities involved are highly correlated: sounds are emitted by the motion of the player performing and in chamber music performances the scores constitute an additional encoding that may be as well leveraged for the automatic analysis of music \cite{li2017see}.

A fundamental problem in music analysis and in general in audio processing is Blind Source Separation (BSS). BSS consists in, given a mixture of signals, recovering the individual signals the mixture is conformed by. Mathematically, a mixture of sounds can be expressed as the sum of individual sources: $s_m = \sum s_i$. Thus, the BSS problem consists in recovering each $s_i$ for a given $s_m$. In speech, it is also known as the Cocktail Party problem, which refers to the task
of recognizing an individual speech in noisy social environments \cite{cherry1953some}.
Single-channel source separation problem can be approached from an audio-only perspective using techniques
such as independent component analysis (ICA) \cite{hyvarinen2000independent}, sparse decomposition \cite{zibulevsky2001blind}, nonnegative
matrix factorization (NMF) \cite{virtanen2007monaural}, computational auditory scene analysis (CASA) \cite{ellis1996prediction}, probabilistic
latent component analysis (PLCA) \cite{smaragdis2006probabilistic} or deep learning techniques (e.g. \cite{chandna2017monoaural, stoller2018wave }).

On the other side, by
visually inspecting the scene we may extract information about the number of sound sources, their
type, spatio-temporal location and also motion, which naturally relates to the emitted sound. Besides, it is possible to carry out self-supervised tasks in which one modality supervises the other one. This entails another research field, the cross-modal correspondence (CMC). We can find pioneering works for both problems BSS and CMC.  \cite{hershey2000audio, kidron2005pixels} make use of audio-visual data for sound localization and \cite{darrell2000audio}, \cite{sodoyer2002separation},\cite{rivet}  for speech separation. In the context of music, visual information has also proven to help model-based methods both in source separation \cite{li2017audiovisual, parekh2017guiding} and localization \cite{li2017see}.
With the flourishing of deep learning techniques many recent works exploit both, audio and video content, to perform music source separation \cite{gao2019co, SoM, xu2019recursive}, source association \cite{li2019online}, localization \cite{arandjelovic2018objects} or both \cite{SoP}. Some CMC works explore features generated from synchronization \cite{owenssync,kobar} and prove these features are reusable for source separation. These works use networks that have been trained in a self-supervised way using pairs of corresponding/non-corresponding audio-visual signals for localization purposes \cite{arandjelovic2018objects} or the mix-and-separate approach for source separation \cite{SoP, gao2019co, SoM, xu2019recursive}. Despite deep learning made possible to solve classical problems in a different way, it also contributed to create new research fields like cross-modal generation, in which the main aim is to generate video from audio\cite{speech2face,rochestergen} or viceversa \cite{visual2sound}. More recent works related to human motion make use of skeleton as an inner representation of the body which can be further converted into video \cite{audio2body,conversationalgestures} which shows the potential of skeletons.
The main contribution of this paper is Solos, a new dataset of musical performance recordings of soloists that can be used to train deep neural networks for any of the aforementioned fields. Compared to a similar dataset of musical instruments presented in \cite{SoP} and its extended version \cite{Musices}, our dataset does contain the same type of chamber orchestra instruments present in the URMP dataset. Solos is a dataset of 755 real-world recordings gathered from YouTube which provides several features missing in the aforementioned datasets:  skeletons and high quality timestamps. Source localization is usually indirectly learned by networks. Thus, providing a practical localization ground-truth is not straightforward. Nevertheless, networks often point to the player hands as if they were the sound source. We expect hands localization can help to provide additional cues to improve audio-visual BSS or can be used as source ground-truth localization. In order to show the benefits of using Solos we trained some popular BSS architectures and compare their results.

\section{Related Work}

The University of Rochester Multi-Modal Music Performance Dataset (URMP) \cite{URPM} is a dataset with 44 multi-instrument video recordings of classical music pieces. Each instrument present in a piece was recorded separately, both with video and high-quality audio with a stand-alone microphone, in order to have ground-truth individual tracks. Although playing separately, the instruments were coordinated by using a conducting video with a pianist playing in order to set the common timing for the different players. After synchronization, the audio of the individual videos was replaced by the high-quality audio of the microphone and then different recordings were assembled to create the mixture: the individual high-quality audio recordings were added up to create the audio mixture and the visual content was composited in a single video with a common background where all players were arranged at the same level from left to right. For each piece, the dataset provides the musical score in MIDI format, the high-quality individual instrument audio recordings and the videos of the assembled pieces. The instruments present in the dataset, shown in Figure \ref{fig: categories}, are common instruments in chamber orchestras. In spite of all its good characteristics, it is a small dataset and thus not appropriate for training  deep learning architectures.

Two other datasets of audio-visual recordings of musical instruments performances have been presented recently:
Music \cite{SoP} and MusicES \cite{Musices}. Music consists of 536 recordings of solos and 149 videos of duets across 11 categories: accordion, acoustic guitar, cello, clarinet, erhu, flute, saxophone, trumpet, tuba, violin and xylophone. This dataset was gathered by querying YouTube. MusicES \cite{Musices} is an extension of MUSIC to around the triple of its original size with approximately  1475 recordings but spread in 9 categories instead: accordion, guitar, cello, flute, saxophone, trumpet, tuba, violin and xylophone.  There are 7  common categories in MUSIC and Solos: violin, cello, flute, clarinet,  saxophone, trumpet and tuba. The common categories between MusicES and Solos are 6 (the former ones except clarinet). Solos and MusicES are complementary.  There is only an small intersection of 5\% between both, which means both datasets can be combined into a bigger one.

We can find in the literature several examples which show the utility of audio-visual datasets. \textit{The Sound of Pixels} \cite{SoP} performs audio source separation generating audio spectral components which are further smartly selected by using visual features coming from the video stream to obtain separated sources. This idea was further extended in \cite{xu2019recursive} in order to separate the different sounds present in the mixture in a recursive way. At each stage, the system separates the most salient source from the ones remaining in the mixture.  \textit{The Sound of Motions} \cite{SoM} uses dense trajectories obtained from optical flow to condition audio source separation, being able even to separate same-instrument mixtures. Visual conditioning is also used in \cite{gao2019co} to separate different instruments; during training, a classification loss is used on the separated sounds to enforce object consistency and a co-separation loss forces the estimated individual sounds to produce the original mixtures once reassembled.
In
\cite{parekh2017guiding}, the authors developed an energy-based method which minimizes a \textit{Non-Negative Matrix Factorization} term with an activation matrix which is forced to be aligned to a matrix containing per-source motion information. This motion matrix
contains the average magnitude velocities of the clustered motion trajectories in each player bounding box.

Recent works show the rising use of skeletons in audiovisual tasks. In  \textit{Audio to body dynamics} \cite{audio2body} authors show it is possible to predict skeletons reproducing the movements of players playing instruments such as piano or violin.
Skeletons have proven to be useful for establishing audio-visual correspondences, such as body or finger motion with note onsets or pitch fluctuations, in chamber music performances  \cite{li2019online}. A recent work
\cite{music2gesture} tackles the source separation problem in a
 similar to \textit{Sound of Motions} \cite{SoM} but  replacing the dense trajectories by skeleton information.

\section{Dataset}
Solos\footnote{\label{github}Dataset available at https://juanfmontesinos.github.io/Solos/}
was designed to have the same categories as the URMP \cite{URPM} dataset, so that URMP can be used as testing dataset in a real-world scenario. This way we aim to establish a standard way of evaluating source separation algorithms' performance avoiding the use of mix-and-separate  in testing.
Solos consists of 755 recordings distributed amongst 13 categories as shown in Figure \ref{fig: categories}, with an average amount of 58 recordings per category and an average duration of 5:16 min. It is interesting to highlight that, for 8 out of 13 categories, the median of resolution is HD, despite being a YouTube-gathered dataset. Per-category statistics can be found in Table \ref{tab:solos_stats}. These recordings were gathered by querying YouTube using the tags solo and auditions in several languages such as English, Spanish, French, Italian, Chinese or Russian.

\begin{figure}[t]
\includegraphics[width=0.5\textwidth]{images/players_instruments.png}
\caption{Solos and URMP instrument categories. Image adapted from \cite{URPM}.}
\label{fig: categories}
\end{figure}

\begin{table}[]
\centering
\begin{tabular}{lccc}
\multicolumn{1}{l}{Category} & \# Recordings & Mean duration & Median resolution \\ \hline
Violin                       & 66            & 6:16          & 1080x720          \\
Viola                        & 55            & 5:31          & 1280x720          \\
Cello                        & 134           & 7:21          & 640x480           \\
DoubleBass                   & 58            & 8:53          & 1280x720          \\ \hdashline

Flute                        & 48            & 4:00          & 640x360           \\
Oboe                         & 53            & 5:45          & 1280x720          \\
Clarinet                     & 49            & 3:23          & 640x360           \\
Bassoon                      & 56            & 5:08          & 1280x720          \\
Saxophone                    & 45            & 2:42          & 1280x720          \\ \hdashline

Trumpet                      & 50            & 1:14          & 640x360           \\
Horn                         & 50            & 5:11          & 1280x720          \\
Trombone                     & 50            & 5:03          & 1280x720          \\
Tuba                         & 41            & 2:49          & 640x360           \\ \hdashline

{\bf TOTAL}                  & {\bf 755}     & {\bf 5:16}    & {\bf 854x480}
\end{tabular}
\caption{Statistics of Solos Dataset}
\label{tab:solos_stats}
\end{table}
\begin{comment}
Solos is not only a set of recordings. We also provide OpenPose body and hand skeletons for each frame of each recording and timestamps indicating useful parts. To do so, video streams are re-sampled to 25 FPS keeping the audio stream intact. An iterative process returns stamps for which there are at least N frames with a detected hand and no more than M consecutive mispredictions. In practice we use N=150 and M=5, thus, a minimum of 6 seconds of video with at most 5 consecutive frames with hand mispredictions. At this point, we have segments of video in which there are hands detected. To refine these results we further applied an energy-based silence detector which allows to discard those segments in which the instrument is not being played, e.g., transitions, music sheet changes, etcetera. Besides, we perform a linear interpolation of the mispredicted keypoints in a relative base of coordinates.
Directly interpolating the absolute coordinates would lead to deformations of the skeleton and inaccuracies.  Since skeletons are tree-like graphs  it is possible to interpolate the relative coordinates of each joint (node in the graph) with respect to its parent node. Then, the absolute coordinates of the joint are recovered with the sum of the absolute coordinates of its parent and the estimated relative coordinates with respect to the parent. Let us denote by $J_i
^t$ the relative coordinates of the $i$th joint with respect to its parent at time $t$. On the other hand, $\hat{J}^t_i$ denotes the estimated value of $J_i^t$ when the $i$th joint is mispredicted. $\hat{J}^t_i$ can be linearly interpolated using the relative coordinates of the closest $i$th detected joint before time $t$ (i.e $J^{t^-}_i$ where $t^- < t$), and analogously with the closest $i$th detected joint after time $t$ (i.e $J^{t^+}_i$ where $t < t
^+$). For example, given the following sequence of detected and misdetected coordinates (that need to be estimated), $J$ and $\hat{J}$ respectively:
\begin{equation*}
\label{eq:set}
\{J_i^{t-n},\hat{J}_i^{t-n+1},...,\hat{J}_i^{t},...,\hat{J}_i^{t+m-1},J_i^{t+m}\}
\end{equation*}
then, the interpolation at time $t$ can be calculated as:
\begin{equation}
\label{eq:interp}
\hat{J}^t_i=\frac{m \, J_i^{t-n}  + n \, J_i^{t+m}}{m+n}.
\end{equation}
OpenPose maps mispredicted joints to the origin of coordinates. We empirically found that such a big jump in the position of a joint induces noise. Using interpolated coordinates helps to address this problem.
\end{comment}
\subsection{OpenPose Skeletons}
Solos is not only a set of recordings. Apart from the videos identificators We also provide: i) body and hand  skeletons estimated by OpenPose \cite{openpose} in each frame of each recording and ii) timestamps indicating useful parts. OpenPose is a system capable to predict body skeleton and hands skeletons making use of two different neural networks. To do so, they predict a confidence map of the belief that a specific body part may be located at any given pixel as well as part affinity fields which encode the degree of association between different body parts. Finally, it predicts 2D skeletons and per-joint confidence via greedy inference. In practice, the body skeleton is estimated with a first network. Then, the position  of the wrists in the body skeleton are used to estimate the position of both hands. A second neural network obtains the skeleton of each hand independently. Note that since each body part is estimated independently, OpenPose makes no assumptions about the limbs to find. It just calculates the most likely skeleton given confidence maps and part affinity fields. The whole process is carried out frame-wise. This leads to a small flickering and mispredictions  between frames.

\subsection{Timestamps estimation and skeleton refinement}
Video streams are re-sampled to 25 FPS keeping the audio stream intact. An iterative process returns stamps for which there are at least N frames with a detected hand and no more than M consecutive mispredictions. In practice we use N=150 and M=5, thus, a minimum of 6 seconds of video with at most 5 consecutive frames with hand mispredictions. At this point, we have segments of video in which there are hands detected. To refine these results we further applied an energy-based silence detector which allows to discard those segments in which the instrument is not being played, e.g., transitions, music sheet changes, etcetera. Besides, we perform a linear interpolation of the mispredicted keypoints in a relative base of coordinates.
Directly interpolating the absolute coordinates would lead to deformations of the skeleton and inaccuracies.  Since skeletons are tree-like graphs  it is possible to interpolate the relative coordinates of each joint (node in the graph) with respect to its parent node. Then, the absolute coordinates of the joint are recovered with the sum of the absolute coordinates of its parent and the estimated relative coordinates with respect to the parent. Let us denote by $J_i
^t$ the relative coordinates of the $i$th joint with respect to its parent at time $t$. On the other hand, $\hat{J}^t_i$ denotes the estimated value of $J_i^t$ when the $i$th joint is mispredicted. $\hat{J}^t_i$ can be linearly interpolated using the relative coordinates of the closest $i$th detected joint before time $t$ (i.e $J^{t^-}_i$ where $t^- < t$), and analogously with the closest $i$th detected joint after time $t$ (i.e $J^{t^+}_i$ where $t < t
^+$). For example, given the following sequence of detected and misdetected coordinates (that need to be estimated), $J$ and $\hat{J}$ respectively:
\begin{equation*}
\label{eq:set}
\{J_i^{t-n},\hat{J}_i^{t-n+1},...,\hat{J}_i^{t},...,\hat{J}_i^{t+m-1},J_i^{t+m}\}
\end{equation*}
then, the interpolation at time $t$ can be calculated as:
\begin{equation}
\label{eq:interp}
\hat{J}^t_i=\frac{m \, J_i^{t-n}  + n \, J_i^{t+m}}{m+n}.
\end{equation}
OpenPose maps mispredicted joints to the origin of coordinates. We empirically found that such a big jump in the position of a joint induces noise. Using interpolated coordinates helps to address this problem.
\vspace{6mm}

\section{Experiments}
In order to show the suitability of Solos, we have focused in the blind source separation problem and have trained  \textit{The Sound of Pixels} (SoP) \cite{SoP} and the Multi-head U-Net (MHU-Net) \cite{multihead} models on the new dataset. We have carried out four experiments: i) we have evaluated the SoP pre-trained model  provided by the authors; ii) we have trained SoP from  scratch; iii) we have fine-tuned SoP on Solos starting from the weights of the pre-trained model on MUSIC and iv) we have trained the Multi-head U-Net from scratch.
MHU-Net has been trained to separate mixtures with the number of sources varied from two to seven following a curriculum learning procedure as it improves the results. SoP has been trained according to the optimal strategy described in \cite{SoP}.
\vspace{4mm}
\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{images/models.PNG}
\caption{Considered architectures. Left, Sound of Pixels: The network takes as input a mixture spectrogram and returns a binary mask given the visual feature vector of the desired source. Right, Multi-Head U-Net: It takes as input a mixture spectrogram and returns 13 ratio masks, one per decoder. }
\label{fig: models}
\end{figure*}

Evaluation is performed on the URMP dataset \cite{URPM} using the real mixtures they provide. URMP tracks are sequentially split in 6s-duration segments. Metrics are obtained from all the resulting splits.

\subsection{Architectures and training details}
We have chosen \textit{The Sound of Pixels} as baseline since its weights are publicly available and the network is trained in a straight-forward way. SoP is composed of three main sub-networks: A dilated ResNet\cite{DResnet} as  video-analysis network, a U-Net\cite{unet} as  audio-processing network and an audio synthesizer network.  We also compare its results against a Multi-head U-Net \cite{multihead}.
\vspace{4mm}

U-Net  \cite{ronneberger2015u} is an encoder-decoder architecture with skip connections in between. Skip connections help to recover the original spatial structure.
MHU-Net is a step forward as it consist of as many decoders as possible sources. Each decoder is specialized in a single source, thus improving performance.

\textit{The Sound of Pixels}\cite{SoP} does not follow the original U-Net architecture proposed for biomedical imaging, but the U-Net described at \cite{unet}, which was tuned for singing voice separation. Instead of having two convolutions per block followed by max-pooling, they use a single convolution with a bigger kernel and striding. The original work  proposes a central block with learnable parameters whereas the central block is a static latent space in \textit{SoP}. U-Net has been widely used as backbone of several architectures for tasks such us image generation \cite{unet-img-synt}, noise suppression and super-resolution \cite{mao2016image}, image-to-image translation\cite{pix2pix}, image segmentation \cite{ronneberger2015u} or audio source separation \cite{unet}.
SoP U-Net consists of 7 blocks with 32, 64, 128, 256, 512, 512 and 512 channels respectively (6 blocks for the MHU-Net). The latent space can be considered as the last output of the encoder.
Dilated ResNet is a ResNet-like architecture which makes use of dilated convolutions to keep the receptive field while increasing the resulting spatial resolution.
The output of the U-Net is a set of 32 spectral components (channels) which are the same size than the input spectrogram, in case of SoP, and a single source per decoder in case of MHU-Net. Given a representative frame, visual features are obtained using the Dilated ResNet. These visual features are nothing but a vector of 32 elements (which corresponds to the number of output channels of U-Net) which are used to select proper spectral components. This selection is performed by the audio analysis network which consist of 32 learnable parameters, $\alpha_k$,  plus a bias, $\beta$. This operation can be mathematically described as follows:

\begin{equation*}
    \beta+\sum_{k=1}^{32} \alpha_k v_{f_k}S_k(t,f),
\end{equation*}
where $S_k(t,f)$ is the $k$-th predicted spectral component at time-frequency bin $(t,f)$.

Figure \ref{fig: models} illustrates the SoP configuration. It is interesting to highlight that making the visual network to  select the spectral components forces it to indirectly learn instrument localization, which can be inferred via activation maps.
\vspace{4mm}

On one hand, MHU-Net has been trained using a curriculum learning strategy that consists of a gradual increment on the amount of sources present in the mixture from two to four. When the loss stays on a plateau for more than 160,000 iterations, the amount of sources is increased by one. We have used mean-square error loss, ADAM  optimizer \cite{Adam}, an initial learning rate (LR)  of $10^{-4}$, weight decay of $10^{-5}$  and dropout of $0.2$ in the decoder. We have also reduced the LR by a half if the loss stays on a plateau for more than 400,000 iterations.

On the other hand, SoP has been trained using a LR of $10^{-3}$ for the U-Net and a LR of $10^{-4}$ for the Dilated ResNet as it was pre-trained on ImageNet.
We have applied a weight on the gradients  based on the magnitude of the mixture spectrogram so that time-frequency points of the predicted source/s contribute to the loss according to the energy of the analogous time-frequency points in the mixture spectrogram.  This reduces overfitting since, given a source,
a time-frequency bin with a low value may be assigned either to one or zero in the ground-truth mask depending on the recorded noise and such weights help reduce its impact on the training. We used different training strategies for SoP and MHU-Net as the optimal training for SoP harms the performance of the MHU-Net.

Ground-truth mask calculation for both SoP and MHU-Net are  described in Eq. \eqref{eq:mask} and Eq. \eqref{eq:mask_ratio}, Sec. \ref{sec: GT}.

\subsection{Data pre-processing}
In order to train the aforementioned architectures, audio is re-sampled to 11025 Hz and 16 bit. Samples fed into the network are 6s duration. We use Short-time Fourier Transform (STFT) to obtain time-frequency representations of waveforms. Following \cite{SoP}, STFT is computed using Hanning window of length 1022 and hop length 256 so that we obtain a spectrogram of size 512$\times$256 for a 6s sample. Later on, we apply a log re-scale on the frequency axis expanding lower frequencies and compressing higher ones. Lastly, we convert  magnitude spectrograms into dB w.r.t. the minimum value of each spectrogram and normalize between -1 and 1.
\vspace{4mm}
\subsection{Ground-truth mask} \label{sec: GT}
Before introducing ground-truth mask computations we would like to point out some considerations. Standard  floating-point audio format imposes a waveform to be bounded between -1 and 1.  At the time of creating artificial mixtures resulting waveforms may be out of these bounds. This can help neural networks to find shortcuts to overfit. To avoid this behaviour spectrograms are clamped according to the equivalent bounds in the time-frequency domain.
\vspace{4mm}

The Discrete Short-time Fourier Transform  can be computed as described in \cite{STFT}:
\begin{equation*}
    S[t,f] = \sum_{k=0}^{L-1}s[k]w[k-t]e^{-j\pi fk/L}.
\end{equation*}
Since $s[k] \in [-1,1]$ it can be easily shown that:
$$
    |S[t,f]|  \leq \sum_{k=0}^{L-1}|w[k-t]|,
$$
i.e., that the  magnitude STFT of an audio signal bounded between [-1,1] is bounded between $[0,\sum |w[k]|]$. Thus,
given  the STFT of N waveforms,  the spectogram of a mixture of sounds is defined the following way:
\begin{equation*}
\label{eq: mix}
     S_{mix}(t,f) = \min\left(\sum_{n=1}^N S_n(t,f),\sum |w[k]|\right),
\end{equation*}
which is equivalent to:
\begin{equation*}
\label{eq: mix_time}
     S_{mix}(t,f) = STFT \left\{\min \left(1,\max\left(\sum_{n=1}^N s_n(t),-1\right)\right)\right\}.
\end{equation*}
For training Sound of Pixels  we have used complementary binary masks as ground-truth masks, defined as:

\begin{equation}
\label{eq:mask}
   M_n(t,f) =
    \begin{cases}
      1, & \text{if}\  |S_n(t,f)|\geq |S_{m}(t,f)| \, \, \, \forall m \in \{1, ..., N\},\\
      0, & \text{otherwise}.
    \end{cases}
\end{equation}
The Multi-head U-Net has been trained with complementary ratio masks, defined as:
\begin{equation}
\label{eq:mask_ratio}
   M_n(t,f) = \frac{|S_n(t,f)|}{ |S_{mix}(t,f)|}.
\end{equation}

\subsection{Results}
Benchmark results for Source to Distortion Ratio (SDR), Source to Interferences Ratio (SIR), Sources to Artifacts Ratio (SAR) proposed in \cite{sdr} are shown in Table \ref{tab:benchmark} in terms of mean and standard deviation.
As it can be observed, Sound of Pixels evaluated using its original weights performs the worst. One possible reason for that could be the absence of some of the URMP categories on the  MUSIC dataset. If we train the network from scratch on Solos, results improve by almost 1 dB.  However, it is possible to achieve an even better result fine-tuning the network, pre-trained with MUSIC, on Solos. We hypothesize that the improvement occurs as the network is exposed to much more training data. Moreover, the table results show how it is possible to reach higher performance by using more powerful architectures like MHU-Net.
\begin{table}[]
    \centering
    \begin{tabular}{c|c|c|c}
          & SDR $\uparrow$ & SIR $\uparrow$ & SAR $\uparrow$ \\ \hline
         SoP \cite{SoP} & $-3.76\pm4.00$ & $-1.45\pm4.68$ & $7.56\pm3.13$ \\
         SoP-Solos & $-2.98\pm5.07$ & $0.46\pm6.76$ & $6.37\pm2.94$ \\
         SoP-ft & $-2.57\pm4.99$ & $0.47\pm6.43$ & $6.89\pm2.48$ \\
         MHU-Net & $ -0.56\pm5.96 $ & $ 1.04\pm7.24 $ & $ 10.37\pm3.48 $ \\
    \end{tabular}
    \caption{Benchmark results (mean $\pm$ standard deviation). SoP:Sound of Pixels original weights, SoP-Solos: Sound of Pixels trained from scratch on Solos. SoP-ft: Sound of Pixels finetuned on Solos. MHU-Net: Multi-head U-Net with 13 decoders. }
    \label{tab:benchmark}
\end{table}

\section{Conclusions}
We have presented Solos, a new audio-visual dataset of music recordings of soloists, suitable for different self-supervised learning tasks such as source separation  using the mix-and-separate strategy, sound localization, cross-modal generation and finding audio-visual correspondences. There are 13 different instruments in the dataset; those are common instruments in chamber orchestras and the ones included in the University  of  Rochester  Multi-Modal  Music  Performance (URMP) dataset \cite{URPM}. The characteristics of URMP -- small dataset of real performances with ground truth individual stems -- make it a suitable dataset for testing purposes but to the best of our knowledge, to date there is no existing large-scale dataset with the same instruments as in URMP. Two different networks for audio-visual source separation based on the U-Net architecture have been trained in the new dataset and further evaluated in URMP, showing the impact of training on the same set of instruments as the test set. Moreover, Solos provides skeletons and timestamps to video intervals where hands are sufficiently visible. This information could be useful for training purposes and also for learning to solve the task of sound localization.
\begin{comment}
We have presented Solos, a new audio-visual dataset of music recordings of soloists, suitable for  training source separation deep neural networks using the mix-and-separate strategy for self-supervised learning. There are 13 different instruments in the dataset; those are common instruments in chamber orchestras and the ones included in the University  of  Rochester  Multi-Modal  Music  Performance (URMP) dataset \cite{URPM}. The characteristics of URMP -- small dataset of real performances with ground truth individual stems -- make it a suitable dataset for testing purposes but to the best of our knowledge, to date there is no existing large-scale dataset with the same instruments as in URMP. Two different networks \gh{for audio-visual source seeparation} based on the U-Net architecture have been trained in the new dataset and further evaluated in URMP, showing the impact of training on the same set of instruments as the test set. Moreover, Solos provides skeletons and timestamps to video intervals where hands are sufficiently visible. This information could be useful for training purposes and also for learning to solve the task of sound localization.
\end{comment}

\bibliographystyle{IEEEtran}


\section{Introduction}
We study the problem of zero-order stochastic optimization, in which we aim to minimize
an unknown strongly convex function via a sequential exploration of its function values,
under measurement error, and a closely related problem of continuous (or continuum-armed) stochastic bandits.  These problems
have received significant attention in the literature, see \cite{agarwal2010,agarwal2011,BP2016,bartlett2018simple,belloni2015escaping,bubeck2012regret,bubeck2017kernel,dvurechensky2018accelerated,hu2016bandit,liang2014zeroth,wang2017stochastic,flaxman2004,Locatelli,malherbe2017global,Shamir13,NS17,Shamir17,duchi2015,rakhlin2011making,saha2011,jamieson2012,Shalev-Shwartz2011}, and are fundamental for many applications in which the derivatives of the function are either too expensive or impossible to compute.
A principal goal of this paper is to exploit higher order smoothness properties of the underlying function in order to improve the performance of search algorithms. We derive upper bounds on the estimation error for a class of projected gradient-like algorithms, as well as close matching lower bounds, that characterize the role played by the number of iterations,
the strong convexity parameter, the smoothness parameter, the number of variables, and the noise level.

Let $f: \mathbb{R}^d \rightarrow \mathbb{R}$ be the function that we wish to minimize over a closed convex subset $\Theta$ of $\mathbb{R}^d$.
Our approach, outlined in Algorithm \ref{algo}, builds upon previous work in which a sequential algorithm queries at each iteration a pair of function values, under a general noise model.
Specifically, at iteration $t$ the current guess $x_t$ for the minimizer of $f$ is used to build two perturbations $x_t+\delta_t$ and $x_t-\delta_t$, where the function values are queried subject to additive measurement errors $\xi_t$ and $\xi_t'$, respectively. The values $\delta_t$ can be chosen in different ways. In this paper, we set $ \delta_t=h_tr_r\zeta_t$ (Line 1), where $h_t>0$ is a suitably chosen small parameter, $r_t$ is random and uniformly distributed on $[-1,1]$, and $\zeta_t$ is uniformly distributed on the unit sphere.  The estimate for the gradient is then computed at Line 2 and used inside a projected gradient method scheme to compute the next exploration point.
We introduce a suitably chosen kernel $K$ that allows us to take advantage of higher order smoothness of~$f$.

The idea of using randomized procedures for derivative-free stochastic optimization can be traced back to Nemirovski and Yudin \cite[Sec.~9.3]{NY1983} who suggested an algorithm with one query per step at point $x_t+h_t\zeta_t$, with $\zeta_t$ uniform on the unit sphere. Its versions with one, two or more queries were studied in several papers including~\cite{agarwal2010,BP2016,flaxman2004,Shamir17}.
Using two queries per step leads to better performance bounds as emphasized in~\cite{PT90, agarwal2010,BP2016,flaxman2004,Shamir17,duchi2015}.
Randomizing sequences other than uniform on the sphere were also explored: $\zeta_t$ uniformly distributed on a cube \cite{PT90}, Gaussian $\zeta_t$~\cite{Nesterov2011, NS17}, $\zeta_t$ uniformly distributed on the vertices of a cube~\cite{Shamir13} or satisfying some general assumptions~\cite{dippon,duchi2015}. Except for \cite{PT90,dippon,BP2016}, these works study settings with low smoothness of $f$ (2-smooth or less) and do not invoke kernels $K$ (i.e. $K(\cdot)\equiv 1$ and $r_t\equiv 1$ in Algorithm \ref{algo}). The use of randomization with smoothing kernels was proposed by Polyak and Tsybakov \cite{PT90} and further developed by Dippon \cite{dippon}, and
Bach and Perchet~\cite{BP2016} to whom the current form of Algorithm~\ref{algo} is due.
\begin{algorithm}[t!]
\caption{Zero-Order Stochastic Projected Gradient} \label{algo}
\begin{algorithmic}
\State
\State {\bfseries Requires} ~ Kernel $K :[-1,1]\rightarrow \mathbb{R}$, step size $\eta_t>0$ and parameter $h_t$, for $t=1,\dots,T$
\vspace{.1cm}
\State {\bfseries Initialization} ~ Generate scalars $r_1,\dots,r_T$ uniformly on the interval $[-1,1]$, vectors $\zeta_1,\dots,\zeta_T$ uniformly distributed on the unit sphere $S_d=\{\zeta\in \mathbb{R}^d:  \norm{\zeta}=1\}$, and choose $x_1\in \Theta$
\vspace{.1cm}
\State {\bfseries For}  $t=1,\dots, T$
\vspace{.1cm}
\State \qquad {1.} ~~~Let $y_t = f(x_t+h_tr_t\zeta_t)+ \xi_t$ and $y_t' = f(x_t-h_tr_t\zeta_t)+ \xi_t',$
\vspace{.05cm}
\State \qquad 2. ~~~Define $\hat g_t = \frac{d}{2h_t} (y_t- y'_t) \zeta_t K(r_t)$
\vspace{.05cm}
\State \qquad 3. ~~~Update $x_{t+1}= \ensuremath{\text{\rm Proj}}_{\Theta}(x_{t}- \eta_t \hat g_t)$
\vspace{.05cm}
\State {\bfseries Return} ~ $(x_t)_{t=1}^T$
\end{algorithmic}
\end{algorithm}

In this paper we consider higher order smooth functions $f$ satisfying the generalized H\"older condition with parameter $\beta \geq 2$, cf. inequality~\eqref{eq:Hclass} below. For integer $\beta$, this parameter can be roughly interpreted as the number of bounded derivatives. Furthermore, we assume that $f$ is $\alpha$-strongly convex. For such functions, we address the following two main questions:
\begin{itemize}
\item[(a)] What is the performance of Algorithm \ref{algo}  in terms of the cumulative regret and optimization error, namely what is the explicit dependency of the rate on the main parameters $d,T,\alpha,\beta$?
\item[(b)] What are the fundamental limits of any sequential search procedure expressed in terms of minimax
optimization error?
\end{itemize}
To handle task (a), we prove upper bounds for Algorithm \ref{algo}, and to handle (b), we prove minimax lower bounds for any sequential search method.

{\bf Contributions.} Our main contributions can be summarized as follows: {\bf i)} Under an adversarial noise assumption (cf. Assumption \ref{ass1} below), we establish for all $\beta\ge 2$ upper bounds of the order $\frac{d^2}{\alpha}T^{-\frac{\beta-1}{\beta}}$ for the optimization risk and $\frac{d^2}{\alpha}T^{\frac{1}{\beta}}$ for the cumulative regret of Algorithm \ref{algo}, both for its constrained and unconstrained versions; {\bf ii)} In the case of independent noise satisfying some natural assumptions (including the Gaussian noise), we prove a minimax lower bound of the order $\frac{d}{\alpha}T^{-\frac{\beta-1}{\beta}}$ for the optimization risk when $\alpha$ is not very small. This shows that to within the factor of $d$ the bound for Algorithm \ref{algo} cannot be improved for all $\beta\ge 2$; {\bf iii)} We show that, when $\alpha$ is too small, below some specified threshold, higher order smoothness does not help to improve the convergence rate. We prove that in this regime the rate cannot be faster than $d/\sqrt{T}$, which is not better (to within the dependency on $d$) than for derivative-free minimization of simply convex functions \cite{agarwal2011,liang2014zeroth}; {\bf iv)} For $\beta=2$, we obtain a bracketing of the optimal rate between $O(d/\sqrt{\alpha T})$ and $\Omega(d/(\max(1,\alpha)\sqrt{T}))$. In a special case when $\alpha$ is a fixed numerical constant, this validates a conjecture in \cite{Shamir13} (claimed there as proved fact) that the optimal rate for $\beta=2$ scales as $d/\sqrt{T}$; {\bf v)} We propose a simple algorithm of estimation of the value $\min_x f(x)$ requiring three queries per step and attaining the optimal rate $1/\sqrt{T}$ for all $\beta\ge 2$. The best previous work on this problem  \cite{belitser} suggested a method with exponential complexity and proved a bound of the order $c(d,\alpha)/\sqrt{T}$ for $\beta> 2$ where $c(d,\alpha)$ is an unspecified constant.

{\bf Notation.} Throughout the paper we use the following notation. We let $\langle \cdot, \cdot \rangle$ and $\|\cdot\|$ be the standard inner product and Euclidean norm on $\mathbb{R}^d$, respectively. For every close convex set $\Theta\subset \mathbb{R}^d$ and $x\in \mathbb{R}^d$ we denote by $\ensuremath{\text{\rm Proj}}_\Theta(x) = {\rm argmin} \{ \|z-x\| : z \in \Theta\}$ the Euclidean projection of $x$ to $\Theta$. We assume everywhere that $T\ge 2$. We denote by ${\cal F_\beta}(L)$ the class of functions with H\"older smoothness $\beta$ (inequality~\eqref{eq:Hclass} below). Recall that $f$ is $\alpha$-strongly convex for some $\alpha>0$ if, for any $x,y \in \mathbb{R}^d$ it holds that $f(y) \geq f(x) + \langle \nabla f(x),y-x\rangle + \frac{\alpha}{2} \|x-y\|^2$.
We further denote by ${\cal F_{\alpha,\beta}}(L)$ the class of all $\alpha$-strongly convex functions belonging to ${\cal F_\beta}(L)$.

{\bf Organization.}
We start in Section~\ref{sec:BM}~with some preliminary results on the gradient estimator.~Section~\ref{sec:UB}  presents our upper bounds for Algorithm \ref{algo}, both in the constrained
and unconstrained case.
In Section \ref{sec:sec4} we observe that a slight modification of Algorithm 1 can be used to estimated the minimum value (rather than the minimizer) of $f$. Section  \ref{sec:sec4} presents improved upper bounds in the case $\beta=2$.
In Section~\ref{sec:LB} we establish minimax lower bounds.
Finally, Section \ref{sec:PW} contrasts our results with previous work in the literature and discusses future directions of research.

\section{Preliminaries}
\label{sec:BM}
In this section, we give the definitions, assumptions and basic facts that will be used throughout the paper.~For $\beta>0$, let  $\ell$ be the greatest integer strictly less than $\beta$. We denote by ${\cal F_\beta}(L)$ the set of all functions $f:\mathbb{R}^d\to \mathbb{R}$ that are  $\ell$ times differentiable and satisfy, for all $x,z \in \Theta$ the H\"older-type condition
\begin{equation}
\bigg|f(z)-\sum_{0\le |m|\leq \ell} \frac{1}{m!}D^{m}f(x)(z-x)^{m} \bigg|\leq L \|z-x\|^{\beta},
\label{eq:Hclass}
\end{equation}
where $L >0$, the sum is over the multi-index $m=(m_{1},...,m_{d}) \in \mathbb{N}^d$, we used the notation $m!=m_{1}! \cdots m_{d}!$, $|m|=m_{1}+ \cdots+m_{d}$, and we defined
\[
D^{m}f(x)\nu^{m} = \frac{\partial ^{|m|}f(x)}{\partial ^{m_{1}}x_{1} \cdots\partial ^{m_{d}}x_{d}}\nu_{1}^{m_{1}} \cdots \nu_{d}^{m_{d}}, \quad \forall \nu=(\nu_1,\dots, \nu_d) \in \mathbb{R}^{d}.
\]
In this paper, we assume that the gradient estimator defined by Algorithm \ref{algo} uses a kernel function $K:[-1,1]\to \mathbb{R}$ satisfying
\begin{equation}
\int K(u) du =0,  \int u K(u) du =1,  \int u^j K(u) du =0, \ j=2,\dots, \ell,~
\int |u|^{\beta} |K(u)| du  <\infty.
\end{equation}
Examples of such kernels obtained as weighted sums of Legendre polynomials are given in \cite{PT90} and further discussed in \cite{BP2016}.

\begin{assumption}
\label{ass1}
It holds, for all $t \in \{1,\dots,T\}$,   that: \text{(i)} the random variables $\xi_t$ and $\xi_t'$ are independent from $\zeta_t$ and from $r_t$, and the random variables $\zeta_t$ and $r_t$ are independent; \text{(ii)} $ \mathbb{E} [\xi_t^2]\le \sigma^2,$ and  $\mathbb{E} [(\xi_t')^2]\le \sigma^2$, where $\sigma\ge 0$.
\end{assumption}
Note that we do not assume $\xi_t$ and $\xi_t'$ to have zero mean. Moreover, they can be non-random and no independence between noises on different steps is required, so that the setting can be considered as adversarial. Having such a relaxed set of assumptions is possible because of randomization that, for example, allows the proofs go through without assuming the zero mean noise.

We will also use the following assumption.

\begin{assumption}\label{ass:lip}
Function $f:\mathbb{R}^d\to \mathbb{R}$ is 2-smooth, that is, differentiable on $\mathbb{R}^d$ and such that $\|\nabla f(x)-\nabla f(x')\|\leq {\bar L} \|x-x'\|$ for all $x,x' \in \mathbb{R}^d$, where $\bar L>0$.
\end{assumption}
It is easy to see that this assumption implies that
$f\in {\cal F}_2(\bar L/2)$.
The following lemma gives a bound on the bias of the gradient estimator.
\begin{restatable}{lemma}{firstlemma}
\label{lem:1}
Let $f \in {\cal F}_\beta(L)$, with $\beta > 1$ and let Assumption \ref{ass1} hold. Let ${\hat g}_t$ and $x_t$ be defined by Algorithm \ref{algo} and let $\kappa_\beta = \int |u|^{\beta} |K(u)| du.$ Then
\begin{equation}
\|\mathbb{E} [{\hat g}_t \hspace{.03truecm}|\hspace{.03truecm} x_t] - \nabla f(x_t)\| \leq  \kappa_\beta L  d h_t^{\beta-1}.
\label{eq:}
\end{equation}
\end{restatable}
{If $K$ be a weighted sum of Legendre polynomials, $\kappa_\beta \le 2\sqrt{2}\beta$, with $\beta \geq 1$ (see e.g.,  \cite[Appendix~A.3]{BP2016}).}

The next lemma provides a bound on the stochastic variability of the estimated gradient
by controlling its second moment.
\begin{restatable}{lemma}{sndlemma}
\label{lem:2unco}
Let Assumption \ref{ass1} hold, let ${\hat g}_t$ and $x_t$ be defined by Algorithm \ref{algo} and set $\kappa = \int K^2(u) du$. Then
\vspace{-2mm}
\begin{itemize}
    \item[(i)] If $\Theta \subseteq \mathbb{R}^d$, $\nabla f(x^*)=0$ and Assumption \ref{ass:lip} holds,
    \[
\mathbb{E} [\|{\hat g}_t \|^2 \hspace{.04truecm}|\hspace{.04truecm} x_t] \leq
9\kappa{\bar L}^2  \left(d\|x_t-x^*\|^2 + \frac{d^2h_t^2}{8} \right) + \frac{3\kappa d^2 \sigma^2}{2h_t^2},
\]
\vspace{-2mm}
    \item[(ii)] If $f\in \mathcal{F}_2(L)$ and $\Theta$ is a closed convex subset of $\mathbb{R}^d$ such that $\displaystyle{\max_{x \in \Theta}} \| \nabla f(x)\|\le G$,  then
    \[
\mathbb{E} [\|{\hat g}_t \|^2 \hspace{.04truecm}|\hspace{.04truecm} x_t] \leq
9\kappa  \left(G^2 d + \frac{{L}^2d^2h_t^2}{2} \right) + \frac{3\kappa d^2 \sigma^2}{2h_t^2}.
\]
\end{itemize}
\end{restatable}
\section{Upper bounds}
\label{sec:UB}
In this section, we provide upper bounds on the cumulative regret and on the optimization error of Algorithm \ref{algo}, {which are defined as
$$\sum_{t=1}^T \mathbb{E} [ f(x_t) - f(x)],$$
and
$$\mathbb{E} [ f({\hat x}_T) - f(x^*)],$$
 respectively, where $x \in \Theta$ and ${\hat x}_T$ is an estimator after $T$ queries. Note that the provided upper bound for cumulative regret is valid for any $x \in \Theta$}.

First we  consider Algorithm \ref{algo} when the convex set $\Theta$ is bounded (constrained case).
\begin{restatable}{theorem}{thdtheorem}{\em (Upper Bound, Constrained Case.)}
\label{thm:1}
	Let $f \in {\cal F}_{\alpha, \beta}(L)$ with $\alpha,L>0$ and $\beta \geq 2$.
	Let Assumptions \ref{ass1} and \ref{ass:lip} hold and let $\Theta$ be a convex compact subset  of $\mathbb{R}^d$. Assume that $\max_{x\in\Theta}\|\nabla f(x)\|\le G $.
	If $\sigma>0$ then the cumulative regret of Algorithm \ref{algo} with \[h_t = \left(\frac{3\kappa\sigma^2}{2(\beta-1)(\kappa_\beta L)^2} \right)^{\frac{1}{2\beta}} t^{-\frac{1}{2\beta}},~~~\eta_t=\frac{2}{\alpha t},~~~t= 1,\dots,T~~~~~~~~\]  satisfies
	\begin{equation}
	\label{eq:bound0}
	\forall x\in \Theta: \ \sum_{t=1}^T \mathbb{E} [ f(x_t) - f(x)]\le \frac{1}{\alpha}
	\left(
d^2 \Big(A_1 T^{{1}/{\beta}}+A_2 \Big)
+  A_3  d \log T
\right),
	\end{equation}
	where $A_1={3\beta}(\kappa\sigma^2)^{\frac{\beta-1}{\beta}}
(\kappa_\beta L)^{\frac{2}{\beta}}$, $A_2=\bar c {\bar L^2} (\sigma/L)^{\frac{2}{\beta}}+{9\kappa G^2}/{d}$ with constant $\bar c>0$ depending only on $\beta$, and $A_3 = 9 \kappa G^2$.
	The optimization error of averaged estimator ${\bar x}_T=\frac{1}{T}\sum_{t=1}^T x_t$ satisfies
	\begin{equation}
	\label{eq:bound1}
	\mathbb{E} [ f({\bar x}_T) - f(x^*)]\leq \frac{1}{\alpha}
\left(  d^2\left(\frac{A_1}{T^{\frac{\beta-1}{\beta}}}+ \frac{A_2}{T} \right)+ A_3 \,\frac{d\log T}{T} \right),
		\end{equation}
		where
	$x^*= \argmin_{x\in\Theta}f(x)$.
			If $\sigma=0$, then the cumulative regret and the optimization error of Algorithm \ref{algo} with any $h_t$ chosen small enough and $\eta_t=\frac{2}{\alpha t}$ satisfy the bounds \eqref{eq:bound0} and \eqref{eq:bound1}, respectively, with $A_1=0, A_2={9\kappa G^2}/{d}$ and $A_3=10\kappa G^2$.
\end{restatable}
\begin{proof}[Proof sketch]
We use the definition of Algorithm \ref{algo} and strong convexity of $f$ to obtain an upper bound for $\sum_{t=1}^{T}\mathbb{E} [ f(x_t) - f(x)  \hspace{.03truecm}|\hspace{.03truecm} x_t]$, which depends on the bias term $\sum_{t=1}^{T}\|\mathbb{E} [{\hat g_t}  \hspace{.03truecm}|\hspace{.03truecm} x_t]-\nabla f(x_t)\|$ and on the stochastic error term $\sum_{t=1}^{T}\mathbb{E} [\|{\hat g}_t\|^2]$.
By substituting $h_{t}$ (that is derived from balancing the two terms) and $\eta_{t}$ in Lemmas \ref{lem:1} and \ref{lem:2unco}  we obtain upper bounds for $\sum_{t=1}^{T}\|\mathbb{E} [{\hat g_t}  \hspace{.03truecm}|\hspace{.03truecm} x_t]-\nabla f(x_t)\|$ and $\sum_{t=1}^{T}\mathbb{E} [\|{\hat g}_t\|^2]$ that imply the desired upper bound for $\sum_{t=1}^{T}\mathbb{E} [ f(x_t) - f(x)  \hspace{.03truecm}|\hspace{.03truecm} x_t]$ due to a recursive argument in the spirit of \cite{bartlett-hazan-rakhlin}.
\end{proof}

In the non-noisy case ($\sigma=0$) we get the rate $\frac{d}{\alpha} \log T$
for the cumulative regret,  and  $\frac{d}{\alpha} \frac{\log T}{T}$ for the optimization error.
In what concerns the optimization error, this rate is not optimal since one can achieve much faster rate under strong convexity \cite{NS17}. However,  for the cumulative regret in our derivative-free setting it remains an open question whether the result of Theorem \ref{thm:1}
can be improved. Previous papers on derivative-free online methods with no noise \citep{agarwal2010,duchi2015,flaxman2004} provide slower rates than $({d}/{\alpha}) \log T$. The best known so far is
$({d^2}/{\alpha})\log T$, cf. \cite[Corollary 5]{agarwal2010}.
We may also notice that the cumulative regret bounds of Theorem \ref{thm:1} trivially extend to the case when we query functions $f_t$ depending on $t$ rather than a single $f$. Another immediate fact is that on the r.h.s. of inequalities \eqref{eq:bound0} and \eqref{eq:bound1} we can take the minimum with $GBT$ and $GB$, respectively, where $B$ is the Euclidean diameter of $\Theta$. Finally, the factor $\log T$ in the bounds for the optimization error can be eliminated by considering averaging from $T/2$ to $T$ rather than from 1 to $T$, in the spirit of \cite{rakhlin2011making}. We refer to Appendix \ref{app:C} for the details and proofs of these facts.

We now study the performance of Algorithm \ref{algo} when $\Theta=\mathbb{R}^d$. In this case we make the following choice for the parameters $h_t$ and $\eta_t$ in Algorithm \ref{algo}:
\begin{equation}
\begin{aligned}
h_t &=T^{-\frac{1}{2\beta}},~~~\eta_t = \frac{1}{\alpha T},~~~t= 1,\dots,T_0, \\
h_t &= t^{-\frac{1}{2\beta}},~~~~~\eta_t = \frac{2}{\alpha t},~~~~t= T_0+ 1,\dots,T,
\end{aligned}
\label{eq:hhh}
\end{equation}
where $T_0 = \max \left\{k \geq 0 : C_1 {\bar L}^2 d > {\alpha^2 k}/{2} \right\}$ and $C_1$ is a positive constant\footnote{If $T_0=0$ the algorithm does not use \eqref{eq:hhh}. Assumptions of Theorem \ref{thm:2} are such that condition $T>T_0$ holds.} depending only on the kernel $K(\cdot)$ (this is defined in the proof of Theorem \ref{thm:2} in Appendix \ref{app:B}) and recall ${\bar L}$ is the Lipschitz constant on the gradient $\nabla f$. Finally, define the estimator
\begin{equation}
\label{eq:ppmean}
{\bar x}_{T_0,T} = \frac{1}{T-T_0} \sum_{t=T_0+1}^T x_t.
\end{equation}
\begin{restatable}{theorem}{fiftheorem}{\em (Upper Bounds, Unconstrained Case.)}
\label{thm:2}
	Let $f\in {\cal F} _{\alpha,\beta}(L)$ with $\alpha, L>0$ and $\beta \geq 2$. Let Assumptions \ref{ass1} and \ref{ass:lip} hold.  Assume also that  $\alpha > \sqrt{{C_*d}/{T}}$, where $C_* > 72 \kappa {\bar L}^2$. Let $x_t$'s be the updates of Algorithm \ref{algo} with $\Theta=\mathbb{R}^d$, $h_t$ and $\eta_t$ as in \eqref{eq:hhh} and a non-random $x_1 \in \mathbb{R}^d$.
	Then the estimator defined by \eqref{eq:ppmean} satisfies
\begin{equation}
\label{eq:pmean}
\mathbb{E}  [f({\bar x}_{T_0,T}) - f(x^*)] \leq
C \kappa {\bar L}^2 \frac{d}{\alpha T} \|x_1- x^*\|^2
 + C
\frac{d^2}{\alpha} \left(
(\kappa_\beta L)^2 + \kappa \big({\bar L}^2
+ \sigma^2\big)
\right)
T^{-\frac{\beta-1}{\beta}}
\end{equation}
where $C>0$ is a constant  depending only on $\beta$
and $x^*=\argmin_{x\in \mathbb{R}^d}f(x)$.
\end{restatable}
\begin{proof}[Proof sketch]
As in the proof of Theorem \ref{thm:1}, we apply Lemmas \ref{lem:1} and \ref{lem:2unco}. But we can only use Lemma \ref{lem:2unco}(i) and not Lemma \ref{lem:2unco}(ii) and thus the bound on the stochastic error now involves  $\norm{x_{t}-x^{*}}^{2}$. So, after taking expectations, we need to control an additional term containing  $r_t=\mathbb{E}[\norm{x_{t}-x^{*}}^{2}]$. However, the issue concerns only small $t$ ($t\le T_0\sim d^2/\alpha$) since for bigger $t$ this term is compensated due to the strong convexity with parameter
$\alpha > \sqrt{{C_*d}/{T}}$.
This motivates the method where we use the first $T_0$ iterations to get a suitably good (but not rate optimal) bound on $r_{T_0+1}$ and then proceed analogously to Theorem \ref{thm:1} for iterations $t\ge T_0+1$.
\end{proof}

\section{Estimation of $f(x^*)$}
\label{sec:sec4}
In this section, we apply the above results to estimation of the minimum value  $f(x^*)=\min_{x\in\Theta}f(x)$ for functions $f$ in the class $\mathcal{F}_{\alpha,\beta}(L)$. The literature related to this problem assumes that $x_t$'s are either i.i.d. with density bounded away from zero on its support \cite{tsy1990} or $x_t$'s are chosen sequentially \cite{mokkadem,belitser}. In the fist case, from the results in \cite{tsy1990} one can deduce that $f(x^*)$ cannot be estimated better than at the slow rate $T^{-\beta/(2\beta+d)}$.
For the second case, which is our setting, the best result so far is obtained in \cite{belitser}.
The estimator of $f(x
^*)$ in \cite{belitser} is defined via a multi-stage procedure whose complexity increases exponentially with the dimension $d$ and it is shown to achieve  (asymptotically, for $T$ greater than an exponent of $d$) the $c(d, \alpha)/\sqrt{T}$ rate for functions in $\mathcal{F}_{\alpha,\beta}(L)$ with $\beta> 2$. Here, $c(d,\alpha)$ is some constant depending on $d$ and $\alpha$ in an unspecified way.

Observe that $f(\bar x_T)$ is not an estimator since it depends on the unknown $f$, so Theorem \ref{thm:1} does not provide a result about estimation of $f(x^*)$.
In this section, we show that using the computationally simple Algorithm \ref{algo} and making one more query per step (that is, having three queries per step in total) allows us to achieve the $1/\sqrt{T}$ rate for all $\beta\ge 2$ with no dependency on the dimension in the main term. Note that the $1/\sqrt{T}$ rate cannot be improved. Indeed, one cannot estimate $f(x^*)$ with  a better rate even using the ideal but non-realizable oracle that makes all queries at point $x^*$. {That is, even if $x^*$ is known and we sample $T$ times $f(x^*) + \xi_t$ with independent centered variables  $\xi_t$, the error is still of the order $1/\sqrt{T}$.}

In order to construct our estimator, at any step $t$ of Algorithm \ref{algo} we make along with $y_t$ and $y_t'$ the third query $y''_t=f(x_t)+\xi''_t$, where $\xi''_t$ is some noise and $x_t$ are the updates of Algorithm \ref{algo}. We estimate $f(x^*)$ by
$
\hat M = \frac{1}{T}\sum_{t=1}^T y''_t.
$
The properties of estimator $\hat M$ are summarized in the next theorem, which is an immediate corollary of Theorem \ref{thm:1}.
\begin{restatable}{theorem}{fththeorem}
\label{thm:min}
	Let the assumptions of Theorem \ref{thm:1} be satisfied. Let $\sigma>0$ and assume that $(\xi''_t)_{t=1}^T$ are independent random variables with $\mathbb{E}[\xi''_t]=0$ and $\mathbb{E}[(\xi''_t)^2]\le \sigma^2$ for $t=1,\dots, T$.
	If $f$ attains its minimum at point
	$x^*\in \Theta$, then
	\begin{equation}
	\label{eq:boundM}
	\mathbb{E} \lvert \hat M - f(x^*) \rvert \leq \frac{\sigma}{T^{\frac{1}{2}}}+ \frac{1}{\alpha}
\left(  d^2\left(\frac{A_1}{T^{\frac{\beta-1}{\beta}}}+ \frac{A_2}{T}\right)+
A_3\frac{d\log T}{T} \right).
		\end{equation}
		\end{restatable}
\begin{remark}\label{remark:min}
With three queries per step, the risk (error) of the oracle that makes all queries at point $x^*$ does not exceed  $\sigma/\sqrt{3T}$. Thus, for $\beta>2$ the estimator $\hat M$ achieves asymptotically as $T\to \infty$ the oracle risk up to a numerical constant factor. We do not obtain such a sharp property for $\beta=2$, in which case the remainder term in Theorem \ref{thm:min} accounting for the accuracy of Algorithm \ref{algo} is of the same order as the main term $\sigma/\sqrt{T}$.
\end{remark}
Note that in Theorem \ref{thm:min} the noises $(\xi''_t)_{t=1}^T$ are assumed to be independent and zero mean random variables, which is essential to obtain the $1/\sqrt{T}$ rate. Nevertheless, we do not require independence between the noises $(\xi''_t)_{t=1}^T$ and the noises in the other two queries $(\xi_t)_{t=1}^T$ and $(\xi'_t)_{t=1}^T$.
Another interesting point is that for $\beta=2$ the third query is not needed and $f(x^*)$ is estimated with the $1/\sqrt{T}$ rate either by
$
\hat M = \frac{1}{T}\sum_{t=1}^T y_t
$
or by
$
\hat M = \frac{1}{T}\sum_{t=1}^T y'_t.
$
This is an easy consequence of the above argument, the property \eqref{lem6.4_2} -- see Lemma \ref{lem6.4} in the appendix -- which is specific for the case $\beta=2$, and the fact that the optimal choice of $h_t$ is of order $t^{-1/4}$ for $\beta=2$.

\section{Improved bounds for $\beta=2$}
\label{sec:beta2}

In this section,  we consider the case $\beta=2$ and obtain improved bounds that scale as $d$ rather than $d^2$ with the dimension in the constrained optimization setting analogous to Theorem \ref{thm:1}. First note that for $\beta=2$ we can simplify the algorithm.  The use of kernel $K$ is redundant when $\beta=2$, and therefore in this section we define the approximate gradient as
\begin{equation}\label{gt}
   \hat g_t = \frac{d}{2h_t} (y_t- y'_t) \zeta_t,
\end{equation}
where $y_t=f(x_t+h_t \zeta_t)+\xi_t$ and $y'_t=f(x_t-h_t  \zeta_t) + \xi_t'$. A well-known observation that goes back to \cite{NY1983} consists in the fact that $\hat g_t$ defined in \eqref{gt} is an unbiased estimator of the gradient at point $x_t$ of the surrogate function $\hat f_t$ defined by
\[
\hat f_t(x)= \mathbb{E} f(x+h_t \tilde \zeta),\quad \forall x\in\mathbb{R}^d,
\]
where the expectation  $\mathbb{E}$ is taken   with respect to the random vector $\tilde \zeta$ uniformly distributed on the unit ball $B_d=\{u\in \mathbb{R}^d:\|u\|\le 1 \}$.  The properties of the surrogate $\hat f_t$  are described in Lemmas \ref{lem6.2} and \ref{lem6.4} presented in the appendix.

The improvement in the rate that we get for $\beta=2$ is due to the fact that we can consider Algorithm \ref{algo} with  $\hat g_t$ defined in \eqref{gt} as the SGD for the surrogate function. Then the bias of approximating $f$ by $\hat f_t$ scales as $h_t^2$, which is smaller than the squared bias of
approximating the gradient arising in the proof of Theorem \ref{thm:1} that scales as $d^2h_t^{2(\beta-1)}=d^2h_t^2$ when $\beta=2$. On the other hand, the stochastic variability terms are the same for both methods of proof. This explains the gain in dependency on $d$. However, this technique does not work for $\beta>2$ since then the error of approximating $f$ by $\hat f_t$, which is of the order $h_t
^\beta$ (with $h_t$ small), becomes too large compared to the bias $d^2h_t^{2(\beta-1)}$ of Theorem \ref{thm:1}.
\begin{restatable}{theorem}{eiththeorem}
\label{th6.3}
Let $f\in \mathcal{F}_{\alpha, 2}(L)$ with $\alpha,L>0$.
	Let Assumption \ref{ass1}  hold and let $\Theta$ be a convex compact subset  of $\mathbb{R}^d$. Assume that $\max_{x\in\Theta}\|\nabla f(x)\|\le G $.
If $\sigma>0$ then for the updates $x_t$ as in item 3 of Algorithm \ref{algo} with $\hat g_t$ defined in \eqref{gt} and parameters $h_t= \left(\frac{3d^2\sigma^2}{4L\alpha t+9L^2d^2}\right)^{1/4}$ and $\eta_t = \frac{1}{\alpha t}$ we have
\begin{equation}\label{eq1:th6.3}
\forall x\in \Theta: \ \ \mathbb{E} \sum_{t=1}^{T} \big(f(x_t) - f(x)\big)
\le  \min\left(GBT, 2\sqrt{3L}\sigma
\frac{d}{\sqrt{\alpha}} \sqrt{T}+ A_4\frac{d^2}{\alpha} \log T\right),
\end{equation}
where $B$ is the Euclidean diameter of $\Theta$ and
$A_4=
6.5 L\sigma+ 22 G^2/d$.
Moreover, if $x^*= \argmin_{x\in\Theta}f(x)$ the optimization error of averaged estimator ${\bar x}_T=\frac{1}{T}\sum_{t=1}^T x_t$ is bounded as
	\begin{equation}
	\label{eq2:th6.3}
	\mathbb{E} [ f({\bar x}_T) - f(x^*)]\leq
	\min\left(GB, 2\sqrt{3L}\sigma
 \frac{d}{\sqrt{\alpha T}}+ A_4 \frac{d^2}{\alpha} \frac{\log T}{T}\right).
		\end{equation}
Finally, if $\sigma=0$, then the cumulative regret of the same procedure $x_t$ with any $h_t>0$ chosen small enough and $\eta_t=\frac{1}{\alpha t}$ and the optimization error of its averaged version are of the order $\frac{d^2}{\alpha} \log T$ and $\frac{d^2}{\alpha} \frac{\log T}{T}$, respectively.                         \end{restatable}
Note that the terms $\frac{d^2}{\alpha} \log T$ and $\frac{d^2}{\alpha} \frac{\log T}{T}$ appearing in these bounds can be improved to $\frac{d}{\alpha} \log T$ and $\frac{d}{\alpha} \frac{\log T}{T}$ at the expense of assuming that the norm $\|\nabla f\|$ is uniformly bounded by $G$ not only on $\Theta$ but also on a large enough Euclidean neighborhood of $\Theta$. Moreover, the $\log T$ factor in the bounds for the optimization error can be eliminated by considering averaging from $T/2$ to $T$ rather than from 1 to $T$ in the spirit of  \cite{rakhlin2011making}. We refer to Appendix \ref{app:C}
for the details and proofs of these facts. A major conclusion is that, when $\sigma>0$ and we consider the optimization error,  those terms are negligible with respect to $d/\sqrt{\alpha T}$ and thus an attainable rate is $\min(1,d/\sqrt{\alpha T})$.

We close this section by noting, in connection with the bandit setting, that the bound \eqref{eq1:th6.3} extends straightforwardly (up to a change in numerical constants) to the cumulative regret of the form $\mathbb{E}\sum_{t=1}^{T} \big(f_t(x_t\pm h_t\zeta_t) - f_t(x)\big)$, where the losses are measured at the query points and $f$ depends on $t$. This fact follows immediately from the proof of Theorem \ref{th6.3} presented in the appendix and the property \eqref{lem6.4_2}, see Lemma \ref{lem6.4} in the appendix.

\section{Lower bound}
\label{sec:LB}
In this section we prove a minimax lower bound on the optimization error over all sequential strategies that allow the query points depend on the past. For $t=1,\dots, T$, we assume that $y_t=f(z_t)+\xi_t$ and we consider strategies of choosing the query points as $z_t= \Phi_t(z_1^{t-1},y_1^{t-1})$ where $\Phi_t$ are Borel functions and $z_1\in \mathbb{R}^d$
is any random variable. We denote by $\Pi_T$ the set of all such strategies. The noises $\xi_1,\dots,\xi_T$ are assumed in this section to be independent with cumulative distribution function $F$ satisfying the condition
\begin{equation}\label{distribution}
\int \log\big(dF(u)/dF(u+v)\big)dF(u)\leq I_{0}v^{2}, \quad |v|< v_{0}
\end{equation}
for some $0<I_{0}<\infty$, $0<v_{0}\leq \infty$. {Using the second order expansion of the logarithm w.r.t. $v$, one can verify that this assumption is satisfied when $F$ has a smooth enough density with finite Fisher information.} For example, for Gaussian distribution $F$
this condition holds with $v_0=\infty$. Note that the class $\Pi_T$ includes the sequential strategy of Algorithm \ref{algo} that corresponds to taking $T$ as an even number, and choosing $z_t=x_t+\zeta_tr_t$ and $z_t=x_t-\zeta_tr_t$  for even $t$ and odd $t$, respectively. The presence of the randomizing sequences $\zeta_t, r_t$ is not crucial for the lower bound. Indeed, Theorem \ref{lb} below is valid conditionally on any randomization, and thus the lower bound remains valid when taking expectation over the randomizing distribution.

\begin{restatable}{theorem}{lowerB}
\label{lb}
Let $\Theta=\{x\in\mathbb{R}^d: \norm{x}\le 1\}$. For $\alpha, L>0$,$ \beta\ge 2$, let  $\mathcal{F}'_{\alpha,\beta}$ denote the set of functions $f$ that attain their minimum over $\mathbb{R}^d$ in $\Theta$ and belong to $\mathcal{F}_{\alpha,\beta}(L)\cap \{f: \max_{x\in\Theta}\|\nabla f(x)\|\le G\}$,  where $G> 2\alpha$. Then for any strategy in the class $\Pi_T$ we have
\begin{equation}\label{eq1:lb}
\sup_{f \in \mathcal{F}'_{\alpha,\beta}}\mathbb{E}\big[f(z_T)-\min_{x}f(x)\big]\geq C\min\Big(\max(\alpha, T^{-1/2+1/\beta}), \frac{d}{\sqrt{T}}, \,\frac{d}{\alpha}T^{-\frac{\beta-1}{\beta}}\Big),
\end{equation}
and
\begin{equation}\label{eq2:lb}
\sup_{f \in \mathcal{F}'_{\alpha,\beta}}\mathbb{E}\big[\norm{z_{T}-x^{*}(f)}^{2}\big]\geq C\min\Big(1, \frac{d}{T^{\frac{1}{\beta}}}, \,\frac{d}{\alpha^{2}}T^{-\frac{\beta-1}{\beta}}\Big),
\end{equation}
 where  $C>0$ is a constant that does not depend  of $T,d$, and $\alpha$, and $x^*(f)$ is the minimizer of $f$ on $\Theta$.
\end{restatable}
The proof is given in  Appendix \ref{app:B}. {It extends the proof technique of Polyak and Tsybakov [28], by applying it to more than two probe functions. The proof takes into account dependency on the dimension $d$, and on $\alpha$. The final result is obtained by applying Assouad’s Lemma, see e.g. \cite{Tsybakov09}.}

We stress that the condition $G> 2\alpha$ in this theorem is necessary. It should always hold if the intersection $\mathcal{F}_{\alpha,\beta}(L)\cap \{f: \max_{x\in\Theta}\|\nabla f(x)\|\le G\}$ is not empty. Notice also that the threshold $T^{-1/2+1/\beta}$ on the strong convexity parameter $\alpha$ plays an important role in bounds \eqref{eq1:lb} and \eqref{eq2:lb}. Indeed, for $\alpha$ below this threshold, the bounds start to be independent of $\alpha$.  Moreover, in this regime, the rate of \eqref{eq1:lb} becomes $\min(T^{1/\beta},d)/\sqrt{T}$, which is asymptotically  $d/\sqrt{T}$ and thus not better as function of $T$ than the rate attained for
zero-order minimization of simply convex functions \cite{agarwal2011,belloni2015escaping}. Intuitively, it seems reasonable that $\alpha$-strong convexity should be of no added value for very small $\alpha$. Theorem \ref{lb} allows us to quantify exactly how small such $\alpha$ should be. Also, quite naturally, the threshold becomes smaller when the smoothness $\beta$ increases.\\
Finally note that for $\beta=2$ the lower bounds \eqref{eq1:lb} and \eqref{eq2:lb} are, in the interesting regime of large enough $T$, of order $d/(\max(\alpha, 1)\sqrt{T})$ and $d/(\max(\alpha^2, 1)\sqrt{T})$, respectively. This highlights the near minimax optimal properties of Algorithm \ref{algo} in the setting of Theorem \ref{th6.3}.

\section{Discussion and related work}
\label{sec:PW}
There is a great deal of attention to zero-order feedback stochastic optimization and convex bandits problems in the recent literature. Several settings are studied: (i) deterministic in the sense that the queries contain no random noise and we query functions $f_t$ depending on $t$ rather than $f$ where $f_t$ are Lipschitz or 2-smooth \cite{flaxman2004,agarwal2010,Nesterov2011,NS17,saha2011,Shamir17}; (ii) stochastic with two-point feedback where the two noisy evaluations are obtained with the same noise and the noisy functions are Lipschitz or 2-smooth \cite{Nesterov2011,NS17,duchi2015} (this setting does not differ much from (i) in terms of the analysis and the results); (iii) stochastic, where the noises $\xi_i$ are independent zero-mean random variables  \cite{Fabian,PT90,dippon,agarwal2011,Shamir13,BP2016, jamieson2012,bartlett2018simple,Locatelli}. In this paper, we considered a setting, which is more general than (iii) by allowing for adversarial noise (no independence or zero-mean assumption in contrast to (iii), no Lipschitz assumption  in contrast to settings (i) and (ii)), which are both covered by our results when the noise is set to zero.

One part of our results are bounds on the cumulative regret, cf. \eqref{eq:bound0} and \eqref{eq1:th6.3}. We emphasize that they remain trivially valid if the queries are from $f_t$ depending on $t$ instead of $f$, and thus cover the setting (i). To the best of our knowledge, there were no such results in this setting previously, except for \cite{BP2016} that gives  bounds with suboptimal dependency on $T$ in the case of classical (non-adversarial) noise. In the non-noisy case, we get bounds on the cumulative regret with faster rates than previously known for the setting (i). It remains an open question whether these bounds can be improved.

The second part of our results dealing with the optimization error $\mathbb{E}[f(\bar x_T)-f(x^*)]$ is closely related to the work on derivative-free stochastic optimization under strong convexity and smoothness assumptions initiated in \cite{Fabian,PT90} and more recently developed in \cite{dippon,jamieson2012,Shamir13,BP2016}.
It was shown in \cite{PT90} that the minimax optimal rate for $f\in\mathcal{F}_{\alpha,\beta}(L)$ scales as $c(\alpha,d)T^{-{(\beta-1)}/{\beta}}$, where $c(\alpha,d)$ is an unspecified function of $\alpha$ and $d$ (for $d=1$ an upper bound of the same order was earlier established in
\cite{Fabian}). The issue of establishing non-asymptotic fundamental limits as function of the main parameters of the problem ($\alpha$, $d$ and $T$) was first addressed in \cite{jamieson2012} giving a lower bound $\Omega(\sqrt{d/T})$ for $\beta=2$. This was improved to $\Omega(d/\sqrt{T})$ when $\alpha\asymp 1$ by Shamir \cite{Shamir13} who conjectured that the rate $d/\sqrt{T}$ is optimal for $\beta=2$, which indeed follows from our Theorem \ref{th6.3} (although \cite{Shamir13} claims the optimality as proved fact by referring to results in \cite{agarwal2010}, such results cannot be applied in setting (iii) because the noise cannot be considered as Lipschitz).
A result similar to Theorem \ref{th6.3} is stated without proof in Bach and Perchet \cite[Proposition 7]{BP2016} but not for the cumilative regret and with  a suboptimal rate in the non-noisy case. For integer $\beta\ge 3$, Bach and Perchet \cite{BP2016} present explicit upper bounds as functions of $\alpha$, $d$ and $T$ with, however, suboptimal dependency on $T$ except for their Proposition 8 that is problematic (see Appendix \ref{app:D} for the details).  Finally, by slightly modifying the proof of Theorem \ref{thm:1} we get that the estimation risk $\mathbb{E}\big[\norm{\bar x_{T}-x^{*}}^{2}\big]$ is $O((d^2/\alpha
^2)T^{-{(\beta-1)}/{\beta}})$, which is to within factor $d$ of the main term in the lower bound \eqref{eq2:lb} (see Appendix \ref{app:C} for details).

The lower bound in Theorem  \ref{lb} is, to the best of our knowledge, the first result providing non-asymptotic fundamental limits
under general configuration of $\alpha$, $d$ and $T$. The known lower bounds \cite{PT90,jamieson2012,Shamir13} either give no explicit dependency on  $\alpha$ and $d$,
or treat the special case $\beta=2$ and $\alpha\asymp 1$. Moreover, as an interesting consequence of our lower bound we find that, for small strong convexity parameter $\alpha$ (namely, below the $T^{-1/2+1/\beta}$ threshold), the best achievable rate cannot be substantially faster
than for simply convex functions, at least for moderate dimensions. Indeed, for such small $\alpha$, our lower bound is asymptotically $\Omega(d/\sqrt{T})$ independently of the smoothness index $\beta$ and on $\alpha$, while the achievable rate for convex functions is shown to be $d^{16}/\sqrt{T}$ in \cite{agarwal2011} and improved to $d^{3.75}/\sqrt{T}$ in  \cite{belloni2015escaping} (both up to log-factors). The gap here is only in the dependency on the dimension. Our results imply that for $\alpha$ above the $T^{-1/2+1/\beta}$ threshold, the gap between upper and lower bounds is much smaller. Thus, our upper bounds in this regime scale as $(d^2/\alpha)T^{-{(\beta-1)}/{\beta}}$ while the lower bound of
Theorem  \ref{lb} is of the order $\Omega\big((d/\alpha)T^{-{(\beta-1)}/{\beta}}\big)$; moreover for $\beta=2$, upper and lower bounds match in the dependency on $d$.

We hope that our work will stimulate further study at the intersection of zero-order optimization and convex bandits in machine learning. An important open problem is to study novel algorithms which match our lower bound simultaneously in all main parameters. For example a class of algorithms worth exploring are those using memory of the gradient in the spirit of Nesterov accelerated method. Yet another important open problem is to study lower bounds for the regret in our setting. Finally, it would be valuable to study extensions of our work to locally strongly convex functions.

\section*{Broader impact}
The present work improves our understanding of zero-order optimization methods in specific scenarios in which the underlying function we wish to optimize has certain regularity properties. We believe that a solid theoretical foundation is beneficial to the development of practical machine learning and statistical methods. We expect no direct or indirect ethical risks from our research.

\begin{ack}
We would like to thank Francis Bach, Vianney Perchet, Saverio Salzo, and Ohad Shamir for helpful discussions. The first and second authors were  partially supported by SAP SE. The research of A.B. Tsybakov is supported by a grant of the French National Research Agency (ANR), “Investissements d'Avenir” (LabEx Ecodec/ANR-11-LABX-0047).
\end{ack}