diff --git a/GPy/examples/GP_regression_demo.py b/GPy/examples/GP_regression_demo.py
index 0fe2fabd..96bd6dc3 100644
--- a/GPy/examples/GP_regression_demo.py
+++ b/GPy/examples/GP_regression_demo.py
@@ -1,7 +1,6 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-
 """
 Simple Gaussian Processes regression with an RBF kernel
 """
@@ -19,8 +18,11 @@ pb.close('all')
 X = np.random.uniform(-3.,3.,(20,1))
 Y = np.sin(X)+np.random.randn(20,1)*0.05
 
+# define kernel
+ker = GPy.kern.rbf(1,ARD=False) + GPy.kern.white(1)
+
 # create simple GP model
-m = GPy.models.GP_regression(X,Y)
+m = GPy.models.GP_regression(X,Y,ker)
 
 # contrain all parameters to be positive
 m.constrain_positive('')
@@ -30,6 +32,7 @@ m.optimize('tnc', max_f_eval = 1000)
 m.plot()
 print(m)
 
+
 ######################################
 ## 2 dimensional example
 
@@ -37,8 +40,11 @@ print(m)
 X = np.random.uniform(-3.,3.,(40,2))
 Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(40,1)*0.05
 
+# define kernel
+ker = GPy.kern.rbf(2,ARD=True) + GPy.kern.white(2)
+
 # create simple GP model
-m = GPy.models.GP_regression(X,Y)
+m = GPy.models.GP_regression(X,Y,ker)
 
 # contrain all parameters to be positive
 m.constrain_positive('')
diff --git a/GPy/examples/sparse_GPLVM_demo.py b/GPy/examples/sparse_GPLVM_demo.py
index 6ca6c941..3f1969fe 100644
--- a/GPy/examples/sparse_GPLVM_demo.py
+++ b/GPy/examples/sparse_GPLVM_demo.py
@@ -10,11 +10,11 @@ print "sparse GPLVM with RBF kernel"
 
 N = 100
 M = 4
-Q = 1
+Q = 2
 D = 2
 #generate GPLVM-like data
 X = np.random.rand(N, Q)
-k = GPy.kern.rbf(Q, 1.0, 2.0) + GPy.kern.white(Q, 0.00001)
+k = GPy.kern.rbf(Q,1.,2*np.ones((1,))) + GPy.kern.white(Q, 0.00001)
 K = k.K(X)
 Y = np.random.multivariate_normal(np.zeros(N),K,D).T
 
diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py
index 0ddc09e3..e1304f11 100644
--- a/GPy/kern/constructors.py
+++ b/GPy/kern/constructors.py
@@ -22,7 +22,7 @@ from Brownian import Brownian as Brownianpart
 #using meta-classes to make the objects construct properly wthout them.
 
 
-def rbf(D,variance=1., lengthscale=1.):
+def rbf(D,variance=1., lengthscale=None,ARD=False):
     """
     Construct an RBF kernel
 
@@ -33,7 +33,7 @@ def rbf(D,variance=1., lengthscale=1.):
     :param lengthscale: the lengthscale of the kernel
     :type lengthscale: float
     """
-    part = rbfpart(D,variance,lengthscale)
+    part = rbfpart(D,variance,lengthscale,ARD)
     return kern(D, [part])
 
 def rbf_ARD(D,variance=1., lengthscales=None):
diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py
index 9e2bb509..775bb797 100644
--- a/GPy/kern/rbf.py
+++ b/GPy/kern/rbf.py
@@ -20,17 +20,33 @@ class rbf(kernpart):
     :type D: int
     :param variance: the variance of the kernel
     :type variance: float
-    :param lengthscale: the lengthscale of the kernel
-    :type lengthscale: float
+    :param lengthscale: the vector of lengthscale of the kernel
+    :type lengthscale: np.ndarray
+    :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension.
+    :type ARD: Boolean
 
-    .. Note: for rbf with different lengthscale on each dimension, see rbf_ARD
     """
 
-    def __init__(self,D,variance=1.,lengthscale=1.):
+    def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
-        self.Nparam = 2
-        self.name = 'rbf'
-        self.set_param(np.hstack((variance,lengthscale)))
+        self.ARD = ARD
+        if ARD == False:
+            self.Nparam = 2
+            self.name = 'rbf'
+            if lengthscale is not None:
+                assert lengthscale.shape == (1,)
+            else:
+                lengthscale = np.ones(1)
+        
+        else:
+            self.Nparam = self.D + 1
+            self.name = 'rbf_ARD'
+            if lengthscale is not None:
+                assert lengthscale.shape == (self.D,)
+            else:
+                lengthscale = np.ones(self.D)
+
+        self.set_param(np.hstack((variance,lengthscale)))        
 
         #initialize cache
         self._Z, self._mu, self._S = np.empty(shape=(3,1))
@@ -40,14 +56,19 @@ class rbf(kernpart):
         return np.hstack((self.variance,self.lengthscale))
 
     def set_param(self,x):
-        self.variance, self.lengthscale = x
+        assert x.size==(self.Nparam)
+        self.variance = x[0]
+        self.lengthscale = x[1:]
         self.lengthscale2 = np.square(self.lengthscale)
         #reset cached results
         self._X, self._X2, self._params = np.empty(shape=(3,1))
         self._Z, self._mu, self._S = np.empty(shape=(3,1)) # cached versions of Z,mu,S
 
     def get_param_names(self):
-        return ['variance','lengthscale']
+        if self.Nparam == 2:
+            return ['variance','lengthscale']
+        else:
+            return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)]        
 
     def K(self,X,X2,target):
         if X2 is None:
@@ -61,7 +82,12 @@ class rbf(kernpart):
     def dK_dtheta(self,partial,X,X2,target):
         self._K_computations(X,X2)
         target[0] += np.sum(self._K_dvar*partial)
-        target[1] += np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
+        if self.ARD == True:
+            dl = self._K_dvar[:,:,None]*self.variance*self._K_dist2/self.lengthscale
+            target[1:] += (dl*partial[:,:,None]).sum(0).sum(0)
+        else:
+            target[1] += np.sum(self._K_dvar*self.variance*(self._K_dist2.sum(-1))/self.lengthscale*partial)
+        #np.sum(self._K_dvar*self.variance*self._K_dist2/self.lengthscale*partial)
 
     def dKdiag_dtheta(self,partial,X,target):
         #NB: derivative of diagonal elements wrt lengthscale is 0
@@ -81,15 +107,13 @@ class rbf(kernpart):
             self._X = X
             self._X2 = X2
             if X2 is None: X2 = X
-            XXT = np.dot(X,X2.T)
-            if X is X2:
-                self._K_dist2 = (-2.*XXT + np.diag(XXT)[:,np.newaxis] + np.diag(XXT)[np.newaxis,:])/self.lengthscale2
-            else:
-                self._K_dist2 = (-2.*XXT + np.sum(np.square(X),1)[:,None] + np.sum(np.square(X2),1)[None,:])/self.lengthscale2
-            # TODO Remove comments if this is fine.
-            # Commented out by Neil as doesn't seem to be used elsewhere.
-            #self._K_exponent = -0.5*self._K_dist2
-            self._K_dvar = np.exp(-0.5*self._K_dist2)
+            self._K_dist = X[:,None,:]-X2[None,:,:] # this can be computationally heavy
+            self._params = np.empty(shape=(1,0))#ensure the next section gets called
+        if not np.all(self._params == self.get_param()):
+            self._params == self.get_param()
+            self._K_dist2 = np.square(self._K_dist/self.lengthscale) 
+            #self._K_exponent = -0.5*self._K_dist2.sum(-1) #ND: commented out because seems not to be used
+            self._K_dvar = np.exp(-0.5*self._K_dist2.sum(-1))
 
     def psi0(self,Z,mu,S,target):
         target += self.variance
@@ -132,7 +156,7 @@ class rbf(kernpart):
         d_length = self._psi2[:,:,:,None]*(0.5*self._psi2_Zdist_sq*self._psi2_denom + 2.*self._psi2_mudist_sq + 2.*S[:,None,None,:]/self.lengthscale2)/(self.lengthscale*self._psi2_denom)
         d_length = d_length.sum(0)
         target[0] += np.sum(partial*d_var)
-        target[1] += np.sum(d_length*partial)
+        target[1:] += (d_length*partial[:,:,None]).sum(0).sum(0)
 
     def dpsi2_dZ(self,partial,Z,mu,S,target):
         """Returns shape N,M,M,Q"""
@@ -175,4 +199,3 @@ class rbf(kernpart):
             self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M
 
             self._Z, self._mu, self._S = Z, mu,S
-