메모, Jacobi Rotation

1. Similarity of Matrices(닮음 행렬)

이상의 닮음행렬의 고유값과 고유벡터에 관한 정리는 다음과 같이 증명한다.



2. Jacobi Rotation( http://en.wikipedia.org/wiki/Jacobi_rotation )

  

   n x n 실대칭 행렬(Real Symmetric Matrix) A에 대해서, 'Jacobi Rotation(Q(k,l))'을 수행한 행렬 A-hat은 비대각 원소(non-diagonal element)이며 대칭 위치에 있는 두 원소, 즉 k행 l열의 값과 l행 k열의 값을 0으로 만드는 닮음 변환(Similarity Transformation)이다. 다시 말해 A-hat은 A의 닮음 행렬이며, A의 고유값을 유지한다.(Jacobi Rotation에 사용되는 행렬 Q는 직교행렬(Orthogonal)이다.)

/* 1. 이 방법을 비대각 원소들에 대해 반복하면,(이전 연산에 의해 0으로 만들었던 원소에 다른값이 들어가긴 하지만... 수렴해가기 때문인가. 자세한 건 항목 3에서 정리한다.) 최초에 주어진 행렬 A를 대각행렬로 바꿀 수 있으며, 대각행렬의 고유값은 대각행렬의 대각원소들과 정확히 일치하므로, 마지막에 얻은 대각행렬의 대각원소를 고유값으로 사용하면 된다!*/

  이 때, 행렬 Q는 다음과 같다.

(행렬 그림 출처: http://en.wikipedia.org/wiki/Jacobi_rotation)  

즉, nxn 단위 행렬에 대해서 네 위치의 원소를 특정 값으로 수정한 행렬이다. 즉 Q(k, k) = cos(theta), Q(l, l) = cos(theta)와 Q(k, l) = sin(theta), Q(l, k) = -sin(theta)를 nxn 단위행렬에 대해 수행하면 Q를 얻을 수 있다. (그림에 나와있는 기호는 s = sin(theta), c = cos(theta))

/* 2. (나처럼 머리만으로는 이해가 잘 되지 않는 사람들은...) 일반적인 n x n 행렬과 임의의 theta 대해 주어진 연산을 손으로 계산해 보면 A-hat이 A에서 k행, l행, k열, l열만 바뀐 행렬이고, 

   A-hat(k, l) = A-hat(l, h) = (cos(theta)^2 - sin(theta)^2)*A(k,j) + sin(theta)*cos(theta)*(A(k,k)-A(l, l))

이 성립함을 볼 수 있다(A(k,j)와 A(j,k)가 같아야 위 식을 얻을 수 있으므로, A는 조건에서와 같이 대칭행렬...인데 그럼 또 대칭행렬이 아닌 일반적인 행렬에 대해 고유값은 어떻게 구해야하는 걸까? ...). 대칭된 위치에 있는 두 값이 0이 되는 연산은 해야하므로, 식을 0으로 두고 theta를 구하면 우리가 필요로 하는 theta는 다음과 같은 식을 이용해서 행렬 Q를 구성하기 직전에 먼저 구할 수 있다.

 A(k,k) = A(l,l)일 경우에는 (cos(theta)^2 - sin(theta)^2)*A(k,j) = 0  이므로, theta는 pi/4를 쓴다. 

*/

3. Jacobi Eigenvalue Algorithm
( http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm :위키 정의와 설명)
Jacobi Rotation을 수행함으로써, 이전의 행렬과 Frobenius Norm은 같고 행렬의 대각원소를 하나의 벡터로 봤을 때, 이전 행렬의 대각 벡터의 크기(norm)보다 이후 행렬의 대각 벡터의 크기가 더 큰 새로운 행렬 A-hat을 얻을 수 있기 때문에 이를 반복하면, 행렬의 Frobenius Norm와 대각 벡터의 크기가 거의 일치하는(그래서 비대각 원소의 값이 매우 작은) 근사 대각행렬을 얻을 수 있다는 점에서 대칭행렬의 고유값 문제에 접근하는 방식이다.(그래서 Frobenius Norm이 보존되는 것과, 대각 벡터의 크기가 항상 커져서 행렬의 Norm로 수렴한다는 것도 보여야한다... 왜 우리학교는 수치해석만 열리고 수치선대는 강의해주지 않는거죠)

4. Rutishauser's Modification of Jacobi Eigenvalue Algorithm with Threshold Pivoting
( http://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_eigenvalue/jacobi_eigenvalue.html 
   :적용한 소스코드 - 은혜로운 GNU licensed open source code)

////////////////// 공부하면서 주석 달아놓는 중.. .

void jacobi_eigenvalue(int matSize, double** a, double** eigVec, double* eigVal){

    double c;
    double g;
  double h;

double s; 
    double t;

double tau; 
    double term;

double termp;
  double termq;

double theta;

double gapq;

identityFunc(matSize, eigVec); //// 아이겐벡터들을 저장하는 저장공간을 단위행렬로 채움

getDiag(matSize, a, eigVal);   //// 아이겐벨유들을 저장하는 저장공간을 행렬a의 대각원소로 채움

double* bw = new double[matSize];

double* zw = new double[matSize];

for(int i = 0; i < matSize; i++){

  bw[i] = eigVal[i]; //임시 저장공간 bw에 아이겐벨유를 복사

  zw[i] = 0.0; //임시 저장공간 zw는 0.0으로 채움

}

int it_num = 0;

while( it_num < 100 ){//ITERATION STOP CONDITION #1) max iteration을 100회로 지정.

    it_num = it_num + 1;

  double thresh = 0.0;

        //lower triangle 부분을 탐색하면서, thresh에 원소의 제곱값을 저장.
   for(int i = 0; i < matSize; i++){
    for(int j = 0; j < i; j++){
      thresh = thresh + a[i][j] * a[i][j];
    }
  }    

  thresh = sqrt(thresh)/(double)(4*matSize);//(non-diagonal)

        //ITERATION STOP CONDITION #2) non-diagonal norm이 0이라면 그대로 while iteration을 중단.
  if( thresh == 0.0 )

    break;

  for(int p = 0; p < matSize; p++){//Visiting Lower Triangular Part of Original Matrix.

    for(int q = p + 1; q < matSize; q++){

      gapq = 10.0 * fabs(a[q][p]);

      termp = gapq + fabs(eigVal[p]);// matrix[p, p]의 절대값을 더함

      termq = gapq + fabs(eigVal[q]);// matrix[q, q]의 절대값을 더함.

      if (4 < it_num && termp == fabs(eigVal[p]) && termq == fabs(eigVal[q]))
                    //작은 값 교정(?)... 아니면 오차 확인 코드인가?;;
                    // iteration count 가 4보다 크고, 동시에 비대각 원소 matrix[q, p]가 0이라면.

        a[q][p] = 0.0;

      else if(thresh <= fabs(a[q][p])){
                    //thresh

        h = eigVal[q] - eigVal[p];

        term = fabs(h) + gapq;

        if(term == fabs(h))

          t = a[q][p]/h;

        else{

          theta = 0.5 * h / a[q][p];

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

          if (theta < 0.0)

            t = - t;

        }

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

        s = t * c;

        tau = s / ( 1.0 + c );

        h = t * a[q][p];

         zw[p] = zw[p] - h;                 

        zw[q] = zw[q] + h;

        eigVal[p] = eigVal[p] - h;

        eigVal[q] = eigVal[q] + h;

        a[q][p] = 0.0;

        for(int j = 0; j < p; j++){

          g = a[p][j];

          h = a[q][j];

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

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

        }

        for(int j = p+1; j < q; j++){

          g = a[j][p];

          h = a[q][j];

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

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

        }

        for(int j = q+1; j < matSize; j++){

          g = a[j][p];

          h = a[j][q];

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

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

        }

        for(int j = 0; j < matSize; j++){

          g = eigVec[p][j];

          h = eigVec[q][j];

          eigVec[p][j] = g - s * ( h + g * tau );

          eigVec[q][j] = h + s * ( g - h * tau );

        }

      }

    }

  }

///bw업데이트, 아이겐벨유 업데이트, zw초기화.

  for(int i = 0; i < matSize; i++){

    bw[i] = bw[i] + zw[i];

    eigVal[i] = bw[i];

    zw[i] = 0.0;

  }

}

//상삼각 행렬 정리.

for(int j = 0; j < matSize; j++){

  for(int i = 0; i < j; i++){

    a[j][i] = a[i][j];

  }

}

///Delete Dynamic Memory Allocation

delete[] bw;

delete[] zw;

for(int i = 0; i < matSize; i++)

  delete[] a[i];

delete[] a;

return;

};

//// 메트릭스 사이즈를 확인하고 행렬 a의 대각원소를 벡터 v에 저장.

void getDiag(int matSize, double** a, double* v){

for(int i = 0; i < matSize; i++)

v[i] = a[i][i];

  return;

};

//// 메트릭스 사이즈를 확인하고 행렬 a를 단위행렬로 만든다.

void identityFunc(int matSize, double** a)

{

for(int i = 0; i < matSize; i++)

for(int j = 0; j < matSize; j++)

if( i == j ) a[i][j] = 1;

else a[i][j] = 0;

  return;

};

//// 행렬의 Norm을 계산한다.

double norm ( int matSize, int** a ){

double value = 0.0;

for (int i = 0; i < matSize; i++ ){

for (int j = 0; j < matSize; j++ ){

value = value + a[i][j]*a[i][j];

}

}

value = sqrt( value );

return value;

};





아... 그냥 궁금해서 정리하기 시작했는데, 참 이거 왜하고 있는건지. 그냥 소스 가져다 쓸거면서...
게다가 자코비는 공부 순서상 배우는 정도이고 실제로 쓰는 더 좋은 Iteration 방법들이 있는 모양.
결론: 우리학교 학부에도 수치선형대수를 개설해주세요.