LAPACK ZHEEVR routine计算特征向量不准确问题求助

我尝试在Cython代码中用C级别的LAPACK ZHEEVR routine替代Scipy的eigh函数，计算维度约50x50的厄米矩阵对应最大、最小特征值的特征向量。制作了3x3的测试示例，确保结果与eigh一致（eigh内部似乎也使用ZHEEVR）。使用Scipy提供的可在Cython中导入的C级LAPACK例程，语法等应该无误，因为其中一个特征向量和所有特征值与eigh完全一致，但ZHEEVR输出的其他特征向量不准确，甚至无法得到正确的分解结果。info=0表明算法执行成功，特征值也正确，但vectors_zheevr_np.T.conj() @ A @ vectors_zheevr_np并非对角矩阵（正确分解应满足对角化），其右上角和左下角元素非零。

测试代码 #

import numpy as np
cimport numpy as cnp

from scipy.linalg import eigh
from scipy.linalg.cython_lapack cimport zheevr
from libc.stdlib cimport malloc, free

dim = 3

A = (np.arange(dim**2) + 2j).reshape(dim, dim) + 2*np.eye(dim)
A = A @ np.transpose(np.conjugate(A))

print(A)
print('')

values_eigh, vectors_eigh = eigh(A)

print(values_eigh)
print('')
print(vectors_eigh)
print('')

cdef:
    char jobz = 'V'
    char rrange = 'A'
    char uplo = 'U'
    int n = dim
    double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
        A.astype(np.complex128))
    int lda = dim
    double vl = 0
    double vu = 0
    int il = 0
    int iu = 0
    double abstol = 1e-12
    int m = n
    double * values_zheever = <double *> malloc(
        n * sizeof(double))
    double complex * vectors_zheever = <double complex *>malloc(
        n * n * sizeof(double complex))
    int ldz = n
    int * isuppz = <int *> malloc(2 * n * sizeof(int))
    double complex * work = <double complex *>malloc(
        1 * sizeof(double complex))
    int lwork = -1
    double * rwork = <double *>malloc(1 * sizeof(double))
    int lrwork = -1
    int * iwork = <int *>malloc(1 * sizeof(int))
    int liwork = -1
    int info

zheevr(&jobz, &rrange, &uplo, &n, A_pointer, &lda, &vl, &vu, &il, &iu, &abstol, &m, values_zheever, vectors_zheever, &ldz, isuppz, work, &lwork, rwork, &lrwork, iwork, &liwork, &info)

lwork = <int> work[0]
lrwork = <int> rwork[0]
liwork = <int> iwork[0]
free(work)
free(rwork)
free(iwork)
work = <double complex *>malloc(lwork * sizeof(double complex))
rwork = <double *>malloc(lrwork * sizeof(double))
iwork = <int *>malloc(liwork * sizeof(int))

zheevr(&jobz, &rrange, &uplo, &n, A_pointer, &lda, &vl, &vu, &il, &iu, &abstol, &m, values_zheever, vectors_zheever, &ldz, isuppz, work, &lwork, rwork, &lrwork, iwork, &liwork, &info)

for i in range(dim):
    print(values_zheever[i])

vectors_zheevr_np = np.zeros((n* n), dtype=np.complex128)
for i in range(n*n):
        vectors_zheevr_np[i] = vectors_zheever[i]
vectors_zheevr_np = vectors_zheevr_np.reshape(n, n).T

print('')
print(vectors_zheevr_np)
print('')
print(vectors_zheevr_np.T.conj() @ A @ vectors_zheevr_np)
print('')
print(info)

运行结果 #

[[ 21. +0.j  34.+18.j  51.+36.j]
 [ 34.-18.j  82. +0.j 122.+18.j]
 [ 51.-36.j 122.-18.j 197. +0.j]]

[  0.82663284   4.         295.17336716]

[[ 0.8755536 -0.00000000e+00j -0.40824829-0.00000000e+00j
  -0.25833936-0.00000000e+00j]
 [ 0.24467905+6.98442317e-02j  0.81649658+1.22124533e-15j
  -0.46103588+2.36713326e-01j]
 [-0.38619551+1.39688463e-01j -0.40824829-6.10622664e-16j
  -0.6637324 +4.73426651e-01j]]

0.8266328441181744
4.000000000000024
295.1733671558813

[[-0.8233493 +2.97808740e-01j -0.40824829-2.22044605e-16j
   0.21031944-1.50016524e-01j]
 [-0.20633357+1.48904370e-01j  0.81649658-1.56819002e-15j
   0.51279727-7.50082620e-02j]
 [ 0.41068216-0.00000000e+00j -0.40824829-0.00000000e+00j
   0.8152751 +0.00000000e+00j]]

[[ 2.72444561e+01+3.55271368e-15j -1.35003120e-13-1.42108547e-14j
   1.95681896e+01-8.18241075e+01j]
 [-1.30784272e-13+2.16125482e-14j  4.00000000e+00-4.07921987e-16j
  -2.66453526e-14+4.85653027e-13j]
 [ 1.95681896e+01+8.18241075e+01j -2.84217094e-14-4.84945417e-13j
   2.68755544e+02+0.00000000e+00j]]

0

问题原因与解决方法 #

问题出在特征向量的存储顺序处理错误：

• LAPACK的ZHEEVR采用列优先存储，输出的特征向量矩阵vectors_zheever中，第k个特征向量是矩阵的第k列，即内存中按[v1_1, v2_1, v3_1, v1_2, v2_2, v3_2, v1_3, v2_3, v3_3]的顺序排列（对应3x3矩阵）。
• 你的代码中将vectors_zheeverreshape为(n,n)后进行了转置（.T），导致特征向量的行列被颠倒，后续的对角化验证自然无法得到正确的对角矩阵。

修正代码 #

将向量转换部分的代码修改为：

vectors_zheevr_np = np.zeros((n* n), dtype=np.complex128)
for i in range(n*n):
        vectors_zheevr_np[i] = vectors_zheever[i]
vectors_zheevr_np = vectors_zheevr_np.reshape(n, n)  # 移除.T转置操作

或者更高效的方式（直接内存拷贝，避免循环）：

vectors_zheevr_np = np.zeros((n, n), dtype=np.complex128)
cnp.memmove(cnp.PyArray_DATA(vectors_zheevr_np), vectors_zheever, n*n*sizeof(double complex))

验证修正效果 #

修正后，执行vectors_zheevr_np.conj().T @ A @ vectors_zheevr_np会得到接近对角的矩阵，特征向量的正确性也会和Scipy eigh的结果一致（允许符号差异，因为特征向量的符号不影响分解正确性）。

额外注意：如果只需要最大/最小的几个特征值和特征向量，可以调整ZHEEVR的rrange、vl/vu或il/iu参数，避免计算全部特征向量，提升效率（比如rrange='V'，设置vl/vu为目标范围，或者rrange='I'设置il/iu为目标索引）。

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