横向场伊辛模型代码优化求助：大矩阵内存溢出问题

核心问题分析 #

横向场伊辛模型的全空间希尔伯特空间维度随自旋数N呈$2^N$指数增长，直接对角化的方法在N≥15时会因内存不足崩溃——即使启用稀疏对角化，也无法突破高维矩阵的内存瓶颈。但该模型是严格可解的，我们可以利用其映射到自由费米子的特性，绕过全空间矩阵对角化，直接计算单个自旋的约化密度矩阵。

优化方案：基于Jordan-Wigner变换的解析+高效数值计算 #

横向场伊辛模型可通过Jordan-Wigner变换转化为无相互作用的费米子模型，单个自旋的约化密度矩阵特征值可由费米子关联函数直接推导，无需处理高维态向量：

• 单个自旋的约化密度矩阵$\rho_1$是2x2矩阵，特征值为$\frac{1\pm\sqrt{1-C^{2}}{2}$，其中$C=\langle\sigma_1}x\rangle$是基态下自旋1的x方向关联函数
• 冯诺依曼熵$S=-\lambda_1\ln\lambda_1-\lambda_2\ln\lambda_2$，$\lambda_1,\lambda_2$为$\rho_1$的特征值

对于周期性边界条件的模型，基态下的$\langle\sigma_j^x\rangle$可通过费米子色散关系求解：
$$C = \frac{1}{N}\sum_{k=0}^{N-1}\frac{\cos k - \lambda}{\sqrt{(\cos k - \lambda)^2 + \sin^2 k}}$$

优化后的代码 #

import numpy as np
import matplotlib.pyplot as plt

def von_neumann_entropy_ising(N, lambd):
    # 生成周期性边界条件下的费米子动量点
    k = 2 * np.pi * np.arange(N) / N
    # 计算每个动量点的关联项
    numerator = np.cos(k) - lambd
    denominator = np.sqrt((np.cos(k) - lambd)**2 + np.sin(k)**2)
    # 计算基态下<σ₁ˣ>的平均值
    C = np.mean(numerator / denominator)
    # 求解约化密度矩阵的特征值
    lambda1 = (1 + np.sqrt(1 - C**2)) / 2
    lambda2 = (1 - np.sqrt(1 - C**2)) / 2
    # 计算冯诺依曼熵（避免0值的对数运算）
    S_vn = 0.0
    if lambda1 > 1e-12:
        S_vn -= lambda1 * np.log(lambda1)
    if lambda2 > 1e-12:
        S_vn -= lambda2 * np.log(lambda2)
    return S_vn

# 计算并绘制结果
Ns = range(2, 31)  # 可轻松计算至N=30甚至更大规模
lambd = 1.0
entropies = []

for N in Ns:
    S_vn = von_neumann_entropy_ising(N, lambd)
    entropies.append(S_vn)

plt.figure(figsize=(10, 6))
plt.plot(Ns, entropies, marker='o')
plt.xlabel('自旋数 (N)')
plt.ylabel('冯诺依曼熵')
plt.title('横向场伊辛模型冯诺依曼熵与自旋数的关系')
plt.grid(True)
plt.show()

效果说明 #

• 该方法时间复杂度为$O(N)$，空间复杂度为$O(N)$，完全摆脱$2^N$维度限制，在Colab环境下可轻松计算N=1000甚至更大规模的系统
• 结果为精确解，完全满足原问题对精确计算的要求，无近似误差

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