text:
\begin{document}

\title{High probability convergence bounds for stochastic gradient descent assuming the Polyak-\L ojasiewicz inequality}

\author{Liam Madden\thanks{\url{liam.madden@colorado.edu}, Dept.\ of Applied Mathematics, University of Colorado Boulder}  \and  Emiliano Dall'Anese\thanks{\url{emiliano.dallanese@colorado.edu}, Dept.\ of Electrical, Computer and Energy Engineering, University of Colorado Boulder} \and Stephen Becker\thanks{\url{stephen.becker@colorado.edu}, Dept.\ of Applied Mathematics, University of Colorado Boulder}   }

\date{\today}

\maketitle

\begin{abstract}
Stochastic gradient descent (with a mini-batch) is one of the most common iterative algorithms used in machine learning. While being computationally cheap to implement, recent literature suggests that it may also have implicit regularization properties that prevent overfitting. This paper analyzes the properties of stochastic gradient descent from a theoretical standpoint to help bridge the gap between theoretical and empirical results. Assuming smoothness, the Polyak-\L ojasiewicz inequality, and sub-gaussian stochasticity, we prove high probability bounds on the convergence rate, matching existing expected bounds. As an application of our convergence results, we combine them with existing uniform stability and generalization bounds to bound the true risk of the iterates of stochastic gradient descent. We find that for a certain number of epochs, the convergence and generalization balance in such a way that the true risk bound goes to the empirical minimum as the number of samples goes to infinity.
\end{abstract}

\section{Introduction}

Mini-batch stochastic gradient descent (SGD) is a popular and effective choice for solving nonconvex optimization problems in the machine learning domain. In its simplest form, SGD minimizes a given function by relying on noisy estimates of the gradients; for example, in a  statistical learning setting, estimates of the gradient can be obtained by processing a sample point from a dataset or a mini-batch of randomly selected samples (as opposed to utilizing all the points, as in the plain-vanilla gradient descent).
This paper proves high probability convergence bounds for SGD in the general setting, then focuses on a parametrized  learning model and addresses the problem of minimizing the empirical risk associated with a given dataset.

To help bridge the gap between theoretical and empirical results, in the analysis of SGD we pay attention to a number of practical aspects and we provide contributions towards relaxing a number of assumptions on the empirical risk. For example, we derive convergence results (in both mean and high-probability) assuming that the cost associated with the  the empirical risk is smooth and satisfies the Polyak-\L ojasiewicz (PL) inequality,
improving on prior works that unrealistically assumed strong convexity.
For example, while the problem of training neural networks is nonconvex, \cite{AZ1,AZ2} show that for sufficiently wide neural nets, the empirical risk is locally semi-smooth and self-bounded and locally satisfies the PL inequality with a $\Omega(1/n^2)$ constant where $n$ is the number of data points (note that self-boundedness is a consequence of nonnegativity and smoothness); further, the PL inequality is satisfied for some neural networks with linear activations, as discussed by~\cite{Charles}. The results are also applicable to machine learning applications where the the empirical risk is convex and satisfies the PL inequality~\citep{Karimi}.
We also analyze SGD as it is often used in practice: without projection steps, and without iterate averaging.

\paragraph{Contributions at-a-glance.}
First, we derive a concentration inequality for a particular stochastic recursion. Our analysis is inspired by Thm. 4.1 of \cite{Harvey1}. Next, we use our concentration inequality to derive high probability convergence bounds for an objective that is smooth and satisfies the PL inequality, and assuming sub-gaussian noise. The first bound matches the lower bound given in \cite{Nemirovski} in the setting where the PL inequality is actually strengthened to strong convexity. The second bound extends the expected bounds of \cite{Orvieto} to high probability. Finally, we apply existing uniform stability and generalization bounds, along with our second convergence bound, to obtain a true risk bound that goes to the empirical minimum as the number of samples increases to infinity. Thus, our theory supports taking multiple passes of the data, but stopping early instead of running to convergence. A more in-depth comparison with prior works will be provided in Section~\ref{sec:Priorwork}.

\paragraph{Notation.} Let $[n]\defeq \{1,2,\ldots,n\}$. We use

$B(n,p)$ to denote the binomial distribution with $n$ trials and probability of success $p$. By convention, we define $\sum_{i=a}^b c_i=0$ if $b<a$ and $\prod_{i=a}^b c_i=1$ if $b<a$.

\section{Preliminaries}

This paper deals with the stochastic optimization problem
\begin{align}
\label{eq:stochprob}
    \min_x f(x)\defeq \E[ f_{\xi}(x)]
\end{align}
where $\xi$ is a random variable. We define the SGD iteration by
\begin{align}
\label{eq:sgd}
    x_{t+1} = x_t-\eta_t g_t
\end{align}
for $g_t=g(x_t,b_t)$ where $g(x,b)=\frac{1}{b}\sum_{i=1}^b \nabla f_{\xi_i}(x)$ and the $\xi_i$ are independently sampled. Thus, $\E [g(x,b)]=\nabla f(x)$. A common assumption on $g$ is the \textit{bounded variance} property: there exists $\sigma<\infty$ such that $\E[\|\nabla f(x)-g(x,1)\|^2]\leq\sigma^2$ for all $x$ \citep{Karimi,Nemirovski,Khaled}. While this is a sufficient assumption for expected convergence, we need a stronger assumption for high probability convergence. In particular, we will assume that $\nabla f(x)-g(x,1)$ is a sub-gaussian random vector with parameter $\sigma/\sqrt{d}$. See Appendix \ref{sec:probability} for details. Note that this assumption implies the bounded variance property with the same constant.

In addition to the sub-gaussian assumption on the noise, we assume $\argmin_{x\in\R^d}f(x)$ is nonempty and let $f\opt$ denote the minimum, and we assume $f$ is smooth and satisfies the PL inequality.

A function $f$ is $L$-smooth if $f$ is $C^1$ and its gradient is $L$-Lipschitz continuous. Further details on standard optimization terminology and results are in Appendix~\ref{subsec:optimizationResults}.

A function $f$ is $\mu$-Polyak-\L ojasiewicz (or $\mu$-PL for short) if $(\forall x \in \R^d)$ $\|\nabla f(x)\|^2\geq 2\mu(f(x)-f\opt)$. As an example, the logistic regression cost function satisfies the PL inequality \citep{Karimi}. Also, deep linear neural networks satisfy the PL inequality in large regions of the parameter space \citep[Thm.\ 4.5]{Charles}. Furthermore, sufficiently wide recurrent neural networks satisfy the PL inequality locally around random initialization \citep{AZ1,AZ2}. While strong convexity implies the PL inequality, a PL function need not even be convex, hence the PL condition is considerably more applicable in the context of neural networks. More details on PL, including a discussion of projections, is in Appendix~\ref{subsec:PL}.

Now, a particular example of Eq.~\eqref{eq:stochprob} that we are concerned with is the supervised statistical learning problem~\citep{Shai,Mohri}. Consider a dataset $S \defeq (s_i)_{i=1}^n\sim \D^n$ consisting of $n$ iid samples, $s_i\sim\D$ for $i\in [n]$, where $\D$ is an arbitrary and unknown distribution. Let $\ell(x,s)$ measure the loss of the model parameterized by $x$ applied to the sample point $s$. Given $S$, let $f_i=\ell(\cdot,s_i)$ and define the \emph{empirical risk}
\begin{equation}\label{eq:ERM}
f(x) \defeq \frac{1}{n}\sum_{i=1}^n f_i(x).
\end{equation}
The ultimate goal is to find an $x$ that minimizes the following \emph{true risk}:
\begin{equation}\label{eq:TR}
F(x) \defeq \E_{s\sim\D}~\ell(x,s).
\end{equation}
This example fits into Eq.~\eqref{eq:stochprob} by defining $\xi_i=s_i$ with probability $1/n$ and setting $f_{s_i}=\ell(\cdot,s_i)$.

Applying SGD with fewer  than $n$ iterations essentially minimizes $F$ directly and can be analyzed via stochastic approximation \citep{Nemirovski}.
However, most realistic scenarios have a finite $n$ and take more than $n$ iterations in order to further reduce the training error. To make generalization guarantees in this scenario, two popular approaches are to either assume a limited model $x\in\mathcal{H}$ and bound the complexity of $\mathcal{H}$, or to show that the algorithm is stable. As it is widely established that successful deep learning models have extremely large VC-dimension, standard bounds based on the first approach are usually vacuous for realistic values of $n$. For this reason, we pursue the second approach, and will analyze the stability of SGD in order to bound the true risk in terms of the empirical risk.

Formally, for $T$ iterations of SGD we define the convergence error as $f(x_T)-f\opt$ and the generalization error as $F(x_T)-f(x_T)$. We will balance the sum of the convergence error and generalization error, which is $F(x_t)-f\opt$. Note that $F(x_t)-\min_x F(x)=f\opt-\min_x F(x)+F(x_t)-f\opt$. The rate at which the $f\opt-\min_x F(x)$ term goes to zero as $n\to\infty$ depends on $\D$, not the algorithm, hence why we focus on bounding $F(x_t)-f\opt$.

For stability, we assume $\ell(\cdot,s)$ is Lipschitz continuous and smooth uniformly over all $s$. For generalization, we further assume that $\ell(x,s)$ is nonnegative and bounded uniformly over all $x$ and $s$. However, these assumptions are interrelated, as the following proposition shows (proof in Appendix \ref{proofs:interrelated}).

\begin{proposition} \label{prop:interrelated}
Define $f$ as in Eq.~\eqref{eq:ERM} and assume it is $\mu$-PL. Assume $\ell(\cdot,s)$ is  $L$-smooth for all $s$. Then $f$ is $L$-smooth. Moreover, if $\nabla f(x)-g(x,1)$ is sub-gaussian and $\ell(x,s)$ is nonnegative for all $x$ and $s$, and assuming interpolation ($f\opt=0$), then for $(x_t)$ generated by SGD with a $\order(1/t)$ step-size sequence there exists a constant $R<\infty$ such that $\P\left(\|x_t-x_0\|\leq R~\forall t\in\N\right)=1$. If $s$ comes from a compact set (e.g., normalized data) and we assume $\ell(x,s)$ is jointly continuous in $x$ and $s$, then there exist constants $\rho$ and $M$ such that, for all $x$ in the ball of radius $R$ around $x_0$ and for all $s$, $\ell(x,s)$ is $\rho$-Lipschitz with respect to $x$ and is bounded above by $M$.
\end{proposition}

\section{Prior work}
\label{sec:Priorwork}

First, we will discuss prior high probability convergence results in the convex non-smooth setting, the convex smooth setting, and the non-convex smooth setting. To the best of our knowledge, there are not high probability upper bounds in the non-convex PL smooth setting.

The lower high probability bound for the final iterate of projected SGD assuming only Lipschitz continuity, convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)\log(t)/\sqrt{t})$. The lower high probability bound for the final iterate of projected SGD assuming only Lipschitz continuity, strong convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)\log(t)/t)$. The lower high probability bound for iterate averaging of projected SGD assuming only Lipschitz continuity, convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)/\sqrt{t})$. The lower high probability bound for iterate averaging of projected SGD assuming only Lipschitz continuity, strong convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)/t)$. Furthermore, upper high probability bounds matching each of the lower bounds have been shown for appropriate step-size sequences and iterate averaging strategies \cite{Harvey1,Harvey2,Nemirovski,Rakhlin}.

The lower high probability bound for the final iterate of projected SGD assuming only smoothness, convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)/\sqrt{t})$. The lower high probability bound for the final iterate of projected SGD assuming only smoothness, strong convexity, and sub-gaussian noise, is $\Omega(\log(1/\delta)/t)$.

On the other hand, in the general non-convex setting, iterating averaging does not make sense and the error has to be in terms of $\min_{t\in[T]}\|\nabla f(x_t)\|^2$. Smoothness and relaxed growth (which is weaker than the bounded variance property) are enough to show almost sure convergence of SGD \cite{Bertsekas,patel2020stopping}. Smoothness and sub-gaussian noise are enough to show the high probability bound $O(\log(t/\delta)\log(t)/t)$ \cite{li2020high}. This is optimal, at least up to $\log(t)$ factors, by the smooth/convex/sub-gaussian lower bound.

If we additionally assume the PL inequality, which does not imply convexity, but is implied by strong convexity, then we can again use the error $f(x_t)-f\opt$. Our first result is the upper high probability bound for SGD assuming smoothness, the PL inequality, and sub-gaussian noise of $O(\log(1/\delta)/t)$. This is optimal by the smooth/strongly-convex/sub-gaussian lower bound.

Note that a function cannot be globally Lipschitz continuous and strongly convex (unless it is essentially a constant function), which is why \emph{projected} SGD is used in such a setting. The same is true of the PL inequality instead of strong convexity, but we no longer have the luxury of projection as Remark \ref{rmk:pl} in Appendix~\ref{subsec:PL} shows. This precludes us from assuming Lipschitz continuity in our setting.

While our high probability convergence bounds seem to be the first in the non-convex PL smooth setting, there are expected convergence bounds in the literature. In particular, \cite{Orvieto} prove a $\order(1/t^a)$ expected convergence bound for $\order(1/t^a)$ step-size sequence where $a\in(0,1)$ \citep[Cor. 3]{Orvieto}. Similarly, \cite{Karimi} prove a $\order(1/t)$ expected convergence bound for a $\order(1/t)$ step-size sequence. Both of these upper bounds assume only bounded variance, not sub-gaussian noise. By assuming sub-gaussian noise, we are able to extend the expected bound that can be derived from Thm. E.3 of \cite{Orvieto} to a high probability bound.

To get our high probability convergence bounds, we prove a concentration inequality with similar analysis to Theorem 4.1 of \cite{Harvey1}. They use Theorem 4.1 to unravel a recursion with a square root of the random variable in it. They are then able to bound an error term that is in terms of the iterates. Finally, they apply a generalized version of Freedman's inequality \citep{Freedman}. On the other hand, since we have the PL inequality instead of strong convexity, we don't have contractive properties at our disposal to work directly with the iterates. But, using smoothness, we do get a recursion with a square root of the random variable in the moment generating function. For concentration inequality references, see \cite{Boucheron} and \cite{Fan}.

Regarding stability, \cite{Hardt} prove an expected stability bound for SGD with batch-size 1 assuming Lipschitz continuity and smoothness. \cite{Mou} also investigate expected stability, but for SGD with Langevin dynamics. Their bound ends up being much better than the one of \cite{Hardt}, but the noise decays too slowly to get fast convergence.

In a different vein, \cite{Charles} prove hypothesis stability results. Their bounds are the sum of a convergence term and a geometry term that comes from the PL inequality. Thus, there is no benefit to early stopping. But, to get from their stability bound to a generalization bound, they have to make a strict assumption (Asm. 1) because the standard tool, Theorem 11 of \cite{Bousquet}, uses hypothesis stability, which is stronger than pointwise hypothesis stability.

The seminal paper of \cite{Bousquet} shows how to turn stability bounds into generalization bounds for deterministic algorithms. \cite{Elisseeff} extend those bounds to randomized algorithms, as do \cite{Shalev}. \cite{Feldman2} prove high probability generalization bounds for SGD, but only in the convex setting. They also balance their generalization result with the convergence result of \cite{Harvey1}.

We extend the result of \cite{Hardt} to mini-batch stochastic gradient descent and apply the result of \cite{Elisseeff} to it, proving a low probability generalization bound. We combine our convergence bound with the generalization bound to get a low probability true risk bound. Assuming the constants in our assumptions do not change with the number of samples (which is an admittedly strong assumption), our true risk bound goes to the empirical minimum as the number of samples goes to infinity. Furthermore, it has the form of taking multiple, but limited, epochs of SGD. We believe this is a first step in finding sharper true risk bounds that support multiple, limited epochs of SGD.

\section{Convergence}
\label{sec:convergence}

In this section, we provide high probability convergence bounds assuming smoothness, the PL inequality, and sub-gaussian gradient noise. First, Proposition \ref{prop:harvey} gives a concentration inequality for a specific stochastic recursion. Theorem \ref{thm:fastconvergence} then leverages Proposition \ref{prop:harvey} to provide a convergence bound for a specific step-size sequence. Furthermore, the bound is optimal. Then, Theorem \ref{thm:slowconvergence} applies Proposition \ref{prop:harvey} to general $\order(1/t)$-decaying step-size sequences to be used for balancing convergence with stability/generalization. We start with the concentration inequality (full proof in Appendix \ref{subsec:convergence}).

\vspace{.1cm}

\begin{proposition}
\label{prop:harvey}
Let $(\Omega,\F,P)$ be a probability space and let $\{\F_t\}_{t\geq 0}$ be a filtration and $\F_{-1}=\emptyset$. For each $t\geq 0$, let $\alpha_t\ge 0$, $\beta_t \geq 0$, and $\gamma_t \ge 0$ be arbitrary scalars; and
$X_t$, $\wh_t$, and $\vh_t$ be random variables on $(\Omega,\F,P)$, with $X_t$ $\F_{t-1}$-measurable and nonnegative almost surely, and $\wh_t$ and $\vh_t$ both $\F_t$-measurable, and $\vh_t$  nonnegative almost surely.
Consider the recursion
\[ X_{t+1}\leq \alpha_t X_t+\wh_t+\vh_t
\]
and assume $\E[\wh_t\mid\F_{t-1}]=0$,  $\E[\exp(\lambda \wh_t) \mid \F_{t-1}]\leq \exp(\frac{\lambda^2}{2}\beta_t^2 X_t)$ for all $\lambda\in\R$, and $\E[\exp(\lambda \vh_t) \mid \F_{t-1}]\leq \exp(\lambda \gamma_t)$ for all $\lambda\in [0,\gamma_t^{-1}]$. Then,
\begin{align} \label{eq:harveyPropExpected}
    \E X_T \leq \left(\prodt \alpha_t\right) X_0+\sumt \left(\prodi \alpha_i\right) \gamma_t.
\end{align}
Furthermore, let $(K_t)$ be a sequence such that $K_0\geq X_0$ and, for all $t\geq 0$, $K_{t+1}^2\geq (\alpha_t K_t+2\gamma_t)K_{t+1}+\beta_t^2 K_t$. Then, for all $T\geq 0$, we have, $\whp$ for all $\delta \in (0,1/e)$,
\begin{align} \label{eq:harveyPropHighProb}
    X_T \leq K_T\log(e/\delta).
\end{align}
\end{proposition}

The result boils down to a high probability bound on the recursion $X_{t+1}\leq\alpha_t X_t+\gamma_t$ with centered noise. The sub-gaussian noise term has the additional property that its sub-gaussian parameter depends on $\sqrt{X_t}$. For the SGD setting, the smoothness assumption implies the $\sqrt{X_t}$ sub-gaussian parameter.

We start with a lemma deriving the recursion for SGD, where we use $e_t$ to denote $\nabla f(x_t)-g_t$.

\begin{lemma}
\label{lem:convergence}
Assume $f$ is $L$-smooth and $\mu$-PL. Let $\eta_t\leq \frac{1}{L}~\forall t\geq 0$. Then,
\begin{align}
\begin{split}
\label{eq:ineqLmu}
    &f(x_{t+1})-f\opt \leq (1-\mu \eta_t) (f(x_t)-f\opt)+ \eta_t(1-L\eta_t)\langle \nabla f(x_t),e_t\rangle+\frac{L\eta_t^2}{2}\|e_t\|^2.
\end{split}
\end{align}
\end{lemma}
\begin{proof}
Observe,
\begin{align*}
    f(x_{t+1})-f\opt&\leq f(x_t)-f\opt -\eta_t(1-L\eta_t/2)\|\nabla f(x_t)\|^2+\eta_t(1-L\eta_t)\langle \nabla f(x_t),e_t\rangle+\frac{L\eta_t^2}{2}\|e_t\|^2\\
    &\leq f(x_t)-f\opt -\frac{\eta_t}{2}\|\nabla f(x_t)\|^2+\eta_t(1-L\eta_t)\langle \nabla f(x_t),e_t\rangle+\frac{L\eta_t^2}{2}\|e_t\|^2\\
    &\leq (1-\mu \eta_t) (f(x_t)-f\opt)+ \eta_t(1-L\eta_t)\langle \nabla f(x_t),e_t\rangle+\frac{L\eta_t^2}{2}\|e_t\|^2
\end{align*}
where the first inequality uses smoothness, the second inequality uses the step-size restriction, and the third inequality uses the PL inequality.
\end{proof}

We will apply Proposition \ref{prop:harvey} to Lemma \ref{lem:convergence} to derive convergence bounds in two cases. For the first case, define the sequence $(\theta_t)$ as
\begin{align}
    \theta_t &= \begin{cases}\frac{1}{L}  &\text{ for }t=0,...,\tau-1 \\ \frac{2t+1}{\mu(t+1)^2} &\text{ for }t=\tau,\tau+1,...\end{cases}
\end{align}
where $\tau=\lfloor 2L/\mu \rfloor$. Notice that $\theta_t\leq \frac{1}{L}~\forall t\geq 0$. The first convergence result is then provided next.

\begin{theorem}
\label{thm:fastconvergence}
Assume $f$ is $L$-smooth and $\mu$-PL and that $\nabla f(x)-g(x,1)$ is $\sigma/\sqrt{d}$-sub-gaussian. Then, SGD with $(\eta_t)=(\theta_t)$ and $b_t=b$ constructs a sequence $(x_t)$ such that
\begin{align}
    \E[f(x_t)-f\opt] &=\order\left(\frac{L\sigma^2}{\mu^2 b t}\right)
\end{align}
and, $\whp$ for all $\delta\in(0,1/e)$,
\begin{align}
    f(x_t)-f\opt
    &= \order\left(\frac{L\sigma^2\log(e/\delta)}{\mu^2 b t}\right).
\end{align}
\end{theorem}

A more detailed version of Theorem \ref{thm:fastconvergence} appears in Appendix~\ref{subsec:convergence} as Theorem \ref{thm:appendixconv}. Theorem \ref{thm:fastconvergence} provides an expected bound of $\order(1/t)$ and a high probability bound of $\order(\log(1/\delta)/t)$. This matches the optimal rates for the smooth and strongly convex setting \citep{Nemirovski}. The expected bound in Theorem \ref{thm:fastconvergence} can also be derived by applying Theorem E.3 of \cite{Orvieto} to $(\theta_t)$. However, to the best of our knowledge, the high probability bound is the first to match the expected bound.

In the following, we consider more general step-size sequences. While smaller constants result in sub-optimal convergence rates, they are also more stable (in the sense of Section \ref{sec:Stability}). The full proof is in Appendix \ref{subsec:convergence}.

\begin{theorem}
\label{thm:slowconvergence}
Assume $f$ is $L$-smooth and $\mu$-PL and that $\nabla f(x)-g(x,1)$ is $\sigma/\sqrt{d}$-sub-gaussian. Then, SGD with $\eta_t =  \frac{2c}{t+2}$, where $c< \frac{1}{L}$, and $b_t=b$ constructs a sequence $(x_t)$ such that
\begin{align}
    \E[f(x_t)-f\opt] &\leq \order\left(\frac{Lc^2\sigma^2}{b(1-\mu c)t^{\mu c}}\right).
\end{align}
Furthermore, $\whp$ for all $\delta\in(0,1/e)$, $\forall t\ge 0$,
\begin{align}
 f(x_t)-f\opt &\leq \order\left(\frac{Lc^2\sigma^2\log(e/\delta)}{b(1-\mu c)t^{\mu c}}\right).
\end{align}
\end{theorem}

In the following, we investigate stability and generalization.

\section{Stability and generalization}
\label{sec:Stability}

The foundations for the interplay between stability and generalization go back to \cite{Devroye}. Most recent results build off of \cite{Bousquet}. In particular, the results of \cite{Elisseeff} can be applied to the expected stability bound of \cite{Hardt} to get a low probability generalization bound. While this bound is not novel (except for the extension of it to mini-batch SGD), its pairing with Theorem \ref{thm:slowconvergence} produces a novel true risk bound that supports running multiple, but limited, epochs of SGD.

\begin{definition}
Let $S$ and $S'$ differ in at most one sample. Let $g$ and $g'$ be the respective gradient estimators. Define $x_0=x_0'$ and, for $t\geq 0$, $x_{t+1}=x_t-\eta_tg(x_t,b_t)$ and $x_{t+1}'=x_t'-\eta_t g'(x_t',b_t)$. We say SGD with $T$ iterations is \underline{$\alpha$-uniformly stable in expectation} if for all $S$ and $S'$, $\sup_s |\E [\ell(x_T,s)-\ell(x_T',s)]| \leq \alpha$.
\end{definition}

Note that
$\P(I_t=I_t')=\left(1-\frac{1}{n}\right)\cdots\left(1-\frac{1}{n-b_t+1}\right)=1-\frac{b_t}{n}$; \cite{Hardt} use this fact to bound the expected stability in the special case $b_t=1$. The proof of Thm.~\ref{thm:stability} below is in Appendix~\ref{subsec:stability}; recall that some of the hypotheses automatically hold under the setting of Prop.~\ref{prop:interrelated}.

\begin{theorem}
\label{thm:stability}
Let $\eta_t\leq c/(t+1)$ and $b_t=b$. Assume $\ell(\cdot,s)$ is $\rho$-Lipschitz and $L$-smooth for all $s$. Then,  SGD with $T$ iterations is $\alpha$-uniformly stable in expectation with
\begin{align}
\label{eq:alpha}
    \alpha \leq \frac{2\rho^2T^{Lc}(c+1/L)b}{n}.
\end{align}
\end{theorem}

Thm.~\ref{thm:stability} is essentially Thm. 3.12 of \cite{Hardt}. In fact, the two theorems consist of the same analysis with only two differences. First, we extend the analysis to mini-batch SGD, though this extension is not difficult. Second, \cite{Hardt} consider the first iteration where instability becomes non-zero. They use this to break up the bound and optimize over it. However, in so doing, they introduce the restriction $(2cL)^{1/(Lc+1)}T^{Lc/(Lc+1)}<n$. Thus, for simplicity, we left out this step. While it seems that makes our bound significantly larger, that is only when the restriction is not met (in which case their bound is vacuous).

Thms. 11 and 12 of \cite{Bousquet} provide generalization bounds from uniform stability bounds. The bounds are for deterministic algorithms, but can be extended to randomized algorithms as shown by \cite{Elisseeff}. We apply Thm. 2 of \cite{Elisseeff} (which we include as Theorem \ref{thm:elisseeff} in Appendix~\ref{subsec:stability}) to Eq.~\eqref{eq:alpha} to prove the following result.

\begin{theorem}
\label{thm:generalization2}
Assume $\ell(x,s)\in[0,M]$ for all $x$ and $s$. Assume $\ell(\cdot,s)$ is $\rho$-Lipschitz and $L$-smooth for all $s$. Then, $T$ iterations of SGD with $\eta_t\leq c/(t+1)$ and $b_t=b$ satisfies, $\whp$ over $S$ and $(I_t)$ for all $\delta\in(0,1)$,
\begin{align}
    \begin{split}&F(x_T)-f(x_T)\\
    &\leq \sqrt{\frac{12M\rho^2T^{Lc}(c+1/L)b+M^2}{2n\delta}}\end{split}\\
    &=\order\left(\sqrt{M\rho^2T^{Lc}(c+1/L)b/(n\delta)}\right).
\end{align}
\end{theorem}

Note that the generalization bound in Thm.~\ref{thm:generalization2} depends on $1/\sqrt{\delta}$, which makes it a ``low'' probability bound. Furthermore, it grows with the number of iterations. Fortunately, the convergence rate in Theorem \ref{thm:slowconvergence} is ``fast enough'' to allow us to balance convergence and generalization, by choosing $T$ to approximately minimize the bound, and obtain true risk bounds that go to the empirical minimum as $n \rightarrow \infty$. However, similarly to \cite{Charles}, one has to assume that $\mu$ is constant for all $n$. The bounds depend on SGD taking a certain number of epochs. Note that the boundedness and Lipschitz assumptions come for free under the conditions of Prop.~\ref{prop:interrelated}.

\begin{theorem}
\label{thm:hprisk}
Assume $\ell(x,s)\in[0,M]$ for all $x$ and $s$. Assume $\ell(\cdot,s)$ is $\rho$-Lipschitz and $L$-smooth for all $s$. Assume $f$ is $\mu$-PL. Let $\kappa=L/\mu$. Assume $\nabla f(x)-g(x,1)$ is $\sigma/\sqrt{d}$-sub-gaussian. Let $b_t=b$, $c=1/(\mu+L)$, and
\begin{align}
    T = \Theta\left(\frac{L^2d^2\sigma^4}{M\rho^2(\mu+L)^4b^2}\frac{n}{b}\right).
\end{align}
Then, $T$ iterations of SGD with $\eta_t= c/(t+1)$ satisfies, $\whp$ over $S$ and $(I_t)$ for all $\delta\in(0,1/e)$,
\begin{align}
\begin{split}
    &F(x_T)-f\opt\\
    &=\order\left( \left(\frac{L\sigma^2}{(\mu+L)^2b\sqrt{\delta}}\right)^{L/(\mu +L)}\left(\sqrt{\frac{Mb}{n\delta}}\rho\right)^{\mu /(\mu +L)}\right)\\
    &= \order\left(\frac{1}{b^{1-.75\kappa^{-1}}n^{\kappa^{-1}}\sqrt{\delta}}\right).
\end{split}
\end{align}
\end{theorem}

The proof of Theorem \ref{thm:hprisk} follows from plugging $c$ and $T$ into Theorems \ref{thm:slowconvergence} and \ref{thm:generalization2}. Note that the right-hand side goes to zero as $n \rightarrow \infty$, under the  assumption that $\mu$ does not change. The form of the bound is expected  for a (mini-batch) SGD: the main message is that one should take multiple passes over the data to converge sufficiently, but then stop early to avoid overfitting. It is also worth noticing that the bound decreases with the batch-size. Consequently, while increasing the batch-size increases the generalization bound, it decreases the convergence bound enough to decrease the true risk bound.

We also notice that while Theorem \ref{thm:hprisk} provides a low probability bound, since the low probability comes from the generalization and not the convergence, we can boost it to high probability by re-running SGD multiple times and selecting the solution with the smallest test error \citep[Exer. 13.1]{Shai}. Each round of SGD can use the same dataset $S$, and selecting the best solution only needs a small amount of test/validation data since this is independent.

\section{Conclusions}

In this paper, we proved high probability convergence bounds for stochastic gradient descent. Our optimization assumptions were smoothness and the PL inequality and our stochasticity assumption was sub-gaussian noise. In this setting, we were able to prove a $\order(\log(1/\delta)/t)$ bound with $\order(1/b)$ dependence on the mini-batch size. This is optimal since $\order(\log(1/\delta)/t)$ is the lower bound in the Lipschitz continuous and strongly convex setting. Also, this matches expected bounds in this setting. In the future, we would like to consider heavier tailed noise.

We do note that by considering the much weaker assumption of the PL inequality, we lose the ability to project, and hence to directly assume Lipschitz continuity. We leave it as a future research direction to determine what extra assumptions we could apply to the constraint set in order for projected SGD or stochastic Franke-Wolfe to converge.

We also note that the application to learning theory is quite preliminary. First of all, we assume that the PL constant depends favorably on the number of samples. Furthermore, the true risk bound has $\order(1/\sqrt{\delta})$ probability dependence instead of $\order(\log(1/\delta))$. This is because the expected uniform stability bound cannot be extended to high probability because the corresponding recursion does not concentrate sufficiently. We leave it as a future research direction to derive a recursion that utilizes the PL inequality in order to concentrate.

\section*{Acknowledgements}

All three authors gratefully acknowledge support from the NSF program ``AMPS-Algorithms for Modern Power Systems'' under award \# 1923298. Thank you to an anonymous reviewer who pointed out that the true risk was converging to the empirical, rather than true, minimum. Also, thank you to Aditya Golatkar for suggesting improvements in the convergence proofs.

\bibliographystyle{spbasic}
\bibliography{main}

\newpage

\appendix

\section{Probability results}
\label{sec:prob}

\subsection{Sub-gaussian properties}
\label{sec:probability}

First, we will review some definitions and results from \cite{Vershynin}. See \cite{Boucheron} as well. A random variable $X$ is called sub-gaussian with parameter $\sigma$ if
\begin{align}
    \E\exp(\lambda(X-\E X)) \leq \exp(\sigma^2\lambda^2)~\forall\lambda\in\R
\end{align}
and sub-exponential with parameter $\sigma$ if
\begin{align}
    \E\exp(\lambda|X|)\leq\exp(\sigma\lambda)~\forall\lambda\in[0,1/\sigma].
\end{align}
For a sub-gaussian random variable $X$, define its sub-gaussian norm to be
\begin{align}
    \|X\|_g=\inf\{\sigma>0\mid \E\exp(X^2/\sigma^2)\leq \exp(1)\}.
\end{align}
For a sub-exponential random variable $X$, define its sub-exponential norm to be
\begin{align}
    \|X\|_e=\inf\{\sigma>0\mid \E\exp(|X|/\sigma)\leq \exp(1)\}.
\end{align}
It is not difficult to verify that these are norms. Furthermore, we have the following facts:
\begin{itemize}[label={$\circ$}]
    \item A sub-gaussian random variable $X$ has parameter $\sigma$ if and only if $\|X\|_g\leq C\sigma$ for an absolute constant $C$ \citep[Prop. 2.5.2]{Vershynin}.
    \item A sub-exponential random variable $X$ has parameter $\sigma$ if and only if $\|X\|_e\leq C\sigma$ for an absolute constant $C$ \citep[Prop. 2.7.1]{Vershynin}.
    \item Hoeffding's inequality: If $X_1,\ldots,X_n$ are independent, centered, sub-gaussian random variables, then $S_n=\sum_{i=1}^n X_i$ is sub-gaussian and $\|S_n\|_g^2\leq C \sum_{i=1}^n \|X_i\|_g^2$ for an absolute constant $C$ \citep[Prop. 2.6.1]{Vershynin}.
    \item A random variable $X$ is sub-gaussian if and only if $X^2$ is sub-exponential. Moreover, $\|X\|_g^2=\|X^2\|_e$ \citep[Lem. 2.7.6]{Vershynin}.
\end{itemize}
Now, we would like to consider random vectors. A random vector $X$ is sub-gaussian with parameter $\sigma$ if $\langle u,X\rangle$ is sub-gaussian with parameter $\sigma$ for all unit vectors $u$ \citep[Def. 3.4.1]{Vershynin}. For a sub-gaussian random vector $X$, define
\begin{align}
    \|X\|_v = \sup_{\|u\|=1}\|\langle u,X\rangle\|_g.
\end{align}

Assume $\nabla f(x)-g(x,1)$ is $\sigma/\sqrt{d}$-sub-gaussian. Then $\nabla f(x)-g(x,b)$ is $C\sigma/\sqrt{bd}$-sub-gaussian. Furthmore, by Lemma 1 of \cite{Jin}, $\|\nabla f(x)-g(x,1)\|$ is $\sigma$-sub-gaussian. So, $\|\nabla f(x)-g(x,1)\|^2$ is $\sigma^2$-sub-exponential. Then, applying Jensen's inequality, we get that $\E\left[\|\nabla f(x)-g(x,1)\|^2\right]\leq \exp(1)$.

\subsection{Convergence of random variables}
\label{subsec:convofrvs}

When the rate of convergence is not an issue, a sequence of random variables $X_t$ converges to $X$ \citep[Sec. 20]{Billingsley}:
\begin{itemize}[label={$\circ$}]
    \item in probability if, for all $\epsilon>0$, $\lim_t \P(|X_t-X|>\epsilon)=0$, denoted $X_t\overset{p}{\to}X$;
    \item in mean if $\lim_t\E|X_t-X| =0$, denoted $X_t\overset{L^1}{\to}X$;
    \item almost surely if $\P(\lim_t X_t=X)=1$, denoted $X_t\overset{a.s.}{\to}X$.
\end{itemize}
On the other hand, when we care about the rate of convergence, we say a sequence of random variables $X_t$ converges to $X$:
\begin{itemize}[label={$\circ$}]
    \item with expected convergence rate $r$ if, for all $t$, $\E|X_t-X|\leq r(t)$;
    \item with $p$-probability convergence rate $r$ if, for all $t$ and $\delta$, $\P\left(|X_t-X|\leq r(t)p(\delta)\,\right)\geq 1-\delta$.
\end{itemize}
Specifically, we use ``high probability'' to refer to $p(\delta)=\order(\log(1/\delta))$ and ``low probability'' otherwise. All five kinds of convergence are interrelated:

\begin{itemize}[label={$\circ$}]
    \item Convergence in mean and convergence almost surely both imply convergence in probability.
    \item Expected convergence with rate $r$ such that $r(t)\to 0$ implies convergence in mean.
    \item By Markov's inequality, expected convergence with rate $r$ implies $p$-probability convergence with rate $r$ where $p(\delta)=1/\delta$.
    \item By the Borel-Cantelli lemma \citep[Thm. 4.3]{Billingsley}, $p$-probability convergence with rate $r$ implies almost sure convergence if $\sum_{t=0}^{\infty}\min\{1,p^{-1}(a/r(t))\}<\infty$ for all $a>0$. If $p(\delta)=1/\delta^c$ for some $c>0$, then $r(t)=o(1/t^c)$ is required. If $p(\delta)=\log(1/\delta)$, then only $r(t)=\order(1/t^c)$ for some $c>0$ is required.
\end{itemize}

The expected convergence rate of the function values for SGD in our setting is $r(t)=\order(1/t)$, so Markov's inequality is insufficient to prove almost sure convergence. However, almost sure convergence can be proved via other techniques and with weaker assumptions. In particular, \cite{Khaled} consider two weaker assumptions than the bounded variance property: the \textit{relaxed growth} property and the \textit{expected smoothness} property, which are equivalent in our setting. Almost sure convergence of SGD can be proved assuming the relaxed growth property \citep{Bertsekas,patel2020stopping}. While we extend the bounded variance property to higher moments in order to get high probability convergence, it is not clear how to do so for the relaxed growth property, so we leave this as a future research direction.

\section{Standard optimization results}
\label{subsec:optimizationResults}

We state a few useful definitions and results from optimization. More details can be found in \cite{Nesterov} and \cite{Bubeck}.  All our definitions are with respect to the Euclidean norm, denoted $\|\cdot\|$, and we let $f:\R^d \rightarrow \R$. We write $C^1(\R^d)$ or just $C^1$ to denote functionals on $\R^d$ that are continuously differentiable.

\begin{definition}[Lipschitz continuity]
We say $f:\R^d\rightarrow \R$ is $\rho$-Lipschitz continuous if $|f(x)-f(y)|\leq \rho\|x-y\|$ for all $x,y\in\R^d$. Likewise, $g:\R^d \rightarrow \R^m$ is $\rho$-Lipschitz continuous if $\|g(x)-g(y)\|\leq \rho\|x-y\|$ for all $x,y\in\R^d$.
\end{definition}
Note that averages of Lipschitz continuous functions are Lipschitz continuous with the same constant. Also, if $f\in C^1$, then $f$ is $\rho$-Lipschitz if and only if $\|\nabla f(x)\|\le \rho$ for all $x\in\R^d$.

\begin{definition}[L-smoothness]
We say $f$ is $L$-smooth if $f\in C^1$ and its gradient is $L$-Lipschitz continuous.
\end{definition}
If $f$ is $L$-smooth, then a standard result is
\begin{equation}
(\forall x,y \in\R^d)\;
f(y)\leq f(x)+\langle \nabla f(x),y-x\rangle+\frac{L}{2}\|x-y\|^2.
\end{equation}
Applying this result to $x-\frac{1}{L}\nabla f(x)$ and $x$ and rearranging terms gives
\begin{equation} \label{eq:LipschitzBound}
\|\nabla f(x)\|^2\leq 2L(f(x)-f\opt)
\end{equation}
and taking $x\in \argmin_{x'}\,f(x')$ gives
\begin{equation} \label{eq:LipschitzBound2}
f(y) - f^\star \le \frac{L}{2}\|y-x^\star\|^2.
\end{equation}
If $f$ is both $\mu$-PL and $L$-smooth, then combining the PL inequality and Eq.~\eqref{eq:LipschitzBound} gives
\[
f(x) - f^\star \le \frac{L}{\mu}\left( f(x) - f^\star \right)
\]
which implies $L \ge \mu$.

The following inequality is known as the \emph{quadratic growth} (QG) condition:
\begin{equation} \label{eq:QG}
\exists \mu > 0 \quad\text{s.t. } (\forall x\in\R^d) \;
\|x - \mathcal{P}_{X\opt}(x) \| \le \frac{2}{\mu}\left( f(x) - f\opt \right)
\end{equation}
where $\mathcal{P}_{X\opt}$ denotes the (set-valued) Euclidean projection onto the set of global minimizers, $X\opt=\argmin_x \, f(x)$. The PL inequality implies the QG condition~\citep{Karimi}.

\section{Proofs of results}

\subsection{PL inequality}\label{subsec:PL}

Strong convexity implies the PL inequality but the PL inequality only implies invexity, not convexity. Furthermore, unlike the convexity and strong convexity properties, the PL property is not preserved under nonnegative sums.  The PL inequality implies the quadratic growth property (Eq.~\eqref{eq:QG}), which can be used to bound how far away the nearest minimizer is in terms of the optimality gap (which can be bounded above by simply evaluating the cost function, if the cost function is nonnegative).

When working with constraints to a set $x\in X \subsetneq \R^d$, if $f$ is strongly convex, then so is $f$ added to the indicator function of $X$, and results for gradient descent methods easily extend to projected gradient descent. On the other hand, if $f$ is PL, then $f$ added to the indicator function is not KL (in particular, it is not even differentiable). This has real impacts on gradient descent algorithms, where gradient descent might converge but projected gradient descent does not.

For example, there is a smooth function, a mollified version of $f(x,y)=\left(a(x)_+^2-b\left(|y|-c\right)_+\right)_+$, such that the PL inequality is satisfied but projected gradient descent does \textit{not} converge to a minimizer;  we formalize this in the remark below.

\begin{remark}
\label{rmk:pl}
Consider $f(x,y)=\left(a(x)_+^2-b\left(|y|-c\right)_+\right)_+$ where $(\cdot)_+$ denotes $\max(\cdot,0)$ and $a,b>0,c\geq 0$. The minimum of $f$ is 0 and $X\opt=\{(x,y)~|~x\leq 0 \text{ or } |y|\geq \frac{a}{b}x^2+c\}$. If we use $\varphi(x)=\X_{B_1(0)}\cdot\exp\left(-1/(1-\|x\|^2)\right)/\Phi$ (where $\Phi$ is the normalization constant) to mollify $f$, then, for $\epsilon<c$, $f_{\epsilon}$ has PL constant $2a$ and smoothness constant $2a$. Consider the starting point $(d,0)$. For $a,b,c,d$ chosen appropriately, the distance from $(d,0)$ to its projection onto $X\opt$ is strictly less than $d$. Thus, if we let $\X$ be the ball centered at $(d,0)$ with radius equal to exactly that distance, then the constrained problem and the unconstrained problem have the same minimum. However, projected gradient flow, starting from $(d,0)$, ends up stuck at a non-minimizer: the point of $\X$ closest to $(0,0)$.
\end{remark}

In order to generalize gradient methods to projected gradient methods under PL-like assumptions, the proper generalization is that the function should satisfy a \emph{proximal} PL inequality~\citep{Karimi}. However, such an assumption is quite restrictive compared to the PL inequality. We leave the problem of convergence with just the PL inequality, via added noise or a Frank-Wolfe type construction, as a future research direction.

\subsection{Proof of Proposition \ref{prop:interrelated}} \label{proofs:interrelated}

\begin{proof}
The smoothness of $f$ easily follows from that of $f_i$ using the definition of $f$ in Eq.~\eqref{eq:ERM}, the definition of smoothness, and the triangle inequality.

The rest of the proof will hinge on bounding the SGD sequence $(x_t)$. Unrolling the SGD updates, we have $x_{T+1} = x_0 + \sum_{t=0}^T \eta_t g_t$. Define $\tilde{f}_t = \frac{1}{b_t}\sum_{i\in I_t}f_i$ so that $g(\cdot,b_t) = \nabla \tilde{f}_t$, and observe $\tilde{f}$ is $L$-smooth for the same reason that $f$ is $L$-smooth. Then,
\begin{align*}
 \| g(x,b_t) \|^2= \|\nabla \tilde{f}_t(x)\|^2 &\leq 2L(\tilde{f}_t(x)-\tilde{f}_t\opt) \\
 &\leq 2L\tilde{f}_t(x) \\
 &\leq 2L\frac{n}{b_t}f(x)
\end{align*}
where the first inequality follows from Eq.~\eqref{eq:LipschitzBound}, the second from nonnegativity, and the third from nonnegativity as well. For simplicity, assume a fixed minibatch size $b_t =b$.

Thus, by the triangle inequality,
$\|x_{T}-x_0\| \le \sqrt{2L\frac{n}{b}} \sum_{t=0}^{T-1}\eta_t\sqrt{f(x_t)}$.
For both Thm.~\ref{thm:fastconvergence} and Thm.~\ref{thm:slowconvergence}, the stepsize is $\eta_t = \order(1/t)$, and, by interpolation, $f(x_t)\le \log(1/\delta_t)\frac{C}{t^c}$ for each $t$ with probability $1-\delta_t$, where $C$ and $c>0$ are constants (e.g., $c=1$ for Thm.~\ref{thm:fastconvergence})\footnote{These theorems and stepsizes are actually in terms of $t+1$ due to the $0$-based indexing convention, so we now use $1$-based indexing to simplify presentation.}. Choose $\delta_t = \order(1/t^2)$, then it holds that
\[
\whp, \quad (\forall T \in\mathbb{N})\quad
\|x_{T}-x_0\| \lesssim C' \sum_{t=1}^T \frac{\log(t)}{t^{1+c/2}} \le C''
\]
for constants $C'$ and $C''$ since the series converges (verified by, e.g., the integral test) and with $\delta = \sum_{t=0}^{T-1} \delta_t \le \sum_{t=0}^\infty \delta_t<1 $
coming from the union bound and the fact $\sum_{t=1}^\infty t^{-2}=\pi^2/6< \infty$.

Let $E_k$ be the event that $\|x_T - x_0\|\le k$ for all $T$. Then the event $\bigcup_{k=1}^\infty E_k$ implies that $(x_t)_{t=0}^\infty$ is a bounded sequence.  Since $E_k$ is nested, by continuity of measure,
$\P(\bigcup_{k=1}^\infty E_k) = \lim_{k\rightarrow\infty} \P(E_k)$, and the argument from the previous paragraph shows that $\P(E_k)$ can be made arbitrarily close to $1$ by choosing a sufficiently large radius $k$. Hence $(x_t)_{t=0}^\infty$ is a bounded sequence with probability $1$.

Now that $(x_t)$ is bounded with probability $1$, we turn to $\rho$ and $M$. First, $\ell(x,s)$ is jointly continuous in $x$ and $s$, and since $(x,s)$ comes from a compact set, we get that $\ell$ attains its maximum value \citep[Thm.\ 1.68]{Hunter}, which we call $M$. The same is true of $\|\nabla_x\ell(x,s)\|$ and we call its maximum value $\rho$. Thus, since $\ell(x,s)$ is $C^1$ with respect to $x$, we have that it is $\rho$-Lipschitz with respect to $x$.
\end{proof}

\subsection{Proofs of results in Section \ref{sec:convergence}}
\label{subsec:convergence}

\paragraph{Proof of Proposition \ref{prop:harvey}}
\begin{proof}
By Jensen's inequality, $\exp\left(\E[\vh_t/\gamma_t\mid\F_{t-1}]\right)\leq\E\left[\exp(\vh_t/\gamma_t)\mid\F_{t-1}\right]\leq\exp(1)$, so $\E[\vh_t~|~\F_{t-1}]\leq \gamma_t$.  Applying the expectation and unraveling the recursion proves Eq.~\eqref{eq:harveyPropExpected}. For the high probability bound, we will use induction over $t\geq 0$ to prove that $\E[\exp(\lambda X_t)]\leq \exp(\lambda K_t)$ for all $\lambda \in (0,K_t^{-1}]$ (i.e. we will show that $X_t$ is $K_t$-sub-exponential). The base case follows by construction. Assume that $\E[\exp(\lambda X_t)]\leq \exp(\lambda K_t)$ for all $\lambda \in (0,K_t^{-1}]$. Let $\lambda \in (0,K_{t+1}^{-1}]$. We will use the following three inequalities that follow from the constraint on $(K_t)$:
\begin{align*}
    &i)~  2\lambda \leq \frac{1}{\gamma_t}, \hspace{.5cm} ii)~ \lambda(\alpha_t+\lambda\beta_t^2)\leq \frac{1}{K_t}, \hspace{.5cm} iii)~ K_{t+1}\geq K_t(\alpha_t +\lambda\beta_t^2)+\gamma_t.
\end{align*}
Then, using the law of total expectation and the Cauchy-Schwarz inequality, we get
\begin{align*}
    \E[\exp(\lambda X_{t+1})] &\leq \E[\exp(\lambda\alpha_t X_t+\lambda \wh_t+\lambda \vh_t)]\\
    &= \E[\exp(\lambda\alpha_t X_t)\E[\exp(\lambda\wh_t)\exp(\lambda\vh_t)~|~\F_{t-1}]]\\
    &\leq \E[\exp(\lambda\alpha_t X_t)\sqrt{\E[\exp(2\lambda\wh_t)~|~\F_{t-1}]\E[\exp(2\lambda\vh_t)~|~\F_{t-1}]}]\\
    &\overset{(i)}{\leq} \E[\exp(\lambda\alpha_t X_t+\lambda^2\beta_t^2X_t+\lambda \gamma_t)]\\
    &= \E[\exp(\lambda X_t(\alpha_t +\lambda\beta_t^2)+\lambda \gamma_t)]\\
    &\overset{(ii)}{\leq} \exp(\lambda K_t(\alpha_t +\lambda\beta_t^2)+\lambda \gamma_t)\\
    &\overset{(iii)}{\leq} \exp(\lambda K_{t+1}),
\end{align*}
thus completing the induction proof. In particular, $\E[\exp(X_t/K_t)]\leq e$, and so, by Markov's inequality, $\P(X_t\ge K_t\log(e/\delta))=\P(\exp(X_t/K_t)\ge e/\delta)\leq \delta$, proving Eq.~\eqref{eq:harveyPropHighProb}.
\end{proof}

\paragraph{Proof of convergence results}

Set
\begin{align*}
    &X_t = f(x_t)-f\opt\\
    &\alpha_t = 1-\mu\eta_t\\
    &\wh_t = \eta_t(1-L\eta_t)\langle \nabla f(x_t),e_t\rangle\\
    &\vh_t = \frac{L\eta_t^2}{2}\|e_t\|^2.
\end{align*}

Then there exists a universal constant $C_1$ such that
\begin{align*}
    \E[\exp(\lambda \wh_t)\mid\F_{t-1}] &= \E[\exp(\lambda\eta_t(1-L\eta_t)\|\nabla f(x_t)\|\langle \nabla f(x_t)/\|\nabla f(x_t)\|,e_t\rangle)\mid\F_{t-1}]\\
    &\leq \exp\left(\frac{\lambda^2\eta_t^2(1-L\eta_t)^2\|\nabla f(x_t)\|^2C_1\sigma^2}{2bd}\right)\\
    &\leq \exp\left(\frac{\lambda^2\eta_t^2(1-L\eta_t)^2 2LX_t C_1\sigma^2}{2bd}\right).
\end{align*}

Thus, we can set
\begin{align*}
    \beta_t^2 = \frac{C_1 2L\sigma^2\eta_t^2(1-L\eta_t)^2}{bd}.
\end{align*}

There also exists a universal constant $C_2$ such that $\|e_t\|$ is $C_2\sigma/\sqrt{b}$ sub-gaussian conditioned on $F_{t-1}$. Thus, we can set
\begin{align*}
    \gamma_t = \frac{C_2L\sigma^2\eta_t^2}{2b}.
\end{align*}

Set $K_0=X_0$. Then, we only need
\begin{align}
\label{eq:ktconstraint}
    K_{t+1}^2\ge (\alpha_tK_t+2\gamma_t)K_{t+1}+\beta_t^2K_t.
\end{align}
For Theorem \ref{thm:fastconvergence}, we have $\eta_t=\theta_t$. For $t=0,\ldots,\tau-1$,
\begin{align*}
    &\alpha_t=1-\frac{\mu}{L}\\
    &\beta_t^2=0\\
    &\gamma_t = \frac{C_2\sigma^2}{2Lb}.
\end{align*}
Thus, Eq.~\eqref{eq:ktconstraint} becomes
\begin{align*}
    K_{t+1} \ge \left(1-\frac{\mu}{L}\right)K_t+\frac{C_2\sigma^2}{Lb}
\end{align*}
and so, for $t=1,\ldots,\tau$, we can set
\begin{align*}
    K_t=\left(1-\frac{\mu}{L}\right)^t X_0+\frac{C_2\sigma^2}{\mu b}.
\end{align*}

For $t=\tau,\tau+1,\ldots$,
\begin{align*}
    &\alpha_t= 1-\frac{2t+1}{(t+1)^2}=\frac{t^2}{(t+1)^2}\\
    &\beta_t^2 = \frac{C_12L\sigma^2\eta_t^2(1-L\eta_t)^2}{bd}\leq \frac{C_12L\sigma^2\eta_t^2}{bd}\leq \frac{C_12L\sigma^2(2t+1)^2}{bd\mu^2(t+1)^4}\leq \frac{18C_1L\sigma^2t^2}{bd\mu^2(t+1)^4}\\
    &\gamma_t = \frac{C_2 L\sigma^2(2t+1)^2}{2b\mu^2(t+1)^4}\leq \frac{2C_2L\sigma^2}{b\mu^2(t+1)^2}.
\end{align*}
Thus, the following equation implies Eq.~\eqref{eq:ktconstraint}:
\begin{align*}
    K_{t+1}^2 &\ge \left(\frac{t^2}{(t+1)^2}K_t+\frac{4C_2L\sigma^2}{b\mu^2(t+1)^2}\right)K_{t+1}+\frac{18C_1L\sigma^2 t^2}{bd\mu^2 (t+1)^4}K_t.
\end{align*}
Multiplying by $(t+1)^4$, this becomes
\begin{align*}
    L_{t+1}^2 \ge \left(L_t+\frac{4C_2L\sigma^2}{b\mu^2}\right)L_{t+1}+\frac{18C_1L\sigma}{bd\mu^2}L_t
\end{align*}
where $L_t=t^2K_t$. This is satisfied, for example, by
\begin{align*}
    L_{t+1} &= L_t + \frac{4C_2L\sigma^2}{b\mu^2}+\frac{18C_1L\sigma^2}{bd\mu^2}=\frac{(18C_1+4C_2d)L\sigma^2}{bd\mu^2}.
\end{align*}
Unraveling this, we get
\begin{align*}
    L_t &= L_{\tau}+\frac{(18C_1+4C_2d)L\sigma^2}{bd\mu^2}(t-\tau)
\end{align*}
and so
\begin{align*}
    K_t &= \frac{K_{\tau}\tau^2}{t^2}+\frac{(18C_1+4C_2d)L\sigma^2(t-\tau)}{bd\mu^2t^2}\\
    &= \frac{\left(\left(1-\frac{\mu}{L}\right)^{\tau}X_0+\frac{C_2\sigma^2}{\mu b}\right)\tau^2}{t^2}+\frac{(18C_1+4C_2d)L\sigma^2(t-\tau)}{bd\mu^2t^2}.
\end{align*}

Thus, we have proved the following explicit version of Theorem \ref{thm:fastconvergence}:

\begin{theorem}
\label{thm:appendixconv}
Assume $f$ is $L$-smooth and $\mu$-PL and that $\nabla f(x)-g(x,1)$ is $\sigma/\sqrt{d}$-sub-gaussian. Then, SGD with $(\eta_t)=(\theta_t)$ and $b_t=b$ constructs a sequence $(x_t)$ such that, for $t=0,...,\tau$,
\begin{align}
    \E[f(x_t)-f\opt] &\leq \left(1-\frac{\mu}{L}\right)^t(f(x_0)-f\opt)+\frac{C_2\sigma^2}{\mu b}
\end{align}
and, $\whp$ for all $\delta\in(0,1/e)$,
\begin{align}
    f(x_t)-f\opt &\leq \left(\left(1-\frac{\mu}{L}\right)^t(f(x_0)-f\opt)+\frac{C_2\sigma^2}{\mu b}\right)\log(e/\delta) \, .
\end{align}
For $t>\tau$, one has that
\begin{align}
    \E[f(x_t)-f\opt] &\leq \frac{\left(\left(1-\frac{\mu}{L}\right)^{\tau}(f(x_0)-f\opt)+\frac{C_2\sigma^2}{\mu b}\right)\tau^2}{t^2}+\frac{2C_2L\sigma^2(t-\tau)}{\mu^2 bt^2}
\end{align}
and, $\whp$ for all $\delta\in(0,1/e)$,
\begin{align}
\begin{split}
    &f(x_t)-f\opt \\
    &\leq\left(\frac{\left(\left(1-\frac{\mu}{L}\right)^{\tau}(f(x_0)-f\opt)+\frac{C_2\sigma^2}{\mu b}\right)\tau^2}{t^2}+\frac{(18C_1+4C_2 d)L\sigma^2(t-\tau)}{\mu^2 bdt^2}\right)\log(e/\delta).
\end{split}
\end{align}
$C_1$ and $C_2$ are universal constants.
\end{theorem}

For Theorem \ref{thm:slowconvergence}, we have $\eta_t=\frac{2c}{t+2}$. So,
\begin{align*}
    &\alpha_t = 1-\frac{2\mu c}{t+2}\leq e^{\frac{-2\mu c}{t+2}}\leq\left(\frac{t+1}{t+2}\right)^{\mu c}\\
    &\beta_t^2 = \frac{C_1 2L\sigma^2\eta_t^2(1-L\eta_t)^2}{bd}\leq \frac{C_12L\sigma^2\eta_t^2}{bd}=\frac{8C_1L\sigma^2c^2}{bd(t+2)^2}\leq \frac{16C_1L\sigma^2c^2(t+1)^{\mu c}}{bd(t+2)^{2+\mu c}}\\
    &\gamma_t = \frac{C_2L\sigma^2\eta_t^2}{2b}=\frac{2C_2L\sigma^2c^2}{b(t+2)^2}.
\end{align*}
Thus, the following equation implies Eq.~\eqref{eq:ktconstraint}:
\begin{align*}
    K_{t+1}^2 &\ge \left(\frac{(t+1)^{\mu c}}{(t+2)^{\mu c}}K_t+\frac{4C_2L\sigma^2c^2}{b(t+2)^2}\right)K_{t+1}+\frac{16C_1L\sigma^2 c^2(t+1)^{\mu c}}{bd (t+2)^{2+\mu c}}K_t.
\end{align*}
Multiplying by $(t+2)^{2\mu c}$, this becomes
\begin{align*}
    L_{t+1}^2 &\ge \left(L_t+\frac{4C_2L\sigma^2c^2}{b(t+2)^{2-\mu c}}\right)L_{t+1}+\frac{16C_1L\sigma^2 c^2}{bd (t+2)^{2-\mu c}}L_t.
\end{align*}
where $L_t=(t+1)^{\mu c}K_t$. This is satisfied, for example, by
\begin{align*}
    L_{t+1} &= L_t + \frac{4C_2L\sigma^2c^2}{b(t+2)^{2-\mu c}}+\frac{16C_1L\sigma^2 c^2}{bd (t+2)^{2-\mu c}}=L_t+\frac{(16C_1+4C_2 d)L\sigma^2 c^2}{bd(t+2)^{2-\mu c}}.
\end{align*}
Unraveling this, we get
\begin{align*}
    L_t &= \frac{(16C_1+4C_2 d)L\sigma^2 c^2}{bd}\sum_{i=0}^{t-1} \frac{1}{(i+2)^{2-\mu c}}\\
    &\leq \frac{(16C_1+4C_2 d)L\sigma^2 c^2}{bd(1-\mu c)}.
\end{align*}
Dividing by $t^{\mu c}$ then proves Theorem \ref{thm:slowconvergence}.

\subsection{Proofs of results in Section \ref{sec:Stability}}
\label{subsec:stability}

Define $\delta_t \defeq \|x_t-x_t'\|$. Note that $\P(I_t=I_t')=\left(1-\frac{1}{n}\right)\cdots\left(1-\frac{1}{n-b_t+1}\right)=1-\frac{b_t}{n}$. Following \cite{Hardt}, a uniform stability bound is sought by  bounding the growth of $\delta_t$; if  $\ell(\cdot,s)$ is $\rho$-Lipschitz [cf. Proposition~\ref{prop:interrelated}], this implies that $\sup_s |\ell(x_t,s)-\ell(x_t',s)|\leq \rho \delta_t$. The next lemma is from \cite{Hardt}.

\vspace{.1cm}

\begin{lemma}
\label{lem:growth}
\citep[Lem. 2.5]{Hardt} Assume $\ell(\cdot,s)$ is $\rho$-Lipschitz and $L$-smooth for all $s$. Then,
\begin{align}
    \delta_{t+1} &\leq \begin{cases}(1+\eta_t L)\delta_t &\text{ if }I_t=I_t'\\\delta_t+2\eta_t\rho &\text{ o.w.}\end{cases}.
\end{align}
\end{lemma}

\paragraph{Proof of Theorem \ref{thm:stability}}
\begin{proof}
Define $a_t=0$ if $I_t=I_t'$ and $a_t=1$ otherwise. Hence, $a_t\sim B(1,p)$ independently with $p=b/n$. By Lemma \ref{lem:growth},
\begin{align*}
    \delta_{t+1} &\leq (1+(1-a_t)\eta_t L)\delta_t+2\eta_t \rho a_t\\
    &\leq (1+\eta_t L)\delta_t+2\eta_t \rho a_t.
\end{align*}
Thus,
\begin{align*}
    \delta_T &\leq 2\rho\sum_{t=0}^{T-1}\prod_{i=t+1}^{T-1}(1+\eta_i L)\eta_t a_t\\
    &\leq 2\rho\sum_{t=0}^{T-1}\exp\left(L\sum_{i=t+1}^{T-1}\eta_i\right)\eta_t a_t\\
    &\leq 2\rho c\sum_{t=0}^{T-1}\exp\left(Lc\log\left(\frac{T}{t+1}\right)\right)\frac{a_t}{t+1}\\
    &\leq 2\rho c T^{Lc}\sum_{t=0}^{T-1}\frac{a_t}{(t+1)^{Lc+1}}.
\end{align*}
Note that
\begin{align*}
    \sumt \frac{1}{(t+1)^{Lc+1}} &\leq 1+\int_0^{T-1}\frac{dt}{(t+1)^{Lc+1}}\\
    &= 1+\frac{1}{Lc}-\frac{1}{LcT^{Lc}}\\
    &\leq 1+\frac{1}{Lc}.
\end{align*}
The result follows by taking expectation and using $\sup_s |\ell(x_t,s)-\ell(x_t',s)|\leq \rho\delta_t$.
\end{proof}

\begin{theorem}
\label{thm:elisseeff}
\citep[Thm. 12]{Elisseeff}
Assume $\ell(x,s)\in[0,M]$ for all $x$ and $s$. Let $\delta\in(0,1)$. If SGD with $T$ steps is $\alpha$-uniformly stable in expectation then, $\whp$ over $S$ and $(i_t)$ for all SGD,
\begin{align}
    F(x_T)-f(x_T) &\leq \sqrt{\frac{6Mn\alpha+M^2}{2n\delta}}.
\end{align}
\end{theorem}

Note that Thm. 12 of \cite{Elisseeff} is actually for pointwise hypothesis stability in expectation, but uniform stability in expectation implies pointwise hypothesis stability in expectation. Furthermore, in the form we have it, Theorem \ref{thm:elisseeff} can be proved using Lem. 9 of \cite{Bousquet} and Chebyshev's inequality.

\end{document}


TITLE: Finding the Smith Normal Form of an integer matrix
QUESTION [1 upvotes]: I want to put the following integer matrix into Smith Normal Form:
$$\begin{pmatrix} -9 & 6 \\ 5 & -2 \\ 6 & 3 \end{pmatrix}$$
I have done this and found the answer to be $$\begin{pmatrix} 1 & 0 \\ 0 & 3 \\ 0 & 0\end{pmatrix}$$
Could someone verify whether this is correct or not?

REPLY [1 votes]: $$\begin{pmatrix}-9 & 6 \\ 5 & -2 \\ 6 & 3\end{pmatrix}\to\begin{pmatrix}2 & -5 \\ 6 & -9 \\ 3 & 6\end{pmatrix}\to \begin{pmatrix}2 & -5 \\ 0 & 6 \\ 1 & 11\end{pmatrix}\to\begin{pmatrix}1 & 11 \\ 0 & 6 \\ 0 & -27\end{pmatrix}\to\begin{pmatrix}1 & 0 \\ 0 & 6 \\ 0 & 3\end{pmatrix}\to\begin{pmatrix}1 & 0 \\ 0 & 3 \\ 0 & 0\end{pmatrix}$$
Rearrange rows and columns, reduce, rearrange again, reduce (twice), rearrange and reduce.


text:
In this section, we will apply the NEBV estimator to construct $100(1-\gamma)\%$ empirical Bayes intervals for the mean parameters assuming the following model:
\begin{eqnarray}\label{eqn:model:mean}
 \begin{cases}
\,\,\,X_i|\theta_i,\sigma_i^2&\overset{ind}\sim\,\,\,\, N(\theta_i,\sigma_i^2),\\
\,\,\,S_i^2|\sigma_i^2         &\overset{ind}\sim\,\,\,\,\sigma_i^2\frac{\chi_k^2}{k},\\
\,\,\,\quad\theta_i             &\overset{iid}\sim\,\,\,\,N(\mu,\tau^2),\\
\,\,\, \quad\sigma_i^2        &\overset{iid}\sim\,\,\,\, \pi(\sigma_i^2).
\end{cases}
\end{eqnarray}

Note that the posterior distribution of $\theta_i$ is
\begin{equation*}
	\,\,\,\theta_i \big|X_i, \sigma_i^2\,\,\,\sim  \,\,\,\,N\left(\,\frac{\tau^2}{\tau^2+\sigma_i^2}X_i+\frac{\sigma_i^2}{\tau^2+\sigma_i^2}\mu,\,\,\, \frac{\tau^2 \sigma_i^2}{\tau^2+\sigma_i^2} \right).
\end{equation*}

When $\sigma_i^2$'s are known, the hyper-parameters $(\mu,\tau^2)$ can be estimated using the method of moments as
\begin{equation*}
                \hat \mu =\bar X \,\textrm{ or }\, \sum_{i=1}^p\frac{X_i/\sigma_i^2}{\sum_{i=1}^p 1/\sigma_i^2}, \,\,\textrm{ and }\,\,	\hat\tau^2=\max\left( \frac{1}{p}\sum_{i=1}^p(X_i-\hat\mu)^2-\frac{1}{p} \sum_{i=1}^p\sigma_i^2,\,\,\, \tau_0^2\right),
\end{equation*}
where $\tau_0^2$, known as the truncated value (\citet{Hwang:Qiu:Zhao:2009}), is defined as, with $\gamma^*\ge\Pr_{\tau^2}(\hat\tau^2=0)$,
\[
    \tau_0^2=\frac{1}{p}\left(z_{\gamma^*}^2+z_{\gamma^*}\sqrt{z_{\gamma^*}^2+2\sum_{i=1}^p\sigma_i^2}\,\right).
\]

For unknown $\sigma_i^2$'s, we replace them by various estimates $\hat{\sigma}^2_{i}$ to get the following $100(1-\gamma)\%$ empirical Bayesian confidence interval
\begin{equation}\label{eq:CI}
	\left(\hat \theta_i - z_{\gamma/2}\sqrt{\frac{\hat\sigma_{i}^2 \hat\tau^2}{\hat\sigma_{i}^2+\hat\tau^2}}\,,\,\,\,\hat \theta_i + z_{\gamma/2}\sqrt{\frac{\hat\sigma_{i}^2 \hat\tau^2}{\hat\sigma_{i}^2+\hat\tau^2}}\,\,\right),
\end{equation}
where
\begin{equation*}
  \hat\theta_i =\frac{\hat\tau^2}{\hat\sigma_{i}^2+\hat\tau^2}X_i+ \frac{\hat\sigma_{i}^2}{\hat\sigma_{i}^2+\hat\tau^2}\hat\mu.
\end{equation*}

In many applications, we might be only interested in $\theta_{(p)}$ for which the corresponding $T$-statistic attains the maximum in magnitude among all the $T$-statistic, $T_i=\frac{X_i}{S_i}$. We then add the parentheses adjustment according to \citet{Hwang:Zhao:2013}.

In the comparison, we applied $\hat{\sigma}_{NEBV}^2$, $\hat{\sigma}_{ELJS}^2$ (\citet{Cui:Hwang:Qiu:Blades:Churchill:2005}) and $\hat{\sigma}_{TW}^2$ (\citet{Tong:Wang:2007}). We wanted to point out that these methods are not designed for estimating the variance post the selection. Actually, to the best of our knowledge, there is no methods designed for estimating the variances post the selection in the exising literature.

\subsection{Independent Case}
We generated $(X_i,S_i^2)$'s according to Model (\ref{eqn:model:mean}) with different choices of the prior for $\sigma_i^2$'s, including
 inverse Gamma prior ($IG(\alpha,\beta)$), mixture of Gamma/inverse Gamma prior, and log-normal prior ($LN(\mu_v,\tau_v^2)$).
The dimension for the numerical study is 3,000.
After obtaining the statistic $(X_i,S_i^2)$, we made the inference for the parameters $(\theta_{(p)},\sigma_{(p)}^2)$ corresponding to the observation with the largest $t$-statistic in magnitude. These steps were replicated 1,000 times.

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Risk_p_2_a_723_b_246.pdf}&
		\includegraphics[width=60mm]{Risk_p_2_a_805_b_472.pdf} \\
		\includegraphics[width=60mm]{Risk_p_1_a_3000_b_1000.pdf} &
		\includegraphics[width=60mm]{Risk_p_1_a_2444_b_722.pdf}
	\end{tabular}
	\caption{Plots of $\ln\textrm{(Risk)}$ for estimating $\sigma^2_{(i)}$ against $\frac{\tau^2}{\tau^2+E(\sigma^2_i)}$ under log-normal (the first row) and inverse gamma prior (the second row) for $\sigma^2_i$.}
	\label{fig:indep_risk}
\end{figure}

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Risk_MixIG.pdf}&
		\includegraphics[width=60mm]{Risk_MixGam.pdf}
	\end{tabular}
	\caption{Plots of $\ln\textrm{(Risk)}$ for estimating $\sigma^2_{(i)}$:  $0.5*IG(6,1)+0.5*IG(6,4)$ for the left panel and  $0.5*Gamma(4,1/16)+0.5*Gamma(6,1/8)$ for the right panel.}\label{fig:risk_mix}
\end{figure}

We first compared the risk of different estimators for estimating the variance $\sigma_{(p)}^2$. The logarithm of the risks based on the loss function (\ref{eq:lossselect}) are reported in Figure \ref{fig:indep_risk}. It is clearly seen that the proposed NEBV has the smallest risk for all the parameter settings.

Next, we considered the confidence interval for $\theta_{(p)}$. In addition to the three intervals by replacing $\sigma_i^2$ with one of the three variance estimators above, we included the second-order corrected (SOC) interval (\citet{Hwang:Zhao:2013}) and Bonferroni's corrected $T$-interval: $X_{(p)}\pm\, t_{k,\frac{\gamma}{2p}}\,S_{(p)}$. The SOC intervals are designed for the log-normal model where $\log S_i^2|\log \sigma_i^2\sim N(\log\sigma_i^2+m_k,\sigma_k^2)$ and $\log\sigma_i^2\sim N(\mu_v,\tau_v^2)$. We set the desired coverage probability as 95\%, namely $\gamma=0.05$.

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Coverage_p_2_a_723_b_246.pdf}&
        \includegraphics[width=60mm]{Coverage_p_2_a_805_b_472.pdf} \\
		\includegraphics[width=60mm]{Length_p_2_a_723_b_246.pdf}&
		\includegraphics[width=60mm]{Length_p_2_a_805_b_472.pdf}
	\end{tabular}
	\caption{Plots of the coverage probabilities (the first row)  and the average length ratios (the second row)  of various intervals for $\theta_{(p)}$ against $\frac{\tau^2}{\tau^2+E(\sigma^2_i)}$ under log-normal prior  with error bar: $0.95-1.96\sqrt{(0.95)(0.05)/1000}$.}
\label{fig:indep_lognormal}
\end{figure}

\begin{figure}
	\begin{tabular}{c c}

          \includegraphics[width=60mm]{Coverage_p_1_a_3000_b_1000.pdf} &
          \includegraphics[width=60mm]{Coverage_p_1_a_2444_b_722.pdf}\\
          \includegraphics[width=60mm]{Length_p_1_a_3000_b_1000.pdf} &
          \includegraphics[width=60mm]{Length_p_1_a_2444_b_722.pdf}
        \end{tabular}
	\caption{Plots of the coverage probabilities (the first row) and the average length ratio (the second row) for various intervals of $\theta_{(p)}$ for the inverse Gamma prior
	}
	\label{fig:indep_gamma}
\end{figure}

We reported the simulated results corresponding to the log-normal prior in Figure \ref{fig:indep_lognormal}. The coverage probabilities were reported in the first row. It is seen that the intervals based on NEBV, the Bonferroni's corrected $t$-interval, and SOC intervals all guarantee good coverage probabilities for all the hyper-parameters. The intervals based on the other two estimators, ELJS and TW, have substantially low coverage probabilities. For instance, for the inverse Gamma prior with $\alpha=3$ and $\beta=1$, the coverage probability of intervals based on TW can be as low as 35\%. This is expected because these method are not designed for the setting post the selection. We also reported the average ratio of the length of all the confidence intervals to the length of $t$-interval based on the Bonferroni's correction in the second row of the figure. It is seen that among the three intervals which guarantee a desired coverage probability, the one based on NEBV is the shortest one.

In Figure \ref{fig:indep_gamma}, we reported the coverage probability and average length ratio of the various approaches for the inverse Gamma prior. The intervals based on NEBV and the Bonferroni's corrected $t$-interval still have good coverage probability. However, the coverage probability of SOC interval dropped below the nominal level. SOC intervals are designed to work well for the log-normal model and does not work well when this assumption does not hold.

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Coverage_MixIG.pdf} &
		\includegraphics[width=60mm]{Coverage_MixGam.pdf} \\
		\includegraphics[width=60mm]{Length_MixIG.pdf} &
		\includegraphics[width=60mm]{Length_MixGam.pdf}
	\end{tabular}
	\caption{Plots of the coverage probabilities and the length ratio of confidence intervals for $\theta_{(p)}$ under Inverse Gamma and Gamma mixture prior for $\sigma^2_i$ :  $0.5*IG(6,1)+0.5*IG(6,4)$ for left panels and  $0.5*Gamma(4,1/16)+0.5*Gamma(6,1/8)$ for right panels.}
	\label{fig:ind_mix}
\end{figure}

In Figure \ref{fig:ind_mix}, the variances $\sigma_i^2$'s were generated from a mixture Gamma/inverse Gamma distribution. Although SOC intervals can be shorter than the NEBV-based intervals, but the coverage probabilities of SOC intervals were much lower than 95\% for some parameter settings. The NEBV-based intervals are robust to the prior of $\sigma_i^2$'s and still work well. We have done extensive simulations and all of which demonstrated the patterns similar to the representative cases reported here.

\subsection{Dependent Case}\label{sec:dependent}
The previous estimators of the variances and the confidence intervals (\ref{eq:CI}) are based on the independence assumption. In this section, we compared their numerical performance when the independence assumption is violated. We generated the data from the multivariate normal distribution. Namely, for $i=1,2,\ldots, n$,
\begin{align}\label{eq:simulmodel2}
\begin{split}
\boldsymbol{X_{i,\cdot}}=(x_{i1},\ldots, x_{ip}) \,|\, \boldsymbol{\theta},\, \Sigma\,\,\,\, &\overset{}{\sim}\,\,\,\, N(\boldsymbol{\theta},\, \Sigma), \qquad \theta_i \,\,\,\, \overset{iid}{\sim}\,\,\,\, N(\mu, \tau^2).

\end{split}
\end{align}
Here, the variance-covariance matrix $\Sigma$ is given as $D^{1/2}RD^{1/2}$ where $D=diag(\sigma_1^2,\ldots,\sigma_p^2)$ with $\sigma_i^2\overset{iid}\sim IG(\alpha,\beta)$. Four types of correlation matrix $R$ were considered:
\begin{enumerate}
\item[(i)] Banded Correlation matrix ($\rho_{ij}=0$ when $|j-i|>d$ and $\rho_{ij}=\rho$ when $0<|j-i|\le d$),
\item[(ii)] AR Correlation matrix ($\rho_{ij}=\rho^{|i-j|}$),
\item[(iii)] Sparse Correlation matrix, and
\item[(iv)] Correlation matrix from real data.
\end{enumerate}
In the banded case, $(d, \rho)$ were set as $(5,0.1), (4,0.2), (3,0.3), (2,0.4)$, and $(1,0.5)$ respectively. For the AR case, $\rho$ was chosen as 0.5, 0.75, 0.9, 0.95, and 0.99 respectively. To get a sparse correlation matrix, we started with a $p\times p$ matrix $A$ by setting 1.5\% of randomly chosen entries as numbers from $uniform(a,b)$ and the rest as zero. Let $R=B^{-\frac{1}{2}}A'AB^{-\frac{1}{2}}$ where $B$ is a diagonal matrix with the diagonal-line elements being the diagonal line of $A'A$. The $R$ is a sparse correlation matrix and about 77\% of the entries are zero under this setting. We considered the following choices of $(a,b)$: $(-1,1),\,\,(0,1),\,\,(0.5,1),\,\,(0.9,1),\,\, $ and  $(0.9,0.9)$.

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Coverage_AR_rho_75_df_2.pdf} &
		\includegraphics[width=60mm]{Coverage_AR_rho_95_df_2.pdf}\\
		\includegraphics[width=60mm]{Coverage_Banded_rho_30_df_2_k_3.pdf} &
		\includegraphics[width=60mm]{Coverage_Banded_rho_40_df_2_k_2.pdf}\\
		\includegraphics[width=60mm]{Coverage_Sparse_a_10_b_20_df_2_spar_77.pdf} &
		\includegraphics[width=60mm]{Coverage_Sparse_a_5_b_5_df_2_spar_77.pdf}
	\end{tabular}
	\caption{Plots of the coverage probabilities of four intervlas with AR, Banded and Sparse correlation matrices.}\label{fig:dependent}
\end{figure}

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Length_AR_rho_75_df_2.pdf} &
		\includegraphics[width=60mm]{Length_AR_rho_95_df_2.pdf}\\
		\includegraphics[width=60mm]{Length_Banded_rho_30_df_2_k_3.pdf} &
		\includegraphics[width=60mm]{Length_Banded_rho_40_df_2_k_2.pdf}\\
		\includegraphics[width=60mm]{Length_Sparse_a_10_b_20_df_2_spar_77.pdf} &
		\includegraphics[width=60mm]{Length_Sparse_a_5_b_5_df_2_spar_77.pdf}
	\end{tabular}
	\caption{Plots of the average length ratio of four intervlas with AR, Banded and Sparse correlation matrices.}\label{fig:dep:length}
\end{figure}

We have done extensive simulation studies and only reported representative results in Figure \ref{fig:dependent}. Under all these dependence cases, the coverage probability of the NEBV-based intervals for $\theta_{(p)}$ guarantee the desired coverage probability. On the other hand, the intervals based on ELJS and TW fail to achieve the satisfactory coverage probability.
In Figure \ref{fig:dep:length}, we reported the average length ratio of the confidence intervals.

\begin{figure}
	\begin{tabular}{c c}
		\includegraphics[width=60mm]{Coverage_Colon.pdf} &
		\includegraphics[width=60mm]{Coverage_Leuk.pdf}\\
		\includegraphics[width=60mm]{Length_Colon.pdf} &
		\includegraphics[width=60mm]{Length_Leuk.pdf}\\

	\end{tabular}
	\caption{Plots of coverage probabilities, length ratio of confidence intervals for $\theta_{(p)}$ based on the Colon Cancer and Leukemia Data.}
	\label{fig:covlength_Colon}
\end{figure}

In the last case, we considered a correlation matrix estimated from two microarray experiments : Leukemia (\citet{Golub:Slonim:Tamayo:1999}) and Colon Cancer (\citet{Alon:Barkai:Notterman:1999}). We downloaded the data from \url{http://www.bioconductor.org/packages/release/data/experiment/html/golubEsets.html} and \url{http://genomics-pubs.princeton.edu/oncology/affydata/index.html}. The leukemia data includes the expressions of genes extracted from patients with two types of leukemia($p=7129$): Acute Lymphoblastic Leukemia(ALL, 47 patients) and Acute Myeloid Leukemia (AML, 25 patients). The colon cancer data contains gene expressions of 2000 genes for 22 normal people and 40 patients. Since leukemia data and colon cancer data have 72 and 62 observations with 7128 and 2000 variables, respectively. We randomly selected 70 and 60 variables respectively from each dataset and calculated the sample covariance matrices which were treated as $\Sigma$ in the simulations.

The results of coverage probabilities, and length of the intervals for $\theta_{(p)}$ are shown in Figure \ref{fig:covlength_Colon}. What we have observed is similar to that for the other settings. The intervals based on NEBV guarantees a desired coverage probability; while the coverage probabilities of the intervals based on ELJS and TW could be much lower than the nominal level.

\begin{table}
  \begin{tabular}{|ccc||ccc|}
    \hline
    Method     & Mean   &Median &Method & Mean & Median\\
    \hline
    \textit{Leukemia}      &    & & \textit{Colon} & & \\
    $S^2$         & 0.211  & 0.209 & $S^2$& 0.238  & 0.232 \\
    \textit{ELJS}  & 0.166  & 0.164 & \textit{ELJS} & 0.176  & 0.170 \\
    \textit{TW}  & 0.169  & 0.167 & \textit{TW}  & 0.178  & 0.172\\
    \textit{NEBV} & 0.149  & 0.146 & \textit{NEBV}& 0.156  & 0.152 \\
    \hline
  \end{tabular}
  \caption{Percentage of discordant genes for various intervals.}\label{tab:leukemia}
\end{table}

\section{Real data Analysis}\label{sec:realdata}

According to a 2016 poll of 1,500 scientists reported in the journal {\it Nature}, 70\% of them had failed to reproduce at least one other scientist's experiment (\citet{Baker:2016}). {\it Replication crisis}, which has long-standing roots, catches growing awareness in the last decades. A natural question is in what degree, the statistical conclusion of a certain method based on the data agrees with the conclusion when provided with different/new data set. In this section, we used the above two microarray data to examine the performance of various statistical methods.

For the Leukemia data set, we first randomly split the subjects into two subgroups such that both subgroups contain the similar numbers of subjects from ALL and AML groups. We then constructed confidence intervals based on one sub-group. We declared a gene to be significant if the corresponding interval does not enclose zero. We did the same for the other sub-group. We call a gene {\it discordant} if the interval based on the first subgroup does not enclose zero while the interval based on the second subgroup encloses zero. Namely, a gene is {\it discordant} implies that the conclusion that it is significant based on the first subgroup cannot be replicated based on the second one. We repeated these steps 1,000 times to calculate the average proportions of discordant genes. We did the same thing for the colon cancer data and reported the results in Table \ref{tab:leukemia} and Figure \ref{fig:boxplot}.

\begin{figure}
	\begin{tabular}{|c|c|}\hline
		\includegraphics[width=55mm,height=60mm]{HistLeuk.pdf} &
		\includegraphics[width=55mm,height=60mm]{HistColon.pdf} \\
		\hline
	\end{tabular}
	\caption{Boxplots of the percentage of discordant genes for leukemia and colon cancer data.}
	\label{fig:boxplot}
\end{figure}

In Figure \ref{fig:boxplot}, we plotted the box-plot based on four interval constructions: intervals based on the sample variance, ELJS, TW, and NEBV. In Table \ref{tab:leukemia}, we summarized the summary statistic based on 1,000 replications. It is seen that the intervals based on NEBV commit the smallest number of {\it discordance} and the difference is statistically significant.

In summary, as shown in both simulations and real data analysis, the performance of the intervals based on NEBV is better than its competitors and is thus recommended.

config: arxiv


TITLE: How does the Cantor set contain more points than the boundary points of those deleted open intervals?
QUESTION [3 upvotes]: I'm trying to get an intuitive grasp of the Cantor set.
Let $$C_0 = [0,1]$$ and $$C_n = \frac{C_{n-1}}{3} \cup (\frac{2}{3} + \frac{C_{n-1}}{3}).$$ Then we call $$\mathcal C = \bigcap_{n \in \mathbb N}$$ the Cantor set.
Here's something that my intuition cannot resolve: The Cantor set is constructed in countably many steps where in each step a finite number of open intervals is deleted.
Now, it's quite clear that all the end points of those deleted intervals will remain in the set. However, these are just countably many since we perform a countable number of steps and each step only adds a finite number of end points.
Now, here's where my intuition escapes me: How is it that any point that is not an end point of one of these deleted open intervals remains in the set as well? Wouldn't any non-end-point be cut out eventually? How can I imagine this process properly?

REPLY [1 votes]: "Each step only adds a finite number of end points". No, each of those endpoints was already there in $C_0$, and each is there in $C_1$, and each is there in $C_2$, and so on.
You are thinking of what each step "adds" when instead you should be thinking of what each step "removes". What each step does is to remove all the points in a certain finite union of open intervals. Step 1 removes the points in $(1/3,2/3)$; step 2 removes the points in $(1/9,2/9) \cup (7/9,8/9)$; and so on. It's the points that are never removed that you should be focussing your attention on.
The point $0$ is never removed, so it is in the Cantor set. Each point which ends up some time down the line being an endpoint of one of the removed intervals is itself never removed, and so each such point in the Cantor set. And there are many other points --- such as $\frac{1}{4}$ mentioned by @user3329719, or $0.202002000200002...$ mentioned by @user --- which are not endpoints and which nonetheless are never removed so each is in the Cantor set.


TITLE: Given a NxN grid and an integer calculate the row and column
QUESTION [0 upvotes]: I know this is most likely elementary school math, but I am having hard time finding a good algorithm. Given a NxN grid and a number which represents the nth tile counting left to right like a typewriter, where at the end of the row you return to column 1... e.g. given a grid of 20x20 and the number is 21, 21 would be located in column 1 row 2. Another example 36 in a 20x20 grid would be column 16 row 2. And so on.
What is the formula I can use to find the row and column.
I tried using modulo to find the column and division for the row, but it just does not work yet.

REPLY [1 votes]: Let's do a similar problem, where you have indices starting from $0$. In other words, you have a grid with $m$ rows and $n$ columns and you need to find the index of some number.
For simplicity let's look at $n=10$ first, letting $m=30$. Then, entry $0$ maps to $(0,0)$, and we follow with similar relationships below:

$0 \mapsto 0,0$
$1 \mapsto 0,1$
\ldots
$9 \mapsto 0,9$
$10 \mapsto 1,0$
\ldots
$19 \mapsto 1,9$

This may be enough to get the pattern. The entry number $N$ maps the second coordinate as $N \pmod{10}$ and the first coordinate is $\lfloor N \div 10 \rfloor$ which we will for convenience denote $N//10$.
In full generality, we map $N \mapsto N//n,N \pmod n$.

The problem you posed is slightly more tricky, since indexing is $1$-based, so given $N$ you translate it to the $N$ from the above problem by subtracting one, then map, and add 1 to compensate, so you get
$$
N \mapsto 1 + [(N-1)//n], 1+[(N-1) \pmod n]
$$

Useful to note that on most computing systems, the first operation is denoted div and the second mod. In addition, some have the faster divmod which does both.


TITLE: Proving that $X_n= 2^n 1_{(0,\frac{1}{n})}, n \geq 1 $ converges to zero in probability
QUESTION [0 upvotes]: Let $\Omega =[0,1],$ the sigma algebra: $F=B(\mathbb{[0,1]})$ and $\mathbb{P}=m$ (Lebesgue measure)
$X_n$ is a random variable: $X_n= 2^n 1_{(0,\frac{1}{n})}, n \geq 1 $
Prove directly, that is,  using the definition, that $X_n \overset{p}\to 0$
The definition of convergence in probability states that If $X, X_n$ are random variables   $\forall \varepsilon \geq 0, \mathbb{P}(|X_n| > \varepsilon) \overset{n \to \infty}\to 0$
My lecture notes give this short proof, that I am struggling to understand:
$\forall \varepsilon \geq 0, \mathbb{P}(|X_n| > \varepsilon)=\mathbb{P}(X_n=2^n)=\mathbb{P}((0,\frac{1}{n}))=\frac{1}{n} \overset{n \to \infty}\to 0$
I don't understand how they go from $\mathbb{P}(|X_n| > \varepsilon)$ to $\mathbb{P}(X_n=2^n)$ to $ \mathbb{P}((0,\frac{1}{n}))$. I trying to separate it into cases according to the value of $\varepsilon$, but I have to consider that the indicator function eventually becomes zero for sufficiently large n and at a fixed $\omega \in \Omega$. Can someone explain it in detail?

REPLY [1 votes]: Notice that for a fixed $n \in \mathbb{N}$,  $X_n(\omega) > 0$ if and only if $\omega \in (0, \frac{1}{n})$. Thus for any $\epsilon > 0$, $P(|X_n| > \epsilon)= P((0,\frac{1}{n})) = \frac{1}{n}$.
Edit: changed $x$ to $\omega$ to make notation a bit clearer


TITLE: An example of a Hilbert's space with application in ODE
QUESTION [1 upvotes]: Problem: I'm trying to prove that the space $H=\{v\in H^{1}(I): v(0)=v(1)\}$ is a Hilbert's space equpped with the $H^{1}$ scalar product.

We know that  and $$\forall u,v \in H^{1}(I):\quad \langle u,v \rangle_{H^{1}(I)}=\int_{I}uv+\int_{I}\nabla u \cdot \nabla v$$
and the induced norm is
$$||v||_{H^{1}(I)}^{2}=||v||_{L^{2}(I)}^{2}+||\nabla v||^{2}_{[L^{2}(I)]^{d}}$$
Let $\{v_{n}\}$ a Cauchy sequence in $H$, so I need to prove that $v_{n}\to v\in H$. Now, we pick a convergent sequence $\{v_n\} \subset H$ such that $v_n \to v$, so I need to show that $v \in H$. Now, how  $\{v_n\}$ converges in the norm $||\cdot ||_{H^1(I)}$  then $v_n \to v$ uniformly,

But how can I continue form here? Perhaps there is an easier way to approach the solution of the problem? or maybe I need to use a more strong result as the Sobolev-embedding?
Apologize me if my solution-approach is not right, I'm a bit confused with this space and your approach with the solution.

REPLY [3 votes]: Since you already know that $H^1((0,1))$ is complete, it suffices to show that $H$ is a closed subspace of $H^1$. For that purpose note that $H$ is the kernel of the continuous map
$$
H^1((0,1))\hookrightarrow C([0,1])\overset{\Phi}{\to}\mathbb{K},
$$
where the first arrow is the continuous embedding guaranteed by the Sobolev embedding theorem and $\Phi(v)=v(1)-v(0)$, which is obviously continuous (as discussed in the comments). The kernel of a continuous linear map is a closed subspace, so we are done.


TITLE: Solving system of equations which contain sin and cos
QUESTION [3 upvotes]: I require some help to push me in the right direction to solve these equations.
$$t_1 = P_1\sin(A)\sin(B) + P_2\cos(A)\cos(B)$$
$$t_2 = P_3\cos(A)\sin(B) + P_4\sin(A)\cos(B)$$
where $t_1, t_2, P_1, P_2, P_3, P_4$ are known coefficients. Solving for $A$ and $B$.
Any help greatly appreciated.
This may seem like it, but this is not homework

REPLY [4 votes]: Using the product to sum formulas we can convert your equations into
$$ e_1 = a \cos x + b \cos y$$
$$ e_2 = c \sin x + d \sin y$$
where $x = A+B$ and $y = A-B$.
This gives us
$$ \left(\frac{a \cos x - e_1}{b}\right)^2 + \left(\frac{c \sin x -e_2}{d}\right)^2 = 1$$
which can be made into a 4th degree polynomial in terms of $\cos x$: by moving $\sin x$ term to the right and using $\sin^2 x = 1 - \cos ^2 x$ and squaring, which is solvable in "closed" form and standard techniques exist to do that.
Once you solve for $\cos x$, it is easy to solve for $\cos y$ and thus get back $A$ and $B$.


TITLE: By using Wallis’s product, prove $\zeta(2)=\frac{\pi^2}{6}$
QUESTION [3 upvotes]: $\sin x = x\left(1-\frac{x^2}{\pi^2}\right)\left(1-\frac{x^2}{(2\pi)^2}\right)\cdots$
This equal is used in both the solution of  Basel’s problem and the proof of Wallis’s product.
So I wonder if I can prove one by using the other.
Help me solve this problem.
Or if you are sure this assumption is wrong, tell me the reason.
I’m sorry if I put wrong tugs.

REPLY [0 votes]: You get the $\,$ Wallis product $\,$, if you set $\,\displaystyle x=\frac{\pi}{2}\,$ .
Basel problem :
The coefficient for $\,z\,$ of the infinite product $\,(1+a_1 z)(1+a_2 z)...\,$ is $\,a_1+a_2+...\,$ .
Therefore with $z:=-x^2$ and using the series of sine we get $\,\displaystyle\frac{1}{3!}=\frac{1}{(1\cdot \pi)^2}+\frac{1}{(2\cdot \pi)^2}+...\,$ .

Both calculations are based on partial fraction expansion or equivalent on the sine product :
$($ Wallis product $\displaystyle )^2 \,\,= \left(\prod\limits_{k=1}^\infty \left(1-\left(\frac{x}{k}\right)^2\right)^{-1}|_{x=\frac{1}{2}}\right)^2=\left(\frac{\pi x}{\sin(\pi x)}|_{x=\frac{1}{2}}\right)^2=\frac{\pi^2}{4}$
$\displaystyle \frac{\pi^2}{6}=\frac{1-\pi x \cot(\pi x)}{2 x^2}|_{x\to 0}=\sum\limits_{k=0}^\infty \frac{1}{k^2-x^2}|_{x=0}=\zeta(2) $
Both calculations can be put together to an equation chain. The core thing is:
$\displaystyle \frac{2}{3}\prod\limits_{k=1}^\infty \left(1-\left(\frac{x+\frac{1}{2}}{k}\right)^2\right)^{-2} \bigg|_{x=0} =\frac{\pi^2}{6} = -\frac{1}{2}\frac{d^2}{dx^2}\ln \prod\limits_{k=1}^\infty \left(1-\left(\frac{x}{k}\right)^2\right)\bigg|_{x=0} $


TITLE: Statements equivalent to the Axiom of Choice
QUESTION [13 upvotes]: The Axiom of Choice reads:

The product of a collection of non-empty sets is non-empty.

As you know well, this axiom is equivalent to many other statements. A few examples (probably the most known) are the following:

Zorn's lemma: Let $(P,\ge)$ be a poset in which every chain has an upper bound. Then $P$ has a maximal element.

This is used in a number of different places to prove very powerful theorems. A couple of examples I can think of right now are: the proof that every ring has a maximal ideal (commutative algebra) and the Hahn-Banach theorem (functional analysis).

Well ordering theorem: Every set can be well-ordered (i.e. it admits a total order such that every subset has a least element).
Tychonoff's theorem: Every product of compact topological spaces is compact.

This is used for example in the proof of Alaoglu's theorem (again functional analysis).
Wikipedia gives various more. Among them:

Tarski's theorem: If $A$ is an infinite set, then there is a bijection $A\to A\times A$.

(As @Wicht dutifully commented, Tarski's theorem is the name given to the implication $|A|=A\times A \Rightarrow$ Axiom of Choice, while the other implication was known before the proof of this fact by Tarski.)

Every vector space has a basis.

My question is: what are other useful statements which are equivalent to the axiom of choice? Where are they used, and to prove what?
P.s.: By "useful statements" I mean either statements which are important theorems by themselves (as Tychonoff's theorem) or that are directly used in (simple versions of) the proof of important theorems (as Zorn's lemma).

REPLY [12 votes]: From the top of my head, some "big" equivalents to the axiom of choice.

Every commutative ring with a unity has a maximal ideal.
Every surjection has an injective inverse.
Every free abelian group is projective.
Every divisible abelian group is injective.
Skolem-Lowenheim theorems in model theory.
Every two cardinals are comparable. This is used often in "hiding" where we have two sets and we know that one necessarily injects into the other.
Every spanning set includes a basis (which is stronger than just the existence of a basis).
Every set is included in $L[A]$ for some $A$. Where $L[A]$ is a model of $\sf ZFC$. This can be used to prove certain things in set theory.

There are many many many more, and one can scrounge through Rubin & Rubin Equivalents of the Axiom of Choice II to find many of them including proofs.