复厄密矩阵的特征分析:EIG 和 EIGH 的不同相位角
Eigenanalysis of complex hermitian matrix: different phase angles for EIG and EIGH
我知道特征向量只定义为一个乘法常数。据我所知,所有 numpy
算法(例如 linalg.eig
、linalg.eigh
、linalg.svd
)都会为 实矩阵 产生相同的特征向量,所以显然他们使用相同的规范化。然而,在复杂矩阵的情况下,算法会产生不同的结果。
也就是说,特征向量相同直到(复数)常数 z
。在对 eig
和 eigh
进行了一些试验后,我意识到 eigh
总是将相位角(定义为 arctan(complex part/real part))设置为 0 每个的第一个分量特征向量而 eig
似乎以一些(任意?)非零相位角开始。
问:有没有办法按照 eig
的方式对 eigh
的特征向量进行归一化(即不强制相角 = 0)?
例子
我有一个复杂的厄密矩阵 G
,我想使用以下两种算法计算其特征向量:
numpy.linalg.eig
对于实数/复数方阵
numpy.linalg.eigh
对于实对称/复厄密矩阵(1 的特例)
检查 G 是否为 hermitian
# check if a matrix is hermitian
def isHermitian(a, rtol=1e-05, atol=1e-08):
return np.allclose(a, a.conjugate().T, rtol=rtol, atol=atol)
print('G is hermitian:', isHermitian(G))
输出:
G is hermitian: True
执行特征分析
# eigenvectors from EIG()
l1,u1 = np.linalg.eig(G)
idx = np.argsort(l1)[::-1]
l1,u1 = l1[idx].real,u1[:,idx]
# eigenvectors from EIGH()
l2,u2 = np.linalg.eigh(G)
idx = np.argsort(l2)[::-1]
l2,u2 = l2[idx],u2[:,idx]
检查特征值
print('Eigenvalues')
print('eig\t:',l1[:3])
print('eigh\t:',l2[:3])
输出:
Eigenvalues
eig : [2.55621629e+03 3.48520440e+00 3.16452447e-02]
eigh : [2.55621629e+03 3.48520440e+00 3.16452447e-02]
两种方法产生相同的特征向量。
检查特征向量
现在查看特征向量(例如 3. 特征向量),它们相差一个常数因子 z
。
multFactors = u1[:,2]/u2[:,2]
if np.count_nonzero(multFactors[0] == multFactors):
print("All multiplication factors are same:", multFactors[0])
else:
print("Multiplication factors are different.")
输出:
All multiplication factors are same: (-0.8916113627685007+0.45280147727156245j)
检查相位角
现在检查 3. 特征向量的第一个分量的相位角:
print('Phase angel (in PI) for first point:')
print('Eig\t:',np.arctan2(u1[0,2].imag,u1[0,2].real)/np.pi)
print('Eigh\t:',np.arctan2(u2[0,2].imag,u2[0,2].real)/np.pi)
输出:
Phase angel (in PI) for first point:
Eig : 0.8504246311627189
Eigh : 0.0
复制图的代码
num = 2
fig = plt.figure()
gs = gridspec.GridSpec(2, 3)
ax0 = plt.subplot(gs[0,0])
ax1 = plt.subplot(gs[1,0])
ax2 = plt.subplot(gs[0,1:])
ax3 = plt.subplot(gs[1,1:])
ax2r= ax2.twinx()
ax3r= ax3.twinx()
ax0.imshow(G.real,vmin=-30,vmax=30,cmap='RdGy')
ax1.imshow(G.imag,vmin=-30,vmax=30,cmap='RdGy')
ax2.plot(u1[:,num].real,label='eig')
ax2.plot((u2[:,num]).real,label='eigh')
ax3.plot(u1[:,num].imag,label='eig')
ax3.plot((u2[:,num]).imag,label='eigh')
for a in [ax0,ax1,ax2,ax3]:
a.set_xticks([])
a.set_yticks([])
ax0.set_title('Re(G)')
ax1.set_title('Im(G)')
ax2.set_title('Re('+str(num+1)+'. Eigenvector)')
ax3.set_title('Im('+str(num+1)+'. Eigenvector)')
ax2.legend(loc=0)
ax3.legend(loc=0)
fig.subplots_adjust(wspace=0, hspace=.2,top=.9)
fig.suptitle('Eigenanalysis of Hermitian Matrix G',size=16)
plt.show()
根据我的经验(这里有很多问题可以支持这一点),当 eigh
是一个选项时,您 永远不会 想使用 eig
- eig
非常慢而且非常不稳定。这样做的相关性是我相信你的问题是落后的——你想将 eig
的特征向量归一化为 eigh
的特征向量,而你知道该怎么做。
正如您所说,特征值问题仅将特征向量固定为标量x
。将特征向量 v
转换为 v = v*x
不会改变其作为特征向量的状态。
有一种“明显”的方式来归一化向量(根据欧氏内积 np.vdot(v1, v1)
),但这只能固定标量的振幅,它可能是复数。
固定角度或“相位”有点武断,没有进一步的上下文。我尝试了 eigh()
,实际上它只是使向量的第一个条目变为实数(带有明显的随机符号!?)。
eig()
而是选择将具有最大实部的向量条目变为实数。例如,这是我得到的随机 Hermitian 矩阵:
n = 10
H = 0.5*(X + X.conj().T)
np.max(la.eig(H)[1], axis=0)
# returns
array([0.57590624+0.j, 0.42672485+0.j, 0.51974879+0.j, 0.54500475+0.j,
0.4644593 +0.j, 0.53492448+0.j, 0.44080532+0.j, 0.50544424+0.j,
0.48589402+0.j, 0.43431733+0.j])
这可以说是更明智的,因为如果第一个条目恰好非常小,那么像 eigh()
那样只选择第一个条目并不是很稳健。选择最大值可以避免这种情况。我不确定 eig()
是否也修复了符号(随机矩阵不是一个很好的测试用例,因为特征向量中的所有条目都具有负实部是非常不寻常的,这是唯一的情况其中会出现一个未固定的标志)。
无论如何,我不会依赖特征求解器使用任何特定的固定相的方法。它没有记录在案,因此原则上将来可能会发生变化。相反,您可以自己修复这些阶段,也许 eig()
现在就是这样做的。
我知道特征向量只定义为一个乘法常数。据我所知,所有 numpy
算法(例如 linalg.eig
、linalg.eigh
、linalg.svd
)都会为 实矩阵 产生相同的特征向量,所以显然他们使用相同的规范化。然而,在复杂矩阵的情况下,算法会产生不同的结果。
也就是说,特征向量相同直到(复数)常数 z
。在对 eig
和 eigh
进行了一些试验后,我意识到 eigh
总是将相位角(定义为 arctan(complex part/real part))设置为 0 每个的第一个分量特征向量而 eig
似乎以一些(任意?)非零相位角开始。
问:有没有办法按照 eig
的方式对 eigh
的特征向量进行归一化(即不强制相角 = 0)?
例子
我有一个复杂的厄密矩阵 G
,我想使用以下两种算法计算其特征向量:
numpy.linalg.eig
对于实数/复数方阵numpy.linalg.eigh
对于实对称/复厄密矩阵(1 的特例)
检查 G 是否为 hermitian
# check if a matrix is hermitian
def isHermitian(a, rtol=1e-05, atol=1e-08):
return np.allclose(a, a.conjugate().T, rtol=rtol, atol=atol)
print('G is hermitian:', isHermitian(G))
输出:
G is hermitian: True
执行特征分析
# eigenvectors from EIG()
l1,u1 = np.linalg.eig(G)
idx = np.argsort(l1)[::-1]
l1,u1 = l1[idx].real,u1[:,idx]
# eigenvectors from EIGH()
l2,u2 = np.linalg.eigh(G)
idx = np.argsort(l2)[::-1]
l2,u2 = l2[idx],u2[:,idx]
检查特征值
print('Eigenvalues')
print('eig\t:',l1[:3])
print('eigh\t:',l2[:3])
输出:
Eigenvalues
eig : [2.55621629e+03 3.48520440e+00 3.16452447e-02]
eigh : [2.55621629e+03 3.48520440e+00 3.16452447e-02]
两种方法产生相同的特征向量。
检查特征向量
现在查看特征向量(例如 3. 特征向量),它们相差一个常数因子 z
。
multFactors = u1[:,2]/u2[:,2]
if np.count_nonzero(multFactors[0] == multFactors):
print("All multiplication factors are same:", multFactors[0])
else:
print("Multiplication factors are different.")
输出:
All multiplication factors are same: (-0.8916113627685007+0.45280147727156245j)
检查相位角
现在检查 3. 特征向量的第一个分量的相位角:
print('Phase angel (in PI) for first point:')
print('Eig\t:',np.arctan2(u1[0,2].imag,u1[0,2].real)/np.pi)
print('Eigh\t:',np.arctan2(u2[0,2].imag,u2[0,2].real)/np.pi)
输出:
Phase angel (in PI) for first point:
Eig : 0.8504246311627189
Eigh : 0.0
复制图的代码
num = 2
fig = plt.figure()
gs = gridspec.GridSpec(2, 3)
ax0 = plt.subplot(gs[0,0])
ax1 = plt.subplot(gs[1,0])
ax2 = plt.subplot(gs[0,1:])
ax3 = plt.subplot(gs[1,1:])
ax2r= ax2.twinx()
ax3r= ax3.twinx()
ax0.imshow(G.real,vmin=-30,vmax=30,cmap='RdGy')
ax1.imshow(G.imag,vmin=-30,vmax=30,cmap='RdGy')
ax2.plot(u1[:,num].real,label='eig')
ax2.plot((u2[:,num]).real,label='eigh')
ax3.plot(u1[:,num].imag,label='eig')
ax3.plot((u2[:,num]).imag,label='eigh')
for a in [ax0,ax1,ax2,ax3]:
a.set_xticks([])
a.set_yticks([])
ax0.set_title('Re(G)')
ax1.set_title('Im(G)')
ax2.set_title('Re('+str(num+1)+'. Eigenvector)')
ax3.set_title('Im('+str(num+1)+'. Eigenvector)')
ax2.legend(loc=0)
ax3.legend(loc=0)
fig.subplots_adjust(wspace=0, hspace=.2,top=.9)
fig.suptitle('Eigenanalysis of Hermitian Matrix G',size=16)
plt.show()
根据我的经验(这里有很多问题可以支持这一点),当 eigh
是一个选项时,您 永远不会 想使用 eig
- eig
非常慢而且非常不稳定。这样做的相关性是我相信你的问题是落后的——你想将 eig
的特征向量归一化为 eigh
的特征向量,而你知道该怎么做。
正如您所说,特征值问题仅将特征向量固定为标量x
。将特征向量 v
转换为 v = v*x
不会改变其作为特征向量的状态。
有一种“明显”的方式来归一化向量(根据欧氏内积 np.vdot(v1, v1)
),但这只能固定标量的振幅,它可能是复数。
固定角度或“相位”有点武断,没有进一步的上下文。我尝试了 eigh()
,实际上它只是使向量的第一个条目变为实数(带有明显的随机符号!?)。
eig()
而是选择将具有最大实部的向量条目变为实数。例如,这是我得到的随机 Hermitian 矩阵:
n = 10
H = 0.5*(X + X.conj().T)
np.max(la.eig(H)[1], axis=0)
# returns
array([0.57590624+0.j, 0.42672485+0.j, 0.51974879+0.j, 0.54500475+0.j,
0.4644593 +0.j, 0.53492448+0.j, 0.44080532+0.j, 0.50544424+0.j,
0.48589402+0.j, 0.43431733+0.j])
这可以说是更明智的,因为如果第一个条目恰好非常小,那么像 eigh()
那样只选择第一个条目并不是很稳健。选择最大值可以避免这种情况。我不确定 eig()
是否也修复了符号(随机矩阵不是一个很好的测试用例,因为特征向量中的所有条目都具有负实部是非常不寻常的,这是唯一的情况其中会出现一个未固定的标志)。
无论如何,我不会依赖特征求解器使用任何特定的固定相的方法。它没有记录在案,因此原则上将来可能会发生变化。相反,您可以自己修复这些阶段,也许 eig()
现在就是这样做的。