diff --git a/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp b/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp
index 43957f1bed8bb60e49f3f8219764ff09618a69de..2d9e7617ba7f4cecc91e00d1c1504e5fa99d77c6 100644
--- a/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp
+++ b/include/eigenpy/solvers/LeastSquaresConjugateGradient.hpp
@@ -1,5 +1,5 @@
 /*
- * Copyright 2017, Justin Carpentier, LAAS-CNRS
+ * Copyright 2017-2018, Justin Carpentier, LAAS-CNRS
  *
  * This file is part of eigenpy.
  * eigenpy is free software: you can redistribute it and/or
@@ -22,67 +22,6 @@
 
 #include "eigenpy/solvers/IterativeSolverBase.hpp"
 
-namespace Eigen
-{
-  template <typename _Scalar>
-  class LeastSquareDiagonalPreconditionerFix
-  : public LeastSquareDiagonalPreconditioner<_Scalar>
-  {
-    typedef _Scalar Scalar;
-    typedef typename NumTraits<Scalar>::Real RealScalar;
-    typedef LeastSquareDiagonalPreconditioner<_Scalar> Base;
-    using DiagonalPreconditioner<_Scalar>::m_invdiag;
-  public:
-    
-    LeastSquareDiagonalPreconditionerFix() : Base() {}
-    
-    template<typename MatType>
-    explicit LeastSquareDiagonalPreconditionerFix(const MatType& mat) : Base()
-    {
-      compute(mat);
-    }
-    
-    template<typename MatType>
-    LeastSquareDiagonalPreconditionerFix& analyzePattern(const MatType& )
-    {
-      return *this;
-    }
-    
-    template<typename MatType>
-    LeastSquareDiagonalPreconditionerFix& factorize(const MatType& mat)
-    {
-      // Compute the inverse squared-norm of each column of mat
-      m_invdiag.resize(mat.cols());
-      if(MatType::IsRowMajor)
-      {
-        m_invdiag.setZero();
-        for(Index j=0; j<mat.outerSize(); ++j)
-        {
-          for(typename MatType::InnerIterator it(mat,j); it; ++it)
-            m_invdiag(it.index()) += numext::abs2(it.value());
-        }
-        for(Index j=0; j<mat.cols(); ++j)
-          if(numext::real(m_invdiag(j))>RealScalar(0))
-            m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j));
-      }
-      else
-      {
-        for(Index j=0; j<mat.outerSize(); ++j)
-        {
-          RealScalar sum = mat.col(j).squaredNorm();
-          if(sum>RealScalar(0))
-            m_invdiag(j) = RealScalar(1)/sum;
-          else
-            m_invdiag(j) = RealScalar(1);
-        }
-      }
-      Base::m_isInitialized = true;
-      return *this;
-    }
-    
-  };
-}
-
 namespace eigenpy
 {
   
diff --git a/src/solvers/solvers.cpp b/src/solvers/solvers.cpp
index 8a30217ddd9fb2854049f9f501afb0e356272cef..1317b2dc2ce40d981b9d60066d180582e71ae2e4 100644
--- a/src/solvers/solvers.cpp
+++ b/src/solvers/solvers.cpp
@@ -32,7 +32,7 @@ namespace eigenpy
     using namespace Eigen;
     ConjugateGradientVisitor< ConjugateGradient<MatrixXd,Lower|Upper> >::expose();
 #if EIGEN_VERSION_AT_LEAST(3,3,5)
-    LeastSquaresConjugateGradientVisitor< LeastSquaresConjugateGradient<MatrixXd, LeastSquareDiagonalPreconditionerFix<MatrixXd::Scalar> > >::expose();
+    LeastSquaresConjugateGradientVisitor< LeastSquaresConjugateGradient<MatrixXd, LeastSquareDiagonalPreconditioner<MatrixXd::Scalar> > >::expose();
 #endif
     
     // Conjugate gradient with limited BFGS preconditioner