Variance Reduced Distributed Non-Convex Optimization Using Matrix Stepsizes

## Abstract

Abstract Matrix-stepsized gradient descent algorithms have been shown to have superior performance in non-convex optimization problems compared to their scalar counterparts.
The det-CGD algorithm, as introduced by Li et al. ( 2024b ) , leverages matrix stepsizes to perform compressed gradient descent for non-convex objectives and matrix-smooth problems in a federated manner.
The authors establish the algorithm’s convergence to a neighborhood of a weighted stationarity point under a convex condition for the symmetric and positive-definite matrix stepsize.
In this paper, we propose two variance-reduced versions of the det-CGD algorithm, incorporating MARINA and DASHA methods.
Notably, we establish theoretically and empirically, that det-MARINA and det-DASHA outperform MARINA , DASHA and the distributed det-CGD algorithms in terms of iteration and communication complexities.

## 1 Introduction

We focus on optimizing the finite sum non-convex objective

$$ $\min_{x\in\mathbb{R}^{d}}\left\{f(x):=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x)\right\}.$ (1) $$

In this context, each function $f_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is differentiable and bounded from below.
This type of objective function finds extensive application in various practical machine learning algorithms, which increase not only in terms of the data size but also in the model size and overall complexity as well.
As a result, most neural network architectures result in highly non-convex empirical losses, which need to be minimized.
In addition, it becomes computationally infeasible to train these models on one device, often excessively large, and one needs to redistribute them amongst different devices/clients.
This redistribution results in a high communication overhead, which often becomes the bottleneck in this framework.

In other words, we have the following setting.
The data is partitioned into $n$ clients, where the $i$-th client has access to the component function $f_{i}$ and its derivatives.
The clients are connected to each other through a central device, called the server.
In this work, we are going to study iterative gradient descent-based algorithms that operate as follows. The clients compute the local gradients in parallel. Then they compress these gradients to reduce the communication cost and send them to the server in parallel.
The server then aggregates these vectors and broadcasts the iterate update back to the clients. This meta-algorithm is called federated learning.
We refer the readers to Konečný et al. (); McMahan et al. (); Kairouz et al. () for a more thorough introduction to federated learning.

### 1.1 Contributions

In this paper, we introduce two novel federated learning algorithms named det-MARINA and det-DASHA.
These algorithms extend a recent method called det-CGD (Li et al., ), which aims to solve problem (1) using matrix stepsized gradient descent.
Under the matrix smoothness assumption proposed by Safaryan et al. (), the authors demonstrate that the matrix stepsized version of the Distributed Compressed Gradient
Descent (Khirirat et al., ) algorithm enhances communication complexity compared to its scalar counterpart.
However, in their analysis, Li et al. () show stationarity only within a certain neighborhood due to stochastic compressors.
Our algorithm addresses this issue by incorporating previously known variance reduction schemes, namely,
MARINA (Gorbunov et al., ) and DASHA (Tyurin & Richtárik, ).
We establish theoretically and empirically, that both algorithms outperform their scalar alternatives, as well as the distributed det-CGD algorithms.
In addition, we describe specific matrix stepsize choices, for which our algorithms beat MARINA, DASHA and distributed det-CGD both in theory and in practice.

## 2 Background and motivation

For a given $\varepsilon>0$, finding an approximately global optimum, that is $x_{\varepsilon}$ such that $f(x_{\varepsilon})-\min_{x}f(x)<\varepsilon$, is known to be NP-hard (Jain et al., ; Danilova et al., ).
However, gradient descent based methods are still useful in this case.
When these methods are applied to non-convex objectives, they treat the function $f$ as locally convex and aim to converge to a local minimum.
Despite this simplification, such methods have gained popularity in practice due to their superior performance compared to other approaches for non-convex optimization, such as convex relaxation-based methods (Tibshirani, ; Cai et al., ).

### 2.1 Stochastic gradient descent

Arguably, one of the most prominent meta-methods for tackling non-convex optimization problems is stochastic gradient descent ( SGD). The formulation of SGD is presented as the following iterative algorithm:
$x^{k+1}=x^{k}-\gamma g^{k}$.
Here, $g^{k}\in\mathbb{R}^{d}$ serves as a stochastic estimator of the gradient $\nabla f(x^{k})$.
SGD essentially mimics the classical gradient descent algorithm, and recovers it when $g^{k}=\nabla f(x^{k})$.
In this scenario, the method approximates the objective function $f$ using a linear function and takes a step of size $\gamma$ in the direction that maximally reduces this approximation.
When the stepsize is sufficiently small, and the function $f$ is suitably smooth, it can be demonstrated that the function value decreases, as discussed in (Bubeck et al., ; Gower et al., ).

However, computing the full gradient can often be computationally expensive.
In such cases, stochastic approximations of the gradient come into play.
Stochastic estimators of the gradient can be employed for various purposes, leading to the development of different methods.
These include stochastic batch gradient descent (Nemirovski et al., ; Johnson & Zhang, ; Defazio et al., ), randomized coordinate descent (Nesterov, ; Wright, ), and compressed gradient descent
(Alistarh et al., ; Khirirat et al., ; Mishchenko et al., ).
The latter, compressed gradient descent, holds particular relevance to this paper, and we will delve into a more detailed discussion of it in subsequent sections.

### 2.2 Second order methods

The stochastic gradient descent is considered as a first-order method as it uses only the first order derivative information.
Although being immensely popular, the first order methods are not always optimal.
Not surprisingly, using higher order derivatives in deciding update direction can yield to faster algorithms.
A simple instance of such algorithms is the Newton Star algorithm (Islamov et al., ):

$$ $x^{k+1}=x^{k}-\left(\nabla^{2}f(x^{\star})\right)^{-1}\nabla f(x^{k}),$ ( NS ) $$

where $x^{\star}$ is the minimum point of the objective function.
The authors establish that under specific conditions, the algorithm’s convergence to the unique solution $x^{\star}$ in the convex scenario occurs at a local quadratic rate.
Nonetheless, its practicality is limited since we do not have prior knowledge of the Hessian matrix at the optimal point.
Despite being proposed recently, the Newton-Star algorithm gives a deeper insight on the generic Newton method (Gragg & Tapia, ; Miel, ; Yamamoto, ):

$$ $x^{k+1}=x^{k}-\gamma\left(\nabla^{2}f(x^{k})\right)^{-1}\nabla f(x^{k}).$ ( NM ) $$

Here, the unknown Hessian of the Newton-Star algorithm, is estimated progressively along the iterations.
The latter causes elevated computational costs, as the inverting a large square matrix is expensive.
As an alternative, quasi-Newton methods replace the inverse of the Hessian at the iterate with a computationally cheaper estimate (Broyden, ; Dennis & Moré, ; Al-Baali & Khalfan, ; Al-Baali et al., ).

### 2.3 Fixed matrix stepsizes

The det-CGD algorithm falls into this framework of the second order methods as well. Proposed by Li et al. ()(^1^11In the original paper, the algorithm is referred to as det-CGD, as there is a variant of the same algorithm named det-CGD2. Since we are going to use only the first one and our framework is applicable to both, we will remove the number in the end for the sake of brevity.), the algorithm suggests using a uniform “upper bound” on the inverse Hessian matrix.
Assuming matrix smoothness of the objective (Safaryan et al., ), they replace the scalar stepsize with a positive definite matrix ${\bm{D}}$.
The algorithm is given as follows:

$$ $x^{k+1}=x^{k}-{\bm{D}}{\bm{S}}^{k}\nabla f(x^{k}).$ ( det-CGD ) $$

##### Matrix 𝑫 𝑫 {\bm{D}} bold_italic_D .

Here, ${\bm{D}}$ plays the role of the stepsize.
Essentially, it uniformly lower bounds the inverse Hessian.
The standard SGD is a particular case of this method, as the scalar stepsize $\gamma$ can be seen as a matrix
$\gamma{\bm{I}}_{d}$, where ${\bm{I}}_{d}$ is the $d$-dimensional identity matrix.
An advantage of using a matrix stepsize is more evident if we take the perspective of the second order methods.
Indeed, the scalar stepsize $\gamma{\bm{I}}_{d}$ uniformly estimates the largest eigenvalue of the Hessian matrix, while ${\bm{D}}$ can capture the Hessian more accurately.
The authors show both theoretical and empirical improvement that comes with matrix stepsizes.

##### Matrix 𝑺 k superscript 𝑺 𝑘 {\bm{S}}^{k} bold_italic_S start_POSTSUPERSCRIPT italic_k end_POSTSUPERSCRIPT .

We assume that ${\bm{S}}^{k}$ is a positive semi-definite, stochastic sketch matrix.
Furthermore, it is unbiased: $\mathbb{E}[{\bm{S}}^{k}]={\bm{I}}_{d}$.
We notice that det-CGD can be seen as a matrix stepsize instance of SGD, with $g^{k}={\bm{S}}^{k}\nabla f(x^{k})$.
The sketch matrix can be seen as a linear compressing operator, hence the name of the algorithm: Compressed Gradient Descent ( CGD) (Alistarh et al., ; Khirirat et al., ).
A commonly used example of such a compressor is the Rand-$\tau$ compressor.
This compressor randomly selects $\tau$ entries from its input and scales them using a scalar multiplier to ensure an unbiased estimation.
By adopting this approach, instead of using all $d$ coordinates of the gradient, only a subset of size $\tau$ is communicated.
Formally, Rand-$\tau$ is defined as follows:

$$ ${\bm{S}}=\frac{d}{\tau}\sum_{j=1}^{\tau}e_{i_{j}}e_{i_{j}}^{\top}.$ (2) $$

Here, $e_{i_{j}}$ denotes the $i_{j}$-th standard basis vector in $\mathbb{R}^{d}$.
For a more comprehensive understanding of compression techniques, we refer to Safaryan et al. ().

### 2.4 The neighborhood of the distributed det-CGD1

The distributed version of det-CGD follows the standard federated learning paradigm (McMahan et al., ).
The pseudocode of the method, as well as the convergence result of Li et al. (), can be found in Appendix F.
Informally, their convergence result can be written as

$$ $\min_{k=1,\ldots,K}\mathbb{E}\left[\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{ 2}\right]\leq\mathcal{O}\left(\frac{(1+\alpha)^{K}}{K}\right)+\mathcal{O}\left (\alpha\right),$ $$

where $\alpha>0$ is a constant that can be controlled.
The crucial insight from this result is that the error bound does not diminish as the number of iterations increases.
In fact, by controlling $\alpha$ and considering a large $K$, it is impossible to make the second term smaller than $\varepsilon$. This implies that the algorithm converges to a certain neighborhood surrounding the (local) optimum.
This phenomenon is a common occurrence in SGD and is primarily attributable to the variance associated with the stochastic gradient estimator. In the case of det-CGD the stochasticity comes from the sketch ${\bm{S}}^{k}$.

### 2.5 Variance reduction

To eliminate this neighborhood, various techniques for reducing variance are employed. One of the simplest techniques applicable to CGD is gradient shifting.
By replacing ${\bm{S}}^{k}\nabla f(x^{k})$ with ${\bm{S}}^{k}(\nabla f(x^{k})-\nabla f(x^{\star}))+\nabla f(x^{\star}),$ the neighborhood effect is removed from the general CGD.
This algorithm is an instance of a more commonly known method called ${\color[rgb]{0.68,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0.68,0,0}
\pgfsys@color@cmyk@stroke{0}{0.87}{0.68}{0.32}\pgfsys@color@cmyk@fill{0}{0.87}
{0.68}{0.32}\text{SGD}}_{\star}$ (Gower et al., ).
However, since the exact optimum $x^{\star}$ is typically unknown, this technique encounters similar challenges as the Newton-Star algorithm mentioned earlier.
Fortunately, akin to quasi-Newton methods, one can employ methods that iteratively learn the optimal shift (Shulgin & Richtárik, ).
A line of research focuses on variance reduction for CGD based algorithms on this insight.

To eliminate the neighborhood in the distributed version of CGD, denoted as det-CGD1, we apply a technique called MARINA (Gorbunov et al., ).
MARINA cleverly combines the general shifting (Shulgin & Richtárik, ) technique with loopless variance reduction techniques (Qian et al., ).
This approach introduces an alternative gradient estimator specifically designed for the federated learning setting.
Thanks to its structure, it allows to establish an upper bound on the stationarity error that diminishes significantly with a large number of iterations.
In this paper, we construct the analog of the this algorithm called det-MARINA, using matrix stepsizes and sketch gradient compressors.
For this new method, we prove a convergence guarantee similar to the results of Li et al. () without a neighborhood term.

Furthermore, we also propose det-DASHA, which is the extension of DASHA in the matrix stepsize setting.
The latter was proposed by Tyurin & Richtárik () and it combines MARINA with momentum variance reduction techniques (Cutkosky & Orabona, ). DASHA offers better practicality compared to MARINA, as it always sends compressed gradients and does not need to synchronize among all the nodes.

### 2.6 Organization of the paper

The rest of the paper is organized as follows. Section 3 discusses the general mathematical framework.
Section 4 and Section 5 present the det-MARINA and det-DASHA algorithms, respectively.
We show the superior theoretical performance of our algorithms compared to the relevant existing algorithms, that is MARINA, DASHA and det-CGD in Section 6.
The experimental results validating our theoretical findings are presented in Section 7, with additional details and setups available in the Appendix.
We conclude the paper by outlining several directions of future work in Section 8.

## 3 Mathematical framework

In this section we present the assumptions that we further require in the analysis.

###### Assumption 1 .

(Lower Bound)
There exists $f^{\star}\in\mathbb{R}$ such that, $f(x)\geq f^{\star}$ for all $x\in\mathbb{R}^{d}$.

This is a standard assumption in optimization, as otherwise the problem of minimizing the objective would not be correct mathematically.
We then need a matrix version of Lipschitz continuity for the gradient.

###### Definition 1 .

(${\bm{L}}$-Lipschitz Gradient)
Assume that $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a continuously differentiable function and matrix ${\bm{L}}\in\mathbb{S}^{d}_{++}$. We say the gradient of $f$ is ${\bm{L}}$-Lipschitz if for all $x,y\in\mathbb{R}^{d}$

$$ $\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}\leq\left\|x-y\right\|_{ {\bm{L}}}.$ (3) $$

In the following, we will assume that (3) is satisfied for component functions
$f_{i}$.

###### Assumption 2 .

Each function $f_{i}$ is ${\bm{L}}_{i}$-gradient Lipschitz, while $f$ is ${\bm{L}}$-gradient Lipschitz.

In fact, the second half of the assumption is a consequence of the first one.
Below, we formalize this claim.

###### Proposition 1 .

If $f_{i}$ is ${\bm{L}}_{i}$-gradient Lipschitz for every $i=1,\ldots,n$, then function $f$ has ${\bm{L}}$-Lipschitz gradient with ${\bm{L}}\in\mathbb{S}^{d}_{++}$ satisfying

$$ $\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}^{-1}\right)\cdot\lambda_ {\max}\left({\bm{L}}_{i}\right)\cdot\lambda_{\max}\left({\bm{L}}_{i}{\bm{L}}^{ -1}\right)=1.$ $$

###### Remark 1 .

In the scalar case, where ${\bm{L}}=L{\bm{I}}_{d}$, ${\bm{L}}_{i}=L_{i}{\bm{I}}_{d}$, the relation becomes
$L^{2}=\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}.$
This corresponds to the statement in Assumption 1.2 in (Gorbunov et al., ).

Nevertheless, the matrix ${\bm{L}}$ found according to Proposition 1 is only an estimate.
In principle, there might exist a better ${\bm{L}}_{f}\preceq{\bm{L}}$ such that $f$ has
${\bm{L}}_{f}$-Lipschitz gradient.

More generally, this condition can be interpreted as follows.
The gradient of $f$ naturally belongs to the dual space of $\mathbb{R}^{d}$,
as it is defined as a linear functional on $\mathbb{R}^{d}$.
In the scalar case, $\ell_{2}$-norm is self-dual, thus (3) reduces to the standard Lipschitz continuity of the gradient. However, with the matrix smoothness assumption, we are using the ${\bm{L}}$-norm for the iterates, which naturally induces the ${\bm{L}}^{-1}$-matrix norm for the gradients in the dual space.
This insight, which is originally presented by Nemirovski & Yudin (), plays a key role in our analysis.

See Appendix C for a more thorough discussion on the properties of 2, as well as its connection to matrix smoothness (Safaryan et al., ).

###### Assumption 1 .

###### Definition 1 .

###### Assumption 2 .

###### Proposition 1 .

###### Remark 1 .

## 4 MARINA-based variance reduction

In this section, we present our algorithm det-MARINA with its convergence result.
We construct a sequence of vectors $g^{k}$ which are stochastic estimators of $\nabla f(x^{k})$.
At each iteration, the server samples a Bernoulli random variable (coin flip) $c_{k}$ and broadcasts it in parallel to the clients, along with the current gradient estimate $g^{k}$.
Each client, then, does a det-CGD-type update with the stepsize ${\bm{D}}$ and a gradient estimate $g^{k}$.
The next gradient estimate $g^{k+1}$ is then computed. With a low probability, that is when $c_{k}=1$, we take the $g^{k+1}$ to be the full gradient $\nabla f(x^{k+1})$. Otherwise, we update it using the compressed gradient differences at each client.
See Algorithm 1 for the pseudocode of det-MARINA.

Figure: Algorithm 1 det-MARINA

### 4.1 Convergence guarantees

In the following theorem, we formulate one of the main results of this paper, which guarantees the convergence of Algorithm 1 under the above-mentioned assumptions.

###### Theorem 1 .

Assume that Assumptions 1 and 2 hold, and the following condition on stepsize matrix ${\bm{D}}\in\mathbb{S}^{d}_{++}$ holds,

$$ ${\bm{D}}^{-1}\succeq\left(\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})}{np}+1 \right){\bm{L}},$ (4) $$

where
$R({\bm{D}},\mathcal{S}):=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}
_{i}\right)\lambda_{\max}\left({\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-
\frac{1}{2}}\right)$ $\times\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{
i}^{k}\right]-{\bm{D}}\right).$
Then, after $K$ iterations of det-MARINA, we have

$$ ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{ \det({\bm{D}})^{1/d}}}\right]\leq\frac{2\left(f(x^{0})-f^{\star}\right)}{\det( {\bm{D}})^{1/d}\cdot K}.$ (5) $$

Here, $\tilde{x}^{K}$ is chosen uniformly randomly from the first $K$ iterates of the algorithm.

The criterion $\left\|\cdot\right\|^{2}_{{\bm{D}}/\det({\bm{D}})^{1/d}}$ is the same as that used in Li et al. (), known as determinant normalization.
The weight matrix of the matrix norm has determinant $1$ after normalization, which makes it comparable to the standard Euclidean norm.

###### Remark 2 .

We notice that the right-hand side of the algorithm vanishes with the number of iterations, thus solving the neighborhood issue of the distributed det-CGD.
Therefore, det-MARINA is indeed the variance reduced version of det-CGD in the distributed setting and has better convergence guarantees.

###### Remark 3 .

Theorem 1 implies the following iteration complexity for the algorithm. In order to get an $\varepsilon^{2}$ stationarity error(^2^22We say a (possibly random) vector $x\in\mathbb{R}^{d}$ is an $\varepsilon$-stationary point of a possibly non-convex function $f:\mathbb{R}^{d}\mapsto\mathbb{R}$, if ${\mathbb{E}}\left[\left\|\nabla f(x)\right\|^{2}\right]\leq\varepsilon^{2}$. The expectation is over the randomness of the algorithm), the algorithm requires $K$ iterations, with

$$ $K\geq\frac{2(f(x^{0})-f^{\star})}{\det({\bm{D}})^{1/d}\cdot\varepsilon^{2}}.$ $$

###### Remark 4 .

In the case where no compression is applied, that is we have ${\bm{S}}_{i}^{k}={\bm{I}}_{d}$, condition (4) reduces to ${\bm{D}}\preceq{\bm{L}}^{-1}$.
The latter is due to ${\mathbb{E}}\left[{\bm{S}}^{k}_{i}{\bm{D}}{\bm{S}}^{k}_{i}\right]={\bm{D}}$, which results in $R({\bm{D}},\mathcal{S})=0$. This is expected, since in the deterministic case det-MARINA reduces to GD with matrix stepsize.

The convergence condition and rate of matrix stepsize GD can be found in (Li et al., ).
Below we do a sanity check to verify that the convergence condition for scalar MARINA can be obtained.

###### Remark 5 .

Let us consider the scalar case.
That is ${\bm{L}}_{i}=L_{i}{\bm{I}}_{d},{\bm{L}}=L{\bm{I}}_{d},{\bm{D}}=\gamma{\bm{I}}_
{d}$
and $\omega=\lambda_{\max}\left({\mathbb{E}}\left[\left({\bm{S}}_{i}^{k}\right)^{
\top}{\bm{S}}_{i}^{k}\right]\right)-1.$
Then, the condition (4) reduces to

$$ $\gamma\leq\left[L\left(1+\sqrt{\frac{(1-p)\omega}{pn}}\right)\right]^{-1}.$ $$

The latter coincides with the stepsize condition of the convergence result of scalar MARINA.

###### Theorem 1 .

###### Remark 2 .

###### Remark 3 .

###### Remark 4 .

###### Remark 5 .

### 4.2 Optimizing the matrix stepsize

Now let us look at the right-hand side of (5).
We notice that it decreases in terms of the determinant of the stepsize matrix.
Therefore, one needs to solve the following optimization problem to find the optimal stepsize:

$$ minimize $\displaystyle\log\det({\bm{D}}^{-1})$ subject to $\displaystyle{\bm{D}}\text{ satisfying }\eqref{eq:stepsize-cond-thm1}.$ $$

The solution of this constrained minimization problem on $\mathbb{S}^{d}_{++}$ is not explicit.
In theory, one may show that the constraint (4) is convex and attempt to solve the problem numerically.
However, as stressed by Li et al. (), the similar stepsize condition for det-CGD is not easily computed using solvers like CVXPY (Diamond & Boyd, ).
Instead, we may relax the problem to certain linear subspaces of $\mathbb{S}^{d}_{++}$.
In particular, we fix a matrix ${\bm{W}}\in\mathbb{S}^{d}_{++}$, and define ${\bm{D}}:=\gamma{\bm{W}}$.
Then, the condition on the matrix ${\bm{D}}$ becomes a condition for the scalar $\gamma$, which is given in the following corollary.

###### Corollary 1 .

Let ${\bm{W}}\in\mathbb{S}^{d}_{++}$, defining ${\bm{D}}:=\gamma\cdot{\bm{W}}$, where $\gamma\in\mathbb{R}_{+}$. then the condition in (4) reduces to the following condition on $\gamma$

$$ $\gamma\leq\frac{2\lambda_{{\bm{W}}}}{1+\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{ W}},\mathcal{S}}\lambda_{{\bm{W}}}}},$ (6) $$

where $\Lambda_{{\bm{W}},\mathcal{S}}=\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_
{i}^{k}{\bm{W}}{\bm{S}}_{i}^{k}\right]-{\bm{W}}\right)$,
$\lambda_{{\bm{W}}}=\lambda_{\max}^{-1}\big{(}{\bm{W}}^{\frac{1}{2}}{\bm{L}}{
\bm{W}}^{\frac{1}{2}}\big{)}$, $\alpha=\frac{1-p}{np}$ and $\beta=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot
\lambda_{\max}\left({\bm{L}}^{-1}{\bm{L}}_{i}\right)$.

This means that for every fixed ${\bm{W}}$, we can find the optimal scaling coefficient $\gamma$.
In section Section 6, we will use this corollary to prove that a suboptimal matrix step size, determined in this efficient way, is already better than the optimal scalar step size.

Extension to det-CGD2.
A variant of det-CGD, called det-CGD2, was also proposed by Li et al. ().
This algorithm, has the same structure as det-CGD with the sketch and stepsize interchanged.
It was shown, that this algorithm has explicit stepsize condition in the single node setting.
In Appendix G, we propose the variance reduced extension of the distributed det-CGD2 following the MARINA scheme.

###### Corollary 1 .

## 5 DASHA-based variance reduction

In this section, we present our second algorithm based on DASHA.
The latter utilizes a different type of variance reduction based on momentum (MVR).
Compared to MARINA, dasha makes simpler optimization steps and does not require periodic synchronization with all the nodes.
Notice that one may further simplify the notations here used in the algorithm.
However, we keep it this way as it is consistent with (Tyurin & Richtárik, ).

Figure: Algorithm 2 det-DASHA

### 5.1 Theoretical guarantees

###### Theorem 2 .

Suppose that Assumptions 1 and 2 hold. Let us initialize $g_{i}^{0}=h_{i}^{0}=\nabla f_{i}(x^{0})$ for all $i\in[n]$ in Algorithm 2, and define

$$ $\displaystyle\Lambda_{{\bm{D}},\mathcal{S}}=\lambda_{\max}\left({\mathbb{E}} \left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right),\quad \omega_{{\bm{D}}}=\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda_{{\bm{D }},\mathcal{S}}.$ $$

If $a=\frac{1}{2\omega_{{\bm{D}}}+1}$, and the following condition on stepsize ${\bm{D}}\in\mathbb{S}^{d}_{++}$ is satisfied

$$ ${\bm{D}}^{-1}\succeq{\bm{L}}-\frac{4\lambda_{\max}\left({\bm{D}}\right)\omega_ {{\bm{D}}}\left(4\omega_{{\bm{D}}}+1\right)}{n^{2}}\sum_{i=1}^{n}\lambda_{\max }\left({\bm{L}}_{i}\right){\bm{L}}_{i},$ $$

then the following inequality holds for the iterates of Algorithm 2

$$ $\displaystyle{\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{{ \bm{D}}/\left(\det({\bm{D}})\right)^{1/d}}\right]\leq\frac{2(f(x^{0})-f^{\star })}{\det({\bm{D}})^{1/d}\cdot K}.$ $$

Here $\tilde{x}^{K}$ is chosen uniformly randomly from the first K iterates of the algorithm.

###### Remark 6 .

The term $\Lambda_{{\bm{D}},\mathcal{S}}$ can be viewed as the matrix version of $\gamma\cdot\omega$, where $\omega$ is associated with the sketch, and $\gamma$ is the scalar stepsize. On the other hand, the $\omega_{{\bm{D}}}$ is the extension of $\omega$ in matrix norm. Similar to Remark 5, plugging in scalar arguments in the algorithm, we recover the result from Tyurin & Richtárik ().

Following the same scheme as in Section 4, we choose ${\bm{D}}=\gamma_{{\bm{W}}}\cdot{\bm{W}}$, where ${\bm{W}}\in\mathbb{S}^{d}_{++}$.
Thus, for a fixed ${\bm{W}}$, we relax the problem of finding the optimal stepsize to the problem of finding the optimal scaling factor $\gamma_{{\bm{W}}}>0$.

###### Corollary 2 .

For a fixed ${\bm{W}}\in\mathbb{S}^{d}_{++}$, the optimal scaling factor $\gamma_{{\bm{W}}}\in\mathbb{R}_{+}$ is given by

$$ $\gamma_{{\bm{W}}}=\frac{2\lambda_{{\bm{W}}}}{1+\sqrt{1+16C_{{\bm{W}}}\lambda_{ \min}\left({\bm{L}}\right)\cdot\lambda_{{\bm{W}}}}},$ $$

where $C_{{\bm{W}}}:={\lambda_{\max}\left({\bm{W}}\right)\cdot\omega_{{\bm{W}}}\left(
4\omega_{{\bm{W}}}+1\right)}/n$,
and $\lambda_{{\bm{W}}}:=\lambda_{\max}^{-1}\left({\bm{L}}^{\frac{1}{2}}{\bm{W}}{
\bm{L}}^{\frac{1}{2}}\right)$.

We observe that the structure of the optimal scaling factor for obtained above is similar to the one obtained in Corollary 1.

###### Theorem 2 .

###### Remark 6 .

###### Corollary 2 .

##### The availability of 𝑳 𝑳 {\bm{L}} bold_italic_L :

For both det-MARINA and det-DASHA, in order to determine the matrix stepsize, the knowledge of ${\bm{L}}$ is needed, if ${\bm{L}}$ is known, better complexities are guaranteed.
When ${\bm{L}}$ is unknown, a closed-form solution can be obtained for generalized linear models.
In more general cases, ${\bm{L}}_{i}$ can be treated as hyperparameters and estimated using first-order information via a gradient-based method (Wang et al., ).
One can think of this as some type of preprocessing step, after which the matrices are learnt.

## 6 Complexities of the algorithms

### 6.1 det-MARINA

The following corollary formulates the iteration complexity for det-MARINA for ${\bm{W}}={\bm{L}}^{-1}$.

###### Corollary 3 .

If we take ${\bm{W}}={\bm{L}}^{-1}$, then the condition (6) on $\gamma$ is given by

$$ $\gamma\leq{2}\left(1+\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{ S}}}\right)^{-1}.$ (7) $$

In order to satisfy $\varepsilon$-stationarity, that is ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{
\det({\bm{D}})^{1/d}}}\right]\leq\varepsilon^{2}$, we require

$$ $K\geq\mathcal{O}\left(\frac{\Delta_{0}\cdot\det({\bm{L}})^{\frac{1}{d}}}{ \varepsilon^{2}}\cdot{\left(1+\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{L}}^{-1}, \mathcal{S}}}\right)}\right),$ $$

where $\Delta_{0}:=f(x^{0})-f(x^{\star})$.
Moreover, this iteration complexity is always better than the one of MARINA.

The proof can be found in the Appendix.
In fact, we can show that in cases where we fix ${\bm{W}}={\bm{I}}_{d}$ and ${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$, the same conclusion also holds, relevant details can be found in Section D.3.
This essentially means that det-MARINA always has a “larger” stepsize compared to MARINA, even if the stepsize is suboptimal for the sake of efficiency, which leads to a better iteration complexity.
In addition, because we are using the same compressor for those two algorithms, the communication complexity of det-MARINA is also provably better than that of MARINA.

In order to compute the communication complexity, we borrow the concept of expected
density from Gorbunov et al. ().

###### Definition 2 .

For a given sketch matrix ${\bm{S}}\in\mathbb{S}^{d}_{+}$, the expected density is defined as

$$ $\zeta_{{\bm{S}}}=\sup_{x\in\mathbb{R}^{d}}{\mathbb{E}}\left[\left\|{\bm{S}}x \right\|_{0}\right],$ $$

where $\left\|x\right\|_{0}$ denotes the number of non-zero components of $x\in\mathbb{R}^{d}$.

In particular, we have $\zeta_{\texttt{Rand-$\tau$}}=\tau$.
Below, we state the communication complexity of det-MARINA with ${\bm{W}}={\bm{L}}^{-1}$ and the Rand-$\tau$ compressor.

###### Corollary 4 .

Assume that we are using sketch ${\bm{S}}\sim\mathcal{S}$ with expected density $\zeta_{\mathcal{S}}$.
Suppose also we are running det-MARINA with probability $p$ and we use the optimal stepsize matrix with respect to ${\bm{W}}={\bm{L}}^{-1}$.
Then the overall communication complexity of the algorithm is given by $\mathcal{O}\big{(}(Kp+1)d+(1-p)K\zeta_{\mathcal{S}}\big{)}$.
Specifically, if we pick $p=\zeta_{\mathcal{S}}/d$,
then the communication complexity is given by

$$ $\mathcal{O}\left(d+\frac{\Delta_{0}\det({\bm{L}})^{\frac{1}{d}}}{\varepsilon^{ 2}}\left(\zeta_{\mathcal{S}}+\sqrt{\frac{\beta}{n}\Lambda_{{\bm{L}}^{-1}, \mathcal{S}}\zeta_{\mathcal{S}}(d-\zeta_{\mathcal{S}})}\right)\right).$ $$

Notice that in case where no compression is applied, the communication complexity reduces to $\mathcal{O}(\nicefrac{{d\Delta_{0}\cdot\det({\bm{L}})^{\frac{1}{d}}}}{{
\varepsilon^{2}}})$. The latter coincides with the rate of matrix stepsize GD (see (Li et al., )).
Therefore, the dependence on $\varepsilon$ is not possible to improve further since GD is optimal among first order methods (Carmon et al., ).

Figure: Figure 1: Comparison of DCGD with optimal scalar stepsize, det-CGD with matrix stepsize ${\bm{D}}^{*}_{3}$, MARINA with optimal scalar stepsize, DASHA with optimal scalar stepsize, det-MARINA with optimal stepsize ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ and det-DASHA with optimal stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$. Throughout the experiment, we are using Rand-$\tau$ sketch with $\tau=60$, and each algorithm is run for a fixed number of iterations $K=10000$. The $G_{K,{\bm{D}}}$ in the y-axis is defined in (63), which is the average squared matrix norm of the gradients.
Refer to caption: x1.png

###### Corollary 3 .

###### Definition 2 .

###### Corollary 4 .

### 6.2 Det-DASHA

The difference of compression mechanisms, does not allow to have a direct comparison of the complexities of these algorithms.
In particular, det-MARINA compresses the gradient difference with some probability $p$, while det-DASHA compresses the gradient difference with momentum in each iteration.

###### Corollary 5 .

If we pick ${\bm{D}}=\gamma_{{\bm{L}}^{-1}}\cdot{\bm{L}}^{-1}$, then in order to reach an $\varepsilon^{2}$ stationary point, det-DASHA needs $K$ iterations with

$$ $K\geq\frac{f(x^{0})-f^{\star}}{\det({\bm{L}})^{-\frac{1}{d}}\varepsilon^{2}} \left(1+\sqrt{1+16C_{{\bm{L}}^{-1}}\lambda_{\min}\left({\bm{L}}\right)}\right).$ $$

The following corollary compares the complexities of DASHA and det-DASHA.
For the sake of brevity, we defer the complexities and other details to the proof of this corollary.

###### Corollary 6 .

Suppose that the conditions in Theorem 2 hold, then compared to DASHA, det-DASHA with ${\bm{W}}={\bm{L}}^{-1}$ always has a better iteration complexity, therefore, communication complexity as well.

The following corollary suggests that the communication complexity of det-DASHA is better than that of det-MARINA,

###### Corollary 7 .

The iteration complexity of det-MARINA with $p=\nicefrac{{1}}{{(\omega_{{\bm{L}}^{-1}}+1)}}$ and det-DASHA with momentum $\nicefrac{{1}}{{(2\omega_{{\bm{L}}^{-1}}+1)}}$ is the same, therefore the communication complexity of det-DASHA is better than the communication complexity of det-MARINA.

This is expected since the same relation occurs between MARINA and DASHA as it is described by Tyurin & Richtárik (, Table 1). We refer the readers to Section E.2.1.

###### Corollary 5 .

###### Corollary 6 .

###### Corollary 7 .

## 7 Experiments

This section contains several plots which confirm our theoretical improvements on the existing methods.
Figure 1 shows that the performance in terms of communication complexity of det-DASHA and det-MARINA is better than their scalar counterpart DASHA and MARINA respectively.
This validates the efficiency of using a matrix stepsize over a scalar stepsize.
Further, we notice that det-DASHA and det-MARINA have better communication complexity in this case, compared to det-CGD.
This demonstrates the effectiveness of applying variance reduction.
Finally, as expected, det-DASHA has better communication complexity than det-MARINA.
We refer the readers to the appendix for more technical details of the experiments.

## 8 Future work

i) In this paper, we have only considered (linear) sketches as the compression operator. However, there exists a variety of compressors which are useful in practice that do not fall into this category. Extending det-CGD and det-MARINA for general unbiased compressors is a promising future work direction.
ii) Additionally, given recent successes with adaptive stepsizes (e.g., (Loizou et al., ; Orvieto et al., ; Schaipp et al., )), designing an adaptive matrix stepsize tailored to our case could be viable.
iii) Finally, recent advances suggest that server step sizes play a key role in accelerating federated learning algorithms (Jhunjhunwala et al., ; Li et al., ; Li & Richtárik, ). Designing a matrix version of the server step size could also be interesting.

## Appendix A Additional details

### A.1 Notations

The standard Euclidean norm on $\mathbb{R}^{d}$ is defined as $\left\|\cdot\right\|$. We use $\mathbb{S}^{d}_{++}$ (resp. $\mathbb{S}^{d}_{+}$) to denote the positive definite (resp. semi-definite) cone of dimension $d$. $\mathbb{S}^{d}$ is used to denote all symmetric matrices of dimension $d$. We use the notation ${\bm{I}}_{d}$ to denote the identity matrix of size $d\times d$, and ${\bm{O}}_{d}$ to denote the zero matrix of size $d\times d$. Given ${\bm{Q}}\in\mathbb{S}^{d}_{++}$ and $x\in\mathbb{R}^{d}$,

$$ $\left\|x\right\|_{{\bm{Q}}}:=\sqrt{x^{\top}{\bm{Q}}x}=\sqrt{\left\langle x,{ \bm{Q}}x\right\rangle},$ $$

where $\left\langle\cdot,\cdot\right\rangle$ is the standard Euclidean inner product on $\mathbb{R}^{d}$. For a matrix ${\bm{A}}\in\mathbb{S}^{d}$, we use $\lambda_{\max}\left({\bm{A}}\right)$ (resp. $\lambda_{\min}\left({\bm{A}}\right)$) to denote the largest (resp. smallest) eigenvalue of the matrix ${\bm{A}}$. For a function $f:\mathbb{R}^{d}\mapsto\mathbb{R}$, its gradient and its Hessian at a point $x\in\mathbb{R}^{d}$ are respectively denoted as $\nabla f(x)$ and $\nabla^{2}f(x)$. For the sketch matrices ${\bm{S}}^{k}_{i}$ used in the algorithm, we use the superscript $k$ to denote the iteration and subscript $i$ to denote the client, the matrix ${\bm{S}}^{k}_{i}$ is thus sampled for client $i$ in the $k$-th iteration from the same distribution $\mathcal{S}$. For any matrix ${\bm{A}}\in\mathbb{S}^{d}$, we use the notation $\operatorname{diag}\left({\bm{A}}\right)\in\mathbb{S}^{d}$ to denote the diagonal of matrix ${\bm{A}}$.

### A.2 Additional prior work

Numerous effective convex optimization techniques have been adapted for application in non-convex scenarios. Here’s a selection of these techniques, although it’s not an exhaustive list: adaptivity (Dvinskikh et al., ; Zhang et al., ), variance reduction (J Reddi et al., ; Li et al., ), and acceleration (Guminov et al., ). Of particular relevance to our work is the paper by Khaled & Richtárik (), which introduces a unified approach for analyzing stochastic gradient descent for non-convex objectives. A comprehensive overview of non-convex optimization can be found in (Jain et al., ; Danilova et al., ).

An illustrative example of a matrix stepsized method is Newton’s method, which has been a long-standing favorite in the optimization community (Gragg & Tapia, ; Miel, ; Yamamoto, ). However, the computational complexity involved in computing the stepsize as the inverse of the Hessian of the current iteration is substantial. Instead, quasi-Newton methods employ a readily computable estimator to replace the inverse Hessian (Broyden, ; Dennis & Moré, ; Al-Baali & Khalfan, ; Al-Baali et al., ).
An important direction of research that is relevant to our work, studies distributed second order methods. Here is a non-exhaustive list of papers in this area: (Wang et al., ; Crane & Roosta, ; Zhang et al., ; Islamov et al., ; Alimisis et al., ; Safaryan et al., ).

The Distributed Compressed Gradient Descent ( DCGD) algorithm, initially proposed by Khirirat et al. (), has seen improvements in various aspects, as documented in works such as (Li et al., ; Horváth et al., ).
Its variance reduced version with gradients shifts was studied by Shulgin & Richtárik () in the (strongly) convex setting.
Additionally, there exists a substantial body of literature on other federated learning algorithms employing unbiased compressors (Alistarh et al., ; Mishchenko et al., ; Gorbunov et al., ; Mishchenko et al., ; Maranjyan et al., ; Horváth et al., ).

Variance reduction techniques have gained significant attention in the context of stochastic batch gradient descent that is prevalent in machine learning.
Numerous algorithms have been developed in this regard, including well-known ones like SVRG (Johnson & Zhang, ), SAG (Schmidt et al., ), SDCA(Richtárik & Takáč, ), SAGA (Defazio et al., ), MISO (Mairal, ), and Katyusha (Allen-Zhu, ).
An overview of more advanced methods can be found in (Gower et al., ).
Notably, SVRG and Katyusha have been extended with loopless variants, namely L-SVRG and L-Katyusha (Kovalev et al., ; Qian et al., ). These loopless versions streamline the algorithms by eliminating the outer loop and introducing a biased coin-flip mechanism at each step.
This simplification eases both the algorithms’ structure and their analyses, while preserving their worst-case complexity bounds.
L-SVRG, in particular, offers the advantage of setting the exit probability from the outer loop independently of the condition number, thus, enhancing both robustness and practical efficiency.

This technique of coin flipping allows to obtain variance reduction for the CGD algorithm.
A relevant example is the DIANA algorithm proposed by Mishchenko et al. (). Its convergence was proved both in the convex and non-convex cases.
Later, MARINA (Gorbunov et al., ) obtained the optimal convergence rate, improving in communication complexity compared to all previous first order methods.
Finally, there is a line of work developing variance reduction in the federated setting using other methods and techniques (Chraibi et al., ; Hanzely & Richtárik, ; Dinh et al., ; Peng et al., ).

Another method to obtain variance reduction is based on momentum.
It was initially studied by Cutkosky & Orabona (), where they propose the STORM algorithm, which is a stochastic gradient descent algorithm with a momentum term for non-convex objectives.
They obtain stationarity guarantees using adaptive stepsizes with optimal convergence rates.
However, they require the variance of the stochastic gradient to be bounded by a constant, which is impractical.
Using momentum for variance reduction has since been widely studied (Liu et al., ; Khanduri et al., ; Tran-Dinh et al., ; Li et al., ).

## Appendix B Basic facts

In this section, we present some basic facts along with their proofs that will be used later in the analysis.

###### Fact 1 .

For two matrices ${\bm{A}},{\bm{B}}\in\mathbb{S}_{+}^{d}$, denote the $i$-th largest eigenvalues of ${\bm{A}},{\bm{B}}$ as $\lambda_{i}({\bm{A}}),\lambda_{i}({\bm{B}})$, if ${\bm{A}}\succeq{\bm{B}}$, then the following holds

$$ $\lambda_{i}({\bm{A}})\geq\lambda_{i}({\bm{B}}).$ (8) $$

###### Proof.

According to the Courant-Fischer theorem, we write

$$ $\displaystyle\lambda_{i}({\bm{B}})$ $\displaystyle=$ $\displaystyle\max_{S:\dim S=i}\min_{x\in S\backslash\{0\}}\frac{x^{\top}{\bm{B }}x}{x^{\top}x}.$ $$

Let $S^{i}_{\max}$ be a subspace of dimension $i$ where the maximum is attained, we then have

$$ $\displaystyle\lambda_{i}({\bm{B}})$ $\displaystyle=$ $\displaystyle\min_{x\in S^{i}_{\max}\backslash\{0\}}\frac{x^{\top}{\bm{B}}x}{x ^{\top}x}$ $\displaystyle\leq$ $\displaystyle\min_{x\in S^{i}_{\max}\backslash\{0\}}\frac{x^{\top}{\bm{A}}x}{x ^{\top}x}\leq\max_{S:\dim S=i}\min_{x\in S\backslash\{0\}}\frac{x^{\top}{\bm{A }}x}{x^{\top}x}=\lambda_{i}({\bm{A}}).$ $$

∎

The following is a generalization of the bias-variance decomposition for the matrix norm.

###### Fact 2 .

(Variance Decomposition)
Given a matrix ${\bm{M}}\in\mathbb{S}^{d}_{++}$, any vector $c\in\mathbb{R}^{d}$, and a random vector $x\in\mathbb{R}^{d}$ such that ${\mathbb{E}}\left[\left\|x\right\|\right]\leq+\infty$, the following bound holds

$$ ${\mathbb{E}}\left[\left\|x-{\mathbb{E}}\left[x\right]\right\|^{2}_{{\bm{M}}} \right]={\mathbb{E}}\left[\left\|x-c\right\|^{2}_{{\bm{M}}}\right]-\left\|{ \mathbb{E}}\left[x\right]-c\right\|^{2}_{{\bm{M}}}.$ (9) $$

###### Proof.

We have

$$ $\displaystyle{\mathbb{E}}\left[\left\|x-c\right\|^{2}_{{\bm{M}}}\right]$ $\displaystyle-\left\|{\mathbb{E}}\left[x\right]-c\right\|^{2}_{{\bm{M}}}$ $\displaystyle={\mathbb{E}}\left[x^{\top}{\bm{M}}x\right]-2{\mathbb{E}}\left[x \right]^{\top}{\bm{M}}c+c^{\top}{\bm{M}}c-{\mathbb{E}}\left[x\right]^{\top}{ \bm{M}}{\mathbb{E}}\left[x\right]+2{\mathbb{E}}\left[x\right]^{\top}{\bm{M}}c- c^{\top}{\bm{M}}c$ $\displaystyle={\mathbb{E}}\left[x^{\top}{\bm{M}}x\right]-{\mathbb{E}}\left[x \right]^{\top}{\bm{M}}{\mathbb{E}}\left[x\right]$ $\displaystyle={\mathbb{E}}\left[x^{\top}{\bm{M}}x\right]-2\cdot{\mathbb{E}} \left[x\right]^{\top}{\bm{M}}{\mathbb{E}}\left[x\right]+{\mathbb{E}}\left[x \right]^{\top}{\bm{M}}{\mathbb{E}}\left[x\right]$ $\displaystyle={\mathbb{E}}\left[\left\|x-{\mathbb{E}}\left[x\right]\right\|^{2 }_{{\bm{M}}}\right].$ $$

This completes the proof.
∎

###### Fact 3 .

The map $\left({\bm{A}},{\bm{B}},{\bm{X}}\right)\mapsto{\bm{A}}-{\bm{X}}{\bm{B}}^{-1}{
\bm{X}}$ is jointly concave on $\mathbb{S}^{d}_{+}\times\mathbb{S}^{d}_{++}\times\mathbb{S}^{d}$. It is also monotone increasing in variables ${\bm{A}}$ and ${\bm{B}}$.

We refer the reader to Corollary 1.5.3 of Bhatia () for the details and the proof.
The following is a result of 1 and 3.

###### Fact 4 .

Suppose ${\bm{L}}_{i}\in\mathbb{S}^{d}_{++}$, for $i=1,\ldots,n$.
Then, for every matrix ${\bm{X}}\in\mathbb{S}^{d}_{++}$, we define the following mapping

$$ $f({\bm{X}},{\bm{L}}_{1},\ldots,{\bm{L}}_{n})=\frac{1}{n}\sum_{i=1}^{n}\lambda_ {\max}({\bm{L}}_{i})\cdot\lambda_{\max}\left({\bm{L}}_{i}{\bm{X}}^{-1}\right) \cdot\lambda_{\max}\left({\bm{X}}^{-1}\right).$ $$

Then the above mapping is monotone decreasing in ${\bm{X}}$.

###### Proof.

First we notice that from 3 the mapping ${\bm{X}}\mapsto{\bm{X}}^{-1}$ is monotone decreasing.
The latter means that if we have any ${\bm{X}}_{1},{\bm{X}}_{2}\in\mathbb{S}^{d}_{++}$ such that ${\bm{X}}_{1}\succeq{\bm{X}}_{2}$, we have

$$ $\displaystyle{\bm{X}}_{1}^{-1}\preceq{\bm{X}}_{2}^{-1}.$ $$

Then it immediately follows, due to 1, that

$$ $\displaystyle 0<\lambda_{\max}({\bm{X}}_{1}^{-1})\leq\lambda_{\max}({\bm{X}}_{ 2}^{-1}).$ $$

We also notice that the relation $\lambda_{\max}\left({\bm{L}}_{i}{\bm{X}}^{-1}\right)=\lambda_{\max}\left({\bm{
L}}_{i}^{\frac{1}{2}}{\bm{X}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)=\lambda_{
\max}\left({\bm{X}}^{-1}{\bm{L}}_{i}\right)$, and that the mapping ${\bm{X}}\mapsto{\bm{L}}_{i}^{\frac{1}{2}}{\bm{X}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}$ is also monotone decreasing for every $i\in[n]$, so we have

$$ $0<\lambda_{\max}\left({\bm{L}}_{i}{\bm{X}}_{1}^{-1}\right)\leq\lambda_{\max} \left({\bm{L}}_{i}{\bm{X}}_{2}^{-1}\right).$ $$

Since we have the coefficient $\lambda_{\max}\left({\bm{L}}_{i}\right)>0$, it follows that,

$$ $\displaystyle f({\bm{X}}_{1},{\bm{L}}_{1},\ldots,{\bm{L}}_{n})\leq f({\bm{X}}_ {2},{\bm{L}}_{1},\ldots,{\bm{L}}_{n}).$ $$

This means that $f({\bm{X}})$ is monotone decreasing in ${\bm{X}}$.
∎

###### Fact 5 .

For any two matrices ${\bm{A}},{\bm{B}}\in\mathbb{S}^{d}_{++}$, the following relation regarding their largest eigenvalue holds

$$ $\lambda_{\max}\left({\bm{A}}{\bm{B}}\right)\leq\lambda_{\max}\left({\bm{A}} \right)\cdot\lambda_{\max}\left({\bm{B}}\right).$ (10) $$

###### Proof.

Using the Courant-Fischer theorem, we can write

$$ $\displaystyle\lambda_{\max}\left({\bm{A}}{\bm{B}}\right)$ $\displaystyle=$ $\displaystyle\min_{S:\dim S=d}\max_{x\in S\backslash\{0\}}\frac{x^{\top}{\bm{A }}{\bm{B}}x}{x^{\top}x}$ $\displaystyle=$ $\displaystyle\max_{x\in\mathbb{R}^{d}\backslash\{0\}}\frac{x^{\top}{\bm{A}}{ \bm{B}}x}{x^{\top}x}$ $\displaystyle\leq$ $\displaystyle\max_{x\in\mathbb{R}^{d}\backslash\{0\}}\frac{x^{\top}{\bm{A}}x}{ x^{\top}x}\cdot\max_{x\in\mathbb{R}^{d}\backslash\{0\}}\frac{x^{\top}{\bm{B}}x }{x^{\top}x}$ $\displaystyle=$ $\displaystyle\lambda_{\max}\left({\bm{A}}\right)\cdot\lambda_{\max}\left({\bm{ B}}\right).$ $$

∎

###### Fact 6 .

Given matrix ${\bm{Q}}\in\mathbb{S}^{d}_{++}$ and its matrix norm $\left\|\cdot\right\|_{{\bm{Q}}}$, its associated dual norm is $\left\|\cdot\right\|_{{\bm{Q}}^{-1}}$.

###### Proof.

Let us first recall the definition of the dual norm $\left\|\cdot\right\|_{*}$.
For any vector $z\in\mathbb{R}^{d}$, it is defined as

$$ $\left\|z\right\|_{*}:=\sup\{z^{\top}x:\left\|x\right\|_{{\bm{Q}}}\leq 1\}.$ $$

Solving this optimization problem is equivalent to solving $\sup\{z^{\top}x:\left\|x\right\|^{2}_{{\bm{Q}}}=1\}$.
The Lagrange function is given as

$$ $f(x,\lambda)=z^{\top}x-\lambda\left(\left\|x\right\|^{2}_{{\bm{Q}}}-1\right)=z ^{\top}x-\lambda\left(x^{\top}{\bm{Q}}x-1\right).$ $$

Computing the derivatives we deduce that

$$ $\frac{\partial f(x,\lambda)}{\partial x}=z-2\lambda\cdot{\bm{Q}}x=0,\qquad \frac{\partial f(x,\lambda)}{\partial\lambda}=\left\|x\right\|^{2}_{{\bm{Q}}}- 1=0.$ $$

This leads to

$$ $\lambda=\frac{\left\|z\right\|_{{\bm{Q}}^{-1}}}{2},\qquad x=\frac{{\bm{Q}}^{-1 }z}{\left\|z\right\|_{{\bm{Q}}^{-1}}}.$ $$

As a result, we have

$$ $\displaystyle\sup\{z^{\top}x:\left\|x\right\|_{{\bm{Q}}}\leq 1\}$ $\displaystyle=$ $\displaystyle\sup\{z^{\top}x:\left\|x\right\|^{2}_{{\bm{Q}}}=1\}$ $\displaystyle=$ $\displaystyle z^{\top}z=\frac{z^{\top}{\bm{Q}}^{-1}z}{\left\|z\right\|_{{\bm{Q }}^{-1}}}=\left\|z\right\|_{{\bm{Q}}^{-1}}.$ $$

∎

###### Fact 1 .

###### Proof.

###### Fact 2 .

###### Proof.

###### Fact 3 .

###### Fact 4 .

###### Proof.

###### Fact 5 .

###### Proof.

###### Fact 6 .

###### Proof.

## Appendix C Properties of matrix smoothness

### C.1 The matrix Lipschitz-continuous gradient

In this section we describe some properties of matrix smoothness, matrix gradient Lipschitzness and their relations. The following proposition describes a sufficient condition for the matrix Lipschitz-continuity of the gradient.

###### Proposition 2 .

Given twice continuously differentiable function $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ with bounded Hessian,

$$ $\nabla^{2}f(x)\preceq{\bm{L}},$ (11) $$

where ${\bm{L}}\in\mathbb{S}^{d}_{++}$ and the generalized inequality holds for any $x\in\mathbb{R}^{d}$. Then $f$ satisfies (3) with the matrix ${\bm{L}}$.

The below proposition is a variant of Proposition 1 and it characterizes the smoothness matrix of the objective function $f$, given the smoothness matrices of the component functions $f_{i}$.

###### Proposition 3 .

Assume that $f_{i}$ has ${\bm{L}}_{i}$-Lipschitz continuous gradient for every $i\in[n]$, then function $f$ has ${\bm{L}}$-Lipschitz gradient with ${\bm{L}}\in\mathbb{S}^{d}_{++}$ satisfying

$$ ${\bm{L}}\cdot\lambda_{\min}\left({\bm{L}}\right)=\frac{1}{n}\sum_{i=1}^{n} \lambda_{\max}\left({\bm{L}}_{i}\right)\cdot{\bm{L}}_{i}.$ (12) $$

###### Proposition 2 .

###### Proposition 3 .

#### C.1.1 Quadratics

Given a matrix ${\bm{A}}\in\mathbb{S}^{d}_{++}$ and a vector $b\in\mathbb{R}^{d}$, consider the function $f(x)=\frac{1}{2}x^{\top}{\bm{A}}x+b^{\top}x+c$.
Then its gradient is computed as $\nabla f(x)={\bm{A}}x+b$ and $\nabla^{2}f(x)={\bm{A}}$.
Inserting gradients formula into (3) we deduce

$$ $\sqrt{(x-y)^{\top}{\bm{A}}{\bm{L}}^{-1}{\bm{A}}(x-y)}\leq\sqrt{(x-y)^{\top}{ \bm{L}}(x-y)},$ $$

for any $x,y\in\mathbb{R}^{d}$. This reduces to

$$ ${\bm{A}}{\bm{L}}^{-1}{\bm{A}}\preceq{\bm{L}}.$ (13) $$

Since ${\bm{A}}\in\mathbb{S}^{d}_{++}$, we can also rewrite (13) as

$$ ${\bm{A}}^{\frac{1}{2}}{\bm{L}}^{-1}{\bm{A}}^{\frac{1}{2}}\preceq{\bm{A}}^{- \frac{1}{2}}{\bm{L}}{\bm{A}}^{-\frac{1}{2}},$ $$

which is equivalent to

$$ ${\bm{A}}\preceq{\bm{L}}.$ (14) $$

Therefore, the “best” ${\bm{L}}\in\mathbb{S}^{d}_{++}$ that satisfies (3) is ${\bm{L}}={\bm{A}}=\nabla^{2}f(x)$, for every $x\in\mathbb{R}^{d}$.
Now, let us look at a more general setting.
Consider $f$ given as follows,

$$ $f(x)=\sum_{i=1}^{s}\phi_{i}({\bm{M}}_{i}x),$ $$

where ${\bm{M}}_{i}\in\mathbb{R}^{q_{i}\times d}$. Here $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ is the sum of functions $\phi_{i}:\mathbb{R}^{q_{i}}\mapsto\mathbb{R}$. We assume that each function $\phi_{i}$ has matrix ${\bm{L}}_{i}$ Lipschitz gradient. We have the following lemma regarding the matrix gradient Lipschitzness of $f$.

###### Proposition 4 .

Assume that functions $f$ and $\{\phi_{i}\}_{i=1}^{s}$ are described above.
Then function $f$ has ${\bm{L}}$-Lipschitz gradient, if the following condition is satisfied:

$$ $\sum_{i=1}^{s}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}{\bm{M}}_{i}{\bm{L }}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}}\right)=1.$ (15) $$

Note that Proposition 4 is a generalization of the previous case of quadratics, if we pick $s=1$, ${\bm{M}}_{i}={\bm{A}}^{\frac{1}{2}}$ and $\phi_{1}(x)=x^{\top}{\bm{I}}_{d}x$, the condition becomes ${\bm{L}}={\bm{A}}$, which is exactly the solution given by (14).
Thus we recover the result for quadratics. The linear term $bx+c$ is ignored in this case. In Proposition 4, we only intend to give a way of finding a matrix ${\bm{L}}\in\mathbb{S}^{d}_{++}$, so that $f$ has ${\bm{L}}$-Lipschitz gradient. This does not mean, however, the ${\bm{L}}$ here is optimal.
The proof is deferred to Section C.3.

###### Proposition 4 .

### C.2 Comparison of the different smoothness conditions

Let us recall the definition of matrix smoothness.

###### Definition 3 .

(${\bm{L}}$-smoothness)
Assume that $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a continuously differentiable function and matrix ${\bm{L}}\in\mathbb{S}^{d}_{++}$. We say that $f$ is ${\bm{L}}$-smooth if for all $x,y\in\mathbb{R}^{d}$

$$ $f(y)\leq f(x)+\left\langle\nabla f(x),x-y\right\rangle+\frac{1}{2}\left\|x-y \right\|_{{\bm{L}}}^{2}.$ (16) $$

We provide a proposition here which describes an equivalent form of stating ${\bm{L}}$-matrix smoothness of a function $f$. This proposition is used to illustrate the relation between matrix smoothness and matrix Lipschitz gradient.

###### Proposition 5 .

Let function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be continuously differentiable. Then the following statements are equivalent.

- (i)
$f$ is ${\bm{L}}$-matrix smooth.
- (ii)
$\left\langle\nabla f(x)-\nabla f(y),x-y\right\rangle\leq\left\|x-y\right\|^{2}
_{{\bm{L}}}$ for all $x,y\in\mathbb{R}^{d}$.

The two propositions, Proposition 6 and Proposition 7, formulated below illustrate the relation between matrix smoothness of $f$ and matrix gradient Lipschitzness of $f$.

###### Proposition 6 .

Assume $f:\mathbb{R}^{d}\mapsto\mathbb{R}$ is a continuously differentiable function, and its gradient is ${\bm{L}}$-Lipschitz continuous with ${\bm{L}}\in\mathbb{S}_{++}^{d}$. Then function $f$ is ${\bm{L}}$-matrix smooth.

###### Proposition 7 .

Assume $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is a continuously differentiable function. Assume also that $f$ is convex and ${\bm{L}}$-matrix smooth. Then $\nabla f$ is ${\bm{L}}$-Lipschitz continuous.

The next proposition shows that standard Lipschitzness of the gradient of a function is an immediate consequence of matrix Lipschitzness.

###### Proposition 8 .

Assume that the gradient of $f$ is ${\bm{L}}$-Lipschitz continuous. Then $\nabla f$ is also $L$-Lipschitz with $L=\lambda_{\max}\left({\bm{L}}\right)$.

###### Definition 3 .

###### Proposition 5 .

###### Proposition 6 .

###### Proposition 7 .

###### Proposition 8 .

### C.3 Proofs of the propositions regarding smoothness

#### C.3.1 Proof of Proposition 1

We start with the definition of ${\bm{L}}$-Lipschitz gradient of function $f$, and pick
two arbitrary points $x,y\in\mathbb{R}^{d}$,

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle=$ $\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}\left(\nabla f_{i}(x)-\nabla f_{i }(y)\right)\right\|^{2}_{{\bm{L}}^{-1}}.$ $$

Applying the convexity of $\left\|\cdot\right\|^{2}_{{\bm{L}}^{-1}}$, we have

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\leq$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\|\nabla f_{i}(x)-\nabla f_{i}(y) \right\|^{2}_{{\bm{L}}^{-1}}.$ $$

For each term within the summation, we use the definition of matrix norms and replace the matrix ${\bm{L}}^{-1}$ with ${\bm{L}}_{i}^{-1/2}{\bm{L}}_{i}^{1/2}{\bm{L}}^{-1}{\bm{L}}_{i}^{1/2}{\bm{L}}_{
i}^{-1/2}$, for every $i=1,\ldots,n$:

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left({\bm{L}}_{i}^{-\frac{1}{2}}(\nabla f _{i}(x)-\nabla f_{i}(y))\right)^{\top}{\bm{L}}_{i}^{\frac{1}{2}}{\bm{L}}^{-1}{ \bm{L}}_{i}^{\frac{1}{2}}\left({\bm{L}}_{i}^{-\frac{1}{2}}(\nabla f_{i}(x)- \nabla f_{i}(y))\right)$ $\displaystyle\leq$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{ 1}{2}}{\bm{L}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|{\bm{L}}_{i}^{- \frac{1}{2}}(\nabla f_{i}(x)-\nabla f_{i}(y))\right\|^{2}$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{ 1}{2}}{\bm{L}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|\nabla f_{i}(x)- \nabla f_{i}(y)\right\|^{2}_{{\bm{L}}_{i}^{-1}}.$ $$

Using the assumption that the gradient of each function $f_{i}$ is ${\bm{L}}_{i}$-Lipschitz, we obtain,

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle{\leq}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{ 1}{2}}{\bm{L}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|x-y\right\|^{2}_{{ \bm{L}}_{i}}.$ $$

Replacing ${\bm{L}}_{i}^{-1}$ with ${\bm{L}}^{-1/2}{\bm{L}}^{1/2}{\bm{L}}_{i}^{-1}{\bm{L}}^{1/2}{\bm{L}}^{-1/2}$

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{ 1}{2}}{\bm{L}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)\cdot\left[({\bm{L}}^{ \frac{1}{2}}(x-y))^{\top}{\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1 }{2}}({\bm{L}}^{\frac{1}{2}}(x-y))\right]$ $\displaystyle\leq$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{ 1}{2}}{\bm{L}}^{-1}{\bm{L}}_{i}^{\frac{1}{2}}\right)\cdot\lambda_{\max}\left({ \bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right)\left\|{\bm{L} }^{\frac{1}{2}}(x-y)\right\|^{2}.$ $$

Using 5, we are deduce the following bound,

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\leq$ $\displaystyle\left(\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}^{-1} \right)\cdot\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot\lambda_{\max}\left({ \bm{L}}_{i}{\bm{L}}^{-1}\right)\right)\cdot\left\|x-y\right\|^{2}_{{\bm{L}}}$ $\displaystyle=$ $\displaystyle\left\|x-y\right\|^{2}_{{\bm{L}}}.$ $$

#### C.3.2 Proof of Proposition 2

We start with picking any two vector $x,y\in\mathbb{R}^{d}$. We have

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\quad=\left\|\int_{0}^{1}\nabla^{2}f(\theta x+(1-\theta)y)(x-y) \operatorname{d}\!{\theta}\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\quad=(x-y)^{\top}\left(\int_{0}^{1}\nabla^{2}f(\theta x+(1- \theta)y)\operatorname{d}\!{\theta}\right)^{\top}{\bm{L}}^{-1}\left(\int_{0}^{ 1}\nabla^{2}f(\theta x+(1-\theta)y)\operatorname{d}\!{\theta}\right)(x-y).$ $$

Denote ${\bm{F}}:=\int_{0}^{1}\nabla^{2}f(\theta x+(1-\theta)y)\operatorname{d}\!{\theta}$, notice that ${\bm{F}}$ is a symmetric matrix.
Then, the previous identity becomes

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle=(x-y)^{\top}{\bm{F}}^{\top}{\bm{L}}^{-1}{\bm{F}}(x-y).$ $$

From the definition of ${\bm{F}}$ and the bounded Hessian assumption, we have ${\bm{F}}\preceq{\bm{L}}$. Let us prove that ${\bm{F}}{\bm{L}}^{-1}{\bm{F}}\preceq{\bm{L}}$:

$$ $\displaystyle{\bm{F}}{\bm{L}}^{-1}{\bm{F}}\preceq{\bm{L}}$ $\displaystyle\iff$ $\displaystyle{\bm{L}}^{-\frac{1}{2}}{\bm{F}}{\bm{L}}{\bm{F}}{\bm{L}}^{-\frac{1 }{2}}\preceq{\bm{I}}_{d}$ $\displaystyle\iff$ $\displaystyle{\bm{L}}^{-\frac{1}{2}}{\bm{F}}{\bm{L}}^{-\frac{1}{2}}\cdot{\bm{L }}^{-\frac{1}{2}}{\bm{F}}{\bm{L}}^{-\frac{1}{2}}\preceq{\bm{I}}_{d}$ $\displaystyle\iff$ $\displaystyle{\bm{L}}^{-\frac{1}{2}}{\bm{F}}{\bm{L}}^{-\frac{1}{2}}\preceq{\bm {I}}_{d}$ $\displaystyle\iff$ $\displaystyle{\bm{F}}\preceq{\bm{L}}.$ $$

This means that

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\leq$ $\displaystyle(x-y)^{\top}{\bm{L}}(x-y)=\left\|x-y\right\|^{2}_{{\bm{L}}},$ $$

which completes the proof.

#### C.3.3 Proof of Proposition 3

Suppose ${\bm{L}}$ is a symmetric positive definite matrix satisfying (12). Let us now show that the function $\nabla f$ is ${\bm{L}}$-Lipschitz continuous.
We start with picking any two points $x,y\in\mathbb{R}^{d}$, and notice that

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}=\left \|\frac{1}{n}\sum_{i=1}^{n}\left(\nabla f_{i}(x)-\nabla f_{i}(y)\right)\right \|^{2}_{{\bm{L}}^{-1}}.$ $$

Applying Jensen’s inequality, we obtain

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}\leq \frac{1}{n}\sum_{i=1}^{n}\left\|\nabla f_{i}(x)-\nabla f_{i}(y)\right\|^{2}_{{ \bm{L}}^{-1}}.$ $$

We then re-weight the norm appears in the summation individually,

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\left(\nabla f_{i}(x)-\nabla f_{i}(y )\right)^{\top}{\bm{L}}_{i}^{-\frac{1}{2}}{\bm{L}}_{i}^{\frac{1}{2}}{\bm{L}}^{ -1}{\bm{L}}_{i}^{\frac{1}{2}}{\bm{L}}_{i}^{-\frac{1}{2}}\left(\nabla f_{i}(x)- \nabla f_{i}(y)\right)$ $\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i} \right)\cdot\lambda_{\max}\left({\bm{L}}^{-1}\right)\cdot\left\|\nabla f_{i}(x )-\nabla f_{i}(y)\right\|^{2}_{{\bm{L}}_{i}^{-1}}.$ $$

Utilizing the assumption that each $f_{i}$ has ${\bm{L}}_{i}$ Lipschitz gradient, we obtain

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}$ $\displaystyle\leq\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i} \right)\cdot\lambda_{\max}\left({\bm{L}}^{-1}\right)\cdot\left\|x-y\right\|^{2 }_{{\bm{L}}_{i}}$ $\displaystyle=\left\|x-y\right\|^{2}_{\lambda_{\max}\left({\bm{L}}^{-1}\right) \cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot{\bm {L}}_{i}}\overset{\eqref{dasha:eq:global-smooth}}{=}\left\|x-y\right\|^{2}_{{ \bm{L}}}.$ $$

#### C.3.4 Proof of Proposition 4

For any $x$ and $y$ from $\mathbb{R}^{d}$, we have

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}$ $\displaystyle\quad=\left\|\sum_{i=1}^{s}{\bm{M}}_{i}^{\top}\nabla\phi_{i}({\bm {M}}_{i}x)-\sum_{i=1}^{s}{\bm{M}}_{i}^{\top}\nabla\phi_{i}({\bm{M}}_{i}y) \right\|_{{\bm{L}}^{-1}}$ $\displaystyle\quad=s\cdot\left\|\frac{1}{s}\sum_{i=1}^{s}{\bm{M}}_{i}^{\top} \left(\nabla\phi_{i}({\bm{M}}_{i}x)-\nabla\phi_{i}\left({\bm{M}}_{i}y\right) \right)\right\|_{{\bm{L}}^{-1}}.$ $$

Applying the convexity of the norm $\left\|\cdot\right\|_{{\bm{L}}^{-1}}$,

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}$ $\displaystyle\leq$ $\displaystyle s\cdot\frac{1}{s}\sum_{i=1}^{s}\left\|{\bm{M}}_{i}^{\top}\left( \nabla\phi_{i}({\bm{M}}_{i}x)-\nabla\phi_{i}({\bm{M}}_{i}y)\right)\right\|_{{ \bm{L}}^{-1}}.$ $$

Expanding the norm and applying the replacement trick for above ${\bm{L}}$ and ${\bm{M}}_{i}$, we obtain

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}$ $\displaystyle=\sum_{i=1}^{s}\sqrt{\left(\nabla\phi_{i}({\bm{M}}_{i}x)-\nabla \phi_{i}({\bm{M}}_{i}y)\right)^{\top}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{ \top}\left(\nabla\phi_{i}({\bm{M}}_{i}x)-\nabla\phi_{i}({\bm{M}}_{i}y)\right)}$ $\displaystyle=\sum_{i=1}^{s}\sqrt{{\bm{B}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2} }{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}}{\bm{B} }_{i}}$ $\displaystyle\leq\sum_{i=1}^{s}\sqrt{\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1 }{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}} \right)}\cdot\left\|\nabla\phi_{i}({\bm{M}}_{i}x)-\nabla\phi_{i}({\bm{M}}_{i}y )\right\|_{{\bm{L}}_{i}^{-1}},$ $$

where ${\bm{B}}_{i}:={\bm{L}}_{i}^{-\frac{1}{2}}\left(\nabla\phi_{i}({\bm{M}}_{i}x)-
\nabla\phi_{i}({\bm{M}}_{i}y)\right)$.
Due to the assumption that the gradient of $\phi_{i}$ is ${\bm{L}}_{i}$-Lipschitz, we have

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}$ $\displaystyle\quad\leq\sum_{i=1}^{s}\sqrt{\lambda_{\max}\left({\bm{L}}_{i}^{ \frac{1}{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1} {2}}\right)}\cdot\left\|{\bm{M}}_{i}(x-y)\right\|_{{\bm{L}}_{i}}$ $\displaystyle\quad=\sum_{i=1}^{s}\sqrt{\lambda_{\max}\left({\bm{L}}_{i}^{\frac {1}{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}} \right)}\cdot\sqrt{\left[{\bm{L}}^{\frac{1}{2}}\left(x-y\right)\right]^{\top}{ \bm{L}}^{-\frac{1}{2}}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}{\bm{M}}_{i}{\bm{L}}^{- \frac{1}{2}}\left[{\bm{L}}^{\frac{1}{2}}(x-y)\right]}$ $\displaystyle\quad\leq\sum_{i=1}^{s}\sqrt{\lambda_{\max}\left({\bm{L}}_{i}^{ \frac{1}{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1} {2}}\right)\cdot\lambda_{\max}\left({\bm{L}}^{-\frac{1}{2}}{\bm{M}}_{i}^{\top} {\bm{L}}_{i}{\bm{M}}_{i}{\bm{L}}^{-\frac{1}{2}}\right)}\cdot\left\|x-y\right\| _{{\bm{L}}}$ $\displaystyle\quad\overset{}{\leq}\sum_{i=1}^{s}\lambda_{\max}\left({\bm{L}}_{ i}^{\frac{1}{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{ \frac{1}{2}}\right)\cdot\left\|x-y\right\|_{{\bm{L}}},$ $$

where the last inequality is due to the fact that,

$$ $\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}{\bm{M}}_{i}{\bm{L}}^{-1}{\bm{M} }_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}}\right)=\lambda_{\max}\left({\bm{L}}^{- \frac{1}{2}}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}{\bm{M}}_{i}{\bm{L}}^{-\frac{1}{2}} \right).$ $$

Recalling the condition of the proposition:

$$ $\sum_{i=1}^{s}\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}{\bm{M}}_{i}{\bm{L }}^{-1}{\bm{M}}_{i}^{\top}{\bm{L}}_{i}^{\frac{1}{2}}\right)=1,$ $$

we deduce

$$ $\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}\leq\left\|x-y\right\|_{ {\bm{L}}}.$ $$

#### C.3.5 Proof of Proposition 5

(i) $\rightarrow$ (ii). If $f$ is ${\bm{L}}$-matrix smooth, then for all $x,y\in\mathbb{R}^{d}$, we have

$$ $f(x)\leq f(y)+\left\langle\nabla f(y),x-y\right\rangle+\frac{1}{2}\left\|x-y \right\|^{2}_{{\bm{L}}},$ $$

and

$$ $f(y)\leq f(x)+\left\langle\nabla f(x),y-x\right\rangle+\frac{1}{2}\left\|x-y \right\|^{2}_{{\bm{L}}}.$ $$

Summing up these two inequalities we get

$$ $\left\langle\nabla f(x)-\nabla f(y),x-y\right\rangle\leq\left\|x-y\right\|^{2} _{{\bm{L}}}.$ $$

(ii) $\rightarrow$ (i). Choose any $x,y\in\mathbb{R}^{d}$, and define $z=x+t(y-x)$, then we have,

$$ $\displaystyle f(y)$ $\displaystyle=$ $\displaystyle f(x)+\int_{0}^{1}\left\langle\nabla f(x+t(y-x)),y-x\right\rangle \operatorname{d}\!{t}$ $\displaystyle=$ $\displaystyle f(x)+\int_{0}^{1}\left\langle\nabla f(z),y-x\right\rangle \operatorname{d}\!{t}$ $\displaystyle=$ $\displaystyle f(x)+\left\langle\nabla f(x),y-x\right\rangle+\int_{0}^{1}\left \langle\nabla f(z)-\nabla f(x),y-x\right\rangle\operatorname{d}\!{t}$ $\displaystyle=$ $\displaystyle f(x)+\left\langle\nabla f(x),y-x\right\rangle+\int_{0}^{1}\left \langle\nabla f(z)-\nabla f(x),z-x\right\rangle\cdot\frac{1}{t}\operatorname{d }\!{t}.$ $$

Using the assumption that for any $x,z\in\mathbb{R}^{d}$, we have

$$ $\left\langle\nabla f(z)-\nabla f(x),z-x\right\rangle\leq\left\|z-x\right\|^{2} _{{\bm{L}}}.$ $$

Plug this back into the previous identity, we obtain

$$ $\displaystyle f(y)$ $\displaystyle\leq$ $\displaystyle f(x)+\left\langle\nabla f(x),y-x\right\rangle+\int_{0}^{1}\left \|z-x\right\|^{2}_{{\bm{L}}}\cdot\frac{1}{t}\operatorname{d}\!{t}$ $\displaystyle=$ $\displaystyle f(x)+\left\langle\nabla f(x),y-x\right\rangle+\int_{0}^{1}\left \|y-x\right\|^{2}_{{\bm{L}}}\cdot t\operatorname{d}\!{t}$ $\displaystyle=$ $\displaystyle f(x)+\left\langle\nabla f(x),y-x\right\rangle+\frac{1}{2}\left\| y-x\right\|_{{\bm{L}}}^{2}.$ $$

#### C.3.6 Proof of Proposition 6

We start with picking any two points $x,y\in\mathbb{R}^{d}$, using the generalized Cauchy-Schwarz inequality for dual norm, we have

$$ $\displaystyle\left\langle\nabla f(x)-\nabla f(y),x-y\right\rangle$ $\displaystyle\leq$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}\cdot\left \|x-y\right\|_{{\bm{L}}}$ $\displaystyle\overset{\eqref{eq:assmp:3}}{\leq}$ $\displaystyle\left\|x-y\right\|_{{\bm{L}}}\cdot\left\|x-y\right\|_{{\bm{L}}}$ $\displaystyle=$ $\displaystyle\left\|x-y\right\|^{2}_{{\bm{L}}}$ $$

According to Proposition 5, this indicates that function $f$ is ${\bm{L}}$-matrix smooth.

#### C.3.7 Proof of Proposition 7

Using Proposition 5, we know that for any $x,y\in\mathbb{R}^{d}$, we have

$$ $\left\langle\nabla f(x)-\nabla f(y),x-y\right\rangle\leq\left\|x-y\right\|^{2} _{{\bm{L}}}.$ (17) $$

Now we pick any three points $x,y,z\in\mathbb{R}^{d}$. Using the ${\bm{L}}$-smoothness of $f$, we have

$$ $f(x+z)\geq f(x)+\left\langle\nabla f(x),z\right\rangle+\frac{1}{2}\left\|z \right\|^{2}_{{\bm{L}}}.$ (18) $$

Using the convexity of $f$ we have

$$ $\left\langle\nabla f(y),x+z-y\right\rangle\leq f(x+z)-f(y).$ (19) $$

Combining (18) and (19), we obtain

$$ $\left\langle\nabla f(y),x+z-y\right\rangle\leq f(x)-f(y)+\left\langle\nabla f( x),z\right\rangle+\frac{1}{2}\left\|z\right\|^{2}_{{\bm{L}}}.$ $$

Rearranging terms we get

$$ $\left\langle\nabla f(y)-\nabla f(x),z\right\rangle-\frac{1}{2}\left\|z\right\| ^{2}_{{\bm{L}}}\leq f(x)-f(y)-\left\langle\nabla f(y),x-y\right\rangle.$ $$

The inequality holds for any $z$ for fixed $x$ and $y$, and the left hand side is maximized (w.r.t. $z$) when $z={\bm{L}}^{-1}\left(\nabla f(y)-\nabla f(x)\right)$. Plugging it in, we get

$$ $\frac{1}{2}\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}\leq f(x) -f(y)-\left\langle\nabla f(y),x-y\right\rangle.$ (20) $$

By symmetry we can also obtain

$$ $\frac{1}{2}\left\|\nabla f(y)-\nabla f(x)\right\|^{2}_{{\bm{L}}^{-1}}\leq f(y) -f(x)-\left\langle\nabla f(x),y-x\right\rangle.$ $$

Adding (20) and its counterpart together, we get

$$ $\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}\leq\left\langle \nabla f(x)-\nabla f(y),x-y\right\rangle.$ (21) $$

Combing (21) and (17), it follows

$$ $\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}\leq\left\|x-y\right \|^{2}_{{\bm{L}}}.$ $$

Note that ${\bm{L}}$ and ${\bm{L}}^{-1}$ are both positive definite matrices, so it is equivalent to

$$ $\left\|\nabla f(x)-\nabla f(y)\right\|_{{\bm{L}}^{-1}}\leq\left\|x-y\right\|_{ {\bm{L}}}.$ $$

This completes the proof.

#### C.3.8 Proof of Proposition 8

Let us start with picking any two points $x,y\in\mathbb{R}^{d}$. With the matrix ${\bm{L}}$-Lipschitzness of the gradient of function $f$, we have

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}_{{\bm{L}}^{-1}}\leq \left\|x-y\right\|^{2}_{{\bm{L}}}.$ $$

This implies

$$ $(x-y)^{\top}{\bm{L}}(x-y)-\left(\nabla f(x)-\nabla f(y)\right)^{\top}{\bm{L}}^ {-1}\left(\nabla f(x)-\nabla f(y)\right)\geq 0.$ $$

Define function $f({\bm{X}}):=a^{\top}{\bm{X}}a-b^{\top}{\bm{X}}^{-1}b$ for ${\bm{X}}\in\mathbb{S}^{d}_{++}$, where $a,b\in\mathbb{R}^{d}$ are fixed vectors. Then $f$ is monotone increasing in ${\bm{X}}$. This can be shown in the following way, picking two matrices ${\bm{X}}_{1},{\bm{X}}_{2}\in\mathbb{S}^{d}_{++}$, where ${\bm{X}}_{1}\succeq{\bm{X}}_{2}$. It is easy to see that $-{\bm{X}}_{1}^{-1}\succeq-{\bm{X}}_{2}^{-1}$, since from 3 the map ${\bm{X}}\mapsto-{\bm{X}}^{-1}$ is monotone increasing for ${\bm{X}}\in\mathbb{S}^{d}_{++}$. Thus,

$$ $\displaystyle f({\bm{X}}_{1})-f({\bm{X}}_{2})$ $\displaystyle=(x-y)^{\top}\left({\bm{X}}_{1}-{\bm{X}}_{2}\right)(x-y)$ $\displaystyle\qquad+\left(\nabla f(x)-\nabla f(y)\right)^{\top}\left(-{\bm{X}} _{1}^{-1}-(-{\bm{X}}_{2}^{-1})\right)\left(\nabla f(x)-\nabla f(y)\right)\geq 0.$ $$

As a result, $f(\lambda_{\max}\left({\bm{L}}\right)\cdot{\bm{I}}_{d})\geq f({\bm{L}})\geq 0$, due to the fact that $\lambda_{\max}\left({\bm{L}}\right)\cdot{\bm{I}}_{d}\succeq{\bm{L}}$. It remains to notice that

$$ $\displaystyle f(\lambda_{\max}\left({\bm{L}}\right)\cdot{\bm{I}}_{d})$ $\displaystyle=\lambda_{\max}\left({\bm{L}}\right)\left\|x-y\right\|^{2}-\frac{ 1}{\lambda_{\max}\left({\bm{L}}\right)}\left\|\nabla f(x)-\nabla f(y)\right\|^ {2}\geq 0,$ $$

which yields

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|^{2}\leq\lambda^{2}_{\max} \left({\bm{L}}\right)\left\|x-y\right\|^{2}.$ $$

Since we are working with ${\bm{L}}\in\mathbb{S}^{d}_{++}$, the above inequality implies

$$ $\displaystyle\left\|\nabla f(x)-\nabla f(y)\right\|\leq\lambda_{\max}\left({ \bm{L}}\right)\left\|x-y\right\|.$ $$

## Appendix D Analysis of det-MARINA

### D.1 Technical lemmas

We first state some technical lemmas.

###### Lemma 1 (Descent lemma) .

Assume that function $f$ is ${\bm{L}}$ smooth, and $x^{k+1}=x^{k}-{\bm{D}}\cdot g^{k}$, where ${\bm{D}}\in\mathbb{S}_{++}^{d}$. Then we will have

$$ $f(x^{k+1})\leq f(x^{k})-\frac{1}{2}\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{ 2}+\frac{1}{2}\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}-\frac{1}{2} \left\|x^{k+1}-x^{k}\right\|_{{\bm{D}}^{-1}-{\bm{L}}}.$ $$

The following lemma is obtained for any sketch matrix ${\bm{S}}\in\mathbb{S}_{+}^{d}$ and any two positive definite matrices ${\bm{D}}$ and ${\bm{L}}$.

###### Lemma 2 (Property of sketch matrix) .

For any sketch matrix ${\bm{S}}\in\mathbb{S}^{d}_{+}$, a vector $t\in\mathbb{R}^{d}$, and matrices ${\bm{D}},{\bm{L}}\in\mathbb{S}^{d}_{++}$, we have

$$ ${\mathbb{E}}\left[\left\|{\bm{S}}t-t\right\|_{{\bm{D}}}^{2}\right]\leq\lambda_ {\max}\left({\bm{L}}^{\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}{\bm{D}}{\bm {S}}\right]-{\bm{D}}\right){\bm{L}}^{\frac{1}{2}}\right)\cdot\left\|t\right\|^ {2}_{{\bm{L}}^{-1}}.$ (22) $$

###### Lemma 3 .

Assume that Definition 1 holds and $h_{i}^{0}=\nabla f_{i}(x^{0})$, then for $h_{i}^{k+1}$ from Algorithm 2, we have for any ${\bm{D}}\in\mathbb{S}^{d}_{++}$

$$ $\displaystyle\left\|h^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}}$ $\displaystyle=\left\|h_{i}^{k+1}-\nabla f_{i}(x^{k+1})\right\|^{2}_{{\bm{D}}}=0,$ $$

and

$$ $\displaystyle\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}$ $\displaystyle\leq\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}_{i}}.$ $$

The following lemmas describe the recurrence applied to terms in the Lyapunov function.

###### Lemma 4 .

Suppose $h^{k+1}$ and $g^{k+1}$ are from Algorithm 2, then the following recurrence relation holds,

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad\leq\frac{2\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{\max} \left({\bm{D}}^{-1}\right)\cdot\lambda_{\max}\left({\bm{D}}\right)}{n^{2}}\sum _{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right){\mathbb{E}}\left[\left\|h_{i }^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}\right]$ $\displaystyle\qquad+\frac{2a^{2}\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{ \max}\left({\bm{D}}^{-1}\right)}{n^{2}}\sum_{i=1}^{n}{\mathbb{E}}\left[\left\| g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]+(1-a)^{2}{\mathbb{E}}\left[ \left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}\right],$ (23) $$

where $\Lambda_{{\bm{D}},\mathcal{S}}=\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_
{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right)$ for ${\bm{D}}\in\mathbb{S}^{d}_{++}$ and ${\bm{S}}_{i}^{k}\sim\mathcal{S}$.

###### Lemma 5 .

Suppose $h_{i}^{k+1}$ and $g_{i}^{k+1}$ for $i\in[n]$ are from Algorithm 2, then the following recurrence holds,

$$ $\displaystyle{\mathbb{E}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1}\right\|^{2}_{{ \bm{D}}}\right]$ $\displaystyle\quad\leq\left(2a^{2}\lambda_{\max}\left({\bm{D}}^{-1}\right) \cdot\Lambda_{{\bm{D}},\mathcal{S}}+(1-a)^{2}\right)\cdot{\mathbb{E}}\left[ \left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+2\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\lambda_{ \max}\left({\bm{D}}\right)\cdot\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{ \max}\left({\bm{L}}_{i}\right)\cdot{\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^ {k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}\right].$ $$

###### Lemma 1 (Descent lemma) .

###### Lemma 2 (Property of sketch matrix) .

###### Lemma 3 .

###### Lemma 4 .

###### Lemma 5 .

### D.2 Proof of Theorem 1

According to Lemma 1, we have

$$ $\displaystyle{\mathbb{E}}\left[f(x^{k+1})\right]$ $\displaystyle\leq{\mathbb{E}}\left[f(x^{k})\right]-{\mathbb{E}}\left[\frac{1}{ 2}\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right]+{\mathbb{E}}\left[\frac {1}{2}\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right]$ $\displaystyle\quad-{\mathbb{E}}\left[\frac{1}{2}\left\|x^{k+1}-x^{k}\right\|^{ 2}_{{\bm{D}}^{-1}-{\bm{L}}}\right].$ (24) $$

We then use the definition of $g^{k+1}$ to derive an upper bound for ${\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}}\right]$. Notice that,

$$ $g^{k+1}=\begin{cases}\nabla f(x^{k+1})&\text{ with probability }p,\\ g^{k}+\frac{1}{n}\sum_{i=1}^{n}{\bm{S}}_{i}^{k}\left(\nabla f_{i}(x^{k+1})- \nabla f_{i}(x^{k})\right)&\text{ with probability }1-p.\end{cases}$ $$

As a result, from the tower property,

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\quad={\mathbb{E}}\left[{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f (x^{k+1})\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k},c_{k}\right]\right]$ $\displaystyle\quad=p\cdot\left\|\nabla f(x^{k+1})-\nabla f(x^{k+1})\right\|^{2 }_{{\bm{D}}}$ $\displaystyle\qquad+(1-p)\cdot{\mathbb{E}}\left[\left\|g^{k}+\frac{1}{n}\sum_{ i=1}^{n}{\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))-\nabla f(x ^{k+1})\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\quad=(1-p)\cdot{\mathbb{E}}\left[\left\|g^{k}+\frac{1}{n}\sum_{i =1}^{n}{\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))-\nabla f(x^ {k+1})\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k}\right].$ $$

Using 2, we have

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle=(1-p)\cdot{\mathbb{E}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n}{\bm {S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))-\left(\nabla f(x^{k+1} )-\nabla f(x^{k})\right)\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\qquad+(1-p)\cdot\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D }}}$ $\displaystyle=(1-p)\cdot{\mathbb{E}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n} \left({\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))-(\nabla f_{i }(x^{k+1})-\nabla f_{i}(x^{k}))\right)\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k }\right]$ $\displaystyle\qquad+(1-p)\cdot\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D }}}.$ $$

Notice that the sketch matrix is unbiased, thus we have

$$ ${\mathbb{E}}\left[{\bm{S}}^{k}_{i}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{ k})\right)\mid x^{k+1},x^{k}\right]=\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}),$ $$

and any two random vectors in the set $\{{\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))\}_{i=1}^{n}$ are independent from each other, if $x^{k+1}$ and $x^{k}$ are fixed.
Therefore, we have

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\quad=\frac{1-p}{n^{2}}\sum_{i=1}^{n}{\mathbb{E}}\left[\left\|{ \bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))-(\nabla f_{i}(x^{k+ 1})-\nabla f_{i}(x^{k}))\right\|^{2}_{{\bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\qquad+(1-p)\cdot\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^ {2}.$ (25) $$

Lemma 2 yields

$$ $\displaystyle{\mathbb{E}}\left[\left\|{\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})- \nabla f_{i}(x^{k}))-(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))\right\|^{2}_{ {\bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\qquad\leq\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}\left({ \mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}} \right){\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|\nabla f_{i}(x^{k+1})-\nabla f_ {i}(x^{k})\right\|^{2}_{{\bm{L}}_{i}^{-1}}.$ (26) $$

2 implies

$$ $\displaystyle{\mathbb{E}}\left[\left\|{\bm{S}}_{i}^{k}(\nabla f_{i}(x^{k+1})- \nabla f_{i}(x^{k}))-(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))\right\|^{2}_{ {\bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\quad\leq\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}\left({ \mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}} \right){\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm {L}}_{i}}.$ (27) $$

Plugging (D.2) into (D.2), we deduce

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{ i}^{\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{ k}\right]-{\bm{D}}\right){\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|x^{k+1}-x^{k} \right\|^{2}_{{\bm{L}}_{i}}+(1-p)\cdot\left\|g^{k}-\nabla f(x^{k})\right\|^{2} _{{\bm{D}}}.$ $$

Replacing ${\bm{L}}_{i}^{-1}$ with ${\bm{L}}^{-1/2}{\bm{L}}^{1/2}{\bm{L}}_{i}^{-1}{\bm{L}}^{1/2}{\bm{L}}^{-1/2}$, we denote that

$$ $\lambda_{i}:=\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}\left({\mathbb{E}} \left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right){\bm{L}}_ {i}^{\frac{1}{2}}\right),$ $$

and rewrite the ${\bm{L}}_{i}$-norm in the first term of RHS by the ${\bm{L}}$-norm:

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\mid x^{k+1},x^{k}\right]$ $\displaystyle=\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{i}\cdot\left({\bm{L}}^{ \frac{1}{2}}(x^{k+1}-x^{k})\right)^{\top}{\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{ \bm{L}}^{-\frac{1}{2}}\left({\bm{L}}^{\frac{1}{2}}(x^{k+1}-x^{k})\right)$ $\displaystyle\quad+(1-p)\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}$ $\displaystyle\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{i}\cdot\lambda_{\max} \left({\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right)\left\| x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}+(1-p)\cdot\left\|g^{k}-\nabla f(x^{k}) \right\|^{2}_{{\bm{D}}}.$ $$

We further use 5 to upper bound $\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}
_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right){\bm{L}}_{i}^{\frac{1}{
2}}\right)$ by the product of $\lambda_{\max}\left({\bm{L}}_{i}\right)$ and $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}
\right]-{\bm{D}}\right)$.
This allows us to simplify the expression since $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}
\right]-{\bm{D}}\right)$ is independent of the index $i$.
Notice that we have already defined

$$ $R({\bm{D}},\mathcal{S})=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\mathbb{ E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right)\cdot \lambda_{\max}\left({\bm{L}}_{i}\right)\cdot\lambda_{\max}\left({\bm{L}}^{- \frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right).$ $$

Taking expectation, using tower property and using the definition above, we deduce

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-\nabla f(x^{k+1})\right\|^{2}_{{ \bm{D}}}\right]$ $\displaystyle\quad\leq\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})}{n}{\mathbb{E}} \left[\left\|x^{k+1}-x^{k}\right\|_{{\bm{L}}}^{2}\right]+(1-p){\mathbb{E}} \left[\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right].$ (28) $$

We construct the following Lyapunov function $\Phi_{k}$,

$$ $\Phi_{k}=f(x^{k})-f^{\star}+\frac{1}{2p}\left\|g^{k}-\nabla f(x^{k})\right\|_{ {\bm{D}}}^{2}.$ (29) $$

Using (D.2) and (D.2), we are able to get

$$ $\displaystyle{\mathbb{E}}\left[\Phi_{k+1}\right]$ $\displaystyle\leq\frac{1}{2p}\left[\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})}{n }{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]+(1-p) \cdot{\mathbb{E}}\left[\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}} \right]\right]$ $\displaystyle\qquad+{\mathbb{E}}\left[f(x^{k})-f^{\star}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]+\frac{1} {2}{\mathbb{E}}\left[\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right]$ $\displaystyle\qquad-\frac{1}{2}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^ {2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]$ $\displaystyle\quad={\mathbb{E}}\left[\Phi_{k}\right]-\frac{1}{2}{\mathbb{E}} \left[\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right]$ $\displaystyle\qquad+\left(\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})}{2np}{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]-\frac{1}{2 }{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|_{{\bm{D}}^{-1}-{\bm{L}}}^{2} \right]\right)$ $\displaystyle\quad={\mathbb{E}}\left[\Phi_{k}\right]-\frac{1}{2}{\mathbb{E}} \left[\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}\right]$ $\displaystyle\qquad+\frac{1}{2}\left(\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})} {np}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]-{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|_{{\bm{D}}^{-1}-{\bm{L}}}^{2} \right]\right).$ $$

We can rewrite the last term as

$$ ${\mathbb{E}}\left[(x^{k+1}-x^{k})^{\top}\left[\frac{(1-p)\cdot R({\bm{D}}, \mathcal{S})}{np}{\bm{L}}+{\bm{L}}-{\bm{D}}^{-1}\right](x^{k+1}-x^{k})\right].$ (30) $$

We require the matrix in between to be negative semi-definite, which is

$$ ${\bm{D}}^{-1}\succeq\left(\frac{(1-p)\cdot R({\bm{D}},\mathcal{S})}{np}+1 \right){\bm{L}}.$ $$

This leads to the result that the expression (30) is always non-positive.
After dropping the last term, the relation between ${\mathbb{E}}\left[\Phi_{k+1}\right]$ and ${\mathbb{E}}\left[\Phi_{k}\right]$ becomes

$$ ${\mathbb{E}}\left[\Phi_{k+1}\right]\leq{\mathbb{E}}\left[\Phi_{k}\right]-\frac {1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right].$ $$

Unrolling this recurrence, we get

$$ $\frac{1}{K}\sum_{k=0}^{K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2 }_{{\bm{D}}}\right]\leq\frac{2\left({\mathbb{E}}\left[\Phi_{0}\right]-{\mathbb {E}}\left[\Phi_{K}\right]\right)}{K}.$ (31) $$

The left hand side can viewed as ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{{\bm{D}}}\right]$,
where $\tilde{x}^{K}$ is drawn uniformly at random from $\{x_{k}\}_{k=0}^{K-1}$.
From $\Phi_{K}>0$, we obtain

$$ $\displaystyle\frac{2\left({\mathbb{E}}\left[\Phi_{0}\right]-{\mathbb{E}}\left[ \Phi_{K}\right]\right)}{K}$ $\displaystyle\leq$ $\displaystyle\frac{2\Phi_{0}}{K}$ $\displaystyle=$ $\displaystyle\frac{2\left(f(x^{0})-f^{\star}+\frac{1}{2p}\left\|g^{0}-\nabla f (x^{0})\right\|^{2}_{{\bm{D}}}\right)}{K}$ $\displaystyle=$ $\displaystyle\frac{2\left(f(x^{0})-f^{\star}\right)}{K}.$ $$

Plugging in the simplified result into (31), and performing determinant normalization, we get

$$ ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{ \det({\bm{D}})^{1/d}}}\right]\leq\frac{2\left(f(x^{0})-f^{\star}\right)}{\det( {\bm{D}})^{1/d}K}.$ (32) $$

###### Remark 7 .

We can achieve a slightly more refined stepsize condition than (4) for det-MARINA, which is given as follows

$$ ${\bm{D}}\succeq\left(\frac{(1-p)\cdot\tilde{R}({\bm{D}},\mathcal{S})}{np}+1 \right){\bm{L}},$ (33) $$

where

$$ $\tilde{R}({\bm{D}},\mathcal{S}):=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left( {\bm{L}}_{i}^{\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm {S}}_{i}^{k}\right]-{\bm{D}}\right){\bm{L}}_{i}^{\frac{1}{2}}\right)\cdot \lambda_{\max}\left({\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}} \right).$ $$

This is obtained if we do not use 5 to upper bound $\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}
_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right){\bm{L}}_{i}^{\frac{1}{
2}}\right)$ by the product of $\lambda_{\max}\left({\bm{L}}_{i}\right)$ and $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{k}
\right]-{\bm{D}}\right)$.
However, (33) results in a condition that is much harder to solve even if we assume ${\bm{D}}=\gamma\cdot{\bm{W}}$.
So instead of using the more refined condition (33), we turn to (4).
Notice that both of the two conditions (33) and (4) reduce to the stepsize condition for MARINA in the scalar setting.

###### Remark 7 .

### D.3 Comparison of different stepsizes

In Corollary 3, we focus on the special stepsize where we fix ${\bm{W}}={\bm{L}}^{-1}$, and show that in this case det-MARINA always beats MARINA in terms of both iteration and communication complexities.
However, other choices for ${\bm{W}}$ are also possible.
Specifically, we consider the cases where
${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$ and ${\bm{W}}={\bm{I}}_{d}$.

#### D.3.1 The diagonal case

We consider ${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$.
The following corollary describes the optimal stepsize and the iteration complexity.

###### Corollary 8 .

If we take ${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$ in Corollary 1, then the optimal stepsize satisfies

$$ ${\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}=\frac{2}{1+\sqrt{1+4\alpha \beta\cdot\Lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right),\mathcal{S}}} }\cdot\operatorname{diag}^{-1}\left({\bm{L}}\right).$ (34) $$

This stepsize results in a better iteration complexity of det-MARINA compared to scalar MARINA.

From this corollary we know that det-MARINA has a better iteration complexity when ${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$.
And since the same sketch is used for MARINA and det-MARINA, the communication complexity is improved as well.
However, in general there is no clear relation between the iteration complexity of ${\bm{W}}={\bm{L}}^{-1}$ case and ${\bm{W}}=\operatorname{diag}^{-1}\left({\bm{L}}\right)$ case.
This is also confirmed by one of our experiments, see Figure 6 to see the comparison of det-MARINA using optimal stepsizes in different cases.

###### Corollary 8 .

#### D.3.2 The identity case

In this setting, ${\bm{W}}$ is the $d$-dimensional identity matrix ${\bm{I}}_{d}$.
Then the stepsize of our algorithm reduces to a scalar $\gamma$, where $\gamma$ is determined through Corollary 1.
Notice that in this case we do not reduce to the standard MARINA case because we are still using the matrix Lipschitz gradient assumption with ${\bm{L}}\in\mathbb{S}^{d}_{++}$.

###### Corollary 9 .

If we take ${\bm{W}}={\bm{I}}_{d}$, the optimal stepsize is given by

$$ ${\bm{D}}^{*}_{{\bm{I}}_{d}}=\frac{2}{1+\sqrt{1+4\alpha\beta\frac{1}{\lambda_{ \max}\left({\bm{L}}\right)}\cdot\omega}}\cdot\frac{{\bm{I}}_{d}}{\lambda_{\max }\left({\bm{L}}\right)}.$ (35) $$

This stepsize results in a better iteration complexity of det-MARINA compared to scalar MARINA.

The result in this corollary tells us that using scalar stepsize with matrix Lipschitz gradient assumption alone can result in acceleration of MARINA.
However, the use of matrix stepsize allows us to also take into consideration the ”structure” of the stepsize, thus allows more flexibility.
When the structure of the stepsize is chosen properly, combining matrix gradient Lipschitzness and matrix stepsize can result in a faster rate, as it can also be observed from the experiments in Figure 6.
The choices of ${\bm{W}}$ we consider here are in some sense inspired by the matrix stepsize GD, where the optimal stepsize is ${\bm{L}}^{-1}$.
In general, how to identify the best structure for the matrix stepsize remains a open problem.

###### Corollary 9 .

### D.4 Proofs of the corollaries

#### D.4.1 Proof of Corollary 1

We start with rewriting (4) as

$$ $\displaystyle\left(\frac{1-p}{np}\cdot R\left({\bm{D}},\mathcal{S}\right)+1 \right){\bm{D}}^{\frac{1}{2}}{\bm{L}}{\bm{D}}^{\frac{1}{2}}\preceq{\bm{I}}_{d}.$ $$

Plugging in the definition of ${R}({\bm{D}},\mathcal{S})$ and ${\bm{D}}=\gamma{\bm{W}}$, we get

$$ $\displaystyle\gamma\left(\frac{1-p}{np}\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{ \max}\left({\bm{L}}_{i}\right)\lambda_{\max}\left({\bm{L}}^{-1}{\bm{L}}_{i} \right)\cdot\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{W}}{\bm {S}}_{i}^{k}\right]-{\bm{W}}\right)\cdot\gamma+1\right){\bm{W}}^{\frac{1}{2}}{ \bm{L}}{\bm{W}}^{\frac{1}{2}}\preceq{\bm{I}}_{d}.$ $$

This generalized inequality is equivalent to the following inequality,

$$ $\displaystyle\gamma\left(\frac{1-p}{np}\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{ \max}\left({\bm{L}}_{i}\right)\lambda_{\max}\left({\bm{L}}^{-1}{\bm{L}}_{i} \right)\cdot\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{W}}{\bm {S}}_{i}^{k}\right]-{\bm{W}}\right)\cdot\gamma+1\right)\cdot\lambda_{\max} \left({\bm{W}}^{\frac{1}{2}}{\bm{L}}{\bm{W}}^{\frac{1}{2}}\right)\leq 1,$ $$

which is a quadratic inequality on $\gamma$.
Notice that we have already defined

$$ $\displaystyle\alpha=\frac{1-p}{np};$ $\displaystyle\beta=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i} \right)\cdot\lambda_{\max}\left({\bm{L}}^{-1}{\bm{L}}_{i}\right);$ $\displaystyle\Lambda_{{\bm{W}},\mathcal{S}}=\lambda_{\max}\left({\mathbb{E}} \left[{\bm{S}}_{i}^{k}{\bm{W}}{\bm{S}}_{i}^{k}\right]-{\bm{W}}\right);$ $\displaystyle\lambda_{{\bm{W}}}=\lambda_{\max}^{-1}\left({\bm{W}}^{\frac{1}{2} }{\bm{L}}{\bm{W}}^{\frac{1}{2}}\right).$ $$

As a result, the above inequality can be written equivalently as

$$ $\alpha\beta\Lambda_{{\bm{W}},\mathcal{S}}\cdot\gamma^{2}+\gamma-\lambda_{{\bm{ W}}}\leq 0,$ $$

which yields the upper bound on $\gamma$

$$ $\displaystyle\gamma\leq\frac{\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{W}}, \mathcal{S}}\lambda_{{\bm{W}}}}-1}{2\alpha\beta\cdot\Lambda_{{\bm{W}},\mathcal {S}}}.$ $$

Since $\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{W}},\mathcal{S}}\lambda_{{\bm{W}}}}+1>0$, we can simplify the result as

$$ $\displaystyle\gamma\leq\frac{2\lambda_{{\bm{W}}}}{1+\sqrt{1+4\alpha\beta\cdot \Lambda_{{\bm{W}},\mathcal{S}}\lambda_{{\bm{W}}}}}.$ $$

#### D.4.2 Proof of Corollary 3

It is obvious that (7) directly follows from plugging ${\bm{W}}={\bm{L}}^{-1}$ into (6).
The optimal stepsize is obtained as the product of $\gamma$ and ${\bm{L}}^{-1}$.
The iteration complexity of MARINA, according to Gorbunov et al. (), is

$$ $K\geq K_{1}=\mathcal{O}\left(\frac{\Delta_{0}L}{\varepsilon^{2}}\left(1+\sqrt{ \frac{(1-p)\omega}{pn}}\right)\right).$ (36) $$

On the other hand,

$$ $\det({\bm{L}})^{\frac{1}{d}}\leq\lambda_{\max}\left({\bm{L}}\right)=L.$ (37) $$

In addition, using the inequality

$$ $\sqrt{1+4t}\leq 1+2\sqrt{t},$ (38) $$

which holds for any $t\geq 0$, we have the following bound

$$ $\displaystyle\frac{\left(1+\sqrt{1+4\alpha\beta\cdot\Lambda_{{\bm{L}}^{-1}, \mathcal{S}}}\right)}{2}\leq 1+\sqrt{\alpha\beta\cdot\Lambda_{{\bm{L}}^{-1}, \mathcal{S}}}.$ $$

Next we prove that

$$ $1+\sqrt{\alpha\beta\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{S}}}\leq 1+\sqrt{\frac {(1-p)}{pn}\cdot\omega},$ (39) $$

which is equivalent to proving

$$ $\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right)\lambda_{\max} \left({\bm{L}}_{i}{\bm{L}}^{-1}\right)\cdot\lambda_{\max}\left({\mathbb{E}} \left[{\bm{S}}_{i}^{k}{\bm{L}}^{-1}{\bm{S}}_{i}^{k}\right]-{\bm{L}}^{-1}\right )\leq\omega.$ $$

The left hand side can be upper bounded by,

$$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right) \lambda_{\max}\left({\bm{L}}^{-1}{\bm{L}}_{i}\right)\cdot\lambda_{\max}\left({ \bm{L}}^{-1}\right)\cdot\frac{\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{ i}^{k}{\bm{L}}^{-1}{\bm{S}}_{i}^{k}\right]-{\bm{L}}^{-1}\right)}{\lambda_{\max }\left({\bm{L}}^{-1}\right)}$ $\displaystyle\quad\leq\frac{\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i} ^{k}{\bm{L}}^{-1}{\bm{S}}_{i}^{k}\right]-{\bm{L}}^{-1}\right)}{\lambda_{\max} \left({\bm{L}}^{-1}\right)},$ $$

where the inequality is a consequence of Proposition 1.
We further bound the last term with

$$ $\displaystyle\frac{\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{ L}}^{-1}{\bm{S}}_{i}^{k}\right]-{\bm{L}}^{-1}\right)}{\lambda_{\max}\left({\bm {L}}^{-1}\right)}$ $\displaystyle=$ $\displaystyle\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}\cdot\frac{ {\bm{L}}^{-1}}{\lambda_{\max}({\bm{L}}^{-1})}\cdot{\bm{S}}_{i}^{k}\right]- \frac{{\bm{L}}^{-1}}{\lambda_{\max}\left({\bm{L}}^{-1}\right)}\right)$ $\displaystyle\leq$ $\displaystyle\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{S}}_{i }^{k}\right]-{\bm{I}}_{d}\right)=:\omega.$ $$

Here, the last inequality is due to the monotonicity of the mapping ${\bm{X}}\mapsto\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{X}}{
\bm{S}}_{i}^{k}\right]-{\bm{X}}\right)$ with ${\bm{X}}\in\mathbb{S}^{d}_{++}$, which can be shown as follows, let us pick any ${\bm{X}}_{1},{\bm{X}}_{2}\in\mathbb{S}^{d}_{++}$ and ${\bm{X}}_{1}\preceq{\bm{X}}_{2}$,

$$ $\displaystyle\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{X}}_{2}{\bm{S}}_{i}^ {k}\right]-{\bm{X}}_{2}\right)-\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{X} }_{1}{\bm{S}}_{i}^{k}\right]-{\bm{X}}_{1}\right)={\mathbb{E}}\left[{\bm{S}}_{i }^{k}\left({\bm{X}}_{2}-{\bm{X}}_{1}\right){\bm{S}}_{i}^{k}\right]-\left({\bm{ X}}_{2}-{\bm{X}}_{1}\right)\succeq{\bm{O}}_{d}.$ $$

The above inequality is due to the convexity of the mapping ${\bm{S}}_{i}^{k}\mapsto{\bm{S}}_{i}^{k}{\bm{X}}{\bm{S}}_{i}^{k}$. As a result, we have

$$ $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{X}}_{2}{\bm{S}}_{i} ^{k}\right]-{\bm{X}}_{2}\right)\geq\lambda_{\max}\left({\mathbb{E}}\left[{\bm{ S}}_{i}^{k}{\bm{X}}_{1}{\bm{S}}_{i}^{k}\right]-{\bm{X}}_{1}\right),$ $$

whenever ${\bm{X}}_{2}\succeq{\bm{X}}_{1}$. Due to the fact that

$$ $\frac{{\bm{L}}^{-1}}{\lambda_{\max}\left({\bm{L}}^{-1}\right)}\preceq{\bm{I}}_ {d},$ $$

we have

$$ $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}\cdot\frac{{\bm{L}}^{-1} }{\lambda_{\max}({\bm{L}}^{-1})}\cdot{\bm{S}}_{i}^{k}\right]-\frac{{\bm{L}}^{- 1}}{\lambda_{\max}\left({\bm{L}}^{-1}\right)}\right)\leq\lambda_{\max}\left({ \mathbb{E}}\left[{\bm{S}}_{i}^{k}\cdot{\bm{I}}_{d}\cdot{\bm{S}}_{i}^{k}\right] -{\bm{I}}_{d}\right)=\omega.$ $$

Combining (37) and (39), we know that the iteration complexity of det-MARINA is always better than that of MARINA.

#### D.4.3 Proof of Corollary 4

The number of bits sent in expectation is

$$ $\mathcal{O}(d+K(pd+(1-p)\zeta_{\mathcal{S}}))=\mathcal{O}((Kp+1)d+(1-p)K\zeta_ {\mathcal{S}}).$ $$

The special case where we choose $p=\zeta_{\mathcal{S}}/d$ indicates that

$$ $\alpha=\frac{1-p}{np}=\frac{1}{n}\left(\frac{d}{\zeta_{\mathcal{S}}}-1\right).$ $$

In order to reach an error of $\varepsilon^{2}$, we need

$$ $\displaystyle K=\mathcal{O}\left(\frac{\Delta_{0}\cdot\det({\bm{L}})^{\frac{1} {d}}}{\varepsilon^{2}}\cdot\left(1+\sqrt{1+\frac{4\beta}{n}\left(\frac{d}{ \zeta_{\mathcal{S}}}-1\right)\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{S}}}\right) \right),$ $$

which is the iteration complexity.
Applying once again (38) and using the fact that $p=\zeta_{\mathcal{S}}/d$, the communication complexity in this case is given by

$$ $\displaystyle\mathcal{O}\left(d+\frac{\Delta_{0}\cdot\det({\bm{L}})^{\frac{1}{ d}}}{\varepsilon^{2}}\cdot\left(1+\sqrt{1+\frac{4\beta}{n}\left(\frac{d}{\zeta _{\mathcal{S}}}-1\right)\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{S}}}\right)\cdot \left(pd+(1-p)\zeta_{\mathcal{S}}\right)\right)$ $\displaystyle\quad\leq\mathcal{O}\left(d+\frac{2\Delta_{0}\cdot\det({\bm{L}})^ {\frac{1}{d}}}{\varepsilon^{2}}\cdot\left(1+\sqrt{\frac{\beta}{n}\left(\frac{d }{\zeta_{\mathcal{S}}}-1\right)\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{S}}}\right )\cdot\left(pd+(1-p)\zeta_{\mathcal{S}}\right)\right)$ $\displaystyle\quad\leq\mathcal{O}\left(d+\frac{4\Delta_{0}\cdot\det({\bm{L}})^ {\frac{1}{d}}}{\varepsilon^{2}}\cdot\left(\zeta_{\mathcal{S}}+\sqrt{\frac{ \beta\cdot\Lambda_{{\bm{L}}^{-1},\mathcal{S}}}{n}\cdot\zeta_{\mathcal{S}}(d- \zeta_{\mathcal{S}})}\right)\right).$ $$

Ignoring the coefficient we get

$$ $\mathcal{O}\left(d+\frac{\Delta_{0}\cdot\det({\bm{L}})^{\frac{1}{d}}}{ \varepsilon^{2}}\cdot\left(\zeta_{\mathcal{S}}+\sqrt{\frac{\beta\cdot\Lambda_{ {\bm{L}}^{-1},\mathcal{S}}}{n}\cdot\zeta_{\mathcal{S}}(d-\zeta_{\mathcal{S}})} \right)\right).$ $$

#### D.4.4 Proof of Corollary 8

Applying Corollary 1, notice that in this case

$$ $\lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}=\lambda_{\max}^{-1} \left(\operatorname{diag}^{-\frac{1}{2}}\left({\bm{L}}\right){\bm{L}} \operatorname{diag}^{-\frac{1}{2}}\left({\bm{L}}\right)\right)=1,$ $$

we obtain ${\bm{D}}^{*}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$.
The iteration complexity is given by

$$ $\mathcal{O}\left(\frac{\det\left(\operatorname{diag}({\bm{L}})\right)^{\frac{1 }{d}}\cdot\Delta_{0}}{\varepsilon^{2}}\cdot\left(\frac{1+\sqrt{1+4\alpha\beta \Lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right),\mathcal{S}}}}{2}\right )\right).$ $$

We now compare it to the iteration complexity of MARINA, which is given in (36).
We know that each diagonal element ${\bm{L}}_{jj}$ satisfies ${\bm{L}}_{jj}\leq\lambda_{\max}\left({\bm{L}}\right)=L$ for $j=1,\ldots,d$.
As a result,

$$ $\det\left(\operatorname{diag}({\bm{L}})\right)^{\frac{1}{d}}\leq L.$ (40) $$

From (38), we deduce

$$ $\displaystyle\frac{1+\sqrt{1+4\alpha\beta\cdot\Lambda_{\operatorname{diag}^{-1 }({\bm{L}}),\mathcal{S}}}}{2}\leq 1+\sqrt{\alpha\beta\cdot\Lambda_{ \operatorname{diag}^{-1}({\bm{L}}),\mathcal{S}}}.$ $$

Now, let us prove the below inequality

$$ $1+\sqrt{\alpha\beta\cdot\Lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right) ,\mathcal{S}}}\leq 1+\sqrt{\frac{(1-p)}{pn}\cdot\omega}.$ (41) $$

The latter is equivalent to

$$ $\beta\cdot\Lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right),\mathcal{S}} \leq\omega.$ $$

Plugging in the definition of $\beta$, $\omega$ and $\Lambda_{\operatorname{diag}^{-1}\left({\bm{L}}\right),\mathcal{S}}$ and using the relation given in Proposition 1, we obtain,

$$ $\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}\frac{\operatorname{diag }^{-1}\left({\bm{L}}\right)}{\lambda_{\max}\left({\bm{L}}^{-1}\right)}{\bm{S}} _{i}^{k}-\frac{\operatorname{diag}^{-1}\left({\bm{L}}\right)}{\lambda_{\max} \left({\bm{L}}^{-1}\right)}\right]\right)\leq\lambda_{\max}\left({\mathbb{E}} \left[{\bm{S}}_{i}^{k}{\bm{I}}_{d}{\bm{S}}_{i}^{k}\right]-{\bm{I}}_{d}\right).$ $$

Thus, it is enough to prove that

$$ $\frac{\operatorname{diag}^{-1}\left({\bm{L}}\right)}{\lambda_{\max}\left({\bm{ L}}^{-1}\right)}\preceq{\bm{I}}_{d}.$ $$

We can further simplify the above inequality as

$$ $\lambda_{\min}\left({\bm{L}}\right)\leq\lambda_{\min}\left(\operatorname{diag} ({\bm{L}})\right),$ $$

which is always true for any ${\bm{L}}\in\mathbb{S}^{d}_{++}$.
Combining (40) and (41) we conclude the proof.

#### D.4.5 Proof of Corollary 9

Using the explicit formula for the optimal stepsize ${\bm{D}}^{*}_{{\bm{I}}_{d}}$,
we deduce the following iteration complexity for

$$ $\mathcal{O}\left(\frac{\lambda_{\max}\left({\bm{L}}\right)\Delta_{0}}{ \varepsilon^{2}}\cdot\left(\frac{1+\sqrt{1+4\alpha\beta\frac{\omega}{\lambda_{ \max}\left({\bm{L}}\right)}}}{2}\right)\right).$ (42) $$

Recall that $\lambda_{\max}\left({\bm{L}}\right)=L$, we obtain using (38) that

$$ $\displaystyle\frac{1+\sqrt{1+4\alpha\beta\frac{\omega}{\lambda_{\max}\left({ \bm{L}}\right)}}}{2}$ $\displaystyle\leq 1+\sqrt{\alpha\beta\frac{\omega}{\lambda_{\max}\left({\bm{L} }\right)}}.$ $$

The comparison of two iteration complexities, given in (42) and (36) reduces to

$$ $1+\sqrt{\alpha\beta\frac{\omega}{\lambda_{\max}\left({\bm{L}}\right)}}\leq 1+ \sqrt{\frac{1-p}{np}\omega}.$ $$

This is equivalent to

$$ $\beta\cdot\frac{1}{\lambda_{\max}\left({\bm{L}}\right)}\leq 1.$ $$

Utilizing Proposition 1, the above inequality can be rewritten as

$$ $\frac{1}{\lambda_{\max}\left({\bm{L}}^{-1}\right)\cdot\lambda_{\max}\left({\bm {L}}\right)}\leq 1,$ $$

which is exactly

$$ $\lambda_{\min}\left({\bm{L}}\right)\leq\lambda_{\max}\left({\bm{L}}\right).$ $$

## Appendix E Analysis of det-DASHA

### E.1 Proof of Theorem 2

Using Lemma 1 and taking expectations, we are able to obtain

$$ $\displaystyle{\mathbb{E}}\left[f(x^{k+1})\right]$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})\right]-\frac{1}{2}{\mathbb{E} }\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}} \right]+\frac{1}{2}{\mathbb{E}}\left[\left\|g^{k}-\nabla f(x^{k})\right\|_{{ \bm{D}}}^{2}\right]$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})\right]-\frac{1}{2}{\mathbb{E} }\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]$ $\displaystyle\qquad+{\mathbb{E}}\left[\frac{1}{2}\left\|g^{k}-h^{k}+h^{k}- \nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})\right]-\frac{1}{2}{\mathbb{E} }\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]$ $\displaystyle\qquad+{\mathbb{E}}\left[\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}} }+\left\|h^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right],$ (43) $$

where the last step is due to the convexity of the norm.
Using Lemma 4, we obtain

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\leq\frac{2\omega_{{\bm{D}}}\cdot\lambda_{\max}\left({\bm{D}} \right)}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right){\mathbb{E }}\left[\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}\right]$ $\displaystyle\qquad+\frac{2a^{2}\omega_{{\bm{D}}}}{n^{2}}\sum_{i=1}^{n}{ \mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]+(1-a )^{2}{\mathbb{E}}\left[\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}\right].$ (44) $$

Using Lemma 5, we get

$$ $\displaystyle{\mathbb{E}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1}\right\|^{2}_{{ \bm{D}}}\right]$ $\displaystyle\leq\left(2a^{2}\omega_{{\bm{D}}}+(1-a)^{2}\right)\cdot{\mathbb{E }}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+2\omega_{{\bm{D}}}\cdot\lambda_{\max}\left({\bm{D}}\right) \cdot\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot{\mathbb{E}}\left[\left\|h_{i }^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}\right].$ (45) $$

Now let us fix $\kappa\in[0,+\infty)$, $\eta\in[0,+\infty)$ which we will determine later, and construct the following Lyapunov function $\Phi_{k}$

$$ $\Phi_{k}={\mathbb{E}}\left[f(x^{k})-f^{\star}\right]+\kappa\cdot{\mathbb{E}} \left[\left\|g^{k}-h^{k}\right\|_{{\bm{D}}}^{2}\right]+\eta\cdot{\mathbb{E}} \left[\frac{1}{n}\sum_{i=1}^{n}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D} }}\right].$ (46) $$

Combining (E.1), (E.1) and (E.1), we get

$$ $\displaystyle\Phi_{k+1}$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})-f^{\star}-\frac{1}{2}\left\| \nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+{\mathbb{E}}\left[-\frac{1}{2}\left\|x^{k+1}-x^{k}\right\| ^{2}_{{\bm{D}}^{-1}-{\bm{L}}}+\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}+\left \|h^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+\kappa(1-a)^{2}{\mathbb{E}}\left[\left\|g^{k}-h^{k}\right \|^{2}_{{\bm{D}}}\right]+\frac{2\kappa\cdot\omega_{{\bm{D}}}\lambda_{\max} \left({\bm{D}}\right)}{n}\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({ \bm{L}}_{i}\right){\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{ {\bm{L}}_{i}^{-1}}\right]$ $\displaystyle\qquad+\frac{2a^{2}\omega_{{\bm{D}}}\cdot\kappa}{n}\cdot\frac{1}{ n}\sum_{i=1}^{n}{\mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm {D}}}\right]+\eta\left(2a^{2}\omega_{{\bm{D}}}+(1-a)^{2}\right)\cdot\frac{1}{n }\sum_{i=1}^{n}{\mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{ D}}}\right]$ $\displaystyle\qquad+2\eta\cdot\omega_{{\bm{D}}}\cdot\lambda_{\max}\left({\bm{D }}\right)\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i}\right) \cdot{\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^ {-1}}\right].$ $$

Rearranging terms, and notice that $\left\|h^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}=0$,

$$ $\displaystyle\Phi_{k+1}$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})-f^{\star}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad-\frac{1}{2}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^ {2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]+\left(1+\kappa(1-a)^{2}\right){\mathbb{E}} \left[\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+\left(\frac{2a^{2}\omega_{{\bm{D}}}\cdot\kappa}{n}+\eta \left(2a^{2}\omega_{{\bm{D}}}+(1-a)^{2}\right)\right)\cdot\frac{1}{n}\sum_{i=1 }^{n}{\mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+\left(\frac{2\kappa\cdot\omega_{{\bm{D}}}\lambda_{\max} \left({\bm{D}}\right)}{n}+2\eta\cdot\omega_{{\bm{D}}}\cdot\lambda_{\max}\left( {\bm{D}}\right)\right)\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L }}_{i}\right)\cdot{\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{ {\bm{L}}_{i}^{-1}}\right].$ $$

In order to proceed, we consider the choice of $\kappa$ and $\eta$, for $\kappa$,

$$ $1+\kappa(1-a)^{2}\leq\kappa.$ (47) $$

It is then clear that the choice of $\kappa=\frac{1}{a}$ satisfies the condition.
On the other hand, we look at the terms involving ${\mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]$, we can rewrite as

$$ $\displaystyle T_{1}$ $\displaystyle:=\left(\frac{2a^{2}\omega_{{\bm{D}}}\cdot\kappa}{n}+\eta\left(2a ^{2}\omega_{{\bm{D}}}+(1-a)^{2}\right)\right)\cdot\frac{1}{n}\sum_{i=1}^{n}{ \mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right].$ $$

Picking $\kappa=\frac{1}{a}$ and $a=\frac{1}{2\omega_{{\bm{D}}}+1}$, the $T_{1}$ can be simplified as

$$ $\displaystyle T_{1}=\left(\frac{2\omega_{{\bm{D}}}}{n\cdot(2\omega_{{\bm{D}}}+ 1)}+\eta\cdot\frac{4\omega_{{\bm{D}}}^{2}+2\omega_{{\bm{D}}}}{\left(2\omega_{{ \bm{D}}}+1\right)^{2}}\right)\cdot\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[ \left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right].$ $$

We pick $\eta$ so that it satisfies

$$ $\displaystyle\left(\frac{2\omega_{{\bm{D}}}}{n\cdot(2\omega_{{\bm{D}}}+1)}+ \eta\cdot\frac{4\omega_{{\bm{D}}}^{2}+2\omega_{{\bm{D}}}}{\left(2\omega_{{\bm{ D}}}+1\right)^{2}}\right)\leq\eta.$ (48) $$

Taking $\eta=\frac{2\omega_{{\bm{D}}}}{n}$, which is the minimum value satisfying (48), we conclude that

$$ $\displaystyle T_{1}\leq\eta\cdot\frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[ \left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right].$ (49) $$

Combining (47) and (49), we are able to conclude that

$$ $\displaystyle\Phi_{k+1}$ $\displaystyle\quad\leq{\mathbb{E}}\left[f(x^{k})-f^{\star}\right]+\kappa\cdot{ \mathbb{E}}\left[\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}\right]+\eta\cdot \frac{1}{n}\sum_{i=1}^{n}{\mathbb{E}}\left[\left\|g_{i}^{k}-h_{i}^{k}\right\|^ {2}_{{\bm{D}}}\right]$ $\displaystyle\qquad-\frac{1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right \|^{2}_{{\bm{D}}}\right]-\frac{1}{2}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k} \right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]$ $\displaystyle\qquad+\left(\frac{2\kappa\cdot\omega_{{\bm{D}}}\lambda_{\max} \left({\bm{D}}\right)}{n}+2\eta\cdot\omega_{{\bm{D}}}\cdot\lambda_{\max}\left( {\bm{D}}\right)\right)\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L }}_{i}\right)\cdot{\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{ {\bm{L}}_{i}^{-1}}\right].$ $$

Using the definition of $\Phi_{k}$ and Lemma 3, we obtain

$$ $\displaystyle\Phi_{k+1}$ $\displaystyle\leq\Phi_{k}-\frac{1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k}) \right\|^{2}_{{\bm{D}}}\right]-\frac{1}{2}{\mathbb{E}}\left[\left\|x^{k+1}-x^{ k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}\right]$ $\displaystyle\qquad\left(\frac{2\kappa\cdot\omega_{{\bm{D}}}\lambda_{\max} \left({\bm{D}}\right)}{n}+2\eta\cdot\omega_{{\bm{D}}}\cdot\lambda_{\max}\left( {\bm{D}}\right)\right)\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L }}_{i}\right)\cdot{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}} _{i}}\right]$ $\displaystyle=\Phi_{k}-\frac{1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k}) \right\|^{2}_{{\bm{D}}}\right]+{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^ {2}_{{\bm{N}}}\right],$ $$

where ${\bm{N}}\in\mathbb{S}^{d}$ is defined as

$$ $\displaystyle{\bm{N}}$ $\displaystyle:=\left(\frac{2\kappa\cdot\omega_{{\bm{D}}}\lambda_{\max}\left({ \bm{D}}\right)}{n}+2\eta\cdot\omega_{{\bm{D}}}\cdot\lambda_{\max}\left({\bm{D} }\right)\right)\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{i} \right)\cdot{\bm{L}}_{i}-\frac{1}{2}{\bm{D}}^{-1}+\frac{1}{2}{\bm{L}}.$ $$

We require ${\bm{N}}\preceq{\bm{O}}_{d}$, which leads to the condition on ${\bm{D}}$:

$$ ${\bm{D}}^{-1}-{\bm{L}}-\frac{4\lambda_{\max}\left({\bm{D}}\right)\cdot\omega_{ {\bm{D}}}\cdot\left(4\omega_{{\bm{D}}}+1\right)}{n}\cdot\frac{1}{n}\sum_{i=1}^ {n}\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot{\bm{L}}_{i}\succeq{\bm{O}}_{d}.$ $$

Given the above condition is satisfied, we have the recurrence

$$ $\frac{1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}} \right]\leq\Phi_{k}-\Phi_{k+1}$ $$

Summing up for $k=0\ldots K-1$, we obtain

$$ $\displaystyle\sum_{k=0}^{K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^ {2}_{{\bm{D}}}\right]\leq 2(\Phi_{0}-\Phi_{k}).$ (50) $$

Notice that we also have

$$ $\displaystyle\Phi_{0}$ $\displaystyle=f(x^{0})-f^{\star}+\left(2\omega_{{\bm{D}}}+1\right)\left\|g^{0} -h^{0}\right\|^{2}_{{\bm{D}}}+frac{2\omega_{{\bm{D}}}}{n}\cdot\frac{1}{n}\sum_ {i=1}^{n}\left\|g_{i}^{0}-h_{i}^{0}\right\|^{2}$ $\displaystyle=f(x^{0})-f^{\star},$ $$

We divide both sides of (50) by $K$, and perform determinant normalization,

$$ $\frac{1}{K}\sum_{k=0}^{K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2 }_{\frac{{\bm{D}}}{\det({\bm{D}})^{1/d}}}\right]\leq\frac{2(f(x^{0})-f^{\star} )}{\det({\bm{D}})^{1/d}\cdot K}.$ $$

This is to say

$$ ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{ \det(D)^{1/d}}}\right]\leq\frac{2(f(x^{0})-f^{\star})}{\det({\bm{D}})^{1/d} \cdot K},$ $$

where $\tilde{x}^{K}$ is chosen uniformly randomly from the first K iterates of the algorithm.

### E.2 Proofs of the corollaries

#### E.2.1 Proof of Corollary 2

Plug ${\bm{D}}=\gamma_{{\bm{W}}}\cdot{\bm{W}}$ into the stepsize condition in Theorem 2, we obtain

$$ $\displaystyle\frac{1}{\gamma_{{\bm{W}}}}\cdot{\bm{W}}^{-1}-{\bm{L}}-\frac{4 \gamma_{{\bm{W}}}\cdot\lambda_{\max}\left({\bm{W}}\right)\cdot\omega_{{\bm{W}} }\left(4\omega_{{\bm{W}}}+1\right)}{n}\cdot\frac{1}{n}\sum_{i=1}^{n}\lambda_{ \max}\left({\bm{L}}_{i}\right)\cdot{\bm{L}}_{i}\succeq{\bm{O}}_{d}.$ $$

We then simplify the above condition as

$$ $\displaystyle\frac{1}{\gamma_{{\bm{W}}}}\cdot{\bm{L}}^{-\frac{1}{2}}{\bm{W}}^{ -1}{\bm{L}}^{-\frac{1}{2}}$ $\displaystyle\quad\succeq{\bm{I}}_{d}+\frac{4\gamma_{{\bm{W}}}\cdot\lambda_{ \max}\left({\bm{W}}\right)\cdot\omega_{{\bm{W}}}\left(4\omega_{{\bm{W}}}+1 \right)}{n}\cdot{\bm{L}}^{-\frac{1}{2}}\left(\frac{1}{n}\sum_{i=1}^{n}\lambda_ {\max}\left({\bm{L}}_{i}\right)\cdot{\bm{L}}_{i}\right){\bm{L}}^{-\frac{1}{2}}.$ $$

Using Proposition 3, we have

$$ $\frac{1}{\gamma_{{\bm{W}}}}\cdot{\bm{L}}^{-\frac{1}{2}}{\bm{W}}^{-1}{\bm{L}}^{ -\frac{1}{2}}-\frac{4\gamma_{{\bm{W}}}\cdot\lambda_{\max}\left({\bm{W}}\right) \cdot\omega_{{\bm{W}}}\left(4\omega_{{\bm{W}}}+1\right)}{n}\cdot\lambda_{\min} \left({\bm{L}}\right)\cdot{\bm{I}}_{d}\succeq{\bm{I}}_{d}.$ $$

Taking the minimum eigenvalue of both sides, we obtain that,

$$ $\displaystyle\frac{1}{\gamma_{{\bm{W}}}}\cdot\lambda_{\min}\left({\bm{L}}^{- \frac{1}{2}}{\bm{W}}^{-1}{\bm{L}}^{-\frac{1}{2}}\right)-\frac{4\gamma_{{\bm{W} }}\cdot\lambda_{\max}\left({\bm{W}}\right)\cdot\omega_{{\bm{W}}}\left(4\omega_ {{\bm{W}}}+1\right)}{n}\cdot\lambda_{\min}\left({\bm{L}}\right)\geq 1,$ $$

If we denote $C_{{\bm{W}}}:=\frac{\lambda_{\max}\left({\bm{W}}\right)\cdot\omega_{{\bm{W}}}
\left(4\omega_{{\bm{W}}}+1\right)}{n}>0$, and $\lambda_{{\bm{W}}}:=\lambda_{\max}^{-1}\left({\bm{L}}^{\frac{1}{2}}{\bm{W}}{
\bm{L}}^{\frac{1}{2}}\right)$, we can write

$$ $\displaystyle 4\cdot C_{{\bm{W}}}\cdot\lambda_{\min}\left({\bm{L}}\right)\cdot \gamma_{{\bm{W}}}^{2}+\gamma_{{\bm{W}}}-\lambda_{{\bm{W}}}\leq 0.$ $$

The solution is given by

$$ $\displaystyle\gamma_{{\bm{W}}}\leq\frac{2\lambda_{{\bm{W}}}}{1+\sqrt{1+16C_{{ \bm{W}}}\lambda_{\min}\left({\bm{L}}\right)\cdot\lambda_{{\bm{W}}}}}.$ $$

#### E.2.2 Proof of Corollary 5

The best scaling factor in this case is given as, according to Corollary 2,

$$ $\gamma_{{\bm{L}}^{-1}}=\frac{2}{1+\sqrt{1+16C_{{\bm{L}}^{-1}}\cdot\lambda_{ \min}\left({\bm{L}}\right)}}.$ $$

In order to reach a $\varepsilon^{2}$ stationary point, we need

$$ $K\geq\frac{\det({\bm{L}})^{\frac{1}{d}}\left(f(x^{0})-f^{\star}\right)}{ \varepsilon^{2}}\cdot\left(1+\sqrt{1+16C_{{\bm{L}}^{-1}}\cdot\lambda_{\min} \left({\bm{L}}\right)}\right).$ $$

#### E.2.3 Proof of Corollary 6

The iteration complexity of det-DASHA is given by, according to, Corollary 5,

$$ $\mathcal{O}\left(\frac{f(x^{0})-f^{\star}}{\epsilon^{2}}\cdot\left(1+\sqrt{1+1 6C_{{\bm{L}}^{-1}}\cdot\lambda_{\min}\left({\bm{L}}\right)}\right)\cdot\det({ \bm{L}})^{\frac{1}{d}}\right).$ $$

Using the inequality $\sqrt{1+t}\leq 1+\sqrt{t}$ for $t>0$, and leaving out the coefficients, we obtain

$$ $\mathcal{O}\left(\frac{f(x^{0})-f^{\star}}{\epsilon^{2}}\cdot\left(1+\sqrt{C_{ {\bm{L}}^{-1}}\cdot\lambda_{\min}\left({\bm{L}}\right)}\right)\cdot\det({\bm{L }})^{\frac{1}{d}}\right).$ $$

Notice that

$$ $C_{{\bm{L}}^{-1}}\cdot\lambda_{\min}\left({\bm{L}}\right)=\lambda_{\max}\left( {\bm{L}}^{-1}\right)\cdot\frac{\omega_{{\bm{L}}^{-1}}\left(4\omega_{{\bm{L}}^{ -1}}+1\right)}{n}\cdot\lambda_{\min}\left({\bm{L}}\right)=\frac{\omega_{{\bm{L }}^{-1}}\left(4\omega_{{\bm{L}}^{-1}+1}\right)}{n}.$ $$

As a result, the iteration complexity can be further simplified as

$$ $\mathcal{O}\left(\frac{f(x^{0})-f^{*}}{\epsilon^{2}}\cdot\left(1+\frac{\omega_ {{\bm{L}}^{-1}}}{\sqrt{n}}\right)\cdot\det({\bm{L}})^{\frac{1}{d}}\right).$ $$

The iteration complexity of DASHA is, according to Tyurin & Richtárik (, Corollary 6.2)

$$ $\mathcal{O}\left(\frac{1}{\epsilon^{2}}\cdot\left(f(x^{0})-f^{\star}\right) \left(L+\frac{\omega}{\sqrt{n}}\widehat{L}\right)\right),$ $$

where $\widehat{L}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}}$. Since $\det({\bm{L}})^{\frac{1}{d}}\leq\lambda_{\max}\left({\bm{L}}\right)=L$, and $L\leq\widehat{L}$, it is easy to see that compared to DASHA, det-DASHA has a better iteration complexity when the momentum is the same. Notice that those two algorithms use the same sketch, thus, it also indicates that the communication complexity of the two algorithms are the same.

#### E.2.4 Proof of Corollary 7

The iteration complexity of det-MARINA is given by

$$ $\mathcal{O}\left(\frac{f(x^{0})-f^{\star}}{\epsilon^{2}}\cdot\det({\bm{L}})^{ \frac{1}{d}}\cdot\left(1+\sqrt{\alpha\beta\Lambda_{{\bm{L}}^{-1},\mathcal{S}}} \right)\right),$ $$

after removing logarithmic factors.
Plugging in the definitions we obtain in the case of $\omega_{{\bm{L}}^{-1}}+1=\frac{1}{p}$, we have

$$ $\mathcal{O}\left(\frac{f(x^{0})-f^{\star}}{\epsilon^{2}}\cdot\det({\bm{L}})^{ \frac{1}{d}}\cdot\left(1+\frac{\omega_{{\bm{L}}^{-1}}}{n}\right)\right).$ $$

From the proof of Corollary 6, we know that the iteration complexity of det-DASHA is

$$ $\mathcal{O}\left(\frac{1}{\epsilon^{2}}\cdot\left(f(x^{0})-f^{\star}\right) \left(L+\frac{\omega}{\sqrt{n}}\widehat{L}\right)\right).$ $$

It is easy to see that in this case the two algorithms have the same iteration complexity asymptotically. Notice that the communication complexity is the product of bytes sent per iteration and the number of iterations. det-DASHA clearly sends less bytes per iteration because it always sent the compressed gradient differences, which means that it has a better communication complexity than det-MARINA.

## Appendix F Distributed det-CGD

This section is a brief summary of the distributed det-CGD algorithm and its theoretical analysis.
The details can be found in (Li et al., ).
The algorithm follows the standard FL paradigm.
See the pseudocode in Algorithm 3.

Figure: Algorithm 3 Distributed det-CGD

Below is the main convergence result for the algorithm.

###### Theorem 3 .

Suppose that $f$ is ${\bm{L}}$-smooth.
Under the Assumptions 1,2, if the stepsize satisfies

$$ ${\bm{D}}{\bm{L}}{\bm{D}}\preceq{\bm{D}},$ (51) $$

then the following convergence bound is true for the iteration of Algorithm 3:

$$ $\min_{0\leq k\leq K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|_{\frac{ {\bm{D}}}{\det({\bm{D}})^{1/d}}}^{2}\right]\leq\frac{2(1+\frac{\lambda_{{\bm{D }}}}{n})^{K}\left(f(x^{0})-f^{\star}\right)}{\det({\bm{D}})^{1/d}\,K}+\frac{2 \lambda_{{\bm{D}}}\Delta^{\star}}{\det({\bm{D}})^{1/d}\,n},$ (52) $$

where $\Delta^{\star}:=f^{\star}-\frac{1}{n}\sum_{i=1}^{n}f_{i}^{\star}$ and

$$ $\displaystyle\lambda_{{\bm{D}}}$ $\displaystyle:=\max_{i}\left\{\lambda_{\max}\left({\mathbb{E}}\left[{\bm{L}}_{ i}^{\frac{1}{2}}\left({\bm{S}}_{i}^{k}-{\bm{I}}_{d}\right){\bm{D}}{\bm{L}}{\bm {D}}\left({\bm{S}}_{i}^{k}-{\bm{I}}_{d}\right){\bm{L}}_{i}^{\frac{1}{2}}\right ]\right)\right\}.$ $$

###### Remark 8 .

On the right hand side of (52) we observe that increasing $K$ will only reduce the first term, that corresponds to the convergence error.
Whereas, the second term, which does not depend on $K$, will remain constant, if the other parameters of the algorithm are fixed.
This testifies to the neighborhood phenomenon which we discussed in Section 2.

###### Remark 9 .

If the stepsize satisfies the below conditions,

$$ ${\bm{D}}{\bm{L}}{\bm{D}}\preceq{\bm{D}},\quad\lambda_{{\bm{D}}}\leq\min\left\{ \frac{n}{K},\frac{n\varepsilon^{2}}{4\Delta^{\star}}\det({\bm{D}})^{1/d}\right \},\quad K\geq\frac{12(f(x^{0})-f^{\star})}{\det({\bm{D}})^{1/d}\,\varepsilon^ {2}},$ (53) $$

then we obtain $\varepsilon$-stationary point.

One can see that in the convergence guarantee of det-CGD in the distributed case, the result (52) is not variance-reduced.
Because of this limitation, in order to reach a $\varepsilon$ stationary point, the stepsize condition in (53) is restrictive.

###### Theorem 3 .

###### Remark 8 .

###### Remark 9 .

## Appendix G Extension of det-CGD2 in MARINA form

In this section we want to extend det-CGD2 into its variance reduced counterpart in MARINA form.

### G.1 Extension of det-CGD2 to its variance reduced counterpart

Figure: Algorithm 4 det-CGD2-VR

We call det-MARINA as the extension of det-CGD1, and Algorithm 4 as the extension of det-CGD2 due to the difference in the order of applying sketches and stepsize matrices.
The key difference between det-CGD1 and det-CGD2 is that in det-CGD1 the gradient is sketched first and then multiplied by the stepsize, while for det-CGD2, the gradient is multiplied by the stepsize first after which the product is sketched. The convergence for Algorithm 4 can be proved in a similar manner as Theorem 1.

###### Theorem 4 .

Let Assumptions 1 and 2 hold, with the gradient of $f$ being ${\bm{L}}$-Lipschitz. If the stepsize matrix ${\bm{D}}\in\mathbb{S}^{d}_{++}$ satisfies

$$ ${\bm{D}}^{-1}\succeq\left(\frac{(1-p)\cdot R^{\prime}({\bm{D}},\mathcal{S})}{ np}+1\right){\bm{L}},$ $$

where

$$ $R^{\prime}({\bm{D}},\mathcal{S})=\frac{1}{n}\sum_{i=1}^{n}\lambda_{\max}\left( {\bm{D}}{\mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{\bm{T}}_{i}^{k}\right] {\bm{D}}{\bm{L}}_{i}^{\frac{1}{2}}-{\bm{L}}_{i}^{\frac{1}{2}}{\bm{D}}\right) \cdot\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot\lambda_{\max}\left({\bm{L}}^ {-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right).$ $$

Then after $K$ iterations of Algorithm 4, we have

$$ ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{ \det({\bm{D}})^{1/d}}}\right]\leq\frac{2\left(f(x^{0})-f^{\star}\right)}{\det( {\bm{D}})^{1/d}\cdot K}.$ $$

This is to say that in order to reach a $\varepsilon$-stationary point, we require

$$ $K\geq\frac{2(f(x^{0})-f^{\star})}{\det({\bm{D}})^{1/d}\cdot\varepsilon^{2}}.$ $$

If we look at the scalar case where ${\bm{D}}=\gamma\cdot{\bm{I}}_{d}$, ${\bm{L}}_{i}=L_{i}\cdot{\bm{I}}_{d}$ and ${\bm{L}}=L\cdot I_{d}$, then the condition in Theorem 4 reduces to

$$ $\frac{(1-p)\omega L^{2}}{np}+L-\frac{1}{\gamma}\leq 0.$ (54) $$

Notice that here $\omega=\lambda_{\max}\left({\mathbb{E}}\left[\left({\bm{T}}_{i}^{k}\right)^{2}
\right]\right)-1$, and we have $L^{2}=\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}$, which is due to the relation given in Proposition 5.
This condition coincides with the condition for convergence of MARINA.
One may also check that, the update rule in Algorithm 4, is the same as MARINA in the scalar case.
However, the condition given in Theorem 4 is not simpler than Theorem 1, contrary to the single-node case.
We emphasize that Algorithm 4 is not suitable for the federated learning setting where the clients have limited resources.
In order to perform the update, each client is required to store the stepsize matrix ${\bm{D}}$ which is of size $d\times d$.
In the over-parameterized regime, the dataset size is $m\times d$ where $m$ is the number of data samples, and we have $d>m$.
This means that the stepsize matrix each client needs to store is even larger than the dataset itself, which is unacceptable given the limited resources each client has.

###### Theorem 4 .

### G.2 Analysis of Algorithm 4

We first present two lemmas which are necessary for the proofs of Theorem 4.

###### Lemma 6 .

Assume that function $f$ is ${\bm{L}}$-smooth, and $x^{k+1}=x^{k}-g^{k}$, and matrix ${\bm{D}}\in\mathbb{S}^{d}_{++}$. Then we will have

$$ $f(x^{k+1})\leq f(x^{k})-\frac{1}{2}\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D} }}+\frac{1}{2}\left\|{\bm{D}}\cdot\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{D}}^ {-1}}-\frac{1}{2}\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}.$ (55) $$

This lemma is formulated in a different way from Lemma 1 on purpose.

###### Lemma 7 .

For any sketch matrix ${\bm{T}}\in\mathbb{S}_{+}^{d}$, vector $t\in\mathbb{R}^{d}$, matrix ${\bm{D}}\in\mathbb{S}^{d}_{++}$ and matrix ${\bm{L}}\in\mathbb{S}^{d}_{++}$, we have

$$ ${\mathbb{E}}\left[\left\|{\bm{T}}{\bm{D}}t-{\bm{D}}t\right\|^{2}_{{\bm{D}}^{-1 }}\right]\leq\lambda_{\max}\left({\bm{L}}^{\frac{1}{2}}{\bm{D}}{\mathbb{E}} \left[{\bm{T}}{\bm{D}}^{-1}{\bm{T}}\right]{\bm{D}}{\bm{L}}^{\frac{1}{2}}-{\bm{ L}}^{\frac{1}{2}}{\bm{D}}{\bm{L}}^{\frac{1}{2}}\right)\left\|t\right\|^{2}_{{ \bm{L}}^{-1}}.$ (56) $$

###### Lemma 6 .

###### Lemma 7 .

### G.3 Proof of Theorem 4

We start with Lemma 6,

$$ $\displaystyle{\mathbb{E}}\left[f(x^{k+1})\right]$ $\displaystyle\leq{\mathbb{E}}\left[f(x^{k})\right]-{\mathbb{E}}\left[\frac{1}{ 2}\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad+{\mathbb{E}}\left[\frac{1}{2}\left\|{\bm{D}}\cdot\nabla f(x ^{k})-g^{k}\right\|^{2}_{{\bm{D}}^{-1}}\right]-{\mathbb{E}}\left[\frac{1}{2} \left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}\right].$ (57) $$

Now we do the same as Theorem 1 and look at the term ${\mathbb{E}}\left[\left\|{\bm{D}}\cdot\nabla f(x^{k+1})-g^{k+1}\right\|^{2}_{{
\bm{D}}^{-1}}\right]$. Recall that $g^{k}$ here is given by

$$ $g^{k+1}=\begin{cases}{\bm{D}}\cdot\nabla f(x^{k+1})&\text{ with probability $p $ }\\ g^{k}+\frac{1}{n}\sum_{i=1}^{n}{\bm{T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k +1})-\nabla f_{i}(x^{k})\right)&\text{ with probability $1-p$ }.\end{cases}$ $$

As a result, we have

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right \|^{2}_{{\bm{D}}^{-1}}\mid x^{k+1},x^{k}\right]$ $\displaystyle={\mathbb{E}}\left[{\mathbb{E}}\left[\left\|g^{k+1}-{\bm{D}} \nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}^{-1}}\mid x^{k+1},x^{k},c_{k}\right]\right]$ $\displaystyle=p\cdot\left\|{\bm{D}}\nabla f(x^{k+1})-{\bm{D}}\nabla f(x^{k+1}) \right\|^{2}_{{\bm{D}}^{-1}}$ $\displaystyle\quad+(1-p)\cdot{\mathbb{E}}\left[\left\|g^{k}+\frac{1}{n}\sum_{i =1}^{n}{\bm{T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}) \right)-{\bm{D}}\nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}^{-1}}\mid x^{k+1},x^{k }\right]$ $\displaystyle=(1-p)\cdot{\mathbb{E}}\left[\left\|g^{k}+\frac{1}{n}\sum_{i=1}^{ n}{\bm{T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}) \right)-{\bm{D}}\nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}^{-1}}\mid x^{k+1},x^{k }\right].$ $$

For the sake of presentation, we use ${\rm E}_{k}\left[\cdot\right]$ to denote the conditional expectation ${\rm E}\left[\left.\cdot\right|x_{k},x_{k+1}\right]$ on $x_{k},x_{k+1}$.
Using 2 with $x=\frac{1}{n}\sum_{i=1}^{n}{\bm{T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})
-\nabla f_{i}(x^{k})\right)$, $c={\bm{D}}\nabla f(x^{k+1})-g^{k}$, we are able to obtain that,

$$ $\displaystyle(1-p){\rm E}_{k}\left[\left\|g^{k}+\frac{1}{n}\sum_{i=1}^{n}{\bm{ T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})\right)-{\bm {D}}\nabla f(x^{k+1})\right\|^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle=(1-p){\rm E}_{k}\left[\left\|\frac{1}{n}\sum_{i=1}^{n}{\bm{T}}_{ i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})\right)-{\bm{D}} \left(\nabla f(x^{k+1})-\nabla f(x^{k})\right)\right\|^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad+(1-p)\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}^{-1}}$ $\displaystyle=(1-p){\rm E}_{k}\left[\left\|\frac{1}{n}\sum_{i=1}^{n}\left[{\bm {T}}_{i}^{k}{\bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})\right)-{ \bm{D}}\left(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})\right)\right]\right\|^{ 2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad+(1-p)\left\|g^{k}-\nabla f(x^{k})\right\|^{2}_{{\bm{D}}^{-1 }}.$ $$

It is not hard to notice that for the sketch matrices we pick, the following identity holds due to the unbiasedness,

$$ ${\rm E}_{k}\left[{\bm{T}}_{i}^{k}{\bm{D}}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x ^{k}))\right]={\bm{D}}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})),$ $$

and any two random vectors in the set $\left\{{\bm{T}}_{i}^{k}{\bm{D}}(\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k}))
\right\}_{i=1}^{n}$ are independent if $x^{k+1},x^{k}$ are fixed. As a result

$$ $\displaystyle{\rm E}_{k}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle=\frac{1-p}{n^{2}}\sum_{i=1}^{n}{\rm E}_{k}\left[\left\|{\bm{T}}_ {i}^{k}\left({\bm{D}}\nabla f_{i}(x^{k+1})-{\bm{D}}\nabla f_{i}(x^{k})\right)- \left({\bm{D}}\nabla f_{i}(x^{k+1})-{\bm{D}}\nabla f_{i}(x^{k})\right)\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad+(1-p)\cdot\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2}_ {{\bm{D}}^{-1}}.$ (58) $$

For each term within the summation, we can further upper bound it
using Lemma 7

$$ $\displaystyle{\rm E}_{k}\left[\left\|{\bm{T}}_{i}^{k}\left({\bm{D}}\nabla f_{i }(x^{k+1})-{\bm{D}}\nabla f_{i}(x^{k})\right)-\left({\bm{D}}\nabla f_{i}(x^{k+ 1})-{\bm{D}}\nabla f_{i}(x^{k})\right)\right\|^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad\leq\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}{\bm{D}}{ \mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{\bm{T}}_{i}^{k}\right]{\bm{D}}{ \bm{L}}_{i}^{\frac{1}{2}}-{\bm{L}}_{i}^{\frac{1}{2}}{\bm{D}}{\bm{L}}_{i}^{ \frac{1}{2}}\right)\left\|\nabla f_{i}(x^{k+1})-\nabla f_{i}(x^{k})\right\|^{2 }_{{\bm{L}}_{i}^{-1}}$ $\displaystyle\quad\leq\lambda_{\max}\left({\bm{L}}_{i}^{\frac{1}{2}}{\bm{D}}{ \mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{\bm{T}}_{i}^{k}\right]{\bm{D}}{ \bm{L}}_{i}^{\frac{1}{2}}-{\bm{L}}_{i}^{\frac{1}{2}}{\bm{D}}{\bm{L}}_{i}^{ \frac{1}{2}}\right)\left\|x^{k+1}-x^{k}\right\|_{{\bm{L}}_{i}}^{2}.$ $$

Where the last inequality is due to 2.
Plugging back into (G.3), we get

$$ $\displaystyle{\rm E}_{k}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{ L}}_{i}^{\frac{1}{2}}{\bm{D}}{\mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{ \bm{T}}_{i}^{k}\right]{\bm{D}}{\bm{L}}_{i}^{\frac{1}{2}}-{\bm{L}}_{i}^{\frac{1 }{2}}{\bm{D}}{\bm{L}}_{i}^{\frac{1}{2}}\right)\left\|x^{k+1}-x^{k}\right\|_{{ \bm{L}}_{i}}^{2}$ $\displaystyle\qquad+(1-p)\cdot\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2} _{{\bm{D}}^{-1}}.$ $$

Applying the replacement trick form the proof of Theorem 1, we obtain

$$ $\displaystyle{\rm E}_{k}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{ i}^{\frac{1}{2}}{\bm{D}}{\mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{\bm{T} }_{i}^{k}\right]{\bm{D}}{\bm{L}}_{i}^{\frac{1}{2}}-{\bm{L}}_{i}^{\frac{1}{2}}{ \bm{D}}{\bm{L}}_{i}^{\frac{1}{2}}\right)$ $\displaystyle\quad\times\left\langle{\bm{L}}^{\frac{1}{2}}\left(x^{k+1}-x^{k} \right),\left({\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right )\cdot{\bm{L}}^{\frac{1}{2}}\left(x^{k+1}-x^{k}\right)\right\rangle+(1-p)\cdot \left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2}_{{\bm{D}}^{-1}}$ $\displaystyle\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{L}}_{ i}^{\frac{1}{2}}\left({\bm{D}}{\mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{ \bm{T}}_{i}^{k}\right]{\bm{D}}-{\bm{D}}\right){\bm{L}}_{i}^{\frac{1}{2}}\right )\cdot\lambda_{\max}\left({\bm{L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{ 1}{2}}\right)\left\|x^{k+1}-x^{k}\right\|_{{\bm{L}}}^{2}$ $\displaystyle\quad+(1-p)\cdot\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2}_ {{\bm{D}}^{-1}}.$ $$

Applying 5, we obtain

$$ $\displaystyle{\rm E}_{k}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\leq\frac{1-p}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{D}}{ \mathbb{E}}\left[{\bm{T}}_{i}^{k}{\bm{D}}^{-1}{\bm{T}}_{i}^{k}\right]{\bm{D}}- {\bm{D}}\right)\lambda_{\max}\left({\bm{L}}_{i}\right)\lambda_{\max}\left({\bm {L}}^{-\frac{1}{2}}{\bm{L}}_{i}{\bm{L}}^{-\frac{1}{2}}\right)\left\|x^{k+1}-x^ {k}\right\|_{{\bm{L}}}^{2}$ $\displaystyle\quad+(1-p)\cdot\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2}_ {{\bm{D}}^{-1}}.$ $$

Recalling the definition of $R^{\prime}({\bm{D}},\mathcal{S})$, we further simplify it to

$$ $\displaystyle{\rm E}_{k}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad\leq\frac{(1-p)\cdot R^{\prime}({\bm{D}},\mathcal{S})}{n} \left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}+(1-p)\cdot\left\|g^{k}-{\bm{D}} \nabla f(x^{k})\right\|^{2}_{{\bm{D}}^{-1}}.$ $$

Taking expectation again and using the tower property, we get

$$ $\displaystyle{\mathbb{E}}\left[\left\|g^{k+1}-{\bm{D}}\nabla f(x^{k+1})\right \|^{2}_{{\bm{D}}^{-1}}\right]$ (59) $\displaystyle\quad\leq\frac{(1-p)\cdot R^{\prime}({\bm{D}},\mathcal{S})}{n}{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]+(1-p)\cdot {\mathbb{E}}\left[\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{2}_{{\bm{D}}^{ -1}}\right].$ (60) $$

Construct the Lyapunov function $\Phi_{k}$ as follows,

$$ $\Phi_{k}=f(x^{k})-f^{\star}+\frac{1}{2p}\left\|g^{k}-{\bm{D}}\nabla f(x^{k}) \right\|^{2}_{{\bm{D}}^{-1}}.$ $$

Utilizing (G.3) and (59), we are able to get

$$ $\displaystyle{\mathbb{E}}\left[\Phi_{k+1}\right]$ $\displaystyle\leq{\mathbb{E}}\left[f(x^{k})-f^{\star}\right]-\frac{1}{2}{ \mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+\frac{1}{2}{\mathbb{E}}\left[\left\|g^{k}-{\bm{D}}\nabla f (x^{k})\right\|^{2}_{{\bm{D}}^{-1}}\right]-\frac{1}{2}{\mathbb{E}}\left[\left \|x^{k+1}-x^{k}\right\|_{{\bm{D}}^{-1}-L}\right]$ $\displaystyle\qquad+\frac{1}{2p}\cdot\frac{(1-p)R^{\prime}({\bm{D}},\mathcal{S })}{n}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]+ \frac{1-p}{2p}{\mathbb{E}}\left[\left\|g^{k}-{\bm{D}}\nabla f(x^{k})\right\|^{ 2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle\quad={\mathbb{E}}\left[\Phi_{k}\right]-\frac{1}{2}{\mathbb{E}} \left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+\frac{1}{2}\left(\frac{(1-p)R^{\prime}({\bm{D}},\mathcal{S })}{np}{\mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{L}}}\right]-{ \mathbb{E}}\left[\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}} \right]\right).$ $$

Now, notice that the last term in the above inequality is non-positive as guaranteed by the condition

$$ ${\bm{D}}^{-1}\succeq\left(\frac{(1-p)R^{\prime}({\bm{D}},\mathcal{S})}{np}+1 \right){\bm{L}}.$ $$

This leads to the recurrence after ignoring the last term,

$$ ${\mathbb{E}}\left[\Phi_{k+1}\right]\leq{\mathbb{E}}\left[\Phi_{k}\right]-\frac {1}{2}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}\right].$ $$

Unrolling this recurrence, we get

$$ $\frac{1}{K}\sum_{k=0}^{K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2 }_{{\bm{D}}}\right]\leq\frac{2\left({\mathbb{E}}\left[\Phi_{0}\right]-{\mathbb {E}}\left[\Phi_{K}\right]\right)}{K}.$ $$

The left hand side can viewed as average over $\tilde{x}^{K}$, which is drawn uniformly at random from $\{x_{k}\}_{k=0}^{K-1}$, while the right hand side can be simplified as

$$ $\displaystyle\frac{2\left({\mathbb{E}}\left[\Phi_{0}\right]-{\mathbb{E}}\left[ \Phi_{K}\right]\right)}{K}\leq\frac{2\Phi_{0}}{K}=\frac{2\left(f(x^{0})-f^{ \star}+\frac{1}{2p}\left\|g^{0}-\nabla f(x^{0})\right\|^{2}_{{\bm{D}}}\right)} {K}.$ $$

Recalling that $g^{0}=\nabla f(x^{0})$ and performing determinant normalization
as Li et al. (), we get

$$ ${\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|^{2}_{\frac{{\bm{D}}}{ \det({\bm{D}})^{1/d}}}\right]\leq\frac{2\left(f(x^{0})-f^{\star}\right)}{\det( {\bm{D}})^{1/d}K}.$ $$

## Appendix H Proofs of the technical lemmas

### H.1 Proof of Lemma 1

Let $\bar{x}^{k+1}:=x^{k}-{\bm{D}}\cdot\nabla f(x^{k})$.
Since $f$ has a matrix ${\bm{L}}$-Lipschitz gradient, $f$ is also ${\bm{L}}$-smooth.
From the ${\bm{L}}$-smoothness of $f$, we have

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k}),x^{k+1}-x^{k}\right\rangle+ \frac{1}{2}\left\langle x^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-g^{k},x^{k+1}-x^{k}\right \rangle+\left\langle g^{k},x^{k+1}-x^{k}\right\rangle+\frac{1}{2}\left\langle x ^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{k})\right\rangle.$ $$

We can merge the last two terms and obtain,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-g^{k},-{\bm{D}}\cdot g^{k} \right\rangle-\left\langle x^{k+1}-x^{k},{\bm{D}}^{-1}(x^{k+1}-x^{k})\right\rangle$ $\displaystyle+\frac{1}{2}\left\langle x^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-g^{k},-{\bm{D}}\cdot g^{k} \right\rangle-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2}{\bm{L }}\right)(x^{k+1}-x^{k})\right\rangle.$ $$

We add and subtract $\left\langle\nabla f(x^{k})-g^{k},{\bm{D}}\cdot g^{k}\right\rangle$,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-g^{k},{\bm{D}}\left(\nabla f (x^{k})-g^{k}\right)\right\rangle-\left\langle\nabla f(x^{k})-g^{k},{\bm{D}} \cdot\nabla f(x^{k})\right\rangle$ $\displaystyle-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2}{\bm{L }}\right)(x^{k+1}-x^{k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\|\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{D}}}- \left\langle x^{k+1}-\bar{x}^{k+1},{\bm{D}}^{-1}\left(x^{k}-\bar{x}^{k+1} \right)\right\rangle$ $\displaystyle-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2}{\bm{L }}\right)(x^{k+1}-x^{k})\right\rangle.$ $$

Decomposing the term $\left\langle x^{k+1}-\bar{x}^{k+1},{\bm{D}}^{-1}\left(x^{k}-\bar{x}^{k+1}
\right)\right\rangle$, we get

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\|\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{D}}}- \left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2}{\bm{L}}\right)(x^{k +1}-x^{k})\right\rangle$ $\displaystyle-\frac{1}{2}\left(\left\|x^{k+1}-\bar{x}^{k+1}\right\|^{2}_{{\bm{ D}}^{-1}}+\left\|x^{k}-\bar{x}^{k+1}\right\|^{2}_{{\bm{D}}^{-1}}-\left\|x^{k+1 }-x^{k}\right\|^{2}_{{\bm{D}}^{-1}}\right).$ $$

Plugging in the definition of $x^{k+1},\bar{x}^{k+1}$, we get

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\|\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{D}}}- \left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-\frac{1}{2}{\bm{L}}}$ $\displaystyle-\frac{1}{2}\left(\left\|{\bm{D}}(\nabla f(x^{k})-g^{k})\right\|^ {2}_{{\bm{D}}^{-1}}+\left\|{\bm{D}}\cdot\nabla f(x^{k})\right\|^{2}_{{\bm{D}}^ {-1}}-\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}}\right)$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\|\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{D}}}- \left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-\frac{1}{2}{\bm{L}}}$ $\displaystyle-\frac{1}{2}\left(\left\|\nabla f(x^{k})-g^{k}\right\|^{2}_{{\bm{ D}}}+\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}-\left\|x^{k+1}-x^{k}\right \|^{2}_{{\bm{D}}^{-1}}\right).$ $$

Rearranging terms we get,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})-\frac{1}{2}\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{2 }+\frac{1}{2}\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}-\left\|x^{k+1 }-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-\frac{1}{2}{\bm{L}}}+\frac{1}{2}\left\|x^{k +1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}}$ $\displaystyle=$ $\displaystyle f(x^{k})-\frac{1}{2}\left\|\nabla f(x^{k})\right\|_{{\bm{D}}}^{2 }+\frac{1}{2}\left\|g^{k}-\nabla f(x^{k})\right\|_{{\bm{D}}}^{2}-\frac{1}{2} \left\|x^{k+1}-x^{k}\right\|_{{\bm{D}}^{-1}-{\bm{L}}}.$ $$

### H.2 Proof of Lemma 2

The definition of the weighted norm yields

$$ $\displaystyle{\mathbb{E}}\left[\left\|{\bm{S}}t-t\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle=$ $\displaystyle{\mathbb{E}}\left[\left\langle t,\left({\bm{S}}-{\bm{I}}_{d} \right){\bm{D}}\left({\bm{S}}-{\bm{I}}_{d}\right)t\right\rangle\right]$ $\displaystyle=$ $\displaystyle\left\langle t,{\mathbb{E}}\left[({\bm{S}}-{\bm{I}}_{d}){\bm{D}}( {\bm{S}}-{\bm{I}}_{d})\right]t\right\rangle$ $\displaystyle=$ $\displaystyle\left\langle t,{\bm{L}}^{-\frac{1}{2}}\cdot{\mathbb{E}}\left[{\bm {L}}^{\frac{1}{2}}({\bm{S}}-{\bm{I}}_{d}){\bm{D}}({\bm{S}}-{\bm{I}}_{d}){\bm{L }}^{\frac{1}{2}}\right]\cdot{\bm{L}}^{-\frac{1}{2}}t\right\rangle$ $\displaystyle=$ $\displaystyle\left\langle{\bm{L}}^{-\frac{1}{2}}t,{\mathbb{E}}\left[{\bm{L}}^{ \frac{1}{2}}({\bm{S}}-{\bm{I}}_{d}){\bm{D}}({\bm{S}}-{\bm{I}}_{d}){\bm{L}}^{ \frac{1}{2}}\right]\cdot{\bm{L}}^{-\frac{1}{2}}t\right\rangle$ $\displaystyle\leq$ $\displaystyle\lambda_{\max}\left({\mathbb{E}}\left[{\bm{L}}^{\frac{1}{2}}({\bm {S}}-{\bm{I}}_{d}){\bm{D}}({\bm{S}}-{\bm{I}}_{d}){\bm{L}}^{\frac{1}{2}}\right] \right)\left\|{\bm{L}}^{-\frac{1}{2}}t\right\|^{2}$ $\displaystyle=$ $\displaystyle\lambda_{\max}\left({\bm{L}}^{\frac{1}{2}}\left({\mathbb{E}}\left [{\bm{S}}{\bm{D}}{\bm{S}}\right]-{\bm{D}}\right){\bm{L}}^{\frac{1}{2}}\right) \cdot\left\|t\right\|^{2}_{{\bm{L}}^{-1}}.$ $$

### H.3 Proof of Lemma 4

Throughout the following proof, we denote ${\mathbb{E}}_{\mathcal{S}}\left[\cdot\right]$ as taking expectation with respect to the randomness contained within the sketch sampled from distribution $\mathcal{S}$. We estimate the term ${\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$ in order to construct the Lyapunov function. For ${\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$, we have

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2 }_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k}+\frac{1}{n}\sum_{i= 1}^{n}m_{i}^{k+1}-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k}+\frac{1}{n}\sum_{i= 1}^{n}{\bm{S}}_{i}^{k}\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right )-h^{k+1}\right\|^{2}_{{\bm{D}}}\right]$ $$

Using 3, we obtain

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2 }_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n} {\bm{S}}_{i}^{k}\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)- \left(h^{k+1}-h^{k}-a(g^{k}-h^{k})\right)\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad+(1-a)^{2}\left\|h^{k}-g^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|\frac{1}{n}\sum_{i=1}^{n} {\bm{S}}_{i}^{k}\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)- \frac{1}{n}\sum_{i=1}^{n}\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k}) \right)\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad+(1-a)^{2}\left\|h^{k}-g^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle=\frac{1}{n^{2}}\sum_{i=1}^{n}{\mathbb{E}}_{\mathcal{S}}\left[ \left\|{\bm{S}}_{i}^{k}\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k}) \right)-\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)\right\|^{2}_ {{\bm{D}}}\right]$ $\displaystyle\quad+(1-a)^{2}\left\|h^{k}-g^{k}\right\|^{2}_{{\bm{D}}}.$ $$

Here, the last identity is obtained from the unbiasedness of the sketches:

$$ ${\mathbb{E}}_{\mathcal{S}}\left[{\bm{S}}_{i}^{k}\left(h_{i}^{k+1}-h_{i}^{k}-a( g_{i}^{k}-h_{i}^{k})\right)\right]=h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k}).$ $$

We can further use Lemma 2, and obtain

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2 }_{{\bm{D}}}\right]$ $\displaystyle\quad\leq\frac{1}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{D} }^{-\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D}}{\bm{S}}_{i}^{ k}\right]-{\bm{D}}\right){\bm{D}}^{-\frac{1}{2}}\right)\left\|h_{i}^{k+1}-h_{i }-a(g_{i}^{k}-h_{i}^{k})\right\|^{2}_{{\bm{D}}}$ $\displaystyle\qquad+(1-a)^{2}\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle\quad\leq\frac{1}{n^{2}}\sum_{i=1}^{n}\lambda_{\max}\left({\bm{D} }^{-1}\right)\cdot\lambda_{\max}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm{D }}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right)\left\|h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^ {k}-h_{i}^{k})\right\|^{2}_{{\bm{D}}}$ $\displaystyle\qquad+(1-a)^{2}\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}.$ $$

We can rewrite the above bound, after applying Jensen’s inequality as

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2 }_{{\bm{D}}}\right]$ $\displaystyle\leq\frac{2\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{\max}\left ({\bm{D}}^{-1}\right)}{n^{2}}\sum_{i=1}^{n}\left\|h_{i}^{k+1}-h_{i}^{k}\right \|^{2}_{{\bm{D}}}$ $\displaystyle\qquad+\frac{2a^{2}\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{ \max}\left({\bm{D}}^{-1}\right)}{n^{2}}\sum_{i=1}^{n}\left\|g_{i}^{k}-h_{i}^{k }\right\|^{2}_{{\bm{D}}}$ $\displaystyle\qquad+(1-a)^{2}\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}.$ $$

Notice that we have

$$ $\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\leq\lambda_{\max}\left({ \bm{D}}\right)\cdot\lambda_{\max}\left({\bm{L}}_{i}\right)\cdot\left\|h_{i}^{k +1}-h_{i}^{k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}.$ $$

Thus, it is not hard to see that

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g^{k+1}-h^{k+1}\right\|^{2 }_{{\bm{D}}}\right]$ $\displaystyle\leq\frac{2\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{\max}\left ({\bm{D}}^{-1}\right)\cdot\lambda_{\max}\left({\bm{D}}\right)}{n^{2}}\sum_{i=1 }^{n}\lambda_{\max}\left({\bm{L}}_{i}\right)\left\|h_{i}^{k+1}-h_{i}^{k}\right \|^{2}_{{\bm{L}}_{i}^{-1}}$ $\displaystyle\qquad+\frac{2a^{2}\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{ \max}\left({\bm{D}}^{-1}\right)}{n^{2}}\sum_{i=1}^{n}\left\|g_{i}^{k}-h_{i}^{k }\right\|^{2}_{{\bm{D}}}$ $\displaystyle\qquad+(1-a)^{2}\left\|g^{k}-h^{k}\right\|^{2}_{{\bm{D}}}.$ $$

We obtain the inequality in the lemma after taking expectation again and applying tower property.

### H.4 Proof of Lemma 5

Similarly, we then try to bound the terms ${\mathbb{E}}_{\mathcal{S}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1}\right\|^{2}_{{
\bm{D}}}\right]$. We start with

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1} \right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|g_{i}^{k}+{\bm{S}}_{i}^{k }\left(h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)-h_{i}^{k+1}\right\| ^{2}_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|{\bm{S}}_{i}^{k}\left(h_{ i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)-\left(h_{i}^{k+1}-h_{i}^{k}-a (g_{i}^{k}-h_{i}^{k})\right)+(1-a)(h_{i}^{k}-g_{i}^{k})\right\|^{2}_{{\bm{D}}} \right].$ $$

Using 3,

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1} \right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle={\mathbb{E}}_{\mathcal{S}}\left[\left\|{\bm{S}}_{i}^{k}\left(h_{ i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right)-\left(h_{i}^{k+1}-h_{i}^{k}-a (g_{i}^{k}-h_{i}^{k})\right)\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\quad+(1-a)^{2}\left\|h^{k}_{i}-g_{i}^{k}\right\|^{2}_{{\bm{D}}}.$ $$

Using Lemma 2

$$ $\displaystyle{\mathbb{E}}_{\mathcal{S}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1} \right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\overset{\eqref{dasha:eq:tech-lemma-sketch-comp}}{\leq}\lambda_{ \max}\left({\bm{D}}^{-\frac{1}{2}}\left({\mathbb{E}}\left[{\bm{S}}_{i}^{k}{\bm {D}}{\bm{S}}_{i}^{k}\right]-{\bm{D}}\right){\bm{D}}^{-\frac{1}{2}}\right)\left \|h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right\|^{2}_{{\bm{D}}}$ $\displaystyle\quad+(1-a)^{2}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle\leq\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda_{{\bm{D} },\mathcal{S}}\left\|h_{i}^{k+1}-h_{i}^{k}-a(g_{i}^{k}-h_{i}^{k})\right\|^{2}_ {{\bm{D}}}+(1-a)^{2}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle\leq 2\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda_{{\bm{ D}},\mathcal{S}}\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}+2a^{2} \lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda_{{\bm{D}},\mathcal{S}} \left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle\quad+(1-a)^{2}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle\leq 2\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\lambda_{\max} \left({\bm{D}}\right)\cdot\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{\max} \left({\bm{L}}_{i}\right)\cdot\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L }}_{i}^{-1}}$ $\displaystyle\quad+2a^{2}\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda_ {{\bm{D}},\mathcal{S}}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}+(1-a)^ {2}\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}$ $\displaystyle=\left(2a^{2}\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\Lambda _{{\bm{D}},\mathcal{S}}+(1-a)^{2}\right)\left\|g_{i}^{k}-h_{i}^{k}\right\|^{2} _{{\bm{D}}}$ $\displaystyle\quad+2\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\lambda_{\max }\left({\bm{D}}\right)\cdot\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{\max} \left({\bm{L}}_{i}\right)\cdot\left\|h_{i}^{k+1}-h_{i}^{k}\right\|^{2}_{{\bm{L }}_{i}^{-1}}.$ $$

Taking expectation again, and using tower property, we are able to obtain,

$$ $\displaystyle{\mathbb{E}}\left[\left\|g_{i}^{k+1}-h_{i}^{k+1}\right\|^{2}_{{ \bm{D}}}\right]$ $\displaystyle\quad\leq\left(2a^{2}\lambda_{\max}\left({\bm{D}}^{-1}\right) \cdot\Lambda_{{\bm{D}},\mathcal{S}}+(1-a)^{2}\right){\mathbb{E}}\left[\left\|g _{i}^{k}-h_{i}^{k}\right\|^{2}_{{\bm{D}}}\right]$ $\displaystyle\qquad+2\lambda_{\max}\left({\bm{D}}^{-1}\right)\cdot\lambda_{ \max}\left({\bm{D}}\right)\cdot\Lambda_{{\bm{D}},\mathcal{S}}\cdot\lambda_{ \max}\left({\bm{L}}_{i}\right)\cdot{\mathbb{E}}\left[\left\|h_{i}^{k+1}-h_{i}^ {k}\right\|^{2}_{{\bm{L}}_{i}^{-1}}\right].$ $$

### H.5 Proof of Lemma 6

From Proposition 5, we know that the objective is ${\bm{L}}$-smooth.
Let $\bar{x}^{k+1}=x^{k}-{\bm{D}}\cdot\nabla f(x^{k})$, then ${\bm{L}}$-smoothness yields

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k}),x^{k+1}-x^{k}\right\rangle+ \frac{1}{2}\left\langle x^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k},x^ {k+1}-x^{k}\right\rangle+\left\langle{\bm{D}}^{-1}\cdot g^{k},x^{k+1}-x^{k}\right\rangle$ $\displaystyle\qquad+\frac{1}{2}\left\langle x^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{ k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k},-g ^{k}\right\rangle-\left\langle x^{k+1}-x^{k},{\bm{D}}^{-1}(x^{k+1}-x^{k})\right\rangle$ $\displaystyle\qquad+\frac{1}{2}\left\langle x^{k+1}-x^{k},{\bm{L}}(x^{k+1}-x^{ k})\right\rangle.$ $$

Simplifying the above inner-products we have,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k},-g ^{k}\right\rangle-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2}{ \bm{L}}\right)(x^{k+1}-x^{k})\right\rangle.$ $$

We then add and subtract $\left\langle\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k},{\bm{D}}\cdot\nabla f(x^{
k})\right\rangle$, which

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\langle\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k},{ \bm{D}}\cdot\nabla f(x^{k})-g^{k}\right\rangle-\left\langle\nabla f(x^{k})-{ \bm{D}}^{-1}\cdot g^{k},{\bm{D}}\cdot\nabla f(x^{k})\right\rangle$ $\displaystyle\qquad-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2} {\bm{L}}\right)(x^{k+1}-x^{k})\right\rangle$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\|\nabla f(x^{k})-{\bm{D}}^{-1}\cdot g^{k}\right\| ^{2}_{{\bm{D}}}-\left\langle{\bm{D}}^{-1}(x^{k+1}-\bar{x}^{k+1}),x^{k}-\bar{x} ^{k+1}\right\rangle$ $\displaystyle\qquad-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^{-1}-\frac{1}{2} {\bm{L}}\right)(x^{k+1}-x^{k})\right\rangle.$ $$

Decomposing the inner product term we deduce,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\left\|{\bm{D}}^{-1}\left({\bm{D}}\cdot\nabla f(x^{k})- g^{k}\right)\right\|^{2}_{{\bm{D}}}-\left\langle x^{k+1}-x^{k},\left({\bm{D}}^ {-1}-\frac{1}{2}{\bm{L}}\right)(x^{k+1}-x^{k})\right\rangle$ $\displaystyle\qquad-\frac{1}{2}\left(\left\|x^{k+1}-\bar{x}^{k+1}\right\|^{2}_ {{\bm{D}}^{-1}}+\left\|x^{k}-\bar{x}^{k+1}\right\|^{2}_{{\bm{D}}^{-1}}-\left\| x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}}\right)$ $\displaystyle=$ $\displaystyle f(x^{k})+\left\|{\bm{D}}\cdot\nabla f(x^{k})-g^{k}\right\|^{2}_{ {\bm{D}}^{-1}}-\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-\frac{1}{2}{\bm {L}}}$ $\displaystyle\qquad-\frac{1}{2}\left(\left\|{\bm{D}}\cdot\nabla f(x^{k})-g^{k} \right\|^{2}_{{\bm{D}}^{-1}}+\left\|{\bm{D}}\cdot\nabla f(x^{k})\right\|^{2}_{ {\bm{D}}^{-1}}-\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}}\right).$ $$

Therefore,

$$ $\displaystyle f(x^{k+1})$ $\displaystyle\leq$ $\displaystyle f(x^{k})+\frac{1}{2}\left\|{\bm{D}}\nabla f(x^{k})-g^{k}\right\| ^{2}_{{\bm{D}}^{-1}}-\frac{1}{2}\left\|\nabla f(x^{k})\right\|^{2}_{{\bm{D}}}- \frac{1}{2}\left\|x^{k+1}-x^{k}\right\|^{2}_{{\bm{D}}^{-1}-{\bm{L}}}.$ $$

### H.6 Proof of Lemma 7

We start with

$$ $\displaystyle{\mathbb{E}}\left[\left\|{\bm{T}}{\bm{D}}t-{\bm{D}}t\right\|^{2}_ {{\bm{D}}^{-1}}\right]$ $\displaystyle=$ $\displaystyle{\mathbb{E}}\left[\left\|({\bm{T}}-{\bm{I}}_{d}){\bm{D}}t\right\| ^{2}_{{\bm{D}}^{-1}}\right]$ $\displaystyle=$ $\displaystyle\left\langle t,{\mathbb{E}}\left[{\bm{D}}({\bm{T}}-{\bm{I}}_{d}){ \bm{D}}^{-1}({\bm{T}}-{\bm{I}}_{d}){\bm{D}}\right]\cdot t\right\rangle$ $\displaystyle=$ $\displaystyle\left\langle t,{\bm{D}}\left({\mathbb{E}}\left[{\bm{T}}{\bm{D}}^{ -1}{\bm{T}}\right]-{\bm{D}}^{-1}\right){\bm{D}}\cdot t\right\rangle$ $\displaystyle=$ $\displaystyle\left\langle{\bm{L}}^{-\frac{1}{2}}t,{\bm{L}}^{\frac{1}{2}}{\bm{D }}\left({\mathbb{E}}\left[{\bm{T}}{\bm{D}}^{-1}{\bm{T}}\right]-{\bm{D}}^{-1} \right){\bm{D}}{\bm{L}}^{\frac{1}{2}}\cdot{\bm{L}}^{-\frac{1}{2}}t\right\rangle$ $\displaystyle\leq$ $\displaystyle\lambda_{\max}\left({\bm{L}}^{\frac{1}{2}}{\bm{D}}{\mathbb{E}} \left[{\bm{T}}{\bm{D}}^{-1}{\bm{T}}\right]{\bm{D}}{\bm{L}}^{\frac{1}{2}}-{\bm{ L}}^{\frac{1}{2}}{\bm{D}}{\bm{L}}^{\frac{1}{2}}\right)\cdot\left\|{\bm{L}}^{- \frac{1}{2}}t\right\|^{2}$ $\displaystyle=$ $\displaystyle\lambda_{\max}\left({\bm{L}}^{\frac{1}{2}}{\bm{D}}{\mathbb{E}} \left[{\bm{T}}{\bm{D}}^{-1}{\bm{T}}\right]{\bm{D}}{\bm{L}}^{\frac{1}{2}}-{\bm{ L}}^{\frac{1}{2}}{\bm{D}}{\bm{L}}^{\frac{1}{2}}\right)\cdot\left\|t\right\|^{2 }_{{\bm{L}}^{-1}}$ $$

This completes the proof.

## Appendix I Experiments

In this section, we conduct numerical experiments to back up the theoretical results for det-MARINA and det-DASHA.
The code for the experiments can be found in [https://anonymous.4open.science/r/detCGD-VR-Code-865B](https://anonymous.4open.science/r/detCGD-VR-Code-865B).
All the codes for the experiments are written in Python 3.11 with NumPy and SciPy package.
The code was run on a machine with AMD Ryzen 9 5900HX Radeon Graphics @ 3.3 GHz and 8 cores 16 threads.
The datasets in LibSVM are typically non-IID real world datasets, and it is randomly distributed across all the clients.

### I.1 The setting

We first state the experiment setting.
We are interested in the following logistic regression problem with a non-convex regularizer. The objective is given as

$$ $f(x)=\frac{1}{n}\sum_{i=1}^{n}f_{i}(x);\qquad f_{i}(x)=\frac{1}{m_{i}}\sum_{j= 1}^{m_{i}}\log\left(1+e^{-b_{i,j}\cdot\left\langle a_{i,j},x\right\rangle} \right)+\lambda\cdot\sum_{t=1}^{d}\frac{x_{t}^{2}}{1+x_{t}^{2}},$ $$

where $x\in\mathbb{R}^{d}$ is the model, $\left(a_{i,j},b_{i,j}\right)\in\mathbb{R}^{d}\times\left\{-1,1\right\}$ is one data point in the dataset of client $i$ whose size is $m_{i}$.
The constant $\lambda>0$ is the coefficient of the regularizer.
Larger $\lambda$ means the model is more regular.
For each function $f_{i}$, its Hessian can be upper bounded by

$$ $\displaystyle{\bm{L}}_{i}=\frac{1}{m_{i}}\sum_{i=1}^{m_{i}}\frac{a_{i}a_{i}^{ \top}}{4}+2\lambda\cdot{\bm{I}}_{d};$ $$

and, therefore, the Hessian of $f$ is bounded by

$$ $\displaystyle{\bm{L}}=\frac{1}{\sum_{i=1}^{n}m_{i}}\sum_{i=1}^{n}\sum_{j=1}^{m _{i}}\frac{a_{i}a_{i}^{\top}}{4}+2\lambda\cdot{\bm{I}}_{d}.$ $$

Due to Proposition 2, it immediately follows that $f_{i}$ and $f$ satisfy Definition 1 with ${\bm{L}}_{i}\in\mathbb{S}^{d}_{++}$ and ${\bm{L}}\in\mathbb{S}^{d}_{++}$, respectively.

In the following subsections, we perform several numerical experiments comparing the performance of DCGD, det-CGD, MARINA, DASHA, det-MARINA and det-DASHA.
The datasets we used are from the LibSVM repository (Chang & Lin, ).

### I.2 Comparison of all the methods

In this section, we present several plots which compare all relevant methods to the det-MARINA and det-DASHA.
The methods are the following: (i) DCGD with scalar stepsize $\gamma_{2}$, *(ii)* det-CGD with matrix stepsize ${\bm{D}}_{3}^{*}$, *(iii)* MARINA with scalar stepsize $\gamma_{1}$, *(iv)* DASHA with scalar stepsize $\gamma_{4}$, *(v)* det-MARINA with ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$, *(vi)* det-DASHA with ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$.
Throughout the experiment, $\varepsilon=0.01$, and $\lambda=0.9$, we are using the same Rand-$\tau$ sketch for all the algorithms, and we run all the algorithms for a fixed number of iteration $K=10000$.

Figure: Figure 2: Comparison of DCGD with optimal scalar stepsize, det-CGD with matrix stepsize ${\bm{D}}^{*}_{3}$, MARINA with optimal scalar stepsize, DASHA with optimal scalar stepsize, det-MARINA with optimal stepsize ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ and det-DASHA with optimal stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$. Throughout the experiment, we are using Rand-$\tau$ sketch with $\tau=60$, and each algorithm is run for a fixed number of iterations $K=10000$. The momentum of DASHA is set as $\nicefrac{{1}}{{2\omega+1}}$ and det-DASHA is $\nicefrac{{1}}{{2\omega_{{\bm{D}}}+1}}$. The notation $n$ in the title stands for the number of clients in each case, and $p$ stands for the probability used by MARINA and det-MARINA.
Refer to caption: x4.png

It can be seen in Figure 2, the performance in terms of communication complexity of det-DASHA and det-MARINA is better than their scalar counterpart DASHA and MARINA respectively.
This validates the efficiency of using a matrix stepsize over a scalar stepsize.
Furthermore, we notice that det-DASHA and det-MARINA have better communication complexity in this case, compared to det-CGD.
In addition, we observe variance reduction.

Notice that the optimal stepsizes of det-CGD and DCGD require information of function value differences at
$x^{\star}$.
Furthermore, the stepsizes are also constrained by the number of iterations $K$ and the error $\varepsilon^{2}$.
Meanwhile, for the variance reduced methods, we do not require such considerations, which is much more practical in general.

### I.3 Improvements over MARINA

The purpose of this experiment is to compare the iteration complexity of MARINA, with det-MARINA using Rand-$\tau$ sketches, thus showing improvements of det-MARINA upon MARINA.
Using Theorem C.1 from (Gorbunov et al., ), we deduce the optimal stepsize for MARINA, is

$$ $\gamma_{1}=\frac{1}{L\left(1+\sqrt{\frac{(1-p)\omega}{pn}}\right)},$ (61) $$

where $\omega$ is the quantization coefficient.
In particular, for the Rand-$\tau$ compressor $\omega=\frac{d}{\tau}-1$.
For the full definition see Section 1.3 of (Gorbunov et al., ).
The stepsize for det-MARINA is determined through Corollary 1.
We use the notation ${\bm{D}}_{{\bm{W}}}^{*}$ to denote the optimal stepsize for each choice of ${\bm{W}}$, here we list some of the optimal stepsizes for different ${\bm{W}}$, which are used in the experiment section.
We have

$$ $\displaystyle{\bm{D}}_{{\bm{I}}_{d}}^{*}$ $\displaystyle=$ $\displaystyle\frac{2}{1+\sqrt{1+4\alpha\beta\frac{1}{\lambda_{\max}\left({\bm{ L}}\right)}\cdot\omega}}\cdot\frac{{\bm{I}}_{d}}{\lambda_{\max}({\bm{L}})},$ $\displaystyle{\bm{D}}^{*}_{{\bm{L}}^{-1}}$ $\displaystyle=$ $\displaystyle\frac{2}{1+\sqrt{1+4\alpha\beta\cdot\lambda_{\max}\left({\mathbb{ E}}\left[{\bm{S}}_{i}^{k}{\bm{L}}^{-1}{\bm{S}}_{i}^{k}\right]-{\bm{L}}^{-1} \right)}}\cdot{\bm{L}}^{-1},$ $\displaystyle{\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}$ $\displaystyle=$ $\displaystyle\frac{2}{1+\sqrt{1+4\alpha\beta\cdot\lambda_{\max}\left({\mathbb{ E}}\left[{\bm{S}}_{i}^{k}\operatorname{diag}^{-1}\left({\bm{L}}\right){\bm{S}} _{i}^{k}\right]-\operatorname{diag}^{-1}\left({\bm{L}}\right)\right)}}\cdot \operatorname{diag}^{-1}\left({\bm{L}}\right).$ (62) $$

In this experiment, we aim to compare det-MARINA with stepsize ${\bm{D}}_{{\bm{L}}^{-1}}^{*}$ to the standard MARINA  with the optimal scalar stepsize.
Rand-$\tau$ compressor is used in the comparison.
Throughout the experiments, $\lambda$ is fixed at $0.3$.
We set the $x$-axis to be the number of iterations, while $y$-axis to be the expectation of the corresponding matrix norm of the gradient of the function, which is defined as

$$ $G_{K,{\bm{D}}}={\mathbb{E}}\left[\left\|\nabla f(\tilde{x}^{K})\right\|_{{\bm{ D}}/\det({\bm{D}})^{1/d}}^{2}\right].$ (63) $$

Notice that this criterion is comparable to the standard Euclidean norm Li et al. (), and for a fixed ${\bm{D}}$, we have

$$ $\displaystyle\lambda_{\min}\left(\frac{{\bm{D}}}{\det({\bm{D}})^{1/d}}\right) \cdot\left\|\nabla f(x)\right\|^{2}\leq\left\|\nabla f(x)\right\|^{2}_{\frac{{ \bm{D}}}{\det({\bm{D}})^{1/d}}}\leq\lambda_{\max}\left(\frac{{\bm{D}}}{\det({ \bm{D}})^{1/d}}\right)\cdot\left\|\nabla f(x)\right\|^{2}.$ $$

Figure: Figure 3: In this experiment, we aim to compare det-MARINA with stepsize ${\bm{D}}_{{\bm{L}}^{-1}}^{*}$ to the standard MARINA with the optimal scalar stepsize. Rand-$\tau$ compressor is used in the comparison. Throughout the experiments, $\lambda$ is fixed at $0.3$. Optimal stepsize is calculated in each case with respect to the sketch used. The $x$-axis denotes the number of iterations while the notation $G_{K,{\bm{D}}}$ for the $y$-axis is defined in (63), which is the averaged matrix norm of the gradient. The notation $p$ in the title denotes the probability used in the two algorithms, $n$ denotes the number of clients in each setting.
Refer to caption: x8.png

As it is illustrated in Figure 3, det-MARINA always has a faster convergence rate compared to MARINA if they use the same sketch, this justifies the result we have in Corollary 3.
Notice that in some cases, det-MARINA with Rand-$1$ sketch even outperforms standard MARINA with Rand-$80$ sketch.
This further demonstrates the superiority of matrix stepsizes and smoothness over the standard scalar setting.

### I.4 Improvements on non variance reduced methods

In this section, we compare two non-variance reduced methods, distributed compressed gradient descent ( DCGD) and distributed det-CGD, with two variance reduced methods, MARINA, and det-MARINA.
Rand-$1$ sketch is used throughout this experiment for all the algorithms, for non variance reduced method $\varepsilon^{2}$ is fixed at $0.01$ in order to determine the optimal stepsize.
The purpose of this experiment is to show the advantages of variance reduced methods over non variance reduced methods.
DCGD was initially proposed in (Khirirat et al., ).
Later on DIANA was proposed in (Mishchenko et al., ) and then combined with variance reduction technique. Recently Shulgin & Richtárik () proposed shifted DCGD, which is a shifted version of DCGD and proved its convergence in the (strongly) convex setting. A general analysis on SGD type methods in the non-convex world is provided by Khaled & Richtárik (), including DCGD and shifted DCGD. In our case, in order to determine the optimal scalar stepsize for DCGD, one can simply use Proposition 4 in (Khaled & Richtárik, ). One can check that in order to satisfy $\min_{0\leq k\leq K-1}{\mathbb{E}}\left[\left\|\nabla f(x^{k})\right\|^{2}
\right]\leq\varepsilon^{2}$ the stepsize condition for DCGD in the non-convex case reduces to

$$ $\gamma_{2}\leq\min\left\{\frac{1}{L},\sqrt{\frac{n}{\omega LL_{\max}K}},\frac{ n\varepsilon^{2}}{4LL_{\max}\omega\cdot\Delta^{\star}}\right\},$ $$

where $L$ is the smoothness constant for $f$, $L_{i}$ is the smoothness constant for $f_{i}$, $L_{\max}=\max_{i}L_{i}$, $K$ is the total number of iterations, $\Delta^{\star}=f(x^{\star})-\frac{1}{n}\sum_{i=1}^{n}f_{i}(x^{\star})$. The constant $\omega$ is associated with the compressor used in the algorithm, for Rand-$\tau$ sketch, it is $\frac{d}{\tau}-1$.
For distributed det-CGD  according to Li et al. (), the stepsize condition in order to satisfy $\min_{0\leq k\leq K-1}{\mathbb{E}}\left[\left\|\nabla f(x)\right\|^{2}_{{\bm{D
}}/\det({\bm{D}})^{1/d}}\right]\leq\varepsilon^{2}$ is

$$ ${\bm{D}}{\bm{L}}{\bm{D}}\preceq{\bm{D}},\qquad\lambda_{{\bm{D}}}\leq\min\left \{\frac{n}{K},\frac{n\varepsilon^{2}}{4\Delta^{\star}}\det({\bm{D}})^{1/d} \right\},$ (64) $$

where $\lambda_{{\bm{D}}}$ is defined as

$$ $\lambda_{{\bm{D}}}=\max_{i}\left\{\lambda_{\max}\left({\mathbb{E}}\left[{\bm{L }}_{i}^{\frac{1}{2}}\left({\bm{S}}_{i}^{k}-{\bm{I}}_{d}\right){\bm{D}}{\bm{L}} {\bm{D}}\left({\bm{S}}_{i}^{k}-{\bm{I}}_{d}\right){\bm{L}}_{i}^{\frac{1}{2}} \right]\right)\right\}.$ (65) $$

In general cases, there is no easy way to find a optimal stepsize matrix ${\bm{D}}$ satisfying (64), alternatively, we choose the optimal diagonal stepsize ${\bm{D}}^{*}_{3}$ similarly to (Li et al., ). The stepsize condition for MARINA has already been described by (61). Note that we only consider MARINA, but not DIANA or shifted DCGD, because DIANA and shifted DCGD offer suboptimal rates compared to MARINA in the non-convex setting. For det-MARINA, we fix ${\bm{W}}={\bm{L}}^{-1}$, and use ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ as the stepsize matrix. In theory, det-MARINA in this case should always out perform MARINA  in terms of iteration complexity.

Figure: Figure 4: Comparison of DCGD with optimal scalar stepsize $\gamma_{2}$, det-CGD with optimal diagonal stepsize ${\bm{D}}_{3}^{*}$, MARINA with optimal scalar stepsize $\gamma_{1}$, and det-MARINA with optimal stepsize ${\bm{D}}_{{\bm{L}}^{-1}}^{*}$ with respect to ${\bm{W}}={\bm{L}}^{-1}$. In each case, probability $p$ is chosen the set $\{0.05,0.2,0.5\}$ for MARINA and det-MARINA. $\lambda=0.3$ is fixed throughout the experiment. The notation $n$ in the title indicates the number of clients in each case.
Refer to caption: x20.png

In Figure 4, in each plot, we observe that det-MARINA outperforms MARINA and the rest of the non-variance reduced methods.
This is expected, since our theory confirms that det-MARINA indeed has a better rate compared to MARINA, and the stepsizes of the non-variance reduced methods are negatively affected by the neighborhood.
When $p$ is reasonably large, the variance reduced methods considered here outperform the non-variance reduced methods.
In this experiment we consider only the comparison involving det-CGD.

### I.5 Improvements over det-CGD

In this section, we compare det-CGD in the distributed case with det-MARINA, which are both algorithms using matrix stepsizes and matrix smoothness.
The purpose of this experiment is to show that det-MARINA improves on the current state of the art matrix stepsize compressed gradient method when the objective function is non-convex.
Throughout the experiment, $\lambda=0.3$ is fixed, and for det-CGD, $\varepsilon^{2}=0.01$ is fixed in order to determine its stepsize.
For a thorough comparison, we select the stepsize for det-CGD in the following way.
Let us denote the stepsize as ${\bm{D}}=\gamma_{{\bm{W}}}\cdot{\bm{W}}$, where $\gamma_{{\bm{W}}}\in\mathbb{R}_{++},{\bm{W}}\in\mathbb{S}^{d}_{++}$.
We first fix a matrix ${\bm{W}}$, in this case, we pick ${\bm{W}}$ from the set $\{{\bm{L}}^{-1},\operatorname{diag}^{-1}({\bm{L}}),{\bm{I}}_{d}\}$, and then we determine the optimal scaling $\gamma_{{\bm{W}}}$ for each case using the condition given in (Li et al., ) (see (64) and (65)).
Then, we denote the matrix stepsizes for det-CGD

$$ ${\bm{D}}_{1}=\gamma_{{\bm{I}}_{d}}\cdot{\bm{I}}_{d},\qquad{\bm{D}}_{2}=\gamma_ {\operatorname{diag}^{-1}({\bm{L}})}\cdot\operatorname{diag}^{-1}\left({\bm{L} }\right),\qquad{\bm{D}}_{3}=\gamma_{{\bm{L}}^{-1}}\cdot{\bm{L}}^{-1}.$ (66) $$

For det-MARINA, we use the stepsize ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$, which is described in (I.3). In this experiment, we compare det-CGD using three stepsizes ${\bm{D}}_{1},{\bm{D}}_{2},{\bm{D}}_{3}$ with det-MARINA using stepsize ${\bm{D}}_{{\bm{L}}^{-1}}^{*}$.

Figure: Figure 5: Comparison of det-CGD with matrix stepsize ${\bm{D}}_{1}$, ${\bm{D}}_{2}$ and ${\bm{D}}_{3}$ and det-MARINA with optimal matrix stepsize with respect to ${\bm{W}}={\bm{L}}^{-1}$. The stepsizes $\{{\bm{D}}_{i}\}_{i=1}^{3}$ are described in (66). Throughout the experiment $\varepsilon^{2}$ is fixed at $0.01$, the notation $p$ in the title refers to the probability for det-MARINA, $n$ denotes the number of clients considered, Rand-$1$ sketch is used in all cases for all the algorithms.
Refer to caption: x26.png

From Figure 5, it is clear that det-MARINA outperforms det-CGD with all matrix optimal stepsizes with respect to a fixed ${\bm{W}}$ considered here.
This is expected, since the convergence rate of non-variance reduced methods are affected by its neighborhood.
This experiment demonstrates the advantages of det-MARINA over det-CGD, and is also supported by our theory.
Notice that though different ${\bm{W}}$ are considered for det-CGD, their convergence rates are similar, which is also mentioned by Li et al. ().

### I.6 Comparing different stepsize choices

This experiment is designed to see the how det-MARINA works under different stepsize choices.
As it is mentioned in Section I.3, for each choice of ${\bm{W}}\in\mathbb{S}^{d}_{++}$, an optimal stepsize ${\bm{D}}^{*}_{{\bm{W}}}$ can be determined.
Here we compare det-MARINA using three different stepsize choices ${\bm{D}}^{*}_{{\bm{L}}^{-1}},{\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}$ and ${\bm{D}}^{*}_{{\bm{I}}_{d}}$.
There stepsizes are explicitly defined in (I.3).
Throughout the experiment, we fix $\lambda=0.3$, Rand-$1$ sketch is used in all cases.

Figure: Figure 6: Comparison of det-MARINA with matrix stepsize ${\bm{D}}^{*}_{{\bm{I}}_{d}}$, ${\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}$ and ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$. The stepsizes are defined in (I.3). Throughout the experiment, $\lambda=0.3$ is fixed, Rand-$1$ sketch is used in all cases. The notation $p$ in the title indicates the probability of sending the true gradient for det-MARINA, $n$ denotes the number of clients considered.
Refer to caption: x38.png

We can observe from Figure 6 that, in almost all cases det-MARINA with stepsize ${\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}$ and ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ outperforms det-MARINA with ${\bm{D}}^{*}_{{\bm{I}}_{d}}$.
As det-MARINA with ${\bm{D}}^{*}_{{\bm{I}}_{d}}$ can be viewed as MARINA using scalar stepsize but under matrix Lipschitz gradient assumption, this demonstrates the effectiveness of using a matrix stepsize over the scalar stepsize.
However, in Figure 6, there are cases where det-MARINA with ${\bm{D}}^{*}_{\operatorname{diag}^{-1}({\bm{L}})}$ outperforms ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$.
This tells us the two stepsizes are perhaps incomparable in general cases.
This is similar to det-CGD, where optimal stepsizes with respect to a subspace associated with a fixed ${\bm{W}}^{-1}$ are incomparable.

### I.7 Comparing communication complexity

In this section, we perform an experiment on how different probabilities $p$ will affect the overall communication complexity of det-MARINA.
We use ${\bm{D}}_{{\bm{L}}^{-1}}^{*}$ as the stepsize, which is determined with respect to the sketch used. Rand-$\tau$ sketches are used in these experiments, and we vary the minibatch size $\tau$ to provide a more comprehensive comparison.
For Rand-$\tau$ sketch ${\bm{S}}$ and any ${\bm{A}}\in\mathbb{S}^{d}_{++}$, one can show that

$$ ${\mathbb{E}}\left[{\bm{S}}{\bm{A}}{\bm{S}}\right]=\frac{d}{\tau}\left(\frac{d- \tau}{\operatorname{d}\!{-}1}\operatorname{diag}({\bm{A}})+\frac{\tau-1}{d-1}{ \bm{A}}\right).$ (67) $$

Combining (67) and (I.3), we can find out the corresponding matrix stepsize easily. In the experiment, a fixed number of iterations ($K=5000$) is performed for each det-MARINA with the corresponding stepsize.

Figure: Figure 7: Comparison of det-MARINA with stepsize ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ using different probability $p$. The probability $p$ here is chosen from the set $\{0.05,0.1,0.2,0.4,0.8\}$. The notation $n$ in the title denote the number of clients considered. The $x$-axis is now the number of bytes sent from a single node to the server. In each case, det-MARINA is run for a fixed number of iterations $K=5000$.
Refer to caption: x47.png

As it can be observed from Figure 7, in each dataset, the communication complexity tends to increase with the increase of probability $p$.
However, when the number of iteration is fixed, a larger $p$ often means a faster rate of convergence.
This difference in communication complexity is more obvious when we are using the Rand-$1$ sketch.
In real federated learning settings, there is often constraints on network bandwidth from clients to the server. Thus, trading off between communication complexity and iteration complexity, i.e. selecting the compression mechanism carefully to guarantee a acceptable speed that satisfies the bandwidth constraints, becomes important.

### I.8 Comparison of DASHA and det-DASHA

In this experiment we plan to compare the performance of original DASHA with det-DASHA. Throughout the experiments, $\lambda$ is fixed at $0.3$. The same Rand-$\tau$ sketch is used in the two algorithms. The stepsize condition on DASHA when the momentum is set as $a=\frac{1}{2\omega+1}$ is given as

$$ $\gamma_{4}\leq\left(L+\sqrt{\frac{16\omega(2\omega+1)}{n}}\widehat{L}\right)^{ -1},$ $$

according to Theorem 6.1 of Tyurin & Richtárik ().
Here the $L$ is the smoothness constant of the function $f$, while $\widehat{L}$ satisfies $\widehat{L}^{2}=\frac{1}{n}\sum_{i=1}^{n}L_{i}^{2}$ where $L_{i}$ is the smoothness constant of local objective $f_{i}$.
In theory we can pick $\widehat{L}=L$.
Similarly, according to Corollary 2, the optimal stepsize matrix ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ is given as

$$ ${\bm{D}}^{**}_{{\bm{L}}^{-1}}=\frac{2}{1+\sqrt{1+16C_{{\bm{L}}^{-1}}\cdot \lambda_{\min}\left({\bm{L}}\right)}}\cdot{\bm{L}}^{-1},$ (68) $$

when the momentum is given as $a=\frac{1}{2\omega_{{\bm{D}}}+1}$. We compare the performance of DASHA with $\omega$ and det-DASHA with ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ using the same sketch where the total number of clients are different.

Figure: Figure 8: Comparison of det-DASHA with matrix stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ and DASHA with optimal scalar stepsize $\gamma$ using different Rand-$\tau$ sketches. $\lambda=0.3$ is fixed throughout the experiments. Optimal stepsize is calculated in each case with respect to the sketch used. The $x$-axis denotes the number of iterations while the notation $G_{K,{\bm{D}}}$ for the $y$-axis denotes the averaged matrix norm of the gradient. The notation $n$ denotes the number of clients in each setting.
Refer to caption: x53.png

As it can be observed in Figure 8, det-DASHA with matrix stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ outperforms DASHA with optimal scalar stepsize using the same sketch in every setting we considered. Note that since the same sketch is used in the two algorithm, the number of bits transferred in each iteration is also the same for the two algorithms. This essentially indicates that det-DASHA has better iteration complexity as well as communication complexity than DASHA given that the same sketch is used for the two algorithm.

### I.9 Comparison of DCGD , det-CGD , DASHA and det-DASHA

In this experiment, we consider the comparison between the two non variance reduced methods DCGD, det-CGD and the two variance reduced method DASHA, det-DASHA. The stepsize choices for DCGD and det-CGD have already been discussed in the previous sections, (for DCGD we use $\gamma_{2}$ and for det-CGD we use ${\bm{D}}^{*}_{3}$) ,for DASHA and det-DASHA, we use the stepsize choices of Section I.8. Note that $\varepsilon^{2}$ is set as $0.01$, and $\lambda$ is fixed at $0.9$ here. Throughout this experiment, we consider the case where Rand-$\tau$ sketch is used in the four algorithms.

Figure: Figure 9: Comparison of DCGD with optimal scalar stepsize $\gamma_{2}$, det-CGD with optimal diagonal stepsize ${\bm{D}}_{3}^{*}$, DASHA with optimal scalar stepsize $\gamma_{1}$, and det-DASHA with optimal stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ with respect to ${\bm{W}}={\bm{L}}^{-1}$. $\lambda=0.9$ is fixed throughout the experiment. The notation $n$ in the title indicates the number of clients in each case. Rand-$\tau$ sketch with $\tau=50$ are used in all four algorithms.
Refer to caption: x62.png

It is easy to observe that in each case of Figure 9, det-DASHA outperforms the rest of the algorithms. It is expected that det-DASHA outperforms DASHA, as it is also illustrated by Figure 8, which is a consequence of using matrix stepsize instead of a scalar stepsize. We also see that det-DASHA and DASHA outperform det-CGD and DCGD respectively, which demonstrate the advantages of the variance reduction technique. Note that in this case, all four algorithms are using the same sketch, which means that the number of bits transferred in each iteration is the same for the four algorithms, as a result, compared to the other algorithms, det-DASHA is better in terms of both iteration complexity and communication complexity.

### I.10 Comparison of det-DASHA and det-CGD with different stepsizes

In this experiment, we try to compare det-DASHA and det-CGD with different matrix stepsizes. Throughout this experiment, we will fix $\varepsilon^{2}=0.01$ and $\lambda=0.9$. The same Rand-$\tau$ sketch is used for the two algorithms. For det-CGD, we use the stepsize ${\bm{D}}_{1}=\gamma_{{\bm{I}}_{d}}\cdot{\bm{I}}_{d},{\bm{D}}_{2}=\gamma_{
\operatorname{diag}^{-1}\left({\bm{L}}\right)}\cdot\operatorname{diag}^{-1}
\left({\bm{L}}\right)$ and ${\bm{D}}_{3}=\gamma_{{\bm{L}}^{-1}}\cdot{\bm{L}}^{-1}$, for det-DASHA we use the stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$.

Figure: Figure 10: Comparison of det-DASHA with stepsize ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ and det-CGD with three different stepsizes ${\bm{D}}_{1}$, ${\bm{D}}_{2}$ and ${\bm{D}}_{3}$. Throughout the experiment, $\lambda$ is fixed at $0.9$, $\varepsilon^{2}$ is fixed at $0.01$, Rand-$\tau$ sketch is used for all the algorithms with $\tau$ selected from $\{20,50,80\}$. The notation $n$ denotes the number of clients in each setting.
Refer to caption: x68.png

It can be observed that in all cases of Figure 10, det-DASHA outperforms det-CGD with different stepsizes. This further corroborates our theory that det-DASHA is variance reduced and thus is better in terms of both iteration complexity, and communication complexity (because in this case the same number of bits are transmitted in each iteration due to the fact that the sketch used is the same).

### I.11 Comparison of different stepsizes of det-DASHA

In this experiment, we try to compare det-DASHA with different matrix stepsizes. Specifically, we fix matrix ${\bm{W}}$ to be three different matrices, ${\bm{I}}_{d}$, $\operatorname{diag}^{-1}\left({\bm{L}}\right)$ and ${\bm{L}}^{-1}$. We denote the optimal stepsizes as ${\bm{D}}^{**}_{{\bm{I}}_{d}}$, ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$ and ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$, respectively. For ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$, it is already given in (68), for ${\bm{D}}^{**}_{{\bm{I}}_{d}}$ and ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$, we use Corollary 2 to compute them. As a result,

$$ $\displaystyle{\bm{D}}^{**}_{{\bm{I}}_{d}}=\frac{2}{1+\sqrt{1+16\cdot\frac{ \omega_{{\bm{I}}_{d}}\left(4\omega_{{\bm{I}}_{d}}+1\right)}{n}}\cdot\frac{ \lambda_{\min}\left({\bm{L}}\right)}{\lambda_{\max}\left({\bm{L}}\right)}} \cdot\frac{{\bm{I}}_{d}}{\lambda_{\max}\left({\bm{L}}\right)},$ (69) $$

and

$$ $\displaystyle{\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}= \frac{2}{1+\sqrt{1+16C_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}\cdot \lambda_{\min}\left({\bm{L}}\right)}}\cdot\operatorname{diag}^{-1}\left({\bm{L }}\right).$ (70) $$

Throughout the experiment, $\lambda$ is fixed at $0.9$, Rand-$\tau$ sketch is used for all the algorithms.

Figure: Figure 11: Comparison of det-DASHA three different stepsizes ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$, ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$ and ${\bm{D}}^{**}_{{\bm{I}}_{d}}$. The definition for those matrix stepsize notation are given in (68), (70) and (69) respectively. Throughout the experiment, $\lambda$ is fixed at $0.9$, Rand-$\tau$ sketch is used for all the algorithms. The notation $n$ denotes the number of clients in each setting.
Refer to caption: x74.png

We can observe from Figure 11, det-DASHA with ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ and ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$ both outperform det-DASHA with ${\bm{D}}^{**}_{{\bm{I}}_{d}}$, which demonstrate the effectiveness of using a matrix stepsize instead of a scalar stepsize. However, depending on the parameters of the problem, it is hard to reach a general conclusion whether ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ is better than ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$ or not.

### I.12 Comparison of det-MARINA and det-DASHA

In this section, we aim to provide a comparison of det-DASHA and det-MARINA. They are similar as they are both variance reduced version of det-CGD. However, the variance reduction techniques that are utilized are different. For det-MARINA, it is based on MARINA, and it requires synchronization from time to time depending on a probability parameter $p$, while for det-DASHA it utilizes the momentum variance reduction technique which was also presented in DASHA, it does not need any synchronization at all. Notice that for a fair comparison, we implement the two algorithms so that they use the same sketch. We mainly focus on the communication complexity, i.e. the convergence with respect to the number of bits transferred. Throughout the experiment, $\lambda=0.9$ is fixed. For det-DASHA we pick $3$ different kinds of stepsizes ${\bm{D}}^{**}_{{\bm{I}}_{d}}$, ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ and ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$. For det-MARINA, we also pick three different kinds of stepsizes correspondingly ${\bm{D}}^{*}_{{\bm{I}}_{d}}$, ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ and ${\bm{D}}^{*}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$. We use the same sketch for all of the algorithms we are trying to compare.

Figure: Figure 12: Comparison of det-DASHA with three different stepsizes ${\bm{D}}^{**}_{{\bm{I}}_{d}}$, ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ and ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$, and det-MARINA with ${\bm{D}}^{*}_{{\bm{I}}_{d}}$, ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ and ${\bm{D}}^{*}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$ in terms of communication complexity. Throughout the experiment, $\lambda$ is fixed at $0.9$, the same Rand-$\tau$ sketch is used for all the algorithms. The notation $n$ denotes the number of clients in each setting. Each algorithm is run for a fixed number of iteration $K=5000$.
Refer to caption: x80.png

It is obvious from Figure 12 that det-DASHA always has a better communication complexity comparing to the det-MARINA counterpart. Notice that here since each algorithm is run for a fixed number of iterations, so $x$-axis actually records the total number of bytes transferred for each algorithm. For det-DASHA, ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ perform similarly to ${\bm{D}}^{**}_{\operatorname{diag}^{-1}\left({\bm{L}}\right)}$, and both are better than ${\bm{D}}^{**}_{{\bm{I}}_{d}}$. This is expected since the same sketch is used, and the number of bytes transferred in each iteration is the same for each variant of det-DASHA. The same relation also holds for det-MARINA.

### I.13 Comparison in terms of function values

In this section, we compare det-MARINA and det-DASHA in terms of function values. The starting points of the two algorithms are set to be the same, and we run the two algorithms for multiple times and we average the function values we obtained in each iteration. For the two algorithms, we use the same sketch, and since we are interested in the performance in terms of communication complexity, we use the number of bytes transferred in the training process as the $x$-axis. We run each of the algorithm for $20$ times, and fix $\lambda=0.9$. The starting point is fixed throughout the experiment. We pick ${\bm{D}}^{**}_{{\bm{L}}^{-1}}$ as the stepsize of det-DASHA, while ${\bm{D}}^{*}_{{\bm{L}}^{-1}}$ as the stepsize of det-MARINA.

Figure: Figure 13: Comparing the performance of det-DASHA with ${\bm{D}}^{**}{{\bm{L}}^{-1}}$ and det-MARINA with ${\bm{D}}^{*}{{\bm{L}}^{-1}}$ in terms of the decreasing function values. The function values for each algorithm represent an average of $20$ runs using different random seeds. Here, $\lambda=0.9$ is fixed throughout the experiment, and the starting point for the two algorithms in different runs is the same. The notation $n$ stands for the number of clients, and $p$ represents the probability used in det-MARINA. The same Rand-$\tau$ sketch is employed for both algorithms.
Refer to caption: x86.png

Observing Figure 13, we can see that the function values continuously decrease as the algorithms progress through more iterations.
However, the stability observed here differs from the case of the average (matrix) norm of gradients.
Our theoretical framework, as presented in this paper, primarily addresses the average norm of gradients in the non-convex case.
Despite this, the experiment reinforces the effectiveness of our algorithms, showcasing consistent decreases in function values.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: I have several questions about Section 4.2.

a) Line 310 says we can relax the problem to certain linear subspace and achieve Corollary 1. Is there any theoretical result to bound the approximation error of such relax?  

b) In the expression of equation (6), we require achieving matrices ${\bf W}^{1/2}$, ${\bf L}^{-1}$, and several eigenvalues. Is the computation of these quantity expensive in practice? 

c) Line 394 says we can estimate ${\bf L}_i$ by gradient method (Wang et al., 2022). I have not detailed checked the paper of Want et al. (2022). Can you summarize their main results and ideas (e.g., the complexity and the algorithm design)?

I have also some questions about the comparisons with related work.

d) I think more detailed introduction for MARINA and DASHA is necessary, which is useful to help the readers understand why the authors focus on the comparisons with these two methods.

e) The iteration with matrix stepsize looks more expensive than classical first-order methods. I think the comparison with classical first-order methods (e.g., distributed GD, SGD, SARAH...) should be involved, including the discussion on the computational cost.
soundness: 3
strengths: The idea of using matrix stepsize and variance-reduction to improve the communication complexity in federated learning is interesting. The thermotical guarantees have been established.
summary: This paper provided two novel federated learning optimization methods based on variance-reduced matrix stepsizes. The theoretical analysis shows the communication complexity of proposed methods are better than existing methods, which is also validated by numerical experiments.
weaknesses: Please see the part of questions.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: See weaknesses.
soundness: 3
strengths: 1) The proposed schemes address an issue that makes the distributed version of det-CGD less desirable. This is achieved via the incorporation of a variance reduction mechanisms that aids in controlling the variance introduced by the compression.
2) The work is well written and the contribution is clear.
summary: The work "Variance reduced distributed nonconvex optimization using matrix stepsizes" proposes decentralized variance reduced schemes with matrix step-sizes. The incorporation of variance reductions allows for convergence to arbitrary precision as opposed to the original work Li et. al (2024b). The authors provide theoretical guarantees for their algorithms and demonstrate their performance empirically.
weaknesses: 1) The novelty of this work is limited, as it seems combine a known algorithm (detCGD) with a given problem with two known solutions to the problem (proposed in MARINA and DASHA). It is unclear to this reviewer what is the main challenge when combining these. I suggest the authors highlight which technical challenges need to be addressed to incorporate variance reduction into detCGD, and why these significantly differ from the case in which the step-size is a scalar.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: (1) Regarding the optimal matrix stepsize in the optimization problem in line 304, the authors minimize the RHS of (5) with respect to D that aims to makes the bound as small as possible. However, the LHS of (5) also contains D. How could you tell that minimizing the RHS yields the optimal improvement?

(2) When comparing det-MARINA and and MARINA, in Corollary 3 for example, they say this iteration complexity of det-MARINA is better than MARINA. Please elaborate this claim more. For example, their gradient norm are not the same, how can we compare them?
soundness: 3
strengths: The paper is generally well-written with sufficient background. The motivation is also clear and natural: matrix stepsize is better than scalar stepsize. Convergence analysis is provided.
summary: The paper proposes matrix-stepsize versions, namely det-MARINA and det-DASHA, of distributed variance-reduction algorithms MARINA and DASHA. The main idea was that the scalar stepsize can be replaced by a matrix stepsize similar to second-order (Newton) methods. The authors derive convergence guarantees to stationarity in the sense that the gradient squared norm vanishes with rate $O(1/k)$ where $k$ is the algorithms' iterate. Numerical experiments show some improvement of the proposed methods compared to MARINA and DASHA.
weaknesses: Since the proposed methods are matrix-stepsize adaptations of MARINA and DASHA, we can expect them to retain the general properties observed in the original MARINA and DASHA schemes. To my understanding, both MARINA and DASHA also converge to stationarity and, unlike det-CGD, do not exhibit a restrictive convergence neighborhood. Thus, these new algorithms would naturally inherit this characteristic, which makes this advantage less novel in my view. From a theoretical standpoint, it would be helpful to directly compare Theorems 1 and 2 with the convergence properties of MARINA and DASHA (see specific questions below). Additionally, the Lipschitz matrix seems impractical in most cases, whereas the Lipschitz constant is typically more accessible; therefore, including the Lipschitz matrix may not be essential. While the experiments do show improvements, the gains appear relatively modest.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. What is the challenge in analysis of combining variance reduction with matrix step size?
2. How does the rate compare with existing literature?
3. Why one should take the trouble of using matrix step size?
soundness: 3
strengths: variance reduced method together with matrix step size shows convergence
summary: the paper considers matrix step size for federated learning, combined with some variance reduction methods to achieve exact convergence
weaknesses: only considered nonconvex case

Meta-review: This paper extends the det-CGD algorithm, which uses matrix stepsizes for non-convex optimization, by introducing two variance-reduced versions: det-MARINA and det-DASHA. These methods improve iteration and communication complexities over the existing methods, with both theoretical and empirical validation of their superior performance.

While the paper's theoretical approach is interesting, the main criticism raised by the reviewers was the lack of novelty in the proposed method. Specifically, the following weaknesses were highlighted:

1. One reviewer noted that the novelty of this work is limited, as it combines a known algorithm (det-CGD) with a given problem and two existing solutions (MARINA and DASHA). The reviewer found it unclear what specific challenges were addressed in combining these components.
2. Another reviewer expressed a similar concern, stating that the results of this paper are unsurprising when compared to MARINA and DASHA.

Decision: Reject