# Asymptotic spectra of large (grid) graphs with a uniform local structure
(part II): numerical applications

Andrea Adriani †,(a), Davide Bianchi †,(b), Paola Ferrari †,(c)
Stefano Serra-Capizzano †,(c,d)

###### Abstract

In the current work we are concerned with sequences of graphs having a grid
geometry, with a uniform local structure in a bounded domain
$\Omega\subset{\mathbb{R}}^{d}$, $d\geq 1$. When $\Omega=[0,1]$, such graphs
include the standard Toeplitz graphs and, for $\Omega=[0,1]^{d}$, the
considered class includes $d$-level Toeplitz graphs. In the general case, the
underlying sequence of adjacency matrices has a canonical eigenvalue
distribution, in the Weyl sense, and it has been shown in the theoretical part
of this work that we can associate to it a symbol $\boldsymbol{\mathfrak{f}}$.
The knowledge of the symbol and of its basic analytical features provides key
information on the eigenvalue structure in terms of localization, spectral
gap, clustering, and global distribution. In the present paper, many different
applications are discussed and various numerical examples are presented in
order to underline the practical use of the developed theory. Tests and
applications are mainly obtained from the approximation of differential
operators via numerical schemes such as Finite Differences (FDs), Finite
Elements (FEs), and Isogeometric Analysis (IgA). Moreover, we show that more
applications can be taken into account, since the results presented here can
be applied as well to study the spectral properties of adjacency matrices and
Laplacian operators of general large graphs and networks, whenever the
involved matrices enjoy a uniform local structure.

$\dagger$ All authors have contributed equally.
(a) Department of Theoretical and Applied Sciences, University of Insubria,
Via Dunant 3, 21100 Varese, Italy.
<EMAIL_ADDRESS>
(b) School of Science, Harbin Institute of Technology (Shenzhen), University
Town of Shenzhen, 518055 Shenzhen, China.
<EMAIL_ADDRESS>
(c) Department of Science and High Technology, University of Insubria, Via
Valleggio 11, 22100 Como, Italy.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
(d) Department of Information Technology, Uppsala University, Uppsala, Sweeden
<EMAIL_ADDRESS>

##  1 Introduction

In [1] a large class of (sequences of) graphs has been defined, according to
the following structural properties:

a)

when we look at them from “far away”, they should reconstruct approximately a
given domain $\Omega\subset[0,1]^{d}$, $d\geq 1$, i.e., the larger is the
number of the nodes the more accurate is the reconstruction of $\Omega$;

b)

when we look at them “locally”, that is from a generic internal node, we want
that the structure is uniform, i.e., we should be unable to understand where
we are in the graphs, except possibly when the considered node is close enough
to the boundaries of $\Omega$.

In particular, the domain $\Omega$ and the internal structure of the graphs in
the sequence are fixed, independently of the index (or multi-index) of the
graph uniquely related to the cardinality of nodes: it follows that the
resulting sequence of graphs has a grid geometry, with a uniform local
structure, in a bounded domain $\Omega\subset{\mathbb{R}}^{d}$, $d\geq 1$. The
domain $\Omega$ is assumed to be Lebesgue measurable with regular boundary,
that is a boundary $\partial\Omega$ of zero Lebesgue measure, and it is
supposed to be contained in the hypercube $[0,1]^{d}$. Such a domain is called
regular. When $\Omega=[0,1]$, it has been proven in [1] that such graphs
include the standard Toeplitz graphs (see [33]) and for $\Omega=[0,1]^{d}$ the
considered class includes $d$-level Toeplitz graphs.

Moreover, given a sequence of graphs having a grid geometry with a uniform
local structure in a domain $\Omega$, the underlying sequence of adjacency
matrices shows a canonical eigenvalue distribution, in the Weyl sense (see
[11, 32] and references therein), and it is possible to associate to it a
symbol function $\boldsymbol{\mathfrak{f}}$. More precisely, when
$\boldsymbol{\mathfrak{f}}$ is smooth enough, if $N$ denotes the size of the
adjacency matrix (i.e. the number of nodes of the graph), then the eigenvalues
of the adjacency matrix are approximately values of a uniform sampling of
$\boldsymbol{\mathfrak{f}}$ in its definition domain, which depends on
$\Omega$ (see Definition 2.5 for the formal definition of eigenvalue
distribution in the Weyl sense; the results in Section 5 of [1] for the
precise characterization of $\boldsymbol{\mathfrak{f}}$ and of its definition
domain; and [4] for the definition of its monotone rearrangement).

The knowledge of the symbol and of some of its basic analytical features
provides a lot of information on the eigenvalue structure, of localization,
spectral gap, clustering, and distribution type.

The mathematical tools are mainly taken from the field of Toeplitz structures
(see [11] and [32, 46, 47]) and of Generalized Locally Toeplitz (GLT) matrix-
sequences (see [41, 42, 45]): for a recent account on the GLT theory and for
several useful applications in the approximation of partial differential and
fractional differential operators see [7, 9, 22, 27, 28, 29, 5, 6, 31] and
references therein.

Since in this paper we are interested in numerical applications of the theory
developed in [1], we first show that many numerical schemes (see e.g. [13, 14,
44]) for approximating partial differential equations (PDEs) and operators
lead to sequences of structured matrices, which can be written as linear
combination of adjacency matrices, associated with sequences of graphs showing
a uniform local structure. More specifically, if the physical domain of the
differential operator is $[0,1]^{d}$ (or any $d$-dimensional rectangle) and
the coefficients are constant, then we encounter $d$-level (weighted) Toeplitz
graphs, when approximating the underlying PDE by using e.g. equispaced Finite
Differences (FDs) or uniform Isogeometric Analysis (IgA). On the other hand,
under the same assumptions on the underlying operator, quadrangular and
triangular Finite Elements (FEs) lead to block $d$-level Toeplitz structures,
where the size of the blocks is related to the degree of the polynomial space
of approximation and to the dimensionality $d$ (see [30]). Finally, in more
generality, the GLT case is encountered by using any of the above numerical
techniques, also with non-equispaced nodes/triangulations, when dealing either
with a general domain $\Omega$ or when the coefficients of the differential
operator are not constant. The given classification of approximated PDE
matrix-sequences is relevant also because the obtained spectral information
can be used for guiding the design of proper iterative solvers (in terms
either of preconditioners or of ad hoc multigrid methods) for the underlying
linear systems with large matrix size. In this direction, we provide a series
of numerical examples where the knowledge of the symbol allows to obtain
preconditioners for the conjugate gradient method, which significantly improve
the computational cost for solving the involved large linear systems, and
guides the choice of optimal projectors for two-grid and multigrid methods.

The paper is organized as follows. In Section 2 we will first review some
basic definitions and notation from graph theory, from the field of Toeplitz
and $d$-level Toeplitz graphs, and we then provide the definitions of
canonical spectral distribution, graph Laplacian etc. We also briefly discuss
some properties and tools useful in the study of multigrid methods. Section 3
will constitute the main core of our work. Here, a number of numerical
applications are studied, highlighting the importance of the subject
introduced in [1] both for the discretization of PDEs and the study of
elliptic problems in the discrete settings of weighted graphs. In particular,
we numerically prove the relevance of the derived spectral information for
building fast iterative solvers and we present various numerical examples to
confirm the validity of our derivations. Finally, in Section 4 we draw
conclusions and discuss possible directions of research.

##  2 Background, notation and classes of locally uniform graphs

In this section we recall some basics on graph theory (see, for example, the
very recent [36] for a modern exposition) and introduce definitions, notation
and several families of graphs with uniform local structure, which will be of
interest for our numerical applications. Before doing so, we briefly introduce
a multi-index notation that will be used hereafter.

For a fixed integer $d\geq 1$, a $d$-index ${\boldsymbol{k}}$ is a vector
${\boldsymbol{k}}=(k_{1},...,k_{d})\in\mathbb{Z}^{d}$. Given two $d$-indices
${\boldsymbol{i}}=(i_{1},...,i_{d})$ and ${\boldsymbol{j}}=(j_{1},...,j_{d})$,
we write ${\boldsymbol{i}}\vartriangleleft{\boldsymbol{j}}$ if $i_{r}<j_{r}$
for the first $r=1,2,\ldots,d$ such that $i_{r}\neq j_{r}$ and we write
${\boldsymbol{i}}<{\boldsymbol{j}}$ if $i_{r}<j_{r}$ is satisfied for every
$r=1,\ldots,d$. The relations
$\trianglelefteq,\vartriangleright,\trianglerighteq$ and $\leq,>,\geq$ are
defined accordingly.

We use bold letters for vectors and vector/matrix-valued functions. The
notation $\bf{0},\bf{1},\bf{2},\ldots$ is used for the $d$-dimensional
constant vectors $\left(0,0,\ldots,0\right)$, $\left(1,1,\ldots,1\right)$,
$\left(2,2,\ldots,2\right),\ldots$, respectively. We intend every operation on
vectors in $\mathbb{Z}^{d}$ computed element-wise, so that, for example, the
notation $\frac{{\boldsymbol{i}}}{{\boldsymbol{n}}}$ is used for the vector
$\frac{{\boldsymbol{i}}}{{\boldsymbol{n}}}=\left(\frac{i_{1}}{n_{1}},\ldots,\frac{i_{d}}{n_{d}}\right)$
and $|{\boldsymbol{i}}|$ is used for the vector
$|{\boldsymbol{i}}|=\left(|i_{1}|,\ldots,|i_{d}|\right)$. Finally, given a
$d$-index ${\boldsymbol{n}}$, we write ${\boldsymbol{n}}\to\infty$ if and only
if $\min_{r=1,\ldots,d}\\{n_{r}\\}\to\infty$.

###  2.1 Graphs: basics, notation and $d$-level diamond Toeplitz graphs

We start this subsection by recalling some notation and basic information
about graphs and subgraphs. After this, we introduce the main structures of
graphs which will be of interest for numerical application, namely Toeplitz
graphs, $d$-level Toeplitz graphs and $d$-level diamond Toeplitz graphs (see
[1]).

A (possily infinite) graph is a quadruple $G=(V,E,w,\kappa)$, where

  * •

$V=\left\\{v_{1},v_{2},\ldots,v_{n},...\right\\}$ denotes a countable set of
nodes, or vertices;

  * •

$E=\\{(v_{i},v_{j})|\,v_{i},v_{j}\in V,\;\exists\text{ an edge from
}v_{i}\text{ to }v_{j}\\}$ is the set of edges between nodes;

  * •

$w:V\times V\to\mathbb{R}$ is the weight function and, given two nodes $v_{i}$
and $v_{j}$, we have $w(v_{i},v_{j})\neq 0$ if and only if there exists an
edge from $v_{i}$ to $v_{j}$;

  * •

a potential term $\kappa:V\to\mathbb{R}$.

We use the notation $G=(V,E,w)$ whenever the potential term $\kappa$ vanishes.
The non-zero values of the weight function $w$ are called weights associated
with the edge $(v_{i},v_{j})$. A walk of length $k$ in $G$ is a set of nodes
$v_{i_{0}},v_{i_{1}},\ldots,v_{i_{k-1}},v_{i_{k}}$ such that, for all $1\leq
r\leq k$, $(v_{i_{r-1}},v_{i_{r}})\in E$. In this paper we work with
undirected and connected graphs, i.e. graphs for which the function $w$ is
symmetric ($w(v_{i},v_{j})=w(v_{j},v_{i})$ for every $v_{i},v_{j}\in V$) and
such that for every couple of vertices $v_{i},v_{j}\in V$ there exists a walk
from $v_{i}$ to $v_{j}$. If $w(v_{i},v_{j})\neq 0$ we write $v_{i}\sim v_{j}$
and we say that $v_{i}$ and $v_{j}$ are neighbors. Finally, the degree of a
node $v_{i}$ of an undirected graph, denoted by $\deg(v_{i})$, is defined as

$\deg(v_{i}):=\sum_{v_{j}\sim v_{i}}w(v_{i},v_{j}).$

We are now ready to introduce the main structures of graphs we are interested
in for our numerical applications. We follow the same notation and the same
increasing level of complexity proposed in [1], hence starting with the
definition of Toeplitz graphs.

######  Definition 2.1 (Toeplitz graph).

Let $n,m,t_{1},\ldots,t_{m}$ be positive integers such that
$0<t_{1}<t_{2}<\ldots<t_{m}<n$, and fix $m$ nonzero real numbers
$w_{t_{1}},\ldots,w_{t_{m}}$. A Toeplitz graph of cardinality $n$ is the graph
$T_{n}\langle(t_{1},w_{t_{1}}),\ldots,(t_{m},w_{t_{m}})\rangle=(V_{n},E_{n},w)$
such that $V_{n}=\\{v_{1},\ldots,v_{n}\\}$ and the weight function $w$
satisfies

$w(v_{i},v_{j})=\begin{cases}w_{t_{k}}&\mbox{if }|i-j|=t_{k},\\\
0&\mbox{otherwise}.\end{cases}$

Note that, by construction, the adjacency matrix of a Toeplitz graph, that is
the matrix $W_{n}=(w(v_{i},v_{j}))_{i,j=1}^{n}=(w_{|i-j|})_{i,j=1}^{n}$, is a
symmetric Toeplitz matrix. Trivially, there exists a 1-1 correspondence
between Toeplitz graphs and (symmetric) Toeplitz matrices.

A generalization of Toeplitz graphs is represented by $d$-level Toeplitz
graphs. In order to introduce them, we first need to recall the concept of
directions associated with a $d$-index.

Given a $d$-index ${\boldsymbol{t}}_{k}=((t_{k})_{1},\ldots,(t_{k})_{d})$ such
that $\mathbf{0}\trianglelefteq{\boldsymbol{t}}_{k}$ and
${\boldsymbol{t}}_{k}\neq\mathbf{0}$, the set

$[{\boldsymbol{t}}_{k}]:={{}^{\textstyle
I_{k}}\big{/}_{\textstyle\sim}},\qquad\mbox{where
}{\boldsymbol{i}}\sim{\boldsymbol{j}}\mbox{ iff
}{\boldsymbol{i}}=\pm{\boldsymbol{j}},$

where

$I_{k}:=\left\\{{\boldsymbol{i}}\in\mathbb{Z}^{d}\,|\,{\boldsymbol{i}}=\left(\pm(t_{k})_{1},\ldots,\pm(t_{k})_{d}\right)\right\\},$

is called the set of directions associated with ${\boldsymbol{t}}_{k}$. For
$\alpha=1,\ldots,\left|[{\boldsymbol{t}}_{k}]\right|$, the elements
$[{\boldsymbol{t}}_{k}]_{\alpha}\in[{\boldsymbol{t}}_{k}]$ are called
directions and clearly $\left|[{\boldsymbol{t}}_{k}]_{\alpha}\right|=2$. The
element in $[{\boldsymbol{t}}_{k}]_{\alpha}$ whose first nonzero component is
positive is denoted by $[{\boldsymbol{t}}_{k}]_{\alpha}^{+}$ and the other one
is denoted by $[{\boldsymbol{t}}_{k}]_{\alpha}^{-}$.

######  Definition 2.2 ($d$-level Toeplitz graph).

Let ${\boldsymbol{n}},{\boldsymbol{t}}_{1},\ldots,{\boldsymbol{t}}_{m}$ be
$d$-indices such that $\bf{0}<{\boldsymbol{n}}$,
$0\vartriangleleft{\boldsymbol{t}}_{1}\triangleleft{\boldsymbol{t}}_{2}\triangleleft\ldots\triangleleft{\boldsymbol{t}}_{m}\triangleleft{\boldsymbol{n}}-\mathbf{1},$
and let ${\boldsymbol{w}}_{1},\ldots,{\boldsymbol{w}}_{m}$ be $m$ nonzero real
vectors such that ${\boldsymbol{w}}_{k}\in\mathbb{R}^{c_{k}}$, with
$c_{k}=\left|[{\boldsymbol{t}}_{k}]\right|$, for every $k=1,\ldots,m$. The
components of the vectors ${\boldsymbol{w}}_{k}$ are denoted as follows:

${\boldsymbol{w}}_{k}=\left(w_{[{\boldsymbol{t}}_{k}]_{1}},w_{[{\boldsymbol{t}}_{k}]_{2}},\ldots,w_{[{\boldsymbol{t}}_{k}]_{c_{k}}}\right).$

A $d$-level Toeplitz graph
$T_{\boldsymbol{n}}\langle\\{[{\boldsymbol{t}}_{1}],{\boldsymbol{w}}_{1}\\},\ldots,\\{[{\boldsymbol{t}}_{m}],{\boldsymbol{w}}_{m}\\}\rangle$
is an undirected graph with node set given by
$V_{\boldsymbol{n}}=\left\\{v_{{\boldsymbol{k}}}\,|\,\bf{1}\trianglelefteq{\boldsymbol{k}}\trianglelefteq{\boldsymbol{n}}\right\\}$
and whose weight function $\omega$ satisfies

$w(v_{\boldsymbol{i}},v_{\boldsymbol{j}})=\begin{cases}w_{[{\boldsymbol{t}}_{k}]_{\alpha}}&\mbox{if
}|{\boldsymbol{i}}-{\boldsymbol{j}}|={\boldsymbol{t}}_{k}\mbox{ and
}({\boldsymbol{i}}-{\boldsymbol{j}})\in[{\boldsymbol{t}}_{k}]_{\alpha}=\left\\{[{\boldsymbol{t}}_{k}]_{\alpha}^{+},[{\boldsymbol{t}}_{k}]_{\alpha}^{-}\right\\}\\\
&\mbox{for some }\alpha=1,\ldots,c_{k};\\\ 0&\mbox{otherwise}.\end{cases}$ (1)

When every node of a $d$-level Toeplitz graph is replaced by a fixed,
undirected graph of dimension $\nu$, which we call diamond, the graph assumes
a block structure. To properly define this concept, recall that a _linking
graph operator_ for the reference node set $[\nu]:=\\{1,\ldots,\nu\\}$ is any
non-zero $\mathbb{R}^{\nu\times\nu}$ matrix, or, equivalently, the adjacency
matrix of a (possibly not undirected) graph. We can then give the following
definition.

######  Definition 2.3 ($d$-level diamond Toeplitz graph).

Let $d,m,\nu$ be fixed integers and let $G\simeq\left([\nu],E,w\right)$ be a
fixed undirected graph which we call mold graph.

Let ${\boldsymbol{n}},{\boldsymbol{t}}_{1},\ldots,{\boldsymbol{t}}_{m}$ be
$d$-indices such that $\bf{0}<{\boldsymbol{n}}$, and
$0\vartriangleleft{\boldsymbol{t}}_{1}\triangleleft{\boldsymbol{t}}_{2}\triangleleft\ldots\triangleleft{\boldsymbol{t}}_{m}\triangleleft{\boldsymbol{n}}-\mathbf{1}.$
For $k=1,\ldots,m$, let $\boldsymbol{L}_{k}$ be a collection of linking graph
operators for the reference set $[\nu]$ such that
$\left|\boldsymbol{L}_{k}\right|=c_{k}$, with
$c_{k}=\left|[{\boldsymbol{t}}_{k}]\right|$ for every $k=1,\ldots,m$. We
indicate the elements of the set $\boldsymbol{L}_{k}$ by the following index
notation,

$\displaystyle\boldsymbol{L}_{k}=\left\\{L_{[{\boldsymbol{t}}_{k}]_{1}},L_{[{\boldsymbol{t}}_{k}]_{2}},\ldots,L_{[{\boldsymbol{t}}_{k}]_{c_{k}}}\right\\},$
$\displaystyle\mathbb{R}^{\nu\times\nu}\ni
L_{[{\boldsymbol{t}}_{k}]_{\alpha}}=\left(l_{[{\boldsymbol{t}}_{k}]_{\alpha}}(r,s)\right)_{r,s=1}^{\nu}\quad\mbox{for
}\alpha=1,\ldots,c_{k}.$

Finally, consider ${\boldsymbol{n}}$ copies $G({\boldsymbol{k}})\simeq G$ of
the mold graph, which we call diamonds.

A $d$-level diamond Toeplitz graph, denoted by

$T_{{\boldsymbol{n}},\nu}^{G}\left\langle\left\\{{\boldsymbol{t}}_{1},\boldsymbol{L}_{1}\right\\},\ldots,\left\\{{\boldsymbol{t}}_{m},\boldsymbol{L}_{m}\right\\}\right\rangle,$

is an undirected graph with nodes

$V_{\boldsymbol{n}}=\left\\{v_{({\boldsymbol{k}},r)}\,|\,({\bf{1}},1)\trianglelefteq({{\boldsymbol{k}}},r)\trianglelefteq({\boldsymbol{n}},\nu)\right\\}$

and weight function $w_{{\boldsymbol{n}}}:V_{\boldsymbol{n}}\times
V_{\boldsymbol{n}}\to\mathbb{R}$ satisfying

$w_{{\boldsymbol{n}}}\left(v_{({\boldsymbol{i}},r)},v_{({\boldsymbol{j}},s)}\right):=\begin{cases}w(r,s)&\mbox{if
}{\boldsymbol{i}}={\boldsymbol{j}},\\\
l_{[{\boldsymbol{t}}_{k}]_{\alpha}}(r,s)&\mbox{if
}|{\boldsymbol{i}}-{\boldsymbol{j}}|={\boldsymbol{t}}_{k}\mbox{ and
}({\boldsymbol{i}}-{\boldsymbol{j}})=[{\boldsymbol{t}}_{k}]_{\alpha}^{+},\\\
l_{[{\boldsymbol{t}}_{k}]_{\alpha}}(s,r)&\mbox{if
}|{\boldsymbol{i}}-{\boldsymbol{j}}|={\boldsymbol{t}}_{k}\mbox{ and
}({\boldsymbol{i}}-{\boldsymbol{j}})=[{\boldsymbol{t}}_{k}]_{\alpha}^{-},\\\
0&\mbox{otherwise}.\end{cases}$

To end this section, we now introduce the main character of our numerical
application, namely the graph Laplacian. First, observe that we use the
notation $C(V)$ to denote the set of real-valued functions on $V$.

######  Definition 2.4 (graph Laplacian).

Let $G=(V,E,w,\kappa)$ be an undirected graph with no self-loops (i.e.
$w(v_{i},v_{i})=0$ for every $v_{i}\in V$). The graph Laplacian is the
symmetric matrix $\Delta_{G}:C(V)\to C(V)$ defined as

$\Delta_{G}:=D+K-W,$

where $D$ is the degree matrix, $K$ is the potential term matrix, that is,

$D:=\mathop{\mathrm{diag}}\left\\{\deg(v_{1}),\ldots,\deg(v_{n})\right\\},\quad
K:=\mathop{\mathrm{diag}}\left\\{\kappa(v_{1}),\ldots,\kappa(v_{n})\right\\},$

and $W$ is the adjacency matrix of the graph $G$, that is,

$W=\left(\begin{array}[]{cccc}0&w(v_{1},v_{2})&\cdots&w(v_{1},v_{n})\\\
w(v_{1},v_{2})&0&\ddots&\vdots\\\ \vdots&\ddots&\ddots&w(v_{n-1},v_{n})\\\
w(v_{1},v_{n})&\cdots&w(v_{n-1},v_{n})&0\\\ \end{array}\right).$

Therefore, $\Delta_{G}$ appears as follows:

$\Delta_{G}=\left(\begin{array}[]{cccc}\deg(v_{1})+\kappa(v_{1})&-w(v_{1},v_{2})&\cdots&-w(v_{1},v_{n})\\\
-w(v_{1},v_{2})&\deg(v_{2})+\kappa(v_{2})&\ddots&\vdots\\\
\vdots&\ddots&\ddots&-w(v_{n-1},v_{n})\\\
-w(v_{1},v_{n})&\cdots&-w(v_{n-1},v_{n})&\deg(v_{n})+\kappa(v_{n})\\\
\end{array}\right).$

###  2.2 Matrix-norms and partial ordering of Hermitian matrices

In this brief subsection, we introduce some notation for matrix-norms which
will be employed in the following sections and we define the partial ordering
on the subset of $\mathbb{C}^{n\times n}$ formed by Hermitian matrices.
First, given a norm $\|\cdot\|$ on $\mathbb{C}^{n}$, it is always possible to
consider the induced norm, which will be still denoted by $\|\cdot\|$, on the
space $\mathbb{C}^{n\times n}$ defined as follows:

$\|A\|=\sup\left\\{\|Ax\|,\,x\in\mathbb{C}^{n}\text{ with
}\|x\|=1\right\\}=\sup_{x\neq\bf{0}}\frac{\|Ax\|}{\|x\|}.$

In our work we will be mainly interested in the case where the considered norm
on $\mathbb{C}^{n}$ is the p-norm, for some $p\in[1,\infty]$, which we denote
by $\|\cdot\|_{p}$, if $p\neq 2$, and by $\|\cdot\|$, if $p=2$. In this last
particular case the norm is also known as spectral norm and the following
equality holds

$\|A\|=\sigma_{\max}(A),$

where $\sigma_{\max}(A)$ is the largest singular value of $A$.

We denote the Frobenius norm of a matrix $A=(a_{i,j})_{i,j=1}^{n}$ by

$\|A\|_{F}=\sqrt{\sum_{i=1}^{n}\sum_{j=1}^{n}|a_{i,j}|^{2}}.$

Let now $A$ be a Hermitian Positive Definite (HPD) matrix. This assumption
will be crucial in what follows, especially in Section 3, for the definition
of our iterative methods. In fact, in this special case, it is possible to
define an inner product on $\mathbb{C}^{n}$ as

$\left\langle u,v\right\rangle_{A}=u^{T}Av,$

which induces the following norm on $\mathbb{C}^{n}$

$\|u\|_{A}=\sqrt{\left\langle u,u\right\rangle_{A}}=\sqrt{u^{T}Au}.$

To conclude this subsection, we define the following partial ordering on the
subset of $\mathbb{C}^{n\times n}$ of Hermitian matrices. Given two Hermitian
matrices $A$ and $B$ we say that $A\geq B$ (respectively, $A>B$) if $A-B$ is
Hermitian Positive Semidefinite (HPSD) (respectively, $A-B$ is HPD).

###  2.3 Spectral symbol and generating function of $d$-level diamond Toeplitz
graphs

A matrix-valued function
$\boldsymbol{\mathfrak{f}}:D\to\mathbb{C}^{\nu\times\nu}$, $\nu\geq 1$,
defined on a measurable set $D\subseteq\mathbb{R}^{m}$, $m\in\mathbb{N}$, is
measurable (resp. continuous, in $\textnormal{L}^{p}(D)$) if its components
$\mathfrak{f}_{i,j}:D\to\mathbb{C},\ i,j=1,\ldots,\nu,$ are measurable (resp.
continuous, in $\textnormal{L}^{p}(D)$). Let $\mu_{m}$ be the Lebesgue measure
on $\mathbb{R}^{m}$ and let $C_{c}(\mathbb{R})$ be the set of continuous
functions with bounded support defined over $\mathbb{R}$. We use the notation
$\\{X_{{\boldsymbol{n}},\nu}\\}$ to denote a sequence of matrices of
increasing dimension $d_{{\boldsymbol{n}}}$, i.e. such that
$d_{{\boldsymbol{n}}}\to\infty$ as ${\boldsymbol{n}}\to\infty$, $\nu$ being a
fixed parameter independent of ${\boldsymbol{n}}$.

We say that $\\{X_{{\boldsymbol{n}},\nu}\\}_{\boldsymbol{n}}$ is zero-
distributed if

$\lim_{{\boldsymbol{n}}\rightarrow\infty}\frac{1}{d_{\boldsymbol{n}}}\sum_{k=1}^{d_{\boldsymbol{n}}}F(\sigma_{k}(X_{{\boldsymbol{n}},\nu}))=F(0)\qquad\forall
F\in C_{c}(\mathbb{R}),$ (2)

and we indicate it by
$\\{X_{{\boldsymbol{n}},\nu}\\}_{\boldsymbol{n}}\sim_{\sigma}0$, where
$\sigma_{1}(X_{{\boldsymbol{n}},\nu}),\ldots,\sigma_{d_{{\boldsymbol{n}}}}(X_{{\boldsymbol{n}},\nu})$
are the singular value of $X_{{\boldsymbol{n}},\nu}$, sorted in non-decreasing
order.

Of great importance for our numerical applications is the knowledge of the
spectral symbol of a sequence of matrices, which we now recall.

######  Definition 2.5 (Spectral symbol).

Let $\\{X_{{\boldsymbol{n}},\nu}\\}_{\boldsymbol{n}}$ be a sequence of
matrices and let $\boldsymbol{\mathfrak{f}}:D\to\mathbb{C}^{\nu\times\nu}$ be
a measurable Hermitian matrix-valued function defined on the measurable set
$D\subset\mathbb{R}^{m}$, with $0<\mu_{m}(D)<\infty$.

We say that $\\{X_{{\boldsymbol{n}},\nu}\\}_{\boldsymbol{n}}$ is distributed
like $\boldsymbol{\mathfrak{f}}$ in the sense of eigenvalues, in symbols
$\\{X_{{\boldsymbol{n}},\nu}\\}_{\boldsymbol{n}}\sim_{\lambda}\boldsymbol{\mathfrak{f}}$,
if

$\lim_{{\boldsymbol{n}}\rightarrow\infty}\frac{1}{d_{\boldsymbol{n}}}\sum_{k=1}^{d_{\boldsymbol{n}}}F(\lambda_{k}(X_{{\boldsymbol{n}},\nu}))=\frac{1}{\mu_{m}(D)}\int_{D}\sum_{k=1}^{\nu}F\left(\lambda_{k}(\boldsymbol{\mathfrak{f}}(\boldsymbol{y}))\right)d\mu_{m}(\boldsymbol{y}),\qquad\forall
F\in C_{c}(\mathbb{R}),$ (3)

where
$\lambda_{1}(\boldsymbol{\mathfrak{f}}(\boldsymbol{y})),\ldots,\lambda_{\nu}(\boldsymbol{\mathfrak{f}}(\boldsymbol{y}))$
are the eigenvalues of $\boldsymbol{\mathfrak{f}}(\boldsymbol{y})$ and
$\lambda_{1}(X_{{\boldsymbol{n}},\nu}),\ldots,\lambda_{d_{{\boldsymbol{n}}}}(X_{{\boldsymbol{n}},\nu})$
are the eigenvalues of $\\{X_{{\boldsymbol{n}},\nu}\\}$, sorted in non-
decreasing order.

It is easy to show that (see [1, Proposition 4.4]), given a $d$-level diamond
Toeplitz graph
$T_{{\boldsymbol{n}},\nu}^{G}\left\langle\left\\{{\boldsymbol{t}}_{1},\boldsymbol{L}_{1}\right\\},\ldots,\left\\{{\boldsymbol{t}}_{m},\boldsymbol{L}_{m}\right\\}\right\rangle$,
there exists a function
${\boldsymbol{f}}:[-\pi,\pi]^{d}\to\mathbb{C}^{\nu\times\nu}$, called the
generating function, such that the adjacency matrix
$W^{G}_{{\boldsymbol{n}},\nu}$ of the graph satisfies

$(W^{G}_{{\boldsymbol{n}},\nu})_{{\boldsymbol{i}},{\boldsymbol{j}}}=\hat{{\boldsymbol{f}}}_{{\boldsymbol{i}}-{\boldsymbol{j}}}$

and, therefore (see [46]),
$\\{W^{G}_{{\boldsymbol{n}},\nu}\\}\sim_{\lambda}\boldsymbol{\mathfrak{f}}\equiv{\boldsymbol{f}}$,
whenever ${\boldsymbol{f}}\in\text{L}^{1}([-\pi,\pi]^{d})$.

###  2.4 Multigrid Methods for Structured Matrices

In the current subsection we summarize the key results on symbol-based
multigrid methods for Toeplitz structures. Throughout the subsection, we
consider the case of an invertible matrix $A_{n}\in\mathbb{C}^{d_{n}\times
d_{n}}$ and a vector $b\in\mathbb{C}^{d_{n}}$.

Multigrid methods (MGM) are iterative procedures for large linear systems that
create a proper sequence of linear systems of decreasing dimensions obtained
by consecutive projections. First, we focus on a general coefficient matrix
$A_{n}\in\mathbb{C}^{d_{n}\times d_{n}}$, then we treat in detail the results
related to Toeplitz matrices $T_{n}(f)$.

Given $x_{n},\,b_{n}\in\mathbb{C}^{d_{n}}$, consider a full-rank matrix
$P_{n,k}\in\mathbb{C}^{d_{n}\times k}$, $k<d_{n}$, and let us consider two
stationary iterative methods: the method $\mathcal{V}_{n,\rm{pre}}$, with
iteration matrix ${V}_{n,\rm{pre}}$, and $\mathcal{V}_{n,\rm{post}}$, with
iteration matrix ${V}_{n,\rm{post}}$.

Given an initial guess $x_{n}^{(0)}\in\mathbb{C}^{d_{n}}$, an iteration of a
two-grid method (TGM) is given by the following steps:

$x_{n}^{(k+1)}=\mathcal{TGM}(A_{n},x_{n}^{(j)},b_{n})$
---
| 0\.
$x_{n}^{\rm{pre}}=\mathcal{V}_{n,\rm{pre}}^{\nu_{\rm{pre}}}(A_{n},{b}_{n},x_{n}^{(j)})$
---
| Pre-smoothing iterations
| 1\. $r_{n}=b_{n}-A_{n}x_{n}^{\rm{pre}}$
---
2\. $r_{k}=P_{n,k}^{H}r_{n}$
3\. $A_{k}=P_{n,k}^{H}A_{n}P_{n,k}$
4\. Solve $A_{k}y_{k}=r_{k}$
5\. $\hat{x}_{n}=x_{n}^{\rm{pre}}+P_{n,k}y_{k}$
| Coarse Grid Correction (CGC)
| 6\.
$x_{n}^{(j+1)}=\mathcal{V}_{n,\rm{post}}^{\nu_{\rm{post}}}(A_{n},{b}_{n},\hat{x}_{n})$
---
| Post-smoothing iterations

The Coarse Grid Correction (CGC) depends on the grid transfer operator
$P_{n,k}$, while step $0.$ and step $6.$ consist, respectively, in applying
$\nu_{\rm{pre}}$ times a pre-smoother and $\nu_{\rm{post}}$ times a post-
smoother of the given iterative methods.

We base the convergence analysis on the Ruge-Stüben theory [38] for TGM. In
particular, we mention the results in [2, Remark 2.2] and [38, Theorem 5.2],
which provide two separate conditions that must be fulfilled in order to have
TGM convergence, namely the “smoothing property” and the “approximation
property”. We highlight that, if such conditions are satisfied, then the
resulting TGM is not only convergent but has also an optimal convergence rate.

The standard V-cycle method is obtained by replacing the direct solution at
step 4. with a recursive call of the TGM applied to the coarser linear system
$A_{k_{\ell}}y_{k_{\ell}}=r_{k_{\ell}}$, where $\ell$ represents the level.
The recursion is usually stopped at level ${\overline{\ell}}$ when
$k_{\overline{\ell}}$ becomes small enough for solving cheaply step 4. with a
direct solver.

In the case where $A_{n}=T_{n}(f)$, an efficient choice of grid transfer
operator is given by $P_{n,k}$ obtained as the product between a Toeplitz
matrix $T_{n}(p)$, with $p$ trigonometric polynomial, and the cutting matrix
$K_{n}$. That is,

$P_{n,k}=T_{n}(p)K_{n},$ (4)

where, if the size $n$ of the coefficient matrix is divisible by a factor
$\mathfrak{g}\geq 2$, the lower level dimension is reduced to
$k=n/\mathfrak{g}$ by the cutting matrix

$K_{n}=[\delta_{i-\mathfrak{g}j}]_{i,j},\quad
i=0,\dots,n-1;\,j=0,\dots,k-1,\quad\quad\delta_{\ell}=\begin{cases}1\,&{\rm
if}\,\ell\equiv 0\,({\rm mod}\,n),\\\ 0\,&{\rm otherwise}\end{cases}.$ (5)

The validation of the approximation property depends on the choice of the
generating function $p$, which has to fulfill specific conditions depending on
the zeros of $f$.

In the case where $f$ is a $d$-variate function, i.e. $T_{\boldsymbol{n}}(f)$
is a $d$-level matrix, we construct $K_{{\boldsymbol{n}}}$ as the tensor
product of the cutting matrices $K_{n_{1}}$, $K_{n_{2}}$, $\dots$,
$K_{n_{d}}$.

In what follows we provide, in the $d$-level setting, conditions on $p$ in
order to obtain a projector which is effective in terms of convergence. For
the generalization to the case where ${\boldsymbol{f}}$ is a $\nu\times\nu$
matrix-valued generating function see [21].

######  Theorem 2.1 ([2, 23, 26]).

Consider a Toeplitz matrix $T_{\boldsymbol{n}}(f)$ generated by a non-negative
trigonometric polynomial $f$. Suppose that $f({\boldsymbol{\theta}})$ vanishes
at exactly one point $\mathbf{{\boldsymbol{\theta}}_{0}}$. Then, the
optimality of the two-grid method applied to $T_{\boldsymbol{n}}(f)$ is
guaranteed if we choose a projection operator of the form (4) associated with
a generating function $p$ such that

$\begin{split}\underset{{{\boldsymbol{\theta}}}\to{{\boldsymbol{\theta}}_{0}}}{\lim\sup}\frac{|p({{\boldsymbol{\eta}}})|^{2}}{f({{\boldsymbol{\theta}}})}&<\infty,\quad{{\boldsymbol{\eta}}}\in\mathcal{M}_{\mathfrak{g}}({{\boldsymbol{\theta}}}),\\\
\sum_{{{\boldsymbol{\eta}}}\in\Omega({{\boldsymbol{\theta}}})}p^{2}({{\boldsymbol{\eta}}})&>0,\end{split}$

where

$\Omega_{\mathfrak{g}}({{\boldsymbol{\theta}}})=\left\\{{\boldsymbol{\eta}}\left|\eta_{j}\in\left\\{\theta_{j}+\frac{2\pi
k}{\mathfrak{g}}(\mbox{\rm mod
}2\pi)\right\\},k=0,\dots,\mathfrak{g}-1,j=1,\dots,d\right.\right\\}$

and
$\mathcal{M}_{\mathfrak{g}}({{\boldsymbol{\theta}}})=\Omega_{\mathfrak{g}}({{\boldsymbol{\theta}}})\setminus\\{{\boldsymbol{\theta}}\\}$
are the sets of corner and mirror points respectively.

For achieving the optimality of the V-cycle method, the first condition needs
to be strengthened, see [2] for details. Moreover, it was proven in [23] that
the coarsening factor $\mathfrak{g}$ needs to be strictly greater than 2 for
an overall optimal method, otherwise each multigrid iteration would require a
computational cost greater than the cost of the matrix vector product.
Finally, we highlight that a coarsening factor greater than 2 permits to
construct less coarser levels before solving the error equation exactly. This
implies that the critical issues (i.e. low rank corrections, increase of the
condition number), that could worsen level after level, are propagated a
smaller number of times.

##  3 Numerical examples

In this section we provide numerical examples which emphasize the
computational help of the theory developed in [1] for studying many different
kinds of discretization of PDEs as well as for studying elliptic problems in
the very general setting of weighted graphs. Using the spectral features of
the considered graphs and in particular the analytical properties of the
underlying symbol, we can define appropriate preconditioners for the
preconditioned conjugate gradient (PCG) (see [39] for some general theory on
preconditioned Krylov methods) and appropriate projectors for the multigrid
method (MGM) (see [37]), with the aim of making such methods optimal. For the
optimality of MGM for unilevel Toeplitz structures corresponding to
approximated differential equations in one dimension see [26]; for multi-level
Toeplitz and non-Toeplitz structures corresponding to approximated PDEs in
dimension greater than one see [2, 40]. Furthermore, we show that without the
proposed preconditioning suggested in the present work, the convergence to the
solution of the problem can be very slow, due to the bad conditioning of the
discrete problems under consideration.

First, we show some numerical experiments which exploit the complete power of
the graph structures recalled in Section 2 and the related spectral properties
developed in [1]. At this purpose, we first perform tests in which graphs are
used as a discretization of the geometry of a sub-domain of $[0,1]^{d}$,
$d\geq 1$. Then, we move into a complete graph setting, where the considered
graph structure can represent a model for various type of networks and real-
world interactions. For the sake of completeness with respect to the first
part of the present work, we also perform numerical tests concerning examples
in [1, Subsections 7.2 and 7.3], which have been already investigated in the
relevant literature (see, for example, [15, 17, 18, 19, 20]).

Before proceeding, we recall some notation concerning (proper) subgraphs and
the concept of immersion of a $d$-level Toeplitz graph in a subdomain of
$[0,1]^{d}$.

Given a graph $\bar{G}=(\bar{V},\bar{E},\bar{w},\bar{\kappa})$ and a (proper)
subset $V\subset\bar{V}$, we say that $G=(V,E,w,\kappa)$ is a (proper)
subgraph of $\bar{G}$, and we write $G\subset\bar{G}$ if

  * •

$V\subset\bar{V}$;

  * •

$E=\\{(v_{i},v_{j})\in\bar{E}\,|\,v_{i},v_{j}\in V\\}\subset\bar{E}$;

  * •

$w=\bar{w}_{|E}$;

  * •

$\kappa=\bar{\kappa}_{|\mathring{V}}$,

where the set of nodes

$\mathring{V}:=\left\\{v_{i}\in
V\,|\,v_{i}\nsim\bar{v}_{j}\,\forall\,\bar{v}_{j}\in\bar{V}\setminus
V\right\\}$

is called interior of $V$ and its element are called interior nodes.
Finally, given a subgraph $G=(V,E,w,\kappa)$ of
$\bar{G}=(\bar{V},\bar{E},\bar{w},\bar{\kappa})$, we can define the boundary
of $V$ as

$\partial V:=\left\\{v_{i}\in V\,|\,v_{i}\sim\bar{v}_{j}\,\mbox{for some
}\,\bar{v}_{j}\in\bar{V}\setminus V\right\\},$

whose elements are called boundary nodes. Note that we do not request that
$\kappa=\bar{\kappa}$ on $\partial V$.

For the concept of immersion, consider a $d$-level Toeplitz graph
$T_{\boldsymbol{n}}\langle\\{[{\boldsymbol{t}}_{1}],{\boldsymbol{w}}_{1}\\},\ldots,\\{[{\boldsymbol{t}}_{m}],{\boldsymbol{w}}_{m}\\}\rangle$
with weight function $w$ and a continuous almost everywhere function
$p:[0,1]^{d}\to\mathbb{R}$. We then consider the immersion map
$\iota:V_{\boldsymbol{n}}\to(0,1)^{d}$ such that

$\iota(v_{\boldsymbol{j}}):={\boldsymbol{j}}\circ{\boldsymbol{h}}=\left(j_{1}h_{1},\ldots,j_{d}h_{d}\right),$

where $\circ$ denotes the Hadamard (component-wise) product and
${\boldsymbol{h}}$ is the $d$-dimensional vector

${\boldsymbol{h}}:=(h_{1},\ldots,h_{d})=\left(\frac{1}{n_{1}+1},\ldots,\frac{1}{n_{d}+1}\right).$

The map $\iota$ induces a grid graph
$G=(V^{\prime}_{{\boldsymbol{n}}},E^{\prime}_{{\boldsymbol{n}}},w^{p})$ in
$[0,1]^{d}$, where

$V^{\prime}_{\boldsymbol{n}}:=\left\\{{\boldsymbol{x}}_{\boldsymbol{k}}=\iota(v_{\boldsymbol{k}})\,|\,\bf{1}\trianglelefteq{\boldsymbol{k}}\trianglelefteq{\boldsymbol{n}}\right\\},\qquad
E^{\prime}_{\boldsymbol{n}}:=\left\\{({\boldsymbol{x}}_{\boldsymbol{i}},{\boldsymbol{x}}_{\boldsymbol{j}})\,|\,w^{p}({\boldsymbol{x}}_{\boldsymbol{i}},{\boldsymbol{x}}_{\boldsymbol{j}})\neq
0\right\\},$
$w^{p}({\boldsymbol{x}}_{\boldsymbol{i}},{\boldsymbol{x}}_{\boldsymbol{j}}):=p\left(\frac{{\boldsymbol{x}}_{\boldsymbol{i}}+{\boldsymbol{x}}_{\boldsymbol{j}}}{2}\right)w(v_{\boldsymbol{i}},v_{\boldsymbol{j}}).$

To extend this concept to any subdomain $\Omega$ of $[0,1]^{d}$ we simply
restrict the set of nodes $V^{\prime}_{{\boldsymbol{n}}}$ to

$V_{n^{\prime}}^{\Omega}:=V^{\prime}_{{\boldsymbol{n}}}\cap\Omega$

and consequently modify the weight function $w^{p}$ to

$w^{\Omega,p}:=w^{p}_{|V_{n^{\prime}}^{\Omega}\times
V_{n^{\prime}}^{\Omega}}.$

It is clearly possible to further extend the same idea to the case of
$d$-level diamond Toeplitz graphs. In this case, however, the choice of the
immersion map $\iota$ is not as natural as in the previous case and different
choices of the immersion map $\iota$ would be able to describe different grid
geometries.

###  3.1 Toeplitz graph with Fourier coefficients immersed in the triangle

In this first numerical example our starting model operator is the classic
(semi-positive definite) Laplacian on a triangular domain with Dirichlet or
Neumann boundary conditions (BCs), that is

$\displaystyle\mathcal{L}:\textnormal{dom}(\mathcal{L})\subset
L^{2}\left(D\right)\to L^{2}\left(D\right),$ (6)
$\displaystyle\mathcal{L}u(x,y):=-\left(\partial^{2}_{xx}u(x,y)+\partial^{2}_{yy}u(x,y)\right)\qquad(x,y)\in
D,$ (7)

where $D$ is an equilateral triangle contained in $[0,1]^{2}$. As usual,
$L^{2}\left(D\right)$ is the space of squares integrable functions and
$\textnormal{dom}(\mathcal{L})$ is an appropriate (Sobolev) subspace of
$L^{2}\left(D\right)$ such that the formal equation (7) is well-defined and
the BCs are satisfied. For example, in the case of Dirichlet BCs,
$\textnormal{dom}(\mathcal{L})=\overline{C_{c}^{\infty}(D)}$ is the closure of
the space of smooth and compactly supported functions in $D$ with respect to
the Sobolev norm. For this numerical test we do not apply a standard
discretization scheme.

The idea is to build a full graph-based discretization scheme looking at the
operator itself rather than at the specific numerical approximation method, by
using the Fourier coefficients of $h(\theta)=\theta^{2}$ defined on $[0,\pi]$.
This is due to the fact that, as it was observed in [7, Corollary C.1], when
discretizing the one-dimensional Laplace operator, the sequence of spectral
symbols associated to the uniform FD scheme with $(2\eta+1)$-points converges
uniformly to $h$ as $\eta\to\infty$. This phenomenon is not restrained to the
FD discretization scheme, but it appears in the IgA framework ([24, Theorem 1,
Theorem 2 and Lemma 1]) and in the Sinc collocation method as well ([35]),
once the class of regularity of the approximating functions is progressively
refined. Let us observe that $h$ is exactly the (inverse and normalized)
asymptotic distribution function, in the Weyl sense, of the one-dimensional
Laplace operator on the unit interval. In some sense, all those methods become
“indistinguishable” in the limit of their regularity parameters, from a
spectral point of view, since their sequences of spectral symbols converge to
the same function, that is, the distribution function of the operator they are
approximating, as one should expect. Therefore, our approach is to exploit
directly $h$ as a weight function.

As a matter of facts, $h(\theta)=\theta^{2}$ is the asymptotic distribution
function of the one-dimensional Laplacian while here we are considering the
two-dimensional Laplacian. This is not an issue, since the space of smooth
functions of the form $g(x,y)=\sum_{k=1}^{n}\alpha_{k}g_{1,k}(x)g_{2,k}(y)$ is
dense in $L^{2}(D)$ (for a reference see [12], for example). Therefore, the
two-dimensional Laplacian can be approximated by “splitting” it into two one-
dimensional Laplacians that act separately along the two different axes.

The full graph-based discretization we implement here is then the following:
first, fix ${\boldsymbol{n}}=(n,n)$ and consider the infinite graph
$\bar{G}_{\boldsymbol{n}}$ given by

$\bar{G}_{{\boldsymbol{n}}}=\left(\bar{T}_{{\boldsymbol{n}}}\left\langle\cup_{k=1}^{\infty}\left\\{[k,0],w_{k}\right\\},\cup_{k=1}^{\infty}\left\\{[0,k],w_{k}\right\\}\right\rangle\right),$

where
$\bar{T}_{{\boldsymbol{n}}}\left\langle\cup_{k=1}^{\infty}\left\\{[k,0],w_{k}\right\\},\cup_{k=1}^{\infty}\left\\{[0,k],w_{k}\right\\}\right\rangle$
is a $2$-level infinite Toeplitz graph characterized by the node set

$\bar{V}_{{\boldsymbol{n}}}=\left\\{(x_{i},y_{j})=\left(\frac{i}{n+1},\frac{j}{n+1}\right)\,:\,i,j\in\mathbb{Z}\right\\},$

and the weight function

$w((x_{i},y_{j}),(x_{r},y_{s})):=w_{k}=(-1)^{k}\frac{2}{k^{2}}\quad\mbox{if
}(|i-r|,|j-s|)\in\\{(k,0),(0,k)\\}.$

We highlight that $(-1)^{k+1}\frac{2}{k^{2}}$ is the $k$-th term of the
Fourier expansion of $h(\theta)=\theta^{2}$ on $[0,\pi]$. Finally, consider
the subgraph $G_{\boldsymbol{n}}\subset\bar{G}_{{\boldsymbol{n}}}$ such that

$G_{{\boldsymbol{n}}}=\left(T^{D}_{{\boldsymbol{n}}}\left\langle\cup_{k=1}^{n}\left\\{[k,0],w_{k}\right\\},\cup_{k=1}^{n}\left\\{[0,k],w_{k}\right\\}\right\rangle,\kappa\right),$

and where the node set of $T^{D}_{{\boldsymbol{n}}}$ is given by

$V_{\boldsymbol{n}}=\bar{V}_{{\boldsymbol{n}}}\cap D.$

The potential term $\kappa$ determines the BCs. If $\kappa=0$ then we call it
the _Neumann potential_ : this is related to the fact that in
$\mathbb{R}^{d}$, the Neumann BCs that characterizes a second-order elliptic
differential operator, as (6)-(7), disappears when passing to the closure of
its associated quadratic form, see [16, Theorem 7.2.1]. Instead, if

$k((x_{i},y_{j}))=\sum_{(x_{r},y_{s})\in\bar{V}_{\boldsymbol{n}}\setminus
V_{\boldsymbol{n}}}w((x_{i},y_{j}),(x_{r},y_{s})),$

then we call it _Dirichlet potential_. It takes into account the edge
deficiency of each node in $G_{\boldsymbol{n}}$ when seen as a node in the
host graph $\bar{G}_{\boldsymbol{n}}$. The Dirichlet potential term arises
naturally for example when discretizing a differential operator with Dirichlet
BCs by means of FD schemes, see [7, pp. 40–42]. For more details about
Dirichlet potentials we refer to [36, Section 2.2] and [10, Definition 2.1 and
Lemma A.1]. Look at Figure 1 for a visual representation of the graph.

$\kappa$$\kappa$$\kappa$$\kappa$$\kappa$$\kappa$$\kappa$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{2}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{2}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{2}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{2}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{2}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{2}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{1}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color@cmyk@fill{0.91}{0}{0.88}{0.12}w_{1}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{1}}$${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}w_{1}}$${\color[rgb]{0,0.88,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.88,0}\pgfsys@color@cmyk@stroke{0.91}{0}{0.88}{0.12}\pgfsys@color<EMAIL_ADDRESS>
Figure 1: Immersion of a $2$-level grid graph inside the equilateral triangle
$D$. The white nodes are the nodes of $V_{\boldsymbol{n}}$ while the gray
nodes belong to $\bar{V}_{\boldsymbol{n}}\setminus V_{\boldsymbol{n}}$. The
green connections represent the weighted edges whose end-nodes are both
interior nodes of $V_{\boldsymbol{n}}$, while the red connections represent
the weighted edges which have at least one end-node that belongs to
$\bar{V}_{\boldsymbol{n}}\setminus V_{\boldsymbol{n}}$. For the seek of
clarity, in the whole figure on the left we explicitly draw with continuous
lines only the edges of distance 1, while a complete representation of the
edges can be seen in the enlargement on the right.

In the end, we are going to solve the equation

$\Delta_{G_{\boldsymbol{n}}}u_{\boldsymbol{n}}=b_{\boldsymbol{n}},$ (8)

where $\Delta_{G_{\boldsymbol{n}}}$ is the graph Laplacian associated to
$G_{\boldsymbol{n}}$. Let us observe that
$(n+1)^{2}\Delta_{G_{\boldsymbol{n}}}$ approximates $\mathcal{L}$, see Table 1
for a comparison of the spectra of the discrete and continuous operators.

$\left\lfloor\frac{k({\boldsymbol{n}})}{d_{\boldsymbol{n}}}\right\rfloor$ | $n=16$ | $n=32$ | $n=64$
---|---|---|---
0.1 | 0.0350 | 0.0429 | 0.0153
0.5 | 0.1177 | 0.0711 | 0.0200
0.8 | 0.1493 | 0.0821 | 0.0196

$\left\lfloor\frac{k({\boldsymbol{n}})}{d_{\boldsymbol{n}}}\right\rfloor$ | $n=16$ | $n=32$ | $n=64$
---|---|---|---
0.1 | 0.0536 | 0.0277 | 0.0153
0.5 | 0.0529 | 0.0297 | 0.0200
0.8 | 0.0819 | 0.0334 | 0.0196

Table 1: Relative errors between the eigenvalue
$\lambda_{k({\boldsymbol{n}})}^{({\boldsymbol{n}})}$ and the corresponding
evaluation of the continuous Laplacian eigenvalue on the equilateral triangle,
with Neumann (left) and Dirichlet (right) BCs. The index $k({\boldsymbol{n}})$
is chosen such that $k({\boldsymbol{n}})/d_{\boldsymbol{n}}$ is constant for
every fixed ${\boldsymbol{n}}=(n,n)$ with $n=16,32,64$.

######  Remark 3.1.

Let us observe that here we allow negative values for the weight function $w$.
This is not uncommon, since in some applications it is natural to admit graphs
with possible negative values associated to their edges. Nevertheless, if we
restrain $w$ to be nonnegative then we can recover again the same operator
$\Delta_{G_{\boldsymbol{n}}}$ as sum of path graph Laplacians ([25]). For such
a construction we refer to [8, Sections II and III].

We have the following preliminary result.

######  Proposition 3.1.

Given the sequence of graph Laplacians $\\{\Delta_{G_{\boldsymbol{n}}}\\}$
defined above, it holds that

$\\{\Delta_{G_{\boldsymbol{n}}}\\}\sim_{\lambda}\mathfrak{f}(\theta_{1},\theta_{2})=\theta_{1}^{2}+\theta_{2}^{2},\quad(\theta_{1},\theta_{2})\in[0,\pi]^{2}.$

###### Proof.

To ease the notation of the proof, we write $v_{r}=(x_{i},y_{j})$ and sort
them by lexicographic ordering. Let us consider the graph

$F_{{\boldsymbol{n}}}=\left(T^{Q}_{{\boldsymbol{n}}}\left\langle\cup_{k=1}^{n}\left\\{[k,0],w_{k}\right\\},\cup_{k=1}^{n}\left\\{[0,k],w_{k}\right\\}\right\rangle,\kappa\right),\qquad
Q=[0,1]^{2},$

where the node set of $T^{D}_{{\boldsymbol{n}}}$ is given by

$V^{Q}_{\boldsymbol{n}}=\bar{V}_{{\boldsymbol{n}}}\cap Q.$

Clearly, it holds that $G_{\boldsymbol{n}}\subset
F_{\boldsymbol{n}}\subset\bar{G}_{\boldsymbol{n}}$. Consider now the graph
Laplacian associated to $F_{\boldsymbol{n}}$ with the Dirichlet potential. We
have that $\Delta_{F_{\boldsymbol{n}}}\in\mathbb{R}^{n^{2}\times n^{2}}$ and,
by Definition 2.4, each element of the main diagonal is of the form

$\displaystyle\left(\Delta_{F_{\boldsymbol{n}}}\right)_{i,i}$
$\displaystyle=\deg(v_{r})+k(v_{r})$
$\displaystyle=\sum_{\begin{subarray}{c}v_{s}\sim v_{r}\\\ v_{s}\in
V^{Q}_{\boldsymbol{n}}\end{subarray}}w(v_{r},v_{s})+\sum_{\begin{subarray}{c}v_{s}\sim
v_{r}\\\ v_{s}\in\bar{V}_{\boldsymbol{n}}\setminus
V^{Q}_{\boldsymbol{n}}\end{subarray}}w(v_{r},v_{s})$
$\displaystyle=\sum_{\begin{subarray}{c}v_{s}\sim v_{r}\\\
v_{s}\in\bar{V}_{\boldsymbol{n}}\end{subarray}}w(v_{r},v_{s})$
$\displaystyle=2\sum_{k\in\mathbb{Z}}(-1)^{k}\frac{2}{k^{2}}$
$\displaystyle=\frac{2\pi^{2}}{3}.$

Noticing that
$\left\\{\frac{2\pi^{2}}{3}\right\\}\cup\left\\{\left((-1)^{k}\frac{2}{k^{2}},(-1)^{j}\frac{2}{j^{2}}\right)\right\\}_{(k,j)\in\mathbb{Z}\times\mathbb{Z}\setminus(0,0)}$
are the Fourier coefficients of
$\mathfrak{f}(\theta_{1},\theta_{2})=\theta_{1}^{2}+\theta_{2}^{2}$ on
$[0,\pi]^{2}$, and that $\Delta_{F_{\boldsymbol{n}}}$ is a $2$-level Toeplitz
matrix whose diagonals of each diagonal blocks are given exactly by the
Fourier coefficients of $\mathfrak{f}$, that is,
$\left(\Delta_{F_{\boldsymbol{n}}}\right)_{{\boldsymbol{i}},{\boldsymbol{j}}}:=\left(\hat{\mathfrak{f}}_{{\boldsymbol{i}}-{\boldsymbol{j}}}\right)_{{\boldsymbol{i}},{\boldsymbol{j}}=\mathbf{1}}^{{\boldsymbol{n}}}$,
then by standard theory it follows easily that

$\\{\Delta_{F_{\boldsymbol{n}}}\\}\sim_{\lambda}\mathfrak{f}.$

Consider now the sequence of matrices
$\\{I_{{\boldsymbol{n}}}(\chi_{D})\\}_{{\boldsymbol{n}}}$ such that

$I_{{\boldsymbol{n}}}(\chi_{D}):=\mathop{\mathrm{diag}}\left(\chi_{D}\left(\frac{{\boldsymbol{i}}}{{\boldsymbol{n}}+\bf{1}}\right)\right)_{{\boldsymbol{i}}=\bf{1},\ldots,{\boldsymbol{n}}}.$

It is known (see for example [3, Lemma 2.3]) that

$\\{I_{{\boldsymbol{n}}(\chi_{D})}\\}_{{\boldsymbol{n}}}\sim_{\lambda}\chi_{D}.$

Moreover, it is not difficult to see that the sequence
$\\{I_{{\boldsymbol{n}}}\\}_{{\boldsymbol{n}}}$ belongs to the wide family of
GLT sequences with symbol $\chi_{D}$. Using the property of $*$-algebra of the
GLT sequences (see [29, GLT properties, pg. 118-119]) it is immediate to
conclude that the sequence
$\\{I_{{\boldsymbol{n}}}(\chi_{D})\Delta_{F_{{\boldsymbol{n}}}}I_{{\boldsymbol{n}}}(\chi_{D})\\}$
satisfies

$\\{I_{{\boldsymbol{n}}}(\chi_{D})\Delta_{F_{{\boldsymbol{n}}}}I_{{\boldsymbol{n}}}(\chi_{D})\\}\sim_{\lambda}\mathfrak{f}\,\chi_{D}.$

Now, $\Delta_{G_{\boldsymbol{n}}}$ is a principal sub-matrix of
$\Delta_{F_{\boldsymbol{n}}}$ (and of
$I_{{\boldsymbol{n}}}(\chi_{D})\Delta_{F_{{\boldsymbol{n}}}}I_{{\boldsymbol{n}}}(\chi_{D})$)
for every ${\boldsymbol{n}}$; therefore, in light of [3, Lemma 3.4], we can
conclude that

$\\{\Delta_{G_{\boldsymbol{n}}}\\}\sim_{\lambda}\mathfrak{f}_{|D}.$

For the case $\kappa=0$ the only difference is that now the main diagonal is
no longer constant. Nevertheless, the main diagonal differs from that of
$\Delta_{F_{{\boldsymbol{n}}}}$ of a term of small norm and this, therefore,
does not affect the resulting symbol function. Proceeding as above we have our
conclusion.

∎

Let us now investigate the convergence rate of the discrete solution
$u_{\boldsymbol{n}}$ to the continuous solution $u$. In the following we
consider the Dirichlet BCs case, but an analogous study can be done for
Neumann BCs. For this purpose, we consider the function $u$ given by the
formula

$u(x,y)=y(y-\sqrt{3}x)(y+\sqrt{3}x-\sqrt{3})$

which vanishes on the boundary of the triangle $D$ and satisfies the equation
$\mathcal{L}[u](x,y)=2\sqrt{3}$. We compute the evaluations $u(x_{i},y_{j})$
for all $(x_{i},y_{j})\in V_{\boldsymbol{n}}$ and we sort them into the vector
$u^{*}_{\boldsymbol{n}}$ in lexicographic ordering. Then, we solve the linear
system (8) with right-hand side $b_{\boldsymbol{n}}$ equal to the vector of
all $2\sqrt{3}$. In Table 2 and Figure 2 we report the relative error of the
discrete solution $u_{\boldsymbol{n}}$ with respect to evaluations in
$u^{*}_{\boldsymbol{n}}$ and we see that it decreases proportionally to the
partial dimension $n$.

$n$ | $d_{\boldsymbol{n}}$ | $\frac{\|u_{\boldsymbol{n}}-u_{\boldsymbol{n}}^{*}\|}{\|u_{\boldsymbol{n}}^{*}\|}$
---|---|---
6 | 18 | 0.3157
14 | 90 | 0.1131
30 | 400 | 0.0638
62 | 1686 | 0.0287
126 | 6920 | 0.0146
254 | 28028 | 0.0071
Table 2: The relative error of the approximation vector $u_{\boldsymbol{n}}$,
solution of (8), with respect to the vector of evaluations
$u_{\boldsymbol{n}}^{*}$ of the exact solution increasing $n$ in the Dirichlet
BCs case.

figureThe log-log scale plot of the 2-norm relative error increasing $n$.

In the following, we propose a multigrid strategy to solve the linear system
(8) based on the spectral distribution provided by Proposition 3.1, taking for
simplicity Gauss-Seidel both as a pre-smoother and post-smoother.

Assume that the partial dimension is of the form $n=2^{t}$, $t\in\mathbb{N}$.
The graph Laplacian $\Delta_{G_{\boldsymbol{n}}}$ has size
$d_{\boldsymbol{n}}$ such that $d_{\boldsymbol{n}}<n^{2}$. Since the spectral
symbol $\mathfrak{f}(\theta_{1},\theta_{2})=\theta_{1}^{2}+\theta_{2}^{2}$ has
a zero of order 2 in $(0,0)$, according to Theorem 2.1 we choose the
trigonometric polynomial
$q(\theta)=4+6\cos(\theta)+4\cos(2\theta)+2\cos(3\theta)$ to construct the
grid transfer operator

$P_{{\boldsymbol{n}},{\boldsymbol{k}}}=T_{\boldsymbol{n}}\left((q(\theta_{1}))(q(\theta_{2}))\right)K_{\boldsymbol{n}}$

with
${\boldsymbol{k}}=\left(\frac{n}{\mathfrak{g}},\frac{n}{\mathfrak{g}}\right)$,
where $\mathfrak{g}$ represents the coarsening factor. Then, we eliminate from
the matrix $P_{{\boldsymbol{n}},{\boldsymbol{k}}}$ the rows
$\mathbf{i}=(i_{1},i_{2})$ such that $(x_{i_{1}},y_{i_{2}})$ does not belong
to $V_{{\boldsymbol{n}}}$ and eliminate the columns $\mathbf{j}=(j_{1},j_{2})$
such that $(x_{j_{1}},y_{j_{2}})$ does not belong to $V_{{\boldsymbol{k}}}$.

In Table 3 we numerically show the validity of our proposed methods, reporting
the number of iterations needed for achieving the tolerance
$\varepsilon=10^{-6}$ when increasing the matrix size. In the first and second
columns, we show the results relative to the two-grid and V-cycle methods for
$\mathfrak{g}=2$, while in the third and fourth columns we consider a more
aggressive coarsening with $\mathfrak{g}=4$. For $\mathfrak{g}=2$, we see that
the number of iterations needed for convergence remains almost constant, when
increasing the size $d_{\boldsymbol{n}}$. The method with $\mathfrak{g}=4$ has
a worse convergence rate, but it can be convenient when considering greater
matix-sizes from the computational complexity point of view, as explained in
Subsection 2.4.

$t$ | $d_{\boldsymbol{n}}$ | Two-grid ($\mathfrak{g}=2$) | V-cycle ($\mathfrak{g}=2$) | Two-grid ($\mathfrak{g}=4$) | V-cycle ($\mathfrak{g}=4)$
---|---|---|---|---|---
3 | 30 | 9 | 9 | 25 | 25
4 | 116 | 10 | 10 | 27 | 27
5 | 454 | 10 | 11 | 33 | 33
6 | 1796 | 10 | 11 | 36 | 37
7 | 6920 | 11 | 12 | 38 | 40
8 | 28028 | 11 | 12 | 39 | 41
Table 3: Comparison of the number of iterations, varying $t$, for the two-grid
and V-cycle methods applied to the linear system in (8) in the Dirichlet BCs
case.

In the Neumann BCs case, the matrix $\Delta_{G_{\boldsymbol{n}}}$ is singular,
since the the column vector $e$ of all ones of dimension $d_{\boldsymbol{n}}$
is an eigenvector for $\Delta_{G_{\boldsymbol{n}}}$ associated to the null
eigenvalue. Then, we consider the matrix
$A_{\boldsymbol{n}}=\Delta_{G_{\boldsymbol{n}}}+\frac{1}{d_{\boldsymbol{n}}}ee^{T}$
and we solve the system
$A_{\boldsymbol{n}}y_{\boldsymbol{n}}=b_{\boldsymbol{n}}$. We have

$b_{\boldsymbol{n}}=A_{\boldsymbol{n}}y_{\boldsymbol{n}}=(\Delta_{G_{\boldsymbol{n}}}+\frac{1}{d_{\boldsymbol{n}}}ee^{T})y_{\boldsymbol{n}}=\Delta_{G_{\boldsymbol{n}}}y_{\boldsymbol{n}}+\frac{1}{d_{\boldsymbol{n}}}e(e^{T}y_{\boldsymbol{n}})$
(9)

and the latter equality splits $b_{\boldsymbol{n}}$ into the sum of a term in
the range of $\Delta_{G_{\boldsymbol{n}}}$ and a term in the null space of
$\Delta_{G_{\boldsymbol{n}}}$. Since the vector $b_{\boldsymbol{n}}$ belongs
to the range of $\Delta_{G_{\boldsymbol{n}}}$, then $ee^{T}y_{\boldsymbol{n}}$
is the vector of all zeros and $y_{\boldsymbol{n}}$ is the solution also of
system (8).

Note that the matrix $ee^{T}$ is a circulant matrix, hence the addition does
not worsen the matrix-vector product computational cost. Since the matrix
$A_{\boldsymbol{n}}$ is SPD, we can apply the PCG for solving the linear
system $A_{\boldsymbol{n}}y_{\boldsymbol{n}}=b_{\boldsymbol{n}}$. In Table 4,
we compare the efficiency, in terms of iteration count increasing the matrix-
size, of two different preconditioning strategies. We see from the first
iteration column that the number of iterations needed by the CG method with no
preconditioning for reaching the desired tolerance $\epsilon=10^{-6}$
increases significantly. The first preconditioning strategy that we use is
inspired by circulant preconditioning. Indeed, we construct the Strang
circulant precondioner (plus a rank 1 correction) and we remove the rows and
columns that correspond to the nodes that do not belong to the considered
equilateral triangle. This procedure is far from giving an optimal convergence
rate. The last strategy that we try consists in applying an MGM as
preconditioner. In particular, we choose a V-cycle method with fixed tolerance
$10^{-1}$ which uses Gauss-Seidel as pre-smoother and post-smoother and the
standard linear interpolation as grid transfer operator. The results are shown
in the last column of Table 4, where the number of iterations needed for the
MGM-PCG to reach the tolerance $\epsilon=10^{-6}$ is almost constant as
$d_{\boldsymbol{n}}$ increases.

$t$ | $d_{\boldsymbol{n}}$ | CG | Circulant–PCG | MGM–PCG
---|---|---|---|---
3 | 30 | 30 | 21 | 6
4 | 116 | 67 | 30 | 8
5 | 454 | $>100$ | 42 | 9
6 | 1796 | $>100$ | 60 | 9
Table 4: Comparison of the number of iterations, varying $t$, for the PGC
method with no preconditioning, a circulant preconditioner and the MGM
preconditioner applied to the linear system
$A_{\boldsymbol{n}}y_{\boldsymbol{n}}=b_{\boldsymbol{n}}$ in the Neumann BCs
case.

###  3.2 A multigrid method for a Toeplitz graph immersed in the disk

In this subsection we develop and study multigrid methods for the example in
[1, Subsection 7.1]. For the reader convenience, we briefly summarize
hereafter the mentioned example.

We consider the operator $\mathcal{L}$,

$\displaystyle\mathcal{L}:W^{1,2}_{0}\left(B_{1/2}\right)\to
L^{2}\left(B_{1/2}\right),$ (10)
$\displaystyle\mathcal{L}u(x,y):=-\textnormal{div}\left[p(x,y)\nabla
u(x,y)\right]+q(x,y)u(x,y)\qquad(x,y)\in B_{\frac{1}{2}},$ (11)

where

$B_{1/2}=\left\\{(x,y)\in\mathbb{R}^{2}\,:\,4\left(x-\frac{1}{2}\right)^{2}+4\left(y-\frac{1}{2}\right)^{2}<1\right\\},$

and $W^{1,2}_{0}\left(B_{1/2}\right)$ is the closure of the space of smooth
and compactly supported functions in the disk with respect to the Sobolev
norm. $\mathcal{L}$ is characterized by Dirichlet boundary conditions on the
boundary of the disk $B_{1/2}\subset[0,1]^{2}$. Fixing the diffusion term
$p(x,y)=1+(x-1/2)^{2}+(y-1/2)^{2}$ and the potential term
$q(x,y)=\textrm{e}^{xy}$, the discretization is made by an equispaced two-
dimensional central FD approximation with $5$-points. Fix
${\boldsymbol{n}}=(n,n)$. Then, the resulting graph $G_{{\boldsymbol{n}}}$
that approximates the underlying geometry is given by

$G_{{\boldsymbol{n}}}=\left(T_{{\boldsymbol{n}}}^{B_{1/2}}\left\langle\left\\{[1,0],w^{p}\right\\},\left\\{[0,1],w^{p}\right\\}\right\rangle,\kappa\right),$

where
$T_{{\boldsymbol{n}}}^{B_{1/2}}\left\langle\left\\{[1,0],w^{p}\right\\},\left\\{[0,1],w^{p}\right\\}\right\rangle$
is a $2$-level Toeplitz graph characterized by the node set

$V^{B_{1/2}}_{{\boldsymbol{n}}}=\left\\{(x_{i},y_{j})\in[0,1]^{2}\,:\,(x_{i},y_{j})\in
B_{\frac{1}{2}}\right\\},\quad(x_{i},y_{j})=\left(\frac{i}{n+1},\frac{j}{n+1}\right)\mbox{
for }i,j=1,\ldots,n,$

and the weight function

$w^{p}((x_{i},y_{j}),(x_{r},y_{s})):=p\left(\frac{x_{i}+x_{r}}{2},\frac{y_{j}+y_{s}}{2}\right)\quad\mbox{if
}(|i-r|,|j-s|)\in\\{(1,0),(0,1)\\}.$

The potential term $\kappa$ is the Dirichlet potential, as we already
introduced it in the previous Subsection 3.1. In particular,
$G_{{\boldsymbol{n}}}$ is a sub-graph of a graph $\bar{G}_{{\boldsymbol{n}}}$,
of the same structure of $G_{{\boldsymbol{n}}}$, but whose node set
$\bar{V}_{{\boldsymbol{n}}}$ is given by all the points $(x_{i},y_{j})$ in
$[0,1]^{2}$ and whose weight function $\bar{w}^{p}$ is obtained by extending
$p$ continuously outside the disk, fixing $p(x,y)\equiv 5/2$ for every
$(x,y)\in[0,1]^{2}\setminus B_{1/2}$. We have that

$\kappa(x_{i},y_{j})=\begin{cases}\kappa_{0}(x_{i},y_{j})=h^{2}q(x_{i},y_{j})&\mbox{if
}(x_{i},y_{j})\in\mathring{V}^{B_{1/2}}_{\boldsymbol{n}},\\\
\kappa_{1}(x_{i},y_{j})=h^{2}q(x_{i},y_{j})+\frac{5}{2},&\mbox{if
}\exists!\mbox{ one neighbor}\\\ &\mbox{in }\bar{V}_{\boldsymbol{n}}\setminus
V^{B_{1/2}}_{\boldsymbol{n}},\\\
\kappa_{2}(x_{i},y_{j})=h^{2}q(x_{i},y_{j})+\frac{10}{2},&\mbox{if
}\exists!\mbox{ two neighbors}\\\ &\mbox{in }\bar{V}_{\boldsymbol{n}}\setminus
V^{B_{1/2}}_{\boldsymbol{n}}.\\\ \end{cases}$

For all the details we refer to [1, Subsection 7.1]. See Figure 2.

$\kappa_{0}$$\kappa_{2}$$\kappa_{2}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{p}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{p}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$$\kappa_{1}$$\kappa_{1}$$\kappa_{1}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{p}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{p}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{\bar{p}}}$

Figure 2: Immersion of a $2$-level grid graph inside the disk
$B_{1/2}\subset[0,1]^{2}$. The white nodes are the nodes of
$V^{B_{1/2}}_{\boldsymbol{n}}$ while the gray nodes belong to
$\bar{V}_{\boldsymbol{n}}\setminus V^{B_{1/2}}_{\boldsymbol{n}}$. The
potential term of the host graph $\bar{G}_{n}$ is determined only by the
potential term $q$ from (10) while the potential term $\kappa$ of the sub-
graph $G_{\boldsymbol{n}}$ is influenced by the nodes in
$\bar{V}_{\boldsymbol{n}}\setminus V^{B_{1/2}}_{\boldsymbol{n}}$ on the
boundary set $\partial V^{B_{1/2}}_{\boldsymbol{n}}$. This influence is due to
the presence of Dirichlet BCs in (11). The green connections represent the
weighted edges whose end-nodes are both interior nodes of
$V^{B_{1/2}}_{\boldsymbol{n}}$, while the red connections represent the
weighted edges which have at least one end-node that belongs to
$\bar{V}_{\boldsymbol{n}}\setminus V^{B_{1/2}}_{\boldsymbol{n}}$. The
potential term $\kappa$ sums the weight of a red edge to every of its end-
nodes which belong to $V^{B_{1/2}}_{\boldsymbol{n}}$.

The graph Laplacian $\Delta_{G_{\boldsymbol{n}}}$, associated to the graph
$G_{\boldsymbol{n}}$, approximates the normalized operator
$(n+1)^{-2}\mathcal{L}$. By [1, Corollary 5.3] it immediately follows that

$\\{\Delta_{G_{\boldsymbol{n}}}\\}\sim_{\lambda}\mathfrak{f}(x,y,\theta_{1},\theta_{2}),\qquad(x,y,\theta_{1},\theta_{2})\in
B_{1/2}\times[-\pi,\pi]^{2}$

where

$\displaystyle\mathfrak{f}(x,y,\theta_{1},\theta_{2})$
$\displaystyle=p(x,y)\left(4-f(\theta_{1},\theta_{2})\right)$
$\displaystyle=\left[1+(x-1/2)^{2}+(y-1/2)^{2}\right]\left(4-2\cos(\theta_{1})-2\cos(\theta_{2})\right),$
(12)

and because of a symmetric argument, we can restrict $\mathfrak{f}$ on
$B_{1/2}\times[0,\pi]^{2}$ without affecting the validity of the identity (3).
Let us define the function $g:[-\pi,\pi]^{2}\rightarrow\mathbb{R}$

$g(\theta_{1},\theta_{2})=4-2\cos(\theta_{1})-2\cos(\theta_{2})$

and the function $a:[0,1]^{2}\rightarrow\mathbb{R}$ as the composition of the
potential term $p(x,y)=1+(x-1/2)^{2}+(y-1/2)^{2}$ with the following mapping
of the square $[0,1]^{2}$ into the disk $B_{1/2}$

$m(x,y)=\left(x\sqrt{1-\frac{1}{2}y^{2}},y\sqrt{1-\frac{1}{2}x^{2}}\right),$

that is,

$a(x,y)=p(m(x,y))=1+\left(x\sqrt{1-\frac{1}{2}y^{2}}-\frac{1}{2}\right)^{2}+\left(y\sqrt{1-\frac{1}{2}x^{2}}-\frac{1}{2}\right)^{2}.$

Let us first consider the multilevel Toeplitz matrices
$T_{\boldsymbol{n}}\left(g\right)$, with ${\boldsymbol{n}}=(n,n)$. According
to the discussion in Subsection 2.4 and Theorem 2.1, in order to obtain a two-
grid method that fulfills the approximation property for the linear systems

$T_{\boldsymbol{n}}\left(g\right)x_{\boldsymbol{n}}=b_{\boldsymbol{n}}$

we can choose the following prolongation operator

$P_{{\boldsymbol{n}},{\boldsymbol{k}}}=T_{\boldsymbol{n}}\left((2+2\cos(\theta_{1}))(2+2\cos(\theta_{2}))\right)K_{\boldsymbol{n}}.$
(13)

By Proposition 3.4 in [43], if we prove that

$T_{\boldsymbol{n}}\left(g\right)\leq
T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}},$
(14)

then $P_{{\boldsymbol{n}},{\boldsymbol{k}}}$ is a suitable prolongation
strategy also for the linear systems

$T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}x_{\boldsymbol{n}}=b_{\boldsymbol{n}}.$

By the Sylvester’s law of inertia, condition (14) is equivalent to

$I_{\boldsymbol{n}}\leq\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)$

and the latter inequality is fulfilled if and only if the function $a(x,y)$ is
greater than or equal to 1 on $[0,1]^{2}$, which is trivially verified.
According to Proposition 3.5 in [43], in order to choose a relaxation
parameter $\omega$ such that the Richardson method satisfies the smoothing
property we need to estimate the following quantity

$\displaystyle\sup_{{\boldsymbol{n}}}\rho\left(T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right)$
$\displaystyle=\sup_{{\boldsymbol{n}}}\left\|T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right\|_{2}$
$\displaystyle\leq\sup_{{\boldsymbol{n}}}\left\|T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right\|_{2}\left\|\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)\right\|_{2}\left\|T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right\|_{2}$
$\displaystyle<\sqrt{8}\cdot\frac{5}{4}\cdot\sqrt{8}=10.$

Then, we can take $\omega$ in the interval $\left(0,\frac{2}{10}\right]$. In
order to study the spectral features of the considered matrix-sequence we
employ the Generalized Locally Toeplitz (GLT) theory tools ([28]). Following
the notation of [1], we use the GLT properties GLT1, GLT3 and GLT4 to show
that the subsequent spectral distribution holds

$\left\\{T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right\\}_{\boldsymbol{n}}\sim_{\lambda}a(x,y)g(\theta_{1},\theta_{2}),$

where, in particular, we exploited the algebra property of GLT sequences,
which include Toeplitz and diagonal-sampling matrix-sequences, provided that
they fulfill specific hypothesis. By the very definition of spectral
distribution, it is immediate to see that the latter is equivalent to

$\left\\{T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\mathop{\mathrm{diag}}_{{\boldsymbol{n}}}\left(a\right)T_{\boldsymbol{n}}\left(g\right)^{\frac{1}{2}}\right\\}_{\boldsymbol{n}}\sim_{\lambda}\mathfrak{f}(x,y,\theta_{1},\theta_{2})$

where $\mathfrak{f}(x,y,\theta_{1},\theta_{2})$ is defined as in (3.2).

Consider now the multi-index ${\boldsymbol{n}}=(n,n)$ such that $n=2^{t}$,
$t\in\mathbb{N}$. From the definition of $V_{{\boldsymbol{n}}}$ in equation
(3.2), we recall that the dimension $d_{\boldsymbol{n}}$ of the graph
Laplacian $\Delta_{G_{\boldsymbol{n}}}$ is such that
$d_{\boldsymbol{n}}<n^{2}$. Taking inspiration from all the previous
considerations, we propose a two-grid method for the linear system

$\Delta_{G_{\boldsymbol{n}}}x_{\boldsymbol{n}}=b_{\boldsymbol{n}}$

having a prolongation operator of the form in equation (13). More precisely,
as a first step construct

$P_{{\boldsymbol{n}},{\boldsymbol{k}}}=T_{\boldsymbol{n}}\left((2+2\cos(\theta_{1}))(2+2\cos(\theta_{2}))\right)K_{\boldsymbol{n}}$

with ${\boldsymbol{k}}=(\frac{n}{2},\frac{n}{2})$. Then, eliminate from the
matrix $P_{{\boldsymbol{n}},{\boldsymbol{k}}}$ the rows
$\mathbf{i}=(i_{1},i_{2})$ such that $(x_{i_{1}},y_{i_{2}})$ does not belong
to $V_{{\boldsymbol{n}}}$ and eliminate the columns $\mathbf{j}=(j_{1},j_{2})$
such that $(x_{j_{1}},y_{j_{2}})$ does not belong to $V_{{\boldsymbol{k}}}$.
With this procedure, the matrix-sizes are consistent and it is also possible
to implement a V-cycle algorithm using the same prolongation/restriction
strategy at all levels.

In Table 5 we numerically show the validity of our proposed methods, reporting
the number of iterations needed for achieving the tolerance
$\varepsilon=10^{-6}$ when increasing the matrix size. In the first column, we
show the results relative to the two-grid method with pre-smoother and post-
smoother one iteration of the Richardson method with relaxation parameters
$1/5$ and $2/15$ respectively. In the second and third columns, we use Gauss-
Seidel for the two-grid and V-cycle algorithms respectively. In all cases, we
see that the number of iterations needed for convergence remains almost
constant, when increasing the size $d_{\boldsymbol{n}}$.

$t$ | $d_{\boldsymbol{n}}$ | Two-grid/Richardson | Two-grid/Gauss Seidel | V-cycle/Gauss Seidel
---|---|---|---|---
3 | 60 | 15 | 7 | 7
4 | 216 | 15 | 9 | 9
5 | 848 | 17 | 10 | 10
6 | 3300 | 17 | 9 | 10
Table 5: Comparison of the number of iterations, varying $t$, for the two-grid
and V-cycle methods applied to the graph Laplacian
$\Delta_{G_{\boldsymbol{n}}}$.

###  3.3 Elliptic problems on diamond Toeplitz graphs

All the examples in the previous subsections were about graphs arising from
PDEs approximations, and the PDEs were of the form of elliptic equations on
bounded subsets of $\mathbb{R}^{d}$. For this numerical example, we move into
a complete graph setting, where the graph can model many other different type
of real-world interactions without being necessarily restricted to be the
physical approximation of a geometric overlying domain.

Given a graph $G$ with node set $V$ such that $V=V_{1}\cup V_{2}$, $V_{1}\cap
V_{2}=\emptyset$, then the model problem we want to solve is

$\begin{cases}\Delta_{G}u(v_{i})=f(v_{i})&\mbox{if }v_{i}\in V_{1},\\\
u(v_{i})=h(v_{i})&\mbox{if }v_{i}\in V_{2}.\end{cases}$ (15)

This kind of problem can be viewed again as an elliptic, nonhomogeneous
Dirichlet problem (NHDP), but on graphs. Indeed, $\Delta_{G}$ is positive
semi-definite, and if $G$ is a subgraph of a graph $\bar{G}$ with node set
$\bar{V}$, then we can fix $V_{1}=\mathring{V}$ and $V_{2}=\partial V$, where
we recall that the (inner) boundary is defined as $\partial
V:=\left\\{v_{i}\in V\,:\,v_{i}\sim\bar{v}_{j}\mbox{ for some
}\bar{v}_{j}\in\bar{V}\setminus V\right\\}$ and
$\mathring{V}=V\setminus\partial V$. With this choice, Problem (15) becomes

$\begin{cases}\Delta_{G}u(v_{i})=f(v_{i})&\mbox{if }v_{i}\in\mathring{V},\\\
u(v_{i})=h(v_{i})&\mbox{if }v_{i}\in\partial V.\end{cases}$ (NHDP)

For simplicity, we will assume now that both $G$ and $\bar{G}$ have zero
potential terms, that is, $\kappa=\bar{\kappa}\equiv 0$. If we explicit now
the action of the graph Laplacian, we can see that, for every
$v_{i}\in\mathring{V}$,

$\displaystyle\Delta_{G}u(v_{i})$
$\displaystyle=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\ v_{j}\in
V\end{subarray}}w(v_{i},v_{j})\left(u(v_{i})-u(v_{j})\right)$
$\displaystyle=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in\mathring{V}\end{subarray}}w(v_{i},v_{j})\left(u(v_{i})-u(v_{j})\right)+\sum_{\begin{subarray}{c}v_{j}\sim
v_{i}\\\ v_{j}\in\partial
V\end{subarray}}w(v_{i},v_{j})\left(u(v_{i})-h(v_{j})\right)$
$\displaystyle=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in\mathring{V}\end{subarray}}w(v_{i},v_{j})\left(u(v_{i})-u(v_{j})\right)+\mathring{\kappa}(v_{i})u(v_{i})-g(v_{i}),$

where

$\mathring{\kappa}(v_{i})=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in\partial V\end{subarray}}w(v_{i},v_{j}),\qquad
g(v_{i})=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\ v_{j}\in\partial
V\end{subarray}}w(v_{i},v_{j})h(v_{j}).$

Observing that

$\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in\mathring{V}\end{subarray}}w(v_{i},v_{j})\left(u(v_{i})-u(v_{j})\right)+\mathring{\kappa}(v_{i})u(v_{i})$

defines the action of the graph Laplacian associated with the graph
$\mathring{G}:=\left(\mathring{V},\mathring{E},\mathring{w},\mathring{\kappa}\right)\subset
G$, where

$\mathring{E}=\left\\{(v_{i},v_{j})\in
E\,|\,v_{i},v_{j}\in\mathring{V}\right\\},\qquad\mathring{w}=w_{|\mathring{E}},$

then Problem (NHDP) is equivalent to solve

$\Delta_{\mathring{G}}u(v_{i})=g(v_{i})+f(v_{i}),\quad v_{i}\in\mathring{V}.$
(NHDP’)

Let us notice that the Dirichlet BCs in Problem (NHDP) have been absorbed by
the forcing term $g$ and the potential term $\mathring{\kappa}$, and that
$g\equiv 0$ if $h\equiv 0$ (zero Dirichlet BCs). Since

$\mathring{\kappa}(v_{i})=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in\partial
V\end{subarray}}w(v_{i},v_{j})=\sum_{\begin{subarray}{c}v_{j}\sim v_{i}\\\
v_{j}\in V\setminus\mathring{V}\end{subarray}}w(v_{i},v_{j}),$

again, we recover the Dirichlet potential as in Subsections 3.1 and 3.2.
$\mathring{\kappa}$ accounts for the edge deficiency of the node $v_{i}$ in
$\mathring{V}$, seen as a node in $V$. It emerges from the computation of the
action of the graph Laplacian associated to $G$, when imposing $u\equiv 0$ on
$\partial V=V\setminus\mathring{V}$.

For the numerical experiments, owing to the equivalence of (NHDP) and (NHDP’),
we work directly with $\mathring{G}$. Let us fix the mold graph
$M=\left([4],E,w\right)$ such that

$w(i,j)=\begin{cases}j-i&\mbox{if }$i=1$\mbox{ and }j=2,3,4,\\\
0&\mbox{otherwise}.\end{cases}$

The diamonds will be denoted by $M(k)=(V(k),E(k),w_{k})\simeq M$. We set

$L=\begin{bmatrix}10&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&1&0&0\end{bmatrix}$

and $G=T_{n,4}^{M}\left\langle(1,L)\right\rangle$ the corresponding $1$-level
diamond Toeplitz graph, see Figure 3. Finally, we set

  * •

$\partial V=V(1)\cup V(n)$, that is, $\partial V$ is given by the nodes that
belong to the node sets of the first and the last diamonds. Remember that
$V(k)=\\{v_{(k,1)},v_{(k,2)},v_{(k,3)},v_{(k,4)}\\}$;

  * •

$h(v_{(1,1)})=0.5,h(v_{(1,2)})=0.25,h(v_{(1,3)})=h(v_{(1,4)})=0$;

  * •

$h(v_{(n,1)})=0.5,h(v_{(n,2)})=0.25,h(v_{(n,3)})=h(v_{(n,4)})=0$;

  * •

$f(v_{(k,i)})=\sin(ki)$.

The symbol associated to the graph Laplacian $\Delta_{\mathring{G}}$ is given
by

${\boldsymbol{f}}(\theta)=D-\left[W+(L+L^{*})\cos(\theta)+(L-L^{*})\textup{i}\sin(\theta)\right]\in\mathbb{R}^{4\times
4},\qquad\theta\in[-\pi,\pi]$

where

$W=\begin{bmatrix}0&1&2&3\\\ 1&0&0&0\\\ 2&0&0&0\\\ 3&0&0&0\end{bmatrix},\qquad
D=\begin{bmatrix}26&0&0&0\\\ 0&2&0&0\\\ 0&0&2&0\\\ 0&0&0&4\end{bmatrix},$

which gives in compact form

${\boldsymbol{f}}(\theta)=\begin{bmatrix}26-20\cos(\theta)&-1&-2&-3\\\
-1&2&0&-\cos(\theta)+\textup{i}\sin(\theta)\\\ -2&0&2&0\\\
-3&-\cos(\theta)-\textup{i}\sin(\theta)&0&4\end{bmatrix}.$

$v_{(2,1)}$
---
$v_{5}$
$v_{(2,4)}$
---
$v_{8}$
$v_{(2,2)}$
---
$v_{6}$
$v_{(2,3)}$
---
$v_{7}$
$v_{(3,4)}$
---
$v_{12}$
$v_{(3,3)}$
---
$v_{11}$
$v_{(3,2)}$
---
$v_{10}$
$v_{(3,1)}$
---
$v_{9}$
$v_{(1,2)}$
---
$v_{2}$
$v_{(1,3)}$
---
$v_{3}$
$v_{(1,4)}$
---
$v_{4}$
$v_{(1,1)}$
---
$v_{1}$
$3$$1$$2$$3$$2$$1$$1$$3$$2$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}10}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}10}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}1}$
Figure 3: Representation of the $1$-level diamond Toeplitz graph
$G=T_{n,4}^{M}\left\langle(1,L)\right\rangle$ for $n=3$.

In order to study the zeros of the four eigenvalue functions of
${\boldsymbol{f}}(\theta)$, we numerically checked that
$0\leq\lambda_{1}({\boldsymbol{f}}(\theta))<\lambda_{2}({\boldsymbol{f}}(\theta))<\lambda_{3}({\boldsymbol{f}}(\theta))<\lambda_{4}({\boldsymbol{f}}(\theta))$
for all $\theta\in[-\pi,\pi]$ and we then computed the determinant, which is
equal to $\det({\boldsymbol{f}}(\theta))=292-292\cos(\theta)$. Hence, we
deduce that both the determinant and $\lambda_{1}({\boldsymbol{f}}(\theta))$
have a zero of order 2 in 0. As a consequence, we define a grid transfer
operator of the form (4), associated to the following matrix-valued
trigonometric polynomial

$\mathbf{p}(\theta)=(4+6\cos(\theta)+4\cos(2\theta)+2\cos(3\theta))\begin{bmatrix}2&1&1&1\\\
1&2&1&1\\\ 1&1&2&1\\\ 1&1&1&2\\\ \end{bmatrix}.$

Indeed, the trigonometric polynomial
$2+6\cos(\theta)+4\cos(2\theta))+2\cos(3\theta)$ has zeros of order 2 in
$\pi/2$ and $\pi$, which are the mirror points of 0 if we consider a
coarsening factor of 4, as in Theorem 2.1. Moreover, the matrix

$\begin{bmatrix}2&1&1&1\\\ 1&2&1&1\\\ 1&1&2&1\\\ 1&1&1&2\\\ \end{bmatrix}$

guarantees a constant convergence rate of the V-cycle method, according to the
analysis in [21].

In Table 6 we numerically show the validity of our proposed methods, reporting
the number of iterations needed for achieving the tolerance
$\varepsilon=10^{-6}$ when increasing the matrix size. In this case, we use
the relative error stopping criterion, where the exact solution is computed
using the LAMG library (v.2.2.1), whose run time and storage were empirically
demonstrated to scale linearly with the number of edges [34].

We use one iteration of Gauss-Seidel both as pre-smoother and post-smoother.
In the first and second columns, we show the results relative to the two-grid
and V-cycle methods for $\mathfrak{g}=2$, while in the third and fourth
columns we consider a coarsening factor $\mathfrak{g}=4$.

In all cases, we see that the number of iterations needed for convergence
remains almost constant, when increasing the size $d_{\boldsymbol{n}}$.

$t$ | $d_{\boldsymbol{n}}$ | Two-grid ($\mathfrak{g}=2$) | V-cycle ($\mathfrak{g}=2$) | Two-grid ($\mathfrak{g}=4$) | V-cycle ($\mathfrak{g}=4)$
---|---|---|---|---|---
4 | 1016 | 5 | 6 | 16 | 16
5 | 4088 | 5 | 6 | 20 | 22
6 | 16376 | 5 | 6 | 19 | 23
7 | 65528 | 5 | 6 | 19 | 24
8 | 262136 | 5 | 6 | 20 | 25
Table 6: Comparison of the number of iterations, varying $t$, for the two-grid
and V-cycle methods applied to the graph Laplacian
$\Delta_{G_{\boldsymbol{n}}}$.

###  3.4 FEM and IgA approximations of elliptic PDEs

In this subsection we perform some numerical tests concerning examples
reported in [1, Subsections 7.2, 7.3]. For the sake of completeness we briefly
recall the setting.

We consider the model problem

$-\Delta u=g\text{ in }\ \Omega,$ (16)

where $\Omega=(0,1)$ and $g\in\text{L}^{2}(\Omega)$. We first approximate (16)
by using quadratic FEs over the uniform mesh with stepsize $\frac{1}{n+1}$,
where we choose as FEs basis the quadratic $C^{0}$ B-spline basis over the
knot sequence
$\\{\frac{1}{n+1},\frac{1}{n+1},\frac{2}{n+1},\frac{2}{n+1},...,\frac{n}{n+1},\frac{n}{n+1}\\}$.
Fixed ${\boldsymbol{n}}=(n,n)$ and proceeding as in [30], it is possible to
show that, in this case, the approximation of the operator $-n^{-1}\Delta$ is
given by the graph Laplacian $\Delta_{G_{{\boldsymbol{n}}}}$ of a $1$-level
diamond Toeplitz graph
$T_{{\boldsymbol{n}},2}^{G}\langle\left(1,L_{1}\right)\rangle$, where

$W=\left(\begin{array}[]{cc}0&2\\\ 2&0\end{array}\right)$

is the adjacency matrix of the mold graph $G$ and

$L_{1}=\left(\begin{array}[]{cc}0&0\\\ 2&2\end{array}\right).$

Note that (see [1]) we have a nonzero potential term which depends on the
choice of the boundary conditions.
It is then easy to see that the sequence $\\{\Delta_{G_{{\boldsymbol{n}}}}\\}$
satisfies

$\\{\Delta_{G_{{\boldsymbol{n}}}}\\}\sim_{\lambda}\boldsymbol{\mathfrak{f}}(\theta),$

where $\boldsymbol{\mathfrak{f}}:[-\pi,\pi]\to\mathbb{C}^{2\times 2}$ is given
by

$\boldsymbol{\mathfrak{f}}(\theta)=\frac{1}{3}\left\\{\left(\begin{array}[]{cc}4&-2\\\
-2&8\end{array}\right)+\left(\begin{array}[]{cc}0&-2\\\
-2&-4\end{array}\right)\cos(\theta)+\left(\begin{array}[]{cc}0&-2\\\
2&0\end{array}\right)\textup{i}\sin(\theta)\right\\}.$

It is now well known from the relevant literature (see [15, 17, 18, 19, 20])
that an optimal preconditioner for the linear system

$\Delta_{G_{{\boldsymbol{n}}}}u_{{\boldsymbol{n}}}=g_{{\boldsymbol{n}}}$ (17)

is given by the matrix
$S_{n}=C_{n}(\boldsymbol{\mathfrak{f}})+\frac{1}{2n}ee^{T}$, $e$ being the
column vector of dimension $2n$ of all ones. This is numerically confirmed in
Table 7.

$t$ |  $d_{\boldsymbol{n}}$($n=2^{t}$) | CG | PCG
---|---|---|---
6 | 64 | 95 | 4
7 | 128 | $>100$ | 4
8 | 256 | $>100$ | 4
9 | 512 | $>100$ | 5
10 | 1024 | $>100$ | 5
11 | 2048 | $>100$ | 5
Table 7: Comparison of the number of iterations, varying $t$, by the CG method
and the PCG method with preconditioner $S_{n}$ with tolerance $10^{-6}$.

If we now consider the discretization of (16) arising from an IgA approach as
in [1, Subsection 7.3] (see also [20]), the discretizing operator
$\Delta_{G_{n}}$ is not the graph Laplacian of one of the structures described
in Section 2; however, the difference is due to the presence of local
perturbations near the boundary nodes which do not affect the symbol function.
In particular, it happens that the sequence $\\{\Delta_{G_{n}}\\}$ has the
same symbol function as the graph Laplacian of the $1$-level Toeplitz graph

$H_{n}=\left(T_{n}\left\langle\left(1,\frac{30}{240}\right),\left(1,\frac{48}{240}\right),\left(1,\frac{2}{240}\right)\right\rangle,\kappa\right),\quad\kappa\equiv
0,$

that is

$\\{\Delta_{G_{n}}\\}\sim_{\lambda}\mathfrak{f}(\theta),$

with
$\mathfrak{f}(\theta)=\frac{160}{240}-\frac{60}{240}\cos(\theta)-\frac{96}{240}\cos(2\theta)-\frac{4}{240}\cos(3\theta)$.
By direct computation, it is easy to see that $\mathfrak{f}$ has a unique zero
of order $2$ at $0$. It follows that, if we choose
$p(\theta)=2-2\cos(\theta)$, we deduce that
$\frac{\mathfrak{f}(\theta)}{p(\theta)}$ is monotonic decreasing in $[0,\pi]$,
globally even, and

$\lim_{\theta\to 0}\frac{\mathfrak{f}(\theta)}{p(\theta)}=1.$

The previous formula suggests as preconditioner for the PCG the matrix
$P_{n}=T_{n}(2-2\cos(\theta))$, for which, following the theory developed in
[17], we have

$\displaystyle\lim_{n\to\infty}\lambda_{\min}(P_{n}^{-1}A_{n})$
$\displaystyle=$
$\displaystyle\min_{\theta\in[0,\pi]}\frac{\mathfrak{f}(\theta)}{p(\theta)}=\frac{\mathfrak{f}(\pi)}{p(\pi)}=\frac{2}{15}$
$\displaystyle\lim_{n\to\infty}\lambda_{\max}(P_{n}^{-1}A_{n})$
$\displaystyle=$
$\displaystyle\sup_{\theta\in[0,\pi]}\frac{\mathfrak{f}(\theta)}{p(\theta)}=\lim_{\theta\to
0}\frac{\mathfrak{f}(\theta)}{p(\theta)}=1.$

As we show in Table 8, the PCG with preconditioner
$P_{n}=T_{n}(2-2\cos(\theta))$ for the solution of the linear system (17) is
in fact optimal.

$t$ |  $d_{\boldsymbol{n}}$($n=2^{t}$) | CG | PCG
---|---|---|---
7 | 128 | 72 | 20
8 | 256 | $>100$ | 20
9 | 512 | $>100$ | 21
10 | 1024 | $>100$ | 21
11 | 2048 | $>100$ | 22
12 | 4096 | $>100$ | 22
Table 8: Comparison of the number of iterations, varying $t$, for the IgA
example by the CG method and the PCG method with preconditioner $P_{n}$ with
tolerance $10^{-6}$ .

In what follows, we show an implementation of the MGM with a projector matrix
whose choice depends on the analytic properties of the symbol function $f$. We
consider now the same examples as above in this subsection. In particular we
emphasize that our analysis allows us to define a projector for which the MGM
turns out to be an optimal method when solving system (17). Further, since the
discrete operator $\Delta_{G_{n}}$ is HPD, we use the Richardson method both
as pre-smoother and post-smoother.
We start with the IgA example and assume that $n=2^{t}-1$, $t\geq 1$.
Following the ideas of Subsection 2.4, we can build our projector as
$P_{n}=T_{n}(p(\theta))K_{n}$, with $p(\theta)=2+2\cos(\theta)$ and $K_{n}$
the cutting matrix defined as in (5). Recall that $f$ has a unique zero in
$[0,\pi]$ attained at $0$. It follows that the trigonometric polynomial $p$
defined above is such that

$\lim_{x\to 0}\frac{p^{2}(\pi-x)}{f(x)}<\infty$

and

$p^{2}(x)+p^{2}(\pi-x)>0\quad\forall x\in[0,\pi],$

where the two conditions above insure the optimality in light of Theorem 2.1.
As noted above, we know that $\Delta_{G_{n}}$ has not the structure of the
graph Laplacian of a Toeplitz graph. We denote by $R_{n}$ the difference
$\Delta_{G_{n}}-\Delta{H_{n}}$, where $\Delta{H_{n}}$ is the graph Laplacian
of the Toeplitz graph $H_{n}$ in the previous subsection, which is a low rank
correction. In order to ensure convergence, in numerical examples we use
$0<\omega<\frac{2}{M_{f}+\rho(R_{n})}$ as the scalar parameter in Richardson
method one iteration of the Richardson method with relaxation parameters
$0.7149$ and $1.4299$ respectively, where $M_{f}\approx 0.5333$ is a maximum
value of $f(\theta)$ and $\rho(R_{n})\approx 0.8654$ is the spectral radius of
$R_{n}$, thus satisfying condition in Proposition 3.5 in [43]. In Table 9 we
numerically show the validity of our proposed methods.

$t$ |  $d_{\boldsymbol{n}}$($n=2^{t}-1$) | Two-grid Method
---|---|---
7 | 127 | 8
8 | 255 | 8
9 | 511 | 8
10 | 1023 | 8
11 | 2047 | 8
12 | 4095 | 8
Table 9: Number of iterations, varying $t$, for the two-grid method applied to
the IgA example with tolerance $10^{-6}$.

An analogous study can be carried for the discretization via the FEM approach.

In this case we assume $n=2^{t}$, with $t$ integer. By direct computation, it
is possible to show that the same projector $P_{n}$ can be used also in this
situation; we only need to pay attention to the dimension of the projector
(which is going to be $2n\times n$ due to the block structure of the matrix
involved) and, as above, we use both as pre-smoother and post-smoother one
iteration of the Richardson method with relaxation parameters $0.25$ and
$0.50$ respectively. In Table 10, we show the optimality of the considered
method.

$t$ |  $d_{\boldsymbol{n}}$($n=2^{t}$) | Two-grid Method
---|---|---
7 | 128 | 9
8 | 256 | 9
9 | 512 | 9
10 | 1024 | 9
11 | 2048 | 9
12 | 4096 | 9
Table 10: Number of iterations, varying $t$, for the two-grid method applied
to the FEs example with tolerance $10^{-6}$.

##  4 Conclusions

In the present work we have treated sequences of graphs having a grid
geometry, with a uniform local structure in a bounded domain
$\Omega\subset{\mathbb{R}}^{d}$, $d\geq 1$. When $\Omega=[0,1]$, such graphs
include the standard Toeplitz graphs and, for $\Omega=[0,1]^{d}$, the
considered class includes $d$-level Toeplitz graphs. In the general case, the
underlying sequence of adjacency matrices has a canonical eigenvalue
distribution, in the Weyl sense, and it has been shown in the theoretical part
of this work that we can associate to it a symbol $\boldsymbol{\mathfrak{f}}$,
also in the case of variable coefficients in connection with the notion of
Generalized Locally Toeplitz sequences.

Here we have given practical evidence that the knowledge of the symbol and of
its basic analytical features provided enough information on the eigenvalue
structure in terms of localization, spectral gap, clustering, and global
distribution, in order to design efficient numerical methods for the
corresponding large linear systems and for approximating in a fast way the
eigenvalue of continuous differential operators. Tests and applications have
been taken from the approximation of differential operators via (local)
numerical schemes such as Finite Differences (FDs), Finite Elements (FEs), and
Isogeometric Analysis (IgA), where the differential operator domains are of
non-cartesian nature.

Nevertheless, more applications can be taken into account, since the results
presented here can be applied as well to study the spectral properties of
adjacency matrices and Laplacian operators of general large graphs and
networks, whenever the involved matrices enjoy a uniform local structure. The
extension of the study concerning this last issue will be a specific direction
to be investigated in future works.

## References

  * [1] A. Adriani, D. Bianchi, S. Serra-Capizzano, Asymptotic spectra of large (grid) graphs with a uniform local structure (part I): theory. Milan J. Math. 88: 409–454 (2020).
  * [2] A. Aricò, M. Donatelli, S. Serra-Capizzano, V-cycle optimal convergence for certain (multilevel) structured linear systems. SIAM J. Matrix Anal. Appl. 26: 186–214 (2004).
  * [3] G. Barbarino, A systematic approach to reduced GLT. BIT: 1–63 (2021).
  * [4] G. Barbarino, D. Bianchi, C. Garoni, Constructive approach to the monotone rearrangement of functions. Expo. Math. (2021), in press.
  * [5] G. Barbarino, C. Garoni, S. Serra-Capizzano, Block Generalized Locally Toeplitz sequences: theory and applications in the unidimensional case. Electron. Trans. Numer. Anal. 53: 28–112 (2020).
  * [6] G. Barbarino, C. Garoni, S. Serra-Capizzano, Block Generalized Locally Toeplitz sequences: theory and applications in the multidimensional case. Electron. Trans. Numer. Anal. 53: 113–216 (2020).
  * [7] D. Bianchi, Analysis of the spectral symbol function for discretization of a linear differential operator and analysis of the relative spectrum, with applications. Calcolo 58: paper 38, 47 pp. (2021).
  * [8] D. Bianchi, M. Donatelli, Graph approximation and generalized Tikhonov regularization for signal deblurring. Proceedings in IEEE ICCSA (2021): The International Conference on Computational Science and its Applications. To appear.
  * [9] D. Bianchi, S. Serra-Capizzano, Spectral analysis of finite-dimensional approximations of 1d waves in non-uniform grids. Calcolo 55: paper 47, 28 pp. (2018).
  * [10] D. Bianchi, A.G. Setti, R.K. Wojciechowski. The generalized porous medium equation on graphs: existence and uniqueness of solutions with $\ell^{1}$ data. Preprint.
  * [11] A. Böttcher, B. Silbermann, Introduction to Large Truncated Toeplitz Matrices. Springer-Verlag, New York (1999).
  * [12] B. Brosowski, F. Deutsch, An elementary proof of the Stone-Weierstrass Theorem. Proc. Amer. Math. Soc. 81(1) (1981).
  * [13] P. Ciarlet, The Finite Element Method for Elliptic Problems. North Holland, Amsterdam (1978).
  * [14] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric analysis: toward integration of CAD and FEA. John Wiley & Sons (2009).
  * [15] D. Cvetkovic, M. Doob, H. Sachs, Spectra of Graphs. Academic Press, New York (1979).
  * [16] E.B. Davies, Spectral theory and differential operators. Cambridge University Press 42 (1996).
  * [17] F. Di Benedetto, G. Fiorentino, S. Serra-Capizzano, CG preconditioning for Toeplitz matrices. Comput. Math. Appl. 25(6): 35–45 (1993).
  * [18] M. Donatelli, C. Garoni, C. Manni, S. Serra-Capizzano, H. Speleers, Robust and optimal multi-iterative techniques for IgA Galerkin linear systems. Comput. Methods Appl. Mech. 284: 230–264 (2015).
  * [19] M. Donatelli, C. Garoni, C. Manni, S. Serra-Capizzano, H. Speleers, Two-grid optimality for Galerkin linear systems based on B-splines. Comput. Vis. Sci. 17: 119–133 (2015).
  * [20] M. Donatelli, C. Garoni, C. Manni, S. Serra-Capizzano, H. Speleers, Symbol-based multigrid methods for Galerkin B-spline isogeometric analysis. SIAM J. Numer. Anal. 55(1): 31–62 (2017).
  * [21] M. Donatelli, P. Ferrari, I. Furci, S. Serra-Capizzano, D. Sesana, Multigrid methods for block-Toeplitz linear systems: convergence analysis and applications. Numer. Linear Algebra Appl. 28(4): e2356 (2021).
  * [22] M. Donatelli, M. Mazza, S. Serra-Capizzano, Spectral analysis and multigrid methods for finite volume approximations of space-fractional diffusion equations. SIAM J. Sci. Comput. 40(6): A4007–A4039 (2018).
  * [23] M. Donatelli, S. Serra-Capizzano, D. Sesana, Multigrid methods for Toeplitz linear systems with different size reduction. BIT 52(2): 305–327 (2012).
  * [24] S.-E. Ekström, I. Furci, C. Garoni, C. Manni, S. Serra-Capizzano, H. Speleers, Are the eigenvalues of the $B$-spline isogeometric analysis approximation of $-\Delta u=\lambda u$ known in almost closed form? Numer. Linear Algebra Appl. 25(5) (2018): e2198.
  * [25] E. Estrada, Path Laplacian matrices: introduction and application to the analysis of consensus in networks. Linear Algebra Appl., 436: 3373–3391 (2012).
  * [26] G. Fiorentino, S. Serra-Capizzano, Multigrid methods for Toeplitz matrices. Calcolo 28: 283–305 (1991).
  * [27] C. Garoni, M. Mazza, S. Serra-Capizzano, Block Generalized Locally Toeplitz sequences: from the theory to the applications. Axioms 7(3): 49 (2018).
  * [28] C. Garoni, S. Serra-Capizzano, The theory of Generalized Locally Toeplitz sequences: theory and applications - Vol I. Springer - Springer Monographs in Mathematics, New York (2017).
  * [29] C. Garoni, S. Serra-Capizzano, The theory of multilevel Generalized Locally Toeplitz sequences: theory and applications - Vol II. Springer - Springer Monographs in Mathematics, New York (2018).
  * [30] C. Garoni, S. Serra-Capizzano, D. Sesana, Spectral analysis and spectral symbol of $d$-variate $\mathbb{Q}_{\boldsymbol{p}}$ Lagrangian FEM stiffness matrices. SIAM J. Matrix Anal. Appl. 36(3): 1100–1128 (2015).
  * [31] C. Garoni, S. Serra-Capizzano, D. Sesana, The Theory of Block Generalized Locally Toeplitz Sequences. Technical Report, N. 1, January 2018, Department of Information Technology, Uppsala University.
  * [32] U. Grenander, G. Szegö, Toeplitz Forms and Their Applications. 2nd ed., Chelsea, New York (1984).
  * [33] S. Hossein Ghorban, Toeplitz graph decomposition. Trans. Combinat. 1(4): 35–41 (2012).
  * [34] O.E. Livne, A. Brandt, Lean algebraic multigrid (LAMG): fast graph Laplacian linear solver. SIAM J. Sci. Comput. 34(4): 499–522 (2012).
  * [35] J. Lund, K.L. Bowers, Sinc methods for quadrature and differential equations. SIAM, (1992).
  * [36] M. Keller, D. Lenz, R.K. Wojciechowski, Graphs and discrete Dirichlet spaces. Grundlehren der mathematischen Wissenschaften. To appear.
  * [37] C. W. Oosterlee, A. Schüller, U. Trottenberg, Multigrid. Academic Press, London (2001).
  * [38] J.W. Ruge, K. Stüben, Multigrid methods. SIAM, Philadelphia, PA (1987).
  * [39] Y. Saad, Iterative methods for sparse linear systems. $2^{\rm nd}$ Edition, SIAM, Philadelphia, PA, (2003).
  * [40] S. Serra-Capizzano, Convergence analysis of two-grid methods for elliptic Toeplitz and PDEs matrix-sequences. Numer. Math. 92(3): 433–465 (2002).
  * [41] S. Serra-Capizzano, Generalized Locally Toeplitz sequences: spectral analysis and applications to discretized partial differential equations. Linear Algebra Appl. 366: 371–402 (2003).
  * [42] S. Serra-Capizzano, The GLT class as a generalized Fourier analysis and applications, Linear Algebra Appl. 419: 180–233 (2006).
  * [43] S. Serra-Capizzano, C. Tablino-Possio, Two-Grid Methods for Hermitian positive definite linear systems connected with an order relation, Calcolo. 51(2): 261–285 (2014).
  * [44] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations. Chapman and Hall, International Thompson Publ., New York (1989).
  * [45] P. Tilli, Locally Toeplitz matrices: spectral theory and applications. Linear Algebra Appl. 278: 91–120 (1998).
  * [46] P. Tilli, A note on the spectral distribution of Toeplitz matrices. Linear Multilin. Algebra 45: 147–159 (1998).
  * [47] E. Tyrtyshnikov, N. Zamarashkin, Spectra of multilevel Toeplitz matrices: advanced theory via simple matrix relationships. Linear Algebra Appl. 270: 15–27 (1998).