Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Nov 2016
Lecture topic: Eigenvalues

Gaussian Elimination Modules & a Norm

Gaussian Elimination with No Output Module

Here is a routine with no intermediate output:

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Return[A];  ] ;

Matrix Inversion Module

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];
A3=Take[A2,{1,n},{n+1,2n}];
Return[A3]];

Norm

norm[v_,p_]:=((((∑_(i=1)^(Length[v]))((Abs[v[[i]]])^p)))^(1/p))

Mathematica Algebra

Mathematica can do neat things with matrices and arrays using the Inner and Outer built-in commands.This nomenclature relates to the fact that dot and cross products in three dimensional space are really special instances of inner and outer products. Here are some vectors defined to show this:

Clear[i,n,u,x,y]
n=4
m=6
vec1=Array[x,n]
vec2=Array[y,n]
vec3=Array[z,m]

Here is how a dot product is defined:

dotproduct=Inner[Times, vec1, vec2, Plus]
Dimensions[dotproduct]

And here are vector products:

vectorproduct=Outer[Times,vec1,vec2];
MatrixForm[vectorproduct]
Dimensions[vectorproduct]

vectorproduct=Outer[Times,vec1,vec3];
MatrixForm[vectorproduct]
Dimensions[vectorproduct]

vectorproduct=Outer[Times,vec3,vec1];
MatrixForm[vectorproduct]
Dimensions[vectorproduct]

You can do other things as well:

Inner[Power, -vec1, vec2, Plus]
Outer[Power, vec1, -vec2];
MatrixForm[%]

Eigensystem Basics

The n x n matrix A will take an n-vector x and map it into an n-vector y,

A x = y

A.vec[1]

Here is a picture of this mapping in two dimensions for a specific case:

n=2;
A={{1,-1},{-1,-1/10}};
vec[1]={1,4};

evecplot[p_]:=Show[Graphics[{Red,Arrow[{{0,0},vec[p]}]}],
Graphics[{Blue,Arrow[{{0,0},A.vec[p]}]}],Axes->True,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},AspectRatio->1]

Show[evecplot[1],PlotRange->{{-5,5},{-5,5}},ImageSize->{500,Automatic},BaseStyle->20]

For some vectors x, this mapping will produce a vector y which is parallel to the vector x, but scaled by some amount λ,

A x = λ x

n=2;
vec[1]={(1/20) (-11+(Sqrt[521])),1}
vec[2]={(1/20) (-11-(Sqrt[521])),1}

Show[{evecplot[1],evecplot[2]},DisplayFunction->$DisplayFunction,PlotRange->{{-3,3},{-4,4}},BaseStyle->20]

The vectors for which this is true are called eigenvectors, and each eigenvector has an associated scaling parameter λ, which is the eigenvalue. We can rewrite our system of equations as follows

A x = λ x
A x - λ x = 0
(A  - λ I)x = 0

For this system to have a nonzero solution x, we require the determinant of the matrix A-λ Ito be zero. We can solve the following equation for (λ^((k))), k=1,2,3,...,n.

Det(A  - (λ^((k))) I) = 0

Once the values of  (λ^((k))) are known, we can solve the equation

(A  - (λ^((k))) I)(x^((k))) = 0

for the eigenvalues (x^((k))).

Gerschgorin Circles

Definition

Let A be an n x n matrix and (Subscript[R, i]) denote the circle in the complex plane with center (Subscript[a, i i]) and radius (Subscript[r, i])=(∑_(j=1
j≠i)^n)|(Subscript[a, i i])|, which we can write as

(Subscript[R, i])={z∈C||z-(Subscript[a, i i])| ≤(Subscript[r, i])}

The eigenvalues of A are contained within R=(⋃_(i=1)^n)(Subscript[R, i]). Also, the union of any k of these circles that do not intersect the remaining n-k contains precisely k (counting multiplicities) of the eigenvalues.

A Module

Here is a module to produce the Gerschgorin circles for a matrix

n=3;
A={{4,1,1},{0,2,1},{-2,0,9}};
MatrixForm[A]

gcircles2[A_]:=Module[{r,i,j,n=Length[A],tab},
r=Table[0,{i,1,n}];
Do[r[[i]]=(∑_(j=1)^(i-1))Abs[A[[i,j]]] +(∑_(j=i+1)^n)Abs[A[[i,j]]],{i,1,n}] ;
tab=Table[{Red,Circle[{Re[A[[i,i]]],Im[A[[i,i]]]},r[[i]]]},{i,1,n}];
Show[Graphics[tab],Axes->True,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,BaseStyle->20]]

gcircles2[A]

There should be two eigenvalues that lie in the region enclosed by the intersecting circles, and one in the disjoint circle. We can modify our code to include the eigenvalues (using Mathematica's intrinsic command to get the eigenvalues).

gcircles[A_]:=Module[{r,i,j,n=Length[A],tab},
r=Table[0,{i,1,n}];
Do[r[[i]]=(∑_(j=1)^(i-1))Abs[A[[i,j]]] +(∑_(j=i+1)^n)Abs[A[[i,j]]],{i,1,n}] ;
tab=Table[{Red,Circle[{Re[A[[i,i]]],Im[A[[i,i]]]},r[[i]]]},{i,1,n}];
lam=N[Eigenvalues[A]];
tab2=Table[{PointSize[0.03],Blue,Point[{Re[lam[[i]]],Im[lam[[i]]]}]},{i,1,n}];
Show[Graphics[tab],Graphics[tab2],Axes->True,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,BaseStyle->20]]

gcircles[A]

We see that this is indeed the case. We are interested in Eigenvalues of the matrix, and the Gerschgorin Circles allow us to get some idea of bounds on the eigenvalues without calculating them.

A More Interesting Example

A={{3,0,-1,-1/4,1/4},{0,5,1/2,0,1},{-1/4,0,6,1/4,1/2},{0,-1,1/2,-3,1/4},{1/6,-1/6,1/3,1/3,4}};
MatrixForm[A]

gcircles[A]

More Examples

A={{4,-1},{2,2}};
MatrixForm[A]
gcircles[A]

SeedRandom[3];
A=Table[If[i==j,3 (i^3),RandomInteger[{-50,50}]],{i,1,7},{j,1,7}];
A
gcircles[A]

A=HilbertMatrix[10];
MatrixForm[A]
gcircles[A]

SeedRandom[3];
A=Table[i (j^7) RandomComplex[],{i,1,5},{j,1,5}];
SetPrecision[A,2]
gcircles[A]

More Examples

This just changes the output range:

gcircles[A_]:=Module[{r,i,j,n=Length[A],tab},
r=Table[0,{i,1,n}];
Do[r[[i]]=(∑_(j=1)^(i-1))Abs[A[[i,j]]] +(∑_(j=i+1)^n)Abs[A[[i,j]]],{i,1,n}] ;
tab=Table[{Red,Circle[{Re[A[[i,i]]],Im[A[[i,i]]]},r[[i]]]},{i,1,n}];
lam=N[Eigenvalues[A]];
tab2=Table[{PointSize[0.02],Blue,Point[{Re[lam[[i]]],Im[lam[[i]]]}]},{i,1,n}];
Show[Graphics[tab],Graphics[tab2],Axes->True,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,PlotRange->{{-6,6},{-6,6}},BaseStyle->20]]

num1=20;
num=30;
n=4;
A=IdentityMatrix[n];
Manipulate[
A[[1,2]]=A[[1,3]]=A[[1,4]]=-(((i/num1))^1);
A[[2,1]]=A[[2,3]]=A[[2,4]]=(((i/num1))^2);
A[[3,1]]=A[[3,2]]=A[[3,4]]=(((i/num1))^3);
A[[4,1]]=A[[4,2]]=A[[4,3]]=Sqrt[(i/num1)];
gcircles[A],{i,0,num}]

The Power Method: Determine the Dominant Eigenvalue

Derivation

The power method is an iterative technique to determine the dominant eigenvalue of a matrix.

We assume an n x n matrix A has the n eigenvalues (Subscript[λ, k])and eigenvectors (v^((k)))for k=1,2,3,...,n. We also assume that the eigenvalues can be orderd in the following manner

|(Subscript[λ, 1])|>|(Subscript[λ, 2])|≥|(Subscript[λ, 3])|≥…≥|(Subscript[λ, n])|≥0,

that is, the largest eigenvalue is of multiplicity one. If the eigenvalues are linearly independent, they will form a basis for vectors in (R^n). Therefore, for a vector x in (R^n) we can write

x = (∑_(j=1)^n)(Subscript[β, j])(v^((j)))
A x = (∑_(j=1)^n)(Subscript[β, j])A (v^((j)))= (∑_(j=1)^n)(Subscript[β, j])(Subscript[λ, j]) (v^((j)))
(A^2) x = (∑_(j=1)^n)(Subscript[β, j])(Subscript[λ, j])A (v^((j)))= (∑_(j=1)^n)(Subscript[β, j])((((Subscript[λ, j])))^2) (v^((j)))

(A^k) x =  (∑_(j=1)^n)(Subscript[β, j])((((Subscript[λ, j])))^k) (v^((j)))
(A^k) x = ((((Subscript[λ, 1])))^k) (∑_(j=1)^n)(Subscript[β, j])(((((Subscript[λ, j])/(Subscript[λ, 1]))))^k) (v^((j)))

Since |(Subscript[λ, 1])|>|(Subscript[λ, j])|, for j=1,2,3,...,n, we have

(lim_(k->∞))(A^k) x =(lim_(k->∞)) ((((Subscript[λ, 1])))^k) (∑_(j=1)^n)(Subscript[β, j])(((((Subscript[λ, j])/(Subscript[λ, 1]))))^k) (v^((j)))= (lim_(k->∞)) ((((Subscript[λ, 1])))^k) (Subscript[β, 1])(v^((1)))
(A^k) x ~ ((((Subscript[λ, 1])))^k) (Subscript[β, 1])(v^((1)))

which converges to zero if |(Subscript[λ, 1])|<1, and diverges if |(Subscript[λ, 1])|>1, if (Subscript[β, 1])≠0. We seek a normalized eigenvector, that is, an eigenvector with ||v||=1. This is to keep terms from becoming too small or too large, and is perfectly valid since eigenvectors are determined up to a constant only (the direction is what is important). These considerations lead to the Power Method to determine the dominant eigenvalue:

Choose (x^((0))) 
 (y^((0)))= (( (x^((0))) )/(||(x^((0)))||))

For  k=1,2,3,...,N

		(x^((k)))=A (x^((k-1))),    (y^((k)))= (( (x^((k))) )/(||(x^((k)))||))

If (y^((N)))~(y^((N-1))), then (y^((N)))is an approximation to an eigenvector of A.

Read eigenvalue from  A (y^((N)))=(Subscript[λ, 1])(y^((N)))

Implementation of Power Method

Here is a sample matrix

n=3;
A={{4,1,1},{0,2,1},{-2,0,9}};

And here is a module to perform the power method (note we need the norm module defined earlier)

powermethod[A_,x0_]:=Module[{x=x0,y,k},
Do[
yold=x/norm[x,2];
x=A.yold;
y = x/norm[x,2];
If[norm[y-yold,2]<(10^(-7)),Return[y]];
Print[y,"   ",norm[y-yold,2],"   ", k]; (*output intermediate steps*)
,{k,1,30}]]

Here we initialize with some vector, and use the power method to determine the largest eqigenvector.

SeedRandom[2]
x=Table[RandomReal[{-3,3}],{i,1,n}]
evec=powermethod[A,x]

We get the estimates of the largest eigenvalue by multiplying the estimated eigenvector by the matrix A.

Table[(A.evec)[[i]]/evec[[i]],{i,1,n}]

Here is Mathematica's answer:

Eigenvalues[N[A]]
Eigenvectors[N[A]]

Power Method: Graphically in (R^2)

n=2;
A={{4,1},{1,2}};
Eigenvalues[A]

powermethod2[A_,x0_]:=Module[{x=x0,y,k},
list={Graphics[{Darker[Green],Arrow[{{0,0}, x/norm[x0,2]}]}]};
Do[
yold=x/norm[x,2];
x=A.yold;
y = x/norm[x,2];
If[norm[y-yold,2]<(10^(-5)),{list = Append[list,Graphics[{Blue,Arrow[{{0,0},y}]}]];Goto[end]}];
list = Append[list,Graphics[{Red,Arrow[{{0,0},y}]}]];
,{k,1,300}];
Label[end];Return[y]]

The blue arrow is the eigenvector associated with the dominant eigenvector, and the green arrow is the initialization used in the power method.

SeedRandom[2]
x=Table[RandomReal[{-1,1}],{i,1,n}]
evec=powermethod2[A,x]
Show[{list,Graphics[Circle[{0,0},1]]},Axes->True,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},AspectRatio->1,PlotRange->{{-1,1},{-1,1}},BaseStyle->20]

evec=powermethod2[A,{-0.4,1.}]
Show[{list,Graphics[Circle[{0,0},1]]},Axes->True,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},AspectRatio->1,PlotRange->{{-1,1},{-1,1}},BaseStyle->20]

Here is Mathematica's result to compare to:

Eigenvectors[N[A]]

Power Method: Graphically in (R^2)

n=2;
A={{4,-1},{2,2}};

evec=powermethod2[A,{-0.4,1.}]
Show[{list,Graphics[Circle[{0,0},1]]},Axes->True,AxesLabel->{"(Subscript[x, 1])","(Subscript[x, 2])"},AspectRatio->1,PlotRange->{{-1,1},{-1,1}},ImageSize->{500,Automatic},BaseStyle->20]

Here is Mathematica's result to compare to:

Eigenvectors[N[A]]
Eigenvalues[A]

We see that the problem is that there is no dominant eigenvalue, so the power method does not converge.

The Inverse Power Method: Determine the Eigenvalue when an Eigenvector is Known

Derivation

The inverse power method is an iterative technique to determine the eigenvalue when an eigenvector is known, or approximately known. Once the eigenvector is known, the eigenvalue can be refined.

We assume the an n x n matrix A has the n eigenvalues (Subscript[λ, k]) and eigenvectors (v^((k))) for k=1,2,3,...,n. Then the n x n matrix (((A-q I))^(-1)) has the n eigenvalues 1/((Subscript[λ, k])-q) and eigenvectors (v^((k)))for k=1,2,3,...,n.

Matrix A: has eigenpairs (Subscript[λ, k]), (v^((k))).

Matrix B=(((A-q I))^(-1)): has eigenpairs (1/((Subscript[λ, k])-q)), (v^((k))).

We then apply the power method to the matrix (((A-q I))^(-1)), where q is chosen to be close to, but not equal to, one of the eigenvalues (Subscript[λ, k]).

Choose (x^((0))) ,q
 (y^((0)))= (( (x^((0))) )/(||(x^((0)))||))

For  k=1,2,3,...,N

		(x^((k)))=(((A-q I))^(-1)) (x^((k-1))),    (y^((k)))= (( (x^((k))) )/(||(x^((k)))||))

If (y^((N)))~(y^((N-1))), then (y^((N)))is an approximation to an eigenvector of (((A-q I))^(-1)), or A.

The Gerschgorin circle method can be used to obtain initial estimates q for the eigenvalues. We could use the techniques we learned earlier to invert a matrix. Typically, the inverse is not found, and instead the solution to the system

(A-q I)(x^((k)))= (x^((k-1)))

is found by using an LU decomposition on the matrix A-q I. Here, I will simply invert the matrix using our MatInv module.

Example

n=3;
A={{4,1,1},{0,2,1},{-2,0,9}};
MatrixForm[A]

Get Initial Approximations using Gerschgorin Circles

gcircles2[A]

Try to find an Eigenvalue around 2:

q=2.1;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

powermethod[A_,x0_]:=Module[{x=x0,y,k},
Do[
yold=x/norm[x,2];
x=A.yold;
y = x/norm[x,2];
If[norm[y-yold,2]<(10^(-7)),Goto[end]];
,{k,1,1000}];
Label[end];Return[y]]

evec=powermethod[B,{1.,1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

We have located an eigenvalue of A at (Subscript[λ, 1])=1.88282, and the associated eigenvector is  (v^((1)))={-0.3826319343811083`,0.9176227407921803`,-0.10752352474141856`} .

Try to find an Eigenvalue around 4:

q=4.0;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

evec=powermethod[B,{1.,1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

We have located an eigenvalue of A at (Subscript[λ, 2])=4.631826645109009`, and the associated eigenvector is  (v^((2)))={0.8980632138691508`,0.1562355043236731`,0.4111847894605775`}.

Try to find an Eigenvalue around 9:

q=9.0;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

evec=powermethod[B,{1.,1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

We have located an eigenvalue of A at (Subscript[λ, 3])=8.485349510620875`, and the associated eigenvector is  (v^((3)))={0.2464738044360242`,0.1476913176132397`,0.9578298066089286`}.

Let's how good our estimates are:

Eigensystem[N[A]]

Pretty darn good indeed!

Example

n=2;
A={{4,-1},{2,2}};
MatrixForm[A]

Get Initial Approximations using Gerschgorin Circles

gcircles2[A]

Try to find an Eigenvalue around 2:

q=2.0;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

evec=powermethod[B,{1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

Well, this didn't work, let's try something else. What happens if we take eigenvalues that are complex instead of real?

Second Attempt  to find an Eigenvalue around 2:

q=2.0+I;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

evec=powermethod[B,{1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

This looks pretty good. We have located an eigenvalue of A at (Subscript[λ, 1])=3+i, and the associated eigenvector is  (v^((1)))={0.5773502470056076`,0.5773502721701564`-0.5773502883931122` i}. Let's check this out for accuracy:

(A.{0.5773502,0.577350272-0.577350288 I})/(3+I)

Yes, we really have found the eigenpairs!

Attempt to find an Eigenvalue around  -i:

q=-I;
B=MatInv[A-q*IdentityMatrix[n]];

We execute the power method on the matrix B:

evec=powermethod[B,{1.,1.}]
Table[(1/((B.evec)[[i]]/evec[[i]]))+q,{i,1,n}]

We have located an eigenvalue of A at (Subscript[λ, 2])=3-i, and the associated eigenvector is  (v^((2)))={0.5773502470056076`,0.5773502721701564`+0.5773502883931122` i}. Since the original matrix was real, this makes sense since the eigenvalues/vecotrs must occur in complex conjugate pairs.

Let's see how good our estimates are:

Eigensystem[N[A]]

The eigenvectors Mathematica reports look different, but they are parallel to the ones we found.