\subsubsection*{\bibname}}

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\usepackage{pifont}
\usepackage{subfigure}
\newcommand{\mathbf{E}}{\mathbf{E}}
\newcommand{\mathbb P}{\mathbf{P}}
\newcommand{\mathbb R}{\mathbb{R}}
\newcommand{\stackrel{\text{def}}{=}}{\stackrel{\text{def}}{=}}
\newcommand{\ve}[2]{\left\langle #1 , #2 \right\rangle}
\def\<#1,#2>{\left\langle #1,#2\right\rangle}

\newcommand{\oldstuff}[1]{ {\color{blue} #1}}
\newcommand{\redstuff}[1]{ {\color{red} #1}}

\theoremstyle{plain}
\newtheorem{lemma}{Lemma}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{definition}{Definition}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{assumption}{Assumption}[section]
\newtheorem{corollary}{Corollary}[section]

\newcommand\tagthis{\addtocounter{equation}{1}\tag{\theequation}}
\newcommand{\mathop{\arg\!\min}}{\mathop{\arg\!\min}}
\newcommand{\mathop{\arg\!\max}}{\mathop{\arg\!\max}}
\newcommand{\text{\ding{172}}}{\text{\ding{172}}}
\newcommand{\text{\ding{173}}}{\text{\ding{173}}}
\newcommand{\text{\ding{174}}}{\text{\ding{174}}}
\newcommand{\text{\ding{175}}}{\text{\ding{175}}}
\newcommand{\text{\ding{176}}}{\text{\ding{176}}}
\newcommand{\text{\ding{177}}}{\text{\ding{177}}}
\newcommand{\text{\ding{178}}}{\text{\ding{178}}}
\newcommand{\text{\ding{179}}}{\text{\ding{179}}}
\newcommand{\text{\ding{180}}}{\text{\ding{180}}}
\newcommand{\text{\ding{181}}}{\text{\ding{181}}}

\usepackage[colorinlistoftodos,bordercolor=orange,backgroundcolor=orange!20,linecolor=orange,textsize=scriptsize]{todonotes}

\newcommand{\marina}[1]{\todo[inline]{{\textbf{Marina:} \emph{#1}}}}
\newcommand{\alexander}[1]{\todo[inline]{{\textbf{Alexander:} \emph{#1}}}}
\newcommand{\eduard}[1]{\todo[inline]{{\textbf{Eduard:} \emph{#1}}}}
\newcommand{\pavel}[1]{\todo[inline]{{\textbf{Pavel:} \emph{#1}}}}
\newcommand{\kesha}[1]{\todo[inline]{{\textbf{Kesha:} \emph{#1}}}}

\newcommand{{\cal A}}{{\cal A}}
\newcommand{{\cal B}}{{\cal B}}
\newcommand{{\cal C}}{{\cal C}}
\newcommand{{\cal D}}{{\cal D}}
\newcommand{{\cal E}}{{\cal E}}
\newcommand{{\cal F}}{{\cal F}}
\newcommand{{\cal G}}{{\cal G}}
\newcommand{{\cal H}}{{\cal H}}
\newcommand{{\cal J}}{{\cal J}}
\newcommand{{\cal K}}{{\cal K}}
\newcommand{{\cal L}}{{\cal L}}
\newcommand{{\cal M}}{{\cal M}}
\newcommand{{\cal N}}{{\cal N}}
\newcommand{{\cal O}}{{\cal O}}
\newcommand{{\cal P}}{{\cal P}}
\newcommand{{\cal Q}}{{\cal Q}}
\newcommand{{\cal R}}{{\cal R}}
\newcommand{{\cal S}}{{\cal S}}
\newcommand{{\cal T}}{{\cal T}}
\newcommand{{\cal U}}{{\cal U}}
\newcommand{{\cal V}}{{\cal V}}
\newcommand{{\cal X}}{{\cal X}}
\newcommand{{\cal Y}}{{\cal Y}}
\newcommand{{\cal W}}{{\cal W}}
\newcommand{{\cal Z}}{{\cal Z}}
\newcommand{\mathrm{Var}}{\mathrm{Var}}

\newcommand{{\bf A}}{{\bf A}}
\newcommand{{\bf B}}{{\bf B}}
\newcommand{{\bf C}}{{\bf C}}
\newcommand{{\bf E}}{{\bf E}}
\newcommand{{\bf F}}{{\bf F}}
\newcommand{{\bf G}}{{\bf G}}
\newcommand{{\bf H}}{{\bf H}}
\newcommand{{\bf I}}{{\bf I}}
\newcommand{{\bf J}}{{\bf J}}
\newcommand{{\bf K}}{{\bf K}}
\newcommand{{\bf L}}{{\bf L}}
\newcommand{{\bf M}}{{\bf M}}
\newcommand{{\bf N}}{{\bf N}}
\newcommand{{\bf O}}{{\bf O}}
\newcommand{{\bf P}}{{\bf P}}
\newcommand{{\bf Q}}{{\bf Q}}
\newcommand{{\bf R}}{{\bf R}}
\newcommand{{\bf S}}{{\bf S}}
\newcommand{{\bf T}}{{\bf T}}
\newcommand{{\bf U}}{{\bf U}}
\newcommand{{\bf V}}{{\bf V}}
\newcommand{{\bf W}}{{\bf W}}
\newcommand{{\bf X}}{{\bf X}}
\newcommand{{\bf Y}}{{\bf Y}}
\newcommand{{\bf Z}}{{\bf Z}}

\newcommand{\mathrm{sign}}{\mathrm{sign}}
\newcommand{\omega}{\omega}
\newcommand{\mathbb E}{\mathbf{E}}
\newcommand{\mathbf{V}}{\mathbf{V}}

\def\mathbb R{\mathbb{R}}
\newcommand{\mathbb E}{{\mathbb E}}
\def\mathcal W{\mathcal W}
\def\mathcal C{\mathcal C}
\def\mathcal X{\mathcal X}
\def\mathcal Y{\mathcal Y}
\def\mathcal M_+^1{\mathcal M_+^1}
\def\mathbb R{\mathbb R}
\def\mathbb E{\mathbb E}
\def\mathbb E{\mathbb E}
\def\mathbb P{\mathbb P}
\def\mathbb P{\mathbb P}
\def\mathbb D{\mathbb D}
\def\mathcal{U}{\mathcal{U}}
\def\mathcal{V}{\mathcal{V}}

\def\varepsilon{\varepsilon}
\def\langle{\langle}
\def\rangle{\rangle}
\def\varphi{\varphi}
\def\mathbf{y}{\mathbf{y}}
\def\mathbf{a}{\mathbf{a}}
\def\mathbf{x}{\mathbf{x}}
\def\mathbf{R}{\mathbf{R}}
\def{\mathbf 1}{{\mathbf 1}}

\defy}{\bar{\lambda}{L}{y}{\bar{\lambda}{L}}
\defy}{\bar{\lambda}{g}{y}{\bar{\lambda}{g}}
\def\widetilde{\nabla}{\widetilde{\nabla}}
\defy}{\bar{\lambda}{x}{y}{\bar{\lambda}{x}}
\defy}{\bar{\lambda}{y}{y}{\bar{\lambda}{y}}

\def\boldsymbol{\lambda}{\boldsymbol{\lambda}}
\def\boldsymbol{\zeta}{\boldsymbol{\zeta}}
\def\boldsymbol{\xi}{\boldsymbol{\xi}}
\def\boldsymbol{\eta}{\boldsymbol{\eta}}
\def\boldsymbol{\bar{\lambda}}{\boldsymbol{\bar{\lambda}}}
\def\boldsymbol{\bar{\zeta}}{\boldsymbol{\bar{\zeta}}}
\def\boldsymbol{\bar{\eta}}{\boldsymbol{\bar{\eta}}}
\def\boldsymbol{x}{\boldsymbol{x}}
\def\boldsymbol{\xi}{\boldsymbol{\xi}}
\def\mathtt{p}{\mathtt{p}}
\def\lambda{\lambda}
\def\lambda_{max}(W){\lambda_{max}(W)}
\def\boldsymbol{\boldsymbol}

\def\overline{x}{\overline{x}}
\def\text{cut}{\text{cut}}
\def\text{clip}{\text{clip}}
\def\overline{\Delta}{\overline{\Delta}}
\def\widetilde{\sigma}{\widetilde{\sigma}}

\def\mathbbm{1}{\mathbbm{1}}

\def\eg#1{{\color{blue}#1}}
\def\ag#1{{\color{red}#1}}
\def\md#1{{\color{cyan}#1}}
\def\pd#1{{\color{orange}#1}}

\newcommand{\mathop{\mathrm{prox}}\nolimits}{\mathop{\mathrm{prox}}\nolimits}
\newcommand{\prox_{\gamma R}}{\mathop{\mathrm{prox}}\nolimits_{\gamma R}}
\newcommand{\prox_{j}^{\gamma}}{\mathop{\mathrm{prox}}\nolimits_{j}^{\gamma}}
\newcommand{\prox_{\gamma^k R}}{\mathop{\mathrm{prox}}\nolimits_{\gamma^k R}}
\newcommand{\overline}{\overline}
\newcommand{\sum_{i=1}^n}{\sum_{i=1}^n}

\usepackage{hyperref}

\usepackage{xspace}
\newcommand{\squeeze}{}

\newcommand{\algname}[1]{{\sf  #1}\xspace}
\newcommand{\algnames}[1]{{\sf \small #1}\xspace}

\usepackage{pifont}
\newcommand{\cmark}{{\color{PineGreen}\ding{51}}}%
\newcommand{\xmark}{{\color{BrickRed}\ding{55}}}%

\usepackage{colortbl}
\definecolor{bgcolor}{rgb}{0.8,1,1}
\definecolor{bgcolor2}{rgb}{0.8,1,0.8}

\usepackage{accents}
\newlength{\dhatheight}
\newcommand{\doublehat}[1]{%
    \settoheight{\dhatheight}{\ensuremath{\hat{#1}}}%
    \addtolength{\dhatheight}{-0.35ex}%
    \hat{\vphantom{\rule{1pt}{\dhatheight}}%
    \smash{\hat{#1}}}}

\begin{document}

\title{Near-Optimal High Probability Complexity Bounds for Non-Smooth Stochastic Optimization with Heavy-Tailed Noise}
\author{Eduard Gorbunov$^{1,2}$ \quad Marina Danilova$^{3}$ \quad Innokentiy Shibaev$^{1,2}$ \\
Pavel Dvurechensky$^{4,2,5}$\quad Alexander Gasnikov$^{1,2,4}$}
\date{$^1$ Moscow Institute of Physics and Technology, Russian Federation\\
$^2$ National Research University Higher School of Economics, Russian Federation \\
$^3$ Institute of Control Sciences RAS, Russian Federation\\
$^4$ Weierstrass Institute for Applied Analysis and Stochastics, Germany\\
$^5$ Institute for Information Transmission Problems RAS Russian Federation}
\maketitle

\begin{abstract}
    Thanks to their practical efficiency and random nature of the data, stochastic first-order methods are standard for training large-scale machine learning models. Random behavior may cause a particular run of an algorithm to result in a highly suboptimal objective value, whereas theoretical guarantees are usually proved for the expectation of the objective value. Thus, it is essential to theoretically guarantee that algorithms provide small objective residual with high probability. Existing methods for non-smooth stochastic convex optimization have complexity bounds with the dependence on the confidence level that is either negative-power or logarithmic but under an additional assumption of sub-Gaussian (light-tailed) noise distribution that may not hold in practice, e.g., in several NLP tasks. In our paper, we resolve this issue and derive the first high-probability convergence results with logarithmic dependence on the confidence level for non-smooth convex stochastic optimization problems with non-sub-Gaussian (heavy-tailed) noise. To derive our results, we propose novel stepsize rules for two stochastic methods with gradient clipping. Moreover, our analysis works for generalized smooth objectives with H\"older-continuous gradients, and for both methods, we provide an extension for strongly convex problems. Finally, our results imply that the first (accelerated) method we consider also has optimal iteration and oracle complexity in all the regimes, and the second one is optimal in the non-smooth setting.

\end{abstract}

\section{Introduction}\label{sec:intro}
Stochastic first-order optimization methods like \algname{SGD} \cite{robbins1951stochastic}, \algname{Adam} \cite{kingma2014adam}, and their various modifications are extremely popular in solving a number of different optimization problems, especially those appearing in statistics \cite{spokoiny2012parametric}, machine learning, and deep learning \cite{Goodfellow-et-al-2016}. The success of these methods in real-world applications motivates the researchers to investigate theoretical properties for the methods and to develop new ones with better convergence guarantees. Typically, stochastic methods are analyzed in terms of the convergence in expectation (see \cite{ghadimi2013stochastic, moulines2011non, gower2019sgd} and references therein), whereas high-probability complexity results are established much rarely. However, as illustrated in \cite{gorbunov2020clipped_sstm}, guarantees in terms of the convergence in expectation have much worse correlation with the real behavior of the methods than high-probability convergence guarantees when the noise in the stochastic gradients has \textit{heavy-tailed distribution}.

Recent studies \cite{simsekli2019tail,csimcsekli2019heavy,zhang2020why} show that in several popular problems such as training {\tt BERT} \cite{vaswani2017attention} on {\tt Wikipedia} dataset the noise in the stochastic gradients is heavy-tailed. Moreover, in \cite{zhang2020why}, the authors justify empirically that in such cases \algname{SGD} works significantly worse than \algname{clipped-SGD} \cite{pascanu2013difficulty} and \algname{Adam}. Therefore, it is important to theoretically study the methods' convergence when the noise is heavy-tailed.

For convex and strongly convex problems  {with Lipschitz continuous gradient, i.e., smooth convex and strongly convex problems,} this question was properly addressed in \cite{nazin2019algorithms,davis2021low,gorbunov2020clipped_sstm} where the first high-probability complexity bounds with logarithmic dependence on the confidence level were derived for the stochastic problems with heavy-tailed noise. However, a number of practically important problems are non-smooth \textit{on the whole space} \cite{zhang2020gradient, mai2021stability}. For example, in deep neural network training, the loss function often grows polynomially fast when the norm of the network's weights goes to infinity.  {Moreover, non-smoothness of the activation functions such as ReLU or loss functions such as hinge loss implies the non-smoothness of the whole problem.} While being well-motivated by practical applications, the existing high-probability convergence guarantees for stochastic first-order methods applied to solve non-smooth convex optimization problems with heavy-tailed noise depend on the negative power of the confidence level that dramatically increases the number of iterations required to obtain high accuracy  {of the solution} with probability close to one. Such a discrepancy in the theory between algorithms for stochastic smooth and non-smooth problems leads us to the natural question: \textit{is it possible to obtain high-probability complexity bounds with logarithmic dependence on the confidence level for \textbf{non-smooth} convex stochastic problems with heavy-tailed noise?} In this paper, we give a positive answer to this question.  {To achieve this we focus on gradient clipping methods \cite{pascanu2013difficulty,gehring2017convolutional,menon2020can,mai2021stability,zhang2020gradient,zhang2020why}.}

\subsection{Preliminaries}\label{sec:prelim}
Before we describe our contributions in detail, we formally state the considered setup.
\paragraph{Stochastic optimization.} We focus on the following problem
\begin{equation}
    \min\limits_{x\in \mathbb R^n}f(x),\quad f(x) = \mathbb E_\xi\left[f(x,\xi)\right],\label{eq:main_problem}
\end{equation}
where $f(x)$ is a convex but possibly non-smooth function. Next, we assume that at each point $x\in \mathbb R^n$ we have an access to the unbiased estimator $\nabla f(x,\xi)$ of $\nabla f(x)$ with uniformly bounded variance
\begin{equation}
    \mathbb E_\xi[\nabla f(x,\xi)] = \nabla f(x),\quad
    \mathbb E_\xi\left[\left\|\nabla f(x,\xi) - \nabla f(x)\right\|_2^2\right] \le \sigma^2,\quad \sigma > 0.\label{eq:bounded_variance_clipped_SSTM}
\end{equation}
This assumption on the stochastic {oracle} is widely used in stochastic optimization literature \cite{ghadimi2012optimal, ghadimi2013stochastic, juditsky2011first, lan2012optimal, nemirovski2009robust}. We emphasize that we do not assume that the stochastic gradients have so-called ``light tails'' \cite{lan2012optimal}, i.e., sub-Gaussian noise distribution meaning that $\mathbb P\{\|\nabla f(x,\xi) - \nabla f(x)\|_2 > b\} \le 2\exp(-\nicefrac{b^2}{(2\sigma^2)})$ for all $b > 0$.

\paragraph{Level of smoothness.} Finally, we assume that function $f$ has $(\nu, M_{\nu})$-H\"older continuous gradients on a compact set $Q \subseteq \mathbb R^n$ for some $\nu \in [0,1]$, $M_\nu > 0$ meaning that
\begin{equation}
    \|\nabla f(x) - \nabla f(y)\|_2 \le M_\nu \|x- y\|_2^\nu\quad \forall x,y \in Q.\label{eq:holder_def}
\end{equation}
When $\nu = 1$ inequality \eqref{eq:holder_def} implies $M_1$-smoothness of $f$, and when $\nu = 0$ we have that $\nabla f(x)$
has bounded variation which is equivalent to being uniformly bounded.
Moreover, when $\nu = 0$ differentiability of $f$ is not needed, and one can assume uniform boundedness of the subgradients of $f$. Linear regression in the case when the noise has generalized Gaussian distribution (Example 4.4 from \cite{Chaux_2007}) serves as a natural example of the situation with $\nu \in (0,1)$. Moreover, when \eqref{eq:holder_def} holds for $\nu = 0$ and $\nu = 1$ simultaneously then it holds for all $\nu\in[0,1]$ with $M_\nu \le M_0^{1-\nu}M_1^\nu$ \cite{nesterov2015universal}. As we show in our results, the set $Q$ should contain the ball centered at the solution $x^*$ of \eqref{eq:main_problem} with radius $\sim R_0 \ge \|x^0 - x^*\|_2$, where $x^0$ is a starting point of the method, i.e., our analysis does not require \eqref{eq:holder_def} to hold on $\mathbb R^n$.

\paragraph{High-probability convergence.} For a given accuracy $\varepsilon > 0$ and confidence level $\beta \in (0,1)$ we are interested in finding $\varepsilon$-solutions of problem \eqref{eq:main_problem} with probability at least $1-\beta$, i.e., such $\widehat x$ that $\mathbb P\{f(\widehat x) - f(x^*) \le \varepsilon\} \ge 1 - \beta$. For brevity, we will call such (in general, random) points $\widehat x$ as $(\varepsilon,\beta)$-solution of \eqref{eq:main_problem}. Moreover, by high-probability complexity of a stochastic method ${\cal M}$ we mean the sufficient number of oracle calls{, i.e.,  number of $\nabla f(x,\xi)$ computations,} needed to guarantee that the output of ${\cal M}$ is an $(\varepsilon,\beta)$-solution of \eqref{eq:main_problem}.

\subsection{Contributions}

\begin{table}[t]
    \centering
    \caption{Summary of known and new high-probability complexity bounds for solving \eqref{eq:main_problem} with $f$ being \textbf{convex} and having $(\nu, M_\nu)$-H\"older continuous gradients. Columns: ``Ref.'' = reference, ``Complexity'' = high-probability complexity ($\varepsilon$ -- accuracy, $\beta$ -- confidence level, numerical constants and logarithmic factors are omitted), ``HT'' = heavy-tailed noise, ``UD'' = unbounded domain, ``HCC'' = H\"older continuity of the gradient is required only on the compact set. The results labeled by $^{\clubsuit}$ are obtained from the convergence guarantees in expectation via Markov's inequality. Negative-power dependencies on the confidence level $\beta$ are colored in red.}
    \label{tab:cvx_case_comparison}
    {\small
    \begin{tabular}{|c|c|c|c|c|c|c|}
        \hline
        Method & Ref. & {\centering Complexity} & $\nu$ & HT? & UD? & HCC?\\
        \hline
        \hline
        \algnames{SGD} & \cite{nemirovski2009robust} & $\max\left\{\frac{M_0^2R_0^{2}}{\varepsilon^2},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $0$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{AC-SA} & \cite{ghadimi2012optimal,lan2012optimal} & $\max\left\{\sqrt{\frac{M_1R_0^{2}}{\varepsilon}},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $1$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{SIGMA} & \cite{dvurechensky2016stochastic} & $\max\left\{\frac{M_\nu^{\frac{2}{1+3\nu}}R_0^{\frac{2(1+\nu)}{1+3\nu}}}{\varepsilon^{\frac{2}{1+3\nu}}},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $[0,1]$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{SGD} & \cite{nemirovski2009robust}$^{\clubsuit}$ & $\max\left\{\frac{M_0^2R_0^{2}}{{\color{BrickRed}\beta^2}\varepsilon^2},\frac{\sigma^2R_0^2}{{\color{BrickRed}\beta^2}\varepsilon^2}\right\}$ & $0$ & \cmark & \xmark & \xmark\\
        \hline
        \algnames{AC-SA} & \cite{ghadimi2012optimal,lan2012optimal}$^{\clubsuit}$ & $\max\left\{\sqrt{\frac{M_1R_0^{2}}{{\color{BrickRed}\beta}\varepsilon}},\frac{\sigma^2R_0^2}{{\color{BrickRed}\beta^2}\varepsilon^2}\right\}$ & $1$ & \cmark & \cmark & \xmark\\
        \hline
        \algnames{SIGMA} & \cite{dvurechensky2016stochastic}$^{\clubsuit}$ & $\max\left\{\frac{M_\nu^{\frac{2}{1+3\nu}}R_0^{\frac{2(1+\nu)}{1+3\nu}}}{{\color{BrickRed}\beta^{\frac{2}{1+3\nu}}}\varepsilon^{\frac{2}{1+3\nu}}},\frac{\sigma^2R_0^2}{{\color{BrickRed}\beta^2}\varepsilon^2}\right\}$ & $[0,1]$ & \cmark & \cmark & \xmark\\
        \hline
        \algnames{clipped-SSTM} & \cite{gorbunov2020clipped_sstm} & $\max\left\{\sqrt{\frac{M_1R_0^{2}}{\varepsilon}},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $1$ & \cmark & \cmark & \xmark \\
        \hline
        \algnames{clipped-SGD} & \cite{gorbunov2020clipped_sstm} & $\max\left\{\frac{M_1R_0^2}{\varepsilon},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $1$ & \cmark & \cmark & \xmark \\
        \hline
        \rowcolor{bgcolor2}\algnames{clipped-SSTM} & Thm.~\ref{thm:main_result_clipped_SSTM_str_cvx_main} & $\max\left\{\frac{M_\nu^{\frac{2}{1+3\nu}}R_0^{\frac{2(1+\nu)}{1+3\nu}}}{\varepsilon^{\frac{2}{1+3\nu}}},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $[0,1]$ & \cmark & \cmark & \cmark \\
        \hline
        \rowcolor{bgcolor2}\algnames{clipped-SGD} & Thm.~\ref{thm:main_result_clipped_SGD_main} & $\max\left\{\frac{M_\nu^{\frac{2}{1+\nu}}R_0^2}{\varepsilon^{\frac{2}{1+\nu}}},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}$ & $[0,1]$ & \cmark & \cmark & \cmark \\
        \hline
    \end{tabular}}
\end{table}

\begin{table}[t]
    \centering
    \caption{Summary of known and new high-probability complexity bounds for solving \eqref{eq:main_problem} with $f$ being \textbf{$\mu$-strongly convex} and having $(\nu, M_\nu)$-H\"older continuous gradients. Columns: ``Ref.'' = reference, ``Complexity'' = high-probability complexity ($\varepsilon$ -- accuracy, $\beta$ -- confidence level, numerical constants and logarithmic factors are omitted), ``HT'' = heavy-tailed noise, ``UD'' = unbounded domain, ``HCC'' = H\"older continuity of the gradient is required only on the compact set. The results labeled by $^{\clubsuit}$ are obtained from the convergence guarantees in expectation via Markov's inequality. Negative-power dependencies on the confidence level $\beta$ are colored in red.}
    \label{tab:str_cvx_case_comparison}
    {\small
    \begin{tabular}{|c|c|c|c|c|c|c|}
        \hline
        Method & Ref. & {\centering Complexity} & $\nu$ & HT? & UD? & HCC?\\
        \hline
        \hline
        \algnames{SGD} & \cite{nemirovski2009robust} & $\max\left\{\frac{M_0^2}{\mu\varepsilon},\frac{\sigma^2}{\mu\varepsilon}\right\}$ & $0$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{AC-SA} & \cite{ghadimi2012optimal,lan2012optimal} & $\max\left\{\sqrt{\frac{M_1}{\mu}},\frac{\sigma^2}{\mu\varepsilon}\right\}$ & $1$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{SIGMA} & \cite{dvurechensky2016stochastic} & \makecell{$\max\left\{\hat{N},\frac{\sigma^2}{\mu\varepsilon}\right\}$,\\ $\hat{N} \!=\! \left(\!\frac{M_\nu}{\mu R_0^{1\!-\!\nu}}\!\right)^{\!\frac{2}{1\!+\!3\nu}}\!+\! \left(\!\frac{M_\nu^2}{\mu^{1\!+\!\nu}\varepsilon^{1\!-\!\nu}}\!\right)^{\frac{1}{1\!+\!3\nu}}$} & $[0,1]$ & \xmark & \cmark & \xmark\\
        \hline
        \algnames{SGD} & \cite{nemirovski2009robust}$^{\clubsuit}$ & $\max\left\{\frac{M_0^2}{\mu{\color{BrickRed}\beta}\varepsilon},\frac{\sigma^2}{\mu{\color{BrickRed}\beta}\varepsilon}\right\}$ & $0$ & \cmark & \xmark & \xmark\\
        \hline
        \algnames{AC-SA} & \cite{ghadimi2012optimal,lan2012optimal}$^{\clubsuit}$ & $\max\left\{\sqrt{\frac{M_1}{\mu}},\frac{\sigma^2}{\mu{\color{BrickRed}\beta}\varepsilon}\right\}$ & $1$ & \cmark & \cmark & \xmark\\
        \hline
        \algnames{SIGMA} & \cite{dvurechensky2016stochastic}$^{\clubsuit}$ & \makecell{$\max\left\{\hat{N},\frac{\sigma^2}{\mu\hat{\varepsilon}}\right\}$, $\hat{\varepsilon} = {\color{BrickRed}\beta}\varepsilon$,\\ $\hat{N} \!=\! \left(\!\frac{M_\nu}{\mu R_0^{1\!-\!\nu}}\!\right)^{\!\frac{2}{1\!+\!3\nu}}\!+\! \left(\!\frac{M_\nu^2}{\mu^{1\!+\!\nu}\hat{\varepsilon}^{1\!-\!\nu}}\!\right)^{\frac{1}{1\!+\!3\nu}}$} & $[0,1]$ & \cmark & \cmark & \xmark\\
        \hline
        \algnames{R-clipped-SSTM} & \cite{gorbunov2020clipped_sstm} & $\max\left\{\sqrt{\frac{M_1}{\mu}},\frac{\sigma^2}{\mu\varepsilon^2}\right\}$ & $1$ & \cmark & \cmark & \xmark \\
        \hline
        \algnames{R-clipped-SGD} & \cite{gorbunov2020clipped_sstm} & $\max\left\{\frac{M_1}{\mu},\frac{\sigma^2}{\mu\varepsilon^2}\right\}$ & $1$ & \cmark & \cmark & \xmark \\
        \hline
        \rowcolor{bgcolor2}\algnames{R-clipped-SSTM} & Thm.~\ref{thm:main_result_clipped_SSTM_main} & \makecell{$\max\left\{\hat{N},\frac{\sigma^2}{\mu\varepsilon}\right\}$,\\ $\hat{N} =\! \left(\!\frac{M_\nu}{\mu R_0^{1\!-\!\nu}}\!\right)^{\!\frac{2}{1\!+\!3\nu}}\!+\! \left(\!\frac{M_\nu^2}{\mu^{1\!+\!\nu}\varepsilon^{1\!-\!\nu}}\!\right)^{\frac{1}{1\!+\!3\nu}}$} & $[0,1]$ & \cmark & \cmark & \cmark \\
        \hline
        \rowcolor{bgcolor2}\algnames{R-clipped-SGD} & Thm.~\ref{thm:main_result_clipped_SGD_str_cvx_main} & $\max\left\{\frac{M_\nu^{\frac{2}{1+\nu}}}{\mu^{\frac{2}{1+\nu}}R_0^{\frac{2(1-\nu)}{1+\nu}}}, \frac{M_\nu^{\frac{2}{1+\nu}}}{\mu \varepsilon^{\frac{1-\nu}{1+\nu}}}, \frac{\sigma^2}{\mu\varepsilon}\right\}$ & $[0,1]$ & \cmark & \cmark & \cmark \\
        \hline
    \end{tabular}}
\end{table}

\begin{itemize}[leftmargin=*]
    \item We propose novel stepsize rules for \algname{clipped-SSTM} \cite{gorbunov2020clipped_sstm} to  {handle} the problems with H\"older continuous gradients and derive high-probability complexity guarantees for convex stochastic optimization problems without using ``light tails'' assumption, i.e., we prove that  our version of \algname{clipped-SSTM}
    \begin{equation*}
        {\cal O}\left(\max\left\{D\ln^{\frac{2(1+\nu)}{1+3\nu}}\frac{D}{\beta},\frac{\sigma^2R_0^2}{\varepsilon^2}\ln\frac{D}{\beta}\right\}\right), \quad D = \frac{M_\nu^{\frac{2}{1+3\nu}}R_0^{\frac{2(1+\nu)}{1+3\nu}}}{\varepsilon^{\frac{2}{1+3\nu}}}
    \end{equation*}
    high-probability complexity. Unlike all previous high-probability complexity results in this setup with $\nu < 1$ (see Tbl.~\ref{tab:cvx_case_comparison}), our result depends only logarithmically on the confidence level $\beta$ that is highly important when $\beta$ is small. Moreover, up to the difference in logarithmic factors the derived complexity guarantees meet the known lower bounds \cite{lan2012optimal, guzman2015lower} obtained for the problems with light-tailed noise.  {In particular, when $\nu = 1$ we recover accelerated convergence rate \cite{nesterov1983method,lan2012optimal}}. That is, neglecting the logarithmic factors our results are unimprovable  {and,  {surprisingly} coincide with the best-known results in the ``light-tailed case''.}
    \item We derive the first high-probability complexity bounds for \algname{clipped-SGD} when the objective functions is convex with $(\nu, M_\nu)$-H\"older continuous gradient and the noise is heavy tailed., i.e., we derive
    \begin{equation*}
        {\cal O}\left(\max\left\{D^2,\max\left\{D^{1+\nu},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}\ln\frac{D^2 + D^{1+\nu}}{\beta}\right\}\right),\quad D = \frac{M_\nu^{\frac{1}{1+\nu}}R_0}{\varepsilon^{\frac{1}{1+\nu}}}
    \end{equation*}
    high-probability complexity bound. Interestingly, when $\nu = 0$ the derived bound for \algname{clipped-SGD} has better dependence on the logarithms than the corresponding one for \algname{clipped-SSTM}. Moreover, neglecting the dependence on $\varepsilon$ under the logarithm, our bound for \algname{clipped-SGD} has the same dependence on the confidence level as the tightest known result in this case under the ``light tails'' assumption \cite{guigues2017non}.
    \item Using restarts technique we extend the obtained results for \algname{clipped-SSTM} and \algname{clipped-SGD} to the strongly convex case (see Tbl.~\ref{tab:str_cvx_case_comparison}). As in the convex case, the obtained results are superior to all previous known results in the general setup we consider.
    \item As one of the key contributions of this work, we emphasize that in our theoretical results it is sufficient to assume H\"older continuity of the gradients of $f$ only on the ball with radius $2R_0 = 2\|x^0 - x^*\|_2$ and centered at a solution of the problem. This makes our results applicable to much larger class of problems than functions with H\"older continuous gradients on $\mathbb R^n$, e.g., our analysis works even for polynomially growing objectives.
    \item To test the performance of the considered methods we conduct several numerical experiments on image classification and NLP tasks, and observe that 1) \algname{clipped-SSTM} and \algname{clipped-SGD} show a comparable performance with \algname{SGD} on the image classification task, when the noise distribution is almost sub-Gaussian, 2) converge much faster than \algname{SGD} on the NLP task, when the noise distribution is heavy-tailed, and 3) \algname{clipped-SSTM} achieves a comparable performance with \algname{Adam} on the NLP task enjoying both the best known theoretical guarantees and good practical performance.
\end{itemize}

\subsection{Related work}
\paragraph{Light-tailed noise.} The theory of high-probability complexity bounds for convex stochastic optimization with light-tailed noise is well-developed. Lower bounds and optimal methods for the problems with $(\nu,M_\nu)$-H\"older continuous gradients are obtained in \cite{nemirovski2009robust} for $\nu = 0$, and in \cite{ghadimi2012optimal} for $\nu = 1$. Up to the logarithmic dependencies these high-probability convergence bounds coincide with the corresponding results for the convergence in expectation (see first two rows of Tbl.~\ref{tab:cvx_case_comparison})
While not being directly derived in the literature, the lower bound for the case when $\nu \in (0,1)$ can be obtained as a combination of lower bounds in the deterministic \cite{nemirovsky1983problem, guzman2015lower} and smooth stochastic settings \cite{ghadimi2012optimal}. The corresponding optimal methods are analyzed in \cite{devolder2013exactness,dvurechensky2016stochastic} through the lens of inexact oracle.

\paragraph{Heavy-tailed noise.} Unlike in the ``light-tailed'' case, the first theoretical guarantees with reasonable dependence on both the accuracy $\varepsilon$ and the confidence level $\beta$ appeared just recently. In \cite{nazin2019algorithms}, the first such results without acceleration \cite{nesterov1983method} were derived for Mirror Descent with special truncation technique for smooth ($\nu = 1$) convex problems on a bounded domain, and then were accelerated and extended in \cite{gorbunov2020clipped_sstm}. For the strongly convex problems the first accelerated high-probability convergence guarantees were obtained in \cite{davis2021low} for the special method called \algname{proxBoost} requiring solving auxiliary nontrivial problems at each iteration. These bounds were tightened in \cite{gorbunov2020clipped_sstm}.

In contrast, for the case when $\nu < 1$ and, in particular, when $\nu = 0$ the best-known high-probability complexity bounds suffer from the negative-power dependence on the confidence level $\beta$, i.e., have a factor $\nicefrac{1}{\beta^{\alpha}}$ for some $\alpha >0$, that affects the convergence rate dramatically for small enough $\beta$. Without additional assumptions on the tails these results are obtained via Markov's inequality $\mathbb P\{f(x) - f(x^*) > \varepsilon\} < \nicefrac{\mathbb E[f(x) - f(x^*)]}{\varepsilon}$ from the guarantees for the convergence in expectation to the accuracy $\varepsilon\beta$, see the results labeled by $^\clubsuit$ in Tbl.~\ref{tab:cvx_case_comparison}. Under an additional assumption on noise tails that $\mathbb P\{\|\nabla f(x,\xi) - \nabla f(x)\|_2^2 > s\sigma^2\} = O(s^{-\alpha})$ for $\alpha > 2$ these results can be tightened \cite{gasnikov2015efficiency} when $\nu = 0$ as
$\sim \max\left\{ \nicefrac{\ln\left(\beta^{-1}\right)}{\varepsilon^2},\left(\nicefrac{1}{\beta\varepsilon^{\alpha}}\right)^{\nicefrac{2}{(3\alpha - 2)}}\right\}$
without removing the negative-power dependence on the confidence level $\beta$. Different stepsize policies allow to change the last term in $\max$ to $\beta^{-\frac{1}{2\alpha - 1}}\varepsilon^{-\frac{2\alpha}{2\alpha - 1}}$ without removing the negative-power dependence on $\beta$.

\paragraph{Gradient clipping.} The methods based on gradient clipping \cite{pascanu2013difficulty} and normalization \cite{hazan2015beyond} are popular in different machine learning and deep learning tasks due to their robustness in practice to the noise in the stochastic gradients and rapid changes of the objective function \cite{Goodfellow-et-al-2016}. In \cite{zhang2020gradient,mai2021stability}, \algname{clipped-GD} and \algname{clipped-SGD} are theoretically studied in applications to non-smooth problems that can grow polynomially fast when $\|x - x^*\|_2 \to \infty$ showing the superiority of gradient clipping methods to the methods without clipping. The results from \cite{zhang2020gradient} are obtained for non-convex problems with almost surely bounded noise, and in \cite{mai2021stability}, the authors derive the stability and expectation convergence guarantees for strongly convex  {under assumption that the central $p$-th moment of the stochastic gradient is bounded for $p \ge 2$. Since the authors of \cite{mai2021stability} do not provide convergence guarantees with explicit dependencies on all important parameters of the problem it complicates direct comparison with our results. Nevertheless, convergence guarantees from \cite{mai2021stability} are sub-linear and are given for the convergence in expectation, and, as a consequence, the corresponding high-probability convergence results obtained via Markov's inequality also suffer from negative-power dependence on the confidence level.} Next, in \cite{zhang2020why}, the authors establish several expectation convergence guarantees for \algname{clipped-SGD} and prove their optimality in the non-convex case under assumption that the central $\alpha$-moment of the stochastic gradient is uniformly bounded, where $\alpha \in (1,2]$. It turns out that \algname{clipped-SGD} is able to converge even when $\alpha < 2$, whereas vanilla \algname{SGD} can diverge in this setting.

\section{Clipped Stochastic Similar Triangles Method}\label{sec:clipped_SSTM_main}
In this section, we propose a novel variation of Clipped Stochastic Similar Triangles Method \citep{gorbunov2020clipped_sstm} adjusted to the class of objectives with H\"older continuous gradients (\algname{clipped-SSTM}, see Alg.~\ref{alg:clipped-SSTM}).
\begin{algorithm}[h]
\caption{Clipped Stochastic Similar Triangles Method (\algname{clipped-SSTM}): case $\nu \in [0,1]$}
\label{alg:clipped-SSTM}
\begin{algorithmic}[1]
\Require starting point $x^0$, number of iterations $N$, batchsizes $\{m_k\}_{k=1}^N $, stepsize parameter $\alpha$, clipping parameter $B$, H\"older exponent $\nu \in [0,1]$.
\State Set $A_0 = \alpha_0 = 0$, $y^0 = z^0 = x^0$
\For{$k=0,\ldots, N-1$}
\State Set $\alpha_{k+1} = \alpha (k+1)^{\frac{2\nu}{1+\nu}}$, $A_{k+1} = A_k + \alpha_{k+1}$, $\lambda_{k+1} = \frac{B}{\alpha_{k+1}}$
\State $x^{k+1} = \nicefrac{(A_k y^k + \alpha_{k+1} z^k)}{A_{k+1}}$
\State Draw mini-batch $m_k$ of fresh i.i.d.\ samples $\xi_1^k,\ldots,\xi_{m_{k}}^k$ and compute $\nabla f(x^{k+1},\boldsymbol{\xi}^k) = \frac{1}{m_k}\sum_{i=1}^{m_k}\nabla f(x^{k+1},\xi_i^k)$
\State Compute $\widetilde{\nabla} f(x^{k+1},\boldsymbol{\xi}^k) = \text{clip}(\nabla f(x^{k+1},\boldsymbol{\xi}^k),\lambda_{k+1})$ using \eqref{eq:clip_operator}
\State $z^{k+1} = z^k - \alpha_{k+1}\widetilde{\nabla} f(x^{k+1},\boldsymbol{\xi}^k)$
\State $y^{k+1} = \nicefrac{(A_k y^k + \alpha_{k+1} z^{k+1})}{A_{k+1}}$
\EndFor
\Ensure $y^N$
\end{algorithmic}
\end{algorithm}

The method is based on the clipping of the stochastic gradients:
\begin{equation}
    \text{clip}(\nabla f(x,\boldsymbol{\xi}), \lambda) = \min\left\{1,\frac{\lambda}{\|\nabla f(x,\boldsymbol{\xi})\|_2}\right\}\nabla f(x,\boldsymbol{\xi}) \label{eq:clip_operator}
\end{equation}
where $\nabla f(x,\boldsymbol{\xi})=\frac{1}{m}\sum_{i=1}^m\nabla f(x,\xi_i)$ is a mini-batched stochastic gradient. Gradient clipping ensures that the resulting vector has a norm bounded by the clipping level $\lambda$. Since the clipped stochastic gradient cannot have arbitrary large norm, the clipping helps to avoid unstable behavior of the method when the noise is heavy-tailed and the clipping level $\lambda$ is properly adjusted.

However, unlike the stochastic gradient, clipped stochastic gradient is a \textit{biased} estimate of $\nabla f(x)$: the smaller the clipping level the larger the bias. The biasedness of the clipped stochastic gradient complicates the analysis of the method. On the other hand, to circumvent the negative effect of the heavy-tailed noise on the high-probability convergence one should choose $\lambda$ to be not too large. Therefore, the question on the appropriate choice of the clipping level is highly non-trivial.

Fortunately, there exists a simple but insightful observation that helps us to obtain the right formula for the clipping level $\lambda_{k}$ in \algname{clipped-SSTM}: if $\lambda_k$ is chosen in such a way that $\|\nabla f(x^{k})\|_2 \le \nicefrac{\lambda_k}{2}$ with high probability, then for the realizations $\nabla f(x^{k+1},\boldsymbol{\xi}^{k})$ of the mini-batched stochastic gradient such that $\|\nabla f(x^{k+1},\boldsymbol{\xi}^{k}) - \nabla f(x^{k+1})\|_2 \le \nicefrac{\lambda_k}{2}$ the clipping is an identity operator. Next, if the probability mass of such realizations is big enough then the bias of the clipped stochastic gradient is properly bounded that helps to derive needed convergence guarantees. It turns out that the choice $\lambda_k \sim \nicefrac{1}{\alpha_k}$ ensures the method convergence with needed rate and high enough probability.

Guided by this observation we derive the precise expressions for all the parameters of \algname{clipped-SSTM} and derive high-probability complexity bounds for the method. Below we provide a simplified version of the main result for \algname{clipped-SSTM} in the convex case. The complete formulation and the full proof of the theorem are deferred to Appendix~\ref{sec:clipped_SSTM_cvx_appendix} (see Thm.~\ref{thm:main_result_clipped_SSTM}).
\begin{theorem}\label{thm:main_result_clipped_SSTM_main}
    Assume that function $f$ is convex and its gradient satisfy \eqref{eq:holder_def} with $\nu \in [0,1]$, $M_\nu > 0$ on $Q = B_{3R_0} = \{x\in\mathbb R^n\mid \|x-x^*\|_2 \le 2R_0\}$, where $R_0 \ge \|x^0 - x^*\|_2$. Then there exist such a choice of parameters that \algname{clipped-SSTM} achieves $f(y^N) - f(x^*) \le \varepsilon$ with probability at least $1-\beta$ after ${\cal O}\left(D\ln^{\frac{2(1+\nu)}{1+3\nu}}\frac{D}{\beta}\right)$ iterations with $D = \frac{M_\nu^{\frac{2}{1+3\nu}}R_0^{\frac{2(1+\nu)}{1+3\nu}}}{\varepsilon^{\frac{2}{1+3\nu}}}$ and requires
    \begin{equation}
        {\cal O}\left(\max\left\{D\ln^{\frac{2(1+\nu)}{1+3\nu}}\frac{D}{\beta},\frac{\sigma^2R_0^2}{\varepsilon^2}\ln\frac{D}{\beta}\right\}\right)\text{ oracle calls.} \label{eq:clipped_SSTM_oracle_complexity_main}
    \end{equation}
\end{theorem}

The obtained result has only logarithmic dependence on the confidence level $\beta$ and optimal dependence on the accuracy $\varepsilon$ up to logarithmic factors \cite{lan2012optimal, guzman2015lower} for all $\nu \in [0,1]$. Moreover, we emphasize that our result does not require $f$ to have $(\nu,M_\nu)$-H\"older continuous gradient on the whole space. This is because we prove that for the proposed choice of parameters the iterates of \algname{clipped-SSTM} stay inside the ball $B_{3R_0} = \{x\in\mathbb R^n\mid \|x-x^*\|_2 \le 3R_0\}$ with probability at least $1 - \beta$, and, as a consequence, H\"older continuity of the gradient is required only inside this ball. In particular, this means that the better starting point leads not only to the reduction of $R_0$, but also it can reduce $M_\nu$. Moreover, our result is applicable to much wider class of functions than the standard class of functions with H\"older continuous gradients in $\mathbb R^n$, e.g., to the problems with polynomial growth.

For the strongly convex problems, we consider restarted version of Alg.~\ref{alg:clipped-SSTM} (\algname{R-clipped-SSTM}, see Alg.~\ref{alg:R-clipped-SSTM}) and derive high-probability complexity result for this version.
\begin{algorithm}[h]
\caption{Restarted \algname{clipped-SSTM} (\algname{R-clipped-SSTM}): case $\nu \in [0,1]$}
\label{alg:R-clipped-SSTM}
\begin{algorithmic}[1]
\Require starting point $x^0$, number of restarts $\tau$, number of steps of \algname{clipped-SSTM} in restarts $\{N_t\}_{t=1}^{\tau}$, batchsizes $\{m_{k}^1\}_{k=1}^{N_1-1}, \{m_{k}^2\}_{k=1}^{N_2-1}, \ldots, \{m_{k}^{\tau}\}_{k=1}^{N_{\tau}-1}$, stepsize parameters $\{\alpha^t\}_{t=1}^\tau$, clipping parameters $\{B_t\}_{t=1}^\tau$, H\"older exponent $\nu \in [0,1]$.
\State $\hat{x}^0 = x^0$
\For{$t=1,\ldots, \tau$}
\State Run \algname{clipped-SSTM} (Alg.~\ref{alg:clipped-SSTM}) for $N_t$ iterations with batchsizes $\{m_k^t\}_{k=1}^{N_t-1}$, stepsize parameter $\alpha_{t}$, clipping parameter $B_t$, and starting point $\hat{x}^{t-1}$. Define the output of \algname{clipped-SSTM} by $\hat{x}^{t}$.
\EndFor
\Ensure $\hat{x}^\tau$
\end{algorithmic}
\end{algorithm}
Below we provide a simplified version of the result. The complete formulation and the full proof of the theorem are deferred to Appendix~\ref{sec:clipped_SSTM_str_cvx_appendix} (see Thm.~\ref{thm:main_result_clipped_SSTM_str_cvx}).
\begin{theorem}\label{thm:main_result_clipped_SSTM_str_cvx_main}
    Assume that function $f$ is $\mu$-strongly convex and its gradient satisfy \eqref{eq:holder_def} with $\nu \in [0,1]$, $M_\nu > 0$ on $Q = B_{3R_0} = \{x\in\mathbb R^n\mid \|x-x^*\|_2 \le 3R_0\}$, where $R_0 \ge \|x^0 - x^*\|_2$. Then there exist such a choice of parameters that \algname{R-clipped-SSTM} achieves $f(\hat{x}^\tau) - f(x^*) \le \varepsilon$ with probability at least $1-\beta$ after
    \begin{equation}
        \hat{N} = O\left(D\ln^{\frac{2(1+\nu)}{1+3\nu}}\frac{D}{\beta}\right),\quad D = \max\left\{\left(\frac{M_\nu}{\mu R_0^{1-\nu}}\right)^{\frac{2}{1+3\nu}}\ln\frac{\mu R_0^2}{\varepsilon}, \left(\frac{M_\nu^2}{\mu^{1+\nu}\varepsilon^{1-\nu}}\right)^{\frac{1}{1+3\nu}}\right\} \label{eq:clipped_SSTM_iter_complexity_str_cvx_main}
    \end{equation}
    iterations of Alg.~\ref{alg:clipped-SSTM} in total and requires
    \begin{equation}
        O\left(\max\left\{D\ln^{\frac{2(1+\nu)}{1+3\nu}}\frac{D}{\beta},\frac{\sigma^2}{\mu\varepsilon}\ln\frac{D}{\beta}\right\}\right)\text{ oracle calls.} \label{eq:clipped_SSTM_oracle_complexity_str_cvx_main}
    \end{equation}
\end{theorem}
Again, the obtained result has only logarithmic dependence on the confidence level $\beta$ and, as our result in the convex case, it has optimal dependence on the accuracy $\varepsilon$ up to logarithmic factors depending on $\beta$ \cite{lan2012optimal, guzman2015lower} for all $\nu \in [0,1]$.

\section{SGD with clipping}\label{sec:clipped_SGD_main}

In this section, we present a new variant of \algname{clipped-SGD} \cite{pascanu2013difficulty} properly adjusted to the class of objectives with $(\nu, M_\nu)$-H\"older continuous gradients (see Alg.~\ref{alg:clipped-SGD}).

\begin{algorithm}[h]
\caption{Clipped Stochastic Gradient Descent (\algname{clipped-SGD}): case $\nu \in [0,1]$}
\label{alg:clipped-SGD}
\begin{algorithmic}[1]
\Require starting point $x^0$, number of iterations $N$, batchsize $m$, stepsize $\gamma$, clipping parameter $B > 0$.
\For{$k=0,\ldots, N-1$}
\State Draw mini-batch of $m$ fresh i.i.d.\ samples $\xi_1^k,\ldots,\xi_{m}^k$ and compute $\nabla f(x^{k+1},\boldsymbol{\xi}^k) = \frac{1}{m}\sum_{i=1}^{m}\nabla f(x^{k+1},\xi_i^k)$
\State Compute $\widetilde{\nabla} f(x^{k},\boldsymbol{\xi}^k) = \text{clip}(\nabla f(x^{k},\boldsymbol{\xi}^k),\lambda)$ using \eqref{eq:clip_operator} with $\lambda = \nicefrac{B}{\gamma}$
\State $x^{k+1} = x^k - \gamma\widetilde{\nabla} f(x^{k},\boldsymbol{\xi}^k)$
\EndFor
\Ensure $\Bar{x}^N = \frac{1}{N}\sum_{k=0}^{N-1}x^k$
\end{algorithmic}
\end{algorithm}

We emphasize that as for \algname{clipped-SSTM} we use clipping level $\lambda$ inversely proportional to the stepsize $\gamma$. Below we provide a simplified version of the main result for \algname{clipped-SGD} in the convex case. The complete formulation and the full proof of the theorem are deferred to Appendix~\ref{sec:clipped_SGD_cvx_appendix} (see Thm.~\ref{thm:main_result_clipped_SGD}).
\begin{theorem}\label{thm:main_result_clipped_SGD_main}
    Assume that function $f$ is convex and its gradient satisfy \eqref{eq:holder_def} with $\nu \in [0,1]$, $M_\nu > 0$ on $Q = B_{7R_0} = \{x\in\mathbb R^n\mid \|x-x^*\|_2 \le 7R_0\}$, where $R_0 \ge \|x^0 - x^*\|_2$. Then there exist such a choice of parameters that \algname{clipped-SGD} achieves $f(\Bar{x}^N) - f(x^*) \le \varepsilon$ with probability at least $1-\beta$ after
    \begin{equation}
        {\cal O}\left(\max\left\{D^2,D^{1+\nu}\ln\frac{D^2 + D^{1+\nu}}{\beta}\right\}\right),\quad D = \frac{M_\nu^{\frac{1}{1+\nu}}R_0}{\varepsilon^{\frac{1}{1+\nu}}} \label{eq:clipped_SGD_iter_complexity_main}
    \end{equation}
    iterations and requires
    \begin{equation}
       {\cal O}\left(\max\left\{D^2,\max\left\{D^{1+\nu},\frac{\sigma^2R_0^2}{\varepsilon^2}\right\}\ln\frac{D^2 + D^{1+\nu}}{\beta}\right\}\right) \text{ oracle calls.} \label{eq:clipped_SGD_oracle_complexity_main}
    \end{equation}
\end{theorem}
As all our results in the paper, this result for \algname{clipped-SGD} has two important features: 1) the dependence on the confidence level $\beta$ is logarithmic and 2) H\"older continuity is required only on the ball $B_{2R_0}$ centered at the solution. Moreover, up to the difference in the expressions under the logarithm the dependence on $\varepsilon$ in the result for \algname{clipped-SGD} is the same as in the tightest known results for non-accelerated \algname{SGD}-type methods \cite{devolder2013exactness, guigues2017non}. Finally, we emphasize that for $\nu < 1$ the logarithmic factors appearing in the complexity bound for \algname{clipped-SSTM} are worse than the corresponding factor in the complexity bound for \algname{clipped-SGD}. Therefore, \algname{clipped-SGD} has the best known high-probability complexity results in the case when $\nu = 0$ and $f$ is convex.

For the strongly convex problems, we consider restarted version of Alg.~\ref{alg:clipped-SGD} (\algname{R-clipped-SGD}, see Alg.~\ref{alg:R-clipped-SGD}) and derive high-probability complexity result for this version.
\begin{algorithm}[h]
\caption{Restarted \algname{clipped-SGD} (\algname{R-clipped-SGD}): case $\nu \in [0,1]$}
\label{alg:R-clipped-SGD}
\begin{algorithmic}[1]
\Require starting point $x^0$, number of restarts $\tau$, number of steps of \algname{clipped-SGD} in restarts $\{N_t\}_{t=1}^{\tau}$, batchsizes $\{m_{t}\}_{k=1}^{\tau}$, stepsizes $\{\gamma_t\}_{t=1}^\tau$, clipping parameters $\{B_t\}_{t=1}^\tau$
\State $\hat{x}^0 = x^0$
\For{$t=1,\ldots, \tau$}
\State Run \algname{clipped-SGD} (Alg.~\ref{alg:clipped-SGD}) for $N_t$ iterations with batchsize $m_t$, stepsize $\gamma_{t}$, clipping parameter $B_t$, and starting point $\hat{x}^{t-1}$. Define the output of \algname{clipped-SGD} by $\hat{x}^{t}$.
\EndFor
\Ensure $\hat{x}^\tau$
\end{algorithmic}
\end{algorithm}
Below we provide a simplified version of the result. The complete formulation and the full proof of the theorem are deferred to Appendix~\ref{sec:clipped_SGD_str_cvx_appendix} (see Thm.~\ref{thm:main_result_clipped_SGD_str_cvx}).
\begin{theorem}\label{thm:main_result_clipped_SGD_str_cvx_main}
    Assume that function $f$ is $\mu$-strongly convex and its gradient satisfy \eqref{eq:holder_def} with $\nu \in [0,1]$, $M_\nu > 0$ on $Q = B_{7R_0} = \{x\in\mathbb R^n\mid \|x-x^*\|_2 \le 7R_0\}$, where $R_0 \ge \|x^0 - x^*\|_2$. Then there exist such a choice of parameters that \algname{R-clipped-SGD} achieves $f(\Bar{x}^N) - f(x^*) \le \varepsilon$ with probability at least $1-\beta$ after
    \begin{equation}
        {\cal O}\left(\max\left\{D_1^{\frac{2}{1+\nu}}\ln\frac{\mu R_0^2}{\varepsilon}, D_2^{\frac{2}{1+\nu}}, \max\left\{D_1\ln\frac{\mu R_0^2}{\varepsilon},D_2\right\}\ln\frac{D}{\beta} \right\}\right) \notag
    \end{equation}
    iterations of Alg.~\ref{alg:clipped-SGD} in total and requires
    \begin{equation*}
        {\cal O}\left(\max\left\{D_1^{\frac{2}{1+\nu}}\ln\frac{\mu R_0^2}{\varepsilon}, D_2^{\frac{2}{1+\nu}}, \max\left\{D_1\ln\frac{\mu R_0^2}{\varepsilon},D_2, \frac{\sigma^2}{\mu\varepsilon}\right\}\ln\frac{D}{\beta} \right\}\right) \; \text{oracle calls, where}
    \end{equation*}

    \begin{equation*}
        D_1 = \frac{M_\nu}{\mu R_0^{1-\nu}},\quad D_2 = \frac{M_\nu}{\mu^{\frac{1+\nu}{2}}\varepsilon^{\frac{1-\nu}{2}}},\quad D = (D_1^{\frac{2}{1+\nu}}+D_1)\ln\frac{\mu R_0^2}{\varepsilon} + D_2 + D_2^{\frac{2}{1+\nu}}.
    \end{equation*}
\end{theorem}
As in the convex case, for $\nu < 1$ the log  factors appearing in the complexity bound for \algname{R-clipped-SSTM} are worse than the corresponding factor in the bound for \algname{R-clipped-SGD}. Thus, \algname{R-clipped-SGD} has the best known high-probability complexity results for strongly convex $f$ and $\nu = 0$.

\section{Numerical experiments}\label{sec:experiments}
We tested the performance of the methods on the following problems:
\begin{itemize}[leftmargin=*]
\setlength{\itemsep}{0em}
    \item {\tt BERT} fine-tuning on {\tt CoLA} dataset \cite{warstadt2018neural}. We use pretrained {\tt BERT} from Transformers library \cite{wolf-etal-2020-transformers} ({\tt bert-base-uncased}) and freeze all layers except the last two linear ones.
    \item {\tt ResNet-18} training on {\tt ImageNet-100} (first 100 classes of {\tt ImageNet} \cite{ILSVRC15}).
\end{itemize}

First, we study the noise distribution for both problem as follows: at the starting point we sample large enough number of batched stochastic gradients $\nabla f(x^0,\boldsymbol{\xi}_1),\ldots, \nabla f(x^0,\boldsymbol{\xi}_K)$ with batchsize $32$ and plot the histograms for $\|\nabla f(x^0,\boldsymbol{\xi}_1) - \nabla f(x^0)\|_2,\ldots, \|\nabla f(x^0,\boldsymbol{\xi}_K) - \nabla f(x^0)\|_2$, see Fig.~\ref{figure:norm_diffs_distribution}. As one can see, the noise distribution for {\tt BERT} + {\tt CoLA} is substantially non-sub-Gaussian, whereas the distribution for {\tt ResNet-18} + {\tt Imagenet-100} is almost Gaussian.

\begin{figure}
  \centering
  \includegraphics[scale=0.35]{figures/ResNet-18_norm_diffs_start.pdf}
  \includegraphics[scale=0.35]{figures/BERT_norm_diffs_start.pdf}
  \caption{Noise distribution of the stochastic gradients for {\tt ResNet-18} on {\tt ImageNet-100} and {\tt BERT} fine-tuning on the {\tt CoLA} dataset before the training. Red lines: probability density functions with means and variances empirically estimated by the samples. Batch count is the total number of samples used to build a histogram.}
  \label{figure:norm_diffs_distribution}
\end{figure}

Next, we compared $4$ different optimizers on these problems: \algname{Adam}, \algname{SGD} (with Momentum), \algname{clipped-SGD} (with Momentum and coordinate-wise clipping) and \algname{clipped-SSTM} (with norm-clipping and $\nu=1$). The results are presented in Fig.~\ref{figure:losses_graph_bert_resnet}. We observed that the noise distributions do not change significantly along the trajectories of the considered methods, see Appendix~\ref{sec:extra_experiments}. During the hyper-parameters search we compared different batchsizes, emulated via gradient accumulation (thus we compare methods with different batchsizes by the number of base batches used). The base batchsize was $32$ for both problems, stepsizes and clipping levels were tuned. One can find additional details regarding our experiments in Appendix~\ref{sec:extra_experiments}.

\begin{figure}[h]
  \centering
  \includegraphics[scale=0.35]{figures/ResNet-18_train_loss.pdf}
  \includegraphics[scale=0.35]{figures/ResNet-18_val_loss.pdf}
  \includegraphics[scale=0.35]{figures/ResNet-18_val_acc.pdf}
  \includegraphics[scale=0.35]{figures/BERT_train_loss.pdf}
  \includegraphics[scale=0.35]{figures/BERT_val_loss.pdf}
  \includegraphics[scale=0.35]{figures/BERT_val_acc.pdf}
  \caption{Train and validation loss + accuracy for different optimizers on both problems. Here, ``batch count'' denotes the total number of used stochastic gradients.}
  \label{figure:losses_graph_bert_resnet}
\end{figure}

\paragraph{Image classification.} On {\tt ResNet-18} + {\tt ImageNet-100} task, \algname{SGD} performs relatively well, and even ties with \algname{Adam} (with batchsize of $4\times 32$) in validation loss. \algname{clipped-SSTM} (with batchsize of $2\times 32$) also ties with \algname{Adam} and \algname{clipped-SGD} is not far from them. The results were averaged from 5 different launches (with different starting points/weight initializations). Since the noise distribution is almost Gaussian even vanilla \algname{SGD} performs well, i.e., gradient clipping is not required. At the same time, the clipping does not slow down the convergence significantly.
\vspace{-0.2cm}
\paragraph{Text classification.} On {\tt BERT} + {\tt CoLA} task, when the noise distribution is heavy-tailed, the methods with clipping outperform \algname{SGD} by a large margin. This result is in good correspondence with the derived high-probability complexity bounds for \algname{clipped-SGD}, \algname{clipped-SSTM} and the best-known ones for \algname{SGD}. Moreover, \algname{clipped-SSTM} (with batchsize of $8\times 32$) achieves the same loss on validation as \algname{Adam}, and has better accuracy. These results were averaged from $5$ different train-val splits and $20$ launches (with different starting points/weight initializations) for each of the splits, $100$ launches in total.

\section*{Acknowledgments}
The work of
Marina Danilova in Sections~\ref{sec:clipped_SGD_main} and \ref{sec:clipped_SGD_appendix} was funded by RFBR, project number 20-31-90073.