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

It is known that the generating function f of a sequence of Toeplitz matrices {Tn(f)}n may not describe the asymptotic distribution of the eigenvalues of Tn(f) if f is not real. In this paper, we assume as a working hypothesis that, if the eigenvalues of Tn(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 f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathfrak {f}$\end{document}, appearing in this expansion is real and describes the asymptotic distribution of the eigenvalues of Tn(f). This eigenvalue symbol f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathfrak {f}$\end{document} 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 f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathfrak {f}$\end{document}. The proposed algorithm, which opposed to previous versions, does not need any information about neither f nor f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathfrak {f}$\end{document} is tested on a wide range of numerical examples; in some cases, we are even able to find the analytical expression of f\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathfrak {f}$\end{document}. Future research directions are outlined at the end of the paper.


Introduction
Given a function f ∈ L 1 ([−π, π]), the n × n Toeplitz matrix generated by f is defined as , where the numbers fk are the Fourier coefficients of f , that is, It is known that the generating function f , also known as the symbol of {T n (f )} n , describes the asymptotic distribution of the singular values of ) 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 [8,15,20] for details and [15, 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 ∼ σ f to indicate that {T n (f )} n has an asymptotic singular value distribution described by f and {T n (f )} n ∼ λ f to indicate that {T n (f )} n has an asymptotic eigenvalue distribution described by f .The cases of interest in this paper are those in which {T n (f )} n ∼ λ f and the eigenvalues of T n (f ) are real for all n.We believe that in these cases there exist a real function g such that {T n (f )} n ∼ λ g and the eigenvalues of T n (f ) admit an asymptotic expansion of the same type as considered in previous works [1,10,12,13].We therefore formulate the following working hypothesis.
Working Hypothesis.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: where: • the eigenvalues of T n (f ) are arranged in non-decreasing order, λ 1 (T n (f )) ≤ . . .≤ λ n (T n (f )); • {g ∶= c 0 , c 1 , c 2 , c 3 , . ..} is a sequence of functions from (0, π) to R which depends only on f ; • h = 1 n+1 and θ j,n = jπ n+1 = jπ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 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 g.The case where the eigenvalues of T n (f ) can be described by a complex-valued or non-monotone function g is out of the scope of this article and warrants further research.

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) fail to give accurate eigenvalues of certain matrices T n (f ); see, e.g., [3,21].High-precision computations, by using packages such as GenericLinearAlgebra.jl [16] in Julia can compute the true eigenvalues, but they are very expensive from the computational point of view.Therefore, approximating g on the grid θ j,n and using matrix-less methods [12] 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 Toplitz matrices having real eigenvalues.
Here is a short description of the four examples we are going to consider.In what follows, we denote by ξ j,n a "perfect" sampling grid, typically not equispaced, such that λ j (T n (f )) = g(ξ j,n ) for j = 1, . . ., n; such grids are discussed in [9].
• Example 1: T n (f ) is non-symmetric tridiagonal, g is known, and the eigenvalues λ j (T n (f )) = g(θ j,n ) are known explicitly; • Example 2: T n (f ) is symmetric pentadiagonal, g = f , and the eigenvalues λ j (T n (f )) = g(ξ j,n ) are not known explicitly; • Example 3: T n (f ) is non-symmetric, g is known, and the eigenvalues λ j (T n (f )) = g(ξ j,n ) are not known explicitly; • Example 4: T n (f ) is non-symmetric, g is not known, and the eigenvalues λ j (T n (f )) = g(ξ j,n ) are not known explicitly.
Example 1.Consider the symbol The matrix T n (f ) is tridiagonal, and there exist a function such that T n (f ) ∼ T n (g), that is, they are similar and hence have the same eigenvalues.The eigenvalues are given explicitly by where θ j,n is defined in the working hypothesis.Now, choose the Fourier coefficients f1 = −1, f0 = 2, and f−1 = −2.In this case, we have and the spectrum of T n (f ) is real, even though the symbol f is complex-valued.The Toeplitz matrices generated by f and g are given by We also note that T n (g) is a symmetrized version of T n (f ), in the sense that there exists a decomposition where D is a diagonal matrix with elements (D) i,i = γ i−1 , and γ = f−1 f1 ; see [17].
Example 2. In this example we consider the symbol which generates a Toeplitz matrix T n (f ) associated with the second order finite difference approximation of the bi-Laplacian, .
The matrices T n (f ) are all Hermitian and sothey have a real spectrum.Moreover, we have f (θ) = g(θ), and {T n (f )} n ∼ σ,λ f .In Figure 3 we represent the symbol g = f and the eigenvalues of T n (f ) for n = 5.The "perfect" sampling grid ξ j,n such that λ j (T n (f )) = g(ξ j,n ) is not equispaced, but can in this case be obtained by either computing , or using the expansion described in [9] for large n.Example 3. In this example we consider the following symbol 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 We note that which is equivalent to (41) in [19,Example 3.] with z = e iθ , a = −1, r = 3, and s = 1.Hence by (43) in the same article we have that and the matrix T n (g) would be full with λ j (T n (f )) ≈ λ j (T n (g)) ∈ (− (r+s) r+s r r s s , 0) = (− 256 27 , 0) for all j.In the left panel of Figure 4 we represent the functions f (dashed black line), g (red line) and the eigenvalues λ j (T n (f )) (green dots) for n = 5.In the right panel of Figure 4 we show the function g (red line) on the interval [0, π] only (since it is even on [−π, π]) and the eigenvalues λ j (T 5 (f )) = g(ξ j,5 ) (green dots).In Figure 5 we )) (red line), and λ j (T n (f )) for n = 5 (green dots).Right: Representation of g and λ j (T 5 (f )) = g(ξ j,5 ).
present the numerically computed spectra Ψ j (T n (f )) (blue dots) and Ψ j (T T n (f )) (beige dots), for n = 1000, using a standard double precision eigenvalue solver.The approximations of the true eigenvalues λ j (T 1000 (f )) = gξ j,1000 (green dots) are also shown, computed with 128 bit precision.
Example 4. In this example we consider a symbol f where the explicit formula for g is not known explicitly.Let which generates the matrix . In the left panel of Figure 6 we represent the symbol f (dashed black line) and the eigenvalues λ 1000 (T n (f )), computed with a 256 bit eigenvalue solver (dashed red line) since g is not known.The eigenvalues λ j (T n (f )) (green dots) for n = 5 are also shown.In the right panel of Figure 6 show again the eigenvalues λ 1000 (T n (f )) arranged in non-decreasing order (dashed red line) since g is not known on the interval [0, π], since it is even on [−π, π].Also the eigenvalues λ j (T 5 (f )) = g(ξ j,5 ) (green dots) are represented.The "perfect" grid ξ j,n is computed using data from Example 8. Numerically we have λ j (T 1000 (f )) ∈ [−22.0912,14.9641] in agreement with [19].
In Figure 7 we present the numerically computed spectra Ψ j (T n (f )) (blue dots) and Ψ j (T T n (f )) (beige dots), for n = 1000, using a standard double precision eigenvalue solver.The true eigenvalues λ j (T 1000 (f )) = g(ξ j,1000 ) (green dots) are approximated using a 256 bit precision computation.

Describing the real-valued eigenvalue distribution
Assuming that g is a real cosine trigonometric (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 that c 0 = g.Subsequently, in Section 3.2 we present procedures to obtain an approximation or even the analytical expression of g.

Approximating the expansion functions c k in grid points θ 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 g = f , was discussed in a series ] Symbol f (θ) (dashed black line), the numerically computed spectra Ψ j (T 1000 (f )) (blue dots), Ψ j (T T 1000 (f )) (beige dots), and λ j (T 1000 (f )) = g(ξ j,n ) (green dots).
of papers [5][6][7]; such expansion can be deduced from where θ j,n , h, and E j,n,α are defined in the working hypothesis.
An algorithm was proposed in [13] to approximate the functions c k (θ), which was subsequently extended to other types of Toeplitz-like matrices A n possessing an asymptotic expansion such as (8); see [1,[10][11][12].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 n1 , . . ., A nα and then they use this approximations to compute the approximate spectrum of A n through the formula Assuming that the eigenvalues of T n (f ) admit an asymptotic expansion in terms of an unknwn function g instead of f , as in our working hypothesis, we can use a slight modification of Algorithm 1 in [14, Section 2.1] in order to find approximations of both g and the eigenvalues of T n (f ) through the following formula, analogous to (9): where the approximation g(θ) ∶= c0 (θ) of g(θ) ∶= c 0 (θ) is obtained from α + 1 small matrices T n0 (f ), . . ., T nα (f ) as mentioned above.
Here follows an implementation in Julia of the algorithm that computes the approximations ck (θ) for k = 0, . . ., α; the algorithm is written for clarity and not performance.All computations in this article are made with Julia 1.1.0[4], using Float64 and BigFloat data types, and the GenericLinearAlgebra.jl package [16].
Algorithm 1. Approximate expansion functions c k (θ) for k = 0, . . ., α on the grid θ j,n0 .Using the output ck (θ j,n0 ), we can employ the interpolation-extrapolation technique described in [12] to efficiently compute very accurate approximations of ck (θ) 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 c0 = g to describe g.

3.
2 Constructing a function g from approximations g(θ j,n 0 ) = c0 (θ j,n 0 ) We here assume, for the sake of simplicity, that the sought function g is real and even, so that it admits a cosine Fourier series of the form As we shall see, if g is a real cosine trigonometric polynomial (RCTP), that is, a function of the form then we will be able to recover the exact expression of g (see Examples 5 and 6); otherwise, we will get a truncated representation of the Fourier expansion of g in (11)

Numerical examples
We now employ the proposed Algorithms 1 and 2 on a the symbols f discussed in Examples 1-4 to highlight the applicability of the approach, in the respective Examples 5-8.
• Example 5: Only c0 is non-zero, since θ j,n gives exact eigenvalues, and the function g is constructed.
Example 5. We return to the non-symmetric symbol f (θ) = −e iθ + 2 − 2e −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 returned eigenvalues are complex-valued (and wrong).In such a circumstance a choice of n 0 and α needs to be such that n α = 2 α (n 0 + 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 0 and α can be chosen.The performance however decreases fast as we increase the computational precision, showing the need for the current proposed algorithms.
In Figure 8 we present the computation of ck (θ j,n0 ), for n 0 = 31, and different precision and α.In the left panel we show the approximated expansion functions ck , k = 0, . . ., α, where α = 4, and 128 bit precision computation.As is seen, the only non-zero ck is c0 , which is expected since the exact eigenvalues are given by g(θ j,n0 ) = c 0 (θ j,n ).In the right panel of Figure 8 we show the absolute error in the approximation of c 0 (θ j,n0 ) for double precision computation, with α = 2, and 128 bit computation, with α = 4.
In Figure 9 we show in the left panel the approximated expansion functions ck (θ j,n0 ) for k = 0, . . ., α, computed using n 0 = 100, α = 4. Computations are made with double precision.In the right panel of Figure 9 we show the absolute error of the approximation of g, that is, log 10 g(θ j,n0 ) − c0 (θ j,n0 ) .
The erratic behavior of c4 (θ) 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 [2,14].Using Algorithm 2 we compute approximations of the Fourier coefficients of the mononically increasing g to be g0 = 6, g1 = −4, g2 = 1, and gk = 0 for k > 2. We thus recover the true symbol g(θ) = f (θ) = ĝ0 + 2ĝ 1 cos(θ) + 2ĝ 2 cos(2θ).If Algorithm 1 was used to compute ck for the monotonically decreasing g instead, the computed Fourier coefficients would be g0 = 6, g1 = 4, and g2 = 1.In fact, for g(θ) = 6 ± 8 cos(θ) + 2 cos(2θ) we have the same eigenvalues for T n (f ) and T n (g).
Example 7. In this example we continue the investigation of f (θ) = e −iθ (6−8 cos(θ)+2 cos(2θ)) from Example 3. In Figure 10 we show in the left panel the approximated expansion functions ck (θ j,n0 ) for n 0 = 100 and α = 4. Computations are made with 256 bit precision and the approximation c0 (θ j,n0 ) overlaps well with g, defined in (7).Note the erratic behavior of c4 close to θ = π.In the right panel of Figure 10 we show the absolute values of the first one houndred approximated Fourier coefficients gk , given by Algorithm 2. In Table 1 we present the first ten true Fourier coefficients, ĝk , computed with g defined in (7) and (1), and the approximations gk from Algorithm 2. Since g 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 g through our Algorithm 2.
Example 8. Finally, we return to the non-symmetric symbol discussed in Example 4, that is, f (θ) = e 3iθ − e 2iθ + 7e iθ + 9e −1iθ − 2e −2iθ + 2e −3iθ − e −4iθ .Again, we employ Algorithms 1 and 2 to study the symbols f and g.In the left panel of Figure 11 we present the approximated expansion functions in the working hypothesis, for n 0 = 100 and α = 4. Computations are made with 512 bit precision.The blue line, c0 (θ j,n0 ) corresponds to the approximation of the unknown symbol g.Recall the curve of λ j (T 1000 (f )) in the right panel of Figure 6, which in principal matches the current c0 .Note how all ck , for k > 0, are zero in apparently the same point θ 0 ∈ [ 55π 101 , 56π 101 ].In the right panel of Figure 11 we see the first one houndred approximated Fourier coefficients of g, by using Algorithm 2. In Table 2 is presented the first ten approximated Fourier coefficients, gk .

Conclusions
The working hypothesis in this article concerns the existence of an asymptotic expansion, such that there exists of a function g describing the eigenvalue distribution of the Toeplitz matrices T n (f ) generated by a symbol f .We have shown numerically that we can recover an approximation of the function g.This is done by a matrixless 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 g to find an approximation of its truncated Fourier series; and if g is an RCTP, we have shown that it we are able to recover the original function g analytically.These approaches can be a valuable tool for the exploration of the spectrum For future research we propose the extension to complex-valued functions g of the results presented herein, and also the study of matrices more general than T n (f ).
(seeExamples 7 and 8).More specifically, what we do is the following: we consider the approximations c0 (θ j,n0 ) provided by Algorithm 1 and we approximate the first n 0 Fourier coefficients ĝ0 , . . ., ĝn0 with the numbers g0 , . . ., gn0 obtained by solving the linear system Compute approximations gk of the Fourier coefficients ĝk of g(θ).