A note on the spectral analysis of matrix-sequences via GLT momentary symbols: from all-at-once solution of parabolic problems to distributed fractional order matrices

Throughout this paper, we will use the following notations. Let ${f}:G\to\mathbb{C}$ be a function belonging to $L^{1}(G)$, with $G\subseteq\mathbb{R}^{\ell}$, $\ell\geq 1$, a measurable set. We denote by $\{A_{n}\}_{n}$ the matrix-sequence, whose elements are given by the matrices $A_{n}$ of dimension $n\times n$. Let $s,d\in\mathbb{N}$. Let $\mathbf{n}=(n_{1},n_{2},\dots,n_{d})$ be a multi-index, we indicate by $\{A_{\mathbf{n}}\}_{\mathbf{n}}$, the $d$-level $s\times s$ block matrix-sequence, whose elements are the matrices $A_{\mathbf{n}}$ of size $d=d(\mathbf{n},s)=sn_{1}n_{2}\cdots n_{d}$.

In the following we report the main background concerning the concepts of Toeplitz and circulant matrices, for simplicity, in the scalar unilevel setting. We only provide the generalization in the block multilevel case of the results that will be exploited for the purpose of the paper.

In the following we focus on the two important sub-classes given by the Toeplitz matrices $T_{n}(f)\in\mathbb{C}^{n\times n}$ and the circulant matrices $C_{n}(f)\in\mathbb{C}^{n\times n}$, associated with a function $f$, called the generating function.

Note that, from the definition follows that if ${f}$ is a trigonometric polynomial of fixed degree less than $n$, the entries of $D_{n}(f)$ are the eigenvalues of $C_{n}(f)$, explicitly given by sampling the generating function $f$ using the grid $\theta_{j,n}^{c}$.

The type of domain (either one-dimensional $[-\pi,\pi]$ or d-dimensional $[-\pi,\pi]^{d}$) and codomain (either the complex field or the space of $s\times s$ complex matrices) of $f$ gives rise to different kinds of Toeplitz matrices, see Table 1 for a complete overview.

In particular, we provide the definition of a $d$-level $s\times s$ block Toeplitz matrices $T_{\mathbf{n}}(\mathbf{f})$ starting from $d$-variate matrix-valued function $\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s\times s}$, with with $\mathbf{f}\in L^{1}([-\pi,\pi]^{d})$.

The following result provides an important relation between tensor products and multilevel Toeplitz matrices.

In this subsection we introduce the definition of asymptotic distribution in the sense of the eigenvalues and of the singular values, first for a generic matrix-sequence $\{A_{n}\}_{n}$, and then we report specific results concerning the distributions of Toeplitz and circulant matrix-sequences. Finally, we recall the notion of GLT algebra and we introduce a general notion of GLT momentary symbols. We remind that in a more specific and limited setting, the notion of momentary symbols is given in [6]: here we generalize the definition in [6].

In the case where the matrix-sequence is a Toeplitz matrix-sequence generated by a function, the singular value distribution and the spectral distribution have been well studied in the past few decades. At the beginning Szegő in [22] showed that the eigenvalues of the Toeplitz matrix $T_{n}(f)$ generated by real-valued $f\in L^{\infty}([-\pi,\pi])$ are asymptotically distributed as $f$. Moreover, under the same assumption on $f$, Avram and Parter [1, 25] proved that the singular values of $T_{n}(f)$ are distributed as $|f|$. This result has been undergone many generalizations and extensions among the years (see [20, 21, 3, 2] and the references therein).

The generalized Szegő theorem that describes the singular value and spectral distribution of Toeplitz sequences generated by a scalar $f\in L^{1}([-\pi,\pi])$ is given as follows [31].

Tilli [29] generalized the proof to the block-Toeplitz setting and we report the extension of the eigenvalue result to the case of multivariate Hermitian matrix-valued generating functions.

Concerning the circulant matrix-sequences, though the eigenvalues of a $C_{\bf n}(\textbf{f})$ are explicitly known, a result like Theorem 1.9 and Theorem 1.10 does not hold for sequences $\left\{C_{\bf n}(\textbf{f})\right\}_{{\bf n}}$ in general. Indeed, the Fourier sum of f converges to f under quite restrictive assumptions (see [32]). In particular, if f belongs to the Dini-Lipschitz class, then $\left\{C_{\bf n}(\textbf{f})\right\}_{{\bf n}}\sim_{\lambda}\textbf{f}$, (see [19] for more relationships between circulant sequences and spectral distribution results).

A part from the circulant algebra, introduced in Section 1.2, we recall that other particular matrix algebras have interesting properties and can be exploited for our purpose. In particular, we mention the well-known $\tau$-algebras, see [7] and references therein. Here, we restrict the analysis to the case of the matrix algebras ${\tau_{\varepsilon,\varphi}}$, introduced in [7], where an element of the algebra is a matrix

We can associate to this matrix a function $g$ of the form $g(\theta)=a+2b\cos\theta$. For some values of $\varepsilon$ and $\varphi$ the exact eigenvalues of $T_{n,\epsilon,\varphi}(g)$ are given by sampling with specific grids; for detailed examples see [6] and [18] for asymptotic results.

In Table 2 we provide the proper grids $\theta^{(\varepsilon,\varphi)}_{j,n}$ and $\Theta_{i,j,n}^{(\varepsilon,\varphi)}$ to give the exact eigenvalues and eigenvectors respectively for $\varepsilon,\varphi\in\{-1,0,1\}$.

Since all grids $\theta_{j,n}^{(\varepsilon,\varphi)}$ associated with $\tau_{\varepsilon,\varphi}$-algebras where $\varepsilon,\varphi\in\{-1,0,1\}$ are uniformly spaced grids, we know that

In this subsection we will introduce the main properties from the theory of Generalized Locally Toeplitz (GLT) sequences and the practical features, which are sufficient for our purposes, see [20, 21, 3, 2].

In particular, we consider the multilevel and block setting with $d$ being the number of levels.

GLT1:

Each GLT sequence has a singular value symbol $\mathbf{f}(\boldsymbol{\theta,x})$ which is measurable according to the Lebesgue measure and according to the second item in Definition 1.7 with $\ell=2d$. In addition, if the sequence is Hermitian, then the distribution also holds in the eigenvalue sense.

We specify that a GLT sequence $\{A_{\mathbf{n}}\}_{\mathbf{n}}$ has GLT symbol $\mathbf{f}(\boldsymbol{\theta,x})$ writing $\{A_{\mathbf{n}}\}_{\mathbf{n}}\sim_{\textsc{glt}}\mathbf{f}(\boldsymbol{\theta,x})$, $(\boldsymbol{\theta,x})\in[-\pi,\pi]^{d}\times[0,1]^{d}$.

GLT2:

The set of GLT sequences form a $*$-algebra, i.e., it is closed under linear combinations, products, inversion (whenever the symbol is singular, at most, in a set of zero Lebesgue measure), and conjugation. Hence, we obtain the GLT symbol of algebraic operations of a finite set of GLT sequences by performing the same algebraic manipulations of the symbols of the considered GLT sequences.

GLT3:

Every Toeplitz sequence $\{T_{\mathbf{n}}(\mathbf{f})\}_{n}$ generated by a function $\mathbf{f}(\boldsymbol{\theta})$ belonging to $L^{1}([-\pi,\pi]^{d})$ is a GLT sequence and with GLT symbol given by $\mathbf{f}$. Every diagonal sampling sequence $\{D_{\mathbf{n}}(\mathbf{a})\}_{n}$ generated by a Riemann integrable function $\mathbf{a}(\boldsymbol{x})$, $\boldsymbol{x}\in[0,1]^{d}$ is a GLT sequence and with GLT symbol given by $\mathbf{a}$.

GLT4:

Every sequence which is distributed as the constant zero in the singular value sense is a GLT sequence with symbol $0$. In particular:

• •:

  every sequence in which the rank divided by the size tends to zero, as the matrix-size tends to infinity;

• •:

  every sequence in which the trace-norm (i.e., sum of the singular values) divided by the size tends to zero, as the matrix-size tends to infinity.

From a practical view-point, on the one hand, one of the main advantages for a sequence of belonging to the GLT class is that, under certain hypotheses, crucial spectral and singular value information can be derived using the concept of GLT symbol. On the other hand, the above properties imply the following important features of the GLT symbol. Given a sequence $\{A_{\mathbf{n}}\}_{\mathbf{n}}$ obtained by algebraic operations of a finite set of GLT sequences, the small-norm and low-rank terms which composes the sequence should be neglected in the computation of the GLT symbol. Consequently, it happens that for small matrix-sizes $n$, the approximations may not be as accurate as it is desirable.

For this reason in [6] it has been introduced and exploited the concept of a (singular value and spectral) “momentary symbols”, starting from a special case of Toeplitz structures. Here we generalize the notion to that of “GLT momentary symbols”: the construction stems from that of the symbol in the GLT sense, but in practice the information of the small norm contributions is kept in the symbol and this may lead to higher accuracy, at least in some emblematic cases, when approximating the singular values and eigenvalues of Toeplitz-like matrices, even for small dimensions.

For clarity in this subsection we consider the matrix-sequences in detail only in the unilevel and scalar setting. We want to avoid a cumbersome notation, but the ideas are extensible in a plain manner to the case where the involved GLT symbols are also matrix-valued and multivariate, as briefly sketched.

As an example, we take the following second-order differential equation with Dirichlet boundary conditions

The well-posedness of the previous diffusion-convection-advection problem holds in the case where $a(x)\in C^{1}(0,1)$. Furthermore, the uniqueness and existence of the solution are guaranteed in the case where $a(x)>0$, $c(x)\geq 0$ and with continuous functions $b(x),c(x)$ on $[0,1]$, with $f(x)\in L^{2}([0,1])$ (see [8]). For a more exhaustive discussion regarding the conditions of existence and uniqueness, even in the multidimensional case, we refer to [27, 28] and references therein.

From a GLT viewpoint, we only require the following much weaker assumptions that are

• •

  $a(x),c(x)$ are real-valued functions, continuous almost everywhere, defined in $[0,1]$,

• •

  $b(x)$ is a real-valued function on $[0,1]$, such that $|b(x)x^{\alpha}|$ is bounded for some $\alpha<3/2$,

while $f(x)$ is a general function.

We employ central second-order finite differences for approximating the given equation. We define the stepsize $h=\frac{1}{n+1}$ and the points $x_{k}=kh$ for $k$ belonging to the interval $[0,n+1]$. Let $a_{k}:=a(x_{\frac{k}{2}})$ for any $k\in[0,2n+2]$ and set $b_{j}:=b(x_{j})$, $c_{j}:=c(x_{j})$, $f_{j}:=f(x_{j})$ for every $j=0,\ldots,n+1$. We compute approximations $u_{j}$ of the values $u(x_{j})$ for $j=1,\ldots,n$ by solving the following linear system

where

In the case where $a(x)\equiv 1$ and $b(x)\equiv 1$, we find the basic Toeplitz structures

which are of importance since $A_{n}=D_{n}(a)K_{n}+E_{n}$ and $B_{n}={\color[rgb]{0,0,0}h}D_{n}(b)H_{n}$ with $\{E_{n}\}_{n}\sim_{\sigma}0$, as an immediate check in [4] can show. Therefore, by using the GLT axioms (as done in detail in [4]) we obtain

As a conclusion $\left\{B_{n}\right\}_{n}\sim_{\textsc{glt}}\text{0}$, $\ \left\{C_{n}\right\}_{n}\sim_{\textsc{glt}}0$ and hence, setting $X_{n}=A_{n}+B_{n}+C_{n}$ the actual coefficient matrix of the linear system in (1.11), again by the $*$-algebra structure of the GLT matrix-sequences, we deduce

Now, following [6], the idea is to consider not only the asymptotic setting, but also the case of moderate sizes. As a consequence, for increasing the precision of the evaluation of eigenvalues and singular values, we can associate to

the specific symbol $f_{n}(x,\theta)=a(x)(2-2\cos\theta)+h\mathbf{i}b(x)\sin\theta+h^{2}c(x)$.

We are now in position to give a formal definition of GLT momentary symbols.

Of course, in line with Section 1.5, the momentary symbol could be matrix-valued with a number of variables equal to $2d$ and domain $[-\pi,\pi]^{d}\times[0,1]^{d}$ if the basic matrix-sequences appearing in Definition 1.11 are, up to proper scaling, matrix-valued and multilevel GLT matrix-sequences. For example in the scalar $d$-variate setting relation (1.13) takes the form

which is a plain multivariate (possibly block) version of (1.13).

Clearly there is a link with the GLT theory stated in the next result.

The given definition of momentary symbols is inspired, as it is clear from the initial example of diffusion-convection-advection equation, by the example of approximated differential equations, where the presence of differential operators of different orders induces, after a possible proper scaling, a structure like that reported in (1.12).

The idea is that the momentary symbol can be used for giving a more precise evaluation either of the spectrum or of the eigenvalues for moderate sizes of the matrices and not only asymptotically. However, we should be aware that, intrinsically, there is no general recipe especially for the eigenvalues. In fact, as already proven in [26], a rank one perturbation of infinitesimal spectral norm actually can change the spectra of matrix-sequences, sharing the same GLT symbol and even sharing the same sequence of momentary symbols.

Example 1:

Take the matrices $T_{n}(e^{\mathbf{i}\theta})$ and $X_{n}=T_{n}(e^{\mathbf{i}\theta})+e_{1}e_{n}^{T}c_{n}^{(1)}$ with $c_{n}^{(1)}=n^{-\alpha}$, $\alpha>0$ any positive number independent of the matrix-size $n$. By direct inspection $\{e_{1}e_{n}^{T}c_{n}^{(1)}\}_{n}\sim_{\sigma}0$ and hence it is a GLT matrix-sequence with zero symbol, independently of the parameter $\alpha$. If we look at the GLT momentary symbols then they coincide with the GLT symbol for both $\{T_{n}(e^{\mathbf{i}\theta})\}_{n}$ and $\{X_{n}\}_{n}$: however while in the first case, the eigenvalues are all equal to zero, in the second case they distribute asymptotically as the GLT symbol $e^{\mathbf{i}\theta}$ (which is also the GLT momentary symbol for any $n$).

Example 2:

Take a positive function $a$ defined on $[0,1]$ and the matrices $D_{n}(a)T_{n}(e^{\mathbf{i}\theta})$ and $X_{n}=D_{n}(a)T_{n}(e^{\mathbf{i}\theta})+e_{1}e_{n}^{T}c_{n}^{(1)}$ with $c_{n}^{(1)}=n^{-\alpha}$, $\alpha>0$ any positive number independent of the matrix-size $n$. Since $\{e_{1}e_{n}^{T}c_{n}^{(1)}\}_{n}$ is a GLT matrix-sequence with zero symbol, independently of the parameter $\alpha$, we deduce that both $\{D_{n}(a)T_{n}(e^{\mathbf{i}\theta})\}_{n}$ and $\{X_{n}\}_{n}$ share the same GLT symbol $a(x)e^{\mathbf{i}\theta}$ (which is also the momentary symbol for any $n$). Again there is dramatic change: while in the first case, the eigenvalues are all equal to zero, in the second case they distribute asymptotically as the function $\hat{a}e^{\mathbf{i}\theta}$, where $\hat{a}$ is the limit (if it exists) of the geometric mean of sampling values present in $D_{n}(a)$, as $n$ tends to infinity: since $n^{-\alpha/n}$ converges to $1$ independently of the parameter $\alpha$ as $n$ tends to infinity, $\hat{a}$ will depend only on the diagonal values of $D_{n}(a)$. As a conclusion the eigenvalue distributions do not coincide with the GLT momentary symbols and this is a message that the present tool could be not effective and even misleading, when very non-normal matrices are considered.

In this setting it must be emphasized that the asymptotic eigenvalue distribution is discontinuous with respect to the standard norms or metrics widely considered in the context of matrix-sequences.

The principal ingredients of the time discretization are:

• •

  Choosing $N_{t}$ equispaced points in $[0,T]$ with stepsize $h_{t}=T/N_{t}$, that is, $t_{j}=jh_{t}$, for $j=1,\ldots,N_{t}$.

• •

  Discretizing in time by standard Euler backwards.

Regarding notations, for the sake of simplicity, since we are considering a 2D problem, the symbols will have as Fourier variable $(\theta,\xi)$ instead of the standard choice $(\theta_{1},\theta_{2})$ indicated in the notations of Section 1.2 (see Definition 1.4).

The resulting matrix is $J_{N_{t}}$, which has the following unilevel scalar Toeplitz structure:

where $f_{J}$ is the generating function of the matrix-sequence $\{T_{n}(f_{J})\}_{n}$ with

The principal elements of the time discretization are:

• •

  Choosing $N_{x}$ equispaced points in $[a,b]$. Since we are considering periodic boundary conditions, we have step size $h_{x}=(a-b)/N_{x}$ and $x_{j}=h_{x}(j-1)$, for $j=1,\ldots,N_{x}$.

• •

  Discretizing in space using second order finite differences.

Consequently, the space discretization matrix will be the circulant matrix $Q_{N_{x}}$ of the form: