1 Introduction↩︎

Throughout this paper, a sequence of matrices is, by definition, a sequence of the form \(\{A_n\}_n\), where \(n\) varies in some infinite subset of \(\mathbb{N}\) and \(A_n\) is a matrix of size \(d_n\times e_n\) such that both \(d_n\) and \(e_n\) tend to \(\infty\) as \(n\to\infty\). The theory of generalized locally Toeplitz (GLT) sequences is an apparatus for computing the spectral and singular value distribution of sequences of matrices that possess a (possibly hidden) Toeplitz-like structure. Sequences of this kind, which are known as GLT sequences, arise in several applications, including the discretization of differential and integral equations. Nowadays, the theory of GLT sequences is a fairly extensive research area with numerous applications. For readers who are new to the subject, we recommend the introduction [1] and the six-page conference paper [2]. For a comprehensive exposition of the topic, we refer to the books [3], [4] and the book-like papers [5]–[8]. Recent noteworthy developments not published in the previous references include the identification between spaces of GLT sequences and function spaces [9], [10], the characterization of diagonal matrix sequences enjoying a spectral distribution in terms of GLT sequences [11], the spectral distribution result for non-Hermitian perturbations of Hermitian GLT sequences [12], the GLT classification of diagonal sampling matrix sequences obtained from quasi-uniform samples of almost everywhere continuous functions [13], the derivation of a “normal form” for GLT sequences [14], and the successful application of matrix-less spectral approximation methods to GLT sequences [15]. It is also worth noting that, as demonstrated by recent research, the spectral distribution of a GLT sequence has significant practical implications. For example, suppose that \(\{A_n\}_n\) is a GLT sequence resulting from the discretization of a differential equation \(\mathscr Au=f\) through a given numerical method. Then, the spectral distribution of \(\{A_n\}_n\) can be used to measure the accuracy of the method in approximating the spectrum of the differential operator \(\mathscr A\) [16], to establish whether the method preserves the so-called average spectral gap [17], or to formulate analytical predictions for the eigenvalues of both \(A_n\) and \(\mathscr A\) [18]. Moreover, the spectral distribution of \(\{A_n\}_n\) can be exploited to design efficient iterative solvers for linear systems with matrix \(A_n\) and to analyze/predict their performance; see [19], [20] for accurate convergence estimates of Krylov methods based on the spectral distribution and [3] for more details on this subject.

Associated with any GLT sequence \(\{A_n\}_n\) is a special function \(\kappa\), the so-called symbol. Given a GLT sequence \(\{A_n\}_n\) with symbol \(\kappa\) and a vector \((n_1,\ldots,n_\nu)\) of fixed length \(\nu\) such that \(n_1+\ldots+n_\nu={\rm size}(A_n)\), we consider two sequences of preconditioners for \(A_n\):

• the sequence of block Jacobi preconditioners for \(A_n\) associated with \((n_1,\ldots,n_\nu)\), i.e., the sequence of matrices \[P_n^{J}=\mathop{\rm blockdiag}_{n_1,\ldots,n_\nu}(A_n)=\begin{bmatrix}A_n^{(11)} & & \\ & \ddots & \\ & & A_n^{(\nu\nu)}\end{bmatrix},\] where \(A_n^{(kk)}\) is the \(k\)th diagonal block of \(A_n\) of size \(n_k\);

• the sequence of block Gauss–Seidel preconditioners for \(A_n\) associated with \((n_1,\ldots,n_\nu)\), i.e., the sequence of matrices \[P_n^{GS}=\mathop{\rm blocktril}_{n_1,\ldots,n_\nu}(A_n)=\begin{bmatrix}A_n^{(11)} & & \\ \vdots & \ddots & \\ A_n^{(\nu1)} & \cdots & A_n^{(\nu\nu)}\end{bmatrix},\] where \(A_n^{(pq)}\) is the block of \(A_n\) of size \(n_p\times n_q\) in position \((p,q)\).

In our main results (Theorems 2–3), we prove that, regardless of \((n_1,\ldots,n_\nu)\), \(\{P_n^J\}_n\) and \(\{P_n^{GS}\}_n\) are GLT sequences with symbol \(\kappa\), just like \(\{A_n\}_n\). This result has two noteworthy implications.

• First, it gives reasons to believe that block Jacobi/Gauss–Seidel preconditioning for matrices \(A_n\) belonging to a GLT sequence is very efficient; see Remark 1. The numerical experiments presented in this paper reveal that this is in fact the case.

• Secondly, it allows us to give an alternative proof of the Fasino–Tilli theorem on the zero distribution of sequences of Hankel matrices generated by \(L^1\) functions (Theorem 4). Actually, it also allows us to extend the Fasino–Tilli theorem to a larger class of matrix sequences that we call generalized locally Hankel (GLH) sequences; see Theorem 5.

The paper is organized as follows. In Section 2, we give an overview of the theory of GLT sequences with a focus on the results that we need in this paper. In Section 3, we extend to the case of GLT sequences with matrix-valued symbols a well-known preconditioning result, which will be used in combination with our first main result (Theorem 2) to theoretically justify the efficiency of block Jacobi preconditioners for GLT sequences. In Section 4, we state and prove our main results described above (Theorems 2–3). In Section 5, as an application of Theorem 2, we give an alternative proof of the Fasino–Tilli theorem on the zero distribution of sequences of Hankel matrices generated by \(L^1\) functions, and we also extend this theorem to GLH sequences. In Section 6, we illustrate through numerical examples the efficiency of block Jacobi and block Gauss–Seidel preconditioners for GLT sequences. In Section 7, we draw conclusions and suggest future lines of research, which include, among others, the application of our main results to the analysis of additive and multiplicative Schwarz methods in the context of domain decomposition techniques.

2 Overview of the theory of GLT sequences↩︎

In this section, we give an overview of the theory of GLT sequences. For conciseness purposes, we only present the results that we need in this paper. For a comprehensive exposition of the topic, see [3]–[8]. For an introduction to the subject, we recommend [1], [2].

Singular value and spectral distribution of a sequence of matrices. Let \(\mu_k\) be the Lebesgue measure in \(\mathbb{R}^k\). Throughout this paper, all terminology from measure theory (such as “measurable set”, “measurable function”, “a.e.”, etc.)always refers to the Lebesgue measure. A matrix-valued function \(f:D\subseteq\mathbb{R}^k\to\mathbb{C}^{s\times t}\) is said to be measurable (respectively, continuous, continuous a.e., in \(L^1(D)\), etc.)if its components \(f_{ij}:D\to\mathbb{C}\), \(i=1,\ldots,s\), \(j=1,\ldots,t\), are measurable (respectively, continuous, continuous a.e., in \(L^1(D)\), etc.). For every \(x,y\in\mathbb{R}\), we define \(x\wedge y=\min(x,y)\). We denote by \(C_c(\mathbb{C})\) (respectively, \(C_c(\mathbb{R})\)) the space of complex-valued functions defined on \(\mathbb{C}\) (respectively, \(\mathbb{R}\)) with bounded support. The singular values of a matrix \(A\in\mathbb{C}^{m\times n}\) are denoted by \(\sigma_i(A)\), \(i=1,\ldots,m\wedge n\), and the eigenvalues of a matrix \(A\in\mathbb{C}^{m\times m}\) are denoted by \(\lambda_i(A)\), \(i=1,\ldots,m\). The spectrum of a matrix \(A\in\mathbb{C}^{m\times m}\) is denoted by \(\Lambda(A)\). Recall that a sequence of matrices is, by definition, a sequence of the form \(\{A_n\}_n\), where \(n\) varies in some infinite subset of \(\mathbb{N}\) and \(A_n\) is a matrix of size \(d_n\times e_n\) such that both \(d_n\) and \(e_n\) tend to \(\infty\) as \(n\to\infty\).

Note that Definition 1 is well-posed by [7], which ensures that the functions \({\boldsymbol{x}}\mapsto\sum_{i=1}^{s\wedge t}F(\sigma_i(f({\boldsymbol{x}})))\) and \({\boldsymbol{x}}\mapsto\sum_{i=1}^{s}F(\lambda_i(f({\boldsymbol{x}})))\) are measurable. We refer the reader to [7] or [21] for the informal meaning behind the singular value and spectral distribution of a sequence of matrices.

The following lemma will be used in Section 3; see [7] for the corresponding proof. Throughout this paper, the closure of a set \(S\) is denoted by \(\overline{S}\).

Zero-distributed sequences. A zero-distributed sequence is a sequence of matrices \(\{Z_n\}_n\) such that \(\{Z_n\}_n\sim_\sigma0\), i.e., \[\lim_{n\to\infty}\frac{1}{d_n\wedge e_n}\sum_{i=1}^{d_n\wedge e_n}F(\sigma_i(Z_n))=F(0),\qquad\forall\,F\in C_c(\mathbb{R}),\] where \(d_n\times e_n\) is the size of \(Z_n\).

Toeplitz sequences. If \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) is in \(L^1([-\pi,\pi])\), its Fourier coefficients are denoted by \(\{f_k\}_{k\in\mathbb{Z}}\) and are defined as follows: \[\label{Fc} f_k=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(\theta){\rm e}^{-{\rm i}k\theta}{\rm d}\theta\:\in\:\mathbb{C}^{s\times t},\qquad k\in\mathbb{Z},\tag{1}\] where the integrals are computed componentwise. The \(n\)th (block) Toeplitz matrix generated by a function \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) belonging to \(L^1([-\pi,\pi])\) is the \(ns\times nt\) matrix denoted by \(T_n(f)\) and defined as follows: \[T_n(f)=[f_{i-j}]_{i,j=1}^n.\] Any sequence of matrices of the form \(\{T_{d_n}(f)\}_n\), with \(d_n\to\infty\) and \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) in \(L^1([-\pi,\pi])\), is referred to as a (block) Toeplitz sequence generated by \(f\).

Sequences of diagonal sampling matrices. If \(a:[0,1]\to\mathbb{C}^{s\times t}\), the \(n\)th (block) diagonal sampling matrix generated by \(a\) is the \(ns\times nt\) (block) diagonal matrix denoted by \(D_n(a)\) and defined as follows: \[D_n(a)=\mathop{\rm diag}_{i=1,\ldots,n}a\Bigl(\frac{i}{n}\Bigr).\] Any sequence of matrices of the form \(\{D_{d_n}(a)\}_n\), with \(d_n\to\infty\) and \(a:[0,1]\to\mathbb{C}^{s\times t}\), is referred to as a sequence of (block) diagonal sampling matrices generated by \(a\).

Approximating classes of sequences. The notion of approximating classes of sequences (a.c.s.)is the cornerstone of an asymptotic approximation theory for sequences of matrices that has been developed since the last years of the XX century; see [3]. The formal definition of a.c.s.is reported in Definition 2. Throughout this paper, we denote by \(\|\cdot\|\) the Euclidean norm (\(2\)-norm) of vectors and the associated operator norm (induced norm) for matrices, which is defined by \(\|A\|=\max\{\|A{\boldsymbol{x}}\|:{\boldsymbol{x}}\in\mathbb{C}^n,\,\|{\boldsymbol{x}}\|=1\}\) for every \(A\in\mathbb{C}^{m\times n}\).

GLT sequences. In this section, we collect the properties of GLT sequences that we need in this paper. We begin with the following definition, which is an equivalent alternative to the usual (complicated) definition of GLT sequences. Throughout this paper, we denote by \(I_s\) the \(s\times s\) identity matrix.

The properties of GLT sequences that we need in this paper are listed below; for the corresponding proofs, see [6], [7]. Before being able to formulate these properties, we need to introduce some notation and terminology. Throughout this paper, we denote by \(O_{s,t}\) the \(s\times t\) zero matrix. The matrix \(O_{s,s}\) is denoted by \(O_s\). If the size is clear from the context, we often write \(O\) instead of \(O_{s,t}\) and \(O_s\). The conjugate transpose of a matrix \(A\) is denoted by \(A^*\). If \(A\) is a matrix, we denote by \(A^\dagger\) the Moore–Penrose pseudoinverse of \(A\). What is relevant for our purposes is that \(A^\dagger=A^{-1}\) whenever \(A\) is invertible. For more details on the pseudoinverse of a matrix, see [22], [23]. If \(A\) is a diagonalizable matrix and \(f\) is a complex-valued function defined at each point of \(\Lambda(A)\), we denote by \(f(A)\) the unique matrix such that \(f(A){\boldsymbol{v}}=f(\lambda){\boldsymbol{v}}\) whenever \(A{\boldsymbol{v}}=\lambda{\boldsymbol{v}}\). Note that a spectral decomposition of \(A\) immediately yields a spectral decomposition of \(f(A)\): \[\label{A-62f40A41} A=V\begin{bmatrix}\lambda_1 & & \\ & \ddots & \\ & & \lambda_n\end{bmatrix}V^{-1}\quad\implies\quad f(A)=V\begin{bmatrix}f(\lambda_1) & & \\ & \ddots & \\ & & f(\lambda_n)\end{bmatrix}V^{-1}.\tag{2}\] For more on matrix functions, we refer the reader to Higham’s book [24]. Hereafter, the composite function \(f\circ g\) is denoted by \(f(g)\).

1. Let \(\{X_n\}_n\) be a sequence of matrices, with \(X_n\) of size \(d_ns\times d_nt\) for some fixed positive integers \(s,t\) and some positive integer sequence \(\{d_n\}_n\) tending to \(\infty\), and let \(\kappa,\xi:[0,1]\times[-\pi,\pi]\to\mathbb{C}^{s\times t}\) be measurable.

   • If \(\{X_n\}_n\sim_{\rm GLT}\kappa\) and \(\kappa=\xi\) a.e., then \(\{X_n\}_n\sim_{\rm GLT}\xi\).

   • If \(\{X_n\}_n\sim_{\rm GLT}\kappa\) and \(\{X_n\}_n\sim_{\rm GLT}\xi\), then \(\kappa=\xi\) a.e.

2. Let \(\{X_n\}_n\) be a sequence of matrices, with \(X_n\) of size \(d_ns\times d_nt\) for some fixed positive integers \(s,t\) and some positive integer sequence \(\{d_n\}_n\) tending to \(\infty\), and let \(\kappa:[0,1]\times[-\pi,\pi]\to\mathbb{C}^{s\times t}\) be measurable.

   • If \(\{X_n\}_n\sim_{\rm GLT}\kappa\), then \(\{X_n\}_n\sim_\sigma\kappa\).

   • If \(\{X_n\}_n\sim_{\rm GLT}\kappa\) and the matrices \(X_n\) are Hermitian, then \(s=t\), \(\kappa\) is Hermitian a.e., and \(\{X_n\}_n\sim_\lambda\kappa\).

3. Let \(s,t\) be positive integers and let \(\{d_n\}_n\) be a positive integer sequence tending to \(\infty\). Then,

   • \(\{T_{d_n}(f)\}_n\sim_{\rm GLT}\kappa(x,\theta)=f(\theta)\) if \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) belongs to \(L^1([-\pi,\pi])\);

   • \(\{D_{d_n}(a)\}_n\sim_{\rm GLT}\kappa(x,\theta)=a(x)\) if \(a:[0,1]\to\mathbb{C}^{s\times t}\) is continuous a.e.;

   • for every sequence of matrices \(\{Z_n\}_n\) with \(Z_n\) of size \(d_ns\times d_nt\), we have \(\{Z_n\}_n\sim_{\rm GLT}\kappa(x,\theta)=O_{s,t}\) if and only if \(\{Z_n\}_n\sim_\sigma0\).

4. Let \(\{X_n\}_n\), \(\{Y_n\}_n\) be sequences of matrices, with \(X_n\) of size \(d_ns\times d_nt\) and \(Y_n\) of size \(d_nu\times d_nv\) for some fixed positive integers \(s,t,u,v\) and some positive integer sequence \(\{d_n\}_n\) tending to \(\infty\), and let \(\kappa:[0,1]\times[-\pi,\pi]\to\mathbb{C}^{s\times t}\) and \(\xi:[0,1]\times[-\pi,\pi]\to\mathbb{C}^{u\times v}\) be measurable. Suppose that \(\{X_n\}_n\sim_{\rm GLT}\kappa\) and \(\{Y_n\}_n\sim_{\rm GLT}\xi\). Then,

   • \(\{X_n^*\}_n\sim_{\rm GLT}\kappa^*\);

   • \(\{\alpha X_n+\beta Y_n\}_n\sim_{\rm GLT}\alpha\kappa+\beta\xi\) for every \(\alpha,\beta\in\mathbb{C}\) if \(\kappa\) and \(\xi\) are summable (i.e., \(s=u\) and \(t=v\));

   • \(\{X_nY_n\}_n\sim_{\rm GLT}\kappa\xi\) if \(\kappa\) and \(\xi\) are multipliable (i.e., \(t=u\));

   • \(\{X_n^\dagger\}_n\sim_{\rm GLT}\kappa^{-1}\) if \(\kappa\) is invertible a.e.;

   • \(\{f(X_n)\}_n\sim_{\rm GLT}f(\kappa)\) if the matrices \(X_n\) are Hermitian and \(f:\mathbb{C}\to\mathbb{C}\) is continuous.

3 GLT preconditioning↩︎

The purpose if this section is to prove Theorem 1, which is used in Section 4 to theoretically justify the efficiency of block Jacobi preconditioners for GLT sequences; see Remark 1. Theorem 1 is the generalization of [1] to the case where the involved GLT sequences have matrix-valued symbols. To prove it, we need the following lemma. Throughout this paper, we use the abbreviations HPD and HPSD for Hermitian positive definite and Hermitian positive semi-definite, respectively.

4 Block Jacobi and block Gauss–Seidel preconditioners for GLT sequences↩︎

This section is devoted to the proof of our main results (Theorems 2–3). Throughout this paper, if \(X\in\mathbb{C}^{m_1\times m_2}\) and \(Y\in\mathbb{C}^{\ell_1\times\ell_2}\), the tensor (Kronecker) product of \(X\) and \(Y\) is the \(m_1\ell_1\times m_2\ell_2\) matrix defined by \[X\otimes Y=[X_{ij}Y]_{i=1,\ldots,m_1}^{j=1,\ldots,m_2}=\begin{bmatrix}X_{11}Y & \cdots & X_{1m_2}Y\\ \vdots & & \vdots\\ X_{m_11}Y & \cdots & X_{m_1m_2}Y\end{bmatrix},\] and the direct sum of \(X\) and \(Y\) is the \((m_1+\ell_1)\times(m_2+\ell_2)\) matrix defined by \[X\oplus Y={\rm diag}(X,Y)=\begin{bmatrix}X & O\\ O & Y\end{bmatrix}.\] A partition of a positive integer \(n\) is a vector of positive integers \((n_1,\ldots,n_\nu)\) such that \(n_1+\ldots+n_\nu=n\).

We are now ready to state and prove Theorems 2–3. They show that any sequence of block Jacobi preconditioners or block Gauss–Seidel preconditioners for a GLT sequence with symbol \(\kappa\) is again a GLT sequence with symbol \(\kappa\), as long as the number of blocks in each preconditioner is fixed. Throughout this paper, the characteristic (indicator) function of a set \(E\) is denoted by \(\chi_E\).

5 GLH sequences and the extended Fasino–Tilli theorem↩︎

In this section, as an application of Theorem 2, we provide an alternative proof of the Fasino–Tilli theorem on the zero distribution of Hankel matrix sequences generated by \(L^1\) functions; see Theorem 4. This alternative proof naturally leads to both the definition of GLH sequences (Definition 5) and the extension of the Fasino–Tilli theorem to such sequences; see Theorem 5. We first recall that the \(n\)th (block) Hankel matrix generated by a function \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) belonging to \(L^1([-\pi,\pi])\) is the \(ns\times nt\) matrix denoted by \(H_n(f)\) and defined as follows: \[H_n(f)=[f_{i+j-1}]_{i,j=1}^n=\left[\begin{array}{ccccc}f_1 & f_2 & f_3 & \cdots & f_n\\ f_2 & f_3 & & \iddots & f_{n+1}\\ f_3 & & \iddots & & f_{n+2}\\ \vdots & \iddots & & & \vdots\\ f_n & \cdots & \cdots & \cdots & f_{2n-1}\end{array}\right],\] where \(\{f_k\}_{k\in\mathbb{Z}}\) are the Fourier coefficients of \(f\) defined in 1 . Any sequence of matrices of the form \(\{H_{d_n}(f)\}_n\), with \(d_n\to\infty\) and \(f:[-\pi,\pi]\to\mathbb{C}^{s\times t}\) in \(L^1([-\pi,\pi])\), is referred to as a (block) Hankel sequence generated by \(f\). The following theorem is due to Fasino and Tilli [25].

The proof of Theorem 4 that we are going to see is based on Theorem 2 and is different from the original one appeared in [25]. In addition to Theorem 2, the proof requires Lemma 3 [21]. Throughout this paper, we denote by \(W_n\) the \(n\times n\) anti-identity matrix (or flip matrix) defined by \[\label{flipm} W_n=\begin{bmatrix*} &&&&1\\ &&&1\\ &&\iddots\\ &1\\ 1 \end{bmatrix*}.\tag{14}\] A matrix \(P\in\mathbb{C}^{m\times n}\) is said to be a compression matrix if \(m\ge n\) and \(P^*P=I_n\).

Our proof of Theorem 4 amounts to applying Theorem 2 and Lemma 3, after observing that the bottom left block \(K_n(f)\) of the matrix \(T_{2n}(f)-\mathop{\rm blockdiag}_{n,n}^{s,t}(T_{2n}(f))\) coincides with the Hankel matrix \(H_n(f)\) multiplied by the permutation \(W_n\otimes I_t\); see 16 . In the case where the Toeplitz sequence \(\{T_{2n}(f)\}_n\) is replaced by a GLT sequence \(\{A_n\}_n\sim_{\rm GLT}\kappa\), with \(A_n\) of size \(2ns\times 2nt\) just like \(T_{2n}(f)\), the bottom left block \(K_n\) of the matrix \(A_n-\mathop{\rm blockdiag}_{n,n}^{s,t}(A_n)\) obviously coincides with some matrix \(H_n\) multiplied by the permutation \(W_n\otimes I_t\); it suffices to take \(H_n=K_n(W_n\otimes I_t)\). In this scenario, \(\{H_n\}_n\) is referred to as a GLH sequence.

Through a simple adaptation of the proof of Theorem 4, one can show that any GLH sequence is zero-distributed (and hence it is a GLT sequence with zero symbol by GLT2). We can therefore state the following result, which is an extension of Theorem 4.

6 Numerical experiments↩︎

In this section, we illustrate through numerical examples the efficiency of block Jacobi and block Gauss–Seidel preconditioners for GLT sequences, as predicted by Remark 1. In each example, we proceed as follows.

1. We fix a GLT sequence \(\{A_n\}_n\sim_{\rm GLT}\kappa\), where \(A_n\) is an invertible square matrix of size \(ns\) for some fixed positive integer \(s\) and \(\kappa:[0,1]\times[-\pi,\pi]\to\mathbb{C}^{s\times s}\) is a measurable function.

2. For every \(\nu\ge2\) and every \(n\ge\nu\), we define \((n_1,\ldots,n_\nu)\) as the partition of \(n\) given by \[n_i=\lfloor n/\nu\rfloor,\qquad i=1,\ldots,\nu-1,\qquad n_\nu=n-(\nu-1)\lfloor n/\nu\rfloor,\] and we consider the corresponding block Jacobi and block Gauss–Seidel preconditioners for \(A_n\), denoted by \(P_n^J=\mathop{\rm blockdiag}_{n_1,\ldots,n_\nu}^{s,s}(A_n)\) and \(P_n^{GS}=\mathop{\rm blocktril}_{n_1,\ldots,n_\nu}^{s,s}(A_n)\), respectively. Both \(P_n^J\) and \(P_n^{GS}\) will turn out to be invertible just like \(A_n\).

3. In the case where the matrices \(A_n\) are HPD, the matrices \(P_n^J\) are HPD as well and we report in a table, for increasing values of \(\nu\) and \(n\), the number of iterations needed by the preconditioned conjugate gradient (PCG) method with preconditioner \(P_n^J\) for solving the linear system \(A_n{\boldsymbol{x}}=\mathbf{1}\) up to a precision of \(10^{-6}\), where \(\mathbf{1}\) is the vector of all ones. The PCG method is started with the zero vector \(\mathbf{0}\). We also report in another table, for the same values of \(\nu\) and \(n\) considered in the previous table, the number of iterations needed by the preconditioned generalized minimal residual (PGMRES) method with preconditioner \(P_n^{GS}\) for solving the linear system \(A_n{\boldsymbol{x}}=\mathbf{1}\) up to a precision of \(10^{-6}\). The PGMRES method is started with the zero vector \(\mathbf{0}\) and applied without restarting.

4. In the case where the matrices \(A_n\) are not HPD, we report in a table, for increasing values of \(\nu\) and \(n\), the number of iterations needed by the PGMRES method with preconditioner \(P_n^J\) for solving the linear system \(A_n{\boldsymbol{x}}=\mathbf{1}\) up to a precision of \(10^{-6}\). We also report in another table, for the same values of \(\nu\) and \(n\) considered in the previous table, the number of iterations needed by the PGMRES method with preconditioner \(P_n^{GS}\) for solving the linear system \(A_n{\boldsymbol{x}}=\mathbf{1}\) up to a precision of \(10^{-6}\). As in item (C), the PGMRES method is started with the zero vector \(\mathbf{0}\) and applied without restarting.

In all the considered examples, for each fixed \(\nu\), the number of PCG/PGMRES iterations is bounded from above by a constant \(C_\nu\) independent of \(n\). This shows an optimal convergence rate of both methods and, consequently, the efficiency of the preconditioners \(P_n^J\) and \(P_n^{GS}\) for \(A_n\).

\footnotesize
\centering
\caption{Example~\ref{e1} --- Number of PCG iterations with preconditioner $P_n^J$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$.}
\label{e1T}
\begin{tabular}{rrrrrrrrr}
\toprule
 & \multicolumn{8}{c}{PCG iterations with preconditioner $P_n^J$ $_{\vphantom{\int}}$}\\
\cline{2-9}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=3$ & $\hphantom{100}\nu=4$ & $\hphantom{100}\nu=5$ & $\hphantom{100}\nu=6$ & $\hphantom{100}\nu=7$ & $\hphantom{100}\nu=8$ & $\hphantom{100}\nu=9^{\vphantom{\int^0}}$\\
\midrule
40 & 5 & 9 & 8 & 9 & 14 & 15 & 12 & 16\\
80 & 6 & 10 & 9 & 10 & 15 & 16 & 13 & 18\\
160 & 6 & 10 & 10 & 11 & 16 & 18 & 14 & 19\\
320 & 6 & 11 & 10 & 12 & 17 & 19 & 16 & 21\\
640 & 7 & 12 & 11 & 12 & 18 & 19 & 17 & 20\\
1280 & 7 & 12 & 12 & 13 & 19 & 20 & 18 & 22\\
2560 & 7 & 13 & 12 & 14 & 20 & 21 & 19 & 23\\
\bottomrule
\end{tabular}
\vspace{20pt}
\caption{Example~\ref{e1} --- Number of PGMRES iterations with preconditioner $P_n^{GS}$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$.}
\label{e1T39}
\begin{tabular}{rrrrrrrrr}
\toprule
 & \multicolumn{8}{c}{PGMRES iterations with  preconditioner $P_n^{GS}$ $_{\vphantom{\int}}$}\\
\cline{2-9}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=3$ & $\hphantom{100}\nu=4$ & $\hphantom{100}\nu=5$ & $\hphantom{100}\nu=6$ & $\hphantom{100}\nu=7$ & $\hphantom{100}\nu=8$ & $\hphantom{100}\nu=9^{\vphantom{\int^0}}$\\
\midrule
40 & 4 & 6 & 7 & 8 & 9 & 9 & 10 & 10\\
80 & 4 & 6 & 7 & 9 & 10 & 10 & 11 & 11\\
160 & 4 & 6 & 8 & 9 & 11 & 11 & 12 & 13\\
320 & 5 & 7 & 8 & 10 & 11 & 12 & 13 & 13\\
640 & 5 & 7 & 9 & 10 & 12 & 12 & 13 & 14\\
1280 & 5 & 7 & 9 & 11 & 12 & 13 & 13 & 15\\
2560 & 5 & 8 & 9 & 11 & 12 & 13 & 14 & 15\\
\bottomrule
\end{tabular}

\footnotesize
\centering
\caption{Example~\ref{e2} --- Number of PGMRES iterations with preconditioner $P_n^J$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$.}
\label{e2T}
\begin{tabular}{rrrrr}
\toprule
 & \multicolumn{4}{c}{PGMRES iterations with preconditioner $P_n^J$ $_{\vphantom{\int}}$}\\
\cline{2-5}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=4$ & $\hphantom{100}\nu=8$ & $\hphantom{100}\nu=16^{\vphantom{\int^0}}$\\
\midrule
40 & 4 & 10 & 22 & 40\\
80 & 4 & 10 & 22 & 43\\
160 & 4 & 10 & 22 & 43\\
320 & 4 & 10 & 22 & 44\\
640 & 4 & 8 & 19 & 44\\
1280 & 4 & 8 & 16 & 35\\
2560 & 4 & 7 & 13 & 27\\
\bottomrule
\end{tabular}
\vspace{20pt}
\centering
\caption{Example~\ref{e2} --- Number of PGMRES iterations with preconditioner $P_n^{GS}$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$.}
\label{e2T39}
\begin{tabular}{rrrrr}
\toprule
 & \multicolumn{4}{c}{PGMRES iterations with preconditioner $P_n^{GS}$ $_{\vphantom{\int}}$}\\
\cline{2-5}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=4$ & $\hphantom{100}\nu=8$ & $\hphantom{100}\nu=16^{\vphantom{\int^0}}$\\
\midrule
40 & 2 & 4 & 8 & 17\\
80 & 2 & 4 & 9 & 17\\
160 & 2 & 4 & 9 & 17\\
320 & 2 & 4 & 9 & 18\\
640 & 2 & 5 & 9 & 18\\
1280 & 2 & 5 & 9 & 19\\
2560 & 2 & 5 & 10 & 19\\
\bottomrule
\end{tabular}

\footnotesize
\centering
\caption{Example~\ref{e3} --- Number of PCG iterations with preconditioner $P_n^J$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$ in the case where $a(x)=1+\sqrt x$.}
\label{e3T}
\begin{tabular}{rrrrr}
\toprule
 & \multicolumn{4}{c}{PCG iterations with preconditioner $P_n^J$ $_{\vphantom{\int}}$}\\
\cline{2-5}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=7$ & $\hphantom{100}\nu=12$ & $\hphantom{100}\nu=17^{\vphantom{\int^0}}$\\
\midrule
40 & 5 & 7 & 8 & 9\\
80 & 4 & 7 & 7 & 8\\
160 & 5 & 7 & 6 & 6\\
320 & 5 & 7 & 6 & 6\\
640 & 5 & 7 & 6 & 6\\
1280 & 5 & 7 & 6 & 6\\
2560 & 5 & 7 & 6 & 6\\
\bottomrule
\end{tabular}
\vspace{20pt}
\centering
\caption{Example~\ref{e3} --- Number of PGMRES iterations with preconditioner $P_n^{GS}$ for solving $A_n\xx=\mathbf1$ up to a precision of $10^{-6}$ in the case where $a(x)=1+\sqrt x$.}
\label{e3T39}
\begin{tabular}{rrrrr}
\toprule
 & \multicolumn{4}{c}{PGMRES iterations with preconditioner $P_n^{GS}$ $_{\vphantom{\int}}$}\\
\cline{2-5}
$n$ & $\hphantom{100}\nu=2$ & $\hphantom{100}\nu=7$ & $\hphantom{100}\nu=12$ & $\hphantom{100}\nu=17^{\vphantom{\int^0}}$\\
\midrule
40 & 3 & 6 & 6 & 7\\
80 & 3 & 6 & 6 & 6\\
160 & 3 & 6 & 6 & 6\\
320 & 3 & 6 & 6 & 6\\
640 & 3 & 6 & 6 & 6\\
1280 & 3 & 5 & 6 & 6\\
2560 & 3 & 5 & 6 & 6\\
\bottomrule
\end{tabular}

7 Conclusions and perspectives↩︎

In our main results (Theorems 2–3), we have proved that any sequence of block Jacobi or block Gauss–Seidel preconditioners \(\{P_n\}_n\) associated with a GLT sequence \(\{A_n\}_n\sim_{\rm GLT}\kappa\) is again a GLT sequence with the same symbol \(\kappa\), as long as the number of blocks appearing in each preconditioner \(P_n\) is a fixed number \(\nu\) independent of \(n\). As a consequence of this result:

1. in Remark 1, we have predicted the efficiency of block Jacobi and block Gauss–Seidel preconditioners for GLT sequences; an efficiency that has been illustrated through numerical experiments in Section 6;

2. in Theorem 5, we have extended the Fasino–Tilli theorem on the zero distribution of Hankel sequences generated by \(L^1\) functions to the larger class of GLH sequences defined in Definition 5.

We conclude this paper by suggesting some possible future lines of research.

• In Theorems 2–3, the considered block Jacobi and block Gauss–Seidel preconditioners have a fixed number \(\nu\) of blocks independent of \(n\). From a theoretical point of view, it could be interesting to investigate whether Theorems 2–3 continue to hold if \(\nu=\nu_n\) depends on \(n\) and satisfies some additional assumptions, such as \(\nu=\nu_n=o(d_n)\) as \(n\to\infty\).

• The block Jacobi preconditioner \(P_n^J\) considered in Theorem 2 for a matrix \(A_n\) belonging to a GLT sequence \(\{A_n\}_n\sim_{\rm GLT}\kappa\) is made up of blocks \(A_n^{(1)},\ldots,A_n^{(\nu)}\) that are diagonal submatrices of \(A_n\). In general, these submatrices are not computationally easy, i.e., solving a linear system associated with them can be computationally expensive. A possible way to improve the efficiency of the preconditioner \(P_n^J\) is to replace \(A_n^{(1)},\ldots,A_n^{(\nu)}\) with computationally easier “GLT blocks” \(\tilde{A}_n^{(1)},\ldots,\tilde{A}_n^{(\nu)}\) such that the new preconditioner \(\tilde{P}_n^J=\tilde{A}_n^{(1)}\oplus\cdots\oplus\tilde{A}_n^{(\nu)}\) satisfies \(\{\tilde{P}_n^J\}_n\sim_{\rm GLT}\tilde{\kappa}\) with either \(\tilde{\kappa}\approx\kappa\) or, even better, \(\tilde{\kappa}=\kappa\). To achieve the desired goal, the “GLT blocks” \(\tilde{A}_n^{(1)},\ldots,\tilde{A}_n^{(\nu)}\) could be defined as appropriate products between diagonal sampling matrices and computationally easy Toeplitz-like matrices such as circulant matrices [26] or \(\tau\) matrices [27].

• In view of possible multidimensional applications of block Jacobi/Gauss–Seidel preconditioning to GLT sequences, an interesting topic for future research is to extend Theorems 2–3 to multilevel GLT sequences, which arise, for instance, in the numerical discretization of partial differential equations [4]–[6], [8].

• In the numerical discretization of differential problems by domain decomposition methods (DDMs), two popular iterative solution methods are the additive and multiplicative Schwarz methods; see [28], [29]. What is relevant to us is that, roughly speaking, the additive Schwarz method is another name for the block Jacobi method and the multiplicative Schwarz method is another name for the block Gauss–Seidel method. Thus, DDMs are a research field where the results of the present paper, along with their multidimensional extensions mentioned in the previous item, could find application. More specifically, Theorems 2–3 and their multidimensional extensions could be used in combination with the basic theoretical tools provided in [30] to analyze the additive and multiplicative Schwarz methods applied to GLT structures arising from DDM discretizations. A first contribution in this direction has recently appeared in [31].

• A successful idea proposed by Pestana and Wathen [32] for solving a real non-symmetric Toeplitz linear system \(T_n(f){\boldsymbol{x}}=\boldsymbol{b}\) consists in pre-multiplying both sides by the flip matrix 14 and solving the resulting real symmetric linear system \(W_nT_n(f){\boldsymbol{x}}=W_n\boldsymbol{b}\) by either the minimal residual (MINRES) method or its preconditioned version. This idea, which we call “flipping strategy”, has gained more and more popularity over time and it has now become a consolidated research area. What is relevant to us is that the spectral distribution of the sequence of system matrices \(W_nT_n(f)\), which plays a significant role in the convergence analysis of MINRES, was established in [33], [34] by leveraging on the Fasino–Tilli theorem; see also the generalization to the multilevel case performed in [35], [36], which again invokes (the multilevel version of) the Fasino–Tilli theorem. Moreover, in the spectral distribution analysis of block matrices with rectangular Toeplitz blocks—which arise in several applications as explained in [21]—the Fasino–Tilli theorem is again a cornerstone; see [21] for more details. Given the central role played by the Fasino–Tilli theorem in the asymptotic spectral analysis of Toeplitz-like matrices such as those mentioned above, we can predict that the same central role is played by the extended Fasino–Tilli theorem (Theorem 5) in the asymptotic spectral analysis of “GLT-like” matrices. This observation opens the way to several future investigations on the possible applications of the extended Fasino–Tilli theorem, including in particular the spectral distribution analysis of: (a) flipped “GLT matrices” resulting from the application of the flipping strategy to “GLT linear systems” \(A_n{\boldsymbol{x}}=\boldsymbol{b}\) with \(\{A_n\}_n\sim_{\rm GLT}\kappa\); (b) block matrices with rectangular “GLT blocks”.

Acknowledgements↩︎

The authors are members of the research group GNCS (Gruppo Nazionale per il Calcolo Scientifico) of INdAM (Istituto Nazionale di Alta Matematica). This work was supported by the INdAM-GNCS project PASTRAMI (sPline And Solver innovaTions foR Adaptive isogeoMetric analysIs, CUP E53C24001950001), by the MUR excellence department project MatMod@TOV awarded to the Department of Mathematics of the University of Rome Tor Vergata (CUP E83C23000330006), by the Department of Mathematics of the University of Rome Tor Vergata through the project METRO (Methods and modEls for arTificial neuRal netwOrks, CUP E83C25000630005), by the Theory, Economics and Systems Laboratory (TESLAB) of the Department of Computer Science of the Athens University of Economics and Business, by the PRIN-PNRR project MATHPROCULT (MATHematical tools for predictive maintenance and PROtection of CULTural heritage, Code P20228HZWR, CUP J53D23003780006), and by the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701. The JU receives support from the European Union’s Horizon 2020 research and innovation programme and Belgium, France, Germany, Switzerland.

References↩︎