Scipy.linalg.eig() 给出与 GNU Octave 的 eig() 不同的特征向量

Question

我想计算具有λ * M * v = K * v的广义特征值问题的特征值，其中λ为特征值，v为特征向量，M和K为矩阵。假设我们有

K =
   1.8000 + 0.0000i  -1.0970 + 0.9550i
  -1.0970 - 0.9550i   1.8000 + 0.0000i
M = 
   209     0
     0   209

在Octave中，如果我选择[V,D]=eig(K,M)，我得到:

V =
   0.53332 - 0.46429i  -0.53332 + 0.46429i
   0.70711 + 0.00000i   0.70711 + 0.00000i
D =
Diagonal Matrix
   0.34555         0
         0   3.25445

但是，如果我执行scipy. linear。eig(K, b=M)在python中使用Scipy，我得到不同的结果:

>>> import numpy as np
>>> import scipy as sp
>>> import scipy.linalg
>>> K = np.mat([[1.8, -1.097+0.995j], [-1.097-0.955j, 1.8]])
>>> M = np.mat([[209., 0.], [0., 209.]])
>>> M
matrix([[ 209.,    0.],
        [   0.,  209.]])
>>> K
matrix([[ 1.800+0.j   , -1.097+0.955j],
        [-1.097-0.955j,  1.800+0.j   ]])
>>> D, V = sp.linalg.eig(K, b=M)
>>> D
array([ 0.00165333 -1.99202696e-19j,  0.01557155 +0.00000000e+00j])
>>> V
array([[ 0.70710678 +0.00000000e+00j, -0.53332494 +4.64289256e-01j],
       [ 0.53332494 +4.64289256e-01j,  0.70710678 -8.38231384e-18j]])

特征值应为D数组中的特征值。

为什么这两个例子的特征值不同?我误解什么了吗?

编辑:更正打字错误和重新计算。

编辑:我使用了Octave 3.8.2。Mac OS X 10.10.3

Answer 1

我想我明白是怎么回事了。Scipy可能提供了正确的特征值。Octave在其eig()函数中接受第二个矩阵，但没有指定它的作用。Matlab的文档确实说它是一个广义的特征值问题，但在八度曲中添加第二个矩阵似乎对特征值没有影响。这看起来像是八度档的bug。

Answer 2

出于好奇，我用Julia中的eig()例程尝试了相同的计算。(这里，为了清楚起见，我使用符号A和B来表示广义特征值问题。)

A = [ 1.8000 + 0.0000im  -1.0970 + 0.9550im ;
     -1.0970 - 0.9550im   1.8000 + 0.0000im ]
B = [  209.0    0.0 ;
         0.0   209.0 ]
D, V = eig( A, B )
err = A * V - B * V * diagm(D)
@show vecnorm( err )

结果与Scipy基本相同。

D = [0.0016533341538531843,0.015571546228921934]
V[:,1] = Complex{Float64}[0.03689085719366286 - 0.03211555936184871im,
                          0.048911598804451846 + 0.0im]
V[:,2] = Complex{Float64}[-0.03689085719366286 + 0.03211555936184871im,
                           0.048911598804451846 + 0.0im]
vecnorm(err) = 3.768784066493578e-17

另一方面，如果我尝试使用具有非对角元素的B矩阵的Octave，

A = [ 1.8000 + 0.0000i  -1.0970 + 0.9550i ;
     -1.0970 - 0.9550i   1.8000 + 0.0000i ]
B = [ 1.0 5.0 ;
      5.0 3.0 ]
[ V, D ] = eig( A, B );
err = A * V - B * V * D ;
printf( "error = %20.10gn", norm( err ) )

它给出正确的结果:

V =
   0.73147 - 0.19955i  -0.96401 - 0.03599i
   0.74058 + 0.25942i   0.72711 - 0.17996i
D =
   Diagonal Matrix
   0.057841 - 0.000000i                      0
   0                     -0.883750 + 0.000000i
error =      1.814738616e-15

如果我将B设置为diag(209.0, 209.0)(与原始帖子一样)，Octave给出

V =
    0.53332 - 0.46429i  -0.53332 + 0.46429i
    0.70711 + 0.00000i   0.70711 + 0.00000i
D =
    Diagonal Matrix
    0.34555         0
    0               3.25445
error =          676.9262577

实际上，上面的V和D被发现是一个标准特征值问题A V = V D的解。因此，当B是一个对角矩阵时，Octave通过假设单位矩阵来忽略它……

我使用Octave3.4.3(在Linux x86_64上)，结果可能与最新版本不同。