added codes

doc/Programs/EigenvaluesDeepLearning/EigRNN.py

@@ -0,0 +1,81 @@
+import numpy as np
+import time
+import sys
+import matplotlib.pyplot as plt
+from hamiltonian import *
+
+class EigRNN:
+    """
+    Finds the lowest eigenvalue of a symmetric matrix
+    by minimizing the Rayleigh quotient with gradient
+    descent
+    """
+    def __init__(self,A,eps=1e-4,maxiter=1000):
+        """
+        A - input matrix
+        eps - stops iterating when eigenvalue is not changing more than eps
+        maxiter - maximum number of iterations
+        """
+        self.A = A
+        self.maxiter = maxiter
+        self.eps = eps
+        self.x = np.random.rand(self.A.shape[0])
+
+    def dR(self):
+        """
+        Returns the gradient of the Rayleigh quotient
+        """
+        x = self.x
+        self.Ax = self.A@x
+        self.xTx = x.T@x
+        self.eig = x.T@self.Ax/self.xTx
+        self.grad = 2*((self.xTx)*self.Ax - (x.T@self.Ax)*x)/((self.xTx)*(self.xTx))
+        return(self.grad)
+
+    def optStep(self):
+        """
+        Returns the optimal step length
+        """
+        dR = self.grad
+        x = self.x
+        xTA = self.Ax.T
+        a = xTA@x
+        b = xTA@dR
+        c = dR.T@self.A@dR
+        e = self.xTx
+        f = x.T@dR
+        g = dR.T@dR
+        return( ( (a*g - e*c ) + np.sqrt( (e*c - a*g)**2 - 4*(a*f - e*b)*(b*g - c*f) ) )/(2*(b*g - c*f)))
+
+    def solve(self):
+        """
+        Returns the lowest eigenvalue of A.
+        To be called after initialization.
+        """
+        convergence = False
+        val = 0
+        for i in range(self.maxiter):
+            grad = self.dR()
+            self.x = self.x - self.optStep()*grad
+            if np.sum(np.abs(grad)) < self.eps*(np.log(self.x.shape[0])):
+                convergence = True
+                break
+
+        if not convergence:
+            print('WARNING: Did not converge. Try increasing maxiter or a smaller eps.')    
+
+        return(self.eig,self.x)
+
+
+
+if __name__ == '__main__':
+    n_pairs = int(sys.argv[1])
+    n_basis = int(sys.argv[2])
+    print('System with {} pairs and {} basis states'.format(n_pairs,n_basis))
+    delta = float(sys.argv[3])
+    g = float(sys.argv[4])
+    epsilon = float(sys.argv[5])
+    H,Eref = hamiltonian(n_pairs,n_basis,delta,g)
+    RNN = EigRNN(H,eps=epsilon)
+    eigval,eigvec = RNN.solve()
+    print('Energy: {}'.format(eigval))

doc/Programs/EigenvaluesDeepLearning/hamiltonian.py

@@ -0,0 +1,44 @@
+import scipy as scipy
+import scipy.special as special
+import numpy as np
+import itertools as it
+from pandas import *
+import matplotlib.pylab as plt
+from scipy.integrate import ode
+import time
+
+def unique_rows(a):
+    a = np.ascontiguousarray(a)
+    unique_a = np.unique(a.view([('', a.dtype)]*a.shape[1]))
+    return unique_a.view(a.dtype).reshape((unique_a.shape[0], a.shape[1]))
+
+def hamiltonian(n_pairs,n_basis,delta,g):
+	"""
+	n_pairs - Number of electron pairs
+	n_basis - Number of spacial basis states
+	Returns - The hamiltonian and the reference energy
+	"""
+	n_SD = int(special.binom(n_basis,n_pairs))
+	H_mat = np.zeros((n_SD,n_SD))
+	S = stateMatrix(n_pairs,n_basis)
+	for row in range(n_SD):
+		bra = S[row,:]
+		for col in range(n_SD):
+			ket = S[col,:]
+			if np.sum(np.equal(bra,ket)) == bra.shape:
+				H_mat[row,col] += 2*delta*np.sum(bra - 1) - 0.5*g*n_pairs
+			if n_pairs - np.intersect1d(bra,ket).shape[0] == 1:
+				H_mat[row,col] += -0.5*g
+	return(H_mat,H_mat[0,0])
+
+
+def stateMatrix(n_pairs,n_basis):
+	L = []
+	states = range(1,n_basis+1)
+	for perm in it.permutations(states,n_pairs):
+		L.append(perm)
+	L = np.array(L)
+	L.sort(axis=1)
+	L = unique_rows(L)
+	return(L)
+