目前Eigen暂时不支持直接对超大稀疏矩阵求最小特征值，需要手写一些步骤。
1. 将超大型稀疏矩阵投影到Krylov子空间，Krylov子空间最大为1W*1W
2. Krylov子空间是通过随机向量迭代的办法生成，思想来源于幂法
3. 该方法将超大型稀疏矩阵转化为三对角矩阵，也可以用稀疏存储。
4. 只需要求解投影后的矩阵的最小特征值，这里C++ Eigen库支持的很完善。
Python demo
def lanczos(hm,length): dimension = 2**length diagElements = [] subDiagElements = [] f0 = np.random.rand(dimension) f0 = f0 / norm(f0) f1 = hm @ f0 lastEnergy = 0 groundEnergy = 0 for truncationIndex in range(0,dimension): hdiag = f0 @ f1 f1 = f1 - hdiag * f0 hsub = norm(f1) f1 = f1 / hsub [f0 , f1] = [f1 , hm @ f1 -hsub*f0] diagElements.append(hdiag) # 常规方法对角化 # [va,ve] = scipy.linalg.eigh(self.fillElements(diagElements,subDiagElements,len(diagElements))) # groundEnergy = va[0] # if abs(groundEnergy - lastEnergy) < 1e-10: # break # 三对角矩阵方法 [va,ve] = eigh_tridiagonal(np.array(diagElements),np.array(subDiagElements)) groundEnergy = va[0] if(abs(groundEnergy - lastEnergy) < 1e-10): break lastEnergy = groundEnergy subDiagElements.append(hsub) return groundEnergyC++ 实现
EigVal GSLanczos(SparseMat &H, const u_int32_t size, const double eps = 1e-4, const int maxStep = 1000, const int K = 30) { std::vector<double> Kdiag; Kdiag.reserve(maxStep); std::vector<double> Ksub; Ksub.reserve(maxStep); double minValue = 0; double lastValue = 0; Eigen::VectorXcd KEigVec; Eigen::VectorXcd f0 = Eigen::VectorXcd::Random(size); f0 /= f0.norm(); Eigen::VectorXcd psi0 = f0; Eigen::VectorXcd f1 = H * f0; int step = 0; for (step = 0; step < maxStep; step++ ) { double a = f0.dot(f1).real(); f1 = f1 - a * f0; double b = f1.norm(); f1 = f1 / b; Eigen::VectorXcd tmp = f0; f0 = f1; f1 = H * f1 - b * tmp; Kdiag.push_back(a); if (step > K) { Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver; eigensolver.computeFromTridiagonal(Eigen::Map<Eigen::VectorXd>(Kdiag.data(), step), Eigen::Map<Eigen::VectorXd>(Ksub.data(), step - 1)); minValue = eigensolver.eigenvalues()[0]; if (abs(minValue - lastValue) < eps) { KEigVec = eigensolver.eigenvectors().col(0); break; } lastValue = minValue; } Ksub.push_back(b); } f0 = psi0; Eigen::VectorXcd groundState = KEigVec[0] * f0; f1 = H * f0; int m = KEigVec.rows(); for (int i = 1; i < m; i++) { double a = f0.dot(f1).real(); f1 = f1 - a * f0; double b = f1.norm(); f1 = f1 / b; groundState += KEigVec[i] * f1; Eigen::VectorXcd tmp = f0; f0 = f1; f1 = H * f1 - b * tmp; } return {minValue, f0};}