It is well known that the class of rotation invariant
algorithms are suboptimal even for learning sparse linear problems
when the number of examples is below the “dimension” of the problem.
This class includes any gradient descent trained
neural net with a fully-connected input layer
(initialized with a rotationally symmetric distribution).
The simplest sparse problem is learning a single feature
out of $d$ features. In that case the classification error or regression loss
grows with $1-\sfrac kd$ where $k$ is the number of examples seen.
These lower bounds become vacuous when the number
of examples $k$ reaches the dimension $d$.

Nevertheless we show that when noise is added to this sparse linear problem, rotation invariant algorithms are still suboptimal after seeing $d$ or more examples.
We prove this via a lower bound for the Bayes optimal algorithm on a rotationally symmetrized problem.
We then prove much lower upper bounds on the same problem for simple non-rotation invariant algorithms.

Finally we analyze the gradient flow trajectories of many
standard optimization algorithms in some simple cases and show how
they veer toward or away from the sparse targets.

We believe that our trajectory categorization will be useful in designing algorithms
that can exploit sparse targets and
our method for proving lower bounds will be crucial for
analyzing other families of algorithms that admit different classes of invariances.

rotation invariance,
feed forward nets,
lower bounds,
multiplicative updates,

Any gradient descent trained neural net with a fully-connected input layer is rotation invariant when initialized with a rotationally symmetric distribution.
The reason is that if the input instances are rotated, then
weights rotate in the same way and therefore the dot products
(logits) feeding into the non-linearities of the first layer remain unchanged.
It is well known that rotation invariant
algorithms are fundamentally inferior when learning certain
sparse linear problems [Warmuth and Vishwanathan, 2005, Ng, 2004, Li et al., 2021, Warmuth et al., 2021].
These bounds are surprising because they hold for such a
general class of algorithms and are complemented by
simple non-rotation invariant algorithms that need
exponentially fewer examples to achieve the same accuracy.
However, these lower bounds are not practically relevant
because they all assume that the number of
examples is less than the dimension of the problem.

The simplest sparse linear problem is learning one out
of $d$ features. When the instances are for example the $d$
orthogonal rows of a Hadamard matrix and a rotation
invariant algorithm is given $k\le d$ training instances
[Learning is harder for sampling with replacement.
However lower bounds are stronger for sampling without replacement
and the boundary of getting a full rank instance set is clearer.
We present bounds for both types of models.]
then its generalization error on all $d$ instances is at
least $1-\sfrac kd$. There are classification [Ng, 2004, Li et al., 2021]
and regression versions of this problem [Warmuth and Vishwanathan, 2005, Warmuth et al., 2021]
but the lower bound is essentially the same, proven with
different techniques.
However, the lower bounds become vacuous when $k\ge d$, which is the case in most relevant settings.
For example linear regression (which is rotation invariant)
finds the unique consistent feature after seeing a
noise-free training set of full rank.
In contrast this paper provides lower bounds on the loss
of any rotation invariant algorithm
when the sample size (without replacement) exceeds the
input dimension of the problem.
Now any weight vector (including sparse ones) can be expressed as a linear combination of the instances.
The key point of our paper is that, when the noise is added to the
sparse target, rotation invariant algorithms still produce poor solutions compared to simple non-invariant algorithms.

Our lower bound technique creates a Bayesian setup
where the learning is presented with a randomly rotated
version of the input instances. We prove
a lower bound on the Bayes-optimal algorithm for this case,
and we prove any rotational invariant algorithm on the
original problem can do no better than the Bayes-optimal
algorithm on the randomly rotated problem.
[This proof technique is interesting in its own right
and is different from the regression techniques used for the
under-constrained case [Warmuth et al., 2021].]

We then show that there are trivial non-rotationally invariant algorithms that can learn
noisy sparse linear much more efficiently: multiplicative updates on a linear neuron or
gradient descent on a `spindly' two-layer linear net
where every linear weight $w_i$ is replaced by a product of parameter $u_iv_i$

for tree=circle,draw,s sep=11pt, scale=.40,l=2

[ [ ] ]
[ [ ] ]
[  ,edge label=node[midway,left=-1mm] $u_i$
[ ,edge label=node[midway,left=-1mm] $v_i$
] ]
[  [ ] ]
[  [ ] ]

(Figure <ref>).

Finally we also prove the same upper bound for a novel “priming method” that was conjectured to
learn sparse linear problems in
[Warmuth and Amid, 2023]: First find the linear least squares solution $\w$
and then after reweighing the $i$th feature with $w_i$,
apply tuned ridge regression.

We are not interested in linear problems per se,
but any fancy non-linear model should at least be able to handle the linear case.
Our lower bound implies that for example
deep neural nets with a complete input layer
that are trained with gradient descent and initialized with a rotation invariant distribution
cannot sample-efficiently learn noisy sparse linear
problems even with more examples than the input dimension.
This lower bound can be circumvented
by pivoting away from gradient descent to the alternative updates listed above,
which are of course not rotation invariant.
We also show that other popular rotationally non-invariant algorithms such as Adagrad and Adam
have an opposite inductive bias such that they learn sparse problems even less efficiently than gradient descent.
We visualize how the performance gap between rotationally invariant
and non-invariant algorithms arises by computing their trajectories in weight space
over the course of learning.
Previous work has focused on how rotational symmetry can be
broken adding a L1 penalty to the loss [Ng, 2004, Li et al., 2021].
This changes the minimum to which the learning algorithms converge to.
Here we focus on symmetry breaking due to the differences
in the iterative learning algorithm,
while keeping the loss function fixed (as in Li et al., 2021).
So all algorithms must converge to the same minimum but follow trajectories
that differ in how close they get to the sparse target.
The optimal early stopping point is the trajectory point that is closest to the target.

We observe that when the number of examples is larger than
the input dimension and the target is sparse linear, then multiplicative updates,
spindly networks, and priming all produce trajectories that pass extremely close to the sparse target before finally veering away
and learning the noise.
Rotation invariant algorithms cannot produce such a bias, even when the input is unbalanced (anisotropic covariance) provided it is drawn from a rotationally symmetric distribution.
In essence, they learn the signal and the noise at the same rate.
Interestingly, adaptive learning rate algorithms such as Adagrad [Duchi et al., 2011] and Adam [Kingma and Ba, 2015] show the opposite bias, producing trajectories that are curved away from sparse solutions and learn the noise before the signal.

We arrive at these observations by deriving closed-form solutions for the
continuous time (gradient flow) weight trajectories of each
algorithm for sparse noisy regression problems. We found
forms for the isotropic covariance case but also visualize the
trajectories of the anisotropic case.
We show how the behavior of the different algorithms
can be understood from a geometric analysis that recasts
mirror descent and other preconditioned gradient methods as
gradient descent under an altered geometry,
where the preconditioner acts as a Riemannian metric.

Finally we make some preliminary experimental observations for nonlinear models
by contrasting gradient descent training on a simple
three layer feed forward neural net when each logit on the first layer is connected to all inputs
(fully connected) versus connected to all inputs via
the spindly network (Figure <ref>).
Our experiments show in case b, gradient descent
produces sparse solutions for the input layer, is much
less confused by additional noisy features and can learn
quickly from informative features.

§.§.§ Additional related work:

There is a long history of contrasting the generalization
ability of additive versus multiplicative updates (see e.g.
[Kivinen et al., 1997, Kivinen and Warmuth, 1997]).
Surprisingly there is a connection between both update families
rooted in the observation that when the weights are products of parameters,
then gradient descent is biased towards sparsity [Gunasekar et al., 2017, Kerekes et al., 2021].
More precisely [Amid and Warmuth, 2020], the gradient flow of multiplicative updates on a
linear neuron equals the gradient flow of
gradient descent on the spindly network of Figure <ref>.

In this paper we prove lower bounds for sparse noisy regression problems.
It might also be possible to obtain lower bounds for the classification setting
when the number of examples exceeds the VC dimension, building on the setup of [Ng, 2004, Li et al., 2021].
However for technical reasons the additional loss in the VC dimension / Rademacher complexity
has a square root term that seems to make this
bounding methodology non-optimal in the noisy case.
Our upper bounds for the spindly network of Figure
<ref> have also been proven using the R.I.P.
assumption [Vaskevicius et al., 2019]. For completeness we
gave a self-contained proof in the appendix that also
works for multiplicative updates and their approximations.

Note that we have hardness results for a fixed noisy sparse linear problem:
First feature plus noise.
We prove that the gap in performance is due to rotation invariance.
So we avoid using staircase, cork screw or cryptographically hard functions [Abbe and Boix-Adsera, 2022]
for proving lower bounds. Also in [Shamir, 2018],
the target class is much more complicated and the hard input distribution is algorithm specific
and therefore no fixed problem is given that is hard for all algorithms
in the invariance class.

Paper [4] has two main contributions: In Section 3 they prove bounds for weak learnability
which means outperforming a random predictor by some small constant.
In contrast we prove bounds away from optimal performance.
In Section 4 they prove hardness results based on another cryptographic assumption:
The hardness of some XOR problems with Gaussian noise.
Their theorem says that certain derived problems cannot be learned
on a fully connected neural net with GD in polynomial time
with access to a sample oracle.

Both hardness results in [3] and [4] don't seem to have contrasting upper bounds
for efficient algorithms that solve the same hard problems as
we do for our simple noisy sparse linear problem. Actually
the algorithm is still gradient descent, but on the spindly network.

THANK YOU for the opportunity to express the uniqueness of our contribution.
We will add as much as possible of this discussion to our final version!

§ THE LOWER BOUND METHOD

§.§ Rotation invariance and problem setup

An example $(\x,y)$ is a $d$-dimensional vector, followed by a real-valued label
$y \in \mathbb{R}$. We specify a training set as a tuple $(\underset{n,
d}{\X},\underset{n}{\y})$ containing $n$ training examples, where the rows of
input matrix $\X$ are the $n$ (transposed) input vectors and the
target $\y$ is a vector of their labels.

A learning algorithm is a mapping which,
given the training set $(\X,\y)$, produces a real-valued prediction function
$\mathbb{R}^d \ni \x \mapsto \yh(\x | \X,\y) \in \mathbb{R}$.
An algorithm is called rotation invariant
[Warmuth and Vishwanathan, 2005, Warmuth et al., 2021] if for any orthogonal matrix $\underset{d,d}{\U}$ and any input $\x\in\mathbb{R}^d$:

\begin{equation}
\yh(\U \x\, | \X \U^\top,\y) = \yh(\x\, | \X ,\y)
\label{eq:def_rot_inv}
\end{equation}

In other words, the prediction $\yh(\x | \X,\y)$ remains the same if we rotate both $\x$ and all examples from $\X$ by the same orthogonal matrix $\U$.
If the algorithm is randomized, based on an internal random variable $Z$,[For example, in neural networks $Z$ would correspond to a random initialization of the parameter vector.]
then $\yh(\x_{\rm{te}} | \X, \y)$ is a random variable given by some function $f_{\x_{\rm{te}},\X,\y}(Z)$,
and the equality sign in equation (<ref>)
should be interpreted as “identically distributed”.

Our lower bounds in sections <ref>-<ref>
hold for any rotation invariant algorithm.
In particular, Warmuth et al., 2021 have shown that any neural network with a fully-connected input layer (and arbitrary remaining layers), in which the weights in the input layer are initialized randomly with a rotation invariant distribution (e.g. i.i.d. Gaussians) and are trained by gradient descent, is rotation invariant and is thus subject to our lower bound.
The reason is that such a network can be written as $f(\w_1\cdot\x,\w_2\cdot\x,\ldots,\w_h\cdot\x,\bm{\theta})$, and
the gradient $\nabla_{\w_i}f$  for unit $i$ in the first hidden layer is equal to the instance $\x$ times a scalar that depends on $\x$ only via $\w_i\cdot\x$, and updating of the later layer weights $\bm{\theta}$ depends on the input only via $\w_i\cdot\x$ (i.e., via the computation in the first layer).
Therefore it is easy to show by induction on $t$,
that rotating all instances results in the same rotation of
all $\w_{i,t}$. Thus the rotation has no effect on the dot
products $\x\cdot\w_i$ computed at the input layer.
In contrast, learning with $f((\u_1\odot\v_1)\cdot\x,
\ldots, (\u_h\odot\v_h)\cdot\x, \bm{\theta})$,
where the parameters $\w_i$ connected to the input $\x$ are replaced by $\u_i\odot\v_i$
(“spindlified”), is not rotation invariant.

§.§ Lower bounds for rotation invariant algorithms

Our method for proving lower bounds builds on the following observation: Given any rotation invariant algorithm and any learning problem, the algorithm will achieve the same loss on all rotated versions of that problem. We can therefore consider a Bayesian setting where the problem is sampled uniformly from all rotated versions, and the optimal solution provides a lower bound on the loss of the algorithm. Intuitively, being rotation invariant forces an algorithm to be agnostic over all possible rotations of the problem, and hedging its bets in this way prevents it from excelling at any specific problem instance. In Section <ref> we apply this reasoning to linear regression to show that a rotation invariant algorithm cannot efficiently learn sparse solutions, because it must be equally efficient at finding any other solution (including rotated, non-sparse ones).

Formally, let the learning problem be defined by
(a) an input distribution $p_{\textrm{in}}(\Xfull)$ with the input matrix $\underset{n+1,d}{\Xfull} = [\X, \Xte^\top]$ consisting of the training matrix $\X$ and the test example $\Xte$
(b) an observation model $q(\yfull | \Xfull)$ which gives the joint conditional distribution over $n$ training outcomes and a test outcome, $\underset{n+1}{\yfull} = [\y, \yte]$,
a loss function $\mathcal{L}(\yh, y)$.
We assume the input distribution is rotationally symmetric, meaning $p_{\textrm{in}}(\Xfull) = p_{\textrm{in}}(\Xfull\U^\top)$ for any orthogonal $\U$.
The task of any algorithm will be to produce predictions $\yh(\Xte | \X, \y)$ to minimize the loss on the test outcomes, $\mathcal{L}(\yh,\yte)$.
Note that this setup allows arbitrary conditional dependencies among observations (not just iid problems), including dependencies between the training and the test sets.

For any orthogonal $\underset{d,d}{\U}$, define the rotated observation model as
$\;q_{\U}(\yfull|\Xfull) = q(\yfull|\Xfull\U^{\top})$.
Now define a new learning problem by first sampling $\U$
uniformly (under the Haar measure $p_{\textrm H}$) and then generating observations according to $q_{\U}$. This is equivalent to a symmetrized observation model  $\mathring{q}$ that is a mixture over all $q_{\U}$:

\begin{align*}
\mathring{q}(\yfull | \Xfull) = \int q_{\U}(\yfull | \Xfull) {\textrm d}p_{\textrm H}(\U)
\end{align*}

The Bayes optimal prediction can be expressed by computing a posterior over $\U$ and integrating expected loss over this posterior:

\begin{align*}
\yh^{\star}(\Xte | \X,\y) = \arg\min_{\yh} \int
\mathbb{E}_{\yte\sim
q_{\U}(\cdot\vert\Xfull, \y)}[\mathcal{L}(\yh,\yte)]
%| \Xfull,\y]
{\textrm d}\U.
\end{align*}

Thus $\mathring{q}$ is difficult, especially for large $d$, because equal prior probability must be given to all possible rotations. We define the optimal expected loss on this problem as

\begin{align*}
L_{\mathcal{B}}(\mathring{q}) = \mathbb{E}_{\Xfull\sim
p_{\textrm{in}}, \yfull\sim\mathring{q}(\cdot|\Xfull)}[\mathcal{L}(\yh^{\star}(\Xte | \X, \y),\yte)].
\end{align*}

Our first result is that the performance of any rotation
invariant algorithm on the original problem, defined by
$q$, is lower bounded by $L_{\mathcal{B}}(\mathring{q})$.

Given a rotationally symmetric $p_{\textrm{in}}(\Xfull)$, an observation model $q(\yfull|\Xfull)$, a loss function $\mathcal{L}$, and a rotationally invariant learning algorithm $\yh_n(\cdot | \X,\y)$, define the expected loss

\begin{align*}
L_{\yh}(q) = \mathbb{E}_{\Xfull \sim p_{\textrm{in}}, \yfull \sim q(\cdot|\Xfull), Z}[\mathcal{L}(\yh(\Xte | \X,\y), \yte)],
\end{align*}

This loss is bounded by $L_{\yh}(q) \ge L_{\mathcal{B}}(\mathring{q})$.

The proof in Appendix <ref> shows any rotationally invariant algorithm will satisfy $L_{\yh}(q_{\U}) = L_{\yh}(q)$ for all $\U$. This implies $L_{\yh}(q) = L_{\yh}(\mathring{q}) \le L_{\mathcal{B}}(\mathring{q})$.

As we show in this paper, a consequence of Theorem <ref> is that rotational invariance prevents efficient learning of problems characterized by properties that are not rotationally invariant, such as sparsity. Algorithms with inductive biases for such properties are necessarily not rotation invariant and can give dramatically better performance.
Although we have stated the theorem in terms of rotational
invariance, it is easily extended to other transformation
groups $\mathcal{T}$ on the input (by replacing $\U$ with
elements of $\mathcal{T}$ and requiring $p_{\textrm{in}}$ to be symmetric under $\mathcal{T}$).
For example, natural gradient descent (NGD) is often touted for being invariant to arbitrary smooth reparameterization [Amari and Douglas, 1998], but this invariance comes at a cost of being unable to efficiently learn in environments that are not invariant in this way.
In particular, since spindly and fully connected networks are related by smooth reparameterization,
NGD performs equivalently on both (see discussion in [Kerekes et al., 2021]).
Thus the lower bounds we prove apply to NGD
on the network in Figure <ref>,
whereas we show vanilla gradient descent on this network breaks the lower bound.

§.§ Lower bound for least-squares regression

We demonstrate how Theorem
<ref> can provide a quantitative lower
bound for a specific class of learning problems. We then show how this bound is easily beaten by non-rotation-invariant algorithms in Section <ref>.
In the specific problem class we consider, the number of training examples is $n=md$ for some integer $m$,
and $\X$ consists of $m$ stacked copies of a matrix
$\H=\sqrt{d}\underset{d,d}{\V}$ (i.e. $\X = [\underset{\times m}{\H; \ldots ;\H}]$),
where $\V$ is a random orthogonal matrix distributed according to the Haar measure.
[We multiply $\V$ by $\sqrt{d}$ in order to keep the lengths of examples $\|\x_i\|$ (rows of $\X$) equal to $\sqrt{d}$, so that their coordinates $x_{ij}$ (individual features) are of order $1$ on average; any other scaling would work as well, resulting only in a relative change of the effective label noise level.]
Thus $p_{\rm in}(\X)$ is rotationally symmetric. The test input is one of the rows of $\H$,  $= _k$, with index $k$ drawn uniformly at random.
We assume the labels are the first feature of $$ plus Gaussian noise, that is
\begin{equation}
q(\yfull|\Xfull) = \mathcal{N}(\yfull | \Xfull \e_1, \sigma^2 \I_{n+1}), \qquad
\text{ where }\e_1=(1,0,\ldots,0)^\top,
\label{eq:model}
\end{equation}
which can be be equivalently written as
$$\y = \X \e_1 + \bm{\xi}, \text{ where } \bm{\xi} \sim N(\bm{0}, \sigma^2 \I_{md}),
\;\text{ and }\;\,
\yte = \h_k^\top \e_1 + \xi_{\rm{te}}, \text{ where } \xi_{\rm{te}} \sim N(0, \sigma^2).
Note that while the inputs are shared in the training and the test parts, the test label
is generated using a `fresh' copy of the noise variable $ξ_te$.
The fixed choice of $_1$ as opposed to any other $_i$ is made w.l.o.g. Indeed, Theorem \ref{t:gen_lower_bound} implies that a rotationally invariant algorithm will have the same loss for $_1$ as it will for $$ for any other unit vector $$.

The accuracy of prediction $= (·| ,)$ on the test set $(_k, )$ %an example $(\x,y)$ is measured by
is measured by the \emph{squared loss} $ℒ(,) = (- )^2$.
Let $$ denote the vector of predictions for all possible test instances, $_k = (_k | , )$, and let $_te$ denote the vector of test labels, $y_te_k = _k^⊤_1 + ξ_te$, which can be
jointly written as $= _1 + ξ_te1$ with $1 = (1,…,1)$.
%(y - \yh(\x))^2$.
%The prediction function is evaluated on a test set
% $(\H, \y_{te})$, where $\y_{te}=\H\e_1+\bm{\xi}_{te}$
% with a \emph{fresh set of noise variables} $\bm{\xi}_{te} \sim N(\bm{0}, \sigma^2 \I_d)$.\footnote{\Red{Here, we deviate from the setup of Section \ref{s:gen_lower_bound}, where the
%test set consisted of just a single example. One can, however, easily reduce the current problem to the previous setup by assuming the test example is just one row of $\H$ chosen uniformly at random: the expected loss will be equal to the average loss over the test set.}}
%Let $\predy = \yh(\H) \in \mathbb{R}^d$ be a vector of predictions obtained by applying the prediction function to all instances (rows) in $\H$.
The expected value of the loss over the random choice of $k ∈{1,…,d}$ and over the independent test label noise is given by:
\begin{align*}
\EE_{k,\xi_{\rm{te}}}[\mathcal{L}(\yh,\yte)] &= \frac{1}{d} \EE_{\xi_{\rm{te}}} \left[ \|\predy - \y_{\rm{te}}\|^2 \right]
= \frac{1}{d} \EE_{\xi_{\rm{te}}} \left[ \|\predy - \H \e_1 + \xi_{\rm{te}}\bm{1} \|^2 \right] \\
&= \frac{1}{d} \|\predy - \H \e_1\|^2
+ \frac{2}{d} \underbrace{\EE_{\xi_{\rm{te}}} \left[\xi_{\rm{te}}\right]}_{=0} (\predy - \H \e_1)^\top  \bm{1}
+ \frac{1}{d} \EE_{\bm{\xi}_{\rm{te}}} \underbrace{\left[ \xi_{\rm{te}}^2 \right]}_{=\sigma^2} \|\bm{1}\|^2 \\
&= \frac{1}{d} \|\predy - \H \e_1\|^2 + \sigma^2.
\end{align*}
Clearly, the expression above is minimized by setting the prediction vector to $^⋆ = _1$, and thus the
smallest achievable expected loss is equal to $σ^2$.
Subtracting this loss, we get the expression for the excess risk of the learning algorithm, which we call the \emph{error} of $$:
\begin{align*}
e(\predy) = \EE_{k,\xi_{\rm{te}}}[\mathcal{L}(\yh,\yte)] -
\EE_{k,\xi_{\rm{te}}}[\mathcal{L}(\yh^{\star},\yte)] =
\frac{1}{d} \|\predy - \H \e_1\|^2.
\end{align*}
When the prediction is \emph{linear}, $= $ for
some weight vector $∈ℝ^d$, we can also refer to the error of $$
as the error of its predictions:
\begin{equation}
e(\predw) = \frac{1}{d} \|\H \predw - \H \e_1\|^2
= \frac{1}{d} (\predw - \e_1)^\top \underbrace{\H^\top \H}_{d \I} (\predw - \e_1)
= \|\predw - \e_1\|^2.
\label{eq:error_w}
\end{equation}
%\ell(\predw) = \frac{1}{d} \|\H \predw - \H \e_1\|^2 + \sigma^2
%= \frac{1}{d} (\predw - \e_1)^\top \underbrace{\H^\top \H}_{d \I} (\predw - \e_1) + \sigma^2  = \|\predw - \e_1\|^2 + \sigma^2.
%and the corresponding error:
%e(\predw) = \ell(\predw) - \ell(\e_1) =  \|\predw - \e_1\|^2.
%\ell(\w) &= \frac{1}{d} \EE_{\bm{\xi}_{te}} \left[ \|\H \w - \y_{te}\|^2 \right]
%= \frac{1}{d} \EE_{\bm{\xi}_{te}} \left[ \|\H \w - \H \e_1 + \bm{\xi}_{te}\|^2 \right] \\
%&= \frac{1}{d} \|\H \w - \H \e_1\|^2
%+ \frac{1}{d} \underbrace{\EE_{\bm{\xi}_{te}} \left[2 (\w - \e_1) \H^\top \bm{\xi}_{te}\right]}_{=0}
%+ \frac{1}{d} \EE_{\bm{\xi}_{te}} \underbrace{\left[ \| \bm{\xi}_{te}\|^2 \right]}_{d \sigma^2} \\
%&= \frac{1}{d} (\w - \e_1)^\top \underbrace{\H^\top \H}_{d \I} (\w - \e_1) + \sigma^2  = \|\w - \e_1\|^2 + \sigma^2.

We prove the following lower bound for rotationally invariant algorithms on this problem:

\begin{theorem}
Let $\underset{d,d}{\V}$ be a random orthogonal matrix, and let $\H = \sqrt{d} \V$. Let $(\X,\y)$ be the training
set with $\underset{md,d}{\X}=[\H; \ldots ;\H]$ and labels $\y$ generated according to \eqref{eq:model}.
Then the expected error (with respect to $\V$) of any rotation-invariant learning algorithm is at least
\[ \EE_{\V}[e(\predy)] \ge \frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}.\]
\label{thm:lower_bound}
\end{theorem}
\vspace{-7mm}

The proof in Appendix~\ref{a:thm2-proof} uses Theorem~\ref{t:gen_lower_bound} and closely follows Section \ref{s:gen_lower_bound}:
We start with a Bayesian setting with the rotated observation model
$q(|^⊤)$ for a random orthogonal matrix $$.
This is equivalent to simply rotating the
target weight vector by $^⊤$, because $^⊤_1 = (^⊤_1)$.  Thus we can equivalently consider a linear model $= + ξ̃$, where $$ is drawn
uniformly from the unit sphere $𝒮^d-1 = {∈ℝ^d =1 }$.
Given the square loss function and linear observation model, the optimal Bayes predictor is based on the posterior mean, $𝔼[| , ]$.

Even though the posterior mean does
not have a simple analytic form for the prior distribution over a unit sphere, we use results
from [Marchand, 1993, Dicker, 2016] to show that the Ridge Regression (RR) predictor (which is
the Bayes predictor for the Gaussian prior) with appropriately chosen regularization constant has expected error at most $1/d σ^2/σ^2 + m$ larger
than that of the Bayes predictor. Thus, it suffices to analyze the RR predictor
which we show to be at least
$σ^2/σ^2 + m$. This implies the
Bayes error is at least $d-1/d σ^2/σ^2 + m$, and no other algorithm can achieve any better error for this problem.
Finally, due to the rotation-symmetric distribution of the inputs, we can now apply Theorem
\ref{t:gen_lower_bound} which implies that every rotation
invariant algorithm has error at least
$d-1/d σ^2/σ^2 + m$ for the original sparse linear regression
problem $= _1 + ξ̃$.
For the sake of completeness, we also give an i.i.d.
version of this lower bound in Appendix \ref{app:lower_bound_iid}.

Previous lower bound proofs for sparse problems for the
case when the number of examples is less than the input
dimension used fixed choices for input matrix $$ such as the $d$ dimensional Hadamard matrix [Kivinen et~al., 1997, Warmuth and Vishwanathan, 2005].
In the over-constrained case, the lower bound does not hold for any
fixed choice of $$ (of full rank).
For any fixed full-rank $$, there exists a row vector $$̌
s.t. $\v \X = \e_1^\top$, and the linear algorithm $\yh(\x | \X, \y) = \v\X\, \x$
achieves minimal loss $\sigma^2$ while being trivially rotationally invariant
because $\v \X \U^\top \,\U\x=\v\X\,\x$.
Thus the assumption of rotationally symmetric input distribution is essential.
[In the Hadamard matrix based lower bounds of [Warmuth and Vishwanathan, 2005],
the ease of learning say $\e_1$ is overcome by averaging over all $n$ targets $\e_i$.]

This counterexample requires knowledge of the order of the
examples (i.e., the algorithm makes different predictions
if the rows of $\X$ are permuted). Now consider $\X =
[\H;\dots;\H]\, {\rm Diag}(\p)$ with $p_1=2,p_{i>1}=1$.
Then the target $\e_1$ is embedded as the first principal
component $\mathbb{PC}_1(\X^\top\X)$, so an algorithm that
used this as its weight vector, $\yh(\x|\X,\y) =
\mathbb{PC}_1(\X^\top\X)^\top\x$, would achieve minimal
loss. This algorithm is rotationally invariant but fails if
the input matrix $\X$ was rotated.

§ UPPER BOUNDS

We now show how to break the above lower bound of
$\frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}$ on the error
(excess risk) of rotation invariant algorithms.
We focus on the upper bounds for EG$^\pm$ and Approximated EGU$^\pm$.
(Similar upper bounds for the spindly network, as well as
the novel priming method [Warmuth and Amid, 2023] are given in
appendices <ref> and <ref>,
We use the same setup as
in the lower bound, i.e. $\X$ consists of $m$ stacked copies of
a matrix $\H = \sqrt{d}\V$, when $\V$ is a rotation matrix
and $\y$ is sparse linear (the first component) plus Gaussian noise with variance $\sigma^2$.
The only difference is that for the lower bound, $\V$ was
randomly chosen, but the upper bounds hold for any rotation
matrix $\V.$ The upper bound on the error that we
[Note that at this point, we do not consider general input
matrices $\X$ in the upper bounds. Allowing arbitrary covariance
structure makes the analysis much more complicated and
the general case is left for future research.]
is essentially $O(\sigma^2 \frac{\log d}{md})$.

We begin with a bound for the normalized version of the
multiplicative update called EG$^\pm$.
The proof relies on the facts
(a) under large learning rate, the EG$^\pm$ weight estimate becomes argmax over the negative gradient coordinates, and
(b) with high probability, the least component of the gradient is the first one, so that the argmax is $\e_1$.
EG$^\pm$ makes use of the fact that the norm of the linear
target $\e_1$ is 1 and this additional knowledge allows the speedup.
[We can also provide an upper bound
on the error of the online version of EG$^{\pm}$ for an arbitrary input matrix $\X$
with fixed feature range via a standard
worst-case regret analysis followed by the online-to-batch
conversion (see e.g. [Kivinen and Warmuth, 1997]).
The bound so obtained would however give a slower rate of
order $O(\sqrt{\log d / (dm)}$, which still has a substantially better dependence on dimension
$d$ than the lower bound from the previous
section.] We then prove our main upper bound for the unnormalized EGU$^\pm$.

§.§ Upper bound for the Exponentiated Gradient update

We consider a batch version of the 2-sided
Exponentiated Update algorithm (EG$^{\pm}$) [Kivinen and Warmuth, 1997].
The batch EG$^{\pm}$ algorithm maintains two vectors, $\v^+_t$ and $\v^-_t$, and its weight estimate is given by $\w_t = \v^+_t - \v^-_t$.
It starts with a set of weights $\v_1^+ = \v_1^- = \frac{1}{2d} \1$ (so that $\|\v_1^+\|_1 + \|\v_1^-\|_1 = 1$)  and updates according to

\[
\v^+_{t+1} \propto \v^+_t \odot e^{-\eta \nabla \emploss(\w_t)}, \qquad \v^-_{t+1} \propto \v^-_t \odot e^{\eta \nabla \emploss(\w_t)},
\]

where $\odot$ is component-wise multiplication, and the normalization enforces
$\|\v^+_{t+1}\|_1 + \|\v^-_{t+1}\|_1 = 1$.

The expected error of the batch EG$^{\pm}$ algorithm after the first iteration is bounded

\[
e(\w_2) \le 2 d e^{-\eta} + 8 d e^{-\frac{md}{32\sigma^2}}.
\]

The gradient of the loss can be written as $2 (\w - \e_1 -
\bm{\zeta})$, where  $\bm{\zeta}$ are i.i.d. Gaussian noise variables which
are combinations of the original noise variables. Using the
deviation bound for Gaussians together with union bound, with probability at
least $1-2d \exp\{-md / (32 \sigma^2)\}$ all noise variables are bounded by $1/4$.
Thus, the first coordinate of the negative gradient is larger than all other
coordinates by at least $1$, so that the first weight exponentially dominates all
the other weights already after one step of the algorithm, and the error drops
to $2 d e^{-\eta}$. If the high probability event does not hold, we bound
the error by its maximum value $4$. See Appendix <ref> for full proof.

§.§ Upper bound for the Approximated Unnormalized Exponentiated Gradient update

We now drop the normalization of EG$^\pm$ and use the approximation $\exp(x)\approx 1+x$.
[This approximated version of EGU was introduced in
[Kivinen and Warmuth, 1997]. It was also used in the normalized update PROD
[Cesa-Bianchi et al., 2007].]
The resulting approximation of the unnormalized EGU$^\pm$ algorithm updates its weights as

\begin{equation}
\w_{t+1} = \w_t - \eta\sqrt{\w_t^2 + 4\beta^2}\; \nabla \emploss(\w_t),
\label{eq:spindly_update_main}
\end{equation}

where $\beta > 0$ is a parameter and all operations
(squaring the weights and taking square root) are done
component-wise (see Appendix <ref> for a derivation). When $\beta=0$ this update is closely related to gradient
descent on the spindly network of Figure <ref>.
We start with $\w_1=\bm{0}$. Note that
unlike EG$^{\pm}$ the update does not constrain
its weights by normalizing. Nevertheless, we show that the algorithm achieves an upper bound
on the error that is essentially $O(\frac{d}{\log d})$ better than the error of
any rotation invariant algorithm:

Assume $d \ge 4$ is such that $\sqrt{d}$ is an integer. Consider the Approximated EGU${}^{\pm}$ algorithm (<ref>) with $\beta=1/(2d)$ and learning rate $\eta=1/4$ and let $m \ge 8 \sigma^2 \ln \frac{2 d}{\delta} = \Omega(\sigma^2 \log (d/\delta))$ . With probability at least $1-\delta$, the algorithm run for $T=4\sqrt{d}$ steps achieves error bounded by:

\[
e(\w_{T+1}) \le \frac{10 \sigma^2 \ln \frac{2d}{\delta}}{md} + 9 e^{-\frac{8}{3} \sqrt{d} + 2 \ln d}
= O(\frac{\sigma^2\log d}{md} + e^{-\frac83\sqrt{d}})
\]

sketch: (full proof in Appendix <ref>)
Similarly as for EG$^{\pm}$, with high probability all noise variables are small. We can interpret the weight update (<ref>) as the gradient descent update on the weights with the
effective learning rates $\eta \sqrt{\w_t^2 + 4 \beta^2}$. Since $\beta=1/(2d)$
and $\w_1 = \bm{0}$, these learning rates are initially small and of order $O(1/d)$. We then show by a careful analysis of the update that the weights $w_{t,i}$ remain small for all coordinates $i$ except $i=1$, and so do the associated learning rates. Therefore, after $T$ steps, the algorithm does not move significantly away from zero on these coordinates. Meanwhile, the weight on the first coordinate $w_{t,1}$ increases towards $1$. While the initial rate of increase is small as well, it accelerates over time as $\eta \sqrt{w_{t,1}^2 + 4 \beta^2}$ increases due to increasing $w_{t,1}$. Eventually $w_{T+1,1}$ gets very close to $1$ after $T$ steps of the algorithm.

Discuss the generality of these bounds: they hold against worst
case data but you need to average over all online weight vectors.

Additional square root terms.

Stronger bounds for sparse linear plus iid noise with mirror descent updates:

First always discuss what the update does if there is no noise
and discuss that the noise is unlike to confuse the $\argmax$
calculation. Proofs in appendices.

Mismatch that we should still try to avoid:
The current proofs are for Hadamard. But the lower bounds are for rows of rotation matrices.
Theorem <ref> holds for $\X=\sqrt{d}[\V;\dots;\V]$ for any $\V$, not just Hadamard: The loss gradient is $2 (\w - \e_1 - \bar{\bm{\xi}})$ where $\bar{\bm{\xi}} = \X^{\top}\bm{\xi}/d \sim \mathcal{N}(0,\frac{\sigma^2 m}{d} \I_d)$. The rest of the proof goes through unmodified.

* EG$^\pm$ exponential convergence bound for single batch update.
In the noise-free case its just an $\argmax$ of the batch
gradients of the initial weight vector.

Caviat? The algorithm must know the L$_1$ norm of the target.

* EGU$^\pm$ exponential convergence bound for iterative batch update.
Discuss that in the noise-free case its already more
tricky, because algorithm does not know norm of target?

Actually we prove this for the spindly reparameterization.

Skip priming:
One advantage of the update family that predicts with
linear combinations of examples is that there are closed
expressions for the weight vector that minimizes the square
loss (I.e. the LLS solution). In <cit.> the question was
asked whether there are closed form expressions for weight
vectors that favor sparsity. Experimentally this was
achieved by first finding the LLS solution $\w^*$. Then the
$i$th feature is “primed” by multiplying it with $w^*_i$.
Finally ridge regression is applied on the transformed examples.
Experimentally this favors sparsity when the
tradeoff parameter for ridge regression optimized.
Here we give a proof for a special case...

§ GRADIENT FLOW TRAJECTORIES

In our setting of overconstrained optimization without explicit regularization, all algorithms operate on the same loss landscape but follow different paths to the minimum.
To understand these differences, we derive exact expressions for the weight trajectories for the continuous time (gradient flow) versions of several representative algorithms.
We then use these results to analyze the algorithms' ability to learn sparse targets without overfitting to noise.

§.§ Preconditioning, mirror descent, reparameterization, and Riemannian descent

Gradient descent (GD) on a given loss function can be generalized in several ways, including preconditioned GD, mirror descent (MD), GD on reparameterized variables, and Riemannian GD. We show how they are all interchangeable in the continuous-time case.

Preconditioned GD premultiplies the gradient with a matrix $\bm{P}$ that can depend on the current weights or the training history. Explicit preconditioning methods include adaptive learning rate algorithms such as Adagrad [Duchi et al., 2011], natural gradient descent which preconditions with the inverse Fisher information matrix [Amari, 1998], and variational Bayesian methods that precondition with the prior or posterior covariance [Lambert et al., 2022, Chang et al., 2023].

\begin{equation}
\dot{\w} = -\bm{P} \nabla_w L
\end{equation}

Reparameterizing refers to defining some injective function $\widehat{\w} = g(\w)$ and performing GD on the new variable:

\begin{equation}
\dot{\widehat{\w}} = -\nabla_{\widehat{\w}} L
\end{equation}

MD [Nemirovsky and Yudin, 1983] involves an invertible mirror map $\widetilde{\w} = f(\w)$ where $f$ is the gradient of a convex function, with the update

\begin{equation}
\dot{\widetilde{\w}} = -\nabla_{\w} L
\end{equation}

Riemannian GD generalizes to parameter spaces with non-Euclidean geometry (specifically, differentiable manifolds) [Bonnabel, 2013]. In Riemannian geometry the update $\dot{\w}$ lies in the tangent space $T_{\w}$ (a direction of motion based at $\w$) whereas the gradient $\nabla_{\w} L$ lies in the cotangent space $T_{\w}^*$ (the dual of $T_{\w}$).
Unlike in Euclidean geometry, $T_{\w}$ and $T_{\w}^*$ cannot be identified and instead mapping between them requires a metric $\bm{\Gamma}_{\w}:T_{\w}\to T_{\w}^{*}$ or equivalently $\bm{\Gamma}_{\w}:T_{\w}\times T_{\w}\to\mathbb{R}$.
This leads to an update that is the direction of steepest descent with respect to the geometry:

\begin{equation}
\dot{\w} = -\bm{\Gamma}_{\w}^{-1} \nabla_{\w} L
\end{equation}

The continuous-time versions of reparameterized GD, MD, and Riemannian GD are all equivalent to preconditioned GD under the relation

\begin{equation}
\bm{P}
= \frac{\partial\w}{\partial\widehat{\w}} \left(\frac{\partial\w}{\partial\widehat{\w}}\right)^\top
= \frac{\partial\w}{\partial\widetilde{\w}}
= \bm{\Gamma}_{\w}^{-1}
\end{equation}

Appendix <ref> gives the proof and then uses this result (a) to derive a close connection between EGU and the spindly network and (b) to show how EGU implicitly operates on a non-Euclidean geometry
(Figure <ref>) that tends to hold it in sparse regions of parameter space.
Rotationally invariant algorithms cannot do this because their implicit geometry must be rotationally symmetric.

§.§ Trajectory solutions

We derive analytic trajectories for the simple linear regression problem from Sections <ref> and <ref>, for the continuous-time versions of EGU, EGU$\pm$, primed GD, and Adagrad [Duchi et al., 2011]. For space reasons we present the expressions for the trajectories here and provide the derivations and full details in Appendix <ref>.

For continuous EGU the trajectory is (with $c_i$ a constant depending on initial conditions):

\begin{align}
\label{e:EGU-trajectory}
w_i(t) &= \tfrac{1}{2}w^{\rm LS}_i \left(1 + \tanh \left( w^{\rm LS}_i t + c_i \right) \right)
\end{align}

For continuous EGU$\pm$ the trajectory is (with $\tau_i$ a transform of $t$):

\begin{align}
\label{e:EGUpm-trajectory}
w_i(t) &=
\frac{w^{\rm LS}_i \sinh\tau_i + 1}
{\sinh\tau_i - w^{\rm LS}_i}
\end{align}

For primed gradient flow, which is similar to the priming method analyzed elsewhere in this paper except the second stage uses gradient flow instead of ridge regression, the trajectory is

\begin{equation}
\label{e:primeGD-trajectory}
w_i(t) = w^{\rm LS}_i + \left(w_i(0) - w^{\rm LS}_i \right) e^{-2t\left(w^{\rm LS}_i\right)^2}
\end{equation}

For continuous-time Adagrad the trajectory is (with constants $k_i,\ell_i$ depending on initial conditions):

\begin{align}
\label{e:Adagrad-trajectory}
w_i(t) &= w^{\rm LS}_i - {\rm sign}(w^{\rm LS}_i - w_i(0)) \sqrt{-\frac{8}{\beta k_i} W\left(-e^{-k_i(t+\ell_i)}\right)}
\end{align}

with Lambert's W function defined by $x = W(x)e^{W(x)}$.

§.§ Experimental visualization

The analytic trajectories derived above are visualized in
Figure <ref>a for a 2d version of our
noisy sparse regression problem. This toy problem is sufficient
to bring out key differences among the algorithms.
Because all considered algorithms optimize the same loss
they all converge to the same least squares solution for
each training set (grey dots), but they follow qualitatively different trajectories.
Priming and EGU$\pm$ pass close to the sparse target (black
dot) while GD goes straight to the LLS solution in each
run. Surprisingly Adagrad veers away from the
sparse target.
[In Theorem <ref> we also prove exponential
convergence of the error of EG$^\pm$. We deemphasize this
algorithm because it normalizes based on the norm of the
true weight vector and this may be unreasonable.
Indeed using learning rate $\eta=200$, EG$^\pm$ achieves error below $10^{-300}$
(not shown). Larger learning rates cause numerical

These differences in trajectory have dramatic implications for excess test loss as shown in Figure <ref>b.

Priming and EGU$\pm$ have deep dips in the curves
corresponding to the fact that their trajectories pass
close by the target while the dips of GD and ridge
regression are too small to be visible.
All algorithms require early stopping to
“catch” bottom dip of their error curves.

Note that the rotation invariant algorithms both reach but do not beat the lower
bound of Theorem <ref> while the non-rotation invariant algorithms greatly outperform the upper
bound of Theorem <ref>, showing that the bound
we were able to prove is conservative.

Finally, the bounds of Sections <ref> and <ref>
and the curves of Figure <ref> assume the instance matrix
$\X$ comprises $m$ copies of a
scaled $d$ dimensional rotation matrix. This unusual setup was needed for technical reasons.
However, if we used a Gaussian instance matrix
$\underset{md,d}{\X}$ together with sparse
linear targets plus Gaussian noise, then the curves
would remain essentially unchanged (not shown).

[Loss curves]

(a): Trajectories for key algorithms on 10 sampled datasets in toy 2d regression. Black dot shows true sparse target. Grey dots show least squares solution for each dataset. All algorithms converge to LS solution but differ in the order of learning the signal versus noise. EGU$\pm$ and priming learn the sparse target efficiently, overfitting to the noise only at the very end of training. Because GD is rotation invariant it cannot produce this behavior and instead learns the signal and noise at equal rates. Intriguingly, Adagrad shows the opposite behavior and is especially inefficient at learning a sparse target, because of its per-parameter adaptive learning rates.
(b): Excess test loss averaged over 100 runs of linear regression with $d=1024$. Dotted lines show the lower bound for rotationally invariant algorithms from Theorem <ref> (applies to GD and Ridge) and the upper bound for Approximated EGU$^{\pm}$ from Theorem <ref> (using $\delta=.001$).

§ NOISE EXPERIMENTS ON FASHION MNIST

The assumption of rotationally symmetric input distribution is too strong for real datasets. Therefore it is important to complement our theoretical results with experimental demonstration that standard rotationally invariant algorithms cannot exploit asymmetries in real data. Here we report experiments on the Fashion MNIST [Xiao et al., 2017] dataset.
One could trivially cheat this task by hard-coding a lookup table of the test set,
and this can even be done with a rotation invariant algorithm:
If $\X^\top\X$ is the (known, asymmetric) covariance of the original training set,
then when presented with rotated training data
$\Z=\X\U^\top$, the algorithm infers the unknown $\U$ to
high precision by solving $\Z^\top\Z=\U\X^\top\X\U^\top$.
It then un-rotates $\bm{z}_{\rm te}$ and uses the lookup table.

Such a bespoke algorithm is of course not of interest; the question is whether standard
rotation invariant algorithms based on gradient descent can do anything similar, or whether they are limited by their rotation invariance just as they are in the idealized symmetric-input setting.

We use a
multilayer feedforward network with two hidden layers of
size 256 each. We consider two cases for the input layer
weights: 1) fully-connected (each 1st layer hidden node is
connected to all inputs) and 2) “spindly” (each 1st layer hidden node
is connected to all inputs via the network of Figure
In the noise-free case, both variants of the network
achieve the same top one accuracy of $85\%$ (although
the spindly version takes longer to train).
Next, we double the number of features of
the examples by augmenting each example with uniformly sampled noise pixels in the range $[-1, 1]$.
With noisy augmentation, the spindly network still achieves $85\%$ test accuracy
while the fully-connected network gets only to $71\%$
(after a much longer training phase).
In addition, the learned weights by
the two networks are significantly different: The
fully-connected network assigns almost the same magnitude
of weights to the noisy features as to the image features,
while the spindly network allocates much larger weights to the image features.
This illustrates that rotation invariant algorithms have a harder time ignoring
the noisy features.

Finally, to further compare the different sensitivity to the
informativeness of features, we augment each example on top of the noise
with its one-hot representation of the 10 target class labels.
This splits the features into three categories
in terms of their informativeness: 1) highly informative label features, 2) less informative image features, and 3) noise features with no (structured) information. The spindly network achieves $100\%$ test accuracy while the fully-connected network gets to $98\%$. We observe that the spindly network assigns much larger weights (in magnitude) to the label features while almost ignoring the rest. This phenomenon is less prominent in the fully-connected network. We defer all the details to Appendix <ref>.

§ OPEN PROBLEMS

We gave a lower bound for rotation invariant neural nets
for learning noisy sparse linear problems.
Our work suggests the following approach: The
structure of the network plus the weight update imply a
certain set of invariances and the invariances lead to
lower bounds for the Bayes optimal algorithm of the model
constructed from the invariances.
For example a certain two-sided invariance characterizes
the matrix version of
multiplicative updates <cit.>. Is it thus possible to prove lower
bounds for these algorithms for dense linear targets that are
beaten by vanilla gradient descent?

A more immediate question is whether the sparse classification
problems that are provably hard to solve for rotation invariant
neural net algorithms but are easy to solve for convolutional neural nets
can be lifted to the regression setting [Li et al., 2021].
Note that the existing lower bounds are for the less interesting case where the
number of examples is less than the VC dimension. The
question is whether with the addition of noise there are
still gaps towards the optimal algorithm when the number of
examples is larger than the VC dimension.
Transformers clearly have an ability to access individual
tokens. The question is which structural feature of
transformers enables them to do that and learn sparse
linear problems. Partial results along this line
recently appeared in <cit.>.

Another interesting question is when transforming the
instances (i.e. a kernelization) “helps” rotation invariant algorithms. Note that
linear neurons that are fed transformed instances when
updated with gradient descent again have the $1-\sfrac kd$
lower bound in the noise-free under constraint case [Warmuth and Vishwanathan, 2005].
Proving similar lower bounds for the noisy over constraint case
as done in this paper would be a first step.

Also, the upper bounds proven in this paper focus (with
the exception of priming) on mirror
descent and its approximations. However there is large body
of work based on L$_1$ regularization that is also biased towards
sparsity (See e.g. <cit.>). How these updates are related to mirror
descent is an interesting research topic.

Recall that multiplicative update are mirror descent using the
$\log$ link functions. By stretching the link in the EG$^\pm$
version it was shown that EG$^\pm$ has gradient descent as
a special case <cit.>. Actually we observe that this stretching method can be used to
enhance any mirror descent update with a hyper parameter
that can realize gradient descent as a special case.
However no practical algorithm have been developed based on this observation.
The “full versions” of many common optimization
algorithms [Abdulkadirov et al., 2023] such as
AdaGrad, Fisher, Adam, RMS Prop are all rotation invariant.
However in practice diagonalized versions of these updates
are used. In a major step forward we were able to solve the
differential equations for all these algorithms
when $\bm{\Sigma}=\lambda\I$ (the arbitrary $\bm{\Sigma}$ cases are
challenging open problems).
This lets us analyze the inductive biases of all these
algorithms on linear neurons based on their gradient flow trajectories.
For example, the diagonalized AdaGrad is biased away from sparsity.
So using sparsity as a yard stick, we show
that already for linear neurons, dramatic differences
can occur compared to the optimal algorithms.
The question is whether these differences remain
when running the same algorithms on neural nets with any number of hidden layers.

[Abbe and Boix-Adsera, 2022]
Emmanuel Abbe and Enric Boix-Adsera.
On the non-universality of deep learning: quantifying the cost of
In Advances in Neural Information Processing Systems,
volume 35, pages 17188–17201. Curran Associates, Inc., 2022.

[Abdulkadirov et al., 2023]
Ruslan Abdulkadirov, Pavel Lyakhov, and Nikolay Nagornov.
Survey of optimization algorithms in modern neural networks, 04 2023.

[Amari and Douglas, 1998]
S. Amari and S.C. Douglas.
Why natural gradient?
In Proceedings of the 1998 IEEE International Conference on
Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. No.98CH36181),
volume 2, pages 1213–1216 vol.2, 1998.

[Amari, 1998]
Shun-Ichi Amari.
Natural gradient works efficiently in learning.
Neural computation, 100 (2):0 251–276, 1998.

[Amid and Warmuth, 2020]
Ehsan Amid and Manfred K. Warmuth.
Reparameterizing mirror descent as gradient descent.
In Proceedings of Advances in Neural Information Processing
Systems, volume 33, pages 8430–8439, 2020.

[Berger, 1985]
James O. Berger.
Statistical decision theory and Bayesian analysis.
Springer, 1985.

[Bonnabel, 2013]
Silvere Bonnabel.
Stochastic gradient descent on riemannian manifolds.
IEEE Transactions on Automatic Control, 580
(9):0 2217–2229, 2013.

[Cesa-Bianchi et al., 2007]
Nicolo Cesa-Bianchi, Yishay Mansour, and Gilles Stoltz.
Improved second-order bounds for prediction with expert advice.
Machine Learning, 660 (2-3):0 321–352, 2007.

[Chang et al., 2023]
Peter Chang, Gerardo Duràn-Martín, Alexander Y Shestopaloff, Matt
Jones, and Kevin Murphy.
Low-rank extended kalman filtering for online learning of neural
networks from streaming data.
arXiv preprint arXiv:2305.19535, 2023.

[Dicker, 2016]
Lee H. Dicker.
Ridge regression and asymptotic minimax estimation over spheres of
growing dimension.
Bernoulli, 220 (1):0 1–37, 2016.

[Duchi et al., 2011]
John Duchi, Elad Hazan, and Yoram Singer.
Adaptive subgradient methods for online learning and stochastic
Journal of Machine Learning Research, 120
(Jul):0 2121–2159, 2011a.

[Duchi et al., 2011]
John Duchi, Elad Hazan, and Yoram Singer.
Adaptive subgradient methods for online learning and stochastic
Journal of machine learning research, 120 (7),

[Gerchinovitz et al., 2020]
Sébastien Gerchinovitz, Pierre Ménard, and Gilles Stoltz.
Fano’s Inequality for Random Variables.
Statistical Science, 350 (2):0 178 – 201,

[Gunasekar et al., 2017]
Suriya Gunasekar, Blake E. Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur,
and Nati Srebro.
Implicit regularization in matrix factorization.
In Proceedings of Advances in Neural Information Processing
Systems (NeurIPS), pages 6151–6159, 2017.

[Jain et al., 2007]
Prateek Jain, Brian J. Kulis, and Inderjit S. Dhillon.
Online linear regression using burg entropy.
Technical Report TR-07-08, The University of Texas at Austin, feb

[Kerekes et al., 2021]
Anna Kerekes, Anna Mészáros, and Ferenc Huszár.
Depth without the magic: Inductive bias of natural gradient descent.
ArXiv, abs/2111.11542, 2021.

[Kingma and Ba, 2015]
Diederik Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
In International Conference on Learning Representations
(ICLR), San Diega, CA, USA, 2015.

[Kivinen and Warmuth, 1997]
Jyrki Kivinen and Manfred K. Warmuth.
Exponentiated gradient versus gradient descent for linear predictors.
Information and Computation, 1320 (1):0 1–63,

[Kivinen et al., 1997]
Jyrki Kivinen, Manfred K. Warmuth, and Peter Auer.
The Perceptron algorithm versus Winnow: linear versus logarithmic
mistake bounds when few input variables are relevant.
Artificial Intelligence, 970 (1-2):0 325–343,

[Lambert et al., 2022]
Marc Lambert, Silvere Bonnabel, and Francis Bach.
The recursive variational gaussian approximation (r-vga).
Statistics and Computing, 320 (1):0 10, 2022.

[Li et al., 2021]
Zhiyuan Li, Yi Zhang, and Sanjeev Arora.
Why are convolutional nets more sample-efficient than fully-connected
In International Conference on Learning Representations, 2021.

[Marchand, 1993]
Eric Marchand.
Estimation of a multivariate mean with constraints on the norm.
The Canadian Journal of Statistics, 210 (4):0
359–366, 1993.

[Nemirovsky and Yudin, 1983]
Arkadi Nemirovsky and David Yudin.
Problem Complexity and Method Efficiency in Optimization.
John Wiley & Sons, 1983.

[Ng, 2004]
Andrew Y. Ng.
Feature selection, $l_1$ vs. $l_2$ regularization, and rotational
In Proceedings of the 21-st International Conference on Machine
Learning, Banff, Canada, July 4-8. ACM New York, NY, 2004.

[Rigollet and Hütter, 2023]
Philippe Rigollet and Jan-Christian Hütter.
High-dimensional statistics, 2023.

[Shamir, 2018]
Ohad Shamir.
Distribution-specific hardness of learning neural networks.
Journal of Machine Learning Research, 190
(32):0 1–29, 2018.

[Vaskevicius et al., 2019]
Tomas Vaskevicius, Varun Kanade, and Patrick Rebeschini.
Implicit regularization for optimal sparse recovery.
In Proceedings of Advances in Neural Information Processing
Systems (NeurIPS), pages 2968–2979, 2019.

[Warmuth and Amid, 2023]
Manfred K. Warmuth and Ehsan Amid.
Open problem: Learning sparse linear concepts by priming the
In Proceedings of Thirty Sixth Conference on Learning Theory,
volume 195 of Proceedings of Machine Learning Research, pages
5937–5942. PMLR, Jul 2023.

[Warmuth and Jagota, 1998]
Manfred K. Warmuth and Arun Jagota.
Continuous and discrete time nonlinear gradient descent: relative
loss bounds and convergence.
In R. Greiner E. Boros, editor, Electronic Proceedings of Fifth
International Symposium on Artificial Intelligence and Mathematics.
Electronic,http://rutcor.rutgers.edu/$\tilde{~}$amai, 1998.

[Warmuth and Vishwanathan, 2005]
Manfred K. Warmuth and S.V.N. Vishwanathan.
Leaving the span.
In Proceedings of the 18th Annual Conference on Learning Theory
(COLT), pages 366–381, 2005.

[Warmuth et al., 2021]
Manfred K. Warmuth, Wojciech Kotłowski, and Ehsan Amid.
A case where a spindly two-layer linear network decisively
outperforms any neural network with a fully connected input layer.
In 32th International Conference on Algorithmic Learning Theory
(ALT), volume 132, pages 1–32. PMLR, 2021.

[Xiao et al., 2017]
Han Xiao, Kashif Rasul, and Roland Vollgraf.
Fashion-mnist: a novel image dataset for benchmarking machine
learning algorithms, 2017.
URL <http://arxiv.org/abs/1708.07747>.
cite arxiv:1708.07747Comment: Dataset is freely available at
https://github.com/zalandoresearch/fashion-mnist Benchmark is available at

Consider a problem where $d$ observations are rows of a $d \times d$ Hadamard matrix $\H$ together with label vector $$ generated according to:
\[
\y = \H \e_1 + \bm{\xi}, \qquad \bm{\xi} \sim N(\0, \I \sigma^2).
\]
Now, the learner is given $m$ copies of this dataset, each with its own set of labels based on a fresh version of noise variables, that is
$_t = _1 + ξ_t$, $t=1,…,m$. The learner is then evaluated on the $(m+1)$-th unseen copy of the dataset using the mean squared error.
Everything which is presented here can be generalized by replacing $$̋ with
$\sqrt{d} \U_t$, where $\U_t$ is an arbitrary orthogonal matrix, possibly different in each copy.

Given a fresh copy of labels in the test sample $\y$, the averaged (over $d$ test examples) expected (over the test label noise) squared loss of weight vector $\w$ is given by:

\begin{align*}
\ell(\w) &= \frac{1}{d} \EE_{\y} \left[ \|\H \w - \y\|^2 \right]
= \frac{1}{d} \EE_{\y} \left[ \|\H \w - \H \e_1 + \bm{\xi}\|^2 \right] \\
&= \frac{1}{d} \EE_{\y} \left[ \|\H \w - \H \e_1\|^2\right]
+ \frac{1}{d} \underbrace{\EE_{\y} \left[2 (\w - \e_1) \H^\top \bm{\xi}\right]}_{=0}
+ \frac{1}{d} \EE_{\y} \underbrace{\left[ \| \bm{\xi}\|^2 \right]}_{d \sigma^2} \\
&= \frac{1}{d} (\w - \e_1)^\top \underbrace{\H^\top \H}_{d \I} (\w - \e_1) + \sigma^2  = \|\w - \e_1\|^2 + \sigma^2.
\end{align*}

Plugging $\w = \e_1$ and subtract, we get the expression for the error of the algorithm, that is the excess risk:

\[
e(\w) = \ell(\w) - \ell(\e_1) =  \|\w - \e_1\|^2
\]

§ PROOF OF THE LOWER BOUND THEOREM <REF>

For any orthogonal $\U$, the algorithm's expected loss is

\begin{align*}
L_{\yh}(q_{\U}) &= \mathbb{E}_{\Xfull\sim p_{\textrm{in}}, \yfull\sim q_{\U}(\cdot|\Xfull),Z}[\mathcal{L}(\yh(\Xte | \X,\y), \yte)] \\
&= \mathbb{E}_{\Xfull\sim p_{\textrm{in}}, \yfull\sim q(\cdot|\Xfull \U^\top),Z}[\mathcal{L}(\yh(\Xte | \X,\y), \yte)] \\
&= \mathbb{E}_{\Xfull'\sim p_{\textrm{in}}, \yfull\sim q(\cdot|\Xfull'),Z}[\mathcal{L}(\yh(\Xte' | \X',\y), \yte)] \\
%&= \mathbb{E}_{\widetilde{\X}\sim p_{\textrm{in}}, \y\sim q(\cdot|\widetilde{\X})}[\mathcal{L}^n(\yh_n(\widetilde{\X}_{1:n},\y_{1:n-1}), y_n)] \\
&= L_{\yh}(q),
\end{align*}

where $\Xfull' = \Xfull \U^\top$

(we also used $\X'=\X\U^{\top}$ and $\Xte' = \U \Xte$),
and where the fourth line uses rotational symmetry of $p_{\textrm
in}$ and rotational invariance of $\yh$. Therefore, presented with the problem $\mathring{q}$, the algorithm will achieve expected loss $L_{\yh}(q)$ regardless of the value of $\U$. This implies that $L_{\yh}(q)$ cannot be less than the optimal value $L_{\mathcal{B}}(\mathring{q})$.

§ PROOF OF LOWER BOUND THEOREM <REF>

Before we give the proof observe that one can equivalently think of the $m$ stacked copies of $\H$ as just a single copy of $$̋ with the variance of label noise reduced by a factor of $m$. Indeed,
if $\y_i$ represents the label vector of size $d$ associated with the $i$-th copy of $\H$, then the training loss for any weight vector $$ is effectively
$$\|\y - \X \predw\|^2 = \sum_{i=1}^m \| \y_i - \H \predw\|^2
= m \|\bar{\y} - \H \predw\|^2 + \sum_{i=1}^m \|\y_i - \bar{\y}\|^2,$$
where $ = m^-1 ∑_i _i ∼N(_1, σ^2/m _d)$ and the last term on the right-hand side does not depend on $$.
Therefore, by varying $m$, we effectively alter the level of noise the algorithms encounter.
We also note that our model is equivalent to the Gaussian Sequence Model (see, for example, Rigollet and Hütter, 2023),
given by $= θ + ξ$ (with $θ = _1$ and $ξ ∼𝒩(0,σ^2/m _d)$), a commonly used model for analyzing nonparametric and high dimensional statistical problems.

\vspace{1mm}
\begin{proof}
We consider the rotated observational model described in the main text, which can be defined as follows. Let $\w \in \mathbb{R}^d$
be a weight vector drawn uniformly from a unit sphere
$\mathcal{S}^{d-1} = \{\w \in \mathbb{R}^d \colon \|\w\|=1 \}$.
The algorithm is given data set $(\X,\y)$ with $\X = [\H ; \ldots ; \H]$ being
$m$ copies of $\H=\sqrt{d} \V$, and $\y = \H \w + \bm{\xi}$, where $\bm{\xi} \sim N(\bm{0},
\sigma^2 \I_{dm})$ is a vector of Gaussian i.i.d. noise variables, each having zero
mean and variance $\sigma^2$. Given $\y$, the algorithm is supposed to produce a vector of predictions
$\predy\in \mathbb{R}^d$ and is evaluated by means of the squared error,
$e(\predy | \w) = \frac{1}{d} \|\predy - \H \w\|^2$.
We first note that without loss of generality, the algorithm produces a weight vector
$\predw$, based on which the predictions are generated, $\predy = \H \predw$; this is due
to the fact that $\H$ is invertible (as multiplicity of an orthogonal matrix), so
for every $\predy$, one can have a corresponding weight vector $\predw = \H^{-1} \predy$.
Thus, using \eqref{eq:error_w} the error can be equivalently written as
\[
e(\predw | \w) = \|\predw - \w\|^2.
\]
It is well-known (see, e.g., Berger, 1985) that the expected squared error,
$\EE_{\w,\bm{\xi}}\left[e(\predw | \w)\right]$ (with expectation with
respect to the prior and the label noise) is minimized by the \emph{posterior mean}
$\predw^{\star} = \EE_{\w|\y} [ \w]$, that is the mean value of $\w$ with respect to the posterior distribution
$q(\w|\X,\y)$. % (since $\X$ is fixed, we drop it in our notation for the sake of clarity).
Even though the posterior mean does not have a nice analytic form, we can still lower bound
its expected squared error using a technique borrowed from [Marchand, 1993, Dicker, 2016]. Let
$\predw_{RR}$ be the ridge regression estimator:
\[
\predw_{RR} = (\X^\top \X + \sigma^2 d \I)^{-1} \X^\top \y,
\]
which is the posterior mean itself (and thus optimal)
when the prior over $\w$ is Gaussian with zero mean and covariance $\frac{1}{d} \I$ (the covariance is multiplied by factor $d^{-1}$ to have $\EE[\|\w\|^2] = \EE[\w \w^\top] = \tr(\frac{1}{d}\I) = 1$ as in the unit sphere prior case). Even though the Gaussian prior differs from the uniform prior over a unit sphere, it turns out that the RR predictor has error only slightly larger than that of the optimal Bayes predictor
\begin{lemma}
\[
\EE_{\w,\bm{\xi}}\left[e(\predw_{RR} | \w)\right] \le
\EE_{\w,\bm{\xi}}\left[e(\predw^{\star} | \w)\right]
+ \frac{1}{d} \frac{\sigma^2}{\sigma^2 + m}.
\]
\label{lem:RR}
\end{lemma}
\begin{proof}
Dicker, 2016 considered a Bayesian setting similar to ours, with $\sigma^2 = 1$ and
$\w$ distributed uniformly over $\tau \mathcal{S}^{d-1} = \{\w \colon \|\w\|=\tau\}$.
To account for this setting, we note that in our setup,
\[
\sigma^{-1} \y = \X^\top (\sigma^{-1} \w) ~+ \underbrace{\sigma^{-1} \bm{\xi}}_{\sim N(\bm{0}, \I_{dm})},
\]
so that we can set $\tau = \sigma^{-1}$ and assume unit variance of the noise. Their `oracle ridge estimator' is thus $\predw_{RR}$. Furthermore, since $\|\predw - \w\|^2 =
\sigma^2 \|\sigma^{-1} \predw - \sigma^{-1} \w\|^2$, we need to multiply their bound by $\sigma^2$. We use their Theorem 2 (adapted to the modifications stated above):
\begin{theorem}{(Theorem 2 by Dicker, 2016, Theorem 3.1 by Marchand, 1993)}
Let $n = md$ and let $s_1 \ge \ldots \ge s_d$ denote the eigenvalues of $n^{-1} \X^\top \X$.
\[
\EE_{\w,\bm{\xi}}\left[e(\predw_{RR} | \w)\right]
\le \EE_{\w,\bm{\xi}}\left[e(\predw^{\star} | \w)\right]
+ \frac{\sigma^2}{d} \frac{s_1}{s_d} \mathrm{tr} \left\{ (\X^\top \X + d \sigma^2 \I_n)^{-1}\right\}
. \]
\end{theorem}
Since $\X^\top \X = n \I$, we have $s_1 = s_d = 1$ and
$(\X^\top \X +  \I_n)^{-1} = (n + d \sigma^2)^{-1} \I_n$, and thus,
\[
\EE_{\w,\bm{\xi}}\left[e(\predw_{RR} | \w)\right]
\le \EE_{\w,\bm{\xi}}\left[e(\predw^{\star} | \w)\right]
+ \frac{\sigma^2}{n + d \sigma^2}
=\EE_{\w,\bm{\xi}}\left[e(\predw^{\star} | \w)\right]
+ \frac{1}{d} \frac{\sigma^2}{m + \sigma^2}.
\]

\vspace{-9mm}
\end{proof}
\vspace{1mm}
Now, we compute the expected error of $\predw_{RR}$. Since $\X^\top \X = md \I$, we get
\[
\predw_{RR} = \frac{1}{md + \sigma^2 d} \X^\top \y
= \frac{1}{md + \sigma^2 d} \X^\top (\X \w + \bm{\xi})
= \frac{md \w + \X^\top \bm{\xi}}{md + \sigma^2 d},
\]
and thus
\begin{align*}
e(\predw_{RR}|\w) &= \|\predw_{RR} - \w\|^2
= \left| \frac{\X^\top \bm{\xi} - \sigma^2 d \w}{md + \sigma^2 d} \right\|^2
\\ &= \frac{\| \X^\top \bm{\xi}\|^2}{(md + \sigma^2 d)^2} - \frac{md \w{}^\top \X^\top \bm{\xi}}{(md + \sigma^2 d)^2} + \frac{\sigma^4 d^2 \|\w\|^2}{(md + \sigma^2 d)^2}.
\end{align*}
We take the expectation over $\w$, under which the middle term in the last line
vanishes (as $\EE[\w] = \bm{0}$ over a unit sphere) and use $\|\w\|=1$ to get
\[
\EE_{\w} [e(\predw_{RR}|\w)] =
\frac{\| \X^\top \bm{\xi}\|^2}{(md + \sigma^2 d)^2}
+ \frac{\sigma^4 d^2}{(md + \sigma^2 d)^2}.
\]
We further take an expectation over $\bm{\xi}$ and use
\[
\EE[\|\X^\top \bm{\xi}\|^2] = \EE[\tr(\X^\top \bm{\xi} \bm{\xi}^\top \X)]
= \tr(\X^\top \EE[\bm{\xi} \bm{\xi}^\top] \X) = \tr(\X^\top \X) = md \tr(\I) = md^2,
\]
to get:
\[
\EE_{\w, \bm{\xi}} [e(\predw_{RR}|\w)] =
\frac{m d^2}{(md + \sigma^2 d)^2}
+ \frac{\sigma^4 d^2}{(md + \sigma^2 d)^2}.
= \frac{\sigma^2 d (md + \sigma^2 d)}{(md + \sigma^2 d)^2}
= \frac{\sigma^2}{m + \sigma^2}.
\]
Using this together with Lemma \ref{lem:RR} gives the lower bound on the Bayes optimal
predictor in the rotated observational model.
\[
\EE_{\w, \bm{\xi}} [e(\predw^{\star}|\w)]
\ge \frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}.
\]
We now use the fact that
that $\H = \sqrt{d} \V$
with orthogonal matrix $\V$ of size $d \times d$, drawn uniformly at random (with respect to Haar measure), so that our input distribution is rotation symmetric. This means that we can apply Theorem 1 and conclude that any rotation invariant algorithm has the expected error
at least $\frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}$ on the original problem,
that is for $\w=\e_1$.
%; all the Bayesian results below are, however, conditioned on a fixed choice of $\V$, and the distribution of $\V$ becomes important only in the last step of the proof.
%Since this lower bounds the optimum error in the Bayesian setting, we conclude that
%\emph{no predictor can achieve expected error lower than}
%$\frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}$.
%So far we kept $\H = \sqrt{d} \V$ fixed. Let us now use the fact that $\V$ is drawn
%uniformly from the Haar measure, and take the expected error with respect to the random choice of $\V$. Since our lower bound on the optimum error does not depend on $\H$, we trivially
%\EE_{\w, \bm{\xi}, \V} [e(\predy|\w, \V)]
%\ge \frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}
%for \emph{any predictor} $\predy = \predy(\y, \V)$, where we explicitly stated the (possible) dependence of error on $\V$. In particular, this applies to every rotation-invariant predictor.
%Following Theorem \ref{t:gen_lower_bound}, a rotation invariant predictor has \emph{the same expected error $\EE_{\bm{\xi}, \V} [e(\predy|\w, \V)]$ for every $\w$}.
%Indeed take any two $\w_1$ and $\w^{\star}_2$. As both lie on a unit sphere, there exists an orthogonal matrix $\U$ such that $\w^{\star}_2 = \U \w^{\star}_1$. Now, we claim
%that for any fixed $\bm{\xi}$, the error of the algorithm on $(\V, \w^{\star}_1)$ (call it Problem 1) is the same as the error of the algorithm on $(\V \U^\top, \w^{\star}_2)$ (call it Problem 2). Indeed, this follows because:
%\item The training input matrix in Problem 1 is $\X = [\H \ldots \H] = \sqrt{d} [\V \ldots \V]$, so it becomes $\X \U^\top$ in Problem 2.
%\item Labels in both problems are the same, as
%$\y = \X \w^{\star}_1 + \bm{\xi} = \X \U^\top \U \w^{\star}_1 + \bm{\xi}
%= (\X \U^\top) \w^{\star}_2 + \bm{\xi}$.
%\item The predictions in both problems are the same, as both the training and the test data in Problem 1 and Problem 2 are rotated versions of each other, and the algorithm is rotation-invariant (see Equation \ref{eq:def_rot_inv}).
%\item The weight vector induced by predictions $\predw = \H^{-1} \predy$ in Problem 1 becomes
%$(\H \U^\top)^{-1} \predy = \U \predw$ in Problem 2
%\item The error of the algorithm in Problem 1, $\|\predw - \w^{\star}_1\|^2$, and in Problem 2,
%$\|\U \predw - \w^{\star}_2\|^2  = \|\U \predw - \U \w^{\star}_1\|^2 = \|\predw - \w^{\star}_1\|^2$ are the same.
%Since $\w^{\star}_1, \w^{\star}_2$ were arbitrary, we can thus state
%e(\predy|\w^{\star}, \V) = e(\predy|\U \w^{\star}, \V \U^\top),
%for any orthogonal $\V, \U$ and any $\w^{\star}$ over the unit sphere. Now, since
%$\V$ is a random unitary matrix draw from the Haar measure, the distribution of $\V$ and
%$\V \U^\top$ for any fixed orthogonal $\U$ are the same. This means that
%\EE_{\bm{\xi},\V} [e(\predy|\w^{\star}, \V)]
%= \EE_{\bm{\xi},\V} [e(\predy|\U \w^{\star}, \V \U^\top)]
%= \EE_{\bm{\xi},\V} [e(\predy|\U \w^{\star}, \V)].
%Since $\U$ were arbitrary, the expected error of the algorithm is the same for any
%$\w^{\star}$, which means that the error is
%at least $\frac{d-1}{d} \frac{\sigma^2}{\sigma^2 + m}$ for any $\w^{\star}$.
%In particular, this holds for $\w^{\star} = \e_1$.
%Let $\H = \sqrt{d} \V$
%with orthogonal matrix $\V$ of size $d \times d$, drawn uniformly at random (with respect to Haar measure); all the Bayesian results below are, however, conditioned on a fixed choice of $\V$, and the distribution of $\V$ becomes important only in the last step of the proof.
\end{proof}
Note that the proof would significantly simplify if we
assumed from the start that the target weight vector
$\w$ is generated from a Gaussian distribution
$N(\bm{0}, \frac{1}{d} \I_d)$ rather than from a unit
sphere $\mathcal{S}^{d-1}$ (both priors give unit squared
norm of $\w$ on expectation), as the Bayes
predictor would be exactly the RR predictor, giving even a
better lower bound of $\frac{\sigma^2}{\sigma^2 + m}$,
without the need to apply results from
Marchand, 1993, Dicker, 2016. This would, however, result in a
random norm of the sparse target vector. In our proof we opted for a bound with a fixed, unit norm of $\w$.

\section{Lower bound for i.i.d. sampling}
\label{app:lower_bound_iid}
Instead of observing copies of the same matrix $\H$, one can instead consider a standard linear model framework, where individual inputs
$\x$ are drawn i.i.d. from some input distribution $p_{\mathrm{in}}(\x)$. The associated labels are given by $y = \x_i^\top \e_1 + \xi$, where $\xi \sim \mathcal{N}(0,\sigma^2)$.
The algorithm observes $n$ such i.i.d. samples $(\X,\y)$ and produces a weight vector $\predw$, which is then evaluated by the squared error $\|\predw - \e_1\|^2$. Assume $p_{\mathrm{in}}(\x)$ is rotationally symmetric with covariance $\EE[\x\x^\top] = \I_d$ (so that the expected input size matches that from our previous setup); for instance, we could have $\x \sim \mathcal{N}(\bm{0}, \I_d)$. Employing our Bayesian argument from section \ref{s:gen_lower_bound},
one can establish a similar lower bound for any rotation-invariant algorithm.
Note that for large $d$, the bound is essentially the same
as that of Theorem \ref{thm:lower_bound}.
The proof makes use of information-theoretic arguments by Gerchinovitz et~al., 2020.
\begin{theorem}
Let the rows of $\X$ be drawn i.i.d. from rotationally symmetric $p_{\mathrm{in}}$ with covariance $\I_d$, and the labels given by $\y = \X \e_1 + \bm{\xi}$,
$\bm{\xi} \sim \mathcal{N}(\bm{0}, \sigma^2 \I_n)$. Assume $d \ge 12$.
Then, the expected error (with a random choice of the sample) of any rotation-invariant learning algorithm is at least
\[
\EE[e(\predw)]  \ge \left(\frac{\sigma^2}{e (\sigma^2 + m / 2)}\right)^{d/(d-1)}.
\]
\label{thm:lower_bound_iid}
\end{theorem}
%\subsection{Lower bounds}
%We always assume the linear model $\y = \x^\top \w + \xi$, with
%$\xi \sim N(0,\sigma^2)$.
%\paragraph{Arbitrary input distribution, Gaussian target}
%Assume inputs $\x$ are generated according to an arbitrary rotation-invariant distribution,
%while the target $\w \sim N(\bm{0}, d^{-1} \I)$. The covariance of the target distribution
%was chosen to have:
%\EE[\|\w\|^2] = \tr(\EE[\w\w^\top]) = \tr(d^{-1} \I_d) = 1,
%that is, the target has unit length on expectation. In the Bayesian framework, the
%ridge regression estimator is the Bayes optimal predictor with its error given by
%e(\w_{RR}) = \sigma^2 \EE_{\X}\left[\tr\left\{(\X^\top \X + d \sigma^2 \I_n)^{-1} \right\} \right],
%where the expectation is with respect to random data $\X$.
%Let us now fix $S$ and try to calculate the minimum over
%$F(\X^\top \X) = \tr\left\{(\X^\top \X + d \sigma^2 \I_n)^{-1}\right\}$ under constraint $\tr(\X^\top \X) = S$.
%Clearly $F$ is a function of only the eigenvalues of $\X^\top \X$, $\lambda_1,\ldots,\lambda_d$ and we have $\sum_i \lambda_i = S$. Explicitly writing:
%F(\X^\top \X) = \sum_{i=1}^d \frac{1}{\lambda_i + d \sigma^2},
%we see that the this $F$ is convex and symmetric in $\lambda_i$'s, and thus is minimized
%for $\lambda_1 = \ldots = \lambda_d$. That is, $\lambda_i = S/d$ for all $i$, and
%$F(\X^\top \X) \ge \frac{d}{(S/d) + d \sigma^2} = \frac{d^2}{S + d^2 \sigma^2}$.
%From \eqref{eq:RR_Bayes_risk} we thus get
%e(\w_{RR}) \ge d^2 \sigma^2 \EE_{S}\left[ \frac{1}{S + d^2 \sigma^2} \right],
%where the expectation is now with respect to a random variable $S = \tr(\X^\top \X)$.
%Using Jensen's inequality for convex function $f(x) = 1/(x + d^2 \sigma^2)$, we have
%e(\w_{RR}) \ge \frac{\sigma^2 d^2}{\EE[\tr(\X^\top \X)] + d^2 \sigma^2}.
%Let $\tau^2 = \EE[\|\x\|^2]$ be the total variance of the input distribution.
%We have:
%\EE[\tr(\X^\top \X)] = \EE\left[\tr\left(\sum_{i=1}^n \x_i \x_i^\top\right)\right]
%= \sum_{i=1}^n \EE\left[\tr\left(\x_i \x_i^\top\right)\right]
%= \sum_{i=1}^n \EE\left[\|\x_i\|^2\right] = n \tau^2 = m d \tau^2,
%which finally gives
%e(\w_{RR}) \ge \frac{\sigma^2 d^2}{m d \tau^2 + d^2 \sigma^2} = \frac{\sigma^2}{m \tau^2/d + \sigma^2}.
%Then, since $\w_{RR}$ was optimal, any other predictor has error at least $\frac{\sigma^2}{m \tau^2/d + \sigma^2}$.
%Now, assume that the distribution of $\x$ is such that $\tau^2 = \EE[\|\x\|^2]=d$, so that
%we match the Hadamard case (as well as our previous proof) and ensure that labels $y = \x^\top \w + \epsilon$ are of order $1+\sigma^2$. We then get a lower bound for any predictor at least $\frac{\sigma^2}{m + \sigma^2}$. We see that our original bound on the $m$ stacked copies of the $\sqrt{d} \U$ was
%actually optimal.
%\paragraph{Arbitrary input distribution, unit target.}
\begin{proof}
We consider a Bayesian setup, in which the target weight vector is drawn
uniformly from a unit sphere,
$\w \sim \mathcal{S}^{d-1}$,
with $\mathcal{S}^{d-1} = \{\w \in \mathbb{R}^d \colon \|\w\|=1\}$.
We will show a lower bound on the error of \emph{any} predictor $\predw$,
$\EE_{\w,\bm{\xi}}\left[\|\predw - \w\|^2\right]$,
by exploiting the fact
that it is not less than the error of the Bayes optimal predictor $\predw^{\star}$ (the posterior mean), followed
by a lower bound on the error of $\predw^{\star}$:
\begin{equation}
\EE_{\w,\bm{\xi}}\left[\|\predw^{\star} - \w\|^2\right] \ge \left(\frac{\sigma^2}{e (\sigma^2 + m / 2)}\right)^{d/(d-1)}.
\label{eq:lower_bound_to_prove_iid}
\end{equation}
Using the rotational symmetry of the input distribution $p_{\mathrm{in}}$, and giving the same arguments as in the proof of Theorem \ref{thm:lower_bound} (or Section \ref{s:gen_lower_bound}), this will imply
the same lower bound \eqref{eq:lower_bound_to_prove_iid} for any rotation-invariant algorithm with a sparse target vector
$\w = \e_1$.

%Let the $n$ rows of $\X$ be drawn i.i.d. from rotationally symmetric distribution $p_{\mathrm{in}}$ with covariance $\I_d$, and the labels given by $\y = \X \e_1 + \bm{\xi}$,
%$\bm{\xi} \sim \mathcal{N}(\bm{0}, \sigma^2 \I_n)$.
%Then the expected error (with a random choice of the sample) of any rotation-invariant learning algorithm is at least
%\EE[e(\predw)]  \ge \left(\frac{\sigma^2}{e (\sigma^2 + m / 2)}\right)^{d/(d-1)}.

%Let $\predw^{\star}$ be the Bayes optimal predictor, which is the posterior mean, so
%that for any other predictor $\predw = \predw(\X,\y)$,
%\EE_{\w,\y|\X}[\|\predw - \w\|^2] \ge \EE_{\w,\y|\X}[\|\predw^{\star} - \w\|^2]
Unfortunately, for prior uniform over a sphere, $\predw^{\star}$ does not have a closed form. We will, however, use a continuous version of Fano's inequality [Gerchinovitz et~al., 2020] to lower bound its Bayes risk.

From now on, we condition our analysis on $\X$ drawn i.i.d. from the input distribution.
We lower-bound the Bayes risk using Markov's inequality: for any $\epsilon > 0$,
\begin{equation}
\EE[\|\predw^{\star} - \w\|^2] \ge \epsilon P(\|\predw^{\star} - \w\|^2 \ge \epsilon)
= \epsilon \left(1 - P(\|\predw^{\star} - \w\|^2 < \epsilon) \right).
\label{eq:bound_sphere}
\end{equation}
where the expectation and the probability is with respect to random $\y, \w$ (conditioned on $\X$). Let $P_{\w} = \mathcal{N}(\X \w, \sigma^2 \I_n)$ be the distribution of $\y$ given $\w$ (and $\X$). Given $\w$, define event
$\mathcal{A}_{\w} = \{ \|\predw^{\star} - \w\|^2 < \epsilon\}$. We can
\[
P\left(\|\predw^{\star} - \w\|^2 < \epsilon \right) = \EE_{\w} \left[P_{\w}(\mathcal{A}_{\w})\right]
\]
Using techniques from Gerchinovitz et~al., 2020,
we have for any distribution $Q$ over $\y$,
\begin{equation}
\EE_{\w} \left[P_{\w}(\mathcal{A}_{\w})\right]
\le \EE_{\w}\left[Q(\mathcal{A}_{\w})\right]
+ \sqrt{\left(\EE_{\w}\left[Q(\mathcal{A}_{\w})\right]\right)
\left(\EE_{\w} \left[\chi^2(P_{\w}, Q) \right] \right)}
\label{eq:Gerchinovitz_etal2020}
\end{equation}
where $\chi^2(P,Q) = \EE_{Q}\left[\left(\frac{\mathrm{d}P}{\mathrm{d}Q}\right)^2\right]$ is the $\chi^2$-divergence between $P$ and $Q$.
We now choose $Q = \mathcal{N}(\bm{0}, \SSigma)$
with $\SSigma = \sigma^2 \I_n + \frac{1}{d} \X\X^\top$.
It follow from the Gaussian integration that
\[
\chi^2(P_{\w},Q) =
\sqrt{\det \left(\sigma^{-2} \SSigma^2 (2 \SSigma - \sigma^2 \I_n)^{-1}\right)}
e^{\w^\top \X^\top \left(2 \SSigma - \sigma^2 \I_n\right)^{-1} \X \w} - 1.
\]
First note that
\[
\sigma^{-2} \SSigma^2 (2 \SSigma - \sigma^2 \I_n)^{-1}
= \left(\I_n + \frac{1}{d \sigma^2} \X\X^\top\right)^2 \left(\I_n + \frac{2}{d \sigma^2} \X \X^\top\right)^{-1} \preceq \I_n + \frac{1}{2 d \sigma^2} \X \X^\top.
\]
Indeed, by taking any eigenvalue $\lambda_i$ of $\frac{1}{2d \sigma^2} \X \X^\top$,
\[
\frac{(1 + \lambda_i)^2}{1 + 2 \lambda_i} \le
1 + \frac{\lambda_i^2}{1 + 2 \lambda_i} \le 1 + \frac{\lambda_i}{2}.
\]
Moreover, from the property of the determinant, that
$\det(\I_m + \bm{A} \bm{B}) = \det(\I_n + \bm{B} \bm{A})$ for $\bm{A}$ of size $m \times n$ and $\bm{B}$ of size $n \times m$, we get
\[
\det\left(\I_n + \frac{1}{2 d \sigma^2} \X \X^\top\right) = \det\left(\I_m + \frac{1}{2 d \sigma^2} \X^\top \X\right).
\]
Finally, note that
\[
\X^\top (2 \SSigma - \sigma^2 \I_n)^{-1} \X
= \X^\top \left(\sigma^2 \I_n + \frac{2}{d} \X \X^\top\right)^{-1} \X
\preceq \frac{d}{2} \X^\top(\X \X^\top)^{\dagger} \X = \frac{d}{2} \I_d,
\]
so that
\[
e^{\w^\top \X^\top \left(2 \SSigma - \sigma^2 \I_n\right)^{-1} \X \w}
\le e^{\frac{d}{2} \|\w\|^2} = e^{\frac{d}{2}}.
\]
Thus, we get:
\[
\chi^2(P_{\w},Q) \le
\sqrt{\det \left(\I_m + \frac{1}{2 d \sigma^2} \X^\top \X\right)} e^{\frac{d}{2}}
\]
For any positive-definite $d \times d$ matrix $\bm{A}$ with eigenvalues $a_1,\ldots,a_d$,
we have $\ln \det \bm{A} = \sum_i \ln a_i$. Using the concavity of negative logarithm,
we can bound:
\[
\ln \det \bm{A} = d \frac{1}{d} \sum_i \ln(a_i)
\le d \ln\left(\frac{1}{d} \sum_i a_i \right)
= d \ln \tr(\bm{A} / d),
\]
which gives $\det \bm{A} \le \tr(\bm{A} / d)^d$ and thus
\[
\chi^2(P_{\w},Q) \le
\left(e \tr\left(\frac{1}{d}\I_m + \frac{1}{2 d^2 \sigma^2} \X^\top \X\right)\right)^{\frac{d}{2}}.
\]
Since this expression does not depend on $\w$, taking expectation on both sides (with respect to $\w$) gives the same bound on $\EE_{\w}\left[\chi^2(P_{\w},Q)\right]$.
%\sqrt{\det \left(\left(\I_n + \frac{1}{d \sigma^2} \X\X^\top\right)^2 \left(\I_n + \frac{2}{d \sigma^2} \X \X^\top\right)^{-1}\right)}
%e^{\w^\top \X^\top \left(\I_n \sigma^2 + \frac{2}{d} \X \X^\top\right)^{-1} \X \w} - 1

%Using Lemma 4.1 of Gerchinovitz et~al., 2020,
%we have for any distribution $Q$ over $\y$,
%\EE_{\w} \left[P_{\w}(\mathcal{A}_{\w})\right]
%\le \frac{\EE_{\w} \left[\mathrm{KL}(P_{\w}, Q) \right] + \ln 2}{-\ln \EE_{\w}\left[Q(\mathcal{A}_{\w})\right]},
%where $\mathrm{KL}(\cdot,\cdot)$ is the Kullback-Leibler (KL) divergence. We
%now choose $Q = \mathcal{N}(\bm{0}, \frac{1}{d} \X \X^\top + \sigma^2 \I_n)$, which coincides with a marginal distribution of $\y$ if $\w$ were drawn from a Gaussian $N(\bm{0}, \I_d / d)$.
%Using the known fact that the KL divergence between two normal distributions in dimension $n$ is given by (see, e.g., \citet{Duchi_derivations})
%\mathrm{KL} \left(\mathcal{N}(\bm{\mu}_0,\SSigma_0), \mathcal{N}(\bm{\mu}_1, \SSigma_1)\right)
%= \frac{1}{2} \left\{ \tr(\SSigma_1^{-1} \SSigma_0) + (\bm{\mu}_1 \!-\! \bm{\mu}_0)^\top \SSigma_1^{-1}(\bm{\mu}_1 \!-\! \bm{\mu}_0) - n + \ln \frac{\det \SSigma_1}{\det \SSigma_0}  \right\},
%and denoting $\SSigma_Q = \frac{1}{d} \X \X^\top + \sigma^2 \I_n$ the covariance matrix of
%$Q$, we get
%\mathrm{KL} (P_{\w}, Q)
%&= \frac{1}{2} \left\{ \tr(\sigma^2 \SSigma_Q^{-1})  + \w^\top \X^\top \SSigma_Q^{-1}\X \w
%- n + \ln \det (\sigma^{-2} \SSigma_Q )\right\} \\
%&= \frac{1}{2} \left\{ \tr(\sigma^2 \SSigma_Q^{-1})  + \tr(\SSigma_Q^{-1}\X \w \w^\top \X^\top)
%- n + \ln \det (\sigma^{-2} \SSigma_Q )\right\}.
%Taking expectation over $\w$:
%\EE_{\w} \left[\mathrm{KL}(P_{\w}, Q) \right]
%&= \frac{1}{2} \bigg\{ \tr(\sigma^2 \SSigma_Q^{-1})  + \tr(\SSigma_Q^{-1}\X \underbrace{\EE_{\w} [\w \w^\top]}_{=\I_d / d} \X^\top)
%- n + \ln \det (\sigma^{-2} \SSigma_Q )\bigg\} \\
%&= \frac{1}{2} \bigg\{
%\tr(\SSigma_Q^{-1} \SSigma_Q) - n + \ln \det (\sigma^{-2} \SSigma_Q )
%= \frac{1}{2} \ln \det \left(\I_n + \frac{1}{\sigma^2 d} \X \X^\top\right) \\
%&= \frac{1}{2} \ln \det \left(\I_d + \frac{1}{\sigma^2 d} \X^\top \X\right),
%where in the last line we used the property of determinant $\det(\I_m + \bm{A} \bm{B}) = \det(\I_n + \bm{B} \bm{A})$ for $\bm{A}$ of size $m \times n$ and $\bm{B}$ of size $n \times m$.
%For any positive-definite $d \times d$ matrix $\bm{A}$ with eigenvalues $a_1,\ldots,a_d$,
%we have $\ln \det \bm{A} = \sum_i \ln a_i$. Using the concavity of negative logarithm,
%we can bound:
%\ln \det \bm{A} = d \frac{1}{d} \sum_i \ln(a_i)
%\le d \ln\left(\frac{1}{d} \sum_i a_i \right)
%= d \ln \tr(\bm{A} / d),
%which gives:
%\EE_{\w} \left[\mathrm{KL}(P_{\w}, Q) \right]
%\le \frac{d}{2} \ln \tr\left(\frac{1}{d}\I_d + \frac{1}{\sigma^2 d^2} \X^\top \X\right)
%= \frac{d}{2} \ln \left(\frac{\sigma^2 d^2 + \tr(\X^\top \X)}{\sigma^2 d^2}\right)
%We thus calculated the numerator of \eqref{eq:Gerchinovitz_etal2020}, and now we will
%upper-bound the expectation in the denominator, $\EE_{\w}[Q(\mathcal{A}_{\w})]$.

We will now upper-bound $\EE_{\w}[Q(\mathcal{A}_{\w})]$.
$\bm{1}\{C\}$ denote the Iverson bracket which is $1$ if $C$ holds and $0$ otherwise.
Using the fact that $\predw^{\star} = \predw^{\star}(\y)$ depends on $\y$ (and $\X$), but
not on $\w$, and that $\y \sim Q$ is independent of $\w$, we have
\begin{align*}
\EE_{\w}[Q(\mathcal{A}_{\w})]
&= \EE_{\w}[Q(\|\predw^{\star} - \w\|^2 < \epsilon)]
= \EE_{\w \sim \mathcal{S}^{d-1},\y \sim Q}\left[\bm{1}\left\{\|\predw^{\star}(\y) - \w\|^2 < \epsilon\right\}\right] \\
&= \EE_{\y \sim Q} \left[\EE_{\w \sim \mathcal{S}^{d-1}}\left[\bm{1}\left\{\|\predw^{\star}(\y) - \w\|^2 < \epsilon\right\}\right]\right] \\
&\le \EE_{\y \sim Q} \left[\max_{\predw} \EE_{\w \sim \mathcal{S}^{d-1}}\left[\bm{1}\left\{\|\predw - \w\|^2 < \epsilon\right\}\right]\right]
= \max_{\predw} \EE_{\w}\left[\bm{1}\left\{\|\predw - \w\|^2 < \epsilon\right\}\right],
\end{align*}
where in the last equality we used the fact that the term inside the expectation over $\y$
does not depend on $\y$. Let us now fix $\predw$ and we will bound
$\EE_{\w}\left[\bm{1}\left\{\|\predw - \w\|^2 < \epsilon\right\}\right]$
Since the distribution
of $\w$ is rotation-invariant (uniform over a sphere), without loss of generality take
$\predw = c \e_1$ with $c > 0$ ($c=0$ would give $0$ for any $\epsilon < 1$). We have
\[
\|\predw^{\star} - \w\|^2 = (c-w_1)^2 + \sum_{j=2}^d w_j^2 =
(c-w_1)^2 + 1 - w_1^2 = 1 + c^2 - 2 c w_1,
\]
and thus
\[
\EE_{\w}\left[\bm{1}\left\{\|\predw - \w\|^2 < \epsilon\right\}\right]
= \EE_{\w} \left[ \bm{1}\left\{1 + c^2 - 2 c w_1 < \epsilon\right\} \right]
= P\left( w_1 > \frac{1+ c^2 - \epsilon}{2c} \right).
\]
Since the probability is nonincreasing in $\frac{1+c^2-\epsilon}{2c}$ we take $c$
which minimizes this quantity to get an upper bound. By inspecting the derivative
\[
\min_c \frac{1+c^2-\epsilon}{2c}
= \frac{1}{2} \min_c c + \frac{1-\epsilon}{c}
~\stackrel{c=\sqrt{1-\epsilon}}{=}~ \sqrt{1-\epsilon},
\]
we have
\[
\EE_{\w}\left[\bm{1}\left\{\|\predw - \w\|^2 < \epsilon\right\}\right]
= P\left(w_1 > \sqrt{1-\epsilon}\right) = \frac{1}{2} P\left(w_1^2 > 1 - \epsilon\right),
\]
where we used the fact that distribution of $w_1$ is symmetric around zero.
To determine the distribution of $w_1^2$ we use the fact that drawing $\w$ uniformly
over a unit sphere amounts to drawing $\u \sim \mathcal{N}(\bm{0}, \I_d)$ and setting
$\w = \frac{\u}{\|\u\|}$. Thus, $w_1^2 = \frac{x}{x + y}$, where
$x \sim \chi^2(1) = \mathrm{Gamma}(\alpha = 1/2, \beta = 1/2)$ and
$y \sim \chi^2(d-1) =  \mathrm{Gamma}(\alpha = (d-1)/2, \beta = 1/2)$ and $x,y$ independent,
$w_1^2$ is distributed according to beta distribution $B(\alpha = 1/2, \beta = (d-1)/2)$.
This gives
\begin{align*}
P\left(w_1^2 > 1 - \epsilon\right)
&= \frac{\Gamma(d/2)}{\Gamma(1/2) \Gamma((d-1)/2)} \int_{1-\epsilon}^1 x^{-1/2} (1-x)^{-(d-1)/2}
\mathrm{d} x \\
&\le \frac{\Gamma(d/2)}{\Gamma(1/2) \Gamma((d-1)/2)} (1-\epsilon)^{-1/2} \int_{1-\epsilon}^1 (1-x)^{(d-1)/2-1} \mathrm{d} x \\
&= \frac{\Gamma(d/2)}{\Gamma(1/2) \Gamma((d-1)/2)} \frac{2(1-\epsilon)^{-1/2} \epsilon^{(d-1)/2}}{d-1}
\end{align*}
Assuming $\epsilon < 1/2$, we have $2(1-\epsilon)^{-1/2} \le 2\sqrt{2}$. We also use
Gautschi's inequality which states that for any $x > 0$ and any $s \in (0,1)$,
$\frac{\Gamma(x+1)}{\Gamma(x+s)} \le (x+1)^{1-s}$. Taking
$x = \frac{d}{2}-1$ and $s=\frac{1}{2}$, we can bound
$\frac{\Gamma(d/2)}{\Gamma((d-1)/2)} \le \sqrt{d/2} \le \sqrt{d-1}$, where we used $d \ge 2$. This, together
with $\Gamma(1/2) = \sqrt{\pi}$ allows us to bound
\[
P\left(w_1^2 > 1 - \epsilon\right)
\le \frac{2 \sqrt{2} \epsilon^{(d-1)/2}}{\sqrt{\pi} \sqrt{d-1}}.
\]
Thus, $\EE_{\w}[Q(\mathcal{A}_{\w})]$ can be bounded by
\[
\EE_{\w}[Q(\mathcal{A}_{\w})] = \frac{1}{2} P\left(w_1^2 > 1 - \epsilon\right)
\le \frac{\sqrt{2} \epsilon^{(d-1)/2}}{\sqrt{\pi} \sqrt{d-1}}
\]
Plugging all bounds to \eqref{eq:Gerchinovitz_etal2020}
\[
\EE_{\w}\left[P_{\w}(\mathcal{A}_{\w})\right]
\le \frac{\sqrt{2} \epsilon^{(d-1)/2}}{\sqrt{\pi(d-1)}}
+ \sqrt{\sqrt{2} \frac{\epsilon^{(d-1)/2}}{\sqrt{(\pi(d-1)}}
\left(e \tr\left(\frac{1}{d}\I_m + \frac{1}{2 d^2 \sigma^2} \X^\top \X\right)\right)^{\frac{d}{2}}}.
\]
Now, take
\[
\epsilon^{-1} = \left(e
\tr\left(\frac{1}{d}\I_m + \frac{1}{2 d^2 \sigma^2} \X^\top \X\right)\right)^{d/(d-1)}.
\]
Note that $\epsilon \le 1/e$.
This gives
\[
\EE_{\w}\left[P_{\w}(\mathcal{A}_{\w})\right]
\le \frac{\sqrt{2} e^{-(d-1)/2}}{\sqrt{\pi(d-1)}}
+ \sqrt{\frac{\sqrt{2}}{\sqrt{\pi(d-1)}}} \le \frac{1}{2},
\]
whenever $d \ge 12$. We thus have shown
\[
P\left(\|\predw^{\star} - \w\|^2 \le \epsilon\right)
\le \frac{1}{2},
\]
which, using \eqref{eq:bound_sphere}, gives:
\[
\EE_{\w,\y|\X}[\|\predw^{\star} - \w\|^2] \ge \frac{\epsilon}{2}.
\]
To get a bound with respect to a random choice of $\X$, note that
function $x \mapsto x^{-d/(d-1)}$ is convex, and thus:
\begin{align*}
\EE[\|\predw^{\star} - \w\|^2] &\ge \frac{1}{2} \EE[\epsilon]
\ge \frac{1}{2}
\left(e
\tr\left(\frac{1}{d}\I_m + \frac{1}{2 d^2 \sigma^2} \EE[\X^\top \X]\right)\right)^{d/(d-1)} \\
&= \left(e \left(1 + \frac{n \tau^2}{2 d \sigma^2}\right)\right)^{-d/(d-1)}
= \left(\frac{\sigma^2}{e (\sigma^2 + n \tau^2 / (2d))}\right)^{d/(d-1)}
\end{align*}
\end{proof}
%Thus, the denominator in \eqref{eq:Gerchinovitz_etal2020} is bounded by
%-\ln \EE_{\w}\left[Q(\mathcal{A}_{\w})\right]
%\ge \frac{d-1}{2} \ln (1/\epsilon) + \frac{1}{2} \ln \frac{(d-1) \pi}{2}
%Given $\X$ and $\y$, let
%\predw = (\X^\top \X)^{\dagger} \X^\top \y
%be the least squares estimator.
%We now compute the optimal Bayesian predictor $\predw^{\star}$,
%which is the posterior mean. As the prior is uniform over $\mathcal{S}^{d-1}$,
%the posterior distribution is given by
%q(\w | \X,\y) \propto \exp\left\{-\frac{1}{2} \|\y - \X \w\|^2 \right\}
%for any $\w \in \mathcal{S}^{d-1}$.
%Doing algebraic transformations:
%\frac{1}{2} \|\y - \X \w\|^2
%&= \frac{1}{2} \|\y - \X \predw + \X \predw - \X \w\|^2 \\
%&= \frac{1}{2} \|\y - \X \predw\|^2 + (\X \predw - \X \w)^\top (\y - \X\predw)
%+ \frac{1}{2} \|\X(\predw - \w)\|^2 \\
%&= \frac{1}{2} \|\y - \X \predw\|^2 + (\predw - \w)^\top \X^\top (\y - \X\predw)
%+ \frac{1}{2} (\predw - \w)^\top \X^\top \X (\predw - \w) \\
%&= \frac{1}{2} \|\y - \X \predw\|^2 + \frac{1}{2}(\predw - \w)^\top \X^\top \X (\predw - \w) .
%The term in the middle disappears, because
%\X^\top (\y - \X\predw)= \X^\top \y - \X^\top \X (\X^\top \X)^{\dagger} \X^\top \y
%= \X^\top \y - \X^\top \y = \bm{0}.
%Thus, we can write down the posterior as
%q(\w | \X,\y) &\propto \exp\left\{-\frac{1}{2} \|\y - \X \predw\|^2 -
%\frac{1}{2}(\predw - \w)^\top \X^\top \X (\predw - \w)  \right\} \\
%&\propto \exp\left\{ -\frac{1}{2} (\predw - \w)^\top \X^\top \X (\predw - \w)  \right\},
%where we used the fact that the first term does not depend on $\w$, so it will disappear
%after normalization of the distribution.
%We conclude that the posterior is a function of $\X^\top \X$ and $\predw$, and thus so
%is the Bayes predictor $\predw^{\star} = \predw^{\star}(\X^\top \X, \predw)$.
%Is it any helpful?

%We now compute the expected Bayes error,
%\EE_{\w, \bm{\xi}, \X}[e(\predw^{\star})] =
%\EE_{\w, \bm{\xi}, \X}\left[\|\predw^{\star} - \w\|^2\right]
%where the expectation is with respect to the prior distribution of $\w$,
%the input distribution of $\X$, and the noise variables $\bm{\xi} \sim N(\bm{0}, \sigma^2 \I_n)$.
%Note that $\bm{\xi}, \X, \w$ are all \emph{independent of each other}: indeed, the input distribution is fixed and does not depend on $\w$, while $\bm{\xi}$ are independent noise
%variables by definition.
%We have $\y = \X \w + \bm{\xi}$, which after multiplication from the left by $(\X^\top \X)^{\dagger} \X^\top$ gives
%\predw = (\X^\top \X)^{\dagger} \X^\top \X \w + (\X^\top \X)^{\dagger} \X^\top \bm{\xi}
%= \bm{P} \w + \bm{\epsilon},
%where $\bm{\epsilon} = (\X^\top \X)^{\dagger} \X^\top \bm{\xi}$ and
%$\bm{P} = (\X^\top \X)^{\dagger} \X^\top \X$ is a projector onto the span of $\X$.
%We now transform
%\|\predw^{\star} - \w\|^2
%=  \|\predw^{\star} - \w\|^2

\section{Upper bound for EG$^\pm$: Proof of Theorem \ref{thm:EG}}
\label{a:EGpm-proof}

Here we prove that the
batch version of EG$^{\pm}$ achieves small error already after one trial, when the learning rate is set to a sufficiently large value.

The (batch) EG$^{\pm}$ algorithm keeps track of two vectors, $\v^+_t$ and $\v^-_t$, and its prediction vector is given by $\w_t = \v^+_t - \v^-_t$.
It starts with a set of weights $\v_1^+ = \v_1^- = \frac{1}{2d} \1$,  and updates according to
\[
\v^+_{t+1} \propto \v^+_t \odot e^{-\eta \nabla \emploss(\w_t)}, \qquad \v^-_{t+1} \propto \v^-_t \odot e^{\eta \nabla \emploss(\w_t)},
\]
where $\odot$ is component-wise multiplication, and the normalization ensures that
$\|\v^+_{t+1}\|_1 + \|\v^-_{t+1}\|_1 = 1$, while $\emploss(\w)$ is the average total loss on the training
\[
\emploss(\w) = \frac{1}{dm} \|\X \w - \y_t\|^2.
\]
We compute the expression for the gradient:
\begin{align*}
\nabla \emploss(\w) &= \frac{2}{dm} \sum_{t=1}^m \sqrt{d} \V^\top (\sqrt{d} \V \w - \y_t)
= \frac{2}{dm} \sum_{t=1}^m \sqrt{d} \V^\top (\sqrt{d} \V (\w - \e_1) - \bm{\xi}_t) \\
&= \frac{2}{dm} \sum_{t=1}^m \left( d(\w - \e_1) + \sqrt{d} \V^\top \bm{\xi}_t \right)
= 2 (\w - \e_1) - \frac{2}{\sqrt{d}} \V^\top \avgxi,
\end{align*}
\[
\avgxi = \frac{1}{m} \sum_{t=1}^m \bm{\xi}_t ~\sim~ N\left(\0, \frac{\sigma^2}{m}\I\right),
\]
and the reduction of variance is due to averaging i.i.d. noise variables. Furthermore,
we rewrite
\begin{equation}
\nabla \emploss(\w) = 2(\w - \e_1 - \bm{\zeta}),
\label{eq:gradient_simplified}
\end{equation}
where the noise vector $\bm{\zeta} =  \frac{1}{\sqrt{d}} \V^\top \avgxi$ has distribution
\[
\bm{\zeta} ~\sim~ N \left(\0, \frac{\sigma^2}{m d} \V^\top \V\right)
~=~ N \left(\0, \frac{\sigma^2}{m d} \I\right).
\]
We also bound the error of the algorithm from above:
\[
e(\w_t) = \|\w_t - \e_1\|^2 = \|\w_t\|^2 - 2 \w_t^\top \e_1 + \|\e_1\|^2 \le 2 - 2 \w_t^\top \e_1
= 2(1-w_{t,1}) = 2(1-v_{t,1}^+ +  v_{t,1}^-).
\]
So it suffices to upper-bound $1-v_{t,1}^+$ and $v_{t,1}^-$.

Consider the weights of the batch EG$^{\pm}$ algorithm after just a \emph{single} trial,
that is $\v_2^+$ and $\v_2^-$
Using \eqref{eq:gradient_simplified} and noting that $\w_1 = \0$, and
$v_{1,i}^+ = v_{1,i}^- = \frac{1}{2d}$ for all $i$, we can concisely write $v_{2,1}^+$ and $v_{2,1}^-$ as:
\begin{equation}
v_{2,1}^+ = \frac{e^{2\eta\left(1 + \zeta_1\right)}}{Z_2},
\qquad
v_{2,1}^- = \frac{e^{-2\eta\left(1 + \zeta_1\right)}}{Z_2},
\qquad
Z_2 = \sum_{i=1}^d e^{2\eta\left(\delta_{1i} + \zeta_i\right)}
+ e^{-2\eta\left(\delta_{1i} + \zeta_i\right)},
\label{eq:v_plus_v_minus}
\end{equation}
We will now lower-bound $v_{2,1}^+$, and later use a relation which directly follows from
\eqref{eq:v_plus_v_minus}:
\begin{equation}
v_{2,1}^- = e^{-4\eta(1+\zeta_1)} v_{2,1}^+.
\label{eq:relation_v21}
\end{equation}

Using the deviation bound for zero-mean Gaussian
$z \sim N(0,\tau^2)$,
$P(|z| \ge \gamma) \le 2\exp\left\{-\frac{\gamma^2}{2\tau^2}\right\}$, we get
that for any $i=1,\ldots,d$,
\[P(|\zeta_i| \ge 1/4) \le 2 \exp \left\{- \frac{md \gamma^2}{32 \sigma^2} \right\}.
\]
Taking the union bound over $i=1,\ldots,i$, we have
\[
P(\exists i \; |\zeta_i| \ge 1/4)
\le 2d \exp \left\{- \frac{md \gamma^2}{32 \sigma^2} \right\}.
\] Denoting the probability on the right-hand side by $\delta$, we conclude that with probability at least $1-\delta$, all noise variables $\zeta_i$ are bounded by $1/4$. Let us call this event $E$, and we condition everything that follows on the fact that
$E$ happened.

Note that for any $i \ge 2$,
\[
\frac{\partial v_{2,1}^+}{\partial \zeta_i}
= - \frac{v_{2,1}^+}{Z_2} 2 \eta \left(e^{2\eta \zeta_i} - e^{-2 \eta \zeta_i} \right),
\]
which is decreasing for $\zeta_i > 0$ and increasing for $\zeta_i < 0$. So,
to lower-bound $v_{2,1}^+$, conditioning on event $E$, we set $\zeta_i = 1/4$ for all $i \ge 2$ in \eqref{eq:v_plus_v_minus} ($\zeta_i=-1/4$ would result in the same value of these weights). This gives:
\begin{align}
v_{2,1}^+ &\ge \frac{e^{2\eta(1+\zeta_1)}}{e^{2\eta(1+\zeta_1)} + e^{-2\eta(1+\zeta_1)}
+ (d-1) (e^{\eta/2} + e^{-\eta / 2})} \nonumber \\
&= \frac{1}{1 + e^{-4\eta(1+\zeta_1)} + e^{-2\eta(1+\zeta_1)}(d-1)(e^{\eta/2} + e^{-\eta / 2})} \nonumber \\
&\ge \frac{1}{1 + e^{-4\eta(1-1/4)} +  e^{-2\eta(1-1/4)}(d-1)(e^{\eta/2} + e^{-\eta / 2})}
\nonumber \\
&\ge \frac{1}{1 + e^{-3\eta} + e^{-3/2 \eta} 2 (d-1)e^{\eta/2}}
\nonumber \\
&\ge \frac{1}{1 + e^{-3\eta} + 2(d-1) e^{-\eta}}
\ge \frac{1}{1+(2d-1) e^{-\eta}}
\label{eq:lower_bound_v21_plus}
\end{align}
This gives:
\[
1 - v_{2,1}^+ \le \frac{(2d-1) e^{-\eta}}{1 + (2d-1) e^{-\eta}}
= \frac{1}{1 + e^{\eta}/ (2d-1)} \le (2d-1) e^{-\eta}.
\]
To upper-bound $v_{2,1}^-$, we use \eqref{eq:relation_v21}. Conditioning on $E$,
\[
v_{2,1}^- =e^{-4\eta(1+\zeta_1)} v_{2,1}^+ \le e^{-4 \eta(1+\zeta_1)}
\le e^{-2 \eta}.
\]
Thus, with probability at least $1-\delta$, the error can be bounded by:
\[
e(\w_2) \le 2(1-v_{2,1}^+ + v_{2,1}^-)
\le (2d-1)e^{-\eta} + e^{-2\eta} \le 2d e^{-\eta}
\]
To get the expected error (with respect to the training data), we bound
\begin{align*}
\EE[e(\w_2)] &= \EE[e(\w_2)|E] P(E) + \EE[e(\w_2)|E'] P(E')
\le \EE[e(\w_2)|E] + \delta \EE[e(\w_2)|E'] \\
&\le \EE[e(\w_2)|E] + 2 \delta
= 2de^{-\eta} + 8 d e^{-\frac{md}{32\sigma^2}},
\end{align*}
where we used the fact that $e(\w_2) \le 2(1-v_{2,1}^+ + v_{2,1}^-) \le 4$ as $v_{2,1}^+, v_{2,1}^- \in [0,1]$, and that maximizing convex
function $e(\w)$ give $\w$
Thus, taking sufficiently large $\eta$, we can drop the error arbitrarily
close to $8d e^{-\frac{md}{32\sigma^2}}$.\hfill $\BlackBox$
\section{Upper bound for Approximated EGU$^{\pm}$: Proof of Theorem \ref{thm:spindly}}

\label{a:ApproxEGU-upper}

We derive the Approximated EGU$^{\pm}$ algorithm defined by \eqref{eq:spindly_update_main}
as a first-order approximation of the
unnormalized Exponentiated Gradient update.

The vanilla EGU$^{\pm}$ algorithm
keeps track of two vectors, $\v^+_t$ and $\v^-_t$, and the prediction vector is given by $\w_t = \v^+_t - \v^-_t$.
Let $\beta$ denote the initial value of weights, that is $\v_1^+ = \v_1^- = \beta \1$. The weights
are updated according to
\begin{equation}
\v^{\pm}_{t+1} = \v^{\pm}_t \odot e^{\mp \eta \nabla \emploss(\w_t)}
= \beta e^{\mp \eta \sum_{j=1}^t \nabla \emploss(\w_j)}.
\label{eq:egupm_update}
\end{equation}
At every timestamp we have $\v^+_t \v^-_t = \beta^2$, which together with
$\w_t = \v^+_t - \v^-_t$, allows us to express $\v^+_t and \v^-_t$ in terms of $\w_t$:
\begin{equation}
\v^+_t = \frac{\w_t + \sqrt{\w_t^2 + 4\beta^2}}{2}, \qquad
\v^-_t = \frac{-\w_t + \sqrt{\w_t^2 + 4\beta^2}}{2}.
\label{eq:v_plus_v_minus_update}
\end{equation}
Expanding the EGU$^{\pm}$ update \eqref{eq:egupm_update} in the learning rate we get
\[
\v_{t+1}^{\pm} = \v_t^{\pm} e^{\mp\eta \nabla L(\w_t)}
= \v_t^{\pm} (1 \mp \eta \nabla L(\w_t)) + O(\eta^2).
\]
Dropping the $O(\eta^2)$ term and using \eqref{eq:v_plus_v_minus_update} gives
\[
\v_{t}^+ + \v_t^- = \sqrt{\w^2 + 4\beta^2},
\]
so that the update becomes \eqref{eq:spindly_update_main}:
\[
\w_{t+1} = \w_t - \eta\sqrt{\w_t^2 + 4\beta^2}\; \nabla L(\w_t).
\]

\noindent{\bf Proof of Theorem \ref{thm:spindly}.}
%We start with derivation of the Spindly update \eqref{eq:spindly_update_main} from
Using \eqref{eq:gradient_simplified}, $\nabla \emploss(\w_t) = 2(\w_t - \e_1 - \bm{\zeta})$,
%expanding the EG$^{\pm}$ update \eqref{eq:egupm_update} in the learning rate we get
%\v_{t+1}^{\pm} = \v_t^{\pm} e^{\mp2\eta(\w_t - \e_1 + \bm{\zeta})}
%= \v_t^{\pm} (1 \pm 2\eta(\w_t - \e_1 + \bm{\zeta})) + O(\eta^2).
%Dropping the $O(\eta^2)$ term and using \eqref{eq:v_plus_v_minus_update} gives
%\v_{t}^+ + \v_t^- = \sqrt{\w^2 + 4\beta^2},
the update \eqref{eq:spindly_update_main} becomes:
\begin{equation}
\w_{t+1} = \w_t - 2\eta\sqrt{\w_t^2 + 4\beta^2}\; (\w_t - \e_1 - \bm{\zeta})
\label{eq:spindly_update}
\end{equation}

We set $\beta=1/(2d)$, $\eta = 1/4$, and the number of steps of the algorithm
$T=4 \sqrt{d}$. Note that $d \ge 4$
Using the deviation bound for zero-mean Gaussian $z \sim N(0,\tau^2)$,
$P(z \ge -\gamma) \le \exp\left\{-\frac{\gamma^2}{2\tau^2}\right\}$, we get
that $P(|\zeta_i| \ge \gamma) \le 2 \exp \left\{- \frac{md \gamma^2}{8 \sigma^2} \right\}$.
Taking the union bound we have
$P(\exists i \; |\zeta_i| \ge \gamma) \le 2d \exp \left\{- \frac{md \gamma^2}{8 \sigma^2} \right\}$. Denoting the probability on the left-hand side by $\delta$, we can solve
for $\gamma$:
\[
\gamma = \sigma \sqrt{\frac{\ln \frac{2d}{\delta}}{md}}.
\]
This means the with probability at least $1-\delta$, $|\zeta_i| \le \gamma$ for all $i=1,\ldots,d$. Let us call this event $E$, and we condition everything what follows on the fact that
$E$ happened. Furthermore, due to our assumption that $m$ grows at least logarithmically with $d$, $m \ge 8 \sigma^2 \ln \frac{2 d}{\delta}$, we have
$\gamma \le \frac{1}{\sqrt{8d}}$.

%We note the following: if $\beta$ is of order 1, we essentially have GD. It will require
%the number of iterations $T$ of order $T\eta = 1$ to reach $1$.
%On the other hand, if $\beta= 1/2d$, we see that for small weights, the learning gets
%multiplied by $1/\beta$, which is exactly what we need (so that small learning rates for
%coordinates $i \ge 2$ remain small).
We rewrite the update \eqref{eq:spindly_update} in terms of $\s_t = \w_t - \e_1 - \bm{\zeta}$
\begin{equation}
\s_{t+1} = (1 - 2\eta\sqrt{(\s_t + \e_1 + \bm{\zeta})^2 + 4\beta^2}) \s_t
\label{eq:update_s}
\end{equation}

%We can rewrite
%\v^+_t = \beta \left( \frac{\w_t}{2\beta} + \sqrt{\left(\frac{\w_t}{2 \beta}\right)^2 + 1} \right)
%= \beta e^{\arcsinh(\w_t / (2\beta))},
%where we used $\arcsinh(x) = \ln\left(x + \sqrt{x^2 +1}\right)$, and similarly
%$\v^-_t = \beta e^{-\arcsinh(\w_t / (2 \beta))}$,
%so that
%\w_{t+1} &= \beta e^{\arcsinh(\w_t / (2\beta)) - \eta \nabla \emploss(\w_t)}
%- \beta e^{-\arcsinh(\w_t / (2\beta)) + \eta \nabla \emploss(\w_t)} \\
%&= 2 \beta \sinh \left(\arcsinh(\w_t / (2\beta)) - \eta \nabla \emploss(\w_t)\right).
%One can also define functions:
%\sinh_{\beta}(x) = 2\beta \sinh(x), \qquad
%\arcsinh_{\beta}(x) = \arcsinh(\w_t / (2\beta)),
%and write
%\w_{t+1} = \sinh_{\beta} \left(\arcsinh_{\beta}(\w_t) - \eta \nabla \emploss(\w_t)\right).

%\subsection{High-probability bound}

%\subsection{Spindly: approximation to EGU$^{\pm}$}.
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi: #1}\else
\providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

[Abbe and Boix-Adsera, 2022]
Emmanuel Abbe and Enric Boix-Adsera.
\newblock On the non-universality of deep learning: quantifying the cost of
\newblock In \emph{Advances in Neural Information Processing Systems},
volume~35, pages 17188--17201. Curran Associates, Inc., 2022.

[Abdulkadirov et~al., 2023]
Ruslan Abdulkadirov, Pavel Lyakhov, and Nikolay Nagornov.
\newblock Survey of optimization algorithms in modern neural networks, 04 2023.

[Amari and Douglas, 1998]
S.~Amari and S.C. Douglas.
\newblock Why natural gradient?
\newblock In \emph{Proceedings of the 1998 IEEE International Conference on
Acoustics, Speech and Signal Processing, ICASSP '98 (Cat. No.98CH36181)},
volume~2, pages 1213--1216 vol.2, 1998.

[Amari, 1998]
Shun-Ichi Amari.
\newblock Natural gradient works efficiently in learning.
\newblock \emph{Neural computation}, 10\penalty0 (2):\penalty0 251--276, 1998.

[Amid and Warmuth, 2020]
Ehsan Amid and Manfred~K. Warmuth.
\newblock Reparameterizing mirror descent as gradient descent.
\newblock In \emph{Proceedings of Advances in Neural Information Processing
Systems}, volume~33, pages 8430--8439, 2020.

[Berger, 1985]
James~O. Berger.
\newblock \emph{Statistical decision theory and Bayesian analysis}.
\newblock Springer, 1985.

[Bonnabel, 2013]
Silvere Bonnabel.
\newblock Stochastic gradient descent on riemannian manifolds.
\newblock \emph{IEEE Transactions on Automatic Control}, 58\penalty0
(9):\penalty0 2217--2229, 2013.

[Cesa-Bianchi et~al., 2007]
Nicolo Cesa-Bianchi, Yishay Mansour, and Gilles Stoltz.
\newblock Improved second-order bounds for prediction with expert advice.
\newblock \emph{Machine Learning}, 66\penalty0 (2-3):\penalty0 321--352, 2007.

[Chang et~al., 2023]
Peter Chang, Gerardo Dur{\`a}n-Mart{\'\i}n, Alexander~Y Shestopaloff, Matt
Jones, and Kevin Murphy.
\newblock Low-rank extended kalman filtering for online learning of neural
networks from streaming data.
\newblock \emph{arXiv preprint arXiv:2305.19535}, 2023.

[Dicker, 2016]
Lee~H. Dicker.
\newblock Ridge regression and asymptotic minimax estimation over spheres of
growing dimension.
\newblock \emph{Bernoulli}, 22\penalty0 (1):\penalty0 1--37, 2016.

[Duchi et~al., 2011]
John Duchi, Elad Hazan, and Yoram Singer.
\newblock Adaptive subgradient methods for online learning and stochastic
\newblock \emph{Journal of Machine Learning Research}, 12\penalty0
(Jul):\penalty0 2121--2159, 2011{\natexlab{a}}.

[Duchi et~al., 2011]
John Duchi, Elad Hazan, and Yoram Singer.
\newblock Adaptive subgradient methods for online learning and stochastic
\newblock \emph{Journal of machine learning research}, 12\penalty0 (7),

[Gerchinovitz et~al., 2020]
S{\'e}bastien Gerchinovitz, Pierre M{\'e}nard, and Gilles Stoltz.
\newblock {Fano’s Inequality for Random Variables}.
\newblock \emph{Statistical Science}, 35\penalty0 (2):\penalty0 178 -- 201,

[Gunasekar et~al., 2017]
Suriya Gunasekar, Blake~E. Woodworth, Srinadh Bhojanapalli, Behnam Neyshabur,
and Nati Srebro.
\newblock Implicit regularization in matrix factorization.
\newblock In \emph{Proceedings of Advances in Neural Information Processing
Systems (NeurIPS)}, pages 6151--6159, 2017.

[Jain et~al., 2007]
Prateek Jain, Brian~J. Kulis, and Inderjit~S. Dhillon.
\newblock Online linear regression using burg entropy.
\newblock Technical Report TR-07-08, The University of Texas at Austin, feb

[Kerekes et~al., 2021]
Anna Kerekes, Anna Mészáros, and Ferenc Huszár.
\newblock Depth without the magic: Inductive bias of natural gradient descent.
\newblock \emph{ArXiv}, abs/2111.11542, 2021.

[Kingma and Ba, 2015]
Diederik Kingma and Jimmy Ba.
\newblock Adam: A method for stochastic optimization.
\newblock In \emph{International Conference on Learning Representations
(ICLR)}, San Diega, CA, USA, 2015.

[Kivinen and Warmuth, 1997]
Jyrki Kivinen and Manfred~K. Warmuth.
\newblock Exponentiated gradient versus gradient descent for linear predictors.
\newblock \emph{Information and Computation}, 132\penalty0 (1):\penalty0 1--63,

[Kivinen et~al., 1997]
Jyrki Kivinen, Manfred~K. Warmuth, and Peter Auer.
\newblock The {P}erceptron algorithm versus {W}innow: linear versus logarithmic
mistake bounds when few input variables are relevant.
\newblock \emph{Artificial Intelligence}, 97\penalty0 (1-2):\penalty0 325--343,

[Lambert et~al., 2022]
Marc Lambert, Silvere Bonnabel, and Francis Bach.
\newblock The recursive variational gaussian approximation (r-vga).
\newblock \emph{Statistics and Computing}, 32\penalty0 (1):\penalty0 10, 2022.

[Li et~al., 2021]
Zhiyuan Li, Yi~Zhang, and Sanjeev Arora.
\newblock Why are convolutional nets more sample-efficient than fully-connected
\newblock In \emph{International Conference on Learning Representations}, 2021.

[Marchand, 1993]
Eric Marchand.
\newblock Estimation of a multivariate mean with constraints on the norm.
\newblock \emph{The Canadian Journal of Statistics}, 21\penalty0 (4):\penalty0
359--366, 1993.

[Nemirovsky and Yudin, 1983]
Arkadi Nemirovsky and David Yudin.
\newblock \emph{Problem Complexity and Method Efficiency in Optimization}.
\newblock John Wiley \& Sons, 1983.

[Ng, 2004]
Andrew~Y. Ng.
\newblock Feature selection, $l_1$ vs. $l_2$ regularization, and rotational
\newblock In \emph{Proceedings of the 21-st International Conference on Machine
Learning, Banff, Canada, July 4-8}. ACM New York, NY, 2004.

[Rigollet and Hütter, 2023]
Philippe Rigollet and Jan-Christian Hütter.
\newblock High-dimensional statistics, 2023.

[Shamir, 2018]
Ohad Shamir.
\newblock Distribution-specific hardness of learning neural networks.
\newblock \emph{Journal of Machine Learning Research}, 19\penalty0
(32):\penalty0 1--29, 2018.

[Vaskevicius et~al., 2019]
Tomas Vaskevicius, Varun Kanade, and Patrick Rebeschini.
\newblock Implicit regularization for optimal sparse recovery.
\newblock In \emph{Proceedings of Advances in Neural Information Processing
Systems (NeurIPS)}, pages 2968--2979, 2019.

[Warmuth and Amid, 2023]
Manfred~K. Warmuth and Ehsan Amid.
\newblock Open problem: Learning sparse linear concepts by priming the
\newblock In \emph{Proceedings of Thirty Sixth Conference on Learning Theory},
volume 195 of \emph{Proceedings of Machine Learning Research}, pages
5937--5942. PMLR, Jul 2023.

[Warmuth and Jagota, 1998]
Manfred~K. Warmuth and Arun Jagota.
\newblock Continuous and discrete time nonlinear gradient descent: relative
loss bounds and convergence.
\newblock In R.~Greiner E.~Boros, editor, \emph{Electronic Proceedings of Fifth
International Symposium on Artificial Intelligence and Mathematics}.
Electronic,http://rutcor.rutgers.edu/$\tilde{~}$amai, 1998.

[Warmuth and Vishwanathan, 2005]
Manfred~K. Warmuth and S.V.N. Vishwanathan.
\newblock Leaving the span.
\newblock In \emph{Proceedings of the 18th Annual Conference on Learning Theory
(COLT)}, pages 366--381, 2005.

[Warmuth et~al., 2021]
Manfred~K. Warmuth, Wojciech Kot{\l}owski, and Ehsan Amid.
\newblock A case where a spindly two-layer linear network decisively
outperforms any neural network with a fully connected input layer.
\newblock In \emph{32th International Conference on Algorithmic Learning Theory
(ALT)}, volume 132, pages 1--32. PMLR, 2021.

[Xiao et~al., 2017]
Han Xiao, Kashif Rasul, and Roland Vollgraf.
\newblock Fashion-mnist: a novel image dataset for benchmarking machine
learning algorithms, 2017.
\newblock URL \url{http://arxiv.org/abs/1708.07747}.
\newblock cite arxiv:1708.07747Comment: Dataset is freely available at
https://github.com/zalandoresearch/fashion-mnist Benchmark is available at

\end{thebibliography}

%\paragraph{Parameter tuning.} For simplicity we assume $\sqrt{d}$ is an integer, and $d \ge 4$.
%The algorithm runs for $T = 4\sqrt{d}$ trials, with learning rate chosen as $\eta=1/2$. However,
%the analysis would also work (with modification) for smaller values of learning rate. We
%set the initial weights to $\beta = 1/(2d)$ (having $\beta$ of order $d^{-1}$ seems necessary as for $\beta = \Theta(1)$ the algorithm essentially becomes GD). Finally note that
%we made an assumption that $\gamma \le \frac{1}{\sqrt{8d}}$ (for which we need $m$ to grow
%at least logarithmically with $d$, which is a very mild assumption).

\paragraph{The analysis for `zero signal' direction.} Since every weights evolves independently of the other weights, we can analyze each coordinate separately. We start with any coordinate
$i \ge 2$ ('zero signal' weights), for which the update becomes
\[
s_{t+1,i} = (1 - 2\eta\sqrt{(s_{t,i} + \zeta_i)^2 + 4\beta^2}) s_{t,i},
\]
with $s_{t,1} = - \zeta_1$. Now, w.l.o.g. assume $\zeta_i < 0$ (the analysis for $\zeta_i > 0$ is analogous). This means that $s_{1,i} > 0$, and $s_{t,i}$ is positive and
monotonically decreasing as long as $2\eta\sqrt{(s_{t,i} + \zeta_i)^2 + 4c^2} < 1$;
this condition is ensured by noticing that
\[
(s_{t,i} + \zeta_i)^2 + 4\beta^2 \le \zeta_i^2 + d^{-2} \le \gamma^2 + \frac{1}{4} \le
\frac{1}{2},
\]
so that using learning rate $\eta < \sqrt{2}/2$ will do the trick (remind that we use $\eta=1/4$).
Since we know that $s_{t,i}$ is positive and monotonically decreasing, we can bound:
\[
s_{t+1,i} \ge (1 - 2\eta\sqrt{(s_{T+1,i} + \zeta_i)^2 + 4\beta^2}) s_{t,i},
\]
so that
\[
s_{T+1,i} \ge (1 - 2\eta\sqrt{(s_{T+1,i} + \zeta_i)^2 + 4\beta^2})^T s_{1,i}
\ge (1 - 2\eta T \sqrt{(s_{T+1,i} + \zeta_i)^2 + 4\beta^2}) \zeta_i,
\]
where we used the Bernoulli inequality $(1+x)^n \ge 1 +xn$.
Returning to the original variable $w_{t,i} = s_{t,i} + \zeta_i$ gives
\[
w_{T+1,i} \ge 2 \eta T \zeta_i \sqrt{w_{T+1,i}^2 + 4\beta^2}.
\]
Since we also know that $w_{T+1,i} < 0$ (because $s_{t,i}=w_{t,i}-\zeta_i$ was
decreasing in $t$ with $s_{1,i}=-\zeta_i$ and $w_{1,i}=0$), we have
\[
w_{T+1,i}^2 \le 4 \eta^2 T^2 \zeta_i^2 (w_{T+1,i}^2 + 4\beta^2),
\]
which can be solved for $w_{T+1,i}^2$:
\[
w_{T+1,i}^2  \le \frac{16 \beta^2 \eta^2 T^2 \zeta_i^2}{1 - 4 \eta^2 T^2 \zeta_i^2}
= \frac{1}{d^2} \frac{4\eta^2 T^2 \gamma^2}{1 - 4\eta^2 T^2 \gamma^2}
\]
This expression is increasing in $\zeta_i^2$ so we can upper-bound it by
\[
w_{T+1,i}^2  \le \frac{1}{d^2} \frac{4\eta^2 T^2 \gamma^2}{1 - 4\eta^2 T^2 \gamma^2}
\]
Since $\eta=1/4$, $T=4\sqrt{d}$ and $\gamma \le \frac{1}{\sqrt{8d}}$, the denominator
is bounded from below
\[
1 - 4\eta^2 T^2 \gamma^2 \ge 1 - \frac{1}{2} \ge \frac{1}{2},
\]
so that, using $\gamma = \sigma \sqrt{\frac{\ln \frac{2d}{\delta}}{md}}$,
\[
w_{T+1,i}^2 \le \frac{1}{d^2} \frac{8d \cdot \sigma^2\ln \frac{2d}{\delta}}{md}
= \frac{8 \sigma^2 \ln \frac{2d}{\delta}}{m d^2}.
\]
Thus, the total error from `zero-signal' coordinates is
\begin{equation}
\sum_{i=2}^d w_{T+1,i}^2 \le \frac{8 \sigma^2 \ln \frac{2d}{\delta}}{m d}
\label{eq:error_bound_2_to_d}
\end{equation}
We also need to show that the same amount of error comes from the first coordinate.
%w_{t+1,i} = (1-\eta \sqrt{w_{t,i}^2+4c^2}) w_{t,i} - \eta \sqrt{w_{t,i}^2 + 4c^2} \zeta_i
%Now, w.l.o.g. assume $\zeta_i < 0$. Since $w_{1,i}=0$, $w_{2,i}$ becomes positive, and
%remains so. We can thus write
%|w_{t+1,i}| =  (1-\eta \sqrt{w_{t,i}^2+4c^2}) |w_{t,i}| + \eta \sqrt{w_{t,i}^2 + 4c^2} |\zeta_i|
%\le |w_{t,i}| + \eta \sqrt{w_{t,i}^2 + 4c^2} |\zeta_i|
%Now, we make the following assumption which will be verified later:
%$m \eta |\zeta_i| \le c/2$. This way, we first show that $|w_{t,i}| \le t \eta c |\zeta_i|$ for all $t=1,\ldots,m$. Indeed, $w_{1,i} = 0$, and by induction if $|w_{t,i}| \le 3 t \eta c |\zeta_i|$, we have
%|w_{t+1,i}| &\le 3 t \eta c |\zeta_i| + \eta \sqrt{( 3 t \eta c |\zeta_i|)^2 + 4c^2} |\zeta_i| \\
%&\le 3 t \eta c |\zeta_i| + \eta \sqrt{9 (m \eta |\zeta_i|)^2 + 4c^2} |\zeta_i| \\
%&\le 3 t \eta c |\zeta_i| + \eta \sqrt{(9/4 + 4)c^2} |\zeta_i|
%\le 3 (t+1) \eta c |\zeta_i|

\paragraph{Analysis for the `signal' coordinate.}
Using \eqref{eq:update_s}, we have for $s_{t,1} = 1 - \zeta_i - w_{t,1}$:
\[
s_{t+1,1} =  (1 - 2\eta \sqrt{(1-s_{t,1}-\zeta_i)^2+4\beta^2}) s_{t,1},
\]
with $s_{1,1} = 1 - \zeta_i$. As before, we note that $s_{t,i}$ is decreasing in $t$,
as long as $2\eta \sqrt{(1-s_{t,1}-\zeta_i)^2+4\beta^2} < 1$. However, this condition is satisfied
for our choice of of $\eta=1/4$, because
\begin{align*}
2\eta \sqrt{(1-s_{t,1}-\zeta_i)^2+4\beta^2} &\le
2\eta \sqrt{(1-\zeta_i)^2+4\beta^2} \le \frac{1}{2} \sqrt{1 + 2|\zeta_i| + \zeta_i^2 + d^{-2}}\\
&\le \frac{1}{2} \sqrt{1 + \frac{1}{\sqrt{2d}} + \frac{9}{8d^2}}
\stackrel{d \ge 4}{\le} \frac{1}{2} \sqrt{1 + \frac{1}{8} + \frac{9}{128}} < 1.
\end{align*}
Now, we need to carefully analyze the update. Initially $s_{t,1}$ decreases slowly, as
the square root term is essentially of order $1/d$. At some point, however,
$s_{t,1}$ falls below a certain constant (say, $s_{t,1} = 1/2$), and the square root
term is of order $O(1)$, and the convergence becomes exponential.

First, to simplify analysis we simply denote $s_{t,i}$ by $s_t$; moreover, define $r=1-\zeta_i$, so that the update becomes
\begin{equation}
s_{t+1}= (1- 2\eta \sqrt{(r-s_t)^2 + 4\beta^2} s_t, \qquad s_1 = r,
\label{eq:s_1_simplified}
\end{equation}

Already after the first iteration,
\[
s_2 = (1- 1/2 d^{-1})s_t = \frac{2d-1}{2d}r,
\]
so that
$(r-s_2)^2 = 1/(4d^2)$ becomes comparable with $4 \beta^2 = 1/d^2$ term. So we can drop
the $4\beta^2$ term from the square root in \eqref{eq:s_1_simplified} and upper bound
\begin{equation}
s_{t+1} \le (1-2\eta(r-s_t))s_t
\label{eq:s_1_upper_bound}
\end{equation}
To get some insight into this expression we solve the corresponding differential equation:
\[
\dot{s} = -2\eta(r-s)s,
\]
which give:
\[
\frac{s_t}{r-s_t} = C e^{-\eta r t} \quad \Longrightarrow \quad s_t = \frac{r}{1 + C e^{2\eta r t}}.
\]
Inspired by this we will bound $\frac{s_t}{r-s_t}$. From \eqref{eq:s_1_upper_bound} we get:
\[
r-s_{t+1} \ge r - (1-2\eta(r-s_t)s_t = (1 + 2\eta s_t) (r-s_t),
\]
so that:
\[
\frac{s_{t+1}}{r-s_{t+1}} \le \underbrace{\left(\frac{1-2\eta(r-s_t)}{1+2\eta s_t}\right)}_{=:A_t} \;
\left(\frac{s_t}{r-s_t}\right).
\]
We will now bound $A_t$ independent of $s_t$. To this end, note that
$A_t$ is maximized when $s_t = r$. Indeed,
\[
A_t = \frac{1+2\eta s_t - 2\eta r}{1 + 2\eta s_t} = 1 - \frac{2\eta r}{1 + 2\eta s_t}
\le 1 - \frac{2\eta r}{1 + 2\eta r} = \frac{1}{1+2\eta r}.
\]
This way, we get an upper bound:
\[
\frac{s_{T+1}}{r-s_{T+1}} \le (1+2\eta r)^{-(T-1)} \frac{s_2}{r-s_2}
= (2d-1) (1+2\eta r)^{-(T-1)},
\]
or by solving for $s_{T+1}$,
\[
s_{T+1} \le \frac{r}{1 + (2d-1)^{-1} (1+2\eta r)^{T-1}} \le r (2d-1) (1 + 2\eta r)^{-(T-1)}.
\]
The expression above is decreasing in $r$ for $T \ge 4$ (can be verified by computing the derivative), so we we will upper-bound it by lower-bounding $r$,
that is $r = 1 - \zeta_i \ge 1 - \frac{1}{\sqrt{8d}}$.
Taking $\eta = 1/4$ and using $1 - \frac{1}{\sqrt{8d}}
\stackrel{d \ge 4}{\ge} 1 - \frac{1}{4 \sqrt{2}}$,
we get $1 + 2\eta r >= \frac{3}{2} - \frac{1}{8 \sqrt{2}} > e^{1/3}$ (checked numerically).
\[
s_{T+1} \le (2d-1) e^{-(T-1)/3}
\le e^{1/3 + \ln 2} e^{-T/3 + \ln d} \le 3 e^{-T/3 + \ln d}.
\]
Using the fact that $T=4\sqrt{d}$, we get
\[
s_{T+1} \le 3 e^{-\frac{4}{3} \sqrt{d} + \ln d}
\]
To bound $(1-w_{T+1,1})^2$ we use
\begin{equation}
(1-w_{T+1,1})^2 = (s_{T+1,1}+\zeta_1)^2 \le
2s_{T+1,1}^2 + 2 \zeta_1^2
\le 9 e^{-\frac{8}{3} \sqrt{d} + 2 \ln d} + \frac{2 \sigma^2 \ln \frac{2d}{\delta}}{md}.
\label{eq:error_bound_1}
\end{equation}

\paragraph{Bound the error of Approximated EGU${}^{\pm}$ algorithm}

The final error of the algorithm is obtained by summing \eqref{eq:error_bound_2_to_d} and \eqref{eq:error_bound_1}:
\begin{align*}
\|\w_{T+1} - \e_1\|^2
&\le 9 e^{-\frac{8}{3} \sqrt{d} + 2 \ln d} + \frac{2 \sigma^2 \ln \frac{2d}{\delta}}{md} + \frac{8 \sigma^2 \ln \frac{2d}{\delta}}{m d} \\
&= \frac{10 \sigma^2 \ln \frac{2d}{\delta}}{m d}
+ 9 e^{-\frac{8}{3} \sqrt{d} + 2 \ln d}.
\end{align*}

\vspace{-.9cm}
\hfill $\BlackBox$

\section{Upper bound for the spindly network} \label{app:spindly}
%\Red{Need a statement of the upper bound Theorem and
%mention the alternate R.I.P. proofs again.}

In this section, we upper-bound the error of the spindly network defined in
Figure \ref{f:spindly}, showing essentially the same bound (up to constants) as for the Approximated
EGU$^{\pm}$. We note that essentially the same algorithm has been analyzed by
Vaskevicius et~al., 2019, giving the same bound $O(\frac{\sigma^2\log d}{md})$ under the
restricted isometry property (RIP) assumption.

\begin{theorem}
Assume $d \ge 4$ is such that $\sqrt{d}$ is an integer. Consider the spindly network given in Figure \ref{f:spindly} trained with gradient descent, with weights initialized as $\u = \sqrt{2/d} \bm{1}$ and $\v = \bm{0}$, and the learning rate set to $\eta=1/4$. Let $m \ge 8 \sigma^2 \ln \frac{2 d}{\delta} = \Omega(\sigma^2 \log (d/\delta))$ . With probability at least $1-\delta$, the algorithm run for $T=4\sqrt{d}$ steps achieves error bounded by:
\[
e(\w_{T+1}) \le
\frac{33 \sigma^2 \ln \frac{2d}{\delta}}{4 m d}
+ 16 e^{-\frac{8}{3} \sqrt{d} + 2 \ln d}
= O\left(\frac{\sigma^2\log d}{md} + e^{-\frac83\sqrt{d}}\right)
\]
\end{theorem}

\begin{proof} The spindly network predicts with $\w_t$ given by $\w_t = \u_t \odot \v_t$, and both
$\u_t$ and $\v_t$ are updated with the gradient descent algorithm:
\[
\u_{t+1} = \u_t - \eta \nabla_{\u_t} L(\w_t), \qquad
\v_{t+1} = \v_t - \eta \nabla_{\v_t} L(\w_t)
\]
Using \eqref{eq:gradient_simplified} and the chain rule, $\nabla_{\u_t} \emploss(\w_t) = 2(\w_t - \e_1 - \bm{\zeta}) \odot \v_t$ and $\nabla_{\v_t} \emploss(\w_t) = 2(\w_t - \e_1 - \bm{\zeta}) \odot \u_t$, so that
\begin{equation}
\u_{t+1} = \u_t - 2\eta (\w_t - \e_1 - \bm{\zeta}) \odot \v_t, \qquad
\v_{t+1} = \v_t - 2\eta (\w_t - \e_1 - \bm{\zeta}) \odot \u_t
\label{eq:uv_update}
\end{equation}
Let us introduce two vectors,
$\v^+_t$ and $\v^-_t$, given by:
\[
\v^+_t = \frac{1}{4} (\u_t + \v_t)^2, \qquad \v^-_t = \frac{1}{4} (\u_t - \v_t)^2,
\]
(square applied coordinatewise) and note that
\[
\w_t = \u_t \cdot \v_t = \v^+_t - \v^-_t.
\]
Subtracting and adding equations in \eqref{eq:uv_update}, followed by squaring both sides
\begin{align}
\v^+_{t+1} &= (1 - 2 \eta (\w_t - \e_1 - \bm{\zeta}))^2 \odot \v^+_t \nonumber \\
\v^-_{t+1} &= (1 + 2 \eta (\w_t - \e_1 - \bm{\zeta}))^2 \odot \v^-_t  \label{eq:v_pm_spindly}
\end{align}
We initialize the algorithm as follows:
\[
\u_1 = \sqrt{\frac{2}{d}}, \qquad \v_1 = \bm{0},
\]
so that $\v_1^+ = \v_1^- = \frac{1}{2d} \bm{1}$ and $\w_1 = \bm{0}$, similarly
as in the Approximated EGU${}^{\pm}$. We set $\eta = \frac{1}{8}$, $T=4\sqrt{d}$,
and use the same assumption as in the previous section, that is $m \ge 8 \sigma^2 \ln \frac{2d}{\delta}$, which implies that with probability
at least $1-\delta$, $|\zeta_i| \le \gamma$ for all $i=1,\ldots,d$, where
\begin{equation}
\gamma = \sigma \sqrt{\frac{\ln \frac{2d}{\delta}}{md}} \le \frac{1}{\sqrt{8d}}.
\label{eq:upper_bound_on_gamma}
\end{equation}
As before, we denote the high probability event above as $E$, and we condition everything what follows on the fact that
$E$ happened. We will also assume that $d \ge 9$.

\paragraph{The analysis for `zero signal' direction.} As before, we can analyze each coordinate separately. We start with any coordinate
$i \ge 2$, for which the update \eqref{eq:v_pm_spindly} becomes
\begin{equation}
v^{\pm}_{t+1,i} = (1 \mp 2\eta(w_{t,i} - \zeta_i))^2 v^{\pm}_{t,i}
\label{eq:v_pm_update_zero_signal}
\end{equation}
W.l.o.g. assume $\zeta_i > 0$ (the analysis for $\zeta_i < 0$ is analogous). We will first
prove by induction on $t$ that $0 \le w_{t,i} \le \zeta_i$ for all $t=1,\ldots,T+1$. Since $w_{1,i}=0$ and $\zeta_i > 0$, it clearly holds for $t=1$. Now, assume that it holds for iterations $1,\ldots,t$ and we prove that it also holds for $t+1$. We have:
\begin{align}
w_{t+1,i} &= v_{t+1,i}^+ - v_{t+1,i}^- \nonumber \\
&= (1 - 2\eta(w_{t,i} - \zeta_i))^2 v_{t,i}^+ -
(1 + 2\eta(w_{t,i} - \zeta_i))^2 v_{t,i}^- \nonumber \\
&= (1+4 \eta^2(w_{t,i} - \zeta_i)^2) w_{t,i}
- 4 \eta (w_{t,i} - \zeta_i) (v^+_{t,i} + v^-_{t,i}) \ge 0, \label{eq:deriv_spindly_induction}
\end{align}
because  $w_{t,i} - \zeta_i < 0$ and $w_{t,i} \ge 0$ from the induction assumption, while
$v^+_{t,i}, v^-_{t,i} \ge 0$ from their definitions.

Now, since $-\zeta_i \le w_{q,i} - \zeta_i \le 0$ for $q=1,\ldots,t$ from the inductive assumption,
$\eta = \frac{1}{8}$ and $|\zeta_i| \le \gamma \le \frac{1}{\sqrt{8d}}$ from \eqref{eq:upper_bound_on_gamma}, we have $0 \ge 2 \eta (w_{q,i} - \zeta_i) \ge - \frac{1}{4 \sqrt{8d}}$,
so that
\[
(1 + 2\eta(w_{q,i} - \zeta_i))^2 < 1, \qquad
(1 - 2\eta(w_{q,i} - \zeta_i))^2 > 1, \qquad q=1,\ldots,t.
\]
Thus, we see from \eqref{eq:v_pm_update_zero_signal} that $v^+_{q,i}$ is monotonically
increasing in $q$, and $v^-_{q,i}$ is monotonically decreasing in $q$. This means that
$w_{q,i} = v^+_{q,i} - v^-_{q,i}$ is monotonically increasing in $q$. Therefore,
\[
(1 + 2\eta(w_{q,i} - \zeta_i))^2 \ge (1 + 2\eta(w_{1,i} - \zeta_i))^2
= (1 - 2 \eta \zeta_i)^2 \qquad
\text{for all~} q=1,\ldots,t,
\]
and, similarly,
\[
(1 - 2\eta(w_{q,i} - \zeta_i))^2 \le (1 + 2\eta(w_{1,i} - \zeta_i))^2 = (1 + 2\eta \zeta_i)^2 \qquad
\text{for all~} q=1,\ldots,t.
\]
This let us upper-bound $v^+_{t+1,i}$ as
\[
v^+_{t+1,i} = \prod_{q=1}^t (1 - 2 \eta(w_{q,i} - \zeta_i))^2 v^+_{1,i}
\le (1 + 2 \eta \zeta_i)^{2t} v^+_{1,i} =
(1 + 2 \eta \zeta_i)^{2t} \frac{1}{2d},
\]
and similarly lower-bound $v^-_{t+1,i}$ as
\[
v^-_{t+1,i} = \prod_{q=1}^t (1 + 2 \eta(w_{q,i} - \zeta_i))^2 v^-_{1,i}
\ge (1 - 2 \eta \zeta_i)^{2t} v^+_{1,i} =
(1 - 2 \eta \zeta_i)^{2t} \frac{1}{2d}.
\]
We have
\[
(1+2 \eta \zeta_i)^{2t} = e^{2t \ln(1+ 2\eta \zeta_i)} \le e^{4 t \eta \zeta_i},
\]
where we used $\ln(1+x) \le x$. Moreover, by Bernoulli's inequality,
\[
(1 - 2 \eta \zeta_i)^{2t} \ge 1 - 4t \eta \zeta_i.
\]
This gives:
\[
w_{t+1,i} =
v^+_{t+1,i} - v^-_{t+1,i}
=  \frac{1}{2d} \left( (1 + 2 \eta \zeta_i)^{2t} - (1 - 2 \eta \zeta_i)^{2t}  \right)
\le \frac{1}{2d} (e^{4 t \eta \zeta_i} - 1 + 4 t \eta \zeta_i).
\]
Now, note that
\[
4 t \eta \zeta_i \le \frac{T}{2} \gamma \le \frac{4 \sqrt{d}}{2 \sqrt{8 d}} \le \frac{1}{\sqrt{2}},
\]
Using the convexity of $e^x$, we have for $x \in [0,a]$ that
\[
e^x = e^{\left(1-\frac{x}{a}\right) \cdot 0 + \frac{x}{a} \cdot a} \le \left(1-\frac{x}{a}\right) e^0 + \frac{x}{a} e^a
= 1 + x \frac{e^a -1}{a}.
\]
Taking $x = 4 t \eta \zeta_i$ and $a = \frac{1}{\sqrt{2}}$, we have
\[
e^{4 t \eta \zeta_i}
\le 1 + 4 t \eta \zeta_i \sqrt{2} (e^{1/\sqrt{2}} - 1) \le 1 + 6 t \eta \zeta_i.
\]
This allows us to bound
\begin{equation}
w_{t+1,i} \le \frac{1}{2d} (e^{4 t \eta \zeta_i} - 1 + 4 t \eta \zeta_i)
\le \frac{1}{2d} 10 t \eta \zeta_i \le \frac{5 T \eta \zeta_i}{d}
= \frac{5}{2  \sqrt{d}} \zeta_i \le \frac{5}{6} \zeta_i \le \zeta_i,
\label{eq:zero_signal_useful_inequality}
\end{equation}
where we used $d \ge 9$. This finishes the inductive proof that $0 \le w_{t,i} \le \zeta_i$
for all $t = 1,\ldots,T+1$. However, applying \eqref{eq:zero_signal_useful_inequality} to $t=T$ gives
\[
w_{T+1,i} \le  \frac{5}{2  \sqrt{d}} \zeta_i,
\]
so that using \eqref{eq:upper_bound_on_gamma}
\[
w^2_{T+1,i} \le  \frac{25}{4 d} \zeta^2_i
\le \frac{25}{4 d} \gamma^2
=\frac{25 \sigma^2 \ln \frac{2d}{\delta}}{4 md^2}.
\]
Thus, the total error from `zero-signal' coordinates is
\begin{equation}
\sum_{i=2}^d w_{T+1,i}^2 \le \frac{25 \sigma^2 \ln \frac{2d}{\delta}}{4 m d}.
\label{eq:spindly_error_bound_for_i_ge_2}
\end{equation}

\paragraph{Analysis for the `signal' coordinate.}
For $i=1$, the update becomes:
\[
v^{\pm}_{t+1,1} = (1 \mp 2\eta(w_{t,i} - 1 - \zeta_1))^2 v^{\pm}_{t,1}.
\]
First, we will show by induction that $0 \le w_{t,1} \le 1 + \zeta_1$ for
all $t=1,\ldots,T+1$.
This is clearly true for $t=1$ as $w_{t,1} = 0$ and $|\zeta_1| \le \frac{1}{\sqrt{8d}}$. Assume now this is true for $q=1,\ldots,t$, and we prove it is also true for $t+1$.
We have
\begin{align}
w_{t+1,1} &= v_{t+1,1}^+ - v_{t+1,1}^- \nonumber \\
&= (1 - 2\eta(w_{t,1} - 1 - \zeta_1))^2 v_{t,1}^+ -
(1 + 2\eta(w_{t,1} - 1 - \zeta_1))^2 v_{t,1}^- \nonumber \\
&= (1+4 \eta^2(w_{t,1} - 1 - \zeta_1)^2) w_{t,1}
- 4 \eta (w_{t,1} - 1 - \zeta_1) (v^+_{t,1} + v^-_{t,1}) \nonumber \\
&= (1+4 \eta^2(w_{t,1} - 1 - \zeta_1)^2) w_{t,1}
- 4 \eta (w_{t,1} - 1 - \zeta_1) (w_{t,1} + 2 v^-_{t,1}) \nonumber \\
&= (1 - 2\eta(w_{t,1} - 1 - \zeta_1))^2 w_{t,1} + 8 \eta (1 + \zeta_i - w_{t,1}) v^-_{t,1}.
\label{eq:expression_for_w_t_plus_1_signal}
\end{align}
It follows from the induction assumption that both terms in the last line are nonnegative,
thus $w_{t+1,1} \ge 0$. To show that $w_{t+1,1} \le 1 + \zeta_1$,
note that by inductive assumption
$1 + 2\eta(w_{t,i} - 1 - \zeta_1)) \ge 1 - 2\eta(1- \zeta_1)
\ge 1 - 2\eta(1 + \gamma) = \frac{3}{4} - \frac{1}{4} \gamma \ge
\frac{3}{4} - \frac{1}{4\sqrt{8}{d}} > 0$, and also
$1 + 2\eta(w_{t,i} - 1 - \zeta_1)) \le 1$. This means that $v^{-}_{q,1}$ is nonincreasing
in $q$, and thus $v^{-}_{t,1} \le v^{-}_{1,1} = \frac{1}{2d} \le \frac{1}{18}$ (due to $d \ge 9$). Plugging this into \eqref{eq:expression_for_w_t_plus_1_signal}, we can bound
\begin{align}
w_{t+1,1} &\le