Sequences of structured matrices of increasing size arise in many scientific applications and especially in the numerical discretization of linear differential problems; for example by using Finite Differences (FD), Finite Elements (FE), Finite Volumes (FV), Discontinuous Galerkin (DG), Isogreometric Analysis (IgA). The eigenvalues \(\lambda _j(A_n)\) of matrices \(A_n\), belonging to such a sequence \(\{A_n\}_n\), can often be approximated by a regular expansion:

\begin{align} \label {eq:expansion} \lambda _j(A_n)= \sum _{k=0}^\alpha c_k(\theta _{j,n})h^k+E_{j,n,\alpha }, \qquad j=1,\ldots ,n\qquad \theta _{j,n}=\frac {j\pi }{n+1}, \end{align} where \(c_k: [-\pi ,\pi ]\to \mathbb {C}\) (\(c_0\) is called the spectral symbol and \(c_k, k>0\) are called higher order symbols) and the errors \(E_{j,n,\alpha }=\mathcal {O}(h^{\alpha +1})\).

Hence, if we know these functions \(c_k(\theta )\), or approximate them since they are often not known analytically, we can accurately (and very fast) approximate some (or all) of the eigenvalues of any matrix \(A_n\) simply by evaluating (1).

It was previously shown (under appropriate assumptions, [4, 5]) [1, 7, 8, 9, 10] that for Hermitian sequences \(\{A_n\}_n\), where \(c_0\) is known, that we can approximate \(c_k(\theta _{j,n_0}), k=1,\ldots ,\alpha \) at specified grid points \(\theta _{j,n_0}\) using so-called matrix-less methods. The name is derived from the fact that the spectrum for any matrix \(A_n\) in the sequence \(\{A_n\}_n\) can be approximated by (1) without ever constructing the matrix; only the spectrum of a few small matrices have to be computed. That is, we have equality in (1), up to machine precision, for some chosen \(n=n_0\) and \(\alpha \). These approximations \(c_k(\theta _{j,n_0})\) can then be used for interpolation-extrapolation to any grid \(\theta _{j,n}\) (for any \(n\)) to approximate \(\lambda _j(A_n)\).

In the current presentation, mainly inspired by [3], but also [6, 11], we extend the previous algorithms with two important features:

• 1. The function \(c_0\) is not needed as an input and is approximated; this is necessary for most non-Hermitian matrix sequences, but also for discretizations of problems with variable coefficients.

• 2. The algorithm can handle discretizations of variable coefficient problems, e.g., \((a(x)u'(x))'\).

We here briefly present these two new features.

1. No knowledge of \(c_0\) necessary.

We begin by presenting two simple but representative pure Toeplitz matrix sequences; one Hermitian \(\{T_n(f_1)\}_n\) and one non-Hermitian \(\{T_n(f_2)\}_n\).

\begin{align*} &\scriptstyle \hspace {1em}f_1(\theta )=6-8\cos (\theta )+2\cos (2\theta )&&\scriptstyle \hspace {3em}f_2(\theta )=-e^{\mathbf {i}\theta }+3-3e^{-\mathbf {i}\theta }+e^{-2\mathbf {i}\theta }\\ \scriptstyle T_n(f_1)&=\scriptstyle \left [\def \arraystretch {0.5} \begin{array}{rrrrrrrrr} \scriptstyle 6&\scriptstyle -4&\scriptstyle 1\\ \scriptstyle -4&\scriptstyle 6&\scriptstyle -4&\scriptstyle 1\\ \scriptstyle 1&\scriptstyle -4&\scriptstyle 6&\scriptstyle -4&\scriptstyle 1\\ &\scriptstyle \ddots &\scriptstyle \ddots &\scriptstyle \ddots &\scriptstyle \ddots &\scriptstyle \ddots \\ &&\scriptstyle 1&\scriptstyle -4&\scriptstyle 6&\scriptstyle -4&\scriptstyle 1\\ &&&\scriptstyle 1&\scriptstyle -4&\scriptstyle 6&\scriptstyle -4\\ &&&&\scriptstyle 1&\scriptstyle -4&\scriptstyle 6\\ \end {array} \right ] &&\scriptstyle T_n(f_2)=\scriptstyle \left [\def \arraystretch {0.5} \begin{array}{rrrrrrrrr} \scriptstyle 3&\scriptstyle -3&\scriptstyle 1\\ \scriptstyle -1&\scriptstyle 3&\scriptstyle -3&\scriptstyle 1\\ &\scriptstyle -1&\scriptstyle 3&\scriptstyle -3&\scriptstyle 1\\ &&\scriptstyle \ddots &\scriptstyle \ddots &\scriptstyle \ddots &\scriptstyle \ddots \\ &&&\scriptstyle -1&\scriptstyle 3&\scriptstyle -3&\scriptstyle 1\\ &&&&\scriptstyle -1&\scriptstyle 3&\scriptstyle -3\\ &&&&&\scriptstyle -1&\scriptstyle 3\\ \end {array} \right ] \end{align*} For matrices in the Hermitian sequence \(\{T_n(f_1)\}_n\), we have that the eigenvalues can be approximated by \(\lambda _j(T_n(f_1))\approx c_0(\theta _{j,n})=f_1(\theta _{j,n})\) where \(\theta _{j,n}=j\pi /(n+1)\); the spectral symbol \(c_0\) is known and equal to \(f_1\).

For matrices in the non-Hermitian sequence \(\{T_n(f_2)\}_n\), we have that \(\lambda _j(T_n(f_1))\not \approx f_1(\theta _{j,n})\), we only know that the eigenvalues lie in the convex hull of the complex valued function \(f_2\); the spectral symbol \(c_0\) is not equal to \(f_2\). For most non-Hermitian matrix sequences the \(c_0\) is not known analytically, and the new matrix-less method presented in [3, 11] does not require it to be know. However, the matrix-less method is more efficient and accurate if it is provided.

For clarity we show a Julia implementation below on how to compute a matrix \(C=[c_k(\theta _{j,n_0})]_{k,j=1}^{\alpha +1,n_0}\); the inputs are \(n_0\) (\(\approx 100\)), \(\alpha \) (\(\approx 3\)) and eigfun (a function that returns an ordered set of eigenvalues \(\lambda _j(A_n)\) for a matrix \(A_n\) in \(\{A_n\}_n\)).

function compute_C(\(n_0\), \(\alpha \), eigfun)

hs = zeros(\(\alpha \)+1)

E = zeros(\(\alpha \)+1,\(n_0\))

for kk = 1:\(\alpha \)+1

nk = 2^(kk-1)*(\(n_0\)+1)-1

jk = 2^(kk-1)*(1:\(n_0\))

hs[kk] = 1/(nk+1)

E[kk,:] = eigfun(nk)[jk]

end

V=[hs[ii]^(jj-1) for ii=1:\(\alpha \)+1, jj=1:\(\alpha \)+1]

return C=V\E

end

As is seen above, the algorithm relies on the computed spectrum for \(\alpha +1\) small matrices (of sizes \(n_k=2^{k-1}(n_0+1)-1\), for \(k=1,\ldots ,\alpha +1\)) to compute the elements of \(C\). Subesequenly \(c_k(\theta _{j,n})\) is approximated, using interpolation-extrapolation, for arbitrary \(n\), and used in (1) to approximate \(\lambda _j(A_n)\).

2. Variable coefficients.

The spectral symbol \(f\) of the 2nd order FD discretization of \((a(x)u'(x))'\) is two-dimensional, namely \(f(x,\theta )=a(x)(2-2\cos (\theta ))\), where \(f: [0,1]\times [-\pi ,\pi ]\to \mathbb {C}\); e.g., see [12].

In [3] we show that we can use the rearranged symbol [2] to compute an expansion (1) for discretizations of variable coefficient problems; i.e., we map the function \(f: [0,1]\times [-\pi ,\pi ]\to \mathbb {R}\) to a rearranged symbol \(g: [0,1]\to \mathbb {R}\). In the new matrix-less method we have \(c_0=g\).

Apart from the two main points mentioned above, we will also discuss the current framework in detail, possible extension and current developments, and possible applications.

References

• [1]  Fayyaz Ahmad, Eman Salem Al-Aidarous, Dina Abdullah Alrehaili, Sven-Erik Ekström, Isabella Furci, and Stefano Serra-Capizzano, Are the eigenvalues of preconditioned banded symmetric Toeplitz matrices known in almost closed form?, Numerical Algorithms, 78(3), pp. 867-893 (2018)

• [2]  Giovanni Barbarino, Davide Bianchi, and Carlo Garoni, Constructive approach to the monotone rearrangement of functions, Expositiones Mathematicae, 40(1), pp. 155–175 (2022)

• [3]  Giovanni Barbarino, Melker Claesson, Sven-Erik Ekström, Carlo Garoni, David Meadon, and Hendrik Speleers, Matrix-less spectral approximation for large structured matrices, BIT Numerical Mathematics, (in press, 2024) DOI:10.1007/s10543-024-01041-w

• [4]  Mauricio Barrera, Albrecht Böttcher, Sergei M. Grudsky, and Egor A. Maximenko, Eigenvalues of even very nice Toeplitz matrices can be unexpectedly erratic, In: Böttcher, A., Potts, D., Stollmann, P., Wenzel, D. (eds) The Diversity and Beauty of Applied Operator Theory. Operator Theory: Advances and Applications, vol 268. Birkhäuser, Cham. 2018

• [5]  Manuel Bogoya, Sven-Erik Ekström, and Stefano Serra-Capizzano, Fast Toeplitz eigenvalue computations, joining interpolation-extrapolation matrix-less algorithms and simple-loop theory, Numerical Algorithms, 91, pp. 1653-1676 (2022)

• [6]  Manuel Bogoya, Sven-Erik Ekström, Stefano Serra-Capizzano, and Paris Vassalos, Matrix-less methods for the spectral approximation of large non-Hermitian Toeplitz matrices: A concise theoretical analysis and a numerical study, Numerical Linear Algebra with Applications, 34(4), pp. e2545 (2024)

• [7]  Sven-Erik Ekström, Isabella Furci, Carlo Garoni, Carla Manni, Stefano Serra-Capizzano, and Hendrik Speleers, Are the eigenvalues of the B-spline isogeometric analysis approximation of \(-\Delta u=\lambda u\) known in almost closed form?, Numerical Linear Algebra with Applications, 25(5), pp. e2198 (2018)

• [8]  Sven-Erik Ekström, Isabella Furci, and Stefano Serra-Capizzano, Exact formulae and matrix-less eigensolvers for block banded symmetric Toeplitz matrices, BIT Numerical Mathematics, 58(4), pp. 937-968 (2018)

• [9]  Sven-Erik Ekström and Carlo Garoni, A matrix-less and parallel interpolation-extrapolation algorithm for computing the eigenvalues of preconditioned banded symmetric Toeplitz matrices, Numerical Algorithms, 80(3), pp. 819-848 (2019)

• [10]  Sven-Erik Ekström, Carlo Garoni, and Stefano Serra-Capizzano, Are the Eigenvalues of Banded Symmetric Toeplitz Matrices Known in Almost Closed Form?, Experimental Mathematics, 27(4), pp. 478-487 (2018)

• [11]  Sven-Erik Ekström and Paris Vassalos, A Matrix-Less Method to Approximate the Spectrum and the Spectral Function of Toeplitz Matrices with Real Eigenvalues, Numerical Algorithms, 89, pp. 701-720 (2022)

• [12]  Carlo Garoni and Stefano Serra-Capizzano, Generalized Locally Toeplitz Sequences: Theory and Applications: Volume I, Volume 1, Springer, 2017.

• [13]  Boris Shapiro, František Štampach, Non-selfajoint Toeplitz matrices whose principal submatrices have real spectrum, Constructive Approximation, 49(2), pp. 191–226 (2019)