여러가지 도구를 이용한 행렬의 고유값(eigenvalue)과 고유벡터(eigenvector) 구하기

여기서는 2x2 행렬

     A = [ 1, 1/2; 1/2, 1/3 ]

의 고유값과 고유벡터를 구하는 여러가지 방법을 소개한다.

정의 (고유값과 고유벡터)
주어진 m by n 행렬 A 에 대하여

   Av = λv   (단, v 는 영아닌 벡터, λ 는 스칼라)

를 만족하는 v 와 λ 가 존재할 때, v 를 A 의 고유벡터, λ 를 A 의 고유값이라고 한다.


* Mathematica를 이용하여 행렬의 고규값과 고유벡터 구하기

* Maxima 를 이용하여 행렬의 고유값과 고유벡터 구하기
** 아래에서 v1 이 고유벡터이고, lambda1 이 고유값이다.
    ( 그런데, v1 이 행벽터로 표현되어 잇다. 
      그러므로 이를 transpose 하여 열벡터로 바꾸어야 한다 .)

** v1 의 transpose 를 취하여 행렬 m 을 좌측에 곱해보기도 하고,
    또 스칼라 lambda1 을 곱해 보기도 한다.
    (그런데, 유리분수식이 간단하지 않다. 유리분수식을 간단히 하려면 어떻게 해야 할까?)

** 유리분수식을 간단하게 하지 않은채 수치적 값을 계산하면 뺄셈한 벡터의 성분이 0이 아니다.
    그러므로 먼저 Maxima 의 ratsimp 함수를 이용하여 유리분수식을 간단히 해야 한다.
    그런 다음 비교하면 고유감과 고유벡터의 정의 대로 계산한 식이 맞게 나온다. 

** m . v1 과 lambda1 * v1 이 같음을 확인하고 그 뺄셈도 계산해 보았다.
   (※ Maxima 의 ratsimp 함수를 써서 유리분수식을 먼자 간단히 하는 과정이 필수적이다. )

* Octave 를 이용하여 행렬의 고유값과 고유벡터 구하기
   ( 아래에서 진하게 된 부분만 입력하면 된다. )

octave.exe:1> A = [1, 1/2; 1/2, 1/3]
A =

   1.00000   0.50000
   0.50000   0.33333

octave.exe:2> [V, LAMBDA] = eig(A)
V =

   0.47186  -0.88167
  -0.88167  -0.47186

LAMBDA =

Diagonal Matrix

   0.065741          0
          0   1.267592

octave.exe:3> v1 = V(:, 1)
v1 =

   0.47186
  -0.88167

octave.exe:4> lambda1 = LAMBDA(1,1)
lambda1 =  0.065741
octave.exe:5> A*v1
ans =

   0.031021
  -0.057963

octave.exe:6> lambda1*v1
ans =

   0.031021
  -0.057963

octave.exe:7> A*v1 - lambda1*v1
ans =

  -2.4286e-017
  -2.0817e-017

* numpy, scipy 를 이용하여 행렬의 고유값과 고유벡터 구하기
   ( 아래에서 진하게 된 부분만 입력하면 된다. )

>>> from numpy import *
>>> from scipy import linalg
>>> A = mat([[1.0, 1/2.0], [1/2.0, 1/3.0]])
>>> print A
[[ 1.          0.5       ]
 [ 0.5         0.33333333]]
>>> A
matrix([[ 1.        ,  0.5       ],
        [ 0.5       ,  0.33333333]])
>>> lam, v = linalg.eig(A)
>>> lam
array([ 1.26759188+0.j,  0.06574145+0.j])
>>> v
array([[ 0.8816746 , -0.47185793],
       [ 0.47185793,  0.8816746 ]])
>>> v2 = v[:,1]
>>> lambda2 = lam[1]
>>> v2
array([-0.47185793,  0.8816746 ])
>>> lambda2
(0.065741454089335127+0j)
>>> v2 = mat(v2)
>>> v2
matrix([[-0.47185793,  0.8816746 ]])
>>> lambda2*v2.T
matrix([[-0.03102063+0.j],
        [ 0.05796257+0.j]])
>>> A*v2.T
matrix([[-0.03102063],
        [ 0.05796257]])
>>> A*v2.T - lambda2*v2.T
matrix([[ -1.73472348e-17+0.j],
        [ -2.08166817e-17+0.j]])

* GSL 라이브러리를 이용하여 행렬의 고유값과 고유벡터를 구하는 C언어 소스

/***********************************************
  Filename: eigenvalues-01.c
 
  Purpose:
              Find eigenvalues and eigenvectors by using the GSL library.
              http://www.gnu.org/software/gsl/
  On Window XP & cygwin:
  On Ubuntu
        Compile: gcc -Wall -I/usr/local/include -c eigenvalues-01.c
               Link: gcc -o eigenvalues-01 -L/usr/local/lib eigenvalues-01.o -lgsl -lgslcblas -lm

 On Mac OS X Lion:
        Compile: gcc -Wall -I/usr/local/include -c eigenvalues-01.c
               Link: gcc -o eigenvalues-01 -L/usr/local/lib eigenvalues-01.o -lgsl -lgslcblas -lm
              
  Date: 2011. 10. 20 (Thu)

  Execute: ./eigenvalues-01

  Execution Result:
A =
[ [   1.00000000  0.50000000 ],
  [   0.50000000  0.33333333 ] ]

v1 =
[  -0.47185793
    0.88167460 ]
lambda1 = 0.0657415
A*v1 =
[  -0.03102063
    0.05796257 ]
lambda1*v1 =
[  -0.03102063
    0.05796257 ]
A*v1 - lambda1*v1 =
[   0.00000000
    0.00000000 ]

v2 =
[   0.88167460
    0.47185793 ]
lambda2 = 1.26759
A*v2 =
[   1.11760356
    0.59812327 ]
lambda2*v2 =
[   1.11760356
    0.59812327 ]
A*v2 - lambda2*v2 =
[   0.00000000
    0.00000000 ]
***************************************************/

#include <stdio.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
    
void printVector(char *msg, gsl_vector *x, const char *fmt) {
   if (msg != NULL && strlen(msg) > 0)
       printf("%s", msg);
    size_t n = x->size;
    size_t i;
    printf("[ ");
    printf (fmt, gsl_vector_get(x, 0));
    for (i = 1; i < n; i++) {
        printf (fmt, gsl_vector_get(x, i));
    }
    printf(" ]\n");
}

void printColVector(char *msg, gsl_vector *x, const char *fmt) {
   if (msg != NULL && strlen(msg) > 0)
       printf("%s", msg);
    size_t n = x->size;
    size_t i;
    printf("[ ");
    printf (fmt, gsl_vector_get(x, 0));
    if (n > 0)
        printf ("\n");
    for (i = 1; i < n; i++) {
        printf ("  ");
        printf (fmt, gsl_vector_get(x, i));
        if (i < n - 1)
            printf ("\n");
    }
    printf(" ]\n");
}

void multiplyVectorByScalar(gsl_vector *c, double x, const gsl_vector *a) {
    size_t n = a->size;
    size_t i;
    double v;
    for (i = 0; i < n; i++) {
        v = x*gsl_vector_get(a, i);
        gsl_vector_set(c, i, v);
    }
}

void addVector(gsl_vector *c, const gsl_vector *a, const gsl_vector *b) {
    size_t n = a->size;
    size_t i;
    double v;
    for (i = 0; i < n; i++) {
        v = gsl_vector_get(a, i) + gsl_vector_get (b, i);
        gsl_vector_set(c, i, v);
    }
}
void subtractVector(gsl_vector *c, const gsl_vector *a, const gsl_vector *b) {
    size_t n = a->size;
    size_t i;
    double v;
    for (i = 0; i < n; i++) {
        v = gsl_vector_get(a, i) - gsl_vector_get (b, i);
        gsl_vector_set(c, i, v);
    }
}

void printMatrix(char *msg, gsl_matrix *a, char *fmt) {
    if (msg != NULL && strlen(msg) > 0)
        printf("%s\n", msg);
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    printf("[ [ ");
    printf (fmt, gsl_matrix_get(a, 0, 0));
    for (j = 1; j < n; j++) {
        printf (fmt, gsl_matrix_get(a, 0, j));
    }
    printf(" ]");
    if (m > 1)
        printf(",\n");
    for (i = 1; i < m; i++) {
        printf("  [ ");
        printf (fmt, gsl_matrix_get(a, i, 0));
        for (j = 1; j < n; j++) {
            printf (fmt, gsl_matrix_get(a, i, j));
    }
    printf(" ]");
    if (i < m -1)
        printf(",\n");
    }
    printf(" ]\n");
}

void product(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2, p = b->size2;
    size_t i, j, k;
    double sum;
    for (i = 0; i < m; i++) {
        for (k = 0; k < p; k++) {
            sum = 0.0;
            for (j = 0; j < n; j++) {
                sum = sum + gsl_matrix_get(a, i, j)*gsl_matrix_get(b, j, k);
            }
            gsl_matrix_set (c, i, k, sum);
        }
    }
}

void productWIthRightVector(gsl_vector *c, const gsl_matrix *a, const gsl_vector *b) {
    size_t m = a->size1, n = b->size;
    size_t i, j;
    double sum;
    for (i = 0; i < m; i++) {
        sum = 0.0;
        for (j = 0; j < n; j++) {
            sum = sum + gsl_matrix_get(a, i, j)*gsl_vector_get(b, j);
        }
        gsl_vector_set(c, i, sum);
    }
}

void productWIthLeftVector(gsl_vector *c, const gsl_vector *a, const gsl_matrix *b) {
    size_t m = a->size, n = b->size2;
    size_t i, j;
    double sum;
    for (j = 0; j < n; j++) {
        sum = 0.0;
        for (i = 0; i < m; i++) {
            sum = sum + gsl_vector_get(a, j)*gsl_matrix_get(b, i, j);
        }
        gsl_vector_set(c, j, sum);
    }
}

void add(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) + gsl_matrix_get (b, i, j);
         gsl_matrix_set (c, i, j, v);
  }
 }
}

void subtract(gsl_matrix *c, const gsl_matrix *a, const gsl_matrix *b) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = gsl_matrix_get (a, i, j) - gsl_matrix_get (b, i, j);
            gsl_matrix_set (c, i, j, v);
        }
    }
}

void multiply(gsl_matrix *c, double x, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    double v;
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            v = x*gsl_matrix_get (a, i, j);
            gsl_matrix_set (c, i, j, v);
        }
    }
}

void set(gsl_matrix *c, const gsl_matrix *a) {
    gsl_matrix_memcpy(c, a);
}

double trace(const gsl_matrix *a) {
 size_t m = a->size1;
 size_t i;
 double sum = 0;
    for (i = 0; i < m; i++) {
        sum = sum + gsl_matrix_get (a, i, i);
 }
 return sum;
}

void transpose(gsl_matrix *c, const gsl_matrix *a) {
    size_t m = a->size1, n = a->size2;
    size_t i, j;
    for (i = 0; i < n; i++) {
        for (j = 0; j < m; j++) {
         gsl_matrix_set (c, i, j, gsl_matrix_get (a, j, i));
  }
 }
}

void getQR(gsl_matrix *Q, gsl_matrix *R, const gsl_matrix *A, const gsl_vector *tau) {
   size_t m = A->size1;
   size_t n = A->size2;
   size_t i, j;
   for (i =0 ; i < m; i++) {
       for (j =0 ; j < n; j++) {
           if (i > j) {
               gsl_matrix_set(R, i, j, 0);
               gsl_matrix_set(Q, i, j, gsl_matrix_get(A, i, j));
           }
           else {
               gsl_matrix_set(R, i, j, gsl_matrix_get(A, i, j));
               gsl_matrix_set(Q, i, j, 0);
           }
       }
   }
}

int main (void) {
    double data[] = {    1.0,  1/2.0,
                       1/2.0,  1/3.0   };
    char msg[40];

    gsl_matrix_view A = gsl_matrix_view_array(data, 2, 2);
    gsl_vector *eval = gsl_vector_alloc(2);
    gsl_matrix *evec = gsl_matrix_alloc(2, 2);
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(4);
    gsl_eigen_symmv(&A.matrix, eval, evec, w);
    gsl_eigen_symmv_free(w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);
    gsl_vector *c = gsl_vector_alloc(2);
    gsl_vector *d = gsl_vector_alloc(2);
    gsl_vector *z = gsl_vector_alloc(2);
    double eval_i;
    gsl_vector_view evec_i;

    printMatrix("A = ", &A.matrix, "%12.8f");
    printf("\n");
    {
         int i;
         for (i = 0; i < 2; i++) {
             eval_i = gsl_vector_get(eval, i);
             evec_i = gsl_matrix_column(evec, i);
             // printf ("v%i = \n", i+1);
             // gsl_vector_fprintf (stdout, &evec_i.vector, "%g");
             sprintf(msg, "v%i = \n", i+1);
             printColVector(msg, &evec_i.vector, "%12.8f");
             printf("lambda%i = %g\n", i+1, eval_i);
   
             productWIthRightVector(c, &A.matrix, &evec_i.vector);
             sprintf(msg, "A*v%i = \n", i+1);
             printColVector(msg, c, "%12.8f");
   
             multiplyVectorByScalar(d, eval_i, &evec_i.vector);
             sprintf(msg, "lambda%i*v%i = \n", i+1, i+1);
             printColVector(msg, d, "%12.8f");
   
             subtractVector(z, c, d);
             sprintf(msg, "A*v%i - lambda%i*v%i = \n", i+1, i+1, i+1);
             printColVector(msg, z, "%12.8f");
             printf("\n");
        }
    }

    
    gsl_vector_free (z);
    gsl_vector_free (d);
    gsl_vector_free (c);
    gsl_vector_free (eval);
    gsl_matrix_free (evec);
    
    return 0;
}