Original articleA note on computation of pseudospectra
Introduction
In applications concerning linear processes represented by a non-normal matrix , 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 . Given ϵ > 0,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:
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 be the relevant left and right singular vectors. Then defines the anticlockwise oriented tangent direction at z0.
For the further reference, given , we will call the minimum singular triplet of the matrix zI − A provided that σ is the smallest singular value of zI − A, while and 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 recursively as zk−1 ↦ zk,
= minimum singular triplet of z0I − A
repeat
% predictor step
= minimum singular triplet of
% corrector step
until the contour is completed
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 . Consider a singular triplet , , σ ≥ 0, , of the matrix B:The vectors u and are called left and right singular vectors related to the singular value σ. Note that . If σ > 0, then . Note also that . Hence, in case σ > 0, the scaling condition (4) reads asObviously, we may replace (5) by or by . As far as the uniqueness is concerned, if 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 , , to be an initial guess of a point on ∂Λϵ. Compute a minimum singular triplet of the matrix ztryI − A. Set
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)
- et al.
Cobra: parallel path following for computing the matrix pseudospectrum
Parallel Comput.
(2001) - et al.
Numerical Continuation Methods
(1990) - et al.
Computing the fields of values and pseudospectra using the Lanczos method with continuation
BIT
(1996) 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...
- 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,...
- et al.
Rubust stability and a criss-cross algorithm for pseudospectra
IMA J. Numer. Anal.
(2003) - et al.
The analytic SVD: on the non-generic points on the path
Electron. Trans. Numer. Anal. (ETNA)
(2010)