diff --git a/src/mlpack/methods/ann/layer/CMakeLists.txt b/src/mlpack/methods/ann/layer/CMakeLists.txt
index 5c8ee2c3b0..1229781cbe 100644
--- a/src/mlpack/methods/ann/layer/CMakeLists.txt
+++ b/src/mlpack/methods/ann/layer/CMakeLists.txt
@@ -22,6 +22,8 @@ set(SOURCES
   dropout_impl.hpp
   elu.hpp
   elu_impl.hpp
+  fast_lstm.hpp
+  fast_lstm_impl.hpp
   glimpse.hpp
   glimpse_impl.hpp
   gru.hpp
diff --git a/src/mlpack/methods/ann/layer/fast_lstm.hpp b/src/mlpack/methods/ann/layer/fast_lstm.hpp
new file mode 100644
index 0000000000..e853bb16d2
--- /dev/null
+++ b/src/mlpack/methods/ann/layer/fast_lstm.hpp
@@ -0,0 +1,309 @@
+/**
+ * @file fast_lstm.hpp
+ * @author Marcus Edel
+ *
+ * Definition of the Fast LSTM class, which implements a Fast LSTM network
+ * layer.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef MLPACK_METHODS_ANN_LAYER_FAST_LSTM_HPP
+#define MLPACK_METHODS_ANN_LAYER_FAST_LSTM_HPP
+
+#include <mlpack/prereqs.hpp>
+#include <limits>
+
+namespace mlpack {
+namespace ann /** Artificial Neural Network. */ {
+
+/**
+ * An implementation of a faster version of the Fast LSTM network layer.
+ * Basically by combining the calculation of the input, forget, output gates
+ * and hidden state in a single step. The standard formula changes as follows:
+ *
+ * @f{eqnarray}{
+ * i &=& sigmoid(W \cdot x + W \cdot h + b) \\
+ * f &=& sigmoid(W  \cdot x + W \cdot h + b) \\
+ * z &=& tanh(W \cdot x + W \cdot h + b) \\
+ * c &=& f \cdot c + i \cdot z \\
+ * o &=& sigmoid(W \cdot x + W \cdot h + b) \\
+ * h &=& o \cdot tanh(c)
+ * @f}
+ *
+ * Note that FastLSTM network layer does not use peephole connections between
+ * the cell and gates.
+ *
+ * For more information, see the following.
+ *
+ * @code
+ * @article{Hochreiter1997,
+ *   author  = {Hochreiter, Sepp and Schmidhuber, J\"{u}rgen},
+ *   title   = {Long Short-term Memory},
+ *   journal = {Neural Comput.},
+ *   year    = {1997}
+ * }
+ * @endcode
+ *
+ * \see LSTM for a standard implementation of the LSTM layer.
+ *
+ * @tparam InputDataType Type of the input data (arma::colvec, arma::mat,
+ *         arma::sp_mat or arma::cube).
+ * @tparam OutputDataType Type of the output data (arma::colvec, arma::mat,
+ *         arma::sp_mat or arma::cube).
+ */
+template <
+    typename InputDataType = arma::mat,
+    typename OutputDataType = arma::mat
+>
+class FastLSTM
+{
+ public:
+  // Convenience typedefs.
+  typedef typename InputDataType::elem_type InputElemType;
+  typedef typename OutputDataType::elem_type ElemType;
+
+  //! Create the Fast LSTM object.
+  FastLSTM();
+
+  /**
+   * Create the Fast LSTM layer object using the specified parameters.
+   *
+   * @param inSize The number of input units.
+   * @param outSize The number of output units.
+   * @param rho Maximum number of steps to backpropagate through time (BPTT).
+   */
+  FastLSTM(const size_t inSize,
+           const size_t outSize,
+           const size_t rho = std::numeric_limits<size_t>::max());
+
+  /**
+   * Ordinary feed forward pass of a neural network, evaluating the function
+   * f(x) by propagating the activity forward through f.
+   *
+   * @param input Input data used for evaluating the specified function.
+   * @param output Resulting output activation.
+   */
+  template<typename InputType, typename OutputType>
+  void Forward(InputType&& input, OutputType&& output);
+
+  /**
+   * Ordinary feed backward pass of a neural network, calculating the function
+   * f(x) by propagating x backwards trough f. Using the results from the feed
+   * forward pass.
+   *
+   * @param input The propagated input activation.
+   * @param gy The backpropagated error.
+   * @param g The calculated gradient.
+   */
+  template<typename InputType, typename ErrorType, typename GradientType>
+  void Backward(const InputType&& input,
+                ErrorType&& gy,
+                GradientType&& g);
+
+  /*
+   * Reset the layer parameter.
+   */
+  void Reset();
+
+  /*
+   * Resets the cell to accept a new input. This breaks the BPTT chain starts a
+   * new one.
+   *
+   * @param size The current maximum number of steps through time.
+   */
+  void ResetCell(const size_t size);
+
+  /*
+   * Calculate the gradient using the output delta and the input activation.
+   *
+   * @param input The input parameter used for calculating the gradient.
+   * @param error The calculated error.
+   * @param gradient The calculated gradient.
+   */
+  template<typename InputType, typename ErrorType, typename GradientType>
+  void Gradient(InputType&& input,
+                ErrorType&& error,
+                GradientType&& gradient);
+
+  //! Get the maximum number of steps to backpropagate through time (BPTT).
+  size_t Rho() const { return rho; }
+  //! Modify the maximum number of steps to backpropagate through time (BPTT).
+  size_t& Rho() { return rho; }
+
+  //! Get the parameters.
+  OutputDataType const& Parameters() const { return weights; }
+  //! Modify the parameters.
+  OutputDataType& Parameters() { return weights; }
+
+  //! Get the input parameter.
+  InputDataType const& InputParameter() const { return inputParameter; }
+  //! Modify the input parameter.
+  InputDataType& InputParameter() { return inputParameter; }
+
+  //! Get the output parameter.
+  OutputDataType const& OutputParameter() const { return outputParameter; }
+  //! Modify the output parameter.
+  OutputDataType& OutputParameter() { return outputParameter; }
+
+  //! Get the delta.
+  OutputDataType const& Delta() const { return delta; }
+  //! Modify the delta.
+  OutputDataType& Delta() { return delta; }
+
+  //! Get the gradient.
+  OutputDataType const& Gradient() const { return grad; }
+  //! Modify the gradient.
+  OutputDataType& Gradient() { return grad; }
+
+  /**
+   * Serialize the layer
+   */
+  template<typename Archive>
+  void serialize(Archive& ar, const unsigned int /* version */);
+
+ private:
+  /**
+   * This speeds up the sigmoid operation by using an approximation.
+   *
+   * @param input The input data.
+   * @param sigmoid The matrix to store the sigmoid approximation into.
+   */
+  template<typename InputType, typename OutputType>
+  void FastSigmoid(InputType&& input, OutputType&& sigmoids)
+  {
+    for (size_t i = 0; i < input.n_elem; ++i)
+      sigmoids(i) = FastSigmoid(input(i));
+  }
+
+  /**
+   * Sigmoid approximation for the given sample.
+   *
+   * @param data The given data sample for the sigmoid approximation.
+   * @tparam The sigmoid approximation.
+   */
+  ElemType FastSigmoid(const InputElemType data)
+  {
+    ElemType x = 0.5 * data;
+    ElemType z;
+    if (x >= 0)
+    {
+      if (x < 1.7)
+        z = (1.5 * x / (1 + x));
+      else if (x < 3)
+        z = (0.935409070603099 + 0.0458812946797165 * (x - 1.7));
+      else
+        z = 0.99505475368673;
+    }
+    else
+    {
+      ElemType xx = -x;
+      if (xx < 1.7)
+        z = -(1.5 * xx / (1 + xx));
+      else if (xx < 3)
+        z = -(0.935409070603099 + 0.0458812946797165 * (xx - 1.7));
+      else
+        z = -0.99505475368673;
+    }
+
+    return 0.5 * (z + 1.0);
+  }
+
+  //! Locally-stored number of input units.
+  size_t inSize;
+
+  //! Locally-stored number of output units.
+  size_t outSize;
+
+  //! Number of steps to backpropagate through time (BPTT).
+  size_t rho;
+
+  //! Locally-stored number of forward steps.
+  size_t forwardStep;
+
+  //! Locally-stored number of backward steps.
+  size_t backwardStep;
+
+  //! Locally-stored number of gradient steps.
+  size_t gradientStep;
+
+  //! Locally-stored weight object.
+  OutputDataType weights;
+
+  //! Locally-stored previous output.
+  OutputDataType prevOutput;
+
+  //! Locally-stored batch size.
+  size_t batchSize;
+
+  //! Current batch step, alias for batchSize - 1.
+  size_t batchStep;
+
+  //! Current gradient step to keep track of the backpropagate through time
+  //! step.
+  size_t gradientStepIdx;
+
+  //! Locally-stored cell activation error.
+  OutputDataType cellActivationError;
+
+  //! Locally-stored delta object.
+  OutputDataType delta;
+
+  //! Locally-stored gradient object.
+  OutputDataType grad;
+
+  //! Locally-stored input parameter object.
+  InputDataType inputParameter;
+
+  //! Locally-stored output parameter object.
+  OutputDataType outputParameter;
+
+  //! Weights between the output and gate.
+  OutputDataType output2GateWeight;
+
+  //! Weights between the input and gate.
+  OutputDataType input2GateWeight;
+
+  //! Bias between the input and gate.
+  OutputDataType input2GateBias;
+
+  //! Locally-stored gate parameter.
+  OutputDataType gate;
+
+  //! Locally-stored gate activation.
+  OutputDataType gateActivation;
+
+  //! Locally-stored state activation.
+  OutputDataType stateActivation;
+
+  //! Locally-stored cell parameter.
+  OutputDataType cell;
+
+  //! Locally-stored cell activation error.
+  OutputDataType cellActivation;
+
+  //! Locally-stored foget gate error.
+  OutputDataType forgetGateError;
+
+  //! Locally-stored previous error.
+  OutputDataType prevError;
+
+  //! Locally-stored output parameters.
+  OutputDataType outParameter;
+
+  //! Locally-stored current rho size.
+  size_t rhoSize;
+
+  //! Current backpropagate through time steps.
+  size_t bpttSteps;
+}; // class FastLSTM
+
+} // namespace ann
+} // namespace mlpack
+
+// Include implementation.
+#include "fast_lstm_impl.hpp"
+
+#endif
diff --git a/src/mlpack/methods/ann/layer/fast_lstm_impl.hpp b/src/mlpack/methods/ann/layer/fast_lstm_impl.hpp
new file mode 100644
index 0000000000..3222275a4d
--- /dev/null
+++ b/src/mlpack/methods/ann/layer/fast_lstm_impl.hpp
@@ -0,0 +1,300 @@
+/**
+ * @file fast_lstm_impl.hpp
+ * @author Marcus Edel
+ *
+ * Implementation of the Fast LSTM class, which implements a fast lstm network
+ * layer.
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef MLPACK_METHODS_ANN_LAYER_FAST_LSTM_IMPL_HPP
+#define MLPACK_METHODS_ANN_LAYER_FAST_LSTM_IMPL_HPP
+
+// In case it hasn't yet been included.
+#include "fast_lstm.hpp"
+
+namespace mlpack {
+namespace ann /** Artificial Neural Network. */ {
+
+template<typename InputDataType, typename OutputDataType>
+FastLSTM<InputDataType, OutputDataType>::FastLSTM()
+{
+  // Nothing to do here.
+}
+
+template <typename InputDataType, typename OutputDataType>
+FastLSTM<InputDataType, OutputDataType>::FastLSTM(
+    const size_t inSize, const size_t outSize, const size_t rho) :
+    inSize(inSize),
+    outSize(outSize),
+    rho(rho),
+    forwardStep(0),
+    backwardStep(0),
+    gradientStep(0),
+    batchSize(0),
+    batchStep(0),
+    gradientStepIdx(0),
+    rhoSize(rho),
+    bpttSteps(0)
+{
+  // Weights for: input to gate layer (4 * outsize * inSize + 4 * outsize)
+  // and output to gate (4 * outSize).
+  weights.set_size(
+      4 * outSize * inSize + 4 * outSize + 4 * outSize * outSize, 1);
+}
+
+template<typename InputDataType, typename OutputDataType>
+void FastLSTM<InputDataType, OutputDataType>::Reset()
+{
+  // Set the weight parameter for the input to gate layer (linear layer) using
+  // the overall layer parameter matrix.
+  input2GateWeight = OutputDataType(weights.memptr(),
+      4 * outSize, inSize, false, false);
+  input2GateBias = OutputDataType(weights.memptr() + input2GateWeight.n_elem,
+      4 * outSize, 1, false, false);
+
+  // Set the weight parameter for the output to gate layer
+  // (linear no bias layer) using the overall layer parameter matrix.
+  output2GateWeight = OutputDataType(weights.memptr() + input2GateWeight.n_elem
+      + input2GateBias.n_elem, 4 * outSize, outSize, false, false);
+}
+
+template<typename InputDataType, typename OutputDataType>
+void FastLSTM<InputDataType, OutputDataType>::ResetCell(const size_t size)
+{
+  if (size == std::numeric_limits<size_t>::max())
+    return;
+
+  rhoSize = size;
+
+  if (batchSize == 0)
+    return;
+
+  bpttSteps = std::min(rho, rhoSize);
+  forwardStep = 0;
+  gradientStepIdx = 0;
+  backwardStep = batchSize * size - 1;
+  gradientStep = batchSize * size - 1;
+
+  const size_t rhoBatchSize = size * batchSize;
+  if (gate.is_empty() || gate.n_cols < rhoBatchSize)
+  {
+    gate.set_size(4 * outSize, rhoBatchSize);
+    gateActivation.set_size(outSize * 3, rhoBatchSize);
+    stateActivation.set_size(outSize, rhoBatchSize);
+    cellActivation.set_size(outSize, rhoBatchSize);
+    prevError.set_size(4 * outSize, batchSize);
+
+    if (prevOutput.is_empty())
+    {
+      prevOutput = arma::zeros<OutputDataType>(outSize, batchSize);
+      cell = arma::zeros(outSize, size * batchSize);
+      cellActivationError = arma::zeros<OutputDataType>(outSize, batchSize);
+      outParameter = arma::zeros<OutputDataType>(
+          outSize, (size + 1) * batchSize);
+    }
+    else
+    {
+      // To preserve the leading zeros, recreate the object according to given
+      // size specifications, while preserving the elements as well as the
+      // layout of the elements.
+      prevOutput.resize(outSize, batchSize);
+      cell.resize(outSize, size * batchSize);
+      cellActivationError.resize(outSize, batchSize);
+      outParameter.resize(outSize, (size + 1) * batchSize);
+    }
+  }
+}
+
+template<typename InputDataType, typename OutputDataType>
+template<typename InputType, typename OutputType>
+void FastLSTM<InputDataType, OutputDataType>::Forward(
+    InputType&& input, OutputType&& output)
+{
+  // Check if the batch size changed, the number of cols is defines the input
+  // batch size.
+  if (input.n_cols != batchSize)
+  {
+    batchSize = input.n_cols;
+    batchStep = batchSize - 1;
+    ResetCell(rhoSize);
+  }
+
+  gate.cols(forwardStep, forwardStep + batchStep) = input2GateWeight * input +
+      output2GateWeight * outParameter.cols(
+      forwardStep, forwardStep + batchStep);
+  gate.cols(forwardStep, forwardStep + batchStep).each_col() += input2GateBias;
+
+  FastSigmoid(std::move(
+      gate.submat(0, forwardStep, 3 * outSize - 1, forwardStep + batchStep)),
+      std::move(gateActivation.cols(forwardStep, forwardStep + batchStep)));
+
+  stateActivation.cols(forwardStep, forwardStep + batchStep) = arma::tanh(
+      gate.submat(3 * outSize, forwardStep, 4 * outSize - 1,
+      forwardStep + batchStep));
+
+  // Update the cell: cmul1 + cmul2
+  // where cmul1 is input gate * hidden state and
+  // cmul2 is forget gate * cell (prevCell).
+  if (forwardStep == 0)
+  {
+    cell.cols(forwardStep, forwardStep + batchStep) =
+        gateActivation.submat(0, forwardStep, outSize - 1,
+        forwardStep + batchStep) %
+        stateActivation.cols(forwardStep, forwardStep + batchStep);
+  }
+  else
+  {
+    cell.cols(forwardStep, forwardStep + batchStep) =
+        gateActivation.submat(0, forwardStep, outSize - 1,
+        forwardStep + batchStep) %
+        stateActivation.cols(forwardStep, forwardStep + batchStep) +
+        gateActivation.submat(2 * outSize, forwardStep, 3 * outSize - 1,
+        forwardStep + batchStep) %
+        cell.cols(forwardStep - batchSize, forwardStep - batchSize + batchStep);
+  }
+
+  cellActivation.cols(forwardStep, forwardStep + batchStep) =
+      arma::tanh(cell.cols(forwardStep, forwardStep + batchStep));
+
+  outParameter.cols(forwardStep + batchSize,
+      forwardStep + batchSize + batchStep) = cellActivation.cols(
+      forwardStep, forwardStep + batchStep) % gateActivation.submat(
+      outSize, forwardStep, 2 * outSize - 1, forwardStep + batchStep);
+
+  output = OutputType(outParameter.memptr() +
+      (forwardStep + batchSize) * outSize, outSize, batchSize, false, false);
+
+  forwardStep += batchSize;
+  if ((forwardStep / batchSize) == bpttSteps)
+  {
+    forwardStep = 0;
+  }
+}
+
+template<typename InputDataType, typename OutputDataType>
+template<typename InputType, typename ErrorType, typename GradientType>
+void FastLSTM<InputDataType, OutputDataType>::Backward(
+  const InputType&& /* input */, ErrorType&& gy, GradientType&& g)
+{
+  if (gradientStepIdx > 0)
+  {
+    gy += output2GateWeight.t() * prevError;
+  }
+
+  cellActivationError = gy % gateActivation.submat(outSize,
+      backwardStep - batchStep, 2 * outSize - 1, backwardStep) %
+      (1 - arma::pow(cellActivation.cols(backwardStep - batchStep,
+      backwardStep), 2));
+
+  if (gradientStepIdx > 0)
+    cellActivationError += forgetGateError;
+
+  forgetGateError = gateActivation.submat(2 * outSize,
+      backwardStep - batchStep, 3 * outSize - 1, backwardStep) %
+      cellActivationError;
+
+  if (backwardStep != 0)
+  {
+    prevError.submat(2 * outSize, 0, 3 * outSize - 1, batchStep) =
+        cell.cols((backwardStep - batchSize) - batchStep,
+        (backwardStep - batchSize)) % cellActivationError %
+        gateActivation.submat(2 * outSize, backwardStep - batchStep,
+        3 * outSize - 1, backwardStep) % (1.0 - gateActivation.submat(
+        2 * outSize, backwardStep - batchStep, 3 * outSize - 1, backwardStep));
+  }
+  else
+  {
+    prevError.submat(2 * outSize, 0, 3 * outSize - 1, batchStep).zeros();
+  }
+
+  prevError.submat(0, 0, outSize - 1, batchStep) =
+      stateActivation.cols(backwardStep - batchStep,
+      backwardStep) % cellActivationError % gateActivation.submat(
+      0, backwardStep - batchStep, outSize - 1, backwardStep) %
+      (1.0 - gateActivation.submat(
+      0, backwardStep - batchStep, outSize - 1, backwardStep));
+
+  prevError.submat(3 * outSize, 0, 4 * outSize - 1, batchStep) =
+      gateActivation.submat(0, backwardStep - batchStep,
+      outSize - 1, backwardStep) % cellActivationError % (1 - arma::pow(
+      stateActivation.cols(backwardStep - batchStep, backwardStep), 2));
+
+  prevError.submat(outSize, 0, 2 * outSize - 1, batchStep) =
+      cellActivation.cols(backwardStep - batchStep,
+      backwardStep) % gy % gateActivation.submat(
+       outSize, backwardStep - batchStep, 2 * outSize - 1, backwardStep) %
+      (1.0 - gateActivation.submat(
+      outSize, backwardStep - batchStep, 2 * outSize - 1, backwardStep));
+
+  g = input2GateWeight.t() * prevError;
+
+  backwardStep -= batchSize;
+  gradientStepIdx++;
+  if (gradientStepIdx == bpttSteps)
+  {
+    backwardStep = bpttSteps - 1;
+    gradientStepIdx = 0;
+  }
+}
+
+template<typename InputDataType, typename OutputDataType>
+template<typename InputType, typename ErrorType, typename GradientType>
+void FastLSTM<InputDataType, OutputDataType>::Gradient(
+    InputType&& input, ErrorType&& /* error */, GradientType&& gradient)
+{
+  // Gradient of the input to gate layer.
+  gradient.submat(0, 0, input2GateWeight.n_elem - 1, 0) =
+      arma::vectorise(prevError * input.t());
+
+  gradient.submat(input2GateWeight.n_elem, 0, input2GateWeight.n_elem +
+      input2GateBias.n_elem - 1, 0) = arma::sum(prevError, 1);
+
+  // Gradient of the output to gate layer.
+  gradient.submat(input2GateWeight.n_elem + input2GateBias.n_elem, 0,
+      gradient.n_elem - 1, 0) = arma::vectorise(prevError *
+      outParameter.cols(gradientStep - batchStep, gradientStep).t());
+
+  if (gradientStep == 0)
+  {
+    gradientStep = batchSize * bpttSteps - 1;
+  }
+  else
+  {
+    gradientStep -= batchSize;
+  }
+}
+
+template<typename InputDataType, typename OutputDataType>
+template<typename Archive>
+void FastLSTM<InputDataType, OutputDataType>::serialize(
+    Archive& ar, const unsigned int /* version */)
+{
+  ar & BOOST_SERIALIZATION_NVP(weights);
+  ar & BOOST_SERIALIZATION_NVP(inSize);
+  ar & BOOST_SERIALIZATION_NVP(outSize);
+  ar & BOOST_SERIALIZATION_NVP(rho);
+  ar & BOOST_SERIALIZATION_NVP(bpttSteps);
+  ar & BOOST_SERIALIZATION_NVP(batchSize);
+  ar & BOOST_SERIALIZATION_NVP(batchStep);
+  ar & BOOST_SERIALIZATION_NVP(forwardStep);
+  ar & BOOST_SERIALIZATION_NVP(backwardStep);
+  ar & BOOST_SERIALIZATION_NVP(gradientStep);
+  ar & BOOST_SERIALIZATION_NVP(gradientStepIdx);
+  ar & BOOST_SERIALIZATION_NVP(cell);
+  ar & BOOST_SERIALIZATION_NVP(stateActivation);
+  ar & BOOST_SERIALIZATION_NVP(gateActivation);
+  ar & BOOST_SERIALIZATION_NVP(gate);
+  ar & BOOST_SERIALIZATION_NVP(cellActivation);
+  ar & BOOST_SERIALIZATION_NVP(forgetGateError);
+  ar & BOOST_SERIALIZATION_NVP(prevError);
+  ar & BOOST_SERIALIZATION_NVP(outParameter);
+}
+
+} // namespace ann
+} // namespace mlpack
+
+#endif
diff --git a/src/mlpack/methods/ann/layer/gru.hpp b/src/mlpack/methods/ann/layer/gru.hpp
index e48587064e..aba74883f6 100644
--- a/src/mlpack/methods/ann/layer/gru.hpp
+++ b/src/mlpack/methods/ann/layer/gru.hpp
@@ -113,10 +113,12 @@ class GRU
                 arma::Mat<eT>&& /* gradient */);
 
   /*
-   * Resets the cell to accept a new input.
-   * This breaks the BPTT chain starts a new one.
+   * Resets the cell to accept a new input. This breaks the BPTT chain starts a
+   * new one.
+   *
+   * @param size The current maximum number of steps through time.
    */
-  void ResetCell();
+  void ResetCell(const size_t size);
 
   //! The value of the deterministic parameter.
   bool Deterministic() const { return deterministic; }
diff --git a/src/mlpack/methods/ann/layer/gru_impl.hpp b/src/mlpack/methods/ann/layer/gru_impl.hpp
index b1c19610cb..c795cef406 100644
--- a/src/mlpack/methods/ann/layer/gru_impl.hpp
+++ b/src/mlpack/methods/ann/layer/gru_impl.hpp
@@ -322,7 +322,7 @@ void GRU<InputDataType, OutputDataType>::Gradient(
 }
 
 template<typename InputDataType, typename OutputDataType>
-void GRU<InputDataType, OutputDataType>::ResetCell()
+void GRU<InputDataType, OutputDataType>::ResetCell(const size_t /* size */)
 {
   outParameter.clear();
   outParameter.push_back(std::move(arma::mat(allZeros.memptr(),
diff --git a/src/mlpack/methods/ann/layer/layer.hpp b/src/mlpack/methods/ann/layer/layer.hpp
index 6f0fbc8c1d..f47eadfd2f 100644
--- a/src/mlpack/methods/ann/layer/layer.hpp
+++ b/src/mlpack/methods/ann/layer/layer.hpp
@@ -22,6 +22,7 @@
 #include "linear_no_bias.hpp"
 #include "lstm.hpp"
 #include "gru.hpp"
+#include "fast_lstm.hpp"
 #include "recurrent.hpp"
 #include "recurrent_attention.hpp"
 #include "sequential.hpp"
diff --git a/src/mlpack/methods/ann/layer/layer_types.hpp b/src/mlpack/methods/ann/layer/layer_types.hpp
index 6ea36e3339..ddf377a84d 100644
--- a/src/mlpack/methods/ann/layer/layer_types.hpp
+++ b/src/mlpack/methods/ann/layer/layer_types.hpp
@@ -51,6 +51,7 @@ template<typename InputDataType, typename OutputDataType> class Linear;
 template<typename InputDataType, typename OutputDataType> class LinearNoBias;
 template<typename InputDataType, typename OutputDataType> class LSTM;
 template<typename InputDataType, typename OutputDataType> class GRU;
+template<typename InputDataType, typename OutputDataType> class FastLSTM;
 template<typename InputDataType, typename OutputDataType> class Recurrent;
 template<typename InputDataType, typename OutputDataType> class Sequential;
 template<typename InputDataType, typename OutputDataType> class VRClassReward;
@@ -105,6 +106,7 @@ using LayerTypes = boost::variant<
     Lookup<arma::mat, arma::mat>*,
     LSTM<arma::mat, arma::mat>*,
     GRU<arma::mat, arma::mat>*,
+    FastLSTM<arma::mat, arma::mat>*,
     MaxPooling<arma::mat, arma::mat>*,
     MeanPooling<arma::mat, arma::mat>*,
     MeanSquaredError<arma::mat, arma::mat>*,
diff --git a/src/mlpack/methods/ann/layer/lstm.hpp b/src/mlpack/methods/ann/layer/lstm.hpp
index fc0071395d..cd0f867136 100644
--- a/src/mlpack/methods/ann/layer/lstm.hpp
+++ b/src/mlpack/methods/ann/layer/lstm.hpp
@@ -2,8 +2,7 @@
  * @file lstm.hpp
  * @author Marcus Edel
  *
- * Definition of the LSTM class, which implements a lstm network
- * layer.
+ * Definition of the LSTM class, which implements a LSTM network layer.
  *
  * mlpack is free software; you may redistribute it and/or modify it under the
  * terms of the 3-clause BSD license.  You should have received a copy of the
@@ -14,16 +13,8 @@
 #define MLPACK_METHODS_ANN_LAYER_LSTM_HPP
 
 #include <mlpack/prereqs.hpp>
-
 #include <limits>
 
-#include "../visitor/delta_visitor.hpp"
-#include "../visitor/output_parameter_visitor.hpp"
-
-#include "layer_types.hpp"
-#include "add_merge.hpp"
-#include "sequential.hpp"
-
 namespace mlpack {
 namespace ann /** Artificial Neural Network. */ {
 
@@ -33,6 +24,32 @@ namespace ann /** Artificial Neural Network. */ {
  * This class allows specification of the type of the activation functions used
  * for the gates and cells and also of the type of the function used to
  * initialize and update the peephole weights.
+
+ * The implementation corresponds to the following algorithm:
+ *
+ * @f{eqnarray}{
+ * i &=& sigmoid(W \cdot x + W \cdot h + W \cdot c + b) \\
+ * f &=& sigmoid(W  \cdot x + W \cdot h + W \cdot c + b) \\
+ * z &=& tanh(W \cdot x + W \cdot h + b) \\
+ * c &=& f \cdot c + i \cdot z \\
+ * o &=& sigmoid(W \cdot x + W \cdot h + W \cdot c + b) \\
+ * h &=& o \cdot tanh(c)
+ * @f}
+ *
+ * For more information, see the following.
+ *
+ * @code
+ * @article{Graves2013,
+ *   author  = {Alex Graves and Abdel{-}rahman Mohamed and Geoffrey E. Hinton},
+ *   title   = {Speech Recognition with Deep Recurrent Neural Networks},
+ *   journal = CoRR},
+ *   year    = {2013},
+ *   url     = {http://arxiv.org/abs/1303.5778},
+ * }
+ * @endcode
+ *
+ * \see FastLSTM for a faster LSTM version which combines the calculation of the
+ * input, forget, output gates and hidden state in a single step.
  *
  * @tparam InputDataType Type of the input data (arma::colvec, arma::mat,
  *         arma::sp_mat or arma::cube).
@@ -60,11 +77,6 @@ class LSTM
        const size_t outSize,
        const size_t rho = std::numeric_limits<size_t>::max());
 
-  /**
-   * Delete the LSTM and the layers it holds.
-   */
-  ~LSTM();
-
   /**
    * Ordinary feed forward pass of a neural network, evaluating the function
    * f(x) by propagating the activity forward through f.
@@ -72,8 +84,8 @@ class LSTM
    * @param input Input data used for evaluating the specified function.
    * @param output Resulting output activation.
    */
-  template<typename eT>
-  void Forward(arma::Mat<eT>&& input, arma::Mat<eT>&& output);
+  template<typename InputType, typename OutputType>
+  void Forward(InputType&& input, OutputType&& output);
 
   /**
    * Ordinary feed backward pass of a neural network, calculating the function
@@ -84,10 +96,23 @@ class LSTM
    * @param gy The backpropagated error.
    * @param g The calculated gradient.
    */
-  template<typename eT>
-  void Backward(const arma::Mat<eT>&& /* input */,
-                arma::Mat<eT>&& gy,
-                arma::Mat<eT>&& g);
+  template<typename InputType, typename ErrorType, typename GradientType>
+  void Backward(const InputType&& input,
+                ErrorType&& gy,
+                GradientType&& g);
+
+  /*
+   * Reset the layer parameter.
+   */
+  void Reset();
+
+  /*
+   * Resets the cell to accept a new input. This breaks the BPTT chain starts a
+   * new one.
+   *
+   * @param size The current maximum number of steps through time.
+   */
+  void ResetCell(const size_t size);
 
   /*
    * Calculate the gradient using the output delta and the input activation.
@@ -96,21 +121,10 @@ class LSTM
    * @param error The calculated error.
    * @param gradient The calculated gradient.
    */
-  template<typename eT>
-  void Gradient(arma::Mat<eT>&& input,
-                arma::Mat<eT>&& /* error */,
-                arma::Mat<eT>&& /* gradient */);
-
-  /*
-   * Resets the cell to accept a new input.
-   * This breaks the BPTT chain starts a new one.
-   */
-  void ResetCell();
-
-  //! The value of the deterministic parameter.
-  bool Deterministic() const { return deterministic; }
-  //! Modify the value of the deterministic parameter.
-  bool& Deterministic() { return deterministic; }
+  template<typename InputType, typename ErrorType, typename GradientType>
+  void Gradient(InputType&& input,
+                ErrorType&& error,
+                GradientType&& gradient);
 
   //! Get the maximum number of steps to backpropagate through time (BPTT).
   size_t Rho() const { return rho; }
@@ -138,12 +152,9 @@ class LSTM
   OutputDataType& Delta() { return delta; }
 
   //! Get the gradient.
-  OutputDataType const& Gradient() const { return gradient; }
+  OutputDataType const& Gradient() const { return grad; }
   //! Modify the gradient.
-  OutputDataType& Gradient() { return gradient; }
-
-  //! Get the model modules.
-  std::vector<LayerTypes>& Model() { return network; }
+  OutputDataType& Gradient() { return grad; }
 
   /**
    * Serialize the layer
@@ -161,51 +172,6 @@ class LSTM
   //! Number of steps to backpropagate through time (BPTT).
   size_t rho;
 
-  //! Current batch size.
-  size_t batchSize;
-
-  //! Locally-stored weight object.
-  OutputDataType weights;
-
-  //! Locally-stored previous output.
-  std::list<arma::mat>::iterator prevOutput;
-
-  //! Locally-stored previous cell state.
-  std::list<arma::mat>::iterator prevCell;
-
-  //! Locally-stored input 2 gate module.
-  LayerTypes input2GateModule;
-
-  //! Locally-stored output 2 gate module.
-  LayerTypes output2GateModule;
-
-  //! Locally-stored input gate module.
-  LayerTypes inputGateModule;
-
-  //! Locally-stored hidden state module.
-  LayerTypes hiddenStateModule;
-
-  //! Locally-stored forget gate module.
-  LayerTypes forgetGateModule;
-
-  //! Locally-stored output gate module.
-  LayerTypes outputGateModule;
-
-  //! Locally-stored cell module.
-  LayerTypes cellModule;
-
-  //! Locally-stored cell activation module.
-  LayerTypes cellActivationModule;
-
-  //! Locally-stored output parameter visitor.
-  OutputParameterVisitor outputParameterVisitor;
-
-  //! Locally-stored delta visitor.
-  DeltaVisitor deltaVisitor;
-
-  //! Locally-stored list of network modules.
-  std::vector<LayerTypes> network;
-
   //! Locally-stored number of forward steps.
   size_t forwardStep;
 
@@ -215,41 +181,138 @@ class LSTM
   //! Locally-stored number of gradient steps.
   size_t gradientStep;
 
-  //! Locally-stored cell parameters.
-  std::list<arma::mat> cellParameter;
+  //! Locally-stored weight object.
+  OutputDataType weights;
 
-  //! Locally-stored output parameters.
-  std::list<arma::mat> outParameter;
+  //! Locally-stored previous output.
+  OutputDataType prevOutput;
 
-  //! Matrix of all zeroes to initialize the output and the cell
-  arma::mat allZeros;
+  //! Locally-stored batch size.
+  size_t batchSize;
 
-  //! Iterator pointed to the last cell output processed by backward
-  std::list<arma::mat>::iterator backIterator;
+  //! Current batch step, alias for batchSize - 1.
+  size_t batchStep;
 
-  //! Iterator pointed to the last output processed by gradient
-  std::list<arma::mat>::iterator gradIterator;
+  //! Current gradient step to keep track of the backpropagate through time
+  //! step.
+  size_t gradientStepIdx;
 
-  //! Locally-stored previous error.
-  arma::mat prevError;
-
-  //! Locally-stored foget gate error.
-  arma::mat forgetGateError;
-
-  //! If true dropout and scaling is disabled, see notes above.
-  bool deterministic;
+  //! Locally-stored cell activation error.
+  OutputDataType cellActivationError;
 
   //! Locally-stored delta object.
   OutputDataType delta;
 
   //! Locally-stored gradient object.
-  OutputDataType gradient;
+  OutputDataType grad;
 
   //! Locally-stored input parameter object.
   InputDataType inputParameter;
 
   //! Locally-stored output parameter object.
   OutputDataType outputParameter;
+
+  //! Weights between the output and input gate.
+  OutputDataType output2GateInputWeight;
+
+  //! Weights between the input and gate.
+  OutputDataType input2GateInputWeight;
+
+  //! Bias between the input and input gate.
+  OutputDataType input2GateInputBias;
+
+  //! Weights between the cell and input gate.
+  OutputDataType cell2GateInputWeight;
+
+  //! Weights between the output and forget gate.
+  OutputDataType output2GateForgetWeight;
+
+  //! Weights between the input and gate.
+  OutputDataType input2GateForgetWeight;
+
+  //! Bias between the input and gate.
+  OutputDataType input2GateForgetBias;
+
+  //! Bias between the input and gate.
+  OutputDataType cell2GateForgetWeight;
+
+  //! Weights between the output and gate.
+  OutputDataType output2GateOutputWeight;
+
+  //! Weights between the input and gate.
+  OutputDataType input2GateOutputWeight;
+
+  //! Bias between the input and gate.
+  OutputDataType input2GateOutputBias;
+
+  //! Weights between cell and output gate.
+  OutputDataType cell2GateOutputWeight;
+
+  //! Locally-stored input gate parameter.
+  OutputDataType inputGate;
+
+  //! Locally-stored forget gate parameter.
+  OutputDataType forgetGate;
+
+  //! Locally-stored hidden layer parameter.
+  OutputDataType hiddenLayer;
+
+  //! Locally-stored output gate parameter.
+  OutputDataType outputGate;
+
+  //! Locally-stored input gate activation.
+  OutputDataType inputGateActivation;
+
+  //! Locally-stored forget gate activation.
+  OutputDataType forgetGateActivation;
+
+  //! Locally-stored output gate activation.
+  OutputDataType outputGateActivation;
+
+  //! Locally-stored hidden layer activation.
+  OutputDataType hiddenLayerActivation;
+
+  //! Locally-stored input to hidden weight.
+  OutputDataType input2HiddenWeight;
+
+  //! Locally-stored input to hidden bias.
+  OutputDataType input2HiddenBias;
+
+  //! Locally-stored output to hidden weight.
+  OutputDataType output2HiddenWeight;
+
+  //! Locally-stored cell parameter.
+  OutputDataType cell;
+
+  //! Locally-stored cell activation error.
+  OutputDataType cellActivation;
+
+  //! Locally-stored forget gate error.
+  OutputDataType forgetGateError;
+
+  //! Locally-stored output gate error.
+  OutputDataType outputGateError;
+
+  //! Locally-stored previous error.
+  OutputDataType prevError;
+
+  //! Locally-stored output parameters.
+  OutputDataType outParameter;
+
+  //! Locally-stored input cell error parameter.
+  OutputDataType inputCellError;
+
+  //! Locally-stored input gate error.
+  OutputDataType inputGateError;
+
+  //! Locally-stored hidden layer error.
+  OutputDataType hiddenError;
+
+  //! Locally-stored current rho size.
+  size_t rhoSize;
+
+  //! Current backpropagate through time steps.
+  size_t bpttSteps;
 }; // class LSTM
 
 } // namespace ann
diff --git a/src/mlpack/methods/ann/layer/lstm_impl.hpp b/src/mlpack/methods/ann/layer/lstm_impl.hpp
index 938363fdb1..5abf58be5d 100644
--- a/src/mlpack/methods/ann/layer/lstm_impl.hpp
+++ b/src/mlpack/methods/ann/layer/lstm_impl.hpp
@@ -2,8 +2,7 @@
  * @file lstm_impl.hpp
  * @author Marcus Edel
  *
- * Implementation of the LSTM class, which implements a lstm network
- * layer.
+ * Implementation of the LSTM class, which implements a lstm network layer.
  *
  * mlpack is free software; you may redistribute it and/or modify it under the
  * terms of the 3-clause BSD license.  You should have received a copy of the
@@ -16,10 +15,6 @@
 // In case it hasn't yet been included.
 #include "lstm.hpp"
 
-#include "../visitor/forward_visitor.hpp"
-#include "../visitor/backward_visitor.hpp"
-#include "../visitor/gradient_visitor.hpp"
-
 namespace mlpack {
 namespace ann /** Artificial Neural Network. */ {
 
@@ -31,319 +26,402 @@ LSTM<InputDataType, OutputDataType>::LSTM()
 
 template <typename InputDataType, typename OutputDataType>
 LSTM<InputDataType, OutputDataType>::LSTM(
-    const size_t inSize,
-    const size_t outSize,
-    const size_t rho) :
+    const size_t inSize, const size_t outSize, const size_t rho) :
     inSize(inSize),
     outSize(outSize),
     rho(rho),
-    batchSize(1),
     forwardStep(0),
     backwardStep(0),
     gradientStep(0),
-    deterministic(false)
+    batchSize(0),
+    batchStep(0),
+    gradientStepIdx(0),
+    rhoSize(rho),
+    bpttSteps(0)
 {
-  input2GateModule = new Linear<>(inSize, 4 * outSize);
-  output2GateModule = new LinearNoBias<>(outSize, 4 * outSize);
-
-  network.push_back(input2GateModule);
-  network.push_back(output2GateModule);
-
-  inputGateModule = new SigmoidLayer<>();
-  hiddenStateModule = new TanHLayer<>();
-  forgetGateModule = new SigmoidLayer<>();
-  outputGateModule = new SigmoidLayer<>();
-
-  network.push_back(inputGateModule);
-  network.push_back(hiddenStateModule);
-  network.push_back(forgetGateModule);
-  network.push_back(outputGateModule);
-
-  cellModule = new IdentityLayer<>();
-  cellActivationModule = new TanHLayer<>();
-
-  network.push_back(cellModule);
-  network.push_back(cellActivationModule);
-
-  prevError = arma::zeros<arma::mat>(4 * outSize, batchSize);
-
-  allZeros = arma::zeros<arma::mat>(outSize, batchSize);
-
-  outParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-    allZeros.n_rows, allZeros.n_cols, false, false)));
-
-  cellParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-    allZeros.n_rows, allZeros.n_cols, false, false)));
-
-  prevOutput = outParameter.begin();
-  prevCell = cellParameter.begin();
-
-  backIterator = cellParameter.end();
-  gradIterator = outParameter.end();
+  weights.set_size(4 * outSize * inSize + 7 * outSize +
+      4 * outSize * outSize, 1);
 }
 
 template<typename InputDataType, typename OutputDataType>
-LSTM<InputDataType, OutputDataType>::~LSTM()
+void LSTM<InputDataType, OutputDataType>::ResetCell(const size_t size)
 {
-  boost::apply_visitor(DeleteVisitor(), input2GateModule);
-  boost::apply_visitor(DeleteVisitor(), output2GateModule);
-  boost::apply_visitor(DeleteVisitor(), inputGateModule);
-  boost::apply_visitor(DeleteVisitor(), hiddenStateModule);
-  boost::apply_visitor(DeleteVisitor(), forgetGateModule);
-  boost::apply_visitor(DeleteVisitor(), outputGateModule);
-  boost::apply_visitor(DeleteVisitor(), cellModule);
-  boost::apply_visitor(DeleteVisitor(), cellActivationModule);
-}
+  if (size == std::numeric_limits<size_t>::max())
+    return;
 
-template<typename InputDataType, typename OutputDataType>
-template<typename eT>
-void LSTM<InputDataType, OutputDataType>::Forward(
-    arma::Mat<eT>&& input, arma::Mat<eT>&& output)
-{
-  if (input.n_cols != batchSize)
+  rhoSize = size;
+
+  if (batchSize == 0)
+    return;
+
+  bpttSteps = std::min(rho, rhoSize);
+  forwardStep = 0;
+  gradientStepIdx = 0;
+  backwardStep = batchSize * size - 1;
+  gradientStep = batchSize * size - 1;
+
+  const size_t rhoBatchSize = size * batchSize;
+  if (inputGate.is_empty() || inputGate.n_cols < rhoBatchSize)
   {
-    batchSize = input.n_cols;
-    prevError.resize(3 * outSize, batchSize);
-  }
+    inputGate.set_size(outSize, rhoBatchSize);
+    forgetGate.set_size(outSize, rhoBatchSize);
+    hiddenLayer.set_size(outSize, rhoBatchSize);
+    outputGate.set_size(outSize, rhoBatchSize);
 
-  boost::apply_visitor(ForwardVisitor(std::move(input), std::move(
-      boost::apply_visitor(outputParameterVisitor, input2GateModule))),
-      input2GateModule);
+    inputGateActivation.set_size(outSize, rhoBatchSize);
+    forgetGateActivation.set_size(outSize, rhoBatchSize);
+    outputGateActivation.set_size(outSize, rhoBatchSize);
+    hiddenLayerActivation.set_size(outSize, rhoBatchSize);
 
-  boost::apply_visitor(ForwardVisitor(std::move(*prevOutput), std::move(
-      boost::apply_visitor(outputParameterVisitor, output2GateModule))),
-      output2GateModule);
+    cellActivation.set_size(outSize, rhoBatchSize);
+    prevError.set_size(4 * outSize, batchSize);
 
-  output = boost::apply_visitor(outputParameterVisitor, input2GateModule) +
-      boost::apply_visitor(outputParameterVisitor, output2GateModule);
-
-  boost::apply_visitor(ForwardVisitor(std::move(output.submat(
-      0, 0, 1 * outSize - 1, batchSize - 1)), std::move(boost::apply_visitor(
-      outputParameterVisitor, inputGateModule))), inputGateModule);
-
-  boost::apply_visitor(ForwardVisitor(std::move(output.submat(
-      1 * outSize, 0, 2 * outSize - 1, batchSize - 1)), std::move(
-      boost::apply_visitor(outputParameterVisitor, hiddenStateModule))),
-      hiddenStateModule);
-
-  boost::apply_visitor(ForwardVisitor(std::move(output.submat(
-      2 * outSize, 0, 3 * outSize - 1, batchSize - 1)), std::move(
-      boost::apply_visitor(outputParameterVisitor, forgetGateModule))),
-      forgetGateModule);
-
-  boost::apply_visitor(ForwardVisitor(std::move(output.submat(
-      3 * outSize, 0, 4 * outSize - 1, batchSize - 1)), std::move(
-      boost::apply_visitor(outputParameterVisitor, outputGateModule))),
-      outputGateModule);
-
-  // Update the cell (nextCell): cmul1 + cmul2
-  // where cmul1 is input gate * hidden state and
-  // cmul2 is forget gate * cell (prevCell).
-  arma::mat tempPrevCell = (boost::apply_visitor(outputParameterVisitor,
-      inputGateModule) % boost::apply_visitor(outputParameterVisitor,
-      hiddenStateModule)) + (boost::apply_visitor(outputParameterVisitor,
-      forgetGateModule) % *prevCell);
-
-  boost::apply_visitor(ForwardVisitor(std::move(tempPrevCell), std::move(
-      boost::apply_visitor(outputParameterVisitor, cellModule))), cellModule);
-
-  boost::apply_visitor(ForwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, cellModule)), std::move(boost::apply_visitor(
-      outputParameterVisitor, cellActivationModule))), cellActivationModule);
-
-  output = boost::apply_visitor(outputParameterVisitor,
-      cellActivationModule) % boost::apply_visitor(outputParameterVisitor,
-      outputGateModule);
-
-  forwardStep++;
-  if (forwardStep == rho)
-  {
-    forwardStep = 0;
-    if (!deterministic)
+    if (cell.is_empty())
     {
-      outParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-        allZeros.n_rows, allZeros.n_cols, false, false)));
-
-      cellParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-        allZeros.n_rows, allZeros.n_cols, false, false)));
-
-      prevOutput = --outParameter.end();
-      prevCell = --cellParameter.end();
+      cell = arma::zeros(outSize, size * batchSize);
+      outParameter = arma::zeros<OutputDataType>(
+          outSize, (size + 1) * batchSize);
     }
     else
     {
-      *prevOutput = std::move(arma::mat(allZeros.memptr(),
-        allZeros.n_rows, allZeros.n_cols, false, false));
-
-      *prevCell = std::move(arma::mat(allZeros.memptr(),
-        allZeros.n_rows, allZeros.n_cols, false, false));
+      // To preserve the leading zeros, recreate the object according to given
+      // size specifications, while preserving the elements as well as the
+      // layout of the elements.
+      cell.resize(outSize, size * batchSize);
+      outParameter.resize(outSize, (size + 1) * batchSize);
     }
   }
-  else if (!deterministic)
+}
+
+template<typename InputDataType, typename OutputDataType>
+void LSTM<InputDataType, OutputDataType>::Reset()
+{
+  // Set the weight parameter for the output gate.
+  input2GateOutputWeight = OutputDataType(weights.memptr(), outSize, inSize,
+      false, false);
+  input2GateOutputBias = OutputDataType(weights.memptr() +
+      input2GateOutputWeight.n_elem, outSize, 1, false, false);
+  size_t offset = input2GateOutputWeight.n_elem + input2GateOutputBias.n_elem;
+
+  // Set the weight parameter for the forget gate.
+  input2GateForgetWeight = OutputDataType(weights.memptr() + offset,
+      outSize, inSize, false, false);
+  input2GateForgetBias = OutputDataType(weights.memptr() +
+      offset + input2GateForgetWeight.n_elem, outSize, 1, false, false);
+  offset += input2GateForgetWeight.n_elem + input2GateForgetBias.n_elem;
+
+  // Set the weight parameter for the input gate.
+  input2GateInputWeight = OutputDataType(weights.memptr() +
+      offset, outSize, inSize, false, false);
+  input2GateInputBias = OutputDataType(weights.memptr() +
+      offset + input2GateInputWeight.n_elem, outSize, 1, false, false);
+  offset += input2GateInputWeight.n_elem + input2GateInputBias.n_elem;
+
+  // Set the weight parameter for the hidden gate.
+  input2HiddenWeight = OutputDataType(weights.memptr() +
+      offset, outSize, inSize, false, false);
+  input2HiddenBias = OutputDataType(weights.memptr() +
+      offset + input2HiddenWeight.n_elem, outSize, 1, false, false);
+  offset += input2HiddenWeight.n_elem + input2HiddenBias.n_elem;
+
+  // Set the weight parameter for the output multiplication.
+  output2GateOutputWeight = OutputDataType(weights.memptr() +
+      offset, outSize, outSize, false, false);
+  offset += output2GateOutputWeight.n_elem;
+
+  // Set the weight parameter for the output multiplication.
+  output2GateForgetWeight = OutputDataType(weights.memptr() +
+      offset, outSize, outSize, false, false);
+  offset += output2GateForgetWeight.n_elem;
+
+  // Set the weight parameter for the input multiplication.
+  output2GateInputWeight = OutputDataType(weights.memptr() +
+      offset, outSize, outSize, false, false);
+  offset += output2GateInputWeight.n_elem;
+
+  // Set the weight parameter for the hidden multiplication.
+  output2HiddenWeight = OutputDataType(weights.memptr() +
+      offset, outSize, outSize, false, false);
+  offset += output2HiddenWeight.n_elem;
+
+  // Set the weight parameter for the cell multiplication.
+  cell2GateOutputWeight = OutputDataType(weights.memptr() +
+      offset, outSize, 1, false, false);
+  offset += cell2GateOutputWeight.n_elem;
+
+  // Set the weight parameter for the cell - forget gate multiplication.
+  cell2GateForgetWeight = OutputDataType(weights.memptr() +
+      offset, outSize, 1, false, false);
+  offset += cell2GateOutputWeight.n_elem;
+
+  // Set the weight parameter for the cell - input gate multiplication.
+  cell2GateInputWeight = OutputDataType(weights.memptr() +
+      offset, outSize, 1, false, false);
+}
+
+template<typename InputDataType, typename OutputDataType>
+template<typename InputType, typename OutputType>
+void LSTM<InputDataType, OutputDataType>::Forward(
+    InputType&& input, OutputType&& output)
+{
+  // Check if the batch size changed, the number of cols is defines the input
+  // batch size.
+  if (input.n_cols != batchSize)
   {
-    outParameter.push_back(output);
-    cellParameter.push_back(std::move(tempPrevCell));
-    prevOutput = --outParameter.end();
-    prevCell = --cellParameter.end();
+    batchSize = input.n_cols;
+    batchStep = batchSize - 1;
+    ResetCell(rhoSize);
+  }
+
+  inputGate.cols(forwardStep, forwardStep + batchStep) = input2GateInputWeight *
+      input + output2GateInputWeight * outParameter.cols(forwardStep,
+      forwardStep + batchStep);
+  inputGate.cols(forwardStep, forwardStep + batchStep).each_col() +=
+      input2GateInputBias;
+
+  forgetGate.cols(forwardStep, forwardStep + batchStep) = input2GateForgetWeight
+      * input + output2GateForgetWeight * outParameter.cols(
+      forwardStep, forwardStep + batchStep);
+  forgetGate.cols(forwardStep, forwardStep + batchStep).each_col() +=
+      input2GateForgetBias;
+
+  if (forwardStep > 0)
+  {
+    inputGate.cols(forwardStep, forwardStep + batchStep) +=
+        cell2GateInputWeight % cell.cols(forwardStep - batchSize,
+        forwardStep - batchSize + batchStep);
+
+    forgetGate.cols(forwardStep, forwardStep + batchStep) +=
+        cell2GateForgetWeight % cell.cols(forwardStep - batchSize,
+        forwardStep - batchSize + batchStep);
+  }
+
+  inputGateActivation.cols(forwardStep, forwardStep + batchStep) = 1.0 /
+      (1 + arma::exp(-inputGate.cols(forwardStep, forwardStep + batchStep)));
+
+  forgetGateActivation.cols(forwardStep, forwardStep + batchStep) = 1.0 /
+      (1 + arma::exp(-forgetGate.cols(forwardStep, forwardStep + batchStep)));
+
+  hiddenLayer.cols(forwardStep, forwardStep + batchStep) = input2HiddenWeight *
+      input + output2HiddenWeight * outParameter.cols(
+      forwardStep, forwardStep + batchStep);
+
+  hiddenLayer.cols(forwardStep, forwardStep + batchStep).each_col() +=
+      input2HiddenBias;
+
+  hiddenLayerActivation.cols(forwardStep, forwardStep + batchStep) =
+      arma::tanh(hiddenLayer.cols(forwardStep, forwardStep + batchStep));
+
+  if (forwardStep == 0)
+  {
+    cell.cols(forwardStep, forwardStep + batchStep) =
+        inputGateActivation.cols(forwardStep, forwardStep + batchStep) %
+        hiddenLayerActivation.cols(forwardStep, forwardStep + batchStep);
   }
   else
   {
-    if (forwardStep == 1)
-    {
-      outParameter.clear();
-      cellParameter.clear();
+    cell.cols(forwardStep, forwardStep + batchStep) =
+        forgetGateActivation.cols(forwardStep, forwardStep + batchStep) %
+        cell.cols(forwardStep - batchSize, forwardStep - batchSize + batchStep)
+        + inputGateActivation.cols(forwardStep, forwardStep + batchStep) %
+        hiddenLayerActivation.cols(forwardStep, forwardStep + batchStep);
+  }
 
-      outParameter.push_back(output);
-      cellParameter.push_back(std::move(tempPrevCell));
+  outputGate.cols(forwardStep, forwardStep + batchStep) = input2GateOutputWeight
+      * input + output2GateOutputWeight * outParameter.cols(
+      forwardStep, forwardStep + batchStep) + cell.cols(forwardStep,
+      forwardStep + batchStep).each_col() % cell2GateOutputWeight;
 
-      prevOutput = outParameter.begin();
-      prevCell = cellParameter.begin();
-    }
-    else
-    {
-      *prevOutput = output;
-      *prevCell = std::move(tempPrevCell);
-    }
+  outputGate.cols(forwardStep, forwardStep + batchStep).each_col() +=
+      input2GateOutputBias;
+
+  outputGateActivation.cols(forwardStep, forwardStep + batchStep) = 1.0 /
+      (1 + arma::exp(-outputGate.cols(forwardStep, forwardStep + batchStep)));
+
+  cellActivation.cols(forwardStep, forwardStep + batchStep) =
+      arma::tanh(cell.cols(forwardStep, forwardStep + batchStep));
+
+  outParameter.cols(forwardStep + batchSize,
+      forwardStep + batchSize + batchStep) =
+      cellActivation.cols(forwardStep, forwardStep + batchStep) %
+      outputGateActivation.cols(forwardStep, forwardStep + batchStep);
+
+  output = OutputType(outParameter.memptr() +
+      (forwardStep + batchSize) * outSize, outSize, batchSize, false, false);
+
+  forwardStep += batchSize;
+  if ((forwardStep / batchSize) == bpttSteps)
+  {
+    forwardStep = 0;
   }
 }
 
 template<typename InputDataType, typename OutputDataType>
-template<typename eT>
+template<typename InputType, typename ErrorType, typename GradientType>
 void LSTM<InputDataType, OutputDataType>::Backward(
-  const arma::Mat<eT>&& input, arma::Mat<eT>&& gy, arma::Mat<eT>&& g)
+  const InputType&& /* input */, ErrorType&& gy, GradientType&& g)
 {
-  if (input.n_cols != batchSize)
+  if (gradientStepIdx > 0)
   {
-    batchSize = input.n_cols;
-    prevError.resize(3 * outSize, batchSize);
+    gy += prevError;
   }
 
-  if ((outParameter.size() - backwardStep  - 1) % rho != 0 && backwardStep != 0)
+  outputGateError =
+      gy % cellActivation.cols(backwardStep - batchStep, backwardStep) %
+      (outputGateActivation.cols(backwardStep - batchStep, backwardStep) %
+      (1.0 - outputGateActivation.cols(backwardStep - batchStep,
+      backwardStep)));
+
+  OutputDataType cellError = gy %
+      outputGateActivation.cols(backwardStep - batchStep, backwardStep) %
+      (1 - arma::pow(cellActivation.cols(backwardStep -
+      batchStep, backwardStep), 2)) + outputGateError.each_col() %
+      cell2GateOutputWeight;
+
+  if (gradientStepIdx > 0)
   {
-    gy += boost::apply_visitor(deltaVisitor, output2GateModule);
+    cellError += inputCellError;
   }
 
-  if (backIterator == cellParameter.end())
+  if (backwardStep != 0)
   {
-    backIterator = --(--cellParameter.end());
+    forgetGateError = cell.cols((backwardStep - batchSize) - batchStep,
+      (backwardStep - batchSize)) % cellError % (forgetGateActivation.cols(
+      backwardStep - batchStep, backwardStep) % (1.0 -
+      forgetGateActivation.cols(backwardStep - batchStep, backwardStep)));
+  }
+  else
+  {
+    forgetGateError.zeros();
   }
 
-  arma::mat g1 = boost::apply_visitor(outputParameterVisitor,
-      cellActivationModule) % gy;
+  inputGateError = hiddenLayerActivation.cols(backwardStep - batchStep,
+      backwardStep) % cellError %
+      (inputGateActivation.cols(backwardStep - batchStep, backwardStep) %
+      (1.0 - inputGateActivation.cols(backwardStep - batchStep, backwardStep)));
 
-  arma::mat g2 = boost::apply_visitor(outputParameterVisitor,
-      outputGateModule) % gy;
+  hiddenError = inputGateActivation.cols(backwardStep - batchStep,
+      backwardStep) % cellError % (1 - arma::pow(hiddenLayerActivation.cols(
+      backwardStep - batchStep, backwardStep), 2));
 
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, cellActivationModule)), std::move(g2),
-      std::move(boost::apply_visitor(deltaVisitor, cellActivationModule))),
-      cellActivationModule);
+  inputCellError = forgetGateActivation.cols(backwardStep - batchStep,
+      backwardStep) %cellError + forgetGateError.each_col() %
+      cell2GateForgetWeight + inputGateError.each_col() % cell2GateInputWeight;
 
-  arma::mat cellActivationError = boost::apply_visitor(deltaVisitor,
-      cellActivationModule);
+  g = input2GateInputWeight.t() * inputGateError +
+      input2HiddenWeight.t() * hiddenError +
+      input2GateForgetWeight.t() * forgetGateError +
+      input2GateOutputWeight.t() * outputGateError;
 
-  if (backwardStep > 0)
+  prevError = output2GateOutputWeight.t() * outputGateError +
+      output2GateForgetWeight.t() * forgetGateError +
+      output2GateInputWeight.t() * inputGateError +
+      output2HiddenWeight.t() * hiddenError;
+
+  backwardStep -= batchSize;
+  gradientStepIdx++;
+  if (gradientStepIdx == bpttSteps)
   {
-    cellActivationError += forgetGateError;
+    backwardStep = bpttSteps - 1;
+    gradientStepIdx = 0;
   }
-
-  arma::mat g4 = boost::apply_visitor(outputParameterVisitor,
-      inputGateModule) % cellActivationError;
-
-  arma::mat g5 = boost::apply_visitor(outputParameterVisitor,
-      hiddenStateModule) % cellActivationError;
-
-  forgetGateError = boost::apply_visitor(outputParameterVisitor,
-      forgetGateModule) % cellActivationError;
-
-  arma::mat g7 = *backIterator % cellActivationError;
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, inputGateModule)), std::move(g5),
-      std::move(boost::apply_visitor(deltaVisitor, inputGateModule))),
-      inputGateModule);
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, hiddenStateModule)), std::move(g4),
-      std::move(boost::apply_visitor(deltaVisitor, hiddenStateModule))),
-      hiddenStateModule);
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, forgetGateModule)), std::move(g7),
-      std::move(boost::apply_visitor(deltaVisitor, forgetGateModule))),
-      forgetGateModule);
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, outputGateModule)), std::move(g1),
-      std::move(boost::apply_visitor(deltaVisitor, outputGateModule))),
-      outputGateModule);
-
-  prevError.submat(0, 0, 1 * outSize - 1, batchSize - 1) = boost::apply_visitor(
-      deltaVisitor, inputGateModule);
-  prevError.submat(1 * outSize, 0, 2 * outSize - 1, batchSize - 1) =
-      boost::apply_visitor(deltaVisitor, hiddenStateModule);
-  prevError.submat(2 * outSize, 0, 3 * outSize - 1, batchSize - 1) =
-      boost::apply_visitor(deltaVisitor, forgetGateModule);
-  prevError.submat(3 * outSize, 0, 4 * outSize - 1, batchSize - 1) =
-      boost::apply_visitor(deltaVisitor, outputGateModule);
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, input2GateModule)), std::move(prevError),
-      std::move(boost::apply_visitor(deltaVisitor, input2GateModule))),
-      input2GateModule);
-
-  boost::apply_visitor(BackwardVisitor(std::move(boost::apply_visitor(
-      outputParameterVisitor, output2GateModule)), std::move(prevError),
-      std::move(boost::apply_visitor(deltaVisitor, output2GateModule))),
-      output2GateModule);
-
-  backwardStep++;
-  backIterator--;
-
-  g = boost::apply_visitor(deltaVisitor, input2GateModule);
 }
 
 template<typename InputDataType, typename OutputDataType>
-template<typename eT>
+template<typename InputType, typename ErrorType, typename GradientType>
 void LSTM<InputDataType, OutputDataType>::Gradient(
-    arma::Mat<eT>&& input,
-    arma::Mat<eT>&& /* error */,
-    arma::Mat<eT>&& /* gradient */)
+    InputType&& input, ErrorType&& /* error */, GradientType&& gradient)
 {
-  if (gradIterator == outParameter.end())
+  // Input2GateOutputWeight and input2GateOutputBias gradients.
+  gradient.submat(0, 0, input2GateOutputWeight.n_elem - 1, 0) =
+      arma::vectorise(outputGateError * input.t());
+  gradient.submat(input2GateOutputWeight.n_elem, 0,
+      input2GateOutputWeight.n_elem + input2GateOutputBias.n_elem - 1, 0) =
+      arma::sum(outputGateError, 1);
+  size_t offset = input2GateOutputWeight.n_elem + input2GateOutputBias.n_elem;
+
+  // Input2GateForgetWeight and input2GateForgetBias gradients.
+  gradient.submat(offset, 0, offset + input2GateForgetWeight.n_elem - 1, 0) =
+      arma::vectorise(forgetGateError * input.t());
+  gradient.submat(offset + input2GateForgetWeight.n_elem, 0,
+      offset + input2GateForgetWeight.n_elem +
+      input2GateForgetBias.n_elem - 1, 0) = arma::sum(forgetGateError, 1);
+  offset += input2GateForgetWeight.n_elem + input2GateForgetBias.n_elem;
+
+  // Input2GateInputWeight and input2GateInputBias gradients.
+  gradient.submat(offset, 0, offset + input2GateInputWeight.n_elem - 1, 0) =
+      arma::vectorise(inputGateError * input.t());
+  gradient.submat(offset + input2GateInputWeight.n_elem, 0,
+      offset + input2GateInputWeight.n_elem +
+      input2GateInputBias.n_elem - 1, 0) = arma::sum(inputGateError, 1);
+  offset += input2GateInputWeight.n_elem + input2GateInputBias.n_elem;
+
+  // Input2HiddenWeight and input2HiddenBias gradients.
+  gradient.submat(offset, 0, offset + input2HiddenWeight.n_elem - 1, 0) =
+      arma::vectorise(hiddenError * input.t());
+  gradient.submat(offset + input2HiddenWeight.n_elem, 0,
+      offset + input2HiddenWeight.n_elem + input2HiddenBias.n_elem - 1, 0) =
+      arma::sum(hiddenError, 1);
+  offset += input2HiddenWeight.n_elem + input2HiddenBias.n_elem;
+
+  // Output2GateOutputWeight gradients.
+  gradient.submat(offset, 0, offset + output2GateOutputWeight.n_elem - 1, 0) =
+      arma::vectorise(outputGateError *
+      outParameter.cols(gradientStep - batchStep, gradientStep).t());
+  offset += output2GateOutputWeight.n_elem;
+
+  // Output2GateForgetWeight gradients.
+  gradient.submat(offset, 0, offset + output2GateForgetWeight.n_elem - 1, 0) =
+      arma::vectorise(forgetGateError *
+      outParameter.cols(gradientStep - batchStep, gradientStep).t());
+  offset += output2GateForgetWeight.n_elem;
+
+  // Output2GateInputWeight gradients.
+  gradient.submat(offset, 0, offset + output2GateInputWeight.n_elem - 1, 0) =
+      arma::vectorise(inputGateError *
+      outParameter.cols(gradientStep - batchStep, gradientStep).t());
+  offset += output2GateInputWeight.n_elem;
+
+  // Output2HiddenWeight gradients.
+  gradient.submat(offset, 0, offset + output2HiddenWeight.n_elem - 1, 0) =
+      arma::vectorise(hiddenError *
+      outParameter.cols(gradientStep - batchStep, gradientStep).t());
+  offset += output2HiddenWeight.n_elem;
+
+  // Cell2GateOutputWeight gradients.
+  gradient.submat(offset, 0, offset + cell2GateOutputWeight.n_elem - 1, 0) =
+      arma::sum(outputGateError %
+      cell.cols(gradientStep - batchStep, gradientStep), 1);
+  offset += cell2GateOutputWeight.n_elem;
+
+  // Cell2GateForgetWeight and cell2GateInputWeight gradients.
+  if (gradientStep != 0)
   {
-    gradIterator = --(--outParameter.end());
+    gradient.submat(offset, 0, offset + cell2GateForgetWeight.n_elem - 1, 0) =
+        arma::sum(forgetGateError % cell.cols(gradientStep - batchStep -
+        batchSize, gradientStep - batchSize), 1);
+    gradient.submat(offset + cell2GateForgetWeight.n_elem, 0, offset +
+        cell2GateForgetWeight.n_elem + cell2GateInputWeight.n_elem - 1, 0) =
+        arma::sum(inputGateError % cell.cols(gradientStep - batchStep -
+        batchSize, gradientStep - batchSize), 1);
+  }
+  else
+  {
+    gradient.submat(offset, 0, offset +
+        cell2GateForgetWeight.n_elem - 1, 0).zeros();
+    gradient.submat(offset + cell2GateForgetWeight.n_elem, 0, offset +
+        cell2GateForgetWeight.n_elem +
+        cell2GateInputWeight.n_elem - 1, 0).zeros();
   }
 
-  boost::apply_visitor(GradientVisitor(std::move(input), std::move(prevError)),
-      input2GateModule);
-
-  boost::apply_visitor(GradientVisitor(
-      std::move(*gradIterator),
-      std::move(prevError)), output2GateModule);
-
-  gradIterator--;
-}
-
-template<typename InputDataType, typename OutputDataType>
-void LSTM<InputDataType, OutputDataType>::ResetCell()
-{
-  outParameter.clear();
-  outParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-    allZeros.n_rows, allZeros.n_cols, false, false)));
-
-  cellParameter.clear();
-  cellParameter.push_back(std::move(arma::mat(allZeros.memptr(),
-    allZeros.n_rows, allZeros.n_cols, false, false)));
-
-  prevOutput = outParameter.begin();
-  prevCell = cellParameter.begin();
-
-  backIterator = cellParameter.end();
-  gradIterator = outParameter.end();
-
-  forwardStep = 0;
-  backwardStep = 0;
+  if (gradientStep == 0)
+  {
+    gradientStep = batchSize * bpttSteps - 1;
+  }
+  else
+  {
+    gradientStep -= batchSize;
+  }
 }
 
 template<typename InputDataType, typename OutputDataType>
@@ -351,31 +429,25 @@ template<typename Archive>
 void LSTM<InputDataType, OutputDataType>::serialize(
     Archive& ar, const unsigned int /* version */)
 {
-  // Clear old layers, if needed.
-  if (Archive::is_loading::value)
-  {
-    boost::apply_visitor(DeleteVisitor(), input2GateModule);
-    boost::apply_visitor(DeleteVisitor(), output2GateModule);
-    boost::apply_visitor(DeleteVisitor(), inputGateModule);
-    boost::apply_visitor(DeleteVisitor(), hiddenStateModule);
-    boost::apply_visitor(DeleteVisitor(), forgetGateModule);
-    boost::apply_visitor(DeleteVisitor(), outputGateModule);
-    boost::apply_visitor(DeleteVisitor(), cellModule);
-    boost::apply_visitor(DeleteVisitor(), cellActivationModule);
-  }
-
+  ar & BOOST_SERIALIZATION_NVP(weights);
   ar & BOOST_SERIALIZATION_NVP(inSize);
   ar & BOOST_SERIALIZATION_NVP(outSize);
   ar & BOOST_SERIALIZATION_NVP(rho);
-
-  ar & BOOST_SERIALIZATION_NVP(input2GateModule);
-  ar & BOOST_SERIALIZATION_NVP(output2GateModule);
-  ar & BOOST_SERIALIZATION_NVP(inputGateModule);
-  ar & BOOST_SERIALIZATION_NVP(hiddenStateModule);
-  ar & BOOST_SERIALIZATION_NVP(forgetGateModule);
-  ar & BOOST_SERIALIZATION_NVP(outputGateModule);
-  ar & BOOST_SERIALIZATION_NVP(cellModule);
-  ar & BOOST_SERIALIZATION_NVP(cellActivationModule);
+  ar & BOOST_SERIALIZATION_NVP(bpttSteps);
+  ar & BOOST_SERIALIZATION_NVP(batchSize);
+  ar & BOOST_SERIALIZATION_NVP(batchStep);
+  ar & BOOST_SERIALIZATION_NVP(forwardStep);
+  ar & BOOST_SERIALIZATION_NVP(backwardStep);
+  ar & BOOST_SERIALIZATION_NVP(gradientStep);
+  ar & BOOST_SERIALIZATION_NVP(gradientStepIdx);
+  ar & BOOST_SERIALIZATION_NVP(cell);
+  ar & BOOST_SERIALIZATION_NVP(inputGateActivation);
+  ar & BOOST_SERIALIZATION_NVP(forgetGateActivation);
+  ar & BOOST_SERIALIZATION_NVP(outputGateActivation);
+  ar & BOOST_SERIALIZATION_NVP(hiddenLayerActivation);
+  ar & BOOST_SERIALIZATION_NVP(cellActivation);
+  ar & BOOST_SERIALIZATION_NVP(prevError);
+  ar & BOOST_SERIALIZATION_NVP(outParameter);
 }
 
 } // namespace ann
diff --git a/src/mlpack/methods/ann/rnn.hpp b/src/mlpack/methods/ann/rnn.hpp
index 184528502a..2616c9ba2c 100644
--- a/src/mlpack/methods/ann/rnn.hpp
+++ b/src/mlpack/methods/ann/rnn.hpp
@@ -221,6 +221,11 @@ class RNN
   //! Modify the maximum length of backpropagation through time.
   size_t& Rho() { return rho; }
 
+  /**
+   * Reset the module infomration (weights/parameters).
+   */
+  void ResetParameters();
+
   //! Serialize the model.
   template<typename Archive>
   void Serialize(Archive& ar, const unsigned int /* version */);
@@ -260,11 +265,6 @@ class RNN
    */
   void SinglePredict(const arma::mat& predictors, arma::mat& results);
 
-  /**
-   * Reset the module infomration (weights/parameters).
-   */
-  void ResetParameters();
-
   /**
    * Reset the module status by setting the current deterministic parameter
    * for all modules that implement the Deterministic function.
@@ -279,9 +279,6 @@ class RNN
   //! Number of steps to backpropagate through time (BPTT).
   size_t rho;
 
-  //! Number of steps to backpropagate through time (BPTT) at the previous step.
-  size_t prevRho;
-
   //! Instantiated outputlayer used to evaluate the network.
   OutputLayerType outputLayer;
 
diff --git a/src/mlpack/methods/ann/rnn_impl.hpp b/src/mlpack/methods/ann/rnn_impl.hpp
index 5281ed846e..a443fe33d7 100644
--- a/src/mlpack/methods/ann/rnn_impl.hpp
+++ b/src/mlpack/methods/ann/rnn_impl.hpp
@@ -37,7 +37,6 @@ RNN<OutputLayerType, InitializationRuleType>::RNN(
     OutputLayerType outputLayer,
     InitializationRuleType initializeRule) :
     rho(rho),
-    prevRho(0),
     outputLayer(std::move(outputLayer)),
     initializeRule(std::move(initializeRule)),
     inputSize(0),
@@ -60,7 +59,6 @@ RNN<OutputLayerType, InitializationRuleType>::RNN(
     OutputLayerType outputLayer,
     InitializationRuleType initializeRule) :
     rho(rho),
-    prevRho(0),
     outputLayer(std::move(outputLayer)),
     initializeRule(std::move(initializeRule)),
     inputSize(0),
@@ -121,7 +119,7 @@ void RNN<OutputLayerType, InitializationRuleType>::ResetCells()
 {
   for (size_t i = 1; i < network.size(); ++i)
   {
-    boost::apply_visitor(ResetCellVisitor(), network[i]);
+    boost::apply_visitor(ResetCellVisitor(rho), network[i]);
   }
 }
 
diff --git a/src/mlpack/methods/ann/visitor/reset_cell_visitor.hpp b/src/mlpack/methods/ann/visitor/reset_cell_visitor.hpp
index f2ccf12dd0..d1fdee57f4 100644
--- a/src/mlpack/methods/ann/visitor/reset_cell_visitor.hpp
+++ b/src/mlpack/methods/ann/visitor/reset_cell_visitor.hpp
@@ -26,23 +26,28 @@ namespace ann {
 class ResetCellVisitor : public boost::static_visitor<void>
 {
  public:
+  //! Reset the cell using the given size.
+  ResetCellVisitor(const size_t size);
+
   //! Execute the ResetCell() function.
   template<typename LayerType>
   void operator()(LayerType* layer) const;
 
  private:
-  //! Execute the ResetCell() function for a module which implements 
+  size_t size;
+
+  //! Execute the ResetCell() function for a module which implements
   //! the ResetCell() function.
   template<typename T>
   typename std::enable_if<
-      HasResetCellCheck<T, void(T::*)()>::value, void>::type
+      HasResetCellCheck<T, void(T::*)(const size_t)>::value, void>::type
   ResetCell(T* layer) const;
 
   //! Do not execute the Reset() function for a module which doesn't implement
   // the Reset() or Model() function.
   template<typename T>
   typename std::enable_if<
-      !HasResetCellCheck<T, void(T::*)()>::value, void>::type
+      !HasResetCellCheck<T, void(T::*)(const size_t)>::value, void>::type
   ResetCell(T* layer) const;
 };
 
diff --git a/src/mlpack/methods/ann/visitor/reset_cell_visitor_impl.hpp b/src/mlpack/methods/ann/visitor/reset_cell_visitor_impl.hpp
index 3d7e063c7d..659b9b5c9d 100644
--- a/src/mlpack/methods/ann/visitor/reset_cell_visitor_impl.hpp
+++ b/src/mlpack/methods/ann/visitor/reset_cell_visitor_impl.hpp
@@ -18,6 +18,12 @@
 namespace mlpack {
 namespace ann {
 
+//! ResetVisitor visitor class.
+inline ResetCellVisitor::ResetCellVisitor(const size_t size) : size(size)
+{
+  /* Nothing to do here. */
+}
+
 //! ResetVisitor visitor class.
 template<typename LayerType>
 inline void ResetCellVisitor::operator()(LayerType* layer) const
@@ -27,15 +33,15 @@ inline void ResetCellVisitor::operator()(LayerType* layer) const
 
 template<typename T>
 inline typename std::enable_if<
-    HasResetCellCheck<T, void(T::*)()>::value, void>::type
+    HasResetCellCheck<T, void(T::*)(const size_t)>::value, void>::type
 ResetCellVisitor::ResetCell(T* layer) const
 {
-  layer->ResetCell();
+  layer->ResetCell(size);
 }
 
 template<typename T>
 inline typename std::enable_if<
-    !HasResetCellCheck<T, void(T::*)()>::value, void>::type
+    !HasResetCellCheck<T, void(T::*)(const size_t)>::value, void>::type
 ResetCellVisitor::ResetCell(T* /* layer */) const
 {
   /* Nothing to do here. */
diff --git a/src/mlpack/tests/ann_layer_test.cpp b/src/mlpack/tests/ann_layer_test.cpp
index 171b958c2f..9c803a3416 100644
--- a/src/mlpack/tests/ann_layer_test.cpp
+++ b/src/mlpack/tests/ann_layer_test.cpp
@@ -152,7 +152,7 @@ double CheckGradient(FunctionType& function, const double eps = 1e-7)
 
   estGradient = arma::zeros(orgGradient.n_rows, orgGradient.n_cols);
 
-  // compute numeric approximations to gradient.
+  // Compute numeric approximations to gradient.
   for (size_t i = 0; i < orgGradient.n_elem; ++i)
   {
     double tmp = function.Parameters()(i);
@@ -736,6 +736,43 @@ BOOST_AUTO_TEST_CASE(SimpleAddMergeLayerTest)
   }
 }
 
+/**
+ * Test the LSTM layer with a user defined rho parameter and without.
+ */
+BOOST_AUTO_TEST_CASE(LSTMRrhoTest)
+{
+  const size_t rho = 5;
+  arma::mat input = arma::randu(5, 1);
+  arma::mat target = arma::mat("1; 1; 1; 1; 1");
+  RandomInitialization init(0.5, 0.5);
+
+  // Create model with user defined rho parameter.
+  RNN<NegativeLogLikelihood<>, RandomInitialization> modelA(
+      input, target, rho, false, NegativeLogLikelihood<>(), init);
+  modelA.Add<IdentityLayer<> >();
+  modelA.Add<Linear<> >(1, 10);
+
+  // Use LSTM layer with rho.
+  modelA.Add<LSTM<> >(10, 3, rho);
+  modelA.Add<LogSoftMax<> >();
+
+  // Create model without user defined rho parameter.
+  RNN<NegativeLogLikelihood<> > modelB(
+      input, target, rho, false, NegativeLogLikelihood<>(), init);
+  modelB.Add<IdentityLayer<> >();
+  modelB.Add<Linear<> >(1, 10);
+
+  // Use LSTM layer with rho = MAXSIZE.
+  modelB.Add<LSTM<> >(10, 3);
+  modelB.Add<LogSoftMax<> >();
+
+  optimization::StandardSGD opt(0.1, 1, 5, -100, false);
+  modelA.Train(input, target, opt);
+  modelB.Train(input, target, opt);
+
+  CheckMatrices(modelB.Parameters(), modelA.Parameters());
+}
+
 /**
  * LSTM layer numerically gradient test.
  */
@@ -779,6 +816,89 @@ BOOST_AUTO_TEST_CASE(GradientLSTMLayerTest)
   BOOST_REQUIRE_LE(CheckGradient(function), 1e-4);
 }
 
+/**
+ * Test the FastLSTM layer with a user defined rho parameter and without.
+ */
+BOOST_AUTO_TEST_CASE(FastLSTMRrhoTest)
+{
+  const size_t rho = 5;
+  arma::mat input = arma::randu(5, 1);
+  arma::mat target = arma::mat("1; 1; 1; 1; 1");
+  RandomInitialization init(0.5, 0.5);
+
+  // Create model with user defined rho parameter.
+  RNN<NegativeLogLikelihood<>, RandomInitialization> modelA(
+      input, target, rho, false, NegativeLogLikelihood<>(), init);
+  modelA.Add<IdentityLayer<> >();
+  modelA.Add<Linear<> >(1, 10);
+
+  // Use FastLSTM layer with rho.
+  modelA.Add<FastLSTM<> >(10, 3, rho);
+  modelA.Add<LogSoftMax<> >();
+
+  // Create model without user defined rho parameter.
+  RNN<NegativeLogLikelihood<> > modelB(
+      input, target, rho, false, NegativeLogLikelihood<>(), init);
+  modelB.Add<IdentityLayer<> >();
+  modelB.Add<Linear<> >(1, 10);
+
+  // Use FastLSTM layer with rho = MAXSIZE.
+  modelB.Add<FastLSTM<> >(10, 3);
+  modelB.Add<LogSoftMax<> >();
+
+  optimization::StandardSGD opt(0.1, 1, 5, -100, false);
+  modelA.Train(input, target, opt);
+  modelB.Train(input, target, opt);
+
+  CheckMatrices(modelB.Parameters(), modelA.Parameters());
+}
+
+/**
+ * FastLSTM layer numerically gradient test.
+ */
+BOOST_AUTO_TEST_CASE(GradientFastLSTMLayerTest)
+{
+  // Fast LSTM function gradient instantiation.
+  struct GradientFunction
+  {
+    GradientFunction()
+    {
+      input = arma::randu(5, 1);
+      target = arma::mat("1; 1; 1; 1; 1");
+      const size_t rho = 5;
+
+      model = new RNN<NegativeLogLikelihood<> >(input, target, rho);
+      model->Add<IdentityLayer<> >();
+      model->Add<Linear<> >(1, 10);
+      model->Add<FastLSTM<> >(10, 3, rho);
+      model->Add<LogSoftMax<> >();
+    }
+
+    ~GradientFunction()
+    {
+      delete model;
+    }
+
+    double Gradient(arma::mat& gradient) const
+    {
+      arma::mat output;
+      double error = model->Evaluate(model->Parameters(), 0, 1);
+      model->Gradient(model->Parameters(), 0, gradient, 1);
+      return error;
+    }
+
+    arma::mat& Parameters() { return model->Parameters(); }
+
+    RNN<NegativeLogLikelihood<> >* model;
+    arma::mat input, target;
+  } function;
+
+  // The threshold should be << 0.1 but since the Fast LSTM layer uses an
+  // approximation of the sigmoid function the estimated gradient is not
+  // correct.
+  BOOST_REQUIRE_LE(CheckGradient(function), 0.2);
+}
+
 /**
  * Check if the gradients computed by GRU cell are close enough to the
  * approximation of the gradients.
@@ -1103,5 +1223,4 @@ BOOST_AUTO_TEST_CASE(SimpleMeanSquaredErrorLayerTest)
   BOOST_REQUIRE_EQUAL(output.n_elem, 1);
 }
 
-
 BOOST_AUTO_TEST_SUITE_END();
diff --git a/src/mlpack/tests/recurrent_network_test.cpp b/src/mlpack/tests/recurrent_network_test.cpp
index 97ade2e60d..38857ef9ef 100644
--- a/src/mlpack/tests/recurrent_network_test.cpp
+++ b/src/mlpack/tests/recurrent_network_test.cpp
@@ -41,8 +41,8 @@ void GenerateNoisySines(arma::mat& data,
                         const size_t sequences,
                         const double noise = 0.3)
 {
-  arma::colvec x =  arma::linspace<arma::Col<double>>(0,
-      points - 1, points) / points * 20.0;
+  arma::colvec x =  arma::linspace<arma::colvec>(0, points - 1, points) /
+      points * 20.0;
   arma::colvec y1 = arma::sin(x + arma::as_scalar(arma::randu(1)) * 3.0);
   arma::colvec y2 = arma::sin(x / 2.0 + arma::as_scalar(arma::randu(1)) * 3.0);
 
@@ -363,11 +363,11 @@ void GenerateNextRecursiveReber(const arma::Mat<char>& transitions,
  * Train the specified network and the construct a Reber grammar dataset.
  */
 template<typename RecurrentLayerType>
-void ReberGrammarTestNetwork(size_t hiddenSize = 4,
-                             bool recursive = false,
-                             size_t averageRecursion = 3,
-                             size_t maxRecursion = 5
-                             )
+void ReberGrammarTestNetwork(const size_t hiddenSize = 4,
+                             const bool recursive = false,
+                             const size_t averageRecursion = 3,
+                             const size_t maxRecursion = 5,
+                             const size_t iterations = 10)
 {
   // Reber state transition matrix. (The last two columns are the indices to the
   // next path).
@@ -379,8 +379,8 @@ void ReberGrammarTestNetwork(size_t hiddenSize = 4,
               << 'P' << 'V' << '3' << '5' << arma::endr
               << 'E' << 'E' << '0' << '0' << arma::endr;
 
-  const size_t trainReberGrammarCount = 800;
-  const size_t testReberGrammarCount = 400;
+  const size_t trainReberGrammarCount = 700;
+  const size_t testReberGrammarCount = 250;
 
   std::string trainReber, testReber;
   arma::field<arma::mat> trainInput(1, trainReberGrammarCount);
@@ -450,16 +450,14 @@ void ReberGrammarTestNetwork(size_t hiddenSize = 4,
     const size_t inputSize = 7;
 
     RNN<MeanSquaredError<> > model(5);
-
     model.Add<Linear<> >(inputSize, hiddenSize);
     model.Add<RecurrentLayerType>(hiddenSize, hiddenSize);
     model.Add<Linear<> >(hiddenSize, outputSize);
     model.Add<SigmoidLayer<> >();
-
-    StandardSGD opt(0.01, 50, 2, -50000);
+    MomentumSGD opt(0.06, 50, 2, -50000);
 
     arma::mat inputTemp, labelsTemp;
-    for (size_t i = 0; i < (10 + offset); i++)
+    for (size_t i = 0; i < (iterations + offset); i++)
     {
       for (size_t j = 0; j < trainReberGrammarCount; j++)
       {
@@ -468,6 +466,7 @@ void ReberGrammarTestNetwork(size_t hiddenSize = 4,
 
         model.Rho() = inputTemp.n_elem / inputSize;
         model.Train(inputTemp, labelsTemp, opt);
+        opt.ResetPolicy() = false;
       }
     }
 
@@ -518,7 +517,7 @@ void ReberGrammarTestNetwork(size_t hiddenSize = 4,
     }
 
     error /= testReberGrammarCount;
-    if (error <= 0.2)
+    if (error <= 0.3)
     {
       ++successes;
       break;
@@ -531,27 +530,19 @@ void ReberGrammarTestNetwork(size_t hiddenSize = 4,
 }
 
 /**
- * Train the specified networks on a Reber grammar dataset.
+ * Train the specified networks on an embedded Reber grammar dataset.
  */
 BOOST_AUTO_TEST_CASE(LSTMReberGrammarTest)
 {
-  ReberGrammarTestNetwork<LSTM<>>(6, false);
+  ReberGrammarTestNetwork<LSTM<> >(10, false);
 }
 
 /**
  * Train the specified networks on an embedded Reber grammar dataset.
  */
-BOOST_AUTO_TEST_CASE(LSTMRecursiveReberGrammarTest)
+BOOST_AUTO_TEST_CASE(FastLSTMReberGrammarTest)
 {
-  ReberGrammarTestNetwork<LSTM<>>(22, true);
-}
-
-/**
- * Train the specified networks on a Reber grammar dataset.
- */
-BOOST_AUTO_TEST_CASE(GRUReberGrammarTest)
-{
-  ReberGrammarTestNetwork<GRU<>>(6, false);
+  ReberGrammarTestNetwork<FastLSTM<> >(8, false);
 }
 
 /**
@@ -559,7 +550,7 @@ BOOST_AUTO_TEST_CASE(GRUReberGrammarTest)
  */
 BOOST_AUTO_TEST_CASE(GRURecursiveReberGrammarTest)
 {
-  ReberGrammarTestNetwork<GRU<>>(20, true);
+  ReberGrammarTestNetwork<GRU<> >(16, true);
 }
 
 /*
@@ -595,8 +586,7 @@ void GenerateDistractedSequence(arma::mat& input, arma::mat& output)
   input = arma::zeros<arma::mat>(10, 10);
   output = arma::zeros<arma::mat>(3, 10);
 
-  arma::Col<size_t> index = arma::shuffle(arma::linspace<arma::Col<size_t> >(
-      0, 7, 8));
+  arma::uvec index = arma::shuffle(arma::linspace<arma::uvec>(0, 7, 8));
 
   // Set the target in the input sequence and the corresponding targets in the
   // output sequence by following the correct order.
@@ -623,10 +613,11 @@ void GenerateDistractedSequence(arma::mat& input, arma::mat& output)
  * dataset.
  */
 template<typename RecurrentLayerType>
-void DistractedSequenceRecallTestNetwork()
+void DistractedSequenceRecallTestNetwork(
+    const size_t cellSize, const size_t hiddenSize)
 {
-  const size_t trainDistractedSequenceCount = 800;
-  const size_t testDistractedSequenceCount = 400;
+  const size_t trainDistractedSequenceCount = 600;
+  const size_t testDistractedSequenceCount = 300;
 
   arma::field<arma::mat> trainInput(1, trainDistractedSequenceCount);
   arma::field<arma::mat> trainLabels(1, trainDistractedSequenceCount);
@@ -646,13 +637,13 @@ void DistractedSequenceRecallTestNetwork()
    * output units. The hidden layer is connected to itself. The network
    * structure looks like:
    *
-   *  Input         Hidden        Output
-   * Layer(10)  Layer(layerSize)   Layer(3)
-   * +-----+       +-----+       +-----+
-   * |     |       |     |       |     |
-   * |     +------>|     +------>|     |
-   * |     |    ..>|     |       |     |
-   * +-----+    .  +--+--+       +-----+
+   *  Input        Recurrent      Hidden       Output
+   * Layer(10)  Layer(cellSize)   Layer(3)     Layer(3)
+   * +-----+       +-----+       +-----+       +-----+
+   * |     |       |     |       |     |       |     |
+   * |     +------>|     +------>|     |------>|     |
+   * |     |    ..>|     |       |     |       |     |
+   * +-----+    .  +--+--+       +-----+       +-----+
    *            .     .
    *            .     .
    *            .......
@@ -671,15 +662,17 @@ void DistractedSequenceRecallTestNetwork()
   {
     RNN<MeanSquaredError<> > model(rho);
     model.Add<IdentityLayer<> >();
-    model.Add<Linear<> >(inputSize, 14);
-    model.Add<RecurrentLayerType>(14, 7);
-    model.Add<Linear<> >(7, outputSize);
+    model.Add<Linear<> >(inputSize, cellSize);
+    model.Add<RecurrentLayerType>(cellSize, hiddenSize);
+    model.Add<Linear<> >(hiddenSize, outputSize);
     model.Add<SigmoidLayer<> >();
 
     StandardSGD opt(0.1, 50, 2, -50000);
 
+    // We increase the number of iterations (training) if the first run didn't
+    // pass.
     arma::mat inputTemp, labelsTemp;
-    for (size_t i = 0; i < (10 + offset); i++)
+    for (size_t i = 0; i < (9 + offset); i++)
     {
       for (size_t j = 0; j < trainDistractedSequenceCount; j++)
       {
@@ -707,7 +700,6 @@ void DistractedSequenceRecallTestNetwork()
     }
 
     error /= testDistractedSequenceCount;
-
     // Can we reproduce the results from the paper. They provide an 95% accuracy
     // on a test set of 1000 randomly selected sequences.
     // Ensure that this is within tolerance, which is at least as good as the
@@ -730,7 +722,16 @@ void DistractedSequenceRecallTestNetwork()
  */
 BOOST_AUTO_TEST_CASE(LSTMDistractedSequenceRecallTest)
 {
-  DistractedSequenceRecallTestNetwork<LSTM<>>();
+  DistractedSequenceRecallTestNetwork<LSTM<> >(4, 8);
+}
+
+/**
+ * Train the specified networks on the Derek D. Monner's distracted sequence
+ * recall task.
+ */
+BOOST_AUTO_TEST_CASE(FastLSTMDistractedSequenceRecallTest)
+{
+  DistractedSequenceRecallTestNetwork<FastLSTM<> >(4, 8);
 }
 
 /**
@@ -739,7 +740,7 @@ BOOST_AUTO_TEST_CASE(LSTMDistractedSequenceRecallTest)
  */
 BOOST_AUTO_TEST_CASE(GRUDistractedSequenceRecallTest)
 {
-  DistractedSequenceRecallTestNetwork<GRU<>>();
+  DistractedSequenceRecallTestNetwork<GRU<> >(4, 8);
 }
 
 /**