diff --git a/src/shogun/mathematics/Statistics.cpp b/src/shogun/mathematics/Statistics.cpp
index 859dabc1878..d79be5e5535 100644
--- a/src/shogun/mathematics/Statistics.cpp
+++ b/src/shogun/mathematics/Statistics.cpp
@@ -305,30 +305,34 @@ SGVector<float64_t> CStatistics::matrix_std_deviation(
 	return var;
 }
 
-#ifdef HAVE_LAPACK
 SGMatrix<float64_t> CStatistics::covariance_matrix(
 		SGMatrix<float64_t> observations, bool in_place)
 {
-	SGMatrix<float64_t> centered=
-			in_place ?
-					observations :
-					SGMatrix<float64_t>(observations.num_rows,
-							observations.num_cols);
+	int32_t N = observations.num_rows;
+	int32_t D = observations.num_cols;
+
+	REQUIRE(N>1, "Number of observations (%d) must be at least 2.\n", N);
+	REQUIRE(D>0, "Number of dimensions (%d) must be at least 1.\n", D);
 
+	/* center observations, potentially in-place */
+	SGMatrix<float64_t> centered;
 	if (!in_place)
 	{
-		memcpy(centered.matrix, observations.matrix,
-				sizeof(float64_t)*observations.num_rows*observations.num_cols);
+		centered = observations.clone();
 	}
-	centered.remove_column_mean();
+	else
+		centered = observations;
 
-	/* compute 1/(m-1) * X' * X */
-	SGMatrix<float64_t> cov=SGMatrix<float64_t>::matrix_multiply(centered,
-			centered, true, false, 1.0/(observations.num_rows-1));
+	Map<MatrixXd> eigen_centered(centered.matrix, N, D);
+	eigen_centered.rowwise() -= eigen_centered.colwise().mean();
+
+	/* compute and store 1/(N-1) * X.T * X */
+	SGMatrix<float64_t> cov(D, D);
+	Map<MatrixXd> eigen_cov(cov.matrix, D, D);
+	eigen_cov = (eigen_centered.adjoint() * eigen_centered) / double(N - 1);
 
 	return cov;
 }
-#endif //HAVE_LAPACK
 
 float64_t CStatistics::confidence_intervals_mean(SGVector<float64_t> values,
 		float64_t alpha, float64_t& conf_int_low, float64_t& conf_int_up)
diff --git a/src/shogun/mathematics/Statistics.h b/src/shogun/mathematics/Statistics.h
index 99bd613a0ab..13084cbc3cc 100644
--- a/src/shogun/mathematics/Statistics.h
+++ b/src/shogun/mathematics/Statistics.h
@@ -163,7 +163,6 @@ class CStatistics: public CSGObject
 	static SGVector<float64_t> matrix_std_deviation(
 			SGMatrix<float64_t> values, bool col_wise=true);
 
-#ifdef HAVE_LAPACK
 	/** Computes the empirical estimate of the covariance matrix of the given
 	 * data which is organized as num_cols variables with num_rows observations.
 	 *
@@ -174,8 +173,6 @@ class CStatistics: public CSGObject
 	 * \f$\bar X\f$. Then \f$\text{cov}(X)=(X-\bar X)^T(X - \bar X)\f$ is
 	 * returned.
 	 *
-	 * Needs SHOGUN to be compiled with LAPACK.
-	 *
 	 * @param observations data matrix organized as one variable per column
 	 * @param in_place optional, if set to true, observations matrix will be
 	 * centered, if false, a copy will be created an centered.
@@ -183,7 +180,6 @@ class CStatistics: public CSGObject
 	 */
 	static SGMatrix<float64_t> covariance_matrix(
 			SGMatrix<float64_t> observations, bool in_place=false);
-#endif //HAVE_LAPACK
 
 	/** Calculates the sample mean of a given set of samples and also computes
 	 * the confidence interval for the actual mean for a given p-value,