python中的递归函数定义

Question

因此，我试图实现一个版本的哈特里福克理论的带系统。基本上，这是一个矩阵收敛问题。我有一个矩阵H0，从它的特征值我可以构造另一个矩阵F。然后程序是定义H1=H0+F，并检查H1的特征值是否接近H0的特征值。如果不是，我从H1的特征值构造一个新的F，并定义H2=H0+F，然后再次检查并迭代。你知道吗

这个问题有点泛化，我的确切代码似乎并不真正相关。所以我只展示这个：

# define the matrix F
def generate_fock(H):
    def fock(k): #k is a 2D array
        matt = some prefactor*outer(eigenvectors of H(k) with itself) #error1
        return matt
    return fock

k0 = linspace(0,5*kt/2,51)
# H0 is considered defined
H = lambda k: H0(k)
evalold = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[the ones I care]
while True:
    fe = generate_fock(H)
    H = lambda k: H0(k)+fe(k) #error2
    evalnew = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[the ones I care]
    if allclose(evalnew, evalold): break
    else: evalold = evalnew

我正在使用内部函数，希望python不会发现我的定义是递归的（我不确定我是否正确地使用了这个词）。但是python知道：有什么建议吗？你知道吗

编辑1： 错误消息将突出显示我标记为error1和error2的行，并显示以下内容：

RecursionError: maximum recursion depth exceeded while calling a Python object

我认为这来自于我定义函数的方式：在循环n中，F（k）依赖于上一个循环的H（k），下一步的H（k）又依赖于F（k）。我的问题是我该怎么解决这个问题？你知道吗

编辑2&3： 让我按照建议为代码添加更多细节。这是我能想出的最简单的方法，能准确地再现我的问题。你知道吗

from numpy import *
from scipy import linalg

# Let's say H0 is any 2m by 2m Hermitian matrix. m = 4 in this case.
# Here are some simplified parameters 
def h(i,k):
    return -6*linalg.norm(k)*array([[0,exp(1j*(angle(k@array([1,1j]))+(-1)**i*0.1/2))],
                                    [exp(-1j*(angle(k@array([1,1j]))+(-1)**i*0.1/2)),0]])
T = array([ones([2,2]),
          [[exp(-1j*2*pi/3),1],[exp(1j*2*pi/3),exp(-1j*2*pi/3)]],
          [[exp(1j*2*pi/3),1],[exp(-1j*2*pi/3),exp(1j*2*pi/3)]]])
g = array([[ 0.27023695,  0.46806412], [-0.27023695,  0.46806412]])
kt = linalg.norm(g[0])
def H0(k):
    "one example"
    matt = linalg.block_diag(h(1,k),h(2,k+g[0]),h(2,k+g[1]),h(2,k+g[0]+g[1]))/2
    for j in range(3): matt[0:2,2*j+2:2*j+4] = T[j]
    return array(matrix(matt).getH())+matt
dim = 4
def bz(x):
    "BZ centered at 0 with (2x)^2 points in it"
    tempList = []
    for i in range(-x,x):
        for j in range(-x,x):
            tempList.append(i*g[0]/2/x+j*g[1]/2/x)
    return tempList
def V(i,G):
    "2D Coulomb interaction"
    if linalg.norm(G)==0: return 0
    if i>=dim: t=1
    else: t=0
    return 2*pi/linalg.norm(G)*exp(0.3*linalg.norm(G)*(-1+(-1)**t)/2)

# define the matrix F for some H
def generate_fock(H):
    def fock(k): #k is a 2D array
        matf = zeros([2*dim,2*dim],dtype=complex128)
        for pt in bz(1): #bz is a list of 2D arrays
            matt = zeros([2*dim,2*dim],dtype=complex128)
            eig_vals1, eig_vecs1 = linalg.eigh(H(pt)) #error1
            idx = eig_vals1.argsort()[::]
            vecs1 = eig_vecs1[:,idx][:dim]
            for vec in vecs1:
                matt = matt + outer(conjugate(vec),vec)
            matt = matt.transpose()/len(bz(1))
            for i in range(2*dim):
                for j in range(2*dim):
                    matf[i,j] = V(j-i,pt-k)*matt[i,j] #V is some prefactor
        return matf
    return fock

k0 = linspace(0,5*kt/2,51)
H = lambda k: H0(k)
evalold = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[dim-1:dim+1]
while True:
    fe = generate_fock(H)
    H = lambda k: H0(k)+fe(k) #error2
    evalnew = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[dim-1:dim+1]
    if allclose(evalnew, evalold): break
    else: evalold = evalnew

Answer 1

问题在于这些线路：

while True:
    fe = generate_fock(H)
    H = lambda k: H0(k)+fe(k) #error2

在每次迭代中，都会生成一个引用旧函数的新函数，而不是该旧函数的最终输出，因此必须将它们全部保留在堆栈中。这也将是非常缓慢的，因为你必须反乘所有矩阵每次迭代。你知道吗

您要做的是保留旧值的输出，可能是从先前迭代的结果中创建一个列表，然后应用该列表中的函数。你知道吗

你甚至可以用一个缓存来实现这一点，尽管它可能会变得很大。保留函数输入的字典并使用它。像这样：

# define the matrix F
def generate_fock(H):
    d = {}
    def fock(k): #k is a 2D array
        if k in d:
            return d[k]
        matt = some prefactor*outer(eigenvectors of H(k) with itself) #error1
        d[k] = matt
        return matt
    return fock

那么它应该只需要引用函数的最后一个版本。你知道吗

编辑：试试这个。除了缓存结果外，还要将索引保存到函数数组中，而不是引用。这应该可以防止递归深度溢出。你知道吗

hList = []

# define the matrix F
def generate_fock(n):
    d = {}
    def fock(k): #k is a 2D array
        if k in d:
            return d[k]
        matt = some prefactor*outer(eigenvectors of hList[n](k) with itself) #error1
        d[k] = matt
        return matt
    return fock

k0 = linspace(0,5*kt/2,51)
# H0 is considered defined
HList.append(lambda k: H0(k))
H = HList[0]
evalold = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[the ones I care]
n = 0
while True:
    fe = generate_fock(n)
    n += 1
    hList.append(lambda k: H0(k)+fe(k)) #error2
    H = hList[-1]
    evalnew = array([sort(linalg.eigvalsh(H(array([0,2*k-kt/2])))) for k in k0[:]])[the ones I care]
    if allclose(evalnew, evalold): break
    else: evalold = evalnew