diff --git a/src/shogun/statistical_testing/MMD.cpp b/src/shogun/statistical_testing/MMD.cpp
index 9c61d0434ae..217ec5ed7fa 100644
--- a/src/shogun/statistical_testing/MMD.cpp
+++ b/src/shogun/statistical_testing/MMD.cpp
@@ -89,6 +89,7 @@ void CMMD::Self::create_statistic_job()
 {
 	const DataManager& dm = owner.get_data_manager();
 	auto Bx = dm.blocksize_at(0);
+	auto By = dm.blocksize_at(1);
 	switch (statistic_type)
 	{
 		case EStatisticType::UNBIASED_FULL:
@@ -102,7 +103,7 @@ void CMMD::Self::create_statistic_job()
 			break;
 		default : break;
 	};
-	permutation_job = mmd::WithinBlockPermutation(Bx, statistic_type);
+	permutation_job = mmd::WithinBlockPermutation(Bx, By, statistic_type);
 }
 
 void CMMD::Self::create_variance_job()
diff --git a/src/shogun/statistical_testing/QuadraticTimeMMD.cpp b/src/shogun/statistical_testing/QuadraticTimeMMD.cpp
index 0d1ae1780d0..e55446d0e43 100644
--- a/src/shogun/statistical_testing/QuadraticTimeMMD.cpp
+++ b/src/shogun/statistical_testing/QuadraticTimeMMD.cpp
@@ -89,6 +89,7 @@ void CQuadraticTimeMMD::Self::create_statistic_job()
 	SG_SDEBUG("Entering\n");
 	const DataManager& dm=owner.get_data_manager();
 	auto Nx=dm.num_samples_at(0);
+	auto Ny=dm.num_samples_at(1);
 	switch (owner.get_statistic_type())
 	{
 		case EStatisticType::UNBIASED_FULL:
@@ -102,7 +103,7 @@ void CQuadraticTimeMMD::Self::create_statistic_job()
 			break;
 		default : break;
 	};
-	permutation_job = mmd::WithinBlockPermutation(Nx, owner.get_statistic_type());
+	permutation_job = mmd::WithinBlockPermutation(Nx, Ny, owner.get_statistic_type());
 	SG_SDEBUG("Leaving\n");
 }
 
diff --git a/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.cpp b/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.cpp
index 1b0d16fb3bb..e21752d81f8 100644
--- a/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.cpp
+++ b/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.cpp
@@ -16,6 +16,7 @@
  * along with this program.  If not, see <http:/www.gnu.org/licenses/>.
  */
 
+#include <unordered_map>
 #include <shogun/io/SGIO.h>
 #include <shogun/lib/SGMatrix.h>
 #include <shogun/lib/SGVector.h>
@@ -32,84 +33,89 @@ using namespace shogun;
 using namespace internal;
 using namespace mmd;
 
-WithinBlockPermutation::WithinBlockPermutation(index_t n, EStatisticType type)
-: n_x(n), stype(type)
+WithinBlockPermutation::WithinBlockPermutation(index_t nx, index_t ny, EStatisticType type)
+: n_x(nx), n_y(ny), stype(type), terms()
 {
+	SG_SDEBUG("number of samples are %d and %d!\n", n_x, n_y);
+	std::fill(&terms.term[0], &terms.term[2]+1, 0);
+	std::fill(&terms.diag[0], &terms.diag[2]+1, 0);
+}
+
+void WithinBlockPermutation::add_term(float64_t val, index_t i, index_t j)
+{
+	if (i<n_x && j<n_x)
+	{
+		SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 1!\n", i, j, val);
+		terms.term[0]+=val;
+		if (i==j)
+			terms.diag[0]+=val;
+	}
+	else if (i>=n_x && j>=n_x)
+	{
+		SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 2!\n", i, j, val);
+		terms.term[1]+=val;
+		if (i==j)
+			terms.diag[1]+=val;
+	}
+	else if (i>=n_x && j<n_x)
+	{
+		SG_SDEBUG("Adding KernelMatrix(%d,%d)=%f to term 3!\n", i, j, val);
+		terms.term[2]+=val;
+		if (i-n_x==j)
+			terms.diag[2]+=val;
+	}
 }
 
 float64_t WithinBlockPermutation::operator()(SGMatrix<float64_t> km)
 {
 	SG_SDEBUG("Entering!\n");
-	inds.resize(km.num_rows);
-	std::iota(inds.data(), inds.data()+inds.size(), 0);
-	SGVector<index_t> permuted_inds(inds.data(), inds.size(), false);
+	SGVector<index_t> permuted_inds(n_x+n_y);
+	std::iota(permuted_inds.vector, permuted_inds.vector+permuted_inds.vlen, 0);
 	CMath::permute(permuted_inds);
 
-	const index_t n_y=km.num_rows-n_x;
-	SG_SDEBUG("number of samples are %d and %d!\n", n_x, n_y);
+	std::unordered_map<index_t, index_t> inds;
+	for (int i=0; i<permuted_inds.vlen; ++i)
+		inds.insert(std::make_pair(permuted_inds[i], i));
 
-	auto term_1=0.0;
-	for (auto i=0; i<n_x; ++i)
-	{
-		for (auto j=0; j<n_x; ++j)
-		{
-			if (i>j)
-				term_1+=km(permuted_inds[i], permuted_inds[j]);
-		}
-	}
-	term_1*=2;
-	SG_SDEBUG("term_1 sum (without diagonal) = %f!\n", term_1);
-	if (stype==EStatisticType::BIASED_FULL)
-	{
-		for (auto i=0; i<n_x; ++i)
-			term_1+=km(permuted_inds[i], permuted_inds[i]);
-		SG_SDEBUG("term_1 sum (with diagonal) = %f!\n", term_1);
-		term_1/=n_x*n_x;
-	}
-	else
-		term_1/=n_x*(n_x-1);
-	SG_SDEBUG("term_1 (normalized) = %f!\n", term_1);
+	std::fill(&terms.term[0], &terms.term[2]+1, 0);
+	std::fill(&terms.diag[0], &terms.diag[2]+1, 0);
 
-	auto term_2=0.0;
-	for (auto i=n_x; i<n_x+n_y; ++i)
+	for (auto j=0; j<n_x+n_y; ++j)
 	{
-		for (auto j=n_x; j<n_x+n_y; ++j)
-		{
-			if (i>j)
-				term_2+=km(permuted_inds[i], permuted_inds[j]);
-		}
+		for (auto i=0; i<n_x+n_y; ++i)
+			add_term(km(i, j), inds.find(i)->second, inds.find(j)->second);
 	}
-	term_2*=2.0;
-	SG_SDEBUG("term_2 sum (without diagonal) = %f!\n", term_2);
-	if (stype==EStatisticType::BIASED_FULL)
+
+	SG_SDEBUG("term_1 sum (with diagonal) = %f!\n", terms.term[0]);
+	SG_SDEBUG("term_2 sum (with diagonal) = %f!\n", terms.term[1]);
+	if (stype!=EStatisticType::BIASED_FULL)
 	{
-		for (auto i=n_x; i<n_x+n_y; ++i)
-			term_2+=km(permuted_inds[i], permuted_inds[i]);
-		SG_SDEBUG("term_2 sum (with diagonal) = %f!\n", term_2);
-		term_2/=n_y*n_y;
+		terms.term[0]-=terms.diag[0];
+		terms.term[1]-=terms.diag[1];
+		SG_SDEBUG("term_1 sum (without diagonal) = %f!\n", terms.term[0]);
+		SG_SDEBUG("term_2 sum (without diagonal) = %f!\n", terms.term[1]);
+		terms.term[0]/=n_x*(n_x-1);
+		terms.term[1]/=n_y*(n_y-1);
 	}
 	else
-		term_2/=n_y*(n_y-1);
-	SG_SDEBUG("term_2 (normalized) = %f!\n", term_2);
-
-	auto term_3=0.0;
-	for (auto i=n_x; i<n_x+n_y; ++i)
 	{
-		for (auto j=0; j<n_x; ++j)
-			term_3+=km(permuted_inds[i], permuted_inds[j]);
+		terms.term[0]/=n_x*n_x;
+		terms.term[1]/=n_y*n_y;
 	}
-	SG_SDEBUG("term_3 sum (with diagonal) = %f!\n", term_3);
+	SG_SDEBUG("term_1 (normalized) = %f!\n", terms.term[0]);
+	SG_SDEBUG("term_2 (normalized) = %f!\n", terms.term[1]);
+
+	SG_SDEBUG("term_3 sum (with diagonal) = %f!\n", terms.term[2]);
 	if (stype==EStatisticType::UNBIASED_INCOMPLETE)
 	{
-		for (auto i=0; i<n_x; ++i)
-			term_3-=km(permuted_inds[i+n_x], permuted_inds[i]);
-		SG_SDEBUG("term_3 sum (without diagonal) = %f!\n", term_3);
-		term_3/=n_x*(n_x-1);
+		terms.term[2]-=terms.diag[2];
+		SG_SDEBUG("term_3 sum (without diagonal) = %f!\n", terms.term[2]);
+		terms.term[2]/=n_x*(n_x-1);
 	}
 	else
-		term_3/=n_x*n_y;
-	SG_SDEBUG("term_3 (normalized) = %f!\n", term_3);
+		terms.term[2]/=n_x*n_y;
+	SG_SDEBUG("term_3 (normalized) = %f!\n", terms.term[2]);
 
 	SG_SDEBUG("Leaving!\n");
-	return term_1+term_2-2*term_3;
+	return terms.term[0]+terms.term[1]-2*terms.term[2];
 }
diff --git a/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.h b/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.h
index 9cc0ee09e99..135c239b6ba 100644
--- a/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.h
+++ b/src/shogun/statistical_testing/internals/mmd/WithinBlockPermutation.h
@@ -34,15 +34,24 @@ namespace internal
 namespace mmd
 {
 
-struct WithinBlockPermutation
+class WithinBlockPermutation
 {
 	using return_type = float64_t;
-	WithinBlockPermutation(index_t, EStatisticType);
+public:
+	WithinBlockPermutation(index_t, index_t, EStatisticType);
 	return_type operator()(SGMatrix<float64_t> kernel_matrix);
 //	return_type operator()(CGPUMatrix<float64_t> kernel_matrix);
+private:
 	const index_t n_x;
+	const index_t n_y;
 	const EStatisticType stype;
-	std::vector<index_t> inds;
+	void add_term(float64_t, index_t, index_t);
+	struct terms_t
+	{
+		float64_t term[3];
+		float64_t diag[3];
+	};
+	terms_t terms;
 };
 
 }
diff --git a/tests/unit/statistical_testing/internals/WithinBlockPermutation_unittest.cc b/tests/unit/statistical_testing/internals/WithinBlockPermutation_unittest.cc
index 6b206ac6269..1e9a791a997 100644
--- a/tests/unit/statistical_testing/internals/WithinBlockPermutation_unittest.cc
+++ b/tests/unit/statistical_testing/internals/WithinBlockPermutation_unittest.cc
@@ -51,16 +51,25 @@ using namespace Eigen;
 TEST(WithinBlockPermutation, biased_full)
 {
 	const index_t dim=2;
-	const index_t n=6;
-	const index_t m=4;
+	const index_t n=13;
+	const index_t m=7;
 
 	using operation=std::function<float64_t(SGMatrix<float64_t>)>;
 
-	SGMatrix<float64_t> data(dim, n+m);
-	std::iota(data.matrix, data.matrix+dim*(n+m), 1);
+	SGMatrix<float64_t> data_p(dim, n);
+	std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
+	std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });
 
-	auto feats=new CDenseFeatures<float64_t>(data);
+	SGMatrix<float64_t> data_q(dim, m);
+	std::iota(data_q.matrix, data_q.matrix+dim*m, n+1);
+	std::for_each(data_q.matrix, data_q.matrix+dim*m, [&m](float64_t& val) { val/=2*m; });
+
+	auto feats_p=new CDenseFeatures<float64_t>(data_p);
+	auto feats_q=new CDenseFeatures<float64_t>(data_q);
+	auto feats=feats_p->create_merged_copy(feats_q);
 	SG_REF(feats);
+	SG_UNREF(feats_p);
+	SG_UNREF(feats_q);
 
 	auto kernel=some<CGaussianKernel>();
 	kernel->set_width(2.0);
@@ -69,7 +78,7 @@ TEST(WithinBlockPermutation, biased_full)
 	auto mat=kernel->get_kernel_matrix();
 
 	// compute using within-block-permutation functor
-    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::BIASED_FULL);
+    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, m, EStatisticType::BIASED_FULL);
 	sg_rand->set_seed(12345);
 	auto result_1=compute(mat);
 
@@ -101,22 +110,31 @@ TEST(WithinBlockPermutation, biased_full)
 	EXPECT_NEAR(result_1, result_2, 1E-15);
 	EXPECT_NEAR(result_1, result_3, 1E-15);
 
-	kernel->remove_lhs_and_rhs();
+	SG_UNREF(feats);
 }
 
 TEST(WithinBlockPermutation, unbiased_full)
 {
 	const index_t dim=2;
-	const index_t n=3;
-	const index_t m=4;
+	const index_t n=13;
+	const index_t m=7;
 
 	using operation=std::function<float64_t(SGMatrix<float64_t>)>;
 
-	SGMatrix<float64_t> data(dim, n+m);
-	std::iota(data.matrix, data.matrix+dim*(n+m), 1);
+	SGMatrix<float64_t> data_p(dim, n);
+	std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
+	std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });
+
+	SGMatrix<float64_t> data_q(dim, m);
+	std::iota(data_q.matrix, data_q.matrix+dim*m, n+1);
+	std::for_each(data_q.matrix, data_q.matrix+dim*m, [&m](float64_t& val) { val/=2*m; });
 
-	auto feats=new CDenseFeatures<float64_t>(data);
+	auto feats_p=new CDenseFeatures<float64_t>(data_p);
+	auto feats_q=new CDenseFeatures<float64_t>(data_q);
+	auto feats=feats_p->create_merged_copy(feats_q);
 	SG_REF(feats);
+	SG_UNREF(feats_p);
+	SG_UNREF(feats_q);
 
 	auto kernel=some<CGaussianKernel>();
 	kernel->set_width(2.0);
@@ -125,7 +143,7 @@ TEST(WithinBlockPermutation, unbiased_full)
 	auto mat=kernel->get_kernel_matrix();
 
 	// compute using within-block-permutation functor
-    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::UNBIASED_FULL);
+    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, m, EStatisticType::UNBIASED_FULL);
 	sg_rand->set_seed(12345);
 	auto result_1=compute(mat);
 
@@ -157,22 +175,30 @@ TEST(WithinBlockPermutation, unbiased_full)
 	EXPECT_NEAR(result_1, result_2, 1E-15);
 	EXPECT_NEAR(result_1, result_3, 1E-15);
 
-	kernel->remove_lhs_and_rhs();
+	SG_UNREF(feats);
 }
 
 TEST(WithinBlockPermutation, unbiased_incomplete)
 {
 	const index_t dim=2;
-	const index_t n=5;
-	const index_t m=5;
+	const index_t n=10;
 
 	using operation=std::function<float64_t(SGMatrix<float64_t>)>;
 
-	SGMatrix<float64_t> data(dim, n+m);
-	std::iota(data.matrix, data.matrix+dim*(n+m), 1);
+	SGMatrix<float64_t> data_p(dim, n);
+	std::iota(data_p.matrix, data_p.matrix+dim*n, 1);
+	std::for_each(data_p.matrix, data_p.matrix+dim*n, [&n](float64_t& val) { val/=n; });
+
+	SGMatrix<float64_t> data_q(dim, n);
+	std::iota(data_q.matrix, data_q.matrix+dim*n, n+1);
+	std::for_each(data_q.matrix, data_q.matrix+dim*n, [&n](float64_t& val) { val/=2*n; });
 
-	auto feats=new CDenseFeatures<float64_t>(data);
+	auto feats_p=new CDenseFeatures<float64_t>(data_p);
+	auto feats_q=new CDenseFeatures<float64_t>(data_q);
+	auto feats=feats_p->create_merged_copy(feats_q);
 	SG_REF(feats);
+	SG_UNREF(feats_p);
+	SG_UNREF(feats_q);
 
 	auto kernel=some<CGaussianKernel>();
 	kernel->set_width(2.0);
@@ -181,7 +207,7 @@ TEST(WithinBlockPermutation, unbiased_incomplete)
 	auto mat=kernel->get_kernel_matrix();
 
 	// compute using within-block-permutation functor
-    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, EStatisticType::UNBIASED_INCOMPLETE);
+    operation compute=shogun::internal::mmd::WithinBlockPermutation(n, n, EStatisticType::UNBIASED_INCOMPLETE);
 	sg_rand->set_seed(12345);
 	auto result_1=compute(mat);
 
@@ -213,5 +239,5 @@ TEST(WithinBlockPermutation, unbiased_incomplete)
 	EXPECT_NEAR(result_1, result_2, 1E-15);
 	EXPECT_NEAR(result_1, result_3, 1E-15);
 
-	kernel->remove_lhs_and_rhs();
+	SG_UNREF(feats);
 }