Answers.Body: <p>This is a guess, as I am not reading all that code!</p>
  
  <blockquote>
    <p>This was previously working when the error I fed back through backprop was just the difference between the correct result and the prediction. But I've since been told that I should be propagating back the Loss Functions error, which I then implemented as Binary Cross Entropy.</p>
  </blockquote>
  
  <p>You may have been right before, and the advice you received wrong. Sort of.</p>
  
  <p>Backpropagation does not feed back the error value directly. It works exclusively with <em>gradients</em> of the error. So you typically start with a gradient value based on the objective function. </p>
  
  <p>If you backprop each function in the network a single step at a time, then you would start with <span class="math-container">$\nabla_{\hat{y}}\mathcal{L}$</span>, the gradient of the loss function with respect to the estimated values (i.e. output of the NN). Then from that you could calculate <span class="math-container">$\nabla_{z}\mathcal{L}$</span>, which is the gradient of the loss function with respect to the pre-activation values of the output layer.</p>
  
  <p>However, it is possible to combine steps analytically, and really common to start with the first gradient being <span class="math-container">$\nabla_{z}\mathcal{L}$</span>. That's because these pre-activation values are what you use to then calculate <span class="math-container">$\nabla_{W}\mathcal{L}$</span> and <span class="math-container">$\nabla_{b}\mathcal{L}$</span> - the gradients with respect to the layer's weights and biases, plus also to calculate <span class="math-container">$\nabla_{z'}\mathcal{L}$</span> - the gradients with respect to pre-activation values of the previous layer.</p>
  
  <p>Either way, once you have <span class="math-container">$\nabla_{z}\mathcal{L}$</span> for the output layer you can just run back though each layer in turn, repeating the same calcuation steps again and again. It is this repetition over each layer that looks like classic backprop in code.</p>
  
  <p>Assuming your neural network outputs a single value, the probablility of survival in the case of this dataset:</p>
  
  <ul>
  <li><p>Your objective/loss function for binary classification <em>should</em> be Binary Cross Entropy.</p></li>
  <li><p>Your output layer should have a sigmoid activation function.</p></li>
  </ul>
  
  <p>If you take the gradient of the loss function and backpropagate it throught that output layer's activation function, you end up with a starting delta <span class="math-container">$\nabla_{z} \mathcal{L}$</span> i.e. the gradient of the loss function with respect to the pre-activation value of the first layer. This is a nice place to start the recursive backprop code.</p>
  
  <p>The value of this gradient is mathematically, per item:</p>
  
  <p><span class="math-container">$$\nabla_{z} \mathcal{L} = \hat{y}_i - y_i$$</span> </p>
  
  <p>i.e. the difference the neural network estimate and the ground truth value. It takes this value because complex terms in the gradient of the loss function and the sigmoid activation function cancel out exactly. This is one reason why you often see sigmoid activation used with binary cross entropy - it is very convenient that the combination simplifies the gradient like this. You need to backpropagate each data point separately, and average the gradients for each batch/minibatch (actually you don't need to take this mean value, but doing so means that you don't need to adjust your learning rate as much to account for the number of items being processed per update step).</p>
  
  <p>It appears to of worked for you.</p>
  
  <p>Whether or not you were accidentally correct depends on whether you were getting an estimate for <span class="math-container">$\nabla_{z}\mathcal{L}$</span> to start your backprop routine, or <span class="math-container">$\nabla_{\hat{y}}\mathcal{L}$</span>. I cannot find what that assumption is in your code, but the fact that it worked before and does not now suggests that you may of been correct, even though you may not of fully understood the maths.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-27T20:52:20.037Z
Answers.Text: Neil, I'm disinclined to close this question since time went into this answer.  Do you think would could reform the question to be more suitable.  (The OP themself voted to close and tried to delete, but was blocked once there was an answer, so I think they'd be ok with editing the question.)
Answers.UserDisplayName: DukeZhou
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-10-19T10:16:15.853Z
Answers.LastActivityDate: 2019-10-19T15:45:29.890Z
Answers.LastEditDate: 2019-10-19T15:45:29.890Z
Answers.UserDisplayName: Neil Slater
Body: <p>My full code is as follows. I have tried to whittle it down to just the code that matters, but the problem I have is that i'm not sure what part of my network code is producing the problem. I've removed my code that loads and sifts through the CSV data because then my code would be too long.</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;array&gt;
  #include &lt;random&gt;
  #include &lt;chrono&gt;
  #include &lt;iomanip&gt;
  #include &lt;fstream&gt;
  #include &lt;algorithm&gt;
  #include &lt;iomanip&gt;
  #include &lt;variant&gt;
  #include &lt;unordered_set&gt;
  
  
  typedef std::variant&lt;std::string,std::uint_fast16_t,bool,float&gt; CSVType;
  
  /* ... functions to load CSV data ... */
  
  typedef float DataType;
  typedef DataType (*ActivationFuncPtr)(const DataType&amp;);
  
  DataType step(const DataType&amp; x, const DataType&amp; threshold)
  {
      return x &gt;= threshold ? 1 : 0;
  }
  
  DataType step0(const DataType&amp; x)
  {
      return step(x,0);
  }
  
  DataType step05(const DataType&amp; x)
  {
      return step(x,0.5);
  }
  
  DataType sigmoid(const DataType&amp; x)
  {
      return DataType(1) / (DataType(1) + std::exp(-x));
  }
  
  DataType sigmoid_derivative(const DataType&amp; x)
  {
      return x * (DataType(1) - x);
  }
  
  template&lt;std::size_t NumInputs&gt;
  class Neuron
  {
  public:
  
      Neuron()
      {
          RandomiseWeights();
      }
  
      void RandomiseWeights()
      {
          std::generate(m_weights.begin(),m_weights.end(),[&amp;]()
          {
              return m_xavierNormalDis(m_mt);
          });
          m_biasWeight = 0;
  
          for(std::size_t i = 0; i &lt; NumInputs+1; ++i)
              m_previousWeightUpdates[i] = 0;
      }
  
      DataType FeedForward(const std::array&lt;DataType,NumInputs&gt; inputValues)
      {
          DataType output = m_biasWeight;
          for(std::size_t i = 0; i &lt; inputValues.size(); ++i)
              output += inputValues[i] * m_weights[i];
  
          m_inputValues = inputValues;
  
          return output;
      }
  
      std::array&lt;DataType,NumInputs&gt; Backpropagate(const DataType&amp; error)
      {
          std::array&lt;DataType,NumInputs&gt; netInputOverWeight;
          for(std::size_t i = 0; i &lt; NumInputs; ++i)
          {
              netInputOverWeight[i] = m_inputValues[i];
          }
  
          DataType netInputOverBias = DataType(1);
  
          std::array&lt;DataType,NumInputs&gt; errorOverWeight;
          for(std::size_t i = 0; i &lt; NumInputs; ++i)
          {
              errorOverWeight[i] = error * netInputOverWeight[i];
          }
  
          DataType errorOverBias = error * netInputOverBias;
  
          for(std::size_t i = 0; i &lt; NumInputs; ++i)
          {
              m_outstandingWeightAdjustments[i] = errorOverWeight[i];
          }
          m_outstandingWeightAdjustments[NumInputs] = errorOverBias;
  
          DataType errorOverNetInput = error;
  
          std::array&lt;DataType,NumInputs&gt; errorWeights;
          for(std::size_t i = 0; i &lt; NumInputs; ++i)
          {
              errorWeights[i] = errorOverNetInput * m_weights[i];
          }
  
          return errorWeights;
      }
  
      void AdjustWeights(const DataType&amp; learningRate, const DataType&amp; momentum)
      {
          for(std::size_t i = 0; i &lt; NumInputs; ++i)
          {
              DataType adjustment = learningRate * m_outstandingWeightAdjustments[i] + momentum * m_previousWeightUpdates[i];
              m_weights[i] = m_weights[i] - adjustment;
              m_previousWeightUpdates[i] = adjustment;
          }
          DataType adjustment = learningRate * m_outstandingWeightAdjustments[NumInputs] + momentum * m_previousWeightUpdates[NumInputs];
          m_biasWeight = m_biasWeight - adjustment;
          m_previousWeightUpdates[NumInputs] = adjustment;
      }
  
      const std::array&lt;DataType,NumInputs&gt;&amp; GetWeights() const
      {
          return m_weights;
      }
  
      const DataType&amp; GetBiasWeight() const
      {
          return m_biasWeight;
      }
  
  protected:
  
      static std::mt19937 m_mt;
      static std::uniform_real_distribution&lt;DataType&gt; m_uniformDisRandom;
      static std::uniform_real_distribution&lt;DataType&gt; m_xavierUniformDis;
      static std::normal_distribution&lt;DataType&gt; m_xavierNormalDis;
  
      std::array&lt;DataType,NumInputs&gt; m_weights;
      DataType m_biasWeight;
  
      std::array&lt;DataType,NumInputs+1&gt; m_previousWeightUpdates;
      std::array&lt;DataType,NumInputs+1&gt; m_outstandingWeightAdjustments;
  
      std::array&lt;DataType,NumInputs&gt; m_inputValues;
  };
  
  template&lt;std::size_t NumInputs&gt;
  std::mt19937 Neuron&lt;NumInputs&gt;::m_mt(std::chrono::duration_cast&lt;std::chrono::milliseconds&gt;(std::chrono::system_clock::now().time_since_epoch()).count());
  
  template&lt;std::size_t NumInputs&gt;
  std::uniform_real_distribution&lt;DataType&gt; Neuron&lt;NumInputs&gt;::m_uniformDisRandom(-1,1);
  
  template&lt;std::size_t NumInputs&gt;
  std::uniform_real_distribution&lt;DataType&gt; Neuron&lt;NumInputs&gt;::m_xavierUniformDis(-std::sqrt(6.f / NumInputs+1),std::sqrt(6.f / NumInputs+1));
  
  template&lt;std::size_t NumInputs&gt;
  std::normal_distribution&lt;DataType&gt; Neuron&lt;NumInputs&gt;::m_xavierNormalDis(0,std::sqrt(2.f / NumInputs+1));
  
  template&lt;std::size_t NumNeurons&gt;
  class ActivationLayer
  {
  public:
  
      ActivationLayer()
      :
          m_outputs({})
      {}
  
      virtual std::array&lt;DataType,NumNeurons&gt; GetOutputs() const final
      {
          return m_outputs;
      }
  
      virtual void CompleteBackprop(const DataType&amp; learningRate, const DataType&amp; momentum) final
      {
      }
  
  protected:
      std::array&lt;DataType,NumNeurons&gt; m_outputs;
  };
  
  template&lt;std::size_t NumNeurons&gt;
  class SigmoidActivation : public ActivationLayer&lt;NumNeurons&gt;
  {
  public:
  
      virtual std::array&lt;DataType,NumNeurons&gt; FeedForward(const std::array&lt;DataType,NumNeurons&gt;&amp; inputValues)
      {
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
              ActivationLayer&lt;NumNeurons&gt;::m_outputs[i] = sigmoid(inputValues[i]);
          return ActivationLayer&lt;NumNeurons&gt;::m_outputs;
      }
  
      virtual std::array&lt;DataType,NumNeurons&gt; Backpropagate(const std::array&lt;DataType,NumNeurons&gt; errors)
      {
          std::array&lt;DataType,NumNeurons&gt; backpropErrors;
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
              backpropErrors[i] = errors[i] * sigmoid_derivative(ActivationLayer&lt;NumNeurons&gt;::m_outputs[i]);
          return backpropErrors;
      }
  };
  
  template&lt;std::size_t NumInputs, std::size_t NumNeurons&gt;
  class FullyConnectedLayer
  {
  public:
  
      FullyConnectedLayer()
      :
          m_neurons([=]()
          {
              std::array&lt;Neuron&lt;NumInputs&gt;,NumNeurons&gt; neurons;
              for(Neuron&lt;NumInputs&gt;&amp; n : neurons)
                  n = Neuron&lt;NumInputs&gt;();
              return neurons;
          }())
      {
      }
  
      virtual std::array&lt;DataType,NumNeurons&gt; FeedForward(const std::array&lt;DataType,NumInputs&gt;&amp; inputValues)
      {
          std::array&lt;DataType,NumNeurons&gt; outputValues;
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
              outputValues[i] = m_neurons[i].FeedForward(inputValues);
          return outputValues;
      }
  
      /** \brief Take a sum of errors for each node and produce the errors for each input node in the previous layer.
       *
       */
  
      virtual std::array&lt;DataType,NumInputs&gt;
      Backpropagate(const std::array&lt;DataType,NumNeurons&gt; errors)
      {
          std::array&lt;std::array&lt;DataType,NumInputs&gt;,NumNeurons&gt; errorValues;
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
          {
              errorValues[i] = m_neurons[i].Backpropagate(errors[i]);
          }
          std::array&lt;DataType,NumInputs&gt; returnErrors;
          std::fill(returnErrors.begin(),returnErrors.end(),0);
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
          {
              for(std::size_t j = 0; j &lt; NumInputs; ++j)
              {
                  returnErrors[j] += errorValues[i][j];
              }
          }
          return returnErrors;
      }
  
      virtual void CompleteBackprop(const DataType&amp; learningRate, const DataType&amp; momentum)
      {
          for(Neuron&lt;NumInputs&gt;&amp; n : m_neurons)
              n.AdjustWeights(learningRate, momentum);
      }
  
      const Neuron&lt;NumInputs&gt;&amp; operator[](const std::size_t&amp; index) const
      {
          return m_neurons[index];
      }
  
      std::array&lt;std::array&lt;DataType,NumInputs&gt;,NumNeurons&gt; GetWeights() const
      {
          std::array&lt;std::array&lt;DataType,NumInputs&gt;,NumNeurons&gt; weights;
          for(std::size_t i = 0; i &lt; NumNeurons; ++i)
          {
              weights[i] = m_neurons[i].GetWeights();
          }
          return weights;
      }
  
  protected:
      std::array&lt;Neuron&lt;NumInputs&gt;,NumNeurons&gt; m_neurons;
  };
  
  template&lt;std::size_t I = 0, typename FuncT, typename... Tp&gt;
  inline typename std::enable_if&lt;I == sizeof...(Tp)&gt;::type for_each(std::tuple&lt;Tp...&gt; &amp;, FuncT)
  {
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename... Tp&gt;
  inline typename std::enable_if&lt;I &lt; sizeof...(Tp)&gt;::type for_each(std::tuple&lt;Tp...&gt;&amp; t, FuncT f)
  {
      f(std::get&lt;I&gt;(t)); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each&lt;I + 1, FuncT, Tp...&gt;(t, f); // process the next element of the tuple with the new output
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename... Tp&gt;
  inline typename std::enable_if&lt;I == sizeof...(Tp)&gt;::type for_each(const std::tuple&lt;Tp...&gt; &amp;, FuncT)
  {
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename... Tp&gt;
  inline typename std::enable_if&lt;I &lt; sizeof...(Tp)&gt;::type for_each(const std::tuple&lt;Tp...&gt;&amp; t, FuncT f)
  {
      f(std::get&lt;I&gt;(t)); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each&lt;I + 1, FuncT, Tp...&gt;(t, f); // process the next element of the tuple with the new output
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename O, typename FinalOutput, typename... Tp&gt;
  inline typename std::enable_if&lt;I == sizeof...(Tp)&gt;::type for_each_get_final_output(std::tuple&lt;Tp...&gt; &amp;, FuncT, O o, FinalOutput&amp; finalOutput)
  {
      finalOutput = o;
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename O, typename FinalOutput, typename... Tp&gt;
  inline typename std::enable_if&lt;I &lt; sizeof...(Tp)&gt;::type for_each_get_final_output(std::tuple&lt;Tp...&gt;&amp; t, FuncT f, O o, FinalOutput&amp; finalOutput)
  {
      auto newO = f(std::get&lt;I&gt;(t),o); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each_get_final_output&lt;I + 1, FuncT, decltype(newO), FinalOutput, Tp...&gt;(t, f, newO, finalOutput); // process the next element of the tuple with the new output
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename O, typename... Tp&gt;
  inline typename std::enable_if&lt;I == 0&gt;::type for_each_reverse_impl(std::tuple&lt;Tp...&gt;&amp; t, FuncT f, O o)
  {
      f(std::get&lt;0&gt;(t),o);
  }
  
  template&lt;std::size_t I = 0, typename FuncT, typename O, typename... Tp&gt;
  inline typename std::enable_if&lt;(I &gt; 0)&gt;::type
  for_each_reverse_impl(std::tuple&lt;Tp...&gt;&amp; t, FuncT f, O o)
  {
      auto newO = f(std::get&lt;I&gt;(t),o); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each_reverse_impl&lt;I - 1, FuncT, decltype(newO), Tp...&gt;(t, f, newO); // process the next element of the tuple with the new output
  }
  
  template&lt;typename FuncT, typename O, typename... Tp&gt;
  inline void for_each_reverse(std::tuple&lt;Tp...&gt;&amp; t, FuncT f, O o)
  {
      for_each_reverse_impl&lt;sizeof...(Tp)-1, FuncT, O, Tp...&gt;(t, f, o);
  }
  
  enum class LOSS_FUNCTION : std::uint_fast8_t
  {
      MEAN_SQUARE_ERROR,
      CROSS_ENTROPY
  };
  
  class ValidationOptions
  {
  public:
      enum class METRIC : std::uint_fast8_t { NONE, ACCURACY, LOSS };
  
      ValidationOptions()
      :
          m_validationSplit(.3f),
          m_enableLoss(false),
          m_lossFunction(LOSS_FUNCTION::MEAN_SQUARE_ERROR),
          m_enableAccuracy(false),
          m_outputFilter([](const DataType&amp; x){ return x; }),
          m_earlyStoppingMetric(METRIC::NONE),
          m_earlyStoppingPatience(1.f),
          m_earlyStoppingDelta(1.f),
          m_earlyStoppingNumEpochsAverage(1)
      {}
  
      ValidationOptions&amp; Loss(const bool enable = true, LOSS_FUNCTION lossFunction = LOSS_FUNCTION::MEAN_SQUARE_ERROR)
      {
          m_enableLoss = enable;
          m_lossFunction = lossFunction;
          return *this;
      }
  
      ValidationOptions&amp; Split(const float dataSplitValidation)
      {
          m_validationSplit = dataSplitValidation;
          return *this;
      }
  
      ValidationOptions&amp; Accuracy(const bool enable = true, ActivationFuncPtr outputFilter = [](const DataType&amp; x){return x;})
      {
          m_enableAccuracy = enable;
          m_outputFilter = outputFilter;
          return *this;
      }
  
      ValidationOptions&amp; EarlyStop(const bool enable = true,
                                   const METRIC metric = METRIC::ACCURACY,
                                   const float patience = .1f,
                                   const DataType delta = .01,
                                   const std::size_t epochNumToAverage = 10)
       {
           if(enable == false)
              m_earlyStoppingMetric = METRIC::NONE;
           else
              m_earlyStoppingMetric = metric;
  
           m_earlyStoppingPatience = patience;
           m_earlyStoppingDelta = delta;
           m_earlyStoppingNumEpochsAverage = epochNumToAverage;
  
           return *this;
       }
  
       float GetValidationSplit() const { return m_validationSplit; }
       bool Loss() const { return m_enableLoss; }
       LOSS_FUNCTION GetLossFunction() const { return m_lossFunction; }
       bool Accuracy() const { return m_enableAccuracy; }
       ActivationFuncPtr GetOutputFilter() const { return m_outputFilter; }
       METRIC GetEarlyStoppingMetric() const { return m_earlyStoppingMetric; }
       float GetEarlyStoppingPatience() const { return m_earlyStoppingPatience; }
       DataType GetEarlyStoppingDelta() const { return m_earlyStoppingDelta; }
       std::size_t GetEarlyStoppingNumEpochsAvg() const { return m_earlyStoppingNumEpochsAverage; }
  
  protected:
  
      float m_validationSplit;        /**&lt; Percentage of the data set aside for validation */
  
      bool m_enableLoss;
      LOSS_FUNCTION m_lossFunction;   /**&lt; Loss function to use */
  
      bool m_enableAccuracy;
      ActivationFuncPtr m_outputFilter;   /**&lt; When measuring accuracy data is passed through this */
  
      METRIC m_earlyStoppingMetric;                   /**&lt; The metric used to stop early */
      float m_earlyStoppingPatience;                  /**&lt; Percentage of total epochs to wait before stopping early */
      DataType m_earlyStoppingDelta;                  /**&lt; The amount that the early stopping metric needs to change in a single step before stopping */
      std::size_t m_earlyStoppingNumEpochsAverage;    /**&lt; The number of epochs averaged over to smooth out the stopping metric */
  };
  
  template&lt;typename... Layers&gt;
  class NeuralNetwork
  {
  public:
  
      NeuralNetwork(Layers... layers)
      :
          m_layers(std::make_tuple(layers...))
      {
  
      }
  
      template&lt;std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumTrainingRows&gt;
      void Fit(const std::size_t&amp; numberEpochs,
               const std::size_t&amp; batchSize,
               DataType learningRate,
               const DataType&amp; momentum,
               std::array&lt;std::array&lt;DataType,NumFeatures&gt;,NumTrainingRows&gt;&amp; trainingData,
               std::array&lt;std::array&lt;DataType,NumOutputs&gt;,NumTrainingRows&gt;&amp; trainingOutput,
               const ValidationOptions validationOptions,
               const bool linearDecayLearningRate = true,
               std::ostream&amp; outputStream = std::cout)
      {
          std::size_t epochNumber = 0;
  
          // need to support more than just MSE to measure loss
          std::vector&lt;DataType&gt; lastEpochLoss(validationOptions.GetEarlyStoppingNumEpochsAvg(),0);
          DataType lastEpochLossAverage = std::numeric_limits&lt;DataType&gt;::max();
  
          std::vector&lt;DataType&gt; lastValidationAccuracys(validationOptions.GetEarlyStoppingNumEpochsAvg(),0);
          DataType lastValidationAccuraryAvg = 0;
  
          std::vector&lt;std::size_t&gt; randomIndices(NumTrainingRows,0);
          for(std::size_t i = 0; i &lt; NumTrainingRows; ++i)
              randomIndices[i] = i;
  
          std::random_shuffle(randomIndices.begin(),randomIndices.end());
          // take some percentage as validation split
          // we do this by taking the first percentage of already shuffled indices and removing them
          // from what is available
          std::size_t numValidationRecords = NumTrainingRows*validationOptions.GetValidationSplit();
          std::size_t numTrainingRecords = NumTrainingRows - numValidationRecords;
          std::vector&lt;std::size_t&gt; validationRecords(numValidationRecords);
          for(std::size_t i = 0; i &lt; numValidationRecords; ++i)
          {
              std::size_t index = randomIndices.back();
              randomIndices.pop_back();
              validationRecords[i] = index;
          }
  
          while(epochNumber &lt; numberEpochs)
          {
              // shuffle the indices so that they are pulled into each batch randomly each time
              std::random_shuffle(randomIndices.begin(),randomIndices.end());
  
              DataType epochLoss = 0;
  
              std::tuple&lt;Layers...&gt; backupLayers = m_layers;
  
              for(std::size_t batchNumber = 0; batchNumber &lt; std::ceil(numTrainingRecords / batchSize); ++batchNumber)
              {
                  std::array&lt;DataType,NumOutputs&gt; propagateError = {0};
  
                  std::size_t startIndex = batchNumber * batchSize;
                  std::size_t endIndex = startIndex + batchSize;
                  if(endIndex &gt; numTrainingRecords)
                      endIndex = numTrainingRecords;
  
                  DataType batchLoss = 0;
  
                  for(std::size_t index = startIndex; index &lt; endIndex; ++index)
                  {
                      std::size_t row = randomIndices[index];
                      const std::array&lt;DataType,NumFeatures&gt;&amp; dataRow = trainingData[row];
                      const std::array&lt;DataType,NumOutputs&gt;&amp; desiredOutputRow = trainingOutput[row];
  
                      // Feed the values through to the output layer
                      // use of "auto" is so this lambda can be used for all layers without
                      // me needing to do any fucking around
                      std::array&lt;DataType,NumOutputs&gt; finalOutput;
                      for_each_get_final_output(m_layers, [](auto&amp; layer, auto o)
                      {
                          return layer.FeedForward(o);
                      }, dataRow, finalOutput);
  
                      DataType totalError = 0;
                      for(std::size_t i = 0; i &lt; NumOutputs; ++i)
                      {
                          if(validationOptions.GetLossFunction() == LOSS_FUNCTION::MEAN_SQUARE_ERROR)
                              totalError += std::pow(desiredOutputRow[i] - finalOutput[i],2.0);
                          else if(validationOptions.GetLossFunction() == LOSS_FUNCTION::CROSS_ENTROPY)
                          {
                              if(NumOutputs == 1)
                              {
                                  // binary cross entropy
                                  totalError += (desiredOutputRow[i] * std::log(1e-15 + finalOutput[i]));
                              }
                              else
                              {
                                  // cross entropy
                              }
                          }
                      }
  
                      batchLoss += totalError;
                  }
  
                  batchLoss *= DataType(1) / (endIndex - startIndex);
  
                  for(std::size_t i = 0; i &lt; NumOutputs; ++i)
                      propagateError[i] = batchLoss;
  
                  // update after every batch
                  for_each_reverse(m_layers, [](auto&amp; layer, auto o)
                  {
                      auto errors = layer.Backpropagate(o);
                      return errors;
                  }, propagateError);
  
                  // once backprop is finished, we can adjust all the weights
                  for_each(m_layers, [&amp;](auto&amp; layer)
                  {
                      layer.CompleteBackprop(learningRate,momentum);
                  });
  
                  epochLoss += batchLoss;
              }
  
              epochLoss *= DataType(1) / numTrainingRecords;
  
              lastEpochLoss.erase(lastEpochLoss.begin());
              lastEpochLoss.push_back(epochLoss);
              DataType avgEpochLoss = 1.f * std::accumulate(lastEpochLoss.begin(),lastEpochLoss.end(),0.f) / (epochNumber &lt; validationOptions.GetEarlyStoppingNumEpochsAvg() ? epochNumber+1 : lastEpochLoss.size());
  
              if(validationOptions.GetEarlyStoppingMetric() == ValidationOptions::METRIC::LOSS
                 &amp;&amp; epochNumber &gt; numberEpochs * validationOptions.GetEarlyStoppingPatience()
                 &amp;&amp; avgEpochLoss &gt; lastEpochLossAverage + validationOptions.GetEarlyStoppingDelta())
              {
                  // the loss average has decreased, so we should go back to the previous run and exit
                  std::cout   &lt;&lt; "Early exit Loss Avg \n"
                              &lt;&lt; "Last Epoch: " &lt;&lt; lastEpochLossAverage &lt;&lt; "\n"
                              &lt;&lt; "This Epoch: " &lt;&lt; avgEpochLoss &lt;&lt; std::endl;
                  m_layers = backupLayers;
                  break;
              }
  
              lastEpochLossAverage = avgEpochLoss;
  
              // check for the error against the reserved validation set
              std::size_t numCorrect = 0;
              for(std::size_t row = 0; row &lt; validationRecords.size(); ++row)
              {
                  const std::array&lt;DataType,NumFeatures&gt;&amp; dataRow = trainingData[row];
                  const std::array&lt;DataType,NumOutputs&gt;&amp; desiredOutputRow = trainingOutput[row];
  
                  std::array&lt;DataType,NumOutputs&gt; finalOutput;
                  for_each_get_final_output(m_layers, [](auto&amp; layer, auto o)
                  {
                      return layer.FeedForward(o);
                  }, dataRow, finalOutput);
  
                  bool correct = true;
                  for(std::size_t i = 0; i &lt; NumOutputs; ++i)
                  {
                      if(validationOptions.GetOutputFilter()(finalOutput[i]) != desiredOutputRow[i])
                          correct = false;
                  }
                  if(correct)
                      ++numCorrect;
              }
  
              DataType validationAccuracy = DataType(numCorrect) / numValidationRecords;
  
              lastValidationAccuracys.erase(lastValidationAccuracys.begin());
              lastValidationAccuracys.push_back(validationAccuracy);
              DataType avgValidationAccuracy = std::accumulate(lastValidationAccuracys.begin(),lastValidationAccuracys.end(),0.f) / (epochNumber &lt; validationOptions.GetEarlyStoppingNumEpochsAvg() ? epochNumber+1 : lastValidationAccuracys.size());
  
              if(validationOptions.GetEarlyStoppingMetric() == ValidationOptions::METRIC::ACCURACY
                 &amp;&amp; epochNumber &gt; numberEpochs * validationOptions.GetEarlyStoppingPatience()
                 &amp;&amp; avgValidationAccuracy &lt; lastValidationAccuraryAvg - validationOptions.GetEarlyStoppingDelta())
              {
                  // the accuracy has decreased, so we should go back to the previous run and exit
                  std::cout   &lt;&lt; "Early exit validation accuracy \n"
                              &lt;&lt; "Last Epoch: " &lt;&lt; lastValidationAccuraryAvg &lt;&lt; "\n"
                              &lt;&lt; "This Epoch: " &lt;&lt; avgValidationAccuracy &lt;&lt; std::endl;
                  m_layers = backupLayers;
                  break;
              }
  
              lastValidationAccuraryAvg = avgValidationAccuracy;
  
              outputStream &lt;&lt; epochNumber &lt;&lt; "," &lt;&lt; epochLoss &lt;&lt; "," &lt;&lt; avgEpochLoss &lt;&lt; "," &lt;&lt; validationAccuracy &lt;&lt; "," &lt;&lt; avgValidationAccuracy &lt;&lt; std::endl;
  
              learningRate -= learningRate / (numberEpochs-epochNumber);
  
              ++epochNumber;
          }
      }
  
      template&lt;std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumEvaluationRows&gt;
      void Evaluate(std::array&lt;std::array&lt;DataType,NumFeatures&gt;,NumEvaluationRows&gt; inputData,
                    std::array&lt;std::array&lt;DataType,NumOutputs&gt;,NumEvaluationRows&gt; correctOutputs,
                    DataType&amp; loss,
                    DataType&amp; accuracy,
                    ActivationFuncPtr outputFilter = [](const DataType&amp; x){return x;})
      {
          loss = 0;
  
          std::size_t numCorrect = 0;
  
          for(std::size_t row = 0; row &lt; NumEvaluationRows; ++row)
          {
              const std::array&lt;DataType,NumFeatures&gt;&amp; dataRow = inputData[row];
              const std::array&lt;DataType,NumOutputs&gt;&amp; outputRow = correctOutputs[row];
  
              // Feed the values through to the output layer
  
              std::array&lt;DataType,NumOutputs&gt; finalOutput;
              for_each_get_final_output(m_layers, [](auto&amp; layer, auto o)
              {
                  layer.FeedForward(o);
                  return layer.GetOutputs();
              }, dataRow, finalOutput);
  
              DataType thisLoss = 0;
              for(std::size_t i = 0; i &lt; NumOutputs; ++i)
                  thisLoss += outputRow[i] - finalOutput[i];
              loss += thisLoss * thisLoss;
  
              bool correct = true;
              for(std::size_t i = 0; i &lt; NumOutputs; ++i)
              {
                  if(outputFilter(finalOutput[i]) != outputRow[i])
                      correct = false;
              }
              if(correct)
                  ++numCorrect;
          }
  
          loss *= DataType(1) / NumEvaluationRows;
          accuracy = DataType(numCorrect) / NumEvaluationRows;
      }
  
      template&lt;std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumRecords&gt;
      void Predict(std::array&lt;std::array&lt;DataType,NumFeatures&gt;,NumRecords&gt; inputData,
                   std::array&lt;std::array&lt;DataType,NumOutputs&gt;,NumRecords&gt;&amp; predictions,
                   ActivationFuncPtr outputFilter = [](const DataType&amp; x){return x;})
      {
          for(std::size_t row = 0; row &lt; NumRecords; ++row)
          {
              const std::array&lt;DataType,NumFeatures&gt;&amp; dataRow = inputData[row];
  
              // Feed the values through to the output layer
  
              std::array&lt;DataType,NumOutputs&gt; finalOutput;
              for_each_get_final_output(m_layers, [](auto&amp; layer, auto o)
              {
                  return layer.FeedForward(o);
              }, dataRow, finalOutput);
  
              for(std::size_t i = 0; i &lt; NumOutputs; ++i)
                  predictions[row][i] = outputFilter(finalOutput[i]);
          }
      }
  
  protected:
      std::tuple&lt;Layers...&gt; m_layers;
  };
  
  main()
  {
      std::vector&lt;std::vector&lt;CSVType&gt;&gt; trainingCSVData;
      /* load training CSV data */
  
      std::vector&lt;std::vector&lt;CSVType&gt;&gt; testCSVData;
      /* load test CSV data */
  
      std::cout &lt;&lt; std::fixed &lt;&lt; std::setprecision(80);
  
      std::ofstream file("error_out.csv", std::ios::out | std::ios::trunc);
      if(!file.is_open())
      {
          std::cout &lt;&lt; "couldn't open file" &lt;&lt; std::endl;
          return 0;
      }
  
      file &lt;&lt; std::fixed &lt;&lt; std::setprecision(80);
  
      /*
          Features
          1   pClass 1
          2   pClass 2
          3   pClass 3
          4   Sex female 1, male 0
          5   Age normalised between 0 and 1  age range 0 to 100
          6   Number siblings between 0 and 1   num range 0 to 8
          7   Number of parents / children        num range 0 to 9
          8   Ticket cost   between 0 and 1       num range 0 to 512.3292
          9   embarked S
          10  embarked Q
          11  embarked C
      */
  
      std::array&lt;std::array&lt;DataType,29&gt;,891&gt; inputData;
      std::array&lt;std::array&lt;DataType,1&gt;,891&gt; desiredOutputs;
  
      /* ... data that loads the titanic data into a series of features. Either class labels or normalised values (like age) */
  
      NeuralNetwork neuralNet{
          FullyConnectedLayer&lt;29,256&gt;(),
          SigmoidActivation&lt;256&gt;(),
          FullyConnectedLayer&lt;256,1&gt;(),
          SigmoidActivation&lt;1&gt;()
      };
  
      neuralNet.Fit(300,
                    1,
                    0.05,
                    0.25f,
                    inputData,
                    desiredOutputs,
                    ValidationOptions().Accuracy(true,step05).Loss(true,LOSS_FUNCTION::CROSS_ENTROPY).Split(0.3),
                    false,
                    file);
  
      file.close();
  
      return 0;
  }
  </code></pre>
  
  <p>The data used is from the titanic problem that you can <a href="https://www.kaggle.com/c/titanic/data" rel="nofollow noreferrer">download from Kaggle here</a>.</p>
  
  <p>The typical output file that's being generated is like this:</p>
  
  <pre><code>0,-4.91843843460083007812500000000000000000000000000000000000000000000000000000000000,-4.91843843460083007812500000000000000000000000000000000000000000000000000000000000,0.65168541669845581054687500000000000000000000000000000000000000000000000000000000,0.65168541669845581054687500000000000000000000000000000000000000000000000000000000
  1,-6.14257431030273437500000000000000000000000000000000000000000000000000000000000000,-6.14257431030273437500000000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  2,-6.43130302429199218750000000000000000000000000000000000000000000000000000000000000,-6.43130302429199218750000000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  3,-6.58864736557006835937500000000000000000000000000000000000000000000000000000000000,-6.58864736557006835937500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  4,-6.70884752273559570312500000000000000000000000000000000000000000000000000000000000,-6.70884752273559570312500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  5,-6.78206682205200195312500000000000000000000000000000000000000000000000000000000000,-6.78206682205200195312500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  6,-6.86832284927368164062500000000000000000000000000000000000000000000000000000000000,-6.86832284927368164062500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  7,-6.92110681533813476562500000000000000000000000000000000000000000000000000000000000,-6.92110681533813476562500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  8,-6.96584081649780273437500000000000000000000000000000000000000000000000000000000000,-6.96584081649780273437500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  9,-7.02414274215698242187500000000000000000000000000000000000000000000000000000000000,-7.02414274215698242187500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  10,-7.06421041488647460937500000000000000000000000000000000000000000000000000000000000,-7.06421041488647460937500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  </code></pre>
  
  <p>(shortened for space)</p>
  
  <p>This was previously working when the error I fed back through backprop was just the difference between the correct result and the prediction. <a href="https://www.reddit.com/r/learnmachinelearning/comments/dixyiy/i_think_ive_got_some_elements_of_neural_networks/" rel="nofollow noreferrer">But I've since been told that I should be propagating back the Loss Functions error</a>, which I then implemented as Binary Cross Entropy.</p>
  
  <p>So either:</p>
  
  <ol>
  <li>I should not feed back the error from the loss function</li>
  <li>I am calculating the loss incorrectly</li>
  <li>My back propagation code is wrong</li>
  <li>I've got something horrible happening in my validation loop</li>
  </ol>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-10-19T06:59:45.257Z
Comments.Text: Do a forward and backrpop by hand and compare to the values you get in your code. Speaking of which, that is a huge slab of code. I would say you would be better posting this to data science exchange as ai is generally for philosophical and theory, but this is such a large amount of code I don't think anyone at data science exchange would want to go over it. I think your best bet is to debug your code in segments, and see how that goes.
Comments.UserDisplayName: Recessive
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-10-19T10:00:07.713Z
Comments.Text: Recessive seems correct to me. I don't usually mind reading a fair amount of code - I've answered similar implementation questions before, and have implemented neural networks in C (not C++ as here). But the size of the challenge to find your bug(s) in the attached code looks like too much effort for me. I will though have one guess.
Comments.UserDisplayName: Neil Slater
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-10-27T20:47:02.387Z
Comments.Text: @NeomerArcana No worries!  *(The overwhelming majority of my questions across stack get close;)*
Comments.UserDisplayName: DukeZhou
Tags: <neural-networks><backpropagation>
Title: Why is my neural network giving me wildly incorrect error and not changing accuracy?
UserDisplayName: NeomerArcana
ThreadText: # Why is my neural network giving me wildly incorrect error and not changing accuracy?
  
  (Asked by: username_0 on 2019-10-19)
  
  My full code is as follows. I have tried to whittle it down to just the code that matters, but the problem I have is that i'm not sure what part of my network code is producing the problem. I've removed my code that loads and sifts through the CSV data because then my code would be too long.
  
  ```
  #include <iostream>
  #include <array>
  #include <random>
  #include <chrono>
  #include <iomanip>
  #include <fstream>
  #include <algorithm>
  #include <iomanip>
  #include <variant>
  #include <unordered_set>
  
  typedef std::variant<std::string,std::uint_fast16_t,bool,float> CSVType;
  
  /* ... functions to load CSV data ... */
  
  typedef float DataType;
  typedef DataType (*ActivationFuncPtr)(const DataType&);
  
  DataType step(const DataType& x, const DataType& threshold)
  {
      return x >= threshold ? 1 : 0;
  }
  
  DataType step0(const DataType& x)
  {
      return step(x,0);
  }
  
  DataType step05(const DataType& x)
  {
      return step(x,0.5);
  }
  
  DataType sigmoid(const DataType& x)
  {
      return DataType(1) / (DataType(1) + std::exp(-x));
  }
  
  DataType sigmoid_derivative(const DataType& x)
  {
      return x * (DataType(1) - x);
  }
  
  template<std::size_t NumInputs>
  class Neuron
  {
  public:
  
      Neuron()
      {
          RandomiseWeights();
      }
  
      void RandomiseWeights()
      {
          std::generate(m_weights.begin(),m_weights.end(),[&]()
          {
              return m_xavierNormalDis(m_mt);
          });
          m_biasWeight = 0;
  
          for(std::size_t i = 0; i < NumInputs+1; ++i)
              m_previousWeightUpdates[i] = 0;
      }
  
      DataType FeedForward(const std::array<DataType,NumInputs> inputValues)
      {
          DataType output = m_biasWeight;
          for(std::size_t i = 0; i < inputValues.size(); ++i)
              output += inputValues[i] * m_weights[i];
  
          m_inputValues = inputValues;
  
          return output;
      }
  
      std::array<DataType,NumInputs> Backpropagate(const DataType& error)
      {
          std::array<DataType,NumInputs> netInputOverWeight;
          for(std::size_t i = 0; i < NumInputs; ++i)
          {
              netInputOverWeight[i] = m_inputValues[i];
          }
  
          DataType netInputOverBias = DataType(1);
  
          std::array<DataType,NumInputs> errorOverWeight;
          for(std::size_t i = 0; i < NumInputs; ++i)
          {
              errorOverWeight[i] = error * netInputOverWeight[i];
          }
  
          DataType errorOverBias = error * netInputOverBias;
  
          for(std::size_t i = 0; i < NumInputs; ++i)
          {
              m_outstandingWeightAdjustments[i] = errorOverWeight[i];
          }
          m_outstandingWeightAdjustments[NumInputs] = errorOverBias;
  
          DataType errorOverNetInput = error;
  
          std::array<DataType,NumInputs> errorWeights;
          for(std::size_t i = 0; i < NumInputs; ++i)
          {
              errorWeights[i] = errorOverNetInput * m_weights[i];
          }
  
          return errorWeights;
      }
  
      void AdjustWeights(const DataType& learningRate, const DataType& momentum)
      {
          for(std::size_t i = 0; i < NumInputs; ++i)
          {
              DataType adjustment = learningRate * m_outstandingWeightAdjustments[i] + momentum * m_previousWeightUpdates[i];
              m_weights[i] = m_weights[i] - adjustment;
              m_previousWeightUpdates[i] = adjustment;
          }
          DataType adjustment = learningRate * m_outstandingWeightAdjustments[NumInputs] + momentum * m_previousWeightUpdates[NumInputs];
          m_biasWeight = m_biasWeight - adjustment;
          m_previousWeightUpdates[NumInputs] = adjustment;
      }
  
      const std::array<DataType,NumInputs>& GetWeights() const
      {
          return m_weights;
      }
  
      const DataType& GetBiasWeight() const
      {
          return m_biasWeight;
      }
  
  protected:
  
      static std::mt19937 m_mt;
      static std::uniform_real_distribution<DataType> m_uniformDisRandom;
      static std::uniform_real_distribution<DataType> m_xavierUniformDis;
      static std::normal_distribution<DataType> m_xavierNormalDis;
  
      std::array<DataType,NumInputs> m_weights;
      DataType m_biasWeight;
  
      std::array<DataType,NumInputs+1> m_previousWeightUpdates;
      std::array<DataType,NumInputs+1> m_outstandingWeightAdjustments;
  
      std::array<DataType,NumInputs> m_inputValues;
  };
  
  template<std::size_t NumInputs>
  std::mt19937 Neuron<NumInputs>::m_mt(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
  
  template<std::size_t NumInputs>
  std::uniform_real_distribution<DataType> Neuron<NumInputs>::m_uniformDisRandom(-1,1);
  
  template<std::size_t NumInputs>
  std::uniform_real_distribution<DataType> Neuron<NumInputs>::m_xavierUniformDis(-std::sqrt(6.f / NumInputs+1),std::sqrt(6.f / NumInputs+1));
  
  template<std::size_t NumInputs>
  std::normal_distribution<DataType> Neuron<NumInputs>::m_xavierNormalDis(0,std::sqrt(2.f / NumInputs+1));
  
  template<std::size_t NumNeurons>
  class ActivationLayer
  {
  public:
  
      ActivationLayer()
      :
          m_outputs({})
      {}
  
      virtual std::array<DataType,NumNeurons> GetOutputs() const final
      {
          return m_outputs;
      }
  
      virtual void CompleteBackprop(const DataType& learningRate, const DataType& momentum) final
      {
      }
  
  protected:
      std::array<DataType,NumNeurons> m_outputs;
  };
  
  template<std::size_t NumNeurons>
  class SigmoidActivation : public ActivationLayer<NumNeurons>
  {
  public:
  
      virtual std::array<DataType,NumNeurons> FeedForward(const std::array<DataType,NumNeurons>& inputValues)
      {
          for(std::size_t i = 0; i < NumNeurons; ++i)
              ActivationLayer<NumNeurons>::m_outputs[i] = sigmoid(inputValues[i]);
          return ActivationLayer<NumNeurons>::m_outputs;
      }
  
      virtual std::array<DataType,NumNeurons> Backpropagate(const std::array<DataType,NumNeurons> errors)
      {
          std::array<DataType,NumNeurons> backpropErrors;
          for(std::size_t i = 0; i < NumNeurons; ++i)
              backpropErrors[i] = errors[i] * sigmoid_derivative(ActivationLayer<NumNeurons>::m_outputs[i]);
          return backpropErrors;
      }
  };
  
  template<std::size_t NumInputs, std::size_t NumNeurons>
  class FullyConnectedLayer
  {
  public:
  
      FullyConnectedLayer()
      :
          m_neurons([=]()
          {
              std::array<Neuron<NumInputs>,NumNeurons> neurons;
              for(Neuron<NumInputs>& n : neurons)
                  n = Neuron<NumInputs>();
              return neurons;
          }())
      {
      }
  
      virtual std::array<DataType,NumNeurons> FeedForward(const std::array<DataType,NumInputs>& inputValues)
      {
          std::array<DataType,NumNeurons> outputValues;
          for(std::size_t i = 0; i < NumNeurons; ++i)
              outputValues[i] = m_neurons[i].FeedForward(inputValues);
          return outputValues;
      }
  
      /** \brief Take a sum of errors for each node and produce the errors for each input node in the previous layer.
       *
       */
  
      virtual std::array<DataType,NumInputs>
      Backpropagate(const std::array<DataType,NumNeurons> errors)
      {
          std::array<std::array<DataType,NumInputs>,NumNeurons> errorValues;
          for(std::size_t i = 0; i < NumNeurons; ++i)
          {
              errorValues[i] = m_neurons[i].Backpropagate(errors[i]);
          }
          std::array<DataType,NumInputs> returnErrors;
          std::fill(returnErrors.begin(),returnErrors.end(),0);
          for(std::size_t i = 0; i < NumNeurons; ++i)
          {
              for(std::size_t j = 0; j < NumInputs; ++j)
              {
                  returnErrors[j] += errorValues[i][j];
              }
          }
          return returnErrors;
      }
  
      virtual void CompleteBackprop(const DataType& learningRate, const DataType& momentum)
      {
          for(Neuron<NumInputs>& n : m_neurons)
              n.AdjustWeights(learningRate, momentum);
      }
  
      const Neuron<NumInputs>& operator[](const std::size_t& index) const
      {
          return m_neurons[index];
      }
  
      std::array<std::array<DataType,NumInputs>,NumNeurons> GetWeights() const
      {
          std::array<std::array<DataType,NumInputs>,NumNeurons> weights;
          for(std::size_t i = 0; i < NumNeurons; ++i)
          {
              weights[i] = m_neurons[i].GetWeights();
          }
          return weights;
      }
  
  protected:
      std::array<Neuron<NumInputs>,NumNeurons> m_neurons;
  };
  
  template<std::size_t I = 0, typename FuncT, typename... Tp>
  inline typename std::enable_if<I == sizeof...(Tp)>::type for_each(std::tuple<Tp...> &, FuncT)
  {
  }
  
  template<std::size_t I = 0, typename FuncT, typename... Tp>
  inline typename std::enable_if<I < sizeof...(Tp)>::type for_each(std::tuple<Tp...>& t, FuncT f)
  {
      f(std::get<I>(t)); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each<I + 1, FuncT, Tp...>(t, f); // process the next element of the tuple with the new output
  }
  
  template<std::size_t I = 0, typename FuncT, typename... Tp>
  inline typename std::enable_if<I == sizeof...(Tp)>::type for_each(const std::tuple<Tp...> &, FuncT)
  {
  }
  
  template<std::size_t I = 0, typename FuncT, typename... Tp>
  inline typename std::enable_if<I < sizeof...(Tp)>::type for_each(const std::tuple<Tp...>& t, FuncT f)
  {
      f(std::get<I>(t)); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each<I + 1, FuncT, Tp...>(t, f); // process the next element of the tuple with the new output
  }
  
  template<std::size_t I = 0, typename FuncT, typename O, typename FinalOutput, typename... Tp>
  inline typename std::enable_if<I == sizeof...(Tp)>::type for_each_get_final_output(std::tuple<Tp...> &, FuncT, O o, FinalOutput& finalOutput)
  {
      finalOutput = o;
  }
  
  template<std::size_t I = 0, typename FuncT, typename O, typename FinalOutput, typename... Tp>
  inline typename std::enable_if<I < sizeof...(Tp)>::type for_each_get_final_output(std::tuple<Tp...>& t, FuncT f, O o, FinalOutput& finalOutput)
  {
      auto newO = f(std::get<I>(t),o); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each_get_final_output<I + 1, FuncT, decltype(newO), FinalOutput, Tp...>(t, f, newO, finalOutput); // process the next element of the tuple with the new output
  }
  
  template<std::size_t I = 0, typename FuncT, typename O, typename... Tp>
  inline typename std::enable_if<I == 0>::type for_each_reverse_impl(std::tuple<Tp...>& t, FuncT f, O o)
  {
      f(std::get<0>(t),o);
  }
  
  template<std::size_t I = 0, typename FuncT, typename O, typename... Tp>
  inline typename std::enable_if<(I > 0)>::type
  for_each_reverse_impl(std::tuple<Tp...>& t, FuncT f, O o)
  {
      auto newO = f(std::get<I>(t),o); // call f, passing the Ith element of the std::tuple t and the existing output O
      for_each_reverse_impl<I - 1, FuncT, decltype(newO), Tp...>(t, f, newO); // process the next element of the tuple with the new output
  }
  
  template<typename FuncT, typename O, typename... Tp>
  inline void for_each_reverse(std::tuple<Tp...>& t, FuncT f, O o)
  {
      for_each_reverse_impl<sizeof...(Tp)-1, FuncT, O, Tp...>(t, f, o);
  }
  
  enum class LOSS_FUNCTION : std::uint_fast8_t
  {
      MEAN_SQUARE_ERROR,
      CROSS_ENTROPY
  };
  
  class ValidationOptions
  {
  public:
      enum class METRIC : std::uint_fast8_t { NONE, ACCURACY, LOSS };
  
      ValidationOptions()
      :
          m_validationSplit(.3f),
          m_enableLoss(false),
          m_lossFunction(LOSS_FUNCTION::MEAN_SQUARE_ERROR),
          m_enableAccuracy(false),
          m_outputFilter([](const DataType& x){ return x; }),
          m_earlyStoppingMetric(METRIC::NONE),
          m_earlyStoppingPatience(1.f),
          m_earlyStoppingDelta(1.f),
          m_earlyStoppingNumEpochsAverage(1)
      {}
  
      ValidationOptions& Loss(const bool enable = true, LOSS_FUNCTION lossFunction = LOSS_FUNCTION::MEAN_SQUARE_ERROR)
      {
          m_enableLoss = enable;
          m_lossFunction = lossFunction;
          return *this;
      }
  
      ValidationOptions& Split(const float dataSplitValidation)
      {
          m_validationSplit = dataSplitValidation;
          return *this;
      }
  
      ValidationOptions& Accuracy(const bool enable = true, ActivationFuncPtr outputFilter = [](const DataType& x){return x;})
      {
          m_enableAccuracy = enable;
          m_outputFilter = outputFilter;
          return *this;
      }
  
      ValidationOptions& EarlyStop(const bool enable = true,
                                   const METRIC metric = METRIC::ACCURACY,
                                   const float patience = .1f,
                                   const DataType delta = .01,
                                   const std::size_t epochNumToAverage = 10)
       {
           if(enable == false)
              m_earlyStoppingMetric = METRIC::NONE;
           else
              m_earlyStoppingMetric = metric;
  
           m_earlyStoppingPatience = patience;
           m_earlyStoppingDelta = delta;
           m_earlyStoppingNumEpochsAverage = epochNumToAverage;
  
           return *this;
       }
  
       float GetValidationSplit() const { return m_validationSplit; }
       bool Loss() const { return m_enableLoss; }
       LOSS_FUNCTION GetLossFunction() const { return m_lossFunction; }
       bool Accuracy() const { return m_enableAccuracy; }
       ActivationFuncPtr GetOutputFilter() const { return m_outputFilter; }
       METRIC GetEarlyStoppingMetric() const { return m_earlyStoppingMetric; }
       float GetEarlyStoppingPatience() const { return m_earlyStoppingPatience; }
       DataType GetEarlyStoppingDelta() const { return m_earlyStoppingDelta; }
       std::size_t GetEarlyStoppingNumEpochsAvg() const { return m_earlyStoppingNumEpochsAverage; }
  
  protected:
  
      float m_validationSplit;        /**< Percentage of the data set aside for validation */
  
      bool m_enableLoss;
      LOSS_FUNCTION m_lossFunction;   /**< Loss function to use */
  
      bool m_enableAccuracy;
      ActivationFuncPtr m_outputFilter;   /**< When measuring accuracy data is passed through this */
  
      METRIC m_earlyStoppingMetric;                   /**< The metric used to stop early */
      float m_earlyStoppingPatience;                  /**< Percentage of total epochs to wait before stopping early */
      DataType m_earlyStoppingDelta;                  /**< The amount that the early stopping metric needs to change in a single step before stopping */
      std::size_t m_earlyStoppingNumEpochsAverage;    /**< The number of epochs averaged over to smooth out the stopping metric */
  };
  
  template<typename... Layers>
  class NeuralNetwork
  {
  public:
  
      NeuralNetwork(Layers... layers)
      :
          m_layers(std::make_tuple(layers...))
      {
  
      }
  
      template<std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumTrainingRows>
      void Fit(const std::size_t& numberEpochs,
               const std::size_t& batchSize,
               DataType learningRate,
               const DataType& momentum,
               std::array<std::array<DataType,NumFeatures>,NumTrainingRows>& trainingData,
               std::array<std::array<DataType,NumOutputs>,NumTrainingRows>& trainingOutput,
               const ValidationOptions validationOptions,
               const bool linearDecayLearningRate = true,
               std::ostream& outputStream = std::cout)
      {
          std::size_t epochNumber = 0;
  
          // need to support more than just MSE to measure loss
          std::vector<DataType> lastEpochLoss(validationOptions.GetEarlyStoppingNumEpochsAvg(),0);
          DataType lastEpochLossAverage = std::numeric_limits<DataType>::max();
  
          std::vector<DataType> lastValidationAccuracys(validationOptions.GetEarlyStoppingNumEpochsAvg(),0);
          DataType lastValidationAccuraryAvg = 0;
  
          std::vector<std::size_t> randomIndices(NumTrainingRows,0);
          for(std::size_t i = 0; i < NumTrainingRows; ++i)
              randomIndices[i] = i;
  
          std::random_shuffle(randomIndices.begin(),randomIndices.end());
          // take some percentage as validation split
          // we do this by taking the first percentage of already shuffled indices and removing them
          // from what is available
          std::size_t numValidationRecords = NumTrainingRows*validationOptions.GetValidationSplit();
          std::size_t numTrainingRecords = NumTrainingRows - numValidationRecords;
          std::vector<std::size_t> validationRecords(numValidationRecords);
          for(std::size_t i = 0; i < numValidationRecords; ++i)
          {
              std::size_t index = randomIndices.back();
              randomIndices.pop_back();
              validationRecords[i] = index;
          }
  
          while(epochNumber < numberEpochs)
          {
              // shuffle the indices so that they are pulled into each batch randomly each time
              std::random_shuffle(randomIndices.begin(),randomIndices.end());
  
              DataType epochLoss = 0;
  
              std::tuple<Layers...> backupLayers = m_layers;
  
              for(std::size_t batchNumber = 0; batchNumber < std::ceil(numTrainingRecords / batchSize); ++batchNumber)
              {
                  std::array<DataType,NumOutputs> propagateError = {0};
  
                  std::size_t startIndex = batchNumber * batchSize;
                  std::size_t endIndex = startIndex + batchSize;
                  if(endIndex > numTrainingRecords)
                      endIndex = numTrainingRecords;
  
                  DataType batchLoss = 0;
  
                  for(std::size_t index = startIndex; index < endIndex; ++index)
                  {
                      std::size_t row = randomIndices[index];
                      const std::array<DataType,NumFeatures>& dataRow = trainingData[row];
                      const std::array<DataType,NumOutputs>& desiredOutputRow = trainingOutput[row];
  
                      // Feed the values through to the output layer
                      // use of "auto" is so this lambda can be used for all layers without
                      // me needing to do any fucking around
                      std::array<DataType,NumOutputs> finalOutput;
                      for_each_get_final_output(m_layers, [](auto& layer, auto o)
                      {
                          return layer.FeedForward(o);
                      }, dataRow, finalOutput);
  
                      DataType totalError = 0;
                      for(std::size_t i = 0; i < NumOutputs; ++i)
                      {
                          if(validationOptions.GetLossFunction() == LOSS_FUNCTION::MEAN_SQUARE_ERROR)
                              totalError += std::pow(desiredOutputRow[i] - finalOutput[i],2.0);
                          else if(validationOptions.GetLossFunction() == LOSS_FUNCTION::CROSS_ENTROPY)
                          {
                              if(NumOutputs == 1)
                              {
                                  // binary cross entropy
                                  totalError += (desiredOutputRow[i] * std::log(1e-15 + finalOutput[i]));
                              }
                              else
                              {
                                  // cross entropy
                              }
                          }
                      }
  
                      batchLoss += totalError;
                  }
  
                  batchLoss *= DataType(1) / (endIndex - startIndex);
  
                  for(std::size_t i = 0; i < NumOutputs; ++i)
                      propagateError[i] = batchLoss;
  
                  // update after every batch
                  for_each_reverse(m_layers, [](auto& layer, auto o)
                  {
                      auto errors = layer.Backpropagate(o);
                      return errors;
                  }, propagateError);
  
                  // once backprop is finished, we can adjust all the weights
                  for_each(m_layers, [&](auto& layer)
                  {
                      layer.CompleteBackprop(learningRate,momentum);
                  });
  
                  epochLoss += batchLoss;
              }
  
              epochLoss *= DataType(1) / numTrainingRecords;
  
              lastEpochLoss.erase(lastEpochLoss.begin());
              lastEpochLoss.push_back(epochLoss);
              DataType avgEpochLoss = 1.f * std::accumulate(lastEpochLoss.begin(),lastEpochLoss.end(),0.f) / (epochNumber < validationOptions.GetEarlyStoppingNumEpochsAvg() ? epochNumber+1 : lastEpochLoss.size());
  
              if(validationOptions.GetEarlyStoppingMetric() == ValidationOptions::METRIC::LOSS
                 && epochNumber > numberEpochs * validationOptions.GetEarlyStoppingPatience()
                 && avgEpochLoss > lastEpochLossAverage + validationOptions.GetEarlyStoppingDelta())
              {
                  // the loss average has decreased, so we should go back to the previous run and exit
                  std::cout   << "Early exit Loss Avg \n"
                              << "Last Epoch: " << lastEpochLossAverage << "\n"
                              << "This Epoch: " << avgEpochLoss << std::endl;
                  m_layers = backupLayers;
                  break;
              }
  
              lastEpochLossAverage = avgEpochLoss;
  
              // check for the error against the reserved validation set
              std::size_t numCorrect = 0;
              for(std::size_t row = 0; row < validationRecords.size(); ++row)
              {
                  const std::array<DataType,NumFeatures>& dataRow = trainingData[row];
                  const std::array<DataType,NumOutputs>& desiredOutputRow = trainingOutput[row];
  
                  std::array<DataType,NumOutputs> finalOutput;
                  for_each_get_final_output(m_layers, [](auto& layer, auto o)
                  {
                      return layer.FeedForward(o);
                  }, dataRow, finalOutput);
  
                  bool correct = true;
                  for(std::size_t i = 0; i < NumOutputs; ++i)
                  {
                      if(validationOptions.GetOutputFilter()(finalOutput[i]) != desiredOutputRow[i])
                          correct = false;
                  }
                  if(correct)
                      ++numCorrect;
              }
  
              DataType validationAccuracy = DataType(numCorrect) / numValidationRecords;
  
              lastValidationAccuracys.erase(lastValidationAccuracys.begin());
              lastValidationAccuracys.push_back(validationAccuracy);
              DataType avgValidationAccuracy = std::accumulate(lastValidationAccuracys.begin(),lastValidationAccuracys.end(),0.f) / (epochNumber < validationOptions.GetEarlyStoppingNumEpochsAvg() ? epochNumber+1 : lastValidationAccuracys.size());
  
              if(validationOptions.GetEarlyStoppingMetric() == ValidationOptions::METRIC::ACCURACY
                 && epochNumber > numberEpochs * validationOptions.GetEarlyStoppingPatience()
                 && avgValidationAccuracy < lastValidationAccuraryAvg - validationOptions.GetEarlyStoppingDelta())
              {
                  // the accuracy has decreased, so we should go back to the previous run and exit
                  std::cout   << "Early exit validation accuracy \n"
                              << "Last Epoch: " << lastValidationAccuraryAvg << "\n"
                              << "This Epoch: " << avgValidationAccuracy << std::endl;
                  m_layers = backupLayers;
                  break;
              }
  
              lastValidationAccuraryAvg = avgValidationAccuracy;
  
              outputStream << epochNumber << "," << epochLoss << "," << avgEpochLoss << "," << validationAccuracy << "," << avgValidationAccuracy << std::endl;
  
              learningRate -= learningRate / (numberEpochs-epochNumber);
  
              ++epochNumber;
          }
      }
  
      template<std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumEvaluationRows>
      void Evaluate(std::array<std::array<DataType,NumFeatures>,NumEvaluationRows> inputData,
                    std::array<std::array<DataType,NumOutputs>,NumEvaluationRows> correctOutputs,
                    DataType& loss,
                    DataType& accuracy,
                    ActivationFuncPtr outputFilter = [](const DataType& x){return x;})
      {
          loss = 0;
  
          std::size_t numCorrect = 0;
  
          for(std::size_t row = 0; row < NumEvaluationRows; ++row)
          {
              const std::array<DataType,NumFeatures>& dataRow = inputData[row];
              const std::array<DataType,NumOutputs>& outputRow = correctOutputs[row];
  
              // Feed the values through to the output layer
  
              std::array<DataType,NumOutputs> finalOutput;
              for_each_get_final_output(m_layers, [](auto& layer, auto o)
              {
                  layer.FeedForward(o);
                  return layer.GetOutputs();
              }, dataRow, finalOutput);
  
              DataType thisLoss = 0;
              for(std::size_t i = 0; i < NumOutputs; ++i)
                  thisLoss += outputRow[i] - finalOutput[i];
              loss += thisLoss * thisLoss;
  
              bool correct = true;
              for(std::size_t i = 0; i < NumOutputs; ++i)
              {
                  if(outputFilter(finalOutput[i]) != outputRow[i])
                      correct = false;
              }
              if(correct)
                  ++numCorrect;
          }
  
          loss *= DataType(1) / NumEvaluationRows;
          accuracy = DataType(numCorrect) / NumEvaluationRows;
      }
  
      template<std::size_t NumFeatures, std::size_t NumOutputs, std::size_t NumRecords>
      void Predict(std::array<std::array<DataType,NumFeatures>,NumRecords> inputData,
                   std::array<std::array<DataType,NumOutputs>,NumRecords>& predictions,
                   ActivationFuncPtr outputFilter = [](const DataType& x){return x;})
      {
          for(std::size_t row = 0; row < NumRecords; ++row)
          {
              const std::array<DataType,NumFeatures>& dataRow = inputData[row];
  
              // Feed the values through to the output layer
  
              std::array<DataType,NumOutputs> finalOutput;
              for_each_get_final_output(m_layers, [](auto& layer, auto o)
              {
                  return layer.FeedForward(o);
              }, dataRow, finalOutput);
  
              for(std::size_t i = 0; i < NumOutputs; ++i)
                  predictions[row][i] = outputFilter(finalOutput[i]);
          }
      }
  
  protected:
      std::tuple<Layers...> m_layers;
  };
  
  main()
  {
      std::vector<std::vector<CSVType>> trainingCSVData;
      /* load training CSV data */
  
      std::vector<std::vector<CSVType>> testCSVData;
      /* load test CSV data */
  
      std::cout << std::fixed << std::setprecision(80);
  
      std::ofstream file("error_out.csv", std::ios::out | std::ios::trunc);
      if(!file.is_open())
      {
          std::cout << "couldn't open file" << std::endl;
          return 0;
      }
  
      file << std::fixed << std::setprecision(80);
  
      /*
          Features
          1   pClass 1
          2   pClass 2
          3   pClass 3
          4   Sex female 1, male 0
          5   Age normalised between 0 and 1  age range 0 to 100
          6   Number siblings between 0 and 1   num range 0 to 8
          7   Number of parents / children        num range 0 to 9
          8   Ticket cost   between 0 and 1       num range 0 to 512.3292
          9   embarked S
          10  embarked Q
          11  embarked C
      */
  
      std::array<std::array<DataType,29>,891> inputData;
      std::array<std::array<DataType,1>,891> desiredOutputs;
  
      /* ... data that loads the titanic data into a series of features. Either class labels or normalised values (like age) */
  
      NeuralNetwork neuralNet{
          FullyConnectedLayer<29,256>(),
          SigmoidActivation<256>(),
          FullyConnectedLayer<256,1>(),
          SigmoidActivation<1>()
      };
  
      neuralNet.Fit(300,
                    1,
                    0.05,
                    0.25f,
                    inputData,
                    desiredOutputs,
                    ValidationOptions().Accuracy(true,step05).Loss(true,LOSS_FUNCTION::CROSS_ENTROPY).Split(0.3),
                    false,
                    file);
  
      file.close();
  
      return 0;
  }
  
  ```
  
  The data used is from the titanic problem that you can [download from Kaggle here](https://www.kaggle.com/c/titanic/data).
  
  The typical output file that's being generated is like this:
  
  ```
  0,-4.91843843460083007812500000000000000000000000000000000000000000000000000000000000,-4.91843843460083007812500000000000000000000000000000000000000000000000000000000000,0.65168541669845581054687500000000000000000000000000000000000000000000000000000000,0.65168541669845581054687500000000000000000000000000000000000000000000000000000000
  1,-6.14257431030273437500000000000000000000000000000000000000000000000000000000000000,-6.14257431030273437500000000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  2,-6.43130302429199218750000000000000000000000000000000000000000000000000000000000000,-6.43130302429199218750000000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  3,-6.58864736557006835937500000000000000000000000000000000000000000000000000000000000,-6.58864736557006835937500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  4,-6.70884752273559570312500000000000000000000000000000000000000000000000000000000000,-6.70884752273559570312500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  5,-6.78206682205200195312500000000000000000000000000000000000000000000000000000000000,-6.78206682205200195312500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  6,-6.86832284927368164062500000000000000000000000000000000000000000000000000000000000,-6.86832284927368164062500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  7,-6.92110681533813476562500000000000000000000000000000000000000000000000000000000000,-6.92110681533813476562500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  8,-6.96584081649780273437500000000000000000000000000000000000000000000000000000000000,-6.96584081649780273437500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  9,-7.02414274215698242187500000000000000000000000000000000000000000000000000000000000,-7.02414274215698242187500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  10,-7.06421041488647460937500000000000000000000000000000000000000000000000000000000000,-7.06421041488647460937500000000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000,0.65543073415756225585937500000000000000000000000000000000000000000000000000000000
  
  ```
  
  (shortened for space)
  
  This was previously working when the error I fed back through backprop was just the difference between the correct result and the prediction. [But I've since been told that I should be propagating back the Loss Functions error](https://www.reddit.com/r/learnmachinelearning/comments/dixyiy/i_think_ive_got_some_elements_of_neural_networks/), which I then implemented as Binary Cross Entropy.
  
  So either:
  
  1. I should not feed back the error from the loss function
  2. I am calculating the loss incorrectly
  3. My back propagation code is wrong
  4. I've got something horrible happening in my validation loop
  
  **Tags:** neural-networks, backpropagation
  
  > **Comments:**
  > **username_1**: Do a forward and backrpop by hand and compare to the values you get in your code. Speaking of which, that is a huge slab of code. I would say you would be better posting this to data science exchange as ai is generally for philosophical and theory, but this is such a large amount of code I don't think anyone at data science exchange would want to go over it. I think your best bet is to debug your code in segments, and see how that goes.  
  > **username_2**: username_1 seems correct to me. I don't usually mind reading a fair amount of code - I've answered similar implementation questions before, and have implemented neural networks in C (not C++ as here). But the size of the challenge to find your bug(s) in the attached code looks like too much effort for me. I will though have one guess.  
  > **username_3**: @username_0 No worries!  *(The overwhelming majority of my questions across stack get close;)*  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2019-10-19. Score: 1
  
  This is a guess, as I am not reading all that code!
  
  > This was previously working when the error I fed back through backprop was just the difference between the correct result and the prediction. But I've since been told that I should be propagating back the Loss Functions error, which I then implemented as Binary Cross Entropy.
  
  You may have been right before, and the advice you received wrong. Sort of.
  
  Backpropagation does not feed back the error value directly. It works exclusively with *gradients* of the error. So you typically start with a gradient value based on the objective function.
  
  If you backprop each function in the network a single step at a time, then you would start with $\nabla\_{\hat{y}}\mathcal{L}$, the gradient of the loss function with respect to the estimated values (i.e. output of the NN). Then from that you could calculate $\nabla\_{z}\mathcal{L}$, which is the gradient of the loss function with respect to the pre-activation values of the output layer.
  
  However, it is possible to combine steps analytically, and really common to start with the first gradient being $\nabla\_{z}\mathcal{L}$. That's because these pre-activation values are what you use to then calculate $\nabla\_{W}\mathcal{L}$ and $\nabla\_{b}\mathcal{L}$ - the gradients with respect to the layer's weights and biases, plus also to calculate $\nabla\_{z'}\mathcal{L}$ - the gradients with respect to pre-activation values of the previous layer.
  
  Either way, once you have $\nabla\_{z}\mathcal{L}$ for the output layer you can just run back though each layer in turn, repeating the same calcuation steps again and again. It is this repetition over each layer that looks like classic backprop in code.
  
  Assuming your neural network outputs a single value, the probablility of survival in the case of this dataset:
  
  * Your objective/loss function for binary classification *should* be Binary Cross Entropy.
  * Your output layer should have a sigmoid activation function.
  
  If you take the gradient of the loss function and backpropagate it throught that output layer's activation function, you end up with a starting delta $\nabla\_{z} \mathcal{L}$ i.e. the gradient of the loss function with respect to the pre-activation value of the first layer. This is a nice place to start the recursive backprop code.
  
  The value of this gradient is mathematically, per item:
  
  $$\nabla\_{z} \mathcal{L} = \hat{y}\_i - y\_i$$
  
  i.e. the difference the neural network estimate and the ground truth value. It takes this value because complex terms in the gradient of the loss function and the sigmoid activation function cancel out exactly. This is one reason why you often see sigmoid activation used with binary cross entropy - it is very convenient that the combination simplifies the gradient like this. You need to backpropagate each data point separately, and average the gradients for each batch/minibatch (actually you don't need to take this mean value, but doing so means that you don't need to adjust your learning rate as much to account for the number of items being processed per update step).
  
  It appears to of worked for you.
  
  Whether or not you were accidentally correct depends on whether you were getting an estimate for $\nabla\_{z}\mathcal{L}$ to start your backprop routine, or $\nabla\_{\hat{y}}\mathcal{L}$. I cannot find what that assumption is in your code, but the fact that it worked before and does not now suggests that you may of been correct, even though you may not of fully understood the maths.
  
  > **Comments:**
  > **username_3**: username_2, I'm disinclined to close this question since time went into this answer.  Do you think would could reform the question to be more suitable.  (The OP themself voted to close and tried to delete, but was blocked once there was an answer, so I think they'd be ok with editing the question.)  
  
  ---
  
  