pyplot (generic function with 1 method)

xxxxxxxxxx
 
begin 
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add("PyPlot")
    Pkg.add("PlutoUI")
    Pkg.add("GraphPlot"); 
    Pkg.add("LightGraphs"); 
    Pkg.add("Colors")
​
    using PlutoUI
    using PyPlot
    using LinearAlgebra
    using SparseArrays
    
    using GraphPlot, LightGraphs,Colors
​
    function pyplot(f;width=3,height=3)
        clf()
        f()
        fig=gcf()
        fig.set_size_inches(width,height)
        fig
    end
​
    
end

8.7 s

Eigenvalue analysis for more general matrices

For 1D heat conduction we had a very special regular structure of the matrix which allowed exact eigenvalue calculations.

We need a generalization to varying coefficients, nonsymmetric problems, unstructured grids $\dots$ \ $\Rightarrow$ what can be done for general matrices ?

6.0 μs

The Gershgorin Circle Theorem

Theorem (Varga, Th. 1.11) Let $A$ be an $n \times n$ (real or complex) matrix. Let $Λ_{i}$ be the sum of the absolute values of the $i$ -th row's off-diagonal entries:

$Λ_{i} = \sum_{\begin{matrix} j = 1 \dots n \\ j \neq i \end{matrix}} | a_{i j} |$

If $λ$ is an eigenvalue of $A$ , then there exists $r$ , $1 \leq r \leq n$ such that $λ$ lies on the disk defined by the circle of radius $Λ_{r}$ around $a_{r r}$ :

$| λ - a_{r r} | \leq Λ_{r}$

5.0 μs

Proof: Assume $λ$ is an eigenvalue, $\vec{x} = (x_{1} \dots x_{n})$ is a corresponding eigenvector. Assume $\vec{x}$ is normalized such that

$max_{i = 1 \dots n} | x_{i} | = | x_{r} | = 1.$

From $A \vec{x} = λ \vec{x}$ it follows that

$\begin{aligned} λ x_{i} & = \sum_{j = 1 \dots n} a_{i j} x_{j} \\ (λ - a_{i i}) x_{i} & = \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq i \end{array}} a_{i j} x_{j} \\ | λ - a_{r r} | & = | \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} a_{r j} x_{j} | \leq \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} | a_{r j} | | x_{j} | \leq \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} | a_{r j} | = Λ_{r} \end{aligned}$

$◻$

4.0 μs

Corollary Any eigenvalue $λ \in σ (A)$ lies in the union of the disks defined by the Gershgorin circles

$\begin{aligned} λ \in ⋃_{i = 1 \dots n} {μ \in C : | μ - a_{i i} | \leq Λ_{i}} \end{aligned}$

Corollary The Gershgorin circle theorem allows to estimate the spectral radius $ρ (A)$ :

$\begin{aligned} ρ (A) \leq max_{i = 1 \dots n} \sum_{j = 1}^{n} | a_{i j} | = | | A | |_{\infty}, \\ ρ (A) \leq max_{j = 1 \dots n} \sum_{i = 1}^{n} | a_{i j} | = | | A | |_{1} . \end{aligned}$

Proof:

$\begin{aligned} | μ - a_{i i} | \leq Λ_{i} \Rightarrow | μ | \leq Λ_{i} + | a_{i i} | = \sum_{j = 1}^{n} | a_{i j} | \end{aligned}$

Furthermore, $σ (A) = σ (A^{T})$ .

$◻$

6.4 μs

This appears to be very easy to use, so let us try:

2.2 μs

gershgorin_circles (generic function with 1 method)

xxxxxxxxxx
 
function gershgorin_circles(A)
    t=0:0.01*π:2π
    
    # α is the trasnparency value.
    circle(x,y,r;α=0.3)=fill(x.+r.*cos.(t),y.+r.*sin.(t),alpha=α)   
    
    n=size(A,1)
    for i=1:n
        Λ=0
        for j=1:n
            if j!=i
                Λ+=abs(A[i,j])
            end
        end
        circle(real(A[i,i]),imag(A[i,i]),Λ,α=0.5/n)
    end
​
    σ=eigvals(Matrix(A))
    scatter(real(σ),imag(σ),sizes=10*ones(n),color=:red)
end

87.8 μs

xxxxxxxxxx
 
n1=5

60.0 ns

5×5 Array{Complex{Float64},2}:
 0.178502+0.0192365im  0.744411+0.093342im   …  0.392995+0.0961931im
 0.285649+0.0786766im  0.196577+0.0252502im     0.564553+0.0803397im
 0.169284+0.0438933im  0.755742+0.0518099im      0.55119+0.0852875im
 0.669057+0.0866883im  0.867834+0.0229806im     0.607358+0.0276727im
 0.480013+0.0704921im  0.581088+0.0177494im     0.247122+0.0167818im

xxxxxxxxxx
 
A1=rand(n1,n1)+0.1*rand(n1,n1)*1im

3.5 μs

Complex{Float64}1

-0.518821-0.0347195im

-0.417636+0.115591im

-0.207874-0.151289im

0.321384+0.0359316im

2.48231+0.220633im

 
eigvals(A1)

32.0 μs

x
 
pyplot(width=5, height=5) do
    gershgorin_circles(A1)
    xlabel("Re")
    ylabel("Im")
    PyPlot.grid(color=:gray)
end

52.7 ms

So this is kind of cool! Let us try this out with our heat example and the Jacobi iteration matrix: $B = I - D^{- 1} A$

2.4 μs

heatmatrix1d (generic function with 1 method)

 
function heatmatrix1d(N;α=100)
    A=zeros(N,N)
    h=1/(N-1)
    A[1,1]=1/h+α
    for i=2:N-1
        A[i,i]=2/h  
    end
    for i=1:N-1
        A[i,i+1]=-1/h
    end
    for i=2:N
        A[i,i-1]=-1/h
    end
    A[N,N]=1/h+α
    A
end

50.7 μs

jacobi_iteration_matrix (generic function with 1 method)

 
jacobi_iteration_matrix(A)=I-inv(Diagonal(A))*A

20.2 μs

xxxxxxxxxx
 
N=10

176 ns

10×10 Tridiagonal{Float64,Array{Float64,1}}:
 109.0  -9.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅ 
  -9.0  18.0  -9.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅ 
    ⋅   -9.0  18.0  -9.0    ⋅     ⋅     ⋅     ⋅     ⋅      ⋅ 
    ⋅     ⋅   -9.0  18.0  -9.0    ⋅     ⋅     ⋅     ⋅      ⋅ 
    ⋅     ⋅     ⋅   -9.0  18.0  -9.0    ⋅     ⋅     ⋅      ⋅ 
    ⋅     ⋅     ⋅     ⋅   -9.0  18.0  -9.0    ⋅     ⋅      ⋅ 
    ⋅     ⋅     ⋅     ⋅     ⋅   -9.0  18.0  -9.0    ⋅      ⋅ 
    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -9.0  18.0  -9.0     ⋅ 
    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -9.0  18.0   -9.0
    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -9.0  109.0

 
A2=Tridiagonal(heatmatrix1d(N))

9.2 μs

10×10 Tridiagonal{Float64,Array{Float64,1}}:
 0.0  0.0825688   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅          ⋅ 
 0.5  0.0        0.5   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅          ⋅ 
  ⋅   0.5        0.0  0.5   ⋅    ⋅    ⋅    ⋅    ⋅          ⋅ 
  ⋅    ⋅         0.5  0.0  0.5   ⋅    ⋅    ⋅    ⋅          ⋅ 
  ⋅    ⋅          ⋅   0.5  0.0  0.5   ⋅    ⋅    ⋅          ⋅ 
  ⋅    ⋅          ⋅    ⋅   0.5  0.0  0.5   ⋅    ⋅          ⋅ 
  ⋅    ⋅          ⋅    ⋅    ⋅   0.5  0.0  0.5   ⋅          ⋅ 
  ⋅    ⋅          ⋅    ⋅    ⋅    ⋅   0.5  0.0  0.5         ⋅ 
  ⋅    ⋅          ⋅    ⋅    ⋅    ⋅    ⋅   0.5  0.0        0.5
  ⋅    ⋅          ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0825688  0.0

 
B2=jacobi_iteration_matrix(A2)

5.2 μs

ρ2

0.9421026930965894

 
ρ2=maximum(abs.(eigvals(Matrix(B2))))

56.6 μs

We have $b_{i i} = 0$ , $Λ_{i} = {\begin{cases} \frac{1}{1 + α h}, & i = 1, n \\ 1 & i = 2 \dots n - 1 \end{cases}$

We see two circles around 0: one with radius 1 and one with radius $\frac{1}{1 + α h}$

$\Rightarrow$ estimate $| λ_{i} | \leq 1$

3.4 μs

We can also caculate the value from the estimate: Gershgorin circles of $B 2$ are centered in the origin, and the spectral radius estimate just consists in the maximum of the sum of the absolute values of the row entries.

2.5 μs

ρ2_gershgorin

1.0

xxxxxxxxxx
 
ρ2_gershgorin=maximum([sum( abs.(B2[i,:])) for i=1:size(B2,1)])

54.7 ms

xxxxxxxxxx
 
pyplot(width=5, height=5) do
    gershgorin_circles(B2)
        xlabel("Re")
    ylabel("Im")
​
    PyPlot.grid(color=:gray)
end

41.1 ms

So the estimate from the Gershgorin Circle theorem is very pessimistic... Can we improve this ?

Matrices and Graphs

3.0 μs

Permutation matrices are matrices which have exactly one non-zero entry in each row and each column which has value $1$ .
There is a one-to-one correspondence permutations $π$ of the the numbers $1 \dots n$ and $n \times n$ permutation matrices $P = (p_{i j})$ such that

$\begin{aligned} p_{i j} = {\begin{cases} 1, & π (i) = j \\ 0, & else \end{cases} \end{aligned}$

Permutation matrices are orthogonal, and we have $P^{- 1} = P^{T}$
$A \to P A$ permutes the rows of $A$
$A \to A P^{T}$ permutes the columns of $A$

7.4 μs

Define a directed graph from the nonzero entries of a matrix $A = (a_{i k})$ :

Nodes: $N = {N_{i}}_{i = 1 \dots n}$
Directed edges: $E = {\vec{N_{k} N_{l}} | a_{k l} \neq 0}$
Matrix entries $\equiv$ weights of directed edges
1:1 equivalence between matrices and weighted directed graphs

6.1 μs

Create a bidirectional graph (digraph) from a matrix in Julia. Create edge labels from off-diagonal entries and node labels combined from diagonal entries and node indices.

2.4 μs

create_graph (generic function with 1 method)

xxxxxxxxxx
 
function create_graph(matrix)
    @assert size(matrix,1)==size(matrix,2) 
    n=size(matrix,1)
    g=LightGraphs.SimpleDiGraph(n)
    elabel=[]
    nlabel=Any[]
    for i in 1:n 
        push!(nlabel,"""$(i) \n $(matrix[i,i])""")
        for j in 1:n
            if i!=j && matrix[i,j]>0
                add_edge!(g,i,j)
                push!(elabel,matrix[i,j])
            end
        end
     end
    g,nlabel,elabel
end

37.0 μs

rndmatrix (generic function with 1 method)

 
# sparse random matrix with entries with limited numbers decimal values
rndmatrix(n,p)=rand(0:0.01:1,n,n).*Matrix(sprand(Bool,n,n,p))

22.4 μs

5×5 Array{Float64,2}:
 0.0   0.0   0.77  0.0   0.0
 0.0   0.0   0.0   0.0   0.0
 0.84  0.9   0.0   0.33  0.0
 0.0   0.0   0.0   0.0   0.32
 0.0   0.41  0.0   0.38  0.0

xxxxxxxxxx
 
A3=rndmatrix(5,0.3)

7.0 μs

{5, 7} directed simple Int64 graph

2Any1

"1 \n 0.0"

"2 \n 0.0"

"3 \n 0.0"

"4 \n 0.0"

"5 \n 0.0"

3Any1

0.77

0.84

0.9

0.33

0.32

0.41

0.38

 
graph3,nlabel3,elabel3=create_graph(A3)

12.1 μs

xxxxxxxxxx
 
GraphPlot.gplot(graph3,
    nodelabel=nlabel3,
    edgelabel=elabel3,
    nodefillc=RGB(1.0,0.6,0.5),
    EDGELABELSIZE=6.0,
    NODESIZE=0.1,
    EDGELINEWIDTH=1
)

77.4 μs

Matrix graph of $A 3$ is strongly connected: false
Matrix graph of $A 3$ is weakly connected: true

19.4 μs

Definition A square matrix $A$ is reducible if there exists a permutation matrix $P$ such that

$\begin{aligned} P A P^{T} = (\begin{array}{c} A_{11} & A_{12} \\ 0 & A_{22} \end{array}) \end{aligned}$

$A$ is irreducible if it is not reducible.

Theorem (Varga, Th. 1.17): $A$ is irreducible $\Leftrightarrow$ the matrix graph is strongly connected, i.e. for each ordered pair $(N_{i}, N_{j})$ there is a path consisting of directed edges, connecting them.

Equivalently, for each $i, j$ there is a sequence of consecutive nonzero matrix entries $a_{i k_{1}}, a_{k_{1} k_{2}}, a_{k_{2} k_{3}} \dots, a_{k_{r - 1} k_{r}} a_{k_{r} j}$ .

8.4 μs

The Taussky theorem

2.1 μs

Theorem (Varga, Th. 1.18) Let $A$ be irreducible. Assume that the eigenvalue $λ$ is a boundary point of the union of all the disks

$\begin{aligned} λ \in \partial ⋃_{i = 1 \dots n} {μ \in C : | μ - a_{i i} | \leq Λ_{i}} \end{aligned}$

Then, all $n$ Gershgorin circles pass through $λ$ , i.e. for $i = 1 \dots n$ ,

$\begin{aligned} | λ - a_{i i} | = Λ_{i} \end{aligned}$

4.9 μs

Proof Assume $λ$ is eigenvalue, $\vec{x}$ a corresponding eigenvector, normalized such that $max_{i = 1 \dots n} | x_{i} | = | x_{r} | = 1$ . From $A \vec{x} = λ \vec{x}$ it follows that

$\begin{aligned} (λ - a_{r r}) x_{r} & = \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} a_{r j} x_{j} \\ | λ - a_{r r} | & \leq \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} | a_{r j} | \cdot | x_{j} | \\ \leq \sum_{\begin{array}{c} j = 1 \dots n \\ j \neq r \end{array}} | a_{r j} | = Λ_{r} (*) \end{aligned}$

$λ$ is boundary point $\Rightarrow$ $| λ - a_{r r} | = \sum_{\begin{matrix} j = 1 \dots n \\ j \neq r \end{matrix}} | a_{r j} | \cdot | x_{j} | = Λ_{r}$

$\Rightarrow$ For all $p \neq r$ with $a_{r p} \neq 0$ , $| x_{p} | = 1$ .

Due to irreducibility there is at least one $p$ with $a_{r p} \neq 0$ . For this $p$ , $| x_{p} | = 1$ and equation (*) is valid (with $p$ in place of $r$ ) $\Rightarrow$ $| λ - a_{p p} | = Λ_{p}$

Due to irreducibility, this is true for all $p = 1 \dots n$ . $◻$

5.7 μs

Apply this to the Jacobi iteration matrix for the heat conduction problem: We know that $| λ_{i} | \leq 1$ , and we can see that the matrix graph is strongly connected.

Assume $| λ_{i} | = 1$ . Then $λ_{i}$ lies on the boundary of the union of the Gershgorin circles. But then it must lie on the boundary of both circles with radius $\frac{1}{1 + α h}$ and $1$ around 0.

Contradiction! $\Rightarrow$ $| λ_{i} | < 1$ , $ρ (B) < 1$ !

3.3 μs

xxxxxxxxxx
 
α=1

150 ns

xxxxxxxxxx
 
N4=5

67.0 ns

5×5 Tridiagonal{Float64,Array{Float64,1}}:
  5.0  -4.0    ⋅     ⋅     ⋅ 
 -4.0   8.0  -4.0    ⋅     ⋅ 
   ⋅   -4.0   8.0  -4.0    ⋅ 
   ⋅     ⋅   -4.0   8.0  -4.0
   ⋅     ⋅     ⋅   -4.0   5.0

xxxxxxxxxx
 
A4=Tridiagonal(heatmatrix1d(N4,α=α))

9.4 μs

5×5 Tridiagonal{Float64,Array{Float64,1}}:
 0.0  0.8   ⋅    ⋅    ⋅ 
 0.5  0.0  0.5   ⋅    ⋅ 
  ⋅   0.5  0.0  0.5   ⋅ 
  ⋅    ⋅   0.5  0.0  0.5
  ⋅    ⋅    ⋅   0.8  0.0

xxxxxxxxxx
 
B4=jacobi_iteration_matrix(A4)

4.0 μs

ρ4

0.9486832980505135

xxxxxxxxxx
 
ρ4=maximum(abs.(eigvals(Matrix(B4))))

27.5 μs

ρ4_gershgorin

1.0

xxxxxxxxxx
 
ρ4_gershgorin=maximum([sum( abs.(B4[i,:])) for i=1:size(B4,1)])

34.6 ms

{5, 8} directed simple Int64 graph

2Any1

"1 \n 0.0"

"2 \n 0.0"

"3 \n 0.0"

"4 \n 0.0"

"5 \n 0.0"

3Any1

0.8

0.5

0.5

0.5

0.5

0.5

0.5

0.8

xxxxxxxxxx
 
graph4,nlabel4,elabel4=create_graph(B4)

13.6 μs

xxxxxxxxxx
 
GraphPlot.gplot(graph4,
    # rand(n),rand(n),
    nodelabel=nlabel4,
    edgelabel=elabel4,
    nodefillc=RGB(1.0,0.6,0.5),
    EDGELABELSIZE=6.0,
    NODESIZE=0.1,
    EDGELINEWIDTH=1
)

85.8 μs

Matrix graph is strongly connected: true
Matrix graph is weakly connected: true

18.7 μs

xxxxxxxxxx
 
pyplot(width=5, height=5) do
    gershgorin_circles(B4)
        xlabel("Re")
    ylabel("Im")
​
    PyPlot.grid(color=:gray)
end

39.8 ms

Unfortunately, we don't get a quantitative estimate here.
Advantage: we don't need to assume symmetry of $A$ or spectral equivalence estimates

4.1 μs