Eigenvalue computation by vector iteration

Vector iteration

The vector iteration is the simplest way to compute one single, real Eigenvalue. The outcome of the vector iteration is an eigenvector to which the eigenvalue can be computed. Usually it finds the eigenvector with the eigenvalue that has the biggest absolute value (but that’s not for sure :-))

Starting from the basic formulation for eigenvalues and their eigenvectors (see Eigenvalues basics)

EigenvalVect

With A as a n*n matrix, x the eigenvector and λ the eigenvalue, the vector iteration (in the common literature) builds a formulation for an iteration that looks a bit strange. Therefore I want to start with a little different formulation. With the above formula an iteration formula would look like:

EigenvalVect

this expression converges to an eigenvector in xk+1 if it is repeated enough times. The only disadvantage is that the eigenvector can run into an overflow or underflow.


If the eigenvalue is taken to the right side the formulation for xk+1 is:

EigenvalVect
and λ can be calculated out of the formulation

EigenvalVect

if we just use one row of A we get for λ

EigenvalVect

In my implementation I use all rows of A, calculate n λ’s and finally take the mean of them all. For this I use the array y.


The formulation for the iteration means there is a multiplication of the matrix A with the vector xk. The result of this is another vector xk+1 and this one is divided by λ and then it is used for the next iteration. And λ is calculated as mentioned above and if it does not change anymore, the vector iteration has converged and we have got the eigenvector in x and the eigenvalue in λ.


The iteration starts with a random eigenvector and looks like:


public void RunVector(int maxOrder)
{
int i, k, count = 0;
double[] x2 = new double[maxOrder];
double[] y = new double[maxOrder];
double delta = 10;
double oldZ = 1E30;
 
while ((delta > 1.0E-50) && (count < 150))
{
count++;
// multiplication A * xn
x2 = MultMatrixRight(a, x1, maxOrder);
 
for (i = 0; i < maxOrder; i++)
{
y[i] = 0;
for (k = 0; k < maxOrder; k++)
{
y[i] = y[i] + a[k, i] * x1[k];
}
}
 
// calculate the eigenvalue
eigenval = 0;
for (i = 0; i < maxOrder; i++)
{
eigenval = eigenval + y[i] / x1[i];
}
eigenval = eigenval / maxOrder;
 
// set xn+1 for the next iteration
for (i = 0; i < maxOrder; i++)
{
x1[i] = x2[i] / eigenval;
}
 
// check the delta between the last iteration and the current
delta = Math.Abs(oldZ - eigenval);
oldZ = eigenval;
}
}




This function gets an initial eigenvector in the array x1 and the array to compute the eigenvalues from in the array a. It uses the function



public double[] MultMatrixRight(double[,] a, double[] z, int size)
{
int i, k;
double[] al = new double[size];
for (i = 0; i < size; i++)
{
al[i] = 0.0;
for (k = 0; k < size; k++)
{
al[i] = al[i] + (a[i, k] * z[k]);
}
}
return al;
}


For the matrix multiplication.



With a sample matrix of

EigenvalVect


I get the Eigenvalue 23.527


And the eigenvector

EigenvalVect



To avoid an overflow or underflow of the eigenvector it is principally a good idea to normalize the eigenvector in each iteration by dividing it by its length. That would be:

EigenvalVect



and

EigenvalVect



Whereas ||wk+1|| is the length of the vector wk+1 But instead of dividing by λ and additionally by ||wk+1|| we can say it’s as good just to divide by the length of the eigenvector. That’s enough and the algorithm converges even better.

EigenvalVect



Or with the auxiliary variable wk+1:

EigenvalVect



and

EigenvalVect


And that’s what the current literature describes as the vector iteration.


According to what’s written in the literature the division by ||wk+1|| is a normalization of wk+1 to avoid an overflow or underflow during iteration.

I say this ||wk+1|| contains λ (even if that’s not too obvious) and the algorithm would not work without this division. ||wk+1|| is used in place of λ * ||wk+1||  to simplify the algorithm and to improve its ability to converge.

Basically the above implementation is rather complicate and does not make too much sense. But I think it is quite interesting for the understanding of the vector iteration as this is straight according to the formulation and the descriptions that can be found in the net are not too clear about this derivation.


Now the real implementation uses just ||A * xk|| for the normalisation of the eigenvector. That’s much more elegant :-) .


And that gives the following procedure:


public void RunVector(int maxOrder)
{
int i, k, count = 0;
double[] x2 = new double[maxOrder];
double[] y = new double[maxOrder];
double delta = 10;
double oldZ = 1E30;
double lenW = 0;
 
while ((delta > 1.0E-50) && (count < 50))
{
count++;
x2 = MultMatrixRight(a, x1, maxOrder);
lenW = 0.0;
 
// calculate the length for normalisation
for (i = 0; i < maxOrder; i++)
{
lenW = lenW + x2[i] * x2[i];
}
lenW = Math.Sqrt(lenW);
 
// normalisation of the eigenvector
for (i = 0; i < maxOrder; i++)
{
x1[i] = x2[i] / lenW;
}
delta = Math.Abs(oldZ - lenW);
oldZ = lenW;
}
 
// calculate the eigenvalue
for (i = 0; i < maxOrder; i++)
{
y[i] = 0;
for (k = 0; k < maxOrder; k++)
{
y[i] = y[i] + a[k, i] * x1[k];
}
}
eigenval = 0;
for (i = 0; i < maxOrder; i++)
{
eigenval = eigenval + y[i] / x1[i];
}
eigenval = eigenval / maxOrder;
}




To check if the calculation is accurate enough I use the length of x2 here. If that does not change anymore the computation is finished.



With a sample matrix of

EigenvalVect


I get the Eigenvalue 23.527

And the eigenvector

EigenvalVect


I use a symmetric matrix to make sure I get a real eigenvalue as this implementation of the Vector Iteration can only find real eigenvalues.

With these sample values here it finds the eigenvalue with the biggest absolute value. But if the matrix is not symmetric that can be different :-) .



Inverse Vector iteration

The vector iteration usually finds the eigenvalue that has the biggest absolute value. But quite often it is more interesting to have the smallest eigenvalue. Therefore the inverse vector iteration should be usable. It runs the opposite way. They say the formulation for the eigenvalues can as well be written as:

EigenvalVect
If λ is an eigenvalue of A, 1/λ is an eigenvalue of A-1.

And with this it finds the biggest inverse value of an eigenvalue and that would be the smallest if it is not inverse. I tried that and it did not work :-)



The above sample has the eigenvalues
 -1,161; 1,173 ; 6,461 ; 23,527

The first two values have absolute values very close to each other. That seems to be too close for the vector iteration and so the algorithm did not work with the inverse matrix. I think that’s the weakness of the vector iteration. It does not converge too well if the absolute values of the eigenvalues are too close to each other.
(but anyhow for the inverse of a matrix see Inverse matrix)


The behaviour of the inverse vector iteration can be improved if a so called shift value is used. If the eigenvalues of the matrix A are “shifted” by a value close to an eigenvalue, it converges quite well. To shift means to subtract a shift-value µ from all main diagonal elements of the matrix A. In a formula this is written as µ*I (I is the identity matrix) subtracted from A.

EigenvalVect
There are two ways to implement this formulation:

We can build the inverse matrix of A - µ*I and multiply this by xk.

But the second approach is more elegant. It says

EigenvalVect

Is the same as solving the matrix equation

EigenvalVect

For x.

And that can be done by a Gaussian elimination algorithm.

So, with x set to yk+1 and y set to xk, the algorithm becomes:

To solve the matrix equation

EigenvalVect

Normalize yk+1 and set this as xk+1

EigenvalVect

The eigenvalue is calculated according to (same as above)

EigenvalVect

And repeat the iteration as long as the eigenvalue changes.


I implemented this into this C# function:


private void RunVector(TMatrix m)
{
int i, k;
double len = 0.0;
double[] y = new double[MaxOrder];
gauss = new TGauss(MaxOrder, m);
for (i = 0; i < MaxOrder; i++)
{
for (k = 0; k < MaxOrder; k++)
{
m.a[k, i] = ori.a[k, i];
}
m.a[i, i] = m.a[i, i] - shift;
m.x[i] = 0;
}
gauss.Eliminate();
gauss.Solve();
 
// calculate the length of x
len = 0;
for (i = 0; i < MaxOrder; i++)
{
len = len + m.x[i] * m.x[i];
}
len = Math.Sqrt(len);
// narmalize the eigenvector and set y for the next loop
for (i = 0; i < MaxOrder; i++)
{
eigenVect[i] = m.x[i] / len; // eigenvector to u
m.y[i] = eigenVect[i]; // set y for the next iteration
}
 
// calculate the eigenvalue
for (i = 0; i < MaxOrder; i++)
{
y[i] = 0;
for (k = 0; k < MaxOrder; k++)
{
y[i] = y[i] + ori.a[k, i] * eigenVect[k];
}
}
eigenvalue = 0;
for (i = 0; i < MaxOrder; i++)
{
eigenvalue = eigenvalue + y[i] / eigenVect[i];
}
eigenvalue = eigenvalue / MaxOrder;
}



Here I use y of my matrix class instead of x1 and solve the matrix equation A*x = y by a Gaussian elimination algorithm. The matrix class contains A, X and Y as arrays.


With this implementation and the same sample as above

EigenvalVect

and a starting value µ = -2

I get the Eigenvalue -1.16095

and the eigenvector

EigenvalVect

and with a starting value µ = 2

I get the Eigenvalue 1.17305

with the eigenvector

EigenvalVect


It is basically possible to find every eigenvalue as long as the starting value of µ is close enough to it. But this is in my opinion the big disadvantage of this algorithm: You need to know the approximate value of the eigenvalue you are looking for. Therefore I think it might be smarter to use a QR decomposition algorithm and compute all eigenvalues at once :-) (see Eigenvalues QR decomposition and Eigenvalues of a symmetric matrix by Jacobi)


C# Demo Projects Vector iteration
  • Eigenvalues_Vectoriteration.zip
  • Eigenvalues_Vector_invers.zip