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);
}
}
};
