HERTA: A High-Efficiency and Rigorous Training Algorithm for Unfolded Graph Neural Networks

## Abstract

Abstract As a variant of Graph Neural Networks (GNNs), Unfolded GNNs offer enhanced interpretability and flexibility over traditional designs.
Nevertheless, they still suffer from scalability challenges when it comes to the training cost. Although many methods have been proposed to address the scalability issues, they mostly focus on per-iteration efficiency, without worst-case convergence guarantees. Moreover, those methods typically add components to or modify the original model, thus possibly breaking the interpretability of Unfolded GNNs. In this paper, we propose HERTA: a High-Efficiency and Rigorous Training Algorithm for Unfolded GNNs that accelerates the whole training process, achieving a nearly-linear time worst-case training guarantee. Crucially, HERTA converges to the optimum of the original model, thus preserving the interpretability of Unfolded GNNs. Additionally, as a byproduct of HERTA, we propose a new spectral sparsification method applicable to normalized and regularized graph Laplacians that ensures tighter bounds for our algorithm than existing spectral sparsifiers do. Experiments on real-world datasets verify the superiority of HERTA as well as its adaptability to various loss functions and optimizers.

## 1 Introduction

Graph Neural Networks (GNNs) have become a powerful modern tool for handling graph data because of their strong representational capabilities and ability to explore the relationships between data points (Kipf & Welling, ; Battaglia et al., ; Xu et al., ). Like many deep learning models, GNNs are generally designed by heuristics and experience, which makes analyzing and understanding them a difficult task. *Unfolded GNNs* (Yang et al., ) are a type of GNNs that are rigorously derived from an optimization problem, which makes their inductive bias and training dynamics easier to interpret.

Despite their enhanced interpretability, training unfolded models can be extremely expensive, especially when the graph is large and (relatively) dense. This issue can come from two aspects:

- 1)
Slow iterations: in every iteration the whole graph must be fed into the model, and the interactive nature of graphs prevents trivially utilizing techniques like mini-batch training;
- 2)
Slow convergence: the training convergence rate is related to the connectivity of the graph and the well-conditionedness of the node features, and the model converges slowly when the data is ill-conditioned.

Many methods have been proposed to address the high cost of training unfolded models (Chen et al., ; Chiang et al., ; Zeng et al., ; Ding et al., ; Fey et al., ; Jiang et al., ), primarily by using graph sampling to enable mini-batch training schemes that reduce the per-iteration cost (Issue 1).
However, these approaches typically require distorting the underlying optimization objective explicitly or implicitly, thus diminishing the rigorous and interpretable nature of Unfolded GNNs. Moreover, existing works have mostly focused on per-iteration efficiency, while the convergence rate of the optimization in training Unfolded GNNs (Issue 2) remains un-addressed, leading to methods that are not robust to ill-conditioned data.

We propose HERTA: a High-Efficiency and Rigorous Training Algorithm for Unfolded GNNs, which is an algorithmic framework designed to address both Issues 1 and 2, while at the same time preserving the rigorous and interpretable nature of Unfolded GNNs.
Unlike many existing methods that require changing the GNN model or defining a new objective, HERTA converges to the optimum of the original training problem, and thus preserves the interpretability of the unfolded models. Moreover, HERTA uses a specialized preconditioner to ensure fast linear convergence to the optimum, requiring only a logarithmic number of passes over the data. We show empirically that HERTA can be used to accelerate training for a variety of GNN objectives, as an extension of popular optimization techniques.

A key novelty of our algorithm lies in the construction of a preconditioner for the GNN objective, which accelerates the convergence rate. To construct this preconditioner, we design a spectral sparsification method for approximating the squared inverse of a normalized and regularized Laplacian matrix, by relying on an extension of the notion of effective resistance to overcome the limitations of existing graph sparsification tools. Below, we present an informal statement of our main result.

###### Theorem 1.1 (Informal version of Theorem 5.1 ) .

HERTA solves the $\lambda$-regularized Unfolded GNN objective (4.2) with $n$ nodes, $m$ edges and $d$-dimensional node features to within accuracy $\epsilon$ in time $\tilde{O}\left((m+nd)\left(\log\frac{1}{\epsilon}\right)^{2}+d^{3}\right)$ as long as the number of large eigenvalues of the graph Laplacian is $O(n/\lambda^{2})$.

In practice, the node feature dimensionality $d$ is generally much smaller than the graph size $m$, in which case, the running time of HERTA is $\tilde{O}\left((m+nd)\left(\log\frac{1}{\epsilon}\right)^{2}\right)$. Notice that to describe a typical graph dataset with $m$ edges, $n$ nodes and $d$-dimensional node features, we need at least $O(m+nd)$ float or integer numbers. This shows the essential optimality of HERTA: its running time is of the same magnitude as reading the input data (up to logarithm factors).

The condition that the graph Laplacian has $O(n/\lambda^{2})$ large eigenvalues is not strictly necessary; it is used here to simplify the time complexity (see Theorem 5.1 for the full statement). This condition is also not particularly restrictive, since in most practical settings the GNN parameter $\lambda$ is chosen as an absolute constant. At the same time, it is natural that the complexity of training a GNN depends on the number of large eigenvalues in the graph Laplacian, since this quantity can be interpreted as the effective dimension in the Laplacian energy function (4.1) used to define the Unfolded GNN (see more discussion in Section 5).

Outline. The paper is organized as follows. In Section 2 we introduce related work. In Section 3 we introduce the mathematical notation and concepts used in this paper. We define our problem setting in Section 4. In Section 5 we present and analyze our algorithm HERTA, and introduce the techniques that are used in proving the main result.
We conduct experiments on real-world datasets and show the results in Section 6. Finally, we conclude with some potential future directions in Section 7.

###### Theorem 1.1 (Informal version of Theorem 5.1 ) .

## 2 Related Work

Unfolded Graph Neural Networks.  Unlike conventional GNN models, which are designed mainly by heuristics, the forward layers of Unfolded GNNs are derived by explicitly optimizing a graph-regularized target function (Zhang et al., ; Yang et al., ; Ma et al., ; Liu et al., ; Yang et al., ; Zhu et al., ; Ahn et al., ). This optimization nature of Unfolded GNNs allows developing models in a more interpretable and controllable way. The presence of an explicit optimization problem allows designers to inject desired properties into the model and better understand it. This approach has been used to overcome known GNN issues such as over-smoothing (Oono & Suzuki, ; Liu et al., ; Yang et al., ) and sensitivity towards spurious edges (Zügner et al., ; Zhu et al., ).

Efficient Graph Neural Network Training.  In order to address the scalability issue of GNNs, many techniques have been proposed. There are two major directions in the exploration of the solutions of this issue. The first direction is to adopt sampling methods. These methods include node-wise sampling (Hamilton et al., ), layer-wise sampling (Chen et al., ) and subgraph-wise sampling (Chiang et al., ; Jiang et al., ). The second direction is to utilize embeddings from past iterations to enhance current iterations and to obtain more accurate representation with mini-batches and fewer forward layers (Fey et al., ; Ding et al., ; Xue et al., ).

As we indicated in Section 1, all of the above
methods aim at addressing Issue 1 (Slow Iterations), and none of them achieves a worst-case convergence guarantee without affecting the underlying GNN model. In contrast, HERTA addresses both Issue 1 and also Issue 2 (Slow Convergence), and possesses a theoretical guarantee on the training convergence rate while preserving the original target model.

Matrix Sketching and Subsampling.
Our techniques and analysis are closely related to matrix sketching and subsampling, which are primarily studied in the area of Randomized Numerical Linear Algebra (RandNLA, Woodruff et al. ; Drineas & Mahoney ; Martinsson & Tropp ; Murray et al. ). Given a large matrix, by sketching or sampling we compress it to a much smaller matrix with certain properties preserved, thus accelerating the solution by conducting expensive computations on the smaller matrix. This type of methods lead to improved randomized algorithms for tasks including low rank approximation (Halko et al., ; Cohen et al., ; Clarkson & Woodruff, ), linear system solving (Peng & Vempala, ; Dereziński & Yang, ), least squares regression (Rokhlin & Tygert, ; Meng et al., ) and so on (Cohen et al., ). Our usage of these methods includes constructing spectral sparsifiers of graphs, obtained by edge subsampling (Spielman & Srivastava, ; Spielman & Teng, ; Koutis et al., ), which are central in designing fast solvers for Laplacian linear systems (Vishnoi et al., ; Kelner et al., ; Peng & Spielman, ).

## 3 Preliminaries

For a vector ${\bm{x}}$, we denote its $\ell_{2}$ norm as $\|{\bm{x}}\|$. We use $\mathrm{diag}(\bm{x})$ to denote diagonalization of $\bm{x}$ into a matrix. For a matrix ${\bm{M}}$, we use $\|{\bm{M}}\|_{\mathcal{F}}$ and $\|\bm{M}\|$ to denote its Frobinius norm and operator norm. We also denote its largest and smallest singular value as $\sigma_{\max}({\bm{M}})$ and $\sigma_{\min}({\bm{M}})$ respectively, and denote its condition number as $\kappa({\bm{M}})=\frac{\sigma_{\max}({\bm{M}})}{\sigma_{\min}({\bm{M}})}$. We use $\mathrm{nnz}(\bm{M})$ to denote the number of non-zero entries of $\bm{M}$. For two positive semidefinite (PSD) matrices ${\bm{\Sigma}}$ and $\widetilde{\bm{\Sigma}}$, if there exists $\epsilon\in(0,1)$ such that

$$ $\displaystyle(1-\epsilon)\widetilde{\bm{\Sigma}}\preceq{\bm{\Sigma}}\preceq(1+ \epsilon)\widetilde{{\bm{\Sigma}}},$ (3.1) $$

then we say ${\bm{\Sigma}}\approx_{\epsilon}\widetilde{{\bm{\Sigma}}}$, where $\preceq$ refers to Loewner order.

We use ${\bm{1}}$ to denote a vector with all entries equal to $1$, and the dimension of which is determined by the context if not specified. We use ${\bm{\delta}}_{u}$ to denote a vector with $u$-th entry equal to $1$ and all other entries equal to 0 0. We also use $\delta_{u,v}$ to represent the $v$-th entry of ${\bm{\delta}}_{u}$, i.e.:

$$ $\delta_{u,v}=\begin{cases}1&u=v\\ 0&u\neq v\end{cases}.$ $$

For a function $\ell:\mathbb{R}^{d}\to\mathbb{R}$ that is bounded from bottom, a point ${\bm{w}}_{0}$ is called a solution of $\ell$ with $\epsilon$ error rate if it satisfies $\ell({\bm{w}}_{0})\leq(1+\epsilon)\ell^{*}$ where $\ell^{*}=\inf_{{\bm{w}}}\ell({\bm{w}})$.

We adopt big-O notations in the time complexity analysis. Moreover, we use the notation $\tilde{O}(\cdot)$ to hide polynomial logarithmic factors of the problem size. As an example, $O(\log(n))$ can be expressed as $\tilde{O}(1)$.

For a matrix ${\bm{\Pi}}\in\mathbb{R}^{s\times n}$, we call it a Gaussian sketch if its entries are i.i.d. $\mathcal{N}\left(0,1/s\right)$ random variables. For a diagonal matrix ${\bm{R}}\in\mathbb{R}^{n\times n}$, we call it a diagonal Rademacher matrix if its entries are i.i.d. Rademacher random variables. See Appendix A for more detailed introductions.

## 4 Problem Setting

Throughout this paper, we consider a graph $\mathcal{G}$ with $n$ nodes and $m$ edges. We denote its adjacency matrix by ${\bm{A}}\in\mathbb{R}^{n\times n}$, degree matrix by ${\bm{D}}=\mathop{\mathrm{diag}}\left({\bm{A}}{\bm{1}}\right)\in\mathbb{R}^{n
\times n}$, and Laplace matrix (or Laplacian) by ${\bm{L}}={\bm{D}}-{\bm{A}}\in\mathbb{R}^{n\times n}$. We also use the normalized adjacency $\hat{\bm{A}}={\bm{D}}^{-1/2}{\bm{A}}{\bm{D}}^{-1/2}$ and normalized Laplacian $\hat{\bm{L}}={\bm{I}}-\hat{\bm{A}}$. Notice that we assume there’s a self-loop for each node to avoid 0 0-degree nodes, which is commonly ensured in GNN implementations (Kipf & Welling, ; Hamilton et al., ), and therefore ${\bm{D}}^{-1/2}$ always exists. We also use ${\bm{B}}\in\mathbb{R}^{m\times n}$ to represent the incidence matrix of graph $\mathcal{G}$, that is, if the $i$-th edge in $\mathcal{G}$ is $(u_{i},v_{i})$, then the $i$-th row of ${\bm{B}}$ is ${\bm{\delta}}_{u_{i}}-{\bm{\delta}}_{v_{i}}$ (the order of $u_{i}$ and $v_{i}$ is arbitrary). It is not hard to see that ${\bm{L}}={\bm{B}}^{\top}{\bm{B}}$.

We also assume that for each node $u$ in the graph, there is a feature vector ${\bm{x}}_{u}\in\mathbb{R}^{d}$, and a label vector(^1^11We allow the label to be a vector to make the concept more general. When labels are one-hot vectors, the task is a classification task, and when labels are scalars, it is a regression task.) ${\bm{y}}_{u}\in\mathbb{R}^{c}$ attached with it. Let ${\bm{X}}\in\mathbb{R}^{n\times d}$ and ${\bm{Y}}\in\mathbb{R}^{n\times c}$ be the stack of ${\bm{x}}_{i}$-s and ${\bm{y}}_{i}$-s. We assume $n\geq\max\{d,c\}$ and ${\bm{X}}$ and ${\bm{Y}}$ are both column full rank.

The Unfolded GNN we consider in this paper is based on TWIRLS (Yang et al., ), which
is defined by optimizing the following objective (called the “energy function”):

$$ $E({\bm{Z}};{\bm{W}}):=\frac{\lambda}{2}\mathrm{Tr}\left({\bm{Z}}^{\top}\hat{ \bm{L}}{\bm{Z}}\right)+\frac{1}{2}\left\|{\bm{Z}}-f({\bm{X}};{\bm{W}})\right\| _{\mathcal{F}}^{2}.$ (4.1) $$

Here ${\bm{Z}}\in\mathbb{R}^{n\times c}$ is the optimization variable that can be viewed as the hidden state of the model, ${\bm{W}}$ is the learnable parameter of the model, and $\lambda>0$ is a hyper-parameter to balance the graph structure information and node feature information. Note that $f:\mathbb{R}^{d}\to\mathbb{R}^{c}$ is a function of $\bm{X}$ and is parameterized by matrix ${\bm{W}}$.

For the model training, given a loss function $\mathcal{\ell}:\mathbb{R}^{n\times c}\times\mathbb{R}^{n\times c}\to\mathbb{R}$, the training of TWIRLS can be represented by the following bi-level optimization problem

$$ $\displaystyle{\bm{W}}^{*}\in\arg\min_{{\bm{W}}\in\mathbb{R}^{d\times c}}\ell \left[{\bm{Z}}^{*}({\bm{W}});{\bm{Y}}\right]$ (4.2) $\displaystyle\text{s.t.}\quad{\bm{Z}}^{*}({\bm{W}})\in\arg\min_{{\bm{Z}}\in \mathbb{R}^{n\times c}}E({\bm{Z}};{\bm{W}}),$ (4.3) $$

where the solution of inner problem eq. 4.3 (i.e. ${\bm{Z}}^{*}({\bm{W}})$) is viewed as the model output. The outer problem eq. 4.2 is viewed as the problem of training.

In (Yang et al., ), the forward process of TWIRLS is obtained by unfolding the gradient descent algorithm for minimizing $E$ given a fixed ${\bm{W}}$:

$$ $\displaystyle{\bm{Z}}^{(t+1)}$ $\displaystyle={\bm{Z}}^{(t)}-\alpha\nabla E\left({\bm{Z}}^{(t)}\right)$ (4.4) $\displaystyle=(1-\alpha\lambda-\alpha){\bm{Z}}^{(t)}+\alpha\hat{\bm{A}}{\bm{Z} }^{(t)}+\alpha f\left({\bm{X}};{\bm{W}}\right),$ (4.5) $$

where $\alpha>0$ is the step size of gradient descent. Notice that we fix ${\bm{W}}$ and view $E$ as a function of ${\bm{Z}}$. The model output is ${\bm{Z}}^{(T)}$ for a fixed iteration number $T$. When $\alpha$ is set to a suitable value and $T$ is large, it is clear that ${\bm{Z}}^{(T)}$ is an approximation of ${\bm{Z}}^{*}$ defined in eq. 4.3.

In this paper, we always assume $T\to\infty$ and thus the output of TWIRLS is exactly the unique global optimum of eq. 4.1. Therefore, we directly consider the explicit optimal solution of eq. 4.1:

$$ $\displaystyle{\bm{Z}}^{*}({\bm{W}})=\arg\min_{{\bm{Z}}}E({\bm{Z}})=\left({\bm{ I}}+\lambda\hat{\bm{L}}\right)^{-1}f({\bm{X}};{\bm{W}}).$ (4.6) $$

In (Yang et al., ), the implementation of $f$ is arbitrary as long as it can be trained (i.e. computable, smooth, etc.). In this paper we consider a simple implementation of $f$ which is a linear function, namely $f({\bm{X}};{\bm{W}})={\bm{X}}{\bm{W}}$ where ${\bm{W}}\in\mathbb{R}^{d\times c}$. Notice that this is also the practical setting used in (Yang et al., ) for many datasets. For loss function we choose the mean squared error (MSE) loss defined as $\ell({\bm{Z}},{\bm{Y}}):=\frac{1}{2}\left\|{\bm{Z}}-{\bm{Y}}\right\|_{\mathcal
{F}}^{2}$.

Figure: Algorithm 1 HERTA: A High-Efficiency and Rigorous Training Algorithm for TWIRLS

By combining the above formulations and assumptions, the overall target function we consider in this paper has the following form:

$$ $\displaystyle\ell\left[{\bm{Z}}^{*}({\bm{W}});{\bm{Y}}\right]$ $\displaystyle=\frac{1}{2}\left\|\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{-1} {\bm{X}}{\bm{W}}-{\bm{Y}}\right\|_{\mathcal{F}}^{2}.$ (4.7) $$

This loss can be viewed as a summation of $c$ independent sub-loss functions. Define $\ell_{i}({\bm{w}}_{i}):=\frac{1}{2}\|({\bm{I}}+\lambda\hat{\bm{L}})^{-1}{\bm{X
}}{\bm{w}}_{i}-{\bm{y}}_{i}\|^{2}$ where ${\bm{w}}_{i}$ and ${\bm{y}}_{i}$ are the $i$-th column of ${\bm{W}}$ and ${\bm{Y}}$ respectively, then we have

$$ $\displaystyle\ell\left[{\bm{Z}}^{*}({\bm{W}});{\bm{Y}}\right]=\sum_{i=1}^{c} \ell_{i}({\bm{w}}_{i}).$ (4.8) $$

Note that we view $\ell_{i}$ as a function of ${\bm{w}}_{i}$. Moreover, by denoting $\hat{\bm{y}}_{i}=({\bm{I}}+\lambda\hat{\bm{L}}){\bm{y}}_{i}$ we can rewrite $\ell_{i}$ as

$$ $\displaystyle\ell_{i}({\bm{w}}_{i})=\frac{1}{2}\left\|\left({\bm{I}}+\lambda \hat{\bm{L}}\right)^{-1}\left({\bm{X}}{\bm{w}}_{i}-\hat{\bm{y}}_{i}\right) \right\|^{2}.$ (4.9) $$

Figure: Algorithm 2 $\mathsf{Sparsify}_{\epsilon}(\hat{\bm{L}})$: Regularized Spectral Sparsifier

## 5 Algorithm and Analysis

In this section we present our main algorithm HERTA, give the convergence result and analyze its time complexity. We also introduce our key techniques and give a proof sketch for the main result. To start we introduce the following notion.

###### Definition 5.1 (Effective Laplacian dimension) .

For normalized Laplacian $\hat{\bm{L}}$ and regularization term $\lambda>0$, define

$$ $\displaystyle n_{\lambda}:=\mathrm{Tr}\left[\hat{\bm{L}}\left(\hat{\bm{L}}+ \lambda^{-1}{\bm{I}}\right)^{-1}\right]$ (5.1) $$

as the effective Laplacian dimension of the graph.

Denote $\{\lambda_{i}\}_{i=1}^{n}$ as the eigenvalues of $\hat{\bm{L}}$. Then the effective Laplacian dimension can be written as $n_{\lambda}=\sum_{i}\frac{\lambda_{i}}{\lambda_{i}+\lambda^{-1}}$. We can see that $n_{\lambda}$ is roughly the number of eigenvalues of $\hat{\bm{L}}$ which are of order $\Omega(\lambda^{-1})$, i.e., the number of “large eigenvalues” mentioned in Theorem 1.1. It is not hard to see that $n_{\lambda}\leq n$ and $n_{\lambda}\to n$ as $\lambda\to\infty$.

Intuitively, $n_{\lambda}$ represents the number of eigen-directions for which the Laplacian regularizer $\frac{\lambda}{2}\mathrm{Tr}({\bm{Z}}^{\top}\hat{\bm{L}}{\bm{Z}})$ is large, thus having a significant impact on the energy function (4.1). For the remaining eigen-directions which are less significant, the effect of the regularizer is minimal. This can also be seen from the loss function (4.9): when we apply matrix $({\bm{I}}+\lambda\hat{\bm{L}})^{-1}$ to a vector, the effect of the Laplacian is dominated by the ${\bm{I}}$ term if this vector is aligned with a “small” Laplacian eigen-direction (i.e., with an eigenvalue significantly smaller than $\lambda^{-1}$), and thus the objective is close to the usual unregularized least squares loss.

Based on the above discussion, it stands to reason that the larger the effective Laplacian dimension $n_{\lambda}$ of the graph, the more work we must do during the training to preserve the graph structure of the GNN. With the above discussion, we formally give our main result for HERTA.

###### Theorem 5.1 (Main result) .

For any $\epsilon>0$, with a proper step size $\eta$, number of iterations $T$ and constant $K>0$, Algorithm 1 finds a solution $\hat{\bm{W}}\in\mathbb{R}^{d\times c}$ in time

$$ $\displaystyle\tilde{O}\left(\left(m+nd\right)c\left(\log\frac{1}{\epsilon} \right)^{2}+n_{\lambda}\lambda^{2}d+d^{3}\right)$ (5.2) $$

such that with probability $1-\frac{1}{n}$,

$$ $\displaystyle\ell\left[{\bm{Z}}^{*}({\bm{\hat{W}}});{\bm{Y}}\right]\leq(1+ \epsilon)\cdot\ell^{*},$ (5.3) $$

where $\ell^{*}:=\min_{{\bm{W}}}\ell\left[{\bm{Z}}^{*}\left({\bm{W}};{\bm{Y}}\right)\right]$.

Notice that if we assume $n_{\lambda}=O(n/\lambda^{2})$, then $\tilde{O}(n_{\lambda}\lambda^{2}d)=\tilde{O}(nd)$. By further assuming the number of classes $c=O(1)$, we recover the bound given by Theorem 1.1.

Before heading into the details of HERTA (Algorithm 1), we first analyze the original (standard) implementation of training TWIRLS and point out the efficiency bottlenecks. We then introduce the key techniques and methods we need to prove our main result. Finally we present our main algorithm HERTA, and give a proof sketch for Theorem 5.1.

###### Definition 5.1 (Effective Laplacian dimension) .

###### Theorem 5.1 (Main result) .

### 5.1 Analysis of TWIRLS Training

In this subsection, we introduce the framework that we use to analyze the complexity of TWIRLS training. We also provide a complexity analysis of the implementation used in (Yang et al., ).

As noted in Section 4, the whole optimization problem can be viewed as a bi-level optimization problem, where the inner-problem eq. 4.3 approximates the linear system solver $\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{-1}({\bm{X}}{\bm{w}}_{i}-\hat{\bm{y
}}_{i})$, and the outer problem eq. 4.2 uses the linear system solution to approximate the optimal solution of the training loss eq. 4.7. When we discuss the solution of outer problems, we treat the solver of the inner problem (linear solver) as a black box. We formally define a linear solver as follows.

###### Definition 5.2 (Linear solver) .

For a positive definite matrix ${\bm{H}}\in\mathbb{R}^{n\times n}$ and a real number $\epsilon>0$, we call $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ a linear solver for ${\bm{H}}$ with $\epsilon$ error rate if it satisfies

$$ $\displaystyle\forall{\bm{u}}\in\mathbb{R}^{n},~{}\left\|f({\bm{u}})-{\bm{H}}^{ -1}{\bm{u}}\right\|_{{\bm{H}}}\leq\epsilon\left\|{\bm{H}}^{-1}{\bm{u}}\right\| _{{\bm{H}}}.$ (5.4) $$

In (Yang et al., ), the outer problem is also solved by standard gradient descent. The gradient of $\ell_{i}\left({\bm{w}}_{i}\right)$ is

$$ $\displaystyle\nabla\ell_{i}({\bm{w}}_{i})={\bm{X}}^{\top}\left({\bm{I}}+ \lambda\hat{\bm{L}}\right)^{-2}\left({\bm{X}}{\bm{w}}_{i}-\hat{\bm{y}}_{i} \right).$ (5.5) $$

In actual implementation, the matrix inverse is calculated by linear solvers. Suppose that we have a linear solver for $({\bm{I}}+\lambda\hat{\bm{L}})$ denoted by $\mathcal{S}$, then by embedding it into eq. 5.5 we can obtain the gradient approximation as

$$ $\displaystyle\widetilde{\nabla\ell_{i}({\bm{w}}_{i})}:={\bm{X}}^{\top}\mathcal {S}\left[\mathcal{S}\left({\bm{X}}{\bm{w}}_{i}-\hat{\bm{y}}_{i}\right)\right].$ (5.6) $$

The following convergence result can be derived for this approximate gradient descent. See Appendix for the proof.

###### Lemma 5.1 (Convergence) .

If we minimize eq. 4.9 using gradient descent with the gradient approximation defined in eq. 5.6, where $\mathcal{S}$ is a linear solver of ${\bm{I}}+\lambda\hat{\bm{L}}$ with $\mu$ error rate satisfying $\mu\leq\min\left\{\frac{\epsilon^{1/2}}{50\kappa({\bm{X}})\lambda^{2}},1\right\}$, then, to obtain an $\epsilon\in(0,1)$ error rate the outer problem eq. 4.9, the number of iterations needed is

$$ $\displaystyle T=O\left(\kappa\log\frac{1}{\epsilon}\right),$ (5.7) $$

where $\kappa:=\kappa({\bm{X}}^{\top}(\bm{I}+\lambda\hat{\bm{L}})^{-2}{\bm{X}})$.

Based on Lemma 5.1, in Appendix B we provide a detailed analysis of the time complexity of the implementation used by (Yang et al., ). From this analysis, one can figure out two bottlenecks of the computational complexity for solving TWIRLS: 1) The number of iterations needed for the inner loop is dependent on $\lambda$; 2) The number of iterations needed for the outer loop is dependent on the condition number of the outer problem. Our acceleration algorithm is based on the following two observations, which correspond to the aforementioned bottlenecks:

- 1.
The data matrix in the inner loop is a graph Laplacian plus identity. It is possible to solve the linear system defined by this matrix faster by exploiting its structure.
- 2.
A good enough preconditioner can be constructed for the outer iteration using techniques from RandNLA.

###### Definition 5.2 (Linear solver) .

###### Lemma 5.1 (Convergence) .

### 5.2 Key Techniques

In this subsection, we introduce the mathematical tools that we will use in the proof of our main result.

#### SDD Solvers.

If a symmetric matrix ${\bm{H}}\in\mathbb{R}^{n\times n}$ satisfies

$$ $\displaystyle\forall 1\leq k\leq n,~{}~{}{\bm{H}}_{k,k}\geq\sum_{i=1}^{n}\left |{\bm{H}}_{k,i}\right|,$ (5.8) $$

it is called a symmetric diagonal dominated (SDD) matrix. It has been shown that the linear system defined by sparse SDD matrices can be solved in an efficient way, via a SDD solver (Kelner et al., ; Peng & Spielman, ).

Since ${\bm{I}}+\lambda\hat{\bm{L}}$ is indeed a SDD matrix in our problem, we can apply off-the-shelf SDD solvers. In the following, we use $\mathsf{SDDSolver}_{\epsilon}({\bm{H}},{\bm{u}})$ to denote the SDD solver that calculates ${\bm{H}}^{-1}{\bm{u}}$ with error rate less or equal to $\epsilon$(^2^22Principally, $\mathsf{SDDSolver}({\bm{H}};\cdot)$ is a vector function, but for convenience we sometimes also apply it to matrices, in which case we mean applying it column-wisely. We also use the same convention for other vector functions.). From Lemma A.1, the time complexity of calculating $\mathsf{SDDSolver}_{\epsilon}({\bm{H}};{\bm{u}})$ is $\tilde{O}\left[\mathrm{nnz}({\bm{H}})\right]$. See Section A.2 for details.

#### Fast Matrix Multiplication.

For a matrix ${\bm{Q}}\in\mathbb{R}^{n\times d}$ where $n>d$, it is known that calculating ${\bm{Q}}^{\top}{\bm{Q}}$ takes $O(nd^{2})$ time if implemented by brute force. In HERTA, we manage to achieve a faster calculation of ${\bm{Q}}^{\top}{\bm{Q}}$ through Subsampled Randomized Hadamard Transformation (SRHT) (Tropp, ; Woodruff et al., ). The key idea of SRHT is to first apply a randomized Hadamard transformation to the matrix to make the “information” evenly distributed in each row, then uniformly sample the rows of ${\bm{Q}}$ to reduce the dimension. In the following, we use $\mathsf{Hadamard}$ to represent the Hadamard transformation, see Section A.3 for the definition of $\mathsf{Hadamard}$ as well as a detailed introduction.

### 5.3 Regularized Spectral Sparsifier

It has been shown that for any connected graph with $n$ nodes, (no matter how dense it is) there is always a sparsified graph with $\tilde{O}\left(\frac{n}{\epsilon^{2}}\right)$ edges whose Laplacian is an $\epsilon$-spectral approximation of the original graph Laplacian (Spielman & Teng, ). Although the existence is guaranteed, it is generally hard to construct such a sparsified graph. An algorithm that finds such a sparsified graph
is called a spectral sparsifier.

While there has been many existing explorations on constructing spectral sparsifier of a given graph with high probability and tolerable time complexity (Spielman & Srivastava, ; Spielman & Teng, ; Koutis et al., ), the setting considered in this paper is different from the standard one: instead of graph Laplacian $\bm{L}$, the matrix we concern is the normalized and regularized Laplacian $\hat{\bm{L}}+\lambda^{-1}{\bm{I}}$.

At first glance, a spectral sparsifier also works for our problem since $\tilde{{\bm{L}}}\approx_{\epsilon}\hat{\bm{L}}$ implies $\tilde{{\bm{L}}}+\lambda^{-1}{\bm{I}}\approx_{\epsilon}\hat{\bm{L}}+\lambda^{-
1}{\bm{I}}$. However, it turns out that it is possible to make the bound even tighter if we take into account the regularization term $\lambda^{-1}{\bm{I}}$ here. A similar argument is also made in (Calandriello et al., ). We note that the setting explored in (Calandriello et al., ) is also different from ours: we consider the normalized Laplacian of a graph, which is not a Laplacian of any graph (even allowing weighted edges).

To this end, we propose a new spectral sparsifier that works for $\hat{\bm{L}}+\lambda^{-1}{\bm{I}}$ which is used in our problem. This spectral sparsifier reduces the number of edges to $O\left(\frac{n_{\lambda}}{\epsilon^{2}}\log n\right)$, which, as we discussed before, is smaller than $O\left(\frac{n}{\epsilon^{2}}\log n\right)$ obtained by standard spectral sparsifier.

###### Lemma 5.2 (Regularized spectral sparsifier) .

Let $\hat{\bm{L}}$ be a (normalized) graph Laplacian with $m$ edges and $n$ nodes. There exists a constant $C$ such that for any $\epsilon>0$ and $\lambda>0$, Algorithm 2 outputs
a (normalized) graph Laplacian $\tilde{\bm{L}}$ with $n$ nodes and $O\left(\frac{n_{\lambda}}{\epsilon^{2}}\log n\right)$ edges in time $\tilde{O}(m)$, such that

$$ $\displaystyle\tilde{\bm{L}}+\lambda^{-1}{\bm{I}}\approx_{\epsilon}\hat{\bm{L}} +\lambda^{-1}{\bm{I}}$ (5.9) $$

holds with probability at least $1-\frac{1}{2n}$.

The basic idea behind Lemma 5.2 is to use ridge leverage score sampling for the normalized incidence matrix ${\bm{B}}{\bm{D}}^{-1/2}$ (Cohen et al., ), which reduces the problem to estimating a regularized version of the effective resistance. We modify the existing methods of estimating effective resistance (Vishnoi et al., ) to make it work for the normalized and regularized Laplacian.

###### Lemma 5.2 (Regularized spectral sparsifier) .

### 5.4 Main Algorithm

Now we are ready to present our main algorithm HERTA, see Algorithm 1. As mentioned before, HERTA is composed of two major components: constructing a preconditioner matrix ${\bm{P}}$ and applying it to the optimization problem. Below we introduce each part.

#### Constructing the Preconditioner.

For the outer problem which can be ill-conditioned, we first precondition it using a preconditioner ${\bm{P}}\in\mathbb{R}^{d\times d}$ which is a constant level spectral approximation of the Hessian, i.e. ${\bm{X}}^{\top}\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{-2}{\bm{X}}$. We claim that the matrix ${\bm{P}}$ constructed inside Algorithm 1 indeed satisfies this property.

###### Lemma 5.3 (Preconditioner) .

Let ${\bm{P}}$ be the matrix constructed in Algorithm 1. With probability at least $1-\frac{1}{n}$,

$$ $\displaystyle{\bm{P}}\approx_{\frac{1}{2}}{\bm{X}}^{\top}\left({\bm{I}}+ \lambda\hat{\bm{L}}\right)^{-2}{\bm{X}}.$ (5.10) $$

Figure: Figure 1: The training loss comparison between HERTA and standard optimizers on MSE loss with $\lambda=1$. Dataset used from left to right: ogbn-arxiv, citeseer, pubmed.
Refer to caption: extracted/2403.18142v1/figures-log/ogbn-arxiv-1-fast,mse.png

###### Lemma 5.3 (Preconditioner) .

#### Solving the Outer Problem.

After obtaining the preconditioner ${\bm{P}}$ which approximates the Hessian by a constant level, we use it to precondition the outer problem and provably reduce the condition number to a constant. With the new problem being well-conditioned, iterative methods (like gradient descent) take much less steps to converge.

###### Lemma 5.4 (Well-conditioned Hessian) .

Let $\ell^{\prime}$ be such that $\ell^{\prime}({\bm{w}}_{i})=\ell({\bm{P}}^{-1/2}{\bm{w}}_{i})$. Suppose for some constant $c_{0}\in(0,1)$ we have ${\bm{P}}\approx_{c_{0}}{\bm{X}}^{\top}\left({\bm{I}}+\lambda\hat{\bm{L}}\right
)^{-2}{\bm{X}}$, then we can bound the condition number of Hessian of $\ell^{\prime}$ as

$$ $\displaystyle\kappa\left(\nabla^{2}\ell^{\prime}\right)\leq(1+c_{0})^{2}.$ (5.11) $$

Moreover, if ${\bm{w}}^{\prime}$ is a solution of $\ell^{\prime}$ with $\epsilon$ error rate, then ${\bm{P}}^{-1/2}{\bm{w}}^{\prime}$ is a solution of $\ell$ with $\epsilon$ error rate.

Notice that, the problem $\ell^{\prime}$ defined in Lemma 5.4 can be viewed as the original problem $\ell$ with ${\bm{X}}$ being replaced by ${\bm{X}}{\bm{P}}^{-1/2}$, therefore we can use Lemma 5.1 and obtain the convergence rate of HERTA. With the convergence result and analysis of the running time for each step, we are able to prove Theorem 5.1. See Section C.7 for the full proof.

Remark on the Unavoidable $\lambda^{2}$ in Runtime.  Notice that $\lambda^{2}$ appears in the time bound given in Theorem 5.1. From the proof in Section C.5 we can see that this term originates from the squared inverse of ${\bm{I}}+\lambda\hat{\bm{L}}$. In Section D.2, we show that in the worst case even if we approximate a matrix to a very high precision, the approximation rate of the corresponding squared version can still be worse by a factor of the condition number of the matrix, which in our case leads to the unavoidable $\lambda^{2}$ in the runtime.

Remark on Graph Sparsifying in Each Iteration.  Note that in Algorithm 1, we use the complete Laplacian $\hat{\bm{L}}$ for gradient iterations (i.e., the for-loop), and the graph sampling only occurs when constructing the preconditioner ${\bm{P}}$. We note that sampling $\hat{\bm{L}}$ in the for-loop can lead to extra loss in the gradient estimation, which forces us to sparsify the graph to a very high precision to ensure an accurate enough gradient for convergence. This could result in a suboptimal running time. Moreover, as discussed in Section 1, the current running time bound for HERTA is optimal up to logarithmic factors, which means that there is little to gain from performing extra sampling in for-loops. See Section D.3 for a more detailed and quantitative discussion.

###### Lemma 5.4 (Well-conditioned Hessian) .

## 6 Experiments

In this section, we verify our theoretical results through experiments on real world datasets. Since for the inner problem we use off-the-shelf SDD solvers, in this section we focus on the outer problem, i.e. the training loss convergence rate. In each setting we compare the training loss convergence rate of TWIRLS trained by our method against that trained by standard gradient descent using exactly the same training hyper-parameters.

Datasets. We conduct experiments on all of the datasets used in the Section 5.1 of (Yang et al., ): Cora, Citeseer and Pubmed collected by (Sen et al., ), as well as the ogbn-arxiv dataset from the Open Large Graph Benchmark (OGB, Hu et al., ).

Notice that for Cora, Citeseer and Pubmed, common practice uses the semi-supervised setting where there is a small training set and relatively large validation set. As we are comparing training convergence rate, we find it more comparative to use a larger training set (which makes solving the optimization problem of training more difficult). Therefore, we randomly select 80% of nodes as training set for Cora, Citeseer and Pubmed. For OGB, we use the standard training split.

Due to limited space, we only include the results with $\lambda=1$ and with datasets Citeseer, Pubmed and ogbn-arxiv in this section, and delay additional results to Appendix F.

### 6.1 Convergence Rate Comparison Under MSE Loss

In this section, we compare the convergence rate for models trained with MSE loss, which is well aligned with the setting used to derive our theory. We also adopt a variation of HERTA that allows using it on other optimizers and apply it on Adam optimizer. See Figure 1 for the results. Notice that for each figure, we shift the curve by the minimum value of the loss(^3^33The minimum value of the loss is approximated by running HERTA for more epochs.) and take logarithm scale to make the comparison clearer.

From the results, it is clear that on all datasets and all optimizers we consider, HERTA converges much faster than standard methods. It generally requires $\leq 10$ iterations to get very close to the smallest training loss. This shows that the guarantee obtained in our theoretical results (Theorem 5.1) not only matches the experiments, but also holds when the setting is slightly changed (using other optimizers). Thus, these experimental results verify the universality of HERTA as an optimization framework.

We also conduct extensive experiments with larger $\lambda$. See Appendix F for the results. These results verify that even when $\lambda$ is relatively large, HERTA converges very fast, which suggests that the dominant term in Theorem 5.1 is the $\tilde{O}(m)$ term instead of the $\tilde{O}(n_{\lambda}\lambda^{2}d)$ term.

### 6.2 Convergence Rate Comparison under Cross Entropy Loss

Note that in the theoretical analysis of Section 5, we focus on the MSE loss. However, for graph node classification tasks, MSE is not the most commonly used loss. Instead, most node classification models use the cross entropy (CE) loss. For example, the original TWIRLS paper (Yang et al., ) uses CE loss.

To this end, we also adopt HERTA on training problems with CE loss. See Appendix E for the details of the implementation. The results are displayed in Figure 2. From these results, it is clear that HERTA also significantly speeds up training with CE loss. The results demonstrate that although originally developed on MSE loss, HERTA is not limited to it and have the flexibility to be applied to other type of loss functions.

Figure: Figure 2: The training loss comparison between HERTA and standard optimizers on cross entropy loss with $\lambda=1$. Dataset used from left to right: ogbn-arxiv, citeseer and pubmed.
Refer to caption: extracted/2403.18142v1/figures-log/ogbn-arxiv-1-fast_ce,ce.png

Remark on the Surprising Effectiveness of HERTA on Cross Entropy Loss.
In Section 6.2, we observe that despite not being guaranteed by the theoretical result, HERTA shows a certain degree of universality in that it also works on CE loss. We claim that this phenomenon might originate from the fact the Hessians of TWIRLS under MSE loss and CE loss are very similar. We provide an analysis of the gradient and Hessian of TWIRLS under MSE and CE losses in Section D.1. The results show that the gradient under CE loss can be viewed as the gradient under MSE loss with one term being normalized by softmax, and the Hessian under CE loss can be viewed as a rescaled version of the Hessian under MSE loss. These comparisons serve as an explanation of why HERTA works so well with CE loss.

## 7 Conclusions

In this paper we present HERTA: a High-Efficiency and Rigorous Training Algorithm that solves the problem of training Unfolded GNNs on a graph within a time essentially of the same magnitude as the time it takes to load the input data. As a component of HERTA, we also propose a new spectral sparsifier that works for normalized and regularized graph Laplacian matrices.

Experimental results on real world datasets show the effectiveness of HERTA. Moreover, it is shown that HERTA works for various loss functions and optimizers, despite being derived from a specific loss function and optimizer. This shows the universality of HERTA and verifies that it is ready to use in practice.

### 7.1 Future Directions

In this subsection, we briefly discuss some potential future directions to explore with HERTA.

Extending to more Complex Models. In this paper, we stick to a rather straightforward implementation of TWIRLS.
Although it is verified in (Yang et al., ) that this implementation is already strong enough in many cases, under more general and subtle scenarios we might want to use more complex implementations of TWIRLS or other Unfolded GNN models. Therefore, it is critical to extend HERTA to those more complex settings. We discuss possible directions for implementing HERTA for more complex $f$ as well as incorporating attention mechanisms in Section D.4.

More Detailed Analysis of CE Loss and other Losses. Although in this paper we provided an analysis explaining why HERTA can work on CE loss, finding a theoretically rigorous bound for CE loss as well as other loss functions is another important direction for future work.

## Appendix A Introduction to the Mathematical Tools

In the main paper, due to space limitations, our introduction of certain mathematical tools and concepts remains brief. In this section, we provide a more comprehensive introduction to the mathematical tools utilized in the main paper.

### A.1 Subsampling

Subsampling is used in multiple places in our Algorithm. When we say subsampling $s$ rows of a matrix ${\bm{M}}\in\mathbb{R}^{n\times d}$ with probability $\{p_{k}\}_{k=1}^{n}$ and obtain a new matrix $\tilde{\bm{M}}\in\mathbb{R}^{s\times d}$, it means each row of $\tilde{\bm{M}}$ is an i.i.d. random vector which follows the same distribution of a random vector ${\bm{\xi}}$ that satisfies the following property:

$$ $\displaystyle\mathbb{P}\left\{{\bm{\xi}}=\sqrt{\frac{1}{sp_{k}}}{\bm{m}}_{k} \right\}=p_{k},$ (A.1) $$

where ${\bm{m}}_{k}$ is the $k$-th row of matrix ${\bm{M}}$. It is clear that constructing such a subsampled matrix $\tilde{\bm{M}}$ takes time $O(sd)$. Moreover, if $p_{k}=\frac{1}{n}$, we say the subsampling is with uniform probability.

The operator of subsampling $s$ rows of a $n$-row matrix is clearly a linear transformation. Therefore it can be represented by a random matrix ${\bm{S}}\in\mathbb{R}^{s\times n}$. Notice that when we apply ${\bm{S}}$ to a matrix (i.e. ${\bm{S}}\cdot{\bm{M}}$) we don’t actually need to construct such a matrix ${\bm{S}}$ and calculate matrix product, since we only carry out sampling.

### A.2 SDD Solvers

As mentioned in the main paper, for a sparse SDD matrix ${\bm{M}}$, it is possible to fast approximate its linear solver. The following Lemma A.1 rigorously states the bound we can obtain for SDD solvers.

###### Lemma A.1 ( Peng & Spielman 2014 ) .

For any SDD matrix ${\bm{M}}\in\mathbb{R}^{n\times n}$ and any real number $\epsilon>0$, there exists an algorithm to construct an operator $\mathsf{SDDSolver}_{\epsilon}\left({\bm{M}};\cdot\right)$, such that $\mathsf{SDDSolver}_{\epsilon}\left({\bm{M}};\cdot\right)$ is a linear solver for ${\bm{M}}$ with $\epsilon$ error rate and for any ${\bm{x}}\in\mathbb{R}^{n}$, calculating $\mathsf{SDDSolver}_{\epsilon}({\bm{M}};{\bm{x}})$ takes time

$$ $\displaystyle O\left(m\left(\log\frac{1}{\epsilon}\right)\mathop{\mathrm{poly} }\log(m)\mathop{\mathrm{poly}}\log\left(\kappa({\bm{M}})\right)\right),$ (A.2) $$

where $m=\mathrm{nnz}({\bm{M}})$.

###### Lemma A.1 ( Peng & Spielman 2014 ) .

### A.3 Fast Matrix Multiplication

As mentioned in the main paper, we use SRHT to achieve fast matrix multiplication. Here we further explain this process. We refer interested readers to (Tropp, ) for more details.
For a positive integer number $n$ that is a power of $2$, the Hadamard transformation of size $n$ is a linear transformation recursively defined as follows:

$$ $\displaystyle{\bm{H}}_{n}:=\begin{bmatrix}{\bm{H}}_{n/2}&{\bm{H}}_{n/2}\\ {\bm{H}}_{n/2}&-{\bm{H}}_{n/2}\end{bmatrix}$ (A.3) $$

and ${\bm{H}}_{1}=1$. For a vector ${\bm{x}}\in\mathbb{R}^{n}$, we define

$$ $\displaystyle\mathsf{Hadamard}({\bm{x}}):=\frac{1}{\sqrt{n}}{\bm{H}}_{n}{\bm{x }}.$ (A.4) $$

Notice that when $n$ is not a power of $2$, we pad the vector ${\bm{x}}$ with 0 0-s beforehand, therefore the requirement of $n$ being a power of $2$ is ignored in the actual usage. From its recursive nature, it is not hard to see that applying a Hadamard transformation to each column of an $n\times d$ matrix only takes $O(d\log n)$ time.

Let ${\bm{R}}\in\mathbb{R}^{n\times n}$ be a diagonal matrix whose diagonal entries are i.i.d. Rademacher variables, and ${\bm{S}}\in\mathbb{R}^{s\times n}$ be a subsampling matrix with uniform probability and $s=O\left(\frac{d}{\beta^{2}}\log n\right)$. For a matrix ${\bm{Q}}\in\mathbb{R}^{n\times d}$, its SRHT with $\beta$ error rate is defined as

$$ $\displaystyle\mathsf{SRHT}_{\beta}({\bm{Q}}):={\bm{S}}\mathsf{Hadamard}\left({ \bm{R}}{\bm{Q}}\right),$ (A.5) $$

which is exactly the matrix $\tilde{\bm{Q}}$ we use in Algorithm 1. A fast matrix multiplication result can be achieved using SRHT.

###### Lemma A.2 ( Tropp 2011 ) .

For matrix ${\bm{Q}}\in\mathbb{R}^{n\times d}$ where $n\geq d$ and $\mathrm{rank}(\bm{Q})=d$, let $\tilde{\bm{Q}}=\mathsf{SRHT}_{\beta}({\bm{Q}})$ where $\beta\in(0,1/4)$, then we have

$$ $\displaystyle\tilde{\bm{Q}}^{\top}\tilde{\bm{Q}}\approx_{\beta}{\bm{Q}}^{\top} {\bm{Q}}$ (A.6) $$

with probability at least $1-\frac{1}{2n}$.

###### Lemma A.2 ( Tropp 2011 ) .

## Appendix B Running Time of the Original Implementation Used in (Yang et al., 2021b )

In this section we analyze the time complexity of the implementation used in (Yang et al., ), which uses gradient descent to solve both inner and outer problem.

#### Inner Problem Analysis.

We first analyze the time complexity of the inner problem solver used in (Yang et al., ). Recall that for the inner problem eq. 4.3, we need to find an approximation of the linear solver of ${\bm{I}}+\lambda\hat{\bm{L}}$. In (Yang et al., ), it is implemented by a standard gradient descent. Here we consider approximately solving the following least square problem by gradient descent:

$$ $\displaystyle{\bm{v}}=\arg\min_{{\bm{v}}\in\mathbb{R}^{n}}\frac{1}{2}\left\| \left({\bm{I}}+\lambda\hat{\bm{L}}\right){\bm{v}}-{\bm{u}}\right\|_{2}^{2}.$ (B.1) $$

The gradient step with step size $\mu$ is:

$$ $\displaystyle{\bm{v}}^{(t+1)}=(1-\mu)\left({\bm{I}}+\lambda\hat{\bm{L}}\right) ^{2}{\bm{v}}^{(t)}+\eta\left({\bm{I}}+\lambda\hat{\bm{L}}\right){\bm{u}}.$ (B.2) $$

###### Theorem B.1 (Inner Analysis) .

If we update ${\bm{v}}^{(t)}$ through eq. B.2 with initialization ${\bm{v}}^{(0)}={\bm{0}}$ and a proper $\eta$, then after $T=O\left(\lambda^{2}\log\frac{{\lambda}}{\epsilon^{2}}\right)$ iterations we can get

$$ $\displaystyle\left\|{\bm{v}}^{(T)}-\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{ -1}{\bm{u}}\right\|_{({\bm{I}}+\lambda\hat{\bm{L}})}\leq\epsilon\left\|\left({ \bm{I}}+\lambda\hat{\bm{L}}\right)^{-1}{\bm{u}}\right\|_{({\bm{I}}+\lambda\hat {\bm{L}})}$ (B.3) $$

for any $\epsilon\in(0,1)$.

###### Proof.

Let ${\bm{H}}={\bm{I}}+\lambda\hat{\bm{L}}$. Notice that the strongly convexity and Lipschitz constant of Problem eq. B.1 is $\mathscr{C}=1$ and $L\leq 9\lambda^{2}$ respectively.
From Lemma C.1, take $\eta=\frac{1}{L}$ and ${\bm{v}}^{(0)}={\bm{0}}$, we have

$$ $\displaystyle\|{\bm{H}}{\bm{v}}^{(T)}-{\bm{u}}\|^{2}\leq(1-9\lambda^{2})^{T}\| {\bm{u}}\|^{2}$ (B.4) $$

for any $T\in\mathbb{N}$.
Therefore,

$$ $\displaystyle\frac{\|{\bm{v}}^{(T)}-{\bm{H}}^{-1}{\bm{u}}\|^{2}_{{\bm{H}}}}{\| {\bm{H}}^{-1}{\bm{u}}\|^{2}_{{\bm{H}}}}$ $\displaystyle=\frac{\|{\bm{H}}{\bm{v}}^{(T)}-{\bm{u}}\|^{2}_{{\bm{H}}^{-1}}}{ \|{\bm{u}}\|^{2}_{{\bm{H}}^{-1}}}$ (B.5) $\displaystyle\leq 3\lambda\frac{\|{\bm{H}}{\bm{v}}^{(T)}-{\bm{u}}\|_{2}^{2}}{ \|{\bm{u}}\|_{2}^{2}}$ (B.6) $\displaystyle\leq 3\lambda(1-9\lambda^{2})^{{T}}.$ (B.7) $$

To obtain an error rate $\epsilon$, the smallest number of iterations $T$ needed is

$$ $\displaystyle T=\left(\log\frac{1}{1-(9\lambda^{2})^{-1}}\right)^{-1}\log\frac {{3\lambda}}{\epsilon^{2}}+1=O\left[\left(\log\frac{1}{1-\lambda^{-2}}\right)^ {-1}\log\frac{\kappa}{\epsilon^{2}}\right]=O\left(\lambda^{2}\log\frac{\lambda }{\epsilon^{2}}\right).$ (B.8) $$

∎

###### Theorem B.1 (Inner Analysis) .

###### Proof.

#### Overall Running Time.

As we proved above in Theorem B.1, the number of iterations needed for solving the inner problem to error rate $\beta$ with the implementation of (Yang et al., ) is $\tilde{O}\left(\lambda^{2}\log\frac{1}{\beta}\right)$, which is related to the hyper-parameter $\lambda$. For each inner iteration (i.e. eq. B.2), we need to compute $({\bm{I}}+\lambda{\bm{L}})({\bm{X}}{\bm{w}}-\hat{\bm{y}})$. Since this is a sparse matrix multiplication, the complexity of this step is $O(m+nd)$. Putting things together, we have the time complexity of calling inner problem solver is $\tilde{O}\left((m+nd)\lambda^{2}\log\frac{1}{\beta}\right)$.

From Lemma 5.1, the number of outer iterations needed is $\tilde{O}\left(\kappa_{o}\log\frac{1}{\epsilon}\right)$, where $\kappa_{o}$ is the condition number of the outer problem. Moreover, Lemma 5.1 also indicates that we require $\beta\leq\frac{\epsilon^{1/2}}{25\kappa({\bm{X}})\lambda}$, so solving inner problem takes $\tilde{O}(\lambda^{2}\log 1/\beta)=\tilde{O}(\lambda^{2}\log 1/\epsilon)$. In each outer iteration, we need to call inner problem solver constant times, and as well
as performing constant matrix vector multiplications whose complexity is $O(nd)$. Therefore the overall running time of solving the training problem is $\tilde{O}\left[\kappa_{o}\left(\lambda^{2}(m+nd)+nd\right)\left(\log\frac{1}{
\epsilon}\right)^{2}\right]$.

## Appendix C Proof of Theoretical Results

We first note that since $\hat{\bm{L}}$ is normalized, all the eigenvalues of $\hat{\bm{L}}$ are in the range $[0,2)$. Therefore, we have $\sigma_{\min}({\bm{I}}+\lambda\hat{\bm{L}})=1$ and $\sigma_{\max}({\bm{I}}+\lambda\hat{\bm{L}})\leq 1+2\lambda$. In the whole paper we view $\lambda$ as a large value (i.e. $\lambda\gg 1$). In practice, it is possible to use a small $\lambda$, in which case the algorithm works no worse than the case where $\lambda=1$. Therefore the $\lambda$ used in the paper should actually be understood as $\max\{\lambda,1\}$. With this assumption, below we assume $\sigma_{\max}({\bm{I}}+\lambda\hat{\bm{L}})\leq 3\lambda$ for convenience.

### C.1 Descent Lemma

In this subsection, we introduce a descent lemma that we will use to analyze the convergence rate. This is a standard result in convex optimization, and we refer interested readers to (Bottou et al., ) for more details.

###### Lemma C.1 .

If $\ell:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-Lipschitz smooth and $c$-strongly convex, and we have a sequence of points $\left\{{\bm{w}}^{(t)}\right\}_{t=1}^{T}$ in $\mathbb{R}^{d}$ such that

$$ $\displaystyle{\bm{w}}^{(t+1)}={\bm{w}}^{(t)}-\eta{\bm{g}}^{(t)},$ (C.1) $$

where $\left\|{\bm{g}}^{(t)}-\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_{2}\leq
\gamma\left\|\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_{2}$ and $\gamma<1$, $\eta=\frac{1-\gamma}{(1+\gamma)^{2}L}$, then we have

$$ $\displaystyle\ell\left({\bm{w}}^{(T)}\right)-\ell^{*}\leq\left[1-\kappa^{-1} \left(\frac{1-\gamma}{1+\gamma}\right)^{2}\right]^{T}\left[\ell\left({\bm{w}}^ {(0)}\right)-\ell^{*}\right],$ (C.2) $$

where $\ell^{*}=\inf_{{\bm{w}}\in\mathbb{R}^{d}}\ell({\bm{w}})$ and $\kappa=\frac{L}{c}$.

###### Proof.

From the condition $\left\|{\bm{g}}^{(t)}-\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_{2}\leq
\gamma\left\|\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_{2}$, we have that

$$ $\displaystyle(1-\gamma)\|{\bm{g}}^{(t)}\|_{2}\leq\left\|\nabla\ell\left({\bm{w }}^{(t)}\right)\right\|_{2}\leq(1+\gamma)\|{\bm{g}}^{(t)}\|_{2}$ (C.3) $$

and

$$ $\displaystyle-2\left<g^{(t)},\nabla\ell\left({\bm{w}}^{(t)}\right)\right>$ $\displaystyle=\left\|{\bm{g}}^{(t)}-\nabla\ell\left({\bm{w}}^{(t)}\right) \right\|_{2}^{2}-\left\|{\bm{g}}^{(t)}\right\|_{2}^{2}-\left\|\nabla\ell\left( {\bm{w}}^{(t)}\right)\right\|_{2}^{2}$ (C.4) $\displaystyle\leq\left[(\gamma^{2}-1)-(1-\gamma)^{2}\right]\left\|\nabla\ell \left({\bm{w}}^{(t)}\right)\right\|_{2}^{2}$ (C.5) $\displaystyle=2(\gamma-1)\left\|\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_ {2}^{2}$ (C.6) $$

From Lipschitz smoothness, we have

$$ $\displaystyle\ell\left({\bm{w}}^{(t+1)}\right)-\ell\left({\bm{w}}^{(t)}\right)$ $\displaystyle\leq-\eta\left<{\bm{g}}^{(t)},\nabla\ell\left({\bm{w}}^{(t)} \right)\right>+\frac{1}{2}L\eta^{2}\|{\bm{g}}^{(t)}\|_{2}^{2}$ (C.7) $\displaystyle\leq\left(\eta(\gamma-1)+\frac{1}{2}L\eta^{2}(1+\gamma)^{2}\right )\left\|\nabla\ell\left({\bm{w}}^{(t)}\right)\right\|_{2}^{2}.$ (C.8) $$

It’s not hard to show that the optimal $\eta$ for eq. C.8 is $\eta=\frac{1-\gamma}{L(1+\gamma)^{2}}$. Substituting $\eta=\frac{1-\gamma}{L(1+\gamma)^{2}}$ to eq. C.8 and use convexity we can get

$$ $\displaystyle\ell\left({\bm{w}}^{(t+1)}\right)-\ell^{*}$ $\displaystyle\leq\left[\ell\left({\bm{w}}^{(t)}\right)-\ell^{*}\right]-\frac{( 1-\gamma)^{2}}{2L(1+\gamma)^{2}}\left\|\nabla\ell\left({\bm{w}}^{(t)}\right) \right\|_{2}^{2}$ (C.9) $\displaystyle\leq\left[\ell\left({\bm{w}}^{(t)}\right)-\ell^{*}\right]-\frac{c (1-\gamma)^{2}}{L(1+\gamma)^{2}}\left[\ell\left({\bm{w}}^{(t)}\right)-\ell^{*}\right]$ (C.10) $$

By induction we have

$$ $\displaystyle\ell\left({\bm{w}}^{(T)}\right)-\ell^{*}$ $\displaystyle\leq\left[1-\kappa^{-1}\left(\frac{1-\gamma}{1+\gamma}\right)^{2} \right]^{T}\left[\ell\left({\bm{w}}^{(0)}\right)-\ell^{*}\right].$ (C.11) $$

∎

###### Lemma C.1 .

###### Proof.

### C.2 Bound of Loss Value by Gradient

In this subsection, we derive an inequality that allows us to bound the value of loss function by the norm of gradient.

###### Lemma C.2 .

If $\ell:\mathbb{R}^{d}\to\mathbb{R}$ is a $\mathscr{C}$-strongly convex and smooth function, and ${\bm{w}}^{*}=\arg\min_{{\bm{w}}\in\mathbb{R}^{d}}\ell({\bm{w}})$ is a global optimal, then for any ${\bm{w}}\in\mathbb{R}^{d}$ and $\epsilon\in(0,1)$, we have

$$ $\displaystyle\ell({\bm{w}})\leq\max\left\{(1+\epsilon)\ell({\bm{w}}^{*}),\frac {1}{\epsilon\mathscr{C}}\|\nabla\ell({\bm{w}})\|^{2}\right\}.$ (C.12) $$

###### Proof.

Using the inequality (4.12) from (Bottou et al., ), we have

$$ $\displaystyle\ell({\bm{w}})\leq\frac{1}{2\mathscr{C}}\|\nabla\ell({\bm{w}})\|^ {2}+\ell({\bm{w}}^{*}).$ (C.13) $$

If $\ell({\bm{w}})\geq(1+\epsilon)\ell({\bm{w}}^{*})$, we have

$$ $\displaystyle\ell({\bm{w}})\leq\frac{1}{2\mathscr{C}}\|\nabla\ell({\bm{w}})\|^ {2}+\frac{1}{1+\epsilon}\ell({\bm{w}}).$ (C.14) $$

Shifting the terms in eq. C.14 gives $\ell({\bm{w}})\leq\frac{1+\epsilon}{2\mathscr{C}\epsilon}\|\nabla\ell({\bm{w}}
)\|^{2}\leq\frac{1}{\epsilon\mathscr{C}}\|\nabla\ell({\bm{w}})\|^{2}$, which proves the claim.

∎

###### Lemma C.2 .

###### Proof.

### C.3 Proof of Lemma 5.1

Given a point ${\bm{w}}\in\mathbb{R}^{d}$, we first consider the error between the estimated gradient $\widetilde{\nabla\ell_{i}({\bm{w}})}$ and the true gradient $\nabla\ell_{i}({\bm{w}})$. Below we denote ${\bm{H}}={\bm{I}}+\lambda\hat{\bm{L}}$ and ${\bm{z}}={\bm{X}}{\bm{w}}-\hat{\bm{y}}_{i}$. It’s not hard to notice that $\ell_{i}({\bm{w}})=\frac{1}{2}\left\|{\bm{H}}^{-1}{\bm{z}}\right\|^{2}$ and $\nabla\ell_{i}({\bm{w}})={\bm{X}}^{\top}{\bm{H}}^{-2}{\bm{z}}$. Moreover, we denote $\ell^{*}_{i}=\inf_{{\bm{w}}}\ell^{*}_{i}$ as the optimal value of $\ell_{i}$.

we have

$$ $\displaystyle\left\|\widetilde{\nabla\ell_{i}({\bm{w}})}-\nabla\ell_{i}({\bm{w }})\right\|=$ $\displaystyle\left\|{\bm{X}}^{\top}\left[\mathcal{S}\left({\bm{H}}^{-1}{\bm{z} }\right)-{\bm{H}}^{-2}{\bm{z}}+\mathcal{S}(\mathcal{S}({\bm{z}}))-\mathcal{S} \left({\bm{H}}^{-1}{\bm{z}}\right)\right]\right\|$ (C.15) $\displaystyle\leq$ $\displaystyle\left\|{\bm{X}}^{\top}\left[\mathcal{S}\left({\bm{H}}^{-1}{\bm{z} }\right)-{\bm{H}}^{-2}{\bm{z}}\right]\right\|+\left\|{\bm{X}}^{\top}\left[ \mathcal{S}(\mathcal{S}({\bm{z}})-{\bm{H}}^{-1}{\bm{z}})\right]\right\|.$ (C.16) $$

As shown, the gradient error can be decomposed into two terms, namely $\left\|\widetilde{\nabla\ell_{i}({\bm{w}})}-\nabla\ell_{i}({\bm{w}})\right\|=E
_{1}+E_{2}$, where $E_{1}=\left\|{\bm{X}}^{\top}\left[\mathcal{S}\left({\bm{H}}^{-1}{\bm{z}}\right
)-{\bm{H}}^{-2}{\bm{z}}\right]\right\|$ and $E_{2}=\left\|{\bm{X}}^{\top}\left[\mathcal{S}(\mathcal{S}({\bm{z}})-{\bm{H}}^{
-1}{\bm{z}})\right]\right\|$. Below we analyze each part separately.

Notice that, as we assumed ${\bm{X}}$ is full-rank, $\ell_{i}$ is a strongly convex function with parameter $\mathscr{C}=\sigma_{\min}\left({\bm{X}}^{\top}{\bm{H}}^{-2}{\bm{X}}\right)\geq
\frac{\sigma_{\min}({\bm{X}})^{2}}{9\lambda^{2}}$. Let $\mathcal{D}=\left\{\left.{\bm{w}}\in\mathbb{R}^{d}\right|\ell_{i}({\bm{w}})
\leq(1+\epsilon)\ell_{i}^{*}\right\}$. If ${\bm{w}}\in\mathcal{D}$, then ${\bm{w}}$ is already a good enough solution. Below we assume ${\bm{w}}\not\in\mathcal{D}$, in which case by Lemma C.2, we have $\ell_{i}({\bm{w}})\leq\frac{1}{\epsilon\mathscr{C}}\|\nabla\ell_{i}({\bm{w}})
\|^{2}$, which means

$$ $\displaystyle\left\|{\bm{H}}^{-1}{\bm{z}}\right\|^{2}\leq\frac{2}{\mathscr{C} \epsilon}\left\|{\bm{X}}^{-1}{\bm{H}}^{-2}{\bm{z}}\right\|^{2}\leq\frac{18 \lambda^{2}}{\sigma_{\min}({\bm{X}})^{2}\epsilon}.$ (C.17) $$

For $E_{1}$, we have

$$ $\displaystyle E_{1}$ $\displaystyle=\left\|{\bm{X}}^{\top}\left[\mathcal{S}\left({\bm{H}}^{-1}{\bm{z }}\right)-{\bm{H}}^{-2}{\bm{z}}\right]\right\|_{2}$ (C.18) $\displaystyle\leq\sigma_{\max}({\bm{X}})\left\|\mathcal{S}\left({\bm{H}}^{-1}{ \bm{z}}\right)-{\bm{H}}^{-2}{\bm{z}}\right\|_{2}$ (C.19) $\displaystyle\leq\sigma_{\max}({{\bm{X}}})\left\|\mathcal{S}\left({\bm{H}}^{-1 }{\bm{z}}\right)-{\bm{H}}^{-2}{\bm{z}}\right\|_{{\bm{H}}}$ (C.20) $\displaystyle\leq\sigma_{\max}({{\bm{X}}})\mu\left\|{\bm{H}}^{-2}{\bm{z}} \right\|_{{\bm{H}}}$ (C.21) $\displaystyle\leq\sigma_{\max}({{\bm{X}}})\mu\left\|{\bm{H}}^{-1}{\bm{z}} \right\|_{2}$ (C.22) $\displaystyle\leq\sqrt{\frac{18}{\epsilon}}\kappa({{\bm{X}}})\lambda\mu\left\| {\bm{X}}^{\top}{\bm{H}}^{-2}{\bm{z}}\right\|_{2}$ (C.23) $\displaystyle\leq 5\epsilon^{-1/2}\kappa({{\bm{X}}})\lambda\mu\|\nabla\ell_{i} ({\bm{w}})\|_{2}.$ (C.24) $$

For $E_{2}$, we have

$$ $\displaystyle E_{2}$ $\displaystyle=\left\|{\bm{X}}^{\top}\left[\mathcal{S}(\mathcal{S}({\bm{z}})-{ \bm{H}}^{-1}{\bm{z}})\right]\right\|_{2}$ (C.25) $\displaystyle\leq{\sigma_{\max}({\bm{X}})}(1+\mu)\|{\bm{H}}^{-1}\left(\mathcal {S}({\bm{z}})-{\bm{H}}^{-1}{\bm{z}}\right)\|_{{\bm{H}}}$ (C.26) $\displaystyle\leq{\sigma_{\max}({\bm{X}})}(1+\mu)\|\mathcal{S}({\bm{z}})-{\bm{ H}}^{-1}{\bm{z}}\|_{{\bm{H}}}$ (C.27) $\displaystyle\leq{\sigma_{\max}({\bm{X}})}(1+\mu)\mu\|{\bm{H}}^{-1}{\bm{z}}\|_ {{\bm{H}}}$ (C.28) $\displaystyle\leq\sqrt{\frac{18}{\epsilon}}\kappa({\bm{X}})(1+\mu)\mu\lambda \sqrt{3\lambda}\|{\bm{X}}^{\top}{\bm{H}}^{-2}{\bm{z}}\|_{2}$ (C.29) $\displaystyle\leq 20\epsilon^{-1/2}\kappa({\bm{X}})\mu\lambda^{2}\|\nabla\ell_ {i}({\bm{w}}_{i})\|_{2}.$ (C.30) $$

Combine the bound of $E_{1}$ and $E_{2}$, we have

$$ $\displaystyle\left\|\widetilde{\nabla\ell_{i}({\bm{w}}_{i})}-\nabla\ell_{i}({ \bm{w}}_{i})\right\|\leq E_{1}+E_{2}\leq 25\epsilon^{-1/2}\kappa({\bm{X}})\mu \lambda^{2}\|\nabla\ell_{i}({\bm{w}}_{i})\|_{2}.$ (C.31) $$

Now, consider the optimization. Let $\left\{{\bm{w}}_{i}^{(t)}\right\}_{t=1}^{T}$ be a sequence such that

$$ $\displaystyle{\bm{w}}_{i}^{(t+1)}={\bm{w}}^{(t)}-\eta\widetilde{\nabla\ell_{i} \left({\bm{w}}_{i}^{(t)}\right)}.$ (C.32) $$

Let $\kappa$ be the condition number of ${\bm{X}}^{\top}{\bm{H}}^{-2}{\bm{X}}$ and let $\ell_{i}^{*}$ be the optimal value of $\ell_{i}$. Let $\gamma=25\epsilon^{-1/2}\kappa({\bm{X}})\lambda^{2}\mu\leq\frac{1}{2}$. From Lemma C.1, we have with a proper value of $\eta$,

$$ $\displaystyle\ell_{i}\left({\bm{w}}_{i}^{(T)}\right)-\ell^{*}$ $\displaystyle\leq\left(1-\left(\frac{1-\gamma}{1+\gamma}\right)^{2}\kappa^{-1} \right)^{T}\left[\ell_{i}\left({\bm{w}}_{i}^{(0)}\right)-\ell_{i}^{*}\right]$ (C.33) $\displaystyle\leq\left(1-(9\kappa)^{-1}\right)^{T}\left[\ell_{i}\left({\bm{w}} _{i}^{(0)}\right)-\ell_{i}^{*}\right].$ (C.34) $$

Therefore, the number of iterations $T$ required to achieve $\epsilon$ error rate is

$$ $\displaystyle T=O\left(\left(\log\frac{1}{1-(9\kappa)^{-1}}\right)^{-1}\log \frac{1}{\epsilon}\right)=O\left(\kappa\log\frac{1}{\epsilon}\right).$ (C.35) $$

### C.4 Proof of Lemma 5.2

The basic idea of Algorithm 2 is to use ridge leverage score sampling methods to obtain a spectral sparsifier.
Let $\hat{\bm{B}}={\bm{B}}{\bm{D}}^{-1/2}$ be the normalized incidence matrix and $\hat{\bm{b}}_{i}$ the $i$-th row of $\hat{\bm{B}}$. Given $\lambda^{-1}>0$, for $i\in\{1,2,\cdots,m\}$, the $i$-th ridge leverage score is defined as

$$ $\displaystyle l_{i}:=$ $\displaystyle~{}\hat{\bm{b}}_{i}^{\top}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1} \hat{\bm{b}}_{i}$ (C.36) $\displaystyle=$ $\displaystyle~{}\hat{\bm{b}}_{i}^{\top}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1} \hat{\bm{L}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_{i}+\lambda^{-1 }\hat{\bm{b}}_{i}^{\top}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-2}\hat{\bm{b}}_{i}$ (C.37) $\displaystyle=$ $\displaystyle~{}\|\hat{\bm{B}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b }}_{i}\|^{2}+\lambda^{-1}\|(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_ {i}\|^{2}.$ (C.38) $$

It is not affordable to compute all the $m$ ridge leverage scores exactly. Therefore, we first use Johnson–Lindenstrauss lemma to reduce the dimension. Recall that in Algorithm 2 we define $\bm{\Pi}_{1}\in\mathbb{R}^{k\times m}$ and $\bm{\Pi}_{2}\in\mathbb{R}^{k\times n}$ to be Gaussian sketches (that is, each entry of the matrices are i.i.d Gaussian random variables $\mathcal{N}(0,1/k)$). We set $k=O(\log m)$, and by using Lemma C.3 we have the following claim holds with probability at least $1-\frac{1}{8n}$:

$$ $\displaystyle\|\bm{\Pi}_{1}\hat{\bm{B}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1} \hat{\bm{b}}_{i}\|\approx_{2^{1/4}-1}\|\hat{\bm{B}}(\hat{\bm{L}}+\lambda^{-1} \bm{I})^{-1}\hat{\bm{b}}_{i}\|~{}~{}~{}\text{for all $i\in\{1,2,\cdots,m\}$}.$ (C.39) $$

Similarly,the following claim also holds with probability at least $1-\frac{1}{8n}$:

$$ $\displaystyle\|\bm{\Pi}_{2}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_ {i}\|\approx_{2^{1/4}-1}\|(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_{ i}\|~{}~{}~{}\text{for all $i\in\{1,2,\cdots,m\}$}.$ (C.40) $$

By summing over above two inequalities (after squared) and taking an union bound, with probability $1-\frac{1}{4n}$ we have

$$ $\displaystyle\|\bm{\Pi}_{1}\hat{\bm{B}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1} \hat{\bm{b}}_{i}\|^{2}+\lambda^{-1}\|\bm{\Pi}_{2}(\hat{\bm{L}}+\lambda^{-1}\bm {I})^{-1}\hat{\bm{b}}_{i}\|^{2}\approx_{\sqrt{2}-1}l_{i}~{}~{}~{}\text{for all }~{}~{}~{}i\in\{1,2,\cdots,m\}.$ (C.41) $$

However it is still too expensive to compute $\bm{\Pi}_{1}\hat{\bm{B}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}$ and $\bm{\Pi}_{2}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}$, since computing $(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}$ itself takes prohibitive $O(n^{3})$ time. Instead, notice that $\hat{\bm{L}}+\lambda^{-1}\bm{I}$ is a SDD matrix, thus we can apply the SDD solver to the $k$ columns of matrix $(\bm{\Pi}_{1}\hat{\bm{B}})^{\top}$ and $(\bm{\Pi}_{2})^{\top}$ respectively.
According to Lemma A.1, there is a linear operator $\mathsf{SDDSolver}_{\epsilon}({\hat{\bm{L}}+\lambda^{-1}\bm{I}};\bm{x})$ that runs in time $\tilde{O}(\mathrm{nnz}({\hat{\bm{L}}+\lambda^{-1}\bm{I}})\cdot\log 1/\epsilon)
=\tilde{O}(m\log 1/\epsilon)$ such that for any $\bm{x}^{\top}\in\mathbb{R}^{n}$, it outputs $\tilde{\bm{x}}$ that satisfies $\|\tilde{\bm{x}}-\bm{x}^{\top}({\hat{\bm{L}}+\lambda^{-1}\bm{I}})^{-1}\|_{{
\hat{\bm{L}}+\lambda^{-1}\bm{I}}}\leq\epsilon\|\bm{x}^{\top}({\hat{\bm{L}}+
\lambda^{-1}})\|_{{\hat{\bm{L}}+\lambda^{-1}\bm{I}}}$. For our purpose we set $\epsilon=2^{1/4}$. Denote $\bm{B}_{\mathcal{S}}$ as the matrix obtained by applying each column of $(\bm{\Pi}_{1}\hat{\bm{B}})^{\top}$ to $\mathsf{SDDSolver}_{\epsilon}({\hat{\bm{L}}+\lambda^{-1}\bm{I}};\bm{x})$, that is, the $j$-th row of $\mathsf{SDDSolver}_{2^{1/4}}({\hat{\bm{L}}+\lambda^{-1}\bm{I}};(\bm{\Pi}_{1}
\hat{\bm{B}})^{\top}_{j})$. Similarly we denote $\bm{\Pi}_{\mathcal{S}}$ as the matrix with $j$-th row equals to $\mathsf{SDDSolver}_{2^{1/4}}({\hat{\bm{L}}+\lambda^{-1}\bm{I}};(\bm{\Pi}_{2})^
{\top}_{j})$. By using Lemma A.1 to both solvers we have

$$ $\displaystyle\|\bm{B}_{\mathcal{S}}\hat{\bm{b}}_{i}\|^{2}+\lambda^{-1}\|\bm{ \Pi}_{\mathcal{S}}\hat{\bm{b}}_{i}\|^{2}\approx_{\sqrt{2}-1}\|\bm{\Pi}_{1}\hat {\bm{B}}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_{i}\|^{2}+\lambda^{ -1}\|\bm{\Pi}_{2}(\hat{\bm{L}}+\lambda^{-1}\bm{I})^{-1}\hat{\bm{b}}_{i}\|^{2}~ {}~{}~{}\text{for all}~{}~{}~{}i\in\{1,2,\cdots,m\}.$ (C.42) $$

By Equation C.41 and Equation C.42, if we set $\tilde{l}_{i}:=\|\bm{B}_{\mathcal{S}}\hat{\bm{b}}_{i}\|^{2}+\lambda^{-1}\|\bm{
\Pi}_{\mathcal{S}}\hat{\bm{b}}_{i}\|^{2}$, then with probability $1-\frac{1}{4n}$, we obtain all the approximation of ridge leverage scores $\{\tilde{l}_{i}\}_{i=1}^{m}$ such that $\tilde{l}_{i}\approx_{1/2}l_{i}$ holds for all $i$.
Notice that since $\bm{B}_{\mathcal{S}},\bm{\Pi}_{\mathcal{S}}\in\mathbb{R}^{k\times n}$, and that each $\hat{\bm{b}}_{i}$ only contains $2$ non-zero entries, thus it takes $\tilde{O}(k\cdot m)$ to pre-compute $\bm{B}_{\mathcal{S}}$ and $\bm{\Pi}_{\mathcal{S}}$, and takes $O(2k\cdot m)$ to compute $\{\bm{B}_{\mathcal{S}}\hat{\bm{b}}_{i},\bm{\Pi}_{\mathcal{S}}\hat{\bm{b}}_{i}
\}_{i=1}^{m}$. To summarize, computing all $\tilde{l}_{i}$ takes $\tilde{O}(km)=\tilde{O}(m\log m)=\tilde{O}(m)$.
With these ridge leverage score approximations, we apply Lemma C.4 to matrix $\hat{\bm{B}}$ with choice $\delta=1/4n$. By setting $\tilde{\bm{L}}:=\hat{\bm{B}}^{\top}\bm{S}^{\top}\bm{S}\hat{\bm{B}}$ we have $\tilde{\bm{L}}+\lambda^{-1}\bm{I}\approx_{\epsilon}\hat{\bm{L}}+\lambda^{-1}
\bm{I}$ holds with probability $1-1/4n$. By applying another union bound we obtain our final result:

$$ $\displaystyle\tilde{\bm{L}}+\lambda^{-1}\bm{I}\approx_{\epsilon}\hat{\bm{L}}+ \lambda^{-1}\bm{I}~{}~{}~{}\text{with probability}~{}~{}~{}1-\frac{1}{2n}.$ (C.43) $$

Finally, according to Lemma C.4, the number of edges of $\tilde{\bm{L}}$ is $s=Cn_{\lambda}\log(n^{2})/\epsilon^{2}=O(n_{\lambda}\log n/\epsilon^{2})$. Since the last step of computing $\tilde{\bm{L}}=(\bm{S}\hat{\bm{B}})^{\top}(\bm{S}\hat{\bm{B}})$ only takes $O(s)=\tilde{O}(n_{\lambda}/\epsilon^{2})$ due to the sparsity of $\hat{\bm{B}}$, the overall time complexity of Algorithm 2 is $\tilde{O}(m+n_{\lambda}/\epsilon^{2})$.

###### Lemma C.3 (Johnson–Lindenstrauss, (Dasgupta & Gupta, 2003 ) ) .

Let $\bm{\Pi}\in\mathbb{R}^{k\times n}$ be Gaussian sketch matrix with each entry independent and equal to $\mathcal{N}(0,1/k)$ where $\mathcal{N}(0,1)$ denotes a standard Gaussian random variable. If we choose $k=O\left(\frac{\log(1/\delta)}{\epsilon^{2}}\right)$, then for any vector $\bm{x}\in\mathbb{R}^{n}$, with probability $1-\delta$ we have

$$ $\displaystyle(1-\epsilon)\|\bm{x}\|\leq\|\bm{\Pi}\bm{x}\|\leq(1+\epsilon)\|\bm {x}\|.$ (C.44) $$

###### Lemma C.4 (Spectral approximation, (Cohen et al., 2017b ) ) .

Let $\bm{S}$ be an $s\times m$ subsampling matrix with probabilities $p_{i}=\tilde{\ell}_{i}/Z$ where $\tilde{\ell}_{i}\approx_{1/2}\ell_{i}$ and $Z$ is the normalization constant. If we have $s\geq Cn_{\lambda}\log(n/\delta)/\epsilon^{2}$ for some constant $C>0$ and $\epsilon,\delta\in(0,1/2]$, then we have

$$ $\displaystyle\hat{\bm{B}}^{\top}\bm{S}^{\top}\bm{S}\hat{\bm{B}}+\lambda^{-1} \bm{I}\approx_{\epsilon}\hat{\bm{B}}^{\top}\hat{\bm{B}}+\lambda^{-1}\bm{I}$ (C.45) $$

holds with probability $1-\delta$. Here $n_{\lambda}=\mathrm{Tr}[\hat{\bm{L}}(\hat{\bm{L}}+\lambda^{-1}{\bm{I}})^{-1}]$.

###### Lemma C.3 (Johnson–Lindenstrauss, (Dasgupta & Gupta, 2003 ) ) .

###### Lemma C.4 (Spectral approximation, (Cohen et al., 2017b ) ) .

### C.5 Proof of Lemma 5.3

We first prove a lemma which is related to the approximation rate of squared matrices.

###### Lemma C.5 .

Suppose that ${\bm{\Sigma}}$ and $\widetilde{\bm{\Sigma}}$ are two $n\times n$ PD matrices, and $\widetilde{\bm{\Sigma}}\approx_{\frac{\beta}{\kappa}}{\bm{\Sigma}}$, where $\kappa$ is the condition number of ${\bm{\Sigma}}$ and $\beta\in\left(0,\frac{1}{8}\right)$, then we have

$$ $\displaystyle\widetilde{\bm{\Sigma}}^{2}\approx_{8\beta}{\bm{\Sigma}}^{2}.$ (C.46) $$

###### Proof.

Let $\epsilon=\frac{\beta}{\kappa}$.
The condition ${\bm{\Sigma}}\approx_{\epsilon}\widetilde{\bm{\Sigma}}$ implies that $\left\|{\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right\|\leq\epsilon\left\|{\bm{
\Sigma}}\right\|$. We have

$$ $\displaystyle\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right){\bm{x}}\right\|$ $\displaystyle=\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right){\bm{ \Sigma}}^{-1}{\bm{\Sigma}}{\bm{x}}\right\|$ (C.47) $\displaystyle\leq\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right){\bm {\Sigma}}^{-1}\right\|\times\left\|{\bm{\Sigma}}{\bm{x}}\right\|$ (C.48) $\displaystyle\leq\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right) \right\|\left\|{\bm{\Sigma}}^{-1}\right\|\times\left\|{\bm{\Sigma}}{\bm{x}}\right\|$ (C.49) $\displaystyle\leq\epsilon\left\|{\bm{\Sigma}}\right\|\left\|{\bm{\Sigma}}^{-1} \right\|\times\left\|{\bm{\Sigma}}{\bm{x}}\right\|$ (C.50) $\displaystyle=\epsilon\kappa\left\|{\bm{\Sigma}}{\bm{x}}\right\|$ (C.51) $\displaystyle=\beta\left\|{\bm{\Sigma}}{\bm{x}}\right\|.$ (C.52) $$

For any ${\bm{x}}\in\mathbb{R}^{n}$, by triangle inequality we have

$$ $\displaystyle\left\|{\bm{\Sigma}}{\bm{x}}\right\|-\left\|\left({\bm{\Sigma}}- \widetilde{\bm{\Sigma}}\right){\bm{x}}\right\|$ $\displaystyle\leq\left\|\widetilde{\bm{\Sigma}}{\bm{x}}\right\|\leq\left\|{\bm {\Sigma}}{\bm{x}}\right\|+\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}} \right){\bm{x}}\right\|.$ (C.53) $$

Subtracting the inequality $\left\|\left({\bm{\Sigma}}-\widetilde{\bm{\Sigma}}\right){\bm{x}}\right\|\leq
\beta\left\|{\bm{\Sigma}}{\bm{x}}\right\|$ we derived before into eq. C.53 and squaring all sides, we have

$$ $\displaystyle(1-\beta)^{2}{\bm{x}}^{\top}{\bm{\Sigma}}^{2}{\bm{x}}\leq{\bm{x}} ^{\top}\widetilde{\bm{\Sigma}}{\bm{x}}\leq(1+\beta)^{2}{\bm{x}}^{\top}{\bm{ \Sigma}}^{2}{\bm{x}}.$ (C.54) $$

Since $0<\beta<\frac{1}{8}$, the claim is proved.
∎

Let ${\bm{T}}={\bm{X}}^{\top}\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{-2}{\bm{X}}$ be the true Hessian. Let $\tilde{\bm{T}}={\bm{X}}^{\top}\left({\bm{I}}+\lambda\tilde{\bm{L}}\right)^{-2}
{\bm{X}}$ be the approximated Hessian using sparsified Laplacian $\tilde{\bm{L}}$. By Lemma 5.2, ${\bm{I}}+\lambda\tilde{\bm{L}}\approx_{\frac{\beta}{3\lambda}}{\bm{I}}+\lambda
{\bm{L}}$ with probability at least $1-\frac{1}{2n}$. From Lemma C.5, we have $\tilde{\bm{T}}\approx_{8\beta}{\bm{T}}$.

Next, we show that ${\bm{Q}}^{\top}{\bm{Q}}\approx_{O(\beta)}{\bm{T}}$. Let $\tilde{\bm{H}}={\bm{I}}+\lambda\tilde{\bm{L}}$. Let $\mathcal{S}$ be the operator defined by $\mathsf{SDDSolver}_{\frac{\beta}{\sqrt{3\lambda}}}\left(\tilde{\bm{H}},\cdot\right)$, i.e. ${\bm{q}}_{j}=\mathcal{S}({\bm{x}}_{j})$ where ${\bm{q}}_{j}$ and ${\bm{x}}_{j}$ are the $j$-th column of ${\bm{Q}}$ and ${\bm{X}}$ respectively. From Lemma A.1 we have for any ${\bm{z}}\in\mathbb{R}^{d}$,

$$ $\displaystyle\left\|{\bm{Q}}{\bm{z}}-\tilde{\bm{H}}^{-1}{\bm{X}}{\bm{z}}\right \|_{2}^{2}$ $\displaystyle=\left\|\sum_{j=1}^{d}z_{j}\left(\mathcal{S}({\bm{x}}_{j})-\tilde {\bm{H}}{\bm{x}}_{j}\right)\right\|_{2}^{2}$ (C.55) $\displaystyle=\left\|\mathcal{S}\left(\sum_{j=1}^{d}z_{j}{\bm{x}}_{j}\right)- \tilde{\bm{H}}\left(\sum_{j=1}^{d}z_{j}{\bm{x}}_{j}\right)\right\|_{2}^{2}$ (C.56) $\displaystyle\leq\frac{\beta^{2}}{3\lambda}\left\|\tilde{\bm{H}}^{-1}{\bm{X}}{ \bm{z}}\right\|_{\tilde{\bm{H}}}^{2}$ (C.57) $\displaystyle\leq\beta^{2}\|\tilde{\bm{H}}^{-1}{\bm{X}}{\bm{z}}\|_{2}^{2}.$ (C.58) $$

Therefore we have

$$ $\displaystyle(1-\beta)\|\tilde{\bm{H}}^{-1}{\bm{X}}{\bm{z}}\|_{2}\leq\|{\bm{Q} }{\bm{z}}\|\leq(1+\beta)\|\tilde{\bm{H}}^{-1}{\bm{X}}{\bm{z}}\|_{2},$ (C.59) $$

and this is equivalent to

$$ $\displaystyle(1-\beta)^{2}{\bm{z}}^{\top}({\bm{X}}\tilde{\bm{H}}^{-2}{\bm{X}}) {\bm{z}}\leq{\bm{z}}^{\top}{\bm{Q}}^{\top}{\bm{Q}}{\bm{z}}\leq(1+\beta)^{2}{ \bm{z}}^{\top}\left({\bm{X}}^{\top}\tilde{\bm{H}}^{-2}{\bm{X}}\right){\bm{z}}.$ (C.60) $$

Notice that ${\bm{X}}^{\top}\tilde{\bm{H}}^{-2}{\bm{X}}=\tilde{\bm{T}}$. We conclude that ${\bm{Q}}^{\top}{\bm{Q}}\approx_{2\beta+\beta^{2}}\tilde{\bm{T}}$. When $\beta<\frac{1}{8}$, we have ${\bm{Q}}^{\top}{\bm{Q}}\approx_{4\beta}\tilde{\bm{T}}$.

Lastly, from Lemma A.2 and the discussions in Section A.3, we have there exists a constant $C^{\prime}\in(0,1/4)$, such that

$$ $\displaystyle{\bm{P}}=\tilde{\bm{Q}}^{\top}\tilde{\bm{Q}}\approx_{C^{\prime}}{ \bm{Q}}^{\top}{\bm{Q}}$ (C.61) $$

with probability at least $1-\frac{1}{2n}$.

Put the results above together and we get

$$ $\displaystyle{\bm{P}}\approx_{12\beta+C^{\prime}}{\bm{T}}.$ (C.62) $$

Notice that $12\beta+C^{\prime}<\frac{1}{2}$. Using a union bound can prove that the fail probability of this whole process is bounded by $\frac{1}{n}$.

###### Lemma C.5 .

###### Proof.

### C.6 Proof of Lemma 5.4

Let ${\bm{F}}=({\bm{I}}+\lambda\hat{\bm{L}}){\bm{X}}$. We have the Hessian of $\ell^{\prime}$ is

$$ $\displaystyle\nabla^{2}\ell^{\prime}={\bm{P}}^{-\frac{1}{2}}{\bm{F}}^{\top}{ \bm{F}}{\bm{P}}^{-\frac{1}{2}}.$ (C.63) $$

From the condition that ${\bm{P}}\approx_{c_{0}}{\bm{F}}^{\top}{\bm{F}}$, we have ${\bm{F}}^{\top}{\bm{F}}\preceq(1+c_{0}){\bm{P}}$, which implies

$$ $\displaystyle{\bm{P}}^{-1/2}{\bm{F}}^{\top}{\bm{F}}{\bm{P}}^{-1/2}\preceq(1+c_ {0}){\bm{I}}.$ (C.64) $$

Similarly, since ${\bm{P}}\preceq(1+c_{0}){\bm{F}}^{\top}{\bm{F}}$, we have

$$ $\displaystyle\left({\bm{F}}^{\top}{\bm{F}}\right)^{-1}\preceq(1+c_{0}){\bm{P}} ^{-1},$ (C.65) $$

which implies

$$ $\displaystyle\left[\lambda_{\min}({\bm{P}}^{-1/2}{\bm{F}}^{\top}{\bm{F}}{\bm{P }}^{-1/2})\right]^{-1}$ $\displaystyle=\lambda_{\max}\left({\bm{P}}^{1/2}\left({\bm{F}}^{\top}{\bm{F}} \right)^{-1}{\bm{P}}^{1/2}\right)$ (C.66) $\displaystyle\leq 1+c_{0}.$ (C.67) $$

.

Notice that $\ell^{\prime}$ and $\ell$ have the same global optimal value, let it be $\ell^{*}$. For any ${\bm{w}}^{\prime}$ satisfies $\ell^{\prime}({\bm{w}}^{\prime})-\ell^{*}\leq\gamma$, we have

$$ $\displaystyle\ell\left({\bm{P}}^{1/2}{\bm{w}}^{\prime}\right)-\ell^{*}=\ell^{ \prime}({\bm{w}}^{\prime})-\ell^{*}\leq\gamma.$ (C.68) $$

Therefore, if ${\bm{w}}^{\prime}$ is a solution of $\ell^{\prime}$ with $\epsilon$ error rate, then ${\bm{P}}^{1/2}{\bm{w}}^{\prime}$ is a solution of $\ell$ with $\epsilon$ error rate.

### C.7 Proof of Theorem 5.1

As we noted in the main paper, Algorithm 1 is composed by two components: constructing the preconditioner ${\bm{P}}$ and applying it to the optimization.

#### Constructing the Preconditioner.

From Lemma 5.2, the first step that applies spectral sparsifier to get $\tilde{\bm{L}}$ requires $\tilde{O}(m)$ time and the number of non-zero entries in $\tilde{\bm{L}}$ is $\tilde{O}\left(n_{\lambda}\lambda^{2}\right)$. By Lemma A.1, running the SDD solver with error rate $\frac{\beta}{3\lambda}$ in the second step requires $\tilde{O}\left({n_{\lambda}\lambda^{2}}\right)$ time. Notice that since ${\bm{X}}$ is an $n\times d$ matrix, we actually need to run SDD solver for $d$ times, and this introduces another $d$ factor in the time complexity.

As we noted in Section A.3, applying the Hadamard transformation to ${\bm{R}}{\bm{Q}}$ requires $\tilde{O}(d)$ time, and applying the subsampling requires $\tilde{O}(n)$ time. Since $\tilde{\bm{Q}}\in\mathbb{R}^{s\times d}$, calculating ${\bm{P}}=\tilde{\bm{Q}}^{\top}{\bm{Q}}$ requires $O(d^{2}s)$ time. Since $s=\tilde{O}(d)$, this step takes $\tilde{O}(d^{3})$ time. We use brute force to calculate ${\bm{P}}^{\prime}={\bm{P}}^{-1/2}$, and this takes $O(d^{3})$ time.

As a summary, constructing the preconditioner ${\bm{P}}^{\prime}$ takes $\tilde{O}(n_{\lambda}\lambda^{2}d+nd+d^{3})$ time. Notice that we only perform this step once during the whole algorithm.

#### Performing the Iterations.

Next we consider the time required for each iteration. We calculate ${\bm{u}}^{(t)}$ from right to left and it takes $O(nd)$ time. Next, we need to perform two SDD solvers with error rate $\mu$. From Lemma A.1, it takes $\tilde{O}\left[m\log\left(\frac{1}{\mu}\right)\right]$, which is $\tilde{O}\left[m\log\left(\frac{1}{\epsilon}\right)\right]$. Calculating ${\bm{g}}^{(t)}$ and ${\bm{w}}^{(t)}$ is straight-forward and takes $\tilde{O}(nd)$ time.

To conclude, performing each iteration requires $\tilde{O}\left[m\log\left(\frac{1}{\epsilon}\right)+nd\right]$ time. By Lemma 5.1 and Lemma 5.3, with a proper step size, the number of iterations needed for solving outer problem is $\tilde{O}\left(\log 1/\epsilon\right)$.

Combining the analysis above together, to overall complexity is

$$ $\displaystyle\tilde{O}\left[n_{\lambda}\lambda^{2}d+nd+d^{3}+\left(m\log\left( \frac{1}{\epsilon}\right)+nd\right)\log\left(\frac{1}{\epsilon}\right)\right]= \tilde{O}\left(n_{\lambda}\lambda^{2}d+d^{3}+(m+nd)\left(\log 1/\epsilon\right )^{2}\right),$ (C.69) $$

which proves the claim.

## Appendix D Further Discussions

In this section, we extend some of the discussions in the main paper.

### D.1 An Analysis of CE Loss v.s. MSE Loss

Although HERTA is derived from MSE loss, the experiment result shows it also works on CE loss. In this subsection we provide an analysis showing the similarity of the gradient and Hessian of TWIRLS on CE loss and MSE loss, to offer a intuitive explanation why a method that is derived from MSE loss can work on CE loss.

In this section, for a node $u\in\{1,2,\cdots,n\}$, we use ${\bm{y}}^{(u)}\in\mathbb{R}^{c}$ to represent the $u$-th row of ${\bm{Y}}$ (notice we use super-script here to distinguish from ${\bm{y}}_{i}$ used before), ${\bm{h}}_{u}\in\mathbb{R}^{n}$ to represent the $u$-th row of $\left({\bm{I}}+\lambda\hat{\bm{L}}\right)^{-1}$. For $p\in\{1,2,\cdots d\}$, we use ${\bm{x}}_{p}\in\mathbb{R}^{n}$ to denote the $p$-th column of ${\bm{X}}$. For $i\in\{1,2,\cdots c\}$, we use $y^{(u)}_{i}$ to denote the $i$-th entry of ${\bm{y}}^{(u)}$ (or in other words the $(u,i)$-th entry of ${\bm{Y}}$), and $w_{p,i}$ to denote the $(p,i)$-th entry of ${\bm{W}}$.

The MSE loss of TWIRLS defined in eq. 4.7 can be decomposed into summation of sub-losses of each node, i.e.

$$ $\displaystyle\ell({\bm{W}})=\sum_{u=1}^{n}\frac{1}{2}\left\|{\bm{h}}_{u}^{\top }{\bm{X}}{\bm{W}}-{\bm{y}}^{(u)}\right\|_{\mathcal{F}}^{2}=\sum_{u=1}^{n}\ell^ {(u)}\left({\bm{W}}\right),$ (D.1) $$

where $\ell^{(u)}({\bm{W}})=\frac{1}{2}\left\|{\bm{W}}^{\top}{\bm{X}}^{\top}{\bm{h}}_
{u}-{\bm{y}}^{(u)}\right\|_{\mathcal{F}}^{2}$. For a specific class $i\in\{1,2,\cdots,c\}$,
we have the gradient of $\ell^{(u)}$ w.r.t. ${\bm{w}}_{i}$ is

$$ $\displaystyle\frac{\partial\ell^{(u)}}{\partial{\bm{w}}_{i}}={\bm{X}}^{\top}{ \bm{h}}_{u}{\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{i}-{\bm{X}}^{\top}{\bm{h}}_{u} \times{{y}}^{(u)}_{i}.$ (D.2) $$

It is not hard to see from eq. D.2 that the Hessian of $\ell^{(u)}$ with respect to ${\bm{w}}_{i}$ is

$$ $\displaystyle\nabla^{2}_{{\bm{w}}_{i}}\ell^{(u)}\left({\bm{W}}\right)={\bm{X}} ^{\top}{\bm{h}}_{u}{\bm{h}}_{u}^{\top}{\bm{X}}.$ (D.3) $$

For classification tasks, the target ${\bm{y}}^{(u)}$-s are one-hot vectors representing the class of the node. Notice when we calculate cross entropy loss, we use $\mathsf{softmax}$ to normalize it before feeding it into the loss function. In the following for a vector ${\bm{v}}\in\mathbb{R}^{c}$ whose $i$-th entry is $v_{i}$, we define

$$ $\displaystyle\mathsf{softmax}({\bm{v}})_{i}=\frac{\exp(v_{i})}{\sum_{j=1}^{c} \exp(v_{j})}.$ (D.4) $$

Suppose the $u$-th node belongs to class $k$, then the cross entropy loss of the $u$-th node is defined as

$$ $\displaystyle\mathsf{CE}^{(u)}({\bm{W}})$ $\displaystyle=-\log\mathsf{softmax}\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{W}} \right)_{k}$ (D.5) $\displaystyle=-{\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{k}+\log\sum_{j=1}^{c}\exp \left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{j}\right),$ (D.6) $$

The gradient of $\mathsf{CE}^{(u)}$ w.r.t. ${\bm{w}}_{i}$ is

$$ $\displaystyle\frac{\partial\mathsf{CE}^{(u)}}{\partial{\bm{w}}_{i}}$ $\displaystyle=-{\bm{\delta}}_{i,k}\times{\bm{X}}^{\top}{\bm{h}_{u}}+\frac{{\bm {X}}^{\top}{\bm{h}}_{u}\exp\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{i}\right )}{\sum_{j=1}^{c}\exp\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{j}\right)}$ (D.7) $\displaystyle={\bm{X}}^{\top}{\bm{h}}_{u}\mathrm{softmax}\left({\bm{h}}_{u}^{ \top}{\bm{X}}{\bm{W}}\right)_{i}-{\bm{X}}^{\top}{\bm{h}}_{u}\times{{y}}_{i}^{( u)}.$ (D.8) $$

By taking another partial differentiation, we have

$$ $\displaystyle\frac{\partial^{2}\mathsf{CE}^{(u)}}{\partial w_{p,i}\partial w_{ q,i}}$ $\displaystyle=\frac{\left({\bm{x}}_{p}^{\top}{\bm{h}}_{u}\right)\left({\bm{x}} _{q}^{\top}{\bm{h}}_{u}\right)\exp\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{i }\right)}{\sum_{j=1}^{c}\exp\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{j} \right)}-\frac{\left({\bm{x}}_{p}^{\top}{\bm{h}}_{u}\right)\left({\bm{x}}_{q}^ {\top}{\bm{h}}_{u}\right)\left[\exp\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{ i}\right)\right]^{2}}{\left[\sum_{j=1}^{c}\exp\left({\bm{h}}_{u}^{\top}{\bm{X} }{\bm{w}}_{j}\right)\right]^{2}}$ (D.9) $\displaystyle=\left({\bm{x}}_{p}^{\top}{\bm{h}}_{u}\right)\left({\bm{x}}_{q}^{ \top}{\bm{h}}_{u}\right)\left[s_{i}-s_{i}^{2}\right],$ (D.10) $$

where $s_{i}=\mathsf{softmax}\left({\bm{h}}_{u}^{\top}{\bm{X}}{\bm{W}}\right)_{i}$. Rewriting eq. D.10 into matrix form, we have

$$ $\displaystyle\nabla^{2}_{{\bm{w}}_{i}}\mathsf{CE}^{(u)}({\bm{W}})=(s_{i}-s_{i} ^{2}){\bm{X}}^{\top}{\bm{h}}_{u}{\bm{h}}_{u}^{\top}{\bm{X}}.$ (D.11) $$

Comparing eq. D.2 and eq. D.8, we can notice that the gradient of MSE loss and CE loss are essentially the same except the term ${\bm{h}}_{u}^{\top}{\bm{X}}{\bm{w}}_{i}$ in eq. D.2 is normalized by softmax in eq. D.8. Additionally, by comparing the Hessian of MSE loss eq. D.3 and the Hessian of CE loss eq. D.11, we have the latter is a rescaled version of the former. These similarities intuitively explains why our method can be also used on CE loss.

### D.2 The Unavoidable λ 2 superscript 𝜆 2 \lambda^{2} italic_λ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT in the Running Time

Lemma C.5 essentially states the following idea: even when two PD matrices are very similar (in terms of spectral approximation), if they are ill-conditioned, their square can be very different. In this subsection, we prove that the bound obtained in Lemma C.5 is strict up to constant factors by explicitly constructing a worst case. Notice that, in this subsection we discuss the issue of squaring spectral approximators in a broader context, thus in this subsection we possibly overload some symbols used before to simplify the notation.

First we consider another definition of approximation rate: for two PD matrices ${\bm{\Sigma}}$ and $\widetilde{\bm{\Sigma}}$, we define the approximation rate as

$$ $\displaystyle\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)=\max\left \{\left\|{\bm{\Sigma}}^{-1/2}\widetilde{\bm{\Sigma}}{\bm{\Sigma}}^{-1/2}\right \|,\left\|\widetilde{\bm{\Sigma}}^{-1/2}{\bm{\Sigma}}\widetilde{\bm{\Sigma}}^{ -1/2}\right\|\right\}.$ (D.12) $$

This definition is easier to calculate in the scenario considered in this subsection, and can be easily translated to the definition we used in Section 3: when $\epsilon=\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)\in(0,1)$, we have $\widetilde{\bm{\Sigma}}\approx_{\epsilon}{\bm{\Sigma}}$, and when $\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)$ is larger than $1$ we don’t have an approximation in the form defined in Section 3.

Let $\gamma>1$ and $\delta\in(0,1)$ be real numbers. Define

$$ $\displaystyle{\bm{\Sigma}}=\begin{bmatrix}1&0\\ 0&1\end{bmatrix}\begin{bmatrix}\gamma&0\\ 0&1\end{bmatrix}\begin{bmatrix}1&0\\ 0&1\end{bmatrix},$ (D.13) $$

and

$$ $\displaystyle\widetilde{\bm{\Sigma}}=\begin{bmatrix}\sqrt{1-\delta^{2}}&- \delta\\ \delta&\sqrt{1-\delta^{2}}\end{bmatrix}\begin{bmatrix}\gamma&0\\ 0&1\end{bmatrix}\begin{bmatrix}\sqrt{1-\delta^{2}}&\delta\\ -\delta&\sqrt{1-\delta^{2}}\end{bmatrix}.$ (D.14) $$

It’s not hard to compute the error rate $\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)=\left\|{\bm{\Sigma}}^{-
1/2}\widetilde{\bm{\Sigma}}{\bm{\Sigma}}^{-1/2}\right\|\approx\Theta(\gamma
\delta^{2})$. While if we consider the squared matrices, we have $\psi\left({\bm{\Sigma}}^{2},\widetilde{\bm{\Sigma}}^{2}\right)=\left\|{\bm{
\Sigma}}^{-1}\widetilde{\bm{\Sigma}}^{2}{\bm{\Sigma}}^{-1}\right\|\approx
\gamma^{2}\delta^{2}$. In this case, $\psi\left({\bm{\Sigma}}^{2},\widetilde{\bm{\Sigma}}^{2}\right)$ is larger than $\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)$ by a factor of $\gamma$, which is the condition number of ${\bm{\Sigma}}$. When $\gamma$ is very large and $\delta$ is very small, we can have a small $\psi\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)$ while large $\psi\left({\bm{\Sigma}}^{2},\widetilde{\bm{\Sigma}}^{2}\right)$, which .

#### A More General Construction.

For a PD matrix ${\bm{\Sigma}}$ and its spectral approximation $\widetilde{\bm{\Sigma}}$, we call $\frac{\psi\left({\bm{\Sigma}}^{2},\widetilde{\bm{\Sigma}}^{2}\right)}{\psi
\left({\bm{\Sigma}},\widetilde{\bm{\Sigma}}\right)}$ the Squared Error Rate. As noted above, in the worst case the squared error rate can be as large as the condition number of ${\bm{\Sigma}}$. However, the construction above is limited to $2\times 2$ matrices. Now we construct a more general worst case of the squared error rate and perform a loose analysis. Although not rigorously proved, the construction and the analysis suggest the origination of large squared error rates: it approaches the upper bound (condition number) when the eigenspace of the two matrices are very well aligned but not exactly the same.

Let ${\bm{A}}$ be an ill-conditioned matrix with all eigenvalues very large except one eigenvalue equals to $1$, and the eigenvalues of ${\bm{B}}$ are all closed to the eigenvalues of ${\bm{A}}$. Specifically, Let the SVD of ${\bm{A}}$ be ${\bm{A}}={\bm{U}}\mathrm{diag}({\bm{\lambda}}){\bm{U}}^{\top}$, where ${\bm{\lambda}}=\begin{bmatrix}\lambda_{1}&\lambda_{2}&\cdots\lambda_{n}\end{bmatrix}$, and suppose $\lambda_{n}=1$ and $\lambda_{k}\gg 1,\forall k\leq n-1$. For simplicity we just let ${\bm{A}}$ and ${\bm{B}}$ have the same eigenvalues. Suppose the SVD of ${\bm{B}}$ is ${\bm{B}}={\bm{V}}\mathrm{diag}({\bm{\lambda}}){\bm{V}}^{\top}$ and let ${\bm{\Lambda}}=\mathrm{diag}({\bm{\lambda}})$. We denote the one sided error rate of ${\bm{A}}$ and ${\bm{B}}$ by $\epsilon$, i.e. $\epsilon=\left\|{\bm{A}}^{-1/2}{\bm{B}}{\bm{A}}^{-1/2}\right\|$. We have

$$ $\displaystyle\epsilon$ $\displaystyle=\left\|{\bm{\Lambda}}^{-1/2}{\bm{U}}^{\top}{\bm{V}}{\bm{\Lambda} }{\bm{V}}^{\top}{\bm{U}}{\bm{\Lambda}}^{-1/2}\right\|$ (D.15) $\displaystyle=\left\|{\bm{\Lambda}}^{-1/2}{\bm{W}}{\bm{\Lambda}}{\bm{W}}^{\top }{\bm{\Lambda}}^{-1/2}\right\|,$ (D.16) $$

where ${\bm{W}}={\bm{U}}^{\top}{\bm{V}}$.

Since $\lambda_{n}=1$ and for all $k\leq n-1$, $\lambda_{k}\gg 1$, we have

$$ $\displaystyle{\bm{\Lambda}}^{-\frac{1}{2}}=\begin{bmatrix}\lambda_{1}^{-1/2}&& &\\ &\lambda_{2}^{-\frac{1}{2}}&&\\ &&\ddots&\\ &&&\lambda_{n}^{-\frac{1}{2}}\end{bmatrix}\approx\begin{bmatrix}0&&&\\ &0&&\\ &&\ddots&\\ &&&1\end{bmatrix},$ (D.17) $$

and therefore

$$ $\displaystyle\epsilon$ $\displaystyle=\left\|{\bm{\Lambda}}^{-1/2}{\bm{W}}{\bm{\Lambda}}{\bm{W}}^{\top }{\bm{\Lambda}}^{-1/2}\right\|\approx({\bm{W}}{\bm{\Lambda}}{\bm{W}}^{\top})_{ n,n}=\sum_{k=1}^{n}{\bm{W}}_{k,n}^{2}\lambda_{k},$ (D.18) $$

where ${\bm{A}}_{i,j}$ represents the $(i,j)$-th entry of a matrix ${\bm{A}}$.

For simplicity, we assume all $\lambda_{k}(k\leq n-1)$ are approximate equal, i.e. $\lambda_{1}\approx\lambda_{2}\cdots\approx\lambda_{n-1}=\gamma$. Then we have

$$ $\displaystyle\epsilon$ $\displaystyle=\sum_{k=1}^{n}{\bm{W}}_{k,n}^{2}\lambda_{k}$ (D.19) $\displaystyle\approx\gamma\left[\sum_{k=1}^{n-1}{\bm{W}}_{k,n}^{2}+\frac{1}{ \gamma}{\bm{W}}_{n,n}^{2}\right]$ (D.20) $\displaystyle\approx\gamma\sum_{k=1}^{n-1}{\bm{W}}_{k,n}^{2}$ (D.21) $\displaystyle=\gamma\left(1-{\bm{W}}_{n,n}^{2}\right).$ (D.22) $$

Recall ${\bm{W}}={\bm{U}}^{\top}{\bm{V}}$, we have ${\bm{W}}_{n,n}={\bm{u}}_{n}^{\top}{\bm{v}}_{n}$, where ${\bm{u}}_{n}$ and ${\bm{v}}_{n}$ are the $n$-th column of ${\bm{U}}$ and ${\bm{V}}$ respectively. Thus we have

$$ $\epsilon\approx\gamma\left[1-\left({\bm{u}}_{n}^{\top}{\bm{v}}_{n}\right)^{2} \right].$ (D.23) $$

Now we can see the spectral approximation error $\epsilon$ is determined by two terms: the condition number $\gamma$ and the matchness of the eigenvectors corresponds to small eigenvalues, which is evaluated by $\left[1-\left({\bm{u}}_{n}^{\top}{\bm{v}}_{n}\right)^{2}\right]$. When $\gamma$ is very large but $\left[1-\left({\bm{u}}_{n}^{\top}{\bm{v}}_{n}\right)^{2}\right]$ is small, ${\bm{B}}$ can still be a spectral approximation of ${\bm{A}}$. For example if $(1-{\bm{u}}_{n}^{\top}{\bm{v}}_{n})\approx\gamma^{-1}$, we can get a spectral approximation error $\epsilon\approx 1$. However, after we square the matrices, the eigenvalues will also get squared, but eigenvectors remains unchanged. That will enlarge the spectral approximation error by a factor of $\gamma$, i.e. $\left\|{\bm{A}}^{-1}{\bm{B}}^{2}{\bm{A}}^{-1}\right\|\approx\gamma^{2}\gamma^{
-1}=\gamma$, which becomes very large. In this case the squared error rate is $\gamma$, the condition number of ${\bm{A}}$.

### D.3 Applying Graph Sparsification in Each Iteration

As mentioned in the main paper, unlike most existing work, in our algorithm we don’t sparsify the graph in each training iteration. The proof of Lemma 5.1 (see Section C.3) suggests the reason why performing graph sparsification in each iteration can lead to suboptimal running time. In this subsection we illustrate this claim in detail. The intuition is that, in order to ensure convergence of the training, we require a small error rate in the gradient estimation, which is of the order $\epsilon^{1/2}$ as we have showed in Section C.3. It is acceptable for the SDD solver because the running time of the SDD solver only logarithmly depends on the error rate. However, if we sparsify the graph, to obtain an $\epsilon^{1/2}$ error rate we will need to sample $O(n_{\lambda}/\epsilon)$ edges, which grows linearly with $1/\epsilon$, and can be large especially when $\epsilon$ is very small. Below is a more detailed analysis.

Consider in the for-loop of Algorithm 1 we replace the $\hat{\bm{L}}$ by a sparsified version ${\bm{L}}^{\prime}=\mathsf{Sparsify}_{\omega}\left(\hat{\bm{L}}\right)$, where $\omega\in(0,1)$. Let ${\bm{H}}^{\prime}={\bm{I}}+\lambda{\bm{L}}^{\prime}$. From Lemma 5.2, we have ${\bm{H}}^{\prime}\approx_{\omega}{\bm{H}}$. Consider $E_{1}$ defined in Section C.3, it now becomes

$$ $\displaystyle E_{1}^{\prime}$ $\displaystyle=\left\|{\bm{X}}^{\top}\left[\mathcal{S}\left({\bm{H}}^{\prime-1} {\bm{z}}\right)-{\bm{H}}^{-2}{\bm{z}}\right]\right\|$ (D.24) $\displaystyle\leq\sigma_{\max}\left({\bm{X}}\right)\left\|\mathcal{S}\left({ \bm{H}}^{\prime-1}{\bm{z}}\right)-{\bm{H}}^{-2}{\bm{z}}\right\|.$ (D.25) $$

Here we can not proceed by using the fact that $\mathcal{S}$ is a linear solver, because now $\mathcal{S}$ is not a linear solver for ${\bm{H}}$ but for ${\bm{H}}^{\prime}$, thus we will have to split the error term again:

$$ $\displaystyle E_{1}^{\prime}\leq\sigma_{\max}\left({\bm{X}}\right)\left\| \mathcal{S}\left({\bm{H}}^{\prime-1}{\bm{z}}\right)-{\bm{H}}^{\prime-2}{\bm{z} }\right\|+\sigma_{\max}({\bm{X}})\left\|{\bm{H}}^{-2}{\bm{z}}+{\bm{H}}^{\prime -2}{\bm{z}}\right\|.$ (D.26) $$

Of the two terms on the right-hand side of eq. D.26, the first one can be bounded with a similar method used in Section C.3, and the second term is bounded by

$$ $\displaystyle\sigma_{\max}({\bm{X}})\left\|{\bm{H}}^{-2}{\bm{z}}-{\bm{H}}^{ \prime-2}{\bm{z}}\right\|$ $\displaystyle\leq\sigma_{\max}\left({\bm{X}}\right)\left\|I-{\bm{H}}^{-1}{\bm{ H}}^{\prime-2}{\bm{H}}^{-1}\right\|\|{\bm{H}}^{-2}{\bm{z}}\|$ (D.27) $\displaystyle\leq\sigma_{\max}\left({\bm{X}}\right)\omega\left\|{\bm{H}}^{-2}{ \bm{z}}\right\|$ (D.28) $\displaystyle\leq\sqrt{\frac{8}{\epsilon}}\kappa({\bm{X}})\lambda\omega\left\| \nabla\ell_{i}\left({\bm{w}}\right)\right\|.$ (D.29) $$

Therefore, in order to obtain $E_{1}^{\prime}\leq\|\nabla\ell_{i}({\bm{w}})\|$, we at least require $\omega\leq\frac{\epsilon^{-1/2}}{\sqrt{8}\kappa({\bm{X}})\lambda}=O\left(
\epsilon^{-1/2}\right)$. From Lemma 5.2, the number of edges in the sparsified graph $\omega$ error rate is $O(n_{\lambda}/\epsilon)$. A similar analysis can be also applied to $E_{2}$, and by repeating the proof in Section C.3, we obtain an overall running time bound

$$ $\displaystyle O(m+n_{\lambda}\lambda^{2}d+d^{3}+n_{\lambda}\epsilon^{-1}\log 1 /\epsilon),$ (D.30) $$

which, although eliminates the $m\left(\log\frac{1}{\epsilon}\right)^{2}$ term, introduces an extra $n_{\lambda}\epsilon^{-1}\log\frac{1}{\epsilon}$ term, and is usually worse than the original bound we derived in Theorem 5.1, especially when requiring a relatively small $\epsilon$.

That being said, although not very likely, when $\epsilon$ and / or $\lambda$ are large, it is possible that directly sparsifying the graph in each iteration is beneficial. Considering this, we can adopt a mixed strategy: when $n_{\lambda}\epsilon^{-1}\leq m\log(1/\epsilon)$, we sparsify the graph in each iteration, elsewise we don’t. This leads to the following overall running time bound:

$$ $\displaystyle O\left[m+n_{\lambda}\lambda^{2}d+d^{3}+\min\left\{m\log\frac{1}{ \epsilon},n_{\lambda}\epsilon^{-1}\right\}\log\frac{1}{\epsilon}\right],$ (D.31) $$

which is slightly better than the one we presented in Theorem 5.1.

### D.4 Possible Directions for Extending HERTA to More Complex Models

The assumption that $f({\bm{X}};{\bm{W}})={\bm{X}}{\bm{W}}$ provides us with a convenience that the Hessian of this model is a constant matrix. Therefore in Algorithm 1 we only need to calculate the preconditioner ${\bm{P}}$ for one time. However, if $f({\bm{X}};{\bm{W}})$ is implemented by a non-linear network, then the Hessian will change by the time and might be hard to calculate, which will be a key challenge to using a more complex $f$. We note that this can be possibly addressed by constructing a linear approximation of $f$ using its Jaccobian. In each iteration, we can use the Jaccobian to replace the ${\bm{X}}$ used in current version of HERTA, and recalculate ${\bm{P}}$ at each iteration. Since the convergence is fast, we only need to recalculate the Jaccobian for a small number of iterations, so should not bring massive change to the running time of the algorithm.

The attention mechanism of TWIRLS in (Yang et al., ) is achieved by adding a concave penalty function to each summand of the the $\mathrm{Tr}\left({\bm{Z}}^{\top}\hat{\bm{L}}{\bm{Z}}\right)$ term in eq. 4.1. For specific penalty functions, it might still possible to find a inner problem solver, as long as the problem stay convex. We note that this depends on concrete implementation of the penalty term used. Investigating how to fast solve the inner problem under various penalty functions should also be an important problem for future study.

## Appendix E Implementation Details

In the experiments, the datasets are loaded and processed using the DGL package (Wang et al., ). We use the original inner problem solver in (Yang et al., ) since it is not computation bottleneck. It can be in principle replaced by any implementation of SDD solvers.

For the calculating the gradient of the preconditioned model, we presented a calculation method in Algorithm 1 which maintains the lowest computational complexity. In preliminary experiments, we tested this calculation method with using the autograd module in pytorch(^4^44https://pytorch.org/) and verified that they have the same output, and similar computational efficiency on real world datasets (again in practice this step is not a bottleneck). Therefore, we simply use pytorch autograd module to compute gradients in experiments.

Using the pytorch autograd module also enables us to apply HERTA on various loss functions and optimizers: we only need to perform the preconditioning and indicate the loss function. The gradient and optimization algorithm with be automatically realized by pytorch.

In order to ensure a fair comparison and prevent confounding factors, we don’t using any common training regularization techniques like weight decay or dropout. For each setting, we repeat the experiment with learning rates in $\{0.001,0.01,0.1,1,10\}$ and choose the trial which the training loss does not explode and with the lowest final training loss to report in the main paper.

Figure: Figure 3: The training loss comparison between HERTA and standard optimizers on MSE loss with $\lambda=20$. Dataset used from left to right: ogbn-arxiv, citeseer, pubmed.
Refer to caption: extracted/2403.18142v1/figures-log/ogbn-arxiv-20-fast,mse.png

Figure: Figure 4: The training loss comparison between HERTA and standard optimizers on CE loss with $\lambda=20$. Dataset used from left to right: ogbn-arxiv, citeseer, pubmed.
Refer to caption: extracted/2403.18142v1/figures-log/ogbn-arxiv-20-fast_ce,ce.png

Figure: Figure 5: The training loss comparison between HERTA and standard optimizers on Cora with $\lambda=1$. Left: CE loss. Right: MSE loss.
Refer to caption: extracted/2403.18142v1/figures-log/cora-1-fast_ce,ce.png

## Appendix F Additional Experiment Results

In this section, we present additional experiment results.

### F.1 Experiments with Larger λ 𝜆 \lambda italic_λ

All the experiments presented in the main paper are with $\lambda=1$. In this subsection, we present results with $\lambda=20$. See Figures 3 and 4 for results with MSE loss and CE loss respectively. The results supports our observation in the main paper that HERTA works consistently well on all settings.

### F.2 Additional Results on Cora

In this subsection we present results on Cora, another citation dataset used in (Yang et al., ). See Figure 5 and Figure 6 for results with $\lambda=1$ and $\lambda=20$ respectively. It is clear that the results on Cora is consistent with our observation on other datasets.

Figure: Figure 6: The training loss comparison between HERTA and standard optimizers on Cora with $\lambda=20$. Left: CE loss. Right: MSE loss.
Refer to caption: extracted/2403.18142v1/figures-log/cora-20-fast_ce,ce.png

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 5
presentation: 4
questions: Despite these issues, the paper has some valuable contributions. I recommend that the authors provide more comprehensive experiments and discuss the limitations and complexity properly. The theoretical contradictions must be resolved and clear practical guidelines on how to use this method should be presented. 
The core ideas are promising, but the current execution has too many significant weaknesses to warrant acceptance, even as a weak accept. The combination of theoretical inconsistencies and experimental limitations suggests the work needs substantial revision before it's ready for publication.
soundness: 3
strengths: 1. Provides a novel spectral sparsification method for normalized/regularized Laplacians
2. Shows promising empirical results on larger datasets
3. Works with different loss functions and optimizers
4. Provides a solid foundation for future work in this direction
summary: This paper introduces HERTA, a novel training algorithm for Unfolded Graph Neural Networks (GNNs) that aims to solve scalability issues while preserving interpretability. The algorithm achieves nearly linear time training guarantees and introduces a new spectral sparsification method for normalized and regularized graph laplacians. The authors demonstrate HERTA's effectiveness through theoretical analysis and empirical validation on real-world datasets.
weaknesses: 1. The main idea (on lines 211-217) behind the effective Laplacian dimension ($n_\lambda$) is that it represents how many eigenvalue directions "matter significantly" in the energy function. Smaller $n_\lambda$ means fewer “important directions” to consider, which in turn gives less computational complexity. This is explained via the closed form solution in Eqn 4.6 . 
When the eigenvalue $\lambda_i$ is much smaller than $1/\lambda$, then the effect of the Laplacian is dominated by the $I$ (identity matrix). Larger $n_\lambda$ means more work during training as there is a need to preserve graph structure in more directions.  

But the problem arises where for their theoretical guarantees, their algorithm requires $O(n/\lambda^2$ “large eigenvalues”, which implies a high effective dimension. This creates a contradiction because the authors require a low effective dimension for efficiency, but their algorithm's guarantees need a high effective dimension. 
Please elaborate on how the algorithm's performance changes with different distributions of eigenvalues and provide examples of real-world graphs where this condition is satisfied.

2. Also, if we were to consider $\lambda$ to be an “absolute constant” (as written in lines 80-81), then the Big-Oh notation would reduce to $O(n)$ and not be $O(n/\lambda^2)$. In the experiments, the authors set $\lambda$ to 1 (and 20 in their appendix), which again means $O(n)$ in essence. So, you would need a linear fraction of your eigenvalues to be large eigenvalues. Is this practical? 
In real-world scenarios, you typically find very few large eigenvalues and the majority of eigenvalues as small (closer to 0). 
Please discuss the typical distribution of eigenvalues in real-world graphs and how this affects the performance of their algorithm.

3. In theorem 5.1, equation 5.1, the authors have a $d^3$ term in the final time complexity analysis equation. In their informal Thm 1.1 (lines 72-73), they say “In practice, the node feature dimensionality d is generally much smaller than the graph size m, in which case, the running time of HERTA is …” And they ignore the $d^3$ term.  

If we take a closer look at this assumption, there arise several problems with it. Modern ML models often use very high-dimensional features. E.g. Word2vec embeddings have dimensions $d$ around 100-300, BERT embeddings have $d=768$, and LLMs with $d=1024$ or larger. 
Let’s say for d=300, $d^3$ is 27M, for d=768, $d^3$ is approx. 453 M and for d=1024, it's around 1B!  

Even if $d \le m$, $d^3$ can dominate. For m = 10k, d=100, $d^3$ is 1M! The cubic growth makes d^3 significant in the analysis and cannot simply be ignored. 
Please justify and provide conditions for when $d^3$ can be ignored. The only time is when $m >> d$, but this would be way too restrictive a condition for modern applications.

In my opinion, this term should remain and be analyzed properly to say when each term dominates and discuss the trade-offs with feature dimensionality.

- Provide a more detailed analysis of when each term in the complexity bound dominates.
- Discuss the trade-offs between feature dimensionality and algorithm performance.
- Consider revising your informal theorem statement to include the $d^3$ term and provide a more nuanced discussion of when it can be ignored.
- Provide empirical results showing how the algorithm's performance scales with feature dimensionality.

**Experimental Section related issues.**

1. The paper doesn't properly justify why these specific datasets (Cora, Citeseer, Pubmed, ogbn-arxiv, and Flickr) were chosen beyond mentioning that some were used in Yang et al. (2021b). There's no discussion of dataset characteristics (like sparsity, feature dimensions, etc.) that would help understand how these choices stress-test different aspects of the algorithm.

2. There needs to be some analysis that shows the distribution of “large eigenvalues” in the graphs that have been chosen. This would help get more insights and better understand the theoretical bounds proposed on large eigenvalues as well. It can be the case that the theoretical bounds aren’t tight enough to explain why your method gets faster convergence. 

3. The experiments only compare against gradient descent and Adam optimizers. Other relevant baselines like SGD with momentum, RMSprop, or specialized GNN optimizers are missing. This limits understanding of the proposed method’s relative performance against the full spectrum of optimization approaches.

4. In Section 6.1, Figures 1 and 2, can you please explain why a log scale was chosen and how it affects the interpretation of the results. There is no justification given in the paper as to why this was necessary.

5. Despite claiming theoretical time complexity benefits, there are no direct runtime comparisons or scaling experiments showing how HERTA's performance changes with graph size, feature dimensionality, or other relevant parameters. It would be nice, even on a synthetic graph dataset, to see how increasing the graph size (and also feature dimensionality) will affect this method’s convergence runtimes. 

6. As mentioned in Section F3 “on small datasets (Cora and Citeseer with <4000 nodes), HERTA's test accuracy advantages are less clear since optimization isn't a major obstacle." This suggests that this method may not have a universal utility.

7. There are no ablation experiments showing the relative importance of different components of HERTA (like the preconditioner construction, spectral sparsification, etc.). This makes it hard to understand which aspects are crucial for performance improvements. Could it be that the graphs considered already have a good gap between their first and second eigenvalues and thus converging well? Please provide more empirical evidence. 

8. The theoretical analysis focuses primarily on a simple linear implementation of f(X;W) = XW, leaving more complex model architectures unexplored. In the empirical studies, it would be interesting to at the least show performance of the method on more complex deep learning architectures for f.

9. There is a limited analysis of hyperparameter effects and sensitivity analysis. A systematic stability analysis would be very useful. Also, how does one go about selecting the parameters? 

These weaknesses suggest that while HERTA shows promising results, more comprehensive experimental evaluation would strengthen the paper's empirical claims and our understanding of when and why the method works best.

- Can you provide a justification for the choice of datasets and discuss their relevant characteristics (sparsity, feature dimensions, etc.) that make them suitable for evaluating HERTA?

- Please include an analysis of the eigenvalue distribution for the chosen graphs and discuss how this relates to the theoretical bounds.

- Consider expanding the comparison to include other relevant optimizers like SGD with momentum, RMSprop, and specialized GNN optimizers.

- Explain the rationale behind using a log scale in Figures 1 and 2, and discuss how this affects the interpretation of the results.

- Include direct runtime comparisons and scaling experiments

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: - In Section 3, $||M||$ for a matrix $M$ is mentioned as its operator norm.
    What does that mean?
- Tiny Questions (as the main purpose is to accelerate) (while the paper claims that HERTA has excellent generalization properties)
  - There is no accuracy or variance to show the effectiveness of the approach.
  - The performance of the approach is not compared with SOTA baselines.
soundness: 2
strengths: HERTA achieves good results on some benchmark datasets like ogbn-arxiv compared with standard optimizers.
summary: The paper introduces HERTA, a new Graph Neural Network that is designed based on Unfolded GNN.

The approach is designed to accelerate the training process and to ensure convergence.

Moreover, a new spectral sparsification method is also proposed.
weaknesses: - In Section 1, the author claims that the running time is $\widetilde{O}((m+nd)(\log{\frac{1}{\epsilon}})^2+d^3)=\widetilde{O}((m+nd)(\log{\frac{1}{\epsilon}})^2)$ as $ d << m$ in practice, where the feature dimensionality is $d$ and the edge number is $m$.

    However, this does not mean that $\widetilde{O}(d^3)<\widetilde{O}(m)$.
    Moreover, in most benchmark datasets, $d^3 >> m$. Take PubMed for example (the paper also conducts experiments on PubMed), $d = 500$ and $m = 44324$, we have $d^3 > 2000m$. Even in a large dataset, ogbn-arxiv, with a relatively small $d = 128$, and a large $m = 1,166,243$, it also holds that $d^3 \approx 1.8m$.

    This shows that the estimation of the running time of the approach is inappropriate.

- Moreover, no experiment on time and space burden is conducted to validate the efficiency of the approach.

- In Theorem 1.1, the number of large eigenvalues of the graph Laplacian is assumed as $O(n/\lambda^2)$.
    While in Definition 5.1, $n_\lambda$ is roughly the number of large eigenvalues, i.e., $O(n/\lambda^2)$.
    It's not true that $n_\lambda\rightarrow n$ as $\lambda\rightarrow\infty$

- In Definition 5.1, $n_\lambda$ roughly represents the number of eigenvalues of $\hat{L}$ which are of order $\Omega(\lambda^{-1})$.
    While the author also claims that $n_\lambda$ represents the number of eigen-directions for which the Laplacian regularizer $\frac{\lambda}{2}\text{Tr}(Z^\top\hat{L}Z)$ is large, where $Z\in\mathbb{R}^{n\times c}$ is the optimization variable.
    What is the relationship between the eigenvalues of $\hat{L}$ and that of $Z^\top\hat{L}Z$?
    I believe that they are not the same.

- The train loss of the approach is displayed in the form of $\log(\text{train loss}-\text{min})$, which enlarges the change of the loss (as the result is negative for most of the time).
    Moreover, the paper only shows the relative change of each method, the exact minimum value is not shown.
    This adds to the concern of the convergence of the method.

- In Appendix C, (C.3) is wrongly derived:
    Firstly, $g^{(t)}$ and $\nabla l(\omega^{(t)})$ in which should be exchanged. (otherwise (C.5) and (C.8) will not hold)
    Secondly, the derivation from the condition to (C.3) only holds for $l_1$-norm.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 1
questions: Please refer to the Weaknesses section.
soundness: 3
strengths: The HERTA algorithm is built on thorough and rigorous theory.

The technical results are detailed and extensive.

The authors conduct experiments on various real-world datasets, and the presented results display convincing performance.
summary: Current training algorithms for unrolled GNNs suffer from slow iterations and slow convergence rates, and solutions for these issues often require the modification of the original model.

The authors bridge this gap by introducing HERTA, an algorithm to efficiently train Unfolded GNNs, without modifying the original unrolled model.

Rigorous convergence and complexity analysis are provided, along with various convincing experiments on real-world datasets.
weaknesses: I will provide here some suggestions on the presentation as well as technical questions.

* I find Theorem 1.1 to be misplaced in the Introduction. It references its formal version and an objective which come much later. We do not have prior exposure to the manipulated objects or the task at hand at that point. More generally speaking, my suggestion would be to limit the level of technicalities in the introduction which hinders the exposure and induces too many references to later sections (lines 66 to 83). Instead of stating results of complexity analysis, I think it would rather be more useful to informally discuss the core ideas behind the HERTA algorithm (which is presented as the core contribution in the introduction) and/or present the Unfolded GNN that is considered in this work (Eq. 4.1).

* It is also somewhat frustrating to go through convergence results and complexity analysis before even being exposed to the algorithm or at least being presented its core ideas. The HERTA algorithm is presented as being the core contribution in the introduction but we have to wait until Section 5.4 (page 8) to finally be exposed to it for the first time. Concretely speaking, I would provide an overview of the algorithm, explaining the general ideas, before discussing Theorem 5.1.

* Overall, in my opinion, the organisation of the paper is somewhat unusual. This is supported, for example, by the numerous references to later sections or paragraphs (for instance in lines 66, 67, 82, 247, 257) which makes the exposure more difficult. As mentioned earlier, I believe the paper would gain in clarity if the HERTA algorithm is introduced (at least informally) earlier on, and before delving into the theoretical analysis.

* I find the formulation "it is not hard to see that [...]" to be overused (lines 161, 216, 803, 1247). Perhaps you should provide brief explanations instead ?

* In line 166, it is not clear whether "optimizing the following objective" refers to its minimization or maximization. In my opinion, it should be specified that the objective is being minimized, as explicited later on.

* The authors do not provide the code for the implementation of their algorithm or their experiments, and it is not clear in the paper if it would be made public upon acceptance. I think the authors should include a statement about code availability.

* HERTA is essentially an optimization procedure for solving the bi-level optimization problem TWIRLS [Yang et al.,2021]. From a graph signal processing point of view, we are then trying to learn a signal representation $Z^*$ on a given graph $\hat{L}$, parametrized by weights $W$, such that it minimizes a distance $l$ to the label signals given by $Y$. The term $Tr(Z^\top \hat{L} Z)$ in the energy function enforces a smooth representation of the signals on the graph. In the context considered in this paper, the graph represented by its Laplacian $\hat{L}$ is considered given. Would it be possible to extend TWIRLS and HERTA to the graph learning paradigm, where in addition to learning a proper signal representation $Z^*$, we are also seeking to learn the graph topology, in other words to learn both $Z^*$ and $\hat{L}$ ? How would this reflect on HERTA ? Perhaps could you discuss this potential extension in the future work section ?

* In most unrolled optimization procedures, the original algorithm's hyperparameters (such as the step-size or the regularization term $\lambda$) are learned in a data-driven fashion. From my understanding, this was not pursued in this work, could you comment on this choice, and how this would reflect on HERTA ?

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - When comparing with standard GNNs, could you list the specific benefits of Unfolded GNNs, such as training speed, accuracy, interpretability?

- In the demonstration of Theorem 5.1, the phrase 'with a proper step size $\eta$ and number of iterations $T$ seems informal. Could you clarify or formalize this expression?
soundness: 3
strengths: - The paper addresses a critical limitation in Unfolded Graph Neural Networks, which is high training costs, by proposing a novel training algorithm, HERTA, that achieves near-linear worst-case convergence without compromising the interpretability of the model.

- The work is theoretically robust and supported by empirical validation on real-world datasets in context of solving the training problem of Unfolded Graph Neural Networks

-  The paper is overall well-organized, with a clear problem statement and motivation for the research
summary: The paper presents HERTA, an algorithm that optimizes the training of Unfolded Graph Neural Networks (GNNs), which are favored for their interpretability and flexibility but face scalability issues due to high training costs. HERTA addresses two key challenges: slow iteration rates and slow convergence. Unlike previous approaches that compromise interpretability for efficiency, HERTA achieves near-linear worst-case convergence while preserving the model's structure. The algorithm introduces a novel spectral sparsification technique that approximates the inverse of normalized and regularized Laplacian matrices, accelerating the training process without altering the underlying model. Experimental results on real-world datasets validate HERTA's ability to improve training speed and adapt to various loss functions and optimizers, showcasing its versatility and effectiveness for practical applications.
weaknesses: - Minor issue with citation format; please ensure correct use of \cite and \citep commands throughout the paper.

-  The definition of the H-norm in Equation 5.2 is unclear. Providing a brief clarification would improve readability.

- Minor wording issue: the third line of Lemma 5.1 appears to be missing a 'for' in its description, which affects clarity.

- The paper’s focus is on Unfolded Graph Neural Networks, yet a comparison with other types of GNNs could further highlight the advantages of Unfolded GNNs over alternative architectures.

- Although HERTA is shown to outperform standard optimizers, it would be beneficial to include additional baseline methods from within the family of Unfolded Graph Neural Networks (papers cited in your work). This would provide a clearer comparison and better highlight HERTA’s specific advantages within this framework.