Why is my neural network giving me wildly incorrect error and not changing accuracy?

Question

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.

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, 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

Answer 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.