LAPACK 的 ZHEEVR 例程的特征向量不准确

Question

我试图用 C 级 LAPACK 替换我的 cython 代码中的 Scipy 的 eigh 例程。本质上，我想计算与尺寸约为 50x50 的埃尔米特矩阵的最大和最小特征值相对应的特征向量。我制作了一个小型 3x3 示例，以确保获得与 eigh 相同的结果，它内部似乎使用了 ZHEEVR 例程。 Scipy 提供了 LAPACK 例程的 C 级版本，可以导入到我正在使用的 cython 中。我认为我的语法等正确，因为 eigh 和 ZHEEVR 的特征向量之一和所有特征值完全相同。然而，ZHEEVR 的其他特征向量不准确。事实上，我没有获得正确的分解。

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

info = 0 告诉我算法执行正确，特征值正确，但是如您所见， Vectors_zheevr_np.T.conj() @ A @ paths_zheevr_np 不是对角矩阵，因为它应该是正确的分解 - 右上和左下元素非零。

有人以前遇到过这个问题，知道如何解决它，或者有任何想法吗？非常感谢！

Answer 1

正如 Nick ODell 指出的，LAPACK 需要 FORTRAN 连续的输入数组。所以，更换

cdef double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
        A.astype(np.complex128))

与

cdef double complex * A_pointer = <double complex *>cnp.PyArray_DATA(
            A.T.conj().astype(np.complex128))

解决问题。此外，请注意，在特征值不同的情况下，埃尔米特矩阵的复特征向量仅在乘法（复数）常数范围内是唯一的。