对于 Hermitian 矩阵，numpy.linalg.eig 给出的结果与 numpy.linalg.eigh 不同

Posted 2023-02-16

【中文标题】对于 Hermitian 矩阵，numpy.linalg.eig 给出的结果与 numpy.linalg.eigh 不同

【英文标题】：Numpy.linalg.eig is giving different results than numpy.linalg.eigh for Hermitian matrices

【发布时间】：2021-04-01 11:33:50

【问题描述】：

我有一个厄米特矩阵（特别是哈密顿矩阵）。尽管单个特征向量的相位可以是任意的，但我正在计算的数量是物理的（我将代码减少了一点，只保留了可重现的部分）。 eig 和 eigh 给出的结果非常不同。

import numpy as np
import numpy.linalg as nlg
import matplotlib.pyplot as plt
   
def Ham(Ny, Nx, t, phi):
    h = np.zeros((Ny,Ny), dtype=complex)
    for ii in range(Ny-1):
        h[ii+1,ii] = t
    h[Ny-1,0] = t    
    h=h+np.transpose(np.conj(h))
    
    u = np.zeros((Ny,Ny), dtype=complex)
    for ii in range(Ny):
        u[ii,ii] = -t*np.exp(-2*np.pi*1j*phi*ii)
    u = u + 1e-10*np.eye(Ny)
    
    H = np.kron(np.eye(Nx,dtype=int),h) + np.kron(np.diag(np.ones(Nx-1), 1),u) + np.kron(np.diag(np.ones(Nx-1), -1),np.transpose(np.conj(u)))
    H[0:Ny,Ny*(Nx-1):Ny*Nx] = np.transpose(np.conj(u))    
    H[Ny*(Nx-1):Ny*Nx,0:Ny] = u                           

    x=[]; y=[]; 
    for jj in range (1,Nx+1):
        for ii in range (1,Ny+1):
            x.append(jj); y.append(ii)        
    x = np.asarray(x)
    y = np.asarray(y)
    return H, x, y

def C_num(Nx, Ny,  E, t, phi):
    H, x, y = Ham(Ny, Nx, t, phi)
    
    ifhermitian = np.allclose(H, np.transpose(np.conj(H)), rtol=1e-5, atol=1e-8)
    assert ifhermitian == True
    
    Hp = H
    V,wf = nlg.eigh(Hp)     ##Check. eig gives different result
    idx = np.argsort(np.real(V))
    wf = wf[:, idx]

    normmat = wf*np.conj(wf)
    norm = np.sqrt(np.sum(normmat, axis=0))
    wf = wf/(norm*np.sqrt(len(H)))
    
    wf = wf[:, V<=E]  ##Chose a subset of eigenvectors              
    
    V01 = wf*np.exp(1j*x)[:,None]; V12 = wf*np.exp(1j*y)[:,None]
    V23 = wf*np.exp(1j*x)[:,None]; V30 = wf*np.exp(1j*y)[:,None]
    
    wff = np.transpose(np.conj(wf))
    C01 = np.dot(wff,V01); C12 = np.dot(wff,V12); C23 = np.dot(wff,V23); C30 = np.dot(wff,V30)

    F = nlg.multi_dot([C01,C12,C23,C30])

    ifhermitian = np.allclose(F, np.transpose(np.conj(F)), rtol=1e-5, atol=1e-8)
    assert ifhermitian == True

    evals, efuns = nlg.eig(F)    ##Check eig gives different result
    C = (1/(2*np.pi))*np.sum(np.angle(evals));
    
    return  C

C = C_num(16, 16, 0, 1, 1/8)
print(C)
            













  

将 nlg.eigh 都更改为 nlg.eig，甚至只更改最后一个，得到截然不同的结果。

【参考方案1】：

正如我提到的elsewhere，特征值和特征向量不是唯一的。

唯一正确的是，对于每个特征值 $A v = lambda v$，eig 和 eigh 返回的两个矩阵描述了这些解决方案，eig 不精确但近似的结果是很自然的。

您可以看到两种解决方案都会以不同的方式三角化您的矩阵

H, x, y = Ham(16, 16, 1, 1./8)
D, V = nlg.eig(H)
Dh, Vh = nlg.eigh(H)

然后

import matplotlib.pyplot as plt
plt.figure(figsize=(14, 7))
plt.subplot(121);
plt.imshow(abs(np.conj(Vh.T) @ H @ Vh))
plt.title('diagonalized with eigh')
plt.subplot(122);
plt.imshow(abs(np.conj(V.T) @ H @ V))
plt.title('diagonalized with eig')

绘制这个

两个对角化都成功了，但特征值的顺序无关紧要。 如果您对特征值进行排序，您会看到它们匹配

plt.plot(np.diag(np.real(np.conj(Vh.T) @ H @ Vh)))
plt.plot(np.diag(np.imag(np.conj(Vh.T) @ H @ Vh)))
plt.plot(np.sort(np.diag(np.real(np.conj(V.T) @ H @ V))))
plt.title('eigenvalues')
plt.legend(['real eigh', 'imag eigh', 'sorted real eig'], loc='upper left')

由于许多特征值是重复的，与给定特征值关联的特征向量也不是唯一的，我们唯一可以保证的是给定特征值的特征向量必须跨越相同的子空间。

我认为对角化测试是最好的。

eigh 总是比 eig 好吗？

如果您在lapack routines 中搜索特征值，您将有很多选择。所以我不能在这里讨论每个可能的实现。常识说我们可以期望对称/厄米特例程表现更好，否则就没有理由再添加一个更有限的例程。但我从未仔细测试过eig 与eigh 的行为。

为了直观地比较对称矩阵的三对角化方程，以及将一般矩阵简化为海森堡形式的方程，请参见 here。

【讨论】：

如果我们知道矩阵是对称矩阵或厄米矩阵，是否可以说使用 eigh 总是更好？