diff --git a/include/crocoddyl/core/action-base.hpp b/include/crocoddyl/core/action-base.hpp
index 7ce25312a279eb568ba3e38931a205a7ab2d14cd..df24d2599227e1f3558a8ff2ec2516d40fc9ef41 100644
--- a/include/crocoddyl/core/action-base.hpp
+++ b/include/crocoddyl/core/action-base.hpp
@@ -90,7 +90,6 @@ struct ActionDataAbstract {
   template <typename Model>
   explicit ActionDataAbstract(Model* const model)
       : cost(0.),
-        costResidual(model->get_nr()),
         xnext(model->get_state().get_nx()),
         r(model->get_nr()),
         Fx(model->get_state().get_ndx(), model->get_state().get_ndx()),
@@ -123,7 +122,6 @@ struct ActionDataAbstract {
   const Eigen::MatrixXd& get_Fu() const { return Fu; }
 
   double cost;
-  Eigen::VectorXd costResidual;
   Eigen::VectorXd xnext;
   Eigen::VectorXd r;
   Eigen::MatrixXd Fx;
diff --git a/include/crocoddyl/core/numdiff/action.hpp b/include/crocoddyl/core/numdiff/action.hpp
index 2253a9811e49298cb6336ec4a48e0267b45cddc2..022a4a29a125802ff25692643018576e5d27a3fc 100644
--- a/include/crocoddyl/core/numdiff/action.hpp
+++ b/include/crocoddyl/core/numdiff/action.hpp
@@ -21,18 +21,18 @@ namespace crocoddyl {
  * It computes the same quantity as a normal model would do but using numerical
  * differentiation.
  * The subtility is in the computation of the Hessian of the cost. Let us
- * concider that the ActionModel owns a cost residual. This means that the cost
+ * consider that the ActionModel owns a cost residual. This means that the cost
  * is the square of a residual \f$ l(x,u) = .5 r(x,u)**2 \f$, with
- * \f$ r(x,u) \f$ being a vector. Therefore the derivatives of the cost
- * \f$ l \f$ can be expressed in function of the derivatives of the residuals
- * (jacobians), denoted by \f$ R_x \f$ and \f$ R_u \f$. Which would be:
+ * \f$ r(x,u) \f$ being the reisual vector. Therefore the derivatives of the
+ * cost \f$ l \f$ can be expressed in function of the derivatives of the
+ * residuals (Jacobians), denoted by \f$ R_x \f$ and \f$ R_u \f$. Which would be:
  * \f{eqnarray*}{
  *     L_x    &=& R_x^T r \\
  *     L_u    &=& R_u^T r \\
  *     L_{xx} &=& R_x^T R_x + R_{xx} r
  * \f}
- * with \f$ R_{xx} \f$ the derivatives of the jacobian (i.e. not a matrix, but a
- * dim-3 tensor). The Gauss approximation boils down to neglecting this terms.
+ * with \f$ R_{xx} \f$ the derivatives of the Jacobian (i.e. not a matrix, but a
+ * dim-3 tensor). The Gauss approximation consists in neglecting these.
  * So \f$ L_{xx} \sim R_x^T R_x \f$. Similarly for \f$ L_{xu} \sim R_x^T R_u \f$
  * and \f$ L_{uu} \sim R_u^T R_u \f$. The above set of equations becomes:
  * \f{eqnarray*}{
@@ -54,26 +54,45 @@ class ActionModelNumDiff : public ActionModelAbstract {
    * @param with_gauss_approx defines if we use the Gauss approximation of the
    * cost hessian or not.
    */
-  explicit ActionModelNumDiff(ActionModelAbstract& model, bool with_gauss_approx = false);
+  explicit ActionModelNumDiff(ActionModelAbstract& model);
+
+  /**
+   * @brief Destroy the ActionModelNumDiff object
+   */
   ~ActionModelNumDiff();
 
+  /**
+   * @brief @copydoc ActionDataAbstract::calc()
+   */
   void calc(const boost::shared_ptr<ActionDataAbstract>& data, const Eigen::Ref<const Eigen::VectorXd>& x,
             const Eigen::Ref<const Eigen::VectorXd>& u);
+  
   /**
-   * @brief calcDiff computes the
-   *
-   * @param data
-   * @param x
-   * @param u
-   * @param recalc
+   * @brief @copydoc ActionDataAbstract::calcDiff()
    */
   void calcDiff(const boost::shared_ptr<ActionDataAbstract>& data, const Eigen::Ref<const Eigen::VectorXd>& x,
                 const Eigen::Ref<const Eigen::VectorXd>& u, const bool& recalc = true);
+
+  /**
+   * @brief Create a Data object from the given model.
+   * 
+   * @return boost::shared_ptr<ActionDataAbstract> 
+   */
   boost::shared_ptr<ActionDataAbstract> createData();
 
+  /**
+   * @brief Get the model_ object
+   * 
+   * @return ActionModelAbstract& 
+   */
   ActionModelAbstract& get_model() const;
+
+  /**
+   * @brief Get the disturbance_ object
+   * 
+   * @return const double& 
+   */
   const double& get_disturbance() const;
-  bool get_with_gauss_approx();
 
  private:
   /**
@@ -89,8 +108,15 @@ class ActionModelNumDiff : public ActionModelAbstract {
    */
   void assertStableStateFD(const Eigen::Ref<const Eigen::VectorXd>& x);
 
+  /**
+   * @brief This is the model to compute the finite differenciation from
+   */
   ActionModelAbstract& model_;
-  bool with_gauss_approx_;
+
+  /**
+   * @brief This is the numerical disturbance value used during the numerical
+   * differenciations
+   */
   double disturbance_;
 };
 
@@ -127,8 +153,8 @@ struct ActionDataNumDiff : public ActionDataAbstract {
     }
   }
 
-  Eigen::MatrixXd Rx;  //!< Cost residual jacobian: d r / dx
-  Eigen::MatrixXd Ru;  //!< Cost residual jacobian: d r / du
+  Eigen::MatrixXd Rx;  //!< Cost residual jacobian: \f$ \frac{d r(x,u)}{dx} \f$
+  Eigen::MatrixXd Ru;  //!< Cost residual jacobian: \f$ \frac{d r(x,u)}{du} \f$
   Eigen::VectorXd dx;  //!< State disturbance
   Eigen::VectorXd du;  //!< Control disturbance
   Eigen::VectorXd xp;  //!< The integrated state from the disturbance on one DoF "\f$ \int x dx_i \f$"
diff --git a/src/core/numdiff/action.cpp b/src/core/numdiff/action.cpp
index de077f566d39f902b16f5c1623ee8d6c6877e586..b53aad49017ff274b8e0d035b6d23ce19b907258 100644
--- a/src/core/numdiff/action.cpp
+++ b/src/core/numdiff/action.cpp
@@ -10,11 +10,9 @@
 
 namespace crocoddyl {
 
-ActionModelNumDiff::ActionModelNumDiff(ActionModelAbstract& model, bool with_gauss_approx)
+ActionModelNumDiff::ActionModelNumDiff(ActionModelAbstract& model)
     : ActionModelAbstract(model.get_state(), model.get_nu(), model.get_nr()), model_(model) {
-  with_gauss_approx_ = with_gauss_approx;
   disturbance_ = std::sqrt(2.0 * std::numeric_limits<double>::epsilon());
-  assert((!with_gauss_approx_ || nr_ > 1) && "No Gauss approximation possible with nr = 1");
 }
 
 ActionModelNumDiff::~ActionModelNumDiff() {}
@@ -58,7 +56,7 @@ void ActionModelNumDiff::calcDiff(const boost::shared_ptr<ActionDataAbstract>& d
     model_.get_state().diff(xn0, xn, data_nd->Fx.col(ix));
 
     data->Lx(ix) = (c - c0) / disturbance_;
-    if (with_gauss_approx_) {
+    if (model_.get_nr() > 0) {
       data_nd->Rx.col(ix) = (data_nd->data_x[ix]->r - data_nd->data_0->r) / disturbance_;
     }
     data_nd->dx(ix) = 0.0;
@@ -76,14 +74,14 @@ void ActionModelNumDiff::calcDiff(const boost::shared_ptr<ActionDataAbstract>& d
     model_.get_state().diff(xn0, xn, data_nd->Fu.col(iu));
 
     data->Lu(iu) = (c - c0) / disturbance_;
-    if (with_gauss_approx_) {
+    if (model_.get_nr() > 0) {
       data_nd->Ru.col(iu) = (data_nd->data_u[iu]->r - data_nd->data_0->r) / disturbance_;
     }
     data_nd->du(iu) = 0.0;
   }
   data->Fu /= disturbance_;
 
-  if (with_gauss_approx_) {
+  if (model_.get_nr() > 0) {
     data->Lxx = data_nd->Rx.transpose() * data_nd->Rx;
     data->Lxu = data_nd->Rx.transpose() * data_nd->Ru;
     data->Luu = data_nd->Ru.transpose() * data_nd->Ru;
@@ -98,7 +96,7 @@ ActionModelAbstract& ActionModelNumDiff::get_model() const { return model_; }
 
 const double& ActionModelNumDiff::get_disturbance() const { return disturbance_; }
 
-bool ActionModelNumDiff::get_with_gauss_approx() { return with_gauss_approx_; }
+bool ActionModelNumDiff::get_with_gauss_approx() { return model_.get_nr() > 0; }
 
 void ActionModelNumDiff::assertStableStateFD(const Eigen::Ref<const Eigen::VectorXd>& /** x */) {
   // do nothing in the general case