Nothing about a neural network is new to you. You already fitted one in Part 10. Logistic regression takes a weighted sum of your features, adds a bias, squashes the result through a sigmoid, and reports a probability. That is a neural network with exactly zero hidden layers, and every piece of vocabulary that gets thrown at deep learning maps onto something in it. Weights are the coefficients. Bias is the intercept. Activation is the sigmoid. Loss is log loss. Training is gradient descent.

Everything else in this part is one idea repeated: put a layer in the middle, make it nonlinear, and push gradients back through it. Once that clicks, PyTorch in Part 17 becomes a matter of syntax rather than concepts, which is exactly how it should feel.

Logistic regression as a network with no hidden layer

Where our project stands: last part we tuned a gradient boosting model on the churn feature matrix and landed on an honest test average precision after accounting for search optimism. This part keeps that exact matrix and asks a different question, which is whether a network can beat it, and what the network costs to run.

I am modelling churn on a 20,000 row matrix with a 3.4 percent positive rate, standing in for the public Telco Customer Churn file this series works from. Average precision remains the scoring metric, for the reasons argued in Part 14. Here is our baseline, and notice the shape of the fitted coefficients.

# tested against: Python 3.10.12, numpy 2.2.6, scipy 1.15.3,
# scikit-learn 1.7.2
import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.metrics import average_precision_score, roc_auc_score

X, y = make_classification(n_samples=20000, n_features=20, n_informative=6,
                           n_redundant=2, weights=[0.97, 0.03], flip_y=0.01,
                           class_sep=0.8, random_state=42)

X_tr, X_te, y_tr, y_te = train_test_split(
    X, y, test_size=0.25, stratify=y, random_state=42)
print('shapes', X_tr.shape, X_te.shape, 'pos rate', round(y_tr.mean(), 4))

lr = make_pipeline(StandardScaler(), LogisticRegression(max_iter=1000))
lr.fit(X_tr, y_tr)
p = lr.predict_proba(X_te)[:, 1]
print('logreg test AP  %.4f  AUC %.4f'
      % (average_precision_score(y_te, p), roc_auc_score(y_te, p)))
print('weights', lr.named_steps['logisticregression'].coef_.shape,
      'bias', lr.named_steps['logisticregression'].intercept_.shape)

shapes (15000, 20) (5000, 20) pos rate 0.0344
logreg test AP  0.3880  AUC 0.8142
weights (1, 20) bias (1,)

Twenty weights, one bias, average precision of 0.3880. Write it as a matrix product and the network reading falls out immediately: a 20 by 1 weight matrix maps our 20 inputs to a single output unit, a bias is added, a sigmoid squashes it. Deep learning papers would call that a single dense layer with sigmoid activation. Same object, different vocabulary.

Gains from adding a hidden layer

Our logistic model can only draw a straight boundary through feature space. Customers who churn when tenure is short and charges are high, but not when either condition holds alone, sit outside what a weighted sum can express. In Part 6 we solved that by hand, building interaction terms one at a time. A hidden layer does the same job, except it learns which interactions matter instead of waiting for you to guess.

Insert a layer of 16 units between input and output and you get a 20 by 16 matrix, then a 16 by 1 matrix. Multiply those two together, though, and you land back at a single 20 by 1 matrix. Stacking linear layers gains nothing at all. Nonlinearity between them is what stops the collapse, and ReLU, which is just max of zero and the input, is the one almost everybody reaches for because its gradient is either zero or one and never vanishes on the positive side.

Sigmoid still sits at the output because we want a probability between zero and one. Hidden layers use ReLU because they are computing intermediate representations, not probabilities, and squashing those into a narrow range throttles the gradient. Mixing the two up is a common beginner error and produces a network that trains agonisingly slowly for no visible reason.

Backpropagation in thirty lines of numpy

Backpropagation gets described as if it were an algorithm you need to study. It is the chain rule, computed once per layer, moving right to left. Write it out for a two layer network and the entire thing fits in a screen. Below is a complete trainable network with no framework anywhere in it.

sc = StandardScaler().fit(X_tr)
Xtr, Xte = sc.transform(X_tr), sc.transform(X_te)
ytr = y_tr.reshape(-1, 1).astype(float)

rng = np.random.default_rng(0)
n_in, n_hid = Xtr.shape[1], 16
# He initialisation, scaled by fan in, keeps early activations alive
W1 = rng.normal(0, np.sqrt(2.0 / n_in), size=(n_in, n_hid))
b1 = np.zeros((1, n_hid))
W2 = rng.normal(0, np.sqrt(2.0 / n_hid), size=(n_hid, 1))
b2 = np.zeros((1, 1))

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

lr, n = 0.5, Xtr.shape[0]
for epoch in range(1, 401):
    Z1 = Xtr @ W1 + b1          # forward
    A1 = np.maximum(0.0, Z1)
    Z2 = A1 @ W2 + b2
    A2 = sigmoid(Z2)
    loss = -np.mean(ytr * np.log(A2 + 1e-9)
                    + (1 - ytr) * np.log(1 - A2 + 1e-9))

    dZ2 = (A2 - ytr) / n        # backward, right to left
    dW2 = A1.T @ dZ2
    db2 = dZ2.sum(0, keepdims=True)
    dA1 = dZ2 @ W2.T
    dZ1 = dA1 * (Z1 > 0)        # ReLU derivative
    dW1 = Xtr.T @ dZ1
    db1 = dZ1.sum(0, keepdims=True)

    W2 -= lr * dW2; b2 -= lr * db2
    W1 -= lr * dW1; b1 -= lr * db1
    if epoch in (1, 50, 100, 200, 400):
        print('epoch %3d  loss %.4f' % (epoch, loss))

A1t = np.maximum(0.0, Xte @ W1 + b1)
pt = sigmoid(A1t @ W2 + b2).ravel()
print('scratch net test AP %.4f' % average_precision_score(y_te, pt))

epoch   1  loss 0.6456
epoch  50  loss 0.1493
epoch 100  loss 0.1354
epoch 200  loss 0.1200
epoch 400  loss 0.1045
scratch net test AP 0.3461

Read the backward block once more, slowly, because three lines carry all the meaning. Line one, dZ2 is the prediction minus the truth, which is the derivative of log loss composed with sigmoid, and it comes out that clean only because those two are paired. Line two, dA1 pushes that error back through W2 by transposing it. Line three, dZ1 masks the error wherever the ReLU output was zero, since a dead unit contributed nothing and therefore deserves no gradient. Every deep learning framework does exactly this, generalised and on a GPU.

Do not miss the last line of output. Loss dropped by 84 percent and test average precision came in at 0.3461, which is worse than the 0.3880 logistic regression managed in one second with twenty parameters. Falling loss is not the same as a better model, particularly on an imbalanced target where a network can drive loss down by getting the easy negatives even more confidently right. Plain gradient descent with a fixed learning rate and no minibatching is also simply a weak optimiser, which is what Adam fixes.

Feature scaling, and a four times collapse in average precision

Real churn tables have columns on wildly different scales. Tenure runs zero to 72 months. Monthly charges run about 18 to 118. Total charges run zero to 8,000 or so. Tree models could not care less, since they only compare a value against a threshold. Networks care enormously, because a gradient step of a fixed size means something completely different for a column with a standard deviation of 3,177 than for one with a standard deviation of 1.

To measure it, I rescaled three of the twenty columns to those realistic ranges and fitted the same MLPClassifier twice, once bare and once inside a pipeline with StandardScaler in front of it.

from sklearn.neural_network import MLPClassifier
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.metrics import accuracy_score

# give three columns realistic ranges, like tenure and TotalCharges
scales = np.array([1,1,72,1,1,1,3200,1,1,1,118,1,1,1,1,1,1,1,1,1],
                  dtype=float)
Xs = X * scales
X_tr, X_te, y_tr, y_te = train_test_split(
    Xs, y, test_size=0.25, stratify=y, random_state=42)
print('majority class accuracy %.4f' % (1 - y_te.mean()))

net = MLPClassifier(hidden_layer_sizes=(128, 64), max_iter=300,
                    early_stopping=True, random_state=42)

bare = net.fit(X_tr, y_tr)                      # no scaling
print('unscaled AP %.4f  acc %.4f'
      % (average_precision_score(y_te, bare.predict_proba(X_te)[:, 1]),
         accuracy_score(y_te, bare.predict(X_te))))

pipe = make_pipeline(StandardScaler(),
                     MLPClassifier(hidden_layer_sizes=(128, 64), max_iter=300,
                                   early_stopping=True, n_iter_no_change=10,
                                   random_state=42)).fit(X_tr, y_tr)
print('scaled   AP %.4f  acc %.4f'
      % (average_precision_score(y_te, pipe.predict_proba(X_te)[:, 1]),
         accuracy_score(y_te, pipe.predict(X_te))))

majority class accuracy 0.9656
unscaled AP 0.1062  acc 0.9660
scaled   AP 0.5445  acc 0.9758

Two failures wait for you here and both are worth meeting on purpose. First, missing values. Gradient boosting in scikit-learn handles NaN natively, so a column you never imputed passes silently. Networks refuse outright.

>>> Xn = X_tr.copy(); Xn[5, 3] = np.nan
>>> make_pipeline(StandardScaler(), MLPClassifier()).fit(Xn, y_tr)
ValueError: Input X contains NaN.
MLPClassifier does not accept missing values encoded as NaN natively.
For supervised learning, you might want to consider
sklearn.ensemble.HistGradientBoostingClassifier and Regressor which
accept missing values encoded as NaNs natively.

# and with max_iter set too low
ConvergenceWarning: Stochastic Optimizer: Maximum iterations (8)
reached and the optimization hasn't converged yet.

ConvergenceWarning is the dangerous one, because a warning does not stop a script and a half trained network still returns predictions that look plausible. Treat it as an error in any pipeline that runs unattended. Handling of missing values, incidentally, is the same discipline covered in cleaning messy data, and moving from trees to networks removes your ability to skip it.

Initialisation, and units that die before they learn

One line in the numpy code above deserves more attention than it usually gets. Weights were drawn from a normal distribution scaled by the square root of two over the number of inputs, which is He initialisation, named for Kaiming He and his co authors who proposed it for exactly this pairing of ReLU with deep stacks. Set every weight to zero instead and the network never trains at all, because all 16 hidden units compute the identical value, receive the identical gradient, and stay identical forever. Symmetry has to be broken before anything can be learned.

Draw weights too large and you get a subtler failure. A hidden unit whose pre activation lands well below zero for every single row in your training set outputs zero for every row, and because the ReLU derivative is zero there, it receives no gradient and cannot recover. That unit is dead for the rest of training. Capacity you paid for and cannot use. Scaling initialisation by fan in keeps early pre activations centred near zero, which keeps roughly half the units firing and all of them reachable.

Diagnosing this takes one line. After a few epochs, compute the fraction of hidden activations that are exactly zero across a batch. Around 40 to 60 percent is healthy sparsity. Anything above about 90 percent means most of your layer has died and no amount of extra epochs will revive it; drop the learning rate or reinitialise. I have seen a network sit at 0.5 average precision for an hour of training because 118 of its 128 hidden units were producing zeros, and the loss curve gave no hint whatsoever, since the surviving ten units were quite happily reducing it. Frameworks default to sensible initialisation now, which is why this is rarely a problem you cause deliberately, but writing your own layer in Part 17 puts it back within reach.

Sizing a network for tabular data

Choosing hidden layer sizes attracts more superstition than any other decision in modelling. I swept five architectures on the identical split, holding everything else fixed, and the results argue against most of the folklore.

Three things stand out. Going from 177 parameters to 11,009 bought 0.0453 average precision, which is real. Going from 11,009 to 46,593 lost 0.0155 while taking 73 percent longer, so depth on a twenty column table is a straightforward waste. And a single wide layer at (32,) beat the deeper (32, 16) outright, which surprises people who assume more layers means more capacity. Twenty features simply do not contain enough structure to justify a hierarchy.

My starting point for tabular work is one hidden layer sized between two and eight times the feature count, early stopping switched on, and Adam left at its defaults. Search width before you search depth. Reach for a second layer only when a wider single layer has clearly plateaued, which on churn data it usually has by around 128 units.

Reading a loss curve without fooling yourself

MLPClassifier exposes loss_curve_ after fitting, and with early_stopping switched on it also exposes validation_scores_. Plot both together. Doing so on our (32, 16) run produced the single most instructive chart in this part.

Look at what those two panels say together. Training loss halved between epoch 16 and epoch 58. Validation accuracy over the same stretch wobbled between 0.9760 and 0.9787 and ended lower than its peak. All that late training bought nothing; it was memorisation. Early stopping caught it, which is exactly why I switch it on by default rather than tuning max_iter by hand.

Note also that scikit-learn scores validation with accuracy, not with your chosen metric, and on a 3.4 percent positive rate that makes early stopping considerably less sensitive than you would like. Best validation accuracy of 0.9787 sits 1.3 points above predicting no churn for everybody. If you want early stopping driven by average precision, you have to leave MLPClassifier and write your own loop, which is one of the better arguments for moving to PyTorch in Part 17.

Network against gradient boosting on the same churn split

Here is the comparison that decides whether any of this ships. Four models, one split, one metric.

My verdict, and it is not the exciting one. Gradient boosting is what I would put in front of a stakeholder for churn on a tabular feature matrix. It cost 0.3 seconds against 2.3, needed no scaler, swallowed missing values without complaint, and gave up 0.0104 average precision, which is inside the fold to fold noise we measured in Part 15. Every operational property favours the tree model, and operational properties are what break at three in the morning. What I would avoid is the unscaled network, obviously, and equally the three layer architecture, which cost 73 percent more compute to score worse than a single wide layer.

Start with boosting, add a network when features stop being tabular

Concretely, for our churn project: keep the tuned HistGradientBoostingClassifier from Part 15 as the production model. Keep the scaled MLP as a benchmark you re run whenever the feature set changes materially, because the day it pulls decisively ahead is a real signal that your features have grown structure a tree cannot see. Do not ship a network that wins by 0.0104 average precision. Do ship one that wins by 0.05.

And remember what the numbers in this part actually taught, which was less about architecture than about measurement. Accuracy separated our four models by one percentage point while average precision separated them by a factor of five. A scaler going missing was invisible on one metric and catastrophic on the other. Choose your metric before you choose your model, a habit that goes all the way back to Part 11.

Part 17 takes the thirty line numpy loop from this part and rewrites it in PyTorch, where autograd computes those three backward lines for you and a readable training loop replaces the fixed step gradient descent we used here. Before you move on, run the loss curve plot on your own model and put the majority class accuracy on the same axis as a dashed line. If your validation curve sits close to that line, your network has learned nothing worth shipping, and you will know it in thirty seconds rather than after a retention campaign.

References

• scikit-learn API reference, MLPClassifier, for hidden_layer_sizes defaulting to (100,), solver defaulting to adam, and the loss_curve_ and validation_scores_ attributes, verified against scikit-learn 1.7.2
• scikit-learn user guide, Neural network models (supervised), which states plainly that multi layer perceptron is sensitive to feature scaling and recommends StandardScaler
• scikit-learn API reference, StandardScaler, for the fit on train and transform on test contract that keeps scaling out of your leakage surface
• NumPy reference, random Generator, for default_rng, the modern replacement for np.random.seed, verified against numpy 2.2.6
• Telco Customer Churn dataset, the public IBM sample file this series works from

About The Author

Dr. Pranay Jha

Dr. Pranay Jha is a Cloud and AI Consultant with 18+ years of experience in hybrid cloud, virtualization, and enterprise infrastructure transformation. He specializes in VMware technologies, multi-cloud strategy, and Generative AI solutions. He holds a PhD in Computer Applications with research focused on Cloud and AI, has published multiple research papers, and has been a VMware vExpert since 2016 and a VMUG Community Leader.

See author's posts