计算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; 
    } 
}