Chapter 43, Projects

Advanced Project: a neural network from scratch

Propagation, backpropagation, stable softmax, momentum optimizer, a gradient check by finite differences, and a spiral that no straight line separates. All just with NumPy, in about 130 lines of code (not counting comments and blank lines). Do this after chapter 37.

The files for this chapter are in projetos/rede_neural/.

The problem

A classifier separates points into classes. When a line is enough, regression solves. The spiral of the project has three intertwined arms, and no line (nor any simple combination of lines) separates them. A neural network solves it because composes linear transformations with a non-linearity (the ReLU), and each layer folds the space a little more. What the project proves is that there is nothing mysterious: it is matrix multiplications and a gradient.

Request	What it does	Chapter
Dense layers	`a @ W + b`, in batch	23
ReLU	`np.maximum(z, 0)`	21
Stable softmax	Probabilities without overflow	32, 36
Cross-entropy	Measures the error, with `log` safe	15, 36
Backpropagation	The chain rule, with transposes	23
Gradient checking	Finite differences against the analytical	36
Initialization	Seeded generator, He scale	16, 35

Mathematics, in four lines

For a batch of n examples, with a₀ = X:

Front: zᵢ = aᵢ₋₁ @ Wᵢ + bᵢ, and aᵢ = ReLU(zᵢ) (the last layer has no ReLU, and it feeds the softmax).
Loss: the average cross-entropy, −log(p do rótulo certo), plus an L2 term that discourages large weights.
The gradient that simplifies everything: at the output, the gradient of the cross-entropy with softmax is (probabilidades − rótulos) / n. A subtraction, without deriving anything.
Back: the weight gradient is aᵢ₋₁ᵀ @ dz, the bias gradient is the sum of dz in the lines, and the gradient that goes back to the previous layer is (dz @ Wᵢᵀ) · (zᵢ₋₁ > 0), where the mask is the derivative of ReLU.

The code

The data: the spiral, the exclusive or (the simplest problem that a line cannot solve) and the split into training and testing. Note that the espiral generates the three classes without a loop, with broadcasting between (1, n) and (classes, 1):

projetos/rede_neural/src/rede_neural/dados.py

import numpy as np

from rede_neural.tipos import Array, Inteiros


def espiral(
    n_por_classe: int, classes: int, ruido: float, rng: np.random.Generator
) -> tuple[Array, Inteiros]:
    """Braços de espiral entrelaçados: um problema que nenhuma reta separa."""
    raio = np.linspace(0.0, 1.0, n_por_classe)
    angulo = np.linspace(0.0, 4.0, n_por_classe)[None, :] + 4.0 * np.arange(classes)[:, None]
    angulo = angulo + rng.normal(0.0, ruido, size=(classes, n_por_classe))
    pontos = np.stack([raio * np.sin(angulo), raio * np.cos(angulo)], axis=-1)
    X: Array = pontos.reshape(-1, 2)
    y: Inteiros = np.repeat(np.arange(classes), n_por_classe)
    return X, y


def ou_exclusivo(n: int, rng: np.random.Generator) -> tuple[Array, Inteiros]:
    """Pontos em [-1, 1]²: a classe é 1 quando os sinais das coordenadas são diferentes."""
    X: Array = rng.uniform(-1.0, 1.0, size=(n, 2))
    y: Inteiros = ((X[:, 0] > 0) != (X[:, 1] > 0)).astype(np.int64)
    return X, y


def dividir(
    X: Array, y: Inteiros, fracao_teste: float, rng: np.random.Generator
) -> tuple[Array, Array, Inteiros, Inteiros]:
    """Separa treino e teste depois de embaralhar (os dados vêm ordenados por classe)."""
    ordem = rng.permutation(len(X))
    corte = int(len(X) * (1 - fracao_teste))
    treino, teste = ordem[:corte], ordem[corte:]
    return X[treino], X[teste], y[treino], y[teste]

The network. The method perda_e_gradientes is the heart: the front holds the inputs and the pre-activations of each layer, and the back reuses them. He initialization scales the weights by the inverse of the square root of the number of inputs, which keeps the variance of the signals stable from layer to layer:

projetos/rede_neural/src/rede_neural/rede.py

import numpy as np

from rede_neural.tipos import Array, Inteiros


def relu(z: Array) -> Array:
    resultado: Array = np.maximum(z, 0.0)
    return resultado


def softmax(z: Array) -> Array:
    """Probabilidades por linha. Subtrair o máximo evita o overflow do exp (capítulo 36)."""
    e = np.exp(z - z.max(axis=1, keepdims=True))
    resultado: Array = e / e.sum(axis=1, keepdims=True)
    return resultado


class RedeDensa:
    """Camadas totalmente conectadas com ReLU e uma saída softmax.

    Convenção de formas: ``X`` é ``(n, entradas)``, o peso da camada i é ``(entradas_i, saidas_i)``
    e a saída é ``(n, classes)``. A forma segue a do scikit-learn: uma linha por exemplo.
    """

    def __init__(self, tamanhos: list[int], rng: np.random.Generator, l2: float = 0.0) -> None:
        if len(tamanhos) < 2:
            raise ValueError("informe ao menos a entrada e a saída")
        self.l2 = l2
        # Inicialização de He: o desvio cresce com o inverso da raiz do número de entradas.
        self.pesos: list[Array] = [
            rng.normal(0.0, np.sqrt(2.0 / entrada), size=(entrada, saida))
            for entrada, saida in zip(tamanhos[:-1], tamanhos[1:], strict=True)
        ]
        self.vieses: list[Array] = [np.zeros(saida) for saida in tamanhos[1:]]

    def parametros(self) -> list[Array]:
        """Os próprios arrays (não cópias): o otimizador e o teste de gradiente os alteram."""
        return [*self.pesos, *self.vieses]

    def _frente(self, X: Array) -> tuple[list[Array], list[Array]]:
        """Devolve a entrada de cada camada e as pré-ativações (antes da ReLU)."""
        entradas = [X]
        pre_ativacoes: list[Array] = []
        a = X
        ultima = len(self.pesos) - 1
        for i, (W, b) in enumerate(zip(self.pesos, self.vieses, strict=True)):
            z = a @ W + b
            pre_ativacoes.append(z)
            a = z if i == ultima else relu(z)
            entradas.append(a)
        return entradas, pre_ativacoes

    def probabilidades(self, X: Array) -> Array:
        _, pre_ativacoes = self._frente(X)
        return softmax(pre_ativacoes[-1])

    def prever(self, X: Array) -> Inteiros:
        previsto: Inteiros = self.probabilidades(X).argmax(axis=1)
        return previsto

    def perda_e_gradientes(self, X: Array, y: Inteiros) -> tuple[float, list[Array]]:
        """Perda de entropia cruzada (mais L2) e o gradiente em relação a cada parâmetro."""
        n = len(X)
        entradas, pre_ativacoes = self._frente(X)
        probs = softmax(pre_ativacoes[-1])
        perda = -np.log(np.maximum(probs[np.arange(n), y], 1e-12)).mean()
        perda += 0.5 * self.l2 * sum(float((W**2).sum()) for W in self.pesos)

        # O gradiente da entropia cruzada com softmax é simplesmente (probabilidade - rótulo) / n.
        dz = probs.copy()
        dz[np.arange(n), y] -= 1.0
        dz /= n

        grad_pesos: list[Array] = []
        grad_vieses: list[Array] = []
        for i in range(len(self.pesos) - 1, -1, -1):
            grad_pesos.insert(0, entradas[i].T @ dz + self.l2 * self.pesos[i])
            grad_vieses.insert(0, dz.sum(axis=0))
            if i > 0:
                dz = (dz @ self.pesos[i].T) * (pre_ativacoes[i - 1] > 0)
        return float(perda), [*grad_pesos, *grad_vieses]

The optimizer, with momentum: the velocity accumulates the previous steps, and this crosses narrow valleys of the loss function faster than pure gradient. It alters the parameters in place, and that is why parametros() returns the arrays themselves, not copies:

projetos/rede_neural/src/rede_neural/otimizador.py

import numpy as np

from rede_neural.tipos import Array


class SGDMomento:
    """Gradiente descendente com momento: a velocidade acumula os passos anteriores."""

    def __init__(self, parametros: list[Array], taxa: float = 0.1, momento: float = 0.9) -> None:
        self.parametros = parametros
        self.taxa = taxa
        self.momento = momento
        self.velocidades = [np.zeros_like(p) for p in parametros]

    def passo(self, gradientes: list[Array]) -> None:
        for p, v, g in zip(self.parametros, self.velocidades, gradientes, strict=True):
            v *= self.momento
            v -= self.taxa * g
            p += v

The training loop:

projetos/rede_neural/src/rede_neural/treino.py

from dataclasses import dataclass, field

import numpy as np

from rede_neural.otimizador import SGDMomento
from rede_neural.rede import RedeDensa
from rede_neural.tipos import Array, Inteiros


@dataclass
class Historico:
    perdas: list[float] = field(default_factory=list)
    acuracias: list[float] = field(default_factory=list)


def acuracia(rede: RedeDensa, X: Array, y: Inteiros) -> float:
    return float((rede.prever(X) == y).mean())


def treinar(
    rede: RedeDensa,
    X: Array,
    y: Inteiros,
    *,
    epocas: int,
    taxa: float = 0.1,
    momento: float = 0.9,
    tamanho_lote: int | None = None,
    rng: np.random.Generator | None = None,
) -> Historico:
    """Treina a rede. Sem ``tamanho_lote``, cada época é um único passo com todos os dados."""
    otimizador = SGDMomento(rede.parametros(), taxa, momento)
    rng = rng or np.random.default_rng(0)
    historico = Historico()
    for _ in range(epocas):
        if tamanho_lote is None:
            _, gradientes = rede.perda_e_gradientes(X, y)
            otimizador.passo(gradientes)
        else:
            ordem = rng.permutation(len(X))
            for inicio in range(0, len(X), tamanho_lote):
                lote = ordem[inicio : inicio + tamanho_lote]
                _, gradientes = rede.perda_e_gradientes(X[lote], y[lote])
                otimizador.passo(gradientes)
        perda, _ = rede.perda_e_gradientes(X, y)
        historico.perdas.append(perda)
        historico.acuracias.append(acuracia(rede, X, y))
    return historico

And the command line:

projetos/rede_neural/src/rede_neural/cli.py

import argparse

import numpy as np

from rede_neural.dados import dividir, espiral, ou_exclusivo
from rede_neural.rede import RedeDensa
from rede_neural.treino import acuracia, treinar


def criar_parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(prog="rede", description="Rede neural só com NumPy")
    parser.add_argument("--semente", type=int, default=0)
    sub = parser.add_subparsers(dest="problema", required=True)
    espiral_p = sub.add_parser("espiral", help="3 braços de espiral entrelaçados")
    espiral_p.add_argument("--epocas", type=int, default=1000)
    espiral_p.add_argument("--oculta", type=int, default=64, help="neurônios por camada oculta")
    espiral_p.add_argument("--taxa", type=float, default=0.1)
    xor_p = sub.add_parser("xor", help="o ou exclusivo, que nenhuma reta resolve")
    xor_p.add_argument("--epocas", type=int, default=500)
    return parser


def main(argv: list[str] | None = None) -> int:
    args = criar_parser().parse_args(argv)
    rng = np.random.default_rng(args.semente)
    if args.problema == "espiral":
        X, y = espiral(100, 3, 0.2, rng)
        tamanhos = [2, args.oculta, args.oculta, 3]
        taxa = args.taxa
    else:
        X, y = ou_exclusivo(200, rng)
        tamanhos = [2, 8, 2]
        taxa = 0.1
    X_treino, X_teste, y_treino, y_teste = dividir(X, y, 0.2, rng)
    rede = RedeDensa(tamanhos, rng, l2=1e-4)

    passo = max(args.epocas // 5, 1)
    historico = treinar(rede, X_treino, y_treino, epocas=args.epocas, taxa=taxa)
    for epoca in range(passo, args.epocas + 1, passo):
        perda, acc = historico.perdas[epoca - 1], historico.acuracias[epoca - 1]
        print(f"época {epoca:>5}  perda {perda:.4f}  acurácia de treino {acc:.3f}")
    print(f"acurácia no treino: {acuracia(rede, X_treino, y_treino):.3f}")
    print(f"acurácia no teste:  {acuracia(rede, X_teste, y_teste):.3f}")
    return 0

The tests

The most important test of this project is the one that checks the gradient. Writing backpropagation by hand is easy to get wrong (one sign, one transpose), and the error does not crash anything: the network just learns poorly. The finite difference check, from chapter 36, perturbs each parameter up and down, measures the change in loss, and compares it with the calculated gradient. If they agree in all the parameters, backpropagation is correct:

projetos/rede_neural/tests/test_rede.py

import numpy as np

from rede_neural.dados import dividir, espiral, ou_exclusivo
from rede_neural.otimizador import SGDMomento
from rede_neural.rede import RedeDensa, softmax
from rede_neural.tipos import Array
from rede_neural.treino import acuracia, treinar


def test_probabilidades_tem_a_forma_certa_e_somam_1() -> None:
    rede = RedeDensa([4, 6, 3], np.random.default_rng(0))
    p = rede.probabilidades(np.random.default_rng(1).normal(size=(5, 4)))
    assert p.shape == (5, 3)
    assert np.allclose(p.sum(axis=1), 1.0)


def test_softmax_e_estavel_com_valores_enormes() -> None:
    p = softmax(np.array([[1000.0, 1001.0, 1002.0]]))
    assert np.isfinite(p).all() and np.isclose(p.sum(), 1.0)


def test_gradiente_confere_com_diferencas_finitas() -> None:
    """O teste mais importante: o gradiente escrito à mão bate com a derivada numérica."""
    rng = np.random.default_rng(3)
    rede = RedeDensa([3, 5, 4, 2], rng, l2=0.01)
    X = rng.normal(size=(6, 3))
    y = rng.integers(0, 2, size=6)
    _, analiticos = rede.perda_e_gradientes(X, y)

    h = 1e-5
    for parametro, analitico in zip(rede.parametros(), analiticos, strict=True):
        numerico = np.zeros_like(parametro)
        for indice in np.ndindex(parametro.shape):
            original = parametro[indice]
            parametro[indice] = original + h
            mais, _ = rede.perda_e_gradientes(X, y)
            parametro[indice] = original - h
            menos, _ = rede.perda_e_gradientes(X, y)
            parametro[indice] = original
            numerico[indice] = (mais - menos) / (2 * h)
        assert np.allclose(analitico, numerico, rtol=1e-4, atol=1e-7)


def test_otimizador_passo_simples_e_momento() -> None:
    p: Array = np.array([1.0])
    sgd = SGDMomento([p], taxa=0.1, momento=0.0)
    sgd.passo([np.array([2.0])])
    assert np.allclose(p, [0.8])

    q: Array = np.array([1.0])
    com_momento = SGDMomento([q], taxa=0.1, momento=0.5)
    com_momento.passo([np.array([2.0])])
    com_momento.passo([np.array([2.0])])
    assert np.allclose(q, [1.0 - 0.2 - 0.3])


def test_a_perda_cai_durante_o_treino() -> None:
    rng = np.random.default_rng(0)
    X, y = espiral(60, 3, 0.2, rng)
    rede = RedeDensa([2, 32, 3], rng)
    historico = treinar(rede, X, y, epocas=200)
    assert historico.perdas[-1] < 0.5 * historico.perdas[0]


def test_aprende_o_ou_exclusivo() -> None:
    rng = np.random.default_rng(0)
    X, y = ou_exclusivo(300, rng)
    Xt, Xv, yt, yv = dividir(X, y, 0.2, rng)
    rede = RedeDensa([2, 8, 2], rng, l2=1e-4)
    treinar(rede, Xt, yt, epocas=500)
    assert acuracia(rede, Xv, yv) >= 0.95


def test_resolve_a_espiral_em_dados_novos() -> None:
    rng = np.random.default_rng(0)
    X, y = espiral(100, 3, 0.2, rng)
    Xt, Xv, yt, yv = dividir(X, y, 0.2, rng)
    rede = RedeDensa([2, 64, 64, 3], rng, l2=1e-4)
    treinar(rede, Xt, yt, epocas=1000)
    assert acuracia(rede, Xt, yt) > 0.97
    assert acuracia(rede, Xv, yv) > 0.95


def test_mesma_semente_mesmo_resultado() -> None:
    def rodar() -> float:
        rng = np.random.default_rng(5)
        X, y = espiral(40, 3, 0.2, rng)
        rede = RedeDensa([2, 16, 3], rng)
        return treinar(rede, X, y, epocas=50).perdas[-1]

    assert rodar() == rodar()

projetos/rede_neural/tests/test_dados.py

import numpy as np

from rede_neural.dados import dividir, espiral, ou_exclusivo


def test_espiral_formas_e_classes() -> None:
    X, y = espiral(50, 3, 0.2, np.random.default_rng(0))
    assert X.shape == (150, 2) and y.shape == (150,)
    assert np.bincount(y).tolist() == [50, 50, 50]


def test_ou_exclusivo_rotulos() -> None:
    X, y = ou_exclusivo(300, np.random.default_rng(1))
    esperado = ((X[:, 0] > 0) != (X[:, 1] > 0)).astype(int)
    assert np.array_equal(y, esperado)
    assert set(np.unique(y)) == {0, 1}


def test_dividir_e_disjunto_e_embaralha() -> None:
    X, y = espiral(50, 3, 0.2, np.random.default_rng(0))
    Xt, Xv, yt, yv = dividir(X, y, 0.2, np.random.default_rng(2))
    assert len(Xt) == 120 and len(Xv) == 30
    assert len(np.unique(yv)) > 1  # o teste não é só de uma classe, porque embaralhou antes

projetos/rede_neural/tests/test_cli.py

import pytest

from rede_neural.cli import main


def test_xor_pela_linha_de_comando(capsys: pytest.CaptureFixture[str]) -> None:
    assert main(["xor", "--epocas", "300"]) == 0
    saida = capsys.readouterr().out
    assert "acurácia no teste" in saida
    assert "época" in saida

Run

Terminal

cd projetos/rede_neural
uv sync
uv run pytest
uv run rede espiral

Output

............                                                             [100%]
12 passed

Treino na espiral (execução real)

época   200  perda 0.0558  acurácia de treino 0.996
época   400  perda 0.0445  acurácia de treino 0.996
época   600  perda 0.0407  acurácia de treino 0.996
época   800  perda 0.0387  acurácia de treino 0.996
época  1000  perda 0.0374  acurácia de treino 0.996
acurácia no treino: 0.996
acurácia no teste:  0.983

And in the exclusive or, which a line can never solve:

Terminal

uv run rede xor

Ou exclusivo (execução real)

época   100  perda 0.0994  acurácia de treino 0.981
época   200  perda 0.0707  acurácia de treino 0.981
época   300  perda 0.0603  acurácia de treino 0.988
época   400  perda 0.0547  acurácia de treino 0.988
época   500  perda 0.0511  acurácia de treino 0.994
acurácia no treino: 0.994
acurácia no teste:  1.000

The accuracy in testing (98.3%) is what matters, because these are points that the network has never seen. Close training and testing indicate that it generalized, instead of memorizing. I did not choose seed 0 by chance: with these parameters (two hidden layers of 64 neurons, rate 0.1, momentum 0.9, 1000 epochs), I tested six different seeds and all passed 98% in testing.

See what the network learned

The fronteira.py classifies all points of a grid of the plane at once (the meshgrid from chapter 27, plus a single vectorized call to prever) and paints each region with the predicted class. It requires matplotlib: uv run --with matplotlib python fronteira.py:

projetos/rede_neural/fronteira.py

"""Desenha a fronteira de decisão que a rede aprendeu na espiral.

Precisa do matplotlib, que não é dependência do projeto:
    uv run --with matplotlib python fronteira.py
"""

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from rede_neural.dados import dividir, espiral
from rede_neural.rede import RedeDensa
from rede_neural.treino import treinar

rng = np.random.default_rng(0)
X, y = espiral(100, 3, 0.2, rng)
X_treino, X_teste, y_treino, y_teste = dividir(X, y, 0.2, rng)
rede = RedeDensa([2, 64, 64, 3], rng, l2=1e-4)
treinar(rede, X_treino, y_treino, epocas=1000)

# Uma malha de pontos cobrindo o plano: a rede classifica todos de uma vez (vetorizado).
xs = np.linspace(-1.2, 1.2, 300)
malha_x, malha_y = np.meshgrid(xs, xs)
pontos = np.column_stack([malha_x.ravel(), malha_y.ravel()])
classes = rede.prever(pontos).reshape(malha_x.shape)

fig, ax = plt.subplots(figsize=(5, 5), dpi=100)
ax.contourf(malha_x, malha_y, classes, alpha=0.35, levels=[-0.5, 0.5, 1.5, 2.5])
ax.scatter(X_treino[:, 0], X_treino[:, 1], c=y_treino, s=10, edgecolors="none")
ax.scatter(X_teste[:, 0], X_teste[:, 1], c=y_teste, s=28, edgecolors="black", linewidths=0.6)
ax.set_title("Fronteira aprendida (pontos com borda: teste)")
ax.set_aspect("equal")
fig.savefig("fronteira.png")

The regions that the network learned for the three classes. The points with a black border are from the test set, which the network has never seen during training.

The boundary follows the wrapping of the arms, and the test points fall, in almost all cases, in the region of their own color.

What this project is not

It is an exercise to understand, and not a substitute for a framework. It only has dense layers, and the backpropagation was written by hand for this architecture. For any other, PyTorch and JAX derive the gradient automatically (and run on GPU), and that is why no one writes this by hand in production. The advantage of having written it once is that when a training does not converge, you know where to look: the gradient, the scale of the initialization, the learning rate.

Challenges

Adam. Replace the optimizer with Adam (moving averages of the gradient and its square) and compare the loss per epoch with momentum.
Mini-batches. The treinar already accepts tamanho_lote. Compare the convergence with batches of 32 and with the full batch, in time and in loss.
Dropout. Add dropout regularization (randomly zeroing a fraction of the neurons during training) and check the gradient again: does the test still pass?
Early stopping. Stop training when the loss on the validation set stops decreasing.
Compare with scikit-learn. Train a MLPClassifier (chapter 40) on the same data and compare the accuracies. Does your network come close?

What I learned from this type of project
Writing backpropagation once demystifies deep learning, and the gradient checking is the habit that most prevents hours of debugging. In any hand-written numerical code, I compare the result with a slow and obviously correct version before trusting it. It is the same principle as the tests from the other two projects.