The calculus and linear algebra behind science, engineering and machine learning: derivatives, integrals, vectors, matrices, gradients and eigenvalues, with Python exercises.

Calculus & Linear Algebra for STEM

## Why this matters

Most quantities depend on more than one input. The stress in a beam depends on load and length, the yield of a reaction on temperature and pressure, an insurance price on age and coverage, and a machine learning model's error on thousands of weights. To ask "how does the output respond if I change one input?", or "which way should I move all the inputs to improve the output fastest?", you need **partial derivatives** and the **gradient**.

## Partial derivatives

A partial derivative measures the rate of change with respect to one variable while holding all the others fixed. It is written with a curly ∂:

```text
f(x, y) = x²y + 3y
∂f/∂x = 2xy          (treat y as a constant)
∂f/∂y = x² + 3       (treat x as a constant)
```

Every rule from Module 1 still applies; you just pretend the other variables are numbers. At the point (1, 2): ∂f/∂x = 4 and ∂f/∂y = 4.

Worked example from physics: the ideal gas pressure is P = nRT / V. Holding n and V fixed, ∂P/∂T = nR / V (pressure rises with temperature). Holding T fixed, ∂P/∂V = −nRT / V² (pressure falls as volume grows).

## The gradient

Collect all the partial derivatives into a vector and you have the **gradient**:

```text
∇f = (∂f/∂x, ∂f/∂y, …)
```

For f(x, y) = x²y + 3y at (1, 2), ∇f = (4, 4). The gradient has two key properties:

1. It points in the direction of **steepest increase** of f.
2. Its length is the rate of increase in that direction.

Picture a hilly landscape where f is the height. Standing anywhere, the gradient is an arrow pointing straight uphill. At the top of a hill, or the bottom of a valley, the ground is flat and the gradient is the zero vector, which is the multivariable version of f′(x) = 0.

## Gradient descent

To minimise a function of many variables, repeatedly step **against** the gradient:

```text
new point = point − learning_rate × ∇f(point)
```

This is **gradient descent**, the algorithm that trains linear regression, logistic regression and neural networks. Example: f(x, y) = (x − 3)² + (y + 1)² has ∇f = (2(x − 3), 2(y + 1)). Starting at (0, 0) with a learning rate of 0.1, the first step goes to (0, 0) − 0.1·(−6, 2) = (0.6, −0.2), heading toward the minimum at (3, −1). Too large a learning rate overshoots and can diverge; too small is slow.

## Computing gradients numerically

When you only have a function, not a formula, estimate each partial with a central difference: nudge one coordinate up and down by h and leave the rest alone.

```python
def partial(f, point, i, h=1e-5):
    up = list(point)
    down = list(point)
    up[i] += h
    down[i] -= h
    return (f(up) - f(down)) / (2 * h)

def gradient(f, point):
    return [partial(f, point, i) for i in range(len(point))]

f = lambda p: p[0] ** 2 * p[1] + 3 * p[1]
print(gradient(f, [1, 2]))   # about [4.0, 4.0]
```

Machine learning frameworks (PyTorch, JAX, TensorFlow) instead use **automatic differentiation**, which applies the chain rule exactly through the code. Numerical gradients are still the standard way to *test* that an analytic or automatic gradient is correct ("gradient checking").

## Common mistakes

- **Changing two variables at once** when estimating a partial. Only the i-th coordinate moves.
- **Forgetting that other variables are constants.** ∂/∂x of x²y is 2xy, not 2x.
- **Stepping up the gradient when minimising.** Descent subtracts the gradient.
- **Mutating the input point** in code. Copy it before nudging, or every later partial is wrong.

## Quick recap

- A partial derivative holds every other variable fixed.
- The gradient ∇f stacks the partials; it points uphill, and its length is the steepness.
- At a minimum or maximum, ∇f = 0.
- Gradient descent steps against the gradient; it is how most models are trained.

## Try it

In the **Do** tab, write the numerical partial derivative and use it to compute gradients and run a few steps of gradient descent.


Take derivatives of functions of several variables, build the gradient, and use it for gradient descent.

Partial derivatives and gradients

Partial derivatives, introduction

Gradient

## Why this matters

Most vectors change direction when a matrix acts on them. A few special vectors do not: the matrix only stretches or shrinks them. Those are **eigenvectors**, and the stretch factors are **eigenvalues**. They describe the natural behaviour of a system: the vibration modes of a bridge, the principal stresses in a material, the stable states of a Markov chain, the energy levels in quantum mechanics, and the principal components in a dataset (PCA). Google's original PageRank was an eigenvector.

## The core idea

A non-zero vector v is an eigenvector of a square matrix A, with eigenvalue λ (lambda), when:

```text
A v = λ v
```

Applying A to v gives a multiple of v: same line, scaled by λ. If λ = 2 the vector doubles, if λ = 0.5 it halves, if λ is negative it flips, and if λ = 1 it is left unchanged.

## A worked example

```text
A = | 2  1 |
    | 1  2 |
```

Try v = (1, 1): A v = (3, 3) = 3·(1, 1). So (1, 1) is an eigenvector with eigenvalue 3.

Try v = (1, −1): A v = (1, −1) = 1·(1, −1). An eigenvector with eigenvalue 1.

Try v = (1, 0): A v = (2, 1), which is not a multiple of (1, 0). Not an eigenvector.

Geometrically, A stretches the plane by 3 along the diagonal line y = x and leaves the line y = −x alone.

## Finding eigenvalues by hand

A v = λ v means (A − λI) v = 0 has a non-zero solution, which happens only when the determinant is zero:

```text
det(A − λI) = 0
```

For the example: (2 − λ)² − 1 = 0, so 2 − λ = ±1, giving λ = 3 or λ = 1. Two useful checks: the eigenvalues add up to the **trace** (sum of the diagonal, 2 + 2 = 4) and multiply to the **determinant** (4 − 1 = 3).

## Finding eigenvalues numerically

For anything larger than 3 × 3, nobody solves the characteristic equation by hand. Two numerical approaches:

**Power iteration** finds the eigenvalue with the largest size. Start with almost any vector, multiply by A, rescale to length 1, and repeat. Each multiplication amplifies the dominant eigenvector more than the others, so the vector lines up with it. The **Rayleigh quotient** (v · A v) / (v · v) then reads off its eigenvalue.

```python
import numpy as np

A = np.array([[2.0, 1.0], [1.0, 2.0]])
v = np.array([1.0, 0.0])
for _ in range(50):
    v = A @ v
    v = v / np.linalg.norm(v)
print(v @ A @ v)   # about 3.0
```

**Library routines** (`np.linalg.eig`, and `np.linalg.eigh` for symmetric matrices) compute all eigenvalues and eigenvectors at once with the QR algorithm:

```python
values, vectors = np.linalg.eigh(A)
print(values)      # [1. 3.]
```

Each column of `vectors` is an eigenvector, scaled to length 1.

## Where eigenvalues show up

- **Stability:** in a system x_(n+1) = A x_n, if every eigenvalue has size below 1, the system settles down; if any is above 1, it grows. Control engineers and population ecologists check exactly this.
- **Vibrations:** a structure's natural frequencies come from the eigenvalues of its stiffness and mass matrices. Resonance happens when forcing matches one.
- **Principal stresses:** the eigenvectors of the stress tensor are the directions with no shear.
- **PCA:** the eigenvectors of a covariance matrix are the directions of greatest variance in the data.

## Common mistakes

- **Allowing v = 0.** The zero vector satisfies A·0 = λ·0 for every λ, so it never counts.
- **Expecting unique eigenvectors.** Any non-zero multiple of an eigenvector is also one. Software returns unit-length versions, and the sign may flip.
- **Using `eig` for symmetric matrices.** `eigh` is faster, more accurate and returns real, sorted values.
- **Forgetting complex eigenvalues.** A rotation matrix moves every real vector off its line, so its eigenvalues are complex.

## Quick recap

- A v = λ v: an eigenvector keeps its direction; λ is its stretch factor.
- Eigenvalues solve det(A − λI) = 0; they sum to the trace and multiply to the determinant.
- Power iteration finds the dominant eigenvalue; `np.linalg.eigh` finds them all.
- They describe stability, vibration modes, principal stresses and PCA.

## Try it

In the **Do** tab, finish power iteration and compare its answer with NumPy's.


Find the vectors a matrix only stretches, compute eigenvalues by hand and by power iteration, and see where they matter.

Eigenvalues

Eigenvectors and eigenvalues | Chapter 14, Essence of linear algebra

## Why this matters

You now have the core of calculus and linear algebra: rates, accumulation, differential equations, vectors, matrices, gradients and eigenvalues. This lesson maps those ideas onto the courses and careers that use them, so you can recognise them when they appear in a statistics paper, a structural model or a machine learning library. There is no new technique here, only connections.

## Probability and statistics

- **Densities are integrated.** For a continuous random variable with density f(x), the probability of landing between a and b is the area ∫ from a to b of f(x) dx. The total area is 1. The mean is ∫ x·f(x) dx, a weighted average computed with an integral.
- **The CDF and PDF are a derivative pair.** The cumulative distribution function F(x) is the integral of the density, and the density is the derivative of F.
- **Maximum likelihood is optimisation.** To fit a model, you choose parameters that maximise the likelihood (usually its logarithm). That means setting derivatives to zero or running gradient ascent, exactly as in Lessons 3 and 10.
- **Covariance matrices** summarise how variables move together. Their eigenvectors give PCA's principal components, and their eigenvalues say how much variance each one explains.
- **Actuarial models** use exponential decay and integrals for survival probabilities and the present value of continuous cash flows, for example ∫ e^(−δt) · cash flow dt.

## Machine learning

- **Data are matrices.** A dataset with n rows and p features is an n × p matrix X. A linear model's predictions for every row at once are X @ w: one matrix-vector product.
- **Training is gradient descent.** The loss (for example mean squared error) is a function of the weights. Its gradient says how to nudge every weight to reduce the loss, and the update w ← w − η∇L repeats until it settles.
- **Backpropagation is the chain rule.** A neural network is a long composition of functions; backprop applies the chain rule layer by layer to get the gradient efficiently.
- **Least squares has a matrix solution.** Linear regression's best weights solve the normal equations (XᵀX) w = Xᵀy, a linear system like those in Lesson 8. Libraries solve it with stable factorisations rather than an inverse.
- **Embeddings and similarity** use dot products and cosine similarity (Lesson 7) to compare words, images and users.

## Engineering and the physical sciences

- **Statics** checks that force vectors sum to zero at every joint, which becomes a matrix system for a whole truss. Finite element software solves systems with millions of unknowns.
- **Dynamics and circuits** are differential equations: Newton's second law, RC and RLC circuits, heat flow. Simulation tools step them forward with methods descended from Euler's.
- **Vibration analysis** finds a structure's natural frequencies as eigenvalues.
- **Thermodynamics** relies on partial derivatives (how energy changes with temperature at constant pressure).
- **Chemistry** uses rate laws (differential equations) for kinetics, integrals for work and heat, and matrices for balancing reactions and fitting spectra.
- **Operations research** optimises costs subject to constraints (linear programming is linear algebra with inequalities), and sensitivity analysis is a set of partial derivatives.

## The tools you will use

Most of this work happens in software, not on paper. NumPy handles vectors, matrices and eigenvalues; SciPy adds integration (`scipy.integrate`), differential equation solvers (`solve_ivp`) and optimisers (`scipy.optimize`); SymPy does symbolic algebra; PyTorch and JAX compute gradients automatically. The `scientific-python`, `probability-and-inference` and `machine-learning-foundations` courses build on this one directly.

The value of knowing the maths is judgement: choosing a sensible step size, spotting a singular matrix, recognising that a model diverged because the learning rate was too large, or that an integral needs more points.

## Quick recap

- Probability: integrals give probabilities and expectations; fitting a model is optimisation.
- Machine learning: data are matrices, training is gradient descent, backprop is the chain rule.
- Engineering and science: systems of equations, differential equations and eigenvalues model structures, circuits and reactions.
- Software does the arithmetic; understanding the maths lets you trust and debug it.

## Try it

This lesson has no code exercise. Use the quiz in the **Do** tab to connect each technique with where it is used.


Connect derivatives, integrals, gradients and matrices to probability, machine learning, engineering and science.

Where this shows up

Gradient descent, how neural networks learn | Deep Learning Chapter 2

Partial derivatives, gradients and eigenvalues, and how all of this course powers probability, machine learning and engineering.