Version 0.0.3

GradieX — C++ Library Reference

The engine keeps v0.0.2’s flat, layer-contiguous memory layout and changes what runs over it: float instead of double, OpenMP simd instead of GCC unroll, accumulated gradients instead of assigned ones.

#Build Status

The suite was ported to the v0.0.3 API: the five-argument layer constructors, float fixtures, and zeroAllGradients / accumulateGradients / applyGradients in place of the removed backwardPass. It runs each test in a forked child, so a crash is reported against one row of the table rather than ending the run.

./tests/run_tests.sh              # 91 checks
./tests/run_tests.sh --asan       # + AddressSanitizer and UBSan
./tests/run_tests.sh Gradients    # only tests whose name contains "Gradients"

Reference values are computed in the test from the definitions, never captured from the implementation, so a failure means the code disagrees with the maths rather than with itself. Because the engine is float32, comparisons use a 1e-6 tolerance for hand-computed references, and snapshot the buffer for “this did not change” — 0.7f != 0.7, so a decimal literal would fail for the wrong reason.

#Building

C++20 is required, as in v0.0.2. The engine is header-only; main.cpp is the only translation unit.

clang++ -std=c++20 -O2 main.cpp -o main
./main

To act on the #pragma omp simd directives the compiler needs OpenMP enabled:

clang++ -std=c++20 -O2 -fopenmp-simd main.cpp -o main   # Apple clang: this works
g++-14  -std=c++20 -O3 -fopenmp      main.cpp -o main   # GCC, full OpenMP
clang++ -std=c++20 -O3 -fopenmp      main.cpp -o main   # LLVM clang, with libomp

core/ uses #pragma omp simd and nothing else — no parallel regions, no runtime calls — so -fopenmp-simd is sufficient and needs no libomp. It is worth up to 2.4× on the widest layers; see SIMD Pragmas.

Under -Wall -Wextra the build emits 18 warnings: 17 copies of 'restrict' type qualifier on return type has no effect and one unused parameter in identityDerivative, which ignores input because the derivative of the identity is 1 regardless. The first group is worth understanding rather than silencing blindly, though tests/run_tests.sh already passes -Wno-ignored-qualifiers.

#Overview

The architecture is v0.0.2's: a chain of HiddenLayer objects followed by an OutputLayer, each storing its weights as one flat vector rather than a vector of vectors, all owned by Network through raw pointers. Activation, initializer and loss behaviour are plain function pointers.

What changed is the training step. v0.0.2 computed a gradient and applied it immediately, one sample at a time. v0.0.3 zeroes the gradient buffers, accumulates over a batch, scales by 1/batchSize, and applies once:

trainBatch(startSample, batchSize, inputs, targets)
    zeroAllGradients()
    for each sample in the batch:
        forwardPass(x)
        accumulateGradients(y)      // deltas, then grad += delta * input
        batchLoss += computeLoss(y)
    scaleAllGradients(1 / batchSize)
    applyGradients()                // w -= learningRate * grad

As in every release, no weights move until every error signal has been computed — and now, until the whole batch has been seen. A hidden layer's error signal reads the next layer's weights, so updating a layer before propagating through it corrupts every layer behind it.

backwardPass and runEpoch are gone. Their work is split across accumulateGradients, zeroAllGradients, scaleAllGradients and applyGradients, which is what makes the batch loop expressible.

#Single Precision

Every double in core/ became a float: the weight, bias, gradient, delta, cache and scratch buffers, the function pointer signatures, the learning rate, and the loss accumulators.

The motivation is bandwidth. A float weight matrix is half the bytes, so twice as much of it fits in L1 and L2, and a 128-bit NEON register holds four lanes instead of two. The project's own TODO.md measures the current engine at about 11.5 GFLOP/s against a measured NEON fp32 FMA peak of 45.3, so the headroom is real — but it is unlocked by vectorisation, and the narrower type is the precondition rather than the win itself.

The cost is accuracy, and it is not small:

Three consequences follow for this engine.

The gradient check has to be re-tuned. A central-difference check compares an analytic gradient against (L(w+h) - L(w-h)) / 2h. v0.0.1 holds that to roughly 10^-11 in double; in float the subtraction cancels away most of the significand, and 10^-3 relative error is a realistic tolerance. The usual remedy is to run the numerical reference in double while the engine runs in float.

Model persistence would need rethinking. v0.0.1 writes 17 significant digits because that round-trips an IEEE-754 double exactly; the equivalent for float is 9. Save/load is not in v0.0.2 or v0.0.3, but it is on the roadmap, and the digit count is not a detail — it is what makes a reloaded model predict identically.

Cross entropy lost headroom, and it mattered. binaryCrossEntropy clamped its probability with eps = 1e-12, carried over unchanged from the double code. In float, 1.0f - 1e-12f rounds to exactly 1.0f, so the upper clamp did not clamp: once a sigmoid saturated, the loss hit log(0) = -inf and the derivative divided by zero, and inf * 0 is NaN — which reaches every weight on the next update and never washes out. A two-moons run diverged at epoch 243 with the loss already down at 0.003. The guard is now 1e-7f, representable either side of 1.0, and five regression tests pin it. v0.0.1 uses 1e-9 in double, which is correct there.

#SIMD Pragmas

v0.0.2's #pragma GCC unroll 4 directives were replaced throughout with #pragma omp simd, with reduction clauses where a loop accumulates:

// core/layer.h, the inner dot product
const float* __restrict__ w = neuronWeights + neuron * inputSize;
float sum = 0.0f;
#pragma omp simd reduction(+:sum)
for (int inputIdx = 0; inputIdx < inputSize; ++inputIdx)
    sum += w[inputIdx] * layerInputs[inputIdx];
neuronOutputsUnactivated[neuron] = sum + neuronBiases[neuron];

Where the two differ matters. #pragma GCC unroll 4 is an unrolling hint that changes nothing about the arithmetic. #pragma omp simd is an assertion that the loop is safe to vectorise, and reduction(+:sum) explicitly permits the compiler to reassociate the sum. Partial sums are then combined in an order the compiler chooses, so the result is no longer reproducible across compilers, flags or vector widths.

v0.0.2's README pinned its four-accumulator tree precisely to keep that order fixed. v0.0.3 gives it up. That is a defensible trade for speed, but it means any test comparing against a stored run has to become a tolerance test.

How much it buys, measured on an Apple M5 with Apple clang 21 at -O2 (microseconds per sample, minimum of nine trials — see bench/run_bench.sh):

The gain tracks layer width — the wider the dot product, the more there is to vectorise — and on the narrowest network it is nothing at all. Against v0.0.2 the same benchmark gives v0.0.3 1.2–2.1× without the flag and 1.2–3.2× with it, rising to 1.7–5.5× at batch 32. -O3 made no measurable difference at any size.

One pragma was also malformed, and enabling the flag proved it. In core/loss.h, softmaxCrossEntropyDerivative carried #pragma omp simd reduction(+:vectorScratchPad[:size]) over a loop that overwrites the array rather than reducing into it. An array-section reduction gives each lane a private copy initialised to the identity and sums them at the end — simply the wrong answer. It looked harmless only because nothing enabled the pragmas. Building the test suite with -fopenmp-simd made OutputLayer_ErrorSignal_FusedSoftmaxCrossEntropyIsProbabilityMinusTarget fail immediately. The clause is gone; the loop needs none. tests/run_tests.sh --simd now builds with the flag so this cannot come back unnoticed.

core/initializer.h deliberately carries no pragmas: the initializers draw from std::normal_distribution, which is stateful and cannot be vectorised. The header says so.

#Gradient Accumulation and Batching

This is the feature v0.0.2 dropped and v0.0.1 had. Two changes make it work.

Gradients accumulate. Layer::computeGradients now uses += where v0.0.2 used =:

// v0.0.2 -- assignment, one sample only
neuronWeightGradients[neuron * inputSize + inputIdx] = delta * savedInputs[inputIdx];

// v0.0.3 -- accumulation across a batch
dW[inputIdx] += delta * savedInputs[inputIdx];
neuronBiasGradients[neuron] += delta;

That makes zeroing mandatory. zeroAllGradients() runs at the top of every batch; skipping it means this batch's update carries every previous batch's gradient.

The batch size is a network property. Network(numHidden, numOutput, learningRate, batches = 1). batches = 1 reproduces v0.0.2's per-sample SGD; a value at or above the sample count is full-batch gradient descent. A non-positive value is treated as full-batch.

The final partial batch is handled: train() clamps each batch to min(batchSize, numSamples - sample) and scales by the actual count, so the last batch of an uneven split is averaged over its own size rather than the nominal one. Epoch loss is the sample-weighted mean of the batch means, which equals the plain mean over samples.

#Derivative Lookup

In v0.0.2 a layer is handed its activation and its derivative as two independent arguments. Nothing checks that they match, so sigmoid paired with reluDerivative compiles and trains to the wrong answer without a diagnostic. It is the reason v0.0.2's Python binding selects matched pairs by name.

v0.0.3 removes the argument instead. The Layer constructor takes five arguments and derives both derivatives itself:

Layer(int width, int inputSize,
      ActivationFunction activation,
      initFunction init,
      LossFunction loss)
{
    this->activationFunction = activation;
    this->activationFunctionDerivative = activations::fetchAssociatedDerivative(activation);
    this->lossFunctionDerivative      = loss::fetchAssociatedDerivative(loss);
    ...
}

Both lookups are chains of function pointer comparisons. The loss side is complete — all six losses map to their derivatives. The activation side is not:

The two lookup functions also carry an unrelated defect that costs more than it looks — see One Translation Unit Only.

#One Translation Unit Only

Both lookup functions are defined in headers without inline:

// core/activations.h:349   -- no 'inline'
ActivationFunctionDerivative fetchAssociatedDerivative(ActivationFunction activation) { ... }

// core/loss.h:288          -- no 'inline'
LossFunctionDerivative fetchAssociatedDerivative(LossFunction lossFunction) { ... }

Every other function in core/ is inline. These two were not, so each translation unit that included the header emitted its own definition and the linker saw duplicates:

$ clang++ -std=c++20 a.cpp b.cpp -o out     # both #include "core/network.h"
duplicate symbol 'core::functions::activations::fetchAssociatedDerivative(...)'
duplicate symbol 'core::functions::loss::fetchAssociatedDerivative(...)'
ld: 2 duplicate symbols

The project has not noticed because main.cpp is the only translation unit. It stops mattering the moment a second one appears — which is exactly what tests/run_tests.sh and the pybind11 build both do. The fix is one keyword on each.

#The __restrict__ Return Qualifier

17 of the build's 18 warnings are this one:

warning: 'restrict' type qualifier on return type has no effect
         [-Wignored-qualifiers]

Every activation is declared inline float* __restrict__ name(...), and the two aliases match:

using ActivationFunction = float* __restrict__ (*)(const float* __restrict__ input,
                                                  float* __restrict__ vectorScratchPad,
                                                  size_t size);

A restrict qualifier on a return type is meaningless — it qualifies the returned pointer value, which the caller immediately copies into its own unqualified variable. The compiler says so and moves on. But it is not merely noise: the qualifier is part of the type for the purpose of matching function pointers, so the alias and the functions have to agree. They now do, which is why the header compiles.

The clean fix is to drop __restrict__ from the return type in both the aliases and all 17 function declarations at once. Dropping it from only one side reintroduces a type mismatch. Suppressing the warning — which tests/run_tests.sh already does with -Wno-ignored-qualifiers — hides the noise without removing the coupling.

#Namespaces and Headers

Unchanged from v0.0.2. Including core/network.h pulls in everything.

#Function Pointer Types

Identical in shape to v0.0.2, with float in place of double.

activations::ActivationFunction

float* __restrict__ (*)(const float* __restrict__ input,
                        float* __restrict__ vectorScratchPad,
                        size_t size)

Applies an activation elementwise over size values, writing into vectorScratchPad and returning it. The pad is sized once by the layer and reused, so the hot path allocates nothing. The return qualifier is inert; see above.

activations::ActivationFunctionDerivative

Same signature. Takes the pre-activation values (z), not the activated outputs.

initializers::initFunction

float* (*)(float* __restrict__ vectorScratchPad,
           size_t size, size_t size2)

Fills size * size2 values. For xavier and he the arguments are fan-in and fan-out.

loss::LossFunction

float (*)(const float* __restrict__ actual,
          const float* __restrict__ predicted,
          size_t size)

Returns a scalar. Note the argument order: target first, prediction second.

loss::LossFunctionDerivative

float* (*)(const float* __restrict__ actual,
           const float* __restrict__ predicted,
           float* __restrict__ vectorScratchPad,
           size_t size)

Writes dL/dpredicted into the scratch pad and returns it.

#class core::Layer

Defined in core/layer.h. The base class for both layer kinds. All members are public.

#Constructors

Layer() noexcept = default

Leaves the object empty.

Warning

The function pointer members were indeterminate after this constructor rather than null — fixed in v0.0.3; all five now carry = nullptr. It was known bug 6. Calling calculateOutputs() on such an object jumps to an arbitrary address. The fix remains = nullptr initialisers on all five.

Layer(int width, int inputSize, ActivationFunction, initFunction, LossFunction)

Five arguments, down from seven. The activation derivative and the loss derivative are looked up rather than passed — see Derivative Lookup. Sizes every internal buffer. Does not initialize the weights; call initializeWeightsAndBiases() for that.

#Data Members

#Member Functions

void initializeWeightsAndBiases()

Calls the stored initializer over the whole weight array with (inputSize, layerWidth) as fan-in and fan-out, then zeroes every bias.

void calculateOutputs(const float* __restrict__ input)

Copies input into layerInputs, computes z[n] = sum_i w[n][i] * layerInputs[i] + b[n], then applies the activation into neuronOutputsActivated. The inner loop carries #pragma omp simd reduction(+:sum), and the loop bounds and base pointers are hoisted into locals first.

Note that the dot product reads layerInputs, the copy, rather than input — v0.0.2 read input directly. The values are identical; the change lets the compiler assume the __restrict__ guarantee over a buffer the layer owns.

input must point to at least inputSize floats. This is not checked.

void zeroGradients()

New in v0.0.3. Sets both gradient buffers to zero. Must run before each batch; skipping it carries the previous batch's gradients into this one.

void scaleGradients(float scale)

New in v0.0.3. Multiplies both gradient buffers in place. Used with 1 / batchSize to turn an accumulated sum into a mean.

void computeGradients()

Requires neuronDeltas to be current. Accumulates dL/dw[n][i] += delta[n] * layerInputs[i] and dL/db[n] += delta[n]. v0.0.2 assigned here, which is why it could not express mini-batching.

Only the inner loop over inputs is vectorised. The outer loop over neurons carries a per-neuron delta, so the header notes it is left scalar.

void updateWeightsAndBiases(float learningRate)

Applies w -= learningRate * dL/dw and the same for biases. Plain SGD; no momentum, no weight decay. Does not zero the gradients afterwards — zeroGradients() is the caller's job at the top of the next batch.

Layer::print() was removed.

#class core::HiddenLayer

class HiddenLayer : public Layer, inheriting the constructors via using Layer::Layer.

void computeErrorSignal(int nextLayerWidth, const float* nextLayerDeltas, const float* nextLayerWeights)

Computes this layer's error signal from the next layer's:

delta[n] = ( sum over m of nextDeltas[m] * nextWeights[m][n] ) * f'(z[n])

The three arguments describe the layer after this one. Because the next layer's fan-in equals this layer's width, nextLayerWeights is indexed as nextWeights[m * this->layerWidth + n] — the stride is this->layerWidth, not this->inputSize.

#class core::OutputLayer

class OutputLayer : public Layer, inheriting the constructors.

void computeErrorSignal(const float* __restrict__ target)

Computes

delta[n] = dL/da[n] * f'(z[n])

by calling the loss derivative into scratchPad1 and the activation derivative into scratchPad2, then multiplying elementwise.

target must point to at least layerWidth floats. This is not checked.

Warning

The elementwise product cannot express softmax, whose Jacobian is a dense matrix. The fused softmaxCrossEntropy path in core/loss.h exists for exactly this, and v0.0.1 routes to it, but v0.0.3 does not: a softmax output layer still multiplies by an activation derivative that fetchAssociatedDerivative returns as null.

#class core::Network

Defined in core/network.h. All members are public.

#Data Members

#Member Functions

Network()

Zero hidden layers, one output slot, learning rate 0.01f, batch size 1. The layer vectors are left empty, so the slots must not be used.

Network(int numHiddenLayers, int numOutputLayers, float learningRate, int batches = 1)

Resizes both layer vectors. Every slot is a null pointer until filled.

void addHiddenLayer(int layerIndex, int width, int inputSize, ActivationFunction, initFunction, LossFunction)

Six arguments, down from eight. Allocates a HiddenLayer with new and stores it at layerIndex.

Warning

layerIndex is not bounds-checked. An out-of-range index writes past the end of the vector. Unchanged from v0.0.2.

void addOutputLayer(int layerIndex, int width, int inputSize, ActivationFunction, initFunction, LossFunction)

As above, for the output layer. Same lack of bounds checking.

void initializeAllWeightsAndBiases()

Calls initializeWeightsAndBiases() on every slot.

Warning

Dereferences every slot without a null check, so this crashes if any slot was never filled. Unchanged from v0.0.2.

void forwardPass(const float* inputs)

Runs each hidden layer in order, then each output layer, feeding each layer's activations to the next. Read the result from outputLayers[0]->neuronOutputsActivated.

void zeroAllGradients(), void scaleAllGradients(float scale), void applyGradients()

New in v0.0.3. Fan the corresponding Layer method out over every layer. applyGradients() uses the network's stored learningRate.

void accumulateGradients(std::vector<float>& targetOutputs)

Replaces backwardPass. Computes error signals and gradients for the output layer, then for the hidden layers from last to first. Does not apply the update — that is applyGradients(), called once per batch. The argument is a non-const reference but is not modified.

float computeLoss(const std::vector<float>& targetOutputs) const

Sums each output layer's configured loss over the corresponding slice of targetOutputs.

Fixed since v0.0.2

v0.0.2's computeLoss always summed squared error regardless of the configured loss, so training on binaryCrossEntropy and reading the loss from here reported a different quantity from the one being minimised. That was v0.0.2 known bug 1. v0.0.3 calls layer->lossFunction.

The comment above the return says the value is averaged over output layers; it is not, it is summed. With the usual single output layer the two agree.

float trainBatch(size_t startSample, int batchSize, std::vector<std::vector<float>>& inputs, std::vector<std::vector<float>>& targetOutputs)

New in v0.0.3. Zeroes the gradients, runs forward and accumulateGradients over batchSize samples starting at startSample, scales by 1 / batchSize, applies the update, and returns the mean loss over the batch. Indices are not bounds-checked against inputs.size().

void train(std::vector<std::vector<float>>& inputs, std::vector<std::vector<float>>& targetOutputs, int numEpochs, float* __restrict__ lossBuffer, int bufferSize, bool verbose = false)

Runs numEpochs passes, splitting each into batches of this->batches and clamping the last to the remainder. Writes the mean loss per epoch into lossBuffer[epoch].

Throws std::invalid_argument if the two containers differ in length, or if lossBuffer is null or shorter than numEpochs; returns immediately if inputs is empty. v0.0.2 checked the buffer inside the epoch loop, after the first epoch had already run; v0.0.3 checks up front.

verbose defaults to false and gates all printing. v0.0.2 printed every epoch unconditionally.

There is still no shuffling, so every batch holds the same samples every epoch.

~Network()

deletes every layer in both vectors.

printAllLayers(), runEpoch() and backwardPass() were removed. predict() is present but commented out — it called an undefined predictSample and returned a pointer into a destroyed temporary. There is currently no prediction API beyond forwardPass() followed by reading outputLayers[0]->neuronOutputsActivated.

#Activation Functions

Namespace core::functions::activations. All are elementwise except softmax. Each derivative takes the pre-activation z.

The In lookup? column is whether fetchAssociatedDerivative recognises the function. The three that do not yield a null derivative pointer, which is stored without complaint and called on the first backward pass.

softmax subtracts the row maximum before exponentiating, so it does not overflow. sigmoid does not branch on the sign of its argument and evaluates exp(-x) directly — v0.0.1 does branch.

Two cosmetic notes on the source. The gelu header comment reads 0.0f44715, a find-and-replace of 0.044715 that caught the comment as well as the code; the constant itself is right. And softmax's doc comment sits between its two parameters rather than above the function — legal, since comments are stripped before parsing, but the parameter list reads as split in half.

#Weight Initializers

Namespace core::functions::initializers. Unchanged from v0.0.2 apart from the scalar type.

float* zeros(float* pad, size_t size, size_t size2)

Sets size * size2 values to zero.

float* xavier(float* pad, size_t fanIn, size_t fanOut)

Normal with standard deviation sqrt(2 / (fanIn + fanOut)). Suited to sigmoid and tanh.

float* he(float* pad, size_t fanIn, size_t fanOut)

Normal with standard deviation sqrt(2 / fanIn). Suited to relu and its relatives.

The shared generator is seeded once from std::random_device and cannot be reseeded, so weight initialization is not reproducible — unlike v0.0.1, which takes a seed and records it. The header explains why these functions carry no SIMD pragmas: std::normal_distribution is stateful and cannot be vectorised.

#Loss Functions

Namespace core::functions::loss. The catalogue is unchanged from v0.0.2; only the scalar type differs. Every function takes target first, prediction second.

All six are covered by loss::fetchAssociatedDerivative, so the loss side of the lookup is complete.

#Weight Storage Layout

Unchanged from v0.0.2, in float. A layer's weights are one flat array of layerWidth * inputSize values, neuron-major:

neuronWeights[neuron * inputSize + input]

  neuron 0: [ w00 w01 w02 ... ]   neuron 1: [ w10 w11 w12 ... ]   ...
  |<------ inputSize ------>|

neuronWeightGradients uses the same indexing, so a weight and its gradient sit at the same offset in their respective buffers. Two strides are easy to confuse: Layer indexes its own weights with stride inputSize, while HiddenLayer::computeErrorSignal indexes the next layer's weights with stride this->layerWidth — which is that layer's fan-in.

At float, a 784 → 128 layer's weights occupy 392 KB rather than 784 KB, which is the whole argument for the change.

#Complete Example

#include <vector>
#include "core/network.h"

namespace act  = core::functions::activations;
namespace init = core::functions::initializers;
namespace loss = core::functions::loss;

int main()
{
    // 1 hidden layer, 1 output layer, lr = 0.3, batch size 3
    core::Network net(1, 1, 0.3f, 3);

    // no derivatives passed -- they are looked up from the forward function
    net.addHiddenLayer(0, /*width=*/4, /*inputSize=*/2,
                       act::sigmoid, init::xavier, loss::mse);
    net.addOutputLayer(0, /*width=*/1, /*inputSize=*/4,
                       act::sigmoid, init::xavier, loss::mse);

    net.initializeAllWeightsAndBiases();

    std::vector<std::vector<float>> X = {{0.1f, 0.2f}, {0.8f, 0.9f},
                                          {0.2f, 0.1f}, {0.9f, 0.8f}};
    std::vector<std::vector<float>> Y = {{0.0f}, {1.0f}, {0.0f}, {1.0f}};

    std::vector<float> lossHistory(1000);
    net.train(X, Y, 1000, lossHistory.data(),
              static_cast<int>(lossHistory.size()), /*verbose=*/true);

    net.forwardPass(X[1].data());
    float p = net.outputLayers[0]->neuronOutputsActivated[0];
    return 0;
}

Three things to watch. The output layer's inputSize must equal the hidden layer's width, and nothing checks it — note that the shipped main.cpp declares a hidden layer of width 2 and an output layer of inputSize 2, which happens to agree. Read predictions through forwardPass and neuronOutputsActivated, since predict() is commented out. And pass verbose=true if you want to see anything.

#Known Bugs and Limitations

Nine of the twelve issues this page listed for v0.0.3 are fixed. They are kept below, struck through with what was done, because several are the kind of defect that reads as correct code — and because the test suite now pins each one.

Still open

1. std::expf is not a standard name. It compiles under libc++ and is not portable; std::exp is the spelling that is. Four call sites in core/activations.h.
2. Only the first output layer back-propagates. accumulateGradients reads outputLayers[0] when propagating into the hidden stack, so with numOutputLayers > 1 the other heads’ gradients are silently discarded. The Python binding refuses to build more than one output layer for this reason; the C++ API still lets you.
3. Network is unsafe to copy. It owns raw pointers and has a destructor but no copy constructor or assignment operator, so a copy double-frees. The Python binding deletes both.

Fixed in v0.0.3

1. The library links only into one translation unit. Both fetchAssociatedDerivative functions are inline, and core/loss.h gained a #pragma once. See One Translation Unit Only.
2. A null derivative pointer is a crash, not a diagnostic. identity and gelu are in the lookup. softmax deliberately is not — its derivative is a Jacobian — and the Python layer rejects it at add_* time, pointing at identity + softmax_cross_entropy, which is the fused form that does work.
3. binaryCrossEntropy’s clamp is inert in float. The guard is now 1e-7f. This one was not theoretical: it took a real run to NaN at epoch 243, with the loss already down at 0.003. 1.0f - 1e-12f is 1.0f, so the upper clamp never bit; a sigmoid saturating to exactly 1.0 gave log(0) and a division by zero, and inf * 0 is NaN, which reaches every weight on the next update and never washes out.
4. No shuffling. train() takes shuffle and an optional seed, and reshuffles a sample-order vector each epoch. The order lives in the engine so that every caller — C++ or Python — gets the same batching and the same update rule.
5. addHiddenLayer and addOutputLayer do not bounds-check layerIndex. Both check before allocating and throw std::out_of_range. Replacing a filled slot now deletes the old layer instead of leaking it.
6. initializeAllWeightsAndBiases dereferences unfilled slots. It skips nulls, as the destructor already did.
7. The default Layer constructor leaves function pointers indeterminate. All five have = nullptr initialisers.
8. An array-section reduction on a non-reduction loop. Removed from softmaxCrossEntropyDerivative. It was not theoretical: with -fopenmp-simd the fused softmax test failed outright.
9. No seed control. initializers::setSeed(unsigned) reseeds the shared generator; gradiex.seed(n) exposes it. Weight initialisation is reproducible, which is what makes a run repeatable end to end.

#Not Supported

Optimizers other than plain SGD (no momentum, no Adam, no learning-rate schedule), L2 weight decay, dropout, train/validation split, early stopping, model save and load, a trainable softmax output, multi-class output, convolution, recurrence, attention, embeddings, skip connections, batch or layer normalization, sparse input, streaming datasets, GPU and threading.

Vectorisation is requested through #pragma omp simd rather than written with intrinsics, so it is present only when the build enables OpenMP — see SIMD Pragmas. Mini-batching is supported in v0.0.3, unlike v0.0.2.

Most of the list above was present in v0.0.1 and has not yet returned. See 0.0.3/TODO.md for measured performance work in priority order.