\section{Introduction}
We consider the stochastic optimization problem:
\begin{equation}
\label{MainProb}
\min_{x\in\mathbb{R}^d} \left[ f(x) = \frac{1}{n} \sum_{i=1}^n f_i(x) \right],
\end{equation}
where each $f_i$ is convex and lower bounded. We denote by $\mathcal{X}^*$ the non-empty set of optimal points $x^*$ of equation~\eqref{MainProb}. We set $f^*:=\min_{x\in\mathbb{R}^d} f(x)$, and $f_i^* := \inf_{x\in\mathbb{R}^d} f_i(x)$.

In this setting, the algorithm of choice is often Stochastic Gradient Descent~(SGD), i.e. $x^{k+1} =  x^k - \gamma_k \nabla f_{\SC_k}(x^k)$, where $\gamma_k>0$ is the stepsize at iteration $k$, $\SC_k\subseteq[n]$ a random subset of datapoints~(minibatch) with cardinality $B$ sampled independently at each iteration $k$, and $\nabla f_{\SC_k}(x):=\frac{1}{B}\sum_{i \in \SC_k} \nabla f_i(x)$ is the minibatch gradient. \looseness=-1

\vspace{-1mm}
A careful choice of $\gamma_k$ is crucial for most applications~\citep{bottou2018optimization,goodfellow2016deep}. The simplest option is to pick $\gamma_k$ to be constant over training, with its value inversely proportional to the  Lipschitz constant of the gradient. While this choice yields fast convergence to the neighborhood of a minimizer, two main problems arise: (a) the optimal $\gamma$ depends on (often unknown) problem parameters --- hence often requires heavy tuning ; and (b) it cannot be guaranteed that $\mathcal{X}^*$ is reached in the limit~\citep{ghadimi2013stochastic,gower2019sgd, gower2021sgd}. A simple fix for the last problem is to allow polynomially decreasing stepsizes~(second option)~\citep{nemirovski2009robust}: this choice for $\gamma_k$ often leads to convergence to $\mathcal{X}^*$, but hurts the overall algorithm speed. The third option, which became very popular with the rise of deep learning, is to implement an \textit{adaptive} stepsize. These methods do not commit to a fixed schedule, but instead use the optimization statistics~(e.g. gradient history, cost history) to tune the value of $\gamma_k$ at each iteration. These stepsizes are known to work very well in deep learning~\citep{zhang2019adaptive}, and include Adam~\citep{kingma2014adam}, Adagrad~\citep{duchi2011adaptive}, and RMSprop~\citep{tieleman2012lecture}.

Ideally, a theoretically grounded adaptive method should yield fast convergence to $\mathcal{X}^*$ without knowledge of problem dependent parameters, such as the gradient Lipshitz constant or the strong convexity constant. As a result, an ideal adaptive method should require very little tuning by the user, while matching the performance of a fine-tuned $\gamma_k$. However, while in practice this is the case for common adaptive methods such as Adam and AdaGrad, the associated convergence rates often rely on strong assumptions --- e.g. that the iterates live on a bounded domain, or that gradients are uniformly bounded in norm~\citep{duchi2011adaptive,ward2019adagrad, vaswani2020adaptive}. While the above assumptions are valid in the constrained setting, they are problematic for problems defined in the whole~$\R^d$.

A promising new direction in the adaptive stepsizes literature is based on the idea of Polyak stepsizes, introduced by~\citet{polyak1987introduction} in the context of deterministic convex optimization. Recently~\citet{loizou2021stochastic} successfully adapted Polyak stepsizes to the stochastic setting, and provided convergence rates matching fine-tuned SGD --- while the algorithm does not require knowledge of the unknown quantities such as the gradient Lipschitz constant. The results especially shines in the overparameterized strongly convex setting, where linear convergence to $x^*$ is shown. This result is especially important since, under the same assumption, no such rate exists for AdaGrad~(see e.g.~\citep{vaswani2020adaptive} for the latest results) or other adaptive stepsizes. Moreover, the method was shown to work surprisingly well on deep learning problems, without requiring heavy tuning~\citep{loizou2021stochastic}.

Even if the stochastic Polyak stepsize~(SPS)~\citet{loizou2021stochastic} comes with strong convergence guarantees, it has two main drawbacks when it is used in non-over-parameterized regimes: (i)  It requires \textit{a priori} knowledge of the optimal mini-batch losses, which are not often available for big batch sizes or regularized objectives~(see discussion in~\S\ref{sec:background}) and (ii) it guarantees convergence only to a neighborhood of the solution.
In this work, we study the dynamics and the convergence properties of SGD equipped with new variants of SPS for solving general convex optimization problems. Our new proposed variants provide solutions to both drawbacks of the original SPS.

 \subsection{Background and Technical Preliminaries}
\label{sec:background}
The stepsize proposed by~\citet{loizou2021stochastic} is
\vspace{-2mm}
\begin{equation}
\label{SPSmax}
\tag{SPS$_{\max}$}
\gamma_k = \min \left\{ \frac{f_{\SC_k}(x^k)-f_{\SC_k}^*}{c\|\nabla f_{\SC_k}(x^k)\|^2}, \gamma_b \right\},
\end{equation}
where $\gamma_b, c>0$ are problem-independent constants, $f_{\SC_k} :=\frac{1}{|\SC_k|}\sum_{i\in\SC_k} f_i$, $f_{\SC_k}^* = \inf_{x\in\mathbb{R}^d} f_{\SC_k}(x)$.

\paragraph{Dependency on $\boldsymbol{f_{\SC_k}}$.} Crucially the algorithm requires knowledge of $f_{\SC_k}^*$ for every realization of the mini-batch $\SC_k$. In the non-regularized overparametrized setting~(e.g. neural networks), $f_{\SC_k}$ is often zero for every subset $\SC$~\citep{zhang2021understanding}. However, this is not the only setting where $f_{\SC}^*$ is computable: e.g., in the regularized logistic loss with batch size 1, it is possible to recover a cheap closed form expression for each $f_{i}^*$~\citep{loizou2021stochastic}. Unfortunately, if the batch-size is bigger than $1$ or the loss becomes more demanding~(e.g. cross-entropy), then \textit{no such closed-form computation is possible}.

\vspace{-3mm}
\paragraph{Rates and comparison with AdaGrad.} In the convex overparametrized setting~(more precisely, under the interpolation condition, i.e. $\exists x^* \in \mathcal{X}^*$ :
$\inf_{x\in\mathbb{R}^d} f_\SC(x) = f_\SC(x^*)$ for all $\SC$, see also~\S\ref{sec:SPS}), SPS$_{\max}$ enjoys a convergence speed of $\mathcal{O}(1/k)$, without requiring knowledge of the gradient Lipschitz constant or other problem parameters. Recently,~\citep{vaswani2020adaptive} showed that the same rate can be achieved for AdaGrad in the same setting. However, there is an important difference: the rate of~\citep{vaswani2020adaptive} is technically $\mathcal{O}(dD^2/k)$, where $d$ is the problem dimension and $D^2$ is a global bound on the squared distance to the minimizer, which is assumed to be finite. Not only does SPS$_{\max}$ not have this dimension dependency, which dates back to crucial arguments in the AdaGrad literature~\citep{duchi2011adaptive,levy2018online}, but also does not require bounded iterates. While this assumption is satisfied in the constrained setting, it has no reason to hold in the unconstrained scenario. Unfortunately, this is a common problem of all AdaGrad variants: with the exception of~\citep{xie2020linear}~(which works in a slightly different scenario), no rate can be provided in the stochastic setting without the bounded iterates/gradients~\citep{ene2020adaptive} assumption --- even after assuming strong convexity. However, in the non-interpolated setting, AdaGrad enjoys a convergence guarantee of $\mathcal{O}(1/\sqrt{k})$~(with the bounded iterates assumption). A similar rate does not yet exist for SPS --- our work aims at filling this gap.

\subsection{Main Contributions}

As we already mentioned, in the non-interpolated setting SPS$_{\max}$ has the following issues:
\begin{enumerate}[leftmargin=52pt]
    \item[\textbf{Issue (1)}:] For $B>1$ (minibatch setting), SPS$_{\max}$ requires the exact knowledge of $f_{\SC}^*$. This is not practical.
    \item[\textbf{Issue (2)}:] SPS$_{\max}$ guarantees convergence to a neighborhood of the solution. It is not clear how to modify it to yield convergence to the exact minimizer.
\end{enumerate}
Having the above two issues in mind, the main contributions of our work (see also Table~\ref{tab:results}, for a summary of the main complexity results obtained in this paper) are summarized as follows:
\vspace{-2mm}
\begin{itemize}[leftmargin=*]
    \item In \S\ref{sec:estimation}, we provide a direct solution for Issue (1). We explain how only a lower bound on $f_{\SC}^*$~(trivial if all $f_i$s are non-negative) is required for convergence to a neighborhood of the solution. While this neighborhood is bigger that the one for SPS$_{\max}$, our modified version provides a practical baseline for the solution to the second issue.
    \item We explain why Issue (2) is highly non-trivial and requires an in-depth study of the bias induced by the interaction between gradients and Polyak stepsizes. Namely, we show that simply multiplying the stepsize of SPS$_{\max}$ by $1/\sqrt{k}$ --- which would work for vanilla SGD~\citep{nemirovski2009robust} --- yields a bias in the solution found by SPS~(\S\ref{sec:bias}), regardless of the estimation of $f_{\SC}^*$.
    \item In \S\ref{sec:convergence_decr}, we provide a solution to the problem (Issue (2)) by introducing additional structure --- as well as the fix to Issue (1) ---  into the stepsize. We call the new algorithm \emph{Decreasing SPS}~(DecSPS), and provide a convergence guarantee under the bounded domain assumption --- matching the standard AdaGrad results.
    \item In \S\ref{sec:no_bound} we go one step further and show that, if strong convexity is assumed, iterates are bounded with probability 1 and hence we can remove the bounded iterates assumption. To the best of our knowledge, DecSPS, is the first stochastic adaptive optimization method that converges to the exact solution without assuming strong assumptions like bounded iterates/gradients.
    \item In \S\ref{sec:non-smooth} we provide extensions of our approach to the non-smooth setting.
    \item In \S\ref{sec:exp}, we corroborate our theoretical results with experimental testing.
\end{itemize}
\vspace{-1mm}

\begin{table}[ht]
    \centering
    \scalebox{0.78}{
    \begin{tabular}{l|c|l|c|c|l}
    \toprule
        \textbf{Stepsize} &  \textbf{Citation}& \textbf{Assumptions} & \begin{tabular}{cc} \textbf{No Knowledge} \\ \textbf{of} $\boldsymbol{f_{\SC}^*}$ \end{tabular} & \begin{tabular}{cc} \textbf{Exact} \\ \textbf{Convergence} \end{tabular} &  \textbf{Theorem}\\
        \midrule
        SPS$_{\max}$ &\citep{loizou2021stochastic} &convex, smooth &\xmark& \xmark&  Thm. \ref{thm:loizou}\\
        \midrule
        SPS$_{\max}^\ell$ & This paper & convex, smooth & \cmark&\xmark  &  Thm. \ref{thm:loizou_bs}, Cor. \ref{cor:loizou_bs_cor}\\
        \midrule
        DecSPS & This paper &  convex, smooth, bounded iterates&\cmark &\cmark & Cor.~\ref{cor:sqrt_smooth}, $\mathcal{O}(1/\sqrt{K})$ \\
         \cline{2-6}
         & This paper &  strongly-convex, smooth &\cmark & \cmark & Thm.~\ref{thm:no_bound}, $\mathcal{O}(1/\sqrt{K})$ \\
        \midrule
        DecSPS-NS & This paper &  convex, bounded iterates/grads&\cmark &\cmark &  Cor.~\ref{cor:sqrt_NS}, $\mathcal{O}(1/\sqrt{K})$ \\
        \bottomrule
    \end{tabular}}
    \caption{\small Summary of the considered stepsizes and the corresponding theoretical results in the non-interpolated setting. The studied quantity in all Theorems, with respect to which all rates are expressed is $\Exp \left[f(\bar{x}^K)-f(x^*)\right]$, where $\bar{x}^K=\frac{1}{K}\sum_{k=0}^{K-1} x^k$. In addition, for all converging methods, we consider the stepsize scaling factor $c_k = \mathcal{O}(\sqrt{k})$, formally defined in the corresponding sections. For the methods without exact convergence, we show in \S\ref{sec:bias} that any different scaling factor cannot make the algorithm convergent.}
    \label{tab:results}
\end{table}

\section{Background on Stochastic Polyak Stepsize}
\label{sec:SPS}
In this section, we provide a concise overview of the results in~\citet{loizou2021stochastic}, and highlight the main assumptions and open questions.

To start, we remind the reader that problem~\eqref{MainProb} is said to be interpolated if there exists a problem solution $x^*\in\mathcal{X}^*$ such that $\inf_{x\in\mathbb{R}^d} f_i(x) = f_i(x^*)$ for all $i\in[n]$. The degree of interpolation at batch size $B$ can be quantified by the following quantity, introduced by~\citet{loizou2021stochastic} and studied also in~\citet{vaswani2020adaptive,d2021stochastic}: fix a batch size $B$, and let $\SC\subseteq [n]$ with $|\SC| = B$.
\begin{equation}
    \sigma^2_B := \Exp_{\SC}[ f_\SC(x^*)-f_\SC^*] = f(x^*)-\Exp_\SC [ f_\SC^*]
\end{equation}
It is easy to realize that as soon as problem~\eqref{MainProb} is interpolated, then $\sigma^2_B=0$ for each $B\le n$. In addition, note that $\sigma^2_B$ is non-increasing as a function of $B$.

We now comment on the main result from~\citet{loizou2021stochastic}.

\begin{restatable}[Main results from~\citet{loizou2021stochastic}]{thm}{loizou}
\label{thm:loizou}
Let each $f_i$ be $L_i$-smooth convex functions. Then SGD with SPS$_{\max}$, mini-batch size $B$, and $c=1$, converges as:  \begin{align*}
\Exp \left[f(\bar{x}^K)-f(x^*)\right]
\leq \frac{\|x^0-x^*\|^2}{\alpha \, K} + \frac{2\gamma_{b}\sigma^2_B}{\alpha},
\end{align*}
where $\alpha=\min \left\{\frac{1}{2cL_{\max}},\gamma_{b}\right\}$ and $\bar{x}^K=\frac{1}{K}\sum_{k=0}^{K-1} x^k$. If in addition $f$ is $\mu$-strongly convex, then, for any $c\geq1/2$, SGD with SPS$_{\max}$ converges as:
\vspace{-2mm}
\begin{align*}
\Exp \|x^{k}-x^*\|^2
\leq \left(1-\mu \alpha \right)^k  \|x^0-x^*\|^2 + \frac{2\gamma_{b} \sigma^2_B }{\mu \alpha},
\end{align*}
where again $\alpha=\min \{\frac{1}{2cL_{\max}},\gamma_{b}\}$ and $L_{\max}=\max \{L_i\}_{i=1}^n$ is the maximum smoothness constant.
\end{restatable}

In the overparametrized setting, the result guarantees convergence to the exact minimizer, without knowledge of the gradient Lipschitz constant~(as vanilla SGD would instead require) and without assuming bounded iterates~(in contrast to~\citet{vaswani2020adaptive}).

As soon as (1) a \textit{regularizer} is applied to the loss~(e.g. $L2$ penalty), or (2) the number of datapoints gets comparable to the dimension, then the problem is not interpolated and SPS$_{\max}$ only converges to a neighborhood and it gets impractical to compute $f_\SC^*$ --- \textit{this is the setting we study in this paper.}

\begin{remark}[What if $\|\nabla f_{\SC_k}\|=0$?]
In the rare case that $\|\nabla f_{S_k}(x^k)\|^2 =0$, there is no need to evaluate the stepsize. In this scenario, the update direction $\nabla f_{S_k}(x^k)=0$ and thus the iterate is not updated irrespective of the choice of step-size. If this happens, the user should simply sample a different minibatch. We note that in our experiments~(see \S\ref{sec:exp}), such event never occurred.
\end{remark}
\vspace{-4mm}
\paragraph{Related work on Polyak stepsize:} The classical Polyak stepsize \citep{polyak1987introduction} has been successfully used in the analysis of deterministic subgradient methods in different settings \citep{boyd2003subgradient, davis2018subgradient, hazan2019revisiting}. First attempts on providing an efficient variant of the stepsize that works well in the stochastic setting were made in \citet{berrada2019training, oberman2019stochastic}. However, as explained in \citet{loizou2021stochastic}, none of these approaches provide a natural stochastic extension with strong theoretical convergence guarantees, and thus \citet{loizou2021stochastic} proposed the stochastic Polyak stepsize SPS$_{\max}$ as a better alternative.\footnote{A variant of SGD with SPS$_{\max}$ was also proposed by~\citet{asi2019stochastic} as a special case of a model-based method called the lower-truncated model. \citet{asi2019stochastic} also proposed a decreasing step-size variant of SPS$_{\max}$ which is closely related but different than the DecSPS that we propose in \S\ref{sec:convergence_decr}. Among some differences, they assume interpolation for their convergence results whereas we do not in~\S\ref{sec:convergence_decr}. We describe the differences between our work and~\citet{asi2019stochastic} in more detail in Appendix~\ref{sec:app_related_work}.} Despite its recent appearance, SPS$_{\max}$ has already been used and analyzed as a stepsize for SGD for solving structured non-convex problems \citep{gower2021sgd}, in combination with other adaptive methods \citep{vaswani2020adaptive}, with a moving target~\citep{gower2021stochastic} and in the update rule of stochastic mirror descent~\citep{d2021stochastic}. These extensions are orthogonal to our approach, and we speculate that our proposed variants can also be used in the above settings. We leave such extensions for future work.

\section[]{Removing $\boldsymbol{f_\SC^*}$ from SPS}
\label{sec:estimation}

\begin{figure*}
    \centering
    \textbf{\scriptsize No Interpolation and $\boldsymbol{f^*>0}$  \quad\quad\quad\quad\quad Interpolation and $\boldsymbol{f^*>0}$\quad\quad\quad\quad\quad\quad Interpolation and $\boldsymbol{f^*=0}$}
    \includegraphics[height = 0.21\textwidth]{fig/ell_i/SPS_ell_i_1.pdf}
    \includegraphics[height = 0.21\textwidth]{fig/ell_i/SPS_ell_i_2.pdf}
    \includegraphics[height = 0.21\textwidth]{fig/ell_i/SPS_ell_i_3.pdf}
    \vspace{-3mm}
    \caption{ \small We consider a $100$ dim problem with $n=100$ datapoints where each $f_i = \frac{1}{2}(x-x_i^*)^\top H_i(x-x_i^*) + f_i^*$, with $f_i^*=1$ for all $i\in[n]$ and $H_i$ a random SPD matrix generated using the standard Gaussian matrix $A_i\in\mathbb{R}^{d\times 3d}$ as $H_i = A_iA_i^\top/(3d)$. If $x_i^*\ne x_j^*$ for $i\ne j$, then the problem does \textbf{not satisfy interpolation (left plot)}. Instead, if all $x_i^*$s are equal, \textbf{then the problem is interpolated (central plot)}. The plot shows the behaviour of SPS$^{\ell}_{\max}$~($\gamma_b=2)$ for different choices of the approximated suboptimality $\ell_i^*$. We plot~(mean and std deviation over 10 runs) the function suboptimality level $f(x)-f(x^*)$ for different values of $\ell_i^*$. Note that, if instead $f_i^*=0$ for all $i$ then all the shown \textbf{algorithms coincide~(right plot)} and converge to the solution.}
    \label{fig:ell_i_plot}
    \vspace{-1.5mm}
\end{figure*}

As motivated in the last sections, computing $f_\SC^*$ in the non-interpolated setting is not practical. In this section, we explore the effect of using a lower bound $\ell_\SC^*\le f_{\SC}^*$ instead in the SPS$_{\max}$ definition.
\begin{tcolorbox}
\vspace{-1mm}
\begin{equation}
\tag{SPS$_{\max}^\ell$}
\gamma_k = \min \left\{ \frac{f_{\SC_k}(x^k)-\ell_{\SC_k}^*}{c\|\nabla f_{\SC_k}(x^k)\|^2}, \gamma_b \right\},
\end{equation}
\end{tcolorbox}
Such a lower bound is easy to get for many problems of interest: indeed, for standard regularized regression and classification tasks, the loss is non-negative hence one can pick $\ell_{\SC}^* = 0$, for any $\SC\subseteq[n]$.

The obvious question is: what is the effect of estimating $\ell_{\SC}^*$ on the convergence rates in Thm.~\ref{thm:loizou}? We found that the proof of~\citep{loizou2021stochastic} is easy to adapt to this case, by using the following fundamental bound~(see also Lemma~\ref{lemma:basic_loizou}):
\vspace{-2mm}
\begin{align}
    \frac{1}{2c L_{\SC_k}}\le \frac{f_{\SC_k}(x^k)-f_{\SC_k}^*}{c\|\nabla f_{\SC_k}(x^k)\|^2} \le\frac{f_{\SC_k}(x^k)-\ell_{\SC_k}^*}{c\|\nabla f_{\SC_k}(x^k)\|^2}.
\end{align}
The following results can be seen as an easy extension of the main result of~\citep{loizou2021stochastic}, under a newly defined suboptimality measure:
\begin{equation}
\label{sigmahat}
    \hat\sigma_B^2 := \Exp_{\SC_k}[f_{\SC_k}(x^*)-\ell_{\SC_k}^*] = f(x^*)-\Exp_{\SC_k}[\ell_{\SC_k}^*].
\end{equation}

\begin{restatable}[]{thm}{loizoubs}
\label{thm:loizou_bs}
Under SPS$_{\max}^\ell$, the same exact rates in Thm.~\ref{thm:loizou} hold~(under the corresponding assumptions), after replacing $\sigma_B^2$ with $\hat\sigma_B^2$.
\end{restatable}
\vspace{-2mm}
And we also have an easy practical corollary.
\begin{restatable}[]{cor}{loizoubscor}
\label{cor:loizou_bs_cor}
In the context of Thm.~\ref{thm:loizou_bs}, assume all $f_i$s are non-negative and estimate $\ell_{\SC}^* = 0$ for all $\SC\subseteq[n]$. Then the same exact rates in Thm.~\ref{thm:loizou} hold for SPS$_{\max}^\ell$, after replacing $\sigma^2_B$ with $f^*=f(x^*)$.
\end{restatable}
\vspace{-2mm}
A numerical illustration of this result can be found in Fig.~\ref{fig:ell_i_plot}. In essence, both theory and experiments confirm that, if interpolation is not satisfied, then we have a rate until a convergence ball, where the size is optimal under exact knowledge of $f^*_{\SC}$. Instead, under interpolation, if all the $f_i$s are non-negative and $f^*=0$, then SPS$_{\max}$=SPS$^{\ell}_{\max}$. Finally, in the less common case in practice where $f^*>0$ but we still have interpolation, then SPS$_{\max}$ converges to the exact solution while SPS$^{\ell}_{\max}$ does not. To conclude SPS$^{\ell}_{\max}$ does not (of course) work better than SPS$_{\max}$, but it is a practical variant which we can use as a baseline in \S\ref{sec:convergence_decr} for an adaptive stochastic Polyak stepsize with convergence to the true $x^*$ in the non-interpolated setting.

\section{Bias in the SPS dynamics.}
\label{sec:bias}

In this section, we study convergence of the standard SPS$_{\max}$ in the non-interpolated regime, under an additional~(decreasing) multiplicative factor, in the most ideal setting: batch size 1, and we have knowledge of each $f_i^*$. That is, we consider $\gamma_k =\min\{\frac{f_{i_k}(x^k)-f_{i_k}^*}{c_k\|\nabla f_{i_k}(x^k)\|^2},\gamma_b\}$ with $c_k\to \infty$, e.g. $c_k=\mathcal{O}(\sqrt{k})$ or $c_k=\mathcal{O}(k)$. We note that, in the SGD case, simply picking e.g. $\gamma_k = \gamma_0/\sqrt{k+1}$ would guarantee convergence of $f(x^k)$ to $f(x^*)$, in expectation and with high probability~\citep{kushner2003stochastic,nemirovski2009robust}. Therefore, it is natural to expect a similar behavior for SPS, if $1/c_k$ safisfies the usual Robbins-Monro conditions~\citep{robbins1951stochastic}: $\sum_{k=0}^\infty 1/c_k =\infty, \sum_{k=0}^\infty 1/c_k^2 <\infty$.

\textit{We show that this is not the case}: quite interestingly, $f(x^k)$ converges to a biased solution due to the \textit{correlation between $\nabla f_{i_k}$ and $\gamma_k$}. we show this formally, in the case of non-interpolation~(otherwise both SGD and SPS do not require a decreasing learning rate).
\paragraph{Counterexample.}
Consider the following finite-sum setting: $f(x) = \frac{1}{2} f_1(x) + \frac{1}{2} f_2(x)$ with $f_1(x)=\frac{a_1}{2}(x-1)^2, f_2(x)=\frac{a_2}{2}(x+1)^2$.
To make the problem interesting, we choose $a_1=2$ and $a_2=1$: this introduces asymmetry in the average landscape with respect to the origin. During optimization, we sample $f_1$ and $f_2$ independently and seek convergence to the unique minimizer $x^*=\frac{a_1-a_2}{a_1+a_2}=1/3$. The first thing we notice is that $x^*$ is not a stationary point for the dynamics under SPS.
Indeed note that since $f_i^*=0$ for $i=1,2$ we have~(assuming $\gamma_b$ large enough): $\gamma_k \nabla f_{i_k}(x) =  \frac{x-1}{2c_k}$, if $i_k = 1$, and $\gamma_k \nabla f_{i_k}(x) = \frac{x+1}{2c_k}$ if $i_k = 2$.

Crucially, note that this update is \textit{curvature-independent}. The expected update is $\Exp_{i_k}[\gamma_k\nabla f_{i_k}(x)] =\frac{x-1}{4c_k}+\frac{x+1}{4c_k}=\frac{1}{2c_k} x$. Hence, the iterates can only converge to $x=0$ --- because this is the only fixed point for the update rule. The proof naturally extends to the multidimensional setting, an illustration can be found in Fig.~\ref{fig:convergence_wrong_SPS}.

\begin{figure}
\centering
    \textbf{\small\ \     No interpolation\ \ \ \ \ \ \ Interpolation\ \ \ \ \ }\\
    \includegraphics[width=0.5\linewidth]{fig/land.png}

    \caption{\small Dynamics of SPS$_{\max}$ with decreasing multiplicative constant~(SGD style) compared with DecSPS. We compared both in the \textbf{interpolated setting~(right)} and in the \textbf{non-interpolated setting~(left)}. In the non-interpolated setting, a simple multiplicative factor introduces a bias in the final solution, as discussed in this section. We consider two dimensional $f_i=\frac{1}{2}(x-x_i^*)^\top H_i (x-x_i^*)$, for $i = 1,2$ and plot the contour lines of the corresponding landscapes, as well as the average landscape $(f_1+f_2)/2$ we seek to minimize. Solution is denoted with a gold star.}

    \label{fig:convergence_wrong_SPS}
\end{figure}

In the same picture, we show how our modified variant of the vanilla stepsize --- we call this new algorithm DecSPS, see \S\ref{sec:convergence_decr} --- instead converges to the correct solution.

\begin{remark}
SGD with~(non-adaptive) stepsize $\gamma_k$ instead keeps the curvature, and therefore is able to correctly estimate the average $\Exp_{i_k}[\gamma_k \nabla f_{i_k}(x)] = \frac{\gamma_k}{2}(a_1+a_2)\left[x-\frac{a_1-a_2}{a_1+a_2}\right]$ --- precisely because $\gamma_k$ is independent from $\nabla f_{i_k}$. From this we can see that SGD can only converge to the correct stationary point $x^* = \frac{a_1-a_2}{a_1+a_2}$ --- because again this is the only fixed point for the update rule.
\end{remark}
\vspace{-2mm}
In the appendix, we go one step further and provide an analysis of the bias of SPS in the one-dimensional quadratic case~(Prop.~\ref{prop:bias}). Yet, we expect the precise characterization of the bias phenomenon in the non-quadratic setting to be particularly challenging. We provide additional insights in \S\ref{app:bias_corr}. Instead, in the next section, we show how to effectively modify $\gamma_k$ to yield convergence to $x^*$ without further assumptions.
\vspace{-1mm}
\section{DecSPS: Convergence to the exact solution}
\label{sec:convergence_decr}
\vspace{-3mm}
We propose the following modification of the vanilla SPS proposed in~\citet{loizou2021stochastic}, designed to yield convergence to the exact minimizer while keeping the main adaptiveness properties\footnote{Similar choices are possible. We found that this leads to the biggest stepsize magnitude, allowing for faster convergence in practice.}. We call it Decreasing SPS~(DecSPS), since it combines a steady stepsize decrease with the adaptiveness of SPS.
\begin{tcolorbox}
\begin{equation}
\tag{DecSPS}
{\small
 \gamma_k:=\frac{1}{c_k}\min\left\{\frac{f_{\SC_k}(x^k)-\ell^*_{\SC_k}}{\|\nabla f_{\SC_k}(x^k)\|^2}, \ c_{k-1}\gamma_{k-1}\right\},}
    \label{eq:decr_step}
\end{equation}
\end{tcolorbox}
for $k\in\mathbb{N}$, where $c_k\ne 0$ for every $k\in\mathbb{N}$. We set $c_{-1}=c_0$ and $\gamma_{-1} = \gamma_b>0$~(stepsize bound, similar to~\citet{loizou2021stochastic}), to get $\gamma_0:=\frac{1}{c_0}\cdot\min\left\{\frac{f_{\SC_k}(x^k)-\ell^*_{\SC_k}}{\|\nabla f_{\SC_k}(x^k)\|^2},\quad c_0\gamma_{b}\right\}$

\begin{restatable}[]{lem}{sandwichbounds}
\label{lemma:easy_bounds}
Let each $f_i$ be $L_i$ smooth and let $(c_k)_{k=0}^\infty$ be any non-decreasing positive sequence of real numbers.  Under DecSPS, we have $\min\left\{\frac{1}{2c_k L_{\max}},\frac{c_0\gamma_b}{c_k}\right\}\le\gamma_k\le \frac{c_0\gamma_b}{c_k}$, and $\gamma_{k-1}\le\gamma_k$
\end{restatable}
\begin{remark}
As stated in the last lemma, under the assumption of $c_k$ non-decreasing, $\gamma_k$ is trivially \textit{non-increasing} since $\gamma_k\le c_{k-1}\gamma_{k-1}/c_k$.
\end{remark}
The proof can be found in the appendix, and is based on a simple induction argument.

\vspace{-2mm}

\subsection{Convergence under bounded iterates}
\vspace{-2mm}
The following result provides a proof of convergence of SGD for the $\gamma_k$ sequence defined above.
\begin{restatable}[]{thm}{decSPScvx}
\label{thm:SPS_bounded_domain_cvx} Consider SGD with DecSPS and let $(c_k)_{k=0}^\infty$ be any non-decreasing sequence such that $c_k\ge1, \forall k\in\mathbb{N}$. Assume that each $f_i$ is convex and $L_i$ smooth. We have:
\begin{equation}
 \Exp[f(\bar{x}^K)-f(x^*)] \leq \frac{2c_{K-1} \tilde L D^2}{K} +  \frac{1}{K}\sum_{k=0}^{K-1} \frac{\hat\sigma^2_{B}}{c_k},
\end{equation}
\vspace{-2mm}
where $D^2 :=\max_{k \in [K-1]}\|x^k-x^*\|^2$, $\tilde L := \max\left \{ \max_i\{L_i\},\frac{1}{2c_0\gamma_b}\right\}$ and $\bar{x}^K=\frac{1}{K}\sum_{k=0}^{K-1} x^k$.
\end{restatable}
If $\hat \sigma^2_B=0$, then $c_k=1$ for all $k\in\mathbb{N}$ leads to a rate $\mathcal{O}(\frac{1}{K})$, well known from~\citet{loizou2021stochastic}. If $\hat\sigma^2_B>0$, as for the standard SGD analysis under decreasing stepsizes, the choice $c_k=\mathcal{O}(\sqrt{k})$ leads to an optimal asymptotic trade-off between the deterministic and the stochastic terms, hence to the asymptotic rate $\mathcal{O}(1/\sqrt{k})$ since $\sum_{k=0}^{K-1}\frac{1}{\sqrt{k+1}}\le 2\sqrt{K}$. Moreover, picking $c_0=1$ minimizes the speed of convergence for the deterministic factor. Under the assumption that $\hat \sigma^2_B\ll \tilde L D^2$~(e.g. reasonable distance initialization-solution and $L_{\max} > 1/\gamma_b$), this factor is dominant compared to the factor involving $\hat \sigma^2_B$. For this setting, the rate simplifies as follows.
\begin{restatable}[]{cor}{CorSqrtSmooth}
\label{cor:sqrt_smooth}
Under the setting of Thm.~\ref{thm:SPS_bounded_domain_cvx}, for $c_k=\sqrt{k+1}$~($c_{-1}=c_0$) we have
\begin{equation}
 \Exp[f(\bar{x}^K)-f(x^*)] \leq \frac{2\tilde L D^2 + 2\hat\sigma^2_{B}}{\sqrt{K}}.
\end{equation}
\end{restatable}

\begin{remark}[Beyond bounded iterates]
The result above crucially relies on the bounded iterates assumption: $D^2<\infty$. To the best of our knowledge, if no further regularity is assumed, modern convergence results for adaptive methods~(e.g. variants of AdaGrad) in convex stochastic programming require\footnote{Perhaps the only exception is the result of~\citep{ward2019adagrad}, where the authors work on a different setting: i.e. they introduce the RUIG inequality.} this assumption, or else require gradients to be globally bounded. To mention a few: \citet{duchi2011adaptive, reddi2019convergence, ward2019adagrad, defossez2020simple,vaswani2020adaptive}. A simple algorithmic fix to this problem is adding a cheap projection step onto a large bounded domain~\citep{levy2018online}. We can of course include this projection step in DecSPS, and the theorem above will hold with no further modification. Yet we found this to be not necessary: the strong guarantees of SPS in the strongly convex setting~\citep{loizou2021stochastic} let us go one step beyond: in \S\ref{sec:no_bound} we show that, if each $f_i$ is strongly convex~(e.g. regularizer is added), then one can bound the iterates globally with probability one, without knowledge of the gradient Lipschitz constant. To the best of our knowledge, no such result exist for AdaGrad --- except~\citet{traore2021sequential}, for the deterministic case.
\end{remark}

\begin{remark}[Dependency on the problem dimension]
In standard results for AdaGrad, a dependency on the problem dimension often appears~(e.g. Thm.~1 in~\citet{vaswani2020adaptive}). This dependency follows from a bound on the AdaGrad preconditioner that can be found e.g. in Thm.~4 in~\citet{levy2018online}. In the SPS case no such dependency appears --- specifically because the stepsize is lower bounded by $1/(2c_kL_{\max})$.
\end{remark}

\subsection[]{Removing the bounded iterates assumption}
\label{sec:no_bound}
We prove that under DecSPS the iterates live in a set of diameter $D_{\max}$ almost surely. This can be done by assuming strong convexity of each $f_i$.

The result uses this alternative definition of \textit{neighborhood}: $\hat\sigma^2_{B,\max} := \max_{\SC\subseteq[n], |\SC| =B}[ f_{\SC}(x^*)-\ell_{\SC}^*]$.

Note that trivially $\hat\sigma^2_{B,\max}<\infty$ under the assumption that all $f_{i}$ are lower bounded and $n<\infty$.

\begin{restatable}{prop}{PropSCUnbound}
\label{prop:strong_convex_bound_2}
Let each $f_i$ be $\mu_i$-strongly convex and $L_i$-smooth. The iterates of SGD with DecSPS with $c_k =\sqrt{k+1}$~(and $c_{-1}=c_0$) are such that $\|x^k-x^*\|^2\le D^2_{\max}$ almost surely $\forall k\in\mathbb{N}$, where
\begin{equation}
    D^2_{\max}:=\max\left\{\|x^{0}-x^*\|^2, \frac{2c_0\gamma_b\hat\sigma^2_{B,\max}}{\min\left\{\frac{\mu_{\min}}{2 L_{\max}},\mu_{\min}\gamma_b\right\}}\right\},
\end{equation}
, with $\mu_{\min} = \min_{i\in[n]}\mu_i$ and $L_{\max} = \max_{i\in[n]}L_i$.
\end{restatable}
\vspace{-1mm}
The proof relies on the variations of constants formula and an induction argument --- it is provided in the appendix. We are now ready to state the main theorem for the unconstrained setting, which follows  from Prop.~\ref{prop:strong_convex_bound_2} and Thm.~\ref{thm:SPS_bounded_domain_cvx}.
\begin{restatable}{thm}{ThmNoBound}
\label{thm:no_bound}
Consider SGD with the DecSPS stepsize $\gamma_k:=\frac{1}{\sqrt{k+1}}\cdot\min\left\{\frac{f_{\SC_k}(x^k)-f^*_{\SC_k}}{\|\nabla f_{\SC_k}(x^k)\|^2},\gamma_{k-1}\sqrt{k}\right\}$, for $k\ge 1$ and $\gamma_0$ defined as at the beginning of this section. Let each $f_i$ be $\mu_i$-strongly convex and $L_i$-smooth:
\vspace{-2mm}
\begin{equation}
 \Exp[f(\bar{x}^K)-f(x^*)] \leq \frac{2\tilde L D^2_{\max}+2\hat\sigma^2_{B}}{\sqrt{K}}.
\end{equation}
\end{restatable}
\begin{remark}
The careful reader might notice that, while we assumed strong convexity, our rate is slower than the optimal $\mathcal{O}(1/K)$. This issue  is due to the adaptive nature of DecSPS. It is indeed notoriously hard to achieve a convergence rate of $\mathcal{O}(1/K)$ for adaptive methods in the strongly convex regime. While further investigations will for sure shed light on this interesting problem, we note that \textit{the result we provide is somewhat unique in the literature}: we are not aware of any adaptive method that enjoys such a convergence rate without either (a) assuming bounded iterates/gradients or (b) assuming knowledge of the gradient Lipschitz constant and or the strong convexity constant.
\end{remark}
\subsection{Extension to the non-smooth setting}
\label{sec:non-smooth}

For any $\SC\subseteq[n]$, we denote in this section by $g_{\SC}(x)$ the subgradient of $f_{\SC}$ evaluated at $x$. We discuss the extension of DecSPS to the non-smooth setting.

A straightforward application of DecSPS leads to a stepsize $\gamma_k$ which is no longer lower bounded~(see Lemma~\ref{lemma:easy_bounds}) by the positive quantity $\min\left\{\frac{1}{2c_k L_{\max}},\frac{c_0\gamma_b}{c_k}\right\}$. Indeed, the gradient Lipschitz constant in the non-smooth case is formally $L_{\max} = \infty$. Hence, $\gamma_k$ prescribed by DecSPS can get arbitrarily small\footnote{Take for instance the deterministic setting one-dimensional setting $f(x) = |x|$. As $x\to 0$, the stepsize prescribed by DecSPS converges to zero. This is not the case e.g. in the quadratic setting.} for finite $k$. One easy solution to the problem is to enforce a lower bound, and adopt a new proof technique. Specifically we propose the following:
\begin{tcolorbox}
\begin{equation}
\tag{DecSPS-NS}
{\small
    \gamma_k:=\frac{1}{c_k}\cdot\min\left\{\max\left\{c_0\gamma_{\ell}, \frac{f_{\SC_k}(x^k)-\ell^*_{\SC_k}}{\|g_{\SC_k}(x^k)\|^2}\right\}, c_{k-1}\gamma_{k-1}\right\},}
\end{equation}
\end{tcolorbox}
where $c_k\ne 0$ for every $k\ge0$, $\gamma_\ell\ge\gamma_b$ is a small positive number and all the other quantities are defined as in DecSPS. In particular, as for DecSPS, we set $c_{-1}=c_0$ and $\gamma_{-1} = \gamma_b$. Intuitively, $\gamma_k$ is selected to live in the interval $[c_0\gamma_\ell/c_k, c_0\gamma_b/c_k]$~(see proof in~\S\ref{sec:app_proof4}, appendix), but has subgradient-dependent adaptive value. In addition, this stepsize is enforced to be monotonically decreasing.

\begin{restatable}[]{thm}{decSPScvxNS}
\label{thm:SPS_bounded_domain_cvx_NS} For any non-decreasing positive sequence $(c_k)_{k=0}^\infty$, consider SGD with DecSPS-NS. Assume that each $f_i$ is convex and lower bounded. We have
\vspace{-2mm}
\begin{equation}
 \Exp[f(\bar{x}^K)-f(x^*)]\leq \frac{c_{K-1} D^2}{\gamma_\ell c_0 K } + \frac{1}{K} \sum_{k=0}^{K-1}\frac{c_0 \gamma_b G^2}{c_k} ,
\end{equation}
where $D^2 :=\max_{k \in [K-1]}\|x^k-x^*\|^2$ and $G^2:=\max_{k \in [K-1]}\|g_{\SC_k}(x^k)\|^2$.
\end{restatable}
\vspace{-2mm}
One can then easily derive an $\mathcal{O}(1/\sqrt{k})$ convergence rate. This is presented in \S\ref{sec:app_proof4}~(appendix).

\section{Numerical Evaluation}
\label{sec:exp}

We evaluate the performance of DecSPS with $c_k = c_0\sqrt{k+1}$ on binary classification tasks, with regularized logistic loss $f(x) = \frac{1}{n}\sum_{i=1}^n \log(1+\exp(y_i\cdot a_i^\top x))+\frac{\lambda}{2} \|x\|^2$, where $a_i\in\mathbb{R}^d$ is the feature vector for the $i$-th datapoint and $y_i\in\{-1,1\}$ is the corresponding binary target.

We study performance on three datasets: (1) a Synthetic
Dataset, (2) The A1A dataset~\citep{chang2011libsvm} and (3) the
Breast Cancer dataset~\citep{Dua:2019}. We choose
different regularization levels and batch sizes bigger than 1.
Details are reported below, and the code is available at~\url{https://github.com/aorvieto/DecSPS}.

\begin{itemize}[leftmargin=*]
  \setlength\itemsep{0.01em}
    \item \textit{Synthetic Dataset} : Following~\citep{gower2019sgd} we generate $n = 500$ datapoints from a standardized Gaussian distribution in $\R^d$, with $d=100$. We sample the corresponding labels at random. We consider a batch size $B=20$ and either $\lambda=0$ or $\lambda = 1e-4$.
    \item \textit{A1A dataset}~(standard normalization) from~\citep{chang2011libsvm}, consisting in $1605$ datapoints in 123 dimensions. We consider again $B=20$ but a substantial regularization with $\lambda = 0.01$.
    \item \textit{Breast Cancer dataset}~(standard normalization)~\citep{Dua:2019}, consisting in $569$ datapoints in 39 dimensions. We consider a small batch size $B=5$ with strong regularization $\lambda = 0.1$.
\end{itemize}

At the batch sizes and regularizer levels we choose, the problems do not satisfy interpolation. Indeed, running full batch gradient descent yields $f^*>0$. While running SPS$_{\max}$ on these problems (1) does not guarantee convergence to $f^*$ and (2) requires full knowledge of the set of optimal function values $\{f^*_\SC\}_{|\SC|=B}$, in DecSPS we can simply pick the lower bound $0=\ell^*_\SC \le f^*_\SC$ for every $\SC$. Supported by Theorems~\ref{thm:SPS_bounded_domain_cvx}~\&~\ref{thm:no_bound}~\&~\ref{thm:SPS_bounded_domain_cvx_NS}, we expect SGD with DecSPS to converge to the minimum $f^*$.

All experiments reported below are repeated 5 times. Shown is mean and 2 standard deviations.

\paragraph{Stability of DecSPS.} DecSPS has two hyperparameters: the upper bound $\gamma_b$ on the first stepsize and the scaling constant $c_0$. While Thm.~\ref{thm:SPS_bounded_domain_cvx_NS} guarantees convergence for any positive value of these hyperparameters, the result of Thm.~\ref{thm:SPS_bounded_domain_cvx} suggests that using $c_0=1$ yields the best performance under the assumption that $\hat \sigma^2_B\ll \tilde L D^2$~(e.g. reasonable distance of initialization from the solution, and $L_{\max} > 1/\gamma_b$). In Fig.~\ref{fig:sensitivity}, we show that on the synthetic dataset that (1) $c_0=1$ yields indeed the best performance at $\gamma_b=10$ and (2) the dynamics of SGD with DecSPS is almost independent of $\gamma_b$ at $c_0=1$. Similar findings hold for the A1A dataset. For the Breast Cancer dataset, however, we see that at $\gamma_b=10$, $c_0=5$ yields the best final suboptimality --- yet $c_0=1$ is clearly the best tradeoff between convergence speed and final accuracy.

\ \\Hence, \textit{for all further experiments}, we choose the hyperparameters $\gamma_b=10, c_0=1$.

\begin{figure}
    \centering

    \textbf{\scriptsize\qquad\quad\quad  Synthetic - Sensitivity $\bf{\gamma_b}$\qquad\quad A1A - Sensitivity $\bf{\gamma_b}$\qquad\quad Breast - Sensitivity $\bf{\gamma_b}$\qquad\quad\qquad }\\

    \includegraphics[height=0.25\textwidth]{fig/sensitivity/Toy_Logistic_regression_SPS_sensitivity_gamma_b_1.png}
    \includegraphics[height = 0.25\textwidth]{fig/real_new/A1A_SPS_sensitivity_gamma_b.png}
    \includegraphics[height = 0.25\textwidth]{fig/real_new/Breast_SPS_sensitivity_gamma_b.png}\\

     \textbf{\scriptsize\qquad\quad\quad  Synthetic - Sensitivity $\bf{c_0}$\qquad\quad A1A - Sensitivity $\bf{c_0}$\qquad\quad Breast - Sensitivity $\bf{c_0}$\qquad\quad\qquad }\\
     \hspace{1.7mm}
    \includegraphics[height=0.25\textwidth]{fig/sensitivity/Toy_Logistic_regression_SPS_sensitivity_c0_1.png}
    \includegraphics[height = 0.25\textwidth]{fig/real_new/A1A_SPS_sensitivity_c0.png}
    \includegraphics[height = 0.25\textwidth]{fig/real_new/Breast_SPS_sensitivity_c0.png}
    \vspace{-2mm}
    \caption{\small DecSPS~($c_k=c_0\sqrt{k+1}$) sensitivity to hyperparameters on the Synthetic, A1A and Breast cancer datasets. Experiments are repeated 3 times and plotted is mean and std.}
    \label{fig:sensitivity}
\end{figure}

\paragraph{Comparison with vanilla SGD with decreasing stepsize.} We compare the performance of DecSPS against the classical decreasing SGD stepsize $\eta/\sqrt{k+1}$, which guarantees convergence to the exact solution at the same asymptotic rate as DecSPS. We show that, while the asymptotics are the same, DecSPS with the theoretically inspired hyperparameters $c_0=1, \gamma_b=10$ performs competitively to a fine-tuned $\eta$ --- where crucially the optimal value of $\eta$ depends on the problem. This behavior is shown on all the considered datasets, and is reported in Figure~\ref{Fig_SGDVsDecSPS}.

\begin{figure}
    \centering
\textbf{\scriptsize Synthetic Logistic Regression -- SGD Vs DecSPS \quad \quad  \quad \quad \quad  \quad  Regularized A1A -- SGD Vs DecSPS}\\
    \includegraphics[height=0.25\linewidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_SGD.pdf}
    \includegraphics[height=0.25\linewidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_SGD_stepsize.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_stepsize.pdf}
    \textbf{\scriptsize \qquad  \qquad Regularized Breast Cancer -- SGD Vs DecSPS \qquad}\\
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_stepsize.pdf}
    \vspace{-2mm}
    \caption{\small DecSPS with hyperparameters $c_0=1, \gamma_b=10$ compared to vanilla SGD with stepsize $\eta/\sqrt{k+1}$. DecSPS performs competitively to a fine-tuned $\eta$ --- where crucially the optimal value of $\eta$ depends on the problem.}
    \label{Fig_SGDVsDecSPS}
    \vspace{-2mm}
\end{figure}

\begin{figure}
    \centering
\textbf{\scriptsize \quad  Synthetic Dataset -- DecSPS Vs AdaGrad-Norm\qquad \quad Regularized A1A -- DecSPS Vs AdaGrad-Norm \qquad}\\
\includegraphics[height=0.25\linewidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada.pdf}
    \includegraphics[height=0.25\linewidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada_stepsize.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_ADA.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_stepsize_ADA.pdf}\\
    \textbf{\scriptsize \qquad \qquad Regularized Breast Cancer -- DecSPS Vs AdaGrad-Norm \qquad }\\
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_ADA_5.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_stepsize_ADA_5.pdf}
    \vspace{-2mm}
    \caption{\small Performance of AdaGrad-Norm compared to DecSPS. After fixing $b_0=0.1$ for Adagrad-Norm as recommended in~\citep{ward2019adagrad}~(see their Figure 3), tuning $\eta$ cannot quite match the performance of DecSPS with theoretically inspired hyperparameters $c_0=1, \gamma_b=10$.}
    \label{Fig_SGDVsADA_sklearn}
\end{figure}

\begin{figure}[ht]
    \centering
    \textbf{\scriptsize\quad \quad Synthetic Logistic Regression, $\boldsymbol{b_0 = 0.05}$\qquad\quad\quad \quad\quad  Synthetic Logistic Regression, $\boldsymbol{b_0 = 0.5}$}\\
    \includegraphics[height = 0.25\textwidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada_tuning1.pdf}
    \includegraphics[height = 0.25\textwidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada_stepsize_tuning1.pdf}
    \includegraphics[height = 0.25\textwidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada_tuning2.pdf}
    \includegraphics[height = 0.25\textwidth]{fig/toy_new/Toy_Logistic_regression_SPS_vs_Ada_stepsize_tuning2.pdf}
    \vspace{-2mm}
    \caption{\small Performance of AdaGrad-Norm compared to DecSPS on the Synthetic dataset, for $b_0\ne 0.1$. This figure is a complement to Figure~\ref{Fig_SGDVsADA_sklearn}.}
    \label{fig:synthetic_ada_tuning}
\end{figure}

\paragraph{DecSPS yields a truly adaptive stepsize.} We inspect the value of $\gamma_k$ returned by DecSPS, shown in Figure~\ref{Fig_SGDVsDecSPS}. Compared to the vanilla SGD stepsize $\eta/\sqrt{k+1}$, a crucial difference appears: $\gamma_k$ decreases faster than $\mathcal{O}(1/\sqrt{k})$. This showcases that, while the factor $\sqrt{k+1}$ can be found in the formula of DecSPS\footnote{We pick $c_k=c_0\sqrt{k+1}$, as suggested by Cor.~\ref{cor:sqrt_smooth}~\&~\ref{cor:sqrt_NS}. }, the algorithm structure provides additional adaptation to curvature. Indeed, in~(regularized) logistic regression, the local gradient Lipschitz constant increases as we approach the solution. Since the optimal stepsize for steadily-decreasing SGD is $1/(L\sqrt{k+1})$, where $L$ is the global Lipschitz constant~\cite{ghadimi2013stochastic}, it is pretty clear that $\eta$ should be decreased over training for optimal converge~(as $L$ effectively increases). This is precisely what DecSPS is doing.

\paragraph{Comparison with AdaGrad stepsizes.} We also compare DecSPS with another adaptive coordinate-independent stepsize with strong theoretical guarantees: the norm version of AdaGrad~(a.k.a. AdaGrad-Norm, AdaNorm), which guarantees the exact solution at the same asymptotic rate as DecSPS~\cite{ward2019adagrad}. AdaGrad-norm at each iteration updates the scalar $b_{k+1}^2=b_{k}^2+\|\nabla f_{\SC_k}(x_k)\|^2$ and then selects the next step as $x_{k+1}=x_{k}-\frac{\eta}{b_{k+1}}\nabla f_i(x_k)$. Hence, it has tuning parameters $b_0$ and $\eta$. In Figure~\ref{Fig_SGDVsADA_sklearn} we show that, after fixing $b_0=0.1$ as recommended in~\citep{ward2019adagrad}~(see their Figure 3), tuning $\eta$ cannot quite match the performance of DecSPS. In Fig.~\ref{fig:synthetic_ada_tuning}, we show that further tuning of $b_0$ likely does not yield a substantial improvement.

\paragraph{Comparison with Adam and AMSgrad.} We provide a comparison with Adam~\citep{kingma2014adam} and AMSgrad~\citep{reddi2019convergence}. For both methods, we set the momentum parameter to zero~(a.k.a RMSprop) for a fair comparison with DecSPS. For $\beta :=\beta_2$, the parameter that controls the moving average of the second moments, we select the value $0.99$ since we found that the standard $0.999$ leads to problematic~(exploding) stepsizes. Findings are pretty similar for both the A1A~(Figure~\ref{Fig_SGDVsAdam_A1A}) and Breast Cancer~(Figure~\ref{Fig_SGDVsAdam_Breast}) datasets: when compared to DecSPS with the usual parameters, fine-tuned Adam with fixed stepsize can reach the same performance after a few tens of thousand iterations --- however, it is slower at the beginning of training. While deriving convergence guarantees for Adam is problematic~\citep{reddi2019convergence}, AMSgrad~\citep{reddi2019convergence} with stepsize $\eta/\sqrt{k+1}$ enjoys a convergence guarantee similar to Adagrad and Adagrad-Norm. This is reflected in the empirical convergence: fine-tuned AMSgrad is able to match the convergence of DecSPS with the usual parameters motivated at the beginning of this section. Yet, we recall that the convergence guarantees of AMSgrad require the iterates to live in a bounded domain, an assumption which is not needed in our DecSPS~(see \S~\ref{sec:no_bound}).

\begin{figure}[ht]
    \centering
\textbf{\scriptsize \qquad \qquad Regularized A1A -- DecSPS Vs Adam \qquad \qquad \quad Regularized A1A -- DecSPS Vs AMSgrad}\\
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_Adam.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_stepsize_Adam.pdf}
\includegraphics[height=0.25\linewidth]{fig/real_new/A1A_AMSgrad.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/A1A_stepsize_AMSgrad.pdf}
    \vspace{-2mm}
\caption{\small Performance of Adam~(with fixed stepsize and no momentum) and AMSgrad~(with sqrt decreasing stepsize and no momentum) compared to DecSPS on the A1A dataset. Plotted is also the average stepsize~(each parameter evolves with a different stepsize).}
    \label{Fig_SGDVsAdam_A1A}
\end{figure}

\begin{figure}[h!]
    \centering
\textbf{\scriptsize \qquad \qquad Regularized Breast Cancer -- DecSPS Vs Adam \qquad  Regularized Breast Cancer -- DecSPS Vs AMSgrad}\\
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_Adam.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_stepsize_Adam.pdf}
\includegraphics[height=0.25\linewidth]{fig/real_new/Breast_AMSgrad.pdf}
    \includegraphics[height=0.25\linewidth]{fig/real_new/Breast_stepsize_AMSgrad.pdf}

    \vspace{-2mm}
\caption{\small Performance of Adam~(with fixed stepsize and no momentum) and AMSgrad~(with sqrt decreasing stepsize and no momentum) compared to DecSPS on the Breast Cancer dataset. Plotted is also the average stepsize~(each parameter evolves with a different stepsize).}    \label{Fig_SGDVsAdam_Breast}
\end{figure}

\section{Conclusions and Future Work}
We provided a practical variant of SPS~\citep{loizou2021stochastic} which converges to the true problem solution without the interpolation assumption in convex stochastic problems --- matching the rate of AdaGrad. If in addition strong convexity as assumed, then we show how, in contrast to current results for AdaGrad, the bounded iterates assumption can be dropped. The main open direction is a proof of a faster rate $\mathcal{O}(1/K)$ under strong convexity. There are also other possible extensions of our work: using the proposed new variants of SPS with accelerated methods, studying further the effect of mini-batching and non-uniform sampling of DecSPS, and extensions to the distributed and decentralized settings.

\section*{Acknowledgements}
This work was partially supported by the Canada CIFAR AI Chair Program. Simon Lacoste-Julien is a CIFAR Associate Fellow in the Learning in Machines \& Brains program.
\newpage