### 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.

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;

}

}

## Comments

## Post a Comment