Autovalores y autovectores: Rotación de Jacobi

En la estadística multivariada existe, en algunos casos, la necesidad de determinar los autovalores y autovectores de una matriz real simétrica; cuyos elementos (en ambos casos) son también números reales. Aunque existen varios algoritmos para esto, se va a emplear el procedimiento de la rotación de Jacobi para determinar la matriz de autovectores y el arreglo de autovalores para la siguiente matriz simétrica (grabada como datos.in en JACOBI/datos):

1 2 3
2 1 4
3 4 1

Se elaboró un pequeño programa en C++ que emplea una subrutina con el procedimiento de la rotación de Jacobi (jacobi) y cuyo código es el siguiente:

    #include <iostream>
    #include <fstream>
    #include <cmath>

    using namespace std;

    void jacobi (int n, double **a, double *d, double **v);

    int main() {

       ifstream label1 ("datos//datos.in"); // Abre el archivo de entrada de datos

       // Definición de variables y asignación dinámica de memoria   

       int i, j, n=3;

       double **a, *d, **v;

       a = new double* [n], d = new double [n], v = new double* [n];

       for(j=0; j<n; j++){

          a[j] = new double [n], v[j] = new double [n];
       
       }   

       // Lectura de la matriz (label1 apunta a datos.in en el subdirectorio datos del home de usuario)

            for(i=0; i<3; i++){

                    for(j=0; j<3; j++){

                             label1 >> a[i][j];
                    }

       }

       cout << "Imprime la matriz de origen\n\n";

       // Matriz de origen; a

       for(i=0; i<3; i++){

                    for(j=0; j<3; j++){

                             cout <<  a[i][j] << " ";
                    }

       cout << endl;

            }

       cout << endl;

       jacobi(n, a, d, v); //Destruye la matriz original

       cout << "Imprime la matriz de autovectores\n\n";

       // Matriz de autovectores; v

       for(i=0; i<3; i++){

                    for(j=0; j<3; j++){

                             cout <<  v[i][j] << " ";
                    }

       cout << endl;

            }

       cout << endl;   

       cout << "Imprime el arreglo de autovalores\n\n";

       //Arreglo de autovalores; d

       for(i=0; i<3; i++){

               cout <<  d[i] << " ";
            }

       cout << endl;   

       return 0;

    }

    void jacobi (int n, double **a, double *d, double **v) {

       // Define las variables de tipo entero
    
       int i, j, ip, iq, nrot;

    // Define las variables de tipo doble

       double *b, *z;

       b = new double [n];

       z = new double [n];

       b = new double [n]; z = new double [n];

       double c, g, h, s, sm, t, tau, theta, tresh;

    // Inicializa a la matriz identidad

       for (ip = 0; ip < n; ip++) { 

          for (iq = 0; iq < n; iq++) {

             v[ip][iq] = 0;

          }
        
          v[ip][ip] = 1;

       }

    // Inicializa b y d a la diagonal de a

       for (ip = 0; ip < n; ip++) {   

          b[ip] = a[ip][ip];

          d[ip] = b[ip];

          z[ip] = 0;

       }

       nrot = 0;

       for (i = 0; i < 50; i++) {
    
          sm = 0;
        
          for (ip = 0; ip < n - 1; ip++) {

             for (iq = ip + 1; iq < n; iq++) {

                sm +=fabs(a[ip][iq]);

             }

          } 

          if (sm == 0) break;
        
          if (i < 4)
        
             tresh = 0.2*sm/(n*n);

          else

             tresh = 0.0;

          for (ip =0; ip < n -1; ip++) {

             for (iq = ip + 1; iq < n; iq++) {

                g = 100.0*fabs(a[ip][iq]);

                if(i>4 && (double)(fabs(d[ip])+g) == (double)fabs(d[ip])

                   && (double)(fabs(d[iq])+g) == (double)fabs(d[iq]))

                   a[ip][iq] = 0.0;

                else if (fabs(a[ip][iq]) > tresh) {

                   h = d[iq] - d[ip];

                   if ((double)(fabs(h)+g) == (double)fabs(h))

                      t = (a[ip][iq])/h;   // t = 1/(2theta)

                   else {

                      theta = 0.5*h/(a[ip][iq]);

                      t = 1.0/(fabs(theta)+sqrt(1.0+theta*theta));

                      if(theta < 0.0) t = -t;

                   }

                      c = 1.0/sqrt(1+t*t);

                      s = t*c;

                      tau = s/(1.0+c);

                      h = t*a[ip][iq];

                      z[ip] -=h;

                      z[iq] +=h;

                      d[ip] -=h;

                      d[iq] +=h;

                      a[ip][iq] = 0.0;

    // Varía desde 0 hasta  ip - 1

                   for (j =0; j < ip; j++) { 

                      g = a[j][ip];

                      h = a[j][iq];

                      a[j][ip] = g - s*(h+g*tau);

                      a[j][iq] = h + s*(g-h*tau);

                   } 

    // Varía desde ip+1 hasta  iq - 1

                   for (j =ip+1; j < iq; j++) { 

                      g = a[ip][j];

                      h = a[j][iq];

                      a[ip][j] = g - s*(h+g*tau);

                      a[j][iq] = h + s*(g-h*tau);

                   } 

                   for (j =iq+1; j < n; j++) { 

                      g = a[ip][j];

                      h = a[iq][j];

                      a[ip][j] = g - s*(h+g*tau);

                      a[iq][j] = h + s*(g-h*tau);

                   } 

                  
                   for (j =0; j < n; j++) { 

                      g = v[j][ip];

                      h = v[j][iq];

                      v[j][ip] = g - s*(h+g*tau);

                      v[j][iq] = h + s*(g-h*tau);


                      } 

                   ++(nrot);

                }

             }

          }

             for (ip = 0; ip < n; ip++) {

                b[ip] = b[ip]+z[ip];

                d[ip] = b[ip];

                z[ip] = 0.0;

             }
          
       }   

       delete [] b, z;

    }

y se graba como jacobi.c++.

Se compila en la carpeta JACOBI con:

g++ jacobi.c++ -o jacobi

La ejecución del mismo (./jacobi [Enter]) produce esta salida:

Imprime la matriz de origen

1 2 3
2 1 4
3 4 1

Imprime la matriz de autovectores

0.824038 0.505785 -0.255232
-0.544925 0.584374 -0.601302
-0.154979 0.634577 0.757161

Imprime el arreglo de autovalores

-0.886791 7.07467 -3.18788 

La verificación de que los arreglos anteriores son válidos se puede obtener en una hoja de cálculo. La matriz de autovectores multiplicada por la matriz de autovalores y el resultado, a su vez, multiplicado por la transpuesta de la matriz de autovectores debe reproducir la matriz original.

Entradas relacionadas

Instalando las librerias Lapackpp en Debian

Multiplicación de matrices con Lapack++

Esta entrada fue publicada en Código C++. Guarda el enlace permanente.

4 respuestas a Autovalores y autovectores: Rotación de Jacobi

  1. Está verificado que funciona!

    Saludos

  2. Pingback: Instalando las librerias Lapack++ en Debian |

  3. Pingback: Instalar Lapackpp (librería estándar para álgebra lineal) en Pardus 2009.2 | Pardus Life

  4. Pingback: Análisis de Componentes Principales en Teledetección: Autovalores y autovectores |

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s