A matrix-less method to approximate the spectrum and the spectral function of Toeplitz matrices with real eigenvalues

Abstract

It is known that the generating function f of a sequence of Toeplitz matrices {T_n(f)}_n may not describe the asymptotic distribution of the eigenvalues of T_n(f) if f is not real. In this paper, we assume as a working hypothesis that, if the eigenvalues of T_n(f) are real for all n, then they admit an asymptotic expansion of the same type as considered in previous works, where the first function, called the eigenvalue symbol $\mathfrak {f}$, appearing in this expansion is real and describes the asymptotic distribution of the eigenvalues of T_n(f). This eigenvalue symbol $\mathfrak {f}$ is in general not known in closed form. After validating this working hypothesis through a number of numerical experiments, we propose a matrix-less algorithm in order to approximate the eigenvalue distribution function $\mathfrak {f}$. The proposed algorithm, which opposed to previous versions, does not need any information about neither f nor $\mathfrak {f}$ is tested on a wide range of numerical examples; in some cases, we are even able to find the analytical expression of $\mathfrak {f}$. Future research directions are outlined at the end of the paper.

1 Introduction

Given a function f ∈ L¹([−π, π]), the n × n Toeplitz matrix generated by f is defined as

$$ \begin{array}{@{}rcl@{}} T_{n}(f)=\left[\hat{f}_{i-j}\right]_{i,j=1}^{n}=\left[ \begin{array}{ccccccccc} \hat{f}_{0}&\hat{f}_{-1}&\cdots&\hat{f}_{1-n}\vphantom{\vdots}\\ \hat{f}_{1}&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\hat{f}_{-1}\\ \hat{f}_{n-1}&\cdots&\hat{f}_{1}&\hat{f}_{0}\vphantom{\vdots}\\ \end{array} \right], \end{array} $$

where the numbers $\hat {f}_k$ are the Fourier coefficients of f, that is,

$$ \hat{f}_{k}=\frac{1}{2\pi}{\int}_{-\pi}^{\pi}f(\theta)e^{-\mathbf{i}k\theta}\mathrm{d}\theta,\qquad f(\theta)=\sum\limits_{k=-\infty}^{\infty}\hat{f}_{k}e^{\mathbf{i}k\theta}, \qquad \mathbf{i}^{2}=-1,\qquad k\in\mathbb{Z}. $$

(1)

It is known that the generating function f, also known as the symbol of the sequence {T_n(f)}_n, describes the asymptotic distribution of the singular values of T_n(f); if f is real or if $f\in L^{\infty }([-\pi ,\pi ])$ and its essential range has empty interior and does not disconnect the complex plane, then f also describes the asymptotic distribution of the eigenvalues of f; see [12, 20, 25] for details and [20, Section 3.1] for the notion of asymptotic singular value and eigenvalue distribution of a sequence of matrices. We will write $\{T_n(f)\}_n\sim _{\lambda } f$ to indicate that {T_n(f)}_n has an asymptotic eigenvalue distribution described by f. The precise definition is the following [20, Section 3.1].

Definition 1 (eigenvalue distribution of a matrix sequence)

Let {A_n}_n be a matrix sequence, with A_n of size n, and let $f: D\subset \mathbb {R}^k\to \mathbb {C}$ be a measurable function defined on a set D with $0<\mu _k(D)<\infty$. We say that it has an asymptotic eigenvalue (or spectral) distribution described by f and we write $\{A_n\}_n\sim _{\lambda } f$, if

$$ \underset{n\to\infty}{\lim}\frac{1}{n}\sum\limits_{j=1}^{n} F(\lambda_{j}(A_{n})=\frac{1}{\mu_{k}(D)}{\int}_{D} F(f(x_{1},\ldots,x_{k}))\mathrm{d}x_{1}\cdots\mathrm{d}x_{k}, \quad \forall F\in C_{c}(\mathbb{C}). $$

In this case, the function f is referred to as the eigenvalue (or spectral) symbol of the matrix-sequence {A_n}_n.

Informally, this means that if we have a continuous a.e. and real f, an equispaced grid θ_{j, n} in its domain, and the eigenvalues suitably ordered, then, λ_j(A_n) ≈ f(θ_{j, n}) for all j = 1, … , n, when n is large enough.

The cases of interest in this paper are those in which A_n = T_n(f) and $\{T_n(f)\}_n\not \sim _{\lambda } f$ (except for Examples 2 and 6, where we show that the same approach works for $\{T_n(f)\}_n\sim _{\lambda } f$) where f is complex-valued (and T_n(f) is non-Hermitian), but the eigenvalues of T_n(f) are real for all n. We believe that in these cases there exist a real function $\mathfrak {f}$ such that $\{T_n(f)\}_n\sim _{\lambda } \mathfrak {f}$ and the eigenvalues of T_n(f) admit an asymptotic expansion of the same type as considered in previous works [2,3,4, 7, 8, 15,16,17,18]. We therefore formulate the following working hypothesis.

Working Hypothesis 1

Suppose that the eigenvalues of T_n(f) are real for all n. Then, for every integer α ≥ 0, every n and every j = 1, … , n, the following asymptotic expansion holds:

$$ \begin{array}{@{}rcl@{}} \lambda_{j}(T_{n}(f))&=&\mathfrak{f}(\theta_{j,n})+\sum\limits_{k=1}^{\alpha} c_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha},\\ &=&\sum\limits_{k=0}^{\alpha} c_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha}, \end{array} $$

(2)

where:

• the eigenvalues of T_n(f) are arranged in non-decreasing order, λ₁(T_n(f)) ≤… ≤ λ_n(T_n(f));

• {c₀, c₁, c₂, c₃, …} is a sequence of functions from (0, π) to $\mathbb R$ which depends only on f, and $\mathfrak {f}:= c_0$;

• $h=\frac {1}{n+1}$ and $\theta _{j,n}=\frac {j\pi }{n+1}=j\pi h$;

• E_{j, n, α} = O(h^α+ 1) is the remainder (the error), which satisfies the inequality |E_{j, n, α}|≤ C_αh^α+ 1 for some constant C_α depending only on α, f.

Remark 1

In item three of the working hypothesis, the grid θ_{j, n} is defined, but for some matrices, it is advantageous to define a different equispaced grid, e.g., by choosing h = 1/n; see, e.g., [15] where the grid θ_{j, n} = jπ/n is almost the “perfect grid” ξ_{j, n} (defined, such that, $\lambda _j(T_n(f))=\mathfrak {f}(\xi _{j,n})$ for j = 1, … , n). The typically not equispaced “perfect grid” ξ_{j, n} is further discussed in [14]. When choosing a different grid θ_{j, n} in item three, the functions c₁, c₂, c₃,… in item two will differ, but $c_0=\mathfrak {f}$ remains the same; when choosing the (typically unknown) grid ξ_{j, n}, then c_k = 0 for k > 0. The functions c_k can also be defined on any interval, not just (0, π), but then h and θ_{j, n} in item three should be modified accordingly. Also, in the working hypothesis, we arrange the eigenvalues of T_n(f) in non-decreasing order; however, using a non-increasing order would result in another function $\mathfrak {f}$. The case where the eigenvalues of T_n(f) can be described by a complex-valued or non-monotone function $\mathfrak {f}$ is out of the scope of this article and warrants further research.

2 Motivation and illustrative examples

In this section, we present four examples in support of our working hypothesis. We also discuss the fact that standard double precision eigenvalue solvers (such as LAPACK, eig in Matlab, and eigvals in Julia [6]) fail to give accurate eigenvalues of certain non-Hermitian matrices T_n(f) (both with real and complex eigenvalues); see, e.g., [5, Section 7], [12, Sections 3 and 5.8], [22], and [26, Chapter II]. A simple illustrative example is the Toeplitz matrix T_n(f) generated by the symbol f(θ) = e^iθ + e^− 2iθ presented in Fig. 1; see, e.g., [5, Section 5.1]. The numerically computed eigenvalues Ψ_j(T_n(f)) (beige) and ${\Psi }_{j}(T_{n}^{\textsc {t}}(f))$ (blue), for n = 999 and using Matlab (double precision), differ and neither is correct. A high-precision computation (256 bit) [1] of the eigenvalues λ_j(T_n(f)) (green) agrees well with the theoretical prediction; see, e.g., [23, Section 7] and [10, Example 11.18]. Along the first arm of the “star,” the eigenvalues are real, and the other arms are simply rotations of these real eigenvalues.

Fig. 1

[Symbol f(θ) = e^iθ + e^− 2iθ, n = 999] The matrix T_n(f) and the numerically (incorrect) computed eigenvalues Ψ_j(T_n(f)) (beige) and ${\Psi }_{j}(T_{n}^{\textsc {t}}(f))$ (blue), using a double precision solver (Matlab). To the right is shown the eigenvalues λ_j(T_n(f)) (green), computed using a high-precision solver [1], agreeing well with the theory

High-precision computations, by using packages such as GenericLinearAlgebra.jl [1] in Julia can compute the true eigenvalues, but they are very expensive from the computational point of view. Therefore, approximating $\mathfrak {f}$ on the grid θ_{j, n} and using matrix-less methods, described in Section 3, to compute the spectrum of T_n(f) can be computationally very advantageous. Also, the presented approaches can be a valuable tool for the analysis of the spectra of non-normal Toeplitz matrices having real eigenvalues.

Here is a short description of the four examples we are going to consider, where ξ_{j, n} is the “perfect grid” defined in Remark 1.

• Example 1: T_n(f) is non-symmetric tridiagonal, $\mathfrak {f}$ is known, and the eigenvalues $\lambda _j(T_n(f))=\mathfrak {f}(\theta _{j,n})$ are known explicitly;

• Example 2: T_n(f) is symmetric pentadiagonal, $\mathfrak {f}=f$, and the eigenvalues $\lambda _j(T_n(f))=\mathfrak {f}(\xi _{j,n})$ are not known explicitly;

• Example 3: T_n(f) is non-symmetric, $\mathfrak {f}$ is known, and the eigenvalues $\lambda _j(T_n(f))=\mathfrak {f}(\xi _{j,n})$ are not known explicitly;

• Example 4: T_n(f) is non-symmetric, $\mathfrak {f}$ is not known, and the eigenvalues $\lambda _j(T_n(f))=\mathfrak {f}(\xi _{j,n})$ are not known explicitly.

Example 1

Consider the symbol

$$ \begin{array}{@{}rcl@{}} f(\theta)=\hat{f}_{-1}e^{-\mathbf{i}\theta}+\hat{f}_{0}+\hat{f}_{1}e^{\mathbf{i}\theta}. \end{array} $$

(3)

The matrix T_n(f) is tridiagonal, and there exists a function

$$ \begin{array}{@{}rcl@{}} \mathfrak{f}(\theta)=\hat{f}_{0}+2\sqrt{\hat{f}_{-1}}\sqrt{\hat{f}_{1}}\cos(\theta), \end{array} $$

(4)

such that $T_n(f)\sim T_n(\mathfrak {f})$, that is, they are similar and hence have the same eigenvalues. The eigenvalues are given explicitly by

$$ \begin{array}{@{}rcl@{}} \lambda_{j}(T_{n}(f))=\mathfrak{f}(\theta_{j,n}), \end{array} $$

(5)

where θ_{j, n} is defined in the working hypothesis. Now, choose the Fourier coefficients $\hat {f}_{-1}=-2$, $\hat {f}_0=2$, and $\hat {f}_{1}=-1$. In this case, we have

$$ f(\theta)=-2e^{-\mathbf{i}\theta} +2 -e^{\mathbf{i}\theta}, $$

and thus,

$$ \begin{array}{@{}rcl@{}} \mathfrak{f}(\theta)=2-2\sqrt{2}\cos(\theta), \end{array} $$

(6)

where the spectrum of T_n(f) is real and is described by $\mathfrak {f}(\theta )$, and the symbol f is complex-valued. The Toeplitz matrices generated by f and $\mathfrak {f}$ are given by

$$ \begin{array}{@{}rcl@{}} {T_{n}(f)=\left[ \begin{array}{rrrrr} 2&-2\\ -1&2&-2\\ &\ddots&\ddots&\ddots\\ &&\ddots&\ddots&-2\\ &&&-1&2 \end{array} \right],} \quad {T_{n}(\mathfrak{f})=\left[ \begin{array}{rrrrr} 2&-\sqrt{2}\\ -\sqrt{2}&2&-\sqrt{2}\\ &\ddots&\ddots&\ddots\\ &&\ddots&\ddots&-\sqrt{2}\\ &&&-\sqrt{2}&2 \end{array} \right]}. \end{array} $$

We also note that $T_n(\mathfrak {f})$ is a symmetrized version of T_n(f), in the sense that there exists a transformation $T_n(\mathfrak {f})=DT_n(f)D^{-1}$ where D is a diagonal matrix with elements (D)_{i, i} = γ^i− 1, and $\gamma =\sqrt {\hat {f}_{-1}}/\sqrt {\hat {f}_1}$; see, e.g., [21].

In the left panel of Fig. 2, we represent the function f (red line), $\mathfrak {f}$ (dashed black line), and the eigenvalues $\lambda _j(T_n(f))=\lambda _j(T_n(\mathfrak {f}))$ (green dots) for n = 5. In the right panel of Fig. 2, we show the function $\mathfrak {f}$ (dashed black line) on the interval [0, π] only (since it is even on [−π, π]) and the eigenvalues $\lambda _j(T_5(f))=\lambda _j(T_5(\mathfrak {f}))=\mathfrak {f}(\theta _{j,5})$ (green dots).

Fig. 2

[Example 1: Symbol f(θ) = − 2e^−iθ + 2 − e^iθ] Left: Plots in the complex plane of f(θ) (red line), $\mathfrak {f}(\theta )=2-2\sqrt {2}\cos \limits (\theta )$ (dashed black line), and $\lambda _j(T_n(f))=\lambda _j(T_n(\mathfrak {f}))$ for n = 5 (green dots). Right: Plots of $\mathfrak {f}$ and $\lambda _j(T_5(f))=\lambda _j(T_5(\mathfrak {f}))=\mathfrak {f}(\theta _{j,n})$

In Fig. 3, we present the symbol f(θ) (red line), the numerically computed spectra Ψ_j(T_n(f)) (beige dots), and ${\Psi }_j(T_n^{\mathrm {T}}(f))$ (blue dots), for n = 1000, using a standard double precision eigenvalue solver. The analytical spectrum, defined by (5) and (6), is also shown (green dots). These numerically computed eigenvalues Ψ_j(A_n), for non-Hermitian A_n, are related to the pseudospectrum, discussed for example in [5, 12, 22, 26].

Fig. 3

[Example 1: Symbol f(θ) = − 2e^−iθ + 2 − e^iθ] Symbol f(θ) (red line), the numerically computed spectra (using a standard double precision eigenvalue solver) Ψ_j(T₁₀₀₀(f)) (beige dots), ${\Psi }_j(T_{1000}^{\mathrm {T}}(f))$ (blue dots), and the analytical spectrum $\lambda _j(T_{1000}(f))=\mathfrak {f}(\theta _{j,1000})$ (green dots)

Example 2

In this example, we consider the symbol

$$ \begin{array}{@{}rcl@{}} f(\theta)=(2-2\cos(\theta))^{2}=6-8\cos(\theta)+2\cos(2\theta), \end{array} $$

which generates a Toeplitz matrix T_n(f) associated with the second order finite difference approximation of the bi-Laplacian,

$$ \begin{array}{@{}rcl@{}} { T_{n}(f)=\left[ \begin{array}{rrrrrrrrrrr} 6&-4&1\vphantom{\ddots}\\ -4&6&-4&1\vphantom{\ddots}\\ 1&-4&6&-4&1\vphantom{\ddots}\\ &\ddots&\ddots&\ddots&\ddots&\ddots\\ &&\ddots&\ddots&\ddots&\ddots&1\\ &&&\ddots&\ddots&\ddots&-4\\ &&&&1&-4&6\vphantom{\ddots} \end{array} \right]}. \end{array} $$

The matrices T_n(f) are all Hermitian and so they have a real spectrum. Moreover, we have $f(\theta )=\mathfrak {f}(\theta )$ and $\{T_n(f)\}_n\sim _{\lambda }f$. In Fig. 4, we represent the symbol $\mathfrak {f}=f$ and the eigenvalues of T_n(f) for n = 5. The “perfect grid” ξ_{j, n} such that $\lambda _j(T_n(f))=\mathfrak {f}(\xi _{j,n})$ is not equispaced but can in this case be obtained by either computing $\xi _{j,n}=2\sin \limits ^{-1}\left ((\lambda _j(T_n(f)))^{1/4}/2\right )$ (since $f(\theta )=16\sin \limits ^4(\theta /2)$), finding the roots in (0, π) of $\mathfrak {f}(\theta )-\lambda _j(T_n(f))$ with eigenvalues given by some standard solver, or using the expansion described in [14] for large n.

Fig. 4

[Example 2: Symbol $f(\theta )=6-8\cos \limits (\theta )+2\cos \limits (2\theta )$] Symbols f(θ) (red line), $\mathfrak {f}(\theta )$ (dashed black line), and $\lambda _j(T_5(f))=\mathfrak {f}(\xi _{j,5})$ (green dots)

Example 3

In this example, we consider the following symbol

$$ \begin{array}{@{}rcl@{}} f(\theta)&=&e^{-3\mathbf{i}\theta}-4e^{-2\mathbf{i}\theta}+6e^{-\mathbf{i}\theta}-4+e^{\mathbf{i}\theta}\\ &=&e^{-\mathbf{i}\theta}\left( 6-8\cos(\theta)+2\cos(2\theta)\right)\\ &=&e^{-\mathbf{i}\theta}\left( 2-2\cos(\theta)\right)^{2}. \end{array} $$

The Toeplitz matrix T_n(f) is a shifted version of the matrix considered in Example 2 (that is, the matrix associated with the second order finite difference approximation of the bi-Laplacian), and it is given by

$$ \begin{array}{@{}rcl@{}} { T_{n}(f)=\left[ \begin{array}{rrrrrrrrrrr} -4&6&-4&1\vphantom{\ddots}\\ 1&-4&6&-4&1\vphantom{\ddots}\\ &\ddots&\ddots&\ddots&\ddots&\ddots\\ &&\ddots&\ddots&\ddots&\ddots&1\\ &&&\ddots&\ddots&\ddots&-4\\ &&&&\ddots&\ddots&6\\ &&&&&1&-4\vphantom{\ddots} \end{array} \right]}. \end{array} $$

We note that,

$$ \begin{array}{@{}rcl@{}} f(\theta)&=&e^{-\mathbf{i}\theta}\left( 2-2\cos(\theta)\right)^{2}\\ &=&e^{-3\mathbf{i}\theta}\left( 1-e^{\mathbf{i}\theta}\right)^{4}, \end{array} $$

which is equivalent to (41) in [24, Example 3.] with z = e^iθ, a = − 1, r = 3, and s = 1. Hence, by (43) in [24], we have,

$$ \begin{array}{@{}rcl@{}} \mathfrak{f}(\theta)=-\frac{\sin^{4}(\theta)}{\sin(\theta/4)\sin^{3}(3\theta/4)}, \end{array} $$

(7)

and the matrix $T_n(\mathfrak {f})$ would be full with $\lambda _j(T_n(f))\approx \lambda _j(T_n(\mathfrak {f}))\in (-\frac {(r+s)^{r+s}}{r^rs^s},0)=(-\frac {256}{27},0)$ for all j.

In the left panel of Fig. 5, we represent the functions f (red line), $\mathfrak {f}$ (dashed black line) and the eigenvalues λ_j(T_n(f)) (green dots) for n = 5. In the right panel of Fig. 5, we show the function $\mathfrak {f}$ (dashed black line) on the interval [0, π] only (since it is even on [−π, π]) and the eigenvalues $\lambda _j(T_5(f))=\mathfrak {f}(\xi _{j,5})$ (green dots).

Fig. 5

[Example 3: Symbol $f(\theta )=e^{-\mathbf {i}\theta }\left (6-8\cos \limits (\theta )+2\cos \limits (2\theta )\right )$] Left: Plots in the complex plane of f(θ) (red line), $\mathfrak {f}(\theta )=-\sin \limits ^4(\theta )/(\sin \limits (\theta /4)\sin \limits ^3(3\theta /4))$ (dashed black line), and λ_j(T_n(f)) for n = 5 (green dots). Right: Plot of $\mathfrak {f}$ and $\lambda _j(T_5(f))=\mathfrak {f}(\xi _{j,5})$

In Fig. 6, we present the symbol f (red line), the numerically computed spectra Ψ_j(T_n(f)) (beige dots) and ${\Psi }_j(T_n^{\mathrm {T}}(f))$ (blue dots), for n = 1000, using a standard double precision eigenvalue solver. The approximations of the true eigenvalues $\lambda _j(T_{1000}(f))=\mathfrak {f}(\xi _{j,1000})$ (green dots) are also shown, computed with 128 bit precision BigFloat.

Fig. 6

[Example 3: Symbol $f(\theta )=e^{-\mathbf {i}\theta }\left (6-8\cos \limits (\theta )+2\cos \limits (2\theta )\right )$] Symbol f(θ) (red line), the numerically computed spectra (using a standard double precision eigenvalue solver in Julia) Ψ_j(T₁₀₀₀(f)) (beige dots), ${\Psi }_j(T_{1000}^{\mathrm {T}}(f))$ (blue dots), and approximated spectrum $\lambda _j(T_{1000}(f))=\mathfrak {f}(\xi _{j,1000})$ (green dots)

Example 4

In this example, we consider a symbol f where the expression for $\mathfrak {f}$ is not known explicitly. Let,

$$ \begin{array}{@{}rcl@{}} f(\theta)= -e^{-4\mathbf{i}\theta} +2e^{-3\mathbf{i}\theta} -2e^{-2\mathbf{i}\theta} +9e^{-\mathbf{i}\theta} +7e^{\mathbf{i}\theta} -e^{2\mathbf{i}\theta} +e^{3\mathbf{i}\theta} , \end{array} $$

which generates the matrix

$$ \begin{array}{@{}rcl@{}} {T_{n}(f)=\left[ \begin{array}{rrrrrrrrrrrrrrr} 0&9&-2&2&-1\vphantom{\ddots}\\ 7&0&9&-2&2&-1\vphantom{\ddots}\\ -1&7&0&9&-2&2&-1\vphantom{\ddots}\\ 1&-1&7&0&9&-2&2&-1\vphantom{\ddots}\\ &\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\ &&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&-1\\ &&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&2\\ &&&&\ddots&\ddots&\ddots&\ddots&\ddots&-2\\ &&&&&\ddots&\ddots&\ddots&\ddots&9\\ &&&&&&1&-1&\hphantom{-}7&0\vphantom{\ddots} \end{array} \right]}. \end{array} $$

From [24, Example 4.] we have strong indications that approximately λ_j(T_n(f)) ∈ [− 22.09,14.96] for all j.

In the left panel of Fig. 7, we represent the symbol f (red line) and the eigenvalues λ₁₀₀₀(T_n(f)), computed with 256 bit precision BigFloat (black dots) since $\mathfrak {f}$ is not known. The eigenvalues λ_j(T_n(f)) (green dots) for n = 5 are also shown. In the right panel of Fig. 7, we again show the eigenvalues λ₁₀₀₀(T_n(f)) arranged in non-decreasing order (black dots) since $\mathfrak {f}$ is not known on the interval [0, π], and it is even on [−π, π]. Also the eigenvalues $\lambda _j(T_5(f))=\mathfrak {f}(\xi _{j,5})$ (green dots) are shown. The “perfect” grid ξ_{j, n} is computed using data from Example 8. Numerically, we have λ_j(T₁₀₀₀(f)) ∈ [− 22.0912, 14.9641] in agreement with [24]. In Fig. 8, we present the numerically computed spectra Ψ_j(T_n(f)) (beige dots) and ${\Psi }_j(T_n^{\mathrm {T}}(f))$ (blue dots), for n = 1000, using a standard double precision eigenvalue solver. The true eigenvalues $\lambda _j(T_{1000}(f))=\mathfrak {f}(\xi _{j,1000})$ (green dots) are approximated using a 256 bit precision computation.

Fig. 7

[Example 4: Symbol f(θ) = −e^− 4iθ + 2e^− 3iθ − 2e^− 2iθ + 9e^−iθ + 7e^iθ − e^2iθ + e^3iθ] Left: Symbol f(θ) (red line), λ_j(T₁₀₀₀(f)) (black dots) (since $\mathfrak {f}(\theta )$ is unknown) and λ_j(T_n(f)) for n = 5 (green dots). Right: Eigenvalues λ_j(T₁₀₀₀(f)) ordered in non-decreasing order and $\lambda _j(T_5(f))=\mathfrak {f}(\xi _{j,5})$. The grid points ξ_j,5 are computed by constructing $\tilde {\mathfrak {f}}(\theta )$ as described in Section 3.2 and finding the zeros of $\tilde {\mathfrak {f}}(\theta )-\lambda _j(T_5(f))$

Fig. 8

[Example 4: Symbol f(θ) = −e^− 4iθ + 2e^− 3iθ − 2e^− 2iθ + 9e^−iθ + 7e^iθ − e^2iθ + e^3iθ] Symbol f(θ) (red line), the numerically computed spectra Ψ_j(T₁₀₀₀(f)) (beige dots), ${\Psi }_j(T_{1000}^{\mathrm {T}}(f))$ (blue dots), and $\lambda _j(T_{1000}(f))=\mathfrak {f}(\xi _{j,n})$ (green dots)

3 Describing the real-valued eigenvalue distribution

Assuming that $\mathfrak {f}$ is a real cosine trigonometric polynomial (RCTP) symbol associated with a symbol f as in the working hypothesis, we introduce in Section 3.1 a new matrix-less method to accurately compute the expansion functions c_k, k = 0, … , α, where we recall from the working hypothesis that $c_0=\mathfrak {f}$. Subsequently, in Section 3.2, we present procedures to obtain an approximation or even the analytical expression of $\mathfrak {f}$.

3.1 Approximating the expansion functions c _k in grid points $\theta _{j,n_0}$

An asymptotic expansion of the eigenvalue errors E_{j, n,0} := E_{j, n} (when sampling the symbol f with the grid θ_{j, n} defined in the working hypothesis, under certain assumptions on f implying that $\mathfrak {f}=f$) was previously discussed in a series of papers [7, 8, 11]; such expansions can be deduced from

$$ \begin{array}{@{}rcl@{}} \lambda_{j}(T_{n}(f))=f(\theta_{j,n})+\underbrace{\sum\limits_{k=1}^{\alpha} c_{k}(\theta_{j,n})h^{k}+E_{j,n,\alpha}}_{E_{j,n}}, \end{array} $$

(8)

where θ_{j, n}, h, and E_{j, n, α} are defined in the working hypothesis.

In [7, pp. 1329], [8], and then [18], an algorithm was proposed to approximate the functions c_k(θ), which was subsequently extended and studied, see [2,3,4, 15,16,17], to cover also other types of Toeplitz-like matrices A_n, possessing an asymptotic expansion such as (8). We call this type of methods matrix-less, since they do not need to construct the large matrix A_n to approximate its eigenvalues; indeed, they approximate the functions c_k(θ) from α small matrices $A_{n_1},\ldots ,A_{n_{\alpha }}$ and then they use this approximations to compute the approximate spectrum of A_n through the formula

$$ \begin{array}{@{}rcl@{}} \lambda_{j}(A_{n})\approx \tilde{\lambda}_{j}(A_{n})= f(\theta_{j,n})+\sum\limits_{k=1}^{\alpha} \tilde{c}_{k}(\theta_{j,n})h^{k}. \end{array} $$

(9)

Assuming that the eigenvalues of T_n(f) admit an asymptotic expansion in terms of a possibly unknown function $\mathfrak {f}$ instead of f, as in our working hypothesis, we can modify Algorithm 1 in [13, Section 2.1] in order to find approximations of both $\mathfrak {f}$ and the eigenvalues of T_n(f) through the following formula, analogous to (9):

$$ \begin{array}{@{}rcl@{}} \lambda_{j}(T_{n}(f))\approx \tilde{\lambda}_{j}(T_{n}(f))&=&\sum\limits_{k=0}^{\alpha} \tilde{c}_{k}(\theta_{j,n})h^{k}\\ &=&\tilde{\mathfrak{f}}(\theta_{j,n})+\sum\limits_{k=1}^{\alpha} \tilde{c}_{k}(\theta_{j,n})h^{k}, \end{array} $$

(10)

where the approximation $\tilde {\mathfrak {f}}(\theta ):=\tilde c_0(\theta )$ of $\mathfrak {f}(\theta ):= c_0(\theta )$ is obtained from α + 1 small matrices $T_{n_0}(f),\ldots ,T_{n_{\alpha }}(f)$ as mentioned above.

Here follows an implementation in Julia of the algorithm that computes the approximations $\tilde c_k(\theta )$ for k = 0, … , α; the algorithm is written for clarity and not performance. All computations in this article are made with Julia 1.6.0 [6], using Float64 and BigFloat data types, and the GenericLinearAlgebra.jl package [1].

Using the output $\tilde {c}_k(\theta _{j,n_0})$, we can employ the interpolation–extrapolation technique described in [17] to efficiently compute very accurate approximations of $\tilde c_k(\theta )$ and, through (10), the eigenvalues of T_n(f) for an arbitrarily large order n. In the next section, we focus on the use of the approximations $\tilde {c}_0=\tilde {\mathfrak {f}}$ to describe $\mathfrak {f}$.

3.2 Constructing a function $\mathfrak {f}$ from approximations $\tilde {\mathfrak {f}}(\theta _{j,n_0})=\tilde {c}_0(\theta _{j,n_0})$

We here assume, for the sake of simplicity, that the sought function $\mathfrak {f}$ is real and even, so that it admits a cosine Fourier series of the form

$$ \begin{array}{@{}rcl@{}} \mathfrak{f}(\theta)=\hat{\mathfrak{f}}_{0}+2\sum\limits_{k=1}^{\infty} \hat{\mathfrak{f}}_{k}\cos(k\theta),\qquad \hat{\mathfrak{f}}_{k}\in\mathbb{R}. \end{array} $$

(11)

As we shall see, if $\mathfrak {f}$ is a real cosine trigonometric polynomial (RCTP), that is, a function of the form

$$ \begin{array}{@{}rcl@{}} \mathfrak{f}(\theta)=\hat{\mathfrak{f}}_{0}+2\sum\limits_{k=1}^{m} \hat{\mathfrak{f}}_{k}\cos(k\theta), \qquad \hat{\mathfrak{f}}_{k}\in\mathbb{R}, \end{array} $$

(12)

then we will be able to recover the exact expression of $\mathfrak {f}$ (see Examples 5 and 6); otherwise, we will get a truncated representation of the Fourier expansion of $\mathfrak {f}$ in (11) (see Examples 7 and 8). More specifically, what we do is the following: we consider the approximations $\tilde c_0(\theta _{j,n_0})$ provided by Algorithm 1 and we approximate the first n₀ Fourier coefficients $\hat {\mathfrak {f}}_0,\ldots ,\hat {\mathfrak {f}}_{n_0}$ with the numbers $\tilde {\hat {\mathfrak {f}}}_0,\ldots ,\tilde {\hat {\mathfrak {f}}}_{n_0}$ obtained by solving the linear system

$$ \begin{array}{@{}rcl@{}} \tilde{\hat{\mathfrak{f}}}_{0}+2\sum\limits_{k=1}^{n_{0}} \tilde{\hat{\mathfrak{f}}}_{k}\cos(k\theta_{j,n_{0}})=\tilde c_{0}(\theta_{j,n_{0}}),\qquad j=1,\ldots,n_{0}. \end{array} $$

(13)

4 Numerical examples

We now employ the proposed Algorithms 1 and 2 on matrices generated by the symbols f discussed in Examples 1–4 to highlight the applicability of the approach, in the respective Examples 5–8.

• Example 5: Only $\tilde {c}_0$ is non-zero, since θ_{j, n} gives exact eigenvalues, and the function $\mathfrak {f}$ is constructed.

• Example 6: Symbol f = c₀, and c_k, k = 1, … , 4, are recovered accurately, and the function $\mathfrak {f}=f$ is constructed.

• Example 7: Symbol $\mathfrak {f}=c_0$, and c_k, k = 1, … , 4, are recovered accurately, and a truncated RCTP representation of of $\mathfrak {f}$ is constructed.

• Example 8: Symbol $\mathfrak {f}=c_0$, and c_k, k = 1, … , 4, are constructed, and a truncated RCTP representation of of $\mathfrak {f}$ is constructed.

Example 5

We return to the non-symmetric symbol f(θ) = − 2e^−iθ + 2 − e^iθ of Example 1, and first use the proposed Algorithm 1. Note that care has to be taken when using for example standard Matlab eig command, since already for n = 160 the computed eigenvalues are wrong and complex-valued. In such a circumstance, a choice of n₀ and α needs to be such that n_α = 2^α(n₀ + 1) − 1 < 160. However, we here also use an arbitrary precision solver, GenericLinearAlgebra.jl in Julia, so we can increase precision such that theoretically any combination of n₀ and α can be chosen. The performance however decreases fast as we increase the computational precision, making obvious the applicability of the current proposed algorithms.

In Fig. 9, we present the computation of $\tilde {c}_k(\theta _{j,n_0})$, for n₀ = 31, and different precision and α. In the left panel, we show the approximated expansion functions $\tilde {c}_k$, k = 0, … , α, where α = 4, and using 128 bit precision computation. As is seen, the only non-zero $\tilde {c}_k$ is $\tilde {c}_0$, which is expected since the exact eigenvalues are given by $\mathfrak {f}(\theta _{j,n_0})=c_0(\theta _{j,n})$. In the right panel of Fig. 9, we show the absolute error in the approximation of $c_0(\theta _{j,n_0})$ for double precision computation, with α = 2, and 128 bit computation, with α = 4.

Fig. 9

[Example 5: Symbol f(θ) = − 2e^−iθ + 2 − e^iθ] The computed $\tilde {c}_k(\theta _{j,n_0})$, k = 0,…α, n₀ = 31 using Algorithm 1. Left: Computation using 128 bit precision with α = 4. Only $\tilde {c}_0(\theta _{j,n_0})=\tilde {\mathfrak {f}}(\theta _{j,n_0})$ are non-zero. Right: The absolute errors of $\tilde {c}_0(\theta _{j,n})$ compared with $\mathfrak {f}(\theta _{j,n})$ for double precision (53 bits Float64) computation, with α = 2, and 128 bit precision BigFloat, with α = 4

Now we employ Algorithm 2 to compute the Fourier coefficients of $\mathfrak {f}$. For illustrative purposes, we only show a subset of the system (13). If we choose n₀ = 2^β − 1, we have three grid points ${\theta }_{j_{\beta },n_0}$ equal to π/4, π/2 and 3π/4, corresponding to indices j_β = 2^β− 2{1,2,3}. We then have

$$ \begin{array}{@{}rcl@{}} \left[ \begin{array}{ccc} 1&\sqrt{2}&0\\ 1&0&-1\\ 1&-\sqrt{2}&0 \end{array} \right] \left[\begin{array}{c} \tilde{\hat{\mathfrak{f}}}_{0}\\ \tilde{\hat{\mathfrak{f}}}_{1}\\ \tilde{\hat{\mathfrak{f}}}_{2} \end{array} \right]\approx \left[\begin{array}{c} \tilde{c}_{0}(\pi/4)\\ \tilde{c}_{0}(\pi/2)\\ \tilde{c}_{0}(3\pi/4) \end{array} \right]. \end{array} $$

We construct the following system with $\tilde {c}_0$ computed with n₀ = 31 (that is, β = 5 and α = 2 above) using double precision, and subsequently we compute $[\tilde {\hat {\mathfrak {f}}}_0, \tilde {\hat {\mathfrak {f}}}_1, \tilde {\hat {\mathfrak {f}}}_2]^{\mathrm {T}}$,

$$ \begin{array}{@{}rcl@{}} \left[ \begin{array}{ccc} 1&\sqrt{2}&0\\ 1&0&-1\\ 1&-\sqrt{2}&0 \end{array} \right] \left[\begin{array}{c} \tilde{\hat{\mathfrak{f}}}_{0}\\ \tilde{\hat{\mathfrak{f}}}_{1}\\ \tilde{\hat{\mathfrak{f}}}_{2} \end{array} \right]\approx \left[\begin{array}{r} 0.000000000006924\\ 1.999999999889880\\ 3.999999989588136 \end{array} \right],\quad \left[\begin{array}{c} \tilde{\hat{\mathfrak{f}}}_{0}\\ \tilde{\hat{\mathfrak{f}}}_{1}\\ \tilde{\hat{\mathfrak{f}}}_{2} \end{array} \right]\approx \left[\begin{array}{r} 1.999999994797530\\ -1.414213558689498\\ -0.000000005092350 \end{array} \right]. \end{array} $$

We conclude from this computation that $\mathfrak {f}(\theta )=\hat {\mathfrak {f}}_0+2\hat {\mathfrak {f}}_1\cos \limits (\theta )+2\hat {\mathfrak {f}}_2\cos \limits (2\theta )=2-2\sqrt {2}\cos \limits (\theta )$, which is the already known analytical expression; see (5) and (6). Note also that, in this simple example, the vector containing $\tilde {c}_0(\theta _{j_{\beta },n_0})$ can be assumed to be equal to [0,2,4]^T, which would yield the exact solution (to machine precision). Using the full system (13) in Algorithm 2 yields the same result. If we now would approximate the monotonically non-increasing $\mathfrak {f}$ (instead of the non-decreasing) in Algorithm 1, the vector containing $\tilde {c}(\theta _{j_{\beta },n_0})$ would be [4,2,0]^T and would yield the symbol $\mathfrak {f}(\theta )=2+2\sqrt {2}\cos \limits (\theta )$. Obviously, the eigenvalues of $T_n(\mathfrak {f})$ are the same for both versions of $\mathfrak {f}$.

Example 6

We here return to the Hermitian symbol $f(\theta )=(2-2\cos \limits (\theta ))^2=6-8\cos \limits (\theta )+2\cos \limits (2\theta )$, as in Example 2. Since we know $\{T_n(f)\}_n\sim _{\lambda }f=\mathfrak {f}$, employing Algorithm 1 will return as $\tilde {c}_0$ an approximation of $\mathfrak {f}$, and as $\tilde {c}_k$, k > 0 the expansion functions previously obtained and studied in [3, 13].

In Fig. 10, we show in the left panel the approximated expansion functions $\tilde {c}_k(\theta _{j,n_0})$ for k = 0, … , α, computed using n₀ = 100, α = 4. Computations are made with double precision. In the right panel of Fig. 10 we show the absolute error of the approximation of $\mathfrak {f}$, that is, $\log _{10}|\mathfrak {f}(\theta _{j,n_0})-\tilde {c}_0(\theta _{j,n_0})|$, for different combinations of n₀ and α, both with high precision (128 bit BigFloat) and standard double precision (Float64). Increasing the precision does not reduce the error, unless the error is close to machine epsilon (black line for Float64). We note that increasing α has larger inpact than increasing n₀ to reduce the errors.

Fig. 10

[Example 6: Symbol $f(\theta )=6-8\cos \limits (\theta )+2\cos \limits (2\theta )$] Left: The approximated expansion functions $\tilde {c}_k(\theta _{j,n_0}), k=0,\ldots ,\alpha$ for n₀ = 100 and α = 4. Right: The absolute error $\log _{10}|\mathfrak {f}(\theta _{j,n_1})-\tilde {c}_0(\theta _{j,n_1})|$, for different combinations of n₀ and α. All computations made with 128 bit BigFloat (lines) and Float64 (dots)

The erratic behavior of $\tilde {c}_4(\theta )$ close to θ = 0 in the left panel and the increased error close to to θ = 0 in the right panel are due to the fact that the symbol f violates the so-called simple-loop conditions, discussed in [3, 13].

Using Algorithm 2, we compute approximations of the Fourier coefficients of the mononically increasing $\mathfrak {f}$ to be $\tilde {\hat {\mathfrak {f}}}_0=6,\tilde {\hat {\mathfrak {f}}}_1=-4$, $\tilde {\hat {\mathfrak {f}}}_2=1$, and $\tilde {\hat {\mathfrak {f}}}_k=0$ for k > 2. We thus recover the true symbol $\mathfrak {f}(\theta )=f(\theta )=\hat {\mathfrak {f}}_0+2\hat {\mathfrak {f}}_1\cos \limits (\theta )+2\hat {\mathfrak {f}}_2\cos \limits (2\theta )$. If Algorithm 1 was used to compute $\tilde {c}_k$ for the monotonically decreasing $\mathfrak {f}$ instead, the computed Fourier coefficients would be $\tilde {\hat {\mathfrak {f}}}_0=6,\tilde {\hat {\mathfrak {f}}}_1=4$, and $\tilde {\hat {\mathfrak {f}}}_2=1$. In fact, for $\mathfrak {f}(\theta )=6\pm 8\cos \limits (\theta )+2\cos \limits (2\theta )$, we have the same eigenvalues for T_n(f) and $T_n(\mathfrak {f})$.

Example 7

In this example, we continue the investigation of $f(\theta )=e^{-\mathbf {i}\theta }(6-8\cos \limits (\theta )+2\cos \limits (2\theta ))$ from Example 3. In Fig. 11, we show in the left panel the approximated expansion functions $\tilde {c}_k(\theta _{j,n_0})$ for n₀ = 100 and α = 4. Computations are made with 256 bit precision and the approximation $\tilde {c}_0(\theta _{j,n_0})$ overlaps well with $\mathfrak {f}$, defined in (7). Note the erratic behavior of $\tilde {c}_4$ close to θ = π. In the right panel of Fig. 11, we show the absolute values of the first one hundred approximated Fourier coefficients $\tilde {\hat {\mathfrak {f}}}_k$, given by Algorithm 2. In Table 1, we present the first ten true Fourier coefficients, $\hat {\mathfrak {f}}_k$, computed with $\mathfrak {f}$ defined in (7) and (1), and the approximations $\tilde {\hat {\mathfrak {f}}}_k$ from Algorithm 2. Since $\mathfrak {f}$ is not an RCTP, we can not recover the original simple expression of the symbol (7), but we can anyway obtain an approximated expression of $\mathfrak {f}$ through our Algorithm 2.

Fig. 11

[Example 7: Symbol $f(\theta )=e^{-\mathbf {i}\theta }(6-8\cos \limits (\theta )+2\cos \limits (2\theta ))$] Left: The approximated expansion functions $\tilde {c}_k(\theta _{j,n_0}), k=0,\ldots ,\alpha$ for n₀ = 100 and α = 4. The approximation $\tilde {c}_0(\theta _{j,n_0})$ overlaps well with $\mathfrak {f}(\theta )=-\sin \limits ^4(\theta )/(\sin \limits (\theta /4)\sin \limits ^3(3\theta /4))$. Note the erratic behavior of $\tilde {c}_4$ close to θ = π. Right: The absolute value of the approximated first one hundred Fourier coefficients, $\log _{10}|\tilde {\hat {\mathfrak {f}}}_k|$

Table 1 [Example 7: Symbol $f(\theta )=e^{-\mathbf {i}\theta }(6-8\cos \limits (\theta )$$+2\cos \limits (2\theta ))$] First ten true ($\hat {\mathfrak {f}}_k$) and computed ($\tilde {\hat {\mathfrak {f}}}_k$) Fourier coefficients of $\mathfrak {f}$. Approximations computed using n₀ = 100 and α = 4, and 256 bit precision

Example 8

Finally, we return to the non-symmetric symbol discussed in Example 4, that is, f(θ) = −e^− 4iθ + 2e^− 3iθ − 2e^− 2iθ + 9e^−iθ + 7e^iθ − e^2iθ + e^3iθ. Again, we employ Algorithms 1 and 2 to study the symbols f and $\mathfrak {f}$. In the left panel of Fig. 12, we present the approximated expansion functions in the working hypothesis, for n₀ = 100 and α = 4. Computations are made with 512 bit precision. The blue line, $\tilde {c}_0(\theta _{j,n_0})$, corresponds to the approximation of the unknown symbol $\mathfrak {f}$. Recall the curve of λ_j(T₁₀₀₀(f)) in the right panel of Fig. 7, which in principal matches the current $\tilde {c}_0$. Note how all $\tilde {c}_k$, for k > 0, are zero in apparently the same point $\theta _0\in [\frac {55\pi }{101},\frac {56\pi }{101}]$. In the right panel of Fig. 12, we see the first one hundred approximated Fourier coefficients of $\mathfrak {f}$, by using Algorithm 2. Table 2 presents the first ten approximated Fourier coefficients, $\tilde {\hat {\mathfrak {f}}}_k$.

Fig. 12

[Example 8: Symbol f(θ) = −e^− 4iθ + 2e^− 3iθ − 2e^− 2iθ + 9e^−iθ + 7e^iθ − e^2iθ + e^3iθ] Left: The approximated expansion functions $\tilde {c}_k(\theta _{j,n_0}), k=0,\ldots ,\alpha$ for n₀ = 100 and α = 4. The approximation $\tilde {c}_0(\theta _{j,n_0})$ corresponds well with λ_j(T₁₀₀₀(f)) in the right panel of Fig. 7 (since $\mathfrak {f}$ is unknown). Right: The absolute value of the approximated first one hundred Fourier coefficients, $\log _{10}|\tilde {\hat {\mathfrak {f}}}_k|$

Table 2 [Example 8: Symbol f(θ) = −e^− 4iθ + 2e^− 3iθ − 2e^− 2iθ + 9e^− 1iθ + 7e^iθ − e^2iθ + e^3iθ] First ten computed ($\tilde {\hat {\mathfrak {f}}}_k$) Fourier coefficients of the unknown $\mathfrak {f}$. Approximations computed using n₀ = 100 and α = 4, and 512 bit precision

5 Conclusions

The working hypothesis in this article concerns the existence of an asymptotic expansion such that there exists a function $\mathfrak {f}$ describing the eigenvalue distribution of the Toeplitz matrices T_n(f) generated by a symbol f. As opposed to previous versions of matrix-less methods no information of f or $\mathfrak {f}$ is needed to approximate the expansion functions $\{\mathfrak {f}:= c_0,c_1,c_2,\ldots \}$. We have shown numerically that we can recover an approximation of the function $\mathfrak {f}$. This is done by a matrix-less method described in Algorithm 1, which in principle can be modified so as to work without any information on f or the way in which the eigenvalues of the smaller versions of T_n(f) are computed. Algorithm 1 can also be used to fast and accurately compute the eigenvalues of T_n(f) for an arbitrarily large order n, as highlighted in (10). However, in this article, we have focused only on using the obtained approximation of $\mathfrak {f}$ to find an approximation of its truncated Fourier series; and if $\mathfrak {f}$ is an RCTP, we have shown that we are able to recover the original function $\mathfrak {f}$ analytically. These approaches can be a valuable tool for the exploration of the spectrum of Toeplitz and Toeplitz-like matrices previously not easily understood because of the high computational cost. For future research, we propose the extension to complex-valued functions $\mathfrak {f}$ of the results presented herein, see, e.g., [19] and the recent [9]. Also, the study of non-banded Toeplitz matrices and matrices more general than T_n(f) is of interest.