diff --git a/lib/AsyncLinUCB_AM.py b/lib/AsyncLinUCB_AM.py
index b49e496..beb8c86 100644
--- a/lib/AsyncLinUCB_AM.py
+++ b/lib/AsyncLinUCB_AM.py
@@ -49,9 +49,9 @@ def getUCB(self, alpha, article_FeatureVector):
         article_FeatureVector_g = article_FeatureVector[:self.d_g]
         article_FeatureVector_l = article_FeatureVector[self.d_g:]
         mean_g = np.dot(self.UserTheta_g, article_FeatureVector_g)
-        var_g = np.sqrt(np.dot(np.dot(article_FeatureVector_g, self.AInv_g), article_FeatureVector_g))
+        var_g = np.sqrt(np.linalg.multi_dot([article_FeatureVector_g, self.AInv_g, article_FeatureVector_g]))
         mean_l = np.dot(self.UserTheta_l, article_FeatureVector_l)
-        var_l = np.sqrt(np.dot(np.dot(article_FeatureVector_l, self.AInv_l), article_FeatureVector_l))
+        var_l = np.sqrt(np.linalg.multi_dot([article_FeatureVector_l, self.AInv_l, article_FeatureVector_l]))
         pta = mean_g + alpha_g * var_g + mean_l + alpha_l * var_l
         return pta
 
diff --git a/lib/CLUB.py b/lib/CLUB.py
index 17dbe73..943f86b 100644
--- a/lib/CLUB.py
+++ b/lib/CLUB.py
@@ -41,7 +41,7 @@ def updateParametersofClusters(self,clusters,userID,Graph,users):
 
 	def getProb(self, alpha, article_FeatureVector,time):
 		mean = np.dot(self.CTheta, article_FeatureVector)
-		var = np.sqrt(np.dot(np.dot(article_FeatureVector, self.CAInv),  article_FeatureVector))
+		var = np.sqrt(np.linalg.multi_dot([article_FeatureVector, self.CAInv,  article_FeatureVector]))
 		pta = mean +  alpha * var*np.sqrt(math.log10(time+1))
 		return pta
 
diff --git a/lib/CoDBand.py b/lib/CoDBand.py
index 9f52b1d..af85ba6 100644
--- a/lib/CoDBand.py
+++ b/lib/CoDBand.py
@@ -147,7 +147,7 @@ def updateLocalUserModel(self, articlePicked_FeatureVector, click):
 			self.lambda_)
 
 	def getCB(self, x):
-		var = np.sqrt(np.dot(np.dot(x, self.AInv), x))
+		var = np.sqrt(np.linalg.multi_dot([x, self.AInv, x]))
 		if self.alpha != -1:
 			return self.alpha * var
 		else:
@@ -382,4 +382,4 @@ def getTheta(self, userID):
 		return self.globalModels[self.userID2globalModelIndex[userID]].Mean
 
 if __name__ == '__main__':
-	pass
\ No newline at end of file
+	pass
diff --git a/lib/DisLinUCB.py b/lib/DisLinUCB.py
index a83dadd..99a5d00 100644
--- a/lib/DisLinUCB.py
+++ b/lib/DisLinUCB.py
@@ -33,7 +33,7 @@ def getUCB(self, alpha, article_FeatureVector):
             alpha = self.alpha_t
 
         mean = np.dot(self.UserTheta, article_FeatureVector)
-        var = np.sqrt(np.dot(np.dot(article_FeatureVector, self.AInv), article_FeatureVector))
+        var = np.sqrt(np.linalg.multi_dot([article_FeatureVector, self.AInv, article_FeatureVector]))
         pta = mean + alpha * var
         return pta
 
diff --git a/lib/DyClu.py b/lib/DyClu.py
index f84afd5..49b83a2 100644
--- a/lib/DyClu.py
+++ b/lib/DyClu.py
@@ -84,7 +84,7 @@ def updateLocalUserModel(self, articlePicked_FeatureVector, click):
             self.lambda_)
 
     def getCB(self, x, useConstantAlpha = False):
-        var = np.sqrt(np.dot(np.dot(x, self.AInv), x))
+        var = np.sqrt(np.linalg.multi_dot([x, self.AInv, x]))
         if useConstantAlpha:
             return self.alpha * var
         else:
@@ -109,7 +109,7 @@ def getInstantaneousBadness(self, articlePicked, click, method="ChiSquare"):
                 else:
                     mean = np.dot(self.UserThetaNoReg, x)
                     rewardEstimationError = (mean - click)**2
-                    rewardEstimationErrorSTD = self.NoiseScale**2 * (1 + np.dot(np.dot(x, np.linalg.pinv(self.A-self.lambda_ * np.identity(n=self.d))), x))
+                    rewardEstimationErrorSTD = self.NoiseScale**2 * (1 + np.linalg.multi_dot([x, np.linalg.pinv(self.A-self.lambda_ * np.identity(n=self.d)), x]))
                     df1 = 1
 
                     chiSquareStatistic = rewardEstimationError / rewardEstimationErrorSTD
@@ -164,7 +164,7 @@ def __init__(self, cluster):
             self.lambda_)
 
     def getCB(self, x, useConstantAlpha = False):
-        var = np.sqrt(np.dot(np.dot(x, self.AInv), x))
+        var = np.sqrt(np.linalg.multi_dot([x, self.AInv, x]))
         if useConstantAlpha:
             return self.alpha * var
         else:
@@ -400,4 +400,4 @@ def resetUserConnectedness(self, userID):
         self.user_graph[0:self.cur_userIndex, self.userID2userIndex[userID]] = 1
 
     def getTheta(self, userID):
-        return self.users[userID].UserTheta
\ No newline at end of file
+        return self.users[userID].UserTheta
diff --git a/lib/LinUCB.py b/lib/LinUCB.py
index a382173..d01ec3c 100644
--- a/lib/LinUCB.py
+++ b/lib/LinUCB.py
@@ -208,10 +208,10 @@ def updateParameters(self, articlePicked_FeatureVector, click, userID):
 	def getProb(self, alpha, article_FeatureVector,userID):
 		x = article_FeatureVector
 		z = vectorize(np.outer(self.users[userID].userFeature, article_FeatureVector))
-		temp =np.dot(np.dot(np.dot( self.A_zInv , np.transpose( self.users[userID].B)) , self.users[userID].AInv), x )
+		temp =np.linalg.multi_dot([self.A_zInv , np.transpose( self.users[userID].B) , self.users[userID].AInv, x ])
 		mean = np.dot(self.users[userID].UserTheta,  x)+ np.dot(self.beta, z)
-		s_t = np.dot(np.dot(z, self.A_zInv),  z) + np.dot(np.dot(x, self.users[userID].AInv),  x)
-		-2* np.dot(z, temp)+ np.dot(np.dot( np.dot(x, self.users[userID].AInv) ,  self.users[userID].B ) ,temp)
+		s_t = np.linalg.multi_dot([z, self.A_zInv,  z]) + np.linalg.multi_dot([x, self.users[userID].AInv,  x])
+		-2* np.dot(z, temp)+ np.linalg.multi_dot([x, self.users[userID].AInv ,  self.users[userID].B  ,temp])
 
 		var = np.sqrt(s_t)
 		pta = mean + alpha * var
diff --git a/lib/dLinUCB.py b/lib/dLinUCB.py
index 77a8b63..0a14ab1 100644
--- a/lib/dLinUCB.py
+++ b/lib/dLinUCB.py
@@ -45,14 +45,14 @@ def getProbInfo(self, alpha, article_FeatureVector):
 		if alpha == -1:
 			alpha = alpha = 0.1*np.sqrt(np.log(self.time+1))
 		mean = np.dot(self.UserTheta,  article_FeatureVector)
-		var = np.sqrt(np.dot(np.dot(article_FeatureVector, self.AInv),  article_FeatureVector))
+		var = np.sqrt(np.linalg.multi_dot([article_FeatureVector, self.AInv,  article_FeatureVector]))
 		self.alpha_t = self.NoiseScale *np.sqrt(np.log(np.linalg.det(self.A)/float(self.sigma * self.lambda_) )) + np.sqrt(self.lambda_)
 		#pta = mean + alpha * var
 		return {'mean':mean, 'var':var, 'alpha':alpha, 'alpha_t':self.alpha_t}
 
 	def getProb_plot(self, alpha, article_FeatureVector):
 		mean = np.dot(self.UserTheta,  article_FeatureVector)
-		var = np.sqrt(np.dot(np.dot(article_FeatureVector, self.AInv),  article_FeatureVector))
+		var = np.sqrt(np.linalg.multi_dot([article_FeatureVector, self.AInv,  article_FeatureVector]))
 		pta = mean + alpha * var
 		return pta, mean, alpha * var
 
@@ -154,7 +154,7 @@ def getSlavePredictionInfo(self, alpha, article_FeatureVector):
         if alpha == -1:
             alpha = alpha = 0.1*np.sqrt(np.log(self.time+1))
         mean = np.dot(self.UserTheta,  article_FeatureVector)
-        var = np.sqrt(np.dot(np.dot(article_FeatureVector, self.AInv),  article_FeatureVector))
+        var = np.sqrt(np.linalg.multi_dot([article_FeatureVector, self.AInv,  article_FeatureVector]))
         # self.alpha_t =  self.NoiseScale*np.sqrt(self.d* np.log( (self.lambda_ + self.update_num)/float(self.sigma * self.lambda_) )) + np.sqrt(self.lambda_)
         self.alpha_t = self.NoiseScale ** 2 * np.sqrt(
             self.d * np.log(1 + self.update_num / (self.d * self.lambda_)) + 2 * np.log(1 / self.delta_1)) + np.sqrt(
diff --git a/lib/factorUCB.py b/lib/factorUCB.py
index 6f0da10..53b08f4 100644
--- a/lib/factorUCB.py
+++ b/lib/factorUCB.py
@@ -82,7 +82,7 @@ def updateParameters(self, articles, clicks, userID):
 
 		self.UserTheta = matrixize(np.dot(self.AInv, self.b), len(articles[0].V)) 
 		self.CoTheta = np.dot(self.UserTheta, self.W)
-		self.CCA = np.dot(np.dot(self.BigW , self.AInv), np.transpose(self.BigW))				
+		self.CCA = np.linalg.multi_dot([self.BigW , self.AInv, np.transpose(self.BigW)])				
 	
 	def getA(self):
 		return self.A
@@ -97,8 +97,8 @@ def getProb(self, alpha, alpha2, article, userID):
 		TempFeatureV = vectorize(TempFeatureM)
 
 		mean = np.dot(self.CoTheta.T[userID], article.V)	
-		var = np.sqrt(np.dot(np.dot(TempFeatureV, self.CCA), TempFeatureV))
-		var2 = np.sqrt(np.dot(np.dot(self.CoTheta.T[userID][self.context_dimension:], article.A2Inv),  self.CoTheta.T[userID][self.context_dimension:]))
+		var = np.sqrt(np.linalg.multi_dot([TempFeatureV, self.CCA, TempFeatureV]))
+		var2 = np.sqrt(np.linalg.multi_dot([self.CoTheta.T[userID][self.context_dimension:], article.A2Inv,  self.CoTheta.T[userID][self.context_dimension:]]))
 		pta = mean + alpha * var + alpha2*var2
 		return pta
 	def getProb_plot(self, alpha, alpha2, article, userID):