How I built a Neural Network from Scratch | Part 1: Dense Layer

Photo by Hal Gatewood on Unsplash

This is 3 Part series, where I will be taking a deep dive into the internal functioning of a Neural Network and implement it with C++ from scratch. Take a look at my Repo for your reference.

Neuronite (https://github.com/Aditya-Dawadikar/Neuronite).

Why Neural Networks from Scratch?

The short answer is, I got bored treating the ML libraries like a black box. I had read about MLPs, CNNs and Transformers — I thought I understood what was happening under the hood. I figured: how hard could it be? A dense layer, a loss function, activation functions like ReLU and Sigmoid…right? RIGHT?

Mathematically I knew the drill. Forward pass = input_signal*weight + base, Backward pass = loss, gradients, updates. On paper.

I thought I could write a Neuron class, plug them into a layer and wire them up. Easy. But soon I realized the math doesn’t map cleanly to the classical Object Oriented Design. I had to go bare metal– use primitives for a smaller memory footprint. Afterall, the language of my choice was C++, and I wanted speed.

The Plan

After some brainstorming, reading up on the internet, and a few logical deductions later, I figured why would they call it TensorFlow? I had my small eureka moment when I realized that the Tensor was not just to address the data, rather those tensors were the building blocks for implementing a Neural Net. If I know how to operate on matrices, I’m halfway there. For my use case, 2D matrices would be enough. But then came the real question — how to compute partial derivatives on 2D matrices?

Now I have a starting point, Matrices will be the foundation. I coded the basic operations on Matrices — Addition, Subtraction, Multiplication with Scalar, Transpose, Dot Product and a helper display utility.

//matrix.hpp

#ifndef MATRIX_HPP
#define MATRIX_HPP

#include <vector>
#include <iostream>

class Matrix {

    public:
        int rows, cols;
        std::vector<std::vector<double>> data;

        Matrix();
        Matrix(int rows, int cols);
        Matrix(const std::vector<std::vector<double>>& values);

        static Matrix dot(const Matrix& A, const Matrix& B);
        Matrix transpose() const;

        Matrix operator+(const Matrix& other) const;
        Matrix operator-(const Matrix& other) const;
        Matrix operator*(const Matrix& other) const;
        Matrix operator*(double scalar) const;

        void print() const;
};

#endif

Next created a base class for all my Layers — Dense Layer and Activation Layers. Each layer consists of dimensions (rows, columns) and implement a forward pass, backward pass and an update function. Also, we will base our Layer class from the Matrix Abstraction.

// layer.hpp

#ifndef LAYER_HPP
#define LAYER_HPP

#include "matrix.hpp"

class Layer{
    public:
        virtual Matrix forward(const Matrix& input) = 0;
        virtual Matrix backward(const Matrix& grad_output) = 0;
        virtual void update(double learning_rate) = 0;
        virtual ~Layer() = default;
};

#endif

For weight initialization, I added a small utility util_random.hpp, capable of generating random numbers and for reproducibility, setting the seed value.

// util_random.hpp

#ifndef UTILS_RANDOM_HPP
#define UTILS_RANDOM_HPP

#include "matrix.hpp"

void set_random_seed(unsigned int seed);
void initialize_random(Matrix& mat, double min=-1.0, double max = 1.0);

#endif

The Math of Neural Networks

Now it was time to dig deeper into the math to implement a Dense Layer. A Dense Layer of a Multi Layer Perceptron consists of 2 things that change over the course of training — The Weights and the Biases. A neural network has a forward and backward phase.

Forward Pass

We take the raw input, pass it throw the layers, till it reach the output layer. This information flows in the form of signal that is influenced by the weight and the bias of a neuron. Mathematically, the layer’s output can be represented as follows:

Backward Pass

During this phase, our signals have already reached the last layer — output layer. This when we compute the loss, i.e how far is our predicted output from the expected output. There are various types of Loss functions, but for our use case, I relied on MSE (Mean Squared Error). Once using this loss, we will compute a gradient. A gradient tells us how much a small change in a parameter, affects the final loss. Mathematically, Gradient is represented as follows:

Update Weights and Bias

Finally, use the gradient to update the weight and bias. The weight will be changed in the direction opposite to the gradient.

Translating the equations to Matrices

Now we discuss the part where we compute the derivatives for real. The above equations boil down to the point where we can represent partial derivates in terms of dot product of two matrices.

Why does the weight gradient become a dot product?

Ps: Apologies for the poor visibility, of the above image.

Now that we know why the partial derivatives can be represented as simple matrix dot product, we will implement this. I designed a class called the DenseLayer that is inherited from the abstract class Layer.

#ifndef DENSELAYER_HPP
#define DENSELAYER_HPP

#include <matrix.hpp>
#include "layer.hpp"

class DenseLayer: public Layer{
    private:

        Matrix input_cache;
        Matrix d_weights;
        Matrix d_bias;

        std::pair<int,int> input_shape;
        std::pair<int,int> output_shape;

    public:
        Matrix weights;
        Matrix bias;

        DenseLayer(int input_dim, int output_dim);

        Matrix forward(const Matrix& input) override;
        Matrix backward(const Matrix& grad_output) override;
        void update(double learning_rate) override;
};

#endif

Then I implemented the functions in the dense_layer.cpp file

// denser_layer.cpp

#include "matrix.hpp"
#include "dense_layer.hpp"
#include "utils_random.hpp"
#include <cmath>

// DenseLayer constructor
// Initializes weights and biases, and allocates memory for gradients
// weight shape: (input_dim × output_dim)
// bias shape:   (1 × output_dim)
// d_weights:    (input_dim × output_dim)
// d_bias:       (1 × output_dim)
DenseLayer:: DenseLayer(int input_dim, int output_dim)
    : weights(input_dim, output_dim),
        bias(1, output_dim),
        d_weights(input_dim, output_dim),
        d_bias(1, output_dim){
    // dont worry about this line, just needed for keeping the random number
    // distribution tight
    double limit = std::sqrt(6.0 / (input_dim + output_dim));
    initialize_random(weights, -1*limit, 1*limit);
}

// Forward pass of the dense layer
// input shape:  (batch_size × input_dim)
// output shape: (batch_size × output_dim)
// Computes: Z = X · W + b
Matrix DenseLayer:: forward(const Matrix& input){
    // Cache input for use in backward pass
    input_cache = input;
    input_shape = {input.rows, input.cols};

    // Matrix multiplication: (batch_size × input_dim) · (input_dim × output_dim)
    Matrix output = Matrix::dot(input, weights);

    // Broadcast and add bias: bias is (1 × output_dim)
    output = output + bias;

    output_shape = {output.rows, output.cols};

    return output;
}

// Backward pass of the dense layer
// grad_output = ∂L/∂Z (gradient of loss w.r.t. layer output)
// grad_output shape: (batch_size × output_dim)
// Returns: ∂L/∂X = grad_input, shape: (batch_size × input_dim)
//
// Computes:
// d_weights = Xᵗ · ∂L/∂Z       (input_dim × output_dim)
// d_bias    = sum_rows(∂L/∂Z)  (1 × output_dim)
// grad_input = ∂L/∂Z · Wᵗ      (batch_size × input_dim)
Matrix DenseLayer:: backward(const Matrix& grad_output){
    // ∂L/∂W = inputᵗ · grad_output
    d_weights = Matrix::dot(input_cache.transpose(),grad_output);

    // ∂L/∂b = row-wise sum of grad_output
    d_bias = grad_output.row_wise_sum();

    // ∂L/∂X = grad_output · weightsᵗ
    Matrix grad_input = Matrix::dot(grad_output, weights.transpose());

    return grad_input;
}

// Update step: performs SGD on weights and bias
// W := W - η ∂L/∂W
// b := b - η ∂L/∂b
void DenseLayer:: update(double learning_rate){
    // apply gradient updates
    Matrix scaled_weights = d_weights*learning_rate;
    Matrix scaled_bias = d_bias*learning_rate;

    weights = weights - scaled_weights;
    bias = bias - scaled_bias;
}

I know this was heavy. But with that, we wrap up the first part of our journey — building a fully functional Dense Layer grounded in the foundational ideas of the McCulloch-Pitts model. From matrix multiplications to backpropagated gradients, we’ve taken a close look at how a single layer of a neural network processes and learns from data.

But this is just the beginning. In the next blog, we’ll peel back the layers even further — exploring activation functions that bring non-linearity to our networks and the loss functions that guide learning. Stay tuned, it only gets more exciting from here. 🚀🚀🚀

PS:

Find the implementation at Neuronite. I named it Neuronite because of two reasons, first it sounds similar to Neural Net, secondly its a combination of two words Neuron + Lite, “lite” because its a bare metal implementation of Neural Networks, stripped down to the essentials. And as a bonus, the suffix “-ite” gives it the feel of being a modular or foundational component of something larger. So yeah, I guess I just came up with a really cool name for my project.

It was very difficult to add equations to the Medium blog, so I had to type it in LaTex and then add the screenshots to this blog. Thanks for your understanding.