Mathematics in Computing: The Maths Behind Modern AI

Every computer program is mathematics in some form. Binary numbers are place-value arithmetic with base two, logic gates are Boolean algebra, a dependency graph is graph theory, and an index lookup is a statement about logarithms. But when people say “the maths of AI”, they usually mean five specific ideas that sit underneath every neural network:

1. Linear algebra: vectors, matrices and tensors, the containers that hold data and weights.
2. Calculus: derivatives, partial derivatives and gradients, which measure how a change in one number affects another.
3. Gradient descent and backpropagation: the method that uses those gradients to make a network learn.
4. Probability distributions: the language for uncertainty, noise and the outputs of models.
5. Bayes’ theorem and statistical inference: how to update beliefs from evidence, and what it means to “fit” a model to data.

This article explains each one from the beginning, in plain words first and symbols second, with short Python snippets you can run. All the code in it was run before publishing, and you only need numpy.

Part 1: Linear algebra, the containers of AI

Scalars, vectors, matrices and tensors

These four words describe the same idea at increasing numbers of dimensions.

• A scalar is a single number: 3.5.
• A vector is an ordered list of numbers. It can be a point in space, but in AI it is more often a description: a house as [area, bedrooms, age], a pixel as [red, green, blue], a word as a list of several hundred numbers capturing its meaning.
• A matrix is a grid of numbers with rows and columns. It can hold a table of data (one row per example), or a set of weights that transform one vector into another.
• A tensor is the general term for an array with any number of dimensions. A scalar is a 0-dimensional tensor, a vector is 1-D, a matrix is 2-D, and a colour image is 3-D (height × width × channels).

The number of dimensions is the tensor’s rank (or ndim), and the size along each dimension is its shape. Shapes are the single most important thing to track when you work with neural networks, because most bugs are shape mismatches.

import numpy as np

t = np.zeros((32, 128, 768))
print(t.ndim, t.shape, t.size)   # 3 (32, 128, 768) 3145728

That shape is typical of a language model: a batch of 32 sentences, each with 128 tokens, each token a vector of 768 numbers. Frameworks such as PyTorch and TensorFlow are named for tensors because everything inside them is one.

Vector operations

Two vectors of the same length can be added element by element, and a vector can be multiplied by a scalar to stretch it.

a = np.array([1.0, 2.0, 3.0])
b = np.array([4.0, 5.0, 6.0])
print(a + b)    # [5. 7. 9.]
print(2 * a)    # [2. 4. 6.]

The most important operation is the dot product: multiply matching elements and add them up. For a and b above that is 1×4 + 2×5 + 3×6 = 32.

print(a @ b)    # 32.0

The dot product measures how much two vectors agree. It is large and positive when they point the same way, zero when they are perpendicular (unrelated), and negative when they oppose. Almost every operation in AI is a dot product in disguise: a neuron computes a dot product between its weights and its inputs.

The length (or norm) of a vector is sqrt(a @ a), written np.linalg.norm(a), here 3.74. Dividing the dot product by both lengths gives cosine similarity, the cosine of the angle between the vectors:

def cosine(u, v):
    return (u @ v) / (np.linalg.norm(u) * np.linalg.norm(v))

print(round(cosine(a, b), 4))    # 0.9746  (nearly the same direction)
print(round(cosine(a, -a), 4))   # -1.0    (exact opposites)

This is how semantic search works. Documents and queries are converted into vectors called embeddings, and “similar meaning” becomes “high cosine similarity”. Recommendation engines, duplicate detection and retrieval-augmented generation all rest on it.

Matrix operations

A matrix multiplication is many dot products at once. To compute A @ B, take each row of A and dot it with each column of B. Two rules matter most:

• Shapes must match: an (m × n) matrix times an (n × p) matrix gives an (m × p) matrix. The inner sizes must be equal.
• Order matters: A @ B is generally not B @ A.

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 1], [1, 0]])
print(A @ B)     # [[2 1] [4 3]]   matrix product
print(A * B)     # [[0 2] [3 0]]   element-wise product (a different thing)
print(A.T)       # [[1 3] [2 4]]   transpose: swap rows and columns

Be careful to keep @ (matrix product) and * (element-wise) apart. A transpose flips a matrix over its diagonal, and you will see .T constantly when shapes need to line up. Broadcasting lets numpy stretch a small array across a larger one, for example adding a bias vector to every row of a batch.

A neural network layer is a matrix multiplication

Here is the central fact. A fully connected layer takes an input vector x, multiplies it by a weight matrix W, and adds a bias vector b:

W = np.array([[0.5, -1.0,  2.0],
              [1.5,  0.3, -0.7]])    # 2 outputs from 3 inputs
b = np.array([0.1, -0.2])
x = np.array([1.0, 2.0, 3.0])

print(W @ x + b)                      # [ 4.6 -0.2]

Row 1 of W is the weights of the first neuron, so its output is 0.5×1 − 1.0×2 + 2.0×3 + 0.1 = 4.6. In practice you process many examples at once by stacking them into a matrix X of shape (batch, features), and one matrix product handles the whole batch:

X = np.random.rand(32, 3)             # 32 examples, 3 features each
print((X @ W.T + b).shape)            # (32, 2)

A large language model is, to a first approximation, billions of such multiplications chained together with simple non-linear functions between them. That is why GPUs, which were built to multiply matrices for graphics, became the engine of AI: rendering a 3D scene and running a neural network are both enormous amounts of linear algebra done in parallel. The cost of running a model is really a count of multiplications, and the gigabytes of matrices that must move through memory.

Part 2: Calculus, the mathematics of change

The derivative

Calculus studies how things change. The derivative of a function at a point is its slope there: how much the output changes for a tiny nudge to the input. If f(x) = x², the derivative is 2x. At x = 3 the slope is 6, which means that increasing x by a small amount h increases f by about 6h.

You can check a derivative numerically, with no calculus at all, by nudging the input both ways:

f = lambda x: x ** 2
h = 1e-6
print((f(3 + h) - f(3 - h)) / (2 * h))    # 6.000000000838668

This numerical derivative is slow, but it is a reliable way to test a derivative you worked out by hand, and Part 3 uses it to test backpropagation. The sign of the derivative tells you the direction: positive means the function rises as x rises, negative means it falls. To minimise a function, move x in the direction that makes the slope smaller, which means stepping against the sign of the derivative.

A few rules cover most of what AI needs:

Function	Derivative
constant c	0
x^n	n·x^(n-1)
e^x	e^x
ln x	1/x
f + g	f' + g'
f(g(x))	f'(g(x)) · g'(x) (the chain rule)

Partial derivatives and gradients

A neural network has not one input but millions of numbers to adjust. When a function depends on several variables, the partial derivative with respect to one variable is its slope when you nudge only that variable and hold the others still. Take g(x, y) = x²y + 3y. The partial derivative with respect to x is 2xy, and with respect to y it is x² + 3.

The gradient is the vector that collects all the partial derivatives, written ∇g. At the point (2, 1) it is (2·2·1, 2² + 3) = (4, 7). Check it numerically:

g = lambda x, y: x**2 * y + 3 * y

def grad(x, y, h=1e-6):
    dx = (g(x + h, y) - g(x - h, y)) / (2 * h)
    dy = (g(x, y + h) - g(x, y - h)) / (2 * h)
    return dx, dy

print(grad(2, 1))    # (4.0000000001, 6.9999999996)

The gradient has a beautiful property: it points in the direction in which the function increases fastest, and its length says how steep that climb is. So the negative gradient points straight downhill. A network with a billion weights has a gradient with a billion entries, one per weight, saying how much each weight affects the error. Training is the art of following that gradient downhill.

The chain rule

A network is a function inside a function inside a function, so we need the derivative of a composition. The chain rule says: multiply the derivatives along the chain. If y = f(u) and u = g(x), then dy/dx = (dy/du) · (du/dx). Intuition: if u changes twice as fast as x, and y changes three times as fast as u, then y changes six times as fast as x. Backpropagation is nothing more than applying this rule systematically, which Part 3 shows.

Part 3: Gradient descent and backpropagation

Loss, and the idea of descent

To train a model you first need a number that says how wrong it is. This is the loss (or cost). For predicting a number, a common loss is the mean squared error, the average of the squared differences between predictions and targets. For classification it is cross-entropy (Part 5). A model’s loss depends on its weights, so you can picture a landscape: each point is one setting of the weights, and its height is the loss. Training means finding the low ground.

Gradient descent does this with one rule, repeated:

weight = weight - learning_rate × (slope of the loss with respect to that weight)

Move each weight a small step downhill, recompute, and repeat. Here it is on a two-variable bowl whose lowest point is at (3, −1):

w = np.array([0.0, 0.0])
for step in range(100):
    gradient = np.array([2 * (w[0] - 3), 2 * (w[1] + 1)])
    w = w - 0.1 * gradient
print(w.round(3))    # [ 3. -1.]

The learning rate (here 0.1) is the size of each step, and choosing it is the most common tuning decision in AI. Minimising x² from x = 5 for 30 steps shows the three behaviours:

for lr in (0.01, 0.1, 1.1):
    x = 5.0
    for _ in range(30):
        x -= lr * 2 * x
    print(lr, x)
# 0.01 2.73       too small: barely moved
# 0.1  0.006      good: nearly at the minimum, 0
# 1.1  1186.9     too large: overshoots and blows up

In real training you rarely compute the gradient on the whole dataset. Stochastic gradient descent uses a random mini-batch (say 32 or 256 examples) each step, which is noisier but far cheaper, and the noise often helps the search. Optimisers such as momentum and Adam add refinements, such as remembering the recent direction and adapting the step size for each weight, but every one of them is gradient descent underneath.

Backpropagation

Gradient descent needs the gradient of the loss with respect to every weight, in a network with many layers. Backpropagation computes them all efficiently using the chain rule, in two passes:

1. Forward pass: run the input through the layers, saving each intermediate result, and compute the loss.
2. Backward pass: start at the loss and work backward layer by layer. At each layer, use the gradient arriving from above and that layer’s local derivative to produce (a) the gradient for the layer’s own weights and (b) the gradient to hand to the layer below.

The saving is enormous. Each local derivative is computed once and reused, so the cost of the backward pass is about the same as the forward pass, not a million times more. Modern frameworks do this automatically (automatic differentiation), but seeing it by hand once removes the mystery.

A complete network, by hand

Here is a tiny network that learns XOR, a function that a single layer provably cannot represent: output 1 when exactly one of two inputs is 1. It has two inputs, four hidden neurons with tanh, and one output with a sigmoid, trained with cross-entropy loss. The derivative formulas used below are worked out in the comments.

import numpy as np

X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], float)
y = np.array([[0], [1], [1], [0]], float)
sigmoid = lambda z: 1 / (1 + np.exp(-z))

rng = np.random.default_rng(1)
W1 = rng.normal(0, 1, (2, 4)); b1 = np.zeros((1, 4))
W2 = rng.normal(0, 1, (4, 1)); b2 = np.zeros((1, 1))

for step in range(5000):
    # ---- forward pass ----
    a1 = np.tanh(X @ W1 + b1)
    p = sigmoid(a1 @ W2 + b2)
    loss = -np.mean(y * np.log(p) + (1 - y) * np.log(1 - p))

    # ---- backward pass (chain rule, from the loss backwards) ----
    dz2 = (p - y) / len(X)          # sigmoid + cross-entropy simplifies to p - y
    dW2 = a1.T @ dz2                # gradient for the output weights
    db2 = dz2.sum(0, keepdims=True)
    dz1 = (dz2 @ W2.T) * (1 - a1**2)   # pass back, times tanh'(z) = 1 - tanh(z)^2
    dW1 = X.T @ dz1
    db1 = dz1.sum(0, keepdims=True)

    # ---- gradient descent update ----
    W1 -= dW1; b1 -= db1
    W2 -= dW2; b2 -= db2

print(round(loss, 4), p.round(2).ravel())    # 0.0005 [0. 1. 1. 0.]

The loss falls from about 0.67 to 0.0005 and the four outputs become [0, 1, 1, 0], exactly XOR. Four things are worth noticing:

• Every gradient is a matrix product of saved forward values with the incoming gradient. Linear algebra and calculus meet here.
• The shapes of gradients always match the shapes of the weights they update (dW1 has shape (2, 4), like W1).
• You can check a backpropagation gradient against the numerical derivative from Part 2. Nudging W1[0, 0] by ±1e-5 and recomputing the loss gives the same value as the analytic dW1[0, 0], 3.2313e-05, to eight digits. This “gradient check” is how you debug a hand-written network.
• It is not guaranteed to succeed. With a different random starting point (seed 0 in the same code), this network gets stuck with outputs [0, 0.5, 1, 0.5] on a flat part of the landscape. Initialisation matters, which is why real frameworks use careful schemes and why the same training run can behave differently with a different random seed.

Deep networks add one more problem: multiplying many small derivatives makes gradients vanish in early layers, and multiplying large ones makes them explode. ReLU activations, residual connections, normalisation layers and gradient clipping all exist to keep these chain-rule products well behaved.

Part 4: Probability distributions

Computers are deterministic, but data is noisy and models are uncertain. Probability is a number from 0 (impossible) to 1 (certain), and a probability distribution says how that probability is spread across the possible outcomes. Two families matter: discrete (countable outcomes, such as a dice roll or a word) and continuous (any value in a range, such as a height or a latency).

Two summaries describe most distributions. The expected value (mean) is the average outcome weighted by probability, and the variance (its square root is the standard deviation) measures how spread out the outcomes are.

The distributions you will meet most

• Bernoulli: one yes/no event with probability p. Did this user click? Is this email spam?
• Binomial: the number of successes in n independent Bernoulli trials. If each request fails with probability 0.2, the chance of exactly 3 failures in 10 requests is C(10,3) · 0.2³ · 0.8⁷:

from math import comb
print(round(comb(10, 3) * 0.2**3 * 0.8**7, 4))    # 0.2013

• Poisson: the number of events in a fixed interval when they occur independently at an average rate λ, such as requests per second or errors per hour. With an average of 4 events per interval, the chance of exactly 2 is e^(-4) · 4² / 2!:

from math import exp, factorial
print(round(exp(-4) * 4**2 / factorial(2), 4))     # 0.1465

• Uniform: every value in a range equally likely, the basis of random number generation.
• Exponential: the waiting time between Poisson events, such as the gap between requests.
• Normal (Gaussian): the bell curve, defined by a mean μ and standard deviation σ. It appears so often because the sum of many small independent effects tends toward it (the central limit theorem). About 68% of values fall within one standard deviation of the mean, 95% within two.

rng = np.random.default_rng(1)
s = rng.normal(100, 15, 100000)
print(s.mean().round(1), s.std().round(1))            # 99.9 14.9
print(((s > 85) & (s < 115)).mean().round(3))         # 0.684  (the 68% rule)
print(np.percentile(s, [50, 95, 99]).round(1))        # [ 99.9 124.6 134.8]

Percentiles are worth a mention. Real latency data is not bell-shaped but has a long tail: if the mean response time is 120 ms but the 99th percentile is 3 seconds, a page that makes fifty backend calls will hit a slow one most of the time. That is why engineers watch p95 and p99 rather than averages.

Where distributions appear in AI

• Weight initialisation draws starting weights from a normal or uniform distribution with a carefully chosen spread.
• Noise and regularisation: dropout randomly switches neurons off, and diffusion image models are trained by adding Gaussian noise and learning to remove it.
• Model outputs are distributions. A classifier ends with the softmax function, which turns raw scores (logits) into a categorical distribution: positive numbers that add up to 1.

def softmax(z, T=1.0):
    z = np.array(z) / T
    e = np.exp(z - z.max())          # subtract the max for numerical stability
    return e / e.sum()

z = [2.0, 1.0, 0.1]
print(softmax(z).round(3))           # [0.659 0.242 0.099]
print(softmax(z, T=0.5).round(3))    # [0.864 0.117 0.019]  sharper
print(softmax(z, T=2.0).round(3))     # [0.502 0.304 0.194]  flatter

A language model produces exactly this kind of distribution over its whole vocabulary at every step, then samples the next token from it. The T above is the temperature setting: low values make the model nearly always pick the most likely token, high values spread the choice and make output more varied.

Part 5: Bayes’ theorem and statistical inference

Bayes’ theorem

Bayes’ theorem tells you how to update a belief when new evidence arrives:

P(H | E) = P(E | H) × P(H) / P(E)

In words: the posterior (belief in hypothesis H after seeing evidence E) equals the likelihood (how probable the evidence is if H is true) times the prior (belief before the evidence), divided by the overall probability of the evidence.

A classic worked example shows why it matters. Suppose 1% of servers have a particular fault. Your monitoring check catches 95% of real faults but also raises a false alarm on 5% of healthy servers. An alert fires. How likely is the server to be faulty?

• Faulty and alert: 0.01 × 0.95 = 0.0095
• Healthy and alert: 0.99 × 0.05 = 0.0495

print(0.0095 / (0.0095 + 0.0495))    # 0.161

Only about 16%, despite a 95%-accurate check. The fault is rare, so most alerts are false alarms. Ignoring the prior, called the base-rate fallacy, is one of the most common errors in reasoning about tests, from spam filters to medical screening to fraud detection, and it is the mathematical root of alert fatigue.

Inference: learning about parameters from data

Statistical inference runs the logic backwards: given observed data, what can we say about the unknown quantities that produced it? Machine learning is inference at scale. Three ideas cover most of it.

1. Maximum likelihood estimation (MLE). Choose the parameters that make the observed data most probable. For a normal distribution, the MLE of the mean is the sample mean and the MLE of the spread is the sample standard deviation:

d = rng.normal(5, 2, 1000)
print(d.mean().round(2), d.std().round(2))    # 5.09 2.08  (true values: 5 and 2)

Here is the link to Part 3. In practice we maximise the log-likelihood, which is equivalent to minimising its negative, the negative log-likelihood, and this is the cross-entropy loss. If a model gives the true class probability 0.9, the loss is −ln 0.9 = 0.105. If it gives only 0.1, the loss is −ln 0.1 = 2.303. So training a classifier by gradient descent on cross-entropy is performing maximum likelihood estimation. Likewise, mean squared error is the negative log-likelihood of a model that assumes Gaussian noise. The loss functions are not arbitrary: each one encodes a probability assumption about the data.

2. Bayesian updating. Instead of one best estimate, keep a whole distribution over the unknown value and update it as data arrives. Suppose you want the success rate p of a new feature. Begin with a flat Beta(1, 1) prior, meaning “no idea”. After seeing 7 successes and 3 failures, the posterior is simply Beta(1 + 7, 1 + 3) = Beta(8, 4), with mean 8/12 ≈ 0.667. Each new observation just adds one to a count, and the distribution narrows as evidence accumulates.

a, b = 1, 1                 # prior: Beta(1, 1)
successes, failures = 7, 3
a += successes; b += failures
print(a / (a + b))          # 0.667, the posterior mean

A prior is where you write down what you already know. It also shows up in deep learning under another name: L2 regularisation (weight decay) is equivalent to placing a Gaussian prior on the weights, which gently prefers small ones. That approach is called maximum a posteriori (MAP) estimation, MLE plus a prior.

3. Uncertainty and confidence intervals. An estimate without an error bar is half an answer. If 248 of 400 users liked a change (62%), the standard error is sqrt(p(1−p)/n), and an approximate 95% confidence interval is p ± 1.96 × standard error:

import math
p, n = 0.62, 400
se = math.sqrt(p * (1 - p) / n)
print(round(p - 1.96 * se, 3), round(p + 1.96 * se, 3))    # 0.572 0.668

So the true rate is plausibly anywhere from about 57% to 67%. The same reasoning drives A/B tests, hypothesis tests and the question every canary release asks: is the difference between old and new version real, or just noise?

Putting it all together

The five topics are one pipeline:

1. Linear algebra stores the data and weights as tensors and computes the forward pass as matrix products.
2. Probability turns the final scores into a distribution, and the loss measures how surprised that distribution was by the truth (negative log-likelihood).
3. Calculus gives the gradient of that loss with respect to every weight, through the chain rule.
4. Backpropagation computes all those gradients in one backward pass, and gradient descent nudges the weights downhill.
5. Bayes and inference tell you how much to trust what the trained model says, and how to compare one model with another.

A study path

1. Vectors, dot products and matrix multiplication, with numpy. Get fluent in reading shapes before anything else.
2. Derivatives and the chain rule, then partial derivatives and gradients. Check every result with a numerical derivative.
3. Gradient descent on small functions, varying the learning rate. Then write the XOR network above from memory.
4. Distributions: normal, binomial, Poisson, softmax. Simulate each one and plot a histogram.
5. Bayes’ theorem and MLE. Work three base-rate problems by hand, then derive why cross-entropy is a negative log-likelihood.

You do not need every theorem. You need the intuition for what a matrix product, a gradient, a distribution and a posterior each do, and the habit of testing small examples in code before trusting a library. Frameworks change every few years, but this mathematics will outlast every one of them.