Eigenvalues of a symmetric matrix by Jacobi
Basically
Carl Gustav Jacob Jacobi’s idea was to transform a symmetric matrix A into a diagonal matrix by the use of orthogonal matrixes (rotation matrixes) U in an iterative loop of two steps. At the end of the transformation the elements in the main diagonal are the Eigenvalues of the matrix and the product of all the transformation matrixes contains the Eigenvectors in its columns.

For the transformation Jacobi first multiplies the transposed matrix UT by A and second the resulting matrix A’ by U and gets A’’:

 A’’ = UT * A * U

Such a transformation should eliminate two elements of the matrix A symmetric to the main diagonal. I write “should” because it does not really. It just brings the elements close to 0. But if the same transformation is done several times it will do so.

So the sequence is to start at a21 and a12, determine the transformation matrix and bring these two elements closer to 0. Then go on to a31 and a13 and bring them closer to 0. Do this for all elements except the ones of the main diagonal. Then starting at a21 and a12 again and bring these elements even closer to 0 in a second loop and so on until all the elements are close enough to 0.

The basic rotation matrix U(p, q, φ) used for these transformations looks like:


Matrix

and we have to compute the angle φ to eliminate the element apq with the transformation by this matrix.

Multiplying UT by A gives for the affected elements:

EigenSymm

for j = 1,2,…n


and multiplying the resulting A’ by U

EigenSymm

for i = 1,2,…n


These equations resolved for a’’pq is:

EigenSymm


and as the matrix A is symmetric apq = aqp and we get:

EigenSymm


Now the element a’’pq shall become 0. That means:

EigenSymm


and therefore

EigenSymm


both sides of this equation are divided by 2 for later on and we substitute

EigenSymm


We have θ as a function of app, aqq and apq. And on the other side we have cos(φ) and sin(φ) which we are looking for.


If we use

EigenSymm


and extend the denominator and enumerator by 1/ cos(φ)2, we get:

EigenSymm


and with t = tan(φ) = sin(φ)/ cos(φ)

EigenSymm


This converted a bit

EigenSymm


and extended by + θ2 – θ2:

EigenSymm


and this can be written as

EigenSymm


or

EigenSymm


or as well as

EigenSymm


From these two results we use to smaller one and say

EigenSymm


for θ <> 0

and

t = 1

for θ = 0



and finally for the angle we are looking for:

EigenSymm


With this we have the elements of the rotation matrix U.


I put this into a function CalcRMatrix() as follows:



public void CalcRMatrix(double[,] a, int p, int q, int size)
{
     int i, j;
     double phi;
     double t;
     phi = (a[q, q] - a[p, p]) / 2 / a[q, p];
     if (phi != 0)
     {
         if (phi > 0)
              t = 1/(phi + Math.Sqrt(phi * phi + 1.0));
         else
              t = 1/(phi - Math.Sqrt(phi * phi + 1.0));
     }
     else
         t = 1.0;
     for (i = 0; i < size; i++)
     {
         for (j = 0; j < size; j++)
         {
              r[i, j] = 0.0;
         }
         r[i, i] = 1.0;
     }
     r[p, p] = 1 / Math.Sqrt(t * t + 1.0);
     r[p, q] = t * r[p, p];
     r[q, q] = r[p, p];
     r[q, p] = -r[p, q];
}
 


This function calculates the sin(φ) and cos(φ) values to eliminate the element apq and writes them into the matrix R. Now we first have to multiply the transposed of this rotation matrix RT by A and the multiply the resulting A’ by R. For these multiplications there are some interesting details to consider:

The rotation matrix has the following form:

EigenSymm


For the transposed matrix only the sign of the two sinus values switches or we can switch the indexing for x and y in the multiplication. That’s what I did in my function MultTranspMatrixLeft().

In the rotation matrix all the rows except p and q only have one element not zero: The diagonal element that means if we multiply this rotation matrix by another matrix from the left side only the rows p and q will affect the other matrix. All the elements in the other rows will remain the same.

If the rotation matrix is multiplied by another matrix from the right side, only the columns p and q will affect the other matrix.

That means we just need to handle these two rows or columns for a multiplication as:



public double[,] MultTranspMatrixLeft(double[,] a, double[,] m, int p, int q, int size)
{
     int j, k;
     double[,] al = new double[2, size];
     for (j = 0; j < size; j++)
     {
         al[0, j] = 0.0;
         al[1, j] = 0.0;
         for (k = 0; k < size; k++)
         {
              al[0, j] = al[0, j] + (m[k, p] * a[k, j]);
              al[1, j] = al[1, j] + (m[k, q] * a[k, j]);
         }
     }
     for (j = 0; j < size; j++)
     {
         a[p, j] = al[0, j];
         a[q, j] = al[1, j];
     }
     return a;
}



This function carries out the multiplication of the transposed matrix m from the left side by the matrix a and returns the result in a.



public double[,] MultMatrixRight(double[,] a, double[,] m, int p, int q, int size)
{
     int i, j, k;
     double[,] al = new double[size, 2];
     for (i = 0; i < size; i++)
     {
         al[i, 0] = 0.0;
         al[i, 1] = 0.0;
         for (k = 0; k < size; k++)
         {
              al[i, 0] = al[i, 0] + (a[i, k] * m[k, p]);
              al[i, 1] = al[i, 1] + (a[i, k] * m[k, q]);
         }
     }
     for (i = 0; i < size; i++)
     {
         a[i, p] = al[i, 0];
         a[i, q] = al[i, 1];
     }
     return a;
}



Multiplies m from the right side by matrix a and returns the result in matrix a.


With these functions one iteration can be run as



for (i = 0; i < maxOrder; i++)
{
     for (k = i + 1; k < maxOrder; k++)
     {
         if (a[k, i] != 0.0)
         {
              CalcRMatrix(a, i, k, maxOrder);
              a = MultTranspMatrixLeft(a, r, i, k, maxOrder);
              a = MultMatrixRight(a, r, i, k, maxOrder);
         }
     }
}



One such iteration brings all elements except the ones in the main diagonal a bit closer to 0 and we have to repeat it until all elements are close enough to 0 that we can regard them as 0. To check whether the elements are close enough to 0 or not it’s enough just to look at the elements right below the main diagonal. If they can be considered as 0 that’s o.k. At the end of these loops we have a diagonal matrix and the elements in the main diagonal are the Eigenvalues.

For the eigenvectors we just have to multiply a unit matrix by all the rotation matrixes and at the end all the eigenvectors are the columns of the resulting matrix. I use the matrix v for this, initialise it as unit matrix and multiply it as


v = MultMatrixRight(v, r, i, k, maxOrder);


in each loop. The entire calculation I put into the function RunJacobi().

With a sample matrix of

EigenSymm

I get the Eigenvalues

EigenSymm


And the eigenvectors to the eigenvalues in the corresponding column

EigenSymm

As symmetric matrixes only have real eigenvalues the eigenvectors have only real elements to.$


But remark: This is valid for symmetric matrixes only.



C# Demo Project Eigenvalues Jakobi
  • Eigenvalues_Sym.zip