7
\$\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\$

Your Answer

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

Browse other questions tagged or ask your own question.