From Python to Native C++: What I Learnt Building ML Systems from the Bottom Up

Modern ML libraries make inference look deceptively simple:

model = SentenceTransformer("some-model")
embedding = model.encode(text)

But underneath, the model is ultimately performing operations on arrays of numbers: memory accesses, matrix multiplications, normalisation, and other numerical kernels.

That led me to a simple question:

What actually happens underneath these high-level ML APIs?

To answer it, I built a native C++ inference implementation and progressively optimised it, from loading model weights to controlling memory and CPU execution.

The process became a journey through eight areas:

1. Understanding Arrays as Memory — how tensors are represented and accessed in contiguous memory.
2. Reading FP16 Weights — loading the model's FP16 weights and converting them to FP32.
3. Transformer with OpenBLAS — implementing QKV projections and attention using optimised matrix multiplication.
4. OpenMP to Parallelise the Heads — using thread-level parallelism to process independent attention heads.
5. AVX2 — using SIMD instructions to process multiple values per CPU instruction.
6. Numerical Stability and Cost of Softmax — implementing numerically stable softmax and understanding its computational cost.
7. Preallocation — allocating temporary memory once and reusing it across inference.
8. Buffer Swapping — avoiding tensor copies by alternating between reusable input and output buffers.

The main lesson was that writing inference code in C++ is not simply about replacing Python with a faster language.

It is about understanding the algorithms, memory layout, numerical kernels, parallelism, and hardware that ultimately determine performance.

1. Understanding Arrays as Memory

One of the first things I had to understand when moving to C++ was that an array is ultimately just a contiguous region of memory.

For example:

float array[6] = {a, b, c, d, e, f};

Memory:

[a][b][c][d][e][f]

If each float occupies 4 bytes, the next element is simply 4 bytes after the previous one. Array indexing is therefore closely related to pointer arithmetic:

array[i]

is effectively accessing the memory at:

base_address + i × sizeof(float)

This becomes important when working with tensors. A 2D tensor does not require a special 2D memory structure. For a row-major matrix:

[ a b c ]
[ d e f ]

the underlying memory is simply:

[a][b][c][d][e][f]

The shape tells us how to interpret that memory; the data itself remains a contiguous block.

This is also how embedding lookup works. An embedding matrix can be viewed as:

[vocabulary_size, embedding_dim]

where each token occupies embedding_dim consecutive values. To locate a token's embedding:

const float *src = embedding_weights.data + token_id * embedding_dim;

The expression:

token_id × embedding_dim

calculates the starting element of that token's memory chunk.

The important lesson is that tensor shapes are an interpretation of memory. Once I understood that, pointer arithmetic, tensor indexing, copying, and later optimisations such as cache-friendly access and buffer reuse became much easier to reason about.

2. Reading FP16 Model Weights

The model weights are stored as FP16 inside the safetensors file. I first access the raw tensor data as a char*, since the file provides the weights as a contiguous block of bytes.

I then interpret every two bytes as one 16-bit value:

const uint16_t *raw = reinterpret_cast<const uint16_t *>(raw_data);

Since FP16 is 16 bits, each uint16_t represents one FP16 value. I use std::bit_cast to interpret those 16 bits as a std::float16_t, then convert it to FP32:

std::float16_t value = std::bit_cast<std::float16_t>(raw[i]);
float weight = static_cast<float>(value);

The important distinction is that reinterpret_cast changes how I interpret the memory address, while std::bit_cast copies the underlying bit pattern into another type without changing the bits. The final static_cast performs the numerical conversion from FP16 to FP32.

I allocate an FP32 array for the converted weights:

tensor.data = new float[elements];

So the overall process is simply:

safetensors
↓
char* raw bytes
↓
uint16_t
↓
bit_cast → FP16
↓
static_cast → FP32
↓
float array

This gives the inference engine a contiguous FP32 representation that can be passed directly to the numerical kernels used later in the Transformer.

3. Implementing Transformer Attention in C++

The Transformer uses 16 attention heads with a hidden size of 1024:

16 × 64 = 1024

For each input tensor X, I first compute the Q, K and V projections:

Q = XWq
K = XWk
V = XWv

Each projection is a matrix multiplication:

X: [sequence_length, 1024]
W: [1024, 1024]
Q: [sequence_length, 1024]

Instead of implementing GEMM with nested loops, I use OpenBLAS:

cblas_sgemm(
    CblasRowMajor,
    CblasNoTrans,
    CblasTrans,
    M, N, K,
    1.0f,
    input.data, K,      // lda
    weight.data, K,     // ldb
    0.0f,
    output.data, N);    // ldc

The interesting part is understanding M, N, K and the leading dimensions lda, ldb, and ldc.

With CblasRowMajor, matrices are stored row by row:

[a b c]
[d e f]

becomes:

a b c d e f

The important part is that the number of columns and the row stride are not necessarily the same.

For a tightly packed 2 × 3 matrix:

[a b c]
[d e f]

the row stride is 3 because the next row starts three elements after the current row.

The starting address of row i is calculated as:

row_address = base_address + i × stride

An individual element at row i, column j is therefore:

element_address = base_address + i × stride + j

For example, with a stride of 3:

row 0 → base + 0 × 3
row 1 → base + 1 × 3
row 2 → base + 2 × 3

This becomes different when the matrix is stored inside a larger buffer.

For example:

[a b X]
[c d X]

contains a 2 × 2 matrix, but the stride is 3:

number of columns = 2
row stride        = 3

The X values are padding or belong to other data.

So when BLAS receives:

input.data, K   // lda

K tells it how many elements to move forward in memory to reach the next row.

This is why lda, ldb, and ldc are called leading dimensions. They describe the physical distance between rows in memory, rather than simply describing the mathematical number of columns.

For the input:

X: [M, K]

the matrix is contiguous, so its row stride is K:

input.data, K

For the stored weight:

W: [N, K]

its physical row stride is also K:

weight.data, K

But I use:

CblasTrans

so BLAS logically uses that stored matrix as:

Wᵀ: [K, N]

without creating a separate transposed copy.

Finally, the output:

Y: [M, N]

is contiguous, so its row stride is N:

output.data, N

This means M, N, and K describe the mathematical dimensions of the operation, while lda, ldb, and ldc describe how the matrices are laid out in memory.

4. OpenMP: Parallelising the Work

After optimising the matrix multiplication with BLAS, I looked at the remaining operations.

Some operations were naturally parallel.

For example, if I have thousands of independent tokens:

token 0
token 1
token 2
...
token N

and each token can be processed independently, there is no reason for one CPU thread to process all of them sequentially.

OpenMP makes this relatively straightforward:

#pragma omp parallel for
for (int i = 0; i < sequence_length; ++i)
{
    process_token(i);
}

This allows multiple threads to process different iterations concurrently.

But there was another opportunity.

Attention contains multiple independent heads:

head 0
head 1
head 2
...
head 15

Each head can perform its attention calculation independently before the heads are combined.

Therefore, attention can also be parallelised at the head level.

Conceptually:

             Attention
                 |
     +--+--+
     |           |           |
   Head 0      Head 1      Head 2 ...
     |           |           |
   compute     compute     compute

This can be especially useful when each head performs a relatively substantial amount of computation.

5. AVX2: Going Below OpenMP

OpenMP solves a different problem from SIMD.

OpenMP allows multiple threads to execute different pieces of work concurrently.

SIMD allows one CPU instruction to operate on multiple values simultaneously.

For example, a scalar implementation might process:

a0 × b0
a1 × b1
a2 × b2
a3 × b3

one operation at a time.

With AVX2, a 256-bit register can contain eight FP32 values:

[a0 a1 a2 a3 a4 a5 a6 a7]

and another:

[b0 b1 b2 b3 b4 b5 b6 b7]

A vector instruction can process all eight lanes together.

For operations such as dot products, this can make a significant difference.

I implemented AVX2-based operations for parts of the attention and softmax computation.

The improvement was particularly noticeable in dot-product-style operations.

6. Numerical Stability and Cost of Softmax

Softmax is mathematically simple:

softmax(xᵢ) = exp(xᵢ) / Σ exp(xⱼ)

For numerical stability, I first subtract the maximum value:

m = max(x)

eᵢ = exp(xᵢ - m)

s = Σ eᵢ

yᵢ = eᵢ / s

The subtraction reduces the risk of numerical overflow when computing the exponentials.

In practice, this requires several passes over each row:

find max
↓
subtract max + exp
↓
sum
↓
divide

The expensive part is the exp() operation, which is considerably more costly than basic addition or multiplication.

I used AVX2 to accelerate the vector operations where possible, but softmax remained a relatively expensive operation in the attention pipeline

7. Preallocating a Transformer Workspace

After optimising the computation itself, I found that memory allocation was becoming another source of overhead. Repeatedly allocating temporary arrays inside a loop is convenient:

for (...)
{
    float *buffer = new float[size];

    // work

    delete[] buffer;
}

but unnecessary when the required size and lifetime are predictable.

I benchmarked this against reusing preallocated memory:

allocate every iteration      14.810 ms
preallocated array[i]          4.083 ms
preallocated pointer           3.869 ms

The exact numbers depend on the workload and machine, but the difference was significant enough to affect inference performance.

I therefore created a reusable TransformerWorkspace that allocates the main temporary buffers once:

TransformerWorkspace

├── query
├── key
├── value
├── context
├── scores
├── intermediate
├── buffer_a
└── buffer_b

The workspace is sized according to the maximum sequence length, and the Transformer layers reuse these buffers throughout inference.

Instead of repeatedly doing:

allocate → compute → free

the runtime becomes:

allocate once

compute → reuse
compute → reuse
compute → reuse

This removes repeated allocation and deallocation from the critical path and makes the memory usage of the Transformer much more predictable.

8. Double Buffering Instead of Copying

I also wanted to avoid copying the entire tensor between Transformer layers. Instead, I use two reusable workspace buffers:

Tensor *input = &workspace.buffer_a;
Tensor *output = &workspace.buffer_b;

After each layer, I simply swap their pointers:

std::swap(input, output);

The buffers therefore alternate:

Layer 0: A → B
Layer 1: B → A
Layer 2: A → B
Layer 3: B → A

No tensor data is copied; only the pointers change.

The principle is simple: move references to data instead of moving the data itself.

Final Notes

After implementing many of these optimisations manually, I started looking more seriously at libraries such as oneDNN.

This was an important reality check.

There is a huge ecosystem of highly optimised numerical libraries that already solve many of these problems.

They can contain specialised implementations for:

• matrix multiplication;
• convolution;
• normalisation;
• attention-related operations;
• CPU-specific kernels;
• vectorised execution;
• thread scheduling;
• cache-aware algorithms.

The lesson is not:

"Implement everything yourself."

It is:

"Understand enough of the implementation to know when you should use a specialised library."

Writing a native inference engine was valuable precisely because it taught me what those libraries are doing for me.

Good ML engineering is not only about knowing how to use models. It is about understanding what the system underneath the model is actually doing.

My github: https://github.com/NgKaiWen7/InferenceEngine
Branch: immintrin