Recycling Givens rotations for the
efficient approximation of pseudospectra
of band-dominated operators

Band-dominated operators. We study bounded linear operators on the space $\ell^{2}:=\ell^{2}({\mathbb{Z}})$ of square-summable bi-infinite complex sequences $x=(x_{k})_{k\in{\mathbb{Z}}}$ with $\|x\|=\sqrt{\sum_{k\in{\mathbb{Z}}}|x_{k}|^{2}}<\infty$. Each linear operator $A$ on $\ell^{2}$ acts via matrix-vector multiplication with a bi-infinite matrix $(a_{ij})_{i,j\in{\mathbb{Z}}}$ – and vice versa. We say that $A$ is a band operator if its matrix $(a_{ij})$ is banded (i.e. supported on only finitely many diagonals) and has uniformly bounded entries, so that $A$ is a bounded linear operator. In that case, $d:=\max\{|i-j|:a_{ij}\neq 0\}$ is called the bandwidth of $A$. Moreover, $A$ is called a band-dominated operator if it is the limit, in the induced operator norm on $\ell^{2}$, of a sequence of band operators; in particular it is a bounded operator, too, and its matrix entries decay with their distance from the main diagonal.

Pseudospectra. Because the spectrum of a non-normal operator $A$ can be highly unstable under small perturbations of $A$, one is interested in the so-called ${\varepsilon}$-pseudospectrum of $A$, that is,

Here we agree upon writing $\|B^{-1}\|=\infty$ if $B$ is not invertible. The second equality sign (see e.g. [28]) shows that ${\rm spec}_{\varepsilon}A$ measures the sensitivity of ${\rm spec}\,A$ w.r.t. additive perturbations of $A$ of norm $<{\varepsilon}$. For normal operators $A$, ${\rm spec}_{\varepsilon}A$ is the ${\varepsilon}$-neigbourhood of ${\rm spec}\,A$; otherwise it is generally larger (but never smaller). The interest in pseudospectra has been increasing over the last two decades. See [28] for many more reasons to study pseudospectra and for more references.

The lower norm. As a counterpart to the operator norm $\|A\|=\sup_{\|x\|=1}\|Ax\|$, we look at the quantity

that is sometimes (by abuse of notation) called the lower norm of $A$. While $\|A\|$ is the largest singular value of $A$, $\nu(A)$ is the smallest – provided maximum/minimum exist, such as in the case of finite matrices. It is well-known (see e.g. [19, p.69f]) that $\nu(A)>0$ holds iff $A$ is injective and has a closed image; moreover, the equality

holds with $1/0:=\infty$ indicating non-invertibility of $A$. In particular, $A$ is invertible iff $\nu(A)$ and $\nu(A^{*})$ are both nonzero, in which case they coincide. Together with the definition of ${\rm spec}_{\varepsilon}A$ it follows that

Approximating the lower norm of band-dominated operators. For $x\in\ell^{2}$, we denote its support by ${\rm supp\,}x:=\{j\in{\mathbb{Z}}:x_{j}\neq 0\}$, and we say that a bounded set $J\subset{\mathbb{Z}}$ has diameter ${\rm diam\,}J:=\max\{|i-j|:i,j\in J\}$. One of the main observations of [11] (also see [13, 21]) is that the lower norm³³3A symmetric result holds for the norm, $\|A\|$, see Proposition 3.4 and inequality (ONL) in [16]. of a band-dominated operator $A$ can be realized, up to a given $\delta>0$, by a unit element $x\in\ell^{2}$ with bounded support, say of diameter less than $n\in{\mathbb{N}}$ (dependent on $\delta$, of course). So one has

for a particular $x\in\ell^{2}$ with $\|x\|=1$ and ${\rm diam\,}({\rm supp\,}x)<n$. If ${\rm supp\,}x$ were known to be contained in the discrete interval $J^{n}_{k}:=\{k+1,...,k+n\}$ with a given $k\in{\mathbb{Z}}$, then the optimal term $\|Ax\|$ in (2) could be practically computed as the lower norm / smallest singular value of the restriction of $A$ to $\ell^{2}(J^{n}_{k})$. Since ${\rm diam\,}({\rm supp\,}x)<n$, the support must be contained in some interval $J^{n}_{k}$ with $k\in{\mathbb{Z}}$. Unfortunately, this $k$ is in general not known. It “remains” to look at – and minimize over – all $k\in{\mathbb{Z}}$:

If $A$ is a band operator then $A|_{\ell^{2}(J^{n}_{k})}$ corresponds to a finite rectangular matrix (containing columns $k+1,...,k+n$ of the infinite matrix, truncated to their joint support that is finite – due to the band structure), so that the smallest singular value, $\nu(A|_{\ell^{2}(J^{n}_{k})})$, can be computed effectively. However, consideration of all $k\in{\mathbb{Z}}$ is, in general, of course practically impossible – unless the set of all restrictions $\{A|_{\ell^{2}(J^{n}_{k})}:k\in{\mathbb{Z}}\}$ is finite, e.g. when $A$ is eventually periodic or otherwise structured.

It is clear that the size $n$ has to be increased in order to decrease the error $\delta$ in (2) and (3). The analysis in [11] (also see §3 and 4 in [13]) shows, for the particular case of tridiagonal (bandwidth $d=1$) bi-infinite matrices $(a_{ij})_{i,j\in{\mathbb{Z}}}$, that $\delta$ is of the order $1/n$; more precisely,

The constant turns out to be optimal. We make use of that result by two simple steps of reduction:

1. $(i)$

   Given an accuracy $\eta>0$, approximate our band-dominated operator $A$ (with a generally full matrix) by a band operator $B$ with $\|A-B\|\leq\eta$ and use the contractivity of $\nu(\cdot)$, so that

   $$|\nu(A)-\nu(B)|\leq\|A-B\|\leq\eta,\quad\textrm{as well as}\quad|\nu(A^{*})-\nu(B^{*})|\leq\|A^{*}-B^{*}\|\leq\eta.$$

   (5)

2. $(ii)$

   Use that the matrix of the band operator $B$ is block-tridiagonal (with block size equal to the band width of $B$, see Figure 1) and that the results of [11, 13] even apply to tridiagonal matrices with operator entries⁴⁴4In that case, $|a_{j+1,j}|$ and $|a_{j-1,j}|$ in (4) are interpreted as operator norms. – hence to block-tridiagonal matrices.

We discuss further details of steps $(i)$ and $(ii)$ in Section 2.

Approximating pseudospectra of band-dominated operators. From (1) and the above approximations and bounds on the lower norm we conclude approximations and bounds on the pseudospectrum:

Inequality (3) and its counterpart for the adjoint, $A^{*}$, lead to

from which we conclude the implications

for all ${\varepsilon}>0$, and consequently

where

Concerning the approximation step $(i)$ above, by (5), we have the implications

and the same holds for the adjoints. Subtracting $\lambda I$ from $A$ and $B$ and using (1), this shows that

so that upper and lower bounds on certain pseudospectra of $B$ yield bounds on ${\rm spec}_{\varepsilon}A$. Moreover, the inclusions (8) are as tight as desired (in the Hausdorff distance) by sending $\eta\to 0$.

Existing results. The probably most natural idea to approximate ${\rm spec}_{\varepsilon}A$ is to look at the pseudospectra ${\rm spec}_{\varepsilon}A_{n}$ of the finite sections $A_{n}=(a_{ij})_{i,j=-n}^{n}$ of $A=(a_{ij})_{i,j\in{\mathbb{Z}}}$ as $n\to\infty$. In some rare cases (Toeplitz operators [24, 8], random Jacobi operators [12]), the sets ${\rm spec}_{\varepsilon}A_{n}$ indeed converge to ${\rm spec}_{\varepsilon}A$ w.r.t. the Hausdorff distance – but in general, the sequence ${\rm spec}_{\varepsilon}A_{n}$ does not converge at all; its cluster points usually contain ${\rm spec}_{\varepsilon}A$ but also further points (see e.g. [27], one speaks of spectral pollution). Even in a simple selfadjoint example such as $A={\rm diag}(...,B,B,B,...)$ with $B={0~1\choose 1~0}$, one has⁵⁵5The ${\varepsilon}$-pseudospectra are the ${\varepsilon}$-neighbourhoods of the spectra in this selfadjoint example. ${\rm spec}\,A=\{-1,1\}$, while ${\rm spec}\,A_{n}$ repeatedly switches between $\{-1,1\}$ and $\{-1,0,1\}$ as $n$ grows. As an alternative that is somewhere between spectra and pseudospectra, [17, 26] study so-called $(N,{\varepsilon})$-pseudospectra, where $2^{N}$-th powers of the resolvent and of $1/{\varepsilon}$ are compared to each other. In [5] the lower norms of rectangular submatrices are suggested for the approximation of the spectrum and the $(N,{\varepsilon})$-pseudospectrum. Needless to say, there is a large amount of literature on the selfadjoint case (see e.g. [1, 14] and the references therein).

One major problem in approximating ${\mathbb{Z}}$ by the intervals $\{-n,...,n\}$ is (besides the potential of spectral pollution) that generally, huge values of $n$ are required to capture spectral properties of $A$ properly. (Think of an infinite diagonal matrix with distinguished entries in very remote locations.) From a computational perspective, such huge sections $\{-n,...,n\}$ are too expensive. The approach of [11] (also see §3 and 4 in [13]) – that is very much in the spirit of Gershgorin and that we adopt here – replaces $\{-n,...,n\}$ with $n\to\infty$ by the family of intervals $J_{k}^{n}=\{k+1,...,k+n\}$ for all $k\in{\mathbb{Z}}$ but with $n$ of moderate size. The price that is obviously paid is the infinite amount of positions $k$ that one has to look at, so that a certain structural simplicity of the infinite matrix is required to make the approach practically feasible. The other major plus of the [11] approach is that it comes with sharp and explicit bounds (4) on the accuracy of the approximation (6), while working for the general non-normal case.

What is new here? The tridiagonal results and the ideas of transferring them to band-dominated operators via $(i)$ and $(ii)$ are from [11], therefore not new. But there are two degrees of freedom in the choice of the blocks in step $(ii)$: Firstly, the size of the blocks, say $b\in{\mathbb{N}}$, could be any number greater than or equal to the bandwidth $d\in{\mathbb{N}}$. Secondly, once this size $b$ is fixed, there are $b$ different choices for the position of the blocks inside the infinite matrix (see Figure 1.1).

We play with that second degree of freedom, arguing that there is usually no best choice (in terms of sharpness of (6)) of block positioning, and instead we consider all $b$ possibilities, thereby improving sharpness of the bounds on ${\rm spec}_{\varepsilon}A$. (We take the union of the $b$ different lower bounds and the intersection of the $b$ upper bounds.) Naively implemented, this increases the computational cost by the factor $b$. However, we present an algorithm that compensates for this increase by reusing much of the effort that was put into the computation of $\nu(A|_{\ell^{2}(J^{n}_{k})})$ for the computation of $\nu(A|_{\ell^{2}(J^{n}_{k+1})})$. This is possible due to the large overlap between the two matrices $A|_{\ell^{2}(J^{n}_{k})}$ and $A|_{\ell^{2}(J^{n}_{k+1})}$. We cannot see a similar idea to work for the $b$-sized step from $\nu(A|_{\ell^{2}(J^{n}_{k})})$ to $\nu(A|_{\ell^{2}(J^{n}_{k+b})})$ in the block matrix, though.

In a nutshell, the smallest singular value of $A|_{\ell^{2}(J^{n}_{k})}$ coincides with that of the upper triangular matrix⁶⁶6For the computation of the smallest singular value of $R_{k}$, one can use an inverse Lanczos method. $R_{k}$ from the factorization $A|_{\ell^{2}(J^{n}_{k})}=Q_{k}R_{k}$ with a unitary $Q_{k}$ that results from a sequence of Givens rotations. The key idea is now to rearrange and reuse most of these Givens rotations for the next step when $k$ is replaced by $k+1$. With this algorithm, the complexity of the computation of $\nu(A|_{\ell^{2}(J^{n}_{k+1})})=\nu(R_{k+1})$ decreases from ${\mathcal{O}}(nd^{2})$ to just ${\mathcal{O}}(nd)$, thereby compensating for the increase by a factor of $b\approx d$ that was mentioned above. The same recycling idea and the same complexity then also apply to the computation of $\nu(A|_{\ell^{2}(J^{n}_{k+2})}),\nu(A|_{\ell^{2}(J^{n}_{k+3})})$, etc.

Contents of the paper. In Section 2 we show the details of both reduction steps $(i)$ and $(ii)$. The heart of the paper is Section 3, where we present the algorithm for the computation of $\nu(A|_{\ell^{2}(J^{n}_{k+1})})$ from $\nu(A|_{\ell^{2}(J^{n}_{k})})$ by appropriately reordering Givens rotations. In Section 4 we illustrate our results in two examples with non-trivial pseudospectra. Moreover, we compare the efficiency of our algorithm with the standard QR decomposition in each step.

Recall that we call an operator $A$ on $\ell^{2}$ band-dominated if it is the limit, in the operator norm, of a sequence of band operators (which are bounded operators with a banded matrix representation). Let us denote the sets of all band and all band-dominated operators on $\ell^{2}$ by ${\rm BO}$ and ${\rm BDO}$, respectively. We make use of the results from [11, 13] for tridiagonal operators by two steps of reduction:

To simplify notation abbreviate, for $k\in{\mathbb{Z}},n\in{\mathbb{N}}$ and $\lambda\in{\mathbb{C}}$,

and treat $A^{k}_{\lambda}$ as a finite rectangular matrix. We define $\overline{A}^{k}_{\lambda}:=(A-\lambda I)^{*}|_{\ell^{2}(J^{n}_{k})}$ analogously.
As has been described in the previous section, we need to approximate $\nu(A^{k}_{\lambda})$ and $\nu(\overline{A}^{k}_{\lambda})$ for different values $\lambda\in{\mathbb{C}}$ and multiple consecutive values of $k$. This can be done by computing the smallest singular values $\sigma_{n}(A^{k}_{\lambda})$ and $\sigma_{n}(\overline{A}^{k}_{\lambda})$ which is strongly related to pseudospectra of rectangular matrices ([31]) and similar computational problems arise.
If the considered matrices $A^{k}_{\lambda}$ were square, we could compute the Schur decomposition of $A^{k}_{0}$ – thus transforming $A^{k}_{0}$ into upper rectangular form – while preserving the shift by $\lambda$. Afterwards we could compute $\sigma_{n}(A^{k}_{\lambda})$ for multiple values of $\lambda\in{\mathbb{C}}$ using a bidiagonalization method [9] on $(A_{\lambda}^{k})^{-1}$. In the rectangular case though, no shift-preserving method to reduce $A_{0}^{k}$ to a simple form appears to be known and an inverse iteration is more difficult to implement for rectangular matrices.
We can however use the fact, that for each $\lambda\in{\mathbb{C}}$ we have two sequences $(A^{k}_{\lambda})_{k}$ and $(\overline{A}^{k}_{\lambda})_{k}$ each of which contains large overlaps between consecutive matrices. We will introduce an algorithm that takes advantage of this property.
We fix $\lambda\in{\mathbb{C}}$ and $n\in{\mathbb{N}}$ and abbreviate $A^{k}:=A^{k}_{\lambda}\in{\mathbb{C}}^{(n+2d)\times n}$ for $k\in{\mathbb{Z}}$.
Let $A^{k_{0}+1},A^{k_{0}+2},\ldots,A^{k_{0}+k_{\max}}$ be a finite sequence of matrices given by (14). W.l.o.g. we consider $k_{0}=0$. We can describe the overlapping property of these matrices by

Since $\nu(A^{k})=\sigma_{n}(A^{k})$, we are interested in computing the set

where $\sigma_{n}$ denotes the smallest singular value, which can be approximated using a QR decomposition

and applying an inverse Golub-Kahan-Lanczos-Bidiagonalization method ([2, 15]), from now on abbreviated as GKLB method, to $\tilde{R}^{k}$ (i.e. applying the GKLB method to ${(\tilde{R}^{k})}^{-1}$). Since this is a unitary transformation, the singular values of $A^{k}$ and $\tilde{R}^{k}$ are the same. The inverse GKLB method requires solving two linear systems of equations in each iteration which can be achieved using backward-substitution, since $\tilde{R}^{k}$ is upper triangular.
Note that unlike convention we write $Q^{k}A^{k}=R^{k}$ instead of $(Q^{k})^{H}A^{k}=R^{k}$ to simplify notation. It is possible to compute a QR decomposition such that the banded structure of $A^{k}$ is preserved in $\tilde{R}^{k}$, i.e. $\tilde{R}^{k}$ has at most $2d+1$ consecutive non-zero diagonals. Therefore solving a linear system of equations involving $\tilde{R}^{k}$ requires only ${\mathcal{O}}(nd)$ flops. The QR decomposition (21) itself however requires ${\mathcal{O}}(nd^{2})$ operations and is therefore the bottleneck of the algorithm for large $d$.
This bottleneck is addressed in the QH-shift-algorithm which we will develop in this section. The idea of the algorithm is to use Givens rotations to compute the factorization $Q^{1}A^{1}=H^{1}$, where $H^{1}\in{\mathbb{C}}^{m\times n}$ is an upper Hessenberg-matrix⁷⁷7We say a matrix $H\in{\mathbb{C}}^{m\times n}$ is an upper Hessenberg-matrix if $H_{i,j}=0$ for all $i>j+1$ with $2d+1$ consecutive non-zero diagonals, and then reuse these rotations to factorize $A^{2},A^{3},\ldots$ the same way.
Having factorized $A^{1},A^{2},\ldots$ into Hessenberg form using unitary transformations, we only need to apply $n$ additional Givens rotations to each matrix to arrive at the QR decomposition (21). The total effort for each QR decomposition of $A^{2},A^{3},\ldots$ is only ${\mathcal{O}}(nd)$ instead of ${\mathcal{O}}(nd^{2})$.

Preliminaries. We will only use Givens rotations acting on consecutive rows and define a rotation on the $i$th and $(i+1)$st row by the mapping

and

where $\overline{{\mathbb{D}}}:=\{z\in{\mathbb{C}}:\,\left|z\right|\leq 1\}$ is the closed complex unit disc. Details on the choice of $c,s$ can be found in standard literature [15, 29]. To simplify notation we will, in most cases, write $G_{i}\equiv G_{i}(c,s)$, if the choice of $c,s$ is clear from the context. This naturally leads to the problem of possibly having multiple rotations on the same row, each having different entries $c,s$ and we hope that it will be clear from the context that these Givens rotations are not the same.
In the interest of readability we will mainly use the arrow-notation introduced by Raf Vandebril et al. in [29, 30]:

The arrows in (28) each depict a Givens rotation operation, acting on the two rows in which the arrow is drawn (see axis of ordinates). The order of application of these rotations is described in the abscissa, i.e. from right to left, so that (28) represents the product $G_{4}G_{3}G_{2}G_{1}$. It is important to note the order of application of the Givens rotations, since they do not commute in general unless they act on disjoint couples of rows:

We say that a product $G_{i_{1}},G_{i_{2}},\ldots,G_{i_{l}}$ is a descending, respectively ascending, sequence of Givens rotations of length $l$, if $i_{p+1}=i_{p}-1$, respectively $i_{p+1}=i_{p}+1$, for $p=1,\ldots,l-1$. (28) is an example of a descending sequence of length $4$.

Step 1: We start the first step by computing a QH factorization of $A^{1}$ using consecutive Givens rotations. The number of subdiagonals⁸⁸8Where we define the main diagonal as the set $\{A^{k}_{i,i}:\,1\leq i\leq n\}$ is $2d$, we therefore require $2d-1$ sequences of Givens rotations to achieve Hessenberg form: