3x3 대칭 매트릭스 스펙트럼 분해를 계산하는 빠른 방법

Question

저는 기본적으로 20-100 점의 세트에서 PCA 수백만 번을 수행하는 프로젝트에서 작업하고 있습니다. 현재 우리는 공분산 행렬에서 SVD를 수행하기 위해 GNU의 GSL 선형 대수 팩을 사용하는 일부 레거시 코드를 사용하고 있습니다. 이것은 효과가 있지만 매우 느립니다.

3x3 대칭 매트릭스에서 고유 분해를 수행하는 간단한 방법이 있는지 궁금해서 GPU에 넣고 병렬로 실행할 수 있습니다.

행렬 자체가 너무 작기 때문에 큰 행렬이나 데이터 집합을 위해 설계된 것처럼 보이기 때문에 어떤 알고리즘을 사용해야할지 확신하지 못했습니다. 또한 데이터 세트에서 직선 SVD를 수행 할 수있는 선택이 있지만 최상의 옵션이 무엇인지 확신 할 수 없습니다.

선형 대수학에서는 알고리즘의 이점을 고려할 때 별의 여지가 없습니다. 어떤 도움이라도 대단히 감사하겠습니다.

Answer 1

특성 다항식을 사용하면 효과적이지만 다소 수치 적으로 불안정한 경향이 있습니다 (또는 적어도 부정확합니다).

대칭 행렬에 대한 고유 시스템을 계산하는 표준 알고리즘이 QR 방법입니다. 3x3 행렬의 경우 회전에서 직교 변환을 만들고 Quaternion으로 표현하면 매우 매끄러운 구현이 가능합니다. 여러분이 3x3 행렬과 Quaternion 클래스를 가지고 있다고 가정 할 때 C++에서이 아이디어를 (아주 짧게!) 구현하면 here을 찾을 수 있습니다. 알고리즘은 반복적 (따라서 자체 수정) 및 사용 가능한 경우 고속 저 차원 벡터 수학 원시 타입을 합리적으로 효과적으로 사용할 수 있고 상당히 적은 수의 벡터 레지스터를 사용하기 때문에 GPU 구현에 상당히 적합해야합니다 (따라서 더 많은 활성 스레드를 허용).

Answer 2

내가 선형 대수학에서 별 아니다 (나는 지금 ++ C에서 일하고 있어요) 중 하나를하지만 머피는 "당신이 무슨 말을하는지 모를 때, 모든 것이 가능하다"고 언급 한 이후, 가능하면 CULA pack might be relevant to your needs입니다. 그들은 SVD와 고유 값 분해를 수행합니다.

Answer 3

대부분의 방법은 큰 행렬에 효율적입니다. 작은 경우 분석 방법이 가장 빠르고 간단하지만 어떤 경우에는 정확하지 않습니다.

요아킴 코프 (Yachim Kopp)는 3x3 대칭 매트릭스에 최적화 된 "하이브리드"를 개발했습니다.이 매트릭스는 분석적 mathod에서 중계되지만 QL 알고리즘으로 되돌아갑니다.

3x3 대칭 행렬에 대한 다른 솔루션은 here (대칭형 3 차원 QL 알고리즘)입니다.

Answer 4

// Slightly modified version of Stan Melax's code for 3x3 matrix diagonalization (Thanks Stan!) 
// source: http://www.melax.com/diag.html?attredirects=0 
void Diagonalize(const Real (&A)[3][3], Real (&Q)[3][3], Real (&D)[3][3]) 
{ 
    // A must be a symmetric matrix. 
    // returns Q and D such that 
    // Diagonal matrix D = QT * A * Q; and A = Q*D*QT 
    const int maxsteps=24; // certainly wont need that many. 
    int k0, k1, k2; 
    Real o[3], m[3]; 
    Real q [4] = {0.0,0.0,0.0,1.0}; 
    Real jr[4]; 
    Real sqw, sqx, sqy, sqz; 
    Real tmp1, tmp2, mq; 
    Real AQ[3][3]; 
    Real thet, sgn, t, c; 
    for(int i=0;i < maxsteps;++i) 
    { 
     // quat to matrix 
     sqx  = q[0]*q[0]; 
     sqy  = q[1]*q[1]; 
     sqz  = q[2]*q[2]; 
     sqw  = q[3]*q[3]; 
     Q[0][0] = (sqx - sqy - sqz + sqw); 
     Q[1][1] = (-sqx + sqy - sqz + sqw); 
     Q[2][2] = (-sqx - sqy + sqz + sqw); 
     tmp1  = q[0]*q[1]; 
     tmp2  = q[2]*q[3]; 
     Q[1][0] = 2.0 * (tmp1 + tmp2); 
     Q[0][1] = 2.0 * (tmp1 - tmp2); 
     tmp1  = q[0]*q[2]; 
     tmp2  = q[1]*q[3]; 
     Q[2][0] = 2.0 * (tmp1 - tmp2); 
     Q[0][2] = 2.0 * (tmp1 + tmp2); 
     tmp1  = q[1]*q[2]; 
     tmp2  = q[0]*q[3]; 
     Q[2][1] = 2.0 * (tmp1 + tmp2); 
     Q[1][2] = 2.0 * (tmp1 - tmp2); 

     // AQ = A * Q 
     AQ[0][0] = Q[0][0]*A[0][0]+Q[1][0]*A[0][1]+Q[2][0]*A[0][2]; 
     AQ[0][1] = Q[0][1]*A[0][0]+Q[1][1]*A[0][1]+Q[2][1]*A[0][2]; 
     AQ[0][2] = Q[0][2]*A[0][0]+Q[1][2]*A[0][1]+Q[2][2]*A[0][2]; 
     AQ[1][0] = Q[0][0]*A[0][1]+Q[1][0]*A[1][1]+Q[2][0]*A[1][2]; 
     AQ[1][1] = Q[0][1]*A[0][1]+Q[1][1]*A[1][1]+Q[2][1]*A[1][2]; 
     AQ[1][2] = Q[0][2]*A[0][1]+Q[1][2]*A[1][1]+Q[2][2]*A[1][2]; 
     AQ[2][0] = Q[0][0]*A[0][2]+Q[1][0]*A[1][2]+Q[2][0]*A[2][2]; 
     AQ[2][1] = Q[0][1]*A[0][2]+Q[1][1]*A[1][2]+Q[2][1]*A[2][2]; 
     AQ[2][2] = Q[0][2]*A[0][2]+Q[1][2]*A[1][2]+Q[2][2]*A[2][2]; 
     // D = Qt * AQ 
     D[0][0] = AQ[0][0]*Q[0][0]+AQ[1][0]*Q[1][0]+AQ[2][0]*Q[2][0]; 
     D[0][1] = AQ[0][0]*Q[0][1]+AQ[1][0]*Q[1][1]+AQ[2][0]*Q[2][1]; 
     D[0][2] = AQ[0][0]*Q[0][2]+AQ[1][0]*Q[1][2]+AQ[2][0]*Q[2][2]; 
     D[1][0] = AQ[0][1]*Q[0][0]+AQ[1][1]*Q[1][0]+AQ[2][1]*Q[2][0]; 
     D[1][1] = AQ[0][1]*Q[0][1]+AQ[1][1]*Q[1][1]+AQ[2][1]*Q[2][1]; 
     D[1][2] = AQ[0][1]*Q[0][2]+AQ[1][1]*Q[1][2]+AQ[2][1]*Q[2][2]; 
     D[2][0] = AQ[0][2]*Q[0][0]+AQ[1][2]*Q[1][0]+AQ[2][2]*Q[2][0]; 
     D[2][1] = AQ[0][2]*Q[0][1]+AQ[1][2]*Q[1][1]+AQ[2][2]*Q[2][1]; 
     D[2][2] = AQ[0][2]*Q[0][2]+AQ[1][2]*Q[1][2]+AQ[2][2]*Q[2][2]; 
     o[0] = D[1][2]; 
     o[1] = D[0][2]; 
     o[2] = D[0][1]; 
     m[0] = fabs(o[0]); 
     m[1] = fabs(o[1]); 
     m[2] = fabs(o[2]); 

     k0  = (m[0] > m[1] && m[0] > m[2])?0: (m[1] > m[2])? 1 : 2; // index of largest element of offdiag 
     k1  = (k0+1)%3; 
     k2  = (k0+2)%3; 
     if (o[k0]==0.0) 
     { 
      break; // diagonal already 
     } 
     thet = (D[k2][k2]-D[k1][k1])/(2.0*o[k0]); 
     sgn  = (thet > 0.0)?1.0:-1.0; 
     thet *= sgn; // make it positive 
     t  = sgn /(thet +((thet < 1.E6)?sqrt(thet*thet+1.0):thet)) ; // sign(T)/(|T|+sqrt(T^2+1)) 
     c  = 1.0/sqrt(t*t+1.0); // c= 1/(t^2+1) , t=s/c 
     if(c==1.0) 
     { 
      break; // no room for improvement - reached machine precision. 
     } 
     jr[0 ] = jr[1] = jr[2] = jr[3] = 0.0; 
     jr[k0] = sgn*sqrt((1.0-c)/2.0); // using 1/2 angle identity sin(a/2) = sqrt((1-cos(a))/2) 
     jr[k0] *= -1.0; // since our quat-to-matrix convention was for v*M instead of M*v 
     jr[3 ] = sqrt(1.0f - jr[k0] * jr[k0]); 
     if(jr[3]==1.0) 
     { 
      break; // reached limits of floating point precision 
     } 
     q[0] = (q[3]*jr[0] + q[0]*jr[3] + q[1]*jr[2] - q[2]*jr[1]); 
     q[1] = (q[3]*jr[1] - q[0]*jr[2] + q[1]*jr[3] + q[2]*jr[0]); 
     q[2] = (q[3]*jr[2] + q[0]*jr[1] - q[1]*jr[0] + q[2]*jr[3]); 
     q[3] = (q[3]*jr[3] - q[0]*jr[0] - q[1]*jr[1] - q[2]*jr[2]); 
     mq  = sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]); 
     q[0] /= mq; 
     q[1] /= mq; 
     q[2] /= mq; 
     q[3] /= mq; 
    } 
}