Backpropagation in Neural Networks

Backpropagation in Neural Networks | CodeSignal Learn

Introduction

Welcome back to "Training Neural Networks: The Backpropagation Algorithm"! You've made excellent progress so far, having learned about loss functions in our first lesson and gradient descent in our second. Today, we're diving into the heart of neural network training: backpropagation . In our previous lesson, we explored how gradient descent updates weights by moving in the direction opposite to the gradient of the loss function . But we left an important question unanswered: How do we actually calculate these gradients in a neural network with multiple layers and thousands or even millions of parameters? That's where backpropagation comes in. Backpropagation (short for "backward propagation of errors") is an efficient algorithm for computing these gradients. Today, we'll focus specifically on implementing the backward pass for a single dense layer, which will form the building block for training complete neural networks. By the end of this lesson, you'll understand how to: Calculate derivatives for different activation functions . Store necessary values during the forward pass. Implement the backward pass to calculate gradients. Connect these gradients to the gradient descent algorithm we learned previously. Let's embark on this crucial step in our neural network journey!

Understanding the Chain Rule for Backpropagation

f(g(x))

Join the 1M+ learners on CodeSignal

Be a part of our community of 1M+ users who develop and demonstrate their skills on CodeSignal

Before diving into code, let's build some intuition about how backpropagation works. The core mathematical principle behind backpropagation is the chain rule from calculus.

The chain rule allows us to calculate the derivative of composite functions. In simple terms, if we have a function $f(g(x))$ , the derivative with respect to x is:

\frac{df}{dx} = \frac{df}{dg} \cdot \frac{dg}{dx}

In neural networks, we have many nested functions. Consider a simple example with one dense layer:

We multiply inputs by weights: $z = x \cdot w + b$ .
We apply an activation function: $a = \sigma(z)$ .
We calculate a loss: $L = \text{loss}(a, y_{\text{true}})$ .

If we want to find how the loss changes with respect to the weights ( $\frac{dL}{dw}$ ), we apply the chain rule:

\frac{dL}{dw} = \frac{dL}{da} \cdot \frac{da}{dz} \cdot \frac{dz}{dw}

Backpropagation gets its name because we start at the output (the loss) and work backward through the network, calculating these gradients layer by layer. This is much more efficient than trying to directly compute the gradient of the loss with respect to each parameter.

For our implementation, we'll use matrix operations to handle batches of data efficiently, but the underlying principle remains the same.

Matrix Operations Helper Functions

Before implementing our activation functions and neural network layer, we need some helper functions for matrix operations in C++. These will replace the convenient matrix operations we would typically get from specialized libraries: C++#include <vector> #include <cmath> #include <algorithm> #include <iostream> #include <random> #include <stdexcept> // Matrix class for 2D operations class Matrix { public: std::vector<std::vector<double>> data; size_t rows, cols; Matrix(size_t r, size_t c) : rows(r), cols(c) { data.resize(rows, std::vector<double>(cols, 0.0)); } Matrix(const std::vector<std::vector<double>>& input_data) { data = input_data; rows = data.size(); cols = rows > 0 ? data[0].size() : 0; } // Matrix multiplication Matrix dot(const Matrix& other) const { if (cols != other.rows) { throw std::invalid_argument("Matrix dimensions don't match for multiplication"); } Matrix result(rows, other.cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < other.cols; ++j) { for (size_t k = 0; k < cols; ++k) { result.data[i][j] += data[i][k] * other.data[k][j]; } } } return result; } // Matrix transpose Matrix transpose() const { Matrix result(cols, rows); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[j][i] = data[i][j]; } } return result; } // Element-wise operations Matrix operator*(const Matrix& other) const { Matrix result(rows, cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][j] = data[i][j] * other.data[i][j]; } } return result; } Matrix operator+(const Matrix& other) const { Matrix result(rows, cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][j] = data[i][j] + other.data[i][j]; } } return result; } // Sum along axis (0 for columns, 1 for rows) Matrix sum(int axis = 0) const { if (axis == 0) { Matrix result(1, cols); for (size_t j = 0; j < cols; ++j) { for (size_t i = 0; i < rows; ++i) { result.data[0][j] += data[i][j]; } } return result; } else { Matrix result(rows, 1); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][0] += data[i][j]; } } return result; } } };#include <vector> #include <cmath> #include <algorithm> #include <iostream> #include <random> #include <stdexcept> // Matrix class for 2D operations class Matrix { public: std::vector<std::vector<double>> data; size_t rows, cols; Matrix(size_t r, size_t c) : rows(r), cols(c) { data.resize(rows, std::vector<double>(cols, 0.0)); } Matrix(const std::vector<std::vector<double>>& input_data) { data = input_data; rows = data.size(); cols = rows > 0 ? data[0].size() : 0; } // Matrix multiplication Matrix dot(const Matrix& other) const { if (cols != other.rows) { throw std::invalid_argument("Matrix dimensions don't match for multiplication"); } Matrix result(rows, other.cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < other.cols; ++j) { for (size_t k = 0; k < cols; ++k) { result.data[i][j] += data[i][k] * other.data[k][j]; } } } return result; } // Matrix transpose Matrix transpose() const { Matrix result(cols, rows); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[j][i] = data[i][j]; } } return result; } // Element-wise operations Matrix operator*(const Matrix& other) const { Matrix result(rows, cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][j] = data[i][j] * other.data[i][j]; } } return result; } Matrix operator+(const Matrix& other) const { Matrix result(rows, cols); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][j] = data[i][j] + other.data[i][j]; } } return result; } // Sum along axis (0 for columns, 1 for rows) Matrix sum(int axis = 0) const { if (axis == 0) { Matrix result(1, cols); for (size_t j = 0; j < cols; ++j) { for (size_t i = 0; i < rows; ++i) { result.data[0][j] += data[i][j]; } } return result; } else { Matrix result(rows, 1); for (size_t i = 0; i < rows; ++i) { for (size_t j = 0; j < cols; ++j) { result.data[i][0] += data[i][j]; } } return result; } } }; This Matrix class provides the essential operations we need for implementing our neural network layer, including matrix multiplication, transpose, element-wise operations, and summation along axes.

Activation Functions and Their Derivatives

\frac{da}{dz}

The DenseLayer Class Structure

Now, let's look at the structure of our DenseLayer class. This class encapsulates both the forward pass (which we've seen in previous lessons) and the backward pass (which we're focusing on today). C++class DenseLayer { private: size_t n_inputs; size_t n_neurons; std::string activation_fn_name; std::string weight_init_strategy; double weight_init_scale; // Function pointers for activation functions Matrix (*activation_fn)(const Matrix&); Matrix (*activation_derivative_fn)(const Matrix&); // Variables to store values needed for backpropagation Matrix inputs; Matrix z; Matrix output; public: Matrix weights; Matrix biases; Matrix d_weights; Matrix d_biases; DenseLayer(size_t n_inputs, size_t n_neurons, const std::string& activation_fn_name = "sigmoid", const std::string& weight_init_strategy = "random_scaled", double weight_init_scale = 0.01) : n_inputs(n_inputs), n_neurons(n_neurons), activation_fn_name(activation_fn_name), weight_init_strategy(weight_init_strategy), weight_init_scale(weight_init_scale), inputs(1, 1), z(1, 1), output(1, 1), weights(n_inputs, n_neurons), biases(1, n_neurons), d_weights(n_inputs, n_neurons), d_biases(1, n_neurons) { // Initialize weights and biases if (weight_init_strategy == "random_scaled") { std::random_device rd; std::mt19937 gen(rd()); std::normal_distribution<double> dist(0.0, weight_init_scale); for (size_t i = 0; i < n_inputs; ++i) { for (size_t j = 0; j < n_neurons; ++j) { weights.data[i][j] = dist(gen); } } } // Biases are initialized to zero (already done by Matrix constructor) // Set activation function and its derivative if (activation_fn_name == "sigmoid") { activation_fn = sigmoid; activation_derivative_fn = sigmoid_derivative; } else if (activation_fn_name == "relu") { activation_fn = relu; activation_derivative_fn = relu_derivative; } else if (activation_fn_name == "linear") { activation_fn = linear; activation_derivative_fn = linear_derivative; } else { throw std::invalid_argument("Unsupported activation: " + activation_fn_name); } } };class DenseLayer { private: size_t n_inputs; size_t n_neurons; std::string activation_fn_name; std::string weight_init_strategy; double weight_init_scale; // Function pointers for activation functions Matrix (*activation_fn)(const Matrix&); Matrix (*activation_derivative_fn)(const Matrix&); // Variables to store values needed for backpropagation Matrix inputs; Matrix z; Matrix output; public: Matrix weights; Matrix biases; Matrix d_weights; Matrix d_biases; DenseLayer(size_t n_inputs, size_t n_neurons, const std::string& activation_fn_name = "sigmoid", const std::string& weight_init_strategy = "random_scaled", double weight_init_scale = 0.01) : n_inputs(n_inputs), n_neurons(n_neurons), activation_fn_name(activation_fn_name), weight_init_strategy(weight_init_strategy), weight_init_scale(weight_init_scale), inputs(1, 1), z(1, 1), output(1, 1), weights(n_inputs, n_neurons), biases(1, n_neurons), d_weights(n_inputs, n_neurons), d_biases(1, n_neurons) { // Initialize weights and biases if (weight_init_strategy == "random_scaled") { std::random_device rd; std::mt19937 gen(rd()); std::normal_distribution<double> dist(0.0, weight_init_scale); for (size_t i = 0; i < n_inputs; ++i) { for (size_t j = 0; j < n_neurons; ++j) { weights.data[i][j] = dist(gen); } } } // Biases are initialized to zero (already done by Matrix constructor) // Set activation function and its derivative if (activation_fn_name == "sigmoid") { activation_fn = sigmoid; activation_derivative_fn = sigmoid_derivative; } else if (activation_fn_name == "relu") { activation_fn = relu; activation_derivative_fn = relu_derivative; } else if (activation_fn_name == "linear") { activation_fn = linear; activation_derivative_fn = linear_derivative; } else { throw std::invalid_argument("Unsupported activation: " + activation_fn_name); } } }; The key components to notice: Activation functions: We store function pointers to both the chosen activation_fn and its activation_derivative_fn. Caching variables: inputs will store the inputs to the layer. z will store the pre-activation outputs. output will store the post-activation outputs. d_weights and d_biases will store the gradients of weights and biases. This caching of intermediate values is crucial for backpropagation. We need to know these values during the backward pass to correctly compute the gradients.

The Forward Pass: Setting Up for Backpropagation

The forward pass not only computes the layer's output but also stores the necessary values for the backward pass. Let's examine the forward method implementation: C++Matrix forward(const Matrix& input_data) { inputs = input_data; z = inputs.dot(weights) + biases; // Pre-activation output output = activation_fn(z); return output; }Matrix forward(const Matrix& input_data) { inputs = input_data; z = inputs.dot(weights) + biases; // Pre-activation output output = activation_fn(z); return output; } This method: Stores the input values in inputs. Calculates and stores the pre-activation outputs z (the weighted sum plus bias). Applies the activation function and stores the results in output. Returns the output for use in subsequent layers. The key insight here is that we're caching all the intermediate values we'll need for the backward pass. This is essential for efficient computation of gradients during backpropagation. As you may recall from our previous lessons, this is how information flows forward through the network. Now, let's see how errors flow backward during the backpropagation process.

The Backward Pass: Calculating Gradients

\frac{dL}{dz}

Backpropagation in Action: A Practical Example

Let's now see how our backpropagation implementation works in practice with a simple example: C++int main() { // Single layer: 2 inputs, 3 neurons, sigmoid activation DenseLayer layer(2, 3, "sigmoid"); Matrix X_sample({{0.5, -0.2}}); // 1 sample, 2 features std::cout << "Input X: "; for (size_t j = 0; j < X_sample.cols; ++j) { std::cout << X_sample.data[0][j] << " "; } std::cout << std::endl; // Forward pass Matrix layer_output = layer.forward(X_sample); std::cout << "Layer output (after sigmoid): "; for (size_t j = 0; j < layer_output.cols; ++j) { std::cout << layer_output.data[0][j] << " "; } std::cout << std::endl; // Assume a dummy gradient from a hypothetical next layer or loss function // This is d(Loss)/d(layer_output) // Shape must match layer_output: (n_samples, n_neurons_in_layer) Matrix dummy_d_loss_wrt_layer_output({{0.1, -0.2, 0.05}}); std::cout << "Dummy dL/d(layer_output): "; for (size_t j = 0; j < dummy_d_loss_wrt_layer_output.cols; ++j) { std::cout << dummy_d_loss_wrt_layer_output.data[0][j] << " "; } std::cout << std::endl; // Backward pass Matrix d_loss_wrt_input = layer.backward(dummy_d_loss_wrt_layer_output); std::cout << "\nCalculated Gradients:" << std::endl; std::cout << " dL/d_weights (shape " << layer.d_weights.rows << "x" << layer.d_weights.cols << "):" << std::endl; for (size_t i = 0; i < layer.d_weights.rows; ++i) { std::cout << " "; for (size_t j = 0; j < layer.d_weights.cols; ++j) { std::cout << layer.d_weights.data[i][j] << " "; } std::cout << std::endl; } std::cout << " dL/d_biases (shape " << layer.d_biases.rows << "x" << layer.d_biases.cols << "):" << std::endl; std::cout << " "; for (size_t j = 0; j < layer.d_biases.cols; ++j) { std::cout << layer.d_biases.data[0][j] << " "; } std::cout << std::endl; std::cout << " dL/d_inputs (to pass to prev layer) (shape " << d_loss_wrt_input.rows << "x" << d_loss_wrt_input.cols << "):" << std::endl; std::cout << " "; for (size_t j = 0; j < d_loss_wrt_input.cols; ++j) { std::cout << d_loss_wrt_input.data[0][j] << " "; } std::cout << std::endl; return 0; }int main() { // Single layer: 2 inputs, 3 neurons, sigmoid activation DenseLayer layer(2, 3, "sigmoid"); Matrix X_sample({{0.5, -0.2}}); // 1 sample, 2 features std::cout << "Input X: "; for (size_t j = 0; j < X_sample.cols; ++j) { std::cout << X_sample.data[0][j] << " "; } std::cout << std::endl; // Forward pass Matrix layer_output = layer.forward(X_sample); std::cout << "Layer output (after sigmoid): "; for (size_t j = 0; j < layer_output.cols; ++j) { std::cout << layer_output.data[0][j] << " "; } std::cout << std::endl; // Assume a dummy gradient from a hypothetical next layer or loss function // This is d(Loss)/d(layer_output) // Shape must match layer_output: (n_samples, n_neurons_in_layer) Matrix dummy_d_loss_wrt_layer_output({{0.1, -0.2, 0.05}}); std::cout << "Dummy dL/d(layer_output): "; for (size_t j = 0; j < dummy_d_loss_wrt_layer_output.cols; ++j) { std::cout << dummy_d_loss_wrt_layer_output.data[0][j] << " "; } std::cout << std::endl; // Backward pass Matrix d_loss_wrt_input = layer.backward(dummy_d_loss_wrt_layer_output); std::cout << "\nCalculated Gradients:" << std::endl; std::cout << " dL/d_weights (shape " << layer.d_weights.rows << "x" << layer.d_weights.cols << "):" << std::endl; for (size_t i = 0; i < layer.d_weights.rows; ++i) { std::cout << " "; for (size_t j = 0; j < layer.d_weights.cols; ++j) { std::cout << layer.d_weights.data[i][j] << " "; } std::cout << std::endl; } std::cout << " dL/d_biases (shape " << layer.d_biases.rows << "x" << layer.d_biases.cols << "):" << std::endl; std::cout << " "; for (size_t j = 0; j < layer.d_biases.cols; ++j) { std::cout << layer.d_biases.data[0][j] << " "; } std::cout << std::endl; std::cout << " dL/d_inputs (to pass to prev layer) (shape " << d_loss_wrt_input.rows << "x" << d_loss_wrt_input.cols << "):" << std::endl; std::cout << " "; for (size_t j = 0; j < d_loss_wrt_input.cols; ++j) { std::cout << d_loss_wrt_input.data[0][j] << " "; } std::cout << std::endl; return 0; } This example: Creates a single DenseLayer with 2 inputs and 3 neurons. Performs a forward pass with a sample input. Simulates receiving gradients from the next layer using a dummy gradient. Performs a backward pass using this gradient. Prints the computed gradients for weights, biases, and inputs.

Output Discussion

When we run this code, we get output similar to the following: Input X: 0.5 -0.2 Layer output (after sigmoid): 0.500012 0.498415 0.501302 Dummy dL/d(layer_output): 0.1 -0.2 0.05 Calculated Gradients: dL/d_weights (shape 2x3): 0.0125 -0.025 0.00625 -0.005 0.01 -0.0025 dL/d_biases (shape 1x3): 0.025 -0.05 0.0125 dL/d_inputs (to pass to prev layer) (shape 1x2): 0.000694 -0.000182 Input X: 0.5 -0.2 Layer output (after sigmoid): 0.500012 0.498415 0.501302 Dummy dL/d(layer_output): 0.1 -0.2 0.05 Calculated Gradients: dL/d_weights (shape 2x3): 0.0125 -0.025 0.00625 -0.005 0.01 -0.0025 dL/d_biases (shape 1x3): 0.025 -0.05 0.0125 dL/d_inputs (to pass to prev layer) (shape 1x2): 0.000694 -0.000182 Looking at this output: Our input is a single sample with two features: [0.5, -0.2] . The forward pass produces outputs around 0.5 (since our weights are initialized close to zero, the sigmoid of values near zero is about 0.5). We provide a dummy gradient [0.1, -0.2, 0.05] representing how the loss would change if each output neuron's value changed slightly. The backward pass calculates: Gradients for each weight (d_weights). Gradients for each bias (d_biases). Gradients to pass to the previous layer (d_loss_wrt_input). This example demonstrates the full cycle of forward and backward passes for a single layer. In a complete neural network, we would perform this process for each layer, starting from the output and working backward (hence the name backpropagation).

Let's start implementing backpropagation by first defining our activation functions and their derivatives. These are crucial because the derivative of the activation function ( $\frac{da}{dz}$ ) is a key component in our chain rule calculations.

C++

// --- Activation Functions & Derivatives ---
Matrix sigmoid(const Matrix& x) {
    Matrix result(x.rows, x.cols);
    for (size_t i = 0; i < x.rows; ++i) {
        for (size_t j = 0; j < x.cols; ++j) {
            result.data[i][j] = 1.0 / (1.0 + std::exp(-x.data[i][j]));
        }
    }
    return result;
}

Matrix sigmoid_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = output.data[i][j] * (1.0 - output.data[i][j]);
        }
    }
    return result;
}

Matrix relu(const Matrix& x) {
    Matrix result(x.rows, x.cols);
    for (size_t i = 0; i < x.rows; ++i) {
        for (size_t j = 0; j < x.cols; ++j) {
            result.data[i][j] = std::max(0.0, x.data[i][j]);
        }
    }
    return result;
}

Matrix relu_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = (output.data[i][j] > 0) ? 1.0 : 0.0;
        }
    }
    return result;
}

Matrix linear(const Matrix& x) {
    return x; // Identity function
}

Matrix linear_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = 1.0;
        }
    }
    return result;
}

// --- Activation Functions & Derivatives ---
Matrix sigmoid(const Matrix& x) {
    Matrix result(x.rows, x.cols);
    for (size_t i = 0; i < x.rows; ++i) {
        for (size_t j = 0; j < x.cols; ++j) {
            result.data[i][j] = 1.0 / (1.0 + std::exp(-x.data[i][j]));
        }
    }
    return result;
}

Matrix sigmoid_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = output.data[i][j] * (1.0 - output.data[i][j]);
        }
    }
    return result;
}

Matrix relu(const Matrix& x) {
    Matrix result(x.rows, x.cols);
    for (size_t i = 0; i < x.rows; ++i) {
        for (size_t j = 0; j < x.cols; ++j) {
            result.data[i][j] = std::max(0.0, x.data[i][j]);
        }
    }
    return result;
}

Matrix relu_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = (output.data[i][j] > 0) ? 1.0 : 0.0;
        }
    }
    return result;
}

Matrix linear(const Matrix& x) {
    return x; // Identity function
}

Matrix linear_derivative(const Matrix& output) {
    Matrix result(output.rows, output.cols);
    for (size_t i = 0; i < output.rows; ++i) {
        for (size_t j = 0; j < output.cols; ++j) {
            result.data[i][j] = 1.0;
        }
    }
    return result;
}

Let's analyze each activation function and its derivative:

Sigmoid:
- The function is $\sigma(x) = \frac{1}{1 + e^{-x}}$ .
- Its derivative is $\sigma'(x) = \sigma(x) \cdot (1 - \sigma(x))$ .
- Note that we pass the output of the sigmoid function to calculate its derivative, as this is computationally efficient.
ReLU (rectified linear unit):
- The function is $\text{ReLU}(x) = \max(0, x)$ .
- Its derivative is 1 if $x > 0$ and 0 otherwise.
- We use C++'s conditional operator to efficiently compute this.
Linear:
- The function is simply $f(x) = x$ .
- Its derivative is always 1.
- We create a matrix of ones with the same shape as the input.

These activation functions and their derivatives are essential building blocks for both the forward and backward passes in our neural network layer.

Now we come to the heart of backpropagation: the backward method. This method calculates the gradients that will be used to update the weights during gradient descent.

C++

Matrix backward(const Matrix& d_loss_wrt_layer_output) {
    // Gradient of loss w.r.t. pre-activation output (dL/dz)
    // dL/dz = dL/dy * dy/dz, where y = activation(z)
    Matrix d_activation_output = d_loss_wrt_layer_output * activation_derivative_fn(output);

    // Gradient of loss w.r.t. weights (dL/dW)
    // dL/dW = dL/dz * dz/dW = d_activation_output * inputs.T
    d_weights = inputs.transpose().dot(d_activation_output);
    
    // Gradient of loss w.r.t. biases (dL/dB)
    // dL/dB = dL/dz * dz/dB = d_activation_output * 1 (summed over samples)
    d_biases = d_activation_output.sum(0);
    
    // Gradient of loss w.r.t. inputs of this layer (dL/dX_current_layer)
    // This will be d_loss_wrt_layer_output for the *previous* layer
    // dL/dX = dL/dz * dz/dX = d_activation_output . weights.T
    Matrix d_loss_wrt_prev_layer_output = d_activation_output.dot(weights.transpose());
    
    return d_loss_wrt_prev_layer_output;
}

Matrix backward(const Matrix& d_loss_wrt_layer_output) {
    // Gradient of loss w.r.t. pre-activation output (dL/dz)
    // dL/dz = dL/dy * dy/dz, where y = activation(z)
    Matrix d_activation_output = d_loss_wrt_layer_output * activation_derivative_fn(output);

    // Gradient of loss w.r.t. weights (dL/dW)
    // dL/dW = dL/dz * dz/dW = d_activation_output * inputs.T
    d_weights = inputs.transpose().dot(d_activation_output);
    
    // Gradient of loss w.r.t. biases (dL/dB)
    // dL/dB = dL/dz * dz/dB = d_activation_output * 1 (summed over samples)
    d_biases = d_activation_output.sum(0);
    
    // Gradient of loss w.r.t. inputs of this layer (dL/dX_current_layer)
    // This will be d_loss_wrt_layer_output for the *previous* layer
    // dL/dX = dL/dz * dz/dX = d_activation_output . weights.T
    Matrix d_loss_wrt_prev_layer_output = d_activation_output.dot(weights.transpose());
    
    return d_loss_wrt_prev_layer_output;
}

Let's break down what's happening:

Input: d_loss_wrt_layer_output is the gradient of the loss with respect to this layer's output. For the output layer, this comes directly from the loss function. For hidden layers, it's passed backward from the next layer.
Gradient with respect to pre-activation output: We first calculate $\frac{dL}{dz}$ by applying the chain rule: $\frac{dL}{dz} = \frac{dL}{dy} \cdot \frac{dy}{dz}$ , where y is the activation output and z is the pre-activation output.
Gradient with respect to weights: We calculate $\frac{dL}{dW}$ using the chain rule again: $\frac{dL}{dW} = \frac{dL}{dz} \cdot \frac{dz}{dW}$ . Since $z = X \cdot W + b$ , we have $\frac{dz}{dW} = X^T$ , leading to $\frac{dL}{dW} = X^T \cdot \frac{dL}{dz}$ .
Gradient with respect to biases: Similarly, $\frac{dL}{db} = \frac{dL}{dz} \cdot \frac{dz}{db}$ . Since $\frac{dz}{db} = 1$ , we sum $\frac{dL}{dz}$ over all samples.
Gradient with respect to inputs: Finally, we calculate $\frac{dL}{dX}$ to pass backward to the previous layer: $\frac{dL}{dX} = \frac{dL}{dz} \cdot \frac{dz}{dX}$ . Since $\frac{dz}{dX} = W^T$ , we get $\frac{dL}{dX} = \frac{dL}{dz} \cdot W^T$ .

Each of these calculations is a direct application of the chain rule, and together they form the core of the backpropagation algorithm.