June 9, 2009

Jacobi eigenvalue algorithm

The Jacobi eigenvalue algorithm is a numerical procedure for the calculation of all eigenvalues and eigenvectors of a real symmetric matrix.

Each Jacobi rotation can be done in n steps when the pivot element p is known. However the search for p requires inspection of all N ≈ ½ n² off-diag elements. We can reduce this to n steps too if we introduce an additional index array with the property that mi is the index of the largest element in row i, (i = 1, … , n - 1) of the current S. Then (k, l) must be one of the pairs (i,mi) . Since only columns k and l change, only mk and ml must be updated, which again can be done in n steps. Thus each rotation has O(n) cost and one sweep has O(n³) cost which is equivalent to one matrix multiplication. Additionally the mi must be initialized before the process starts, this can be done in n² steps.
Typically the Jacobi method converges within numerical precision after a small number of sweeps. Note that multiple eigenvalues reduce the number of iterations since NS < N. -- WIKIPEDIA

See also:
- Jacobi eigenvalue algorithm, Wikipedia
- A parallel algorithm for the eigenvalues and eigenvectors of a general complex matrix

Here's my test code.

/**
* Jacobi eigenvalue algorithm
*/
public class Jacobi {
public static void main(String[] args) {
double[][] A = new double[][] { { 4, 3, 2, 1, 0, 0 }, { 3, 4, 3, 2, 0, 0 },
{ 2, 3, 4, 3, 0, 0 }, { 1, 2, 3, 4, 0, 0 }, { 0, 0, 0, 0, 1, 0 },
{ 0, 0, 0, 0, 0, 1 } };

double t, c, s;
int p, q, icount, state, size = 6;
double tol = 1.e-5; // the tolerance level of convergence
int icmax = 100; // the maximum iterations number

int[] colRowOfElMax = new int[size], rowOfElMax = new int[1];
double[][] temp = new double[size][size], D = new double[size][size];
double[][] V, diagD;

double[] maxElColRow = new double[size], maxElRow = new double[1];
double[][] dMinusDiagD = new double[size][size], absDminusDiagD = new double[size][size];
double[][] rot = new double[2][2], rotT = new double[2][2];

// makes V into a unit matrix
V = new double[size][size];
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
V[i][j] = 0;
}
V[i][i] = 1.0;
}

D = A; // copies A to D
diagD = diag(D, size);// outputs DiagD=diagonal of D
dMinusDiagD = minus(D, diagD, size); // does D-DiagD
abs(dMinusDiagD, absDminusDiagD, size);// does abs(D-DiagD)
maxMatrix(absDminusDiagD, size, colRowOfElMax, maxElColRow);
maxVector(maxElColRow, size, rowOfElMax, maxElRow);
q = rowOfElMax[0];
p = colRowOfElMax[q];
icount = 0;
state = 1;

// Iterations
while (state == 1 && icount < icmax) {
icount = icount + 1;
if (D[q][q] == D[p][p]) { // check to prevent t from diverging
D[q][q] = D[p][p] + 1.e-10;
}
t = D[p][q] / (D[q][q] - D[p][p]);
c = 1 / Math.sqrt(t * t + 1);
s = c * t;
rot[0][0] = c;
rot[0][1] = s;
rot[1][0] = -s;
rot[1][1] = c;
transpose(rot, rotT, 2);// rotT=transpose(Rot)

for (int i = 0; i < size; i++) {
temp[p][i] = rotT[0][0] * D[p][i] + rotT[0][1] * D[q][i];
temp[q][i] = rotT[1][0] * D[p][i] + rotT[1][1] * D[q][i];
D[p][i] = temp[p][i];
D[q][i] = temp[q][i];
}
for (int i = 0; i < size; i++) {
temp[i][p] = D[i][p] * rot[0][0] + D[i][q] * rot[1][0];
temp[i][q] = D[i][p] * rot[0][1] + D[i][q] * rot[1][1];
D[i][p] = temp[i][p];
D[i][q] = temp[i][q];
}
for (int i = 0; i < size; i++) {
temp[i][p] = V[i][p] * rot[0][0] + V[i][q] * rot[1][0];
temp[i][q] = V[i][p] * rot[0][1] + V[i][q] * rot[1][1];
V[i][p] = temp[i][p];
V[i][q] = temp[i][q];
}

// find the new q, p element array values that need to be changed
diagD = diag(D, size); // outputs diagD=diagonal of D
dMinusDiagD = minus(D, diagD, size); // does D-DiagD
abs(dMinusDiagD, absDminusDiagD, size); // does abs(D-DiagD)
maxMatrix(absDminusDiagD, size, colRowOfElMax, maxElColRow);
maxVector(maxElColRow, size, rowOfElMax, maxElRow);
q = rowOfElMax[0];
p = colRowOfElMax[q];
if (Math.abs(D[p][q]) < tol * Math.sqrt(sumDiagElSq(diagD, size)) / size) {
state = 0;
}
}

// V is the eigen vectors
System.out.println("Jacobi Eigenvalues");
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
System.out.print(diagD[i][j] + "\t");
}
System.out.println(" ");
}
}

/**
* finds the diagonal elements of A and puts them into B
*
* @param A
* @param n
* @return
*/
public static double[][] diag(double A[][], int n) {
double[][] B = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
B[i][j] = 0;
}
B[i][i] = A[i][i];
}

return B;
}

/**
* C = A - B
*
* @param A
* @param B
* @param n
* @return C
*/
public static double[][] minus(double A[][], double B[][], int n) {
double[][] C = new double[n][n];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
C[i][j] = A[i][j] - B[i][j];
}
}
return C;
}

/**
* Finds the absolute value of a matrix
*
* @param A
* @param B
* @param n
*/
public static void abs(double A[][], double B[][], int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
B[i][j] = Math.abs(A[i][j]);
}
}
}

/**
* finds the maximum elements of each column; returns the maximums in Max and
* their array positions in Row
*
* @param A
* @param n
* @param Row
* @param Max
*/
public static void maxMatrix(double A[][], int n, int Row[], double Max[]) {
for (int i = 0; i < n; i++) {
int k = 0;
Max[i] = A[k][i];
Row[i] = k;
for (int j = 0; j < n; j++) {
if (A[j][i] > Max[i]) {
Max[i] = A[j][i];
Row[i] = j;
}
}
k = k + 1;
}
}

/**
* finds the maximum elements of a column of A; returns the maximum of a
* column as Max and its array position as Row
*
* @param A
* @param n
* @param Row
* @param Max
*/
public static void maxVector(double A[], int n, int Row[], double Max[]) {
Max[0] = A[0];
Row[0] = 0;
for (int i = 0; i < n; i++) {
if (A[i] > Max[0]) {
Max[0] = A[i];
Row[0] = i;
}
}
}

/**
* finds the transpose of A and puts it into B
*
* @param A
* @param B
* @param n
*/
public static void transpose(double A[][], double B[][], int n) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
B[i][j] = A[j][i];
}
}
}

/**
* finds the sums of the squared of the diagonal elements of A
*
* @param A
* @param n
* @return
*/
public static double sumDiagElSq(double A[][], int n) {
double sum = 0;
for (int i = 0; i < n; i++) {
sum = A[i][i] * A[i][i] + sum;
}
return sum;
}
}

No comments:

Post a Comment