計算3×3對稱矩陣譜分解的快速方法

Question

我正在開發一個項目，我基本上在20-100個點的集合上預製PCA數百萬次。目前，我們正在使用一些使用GNU的GSL線性代數包的遺留代碼來對協方差矩陣進行SVD​​。這工作，但非常緩慢。

我想知道是否有任何簡單的方法在3x3對稱矩陣上進行特徵分解，以便我可以將它放在GPU上並讓它並行運行。

由於矩陣本身非常小，我不確定使用什麼樣的算法，因爲它看起來像是爲大型矩陣或數據集設計的。還有對數據集進行直SVD的選擇，但我不確定什麼是最好的選擇。

我不得不承認，我對線性代數並不擅長，特別是在考慮算法優勢時。任何幫助將不勝感激。

（我用C++工作現在）

Answer 1

使用特徵多項式工作，但它傾向於在數值上有些不穩定（或至少不準確）。

用於計算對稱矩陣的特徵系統的標準算法是QR方法。對於3x3矩陣，可以通過構建旋轉外的正交變換並將它們表示爲四元數來實現非常靈活的實現。在C++中實現這個想法（非常簡短！），假設您有一個3x3矩陣和一個Quaternion類，可以找到here。該算法應該非常適合GPU實現，因爲它是迭代的（並且因此是自校正的），當它們可用並且使用相當少量的向量寄存器時可以合理良好地使用快速低維矢量數學基元（所以它允許更活躍的線程）。

Answer 2

我不是恆星的線性代數要麼但由於墨菲說，「當你不知道你在說什麼，一切皆有可能」，有可能是CULA pack might be relevant to your needs。他們做SVD和特徵值分解

Answer 3

大多數方法對較大的矩陣都是有效的。對於小分析來說，分析方法是最快最簡單的，但在某些情況下是不準確的。

Joachim Kopp針對3x3對稱矩陣開發了一個優化的「混合型」method，該矩陣在分析方法上進行中繼，但又回退到QL算法。

可以找到3×3對稱矩陣的其他解決方案here（對稱三對角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; 
    } 
}