如何求矩阵的最大特征向量?

2024-10-04 05:31:52 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图实现这一点而不成功,我必须这样做,而不使用外部模块NUMPY等。在APP中有3个模块,我正在编码这个,Python和C**,C++,但是没有其他的标准库。在

在一个单独的应用程序中,我使用了numpy的svd,它非常准确地工作。所以我用它来匹配我的结果。我的方法是PCA,到目前为止一切都很好。但是在计算了对称协方差矩阵之后,我不知道如何找到最大的特征向量。在

数据集始终是三维点,x、y和z

vector center;
for(point p in points):
    center += p;
center /= points.count;

sumxx = 0;
sumxy = 0;
sumxz = 0;
sumyy = 0;
sumyz = 0;
sumzz = 0;
for(point p in points):
    vector s = p - center;
    sumxx += s.x * s.x;
    sumxy += s.x * s.y;
    sumxz += s.x * s.z;
    sumyy += s.y * s.y;
    sumyz += s.y * s.z;
    sumzz += s.z * s.z;

matrix3 mat = invert(matrix3(sumxx, sumxy, sumxz, sumxy, sumyy, sumyz, sumxz, sumyz, sumzz));
vector n;
if (determinant(mat) > 0)
    normal = find_largest_eigenvalue

Tags: 模块inforpointspointcentervectormat
2条回答

求特征值有不同的算法。有些是从最小到最大的,比如QR;另一些是从大到小,比如幂迭代或Jacobi-Davidson。在

也许你想换个算法。试试能量法,看看是否有用。在

{a1}

让我们回顾一下您的要求,以澄清:

  1. 求矩阵的特征向量mat
  2. 该特征向量应与矩阵的最大特征值相关联
  3. 该矩阵是principal component analysis的对称协方差矩阵。尤其是对称的。
  4. 您的矩阵是3乘3的正方形,如代码中所示,并在(现在已删除)注释中确认。在

在这些非常具体的情况下,以下答案适用。然而,tmyklebu警告说,对于某些病理性矩阵,这种方法的数值不稳定性,特别是当r接近-1时。在


好吧,让我们从wikipedia's page on Characteristic polynomials开始读一点

In linear algebra, the characteristic polynomial of a square matrix is a polynomial, which is invariant under matrix similarity and has the eigenvalues as roots. det(zI-A)

等等,让我们直接跳到3x3 matrix section in the page on Eigenvalue algorithms。在

If A is a 3×3 matrix, then its characteristic equation can be expressed as: det(zI-A)

后面几行后面是(或多或少)这个伪代码,对于对称矩阵(如果我没弄错的话,您可能有复特征值):

p1 = A(1,2)^2 + A(1,3)^2 + A(2,3)^2
if (p1 == 0) 
   % A is diagonal.
   eig1 = A(1,1)
   eig2 = A(2,2)
   eig3 = A(3,3)
else
   q = (A(1,1) + A(2,2) + A(3,3)) / 3
   p2 = (A(1,1) - q)^2 + (A(2,2) - q)^2 + (A(3,3) - q)^2 + 2 * p1
   p = sqrt(p2 / 6)
   B = (1 / p) * (A - q * I)       % I is the identity matrix
   r = determinant(B) / 2

   % In exact arithmetic for a symmetric matrix  -1 <= r <= 1
   % but computation error can leave it slightly outside this range.
   if (r <= -1) 
      phi = pi / 3
   elseif (r >= 1)
      phi = 0
   else
      phi = acos(r) / 3
   end

   % the eigenvalues satisfy eig3 <= eig2 <= eig1
   eig1 = q + 2 * p * cos(phi)
   eig3 = q + 2 * p * cos(phi + (2*pi/3))
   eig2 = 3 * q - eig1 - eig3     % since trace(A) = eig1 + eig2 + eig3
end

所以第一种情况是max(eig1,eig2,eig3),第二种情况是eig1。让我们称e这个最大的特征值。在

对于特征向量,您现在只需求解(A-e*I)x=0

相关问题 更多 >