SFPN Report : Computation of Pseudospectra

AIT KHELIFA Tanina and BAILLARGUET-GAJIC Sacha

March 2023

1
Contents
1 Introduction 3

2 Definitions 3

3 GRID Algorithm 4

4 Prediction-Correction Algorithm 5

5 Parallelisation and Performance 8

6 Component-Wise Pseudospectrum 10

7 Pseudospectral Abscissa and Radius 11

8 Proof 19

References 22

2
1 Introduction
The pseudospectrum of a matrix is a set of complex numbers that indicate how
far the eigenvalues of the matrix can move under small perturbations to the
matrix. In linear algebra, the eigenvalues of a matrix are important in under-
standing the properties of the matrix, such as its invertibility, stability, and the
behavior of solutions to linear systems of equations. Pseudospectra are useful
in many areas of mathematics and science, including control theory, numeri-
cal analysis, and quantum mechanics. They can help identify potential sources
of instability in systems and provide insights into the behavior of numerical
algorithms used to solve linear systems.
Under the supervivision of Steph GRAILLAT, our work will first consist in
understanding the underlying theory of pseudospectra by studying some defini-
tions and proofs. We will then implement different algorithms that can compute
pseudospectra. After that, we will compare them and try to make them more
efficient by using parallelization. Finally we will introduce algorithms which
give information on the pseudospectrum without computing it completely.

2 Definitions
Here are some definitions of the pseudospectrum. Some are more useful for
practical matters while others help understanding the concept or demonstrating
some theoretical proofs.
Definition 1 Λε (A) := {z ∈ C|z ∈ Sp{A + E}, ||E||2 ≤ ε}

Definition 2 Λε (A) := {z ∈ C : ||(zIn − A)−1 ||2 > ε−1 }

Definition 3 Λε (A) := {z ∈ C : σmin (zIn − A) < ε}

Figure 1: Pseudospectrum of a randomly generated 4 × 4 matrix, the blue dots

being its eigenvalues

3
3 GRID Algorithm
The grid algorithm is a numerical method used to compute the pseudospectrum
of a matrix[5]
The grid algorithm works by discretizing the complex plane and computing a
singular value decomposition of (zIn − A) and then checking if the minimal
singular value is smaller than ε according to Definition 3, at each grid point z.
As we want to have the smallest grid possible for evident computing time issues,
an important question is to find a minimal set so that the pseudospectrum is
always contained. √
We demonstrated in Section 8: Proof that Λε ⊆ S, S := ∪i B(Ai,i , ri (A) + nε)
Thus we will construct a rectangular grid based on the largest and smallest
imaginary part of the points of S, as well as the largest and smallest real part.
the precision p sets the distance between 2 adjacent points of the grid.
The code in Julia below :

using LinearAlgebra
using PyPlot

function svdgrid(A, epsilon, p)

n = size(A)[1]
eigenvalues = eigvals(A)

#zone du quadrillage où on calcule le svd

r = zeros(n)
d = diag(A)
e = sqrt(n)*epsilon
for i in 1:n
r[i] = sum(abs.(A[i, :])) - abs(A[i, i])
end
gleft = minimum(real.(d) - (r .+ e))
gright = maximum(real.(d) + (r .+ e))
gbottom = minimum(imag.(d) - (r .+ e))
gtop = maximum(imag.(d) + (r .+ e))

npointsx = Int(floor((gright - gleft)/p))

npointsy = Int(floor((gtop - gbottom)/p))

x = range(gleft, stop = gright, length = npointsx)

y = range(gbottom, stop = gtop, length = npointsy)
s = complex.(zeros(npointsy,npointsx))
for j in 1:npointsx
for k in 1:npointsy
s[k,j] = (minimum(svdvals(((x[j]+y[k]*im)*I-A)::AbstractMatrix)))
end
end

4
 contour(x,y, s, levels = [epsilon])

# Tracer les valeurs propres

scatter(real(eigenvalues), imag(eigenvalues))

show()
#return s
end

The grid algorithm is relatively simple and easy to implement, but it can be
computationally expensive for large matrices and high accuracy requirements.
It is easy to parallelize as the main loop can be distributed evenly.

4 Prediction-Correction Algorithm
As its name yields, this second algorithm uses predictions and corrections in
order to compute the points that make up the pseudospectrum of a matrix A,
one component at a time. It comes from the definition 3 of the pseudospectrum.
We now consider it as a curve in the complex plane denoted by g and defined
by g(z) = ϵ = σmin (zI − A) [3]. Such an interpretation of the problem enables
the use of continuation methods. The algorithm can be broken down into three
main steps.
The first step would be to find a starting point for the very first prediction.
Let’s recall that all the eigenvalues of A are included in its pseudospectrum. One
can search along a straight line that goes through an eigenvalue, until they find
a point where it intersects the pseudospectrum of the matrix. This intersection
point z1 is the starting point for our algorithm.
The next two steps are nested in a loop which stops when convergence is
reached, that is to say when the last point computed zk+1 is close enough to the
starting point z1 . At the beginning of each iteration, we compute the singular
triplet (σmin , umin , vmin ) for the matrix zk I − A.
First, we need to try and predict a value for the next point of the pseu-
dospectrum. Instinctively, one will look for the next point along the tangent
line at the last computed point zk . Let rk be the direction of the prediction.
One can make sure that rk follows the desired direction by defining it as or-
thogonal to the gradient of g at zk , that is : rk = ivmin umin /|vmin ∗
umin | [3].
As for the length of the step we take towards zk+1 , we arbitrarily defined it as
τk = min(0.1, 10ϵ) as it had the best results. The predicted value is obtained
by z̃k+1 = zk + τk rk [3].
In order to correct the predicted point, we compute the singular triplet for
the matrix z̃k+1 I −A, and apply the following correction : zk+1 = z̃k+1 −(σmin −
∗
ϵ)/vmin umin [3].
As we have no way of knowing beforehand whether multiple eigenvalues are
included in a single component of the pseudospectrum, we have to repeat the
fore-mentioned steps for each eigenvalue of the input matrix.

5
Figure 2: Breaking down the Newton steps. Left: Prediction of z̃k along the
tangent rk of the pseudospectrum of A at the point zk−1 . Right: Correction of
the predicted point using dk , resulting in zk .

Our implementation of this algorithm in Julia :

using LinearAlgebra
using SpecialMatrices
using PyPlot

function iterate_spectrum(A, epsilon, K, d0, tol)

vp = eigvals(A)
n_vp = length(vp)
m, n = size(A)
to_plot = complex.(zeros(0))
zk = complex.(zeros(K))

for eigen in eachindex(vp)

lambda = vp[eigen]

#____STEP 0 : COMPUTE THE FIRST POINT z1____

teta = epsilon
z1_new = lambda + teta*d0
U, S, V = svd((z1_new*I - A)::AbstractMatrix)
smin = S[m]

while (abs(smin-epsilon)/epsilon > tol)

z1 = z1_new
#Computing the singular triplet for z1I - A
umin = U[:, m]
vmin = V[:, m]
z1_new = z1_new - d0*(smin - epsilon)/real(dot((conj(d0)*vmin), umin))
Ap = (z1_new*I - A)::AbstractMatrix
U, S, V = svd(Ap)

6
 smin = S[m]
end
#_________________END OF STEP 0 _______________

z1 = z1_new
fin = 1
zk[1] = z1

# _________DRAWING THE PSEUDO SPECTRUM_________

for k = 2:K
U, S, V = svd(zk[k - 1]*I - A)
i2s = CartesianIndices(A)
umin = U[:, m]
vmin = V[:, m]

#__________STEP 1 : PREDICITON_____________
rk = im*dot(vmin, umin)/abs(dot(vmin, umin))
tk = min(0.1, 10*epsilon)
z_pred = zk[k - 1] + tk*rk
#______________END OF STEP 1_______________

#___________STEP 2 : CORRECTION____________
U, S, V = svd(z_pred*I - A)
smin = minimum(S);
i2s = CartesianIndices(A)
smin = S[m]
umin = U[:, m]
vmin = V[:, m]
# Equation 2.3
z_pred = z_pred - (smin - epsilon)/dot(umin, vmin)
zk[k] = z_pred
#______________END OF STEP 2_______________

#Fin de boucle
if (k > 20) && abs(z1 - z_pred) <= 10*epsilon
fin = k
break
end
end
append!(to_plot, zk[1:fin])
end
scatter(real(to_plot), imag(to_plot))
end

7
Figure 3: Left: pseudospectrum of a randomly generated 4×4 matrix computed
with the prediction-correction algorithm. Right: Superposition of the first image
with the boundary computed with the GRID algrorithm on the same matrix

5 Parallelisation and Performance

Both the GRID and the prediction-correction algorithm can be parallelized eas-
ily and efficiently. Indeed for the first one it is sufficient to split the GRID into
equal parts while for the second one we split on the eigenvalues (from which we
compute the components).
In practice we just need to add this line of code before the main loop when
coding in Julia :

Threads.@threads for ...

We compared the efficiency of each algorithm for a random matrix of size 4

using the benchmark function "@btime" from Julia, which executes the code
many times and performs additional operations to minimize the effects of system
perturbations on performance measurements. We can see that the prediction-
correction algorithm is way faster than GRID, even for a very small matrix.
Here are the results using 4 threads for the parralelized versions:
GRID:
84.111 ms (161866 allocations: 101.16 MiB)
GRID parralelized:
31.996 ms (176810 allocations: 109.801 MiB)
pred-corr:
3.508 ms (4487 allocations: 2.91 MiB)
pred-corr parallelized:
1.091 ms (4525 allocations: 2.97 MiB)
We did other tests to measure the efficiency of the parralellized algorithms,
which is the ratio of the speedup achieved by using p threads compared to
using a single thread, normalized by the number of threads used. Using 2 to

8
5 threads on different sizes of matrices we obtained efficiencies between 60%
and 100% most of the time. In theory, this should be close to 100% because
most of the work can be distributed evenly. However in practice it also depends
on the hardware and some random factors such has the size of each connected
component of the pseudospectrum.

Figure 4: Computation time on a 4x4 random matrix using GRID depending

on the number of threads

Above using GRID the efficiency is roughly 1 because the grid can be split
evenly whatever its size and the number of threads.
On the contrary the graph below shows that the efficiency is worth for the
prediction-correction algorithm because the number of component to compute
per threads will not be equal if the number of threads don’t divide the size of
the matrix. Furthermore, some components can take more time to compute
than others.

9
Figure 5: Computation time on a 8x8 random matrix using the prediction-
correction algorithm depending on the number of threads

6 Component-Wise Pseudospectrum
Until now, we’ve considered only normwise pseudospectrum. In this section, we
will introduce componentwise pseudospectrum. It can be defined as follow for
a given ε and a perturbation matrix E:
Definition 4 Λε (A, E) := {z ∈ C : ∃∆ ∈ Cn∗n such that |∆| < εE and
A + ∆ − zI is singular } with |∆| < εE ⇔ ∀i, j, |∆i,j | < εEi,j

Here the ε-related condition is not on the norm of the perturbed matrix but
directly on its coefficients and also depends on a perturbation matrix E which
represents how we allow each individual coefficient to get perturbed.
Computing the component-wise pseudospectrum of a matrix is NP-hard[4].
Therefore, we use a satisfactory enough approximation Γε (A, E) to compute
it at a lesser cost:
Definition 5 Γε (A, E) := {z ∈ C : ρ(|(zIn − A)−1 |E) > ε−1 } with ρ(A) :=
max{|λ| : λ ∈ Sp{A}}

That way, we can use almost the same algorithm than GRID except we have
to compute a spectral radius instead of a singuar value decomposition for each
point of the grid.

10
We also have to take a slightly larger grid as we can demonstrate that Γε (A, E) ⊆
i (A) + nε). The only difference with GRID for S is that we
S, S := ∪i B(Ai,i , r√
have n instead of n.
Instead of the GRID loop, we have in Julia:

for j in 1:npointsx
for k in 1:npointsy
s[k,j] = maximum(abs.(eigvals(((abs.(inv((x[j]+y[k]*im)*I-A)))*E))))

This algorithm can be parallelized efficiently as the others with a good efficiency.

Figure 6: In pink : Normwise-pseudospectrum (with epsilon = 0.5) of a

randomly generated 4 × 4 matrix; in blue (also includes the pink surface:
componentwise-pseudospectrum (with epsilon = 0.5 and E having all its co-
efficients equal to 1) of the same matrix

7 Pseudospectral Abscissa and Radius

In some cases, one does not need to know the entire pseudospectrum of a ma-
trix. Instead, we can focus on scalar quantities derived from them, such as
the pseudospectral abscissa and the pseudospectral radius. In this section, we
will describe two criss-cross algorithms that are used to compute the previous
quantities.
The pseudospectral abscissa is the point with the largest real part on the
pseudospectrum. We can compute it with the help of a four step method.
Firstly, we need to find the rightmost eigenvalue of the matrix, λ. From this
point, the second step consists in finding the rightmost point on the pseudospec-
trum that intersects the horizontal line y = Im(λ). We will denote this point
z1 .
In the third step of the algorithm, we will look for all the points in which
the vertical line x = Re(z1 ) intersects the pseudospectrum of the matrix. We

11
will interpret this set of points as a set of intervals : one point will be the lower
bound, the other will be the higher bound.
From there, the fourth step begins with the computation of the midpoints
of these intervals. For each midpoint, we search horizontally for the rightmost
intersection with the pseudospectrum. We only keep the rightmost point found
through this search, and call it z2 .
We build the sequence (zk−1 ) by repeating step three and four until conver-
gence. When convergence is reached, the real part of the last point computed
will be the pseudosprectal abscissa.

The vertical search (which the horizontal search is easily derived from) can
be deduced from the following lemma [6] : the matrix (x + iy) − A has a singular
value ϵ if and only if iy is an eigenvalue of
x − A∗ −ϵI
 
ϵI A−x
Let iŷ be a purely imaginary eigenvalue of the previous matrix. If ϵ is the
minimal singular value of (x + iŷ) − A, then (x + iŷ) is a point at which the
vertical line though x intersects the pseudospectrum of A.
In order to perform a horizontal search, we simply multiply (x + iy) − A by
i and we work with the following matrix [6] :
−y + iA∗
 
−ϵI
ϵI iA + y
We implemented this algorithm using Julia as follows :

using LinearAlgebra
using SpecialMatrices
using PyPlot
using BlockArrays
using SymPy

function abscissa(A, epsilon, tol, iter_max=300)

#Step 1 : find the rightmost eigenvalue of A
vp = eigvals(A)
lambda = vp[argmax(real(vp))]
m, n = size(A)

#Step 2 : Find the rightmost point on the pseudospectrum z1

# s.t Im(z1) = Im(lambda)
#Going "right" with a step of epsilon
z1_new = lambda + epsilon
U, S, V = svd((z1_new*I - A)::AbstractMatrix)
smin = S[m]
while (abs(smin-epsilon)/epsilon > tol)

12
 z1 = z1_new
umin = U[:, m]
vmin = V[:, m]
z1_new = z1_new - (smin - epsilon)/real(dot(vmin, umin))
U, S, V = svd((z1_new*I - A)::AbstractMatrix)
smin = S[m]
end

#Rightmost point z1 & its abscissa x

z1 = z1_new
x = real(z1)
print("start is :", z1, "\n\n") #debug puposes
cpt = 0
scatter(x, imag(z1))

#Until convergence
while(cpt < iter_max)
#Step 3 : vertical search

#Building M by blocks (cf 42.1 & 42.2)

UL = (x*I - A')::AbstractMatrix
UR = -epsilon*I
DL = epsilon*I
DR = (A - x*I)::AbstractMatrix
M = [UL UR; DL DR]

#Getting purely imanigary eigenvalues of M

eig_mul = Dict{Any, Any}(sympy.Matrix(M).eigenvals())
eigv = (collect(keys(eig_mul)))
im_eigM = convert.(Complex{Float64}, N.(eigv))
im_eigM2 = deepcopy(im_eigM)
im_eigM = im_eigM[abs.(real.(im_eigM)) .< 1e-15]
tmp = deepcopy(im_eigM)
im_eigM = []

#Only keeping the ones s.t ((x + iŷ) - A)

#has minimum singular value epsilon
for ima in tmp
if(abs(minimum(svdvals(((x + im*imag(ima))*I
- A)::AbstractMatrix)) - epsilon) < 1e-15)
append!(im_eigM, ima)
end
end

#Step 4 : computing the midpoints + looking for the

13
#rightmost one
i = 1;

while(i <= length(im_eigM))

val = im_eigM[i]
#Index of the current eigenvalue
tmp = findfirst(x -> x==val, im_eigM2)
#Current eigenvalue has x as real part, so we really
#only need its imaginary part
val = imag(val)

#Noncrossing situation <=> eigenvalue has even

#algebraic multiplicity
if iseven(eig_mul[eigv[tmp]])
midpoint = x + val;
print("EVEN :", midpoint, "'\n") #debug purposes
else
#Going to the next value as it is the top of
#current interval
i += 1
val2 = im_eigM[i]

#Computing the midpoint

midpoint = x + im*(val + imag(val2))/2
end

#y is the imaginary part of our "new" rightest point

y = imag(midpoint)

#Getting the intersection by building M by blocks

#(cf 42.1 & 42.2)
UL = (-y*I + im*(A'))::AbstractMatrix
UR = -epsilon*I
DL = epsilon*I
DR = (y*I + im*A)::AbstractMatrix
M = [UL UR; DL DR]

#Getting purely imaginary eigenvalues of M

e = eigvals(M)
e = e[real(e) .< 1e-15]
#Here, the intersection is always the largest purely
#imaginary eigval
tmp = maximum(imag.(e))
midpoint = e[findfirst(x -> imag(x)==tmp, e)]
midpoint = imag(midpoint) + im*y

14
 #scatter(real(midpoint), imag(midpoint))

#Updating the value of the new rightest point

if(real(midpoint) > real(z1_new))
z1_new = midpoint
end

#Next iteration
i += 1;
end

#if convergence is reached => end

if(abs(z1 - z1_new) < 1e-20)
print("old : ", z1, "\nnew : ", z1_new, "\n\n")
break
end
z1 = z1_new
x = real(z1_new)
cpt += 1
end
#Return the pseudospectral abscissa
return x
end

Here are the results we obtained for the abscissa using our Julia program.

Figure 7: Pseudospectrum of a randomly generated 4×4 matrix and its abscissa

The algorithm used to compute the pseudospectral radius of a matrix A

15
is really similar to the previous one. We replace the horizontal searches with
searches along a radial line through the origin, and the vertical intervals with
arc-intervals of a circle centered at the origin
Here, the first step would be to find the point on the pseudospectrum with
largest magnitude, that is the point furthest from the origin. In order to do that,
we look for the eigenvalue with largest magnitude. Now, we need to find the
largest magnitude point on the pseudospectrum that intersects the line through
the origin and the fore-mentioned eigenvalue. We will denote this point z1 .
The third step consists in finding the points at the intersections of the circle
of radius |z1 | and the pseudospectrum. Using these points, we define the arcs
of the circle that intersect Λϵ (A), with Λϵ (A) being the peudospectra of A.
The last step of the algorithm would be to compute the midpoints of the
previously defined arcs. From each midpoint, we start a radial search for the
largest magnitude intersection with the pseudospectrum. We call z2 the largest
point found, and repeat step three and four until convergence is reached.

Once again, the radial search can be derived from the vertical search. This
time we need to multiply reiθ − A by ie−iθ and apply the very same lemma to
the following matrix [6]:
 iθ ∗ 
ie A −ϵI
ϵI e−iθ A
The only new element in this algorithm would be the search for the inter-
sections of Λϵ (A) with a circle of radius r. These intersections only occur for
values of θ for which eiθ solves the generalized eigenvalue problem [6]
   
−ϵI A iθ 0 rI
v=e v
rI 0 A∗ −ϵI

Thankfully, the Julia language comes with a lot of tools for matrix computation
[2], and these values can be easily computed using a function from the Linear-
Algebra package.

The following Julia code implements this algorithm :

using LinearAlgebra
using SpecialMatrices
using PyPlot
using BlockArrays
using SymPy

# This function returns lambda s.t imagin = lambda*exp(im*angle)

function check_eig(imagin)
teta = angle(imagin)
lambda = imagin/exp(teta*im)
return lambda

16
end

function radius(A, epsilon, tol, iter_max=300)

#Step 1 : find the eigenvalue with largest magnitude
vp = eigvals(A)
lambda = vp[argmax(abs.(vp))]
m, n = size(A)

#Step 2 : Find largest magnitude point on the pseudospectrum

# z1 st arg(z1) = arg(lambda) cf 42.1 modified
UL = (im*exp(im*angle(lambda))*A')::AbstractMatrix
UR = -epsilon*I
DL = epsilon*I
DR = (im*exp(-im*angle(lambda))*A)::AbstractMatrix

M = [UL UR; DL DR]

eigM = eigvals(M)
#Getting the largest purely imaginary eigval of M
eigM = eigM[abs.(real.(eigM)) .< 1e-10]
z1 = maximum(imag.(eigM))*exp(im*angle(lambda))

#Checking first point

scatter(real(z1), imag(z1), color="blue")

z1_new = z1
cpt = 0

while(cpt < iter_max)

#Step 3 : find points at which the pseudospectrum

# intersects with the circle of radius zk
r = check_eig(z1_new)
zero = complex.(zeros(m, n))
rI = r*I
epsI = -epsilon*I

L = [epsI A; rI zero]
R = [zero rI; A' epsI]

res = eigen(L, R)
e = res.values

#Here, we only keep the eigenvalues e = 1*exp(im*teta)

test = check_eig.(e)

17
 val_kept = test[abs.(real(test) .-1) .< 1e-16]

#Step 4 : computing the midpoints + intersection with

# the pseudosprectra
i = 1
while(i <= length(val_kept))
midpoint = (val_kept[i] + val_kept[i+1])/2

#Actually the same search as step 2, except we are

#using the midpoints instead of lambda
UL = (im*exp(im*angle(midpoint))*A')::AbstractMatrix
UR = -epsilon*I
DL = epsilon*I
DR = (im*exp(-im*angle(midpoint))*A)::AbstractMatrix
M = [UL UR; DL DR]

eigM = eigvals(M)
print(eigM, "\n")
eigM = eigM[abs.(real.(eigM)) .< 1e-10]
midpoint = maximum(imag.(eigM))*exp(im*angle(midpoint))

#If current midpoint magnitude > z1 magnitude => z1

if(abs(midpoint) > abs(z1_new))
z1_new = midpoint
end

#Onto the next arc midpoint

i += 2
end

#Checking for convergence

if(abs(z1 - z1_new) < 1e-20)
break
end

#Onto the next iteration

z1 = z1_new
cpt += 1
print("Iteration : ", cpt, "\n\n")
end
#Ploting the circle of radius |z1| centered at the origin
plt.gcf().gca().add_artist(plt.Circle((0,0), abs(z1),
fill=false))
end

Here are the results we obtained for the radius using our Julia program.

18
Figure 8: Pseudospectrum of a randomly generated 4 × 4 matrix and its radius

8 Proof
Let A ∈ Mn (C). We denote:
Pn
• ri (A) := j̸=i |Ai,j |
• Di (A) := B(Ai,i , ri (A))

• ||A||2 := σmax (A)

• Λε := {z ∈ C|z ∈ Sp{A + E}, ||E||2 ≤ ε}
Some properties that we will use later:
1. z ∈ Sp{A} ⇒ z ∈ ∪i Di [1]

2. B(z1 + z2 , r) ⊆ B(z1 , r + |z2 − z1 |)

We want to show that :
√
Λε ⊆ ∪i B(Ai,i , ri (A) + nε)

Let E ∈ Mn (C) such that ||E||2 ≤ ε.

19
 n
X
ri (A + E) = |(A + E)i,j |
j̸=i
X n n
X
ri (A + E) ≤ |Ei,j | + |Ai,j |
j̸=i j̸=i
n
X
ri (A + E) ≤ ri (A) + |Ei,j |
j̸=i

So ∀i, 1 ≤ i ≤ n :
n
X
Di (A + E) ⊆ B(Ai,i + Ei,i , ri (A) + |Ei,j |) (1)
j̸=i

By property 2.:
n
X
Di (A + E) ⊆ B(Ai,i , ri (A) + |Ei,j | + |Ei,i |)
j̸=i

And by property 1.:

n
X
Sp{A + E} ⊆ ∪i B(Ai,i , ri (A) + |Ei,j |)
j

Since this holds for any E ∈ Mn (C) such that ||E||2 ≤ ε:

n
X
Λε ⊆ ∪i B(Ai,i , ri (A) + |Ei,j |)
j

Now we will show that:

n
X √
∀i, |Ei,j | ≤ nε
j

By Cauchy-Schwartz inequality:
v
n
u n √
X uX
|Ei,j | ≤ t |Ei,j |2 n
j j

20
 qP
n
Let us suppose that j |Ei,j |2 > ε. Then:

n
X
|Ei,j |2 > ε2
j
n X
X n
|Ei,j |2 > nε2
i j

T r(E ∗ E) > nε2

Xn
σi2 > nε2
i
n
X
σi2 > nσmax
2

This is absurd, so
v
u n
uX
∀i, t |Ei,j |2 ≤ ε
j

Thus we have shown that

n
X √
∀i, |Ei,j | ≤ nε
j

To conclude:
√
Λε ⊆ ∪i B(Ai,i , ri (A) + nε)

21
References
[1] Disques de Gergshgorin. https://www.bibmath.net/dico/index.php?action
=affiche&quoi=./g/gershgorin.html.

[2] Julia Documentation · The Julia Language.

[3] Brühl, M. A curve tracing algorithm for computing the pseudospectrum.
BIT Numerical Mathematics 36, 3 (Sept. 1996), 441–454.
[4] Malyshev, A., and Sadkane, M. Componentwise pseudospectrum of a
matrix. Linear Algebra and its Applications 378 (Feb. 2004), 283–288.
[5] Trefethen, L. N., and Embree, M. IX. Computation of Pseudospectra.
In IX. Computation of Pseudospectra. Princeton University Press, May 2020,
pp. 369–420.
[6] Trefethen, L. N., and Embree, M. Spectra and Pseudospectra: The
Behavior of Nonnormal Matrices and Operators. Princeton University Press,
May 2020. Publication Title: Spectra and Pseudospectra.

22