Original article
A note on computation of pseudospectra

https://doi.org/10.1016/j.matcom.2012.03.002Get rights and content

Abstract

The aim is to contribute to pseudospectra computation via a path following technique. Given a matrix An×n, we compute the branch consisting of a fixed singular value ϵ and corresponding left and right singular vectors of the parameter dependent matrix (x + iy)I  A. The fact that the branch corresponds to the smallest singular value σmin((x + iy)I  A) = ϵ is sufficient to verify at just one point of the branch due to the continuity argument. We can exploit a standard ready-made software.

Introduction

In applications concerning linear processes represented by a non-normal matrix An×n, the classical spectral information may be misleading. Instead, the information concerning ϵ-pseudospectra yields a deeper insight. For the motivation, see e.g. the monograph [16].

Out of equivalent definitions of the pseudospectra we stick to the following one: Let An×n. Given ϵ > 0,Λϵz:σmin(zIA)<ϵ,where σmin(zI  A) denotes the smallest singular value of zI  A and I is the identity matrix.

We discuss computing pseudospectra. As the basic computational tool, variants of grid techniques are used:

  • 1.

    Construct a mesh D in the complex plane which envelopes the required part of Λϵ for selected values of ϵ.

  • 2.

    Compute σmin(zI  A) at each grid point z  D.

  • 3.

    Consider the level sets (2) for the selected values of ϵ. Visualize them as the contour plots on the grid:

Λϵ=z:σmin(zIA)=ϵ.

As step 2 is concerned, iterative techniques are recommended: inverse iterations to compute the smallest eigenvalue of (zI  A)*(zI  A) see e.g. [11], or inverse Lanczos iterations which approximate the minimal singular value of (zI  A), see e.g. [3]. We point out the algorithms which take advantage of the structure of (zI  A) namely, A being shifted by zI.

In the sequel, we will use Eigtool, see [17], which is the Matlab toolbox for computing pseudospectra, as a benchmark. Eigtool uses a grid technique based on inverse Lanczos iterations, see [15].

In this contribution, we concentrate on curve-tracing techniques, see [16], § 41, p. 394. As observed in [4], the following statement is a consequence of [13]:

Theorem 1

Consider a point z0  ∂Λϵ. Let σmin be a simple singular value of z0I  A. Then ∂Λϵ is a real analytic curve in a neighborhood of z0. Let umin and vmin be the relevant left and right singular vectors. Then ivminumin/vminumin defines the anticlockwise oriented tangent direction at z0.

For the further reference, given z, we will call (σ,u,v) the minimum singular triplet of the matrix zI  A provided that σ is the smallest singular value of zI  A, while un and vn are the relevant left and right singular vectors.

Brühl [4] proposed the following algorithm to continue a connected component of ∂Λϵ from z0 which satisfies the assumptions of the above theorem. The presentation of the algorithm is taken over from [16], p. 395: Let z0  ∂Λϵ. Given a fixed step size τ, generate the sequence zk recursively as zk−1  zk,

  • (σ,u,v) = minimum singular triplet of z0I  A

  • repeat

    • zˆk=zk1+τi(v*u)/v*u     % predictor step

    • (σ,u,v) = minimum singular triplet of zˆkIA

    • zk=zˆk(σϵ)/(u*v)     % corrector step

  • until the contour is completed

It is reported that just one corrector step is sufficient to guarantee a small relative error (σ  ϵ)/ϵ.

Brühl [4] also proposed a Newton-like technique to find a point z0  ∂Λϵ on the contour from an initial guess. As an alternative, one can use Burke, Lewis, Overton [8].

The algorithm Cobra, see Bekas, Gallopoulos [2], is a parallel version of the above curve-tracing algorithm. Again, it is reported that the relevant corrector step is sufficient to be performed only once. For parallel pathfollowing in the context of the pseudospectra computation, see also Mezher, Philippe [12].

In order to treat large scale matrices, projection techniques had been proposed, see e.g. [14], [18]. For a general reference, see [16], § 40, p. 381. The projected version of Cobra is analyzed in Bekas, Gallopoulos, Simoncini [5] namely, the transfer function methodology is implemented.

Section snippets

Continuation idea

Let Bn×n. Consider a singular triplet (σ,u,v), σ, σ  0, un, vn of the matrix B:Bv=σu,B*u=σv,v*v=1,u*u=1.The vectors u and v are called left and right singular vectors related to the singular value σ. Note that σu*u=u*Bv=(u*Bv)*=v*B*u=σv*v. If σ > 0, then u*u=v*v. Note also that (v*v)=(u*u)=0. Hence, in case σ > 0, the scaling condition (4) reads asR(v*v)=1.Obviously, we may replace (5) by R(u*u)=1 or by R(v*v)+R(u*u)=2. As far as the uniqueness is concerned, if (σ,u,v) is a singular

Implementation

In order to initialize the continuation, we need to supply a starting point. Hence, we have to solve the following problem: Given ϵ > 0, find a point z0 = x0 + iy0, z0  ∂Λϵ. We proceed as follows:

Consider (xtry;ytry)2, ztry=xtry+iytry, to be an initial guess of a point on ∂Λϵ. Compute a minimum singular triplet (σmin,umin,vmin) of the matrix ztryI  A. Set

ξtry(xtry,ytry,Rumin,umin,Rvmin,vmin)2+4n.

Generically, σmin > 0 . Define f as in (6), where in (7), (8), (9) we set ϵ  σmin. In that case, f(ξtry)

Experiments

We give examples of computations on a larger scale. No projection technique, see [16], has been used.

Concluding remarks

We applied the standard pathfollowing code to trace a selected level set ∂Λϵ of the pseudospectrum. In particular, we computed a branch consisting of the fixed singular value ϵ and the corresponding left and right singular vectors of the parameter dependent matrix (x + iy)I  A. The fact that the branch corresponds to the smallest singular value σmin((x + iy)I  A) = ϵ is sufficient to verify at just one point of the branch due to the continuity argument.

Pathfollowing the curve (10) formulates the

Acknowledgements

The research was supported by the Grant Agency of the Czech Republic (Grant No. 201/06/0356) and also by the research projects MSM 6046137306 and MSM 0021620839 of The Ministry of Education, Youth and Sports, Czech Republic.

References (18)

  • BekasC. et al.

    Cobra: parallel path following for computing the matrix pseudospectrum

    Parallel Comput.

    (2001)
  • AllgowerE.L. et al.

    Numerical Continuation Methods

    (1990)
  • BraconnierT. et al.

    Computing the fields of values and pseudospectra using the Lanczos method with continuation

    BIT

    (1996)
  • BrühlM.

    A curve tracing algorithm for computing the pseudospectrum

    BIT

    (1996)
  • C.N. Bekas, E. Gallopoulos, V. Simoncini, Pseudospectra computation of large matrices, Preprint, University of...
  • DhoogeA. et al.

    MATCONT: a Matlab package for numerical bifurcation analysis of ODEs

    ACM Trans. Math. Software

    (2003)
  • W. Govaerts, Y.A. Kuznetsov, MATCONT: Continuation Software in Matlab,...
  • BurkeJ.V. et al.

    Rubust stability and a criss-cross algorithm for pseudospectra

    IMA J. Numer. Anal.

    (2003)
  • JanovskáD. et al.

    The analytic SVD: on the non-generic points on the path

    Electron. Trans. Numer. Anal. (ETNA)

    (2010)
There are more references available in the full text version of this article.

Cited by (0)

View full text