diff --git a/fastlib2/mlpack/kde/dataset_scaler.h b/fastlib2/mlpack/kde/dataset_scaler.h
index 912a9bc998..dc6cbdccd2 100644
--- a/fastlib2/mlpack/kde/dataset_scaler.h
+++ b/fastlib2/mlpack/kde/dataset_scaler.h
@@ -86,29 +86,28 @@ class DatasetScaler {
   static void ScaleDataByMinMax(Matrix &qset, Matrix &rset,
 				bool queries_equal_references) {
     
-    int num_dims = rset.n_rows();
-    DHrectBound<2> qset_bound;
-    DHrectBound<2> rset_bound;
-    qset_bound.Init(qset.n_rows());
-    rset_bound.Init(qset.n_rows());
+    index_t num_dims = qset.n_rows();
+    DHrectBound<2> total_bound;
+    total_bound.Init(qset.n_rows());
 
     // go through each query/reference point to find out the bounds
     for(index_t r = 0; r < rset.n_cols(); r++) {
       Vector ref_vector;
       rset.MakeColumnVector(r, &ref_vector);
-      rset_bound |= ref_vector;
+      total_bound |= ref_vector;
     }
-    for(index_t q = 0; q < qset.n_cols(); q++) {
-      Vector query_vector;
-      qset.MakeColumnVector(q, &query_vector);
-      qset_bound |= query_vector;
+    if(!queries_equal_references) {
+      for(index_t q = 0; q < qset.n_cols(); q++) {
+	Vector query_vector;
+	qset.MakeColumnVector(q, &query_vector);
+	total_bound |= query_vector;
+      }
     }
 
     for(index_t i = 0; i < num_dims; i++) {
-      DRange qset_range = qset_bound.get(i);
-      DRange rset_range = rset_bound.get(i);
-      double min_coord = min(qset_range.lo, rset_range.lo);
-      double max_coord = max(qset_range.hi, rset_range.hi);
+      DRange total_range = total_bound.get(i);
+      double min_coord = total_range.lo;
+      double max_coord = total_range.hi;
       double width = max_coord - min_coord;
 
       printf("Dimension %d range: [%g, %g]\n", i, min_coord, max_coord);
@@ -135,6 +134,86 @@ class DatasetScaler {
     }
   }
 
+  /** @brief Standardize the given query and the reference datasets in
+   *         each dimension to have zero mean and at most unit
+   *         variance.
+   *
+   *         Assumes that the query and the reference together contain
+   *         more than one instance.
+   *
+   *  @param qset The column-oriented query set.
+   *  @param rset The column-oriented reference set.
+   *  @param queries_equal_references The boolean flag that tells whether
+   *                                  the queries equal the references.
+   */
+  static void StandardizeData(Matrix &qset, Matrix &rset,
+			      bool queries_equal_references) {
+
+    Vector mean_vector, standard_deviation_vector;
+
+    mean_vector.Init(qset.n_rows());
+    mean_vector.SetZero();
+    standard_deviation_vector.Init(qset.n_rows());
+    standard_deviation_vector.SetZero();
+
+    // Go through each query/reference point to find out the mean
+    // vectors.
+    for(index_t r = 0; r < rset.n_cols(); r++) {
+      la::AddTo(rset.n_rows(), rset.GetColumnPtr(r), mean_vector.ptr());
+    }
+    if(!queries_equal_references) {
+      for(index_t q = 0; q < qset.n_cols(); q++) {
+	la::AddTo(qset.n_rows(), qset.GetColumnPtr(q), mean_vector.ptr());
+      }
+      la::Scale(qset.n_rows(), 1.0 / ((double) qset.n_cols() + rset.n_cols()),
+		mean_vector.ptr());
+    }
+    else {
+      la::Scale(qset.n_rows(), 1.0 / ((double) qset.n_cols()),
+		mean_vector.ptr());
+    }
+
+    // Now find out the standard deviation along each dimension.
+    for(index_t r = 0; r < rset.n_cols(); r++) {
+      for(index_t i = 0; i < rset.n_rows(); i++) {
+	standard_deviation_vector[i] += 
+	  math::Sqr(rset.get(i, r) - mean_vector[i]);
+      }
+    }
+    if(!queries_equal_references) {
+      for(index_t q = 0; q < qset.n_cols(); q++) {
+	for(index_t i = 0; i < qset.n_rows(); i++) {
+	  standard_deviation_vector[i] +=
+	    math::Sqr(qset.get(i, q) - mean_vector[i]);
+	}
+      }
+      la::Scale(qset.n_rows(), 
+		1.0 / ((double) qset.n_cols() + rset.n_cols() - 1),
+		standard_deviation_vector.ptr());
+    }
+    else {
+      la::Scale(rset.n_rows(), 1.0 / ((double) rset.n_cols()),
+		standard_deviation_vector.ptr());
+    }
+
+    // Now scale the datasets using the computed mean and the standard
+    // deviation.
+    for(index_t r = 0; r < rset.n_cols(); r++) {
+      for(index_t d = 0; d < rset.n_rows(); d++) {
+	rset.set(d, r, (rset.get(d, r) - mean_vector[d]) / 
+		 standard_deviation_vector[d]);
+      }
+    }
+    if(!queries_equal_references) {
+      for(index_t q = 0; q < qset.n_cols(); q++) {
+	for(index_t d = 0; d < qset.n_rows(); d++) {
+	  qset.set(d, q, (qset.get(d, q) - mean_vector[d]) /
+		   standard_deviation_vector[d]);
+	}
+      }
+    }
+  }
+
 };
 
 #endif
diff --git a/fastlib2/mlpack/kde/dualtree_kde_main.cc b/fastlib2/mlpack/kde/dualtree_kde_main.cc
index 98df5352b1..b49a35ac48 100644
--- a/fastlib2/mlpack/kde/dualtree_kde_main.cc
+++ b/fastlib2/mlpack/kde/dualtree_kde_main.cc
@@ -47,8 +47,11 @@
  * has to be positive.
  *
  * 5. kde/scaling (optional): whether to prescale the dataset 
- * - range: scales both the query and the reference sets to be within the 
- *   unit hypercube [0, 1]^D where D is the dimensionality. 
+ *
+ * - range: scales both the query and the reference sets to be within
+ * the unit hypercube [0, 1]^D where D is the dimensionality.
+ * - standardize: scales both the query and the reference set to have
+ * zero mean and unit variance.
  * - none: default value; no scaling
  *
  * 6. kde/multiplicative_expansion (optional): If this flag is
@@ -117,6 +120,11 @@ int main(int argc, char *argv[]) {
     DatasetScaler::ScaleDataByMinMax(queries, references,
                                      queries_equal_references);
   }
+  else if(!strcmp(fx_param_str(kde_module, "scaling", "none"), 
+		  "standardize")) {
+    DatasetScaler::StandardizeData(queries, references, 
+				   queries_equal_references);
+  }
 
   if(!strcmp(fx_param_str(kde_module, "kernel", "gaussian"), "gaussian")) {