Coding a neural network

1. Know the tools before using them

We choose Python because it provides a short path from array mathematics to PyTorch. Python is a general-purpose programming language, not a special neural-network environment. The language does not make the model learn; the calculations in our program do.

A virtual environment is a directory, here called .venv, with a Python entry point and its own installed packages. It separates this project's dependencies from another project's dependencies. It is not a virtual machine and does not create a GPU.

The library called PyTorch is installed and imported under the name torch. NumPy uses numpy. Package names and product names are not always identical.

You may have seen Jupyter notebooks: documents containing text and executable code cells. They can be useful, but execution order can hide state. We will start with ordinary scripts that run from top to bottom. An editor and terminal are enough. CPU execution is deliberate; no notebook, cloud account, or GPU is required.

Make a small, reproducible working folder

Download the Python project, extract it, and open a terminal in the extracted coding_lesson_code directory—the folder containing experiment.json. Python 3.13 is the tested interpreter. If needed, obtain it from python.org.

-m asks that interpreter to run a module such as venv or pip. -r tells pip to read package versions from a requirements file. We call the environment's Python by its path, so activation is not required. If you use an editor's Run button, configure it to use that same interpreter.

Create a file named hello.py containing these Python statements:

import numpy as np
print(np.__version__)

Then launch it from the terminal:

If the version prints, your interpreter can import NumPy. This establishes that the environment works; it says nothing yet about a network. Installation normally needs internet access. Once packages are installed, this project's images and calculations are local.

If installation succeeds but you see ModuleNotFoundError, check which interpreter ran the file. Repeatedly installing into the wrong environment will not solve that mismatch.

2. Put data into arrays

We have an interpreter and NumPy. Before implementing a layer, we need one consistent answer to “what does each dimension mean?” Many apparent calculus bugs are actually shape or axis mistakes.

A NumPy ndarray is an N-dimensional numerical array. A vector, matrix, and stack of images are different shapes of this basic object. shape gives the size of each dimension; dtype describes the numerical type stored in its entries.

Our first image is a handwritten 0. It contains 784 brightness values, originally integers from 0 to 255. We flatten the 28 × 28 grid into one row and divide those intensities by 255, putting brightness between zero and one. A batch is a group of examples processed together. Stack B images as B rows:

Only W1, b1, W2, and b2 are trainable parameters. X and y are fixed data. Intermediate arrays, such as activations and predictions, are recalculated from those inputs and parameters; they are not additional independent parameters.

The comma in (B,) denotes a one-element Python tuple. This is a one-dimensional array, not a two-dimensional column (B, 1). The value y[0] = 0 identifies the first image's class. Its one-hot row has Y[0, 0] = 1 and zeros in the other nine positions. The class label and the one-hot target are two representations of the same known answer, not predictions.

We store weights as (output, input). Each row of W1 holds one hidden neuron's 784 weights. Some mathematical explanations store examples as columns and write W @ X. Here examples are rows and the equivalent operation is X @ W.T. This changes the representation, not the network.

Just enough NumPy syntax

import numpy as np       # give the imported module a short name
X.shape                 # a tuple, not a function call
X[0]                    # first image: shape (784,)
X[0:1]                  # a batch containing it: shape (1, 784)
W1.T                    # transpose this two-dimensional array
X @ W1.T                # matrix multiplication
X * X                   # element-by-element multiplication

Indices start at zero. The slice start:stop includes start but excludes stop. Keeping the batch dimension with X[0:1] lets the same code handle one or many images. Later, params['W1'] accesses a named array in a Python dictionary; the dictionary simply groups the four parameter arrays.

Broadcasting and axes are part of the algorithm

Adding b1, shape (16,), to the result (B, 16) uses broadcasting: the same bias vector is applied to every row. There is not a separate learned bias for every image. NumPy aligns trailing dimensions; corresponding sizes must be equal or one must be 1.

An axis is a dimension. Axis 0 is the batch; axis 1 is the features or classes. A reduction such as sum(axis=0) removes that axis by summing over it. keepdims=True retains a length-one dimension, which can be important when broadcasting the result back. We will use this in softmax and in bias gradients.

The supplied load_data helper decodes a compact local file and returns X, y in these formats. File decoding is not a hidden network operation. The file contains 200 real MNIST training images and 100 real test images, with original split indices—not a substitute synthetic dataset.

3. Write the forward calculation

We have chosen batch-first arrays and know which dimensions must match. NumPy can perform the arithmetic, but we still decide which operations make up the network.

“From scratch” here means we specify the model, loss, derivatives, and update rule ourselves. NumPy supplies array storage and efficient arithmetic. Reimplementing an optimized matrix-multiplication library would be a different lesson.

One entry, with actual numbers

For hidden neuron 0, multiply each pixel by its corresponding weight, add all 784 products, and add its bias. Two terms from our first image are:

A black pixel contributes zero to this sum. These are only two of the 784 terms, not a reduced two-input network. Including all products and the zero bias gives a weighted sum of -0.4108369. ReLU returns the larger of that value and zero, so this neuron's output is 0 for this image.

Here is the weighted sum in ordinary Python. def defines a function, indentation groups its body, and return sends a value back to its caller.

def one_neuron(x, weights, bias):
    total = float(bias)
    for i in range(len(x)):
        total += float(x[i]) * float(weights[i])
    return total

This is the actual function in network_numpy.py. Calling it with X[0], params['W1'][j], and params['b1'][j] calculates one entry of Z1. The expression X @ W1.T + b1 performs the collection of dot products for every image and hidden neuron. Vectorizing the arithmetic does not create a different learning algorithm.

Initialize once, then reuse the parameters

Weights begin as small random floating-point values. Biases begin at zero. Different initial weights keep the hidden units from being identical copies. The scale factors and seed are explicit settings in experiment.json; they are initialization choices, not learned values.

def initialize(cfg):
    D, H, K = cfg['input_size'], cfg['hidden_size'], cfg['classes']
    rng = np.random.default_rng(cfg['seed'])
    return {
        'W1': (rng.standard_normal((H, D)) * np.sqrt(cfg['hidden_variance'] / D)).astype(cfg['dtype']),
        'b1': np.zeros(H, dtype=cfg['dtype']),
        'W2': (rng.standard_normal((K, H)) * np.sqrt(cfg['output_variance'] / H)).astype(cfg['dtype']),
        'b2': np.zeros(K, dtype=cfg['dtype']),
    }

Initialization happens before the training loop. Reinitializing at every step would erase the previous step's progress. A forward pass reuses the current parameter values; it does not draw a new random model.

The full forward function

def forward(params, X):
    Z1 = X @ params['W1'].T + params['b1']
    A1 = np.maximum(Z1, 0)
    logits = A1 @ params['W2'].T + params['b2']
    shifted = logits - logits.max(axis=1, keepdims=True)
    logP = shifted - np.log(np.exp(shifted).sum(axis=1, keepdims=True))
    P = np.exp(logP)
    return {'X': X, 'Z1': Z1, 'A1': A1, 'logits': logits, 'P': P, 'logP': logP}

The first two lines implement the hidden layer. ReLU is elementwise. Without a hidden nonlinearity, the two affine layers could be collapsed into one affine mapping to the output scores.

logits contains ten unconstrained scores per image, not probabilities. Softmax normalizes each row. Subtracting the row's largest score does not change its softmax probabilities, but avoids exponentiating a large positive number. The row maximum and row sum retain shape (B, 1) so they broadcast along the class dimension, not across unrelated images.

logP stores log-probabilities directly. This is useful because computing the logarithm of a probability that has underflowed to zero would fail numerically. The stable log-probability calculation avoids that problem for our cross-entropy loss.

The returned dictionary is a cache. Backpropagation will need X, the ReLU inputs Z1, and hidden activations A1 from this exact forward pass. Saving them is bookkeeping, not another learned layer. Notice that y does not appear: predictions must be calculated without knowing the answer.

For our starting image, the probability of digit 0 is 0.1351491. It happens to be the largest of the ten probabilities, so this randomly initialized model already guesses this image's class correctly. That coincidence is not evidence that training has occurred, and it does not make the loss zero.

4. Write the derivatives

Forward propagation produced predictions and a cache. It did not measure their mismatch with the labels or change a parameter. We now add a scalar objective and its derivatives.

Begin with one image and one output bias

For the first image, the label is 0 and the model assigns it probability 0.1351491. Cross entropy for this one-hot target is the negative natural logarithm of that probability:

This is a scalar loss. It is different from the derivative telling us how the loss changes when a particular value changes. For softmax followed by this cross entropy, the derivative with respect to logit 0 is:

For one image, increasing output bias 0 increases logit 0 by the same amount. Its gradient is therefore also -0.8648509. This is one component of the whole parameter gradient, not a new one-parameter model. We will use this component to inspect an update in section 5.

For a batch, we average the per-image losses. The initial logit derivatives inherit that average. In code, dZ2 means the derivative of the scalar loss with respect to the output-logit array; it has shape (B, 10). In the formulas, $D_2$ denotes this same array, and $G_P$ denotes the gradient with respect to the probability array.

Two loss choices, one backward interface

With one-hot target array Y, the mean cross-entropy loss and its logit signal are:

The compact difference in the second expression is a derivative, not a linear cost function. The scalar objective still uses logarithms.

For MSE on probabilities, average the squared differences over B × K entries, with K = 10. Differentiating each square supplies the factor 2; the mean supplies division by B × K. Then pass that probability derivative through softmax. Softmax couples the classes, so the logit gradient is not just ten independent sigmoid derivatives:

The row-wise dot-product subtraction in the MSE branch below accounts for that coupling. Both branches return the same kind of object: a scalar loss and the gradient with respect to logits.

def loss_signal(cache, y, kind):
    P = cache['P']
    B, K = P.shape
    Y = np.eye(K, dtype=P.dtype)[y]
    if kind == 'ce':
        loss = -cache['logP'][np.arange(B), y].mean()
        dZ2 = (P - Y) / B
    elif kind == 'mse':
        error = P - Y
        loss = (error ** 2).mean()
        dP = 2 * error / (B * K)
        dZ2 = P * (dP - (dP * P).sum(axis=1, keepdims=True))
    else:
        raise ValueError('Loss must be ce or mse')
    return float(loss), dZ2

np.eye(K)[y] selects rows of an identity matrix to construct one-hot targets. The paired indices np.arange(B), y select the true class's log-probability in each row. ** 2 squares each error. The code returns an ordinary floating-point loss because NumPy is not tracking an automatic-differentiation graph here.

This choice affects the later derivative previews and resets the step walkthrough. It does not edit the terminal project's experiment.json. There, loss is the setting to change when you want to run the other objective.

Continue through the two layers

Consider the weight linking hidden unit 1 to output logit 0. For our fixed one-image CE example, its input activation is 0.1643927. A small increase in this weight changes the logit in proportion to that activation, so the chain rule multiplies the logit sensitivity by the activation:

For a batch, sum one such product per image. Matrix multiplication collects the products for every second-layer weight: (10, B) @ (B, 16) produces (10, 16), exactly W2's shape. Biases use the corresponding sum without an activation multiplier, because each bias enters its weighted sum with coefficient one.

def backward(params, cache, dZ2):
    dW2 = dZ2.T @ cache['A1']
    db2 = dZ2.sum(axis=0)
    dA1 = dZ2 @ params['W2']
    dZ1 = dA1 * (cache['Z1'] > 0)
    dW1 = dZ1.T @ cache['X']
    db1 = dZ1.sum(axis=0)
    return {'W1': dW1, 'b1': db1, 'W2': dW2, 'b2': db2}

Each output bias was reused across all B rows, so its derivative adds those image contributions. Summing axis 0 leaves one entry per class. To reach the hidden activations, (B, 10) @ (10, 16) produces (B, 16). Then the elementwise ReLU mask passes a derivative where the cached pre-activation was positive and blocks it otherwise. We choose derivative zero at exactly zero.

Finally, dZ1.T @ X applies the same dense-layer rule to the first layer. The gradient has shape (16, 784), matching W1. We did not derive a new kind of calculus for “layer 1” versus “layer 2”; the operations repeat with different arrays.

There is no second division by B. The loss's averaging factor is already in dZ2. Later products and sums carry it through the graph. Dividing the weight gradients again would shrink them by an unintended extra factor.

5. Assemble the learning step

Forward propagation gave us predictions. Backpropagation gave us derivatives. We now add the operation that changes the stored parameters, using gradients calculated from the old forward pass.

Plain stochastic gradient descent, or SGD, subtracts the learning rate times each parameter's gradient. For the one-image CE example, the inspected bias starts at zero and its derivative is -0.8648509. With learning rate 0.2:

The learning rate controls the step size; backpropagation does not choose it. An excessively large step can raise the loss even when the derivative is correct.

In network_numpy.py, PARAMS lists the four dictionary keys: W1, b1, W2, and b2. The update modifies their arrays:

def step(params, grads, learning_rate):
    for name in PARAMS:
        params[name] -= learning_rate * grads[name]

Wait until all derivatives have been calculated. The first-layer backward computation still needs the old W2. Updating it early would mix old and new parameters instead of differentiating one consistent forward pass.

The actual step function in train_numpy.py joins these pieces:

def train_step(params, X_batch, y_batch, cfg):
    cache = forward(params, X_batch)
    loss, dZ2 = loss_signal(cache, y_batch, cfg['loss'])
    grads = backward(params, cache, dZ2)
    step(params, grads, cfg['learning_rate'])
    return loss

The returned loss is the value before the update. Its Python variable does not recompute itself when parameters change. To measure the new cost, run forward and score the new predictions.

Repeat the step without reinitializing

With 200 training images and a batch size of 8, one pass supplies 25 updates. That full pass is an epoch. Eight epochs give 200 updates, not eight. If the size is not divisible by B, the last batch is smaller; our loss uses its actual size.

for epoch in range(cfg['epochs']):
    order = rng.permutation(len(X)) if cfg['shuffle'] else np.arange(len(X))
    for start in range(0, len(X), B):
        indices = order[start:start + B]
        train_step(params, X[indices], y[indices], cfg)
        updates += 1

The same indices select X and y, preserving each image–label pair. Parameters are initialized before these loops and reused after each update. The full file contains the setup and imports; the extracts show the operations we are studying, not independent scripts to paste without their definitions.

Before using the entire dataset, repeatedly train on one fixed small batch. This deliberate overfitting test checks whether the pipeline can learn a tiny problem. It is not a generalization test:

In the checked default CE run, both implementations take the same eight-image batch from loss about 2.28867 to 0.001583 in 300 updates, classifying all eight correctly. These are measured results for these settings, not promised outcomes for arbitrary data.

If a small batch cannot be learned, inspect shapes, normalization, label pairing, the learning rate, and gradients before increasing the model's size.

6. Let autograd calculate the derivatives

The NumPy implementation exposed every derivative. We can replace that manual differentiation without immediately hiding the forward operations inside layer objects.

A PyTorch tensor is an array-like numerical object with a shape, dtype, and device. The file starts with import torch. Here we use CPU float64 tensors to match NumPy closely. Many training projects use float32; GPU execution and lower precision are separate later topics.

For a parameter marked requires_grad=True, PyTorch records relevant tensor operations during forward computation. This record is a computation graph. Each supported operation has a local backward rule. Calling loss.backward() combines those rules in reverse with the chain rule and stores derivatives in the parameters' .grad fields.

Autograd is not an AI guessing derivatives. It is not perturbing each of 12,730 parameters with finite differences either. It executes implemented differentiation rules for the actual operations in the graph. You still choose the forward program and the scalar objective.

Keep the explicit matrix operations

The one-step bridge in torch_tensors.py keeps our matrix expressions. objective is the PyTorch loss helper explained in section 7; on the CE route it calls F.cross_entropy(logits, y_tensor). For now, its needed contract is simple: it returns one differentiable scalar tensor.

def tensor_step(params, X, y, cfg):
    weights = {name: torch.tensor(value, requires_grad=True)
               for name, value in params.items()}
    X_tensor, y_tensor = torch.from_numpy(X), torch.from_numpy(y)
    Z1 = X_tensor @ weights['W1'].T + weights['b1']
    A1 = torch.relu(Z1)
    logits = A1 @ weights['W2'].T + weights['b2']
    loss = objective(logits, y_tensor, cfg['loss'])
    loss.backward()
    gradients = {name: p.grad.detach().numpy().copy() for name, p in weights.items()}
    with torch.no_grad():
        for p in weights.values():
            p -= cfg['learning_rate'] * p.grad
    after = {name: p.detach().numpy().copy() for name, p in weights.items()}
    return loss.item(), gradients, after

The opening dictionary comprehension is a compact loop: for each name–array pair, create a tensor under the same name. torch.tensor makes a copy of each parameter array. torch.from_numpy can share the input array's memory; we do not modify those input data.

torch.relu provides a PyTorch operation with a known backward rule. Keep NumPy conversions out of the differentiable forward path. The NumPy library is not automatically added to PyTorch's graph merely because the values look similar.

loss.backward() replaces the manual starting derivative and backward function. It does not change the weights. We still implement the subtraction ourselves here. torch.no_grad() keeps that housekeeping update from becoming part of the graph being differentiated.

The detach().numpy().copy() expressions export independent arrays for verification, not training. loss.item() converts a scalar tensor into an ordinary Python number for reporting, after backward has used the tensor.

This demonstration creates fresh parameter tensors for one step. A repeated training loop should retain its tensors and clear old gradient buffers before accumulating the next batch's derivatives.

Do not repeatedly reuse an already-consumed backward graph to implement the ordinary training loop. Compute a fresh forward pass after each update; saved backward intermediates are normally released after backward. This loop does not require retain_graph=True.

7. Use layers and an optimizer

Tensors and autograd already implement the whole algorithm. PyTorch's higher-level building blocks now remove repetitive storage and update code—not the need to understand the calculation.

A layer owns its parameters

from torch import nn gives us PyTorch's neural-network module. nn.Linear(D, H) owns weights shaped (H, D) and a bias shaped (H,). Calling it computes X @ weight.T + bias. nn.ReLU() performs the elementwise gate and has no trainable parameters.

def build_model(cfg):
    D, H, K = cfg['input_size'], cfg['hidden_size'], cfg['classes']
    model = nn.Sequential(
        nn.Linear(D, H),
        nn.ReLU(),
        nn.Linear(H, K),
    )
    return model.to(device=cfg['device'], dtype=getattr(torch, cfg['dtype']))

Sequential runs these modules in order. model(X) returns logits. model.parameters() exposes the registered trainable tensors to the optimizer. The final .to selects the configured device and dtype; getattr(torch, 'float64') means torch.float64.

No custom model class is necessary for this stack. A more complicated forward computation can use a subclass of nn.Module, but that is not a prerequisite for the present network.

The loss API has a contract

Why is there no Softmax module at the end? The common cross-entropy API expects raw logits and combines stable normalization with the logarithmic loss. For display, use logits.softmax(dim=1). To pick the winning class, logits.argmax(dim=1) already has the same result.

In network_torch.py, from torch.nn import functional as F gives the functional API a short name. F.cross_entropy(logits, y) provides the same standard calculation as calling an nn.CrossEntropyLoss() object. Our CE case uses (B, 10) logits and (B,) integer torch.long class indices. Probability targets are also supported by that API, but are not our CE representation here.

def objective(logits, y, kind):
    if kind == 'ce':
        return F.cross_entropy(logits, y, reduction='mean')
    if kind == 'mse':
        P = logits.softmax(dim=1)
        Y = F.one_hot(y, num_classes=logits.shape[1]).to(dtype=P.dtype)
        return F.mse_loss(P, Y, reduction='mean')
    raise ValueError('Loss must be ce or mse')

MSE here is deliberately on probabilities: apply softmax and construct floating-point one-hot targets. Mean reduction averages all B × 10 squared errors, matching NumPy. Applying MSE directly to the logits would define a different objective.

Give the optimizer its own job

optimizer = torch.optim.SGD(model.parameters(), lr=cfg['learning_rate'])

Then the actual function from train_torch.py performs our familiar sequence:

def train_step(model, optimizer, X_batch, y_batch, kind):
    optimizer.zero_grad(set_to_none=True)
    logits = model(X_batch)
    loss = objective(logits, y_batch, kind)
    loss.backward()
    optimizer.step()
    return loss.item()

zero_grad(set_to_none=True) clears previous gradient buffers by setting them to None. It does not zero parameters. backward() calculates and accumulates derivatives. step() applies the specified plain SGD update. The optimizer does not read labels or compute the loss for us, and it does not automatically clear gradients afterward.

Compare actual numbers, not matching seed labels

nn.Linear has its own default initialization. A shared seed does not impose identical random draws across NumPy and PyTorch. Our comparison therefore copies the NumPy arrays into the layers:

def copy_parameters(model, params):
    with torch.no_grad():
        model[0].weight.copy_(torch.from_numpy(params['W1']))
        model[0].bias.copy_(torch.from_numpy(params['b1']))
        model[2].weight.copy_(torch.from_numpy(params['W2']))
        model[2].bias.copy_(torch.from_numpy(params['b2']))

model[0] is the first Linear layer; model[2] is the second. copy_ writes into their existing storage; the underscore denotes an in-place operation. This copy is initialization, not learning. The models now start from the same actual values.

8. Run, verify, and extend the program

We can name each library's responsibility. The final step is not trusting the shortest implementation: it is checking that each implementation calculates the intended thing.

Use the complete project, not disconnected fragments

Download the complete Python project. Work in its extracted copy so your experiments do not alter this lesson's fixed browser examples. The command format below follows your earlier OS selection.

experiment.json           # architecture, objective, and training settings
requirements-numpy.txt    # install NumPy first
requirements-torch.txt    # add PyTorch afterward
data/mnist.json           # local real MNIST images
python/
  common.py               # configuration and file decoding
  network_numpy.py        # forward, loss, backward, and SGD
  train_numpy.py          # complete NumPy training program
  torch_tensors.py        # explicit tensor math plus autograd
  network_torch.py        # the same model as layers
  train_torch.py          # complete PyTorch training program
  verify.py               # controlled numerical comparisons

A useful first implementation sequence is:

1. Run the supplied NumPy overfit check once to verify the environment.
2. In a separate file, type your own forward calculation, then the loss and backward functions. Reuse the data loader. Check shapes after each operation.
3. Check selected derivatives numerically before trusting a long run.
4. Read the tensor-only bridge and identify exactly what backward work disappeared.
5. Compare every parameter gradient with PyTorch, then compare one update.
6. Only then run multiple batches and evaluate unseen images.

For one selected parameter, a centered finite difference compares the analytic derivative with two nearby losses:

Use copies or restore the parameter afterward, and avoid crossing a ReLU kink. This is an approximation for checking a derivative, not the algorithm we use to train all parameters. The checker warns if its perturbation crosses a gate boundary.

These programs were executed with Python 3.13.5, NumPy 2.4.4, and PyTorch 2.14.0, on an Apple Silicon CPU using float64. Windows command formatting is provided, but Windows and Linux execution were not directly tested. The Python project needs no notebook, Node.js, or running webpage.

Let a DataLoader deliver batches

A TensorDataset pairs input and target tensors along their first dimension. A DataLoader iterates over those pairs in batches and can shuffle them. It delivers data; it does not predict, differentiate, or train.

dataset = TensorDataset(X, y)
generator = torch.Generator().manual_seed(cfg['shuffle_seed'])
loader = DataLoader(dataset, batch_size=B, shuffle=cfg['shuffle'], generator=generator)
model.train()
updates = 0
for epoch in range(cfg['epochs']):
    for X_batch, y_batch in loader:
        train_step(model, optimizer, X_batch, y_batch, cfg['loss'])
        updates += 1

The two full-data programs use their libraries' own seeded shuffles. Those orders need not match across libraries, so final training trajectories can differ despite matching same-batch gradients. To compare a whole trajectory, also supply the exact same batch sequence.

Keep prediction separate from learning

For held-out images, run forward and score the results without backward or optimizer updates. The supplied PyTorch evaluation helper makes that boundary explicit:

def scores(model, X, y, kind):
    model.eval()
    with torch.no_grad():
        logits = model(X)
        loss = objective(logits, y, kind).item()
        accuracy = (logits.argmax(dim=1) == y).double().mean().item()
    return loss, accuracy

model.eval() selects evaluation behavior for modules such as dropout and batch normalization. torch.no_grad() disables gradient tracking for the enclosed calculation. They are different things. Our Linear/ReLU model has no mode-dependent layers, but the distinction matters when extending it. Conversely, model.train() selects a mode; calling it does not train the model by itself.

Our 200-training / 100-test image subset is a reproducible exercise, not a full-MNIST benchmark. Repeatedly using held-out scores to choose settings makes that set validation data rather than an untouched final test. The scripts print results and exit; trained weights are not saved.

Synthesis: what did each tool take over?

You do not need another neural-network framework immediately. First explain the PyTorch loop in terms of the manual program. Full MNIST, a convolutional network, checkpoints, and GPU execution are later extensions—not prerequisites for understanding this calculation.

Sources and execution notes

The following official sources were consulted to check the explanation and API contracts. They are optional further reading; the core lesson is local. Research checked on 21 September 2026.

• Python Packaging User Guide: pip and virtual environments — interpreters, packages, requirements files, and terminal commands.
• Python 3.13: venv — explicit interpreter paths work without activation.
• NumPy: the absolute basics and broadcasting — arrays, shapes, indexing, axes, and compatible dimensions.
• PyTorch: tensors — array-like operations, dtypes, devices, and memory sharing.
• PyTorch: automatic differentiation — graphs, backward, accumulation, and detaching.
• PyTorch 2.14: Linear, CrossEntropyLoss, and MSELoss — parameter layout and objective contracts.
• PyTorch: optimizing parameters and datasets and loaders — optimization and batch iteration.
• Autograd mechanics and reproducibility — evaluation versus no-grad and limits of matching seeds.

The images are genuine MNIST by Yann LeCun, Corinna Cortes, and Christopher J. C. Burges, retained from the Keras-hosted archive. Source checksums and original split indices are retained. No Fashion-MNIST images or synthetic substitutes are used.

See the local README, validation notes, and this complete Markdown manuscript. Numerical tests are evidence about the implementation, not a guarantee that every reader will find every explanation sufficient.