# Foundations of Artificial Neural Networks

## 1. From the Biological Neuron to the Artificial Neuron

The artificial neuron is a radical abstraction of its biological counterpart. A biological neuron receives electrochemical signals through its **dendrites**, integrates them in the **soma** (cell body), and — if the aggregated membrane depolarization crosses a threshold — emits an **action potential** that propagates along the **axon** to **synapses** connecting to downstream neurons. Two properties of this system are retained in the mathematical abstraction: (i) synaptic transmission is *weighted* (a synapse may be excitatory or inhibitory, strong or weak), and (ii) firing is *nonlinear and threshold-based* (all-or-none).

Everything else is discarded: spike timing, refractory periods, dendritic compartmentalization, neuromodulation, and the fact that biological learning is local and largely unsupervised. The artificial neuron is therefore best understood as a *computational primitive inspired by* — not a *model of* — neurophysiology.

### The McCulloch–Pitts neuron (1943)

The foundational text is Warren S. McCulloch and Walter Pitts, *"A Logical Calculus of the Ideas Immanent in Nervous Activity"*, **Bulletin of Mathematical Biophysics**, 5(4):115–133, 1943. It is widely credited as a seminal contribution to neural network theory, automata theory, the theory of computation, and cybernetics.

The MP neuron takes Boolean inputs $x_i \in \{0,1\}$ and produces a Boolean output via a threshold (Heaviside) function:

$$y = \Theta\!\left(\sum_{i=1}^{n} w_i x_i - \theta\right), \qquad \Theta(u) = \begin{cases} 1 & \text{if } u \geq 0 \\ 0 & \text{if } u < 0 \end{cases}$$

In the original formulation, weights are fixed at $w_i = +1$ for excitatory inputs, and inhibitory inputs are *absolute*: a single active inhibitory input vetoes firing regardless of the excitatory sum. Writing $x_1,\dots,x_n$ for excitatory and $z_1,\dots,z_m$ for inhibitory inputs:

$$y = \Theta\!\left(\sum_{i=1}^{n} x_i - \theta\right)\prod_{j=1}^{m}(1 - z_j)$$

McCulloch and Pitts showed that networks of such units can implement any Boolean function — AND ($\theta = n$), OR ($\theta = 1$), NOT (via inhibition) — and with cycles, any finite-state automaton. The decisive limitation is that **the MP neuron does not learn**: $w_i$ and $\theta$ are set by the designer.

---

## 2. The Perceptron (Rosenblatt, 1958)

Frank Rosenblatt introduced the perceptron in *"The Perceptron: A Probabilistic Model for Information Storage and Organization in the Brain"*, **Psychological Review**, 65(6):386–408, 1958, developed at length in *Principles of Neurodynamics* (1962). The critical advance over McCulloch–Pitts is a **learning procedure**: weights are adjusted from examples rather than hand-designed.

### Output equation

With real-valued inputs $\mathbf{x} \in \mathbb{R}^n$, weights $\mathbf{w} \in \mathbb{R}^n$, and bias $b$ (equivalently, a negative threshold):

$$z = \mathbf{w}^\top \mathbf{x} + b = \sum_{i=1}^{n} w_i x_i + b, \qquad \hat{y} = \operatorname{sign}(z) = \begin{cases} +1 & z \geq 0 \\ -1 & z < 0\end{cases}$$

The bias is conventionally absorbed by augmenting $\mathbf{x} \leftarrow (\mathbf{x}, 1)$ and $\mathbf{w} \leftarrow (\mathbf{w}, b)$, giving $\hat{y} = \operatorname{sign}(\mathbf{w}^\top \mathbf{x})$. Geometrically, the perceptron defines a **hyperplane** $\mathbf{w}^\top \mathbf{x} + b = 0$ splitting input space into two half-spaces.

### The perceptron learning rule

For each misclassified example $(\mathbf{x}^{(k)}, y^{(k)})$, with learning rate $\eta > 0$:

$$\mathbf{w} \leftarrow \mathbf{w} + \eta\left(y^{(k)} - \hat{y}^{(k)}\right)\mathbf{x}^{(k)}, \qquad b \leftarrow b + \eta\left(y^{(k)} - \hat{y}^{(k)}\right)$$

In the $\pm1$ convention this simplifies: correctly classified points produce no update, and a misclassified point gives $\mathbf{w} \leftarrow \mathbf{w} + \eta\, y^{(k)} \mathbf{x}^{(k)}$. This is error-driven, online, and requires no differentiability — the step function's derivative is zero almost everywhere, so this is *not* gradient descent on the 0-1 loss. It is, however, equivalent to stochastic subgradient descent on the **perceptron criterion** $L = \max(0, -y\,\mathbf{w}^\top\mathbf{x})$.

### The perceptron convergence theorem

Formalized by Novikoff (*"On Convergence Proofs for Perceptrons"*, 1962): suppose a unit vector $\mathbf{w}^\star$ with $\|\mathbf{w}^\star\| = 1$ separates the data with margin $\gamma > 0$, i.e. $y_i(\mathbf{w}^{\star\top}\mathbf{x}_i) \geq \gamma$ for all $i$, and let $R = \max_i \|\mathbf{x}_i\|$. Then the perceptron algorithm makes at most

$$T \leq \left(\frac{R}{\gamma}\right)^2$$

updates before converging. Note the bound is independent of dimension and of the number of samples — it depends only on the normalized margin $\gamma/R$.

### Limits: Minsky & Papert and the XOR problem

Marvin Minsky and Seymour Papert, *Perceptrons: An Introduction to Computational Geometry* (MIT Press, 1969), gave a rigorous analysis of what single-layer perceptrons cannot represent. The canonical counterexample is **XOR**:

| $x_1$ | $x_2$ | XOR |
|---|---|---|
| 0 | 0 | 0 |
| 0 | 1 | 1 |
| 1 | 0 | 1 |
| 1 | 1 | 0 |

Suppose a separating line existed. Then $b < 0$ (from $(0,0)\mapsto 0$), $w_1 + b \geq 0$ and $w_2 + b \geq 0$ (from the positive cases), hence $w_1 + w_2 + 2b \geq 0$, so $w_1 + w_2 + b \geq -b > 0$ — which forces $(1,1)\mapsto 1$, a contradiction. The positive vertices lie on one diagonal of the unit square, the negative on the other; no line separates two diagonals.

Minsky and Papert also proved deeper results — e.g. that the **parity** and **connectedness** predicates require perceptrons whose order (number of inputs per feature detector) grows with the input size. Their pessimism about multilayer extensions, combined with a lack of a training algorithm for hidden layers, is commonly cited as a trigger for the first "AI winter" in connectionism, lasting roughly until 1986.

Crucially, XOR *is* solvable by two layers: $\text{XOR}(x_1,x_2) = \text{OR}(x_1,x_2) \wedge \neg\text{AND}(x_1,x_2)$, i.e. one hidden layer of two threshold units plus an output unit.

---

## 3. ADALINE and MADALINE (Widrow & Hoff, 1960)

Bernard Widrow and Marcian E. Hoff, *"Adaptive Switching Circuits"*, **IRE WESCON Convention Record**, Part 4, pp. 96–104, 1960, introduced **ADALINE** (ADAptive LINear Element / Adaptive Linear Neuron).

The key difference from the perceptron: the **error is computed on the linear pre-activation**, not on the thresholded output. Let $z = \mathbf{w}^\top\mathbf{x} + b$ and target $d$. The objective is squared error:

$$L(\mathbf{w}) = \tfrac{1}{2}(d - z)^2$$

Differentiating, $\partial L/\partial w_i = -(d-z)x_i$, giving the **delta rule** (also known as the Widrow–Hoff rule or the **LMS**, least-mean-squares, algorithm):

$$\Delta w_i = \eta\,(d - z)\,x_i, \qquad \mathbf{w} \leftarrow \mathbf{w} + \eta\,(d - z)\,\mathbf{x}$$

Because $z$ is differentiable in $\mathbf{w}$, this *is* genuine stochastic gradient descent, and it converges (in the mean) for $0 < \eta < 2/\lambda_{\max}$ where $\lambda_{\max}$ is the largest eigenvalue of the input autocorrelation matrix $R = \mathbb{E}[\mathbf{x}\mathbf{x}^\top]$. Unlike the perceptron rule, LMS keeps improving even when the data are already separated, driving the solution toward the minimum-MSE hyperplane — and it does not diverge on non-separable data.

**MADALINE** (Many ADALINEs) stacks ADALINE units into a layered network with a fixed voting/majority output unit. MADALINE Rule I (MRI, 1962) adapts the hidden ADALINE whose pre-activation is closest to zero — the "minimum disturbance" principle — flipping it if this reduces error. MADALINE Rule II (MRII, Winter & Widrow, 1988) generalizes the trial-and-adopt search to multiple layers. MADALINE III later approximated backpropagation with sigmoidal units. LMS/ADALINE became one of the first industrially deployed neural methods, notably in adaptive echo cancellation and channel equalization in telephony, and remains the workhorse of adaptive signal processing.

---

## 4. The Multilayer Perceptron (MLP)

An MLP is a **feedforward** network of $L$ layers where each layer applies an affine map followed by a pointwise nonlinearity.

### Forward propagation, layer by layer

Let $\mathbf{a}^{[0]} = \mathbf{x} \in \mathbb{R}^{n_0}$. For $\ell = 1, \dots, L$ with $\mathbf{W}^{[\ell]} \in \mathbb{R}^{n_\ell \times n_{\ell-1}}$ and $\mathbf{b}^{[\ell]} \in \mathbb{R}^{n_\ell}$:

$$\mathbf{z}^{[\ell]} = \mathbf{W}^{[\ell]}\mathbf{a}^{[\ell-1]} + \mathbf{b}^{[\ell]}$$
$$\mathbf{a}^{[\ell]} = \sigma^{[\ell]}\!\left(\mathbf{z}^{[\ell]}\right)$$

with the network output $\hat{\mathbf{y}} = \mathbf{a}^{[L]}$. Componentwise, $z_j^{[\ell]} = \sum_{i} W_{ji}^{[\ell]} a_i^{[\ell-1]} + b_j^{[\ell]}$.

In **batched** form with a design matrix $\mathbf{X} \in \mathbb{R}^{m \times n_0}$ ($m$ examples as rows), $\mathbf{Z}^{[\ell]} = \mathbf{A}^{[\ell-1]}\mathbf{W}^{[\ell]\top} + \mathbf{1}_m \mathbf{b}^{[\ell]\top}$, which maps directly onto GEMM primitives — the reason MLPs are hardware-efficient.

The nonlinearity is essential: composing affine maps yields an affine map, so a network of any depth with $\sigma = \text{id}$ collapses to a single linear layer.

### The Universal Approximation Theorem

**Cybenko (1989)** — George Cybenko, *"Approximation by Superpositions of a Sigmoidal Function"*, **Mathematics of Control, Signals and Systems**, 2(4):303–314, 1989 — proved that finite sums of the form

$$G(\mathbf{x}) = \sum_{j=1}^{N} \alpha_j\, \sigma\!\left(\mathbf{w}_j^\top\mathbf{x} + \theta_j\right)$$

are **dense** in $C(I_n)$, the continuous functions on the unit hypercube $I_n = [0,1]^n$, under the uniform norm, whenever $\sigma$ is any continuous **sigmoidal** function ($\sigma(t)\to 1$ as $t\to+\infty$, $\sigma(t)\to 0$ as $t\to-\infty$). Formally: for any $f \in C(I_n)$ and $\varepsilon > 0$ there exists such a $G$ with $|G(\mathbf{x}) - f(\mathbf{x})| < \varepsilon$ for all $\mathbf{x} \in I_n$. Cybenko's proof is non-constructive, relying on the Hahn–Banach theorem and the Riesz representation theorem to show that the closure of the span of these functions cannot be a proper subspace.

**Hornik, Stinchcombe & White (1989)** — *"Multilayer Feedforward Networks are Universal Approximators"*, **Neural Networks**, 2(5):359–366 — obtained the result independently and more generally, showing that single-hidden-layer networks with any *squashing* activation are universal approximators for Borel measurable functions, in $L^p(\mu)$ for arbitrary finite measures $\mu$, and with derivatives (Hornik, 1991, *"Approximation Capabilities of Multilayer Feedforward Networks"*, **Neural Networks** 4(2):251–257).

**Leshno, Lin, Pinkus & Schocken (1993)** — *"Multilayer Feedforward Networks with a Nonpolynomial Activation Function Can Approximate Any Function"*, **Neural Networks**, 6(6):861–867 — gave the sharpest classical statement: a network with a locally bounded, piecewise-continuous activation is a universal approximator **if and only if the activation is not a polynomial**. This is why ReLU, despite not being sigmoidal, is universal.

Two caveats matter in practice. First, these are **existence** results: they say nothing about how many hidden units are needed (the width $N$ may be exponential in $n$), nor whether gradient descent will *find* the approximating weights. Second, they concern shallow networks; **depth-separation** results (e.g. Telgarsky 2016, Eldan & Shamir 2016) show functions representable by a deep network with polynomially many units that require exponentially many units at shallower depth — the modern justification for depth.

---

## 5. Backpropagation (Rumelhart, Hinton & Williams, 1986)

David E. Rumelhart, Geoffrey E. Hinton and Ronald J. Williams, *"Learning Representations by Back-Propagating Errors"*, **Nature**, 323:533–536, 1986 (DOI: 10.1038/323533a0), popularized the algorithm that made hidden-layer training practical. The paper's stated contribution is a procedure that "repeatedly adjusts the weights of the connections in the network so as to minimize a measure of the difference between the actual output vector of the net and the desired output vector," with the consequence that "hidden units come to represent important features of the task domain" — the ability to *create useful new features* is exactly what distinguishes it from the perceptron convergence procedure.

Historically the method is older: reverse-mode automatic differentiation was described by Seppo Linnainmaa (1970), applied to networks by Paul Werbos in his 1974 Harvard PhD thesis (*Beyond Regression*), and independently derived by Parker (1985) and LeCun (1985).

### Full derivation

Let $L$ be the loss on a single example, with layers indexed $\ell = 1,\dots,L$. Define the **error term** (or "delta") of layer $\ell$ as the gradient of the loss with respect to the pre-activation:

$$\boldsymbol{\delta}^{[\ell]} \;\equiv\; \frac{\partial L}{\partial \mathbf{z}^{[\ell]}} \in \mathbb{R}^{n_\ell}$$

**Output layer.** By the chain rule through $\mathbf{a}^{[L]} = \sigma^{[L]}(\mathbf{z}^{[L]})$:

$$\boldsymbol{\delta}^{[L]} = \nabla_{\mathbf{a}^{[L]}} L \;\odot\; \sigma^{[L]\prime}\!\left(\mathbf{z}^{[L]}\right)$$

where $\odot$ is the Hadamard (elementwise) product.

**Recursive backward pass.** Since $\mathbf{z}^{[\ell+1]} = \mathbf{W}^{[\ell+1]}\sigma^{[\ell]}(\mathbf{z}^{[\ell]}) + \mathbf{b}^{[\ell+1]}$, each $z_k^{[\ell+1]}$ depends on $z_j^{[\ell]}$ through $W^{[\ell+1]}_{kj}\sigma^{[\ell]\prime}(z_j^{[\ell]})$. Summing over all downstream paths:

$$\delta_j^{[\ell]} = \sum_{k} \frac{\partial L}{\partial z_k^{[\ell+1]}}\frac{\partial z_k^{[\ell+1]}}{\partial z_j^{[\ell]}} = \left(\sum_k \delta_k^{[\ell+1]} W_{kj}^{[\ell+1]}\right)\sigma^{[\ell]\prime}\!\left(z_j^{[\ell]}\right)$$

In matrix form:

$$\boxed{\;\boldsymbol{\delta}^{[\ell]} = \left(\mathbf{W}^{[\ell+1]\top}\boldsymbol{\delta}^{[\ell+1]}\right)\odot\sigma^{[\ell]\prime}\!\left(\mathbf{z}^{[\ell]}\right)\;}$$

This is where the name comes from: the error is *propagated backwards* through the transpose of the forward weight matrices.

**Parameter gradients.** Since $z_j^{[\ell]} = \sum_i W_{ji}^{[\ell]}a_i^{[\ell-1]} + b_j^{[\ell]}$, we have $\partial z_j^{[\ell]}/\partial W_{ji}^{[\ell]} = a_i^{[\ell-1]}$ and $\partial z_j^{[\ell]}/\partial b_j^{[\ell]} = 1$, hence

$$\frac{\partial L}{\partial W_{ji}^{[\ell]}} = \delta_j^{[\ell]}\,a_i^{[\ell-1]} \quad\Longleftrightarrow\quad \frac{\partial L}{\partial \mathbf{W}^{[\ell]}} = \boldsymbol{\delta}^{[\ell]}\,\mathbf{a}^{[\ell-1]\top}$$

$$\frac{\partial L}{\partial \mathbf{b}^{[\ell]}} = \boldsymbol{\delta}^{[\ell]}$$

For a mini-batch of size $m$, gradients are averaged: $\partial L/\partial\mathbf{W}^{[\ell]} = \frac{1}{m}\boldsymbol{\Delta}^{[\ell]\top}\mathbf{A}^{[\ell-1]}$ and $\partial L/\partial\mathbf{b}^{[\ell]} = \frac{1}{m}\sum_{k=1}^m \boldsymbol{\delta}^{[\ell](k)}$.

**Complexity.** One backward pass costs the same order as one forward pass, $O(\sum_\ell n_\ell n_{\ell-1})$ — the fundamental efficiency result of reverse-mode automatic differentiation: the full gradient of a scalar with respect to $P$ parameters costs $O(1)$ forward passes, not $O(P)$.

**A useful special case.** With softmax output and categorical cross-entropy, the two Jacobians cancel and the output delta collapses to $\boldsymbol{\delta}^{[L]} = \hat{\mathbf{y}} - \mathbf{y}$. The same holds for sigmoid + binary cross-entropy and for linear output + MSE. This is not a coincidence: it holds for any matched pair of a canonical link function and its exponential-family negative log-likelihood.

**Vanishing/exploding gradients.** The recursion multiplies $\sigma'$ terms at every layer. With sigmoid, $\sigma'(z) \leq 1/4$, so gradients shrink by at least $4^{-L}$ across $L$ layers — the vanishing gradient problem identified by Hochreiter (1991) and Bengio, Simard & Frasconi (1994). This motivates ReLU-family activations, careful initialization, normalization layers, and residual connections.

---

## 6. Activation Functions

| Function | Definition | Derivative | Range |
|---|---|---|---|
| Sigmoid | $\sigma(x) = \dfrac{1}{1+e^{-x}}$ | $\sigma(x)\left(1-\sigma(x)\right)$ | $(0,1)$ |
| Tanh | $\tanh(x) = \dfrac{e^{x}-e^{-x}}{e^{x}+e^{-x}}$ | $1-\tanh^2(x)$ | $(-1,1)$ |
| ReLU | $\max(0,x)$ | $\mathbb{1}[x>0]$ | $[0,\infty)$ |
| Leaky ReLU | $\max(\alpha x, x),\ \alpha{=}0.01$ | $\alpha$ if $x<0$, else $1$ | $(-\infty,\infty)$ |
| PReLU | same, $\alpha$ learned | idem, plus $\partial f/\partial\alpha = \min(0,x)$ | $(-\infty,\infty)$ |
| ELU | $x$ if $x>0$; $\alpha(e^x-1)$ else | $1$ if $x>0$; $\alpha e^x$ else | $(-\alpha,\infty)$ |
| SELU | $\lambda\cdot\text{ELU}_\alpha(x)$ | $\lambda$ if $x>0$; $\lambda\alpha e^x$ else | scaled |
| Softplus | $\ln(1+e^x)$ | $\sigma(x)$ | $(0,\infty)$ |

**Sigmoid.** Historically dominant, now largely confined to gates and binary outputs. Two defects: outputs are not zero-centred (inducing correlated gradient signs across a layer, causing zig-zag descent), and it **saturates** — $\sigma'(x)\to 0$ for $|x|\gtrsim 5$, killing gradient flow. Note $\tanh(x) = 2\sigma(2x)-1$.

**Tanh.** Zero-centred, with max derivative $1$ at the origin. Still saturating, but empirically better-conditioned than sigmoid for hidden layers; still standard in LSTM/GRU cell candidates.

**ReLU.** Used by Fukushima (1969) for visual feature extraction; repopularized by **Nair & Hinton (2010)**, *"Rectified Linear Units Improve Restricted Boltzmann Machines"* (ICML), and **Glorot, Bordes & Bengio (2011)**, *"Deep Sparse Rectifier Neural Networks"* (AISTATS); cemented by AlexNet (Krizhevsky et al., 2012). Advantages: no saturation for $x>0$, gradient exactly $1$ there, trivial to compute, induces sparse activations. Drawback: the **dying ReLU** problem — a unit pushed into the negative regime for all inputs receives zero gradient forever. Non-differentiable at $0$; frameworks conventionally set $f'(0)=0$.

**Leaky ReLU** (Maas, Hannun & Ng, 2013, *"Rectifier Nonlinearities Improve Neural Network Acoustic Models"*, ICML workshop) fixes dying units with a small negative slope. **PReLU** (He, Zhang, Ren & Sun, 2015, *"Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification"*, ICCV) makes $\alpha$ a learned per-channel parameter at negligible cost; the same paper introduced **He/Kaiming initialization**, $\mathrm{Var}(W) = 2/n_{\text{in}}$, which corrects Xavier initialization for the fact that ReLU zeroes half the activations.

**ELU** (Clevert, Unterthiner & Hochreiter, 2015, *"Fast and Accurate Deep Network Learning by Exponential Linear Units"*, ICLR 2016) saturates smoothly to $-\alpha$ for large negative inputs, pushing mean activations toward zero (a batch-norm-like effect) while remaining robust to noise.

**SELU** (Klambauer, Unterthiner, Mayr & Hochreiter, 2017, *"Self-Normalizing Neural Networks"*, NIPS) fixes

$$\lambda \approx 1.0507009873554804934193349852946, \qquad \alpha \approx 1.6732632423543772848170429916717$$

$$\text{SELU}(x) = \lambda\begin{cases} x & x > 0\\ \alpha(e^x - 1) & x \leq 0\end{cases}$$

These constants are derived so that, under LeCun-normal initialization ($\mathrm{Var}(W)=1/n_{\text{in}}$) and with the "alpha-dropout" variant, activation mean and variance have an **attracting fixed point at $(0,1)$** — the network self-normalizes without batch norm. The guarantee is fragile: it requires the specific initialization, fully-connected architecture, and alpha-dropout.

**GELU** (Hendrycks & Gimpel, 2016, *"Gaussian Error Linear Units (GELUs)"*, arXiv:1606.08415) weights the input by the probability that a standard Gaussian falls below it:

$$\text{GELU}(x) = x\,\Phi(x) = \frac{x}{2}\left[1 + \operatorname{erf}\!\left(\frac{x}{\sqrt{2}}\right)\right], \qquad \frac{d}{dx}\text{GELU}(x) = \Phi(x) + x\,\phi(x)$$

where $\phi$ is the standard normal density. The widely-used tanh approximation (used in BERT and GPT-2) is

$$\text{GELU}(x) \approx 0.5\,x\left(1 + \tanh\!\left[\sqrt{\tfrac{2}{\pi}}\left(x + 0.044715\,x^3\right)\right]\right)$$

GELU can be read as a *deterministic* version of stochastic-regularizer gating: instead of dropping a unit with probability $1-\Phi(x)$, it scales by the expected mask. It is the default in most Transformers.

**Swish / SiLU** (Ramachandran, Zoph & Le, 2017, *"Searching for Activation Functions"*, discovered by automated search; the $\beta=1$ case, SiLU, appears earlier in Hendrycks & Gimpel 2016 and Elfwing et al. 2017):

$$\text{Swish}_\beta(x) = x\,\sigma(\beta x), \qquad \text{SiLU}(x) = \frac{x}{1+e^{-x}}$$
$$\text{SiLU}'(x) = \sigma(x)\left(1 + x\left(1 - \sigma(x)\right)\right) = \frac{1 + e^{-x} + x e^{-x}}{(1+e^{-x})^2}$$

Smooth, non-monotonic (a small negative dip near $x\approx-1.28$), unbounded above, bounded below; $\beta$ may be learned, with $\beta\to 0$ recovering a linear unit and $\beta\to\infty$ recovering ReLU. **SwiGLU** — a gated variant $\text{Swish}(xW)\odot(xV)$ (Shazeer, 2020) — is now standard in LLM feed-forward blocks.

**Mish** (Misra, 2019, *"Mish: A Self Regularized Non-Monotonic Activation Function"*, BMVC 2020):

$$\text{Mish}(x) = x\tanh\left(\text{softplus}(x)\right) = x\tanh\left(\ln(1+e^x)\right)$$

Similar in shape to Swish, with a smoother profile; adopted in several YOLO variants.

**Softmax** (Bridle, 1990) converts a logit vector to a probability simplex:

$$\text{softmax}(\mathbf{z})_i = \frac{e^{z_i}}{\sum_{j=1}^{K} e^{z_j}}, \qquad \frac{\partial\,\text{softmax}_i}{\partial z_j} = \text{softmax}_i\left(\delta_{ij} - \text{softmax}_j\right)$$

with $\delta_{ij}$ the Kronecker delta. It is shift-invariant, so implementations subtract $\max_j z_j$ before exponentiating for numerical stability. With a temperature $T$, $\text{softmax}(\mathbf{z}/T)$ interpolates between argmax ($T\to0$) and uniform ($T\to\infty$).

**Softplus** $= \ln(1+e^x)$ is the smooth ReLU; its derivative is exactly the sigmoid. Implemented stably as $\max(x,0) + \ln(1+e^{-|x|})$.

---

## 7. Loss Functions

**Mean Squared Error (L2).** For regression, with $n$ examples:

$$L_{\text{MSE}} = \frac{1}{n}\sum_{i=1}^{n}\left(y_i - \hat{y}_i\right)^2, \qquad \frac{\partial L}{\partial \hat{y}_i} = -\frac{2}{n}(y_i - \hat{y}_i)$$

Corresponds to Gaussian-noise maximum likelihood; strongly penalizes outliers (quadratic growth).

**Mean Absolute Error (L1).**

$$L_{\text{MAE}} = \frac{1}{n}\sum_{i=1}^{n}\left|y_i - \hat{y}_i\right|, \qquad \frac{\partial L}{\partial\hat{y}_i} = -\frac{1}{n}\operatorname{sign}(y_i - \hat{y}_i)$$

Robust to outliers (Laplace-noise MLE; the minimizer is the conditional median rather than the mean), but non-differentiable at zero and with constant gradient magnitude, which impedes fine convergence.

**Huber loss** (Peter J. Huber, *"Robust Estimation of a Location Parameter"*, **Annals of Mathematical Statistics**, 35(1):73–101, 1964), the standard compromise. With $r = y - \hat{y}$ and threshold $\delta$:

$$L_\delta(r) = \begin{cases} \frac{1}{2}r^2 & |r| \leq \delta\\[4pt] \delta\left(|r| - \frac{1}{2}\delta\right) & |r| > \delta\end{cases}, \qquad \frac{dL_\delta}{dr} = \begin{cases} r & |r|\leq\delta\\ \delta\operatorname{sign}(r) & |r| > \delta\end{cases}$$

Quadratic near zero (fast convergence on small residuals), linear in the tails (bounded gradient, outlier-resistant), and $C^1$ everywhere. The *smooth L1* loss used in object detection is $L_\delta/\delta$ with $\delta=1$.

**Binary cross-entropy (log loss).** For $y \in \{0,1\}$ and $\hat{y} = \sigma(z) \in (0,1)$:

$$L_{\text{BCE}} = -\frac{1}{n}\sum_{i=1}^{n}\left[y_i \ln \hat{y}_i + (1-y_i)\ln(1-\hat{y}_i)\right]$$

$$\frac{\partial L}{\partial \hat{y}} = \frac{\hat{y}-y}{\hat{y}(1-\hat{y})}, \qquad \frac{\partial L}{\partial z} = \hat{y} - y$$

The clean gradient w.r.t. the logit is why BCE is paired with sigmoid: the $\sigma'$ factor cancels the denominator, avoiding the learning slowdown that MSE-plus-sigmoid suffers when the output saturates on a wrong prediction.

**Categorical cross-entropy.** With one-hot $\mathbf{y}$ and $\hat{\mathbf{y}} = \text{softmax}(\mathbf{z})$ over $K$ classes:

$$L_{\text{CE}} = -\sum_{k=1}^{K} y_k \ln\hat{y}_k, \qquad \nabla_{\mathbf{z}} L_{\text{CE}} = \hat{\mathbf{y}} - \mathbf{y}$$

Equivalently the KL divergence $D_{\text{KL}}(\mathbf{y}\,\|\,\hat{\mathbf{y}})$ up to the constant entropy of $\mathbf{y}$, and equivalently the negative log-likelihood of a categorical model.

**Hinge loss** (the SVM loss; Cortes & Vapnik, 1995). For $y \in \{-1,+1\}$ and raw score $\hat{y}$:

$$L_{\text{hinge}} = \max\left(0,\; 1 - y\hat{y}\right), \qquad \frac{\partial L}{\partial\hat{y}} = \begin{cases}-y & y\hat{y} < 1\\ 0 & \text{otherwise}\end{cases}$$

Zero loss once the example is correctly classified *with margin at least 1* — unlike cross-entropy, which never reaches exactly zero and keeps pushing confident predictions. The squared hinge $\max(0,1-y\hat{y})^2$ is differentiable everywhere. The multiclass version (Crammer–Singer / Weston–Watkins) is $\sum_{k \neq y}\max(0, \hat{y}_k - \hat{y}_y + 1)$.

---

## 8. Optimizers

Throughout, $\theta_t$ are parameters at step $t$, $g_t = \nabla_\theta L(\theta_t)$ the (mini-batch) gradient, and $\eta$ the learning rate. Operations are elementwise.

**SGD** (Robbins & Monro, 1951, *"A Stochastic Approximation Method"*):

$$\theta_{t+1} = \theta_t - \eta\, g_t$$

Robbins–Monro convergence requires $\sum_t \eta_t = \infty$ and $\sum_t \eta_t^2 < \infty$.

**Momentum / heavy ball** (Boris Polyak, *"Some Methods of Speeding Up the Convergence of Iteration Methods"*, **USSR Comp. Math. and Math. Physics**, 4(5):1–17, 1964):

$$v_{t} = \beta v_{t-1} + g_t, \qquad \theta_{t+1} = \theta_t - \eta\, v_t$$

(equivalently $v_t = \beta v_{t-1} + (1-\beta)g_t$ in the EMA convention, $\beta$ typically $0.9$). The velocity accumulates consistent gradient directions and cancels oscillatory ones; effective step size in a consistent direction is amplified by $1/(1-\beta)$.

**Nesterov Accelerated Gradient** (Yurii Nesterov, 1983, *"A method for solving the convex programming problem with convergence rate $O(1/k^2)$"*). The gradient is evaluated at the *look-ahead* point:

$$v_t = \beta v_{t-1} + \nabla_\theta L\left(\theta_t - \eta\beta v_{t-1}\right), \qquad \theta_{t+1} = \theta_t - \eta v_t$$

The essential difference from Polyak: momentum is applied *before* the gradient is measured, letting the update "see" where it is heading and correct in advance. For smooth convex objectives NAG attains the optimal $O(1/k^2)$ rate versus $O(1/k)$ for plain gradient descent. Deep learning frameworks implement the Sutskever et al. (2013) reparameterization.

**AdaGrad** (Duchi, Hazan & Singer, *"Adaptive Subgradient Methods for Online Learning and Stochastic Optimization"*, **JMLR** 12:2121–2159, 2011):

$$G_t = G_{t-1} + g_t^2, \qquad \theta_{t+1} = \theta_t - \frac{\eta}{\sqrt{G_t} + \epsilon}\odot g_t$$

Per-coordinate learning rates inversely proportional to the accumulated gradient magnitude: rare (sparse) features get large steps, frequent ones small steps. Excellent for sparse, convex problems (NLP with bag-of-words); its flaw in deep learning is that $G_t$ grows monotonically, so the effective learning rate decays to zero and learning stalls.

**RMSProp** (Tieleman & Hinton, Coursera *Neural Networks for Machine Learning*, Lecture 6.5, 2012 — never formally published) replaces the sum with an exponential moving average:

$$E[g^2]_t = \rho\,E[g^2]_{t-1} + (1-\rho)\,g_t^2, \qquad \theta_{t+1} = \theta_t - \frac{\eta}{\sqrt{E[g^2]_t}+\epsilon}\odot g_t$$

with $\rho \approx 0.9$. Old gradients are forgotten, so the effective learning rate no longer collapses. AdaDelta (Zeiler, 2012) is a closely related variant that additionally eliminates $\eta$ via a second accumulator of parameter updates.

**Adam** (Diederik P. Kingma & Jimmy Ba, *"Adam: A Method for Stochastic Optimization"*, ICLR 2015, arXiv:1412.6980) — "adaptive moment estimation," combining momentum (first moment) and RMSProp (second moment) with bias correction:

$$m_t = \beta_1 m_{t-1} + (1-\beta_1)\,g_t$$
$$v_t = \beta_2 v_{t-1} + (1-\beta_2)\,g_t^2$$
$$\hat{m}_t = \frac{m_t}{1-\beta_1^{\,t}}, \qquad \hat{v}_t = \frac{v_t}{1-\beta_2^{\,t}}$$
$$\theta_{t+1} = \theta_t - \eta\,\frac{\hat{m}_t}{\sqrt{\hat{v}_t}+\epsilon}$$

Defaults: $\eta = 0.001$, $\beta_1 = 0.9$, $\beta_2 = 0.999$, $\epsilon = 10^{-8}$. The bias correction matters because $m_0 = v_0 = 0$ biases the early EMAs toward zero; since $\mathbb{E}[m_t] \approx (1-\beta_1^t)\mathbb{E}[g_t]$, dividing by $(1-\beta_1^t)$ removes the bias. Without it, the first steps would be drastically too small — especially for $v$ with $\beta_2 = 0.999$. The ratio $\hat{m}/\sqrt{\hat{v}}$ is a signal-to-noise estimate, and the effective step is bounded by roughly $\eta$ regardless of gradient scale, making Adam invariant to gradient rescaling. Reddi et al. (2018) showed a flaw in the original convergence proof and proposed AMSGrad (using $\max$ of past $v_t$).

**AdamW** (Ilya Loshchilov & Frank Hutter, *"Decoupled Weight Decay Regularization"*, ICLR 2019, arXiv:1711.05101). The key observation: L2 regularization and weight decay are equivalent for plain SGD but **not** for adaptive methods. Adding $\lambda\theta$ to the gradient makes the decay be divided by $\sqrt{\hat v_t}$, so parameters with large historical gradients are decayed *less* — the opposite of the intent. AdamW decouples the two:

$$\theta_{t+1} = \theta_t - \eta_t\left(\frac{\hat{m}_t}{\sqrt{\hat{v}_t}+\epsilon} + \lambda\,\theta_t\right)$$

where $\eta_t$ is a schedule multiplier applied to *both* terms simultaneously (so that the tuned $\lambda$ remains valid under learning-rate schedules). This substantially improves Adam's generalization and decouples the optimal $\lambda$ from the optimal $\eta$; AdamW is the default optimizer for essentially all modern Transformers.

---

## 9. Regularization

**L2 regularization / weight decay / ridge.**

$$L_{\text{total}} = L_{\text{data}} + \frac{\lambda}{2}\sum_{\ell}\left\|\mathbf{W}^{[\ell]}\right\|_F^2, \qquad \frac{\partial L_{\text{total}}}{\partial \mathbf{W}} = \frac{\partial L_{\text{data}}}{\partial\mathbf{W}} + \lambda\mathbf{W}$$

The SGD update becomes $\mathbf{W} \leftarrow (1-\eta\lambda)\mathbf{W} - \eta\,\partial L_{\text{data}}/\partial\mathbf{W}$ — an explicit multiplicative shrinkage, hence "weight decay." From a Bayesian standpoint this is a zero-mean Gaussian prior on the weights (MAP estimation). Biases are conventionally excluded.

**L1 regularization / lasso.**

$$L_{\text{total}} = L_{\text{data}} + \lambda\sum_{\ell}\left\|\mathbf{W}^{[\ell]}\right\|_1, \qquad \frac{\partial}{\partial\mathbf{W}} = \frac{\partial L_{\text{data}}}{\partial\mathbf{W}} + \lambda\operatorname{sign}(\mathbf{W})$$

The constant-magnitude gradient drives small weights exactly to zero, producing **sparse** solutions (feature selection). Corresponds to a Laplace prior. **Elastic net** combines both: $\lambda_1\|\mathbf{W}\|_1 + \frac{\lambda_2}{2}\|\mathbf{W}\|_2^2$.

**Dropout** (Srivastava, Hinton, Krizhevsky, Sutskever & Salakhutdinov, *"Dropout: A Simple Way to Prevent Neural Networks from Overfitting"*, **JMLR** 15:1929–1958, 2014; building on Hinton et al., 2012). During training, each unit is deleted independently with probability $1-p$:

$$r_j^{[\ell]} \sim \text{Bernoulli}(p), \qquad \tilde{\mathbf{a}}^{[\ell]} = \mathbf{r}^{[\ell]}\odot\mathbf{a}^{[\ell]}, \qquad \mathbf{z}^{[\ell+1]} = \mathbf{W}^{[\ell+1]}\tilde{\mathbf{a}}^{[\ell]} + \mathbf{b}^{[\ell+1]}$$

At test time all units are kept and weights are scaled: $\mathbf{W}_{\text{test}} = p\,\mathbf{W}$, so that the expected input to each unit matches training. The practical implementation is **inverted dropout**, which divides by $p$ during training instead:

$$\tilde{\mathbf{a}}^{[\ell]} = \frac{1}{p}\,\mathbf{r}^{[\ell]}\odot\mathbf{a}^{[\ell]} \quad\text{(train)}, \qquad \tilde{\mathbf{a}}^{[\ell]} = \mathbf{a}^{[\ell]}\quad\text{(test)}$$

leaving inference untouched. Typical $p = 0.5$ for hidden layers, $0.8$ for inputs. The mechanism prevents **co-adaptation** — no unit can rely on the presence of any particular other unit — and can be interpreted as training an exponential ensemble of $2^N$ thinned subnetworks with shared weights, approximately averaged at test time by the geometric mean.

**Batch Normalization** (Sergey Ioffe & Christian Szegedy, *"Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift"*, ICML 2015, arXiv:1502.03167). For each feature, over a mini-batch $\mathcal{B} = \{x_1,\dots,x_m\}$:

$$\mu_{\mathcal{B}} = \frac{1}{m}\sum_{i=1}^{m} x_i \qquad\text{(mini-batch mean)}$$
$$\sigma^2_{\mathcal{B}} = \frac{1}{m}\sum_{i=1}^{m}\left(x_i - \mu_{\mathcal{B}}\right)^2 \qquad\text{(mini-batch variance)}$$
$$\hat{x}_i = \frac{x_i - \mu_{\mathcal{B}}}{\sqrt{\sigma^2_{\mathcal{B}}+\epsilon}} \qquad\text{(normalize)}$$
$$y_i = \gamma\hat{x}_i + \beta \equiv \text{BN}_{\gamma,\beta}(x_i) \qquad\text{(scale and shift)}$$

The learned $\gamma,\beta$ restore representational capacity — setting $\gamma = \sqrt{\sigma^2_\mathcal{B}+\epsilon}$, $\beta=\mu_\mathcal{B}$ recovers the identity, so BN never *removes* expressiveness. Backpropagation must flow through $\mu_\mathcal{B}$ and $\sigma^2_\mathcal{B}$ as well:

$$\frac{\partial L}{\partial \hat x_i} = \frac{\partial L}{\partial y_i}\gamma, \quad \frac{\partial L}{\partial\sigma^2_\mathcal{B}} = \sum_i \frac{\partial L}{\partial \hat x_i}(x_i-\mu_\mathcal{B})\cdot\frac{-1}{2}(\sigma^2_\mathcal{B}+\epsilon)^{-3/2}$$
$$\frac{\partial L}{\partial\mu_\mathcal{B}} = \sum_i\frac{\partial L}{\partial\hat x_i}\cdot\frac{-1}{\sqrt{\sigma^2_\mathcal{B}+\epsilon}} + \frac{\partial L}{\partial\sigma^2_\mathcal{B}}\cdot\frac{-2\sum_i(x_i-\mu_\mathcal{B})}{m}$$
$$\frac{\partial L}{\partial x_i} = \frac{\partial L}{\partial\hat x_i}\frac{1}{\sqrt{\sigma^2_\mathcal{B}+\epsilon}} + \frac{\partial L}{\partial\sigma^2_\mathcal{B}}\frac{2(x_i-\mu_\mathcal{B})}{m} + \frac{\partial L}{\partial\mu_\mathcal{B}}\frac{1}{m}$$
$$\frac{\partial L}{\partial\gamma} = \sum_i\frac{\partial L}{\partial y_i}\hat x_i, \qquad \frac{\partial L}{\partial\beta} = \sum_i\frac{\partial L}{\partial y_i}$$

At **inference**, batch statistics are unavailable, so running (EMA) estimates collected during training are used: $\mathbb{E}[x] \leftarrow \text{EMA}(\mu_\mathcal{B})$, $\mathrm{Var}[x] \leftarrow \frac{m}{m-1}\text{EMA}(\sigma^2_\mathcal{B})$, and the whole transform becomes a fixed affine map foldable into the preceding convolution. BN permits much larger learning rates, reduces sensitivity to initialization, and has a mild regularizing effect from mini-batch noise. Its explanation via "internal covariate shift" was later challenged — Santurkar et al. (2018) argue the real benefit is a smoother loss landscape. Its weaknesses: dependence on batch size (poor for $m \lesssim 8$) and awkwardness in recurrent and online settings.

**Layer Normalization** (Jimmy Lei Ba, Jamie Ryan Kiros & Geoffrey E. Hinton, *"Layer Normalization"*, arXiv:1607.06450, 2016) transposes the computation: statistics are taken over the **features of a single example**, not over the batch. For $\mathbf{x}\in\mathbb{R}^d$:

$$\mu = \frac{1}{d}\sum_{i=1}^{d}x_i, \qquad \sigma^2 = \frac{1}{d}\sum_{i=1}^{d}(x_i-\mu)^2$$
$$\text{LN}(\mathbf{x}) = \boldsymbol{\gamma}\odot\frac{\mathbf{x}-\mu}{\sqrt{\sigma^2+\epsilon}} + \boldsymbol{\beta}$$

Because it is batch-independent, LN behaves identically at train and test time, works with batch size 1, and handles variable-length sequences — which is why it, not BN, is used in RNNs and Transformers. Related variants: **GroupNorm** (Wu & He, 2018) normalizes over feature groups; **RMSNorm** (Zhang & Sennrich, 2019) drops the mean-centring, $\text{RMSNorm}(\mathbf{x}) = \boldsymbol{\gamma}\odot\mathbf{x}/\sqrt{\frac{1}{d}\sum_i x_i^2 + \epsilon}$, and is now standard in LLaMA-class models.

Other standard regularizers: **early stopping** (implicit L2 for quadratic objectives), **data augmentation**, **label smoothing** ($y_k \leftarrow (1-\varepsilon)y_k + \varepsilon/K$), and **max-norm constraints** $\|\mathbf{w}_j\|_2 \leq c$, which the dropout paper recommends alongside dropout.

---

## 10. Radial Basis Function (RBF) Networks

RBF networks were introduced by **D. S. Broomhead & D. Lowe**, *"Multivariable Functional Interpolation and Adaptive Networks"*, **Complex Systems**, 2:321–355, 1988 (RSRE Memorandum 4148), with independent contemporaneous work by Lee & Kil (1988) and Niranjan & Fallside (1988), and refined by **J. Moody & C. Darken**, *"Fast Learning in Networks of Locally-Tuned Processing Units"*, **Neural Computation**, 1(2):281–294, 1989. **Michael Powell**'s earlier work on radial basis function interpolation (1985–1987) supplied the mathematical basis, and **Park & Sandberg** (1991) proved universal approximation for RBF networks.

### Architecture

Strictly three layers, with only one hidden layer:

1. **Input layer** — $n$ nodes, pure fan-out.
2. **Hidden layer** — $J$ *locally tuned* units, each with a centre $\boldsymbol{\mu}_j \in \mathbb{R}^n$ and a width $\sigma_j$, computing a radial function of the distance between input and centre.
3. **Output layer** — $L$ **linear** units.

The output is a weighted sum of basis functions:

$$f_l(\mathbf{x}) = \sum_{j=1}^{J} w_{lj}\,\varphi\!\left(\left\|\mathbf{x} - \boldsymbol{\mu}_j\right\|\right) + b_l$$

The most common kernel is the **Gaussian**:

$$\varphi_j(\mathbf{x}) = \exp\!\left(-\frac{\left\|\mathbf{x}-\boldsymbol{\mu}_j\right\|^2}{2\sigma_j^2}\right)$$

often written $\exp(-\beta_j\|\mathbf{x}-\boldsymbol{\mu}_j\|^2)$ with $\beta_j = 1/(2\sigma_j^2)$. Anisotropic units generalize to a Mahalanobis form $\exp\!\left(-\tfrac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_j)^\top\boldsymbol{\Sigma}_j^{-1}(\mathbf{x}-\boldsymbol{\mu}_j)\right)$. Other classical kernels: multiquadric $\sqrt{r^2+c^2}$, inverse multiquadric $1/\sqrt{r^2+c^2}$, and thin-plate spline $r^2\ln r$.

### The essential contrast with the MLP

An MLP hidden unit computes an **inner product** $\mathbf{w}^\top\mathbf{x}$ and responds along a hyperplane — a *global*, distributed representation. An RBF unit computes a **distance** $\|\mathbf{x}-\boldsymbol{\mu}\|$ and responds only in a localized neighbourhood — a *local*, hypersphere-shaped receptive field. Consequences: RBF networks train much faster (the output layer is linear), interpolate cleanly, and degrade gracefully outside the data (activations vanish, so the output tends to the bias); but they suffer the curse of dimensionality, since covering a high-dimensional space with local bumps requires exponentially many centres.

### Training

The standard procedure is **two-stage / hybrid**:

1. **Unsupervised** placement of centres $\boldsymbol{\mu}_j$: random subsampling of the training data or a coarse lattice (Broomhead & Lowe, 1988), or **$k$-means clustering** (Moody & Darken, 1989). Widths are then set heuristically, e.g. $\sigma_j = d_{\max}/\sqrt{2J}$ with $d_{\max}$ the maximum inter-centre distance, or by the $p$-nearest-neighbour rule $\sigma_j = \left(\frac{1}{p}\sum_{k=1}^{p}\|\boldsymbol{\mu}_j-\boldsymbol{\mu}_k\|^2\right)^{1/2}$ (Moody & Darken).

2. **Supervised** solution of the output weights. Since the model is linear in $\mathbf{W}$, the least-squares solution is closed-form. With the design (interpolation) matrix $\Phi_{ij} = \varphi_j(\mathbf{x}_i)$:

$$\mathbf{W} = \left(\boldsymbol{\Phi}^\top\boldsymbol{\Phi} + \lambda\mathbf{I}\right)^{-1}\boldsymbol{\Phi}^\top\mathbf{Y} = \boldsymbol{\Phi}^{+}\mathbf{Y}$$

using the Moore–Penrose pseudoinverse (Broomhead & Lowe's approach), with $\lambda$ a ridge term. In the *exact interpolation* case $J = m$ (one centre per data point), Michelli's theorem guarantees $\boldsymbol{\Phi}$ is nonsingular for Gaussian kernels and distinct points, giving $\mathbf{W}=\boldsymbol{\Phi}^{-1}\mathbf{Y}$ — but this overfits, hence the use of $J \ll m$.

All parameters $\{w_{lj}, \boldsymbol{\mu}_j, \sigma_j\}$ can alternatively be trained jointly by gradient descent, e.g. $\partial L/\partial\boldsymbol{\mu}_j = \sum_l \delta_l w_{lj}\varphi_j(\mathbf{x})\frac{\mathbf{x}-\boldsymbol{\mu}_j}{\sigma_j^2}$, at the cost of the convexity that makes the hybrid scheme attractive.

RBF networks are close relatives of kernel methods (an SVM with a Gaussian kernel is an RBF network whose centres are the support vectors and whose weights come from the dual QP), of Gaussian mixture models, and of normalized-RBF / Nadaraya–Watson regression. They remain in use for function interpolation, meshless PDE solvers, time-series prediction, and control.

---

## Chronological Summary of Founding Papers

| Year | Authors | Contribution |
|---|---|---|
| 1943 | McCulloch & Pitts | *A Logical Calculus of the Ideas Immanent in Nervous Activity* — threshold neuron |
| 1949 | Hebb | *The Organization of Behavior* — Hebbian learning |
| 1958 | Rosenblatt | *The Perceptron: A Probabilistic Model…* — first learning rule |
| 1960 | Widrow & Hoff | *Adaptive Switching Circuits* — ADALINE, LMS/delta rule |
| 1962 | Novikoff | *On Convergence Proofs for Perceptrons* — mistake bound |
| 1964 | Polyak / Huber | Heavy-ball momentum / robust loss |
| 1969 | Minsky & Papert | *Perceptrons* — XOR and the limits of linear separability |
| 1974 | Werbos | PhD thesis — backpropagation (reverse-mode AD) |
| 1983 | Nesterov | Accelerated gradient, $O(1/k^2)$ |
| 1986 | Rumelhart, Hinton & Williams | *Learning Representations by Back-Propagating Errors*, **Nature** |
| 1988–89 | Broomhead & Lowe; Moody & Darken | RBF networks |
| 1989 | Cybenko; Hornik, Stinchcombe & White | Universal approximation |
| 1993 | Leshno, Lin, Pinkus & Schocken | Universal approximation iff non-polynomial |
| 2010–11 | Nair & Hinton; Glorot, Bordes & Bengio | ReLU for deep networks |
| 2011 | Duchi, Hazan & Singer | AdaGrad |
| 2014 | Kingma & Ba; Srivastava et al. | Adam; Dropout |
| 2015 | Ioffe & Szegedy; He et al.; Clevert et al. | BatchNorm; PReLU + He init; ELU |
| 2016 | Ba, Kiros & Hinton; Hendrycks & Gimpel | LayerNorm; GELU |
| 2017 | Klambauer et al.; Ramachandran et al. | SELU; Swish |
| 2019 | Loshchilov & Hutter; Misra | AdamW; Mish |

Sources: [A Logical Calculus (Wikipedia)](https://en.wikipedia.org/wiki/A_Logical_Calculus_of_the_Ideas_Immanent_in_Nervous_Activity), [Perceptrons (book)](https://en.wikipedia.org/wiki/Perceptrons_(book)), [Widrow & Hoff / ADALINE](https://en.wikipedia.org/wiki/Bernard_Widrow), [Hornik, Stinchcombe & White 1989](https://www.cs.cmu.edu/~epxing/Class/10715/reading/Kornick_et_al.pdf), [Hornik 1991](https://web.njit.edu/~usman/courses/cs677/hornik-nn-1991.pdf), [Note on Cybenko's UAT](https://arxiv.org/html/2508.18893v1), [Rumelhart, Hinton & Williams 1986 (Nature)](https://www.nature.com/articles/323533a0), [Activation function reference](https://en.wikipedia.org/wiki/Activation_function), [Self-Normalizing Neural Networks (SELU)](https://arxiv.org/pdf/1706.02515), [Delving Deep into Rectifiers (PReLU)](https://arxiv.org/pdf/1502.01852), [Rectifier (ReLU) history](https://en.wikipedia.org/wiki/ReLU), [Swish function](https://en.wikipedia.org/wiki/Swish_function), [Batch Normalization (Ioffe & Szegedy)](https://arxiv.org/abs/1502.03167), [Dropout (Srivastava et al. 2014)](https://nitishsrivastava.github.io/publication/2014-01-01), [Layer Normalization (Ba, Kiros & Hinton)](https://www.semanticscholar.org/paper/Layer-Normalization-Ba-Kiros/97fb4e3d45bb098e27e0071448b6152217bd35a5), [Decoupled Weight Decay Regularization (AdamW)](https://arxiv.org/pdf/1711.05101), [AdaGrad (Cornell Optimization Wiki)](https://optimization.cbe.cornell.edu/index.php?title=AdaGrad), [Novikoff perceptron convergence proof](https://apps.dtic.mil/sti/tr/pdf/AD0298258.pdf), [Leshno et al. 1993](https://www.sciencedirect.com/science/article/abs/pii/S0893608005801315), [RBF networks (MIT book, ch. 6)](https://neuron.eng.wayne.edu/tarek/MITbook/chap6/6_1.html), [Broomhead & Lowe 1988](https://www.sciepub.com/reference/93721)