Building a Neural Network Library

Welcome to the first lesson of "Building and Applying Your Neural Network Library", the fourth and final course in our "Neural Networks from Scratch using C++" path!

Throughout our journey so far, we've built a solid foundation in neural network principles. In the first course, we explored the fundamentals of neural networks, including perceptrons and the theory behind them. In the second course, we implemented forward propagation and activation functions. Most recently, in our third course, we mastered backpropagation and stochastic gradient descent, culminating in training a neural network on the diabetes dataset.

Now that we understand the core algorithms and mathematics, we're ready for the final stage: transforming our code into a proper, reusable neural network library. In this course, we'll take all the code we've produced in previous courses and restructure it into a more organized, modular framework — similar in spirit to popular libraries like TensorFlow, but built from scratch by us using C++!

Our first task is to modularize the core components we've already built: dense layers and activation functions. By the end of this lesson, you'll have created a well-structured C++ library that separates concerns using header files, namespaces, and proper compilation units, making your neural network code more maintainable and extensible.

Before we dive into implementation details, let's talk about why we're actually restructuring our code. So far, we've focused primarily on understanding the algorithms that power neural networks — the math, the theory, and the implementation of key concepts. While this understanding is crucial, there's another dimension to building effective ML systems: software engineering.

Software engineering principles are vital when building machine learning systems for several key reasons:

• Maintainability: As models grow in complexity, well-structured code becomes easier to debug and update.
• Reusability: Modular components can be reused across different projects.
• Testability: Isolated components with clear interfaces are easier to test.
• Collaboration: Well-organized code enables multiple people to work on different parts simultaneously.
• Extensibility: Adding new features becomes simpler when code is properly modularized.

In the industry, ML practitioners rarely write monolithic programs. Instead, they organize code into libraries and modules with clearly defined responsibilities. This is the approach we'll take as we build our neural network library.

Our library will be called neuralnets, and we'll structure it with separate header files and namespaces for different components. This structure separates concerns: activation functions live in their own namespace and files, layer implementations in another, and so on. As we continue through this course, we'll expand this structure to include losses, optimizers, and model classes.

Let's take a look at the complete directory structure we'll be building throughout this course:

neuralnets/
├── CMakeLists.txt
├── main.cpp
├── include/
│   └── neuralnets/
│       ├── neuralnets.hpp
│       ├── activations/
│       │   └── functions.hpp
│       ├── layers/
│       │   └── dense.hpp
│       ├── losses/
│       │   └── mse.hpp
│       ├── models/
│       │   ├── model.hpp
│       │   └── sequential.hpp
│       └── optimizers/
│           └── sgd.hpp
└── src/
    ├── activations/
    │   └── functions.cpp
    ├── layers/
    │   └── dense.cpp
    ├── losses/
    │   └── mse.cpp
    ├── models/
    │   ├── model.cpp
    │   └── sequential.cpp
    └── optimizers/
        └── sgd.cpp

This structure follows C++ library conventions and separates different concerns into their own compilation units:

• include/neuralnets/: Contains header files (.hpp) with class declarations and function prototypes.
• src/: Contains source files (.cpp) with the actual implementations.
• activations/: Contains activation functions and their derivatives.
• layers/: Houses different layer types (starting with our dense layer).
• losses/: Will contain loss functions for training (coming in later lessons).
• models/: Will include high-level model classes (coming in later lessons).
• optimizers/: Will contain optimization algorithms (coming in later lessons).
• CMakeLists.txt: Build configuration file for compiling our library.

In this lesson, we'll focus on implementing the activations/ and layers/ modules, along with the main library structure. The other modules will be added as we progress through the course, building up our complete neural network library step by step.

Let's begin by setting up the basic structure of our library. In C++, we organize code using header files for declarations and source files for implementations, along with namespaces to avoid naming conflicts.

First, let's create our build configuration file CMakeLists.txt:

cmake_minimum_required(VERSION 3.10)
project(neuralnets)

set(CMAKE_CXX_STANDARD 17)

# Find Eigen library
find_package(Eigen3 REQUIRED)

# Include directories
include_directories(include)

# Source files
set(SOURCES
    src/activations/functions.cpp
    src/layers/dense.cpp
)

# Create library
add_library(neuralnets ${SOURCES})
target_link_libraries(neuralnets Eigen3::Eigen)

# Create executable for testing
add_executable(neuralnets_test main.cpp)
target_link_libraries(neuralnets_test neuralnets Eigen3::Eigen)

This CMake configuration sets up our project to use C++17, finds the Eigen library for matrix operations, specifies our source files, and creates both a library and a test executable.

Next, we need to create our main namespace header. Let's create include/neuralnets/neuralnets.hpp:

#ifndef NEURALNETS_HPP
#define NEURALNETS_HPP

#include "neuralnets/activations/functions.hpp"
#include "neuralnets/layers/dense.hpp"

namespace neuralnets {
    // Main namespace that includes all subcomponents
    using namespace neuralnets::activations;
    using namespace neuralnets::layers;
}

#endif // NEURALNETS_HPP

This header file serves as the main entry point for our library, providing access to all components through include guards and namespace organization. The include guards prevent multiple inclusions of the same header, and the namespace structure keeps our code organized and prevents naming conflicts.

Now that we've set up our library structure, let's take a moment to understand some important C++ library concepts that will help you work more effectively with C++ projects and avoid common pitfalls.

• Header Guards and Include Management: The #ifndef, #define, and #endif directives we use create include guards that prevent the same header from being included multiple times in a single compilation unit. This prevents redefinition errors and speeds up compilation. Modern C++ also supports #pragma once as an alternative, but traditional include guards are more portable. Proper include management is crucial in C++ because, unlike higher-level languages, C++ uses a preprocessor that literally copies header content into source files.
• Compilation Units and Linking: In C++, each .cpp file is compiled separately into an object file, then all object files are linked together to create the final executable or library. This is why we separate declarations (in .hpp files) from definitions (in .cpp files). The linker resolves references between compilation units, which is why we can declare a function in a header and define it in a separate source file. Understanding this process helps explain why forward declarations work and why certain linking errors occur.
• Namespaces and Scope Resolution: We use namespaces like neuralnets::activations to organize our code and prevent naming conflicts. Unlike some languages that use file-based modules, C++ namespaces can span multiple files and can be nested. The using namespace directive brings names from one namespace into another, but should be used carefully to avoid polluting the global namespace. Our library structure uses nested namespaces to create a clear hierarchy of components.

Now let's move our activation functions to their own module. We'll create a header file for declarations (include/neuralnets/activations/functions.hpp):

#ifndef NEURALNETS_ACTIVATIONS_FUNCTIONS_HPP
#define NEURALNETS_ACTIVATIONS_FUNCTIONS_HPP

#include <Eigen/Dense>

namespace neuralnets {
namespace activations {

// Activation function declarations
Eigen::MatrixXd sigmoid(const Eigen::MatrixXd& x);
Eigen::MatrixXd sigmoid_derivative(const Eigen::MatrixXd& output);

Eigen::MatrixXd relu(const Eigen::MatrixXd& x);
Eigen::MatrixXd relu_derivative(const Eigen::MatrixXd& output);

Eigen::MatrixXd linear(const Eigen::MatrixXd& x);
Eigen::MatrixXd linear_derivative(const Eigen::MatrixXd& output);

} // namespace activations
} // namespace neuralnets

#endif // NEURALNETS_ACTIVATIONS_FUNCTIONS_HPP

And the corresponding implementation file (src/activations/functions.cpp):

#include "neuralnets/activations/functions.hpp"
#include <cmath>

namespace neuralnets {
namespace activations {

Eigen::MatrixXd sigmoid(const Eigen::MatrixXd& x) {
    return (1.0 / (1.0 + (-x.array()).exp())).matrix();
}

Eigen::MatrixXd sigmoid_derivative(const Eigen::MatrixXd& output) {
    return (output.array() * (1.0 - output.array())).matrix();
}

Eigen::MatrixXd relu(const Eigen::MatrixXd& x) {
    return x.cwiseMax(0.0);
}

Eigen::MatrixXd relu_derivative(const Eigen::MatrixXd& output) {
    return (output.array() > 0.0).cast<double>().matrix();
}

Eigen::MatrixXd linear(const Eigen::MatrixXd& x) {
    return x;
}

Eigen::MatrixXd linear_derivative(const Eigen::MatrixXd& output) {
    return Eigen::MatrixXd::Ones(output.rows(), output.cols());
}

} // namespace activations
} // namespace neuralnets

As you may recall from our previous courses, these functions implement three common activation functions and their derivatives:

1. Sigmoid: A smooth, S-shaped function that maps any input to a value between 0 and 1.
2. ReLU (rectified linear unit): Returns the input if positive; otherwise, returns 0.
3. Linear: Simply returns the input unchanged (used in regression tasks).

Each activation function has a corresponding derivative function used during backpropagation. Note that our derivative functions expect the output of the activation function rather than the original input, which is a common optimization.

The C++ implementation uses Eigen matrices for efficient linear algebra operations. The .array() method converts matrices to element-wise arrays for operations like exponentiation and comparison, while .matrix() converts back to matrix form.

Now let's implement our DenseLayer class in its own module. Since we've already built and understood the implementation details of this fully connected layer in previous courses, we'll focus on how it fits into our new modular structure and what interface it provides.

First, let's create the header file (include/neuralnets/layers/dense.hpp):

#ifndef NEURALNETS_LAYERS_DENSE_HPP
#define NEURALNETS_LAYERS_DENSE_HPP

#include <Eigen/Dense>
#include <string>
#include <functional>

namespace neuralnets {
namespace layers {

class DenseLayer {
public:
    // Constructor
    DenseLayer(int n_inputs, int n_neurons, const std::string& activation_fn_name = "sigmoid",
               const std::string& weight_init_strategy = "random_scaled", double weight_init_scale = 0.01);

    // Forward and backward pass methods
    Eigen::MatrixXd forward(const Eigen::MatrixXd& inputs);
    Eigen::MatrixXd backward(const Eigen::MatrixXd& d_loss_wrt_layer_output);

    // Getters for layer properties
    int get_n_inputs() const { return n_inputs_; }
    int get_n_neurons() const { return n_neurons_; }
    const std::string& get_activation_fn_name() const { return activation_fn_name_; }
    const Eigen::MatrixXd& get_output() const { return output_; }

private:
    // Layer parameters
    int n_inputs_;
    int n_neurons_;
    std::string activation_fn_name_;
    
    // Weights and biases
    Eigen::MatrixXd weights_;
    Eigen::VectorXd biases_;
    
    // Cached values for backpropagation
    Eigen::MatrixXd inputs_;
    Eigen::MatrixXd z_;
    Eigen::MatrixXd output_;
    
    // Function pointers for activation functions
    std::function<Eigen::MatrixXd(const Eigen::MatrixXd&)> activation_fn_;
    std::function<Eigen::MatrixXd(const Eigen::MatrixXd&)> activation_derivative_fn_;
    
    // Helper methods
    void initialize_weights(const std::string& strategy, double scale);
    void set_activation_functions(const std::string& activation_name);
};

} // namespace layers
} // namespace neuralnets

#endif // NEURALNETS_LAYERS_DENSE_HPP

And the implementation file (src/layers/dense.cpp):

#include "neuralnets/layers/dense.hpp"
#include "neuralnets/activations/functions.hpp"
#include <random>
#include <stdexcept>

namespace neuralnets {
namespace layers {

DenseLayer::DenseLayer(int n_inputs, int n_neurons, const std::string& activation_fn_name,
                       const std::string& weight_init_strategy, double weight_init_scale)
    : n_inputs_(n_inputs), n_neurons_(n_neurons), activation_fn_name_(activation_fn_name) {
    
    // Initialize weights and biases
    initialize_weights(weight_init_strategy, weight_init_scale);
    
    // Set activation functions
    set_activation_functions(activation_fn_name);
}

Eigen::MatrixXd DenseLayer::forward(const Eigen::MatrixXd& inputs) {
    inputs_ = inputs;
    z_ = inputs * weights_ + biases_.transpose().replicate(inputs.rows(), 1);
    output_ = activation_fn_(z_);
    return output_;
}

Eigen::MatrixXd DenseLayer::backward(const Eigen::MatrixXd& d_loss_wrt_layer_output) {
    Eigen::MatrixXd d_activation = activation_derivative_fn_(output_);
    Eigen::MatrixXd d_z = d_loss_wrt_layer_output.cwiseProduct(d_activation);
    
    // Compute gradients (not used for weight updates in this lesson, but computed for completeness)
    Eigen::MatrixXd d_weights = inputs_.transpose() * d_z;
    Eigen::VectorXd d_biases = d_z.colwise().sum();
    
    // Compute gradient with respect to inputs (for backpropagation to previous layer)
    Eigen::MatrixXd d_inputs = d_z * weights_.transpose();
    
    return d_inputs;
}

void DenseLayer::initialize_weights(const std::string& strategy, double scale) {
    std::random_device rd;
    std::mt19937 gen(rd());
    std::normal_distribution<double> dist(0.0, scale);
    
    weights_ = Eigen::MatrixXd(n_inputs_, n_neurons_);
    biases_ = Eigen::VectorXd::Zero(n_neurons_);
    
    if (strategy == "random_scaled") {
        for (int i = 0; i < weights_.rows(); ++i) {
            for (int j = 0; j < weights_.cols(); ++j) {
                weights_(i, j) = dist(gen);
            }
        }
    } else {
        throw std::invalid_argument("Unknown weight initialization strategy: " + strategy);
    }
}

void DenseLayer::set_activation_functions(const std::string& activation_name) {
    if (activation_name == "sigmoid") {
        activation_fn_ = activations::sigmoid;
        activation_derivative_fn_ = activations::sigmoid_derivative;
    } else if (activation_name == "relu") {
        activation_fn_ = activations::relu;
        activation_derivative_fn_ = activations::relu_derivative;
    } else if (activation_name == "linear") {
        activation_fn_ = activations::linear;
        activation_derivative_fn_ = activations::linear_derivative;
    } else {
        throw std::invalid_argument("Unknown activation function: " + activation_name);
    }
}

} // namespace layers
} // namespace neuralnets

This implementation demonstrates the power of our modular approach — the layer can now cleanly access activation functions through proper C++ includes and namespaces. The class uses function pointers (std::function) to store references to the appropriate activation functions, making the code flexible and efficient.

By organizing our layer implementation this way, we've created a self-contained component with a clear interface. Users of our library don't need to understand the internal mathematics — they simply create layer instances and call forward and backward methods. This encapsulation is a fundamental principle of good software design and makes our neural network library much more user-friendly.

Now that we have our core components in place, let's create a main program to test everything. Let's start by creating our main program to include our headers and set up a simple two-layer network in main.cpp:

#include <iostream>
#include <iomanip>
#include "neuralnets/neuralnets.hpp"

int main() {
    // Sample data - 2 samples, 3 features each
    Eigen::MatrixXd X_sample(2, 3);
    X_sample << 0.1, 0.2, -0.1,
                0.5, -0.3, 0.8;

    // Define layers directly
    neuralnets::layers::DenseLayer layer1(3, 5, "relu");
    neuralnets::layers::DenseLayer layer2(5, 2, "sigmoid");

    std::cout << "Network Architecture (Manual Sequential Execution):" << std::endl;
    std::cout << "  Layer 1: " << layer1.get_n_inputs() << " inputs -> " 
              << layer1.get_n_neurons() << " neurons, Activation: " 
              << layer1.get_activation_fn_name() << std::endl;
    std::cout << "  Layer 2: " << layer2.get_n_inputs() << " inputs -> " 
              << layer2.get_n_neurons() << " neurons, Activation: " 
              << layer2.get_activation_fn_name() << std::endl;

    return 0;
}

In this first part, we include our main library header, create a sample input matrix with 2 samples and 3 features each, define two layers (a hidden layer with ReLU activation and an output layer with sigmoid activation), and print information about the network architecture.

To compile and run this code, use the following commands:

mkdir build
cd build
cmake ..
make
./neuralnets_test

When we run this code, we get:

Network Architecture (Manual Sequential Execution):
  Layer 1: 3 inputs -> 5 neurons, Activation: relu
  Layer 2: 5 inputs -> 2 neurons, Activation: sigmoid

Now let's extend our program to perform a forward pass through the network and examine the results:

#include <iostream>
#include <iomanip>
#include "neuralnets/neuralnets.hpp"

int main() {
    // Sample data - 2 samples, 3 features each
    Eigen::MatrixXd X_sample(2, 3);
    X_sample << 0.1, 0.2, -0.1,
                0.5, -0.3, 0.8;

    // Define layers directly
    neuralnets::layers::DenseLayer layer1(3, 5, "relu");
    neuralnets::layers::DenseLayer layer2(5, 2, "sigmoid");

    std::cout << "Network Architecture (Manual Sequential Execution):" << std::endl;
    std::cout << "  Layer 1: " << layer1.get_n_inputs() << " inputs -> " 
              << layer1.get_n_neurons() << " neurons, Activation: " 
              << layer1.get_activation_fn_name() << std::endl;
    std::cout << "  Layer 2: " << layer2.get_n_inputs() << " inputs -> " 
              << layer2.get_n_neurons() << " neurons, Activation: " 
              << layer2.get_activation_fn_name() << std::endl;

    // Perform a forward pass manually
    Eigen::MatrixXd output_layer1 = layer1.forward(X_sample);
    Eigen::MatrixXd predictions = layer2.forward(output_layer1);
    
    std::cout << std::fixed << std::setprecision(8);
    std::cout << "\nInput X (shape " << X_sample.rows() << "x" << X_sample.cols() << "):" << std::endl;
    std::cout << X_sample << std::endl;
    
    std::cout << "\nPredictions from network (shape " << predictions.rows() << "x" << predictions.cols() << "):" << std::endl;
    std::cout << predictions << std::endl;

    std::cout << "\n--- Verifying individual layer outputs (after forward pass) ---" << std::endl;
    std::cout << "Layer 1 (ReLU) output (shape " << layer1.get_output().rows() << "x" << layer1.get_output().cols() << "):" << std::endl;
    std::cout << layer1.get_output() << std::endl;
    
    std::cout << "Layer 2 (Sigmoid) output (shape " << layer2.get_output().rows() << "x" << layer2.get_output().cols() << "):" << std::endl;
    std::cout << layer2.get_output() << std::endl;

    return 0;
}

Here, we perform a forward pass through both layers and print the results. This produces output similar to:

Network Architecture (Manual Sequential Execution):
  Layer 1: 3 inputs -> 5 neurons, Activation: relu
  Layer 2: 5 inputs -> 2 neurons, Activation: sigmoid

Input X (shape 2x3):
 0.10000000  0.20000000 -0.10000000
 0.50000000 -0.30000000  0.80000000

Predictions from network (shape 2x2):
0.49999341 0.49998954
0.50000289 0.49999404

--- Verifying individual layer outputs (after forward pass) ---
Layer 1 (ReLU) output (shape 2x5):
0.00000000 0.00000000 0.00090953 0.00315106 0.00060923
0.00000000 0.00000000 0.00000000 0.00000000 0.00421654

Layer 2 (Sigmoid) output (shape 2x2):
0.49999341 0.49998954
0.50000289 0.49999404

The outputs confirm our network is working as expected: the first layer's ReLU activation produces small positive values or zeros, and the second layer's sigmoid activation maps these values to numbers close to 0.5, which is expected for a randomly initialized network.

Note that we compile and run our program using CMake, which handles the compilation of our library and linking with Eigen. This build system approach is standard practice in professional C++ development and ensures that all dependencies are properly managed.

Congratulations! You've successfully transformed our neural network code into a well-structured C++ library. By separating our code into distinct header and source files with clear namespaces and responsibilities, we've taken a big step toward building a robust, reusable neural network library. This modular design provides the foundation for the rest of this course, where we'll continue expanding our library by adding modules for loss functions, optimizers, and a high-level model class.

The journey from understanding neural network principles to building a complete, well-structured library mirrors the path many practitioners take in the field. As you progress through this course, you'll not only deepen your understanding of neural networks but also develop valuable software engineering skills that are essential for real-world machine learning applications. The C++ approach gives you fine-grained control over memory management and performance, skills that are highly valued in production machine learning systems.

Now, it's time to get ready for some practice. Happy coding!