diff --git a/bob/math/cpp/eig.cpp b/bob/math/cpp/eig.cpp
index 2adae9941c1f7e9c50ad6a3484d0c9d42ec6107c..93a49ef611ecc886a54f2255e33abf12f9a776c4 100644
--- a/bob/math/cpp/eig.cpp
+++ b/bob/math/cpp/eig.cpp
@@ -57,16 +57,6 @@ void bob::math::eig(const blitz::Array<double,2>& A,
   bob::core::array::assertSameShape(V,shape2);
   bob::core::array::assertSameShape(D,shape1);
 
-  bob::math::eig_(A, V, D);
-}
-
-void bob::math::eig_(const blitz::Array<double,2>& A,
-  blitz::Array<std::complex<double>,2>& V,
-  blitz::Array<std::complex<double>,1>& D)
-{
-  // Size variable
-  const int N = A.extent(0);
-
   // Prepares to call LAPACK function
   // Initialises LAPACK variables
   const char jobvl = 'N'; // Do NOT compute left eigen-vectors
@@ -147,15 +137,6 @@ void bob::math::eigSym(const blitz::Array<double,2>& A,
   bob::core::array::assertSameShape(V,shape2);
   bob::core::array::assertSameShape(D,shape1);
 
-  bob::math::eigSym_(A, V, D);
-}
-
-void bob::math::eigSym_(const blitz::Array<double,2>& A,
-  blitz::Array<double,2>& V, blitz::Array<double,1>& D)
-{
-  // Size variable
-  const int N = A.extent(0);
-
   // Prepares to call LAPACK function
   // Initialises LAPACK variables
   const char jobz = 'V'; // Get both the eigenvalues and the eigenvectors
@@ -234,15 +215,6 @@ void bob::math::eigSym(const blitz::Array<double,2>& A, const blitz::Array<doubl
   bob::core::array::assertSameShape(V,shape2);
   bob::core::array::assertSameShape(D,shape1);
 
-  bob::math::eigSym_(A, B, V, D);
-}
-
-void bob::math::eigSym_(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
-  blitz::Array<double,2>& V, blitz::Array<double,1>& D)
-{
-  // Size variable
-  const int N = A.extent(0);
-
   // Prepares to call LAPACK function
   // Initialises LAPACK variables
   const int itype = 1;
diff --git a/bob/math/cpp/pinv.cpp b/bob/math/cpp/pinv.cpp
index 10f6d9d0453e6f823854894dc5f606e0f7f934d3..4e6ec1f1ef29e2eb29d665264d30d32967f0ef74 100644
--- a/bob/math/cpp/pinv.cpp
+++ b/bob/math/cpp/pinv.cpp
@@ -18,6 +18,7 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
   // Size variables
   const int M = A.extent(0);
   const int N = A.extent(1);
+  const int nb_singular = std::min(M,N);
 
   // Checks zero base
   bob::core::array::assertZeroBase(A);
@@ -26,18 +27,6 @@ void bob::math::pinv(const blitz::Array<double,2>& A,
   bob::core::array::assertSameDimensionLength(B.extent(0), N);
   bob::core::array::assertSameDimensionLength(B.extent(1), M);
 
-  // Calls pinv_()
-  bob::math::pinv_(A, B, rcond);
-}
-
-void bob::math::pinv_(const blitz::Array<double,2>& A,
-  blitz::Array<double,2>& B, const double rcond)
-{
-  // Size variables
-  const int M = A.extent(0);
-  const int N = A.extent(1);
-  const int nb_singular = std::min(M,N);
-
   // Allocates arrays for the SVD
   blitz::Array<double,2> U(M,M);
   blitz::Array<double,1> sigma(nb_singular);
@@ -46,7 +35,7 @@ void bob::math::pinv_(const blitz::Array<double,2>& A,
   blitz::Array<double,2> sigmai_Ut(N, M);
   blitz::Array<double,1> Ut_sum(M);
   // Computes the SVD
-  bob::math::svd_(A, U, sigma, Vt);
+  bob::math::svd(A, U, sigma, Vt);
 
   // Cuts off small sigmas for the inverse similar to Numpy:
   // cf. https://github.com/numpy/numpy/blob/maintenance/1.7.x/numpy/linalg/linalg.py#L1577
@@ -68,4 +57,3 @@ void bob::math::pinv_(const blitz::Array<double,2>& A,
   blitz::Array<double,2> V = Vt.transpose(1,0);
   bob::math::prod(V, sigmai_Ut, B);
 }
-
diff --git a/bob/math/cpp/sqrtm.cpp b/bob/math/cpp/sqrtm.cpp
index c02bcff6e1b7fee2e48fe59e8b42ea85eb0b41bc..09b329a277cba338b8a91e4f3d107242f7f9a995 100644
--- a/bob/math/cpp/sqrtm.cpp
+++ b/bob/math/cpp/sqrtm.cpp
@@ -34,7 +34,7 @@ void bob::math::sqrtSymReal(const blitz::Array<double,2>& A,
   blitz::Array<double,2> Vt = V.transpose(1,0);
   blitz::Array<double,1> D(N);
   blitz::Array<double,2> tmp(N,N); // Cache for multiplication
-  bob::math::eigSym_(A,V,D);
+  bob::math::eigSym(A,V,D);
 
   // 2/ Updates the diagonal matrix D, such that D=sqrt(|D|)
   //    |.| is used to deal with values close to zero (-epsilon)
diff --git a/bob/math/cpp/svd.cpp b/bob/math/cpp/svd.cpp
index 854cd47328f134cabab3c4f28040c879243ab409..9bde71ee4c47a66d747d732f7b2ef9d644e60332 100644
--- a/bob/math/cpp/svd.cpp
+++ b/bob/math/cpp/svd.cpp
@@ -79,10 +79,10 @@ static void svd_lapack( const char jobz, const int M, const int N,
     if (info != 0)
       throw std::runtime_error("The LAPACK dgesdd function returned a non-zero value. You may consider using LAPACK dgsevd instead (see #171) by enabling the 'safe' option.");
   }
-  
+
   // Defining the sign of the eigenvectors
-  // Approch extracted from page 8 - http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf  
-  if(U[0] < 0){    
+  // Approch extracted from page 8 - http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf
+  if(U[0] < 0){
     int ucol=0; ucol= (jobz=='A')? M : std::min(M,N);
     for (int i=0; i<ldu*ucol; i++){
       U[i] = -1*U[i];
@@ -114,17 +114,6 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
   bob::core::array::assertSameDimensionLength(Vt.extent(0), N);
   bob::core::array::assertSameDimensionLength(Vt.extent(1), N);
 
-  bob::math::svd_(A, U, sigma, Vt, safe);
-}
-
-void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
-  blitz::Array<double,1>& sigma, blitz::Array<double,2>& Vt, bool safe)
-{
-  // Size variables
-  const int M = A.extent(0);
-  const int N = A.extent(1);
-  const int nb_singular = std::min(M,N);
-
   // Prepares to call LAPACK function:
   // We will decompose A^T rather than A to reduce the required number of copy
   // We recall that FORTRAN/LAPACK is column-major order whereas blitz arrays
@@ -161,7 +150,7 @@ void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
 
   // Call the LAPACK function
   svd_lapack(jobz, N, M, A_lapack, lda, S_lapack, U_lapack, ldu,
-    VT_lapack, ldvt, safe);  
+    VT_lapack, ldvt, safe);
 
 
   // Copy singular vectors back to U, V and sigma if required
@@ -188,17 +177,6 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
   bob::core::array::assertSameDimensionLength(U.extent(1), nb_singular);
   bob::core::array::assertSameDimensionLength(sigma.extent(0), nb_singular);
 
-  bob::math::svd_(A, U, sigma, safe);
-}
-
-void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
-  blitz::Array<double,1>& sigma, bool safe)
-{
-  // Size variables
-  const int M = A.extent(0);
-  const int N = A.extent(1);
-  const int nb_singular = std::min(M,N);
-
   // Prepares to call LAPACK function
 
   // Initialises LAPACK variables
@@ -249,16 +227,6 @@ void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,1>& sig
   // Checks and resizes if required
   bob::core::array::assertSameDimensionLength(sigma.extent(0), nb_singular);
 
-  bob::math::svd_(A, sigma, safe);
-}
-
-void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,1>& sigma, bool safe)
-{
-  // Size variables
-  const int M = A.extent(0);
-  const int N = A.extent(1);
-  const int nb_singular = std::min(M,N);
-
   // Prepares to call LAPACK function
 
   // Initialises LAPACK variables
diff --git a/bob/math/include/bob.math/eig.h b/bob/math/include/bob.math/eig.h
index f225de09c6bd10111fe3c8a53c04b6df360f9064..f269242fef506826dcc0c8dd4529593d7bfb21ef 100644
--- a/bob/math/include/bob.math/eig.h
+++ b/bob/math/include/bob.math/eig.h
@@ -32,23 +32,6 @@ namespace bob { namespace math {
       blitz::Array<std::complex<double>,2>& V,
       blitz::Array<std::complex<double>,1>& D);
 
-  /**
-   * @brief Function which performs an eigenvalue decomposition of a real
-   * matrix, using the LAPACK function <code>dgeev</code>. This version does
-   * <b>NOT</b> perform any checks on the input data and should be, therefore,
-   * faster.
-   *
-   * @param A The A matrix to decompose (size NxN)
-   *
-   * @param V The V matrix of eigenvectors (size NxN) stored in columns
-   *
-   * @param D The vector of eigenvalues (size N) (in ascending order). Note
-   * that this is a 1D array rather than a 2D diagonal matrix!
-   */
-  void eig_(const blitz::Array<double,2>& A,
-      blitz::Array<std::complex<double>,2>& V,
-      blitz::Array<std::complex<double>,1>& D);
-
   /**
    * @brief Function which performs an eigenvalue decomposition of a real
    * symmetric matrix, using the LAPACK function <code>dsyevd</code>.
@@ -65,24 +48,6 @@ namespace bob { namespace math {
   void eigSym(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
       blitz::Array<double,1>& D);
 
-  /**
-   * @brief Function which performs an eigenvalue decomposition of a real
-   * symmetric matrix, using the LAPACK function <code>dsyevd</code>. This
-   * version does <b>NOT</b> perform any checks on the input data and should
-   * be, therefore, faster.
-   *
-   * @warning The input matrix should be symmetric.
-   *
-   * @param A The A matrix to decompose (size NxN)
-   *
-   * @param V The V matrix of eigenvectors (size NxN) stored in columns
-   *
-   * @param D The vector of eigenvalues (size N) (in ascending order). Note
-   * that this is a 1D array rather than a 2D diagonal matrix!
-   */
-  void eigSym_(const blitz::Array<double,2>& A, blitz::Array<double,2>& V,
-      blitz::Array<double,1>& D);
-
   /**
    * @brief Computes all the eigenvalues and the eigenvectors of a real
    * generalized symmetric-definite eigenproblem, of the form:
@@ -103,28 +68,6 @@ namespace bob { namespace math {
   void eigSym(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
       blitz::Array<double,2>& V, blitz::Array<double,1>& D);
 
-  /**
-   * @brief Computes all the eigenvalues and the eigenvectors of a real
-   * generalized symmetric-definite eigenproblem, of the form:
-   * A*x=(lambda)*B*x, using the LAPACK function <code>dsygvd</code>. This
-   * version does <b>NOT</b> perform any checks on the input data and should
-   * be, therefore, faster.
-   *
-   * @warning The input matrices A and B are assumed to be symmetric and B is
-   * also positive definite.
-   *
-   * @param A The A input matrix (size NxN) of the problem
-   *
-   * @param B The B input matrix (size NxN) of the problem
-   *
-   * @param V The V matrix of eigenvectors (size NxN) stored in columns
-   *
-   * @param D The vector of eigenvalues (size N) (in ascending order). Note
-   * that this is a 1D array rather than a 2D diagonal matrix!
-   */
-  void eigSym_(const blitz::Array<double,2>& A, const blitz::Array<double,2>& B,
-      blitz::Array<double,2>& V, blitz::Array<double,1>& D);
-
 }}
 
 #endif /* BOB_MATH_EIG_H */
diff --git a/bob/math/include/bob.math/pinv.h b/bob/math/include/bob.math/pinv.h
index 3d55a882c7c8d12bc2d14b13114b2458001a52c0..1a8230fe45a2eb13df85ed208bd66aec931a591a 100644
--- a/bob/math/include/bob.math/pinv.h
+++ b/bob/math/include/bob.math/pinv.h
@@ -27,20 +27,7 @@ namespace bob { namespace math {
    */
   void pinv(const blitz::Array<double,2>& A, blitz::Array<double,2>& B,
       const double rcond=1e-15);
-  /**
-   * @brief Function which computes the pseudo-inverse using the SVD method.
-   * @warning The output blitz::array B should have the correct
-   *   size, with zero base index. Checks are NOT performed.
-   * @param A The A matrix to decompose (size MxN)
-   * @param B The pseudo-inverse of the matrix A (size NxM)
-   * @param rcond Cutoff for small singular values. Singular values smaller
-   *   (in modulus) than rcond * largest_singular_value (again, in modulus)
-   *   are set to zero.
-   */
-  void pinv_(const blitz::Array<double,2>& A, blitz::Array<double,2>& B,
-      const double rcond=1e-15);
 
 }}
 
 #endif /* BOB_MATH_PINV_H */
-
diff --git a/bob/math/include/bob.math/stats.h b/bob/math/include/bob.math/stats.h
index bc6ab6b224506631d8b74bd4549654be4fd00c6a..6d15c1d8d4f33a1192f8cc06b46a1b720b309b79 100644
--- a/bob/math/include/bob.math/stats.h
+++ b/bob/math/include/bob.math/stats.h
@@ -23,15 +23,20 @@ namespace bob { namespace math {
    * organized row-wise (each sample is a row, each feature is a column).
    * Outputs the sample mean M and the scatter matrix S.
    *
-   * @warning No checks are performed on the array sizes and is recommended
-   * only in scenarios where you have previously checked conformity and is
-   * focused only on speed.
+   * The input and output data have their sizes checked and this method will
+   * raise an appropriate exception if that is not cased.
    *
    * This version of the method also returns the sample mean of the array.
    */
   template<typename T>
-    void scatter_(const blitz::Array<T,2>& A, blitz::Array<T,2>& S,
+    void scatter(const blitz::Array<T,2>& A, blitz::Array<T,2>& S,
         blitz::Array<T,1>& M) {
+
+      // Check output
+      bob::core::array::assertSameDimensionLength(A.extent(1), M.extent(0));
+      bob::core::array::assertSameDimensionLength(A.extent(1), S.extent(0));
+      bob::core::array::assertSameDimensionLength(A.extent(1), S.extent(1));
+
       blitz::firstIndex i;
       blitz::secondIndex j;
       blitz::Range a = blitz::Range::all();
@@ -46,52 +51,13 @@ namespace bob { namespace math {
       }
     }
 
-  /**
-   * @brief Computes the scatter matrix of a 2D array considering data is
-   * organized row-wise (each sample is a row, each feature is a column).
-   * Outputs the sample mean M and the scatter matrix S.
-   *
-   * The input and output data have their sizes checked and this method will
-   * raise an appropriate exception if that is not cased. If you know that
-   * the input and output matrices conform, use the scatter_() variant.
-   *
-   * This version of the method also returns the sample mean of the array.
-   */
-  template<typename T>
-    void scatter(const blitz::Array<T,2>& A, blitz::Array<T,2>& S,
-        blitz::Array<T,1>& M) {
-
-      // Check output
-      bob::core::array::assertSameDimensionLength(A.extent(1), M.extent(0));
-      bob::core::array::assertSameDimensionLength(A.extent(1), S.extent(0));
-      bob::core::array::assertSameDimensionLength(A.extent(1), S.extent(1));
-
-      scatter_<T>(A, S, M);
-    }
-
-  /**
-   * @brief Computes the scatter matrix of a 2D array considering data is
-   * organized row-wise (each sample is a row, each feature is a column).
-   * Outputs the sample scatter matrix S.
-   *
-   * @warning No checks are performed on the array sizes and is recommended
-   * only in scenarios where you have previously checked conformity and is
-   * focused only on speed.
-   */
-  template<typename T>
-    void scatter_(const blitz::Array<T,2>& A, blitz::Array<T,2>& S) {
-      blitz::Array<T,1> M;
-      scatter_<T>(A, S, M);
-    }
-
   /**
    * @brief Computes the scatter matrix of a 2D array considering data is
    * organized row-wise (each sample is a row, each feature is a column).
    * Outputs the sample scatter matrix S.
    *
    * The input and output data have their sizes checked and this method will
-   * raise an appropriate exception if that is not cased. If you know that
-   * the input and output matrices conform, use the scatter_() variant.
+   * raise an appropriate exception if that is not cased.
    */
   template<typename T>
     void scatter(const blitz::Array<T,2>& A, blitz::Array<T,2>& S) {
@@ -149,16 +115,20 @@ namespace bob { namespace math {
    *
    * This method was designed based on the previous design at
    * torch3Vision 2.1, by SM.
-   *
-   * @warning No checks are performed on the array sizes and is recommended
-   * only in scenarios where you have previously checked conformity and is
-   * focused only on speed.
    */
   template <typename T>
-    void scatters_(const std::vector<blitz::Array<T,2> >& data,
+    void scatters(const std::vector<blitz::Array<T,2> >& data,
         blitz::Array<T,2>& Sw, blitz::Array<T,2>& Sb,
         blitz::Array<T,1>& m)
     {
+      // Check output
+      for (size_t i=0; i<data.size(); ++i)
+        bob::core::array::assertSameDimensionLength(data[i].extent(1), m.extent(0));
+      bob::core::array::assertSameDimensionLength(m.extent(0), Sw.extent(0));
+      bob::core::array::assertSameDimensionLength(m.extent(0), Sw.extent(1));
+      bob::core::array::assertSameDimensionLength(m.extent(0), Sb.extent(0));
+      bob::core::array::assertSameDimensionLength(m.extent(0), Sb.extent(1));
+
       // checks for data shape should have been done before...
       const int n_features = data[0].extent(1);
 
@@ -193,74 +163,6 @@ namespace bob { namespace math {
       }
     }
 
-  /**
-   * @brief Calculates the within and between class scatter matrices Sw and
-   * Sb. Returns those matrices and the overall means vector (m).
-   *
-   * Strategy implemented:
-   * 1. Evaluate the overall mean (m), class means (m_k) and the total class
-   *    counts (N).
-   * 2. Evaluate Sw and Sb using normal loops.
-   *
-   * Note that Sw and Sb, in this implementation, will be normalized by N-1
-   * (number of samples) and K (number of classes). This procedure makes
-   * the eigen values scaled by (N-1)/K, effectively increasing their
-   * values. The main motivation for this normalization are numerical
-   * precision concerns with the increasing number of samples causing a
-   * rather large Sw matrix. A normalization strategy mitigates this
-   * problem. The eigen vectors will see no effect on this normalization as
-   * they are normalized in the euclidean sense (||a|| = 1) so that does
-   * not change those.
-   *
-   * This method was designed based on the previous design at
-   * torch3Vision 2.1, by SM.
-   */
-  template <typename T>
-    void scatters(const std::vector<blitz::Array<T,2> >& data,
-        blitz::Array<T,2>& Sw, blitz::Array<T,2>& Sb,
-        blitz::Array<T,1>& m)
-    {
-      // Check output
-      for (size_t i=0; i<data.size(); ++i)
-        bob::core::array::assertSameDimensionLength(data[i].extent(1), m.extent(0));
-      bob::core::array::assertSameDimensionLength(m.extent(0), Sw.extent(0));
-      bob::core::array::assertSameDimensionLength(m.extent(0), Sw.extent(1));
-      bob::core::array::assertSameDimensionLength(m.extent(0), Sb.extent(0));
-      bob::core::array::assertSameDimensionLength(m.extent(0), Sb.extent(1));
-
-      scatters_<T>(data, Sw, Sb, m);
-    }
-
-  /**
-   * @brief Calculates the within and between class scatter matrices Sw and
-   * Sb. Returns those matrices.
-   *
-   * Strategy implemented:
-   * 1. Evaluate the overall mean (m), class means (m_k) and the total class
-   *    counts (N).
-   * 2. Evaluate Sw and Sb using normal loops.
-   *
-   * Note that Sw and Sb, in this implementation, will be normalized by N-1
-   * (number of samples) and K (number of classes). This procedure makes
-   * the eigen values scaled by (N-1)/K, effectively increasing their
-   * values. The main motivation for this normalization are numerical
-   * precision concerns with the increasing number of samples causing a
-   * rather large Sw matrix. A normalization strategy mitigates this
-   * problem. The eigen vectors will see no effect on this normalization as
-   * they are normalized in the euclidean sense (||a|| = 1) so that does
-   * not change those.
-   *
-   * This method was designed based on the previous design at
-   * torch3Vision 2.1, by SM.
-   */
-  template<typename T>
-    void scatters_(const std::vector<blitz::Array<T,2> >& data,
-        blitz::Array<T,2>& Sw, blitz::Array<T,2>& Sb)
-    {
-      blitz::Array<T,1> M(data[0].extent(1));
-      scatters_<T>(data, Sw, Sb, M);
-    }
-
   /**
    * @brief Calculates the within and between class scatter matrices Sw and
    * Sb. Returns those matrices.
diff --git a/bob/math/include/bob.math/svd.h b/bob/math/include/bob.math/svd.h
index cf90156bbc6435e6fca4aa64001bd50407008a76..dd1ebb8f6b5ff8884a1b832ced9d5e965c29130e 100644
--- a/bob/math/include/bob.math/svd.h
+++ b/bob/math/include/bob.math/svd.h
@@ -22,7 +22,7 @@ namespace bob { namespace math {
 /**
  * @brief Function which performs a 'full' Singular Value Decomposition
  *   using the divide and conquer routine dgesdd of LAPACK.
- * @warning The output blitz::array U, sigma and Vt should have the correct 
+ * @warning The output blitz::array U, sigma and Vt should have the correct
  *   size, with zero base index. Checks are performed.
  * @param A The A matrix to decompose (size MxN)
  * @param U The U matrix of left singular vectors (size MxM)
@@ -33,28 +33,13 @@ namespace bob { namespace math {
  */
 void svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
   blitz::Array<double,1>& sigma, blitz::Array<double,2>& Vt, bool safe=false);
-/**
- * @brief Function which performs a 'full' Singular Value Decomposition
- *   using the divide and conquer routine dgesdd of LAPACK.
- * @warning The output blitz::array U, sigma and Vt should have the correct 
- *   size, with zero base index. Checks are NOT performed.
- * @param A The A matrix to decompose (size MxN)
- * @param U The U matrix of left singular vectors (size MxM)
- * @param sigma The vector of singular values (size min(M,N))
- *    Please note that this is a 1D array rather than a 2D diagonal matrix!
- * @param Vt The V^T matrix of right singular vectors (size NxN)
- * @param safe If enabled, use LAPACK dgesvd instead of dgesdd
- */
-void svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
-  blitz::Array<double,1>& sigma, blitz::Array<double,2>& Vt, bool safe=false);
-
 
 /**
  * @brief Function which performs a 'partial' Singular Value Decomposition
- *   using the 'simple' driver routine dgesvd of LAPACK. It only returns 
- *   the first min(M,N) columns of U, and is somehow similar to the 
+ *   using the 'simple' driver routine dgesvd of LAPACK. It only returns
+ *   the first min(M,N) columns of U, and is somehow similar to the
  *   'economical' SVD variant of matlab (except that it does not return V).
- * @warning The output blitz::array U and sigma should have the correct 
+ * @warning The output blitz::array U and sigma should have the correct
  *   size, with zero base index. Checks are performed.
  * @param A The A matrix to decompose (size MxN)
  * @param U The U matrix of left singular vectors (size M x min(M,N))
@@ -64,28 +49,12 @@ void svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
  */
 void svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
   blitz::Array<double,1>& sigma, bool safe=false);
-/**
- * @brief Function which performs a 'partial' Singular Value Decomposition
- *   using the 'simple' driver routine dgesvd of LAPACK. It only returns 
- *   the first min(M,N) columns of U, and is somehow similar to the 
- *   'economical' SVD variant of matlab (except that it does not return V).
- * @warning The output blitz::array U and sigma should have the correct 
- *   size, with zero base index. Checks are NOT performed.
- * @param A The A matrix to decompose (size MxN)
- * @param U The U matrix of left singular vectors (size M x min(M,N))
- * @param sigma The vector of singular values (size min(M,N))
- *    Please note that this is a 1D array rather than a 2D diagonal matrix!
- * @param safe If enabled, use LAPACK dgesvd instead of dgesdd
- */
-void svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
-  blitz::Array<double,1>& sigma, bool safe=false);
-
 
 /**
  * @brief Function which performs a 'partial' Singular Value Decomposition
- *   using the 'simple' driver routine dgesvd of LAPACK. It only returns 
+ *   using the 'simple' driver routine dgesvd of LAPACK. It only returns
  *   the singular values.
- * @warning The output blitz::array sigma should have the correct 
+ * @warning The output blitz::array sigma should have the correct
  *   size, with zero base index. Checks are performed.
  * @param A The A matrix to decompose (size MxN)
  * @param sigma The vector of singular values (size min(M,N))
@@ -94,19 +63,6 @@ void svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
  */
 void svd(const blitz::Array<double,2>& A, blitz::Array<double,1>& sigma,
   bool safe=false);
-/**
- * @brief Function which performs a 'partial' Singular Value Decomposition
- *   using the 'simple' driver routine dgesvd of LAPACK. It only returns 
- *   the singular values.
- * @warning The output blitz::array sigma should have the correct 
- *   size, with zero base index. Checks are NOT performed.
- * @param A The A matrix to decompose (size MxN)
- * @param sigma The vector of singular values (size min(M,N))
- *    Please note that this is a 1D array rather than a 2D diagonal matrix!
- * @param safe If enabled, use LAPACK dgesvd instead of dgesdd
- */
-void svd_(const blitz::Array<double,2>& A, blitz::Array<double,1>& sigma,
-  bool safe=false);
 
 /**
  * @}
diff --git a/bob/math/main.cpp b/bob/math/main.cpp
index 4ec08bac052f7371cef8c3af99a090d9107810d4..8fb8d17c87a665f08a830b5048860c06f28eb1f7 100644
--- a/bob/math/main.cpp
+++ b/bob/math/main.cpp
@@ -169,23 +169,6 @@ static bob::extension::FunctionDoc s_scatter = bob::extension::FunctionDoc(
   .add_return("m", "array_like (float, 1D)", "The mean matrix, with with the row means of ``a``")
 ;
 
-static bob::extension::FunctionDoc s_scatter_nocheck = bob::extension::FunctionDoc(
-  "scatter_",
-  "Computes scatter matrix of a 2D array.",
-  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :py:func:`scatter`."
-  "Use it when you are sure your input matrices sizes match.\n\n"
-  "Computes the scatter matrix of a 2D array *considering data is organized row-wise* (each sample is a row, each feature is a column). "
-  "The resulting array ``s`` is squared with extents equal to the number of columns in ``a``. "
-  "The resulting array ``m`` is a 1D array with the row means of ``a``. "
-  "This function supports many calling modes, but you should provide, at least, the input data matrix ``a``. "
-  "All non-provided arguments will be allocated internally and returned."
-  )
-  .add_prototype("a, s, m")
-  .add_parameter("a", "array_like (float, 2D)", "The sample matrix, *considering data is organized row-wise* (each sample is a row, each feature is a column)")
-  .add_parameter("s", "array_like (float, 2D)", "The scatter matrix, squared with extents equal to the number of columns in ``a``")
-  .add_parameter("m", "array_like (float,1D)", "The mean matrix, with with the row means of ``a``")
-;
-
 
 static bob::extension::FunctionDoc s_scatters = bob::extension::FunctionDoc(
   "scatters",
@@ -218,24 +201,6 @@ static bob::extension::FunctionDoc s_scatters = bob::extension::FunctionDoc(
   .add_return("m", "array_like (float, 1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
 ;
 
-static bob::extension::FunctionDoc s_scatters_nocheck = bob::extension::FunctionDoc(
-  "scatters_",
-  "Computes :math:`S_w` and :math:`S_b` scatter matrices of a set of 2D arrays.",
-  ".. warning:: This variant does not perform any checks on the input matrices and is faster then :py:func:`scatters`. "
-  "Use it when you are sure your input matrices sizes match.\n\n"
-  "For a detailed description of the function, please see :func:`scatters`."
-  )
-  .add_prototype("data, sw, sb, m")
-  .add_prototype("data, sw, sb")
-  .add_parameter("data", "[array_like (float, 2D)]", "The list of sample matrices. "
-      "In each sample matrix the data is organized row-wise (each sample is a row, each feature is a column). "
-      "Each matrix stores the data of a particular class. "
-      "**Every matrix in ``data`` must have exactly the same number of columns.**")
-  .add_parameter("sw", "array_like (float, 2D)", "The within-class scatter matrix :math:`S_w`, squared with extents equal to the number of columns in ``data``")
-  .add_parameter("sb", "array_like (float, 2D)", "The between-class scatter matrix :math:`S_b`, squared with extents equal to the number of columns in ``data``")
-  .add_parameter("m", "array_like (float,1D)", "The mean matrix, representing the ensemble mean with no prior (i.e., biased towards classes with more samples)")
-;
-
 
 static bob::extension::FunctionDoc s_gsvd = bob::extension::FunctionDoc(
   "gsvd",
@@ -346,24 +311,12 @@ static PyMethodDef module_methods[] = {
       METH_VARARGS|METH_KEYWORDS,
       s_scatter.doc()
     },
-    {
-      s_scatter_nocheck.name(),
-      (PyCFunction)py_scatter_nocheck,
-      METH_VARARGS|METH_KEYWORDS,
-      s_scatter_nocheck.doc()
-    },
     {
       s_scatters.name(),
       (PyCFunction)py_scatters,
       METH_VARARGS|METH_KEYWORDS,
       s_scatters.doc()
     },
-    {
-      s_scatters_nocheck.name(),
-      (PyCFunction)py_scatters_nocheck,
-      METH_VARARGS|METH_KEYWORDS,
-      s_scatters_nocheck.doc()
-    },
     {
       s_gsvd.name(),
       (PyCFunction)py_gsvd,
diff --git a/bob/math/scatter.cpp b/bob/math/scatter.cpp
index 44f9154e8cba3bdba633ef35dae5f4ba201b8dcf..d848d9ea27f648db0e2ec08b7493d5d3b5179f67 100644
--- a/bob/math/scatter.cpp
+++ b/bob/math/scatter.cpp
@@ -111,79 +111,6 @@ PyObject* py_scatter (PyObject*, PyObject* args, PyObject* kwds) {
 
 }
 
-PyObject* py_scatter_nocheck (PyObject*, PyObject* args, PyObject* kwds) {
-
-  /* Parses input arguments in a single shot */
-  static const char* const_kwlist[] = { "a", "s", "m", 0 /* Sentinel */ };
-  static char** kwlist = const_cast<char**>(const_kwlist);
-
-  PyBlitzArrayObject* a = 0;
-  PyBlitzArrayObject* s = 0;
-  PyBlitzArrayObject* m = 0;
-
-  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&O&",
-        kwlist,
-        &PyBlitzArray_Converter, &a,
-        &PyBlitzArray_OutputConverter, &s,
-        &PyBlitzArray_OutputConverter, &m
-        )) return 0;
-
-  //protects acquired resources through this scope
-  auto a_ = make_safe(a);
-  auto s_ = make_safe(s);
-  auto m_ = make_safe(m);
-
-  // basic checks
-  if (a->ndim != 2 || (a->type_num != NPY_FLOAT32 && a->type_num != NPY_FLOAT64)) {
-    PyErr_SetString(PyExc_TypeError, "input data matrix `a' should be either a 32 or 64-bit float 2D array");
-    return 0;
-  }
-
-  if (s->ndim != 2 || (s->type_num != a->type_num)) {
-    PyErr_SetString(PyExc_TypeError, "output data matrix `s' should be either a 32 or 64-bit float 2D array, matching the data type of `a'");
-    return 0;
-  }
-
-  if (m->ndim != 1 || (m->type_num != a->type_num)) {
-    PyErr_SetString(PyExc_TypeError, "output data vector `m' should be either a 32 or 64-bit float 1D array, matching the data type of `a'");
-    return 0;
-  }
-
-  try {
-    switch (a->type_num) {
-      case NPY_FLOAT32:
-        bob::math::scatter(
-            *PyBlitzArrayCxx_AsBlitz<float,2>(a),
-            *PyBlitzArrayCxx_AsBlitz<float,2>(s),
-            *PyBlitzArrayCxx_AsBlitz<float,1>(m)
-            );
-        break;
-
-      case NPY_FLOAT64:
-        bob::math::scatter(
-            *PyBlitzArrayCxx_AsBlitz<double,2>(a),
-            *PyBlitzArrayCxx_AsBlitz<double,2>(s),
-            *PyBlitzArrayCxx_AsBlitz<double,1>(m)
-            );
-        break;
-
-      default:
-        PyErr_Format(PyExc_TypeError, "(no-check) scatter calculation currently not implemented for type '%s'", PyBlitzArray_TypenumAsString(a->type_num));
-        return 0;
-    }
-  }
-  catch (std::exception& e) {
-    PyErr_SetString(PyExc_RuntimeError, e.what());
-    return 0;
-  }
-  catch (...) {
-    PyErr_SetString(PyExc_RuntimeError, "(no-check) scatter calculation failed: unknown exception caught");
-    return 0;
-  }
-
-  Py_RETURN_NONE;
-
-}
 
 /**
  * Converts the input iterable d into a tuple of PyBlitzArrayObject's. Checks
@@ -373,94 +300,3 @@ PyObject* py_scatters (PyObject*, PyObject* args, PyObject* kwds) {
   return retval;
 
 }
-
-PyObject* py_scatters_nocheck (PyObject*, PyObject* args, PyObject* kwds) {
-
-  /* Parses input arguments in a single shot */
-  static const char* const_kwlist[] = { "data", "sw", "sb", "m", 0 /* Sentinel */ };
-  static char** kwlist = const_cast<char**>(const_kwlist);
-
-  PyObject* data = 0;
-  PyBlitzArrayObject* sw = 0;
-  PyBlitzArrayObject* sb = 0;
-  PyBlitzArrayObject* m = 0;
-
-  if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&O&O&O&", kwlist,
-        &BzTuple_Converter, &data,
-        &PyBlitzArray_OutputConverter, &sw,
-        &PyBlitzArray_OutputConverter, &sb,
-        &PyBlitzArray_OutputConverter, &m
-        )) return 0;
-
-  //protects acquired resources through this scope
-  auto data_ = make_safe(data);
-  auto sw_ = make_safe(sw);
-  auto sb_ = make_safe(sb);
-  auto m_ = make_safe(m);
-
-  PyBlitzArrayObject* first = (PyBlitzArrayObject*)PyTuple_GET_ITEM(data, 0);
-
-  if (sw->ndim != 2 || (sw->type_num != first->type_num)) {
-    PyErr_SetString(PyExc_TypeError, "output data matrix `sw' should be either a 32 or 64-bit float 2D array, matching the data type of `data'");
-    return 0;
-  }
-
-  if (sb->ndim != 2 || (sb->type_num != first->type_num)) {
-    PyErr_SetString(PyExc_TypeError, "output data matrix `sb' should be either a 32 or 64-bit float 2D array, matching the data type of `data'");
-    return 0;
-  }
-
-  if (m->ndim != 1 || (m->type_num != first->type_num)) {
-    PyErr_SetString(PyExc_TypeError, "output data vector `m' should be either a 32 or 64-bit float 1D array, matching the data type of `data'");
-    return 0;
-  }
-
-  try {
-    switch (first->type_num) {
-      case NPY_FLOAT32:
-        {
-          std::vector<blitz::Array<float,2>> cxxdata;
-          for (Py_ssize_t i=0; i<PyTuple_GET_SIZE(data); ++i) {
-            cxxdata.push_back(*PyBlitzArrayCxx_AsBlitz<float,2>
-                ((PyBlitzArrayObject*)PyTuple_GET_ITEM(data,i)));
-            bob::math::scatters_(cxxdata,
-                *PyBlitzArrayCxx_AsBlitz<float,2>(sw),
-                *PyBlitzArrayCxx_AsBlitz<float,2>(sb),
-                *PyBlitzArrayCxx_AsBlitz<float,1>(m)
-                );
-          }
-        }
-        break;
-
-      case NPY_FLOAT64:
-        {
-          std::vector<blitz::Array<double,2>> cxxdata;
-          for (Py_ssize_t i=0; i<PyTuple_GET_SIZE(data); ++i) {
-            cxxdata.push_back(*PyBlitzArrayCxx_AsBlitz<double,2>
-                ((PyBlitzArrayObject*)PyTuple_GET_ITEM(data,i)));
-            bob::math::scatters_(cxxdata,
-                *PyBlitzArrayCxx_AsBlitz<double,2>(sw),
-                *PyBlitzArrayCxx_AsBlitz<double,2>(sb),
-                *PyBlitzArrayCxx_AsBlitz<double,1>(m)
-                );
-          }
-        }
-        break;
-
-      default:
-        PyErr_Format(PyExc_TypeError, "(no-check) scatters calculation currently not implemented for type '%s'", PyBlitzArray_TypenumAsString(first->type_num));
-        return 0;
-    }
-  }
-  catch (std::exception& e) {
-    PyErr_SetString(PyExc_RuntimeError, e.what());
-    return 0;
-  }
-  catch (...) {
-    PyErr_SetString(PyExc_RuntimeError, "(no-check) scatters calculation failed: unknown exception caught");
-    return 0;
-  }
-
-  Py_RETURN_NONE;
-
-}
diff --git a/bob/math/scatter.h b/bob/math/scatter.h
index e45a389609a956622880397dc198bd8c4e621d3e..2b22c616d3e98d2beb84fe512476baad669b7422 100644
--- a/bob/math/scatter.h
+++ b/bob/math/scatter.h
@@ -1,6 +1,6 @@
 /**
  * @author Andre Anjos <andre.anjos@idiap.ch>
- * @date Thu  5 Dec 12:10:18 2013 
+ * @date Thu  5 Dec 12:10:18 2013
  *
  * @brief Declaration of components relevant for main.cpp
  */
@@ -8,6 +8,4 @@
 #include <Python.h>
 
 PyObject* py_scatter(PyObject*, PyObject* args, PyObject* kwds);
-PyObject* py_scatter_nocheck(PyObject*, PyObject* args, PyObject* kwds);
 PyObject* py_scatters(PyObject*, PyObject* args, PyObject* kwds);
-PyObject* py_scatters_nocheck(PyObject*, PyObject* args, PyObject* kwds);
diff --git a/bob/math/test_scatter.py b/bob/math/test_scatter.py
index 036acc4bbf066560ec9c2cd4d394136ce4f83305..6dd99f193e74f8758d025e2760decb6e12e83ef7 100644
--- a/bob/math/test_scatter.py
+++ b/bob/math/test_scatter.py
@@ -9,7 +9,7 @@
 """
 
 import os, sys
-from bob.math import scatter, scatter_, scatters, scatters_
+from bob.math import scatter, scatters
 import numpy
 import nose.tools
 
@@ -112,7 +112,7 @@ def test_fast_scatter():
   # This test demonstrates how to use the scatter matrix function of bob.
   S = numpy.ndarray((data.shape[1], data.shape[1]), dtype=float)
   M = numpy.ndarray((data.shape[1],), dtype=float)
-  scatter_(data, S, M)
+  scatter(data, S, M)
   S /= (data.shape[0]-1)
 
   # Do the same with numpy and compare. Note that with numpy we are computing
@@ -185,7 +185,7 @@ def test_fast_scatters():
   Sw = numpy.empty_like(Sw_)
   Sb = numpy.empty_like(Sb_)
   m = numpy.empty_like(m_)
-  scatters_(data, Sw, Sb, m)
+  scatters(data, Sw, Sb, m)
   assert numpy.allclose(Sw, Sw_)
   assert numpy.allclose(Sb, Sb_)
   assert numpy.allclose(m, m_)