diff --git a/graph_kernel.py b/graph_kernel.py
new file mode 100644
index 0000000000000000000000000000000000000000..7478ce6ad272bb4c36abb56bef6cfe9a1457936b
--- /dev/null
+++ b/graph_kernel.py
@@ -0,0 +1,69 @@
+import numpy as np
+from numpy import linalg as lina
+import math
+from scipy import linalg
+
+from adjacency import Beta_adjacency_matrix, Coulomb_matrix
+
+
+
+#-------------------------------------------------------------#
+#--------------- Exponential graph kernel --------------------#
+#-------------------------------------------------------------#
+
+#--- compute the exponential kernel from decompositions
+def Graph_kernel (D1,L1,invL1,D2,L2,invL2,gamma):
+
+
+    #Dimensions
+    (n1,a)=D1.shape
+    (n2,a)=D2.shape
+
+    #--- Direct product graphs
+    #--- (keep in mind that we procede to a normalization,
+    #--- that's why we have three matrices)
+    Dx=np.kron(D1,D2)
+    Dx1=np.kron(D1,D1)
+    Dx2=np.kron(D2,D2)
+    
+    Left=np.kron(L1,L2)
+    Right=np.kron(invL1,invL2)
+
+    Left1=np.kron(L1,L1)
+    Right1=np.kron(invL1,invL1)
+
+    Left2=np.kron(L2,L2)
+    Right2=np.kron(invL2,invL2)
+    
+    #--- Diagonal matrices
+    U=np.zeros(n1*n2)
+    U1=np.zeros(n1*n1)
+    U2=np.zeros(n2*n2)
+
+
+    #choose any "reasonable" function ...
+    for i in range(0,n1*n2):
+#        U[i]=1/(1-gamma*Dx[i,i])
+        U[i]=math.exp(gamma*Dx[i,i])
+
+    for i in range(0,n1*n1):
+        #U1[i]=1/(1-gamma*Dx1[i,i])
+        U1[i]=math.exp(gamma*Dx1[i,i])
+
+    for i in range(0,n2*n2):
+        #U2[i]=1/(1-gamma*Dx2[i,i])
+        U2[i]=math.exp(gamma*Dx2[i,i])
+
+
+    #normalization step
+    K12= Left.dot(np.multiply(U,Right))
+    K1= Left1.dot(np.multiply(U1,Right1))
+    K2= Left2.dot(np.multiply(U2,Right2))
+
+
+    K=K12/(math.sqrt(K1)*math.sqrt(K2))
+
+    
+    return K
+
+