Neural Network in C Part 4: Hidden Layer Analysis

I want to look under the hood today. Obviously we can see the input layer and output layer. Now it would be cool to take a look under the hood and try to infer how this black box is actually working. Lets take a look at that hidden layer.

Remember the hidden layer consists of 128 neurons. We defined this earlier.

#define HIDDEN_SIZE 128        

What if we could visualize these neurons? We could train the network rigorously with say 100 epochs and then take a look at each neuron after training.

Now if you read the last tutorial, you realized that converting C code to assembly makes it faster, we will do that as well to save a lot of training time.

So the organized way to do this, and this is long overdue, is two split our program into two separate programs.

1. train.c will train our model.
2. main.c runs the model.

Now in order to take a peek into the hidden layer, a Python program could open up the model and run some image analysis. This could be done with matplotlib and numpy.

I spent some time adding some error handling elements and separating the training related code. The repository is also updated.

Here is train.c. Create a new file for this.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>
#include <time.h>

#define INPUT_SIZE 784
#define HIDDEN_SIZE 128
#define OUTPUT_SIZE 10
#define LEARNING_RATE 0.01
#define EPOCHS 1000

typedef struct {
    double weights[INPUT_SIZE][HIDDEN_SIZE];
    double biases[HIDDEN_SIZE];
} HiddenLayer;

typedef struct {
    double weights[HIDDEN_SIZE][OUTPUT_SIZE];
    double biases[OUTPUT_SIZE];
} OutputLayer;

void initialize_layer(HiddenLayer *hidden, OutputLayer *output) {
    srand(time(NULL));
    for (int i = 0; i < INPUT_SIZE; i++) {
        for (int j = 0; j < HIDDEN_SIZE; j++) {
            hidden->weights[i][j] = ((double)rand() / RAND_MAX) * 2 - 1;
        }
    }
    for (int i = 0; i < HIDDEN_SIZE; i++) {
        hidden->biases[i] = ((double)rand() / RAND_MAX) * 2 - 1;
        for (int j = 0; j < OUTPUT_SIZE; j++) {
            output->weights[i][j] = ((double)rand() / RAND_MAX) * 2 - 1;
        }
    }
    for (int i = 0; i < OUTPUT_SIZE; i++) {
        output->biases[i] = ((double)rand() / RAND_MAX) * 2 - 1;
    }
}

double sigmoid(double x) {
    return 1.0 / (1.0 + exp(-x));
}

double sigmoid_derivative(double x) {
    return x * (1.0 - x);
}

void forward_pass(HiddenLayer *hidden, OutputLayer *output, uint8_t *input, double *hidden_output, double *output_output) {
    for (int i = 0; i < HIDDEN_SIZE; i++) {
        double sum = hidden->biases[i];
        for (int j = 0; j < INPUT_SIZE; j++) {
            sum += input[j] / 255.0 * hidden->weights[j][i];
        }
        hidden_output[i] = sigmoid(sum);
    }

    for (int i = 0; i < OUTPUT_SIZE; i++) {
        double sum = output->biases[i];
        for (int j = 0; j < HIDDEN_SIZE; j++) {
            sum += hidden_output[j] * output->weights[j][i];
        }
        output_output[i] = sigmoid(sum);
    }
}

void backward_pass(HiddenLayer *hidden, OutputLayer *output, uint8_t *input, double *hidden_output, double *output_output, uint8_t label) {
    double output_errors[OUTPUT_SIZE];
    double hidden_errors[HIDDEN_SIZE];

    for (int i = 0; i < OUTPUT_SIZE; i++) {
        double target = (i == label) ? 1.0 : 0.0;
        output_errors[i] = (target - output_output[i]) * sigmoid_derivative(output_output[i]);
    }

    for (int i = 0; i < HIDDEN_SIZE; i++) {
        hidden_errors[i] = 0.0;
        for (int j = 0; j < OUTPUT_SIZE; j++) {
            hidden_errors[i] += output_errors[j] * output->weights[i][j];
        }
        hidden_errors[i] *= sigmoid_derivative(hidden_output[i]);
    }

    for (int i = 0; i < OUTPUT_SIZE; i++) {
        output->biases[i] += LEARNING_RATE * output_errors[i];
        for (int j = 0; j < HIDDEN_SIZE; j++) {
            output->weights[j][i] += LEARNING_RATE * output_errors[i] * hidden_output[j];
        }
    }

    for (int i = 0; i < HIDDEN_SIZE; i++) {
        hidden->biases[i] += LEARNING_RATE * hidden_errors[i];
        for (int j = 0; j < INPUT_SIZE; j++) {
            hidden->weights[j][i] += LEARNING_RATE * hidden_errors[i] * (input[j] / 255.0);
        }
    }
}

void train(HiddenLayer *hidden, OutputLayer *output, uint8_t *images, uint8_t *labels, int num_images, int num_epochs) {
    double hidden_output[HIDDEN_SIZE];
    double output_output[OUTPUT_SIZE];

    for (int epoch = 0; epoch < num_epochs; epoch++) {
        for (int i = 0; i < num_images; i++) {
            forward_pass(hidden, output, &images[i * INPUT_SIZE], hidden_output, output_output);
            backward_pass(hidden, output, &images[i * INPUT_SIZE], hidden_output, output_output, labels[i]);
        }
        printf("Epoch %d/%d completed\n", epoch + 1, num_epochs);
    }
}

uint32_t read_uint32(FILE *f) {
    uint32_t result;
    if (fread(&result, sizeof(result), 1, f) != 1) {
        perror("Failed to read uint32_t");
        exit(1);
    }
    return __builtin_bswap32(result);
}

uint8_t* load_mnist_images(const char *filename, int *num_images, int *image_size) {
    FILE *f = fopen(filename, "rb");
    if (!f) {
        perror("Failed to open file");
        exit(1);
    }

    uint32_t magic = read_uint32(f);
    *num_images = read_uint32(f);
    uint32_t rows = read_uint32(f);
    uint32_t cols = read_uint32(f);
    *image_size = rows * cols;

    uint8_t *images = (uint8_t*)malloc((*num_images) * (*image_size));
    if (!images) {
        perror("Failed to allocate memory for images");
        exit(1);
    }

    if (fread(images, *image_size, *num_images, f) != (size_t)(*num_images)) {
        perror("Failed to read images");
        free(images);
        exit(1);
    }

    fclose(f);
    return images;
}

uint8_t* load_mnist_labels(const char *filename, int *num_labels) {
    FILE *f = fopen(filename, "rb");
    if (!f) {
        perror("Failed to open file");
        exit(1);
    }

    uint32_t magic = read_uint32(f);
    *num_labels = read_uint32(f);

    uint8_t *labels = (uint8_t*)malloc(*num_labels);
    if (!labels) {
        perror("Failed to allocate memory for labels");
        exit(1);
    }

    if (fread(labels, 1, *num_labels, f) != (size_t)(*num_labels)) {
        perror("Failed to read labels");
        free(labels);
        exit(1);
    }

    fclose(f);
    return labels;
}

void save_model(const char *hidden_file, const char *output_file, HiddenLayer *hidden, OutputLayer *output) {
    FILE *f_hidden = fopen(hidden_file, "wb");
    if (!f_hidden) {
        perror("Failed to open hidden layer file");
        exit(1);
    }

    if (fwrite(hidden, sizeof(HiddenLayer), 1, f_hidden) != 1) {
        perror("Failed to write hidden layer");
        fclose(f_hidden);
        exit(1);
    }
    fclose(f_hidden);

    FILE *f_output = fopen(output_file, "wb");
    if (!f_output) {
        perror("Failed to open output layer file");
        exit(1);
    }

    if (fwrite(output, sizeof(OutputLayer), 1, f_output) != 1) {
        perror("Failed to write output layer");
        fclose(f_output);
        exit(1);
    }
    fclose(f_output);
}

int main() {
    int num_images, image_size, num_labels;
    uint8_t *images = load_mnist_images("mnist_data/train-images-idx3-ubyte", &num_images, &image_size);
    uint8_t *labels = load_mnist_labels("mnist_data/train-labels-idx1-ubyte", &num_labels);

    if (num_images != num_labels) {
        fprintf(stderr, "Number of images and labels do not match\n");
        exit(1);
    }

    HiddenLayer hidden;
    OutputLayer output;
    initialize_layer(&hidden, &output);

    int num_epochs = EPOCHS;
    train(&hidden, &output, images, labels, num_images, num_epochs);

    save_model("hidden_layer.bin", "output_layer.bin", &hidden, &output);

    free(images);
    free(labels);
    return 0;
}        

The review this code defines a simple neural network in C with an input layer (784 neurons), one hidden layer (128 neurons), and an output layer (10 neurons) for classifying MNIST digits.

It initializes weights and biases randomly, and uses a hand crafted sigmoid activation function. It also includes functions for forward and backward passes.

The network trains on the wonderful MNIST dataset over a specified number of epochs, adjusting weights using backpropagation and saves the trained model to a new binary file.

Here is our main.c code.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define INPUT_SIZE 784
#define HIDDEN_SIZE 128
#define OUTPUT_SIZE 10
#define TEST_IMAGES 10

typedef struct {
    double weights[INPUT_SIZE][HIDDEN_SIZE];
    double biases[HIDDEN_SIZE];
} HiddenLayer;

typedef struct {
    double weights[HIDDEN_SIZE][OUTPUT_SIZE];
    double biases[OUTPUT_SIZE];
} OutputLayer;

double sigmoid(double x) {
    return 1.0 / (1.0 + exp(-x));
}

void forward_pass(HiddenLayer *hidden, OutputLayer *output, uint8_t *input, double *hidden_output, double *output_output) {
    for (int i = 0; i < HIDDEN_SIZE; i++) {
        double sum = hidden->biases[i];
        for (int j = 0; j < INPUT_SIZE; j++) {
            sum += input[j] / 255.0 * hidden->weights[j][i];
        }
        hidden_output[i] = sigmoid(sum);
    }

    for (int i = 0; i < OUTPUT_SIZE; i++) {
        double sum = output->biases[i];
        for (int j = 0; j < HIDDEN_SIZE; j++) {
            sum += hidden_output[j] * output->weights[j][i];
        }
        output_output[i] = sigmoid(sum);
    }
}

uint8_t recognize(HiddenLayer *hidden, OutputLayer *output, uint8_t *image) {
    double hidden_output[HIDDEN_SIZE];
    double output_output[OUTPUT_SIZE];

    forward_pass(hidden, output, image, hidden_output, output_output);

    uint8_t recognized_digit = 0;
    double max_output = output_output[0];
    for (int i = 1; i < OUTPUT_SIZE; i++) {
        if (output_output[i] > max_output) {
            max_output = output_output[i];
            recognized_digit = i;
        }
    }
    return recognized_digit;
}

uint32_t read_uint32(FILE *f) {
    uint32_t result;
    if (fread(&result, sizeof(result), 1, f) != 1) {
        perror("Failed to read uint32_t");
        exit(1);
    }
    return __builtin_bswap32(result);
}

uint8_t* load_mnist_images(const char *filename, int *num_images, int *image_size) {
    FILE *f = fopen(filename, "rb");
    if (!f) {
        perror("Failed to open images file");
        exit(1);
    }
    printf("Opened image file: %s\n", filename);

    uint32_t magic = read_uint32(f);
    *num_images = read_uint32(f);
    uint32_t rows = read_uint32(f);
    uint32_t cols = read_uint32(f);
    *image_size = rows * cols;

    uint8_t *images = (uint8_t*)malloc((*num_images) * (*image_size));
    if (!images) {
        perror("Failed to allocate memory for images");
        exit(1);
    }

    if (fread(images, *image_size, *num_images, f) != (size_t)(*num_images)) {
        perror("Failed to read images");
        free(images);
        exit(1);
    }

    fclose(f);
    return images;
}

void load_model(const char *hidden_file, const char *output_file, HiddenLayer *hidden, OutputLayer *output) {
    printf("Loading hidden layer model from %s\n", hidden_file);
    FILE *f_hidden = fopen(hidden_file, "rb");
    if (!f_hidden) {
        perror("Failed to open hidden layer file");
        exit(1);
    }

    if (fread(hidden, sizeof(HiddenLayer), 1, f_hidden) != 1) {
        perror("Failed to read hidden layer");
        fclose(f_hidden);
        exit(1);
    }
    fclose(f_hidden);

    printf("Loading output layer model from %s\n", output_file);
    FILE *f_output = fopen(output_file, "rb");
    if (!f_output) {
        perror("Failed to open output layer file");
        exit(1);
    }

    if (fread(output, sizeof(OutputLayer), 1, f_output) != 1) {
        perror("Failed to read output layer");
        fclose(f_output);
        exit(1);
    }
    fclose(f_output);
}

double calculate_accuracy(HiddenLayer *hidden, OutputLayer *output, uint8_t *images, uint8_t *labels, int num_images) {
    int correct_predictions = 0;
    for (int i = 0; i < num_images; i++) {
        uint8_t recognized_digit = recognize(hidden, output, &images[i * INPUT_SIZE]);
        if (recognized_digit == labels[i]) {
            correct_predictions++;
        }
    }
    return (double)correct_predictions / num_images;
}

uint8_t* load_mnist_labels(const char *filename, int *num_labels) {
    FILE *f = fopen(filename, "rb");
    if (!f) {
        perror("Failed to open labels file");
        exit(1);
    }
    printf("Opened label file: %s\n", filename);

    uint32_t magic = read_uint32(f);
    *num_labels = read_uint32(f);

    uint8_t *labels = (uint8_t*)malloc(*num_labels);
    if (!labels) {
        perror("Failed to allocate memory for labels");
        exit(1);
    }

    if (fread(labels, 1, *num_labels, f) != (size_t)(*num_labels)) {
        perror("Failed to read labels");
        free(labels);
        exit(1);
    }

    fclose(f);
    return labels;
}

int main() {
    HiddenLayer hidden;
    OutputLayer output;

    load_model("hidden_layer.bin", "output_layer.bin", &hidden, &output);

    int num_images, image_size, num_labels;
    uint8_t *images = load_mnist_images("mnist_data/train-images-idx3-ubyte", &num_images, &image_size);
    uint8_t *labels = load_mnist_labels("mnist_data/train-labels-idx1-ubyte", &num_labels);

    if (num_images != num_labels) {
        fprintf(stderr, "Number of images and labels do not match\n");
        exit(1);
    }

    double accuracy = calculate_accuracy(&hidden, &output, images, labels, num_images);
    printf("Accuracy: %.2f%%\n", accuracy * 100);

    for (int i = 0; i < TEST_IMAGES; i++) {
        uint8_t recognized_digit = recognize(&hidden, &output, &images[i * INPUT_SIZE]);
        printf("Image %d: Recognized as %d, Actual %d\n", i + 1, recognized_digit, labels[i]);
    }

    free(images);
    free(labels);
    return 0;
}        

The main.c code loads pre-trained weights and biases for the hidden and output layers from our binary file.

The network uses a forward pass to compute activations using the custom sigmoid function.

The recognize function determines the predicted digit for an input image.

The program calculates the overall accuracy on the dataset and prints the recognized and actual digits for a few test images.

So lets convert train.c to assembly code. I want to run it all night long. 1000 epochs!

gcc -S -O3 -o train.s train.c
as -o train.o train.s
gcc -o train train.o -lm
./train        

Now, after about 11 hours of training our highly inefficient neural network, we have a .bin file. Lets bring in our old friend Python for a minute because we can visualize the neurons with matplotlib.

import numpy as np
import matplotlib.pyplot as plt
import os

INPUT_SIZE = 784
HIDDEN_SIZE = 128
OUTPUT_SIZE = 10
def load_hidden_layer(filename):
    if not os.path.exists(filename):
        raise FileNotFoundError(f"The file {filename} does not exist.")
    
    hidden_layer = {}
    with open(filename, 'rb') as f:
        hidden_layer['weights'] = np.fromfile(f, dtype=np.float64, count=INPUT_SIZE * HIDDEN_SIZE).reshape((INPUT_SIZE, HIDDEN_SIZE))
        hidden_layer['biases'] = np.fromfile(f, dtype=np.float64, count=HIDDEN_SIZE)
    
    if hidden_layer['weights'].shape != (INPUT_SIZE, HIDDEN_SIZE):
        raise ValueError("The shape of the loaded weights does not match the expected shape.")
    if hidden_layer['biases'].shape != (HIDDEN_SIZE,):
        raise ValueError("The shape of the loaded biases does not match the expected shape.")
    
    return hidden_layer

def save_weights_image(weights, layer_name, num_neurons, output_file):
    num_rows = (num_neurons + 7) // 8  # Ensure at most 8 neurons per row
    fig, axes = plt.subplots(num_rows, 8, figsize=(20, num_rows * 2))
    for i in range(num_neurons):
        row = i // 8
        col = i % 8
        axes[row, col].imshow(weights[:, i].reshape(28, 28), cmap='viridis')
        axes[row, col].set_title(f'Neuron {i+1}')
        axes[row, col].axis('off')
    for j in range(num_neurons, num_rows * 8):
        row = j // 8
        col = j % 8
        fig.delaxes(axes[row, col])  # Remove empty subplots
    plt.suptitle(f'Weights Visualization for {layer_name}')
    plt.savefig(output_file, bbox_inches='tight')
    plt.close()

hidden_layer = load_hidden_layer('hidden_layer.bin')
save_weights_image(hidden_layer['weights'], 'Hidden Layer', HIDDEN_SIZE, 'hidden_layer_weights.png')
        

The code defines two functions: load_hidden_layer and save_weights_image. The load_hidden_layer function reads the weights and biases for a hidden layer from a binary file, ensuring the file exists and the data matches the expected shapes. The weights are reshaped into a matrix, and biases are loaded as a vector of length HIDDEN_SIZE. The save_weights_image function visualizes the weights of each neuron as 28x28 pixel images.

Then it creates a matplotlib figure, where each subplot corresponds to a neuron's weights, and saves the figure to an output file.

Looking at these neurons you can kind of see various shapes and strange figures. I see eights and some triangles. Certain neurons seemed fine tuned to recognize a part of a digit like the middle of a zero, a "2" or "6". Let me know how I can improve this -Henry