SymPy生成的简并子空间特征向量无法构造有效投影算子的问题排查

我尝试用SymPy符号化计算列表An中各矩阵对应特征子空间的投影算子集合，生成An的代码如下：

import numpy as np
import sympy as sp
qutritketsn = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
qutritketss = sp.Matrix(qutritketsn)
an = [(1/sp.sqrt(1+sp.cos(sp.pi/5)))*(sp.cos(4*sp.pi*n/5)*qutritketss.col(0)+sp.sin(4*sp.pi*n/5)*qutritketss.col(1)+sp.sqrt(sp.cos(sp.pi/5))*qutritketss.col(2)) for n in range(5)]
An = [sp.eye(3)-2*(v*v.T.conjugate()) for v in an]

注：An中的矩阵特征值均应为-1、1、1（1为二重简并）。

尝试的实现代码 #

我用以下代码实现投影算子计算：

eigproj = []
for num in range(len(An)):
    eig_stuff = An[num].eigenvects(simplify=True)
    eigproj_temp = []
    for i in range(len(eig_stuff)):
        if eig_stuff[i][1]<2:
            evnorm = eig_stuff[i][2][0].normalized()
            proj = evnorm*evnorm.T.conjugate()
            eigproj_temp.append(proj)
        else: # 对简并子空间的投影算子求和
            shp = eig_stuff[i][2][0].shape
            init = sp.zeros(shp[0], shp[0])
            for j in range(len(eig_stuff[i][2])):
                evnorm1 = eig_stuff[i][2][j].normalized()
                init += evnorm1*evnorm1.T.conjugate()
            eigproj_temp.append(init)

遇到的问题 #

运行后发现，eig_stuff中返回的简并特征值对应的特征向量，经归一化后逐个计算投影算子再求和，无法得到有效投影算子：

• 生成的简并子空间投影算子不满足幂等性
• 与非简并特征值对应的投影算子不正交
• 两者之和不等于单位矩阵

而非简并特征值对应的投影算子是满足幂等性的，问题应该出在简并特征值的处理上。

NumPy对比实现 #

用等价的NumPy代码处理An的NumPy版本时，返回的特征向量可以构造有效投影算子，参考代码如下：

from collections import Counter
eigproj = []
for num in range(len(An)):
    eigs_temp, eigvs_temp = np.linalg.eigh(An[num])
    # 将接近整数的特征值取整
    for i in range(len(eigs_temp)):
        if np.abs(eigs_temp[i].real-np.round(eigs_temp[i].real)) < 10**(-13):
            temps = np.round(eigs_temp[i].real)
            eigs_temp.real[i] = temps    
    eig_mult = Counter(eigs_temp) # 统计特征值重数
    eig_mult = list(eig_mult.values())
    # 根据重数获取特征值索引
    dup_inds = [[k for k in range(int(int(i>0)*sum([eig_mult[j] for j in range(max(i, 0))])), sum([eig_mult[j] for j in range(i+1)]))] for i in range(len(eig_mult))]
    eigproj_temp = [np.outer(eigvs_temp[:, i], np.conj(eigvs_temp[:, i])) for i in range(eigvs_temp.shape[0])]
    eigproj_temp1 = []
    for i in dup_inds:
        # 对简并子空间的投影算子求和
        if len(i)>1:
            og = np.zeros(eigproj_temp[0].shape)
            for j in i:
                og += eigproj_temp[j]
            eigproj_temp1.append(og)
        else:
            eigproj_temp1.append(eigproj_temp[i[0]])

    eigproj_temp1 = np.array(eigproj_temp1)
    eigproj.append(eigproj_temp1)

示例结果对比 #

以如下矩阵为例：

Matrix([
[0.27639320225002103036, 0.52573111211913360603, 0.80449586419071039205],
[0.52573111211913360603, 0.61803398874989484820, -0.58450045893897621288],
[0.80449586419071039205, -0.58450045893897621288, 0.10557280900008412144]
])

SymPy下20位精度的投影算子结果：

# 特征值-1对应的投影算子
Matrix([
[0.36180339887498948482, -0.26286555605956680301, -0.40224793209535519602],
[-0.26286555605956680301, 0.19098300562505257590, 0.29225022946948810644],
[-0.40224793209535519602, 0.29225022946948810644, 0.44721359549995793928]
])

# 简并特征值1对应的投影算子
Matrix([
[0.89827790731256834867, 0.47552825814757678606, 0.49720578787857844011],
[0.47552825814757678606, 0.65450849718747371205, 0],
[0.49720578787857844011, 0, 0.44721359549995793928]
])

NumPy下的输出结果：

# 特征值-1对应的投影算子
array([
[ 0.3618034 , -0.26286556, -0.40224793],
[-0.26286556,  0.19098301,  0.29225023],
[-0.40224793,  0.29225023,  0.4472136 ]
])

# 简并特征值1对应的投影算子
array([
[ 0.6381966 ,  0.26286556,  0.40224793],
[ 0.26286556,  0.80901699, -0.29225023],
[ 0.40224793, -0.29225023,  0.5527864 ]
])

核心疑问 #

SymPy为何返回这类特征向量？问题出在我的代码实现还是SymPy本身？如果是我的问题，哪里出错了？

内容的提问来源于stack exchange，提问作者Riemann