# Implicit Regularization Properties of Variance Reduced Stochastic Mirror
Descent

Yiling Luo, Xiaoming Huo, Yajun Mei

(April 2022)

###### Abstract

In machine learning and statistical data analysis, we often run into objective
function that is a summation: the number of terms in the summation possibly is
equal to the sample size, which can be enormous. In such a setting, the
stochastic mirror descent (SMD) algorithm is a numerically efficient
method—each iteration involving a very small subset of the data. The variance
reduction version of SMD (VRSMD) can further improve SMD by inducing faster
convergence. On the other hand, algorithms such as gradient descent and
stochastic gradient descent have the implicit regularization property that
leads to better performance in terms of the generalization errors. Little is
known on whether such a property holds for VRSMD. We prove here that the
discrete VRSMD estimator sequence converges to the minimum mirror interpolant
in the linear regression. This establishes the implicit regularization
property for VRSMD. As an application of the above result, we derive a model
estimation accuracy result in the setting when the true model is sparse. We
use numerical examples to illustrate the empirical power of VRSMD.

##  1 Introduction

In statistics and machine learning, it is common to optimize an objective
function that is a finite-sum. SMD efficiently optimizes such an objective by
using a subset of data to do one step update of the variable/parameter.
Further adopting the variance reduction technique to SMD, we get the VRSMD
algorithm that enjoys fast convergence [1, 2].

The implicit regularization is a relatively new concept [3] that explains why
a result of an algorithm generalizes well in some overparameterized models [3,
4]. It refers to the fact that an algorithm can automatically select a minimum
norm solution, which is not explicitly induced by the objective function.
There are works on implicit regularization for Gradient Descent [5, 6, 7, 8],
Stochastic Gradient Descent [9, 10, 11, 12], and Stochastic Mirror Descent
[13]. Considering the computational advantage of VRSMD compared to all the
algorithms above, it would be even better if VRSMD also has the useful
implicit regularization property.

From technical point of view, our work contains the following two results:

  * •

In linear regression (including underfitting and overfitting), we show that
the solution sequence of VRSMD converges to the minimum mirror interpolant,
which is the implicit regularization property of VRSMD, and we also specify
the convergence rate.

  * •

In sparse regression, by choosing a proper mirror map, we show that the
implicit regularization estimator finds the sparse true parameter with a small
error. Moreover, compared with the deterministic algorithms in [5, 6], our
algorithm is equally good in estimating a sparse truth while being
computationally faster, as supported by our experiments.

From the application point of view, our result shows that the Mirror Descent
and its variants are useful to explore the low dimensional geometric structure
from high dimensional data, which leads to nice generalization properties.

Notation. The following notations are used throughout this paper. For a matrix
$X\in\mathbbm{R}^{n\times p}$, we denote by
$\mathrm{col}(X):=\\{\mathbf{u}\in\mathbbm{R}^{n}:\exists\mathbf{v}\in\mathbbm{R}^{p},\mathbf{u}=X\mathbf{v}\\}$
the column space of $X$, and we denote by
$\mathcal{N}(X):=\\{\mathbf{v}\in\mathbbm{R}^{p}:X\mathbf{v}=\mathbf{0}\\}$
the null space of $X$. For a vector $\mathbf{v}\in\mathbbm{R}^{p}$, we use the
definition of $\ell_{p}$ norm of $\mathbf{v}$ that
$\|\mathbf{v}\|_{p}=(\sum_{i}|v_{i}|^{p})^{1/p}$ for $p\geq 1$ and we denote
the number of non-zero elements in $\mathbf{v}$ as $\|\mathbf{v}\|_{0}$. For a
subset of indexes $I\subset\\{1,\ldots,p\\}$, we define
$\mathbf{v}_{I}:=(v_{i})_{i\in I}$, and denote the cardinality of $I$ as
$|I|$. For a set $\mathcal{X}\subset\mathbbm{R}^{p}$, define
$P_{\mathcal{X}}\mathbf{v}=\operatorname{argmin}_{\mathbf{u}\in\mathcal{X}}\|\mathbf{u}-\mathbf{v}\|_{2}$.
For two non-negative-valued functions $a(x)$ and $b(x)$, we denote
$a(x)\sim\mathcal{O}(b(x))$ if there exists an absolute constant $C$ such that
$a(x)\leq Cb(x)$; and we denote $a(x)\sim\Theta(b(x))$ if there are absolute
constants $c,C$ such that $cb(x)\leq a\leq Cb(x)$.

Organization. The rest of the paper is as follows. In Section 2, we describe
our problem formulation and algorithm. Section 3 states the main theory on the
implicit regularization. Section 4 develops insight into the implicit
regularization and establishes the sparse recovery property. Section 5
supports the theory on implicit regularization by simulations and experiments.
In Section 6, we discuss the finding of our work and some future directions.
All proofs are deferred to the appendix.

##  2 Formulation and Algorithm

To better present the material, we split this section into two subsections. In
Subsection 2.1 we formulate the optimization problem motivated from linear
regression. In Subsection 2.2, we present the Variance Reduction Stochastic
Mirror Descent (VRSMD) algorithm for solving such an optimization problem.

###  2.1 Formulation

Assume we observe data pairs
$\\{(\mathbf{x}_{i},y_{i})\in\mathbbm{R}^{p}\times\mathbbm{R}\\}_{i=1}^{n}$,
the goal is to predict the response $y$ based on $\mathbf{x}$. Under the
empirical risk minimization framework, we consider the general optimization
problem of the form

$\min_{\bm{\beta}}F(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\bm{\beta};(\mathbf{x}_{i},y_{i})),$
(1)

and we shorten $f_{i}(\bm{\beta};(\mathbf{x}_{i},y_{i}))$ as
$f_{i}(\bm{\beta})$ to simplify the notation.

As a concrete example, for the linear regression model, the classical least
squares method is to find coefficient $\bm{\beta}$ that minimizes the
objective function

$\min_{\bm{\beta}}F(\bm{\beta})=\frac{1}{2n}\sum_{i=1}^{n}\left(\mathbf{x}_{i}^{T}\bm{\beta}-y_{i}\right)^{2}=\frac{1}{2n}\|X\bm{\beta}-y\|_{2}^{2},$
(2)

where we denote
$X=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]^{T}\in\mathbbm{R}^{n\times p}$, and
$\mathbf{y}=[y_{1},\ldots,y_{n}]^{T}\in\mathbbm{R}^{n}$.

When problem (2) has non-unique solutions, it is important yet nontrivial to
find a solution that has nice generalization property. It is well known that
when running Gradient Descent algorithm with initialization
$\bm{\beta}_{0}=\mathbf{0}$ on (2), the corresponding solution is the minimal
$\ell_{2}$ norm solution among all solutions of (2), see [14]. Also, the
properties of SMD are studied in [13]. It is unknown what happens to the
solution if we run variants of SMD, for example, variance reduced SMD.

###  2.2 VRSMD Algorithm

Let us now present the main idea of the variance reduced stochastic mirror
descent (VRSMD) algorithm as follows.

To understand why we need variance reduction, consider the Stochastic Mirror
Descent(SMD) algorithm using a strictly convex and differentiable mirror map
$\psi(\cdot)$. At step $t$, the SMD updates $\bm{\beta}_{t+1}$ such that

$\nabla\psi(\bm{\beta}_{t+1})=\nabla\psi(\bm{\beta}_{t})-\eta_{t}\nabla
f_{i_{t}}(\bm{\beta}_{t}),$

where $i_{t}$ is randomly sampled from $\\{1,\ldots,n\\}$. Now the term
$\nabla f_{i_{t}}(\bm{\beta}_{t})$ has $\mathbbm{E}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]=\nabla F(\bm{\beta}_{t})$, so SMD has unbiased
update compared to Mirror Descent, where the update is
$\nabla\psi(\bm{\beta}_{t+1})=\nabla\psi(\bm{\beta}_{t})-\eta_{t}\nabla
F(\bm{\beta}_{t})$. However, in general $\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]\neq 0$ for any $\bm{\beta}_{t}$, so we need
$\eta_{t}\to 0$, which may lead to slow convergence.

Variance reduction addresses the issues above by replacing $\nabla
f_{i_{t}}(\bm{\beta}_{t})$ with term $A_{t}$ such that

  * •

$\mathbbm{E}[A_{t}]=\mathbbm{E}[\nabla f_{i_{t}}(\bm{\beta}_{t})]$ to keep
unbiased update;

  * •

$\operatorname{Var}[A_{t}]<\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]$ to control variance.

One choice of $A_{t}$ is $A_{t}=\nabla
f_{i_{t}}(\bm{\beta}_{t})-B_{t}+\mathbbm{E}[B_{t}]$, where $B_{t}$ and $\nabla
f_{i_{t}}(\bm{\beta}_{t})$ are positively correlated with correlation
coefficient $r>0.5$ and
$\operatorname{Var}[B_{t}]\approx\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]$. For this $A_{t}$ one can check that
$\mathbbm{E}[A_{t}]=\mathbbm{E}[\nabla f_{i_{t}}(\bm{\beta}_{t})]$ and
$\operatorname{Var}[A_{t}]=\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})-B_{t}]=\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]-2r\sqrt{\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]\operatorname{Var}[B_{t}]}+\operatorname{Var}[B_{t}]<\operatorname{Var}[\nabla
f_{i_{t}}(\bm{\beta}_{t})]$. For a proper $B_{t}$ such that
$\operatorname{Var}(A_{t})\stackrel{{\scriptstyle
t\to\infty}}{{\longrightarrow}}0$, the algorithm converges for a fixed $\eta$.

For illustration purpose, we use the variance reduction technique in [15] to
get the VRSMD Algorithm 1. However, one should note that this framework
applies to other variance reduction methods such as SARAH [16] and SPIDER
[17].

Input: An objective function $F(\cdot)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\cdot)$,
and a strictly convex and differentiable mirror map $\psi(\cdot)$;

Initialization: Initialize $\tilde{\bm{\beta}}^{0}$. Choose the step-size
$\eta$, outer iteration number $S$, inner iteration number $m$. Denote the
estimator at $t$th inner iteration of $s$th outer iteration as
$\bm{\beta}_{t}^{s}$. Set
$\bm{\beta}_{1}^{1}=\bm{\beta}_{m+1}^{0}=\tilde{\bm{\beta}}^{0}$;

for  _Outer iteration s = 1,…,S_ do

Calculate $\nabla F(\tilde{\bm{\beta}}^{s-1})$;

for  _Inner iteration t = 1,…,m_ do

Randomly sample $i_{t}$ from $\\{1,\ldots,n\\}$, calculate
$\mathbf{v}_{t}^{s}=\nabla f_{i_{t}}(\bm{\beta}_{t}^{s})-\nabla
f_{i_{t}}(\tilde{\bm{\beta}}^{s-1})+\nabla F(\tilde{\bm{\beta}}^{s-1}),$ (3)
and update $\bm{\beta}_{t+1}^{s}$ such that
$\nabla\psi(\bm{\beta}_{t+1}^{s})=\nabla\psi(\bm{\beta}_{t}^{s})-\eta\mathbf{v}_{t}^{s}.$
(4)

Set $\tilde{\bm{\beta}}^{s}$ to be a uniform random sample from
$\\{\bm{\beta}_{1}^{s},\ldots,\bm{\beta}_{m}^{s}\\}$;

Option I: Set $\bm{\beta}_{1}^{s+1}=\bm{\beta}_{m+1}^{s}$;

Option II: Set $\bm{\beta}_{1}^{s+1}=\tilde{\bm{\beta}}^{s}$;

Option I: Output $\bm{\beta}_{a}$ chosen uniformly random from
$\\{\\{\bm{\beta}_{t}^{s}\\}_{t=1}^{m}\\}_{s=1}^{S}$;

Option II: Output $\bm{\beta}_{a}=\tilde{\bm{\beta}}^{S}$.

Algorithm 1 Variance Reduced Stochastic Mirror Descent (VRSMD)

######  Remark 1.

We note the complexity of VRSMD as follows: The total number of stochastic-
first-order calls (i.e. SFO complexity) of Algorithm 1 is
$\mathcal{O}(nS+mS)$. A popular choice of the inner loop number $m$ is
$\Theta(n)$, which leads to $\mathcal{O}(1)$ SFO complexity per inner loop.

The variance reduction component of VRSMD is $\mathbf{v}_{t}^{s}$ in (3). Note
that it is equivalent to the variance reduction scheme we talked above by
taking $B_{t}=\nabla f_{i_{t}}(\tilde{\bm{\beta}}^{s-1})$. The conditions we
list there on $B_{t}$ hold when $\bm{\beta}_{t}^{s}$ and
$\tilde{\bm{\beta}}^{s-1}$ are close, which happens by taking moderate values
of the inner iteration number $m$ and the step-size $\eta$.

We show that the VRSMD algorithm is a generalization of the SVRG algorithm in
[15, 18]: For the special case of $\psi(\cdot)=\frac{1}{2}\|\cdot\|_{2}^{2}$,
we have $\nabla\psi(\bm{\beta})=\bm{\beta}$, then (3) updates
$\bm{\beta}_{t+1}^{s}$ as:

$\bm{\beta}_{t+1}^{s}=\bm{\beta}_{t}^{s}-\eta\mathbf{v}_{t}^{s},$

which is the SVRG update, and VRSMD reduces to SVRG.

##  3 Implicit Regularization

In this section, we present the implicit regularization property of the VRSMD
solution. To do so, it is necessary to first show that the VRSMD converges.

To begin with, we introduce some definitions that will be useful in our
theory.

######  Definition 1 (L-smoothness).

$f$ is $L$-smooth with respect to $\|\cdot\|$ norm if there exists a constant
$L>0$ such that

$\|\nabla f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}\leq
L\|\mathbf{u}-\mathbf{w}\|,\forall\mathbf{u},\mathbf{w},$

where $\|\cdot\|_{*}:=\max_{y:\|y\|=1}\langle y,\cdot\rangle$ is the dual norm
of $\|\cdot\|$.

######  Definition 2 ($\alpha$-strongly convex).

$f$ is $\alpha$-strongly convex with respect to $\|\cdot\|$ norm if there
exists a constant $\alpha>0$ such that

$f(\mathbf{u})\geq f(\mathbf{w})+\nabla
f(\mathbf{w})^{T}(\mathbf{u}-\mathbf{w})+\frac{\alpha}{2}\|\mathbf{u}-\mathbf{w}\|^{2},\forall\mathbf{u},\mathbf{w}.$

######  Definition 3 (Quadratic growth (QG)).

Let $\mathcal{X}$ be the set of all minimizers of $f$. $f$ satisfies QG
condition w.r.t. $\|\cdot\|$ if

$\frac{\mu}{2}\|\mathbf{u}-P_{\mathcal{X}}\mathbf{u}\|^{2}\leq
f(\mathbf{u})-f(P_{\mathcal{X}}\mathbf{u}),\forall\mathbf{u}.$ (5)

######  Definition 4 ($\epsilon$-solution).

For the optimization problem

$\displaystyle\begin{split}Opt=\min_{x\in B}\\{f(x):g_{i}(x)\leq 0,1\leq i\leq
m\\}.\end{split}$ (6)

$x_{\epsilon}\in B$ is called an $\epsilon$-solution to (6) if

$\displaystyle\begin{split}&f(x_{\epsilon})-Opt\leq\epsilon,\\\
&g_{i}(x_{\epsilon})\leq\epsilon,1\leq i\leq m.\end{split}$

######  Definition 5 (Restricted eigenvalue (RE)).

$X$ satisfies $(s,\gamma)$-RE condition if for any $\bm{\beta}$ such that
$\|\bm{\beta}\|_{0}\leq s$ we have

$\frac{\frac{1}{n}\|X\bm{\beta}\|_{2}^{2}}{\|\bm{\beta}\|_{2}^{2}}\geq\gamma.$

######  Definition 6 ($s$-good).

A matrix $X_{n\times p}$ is $s$-good if $\exists\kappa<\frac{1}{2}$ such that
$\forall\mathbf{u}\in\mathcal{N}(X)\subset\mathbbm{R}^{p}$ and $\forall
I\subset\\{1,\ldots,p\\}$ with $|I|\leq s$, we have

$\|\mathbf{u}_{I}\|_{1}\leq\kappa\|\mathbf{u}\|_{1}$.

Next, let us present the convergence result of VRSMD:

######  Proposition 1.

Assume $F(\cdot)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\cdot)$ has every
$f_{i}(\cdot)$ convex and $L$-smooth w.r.t. an arbitrary norm $\|\cdot\|$, and
$\psi(\cdot)$ is $\alpha$-strongly convex w.r.t. $\|\cdot\|$. Denote
$\bm{\beta}^{*}=\arg\min F(\cdot)$.

(a) Run Option I of Algorithm 1 on $F$ with $\eta<\frac{\alpha}{24L}$, then we
have

$\displaystyle\begin{split}\mathbbm{E}[F(\bm{\beta}_{a})-F(\bm{\beta}^{*})]\leq\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\times\\\
\quad\left[D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}(F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{*}))\right],\end{split}$
(7)

where $T=m\cdot S$, and
$D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{0}):=\psi(\bm{\beta}^{*})-\psi(\tilde{\bm{\beta}}^{0})-\langle\nabla\psi(\tilde{\bm{\beta}}^{0}),\bm{\beta}^{*}-\tilde{\bm{\beta}}^{0}\rangle$
is the Bregman divergence.

(b) If we further assume that $F(\cdot)$ satisfies the QG condition in (5)
with constant $\mu$, and that $\psi(\cdot)$ is $\ell$-smooth, all w.r.t.
$\|\cdot\|$, and also suppose that we run Option II of Algorithm 1 with a
large enough $m$ such that

$\tau:=\frac{12L\eta^{2}/\alpha+\ell/(m\mu)}{\eta-12L\eta^{2}/\alpha}<1,$ (8)

then the VRSMD has a stronger linear convergence rate:

$\mathbbm{E}[F(\bm{\beta}_{a})-F(\bm{\beta}^{*})]\leq\tau^{S}[F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{*})].$
(9)

######  Remark 2.

We analyze the computational complexity implied by Proposition 1: In (a), let
$m=n$ and take $\eta=\frac{\alpha}{48L}$, then we have

$\displaystyle\mathbbm{E}$
$\displaystyle[F(\bm{\beta}_{a})-F(\bm{\beta}^{*})]\leq\frac{96L}{\alpha
T}\times$
$\displaystyle\left[D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{0})+\frac{\alpha
n}{192L}(F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{*}))\right].$

In this case, the number of gradient computations for achieving an
$\epsilon$-solution is $\mathcal{O}(\frac{L}{\epsilon}+\frac{n}{\epsilon})$.

######  Remark 3.

The assumption in (b) is moderate. Take $m=\frac{110L\ell}{\alpha\mu}$ and
$\eta=\frac{\alpha}{36L}$, we have $\tau<1$. This choice of $m$ does not
violate the $\mathcal{O}(1)$ SFO comlexity per iteration – take a good mirror
map so that $\ell/\alpha=\mathcal{O}(1)$, and consider the most indicative
case [15] where the condition number $L/\mu=n$, we have $m=\Theta(n)$.

Our results in Proposition 1 are consistent with those for SVRG. In part (a),
the $\mathcal{O}(1/T)$ convergence rate matches the rate in [18]; in part (b),
the linear rate matches the rate in [15], while we reduce their strong
convexity assumption to quadratic growth.

Finally, we are ready to present our main result on implicit regularization.
In the following theorem, we show that VRSMD finds an $\epsilon$-solution of
the minimum mirror interpolation problem.

######  Theorem 1.

For the objective function in (2), assume $\psi(\bm{\beta})$ is
$\alpha$-strongly convex w.r.t. $\|\cdot\|_{2}$, denote
$L=\max_{i}\|\mathbf{x}_{i}\|_{2}^{2}$ and let $s_{m}$ be the smallest nonzero
singular value of $X$. The VRSMD algorithm converges to the minimum mirror map
interpolant

$\displaystyle\begin{split}\bm{\beta}^{\psi}:=\operatorname{argmin}_{\bm{\beta}}&\quad\psi(\bm{\beta})\\\
s.t.&\quad
F(\bm{\beta})=\min_{\bm{\beta}^{\prime}}F(\bm{\beta}^{\prime}).\end{split}$
(10)

We describe the convergence by the following $\epsilon$-solution:

(a) Run Option I of Algorithm 1 with choice $\eta<\frac{\alpha}{24L}$ and
initialization $\tilde{\bm{\beta}}^{0}$ such that
$\nabla\psi(\tilde{\bm{\beta}}^{0})\in\mathrm{col}(X^{T})$, assume that the
output $\bm{\beta}^{a}$ satisfies $\|\nabla\psi(\bm{\beta}^{a})\|_{2}\leq B$,
then we will have

$\displaystyle\begin{split}&\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]\leq\frac{B}{s_{m}}\sqrt{\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}}\times\\\
&\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{24nL\eta^{2}m}{\alpha}\left(F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})\right)\right]^{.5},\end{split}$
(11)

which describes how far is the objective value at $\bm{\beta}^{a}$ away from
the optimal solution to (10). Moreover, we have

$\displaystyle\begin{split}&\mathbbm{E}[F(\bm{\beta}^{a})-F(\bm{\beta}^{\psi})]\leq\frac{\alpha}{(\alpha\eta-8L\eta^{2})T}\times\\\
&\left[D_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})\right)\right],\end{split}$
(12)

which characterizes how much does $\bm{\beta}^{a}$ violates the constraints of
(10). They together show that the VRSMD algorithm finds an $\epsilon$-solution
to (10) for $T=\mathcal{O}(\frac{1}{\epsilon}+\frac{1}{\epsilon^{2}})$.

(b) Further assume that $\psi(\cdot)$ is $\ell$-smooth w.r.t. $\|\cdot\|_{2}$
and

$\tau^{\prime}=\frac{12L\eta^{2}/\alpha+\ell
n/(ms_{m}^{2})}{\eta-12L\eta^{2}/\alpha}<1.$ (13)

Run Option II of Algorithm 1, we can show that:

$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]$
$\displaystyle\leq\frac{B(\tau^{\prime})^{S/2}\sqrt{2n}}{s_{m}}\sqrt{F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})},$
$\displaystyle\mathbbm{E}[F(\bm{\beta}^{a})-F(\bm{\beta}^{\psi})]$
$\displaystyle\leq(\tau^{\prime})^{S}\left(F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})\right).$
(14)

######  Remark 4.

We need to point out that the assumptions in Theorem 1 are moderate. For
instance, for the assumption that
$\nabla\psi(\tilde{\bm{\beta}}^{0})\in\mathrm{col}(X^{T})$, we can take
$\tilde{\bm{\beta}}^{0}=(\nabla\psi)^{-1}(X^{T}\mathbf{a})$ for any
$\mathbf{a}$. One feasible choice is $\mathbf{a}=\mathbf{0}$, resulting in
$\tilde{\bm{\beta}}_{0}=\arg\min_{\bm{\beta}\in R^{p}}\psi(\bm{\beta})$. Since
$\psi$ is strongly convex, this minimizer is not hard to calculate, for
example: we have
$\cdot\psi(\cdot)=\|\cdot\|_{q}^{q}$ or $\psi(\cdot)=\|\cdot\|_{q}^{2}$ for
$q>1\Rightarrow$ $\arg\min\psi(\cdot)=\mathbf{0}$;
$\cdot\psi(\bm{\beta})=\bm{\beta}^{T}H\bm{\beta}$ for a positive definite $H$
$\Rightarrow\arg\min\psi(\cdot)=\mathbf{0}$;
$\cdot\psi(\bm{\beta})=\sum_{i=1}^{p}{\beta}_{i}\log({\beta}_{i})-{\beta}_{i}$
$\Rightarrow\arg\min\psi(\cdot)=\mathbf{1}$.

######  Remark 5.

As for the assumption in (b), we can take $m=\frac{110L\ell n}{\alpha
s_{m}^{2}}$ and $\eta=\frac{\alpha}{36L}$ to get
$\tau^{\prime}=(1+108/110)/2<1$. Take a good mirror map such that
$\ell/\alpha=\mathcal{O}(1)$ and assume $L/s_{m}^{2}=\mathcal{O}(1)$, we
further have $m=\Theta(n)$, so the algorithm can be implemented efficiently.

It is useful to provide a high level understanding of Theorem 1. It implies
that the discrete update of the VRSMD Algorithm on the unregularized objective
(2) is an $\epsilon$-solution of the regularized optimization problem (10), so
it is an implicit regularization result. Furthermore, since (10) minimizes a
strictly convex function over a convex set, the solution will be unique, thus
the estimator from VRSMD must converge to this unique solution. The finding in
Theorem 1 can be further extended to other variants of SMD.

##  4 Further Understanding of Implicit Regularization

In this section, we provide a deeper understanding of our theoretical results
in the previous section by analyzing two special mirror maps for linear
regression model: one is $\psi(\bm{\beta})=\|\bm{\beta}\|_{2}^{2}/2$, the
other is $\psi(\bm{\beta})=\|\bm{\beta}\|_{1+\delta}^{1+\delta}$ for small
$\delta>0$.

Let us first consider $\psi(\cdot)=\|\cdot\|_{2}^{2}/2$, where VRSMD reduces
to SVRG algorithm in [15]. In this case, we further show that
$\|\bm{\beta}^{a}-\bm{\beta}^{\psi}\|_{2}^{2}$ linearly converges to $0$,
where $\bm{\beta}^{\psi}$ is the minimum $\ell_{2}$ norm solution
$\bm{\beta}^{\psi}=(X^{T}X)^{+}X^{T}\mathbf{y}=X^{+}\mathbf{y}$:

######  Corollary 1.

Denote $L=\max_{i}\|\mathbf{x}_{i}\|_{2}^{2}$, take
$\psi(\cdot)=\|\cdot\|_{2}^{2}/2$, let
$\tilde{\bm{\beta}}^{0}=\bm{\beta}_{m}^{0}=\mathbf{0}$,
$\eta_{t}=\eta<1/(24L)$, and assume

$\tau^{\prime\prime}=\frac{12L\eta^{2}+n/(ms_{m}^{2})}{\eta-12L\eta^{2}}<1.$
(15)

Run Option II of Algorithm 1 on (2), we have

$\mathbbm{E}\|\bm{\beta}^{a}-X^{+}\mathbf{y}\|_{2}^{2}\leq\frac{(\tau^{\prime\prime})^{S}}{s_{m}^{2}}\|P_{\mathrm{col}(X)}\mathbf{y}\|_{2}^{2}.$
(16)

Next, let us consider the mirror map
$\psi(\bm{\beta})=\|\bm{\beta}\|_{1+\delta}^{1+\delta}$ for a small $\delta>0$
in VRSMD, which leads to a sparse solution. Assume that
$\|\bm{\beta}_{t}^{s}\|_{\infty}\leq K$ for a large enough $K$ throughout the
updates, we then have $\psi(\bm{\beta}_{t}^{s})$ is
$\frac{(1+\delta)\delta}{K^{1-\delta}}$-strongly convex and
$\|\nabla\psi(\bm{\beta}^{a})\|_{2}\leq\sqrt{p}(1+\delta)K^{\delta}$. By
Theorem 1 we have the estimator sequence converges to the penalized solution

$\displaystyle\bm{\beta}^{(\delta)}=\operatorname{argmin}_{\bm{\beta}}\left\\{\|\bm{\beta}\|_{1+\delta}^{1+\delta}:X\bm{\beta}=P_{\mathrm{col}(X)}\mathbf{y}\right\\}.$
(17)

This choice of mirror map has $\lim_{\delta\to
0}\|\bm{\beta}\|_{1+\delta}^{1+\delta}=\|\bm{\beta}\|_{1}$ and the solution
$\bm{\beta}^{(0)}:=\operatorname{argmin}_{\bm{\beta}}\\{\|\bm{\beta}\|_{1}:X\bm{\beta}=P_{\mathrm{col}(X)}\mathbf{y}\\}$
is sparse. Thus we expect that for small $\delta$, $\bm{\beta}^{(\delta)}$ is
close to $\bm{\beta}^{(0)}$ and recovers a sparse true parameter.

We now provide a rigorous argument for the sparse recovery. Assume that the
data is generated by $\mathbf{y}=X\bm{\beta}^{o}$ for a sparse
$\bm{\beta}^{o}$. In the following theorem we have $\bm{\beta}^{(\delta)}$
accurately recovers $\bm{\beta}^{o}$ when the design matrix $X$ satisfies some
proper conditions.

######  Theorem 2 (Sparse Recovery).

Under the sparse setting defined above, denote $s=\|\bm{\beta}^{o}\|_{0}$.
Assume that the design matrix $X$ satisfies $(s,\gamma)$-RE condition and is
$s$-good with constant $\kappa<\frac{1}{2}$. For any $\xi>0$, if we choose

$\delta\leq\frac{\log\left(1+\frac{(1-2\kappa)\sqrt{n\gamma}}{\sqrt{s}\|\mathbf{y}\|_{2}}\xi\right)}{\log
p-\log\left(1+\frac{(1-2\kappa)\sqrt{n\gamma}}{\sqrt{s}\|\mathbf{y}\|_{2}}\xi\right)},$
(18)

we have

$\|\bm{\beta}^{(\delta)}-\bm{\beta}^{o}\|_{1}\leq\xi.$ (19)

By Theorem 2, the estimator $\bm{\beta}^{(\delta)}$ estimates the sparse truth
with a small error $\xi$. In this way, VRSMD algorithm 1 achieves near sparse
recovery via implicit regularization.

Figure 1: Run VRSMD Algorithm on simulated noiseless data. (a) the squared
error objective function converges quickly to $0$. In (b), the $\ell_{1}$
estimation error converges to smaller value for smaller $\delta$, indicating
the nearly exact recovery of the sparse signal, which supports Theorem 2. In
(c), for $\psi(\cdot)=\|\cdot\|_{1.05}^{1.05}$, the VRSMD estimator converges
to the true parameter values. In (d), for a smaller $\delta$, the convergence
value of VRSMD estimator is closer to ground truth.

##  5 Numerical Experiment

Simulation. We generate data by a sparse model as follows: Set
$n=1000,p=5000>n$ for the design matrix X. Simulate $X=\Sigma^{1/2}W$ where
the entries of $W$ are i.i.d. $N(0,1)$ and
$\Sigma=0.5*\operatorname{diag}(\mathbf{1}_{n})+0.5*\mathbf{1}_{n\times n}$.
The true parameter $\bm{\beta}^{o}\in\mathbbm{R}^{p}$ has its first $30$
entries sampled from i.i.d. $N(0,1)$ and the rest entries set to $0$. Compute
responses $\mathbf{y}=X\bm{\beta}^{o}$.

We then run VRSMD on objective function (2) for this simulated $X$ and
$\mathbf{y}$. For a range of $\delta$, set the mirror map as
$\psi(\cdot)=\|\cdot\|_{1+\delta}^{1+\delta}$, and run VRSMD with
initialization $\tilde{\bm{\beta}}^{0}=\mathbf{0}$, step-size $\eta=0.0002$,
outer iteration number $50$ and inner iteration number $1000=n$. The result is
in Fig. 1.

Experiment on RNA dataset. We use the gene expression cancer RNA-Seq data
set111https://archive.ics.uci.edu/ml/datasets/gene+expression+cancer+RNA-Seq
for experiment. The data consists of $801$ observations, each of dimension
$20,531$. Randomly split the data into $600$ training data and $201$ testing
data. Run VRSMD algorithm on training data using mirror function
$\psi=\|\cdot\|_{1.1}^{1.1}$ where initialization
$\tilde{\bm{\beta}}^{0}=\mathbf{0}$, step-size $\eta=0.015$, inner iteration
number $400$ and outer iteration number determined by $5$-fold cross
validation (i.e. early stopping). We also compare VRSMD with the Hadamard GD
[6, 5], which also has implicit regularization for sparsity. The result is in
Fig. 2.

Figure 2: Run VRSMD Algorithm on RNA dataset. In (a), we plot the solution
path of the first $100$ entries of $\bm{\beta}$, and it shows that the early
stopped estimator is sparse. In (b) and (c), we compare the performance of
VRSMD with Hadamard GD. Plot (b) shows that VRSMD trains faster than Hadamard
GD, which is tested significant by one-sided Wilcoxon signed-rank test
(p-value = $9.77\times 10^{-4}$). Plot (c) shows that the two algorithms have
same prediction error on testing data, for which we test by two-sided Wilcoxon
signed-rank test (p-value = $0.56$).

##  6 Conclusions and Future Research

Our work analyzes the implicit regularization property of VRSMD, which covers
both underfitting and overfitting cases in linear regression. In particular,
our theorem shows that the implicit regularization property can help VRSMD
find a sparse ground truth with a small error. Our experiments illustrate that
the VRSMD is computationally efficient compared to the Hadamard GD algorithm
that has implicit regularization for sparsity [6], thanks to the stochastic
nature of VRSMD.

We discuss some future directions of our research. First, it is useful to
study the implicit regularization properties of the VRSMD in the nonEuclidean
setup. For example, one can consider the generalized linear model (GLM) where
the data lies in a Riemannian manifold and mirror descent and/or natural
gradient descent are efficient. Our analysis can be extended from the linear
regression model to the GLM case.

Second, it is interesting to investigate the minimax property of the VRSMD
estimator. We see from our experiment that VRSMD with early stopping has
comparable prediction performance to Hadamard GD that is minimax optimal for
sparse regression [6, 5]. This leads to the open question of how to select a
good mirror map and an optimal stopping time in VRSMD (or variants of VRSMD)
such that the resulting estimator is minimax optimal and thus generalizes
well.

Third, one can explore the implicit regularization properties of other
variants of SMD. For example, one can consider the accelerated version of
VRSMD as the Katyusha algorithm in [19]. Such an algorithm is well studies in
optimization literature, and it has a better convergence rate that can match
the theoretical optimum. We hypothesize that it also enjoys the implicit
regularization property, but the proof is out of the scope of this work.

Finally, our results might provide a better understanding of deep neural
networks. By [20], Gradient Descent on Hadamard reparameterized linear
regression, which is related to a neural network with multiple layers, can be
approximated by Mirror Descent on original parameters. This point of view
allows us to study a neural network from the VRSMD perspective and helps to
explain why the Gradient Descent gives a sparse estimator in some deep
learning models.

## References

  * [1] Wenjie Li, Zhanyu Wang, Yichen Zhang and Guang Cheng  “Variance Reduction on Adaptive Stochastic Mirror Descent”  In _arXiv preprint arXiv:2012.13760_ , 2021
  * [2] Zeyuan Allen-Zhu  “Katyusha: The first direct acceleration of stochastic gradient methods”  In _The Journal of Machine Learning Research_ 18.1 JMLR. org, 2017, pp. 8194–8244
  * [3] Peter L Bartlett, Andrea Montanari and Alexander Rakhlin  “Deep learning: a statistical viewpoint”  In _arXiv preprint arXiv:2103.09177_ , 2021
  * [4] Michal Derezinski, Feynman T Liang and Michael W Mahoney  “Exact expressions for double descent and implicit regularization via surrogate random design”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 5152–5164  URL: https://proceedings.neurips.cc/paper/2020/file/37740d59bb0eb7b4493725b2e0e5289b-Paper.pdf
  * [5] Tomas Vaskevicius, Varun Kanade and Patrick Rebeschini  “Implicit Regularization for Optimal Sparse Recovery”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper/2019/file/5cf21ce30208cfffaa832c6e44bb567d-Paper.pdf
  * [6] Peng Zhao, Yun Yang and Qiao-Chu He  “Implicit regularization via Hadamard product over-parametrization in high-dimensional linear regression”  In _arXiv preprint arXiv:1903.09367_ , 2019
  * [7] Jianqing Fan, Zhuoran Yang and Mengxin Yu  “Understanding Implicit Regularization in Over-Parameterized Nonlinear Statistical Model”, 2021  arXiv:2007.08322 [stat.ML]
  * [8] Sanjeev Arora, Nadav Cohen, Wei Hu and Yuping Luo  “Implicit Regularization in Deep Matrix Factorization”  In _Advances in Neural Information Processing Systems_ 32 Curran Associates, Inc., 2019  URL: https://proceedings.neurips.cc/paper/2019/file/c0c783b5fc0d7d808f1d14a6e9c8280d-Paper.pdf
  * [9] Alnur Ali, Edgar Dobriban and Ryan Tibshirani  “The implicit regularization of stochastic gradient flow for least squares”  In _International Conference on Machine Learning_ , 2020, pp. 233–244  PMLR
  * [10] Shahar Azulay et al.  “On the Implicit Bias of Initialization Shape: Beyond Infinitesimal Mirror Descent”  In _arXiv preprint arXiv:2102.09769_ , 2021
  * [11] Samuel L Smith, Benoit Dherin, David Barrett and Soham De  “On the Origin of Implicit Regularization in Stochastic Gradient Descent”  In _International Conference on Learning Representations_ , 2021  URL: https://openreview.net/forum?id=rq_Qr0c1Hyo
  * [12] Guy Blanc, Neha Gupta, Gregory Valiant and Paul Valiant  “Implicit regularization for deep neural networks driven by an Ornstein-Uhlenbeck like process”  In _Proceedings of Thirty Third Conference on Learning Theory_ 125, Proceedings of Machine Learning Research  PMLR, 2020, pp. 483–513  URL: http://proceedings.mlr.press/v125/blanc20a.html
  * [13] Navid Azizan and Babak Hassibi  “Stochastic Gradient/Mirror Descent: Minimax Optimality and Implicit Regularization”  In _International Conference on Learning Representations_ , 2019  URL: https://openreview.net/forum?id=HJf9ZhC9FX
  * [14] Jingfeng Wu, Difan Zou, Vladimir Braverman and Quanquan Gu  “Direction Matters: On the Implicit Bias of Stochastic Gradient Descent with Moderate Learning Rate”  In _International Conference on Learning Representations_ , 2021  URL: https://openreview.net/forum?id=3X64RLgzY6O
  * [15] Rie Johnson and Tong Zhang  “Accelerating stochastic gradient descent using predictive variance reduction”  In _Advances in Neural Information Processing Systems_ 26, 2013, pp. 315–323
  * [16] Lam M. Nguyen, Jie Liu, Katya Scheinberg and Martin Takáč  “SARAH: A Novel Method for Machine Learning Problems Using Stochastic Recursive Gradient”  In _Proceedings of the 34th International Conference on Machine Learning_ 70, Proceedings of Machine Learning Research  PMLR, 2017, pp. 2613–2621  URL: http://proceedings.mlr.press/v70/nguyen17b.html
  * [17] Cong Fang, Chris Junchi Li, Zhouchen Lin and Tong Zhang  “SPIDER: Near-Optimal Non-Convex Optimization via Stochastic Path-Integrated Differential Estimator”  In _Advances in Neural Information Processing Systems_ 31 Curran Associates, Inc., 2018  URL: https://proceedings.neurips.cc/paper/2018/file/1543843a4723ed2ab08e18053ae6dc5b-Paper.pdf
  * [18] Sashank J Reddi et al.  “Stochastic variance reduction for nonconvex optimization”  In _International Conference on Machine Learning_ , 2016, pp. 314–323
  * [19] Zeyuan Allen-Zhu  “Katyusha: The First Direct Acceleration of Stochastic Gradient Methods”, 2018  arXiv:1603.05953 [math.OC]
  * [20] Tomas Vaskevicius, Varun Kanade and Patrick Rebeschini  “The Statistical Complexity of Early-Stopped Mirror Descent”  In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 253–264  URL: https://proceedings.neurips.cc/paper/2020/file/024d2d699e6c1a82c9ba986386f4d824-Paper.pdf
  * [21] Hamed Karimi, Julie Nutini and Mark Schmidt  “Linear convergence of gradient and proximal-gradient methods under the Polyak-Łojasiewicz condition”  In _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ , 2016, pp. 795–811  Springer
  * [22] Aaron Sidford  “Introduction to Optimization Theory Lecture NOtes Chapter 5”, 2019  URL: https://web.stanford.edu/~sidford/courses/19fa_opt_theory/sidford_mse213_2019fa_chap_5_smooth_convex_gen.pdf
  * [23] Yurii Nesterov  “Introductory lectures on convex programming volume i: Basic course”  In _Lecture notes_ 3.4, 1998, pp. 5
  * [24] Aharon Ben-Tal and Arkadi Nemirovski  “Lectures on modern convex optimization (2012)”  In _SIAM, Philadelphia, PA. Google Scholar Google Scholar Digital Library Digital Library_ , 2011

##  Appendix A Comparison with the VRSMD work [1]

Besides the implicit regularization property that is emphasised in this work,
the variance reduced stochastic mirror descent (VRSMD) itself is an algorithm
of independent interest to the optimization community. We notice a concurrent
work that also applies variance reduction to SMD [1]. Our work is
independently done without been affected by their work, and we have different
focuses. The high level difference of their work and this one is that: their
main goal is the convergence analysis, while this work focus on the
interesting implicit regularization property and its implication of VRSMD.
Nevertheless, we can compare the overlapping part of our work with theirs,
which is the convergence analysis. Here we do it from three perspectives:
assumptions, measurement of convergence, and the the number of stochastic
gradient $\nabla f_{i}(x)$ (i.e. SFO calls) required to achieve a $\epsilon$
convergence:

Difference in assumption: There are 4 main difference in assumptions:

1\. They assume their objective function to be
$\min\frac{1}{n}\sum_{i}f_{i}(x)+g(x)$, where $g$ is possibly non-smooth but
convex and $f_{i}$ is possibly non-convex but smooth. Our formulation does not
have the non-smooth component, and we assume $f_{i}$ to be convex. In their
work, they assume the proximal operator in $g$ can be efficiently solved to
handle the non-smoothness. And they use different measurement for convergence
compared to our work due to non-convexity.

2\. Their assumptions on the strong convexity of mirror map and the smoothness
of the objective $f_{i}$s are with respect to $\ell_{2}$ norm. Our assumptions
of the strong convexity and smoothness are for an arbitrary norm. Considering
that the mirror descent is a nonEuclidean variant of gradient descent, assume
a nonEuclidean norm would be more suitable.

3\. They assume that the $f_{i}$ has bounded variance as $E_{i}[\|\nabla
f_{i}(x)-\nabla f(x)\|_{2}^{2}]\leq\sigma^{2}$. This assumption is not
required in our work.

4\. They assume the PL inequality (w.r.t. the $\ell_{2}$ norm) holds to
achieve a linear convergence rate. We prove the same rate, but under Quadratic
Growth (QC) condition (w.r.t. an arbitrary norm). When considering the
$\ell_{2}$ norm, the PL condition implies the QC condition [21], thus their
assumption is stronger than our assumption here.

Difference in measuring convergence: Since the objective in [1] is nonconvex,
the convergence is evaluated for the stability gap $E[\|\nabla
f(x)\|_{2}^{2}]$. Our works applies to a convex objective, so the convergence
is evaluated for the optimization gap $E[f(x)]-\min_{x}f(x)$.

Comparing the SFO complexity: Unify all other notations and conditions in our
work and theirs, we compare the SFO complexity required to achieve a
$\epsilon$ stability gap, i.e. $E[\|\nabla f(x)\|_{2}^{2}]\leq\epsilon$:

For the general case, our SFO is
$\mathcal{O}(\frac{1}{\epsilon}+\frac{n}{\epsilon})$, and theirs is
$\mathcal{O}(\min(\frac{n}{\epsilon},\frac{\sigma^{2}}{\epsilon^{2}})+\frac{1}{\epsilon})$.
Ignoring the bounded variance assumption by taking $\sigma=\infty$, then their
rate is the same as ours;

For the special case where we have linear convergence rate, our SFO is
$\mathcal{O}(n\log\frac{1}{\epsilon}+\frac{1}{\mu}\log\frac{1}{\epsilon})$,
and theirs is
$\mathcal{O}(\min(n,\frac{\sigma^{2}}{\mu\epsilon})\frac{1}{\mu}\log\frac{1}{\epsilon}+\frac{1}{\mu}\log\frac{1}{\epsilon})$.
Again taking $\sigma=\infty$, our SFO complexity is better than theirs,
possibly because the difference in problem setting.

##  Appendix B Lemmas

We first list some lemmas that will be used in the proof. For each lemma, we
either give its reference or show the proof.

######  Lemma 1 ([22]).

For a $L$-smooth function $f$ with respect to $\|\cdot\|$, we have:

$f(\mathbf{u})-f(\mathbf{w})-\nabla
f(\mathbf{w})^{T}(\mathbf{u}-\mathbf{w})\leq\frac{L}{2}\|\mathbf{u}-\mathbf{w}\|^{2}.$
(20)

######  Lemma 2 (Theorem 2.1.5 in [23]).

For a $L$-smooth function $f$ with respect to $\|\cdot\|$, we have:

$\|\nabla f(\mathbf{u})-\nabla f(\mathbf{w})\|^{2}_{*}\leq
2L[f(\mathbf{u})-f(\mathbf{w})-\langle\nabla
f(\mathbf{w}),\mathbf{u}-\mathbf{w}\rangle].$ (21)

And especially, when $\nabla f(\mathbf{w})=0$:

$\|\nabla f(\mathbf{u})\|^{2}_{*}\leq 2L[f(\mathbf{u})-f(\mathbf{w})].$ (22)

######  Lemma 3.

For the problem
$\min_{\bm{\beta}}F(\bm{\beta})=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\bm{\beta})$,
suppose all $f_{i}$s are convex and $L-$smooth w.r.t. norm $\|\cdot\|$, denote
$\bm{\beta}^{*}=\operatorname{argmin}_{\bm{\beta}}F(\bm{\beta})$, then for $i$
chosen uniformly random from $\\{1,\ldots,n\\}$, the term $\mathbf{v}=\nabla
f_{i}(\bm{\beta})-\nabla f_{i}(\tilde{\bm{\beta}})+\nabla
F(\tilde{\bm{\beta}})$ have:

$\mathbbm{E}\|\mathbf{v}\|_{*}^{2}\leq
12L[F(\bm{\beta})-F(\bm{\beta}^{*})+F(\tilde{\bm{\beta}})-F(\bm{\beta}^{*})].$
(23)

###### Proof.

$\displaystyle\mathbbm{E}\|\mathbf{v}\|_{*}^{2}$
$\displaystyle=\mathbbm{E}\|\nabla f_{i}(\bm{\beta})-\nabla
f_{i}(\tilde{\bm{\beta}})+\nabla F(\tilde{\bm{\beta}})\|_{*}^{2}$
$\displaystyle\leq 3\mathbbm{E}\|\nabla f_{i}(\bm{\beta})-\nabla
f_{i}(\bm{\beta}^{*})\|_{*}^{2}+3\mathbbm{E}\|\nabla
f_{i}(\bm{\beta}^{*})-\nabla f_{i}(\tilde{\bm{\beta}})\|_{*}^{2}+3\|\nabla
F(\tilde{\bm{\beta}})\|_{*}^{2}$ $\displaystyle=3\mathbbm{E}\|\nabla
f_{i}(\bm{\beta})-\nabla f_{i}(\bm{\beta}^{*})\|_{*}^{2}+3\mathbbm{E}\|\nabla
f_{i}(\bm{\beta}^{*})-\nabla f_{i}(\tilde{\bm{\beta}})\|_{*}^{2}+3\|\nabla
F(\tilde{\bm{\beta}})-\nabla F(\bm{\beta}^{*})\|_{*}^{2}$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem0_0}}}{{\leq}}6\mathbbm{E}\\{L[f_{i}(\bm{\beta})-f_{i}(\bm{\beta}^{*})-\langle\nabla
f_{i}(\bm{\beta}^{*}),\bm{\beta}-\bm{\beta}^{*}\rangle]\\}+6\mathbbm{E}\\{L[f_{i}(\tilde{\bm{\beta}})-f_{i}(\bm{\beta}^{*})-\langle\nabla
f_{i}(\bm{\beta}^{*}),\tilde{\bm{\beta}}-\bm{\beta}^{*}\rangle]\\}$
$\displaystyle\quad+6\\{L[F(\tilde{\bm{\beta}})-F(\bm{\beta}^{*})-\langle\nabla
F(\bm{\beta}^{*}),\tilde{\bm{\beta}}-\bm{\beta}^{*}\rangle]\\}$
$\displaystyle=6L[F({\bm{\beta}})-F(\bm{\beta}^{*})]+12L[F(\tilde{\bm{\beta}})-F(\bm{\beta}^{*})]$
$\displaystyle\leq
12L[F({\bm{\beta}})-F(\bm{\beta}^{*})+F(\tilde{\bm{\beta}})-F(\bm{\beta}^{*})]$

∎

######  Lemma 4.

For the linear regression objective function
$F(\bm{\beta})=\|X\bm{\beta}-y\|_{2}^{2}/(2n)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\bm{\beta})$
where
$f_{i}(\bm{\beta})=\left(\mathbf{x}_{i}^{T}\bm{\beta}-y_{i}\right)^{2}/2$, the
term $\mathbf{v}=\nabla f_{i}(\bm{\beta})-\nabla
f_{i}(\tilde{\bm{\beta}})+\nabla F(\tilde{\bm{\beta}})$ is in
$\mathrm{col}(X^{T})$.

###### Proof.

Check

$\displaystyle\nabla
f_{i}(\bm{\beta})=\mathbf{x}_{i}(\mathbf{x}_{i}^{T}\bm{\beta}-y_{i})=X^{T}\begin{bmatrix}\mathbf{0}_{i-1}\\\
1\\\
\mathbf{0}_{n-i}\end{bmatrix}(\mathbf{x}_{i}^{T}\bm{\beta}-y_{i})\in\mathrm{col}(X^{T}).$

And $\nabla
F(\tilde{\bm{\beta}})=X^{T}\left[\frac{1}{n}(X\tilde{\bm{\beta}}-\mathbf{y})\right]\in\mathrm{col}(X^{T})$.
Thus $\mathbf{v}=\nabla f_{i}(\bm{\beta})-\nabla
f_{i}(\tilde{\bm{\beta}})+\nabla F(\tilde{\bm{\beta}})\in\mathrm{col}(X^{T})$.
∎

######  Lemma 5.

For a matrix $X\in\mathbbm{R}^{n\times p}$, for any
$\mathbf{u}\in\mathrm{col}(X^{T})$, we have

$\|X\mathbf{u}\|_{2}^{2}\geq s_{m}^{2}\|\mathbf{u}\|_{2}^{2}$ (24)

where $s_{m}^{2}$ is the smallest non-zero eigenvalue of $X^{T}X$.

###### Proof.

Consider the SVD of $X^{T}=Q\Sigma V$ where $Q^{T}Q=QQ^{T}=I_{p}$ and
$V^{T}V=VV^{T}=I_{n}$

$\Sigma=\begin{bmatrix}\operatorname{diag}(\mathbf{s})_{m}&\mathbf{0}_{m\times(n-m)}\\\
\mathbf{0}_{(p-m)\times m}&\mathbf{0}_{(p-m)\times(n-m)}\end{bmatrix}_{p\times
n}$

with $\mathbf{s}=[s_{1},\ldots,s_{m}]$, $s_{1}>\ldots>s_{m}>0$. And write
$Q=[\mathbf{q}_{1},\ldots,\mathbf{q}_{p}]$, then for any
$\mathbf{u}\in\mathrm{col}(X^{T})$, we can write
$\mathbf{u}=\sum_{i=1}^{m}w_{i}\mathbf{q}_{i}$ for some $w_{1},\ldots,w_{m}$,
and $\|\mathbf{u}\|_{2}^{2}=\sum_{i=1}^{m}w_{i}^{2}$.

Then

$\displaystyle\begin{split}X\mathbf{u}&=V^{T}\Sigma^{T}Q^{T}(\sum_{i=1}^{m}w_{i}\mathbf{q}_{i})\\\
&=V^{T}\Sigma^{T}\begin{bmatrix}w_{1}\\\ \vdots\\\ w_{m}\\\ 0\\\ \vdots\\\
0\end{bmatrix}\\\ &=V^{T}\begin{bmatrix}s_{1}w_{1}\\\ \vdots\\\ s_{m}w_{m}\\\
0\\\ \vdots\\\ 0\end{bmatrix}.\end{split}$

Thus

$\displaystyle\begin{split}\|X\mathbf{u}\|_{2}^{2}=\|V^{T}\begin{bmatrix}s_{1}w_{1}\\\
\vdots\\\ s_{m}w_{m}\\\ 0\\\ \vdots\\\
0\end{bmatrix}\|_{2}^{2}=\sum_{i=1}^{m}s_{i}^{2}w_{i}^{2}\geq
s_{m}^{2}\sum_{i=1}^{m}w_{i}^{2}=s_{m}^{2}\|\mathbf{u}\|_{2}^{2}.\end{split}$

∎

######  Lemma 6.

If $f$ is $\alpha-$strongly convex w.r.t. $\|\cdot\|$, then

$f(\mathbf{u})\leq f(\mathbf{w})+\nabla
f(\mathbf{w})^{T}(\mathbf{u}-\mathbf{w})+\frac{1}{2\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}.$ (25)

###### Proof.

Consider $g(\mathbf{v})=f(\mathbf{v})-\mathbf{v}^{T}\nabla f(\mathbf{w})$,
then $g$ is also $\alpha-$strongly convex, and it’s easy to check that $g$ has
a unique minimizer $\mathbf{w}$.

Next, by $\alpha-$strong convexity of $g$, for any $\mathbf{u}$ we have

$g(\mathbf{v})\geq g(\mathbf{u})+\nabla
g(\mathbf{u})^{T}(\mathbf{v}-\mathbf{u})+\frac{\alpha}{2}\|\mathbf{v}-\mathbf{u}\|^{2}.$
(26)

Minimizing both side of (26) w.r.t. $\mathbf{v}$, we have

$\displaystyle\begin{split}g(\mathbf{w})=\min_{\mathbf{v}}g(\mathbf{v})&\geq\min_{\mathbf{v}}g(\mathbf{u})+\nabla
g(\mathbf{u})^{T}(\mathbf{v}-\mathbf{u})+\frac{\alpha}{2}\|\mathbf{v}-\mathbf{u}\|^{2}\\\
&\geq\min_{\mathbf{v}}g(\mathbf{u})-\|\nabla
g(\mathbf{u})\|_{*}\|\mathbf{v}-\mathbf{u}\|+\frac{\alpha}{2}\|\mathbf{v}-\mathbf{u}\|^{2}\\\
&\geq g(\mathbf{u})-\frac{1}{2\alpha}\|\nabla
g(\mathbf{u})\|_{*}^{2}.\end{split}$

That is,

$f(\mathbf{w})-\mathbf{w}^{T}\nabla f(\mathbf{w})\geq
f(\mathbf{u})-\mathbf{u}^{T}\nabla f(\mathbf{w})-\frac{1}{2\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}.$

Rearrange the terms we get (25). ∎

######  Lemma 7.

If $f$ is $\alpha-$strongly convex w.r.t. $\|\cdot\|$, then

$[\nabla f(\mathbf{u})-\nabla
f(\mathbf{w})]^{T}[\mathbf{u}-\mathbf{w}]\leq\frac{1}{\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}.$ (27)

###### Proof.

Apply (25) twice to get:

$\displaystyle f(\mathbf{u})$ $\displaystyle\leq f(\mathbf{w})+\nabla
f(\mathbf{w})^{T}(\mathbf{u}-\mathbf{w})+\frac{1}{2\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2},$ (28) $\displaystyle
f(\mathbf{w})$ $\displaystyle\leq f(\mathbf{u})+\nabla
f(\mathbf{u})^{T}(\mathbf{w}-\mathbf{u})+\frac{1}{2\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}.$ (29)

Sum (28) and (29) we have:

$\displaystyle 0\leq(\nabla
f(\mathbf{w})-f(\mathbf{u}))^{T}(\mathbf{u}-\mathbf{w})+\frac{1}{\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}$ $\displaystyle\Longrightarrow$
$\displaystyle[\nabla f(\mathbf{u})-\nabla
f(\mathbf{w})]^{T}[\mathbf{u}-\mathbf{w}]\leq\frac{1}{\alpha}\|\nabla
f(\mathbf{u})-\nabla f(\mathbf{w})\|_{*}^{2}.$

∎

######  Lemma 8 ([21]).

If $f$ satisfies $PL$ inequality with respect to $\ell_{2}$ norm
$\|\cdot\|_{2}$, then $f$ satisfies the Quadratic Growth (QG) condition w.r.t.
$\|\cdot\|_{2}$ for the same constant $\mu$.

##  Appendix C Proof for Proposition 1

###### Proof.

Consider the Bregman divergence

$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{*})-\psi(\bm{\beta}_{t+1}^{s})-\nabla\psi(\bm{\beta}_{t+1}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{*})-\psi(\bm{\beta}_{t+1}^{s})-[\nabla\psi(\bm{\beta}_{t}^{s})-\eta\mathbf{v}_{t}^{s}]^{T}(\bm{\beta}^{*}-\bm{\beta}_{t}^{s}+\bm{\beta}_{t}^{s}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{*})-\psi(\bm{\beta}_{t}^{s})+\psi(\bm{\beta}_{t}^{s})-\psi(\bm{\beta}_{t+1}^{s})-\nabla\psi(\bm{\beta}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t}^{s})$
$\displaystyle+\nabla\psi(\bm{\beta}_{t}^{s})^{T}(\bm{\beta}_{t+1}^{s}-\bm{\beta}_{t}^{s})+(\eta\mathbf{v}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t}^{s})-D_{\psi}(\bm{\beta}_{t+1}^{s},\bm{\beta}_{t}^{s})+(\eta\mathbf{v}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t+1}^{s})].$
(30)

The last term in (30):

$\displaystyle\mathbbm{E}[(\mathbf{v}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[(\mathbf{v}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t}^{s}+\bm{\beta}_{t}^{s}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[\mathbbm{E}[(\mathbf{v}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t}^{s})|\bm{\beta}_{t}^{s}]]+\mathbbm{E}[(\mathbf{v}_{t}^{s})^{T}(\bm{\beta}_{t}^{s}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle=$ $\displaystyle\mathbbm{E}[(\nabla
F(\bm{\beta}_{t}^{s})^{T}(\bm{\beta}^{*}-\bm{\beta}_{t}^{s})]+\frac{1}{\eta}\mathbbm{E}[(\nabla\psi(\bm{\beta}_{t}^{s})-\nabla\psi(\bm{\beta}_{t+1}^{s}))^{T}(\bm{\beta}_{t}^{s}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}[(F(\bm{\beta}^{*})-F(\bm{\beta}_{t}^{s})]+\frac{1}{\eta}\mathbbm{E}[(\nabla\psi(\bm{\beta}_{t}^{s})-\nabla\psi(\bm{\beta}_{t+1}^{s}))^{T}(\bm{\beta}_{t}^{s}-\bm{\beta}_{t+1}^{s})]$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem6}}}{{\leq}}$
$\displaystyle\mathbbm{E}[(F(\bm{\beta}^{*})-F(\bm{\beta}_{t}^{s})]+\frac{1}{\eta\alpha}\mathbbm{E}[\|\nabla\psi(\bm{\beta}_{t}^{s})-\nabla\psi(\bm{\beta}_{t+1}^{s})\|_{*}^{2}]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[(F(\bm{\beta}^{*})-F(\bm{\beta}_{t}^{s})]+\frac{\eta}{\alpha}\mathbbm{E}[\|\mathbf{v}_{t}^{s}\|_{*}^{2}].$

Thus we have

$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t+1}^{s})]$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t}^{s})]-\mathbbm{E}[D_{\psi}(\bm{\beta}_{t+1}^{s},\bm{\beta}_{t}^{s})]+\eta\mathbbm{E}[F(\bm{\beta}^{*})-F(\bm{\beta}_{t}^{s})]+\frac{\eta^{2}}{\alpha}\mathbbm{E}[\|\mathbf{v}_{t}^{s}\|_{*}^{2}]$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem1}}}{{\leq}}$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t}^{s})]-\mathbbm{E}[D_{\psi}(\bm{\beta}_{t+1}^{s},\bm{\beta}_{t}^{s})]$
$\displaystyle\quad+(\eta-\frac{12L\eta^{2}}{\alpha})\mathbbm{E}[F(\bm{\beta}^{*})-F(\bm{\beta}_{t}^{s})]+\frac{12L\eta^{2}}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{t}^{s})]+\left(\frac{24L\eta^{2}}{\alpha}-\eta\right)\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\frac{12L\eta^{2}}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]-\frac{12L\eta^{2}}{\alpha}\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})].$
(31)

Sum (31) for $t=1,\ldots,m$, we have

$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{1}^{s})]+\sum_{t=1}^{m}\left(\frac{24L\eta^{2}}{\alpha}-\eta\right)\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]-\frac{12L\eta^{2}}{\alpha}\sum_{t=1}^{m}\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})].$
(32)

For option I, $\bm{\beta}_{1}^{s}=\bm{\beta}_{m+1}^{s-1}$,
$\tilde{\bm{\beta}}^{s}$ is a uniform random sample from
$\\{\bm{\beta}_{1}^{s},\ldots,\bm{\beta}_{m}^{s}\\}$, (32) becomes:

$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s-1})]+\sum_{t=1}^{m}\left(\frac{24L\eta^{2}}{\alpha}-\eta\right)\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]-\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})].$

That is,

$\displaystyle\left(\eta-\frac{24L\eta^{2}}{\alpha}\right)\sum_{t=1}^{m}\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s-1})]+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]$
$\displaystyle\quad-\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]-\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})].$
(33)

Thus we can define Lyapunov function

$\displaystyle
P^{s}:=\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})]\geq
0.$

Then

$\displaystyle\mathbbm{E}[F(\bm{\beta}_{a})-F(\bm{\beta}^{*})]$
$\displaystyle=\frac{1}{T}\sum_{s=1}^{S}\sum_{t=1}^{m}\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:thm2:eq5}}}{{\leq}}\frac{1}{T}\sum_{s=1}^{S}\frac{1}{\eta-\frac{24L\eta^{2}}{\alpha}}[P^{s-1}-P^{s}]$
$\displaystyle\leq\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}P^{0}$

For Option II, $\bm{\beta}_{1}^{s}=\tilde{\bm{\beta}}^{s-1}$, (32) becomes:

$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]\leq$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{s-1})]+\sum_{t=1}^{m}\left(\frac{24L\eta^{2}}{\alpha}-\eta\right)\mathbbm{E}[F(\bm{\beta}_{t}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]-\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle=$
$\displaystyle\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{s-1})]+m\left(\frac{12L\eta^{2}}{\alpha}-\eta\right)\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\frac{12L\eta^{2}m}{\alpha}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})].$
(34)

Since $\psi$ is $\ell-$smooth and $F$ satisfies QG condition, $\bm{\beta}^{*}$
is any minimum point of $F$, we can take $\bm{\beta}^{*}$ as required by Lemma
8 to get

$\displaystyle D_{\psi}(\bm{\beta}^{*},\tilde{\bm{\beta}}^{s-1})$
$\displaystyle=\psi(\bm{\beta}^{*})-\psi(\tilde{\bm{\beta}}^{s-1})-\nabla\psi(\tilde{\bm{\beta}}^{s-1})^{T}(\bm{\beta}^{*}-\tilde{\bm{\beta}}^{s-1})$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lsmooth_quad_upper}}}{{\leq}}\frac{\ell}{2}\left\|\bm{\beta}^{*}-\tilde{\bm{\beta}}^{s-1}\right\|^{2}$
$\displaystyle\stackrel{{\scriptstyle QG\text{ of
}f}}{{\leq}}\frac{\ell}{\mu}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})].$
(35)

Plug (35) into (34) we have

$\displaystyle
0\leq\mathbbm{E}[D_{\psi}(\bm{\beta}^{*},\bm{\beta}_{m+1}^{s})]\leq$
$\displaystyle\left(\frac{12L\eta^{2}m}{\alpha}-m\eta\right)\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\quad+\left(\frac{12L\eta^{2}m}{\alpha}+\frac{\ell}{\mu}\right)\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})].$

That is,

$\displaystyle\mathbbm{E}[F(\tilde{\bm{\beta}}^{s})-F(\bm{\beta}^{*})]$
$\displaystyle\leq\frac{\frac{12L\eta^{2}}{\alpha}+\frac{\ell}{m\mu}}{\eta-\frac{12L\eta^{2}}{\alpha}}\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})]$
$\displaystyle=\tau\mathbbm{E}[F(\tilde{\bm{\beta}}^{s-1})-F(\bm{\beta}^{*})].$

Thus,

$\displaystyle\mathbbm{E}[F(\tilde{\bm{\beta}}^{S})-F(\bm{\beta}^{*})]$
$\displaystyle\leq\tau^{S}\mathbbm{E}[F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{*})]$

and by $\bm{\beta}^{a}=\tilde{\bm{\beta}}^{S}$ we get the desired bound.

∎

##  Appendix D Proof for Theorem 1

###### Proof.

We first check that every $f_{i}$ is $L-$smooth w.r.t. $\|\cdot\|_{2}$.

$\displaystyle\|\nabla f_{i}(\mathbf{u})-\nabla f_{i}(\mathbf{w})\|_{2}$
$\displaystyle=\|(\mathbf{x}_{i}^{T}\mathbf{u}-y_{i})\mathbf{x}_{i}-(\mathbf{x}_{i}^{T}\mathbf{w}-y_{i})\mathbf{x}_{i}\|_{2}$
$\displaystyle=\|\mathbf{x}_{i}^{T}(\mathbf{u}-\mathbf{w})\mathbf{x}_{i}\|_{2}$
$\displaystyle=|\mathbf{x}_{i}^{T}(\mathbf{u}-\mathbf{w})|*\|\mathbf{x}_{i}\|_{2}$
$\displaystyle\leq\|\mathbf{x}_{i}\|_{2}\|\mathbf{u}-\mathbf{w}\|_{2}\|\mathbf{x}_{i}\|_{2}$
$\displaystyle\leq(\max_{i}\|x_{i}\|_{2}^{2})\|\mathbf{u}-\mathbf{w}\|_{2}.$

First, we prove part (a) of Theorem 1. When $\psi(\bm{\beta})$ is
$\alpha-$strongly convex w.r.t. $\|\cdot\|_{2}$, we can apply part (a) of
Theorem 1, take $\bm{\beta}^{*}=\bm{\beta}^{\psi}$ to get

$\displaystyle\mathbbm{E}\left[\frac{1}{2n}\|X\bm{\beta}^{a}-\mathbf{y}\|_{2}^{2}-\frac{1}{2n}\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right]$
$\displaystyle\leq$
$\displaystyle\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\left[D_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\frac{1}{2n}\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\frac{1}{2n}\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right].$
(36)

For the L.H.S. of (36) we have

$\displaystyle\|X\bm{\beta}^{a}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}+X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}+\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}+2\langle
X\bm{\beta}^{a}-X\bm{\beta}^{\psi},P_{\mathrm{col}(X)}\mathbf{y}-\mathbf{y}\rangle-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}+2\langle
X(\bm{\beta}^{a}-\bm{\beta}^{\psi}),-P_{\mathcal{N}(X^{T})}\mathbf{y}\rangle$
$\displaystyle=$
$\displaystyle\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}=\|X\bm{\beta}^{a}-P_{\mathrm{col}(X)}\mathbf{y}\|_{2}^{2}.$

Thus (36) becomes

$\displaystyle\mathbbm{E}[\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}]=\mathbbm{E}[\|X\bm{\beta}^{a}-P_{\mathrm{col}(X)}\mathbf{y}\|_{2}^{2}]$
$\displaystyle\leq$
$\displaystyle\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right].$
(37)

Since
$\nabla\psi(\tilde{\bm{\beta}}^{0})=\nabla\psi(\bm{\beta}_{m}^{0})\in\mathrm{col}(X^{T})$,
by Lemma 4 we will have $\nabla\psi(\bm{\beta}_{t}^{s})\in\mathrm{col}(X^{T})$
for all $s$ and $t$, thus

$\nabla\psi(\bm{\beta}^{a})\in\mathrm{col}(X^{T}).$ (38)

Then

$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}\langle\nabla\psi(\bm{\beta}^{a}),\bm{\beta}^{a}-\bm{\beta}^{\psi}\rangle$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:16}}}{{=}}$
$\displaystyle\mathbbm{E}\langle\nabla\psi(\bm{\beta}^{a}),P_{\mathrm{col}(X^{T})}(\bm{\beta}^{a}-\bm{\beta}^{\psi})\rangle$
$\displaystyle\leq$
$\displaystyle\mathbbm{E}[\|\nabla\psi(\bm{\beta}^{a})\|_{2}\|P_{\mathrm{col}(X^{T})}(\bm{\beta}^{a}-\bm{\beta}^{\psi})\|_{2}]$
$\displaystyle\stackrel{{\scriptstyle}}{{\leq}}$ $\displaystyle
B*\mathbbm{E}\|P_{\mathrm{col}(X^{T})}(\bm{\beta}^{a}-\bm{\beta}^{\psi})\|_{2}$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem4}}}{{\leq}}$
$\displaystyle\frac{B}{s_{m}}\mathbbm{E}\|XP_{\mathrm{col}(X^{T})}(\bm{\beta}^{a}-\bm{\beta}^{\psi})\|_{2}$
$\displaystyle=$
$\displaystyle\frac{B}{s_{m}}\mathbbm{E}\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}$
$\displaystyle\leq$
$\displaystyle\frac{B}{s_{m}}(\mathbbm{E}\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2})^{1/2}$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:15}}}{{\leq}}$
$\displaystyle\frac{B}{s_{m}}\sqrt{\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right]}.$

Thus the $\bm{\beta}^{a}$ from SMD will have

$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]$
$\displaystyle\quad\leq\frac{B}{s_{m}}\sqrt{\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right]}$
(39)
$\displaystyle\mathbbm{E}[\|X\bm{\beta}^{a}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}]$
$\displaystyle\quad\leq\frac{\alpha}{(\alpha\eta-24L\eta^{2})T}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right].$
(40)

Set the R.H.S. of (39) and (40) to $\epsilon$ to get an $\epsilon-$solution,
we can solve

$\displaystyle T$
$\displaystyle\geq\frac{B^{2}}{s_{m}^{2}\epsilon^{2}}\frac{\alpha}{(\alpha\eta-24L\eta^{2})}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right]$
$\displaystyle T$
$\displaystyle\geq\frac{1}{\epsilon}\frac{\alpha}{(\alpha\eta-24L\eta^{2})}\left[2nD_{\psi}(\bm{\beta}^{\psi},\tilde{\bm{\beta}}^{0})+\frac{12L\eta^{2}m}{\alpha}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)\right].$

That is, $T\sim\mathcal{O}(\frac{1}{\epsilon}+\frac{1}{\epsilon^{2}})$.

For part (b) of Theorem 1, we first need to show that $F(\cdot)$ satisfies PL
inequality w.r.t. $\|\cdot\|_{2}$. Remember
$F(\bm{\beta})=\frac{1}{2n}\|X\bm{\beta}-y\|_{2}^{2}$, then

$\displaystyle\frac{1}{2}\|\nabla F(\bm{\beta})\|_{2}^{2}$
$\displaystyle=\frac{1}{2}\|\frac{1}{n}X^{T}(X\bm{\beta}-y)\|_{2}^{2}$
$\displaystyle=\frac{1}{2n^{2}}\|X^{T}(X\bm{\beta}-P_{\mathrm{col}(X)}y)\|_{2}^{2}$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem4}}}{{\geq}}\frac{s_{m}^{2}}{2n^{2}}\|X\bm{\beta}-P_{\mathrm{col}(X)}y\|_{2}^{2}$
$\displaystyle=\frac{s_{m}^{2}}{2n^{2}}\|X\bm{\beta}-y+P_{\mathcal{N}(X^{T})}y\|_{2}^{2}$
$\displaystyle=\frac{s_{m}^{2}}{2n^{2}}[\|X\bm{\beta}-y\|_{2}^{2}-\|P_{\mathcal{N}(X^{T})}y\|_{2}^{2}]$
$\displaystyle=\frac{s_{m}^{2}}{n}[F(\bm{\beta})-F(\bm{\beta}^{\psi})].$

Thus $F(\cdot)$ satisfies PL inequality with constant $\frac{s_{m}^{2}}{n}$.
By Lemma 8, $F(\cdot)$ satisfies QG with the same constant, so part (b) of
Theorem 1 applies here, and we get

$\mathbbm{E}[F(\bm{\beta}_{a})-F(\bm{\beta}^{\psi})]\leq(\tau^{\prime})^{S}[F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})].$

where

$\tau^{\prime}=\frac{\frac{12L\eta^{2}}{\alpha}+\frac{\ell
n}{ms_{m}^{2}}}{\eta-\frac{12L\eta^{2}}{\alpha}},$

By the similar statement as the first part of the proof, we have

$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]$
$\displaystyle\leq$
$\displaystyle\frac{B}{s_{m}}[2n\mathbbm{E}[F(\bm{\beta}^{a})-F(\bm{\beta}^{\psi})]]^{.5}$
$\displaystyle\leq$
$\displaystyle\frac{B\sqrt{2n}(\tau^{\prime})^{S/2}}{s_{m}}[F(\tilde{\bm{\beta}}^{0})-F(\bm{\beta}^{\psi})]^{.5}.$

Thus

$\displaystyle\mathbbm{E}[\psi(\bm{\beta}^{a})-\psi(\bm{\beta}^{\psi})]$
$\displaystyle\leq\frac{B(\tau^{\prime})^{S/2}}{s_{m}}\sqrt{\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}}$
$\displaystyle\mathbbm{E}[\|X\bm{\beta}^{a}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}]$
$\displaystyle\leq(\tau^{\prime})^{S}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right).$

∎

##  Appendix E Proof of Corollary 1

###### Proof.

For the specific choice of $\psi(\cdot)=\frac{1}{2}\|\cdot\|_{2}^{2}$, we have

$\nabla^{2}\psi(\cdot)=I.$

That is, $\psi(\cdot)$ is $1-$smooth and $1-$strongly convex, for

$\tau^{\prime\prime}=\frac{12L\eta^{2}+\frac{n}{ms_{m}^{2}}}{\eta-12L\eta^{2}}<1,$

we can apply part (b) of Theorem 2 to get

$\displaystyle\mathbbm{E}[\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}$
$\displaystyle\leq(\tau^{\prime\prime})^{S}\left(\|X\tilde{\bm{\beta}}^{0}-\mathbf{y}\|_{2}^{2}-\|X\bm{\beta}^{\psi}-\mathbf{y}\|_{2}^{2}\right)$
$\displaystyle=(\tau^{\prime\prime})^{S}\left(\|\mathbf{y}\|_{2}^{2}-\|P_{\mathcal{N}(X^{T})}\mathbf{y}\|_{2}^{2}\right)=(\tau^{\prime\prime})^{S}\|P_{\mathrm{col}(X)}\mathbf{y}\|_{2}^{2}.$

Then since $\nabla\psi(\bm{\beta}^{a})=\bm{\beta}^{a}\in\mathrm{col}(X^{T})$
and
$\bm{\beta}^{\psi}=(X^{T}X)^{+}X^{T}\mathbf{y}=X^{T}(XX^{T})^{+}\mathbf{y}\in\mathrm{col}(X^{T})$,
we have

$\mathbbm{E}\|\bm{\beta}_{a}-\bm{\beta}^{\psi}\|_{2}^{2}\stackrel{{\scriptstyle\eqref{ch3:eq:lem4}}}{{\leq}}\frac{1}{s_{m}^{2}}\mathbbm{E}\|X\bm{\beta}^{a}-X\bm{\beta}^{\psi}\|_{2}^{2}\leq\frac{\tau^{S}}{s_{m}^{2}}\|P_{\mathrm{col}(X)}\mathbf{y}\|_{2}^{2}.$

∎

##  Appendix F Proof of Theorem 2

To prove Theorem 2, we will need a few preliminaries.

######  Lemma 9.

Denote
$\bm{\beta}^{(0)}=\operatorname{argmin}_{\bm{\beta}}\\{\|\bm{\beta}\|_{1}:X\bm{\beta}=P_{\mathrm{col}(X)}\mathbf{y}\\}$
and
$\bm{\beta}^{(\delta)}=\operatorname{argmin}_{\bm{\beta}}\\{\|\bm{\beta}\|_{1+\delta}^{1+\delta}:X\bm{\beta}=P_{\mathrm{col}(X)}\mathbf{y}\\}$,
we have

$\|\bm{\beta}^{(\delta)}\|_{1}\leq
p^{\frac{\delta}{1+\delta}}\|\bm{\beta}^{(0)}\|_{1}$ (41)

###### Proof.

$\displaystyle\|\bm{\beta}^{(\delta)}\|_{1}$
$\displaystyle=\sum_{i=1}^{p}1*|\beta_{i}^{(\delta)}|$
$\displaystyle\leq\left(\sum_{i=1}^{p}1^{\frac{1+\delta}{\delta}}\right)^{\frac{\delta}{1+\delta}}*\left(\sum_{i=1}^{p}|\beta_{i}^{(\delta)}|^{1+\delta}\right)^{\frac{1}{1+\delta}}$
$\displaystyle=p^{\frac{\delta}{1+\delta}}\|\bm{\beta}^{(\delta)}\|_{1+\delta}$
$\displaystyle\leq p^{\frac{\delta}{1+\delta}}\|\bm{\beta}^{(0)}\|_{1+\delta}$
$\displaystyle\leq
p^{\frac{\delta}{1+\delta}}\left(\sum_{i=1}^{p}\|\bm{\beta}^{(0)}\circ\mathbf{e}_{i}\|_{1+\delta}\right)$
$\displaystyle=p^{\frac{\delta}{1+\delta}}\left(\sum_{i=1}^{p}|\beta^{(0)}_{i}|\right)$
$\displaystyle=p^{\frac{\delta}{1+\delta}}\|\bm{\beta}^{(0)}\|_{1}$

where the first inequality holds by Hölder’s inequality, the second inequality
follows from $\bm{\beta}^{(\delta)}$ minimizes
$\|\bm{\beta}\|_{1+\delta}^{1+\delta}$ thus minimizes
$\|\bm{\beta}\|_{1+\delta}$ among all
$\bm{\beta}\in\mathcal{B}:=\\{\bm{\beta}:X\bm{\beta}=P_{\mathrm{col(X)}}\mathbf{y}\\}$,
and the third inequality follows from triangle inequality. ∎

######  Lemma 10.

Assume $\mathbf{y}=X\bm{\beta}^{o}$ where $\bm{\beta}^{o}\leq s$ and $X$ is
$s-$good with parameter $\kappa<1/2$. Use the notation in Lemma 9 that
$\bm{\beta}^{(0)}=\operatorname{argmin}_{\bm{\beta}}\\{\|\bm{\beta}\|_{1}:X\bm{\beta}=\mathbf{y}\\}$,
we have $\bm{\beta}^{(0)}=\bm{\beta}^{o}$.

###### Proof.

$\bm{\beta}^{o}$ is a feasible solution to
$\min_{\bm{\beta}}\\{\|\bm{\beta}\|_{1}:X\bm{\beta}=\mathbf{y}\\}$. Denote
$\hat{\mathbf{u}}=\bm{\beta}^{o}-\bm{\beta}^{(0)}$, we have

$\displaystyle\left(\|\bm{\beta}^{(0)}_{S}\|_{1}+\|\bm{\beta}^{(0)}_{S^{c}}\|_{1}\right)=\|\bm{\beta}^{(0)}\|_{1}\leq\|\bm{\beta}^{o}\|_{1}=\|\bm{\beta}^{o}_{S}\|_{1}$
(42) $\displaystyle\Longrightarrow$
$\displaystyle\|\bm{\beta}^{(0)}_{S^{c}}\|_{1}\leq\|\bm{\beta}^{o}_{S}\|_{1}-\|\bm{\beta}^{(0)}_{S}\|_{1}\leq\|\bm{\beta}^{o}_{S}-\bm{\beta}^{(0)}_{S}\|_{1}$
(43) $\displaystyle\Longleftrightarrow$
$\displaystyle\|\hat{\mathbf{u}}_{S^{c}}\|_{1}\leq\|\hat{\mathbf{u}}_{S}\|_{1}.$
(44)

On the other hand,

$\displaystyle
X\hat{\mathbf{u}}=X\bm{\beta}^{o}-X\bm{\beta}^{(0)}=\mathbf{y}-\mathbf{y}=\mathbf{0}.$

That is, $\hat{\mathbf{u}}\in\mathcal{N}(X)$, then by $X$ is $s-$good, we have
for $\kappa<\frac{1}{2}$

$\displaystyle\|\hat{\mathbf{u}}_{S}\|_{1}\leq\kappa\|\hat{\mathbf{u}}\|_{1}$
$\displaystyle\Longrightarrow$
$\displaystyle\|\hat{\mathbf{u}}_{S}\|_{1}\leq\frac{\kappa}{1-\kappa}\|\hat{\mathbf{u}}_{S^{c}}\|_{1}.$
(45)

Combine (42) and (F) we get

$\displaystyle\|\hat{\mathbf{u}}_{S^{c}}\|_{1}\leq\frac{\kappa}{1-\kappa}\|\hat{\mathbf{u}}_{S^{c}}\|_{1}$
$\displaystyle\Longrightarrow$
$\displaystyle\|\hat{\mathbf{u}}_{S^{c}}\|_{1}=0$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem10eq3}}}{{\Longrightarrow}}$
$\displaystyle\|\hat{\mathbf{u}}_{S}\|_{1}=0.$

Thus $\hat{\mathbf{u}}=\mathbf{0}$, $\bm{\beta}^{(0)}=\bm{\beta}^{o}$. ∎

######  Lemma 11.

Assume $X\in\mathbbm{R}^{n\times p}$ is $s-good$ with parameter $\kappa<1/2$,
then $\forall\bm{\beta}\in\mathbbm{R}^{p}$ and $I\subset\\{1,\ldots,p\\}$ with
$|I|\leq s$ we have

$\|\bm{\beta}_{I}\|_{1}\leq\frac{\sqrt{s}+\kappa\sqrt{p}}{s_{m}}\|X\bm{\beta}\|_{2}+\kappa\|\bm{\beta}\|_{1}$

where $s_{m}^{2}$ is the smallest nonzero eigenvalue of $X^{T}X$.

###### Proof.

$\displaystyle\|\bm{\beta}_{I}\|_{1}$
$\displaystyle\leq\|(P_{\mathrm{col}(X^{T})}\bm{\beta})_{I}\|_{1}+\|(P_{\mathcal{N}(X)}\bm{\beta})_{I}\|_{1}$
$\displaystyle\leq\sqrt{s}\|(P_{\mathrm{col}(X^{T})}\bm{\beta})_{I}\|_{2}+\kappa\|P_{\mathcal{N}(X)}\bm{\beta}\|_{1}$
$\displaystyle\leq\sqrt{s}\|P_{\mathrm{col}(X^{T})}\bm{\beta}\|_{2}+\kappa(\|P_{\mathcal{N}(X)}\bm{\beta}+P_{\mathrm{col}(X^{T})}\bm{\beta}\|_{1}+\|P_{\mathrm{col}(X^{T})}\bm{\beta}\|_{1})$
$\displaystyle\stackrel{{\scriptstyle\eqref{ch3:eq:lem4}}}{{\leq}}\frac{\sqrt{s}}{s_{m}}\|XP_{\mathrm{col}(X^{T})}\bm{\beta}\|_{2}+\kappa\left(\|\bm{\beta}\|_{1}+\frac{\sqrt{p}}{s_{m}}\|XP_{\mathrm{col}(X^{T})}\bm{\beta}\|_{2}\right)$
$\displaystyle=\frac{\sqrt{s}+\kappa\sqrt{p}}{s_{m}}\|X\bm{\beta}\|_{2}+\kappa\|\bm{\beta}\|_{1}$

where the first inequality holds by triangle inequality, the second inequality
holds by Cauchy-Schwartz inequality and s-goodness of $X$, the third
inequality again uses the triangle inequality. ∎

######  Theorem 3.

[Theorem 1.3.1 in [24]] Consider the following optimization problem with any
norm $\|\cdot\|$

$\displaystyle\begin{split}\hat{\bm{\beta}}=\arg\min_{\bm{\beta}}&\quad\|\bm{\beta}\|_{1}\\\
s.t.&\quad\|X\bm{\beta}-\mathbf{y}\|\leq\Delta\end{split}$ (46)

Assume $\exists\alpha\geq 0$ and $\kappa^{\prime}<1/2$ that
$\forall\bm{\beta}\in\mathbbm{R}^{p}$ and $\forall I\subset\\{1,\ldots,p\\}$
with $|I|\leq s$:

$\|\bm{\beta}_{I}\|_{1}\leq\alpha\|X\bm{\beta}\|+\kappa^{\prime}\|\bm{\beta}\|_{1}$

Suppose $\tilde{\bm{\beta}}$ is an approximate solution to (46) such that

$\displaystyle\|\tilde{\bm{\beta}}\|_{1}$
$\displaystyle\leq\|\hat{\bm{\beta}}\|_{1}+\nu$
$\displaystyle\|X\tilde{\bm{\beta}}-\mathbf{y}\|$
$\displaystyle\leq\Delta+\epsilon.$

Then for any nearly $s-$sparse and feasible $\bm{\beta}^{\prime}$ for (46),
where near $s-$sparsity means $\exists$ $s-$sparse $\bm{\beta^{s}}$ such that
$\|\bm{\beta}^{\prime}-\bm{\beta^{s}}\|_{1}\leq v$, we have

$\|\tilde{\bm{\beta}}-\bm{\beta}^{\prime}\|_{1}\leq\frac{2\alpha(2\Delta+\epsilon)+2v+\nu}{1-2\kappa^{\prime}}.$

We are ready to prove Theorem 2.

###### Proof for Theorem 2.

By Lemma 9 we have

$\displaystyle\|\bm{\beta}^{(\delta)}\|_{1}\leq\|\bm{\beta}^{(0)}\|_{1}+\left(p^{\frac{\delta}{1+\delta}}-1\right)\|\bm{\beta}^{(0)}\|_{1}$

And by Lemma 11

$\|\bm{\beta}_{I}\|_{1}\leq\frac{\sqrt{s}+\kappa\sqrt{p}}{s_{m}}\|X\bm{\beta}\|_{2}+\kappa\|\bm{\beta}\|_{1}$

Thus we can apply Theorem 3 by setting $\Delta=0$, and
$\|\cdot\|=\|\cdot\|_{2}$, then $\bm{\hat{\beta}}=\bm{\beta}^{(0)}$. Take
$\tilde{\bm{\beta}}=\bm{\beta}^{(\delta)}$,
$\bm{\beta}^{\prime}=\bm{\beta}^{o}$, then
$\nu=\left(p^{\frac{\delta}{1+\delta}}-1\right)\|\bm{\beta}^{(0)}\|_{1}$,
$\epsilon=0$ , $v=0$, $\alpha=\frac{\sqrt{s}+\kappa\sqrt{p}}{s_{m}}$,
$\kappa^{\prime}=\kappa$, then we have

$\|\bm{\beta}^{(\delta)}-\bm{\beta}^{o}\|_{1}\leq\frac{\left(p^{\frac{\delta}{1+\delta}}-1\right)\|\bm{\beta}^{(0)}\|_{1}}{1-2\kappa}=\frac{\left(p^{\frac{\delta}{1+\delta}}-1\right)\|\bm{\beta}^{o}\|_{1}}{1-2\kappa}$
(47)

where the equality is from Lemma 10.

We can further bound $\|\bm{\beta}^{o}\|_{1}$. By $(s,\gamma)-RE$ condition we
have

$\displaystyle\|\bm{\beta}^{o}\|^{2}_{2}\leq\frac{\|X\bm{\beta}^{o}\|^{2}_{2}}{n\gamma}=\frac{\|\mathbf{y}\|_{2}^{2}}{n\gamma}.$

Thus

$\displaystyle\|\bm{\beta}^{o}\|_{1}\leq\sqrt{s}\|\bm{\beta}^{o}\|_{2}\leq\sqrt{\frac{s}{n\gamma}}\|\bm{y}\|_{2}.$
(48)

Plug (48) into (47), we have

$\displaystyle\|\bm{\beta}^{(\delta)}-\bm{\beta}^{o}\|_{1}$
$\displaystyle\leq\frac{\left(p^{\frac{\delta}{1+\delta}}-1\right)}{1-2\kappa}\sqrt{\frac{s}{n\gamma}}\|\bm{y}\|_{2}$
$\displaystyle\leq\frac{\left(p^{\frac{\log\left(1+\frac{(1-2\kappa)\sqrt{n\gamma}}{\sqrt{s}\|\mathbf{y}\|_{2}}\xi\right)}{\log
p}}-1\right)}{1-2\kappa}\sqrt{\frac{s}{n\gamma}}\|\bm{y}\|_{2}$
$\displaystyle=\frac{\left(1+\frac{(1-2\kappa)\sqrt{n\gamma}}{\sqrt{s}\|\mathbf{y}\|_{2}}\xi-1\right)}{1-2\kappa}\sqrt{\frac{s}{n\gamma}}\|\bm{y}\|_{2}=\xi.$

∎


# Gravitational Radiation from Primordial Helical MHD Turbulence

Tina Kahniashvili  CCPP, New York University, 4 Washington Plaza, New York, NY
10003, USA  Department of Physics, Kansas State University, 116 Cardwell Hall,
Manhattan, KS 66506, USA  Department of Physics, Laurentian University, Ramsey
Lake Road, Sudbury, ON P3E 2C6, Canada  National Abastumani Astrophysical
Observatory, 2A Kazbegi Ave, Tbilisi, GE-0160, Georgia  Grigol Gogoberidze
Centre for Plasma Astrophysics, K.U. Leuven, Celestijnenlaan 200B, 3001
Leuven, Belgium  Department of Physics, Kansas State University, 116 Cardwell
Hall, Manhattan, KS 66506, USA  National Abastumani Astrophysical Observatory,
2A Kazbegi Ave, Tbilisi, GE-0160, Georgia  Bharat Ratra  Department of
Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, USA

( KSUPT-08/1)

###### Abstract

We consider gravitational waves (GWs) generated by primordial inverse-cascade
helical magneto-hydrodynamical (MHD) turbulence produced by bubble collisions
at the electroweak phase transitions (EWPT). Compared to the unmagnetized EWPT
case, the spectrum of MHD-turbulence-generated GWs peaks at lower frequency
with larger amplitude and can be detected by the proposed Laser Interferometer
Space Antenna (LISA).

###### pacs:

98.70.Vc, 98.80.-k

When detected, primordial cosmological GWs will provide a very valuable probe
of the very early Universe sources . Various mechanisms that generate such GWs
have been discussed: quantum fluctuations inflation ; bubble wall motion and
collisions during phase transitions bubble ; kos1 ; cosmic strings strings ;
cosmological magnetic fields magnet ; cdk04 ; cd06 ; and plasma turbulence
kmk02 ; dolgov ; kgr05 ; gkk07 . From the direct detection point of view, GWs
generated during the EWPT are promising since their peak frequency lies in or
near the LISA lisa  frequency band m00 , however, to produce a detectable
signal the EWPT must be strong enough detection ; gs06 ; grojean . Currently
discussed EWPT models do not predict an observable GW signal from bubble
collisions grojean , nor for GWs produced by unmagnetized turbulence gkk07 .

Here we study the generation of GWs during a first-order EWPT assuming that
bubble collisions produce helical MHD turbulence.111Kinetic or magnetic
helicity generation at the EWPT is studied in Refs. helicity . Previously we
studied generation of GWs by direct-cascade turbulence and found that, due to
parity violation in the early Universe, the induced GWs are circularly
polarized kgr05 . Polarized GWs are present in other models gw-pol , and the
polarization of the GW background is in principle observable, either directly
seto  or through the CMB cdk04 ; a . In the case of unmagnetized
hydrodynamical turbulence the peak frequency of the GW power spectrum is
determined by the inverse turn-over time of the largest eddy and the energy-
scale when the GW is generated. Recently discussed modifications of the
standard EWPT model place the transition at a higher energy-scale gs06 . As a
result, the GW power spectrum peak frequency is shifted to higher frequency
which, since the GW spectrum is sharply peaked, reduces the possibility of
detection by LISA. On the other hand, in the case of MHD turbulence the
presence of an energy inverse-cascade leads to an increase in the effective
size of the largest eddy (now associated with an helical magnetic field), and
can result in the GW power spectrum peaking in the LISA band, with amplitude
large enough to be detected by LISA. We adapt the technique developed in Ref.
gkk07  to study this case here. We model MHD turbulence and obtain the GW
spectrum by using an analogy with the theory of sound wave production by
hydrodynamical turbulence L52 ; P52 ; G ; my75 .

Since the turbulent fluctuations are stochastic, so are the generated GWs. The
GW energy density is m00

$\displaystyle\rho_{\rm GW}({\bf x})$ $\displaystyle=$
$\displaystyle\frac{1}{32\pi
G}\langle\partial_{t}h_{ij}({\mathbf{x}},t)\partial_{t}h_{ij}({\mathbf{x}},t)\rangle$
$\displaystyle=$ $\displaystyle\frac{G}{2\pi}\int{\rm d}^{3}{\bf
x}^{\prime}{\rm d}^{3}{\bf
x}^{\prime\prime}\frac{\langle\partial_{t}S_{ij}({\mathbf{x}}^{\prime},t^{\prime})\partial_{t}S_{ij}({\mathbf{x}}^{\prime\prime},t^{\prime\prime})\rangle}{|{\bf
x}-{\bf x}^{\prime}||{\bf x}-{\bf x}^{\prime\prime}|}.$

Here the times $t^{\prime(\prime\prime)}=t-|{\bf x}-{\bf
x}^{\prime(\prime\prime)}|$, $i$ and $j$ are spatial indices (repeated indices
are summed), the source
$S_{ij}({\mathbf{x}},t)=T_{ij}({\mathbf{x}},t)-\delta_{ij}T^{k}_{k}({\mathbf{x}},t)/3$
is the traceless part of the stress-energy tensor $T_{ij}$, $G$ is the
gravitational constant, and we use natural units with $\hbar=1=c$. We assume
that the turbulence exists for a time short enough to neglect the cosmological
expansion during GW production. We consider metric perturbations in the far-
field limit (i.e. for $x\gg d$, where $d$ is a characteristic length-scale of
the source region), where GWs are the only metric perturbations W , and
replace $|{\bf x}-{\bf x}^{\prime}|$ by $|\bf x|$ in Eq. (Gravitational
Radiation from Primordial Helical MHD Turbulence). If the turbulence is
stationary then the GW spectral energy density $I({\bf x},\omega)$ ($\rho_{\rm
GW}({\bf x})=\int{\rm d}\omega I({\bf x},\omega)$ where $\omega$ is the
angular frequency) is gkk07

$I({\bf x},\omega)=\frac{4\pi^{2}\omega^{2}G{\rm w}^{2}}{|{\bf
x}|^{2}}\int{\rm d}^{3}{\bf x}^{\prime}H_{ijij}\left({\bf
x}^{\prime},\frac{{\bf x}}{|{\bf x}|}\omega,\omega\right).$ (2)

Here $H_{ijij}({\bf x}^{\prime},{\bf k},\omega)$ (where ${\bf k}$ is a proper
wavevector) is the (double traced) four-dimensional Fourier transform of the
two-point time-delayed fourth-order correlation tensor, $\langle S_{ij}({\bf
x}^{\prime},t)S_{lm}({\bf x}^{\prime\prime},t+\tau)\rangle/{\rm w}^{2}$, with
respect to ${\bf x}^{\prime\prime}-{\bf x}^{\prime}$ and $\tau$, where ${\rm
w}=\rho+p$ is the enthalpy density and $p$ and $\rho$ the pressure and energy
density of the plasma.

We assume that primordial MHD turbulence is generated at time $t_{\star}$ at
proper length-scale $l_{0}=2\pi/k_{0}$ with characteristic velocity
perturbation $v_{0}$.222We assume that the usual and magnetic Reynolds numbers
are much greater than unity on scales $\sim l_{0}$, otherwise there is no
turbulence. Throughout this paper the symbol $\sim$ represents equality to the
accuracy of a dimensionless multiplicative constant of order unity. The
dynamics of MHD turbulence is dominated by Alfvén waves for which the magnetic
and kinetic energy densities are in approximate equipartition B03 . In this
case $v_{0}\sim b_{0}$, where $b_{0}=B_{0}/\sqrt{4\pi{\rm w}}$ is the
characteristic magnetic field perturbation expressed in velocity units. While
MHD turbulence is isotropic on large scales, it is locally anisotropic on
small scales s83 , resulting in small-scale anisotropy in the generated GW
background. However, GWs are generated mainly by the largest eddies gkk07  so
we adopt an isotropic turbulence model, and thus the magnetic field two-point
correlation function is $\langle b_{i}^{\star}({\bf k},t)b_{j}({\bf
k^{\prime}},t+\tau)\rangle=F_{ij}^{M}\\!({\bf k},t)f(\eta(k),\tau)\delta({\bf
k}-{\bf k^{\prime}})$, with my75

$F_{ij}^{M}\\!({\bf k},\tau)=P_{ij}({\bf k})\frac{E^{M}\\!(k,t)}{4\pi
k^{2}}+i\varepsilon_{ijl}{k_{l}}\frac{H^{M}\\!(k,t)}{8\pi k^{2}}.$ (3)

Here $P_{ij}({\bf k})=\delta_{ij}-{k_{i}k_{j}}/{k^{2}}$, and $E^{M}\\!(k,t)$
and $H^{M}\\!(k,t)$ are the magnetic field energy and helicity densities. The
Schawarz inequality implies $|H^{M}\\!(k,t)|\leq 2E^{M}\\!(k,t)/k$ B03 . For
the total magnetic energy ${\mathcal{E}}_{M}(t)=\int E^{M}\\!(k,t)dk$ and
helicity ${\mathcal{H}}_{M}(t)=\int H^{M}\\!(k,t)dk$ we get
${\mathcal{H}}_{M}(t)\leq 2\xi_{M}(t){\mathcal{E}}_{M}(t)$, where
$\xi_{M}(t)\equiv\int E^{M}(k,t)k^{-1}dk/{\mathcal{E}}_{M}(t)$ is the
magnetic-eddy correlation length. $\eta(k)$ is an autocorrelation function
that determines the characteristic function $f(\eta(k),\tau)$ that describes
the temporal decorrelation of turbulent fluctuations. In the following we use
$f(\eta(k),\tau)=\exp\left(-{\pi}\eta^{2}(k)\tau^{2}/4\right)$ K64 .

After generation primordial turbulence freely decays. We adopt the decaying
MHD turbulence model of Refs. BM99 ; CHB05 . For non-zero initial magnetic
helicity turbulence decay is a two stage process. First decay stage dynamics
is governed by a direct cascade of energy density lasting for a time
$\tau_{s0}=s_{0}\tau_{0}$, several times ($s_{0}\sim 3-5$) longer then the
characteristic largest-eddy turn-over time
$\tau_{0}=l_{0}/v_{0}=2\pi/k_{0}v_{0}$. During the first stage energy density
flows from large to small scales and finally dissipates on scales $\sim
l_{d}=2\pi/k_{d}$ ($k_{d}\gg k_{0}$) where one of the Reynolds numbers becomes
of order unity. Due to the selective decay effect B03  magnetic helicity is
nearly conserved during this stage CHB05 . To compute the GWs generated by
decaying MHD turbulence, we assume that decaying turbulence lasting for time
$\tau_{s0}$ is equivalent to stationary turbulence lasting for time
$\tau_{s0}/2$. This can be justified using the Proudman P52 ; my75  argument
for (unmagnetized) hydrodynamical turbulence. Consequently, when computing the
emitted GWs we ignore the time dependence of $E^{M}(k,t)$ and $H^{M}(k,t)$. We
also assume small initial magnetic helicity,
$\alpha_{\star}\equiv{\mathcal{H}}_{M}(t_{\star})/[2\xi_{M}(t_{\star}){\mathcal{E}}_{M}(t_{\star})]\ll
1$. For $E^{M}(k,t)$ and $\eta(k)$ we use the Kolmogorov model,

$E^{M}\\!(k,t)\sim\varepsilon^{2/3}k^{-5/3},~{}~{}~{}\eta(k)=\varepsilon^{1/3}k^{2/3}/\sqrt{2\pi},$
(4)

for $k_{0}<k<k_{d}$. Here $\varepsilon\sim k_{0}v_{0}^{3}$ is the energy
dissipation rate per unit enthalpy.

At the end of the first stage turbulence relaxes to a maximally helical state,
$\alpha_{s0}\sim 1$ CHB05 ; jedamzik . Accounting for conservation of magnetic
helicity, the characteristic velocity and magnetic field perturbations at this
stage are $v_{1}\sim\alpha_{\star}^{1/2}v_{0}$ and
$b_{1}\sim\alpha_{\star}^{1/2}b_{0}$. Second stage dynamics is governed by a
magnetic helicity inverse cascade. If both Reynolds numbers are large at the
end of the first stage, magnetic helicity is conserved during the second
stage. The magnetic eddy correlation length evolves as $\xi_{M}(t)\sim
l_{0}\sqrt{1+t/\tau_{1}}$ BM99 ; CHB05  where $\tau_{1}\sim
l_{0}/v_{1}=\tau_{0}/\sqrt{\alpha_{\star}}$ is the characteristic energy
containing eddy turn-over time at the beginning of the second stage. The
magnetic ${\mathcal{E}}_{M}(t)$ and kinetic ${\mathcal{E}}_{K}(t)$ energy
densities evolve as BM99 ; CHB05

$\displaystyle{\mathcal{E}}_{M}(t)$ $\displaystyle\propto$
$\displaystyle(1+t/\tau_{1})^{-1/2}{\rm w}b_{1}^{2},$
$\displaystyle{\mathcal{E}}_{K}(t)$ $\displaystyle\propto$
$\displaystyle(1+t/\tau_{1})^{-1}{\rm w}v_{1}^{2}.$ (5)

These imply that the characteristic turn-over ($\tau_{\rm to}$) and cascade
($\tau_{\rm cas}$) timescales evolve as

$\displaystyle\tau_{\rm to}\sim\tau_{\rm cas}\sim\tau_{1}(1+t/\tau_{1}).$ (6)

To compute the GWs emitted during the second stage we use the stationary
turbulence model that has the same GW output. Introducing the characteristic
wavenumber $k_{\xi}(t)=2\pi/\xi_{M}(t)$ and using Eqs. (5) we find
${\mathcal{E}}_{M}\sim{\rm w}v_{1}^{2}k_{\xi}(t)/k_{0}$ and
${\mathcal{E}}_{K}\sim{\rm w}v_{1}^{2}[k_{\xi}(t)/k_{0}]^{2}$ since $b_{1}\sim
v_{1}$. The time when turbulence is present on scale $\xi_{M}(t)$ is
determined by Eq. (6) which can be rewritten as $\tau_{\rm
cas}\sim\tau_{1}[k_{0}/k_{\xi}(t)]^{2}$. So instead of considering decaying
turbulence, we consider stationary turbulence with a scale-dependent duration
time (time during which the magnetic energy is present at the scale),
$\tau_{s1}\sim\tau_{1}[k_{0}/k]^{2}$ (for $k=k_{\xi}$ this coincides with
$\tau_{\rm cas}$).

The expression for ${\mathcal{E}}_{M}$ yields the time-independent

$E^{M}\\!(k,t)=C_{1}v_{1}^{2}/k_{0}=kH^{M}\\!(k,t)/2,~{}~{}~{}k_{S}<k<k_{0}.$
(7)

Here $C_{1}$ is a constant of order unity, $k_{S}$ is the smallest wavenumber
where the inverse cascade stops, and the second equation follows from
saturating the causality condition. For the second stage autocorrelation
function, which is inversely proportional to the turn-over time (6), we assume
$\eta(k)=(k/k_{0})^{2}/\sqrt{2\pi}\tau_{1}$. At the largest scales there is no
efficient dissipation mechanism, so the inverse cascade will be stopped at
scale $l_{S}(t)=2\pi/k_{S}$ where either the cascade timescale $\tau_{\rm
cas}$ reaches the expansion timescale $H^{-1}_{\star}=H^{-1}(t_{\star})$, or
when the characteristic length scale $\xi_{M}(t)\sim l_{S}$ reaches the Hubble
radius. These conditions are $\alpha_{\star}^{-1/2}l_{S}^{2}/v_{0}l_{0}\leq
H^{-1}_{\star}$ or $l_{S}\leq H^{-1}_{\star}$ (the cascade time is scale
dependent and maximal at $k=k_{S}$). Defining $\gamma=l_{0}/H^{-1}_{\star}$ (
$\gamma\leq 1$), it is easy to see that the first condition is fulfilled first
and consequently $k_{0}/k_{S}\leq(v_{0}/\gamma)^{1/2}\alpha_{\star}^{1/4}$. To
have an inverse cascade requires $k_{0}/k_{S}\geq 1$, leading to a constraint
on initial helicity, $\gamma\leq M\alpha_{\star}^{1/2}$ (where $M=v_{0}$ is
the turbulence Mach number).

The magnetic field perturbation stress-energy tensor is
$T_{ij}^{M}({\mathbf{x}},t)={\rm
w}b_{i}({\mathbf{x}},t)b_{j}({\mathbf{x}},t)$. For the first decay stage we
compute for this magnetic part and then double the result to account for
approximate magnetic and kinetic energy equipartition for Alfvén waves. During
the second stage, according to Eqs. (5), kinetic energy can be neglected
compared to magnetic energy. To compute $H_{ijij}({\bf k},\omega)$ we assume
Millionshchikov quasi-normality my75  and use the convolution theorem (for
details see Sec. III of Ref. gkk07 ). Using the ($k\rightarrow 0$) aero-
acoustic approximation, which is accurate for low Mach number ($M\leq 1$, (and
slightly overestimates GWs amplitude for the Mach number approaching unity
($M\rightarrow 1$) gkk07 ), we find

$\displaystyle\\!\\!\\!\\!\\!\\!\\!H_{ijij}({\bf k},\omega)\approx
H_{ijij}(0,\omega)=~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (8)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\frac{7C_{1}^{2}M^{3}\alpha_{\star}^{3/2}\\!}{6\pi^{3/2}k_{0}}\\!\\!\int_{k_{S}}^{k_{0}}\\!\\!\frac{{\rm
d}k}{k^{4}}\exp\\!\left(-\frac{\omega^{2}k_{0}^{2}}{\alpha_{\star}M^{2}k^{4}}\right)\\!{\rm
erfc}\\!\left(-\frac{\omega k_{0}}{\sqrt{\alpha_{\star}}Mk^{2}}\right).$

The integral is dominated by the contribution of large scale ($k\simeq k_{S}$)
perturbations and is maximal at $\omega_{\rm
max}^{(II)}\sim\alpha_{\star}^{1/2}Mk_{S}^{2}/k_{0}=2\pi H_{\star}$. For the
first-stage direct-cascade turbulence the peak frequency is $\omega_{\rm
max}^{(I)}\sim k_{0}M$ gkk07 . To determine the peak frequency at the current
epoch we need to account for the cosmological expansion which decreases the GW
amplitude and frequency by the factor $a_{\star}/a_{0}$, where $a_{\star}$ and
$a_{0}$ are the values of the cosmological scale factor at the GW generation
and current epochs

The total GW energy spectrum at a given space-time event is obtained by
integrating over all source regions with a light-like separation from that
event, and includes contributions from GW generated during the first and
second stages. For the first stage (with duration time
$\tau^{(I)}_{T}=s_{0}\tau_{0}$) $\rho_{GW}^{(I)}(\omega)$ is given by Eqs.
(21) and (A3) of Ref. gkk07 . For the second stage contribution we must
account for the scale dependence of the cascade time. The total GW fractional
energy density parameter at the moment of emission $\Omega_{{\rm GW},\star}$
is
$105~{}H_{\star}^{4}\omega^{3}\sum_{m}\tau_{T}^{(m)}H_{ijij}^{(m)}(0,\omega_{\star})/H_{0}^{2}$
gkk07 . Here the index $m$ runs over $I$ and $II$ for the first and second
decay stages, $H_{0}$ $\omega_{\star}$ is an angular frequency at the moment
of emission. The current GW amplitude is related to the current fractional
energy density parameter through $h_{C}(f)=1.26\times 10^{-18}\left({1\,{\rm
Hz}}/{f}\right)\left[h_{0}^{2}\Omega_{GW}(f)\right]^{1/2}$ (where $h_{0}$ is
the current Hubble parameter $H_{0}$ in units of 100 ${\rm km}\,{\rm
sec}^{-1}{\rm Mpc}^{-1}$) m00 , and

$\displaystyle h_{C}(f)$ $\displaystyle\simeq$ $\displaystyle 2\times
10^{-14}\left(\frac{100\,{\rm
GeV}}{T_{*}}\right)\left(\frac{100}{g_{*}}\right)^{1/3}$ (9)
$\displaystyle\times\sum_{m}\left[\tau_{T}^{(m)}\omega_{\star}H_{\star}^{4}H_{ijij}^{(m)}(0,\omega_{\star})\right]^{1/2}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

Here the linear frequency $f=(a_{\star}/a_{0})f_{\star}$ with
$f_{\star}=\omega_{\star}/2\pi$, $T_{\star}$ and $g_{\star}$ are the
temperature and effective number of relativistic (all fields) degrees of
freedom at scale factor $a_{\star}$.

Figure 1: The spectrum of gravitational radiation from MHD turbulence for
$g_{\ast}=100,~{}T_{\ast}=100~{}{\rm GeV},~{}\gamma=0.01$, and $M=1/\sqrt{3}$,
for four different initial magnetic helicity values, $\alpha_{\star}=0$ (solid
line), $\alpha_{\star}=0.02$ (dashed line), $\alpha_{\star}=0.05$ (dash-dotted
line), and $\alpha_{\star}=0.1$ (dotted line). The bold line is the LISA
design sensitivity curve.

Figure 1 shows $h_{C}(f)$ for a few initial magnetic helicity values. GWs
emitted during direct-cascade unmagnetized turbulence peak at current $f_{\rm
max}^{(I)}\simeq M\nu_{\star}$ gkk07 . We find that the MHD-inverse-cascade
generated GW (current epoch) peak frequency is determined by cosmology
parameters, $f_{\rm max}^{(II)}=H_{\star}a_{0}/a_{\star}=1.6\times
10^{-5}\,{\rm Hz}\,\left({g_{*}}/{100}\right)^{1/6}\left({T_{*}}/{100\,{\rm
GeV}}\right)$ and is independent of turbulence parameters. On the other hand,
$f_{\rm max}^{(II)}=\gamma f_{\rm max}^{(I)}/M$ is shifted to lower frequency
compared to the unmagnetized case. From Eq. (9), the amplitude of MHD-
turbulence-generated GWs at the peak is a factor
$\sim\alpha_{\star}^{9/8}\gamma^{-3/4}M^{3/4}$ larger than that in the
unmagnetized case.

When modeling turbulence we used the Biskamp and Muller model, BM99 ; CHB05 .
If we adopt the helical MHD turbulence model of Banerjee and Jedamzik jedamzik
(also see Refs. campanelli ) the GW peak frequency remains the same while the
amplitude of the signal doubles.

Figure 1 shows that even for small values of magnetic helicity the main
contribution to the GW energy density is from the second, inverse-cascade
stage. The GWs will be strongly polarized since magnetic helicity is maximal
at the end of the first stage kgr05 . LISA should be able to detect such GW
polarization seto . Unlike the unmagnetized case due to the second (inverse-
cascade) stage contribution the GW amplitude is large enough at $10^{-4}$ Hz
to be detectable by LISA. If the EWPT occurs at higher energies
($T_{\star}>100$ Gev) the peak is shifted to higher frequency, closer to LISA
sensitivity peak, which leads to a stronger signal. Our formalism is
applicable for GW production at an earlier QCD phase transition, assuming the
presence of colored magnetic fields QCD , or for any other phase transitions
bubble ; the peak frequency will be shifted according to the changes in
$T_{\star}$ and $g_{\star}$. The GW signal estimated here exceeds that from
bubble collisions kos1 ; detection ; gs06  or from hydrodynamical,
unmagnetized turbulence kmk02 ; dolgov ; gkk07 . Of course, this strong signal
assumes initial non-zero (although small) magnetic helicity, so detection of
polarized GWs by LISA will indicate parity violation during the EWPT as
proposed in Refs. helicity .

We greatly appreciate useful comments from A. Gruzinov and the referees. We
acknowledge helpful discussions with R. Durrer, G. Gabadadze, and A. Kosowsky.
T. K. acknowledges the hospitality of the Abdus Salam International Center for
Theoretical Physics. G. G. and T. K. acknowledge partial support from INTAS
061000017-9258 and Georgian NSF ST06/4-096 grants. B. R. acknowledges US DOE
grant DE-FG03-99EP4103.

## References

  * (1) S. Chongchitnan and G. Efstathiou, Phys. Rev. D 73, 083511 (2006); C. J. Hogan, AIP Conf. Proc. 873, 30 (2006); A. Buonanno, arXiv:0709.4682 [qr-qc].
  * (2) L. P. Grishchuk, Sov. Phys. JETP 40, 409 (1975) [Zh. Eksp. Teor. Fiz. 67, 825 (1974)]; A. Starobinsky, Sov. Phys. JETP Lett. 30, 682 (1979) [Pisma Zh. Eksp. Teor. Fiz. 30, 719 (1979)]; V. A. Rubakov, M.V. Sazhin, and A. V. Veryaskin, Phys. Lett. B 115, 189 (1982); B. Allen, Phys. Rev. D 37, 2078 (1988); B. Ratra, Phys. Rev. D 45 1913 (1992); M. Giovannini, Phys. Rev. D 60, 123511 (1999).
  * (3) E. Witten, Phys. Rev. D 30, 272 (1984); C. J. Hogan, Mon. Not. R. Astron. Soc. 218, 629 (1986); M. S. Turner and F. Wilczek, Phys. Rev. Lett. 65, 3080 (1990); A. Kosowsky, M. S. Turner, and R. Watkins, Phys. Rev. D 45, 4514 (1992).
  * (4) M. Kamionkowski, A. Kosowsky, and M. S. Turner, Phys. Rev. D 49, 2837 (1994); R. Apreda, et al., Class. Quant. Grav. 18, L155 (2001); J. F. Dufaux, et al., Phys. Rev. D 76, 123517 (2007); S. J. Huber and T. Konstandin, arXiv:0709.2091 [hep-ph]; C. Caprini, R. Durrer, and G. Servant, arXiv:0711.2593 [astro-ph]; R. Easther, J. T. Giblin, and E. A. Lim, arXiv:0712.2991 [astro-ph]; J. Garcia-Bellido, D. G. Figueroa and A. Sastre, Phys. Rev. D 77, 043517 (2008); A. Megevand, arXiv:0804.0391 [astro-ph]. .
  * (5) T. Vachaspati and A. Vilenkin, Phys. Rev. D 31, 3052 (1985); M. P. Infante and N. Sánchez, Phys. Rev. D 61, 083515 (2000).
  * (6) D. V. Deriagin, et al., Mon. Not. R. Astron. Soc. 229, 357 (1987); R. Durrer, P. Ferreira, and T. Kahniashvili, Phys. Rev. D 61, 043001 (2000); A. Lewis, Phys. Rev. D. 70, 043011 (2004).
  * (7) C. Caprini, R. Durrer, and T. Kahniashvili, Phys. Rev. D 69, 063006 (2004).
  * (8) C. Caprini and R. Durrer, Phys. Rev. D 74, 063521 (2006).
  * (9) A. Kosowsky, A. Mack, and T. Kahniashvili, Phys. Rev. D 66, 024030 (2002).
  * (10) A. D. Dolgov, D. Grasso, and A. Nicolis, Phys. Rev. D 66, 103505 (2002).
  * (11) T. Kahniashvili, G. Gogoberidze, and B. Ratra, Phys. Rev. Lett. 95, 151301 (2005).
  * (12) G. Gogoberidze, T. Kahniashvili, and A. Kosowsky, Phys. Rev. D 76, 083002 (2007).
  * (13) http://lisa.nasa.gov/
  * (14) M. Maggiore, Phys. Rept. 331, 28 (2000).
  * (15) R. Apreda, et al., Nucl. Phys. B 631, 342 (2002); A. Nicolis, Class. Quantum Grav. 21, L27 (2004).
  * (16) C. Grojean and G. Servant, Phys. Rev. D 75 043507 (2007).
  * (17) C. Delaunay, C. Grojean, and J. D. Wells, arXiv:0711.2511 [hep-ph].
  * (18) J. Cornwall, Phys. Rev. D 56, 6146 (1997); M. Giovannini and M. E. Shaposhnikov, Phys. Rev. D 57, 2186 (1998); G. B. Field and S. M. Carroll, Phys. Rev. D 62, 103008 (2000); T. Vachaspati, Phys. Rev. Lett. 87, 251302 (2001); G. Sigl, Phys. Rev. D 66, 123002 (2002); K. Subramanian and A. Brandenburg, Phys. Rev. Lett. 93, 205001 (2004); L. Campanelli and M. Giannotti, Phys. Rev. D 72, 123001 (2005); V. B. Semikoz and D. D. Sokoloff, Astron. Astrophys. 413, L53 (2005); A. Diaz-Gil, et al., arXiv:0712.4263 [hep-ph]; L. Campanelli, arXiv:0805.0575 [astro-ph].
  * (19) A. Lue, L. M. Wang, and M. Kamionkowski, Phys. Rev. Lett. 83, 1506 (1999); D. Lyth, C. Quimbay, and Y. Rodriguez, JHEP 0503, 016 (2005); M. Satoh, S. Kanno, and J. Soda, Phys. Rev. D 77, 023526 (2008).
  * (20) N. Seto, Phys. Rev. Lett. 97, 151101 (2006).
  * (21) S. H. S. Alexander, M. E. Peskin, and M. M. Sheikh-Jabbari, Phys. Rev. Lett. 96, 081301 (2006); S. Saito, K. Ichiki, and A. Taruya, J. Cosmol. Astropart. Phys. 0709, 002 (2007).
  * (22) A. S. Monin and A. M. Yaglom, Statistical Fluid Mechanics (MIT Press, Cambridge, MA, 1975).
  * (23) M. E. Goldstein, Aeroacoustics (McGraw-Hill, New York, 1976).
  * (24) M. J. Lighthill, Proc. R. Soc. London Ser. A 211, 564 (1952), 222, 1 (1954).
  * (25) I. Proudman, Proc. R. Soc. London Ser. A 214, 119 (1952).
  * (26) S. Weinberg, Gravitation and Cosmology (Wiley & Sons, New York, 1972).
  * (27) D. Biskamp, Magnetohydrodynamic Turbulence (Cambridge University, Cambridge, 2003); M. K. Verma, Phys. Rept. 401, 229 (2004).
  * (28) J. V. Shebalin, W. H. Matthaeus, and D. Montgomery, J. Plasma Phys. 29, 525 (1983); P. Goldreich and S. Sridhar, Astrophys. J. 438, 763 (1995).
  * (29) R. H. Kraichnan, Phys. Fluids 7, 1163 (1964).
  * (30) D. Biskamp and W.-C. Mueller, Phys. Rev. Lett. 83, 2195 (1999), Phys. Plasma 7, 4889 (2000).
  * (31) M. Christensson, M. Hindmarsh, and A. Brandenburg, Phys. Rev. E 64, 056405 (2001), Astron. Nachrichten 326, 393 (2005).
  * (32) R. Banerjee and K. Jedamzik, Phys. Rev. D 70, 123003 (2004).
  * (33) D. T. Son, Phys. Rev. D 59, 063008 (1999); L. Campanelli, Phys. Rev. Lett. 98, 251302 (2007).
  * (34) M. H. P. van Putten, Phys. Rev. D 50, 6640 (1994).