9
\$\begingroup\$

I was implementing a simple neural network, and I noticed that, if I ever wanted to change the layers' activation functions, i would have had to completely rewrite some parts of the code, so I tried adding some flexibility to the network, using template meta-programming, because I thought that having most things resolved at compile time would have made everything faster.

The network is a feedforward neural-network with stochastic gradient descent, and for the matrix calculations I'm using the Eigen library.

I only implemented the Sigmoid and ReLU neuron types, but, with this architecture, it is trivial to add other ones; the same holds for the cost functions.

These neurons are implemented as structs, containing two functions, one computing the output of the neuron and the other computing its first derivative; these functions are static, to access them directly from the neuron type.

The neuron types are stored in the template parameters of the Network class, and to access them I wrapped them in a std::tuple:

using ActivationFuncsTuple = std::tuple<ActivationFuncs...>;

To access the type of neuron of a layer I used std::tuple_element_t, like this:

std::tuple_element_t<Layer, ActivationFuncsTuple>

Actually the neuron type always refers to Layer + 1, because the input layer doesn't have an associated neuron type

To iterate over the tuple, I took a static-for implementation from this thread Is it possible to develop static for loop in c++?, but I modified it a bit to allow going backwards.

Declaring a Network object is pretty straightforward, for example:

Network<Sigmoid, Sigmoid> net({784, 30, 10});
Network<Sigmoid, ReLU, Sigmoid> net({1000, 70, 40, 20});

The number of activation function is one less the the number of layers because the first layer is the input layer.

Training the network is just as easy:

net.train<CrossEntropyCost>(training_dataset, epochs, mini_batch_size, learning rate, regularization parameter, test_dataset);

Dataset is just a class I created to keep together the training or test inputs and the expected outputs, the make the code more compact.

What I'm most interested in is the correctness of the template-metaprogramming, as it is the first time I employed it in a "big" project. And can it be done in a cleaner and more elegant (maybe faster) way?

Network.h

#pragma once

#include <array>
#include <cmath>
#include <random>
#include <tuple>

#include <Eigen/Eigen>

#include "Functions.h"
#include "Helpers.h"

struct Dataset
{
    Eigen::MatrixXf samples;
    Eigen::MatrixXf expected_outputs;

    Dataset slice(int startIndex, int n) const
    {
        return { samples.block(0, startIndex, samples.rows(), n), expected_outputs.block(0, startIndex, expected_outputs.rows(), n) };
    }

    Eigen::Index size() const
    {
        return samples.cols();
    }
};

template <typename ...ActivationFuncs>
class CustomNeuralNetwork
{
private:
    static constexpr std::size_t n_layers = sizeof...(ActivationFuncs) + 1;

    using ActivationFuncsTuple = std::tuple<ActivationFuncs...>;

    std::array<Eigen::VectorXf, n_layers - 1> biases;
    std::array<Eigen::MatrixXf, n_layers - 1> weights;
public:
    CustomNeuralNetwork(const std::array<std::size_t, n_layers>& layersSizes)
    {
        std::mt19937 gen(std::random_device{}());
        std::normal_distribution<float> normalDist;

        for (int l = 0; l < layersSizes.size() - 1; l++)
        {
            biases[l] = Eigen::VectorXf(layersSizes[l + 1]);
            std::generate_n(biases[l].data(), biases[l].size(), std::bind(normalDist, gen));

            weights[l] = Eigen::MatrixXf(layersSizes[l + 1], layersSizes[l]);
            std::generate_n(weights[l].data(), weights[l].size(), std::bind(normalDist, gen));
            weights[l] /= std::sqrt(static_cast<float>(layersSizes[l + 1]));
        }
    }

    template<typename CostFunc>
    void train(Dataset training_dataset, int epochs, int mini_batch_size, float eta, float lambda, const Dataset& test_dataset = Dataset())
    {
        Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> permMatrix(training_dataset.size());
        permMatrix.setIdentity();
        std::mt19937 gen(std::random_device{}());

        for (int epoch = 0; epoch < epochs; epoch++)
        {
            std::shuffle(permMatrix.indices().data(), permMatrix.indices().data() + permMatrix.indices().size(), gen);
            training_dataset.samples *= permMatrix;
            training_dataset.expected_outputs *= permMatrix;

            for (int i = 0; i < training_dataset.size(); i += mini_batch_size)
            {
                update_mini_batch<CostFunc>(training_dataset.slice(i, mini_batch_size), eta, lambda, training_dataset.size());
            }

            if (test_dataset.size() > 0)
                std::cout << "Epoch " << epoch + 1 << " : " << evaluate(test_dataset) << " / " << test_dataset.size() << "\n";
        }
    }

    template<typename CostFunc>
    void update_mini_batch(const Dataset& mini_batch, float eta, float lambda, int n)
    {
        std::array<Eigen::MatrixXf, n_layers - 1> z;
        std::array<Eigen::MatrixXf, n_layers> a;

        a[0] = mini_batch.samples;

        static_for<0, n_layers - 1>::apply([&](auto l)
            {
                z[l] = biases[l].replicate(1, mini_batch.size());
                z[l].noalias() += weights[l] * a[l];
                a[l + 1] = std::tuple_element_t<l, ActivationFuncsTuple>::value(z[l]);
            });

        Eigen::MatrixXf delta = compute_delta<CostFunc>(z.back(), a.back(), mini_batch.expected_outputs);  //BACKPROPROPAGATION BEGINS HERE
        Eigen::VectorXf delta_b = delta.rowwise().sum();
        Eigen::MatrixXf delta_w = delta * a[a.size() - 2].transpose();

        static_for<n_layers - 2, 0>::apply([&](auto l)
            {
                delta = (weights[l].transpose() * delta).cwiseProduct(std::tuple_element_t<l - 1, ActivationFuncsTuple>::prime(a[l]));

                biases[l] -= delta_b * (eta / mini_batch.size());
                weights[l] = weights[l] * (1 - eta * lambda / n) - delta_w * (eta / mini_batch.size());

                delta_b.noalias() = delta.rowwise().sum();
                delta_w.noalias() = delta * a[l - 1].transpose();
            });

        biases[0] -= delta_b * (eta / mini_batch.size());
        weights[0] = weights[0] * (1 - eta * lambda / n) - delta_w * (eta / mini_batch.size());
    }

    template<typename CostFn>
    Eigen::MatrixXf compute_delta(const Eigen::MatrixXf& z, const Eigen::MatrixXf& a, const Eigen::MatrixXf& y)
    {
        if constexpr (std::is_same<std::tuple_element_t<n_layers - 2, ActivationFuncsTuple>, Sigmoid>::value && std::is_same<CostFn, CrossEntropyCost>::value)
        {
            return a - y;      //Specialization for Sigmoid neuron and Cross-Entropy cost function
        }
        else return CostFn::gradient(a, y).cwiseProduct(std::tuple_element_t<n_layers - 2, ActivationFuncsTuple>::prime(z, a));
    }

    Eigen::MatrixXf feedforward(Eigen::MatrixXf input)
    {
        static_for<0, n_layers - 1>::apply([&](auto l)
            {
                input = std::tuple_element_t<l, ActivationFuncsTuple>::value((weights[l] * input).colwise() + biases[l]);
            });

        return input;
    }

    int evaluate(const Dataset& test_dataset)
    {
        int accuracy = 0;
        Eigen::MatrixXf outputs = feedforward(test_dataset.samples);

        for (int i = 0; i < test_dataset.size(); i++)
        {
            int prediction;
            outputs.col(i).maxCoeff(&prediction);

            int correct;
            test_dataset.expected_outputs.col(i).maxCoeff(&correct);

            accuracy += (prediction == correct);
        }

        return accuracy;
    }
};

Functions.h

#pragma once

#include <Eigen/Eigen>

struct Sigmoid
{
    static Eigen::MatrixXf value(const Eigen::MatrixXf& z)
    {
        return 1.f / (1.f + z.array().exp().inverse());
    }

    static Eigen::MatrixXf prime(const Eigen::MatrixXf& a)
    {
        return a.array() * (1 - a.array());
    }
};

struct ReLU
{
    static Eigen::MatrixXf value(const Eigen::MatrixXf& z)
    {
        return (z.array().max(0.f));
    }

    static Eigen::MatrixXf prime(const Eigen::MatrixXf& z)
    {
        return (z.array() > 0.f).cast<float>();
    }
};


struct QuadraticCost
{
    static float cost(const Eigen::MatrixXf& output, const Eigen::MatrixXf& expected_outputs)
    {
        return 0.5f * (output - expected_outputs).colwise().squaredNorm().sum() / output.cols();
    }

    static Eigen::MatrixXf gradient(const Eigen::MatrixXf& outputs, const Eigen::MatrixXf& expected_outputs)
    {
        return (outputs - expected_outputs).cwiseProduct(Sigmoid::prime(outputs));
    }
};

struct CrossEntropyCost
{
    static float cost(const Eigen::MatrixXf& outputs, const Eigen::MatrixXf& expected_outputs)
    {
        return -(expected_outputs.array() * outputs.array().log() + expected_outputs.array() * (1.f - outputs.array()).log()).colwise().sum().sum() / outputs.cols();
    }

    static Eigen::MatrixXf gradient(const Eigen::MatrixXf& outputs, const Eigen::MatrixXf& expected_outputs)
    {
        return (outputs - expected_outputs).cwiseQuotient(Sigmoid::prime(outputs));
    }
};

Helpers.h

#pragma once

template <int First, int Last>
struct static_for
{
    template <typename Func>
    static constexpr void apply(Func&& f)
    {
        if constexpr(First < Last)
        {
            f(std::integral_constant<int, First>{});
            static_for<First + 1, Last>::apply(f);
        }
        else if constexpr (First > Last)
        {
            f(std::integral_constant<int, First>{});
            static_for<First - 1, Last>::apply(f);
        }
    }
};
\$\endgroup\$

1 Answer 1

1
\$\begingroup\$

Awesome work. Some further suggestions are mentioned as below.

Consider to add some unit tests with a testing framework

What I'm most interested in is the correctness of the template-metaprogramming

You can add some tests for verify the correctness. No matter the single layer cases or multiple layer cases, and the usage of various types of activation functions.

Consider to use size_t for sizes and epochs

The line for (int l = 0; l < layersSizes.size() - 1; l++) in CustomNeuralNetwork class constructor, you use int for iteration. The type size_t can be used here to represent all sizes of things.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.