Recently, my “Calculus for Computer Science” teacher assigned the following problem.

So to complete this assignment, I decided to use Actionscript 3 and Flash. A friend of mine and I at Georgia Tech are developing a complex open source matrix library called as3matrix, and we were planning on implementing Jacobi to find eigenvectors of a MxM matrix anyway, so we decided to just apply it to our library.

The Jacobi algorithm is pretty straight forward. It’s impossible to use a formula to find the eigenvalues of a matrix larger than 5×5 because there is no equation solver for equations to that degree. So let’s say you have an 6×6 matrix. To find the eigenvectors, the Jacobi algorithm creates a smaller 2×2 matrix inside that matrix, diagonlizes that, then reapplies it to the matrix.

So what the program I turned in does, is generates a random 5×5 matrix. The program then begins to diagnolize the Matrix using the Jacobi algorithm. Then the program attempts to try diagnolizing the Matrix using the Jacobi method but this time ignoring sorting to solve the 2×2. Obviously, sorting is much faster since it ensures the 2×2 can be diagnolized (since the corner entries will be the largest absolute value of the matrix).

So anyway, here’s how the process works in our as3Matrix library. The jacobi() method computes one interation, while diagonalize() continues the jacobi method until the Off (the sum of the square of off-diagonal elements) is less than 1e-10.

First, take the matrix A. Find the i,j element in the matrix that have the largest absolute value. Create a 2×2 matrix from the i,j elements where a = i,i ; b = i,j ; c = j,i ; d = j,j . This step was probably the hardest part because I kept mixing up the i’s and j’s! Quite annoying when you accidently flip them…

Next, take that 2×2 matrix and diagonalize it. The formula for the eigenvalues that the library uses for 2×2 matrices is:

var L1:Number = ( (a+d)/2 ) + Math.sqrt( 4*b*c + ((a-d)*(a-d)))/2;
var L2:Number = ( (a+d)/2 ) – Math.sqrt( 4*b*c + ((a-d)*(a-d)))/2;

For the eigenvectors, I use a nice trick found by Harvard professor Oliver Knill. I then normalize (which is something Oliver’s page fails to mention) the eigenvectors. Combining the eigenvectors to {u1,u2}, I now have my matrix U.I take that matrix and embed it into the identity (of size of the original, original matrix).  I call that matrix G. Then D is Transpose(G)*A*G.

Then outside of the method I check if the Off(D) is < 1e-10. If so, then I consider the Matrix diagonalized!

Here are the results of Jacobi (with sorting) vs Theoretical Bound and Jacobi (without sorting) vs Theoretical Bound. Since AS3 doesn’t have an LN function, I just used the change of base formula (log(X)/log(2)). I hard coded log(2) to optimize the code.

A couple of random 5×5 matrix sample:

After running 100 random 5×5 symmetric matrics through the Jacobi algorithms, these were the average number of iterations for each:

Average Sorting = 25.11
Average no Sorting = 102.94

Sorting is clearly the best method.

Anyway, you can browse/download the as3matrix library here. Check out the TestJacobi.as in the trunk.