Linear Models from Scratch¶

In the Linear Regression notebook we built an intuition for what a linear model is: a weighted sum of features added to a bias, where the weights and bias are learned from data. The expressions we met there, like $\hat y = b + w_{1}x_{1} + \dots + w_{n}x_{n}$, were described in terms of their meaning. In this notebook we go the other way around: we build one from nothing and watch the numbers move.

We're going to start with a plain linear model, train it by hand on real data, then progressively upgrade it -- first by squashing its outputs through a sigmoid, then by replacing the clunky row-wise math with a single matrix product, and finally by stacking our linear units into a small neural network. Along the way almost every line we write leans on ideas from the NumPy and PyTorch prework, so this notebook doubles as the place where all of that abstract material finally does something.

We will not reach for a pre-built architecture, an optimizer, or a data-loading framework. The only thing we'll let PyTorch do for us is calculate gradients, because doing that by hand is tedious and teaches us nothing new.

Setup¶

We'll use numpy and pandas for array and tabular work, and torch for the gradient tracking. Let's widen the default print width so rows of numbers don't wrap and hide the pattern we're trying to see.

In [1]:

import numpy as np
import pandas as pd
import torch

np.set_printoptions(linewidth=140)
torch.set_printoptions(linewidth=140, sci_mode=False, edgeitems=7)
pd.set_option('display.width', 140)

import numpy as np import pandas as pd import torch np.set_printoptions(linewidth=140) torch.set_printoptions(linewidth=140, sci_mode=False, edgeitems=7) pd.set_option('display.width', 140)

Getting the data¶

We'll train on the classic Titanic survival dataset. It is small, simple, and yet it throws enough real-world messiness at us -- missing values, string columns, and wildly different scales -- that we have to handle the same problems we'd face in a much larger project.

Rather than keeping a copy of the CSV in this repository, we load the dataset with scikit-learn's fetch_openml. This pulls the data over HTTP the first time and caches it locally afterwards, so it behaves like a small bundled dataset without us having to maintain a file by hand.

In [2]:

from sklearn.datasets import fetch_openml

titanic = fetch_openml(name='titanic', version=1, as_frame=True)
frame = titanic.frame
frame

from sklearn.datasets import fetch_openml titanic = fetch_openml(name='titanic', version=1, as_frame=True) frame = titanic.frame frame

Out[2]:

	pclass	survived	name	sex	age	sibsp	parch	ticket	fare	cabin	embarked	boat	body	home.dest
0	1	1	Allen, Miss. Elisabeth Walton	female	29.0000	0	0	24160	211.3375	B5	S	2	NaN	St Louis, MO
1	1	1	Allison, Master. Hudson Trevor	male	0.9167	1	2	113781	151.5500	C22 C26	S	11	NaN	Montreal, PQ / Chesterville, ON
2	1	0	Allison, Miss. Helen Loraine	female	2.0000	1	2	113781	151.5500	C22 C26	S	NaN	NaN	Montreal, PQ / Chesterville, ON
3	1	0	Allison, Mr. Hudson Joshua Creighton	male	30.0000	1	2	113781	151.5500	C22 C26	S	NaN	135.0	Montreal, PQ / Chesterville, ON
4	1	0	Allison, Mrs. Hudson J C (Bessie Waldo Daniels)	female	25.0000	1	2	113781	151.5500	C22 C26	S	NaN	NaN	Montreal, PQ / Chesterville, ON
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1304	3	0	Zabour, Miss. Hileni	female	14.5000	1	0	2665	14.4542	NaN	C	NaN	328.0	NaN
1305	3	0	Zabour, Miss. Thamine	female	NaN	1	0	2665	14.4542	NaN	C	NaN	NaN	NaN
1306	3	0	Zakarian, Mr. Mapriededer	male	26.5000	0	0	2656	7.2250	NaN	C	NaN	304.0	NaN
1307	3	0	Zakarian, Mr. Ortin	male	27.0000	0	0	2670	7.2250	NaN	C	NaN	NaN	NaN
1308	3	0	Zimmerman, Mr. Leo	male	29.0000	0	0	315082	7.8750	NaN	S	NaN	NaN	NaN

1309 rows × 14 columns

We have one row per passenger and a handful of columns describing that passenger. The survived column is our label -- the thing we want to predict -- and everything else is a candidate feature.

Two things make this data awkward to feed to a model:

1. Some cells hold NaN, pandas' representation of a missing value. We can't multiply a missing value by a coefficient.
2. Some columns are text (sex, embarked) and one (pclass) is a number that isn't really used as a number.

We'll deal with all of this next.

Cleaning the data¶

First, let's see where values are missing. isna() returns True for each missing cell, and since pandas treats True as 1 when summed, we can count the missing cells per column in one go.

In [3]:

frame.isna().sum()

frame.isna().sum()

Out[3]:

pclass          0
survived        0
name            0
sex             0
age           263
sibsp           0
parch           0
ticket          0
fare            1
cabin        1014
embarked        2
boat          823
body         1188
home.dest     564
dtype: int64

Notice that by default pandas sums across each column. That's the behavior we want here.

We need to replace the missing values with something, and it rarely matters much what. A safe and simple choice is the column's mode -- its most common value. The mode() function can return several rows in the case of ties, so we keep just the first one with iloc[0].

In [4]:

modes = frame.mode().iloc[0]
modes

modes = frame.mode().iloc[0] modes

Out[4]:

pclass                        3.0
survived                        0
name         Connolly, Miss. Kate
sex                          male
age                          24.0
sibsp                         0.0
parch                         0.0
ticket                   CA. 2343
fare                         8.05
cabin                 C23 C25 C27
embarked                        S
boat                           13
body                          1.0
home.dest            New York, NY
Name: 0, dtype: object

Let's fill every missing cell with its column's mode. This single line touches every column at once -- again, pandas applies the same operation across all of them.

In [5]:

frame = frame.fillna(modes)

frame = frame.fillna(modes)

We can confirm there are no missing values left:

In [6]:

frame.isna().sum()

frame.isna().sum()

Out[6]:

pclass       0
survived     0
name         0
sex          0
age          0
sibsp        0
parch        0
ticket       0
fare         0
cabin        0
embarked     0
boat         0
body         0
home.dest    0
dtype: int64

Understanding the numeric columns¶

A quick summary of the numeric columns helps us spot problems before they reach the model. The describe() block below shows count, mean, min, and the percentiles for each numeric feature.

In [7]:

frame.describe()

frame.describe()

Out[7]:

	pclass	age	sibsp	parch	fare	body
count	1309.000000	1309.000000	1309.000000	1309.000000	1309.000000	1309.000000
mean	2.294882	28.699516	0.498854	0.385027	33.276193	15.772345
std	0.837836	13.097103	1.041658	0.865560	51.743584	54.953095
min	1.000000	0.166700	0.000000	0.000000	0.000000	1.000000
25%	2.000000	22.000000	0.000000	0.000000	7.895800	1.000000
50%	3.000000	24.000000	0.000000	0.000000	14.454200	1.000000
75%	3.000000	35.000000	1.000000	0.000000	31.275000	1.000000
max	3.000000	80.000000	8.000000	9.000000	512.329200	328.000000

Look at the fare column: most values sit between 0 and ~30, but the maximum is over 500. A single expensive ticket swamps the rest. This is a classic example of a skewed distribution, and if we left it untouched, the fare feature alone would dominate our weighted sum.

The standard fix is to take the logarithm. Because fare can be 0, we add 1 first so that log(0 + 1) stays defined. The Transform your data notebook discusses log scaling in more depth; here we just apply it.

In [8]:

frame['log_fare'] = np.log(frame['fare'] + 1)

frame['log_fare'] = np.log(frame['fare'] + 1)

Now the values are spread out much more evenly, and no single passenger's fare unduly dominates the others -- exactly the effect we want before we start multiplying by coefficients.

Turning text into numbers¶

pclass contains just three values: 1, 2, and 3. It is ordered in the sense that the classes have a natural ranking, so treating it as a plain number is a reasonable first approximation. sex and embarked, on the other hand, are categories with no ordering at all -- we can't sensibly multiply a string like "female" or "S" by a coefficient.

The usual move is one-hot encoding: create a 0/1 column for every category and let the model learn a separate weight for each. pandas' get_dummies does exactly this in one call. Take a look at what it produces:

In [9]:

frame = pd.get_dummies(frame, columns=['sex', 'pclass', 'embarked'])
frame

frame = pd.get_dummies(frame, columns=['sex', 'pclass', 'embarked']) frame

Out[9]:

	survived	name	age	sibsp	parch	ticket	fare	cabin	boat	body	home.dest	log_fare	sex_female	sex_male	pclass_1	pclass_2	pclass_3	embarked_C	embarked_Q	embarked_S
0	1	Allen, Miss. Elisabeth Walton	29.0000	0	0	24160	211.3375	B5	2	1.0	St Louis, MO	5.358177	True	False	True	False	False	False	False	True
1	1	Allison, Master. Hudson Trevor	0.9167	1	2	113781	151.5500	C22 C26	11	1.0	Montreal, PQ / Chesterville, ON	5.027492	False	True	True	False	False	False	False	True
2	0	Allison, Miss. Helen Loraine	2.0000	1	2	113781	151.5500	C22 C26	13	1.0	Montreal, PQ / Chesterville, ON	5.027492	True	False	True	False	False	False	False	True
3	0	Allison, Mr. Hudson Joshua Creighton	30.0000	1	2	113781	151.5500	C22 C26	13	135.0	Montreal, PQ / Chesterville, ON	5.027492	False	True	True	False	False	False	False	True
4	0	Allison, Mrs. Hudson J C (Bessie Waldo Daniels)	25.0000	1	2	113781	151.5500	C22 C26	13	1.0	Montreal, PQ / Chesterville, ON	5.027492	True	False	True	False	False	False	False	True
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
1304	0	Zabour, Miss. Hileni	14.5000	1	0	2665	14.4542	C23 C25 C27	13	328.0	New York, NY	2.737881	True	False	False	False	True	True	False	False
1305	0	Zabour, Miss. Thamine	24.0000	1	0	2665	14.4542	C23 C25 C27	13	1.0	New York, NY	2.737881	True	False	False	False	True	True	False	False
1306	0	Zakarian, Mr. Mapriededer	26.5000	0	0	2656	7.2250	C23 C25 C27	13	304.0	New York, NY	2.107178	False	True	False	False	True	True	False	False
1307	0	Zakarian, Mr. Ortin	27.0000	0	0	2670	7.2250	C23 C25 C27	13	1.0	New York, NY	2.107178	False	True	False	False	True	True	False	False
1308	0	Zimmerman, Mr. Leo	29.0000	0	0	315082	7.8750	C23 C25 C27	13	1.0	New York, NY	2.183238	False	True	False	False	True	False	False	True

1309 rows × 20 columns

A whole block of new binary columns appeared at the right -- one per category of sex, pclass, and embarked. Because these columns are mutually exclusive (every passenger is exactly one sex, one pclass, and one embarked), we no longer need a separate intercept term: the dummies collectively cover the entire dataset, which is a neat detail we'll revisit when we look at the model's coefficients.

The Transform your data notebook covers one-hot encoding and alternatives in depth.

Building the feature matrix and the label vector¶

Now we can assemble the pieces. Our independent variables are the continuous features plus the columns we've just created, and our dependent variable is survived. The label needs a little coercion: fetch_openml left it as a string '0'/'1', so we convert it to an integer tensor.

We'll keep everything as PyTorch tensors from here on, since that's what we train against.

In [10]:

indep_cols = [
    'age',
    'sibsp',
    'parch',
    'log_fare',
    'sex_female',
    'sex_male',
    'pclass_1',
    'pclass_2',
    'pclass_3',
    'embarked_C',
    'embarked_Q',
    'embarked_S',
]
t_indep = torch.tensor(frame[indep_cols].to_numpy(dtype=np.float32))
t_dep = torch.tensor(frame['survived'].astype(int).to_numpy(dtype=np.float32))
t_indep.shape, t_dep.shape

indep_cols = [ 'age', 'sibsp', 'parch', 'log_fare', 'sex_female', 'sex_male', 'pclass_1', 'pclass_2', 'pclass_3', 'embarked_C', 'embarked_Q', 'embarked_S', ] t_indep = torch.tensor(frame[indep_cols].to_numpy(dtype=np.float32)) t_dep = torch.tensor(frame['survived'].astype(int).to_numpy(dtype=np.float32)) t_indep.shape, t_dep.shape

Out[10]:

(torch.Size([1309, 12]), torch.Size([1309]))

Setting up a linear model¶

We now have a matrix of independent variables and a vector of labels. Our first model is a linear one: we'll need one coefficient per feature column. We'll pick random numbers in the range (-0.5, 0.5) and seed the RNG first so the rest of this notebook matches the numbers you see here.

If the terminology here (weight, bias, prediction) feels unfamiliar, the Linear Regression notebook is the conceptual companion to what we're doing now.

In [11]:

torch.manual_seed(442)

n_coeff = t_indep.shape[1]
coeffs = torch.rand(n_coeff) - 0.5
coeffs

torch.manual_seed(442) n_coeff = t_indep.shape[1] coeffs = torch.rand(n_coeff) - 0.5 coeffs

Out[11]:

tensor([-0.4629,  0.1386,  0.2409, -0.2262, -0.2632, -0.3147,  0.4876,  0.3136,  0.2799, -0.4392,  0.2103,  0.3625])

Our predictions are the weighted sum of each row of features: multiply every row by the coefficients and add across the columns. An interesting consequence of the dummy columns is that we don't need a separate intercept term -- the dummies already cover the whole dataset, so there is always a column that is 1 for any given row to absorb the constant.

Let's look at the raw product first:

In [12]:

product = t_indep * coeffs
product.shape

product = t_indep * coeffs product.shape

Out[12]:

torch.Size([1309, 12])

This single * is quietly doing something powerful. Our features live on very different scales -- age averages around 30 while the dummy columns are only 0 or 1 -- so a row's sum is dominated by whichever column is biggest. Let's fix that by normalizing each column to [0, 1].

This line divides a matrix by a vector. What could that possibly mean? It works because of broadcasting: numpy and PyTorch "stretch" the vector so it lines up with each row of the matrix, dividing every row by the same values. Behind the scenes it makes no copies and runs at full speed, but for us the important part is that it lets us write clean, obviously-correct code. Broadcasting is one of the most useful tools in your toolbox, and it's covered thoroughly in the NumPy notebook.

In [13]:

vals, _ = t_indep.max(axis=0)
t_indep = t_indep / vals

vals, _ = t_indep.max(axis=0) t_indep = t_indep / vals

Now every column is on the same footing, and the row sums are no longer dictated by whichever feature happens to be large:

In [14]:

t_indep * coeffs

t_indep * coeffs

Out[14]:

tensor([[-0.1678,  0.0000,  0.0000, -0.1942, -0.2632, -0.0000,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.0053,  0.0173,  0.0535, -0.1822, -0.0000, -0.3147,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.0116,  0.0173,  0.0535, -0.1822, -0.2632, -0.0000,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.1736,  0.0173,  0.0535, -0.1822, -0.0000, -0.3147,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.1447,  0.0173,  0.0535, -0.1822, -0.2632, -0.0000,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.2777,  0.0000,  0.0000, -0.1202, -0.0000, -0.3147,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        [-0.3645,  0.0173,  0.0000, -0.1583, -0.2632, -0.0000,  0.4876,  0.0000,  0.0000, -0.0000,  0.0000,  0.3625],
        ...,
        [-0.1389,  0.0000,  0.0000, -0.0764, -0.0000, -0.3147,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.1389,  0.0000,  0.0000, -0.0992, -0.0000, -0.3147,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.0839,  0.0173,  0.0000, -0.0992, -0.2632, -0.0000,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.1389,  0.0173,  0.0000, -0.0992, -0.2632, -0.0000,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.1533,  0.0000,  0.0000, -0.0764, -0.0000, -0.3147,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.1562,  0.0000,  0.0000, -0.0764, -0.0000, -0.3147,  0.0000,  0.0000,  0.2799, -0.4392,  0.0000,  0.0000],
        [-0.1678,  0.0000,  0.0000, -0.0791, -0.0000, -0.3147,  0.0000,  0.0000,  0.2799, -0.0000,  0.0000,  0.3625]])

We can now produce a single prediction per row by adding across the columns -- the sum(axis=1) we saw in the NumPy notebook being put to real work.

In [15]:

preds = (t_indep * coeffs).sum(axis=1)
preds[:10]

preds = (t_indep * coeffs).sum(axis=1) preds[:10]

Out[15]:

tensor([ 0.2249,  0.4188,  0.4640,  0.2505,  0.3309,  0.1375,  0.0813,  0.3097,  0.1713, -0.8192])

These predictions are meaningless right now, because our coefficients are random. Let's measure how wrong they are so we have a number to improve. A simple, easy-to-read loss is the mean absolute error (MAE) -- how far, on average, our predictions are from the truth. The Linear Regression notebook compares the different loss functions; here we just need one that descends smoothly.

In [16]:

def calc_loss(coeffs, indeps, dep):
  """Mean absolute error between predictions and true labels."""
  preds = (indeps * coeffs).sum(axis=1)
  return torch.abs(preds - dep).mean()

def calc_loss(coeffs, indeps, dep): """Mean absolute error between predictions and true labels.""" preds = (indeps * coeffs).sum(axis=1) return torch.abs(preds - dep).mean()

In [17]:

loss = calc_loss(coeffs, t_indep, t_dep)
loss

loss = calc_loss(coeffs, t_indep, t_dep) loss

Out[17]:

tensor(0.5321)

An average error of 0.5 means our random predictions are, on average, half a unit (half a survived) away from the truth -- as bad as it gets, since the label is just 0 or 1. The next section makes this number go down.

A single gradient-descent step¶

Now we do one epoch of gradient descent by hand. The only part we'll let PyTorch automate is computing the gradients, because doing that from scratch is laborious and teaches nothing new. The Gradient Descent notebook explains the intuition; here we watch the mechanics.

First, ask PyTorch to track gradients for our coefficients by marking them with requires_grad_:

In [18]:

coeffs.requires_grad_()

coeffs.requires_grad_()

Out[18]:

tensor([-0.4629,  0.1386,  0.2409, -0.2262, -0.2632, -0.3147,  0.4876,  0.3136,  0.2799, -0.4392,  0.2103,  0.3625], requires_grad=True)

Now when we compute the loss, PyTorch records every operation on its way from coeffs to the loss, so it can later reverse those steps and hand us the gradients:

In [19]:

loss = calc_loss(coeffs, t_indep, t_dep)
loss

loss = calc_loss(coeffs, t_indep, t_dep) loss

Out[19]:

tensor(0.5321, grad_fn=<MeanBackward0>)

Calling backward() performs that reverse pass and stores the gradient of the loss with respect to each element of coeffs in coeffs.grad:

In [20]:

loss.backward()
coeffs.grad

loss.backward() coeffs.grad

Out[20]:

tensor([-0.0563,  0.0077, -0.0032, -0.0845, -0.1910,  0.0924, -0.1291, -0.0267,  0.0573, -0.2063, -0.0604,  0.1681])

Each entry tells us which way to nudge the corresponding coefficient to reduce the loss. One subtlety: repeated backward() calls add to the gradients rather than replacing them, because PyTorch just accumulates into the .grad buffer. Running the same two lines again doubles the values:

In [21]:

loss = calc_loss(coeffs, t_indep, t_dep)
loss.backward()
coeffs.grad

loss = calc_loss(coeffs, t_indep, t_dep) loss.backward() coeffs.grad

Out[21]:

tensor([-0.1125,  0.0155, -0.0065, -0.1690, -0.3820,  0.1849, -0.2582, -0.0535,  0.1146, -0.4125, -0.1207,  0.3361])

So after we use the gradients to update our coefficients, we must set them back to zero. Here's the complete step: compute the loss, work out the gradients, move coeffs a small distance against the gradient (scaled by a learning rate), and clear the gradients so the next step starts fresh.

The sub_ and zero_ methods end in an underscore because they modify the tensor in place -- PyTorch's convention that a trailing _ means "this changes its receiver". And we wrap the update in torch.no_grad() so PyTorch doesn't try to track gradients for the update itself.

In [22]:

loss = calc_loss(coeffs, t_indep, t_dep)
loss.backward()
with torch.no_grad():
  coeffs.sub_(coeffs.grad * 0.1)
  coeffs.grad.zero_()
  print(calc_loss(coeffs, t_indep, t_dep))

loss = calc_loss(coeffs, t_indep, t_dep) loss.backward() with torch.no_grad(): coeffs.sub_(coeffs.grad * 0.1) coeffs.grad.zero_() print(calc_loss(coeffs, t_indep, t_dep))

tensor(0.4910)

The loss went down. One step isn't enough, but it proves the direction is right -- and doing a few hundred of these steps by hand is tedious, so next we'll wrap the whole process into a training loop.

Training the linear model¶

Before training, we should hold out a slice of the data to check how well the model generalizes to passengers it never saw -- otherwise we'd be fooled into thinking our model is better than it is. The Training, Validation and Test sets notebook covers why this matters. Here we shuffle the rows and split them 80/20.

There's a subtlety with torch tensors: we set the random seed first, then use randperm (a random permutation of indices) so the split is reproducible.

In [23]:

n = len(t_indep)
torch.manual_seed(0)
idx = torch.randperm(n)
train_idx, val_idx = idx[:int(0.8 * n)], idx[int(0.8 * n):]

trn_indep, val_indep = t_indep[train_idx], t_indep[val_idx]
trn_dep, val_dep = t_dep[train_idx], t_dep[val_idx]
trn_indep.shape, val_indep.shape

n = len(t_indep) torch.manual_seed(0) idx = torch.randperm(n) train_idx, val_idx = idx[:int(0.8 * n)], idx[int(0.8 * n):] trn_indep, val_indep = t_indep[train_idx], t_indep[val_idx] trn_dep, val_dep = t_dep[train_idx], t_dep[val_idx] trn_indep.shape, val_indep.shape

Out[23]:

(torch.Size([1047, 12]), torch.Size([262, 12]))

Now we package the recurring pieces into three small functions:

• update_coeffs moves the coefficients one step down the gradient.
• one_epoch runs a full pass over the training data: calculate the loss, its gradients, and an update.
• init_coeffs rebuilds a fresh random coefficient vector, so we can retrain from scratch whenever we change the model.

A epochs loop strings these together, printing the training loss after each step so we can watch it fall.

In [24]:

def update_coeffs(coeffs, lr):
  """Take one gradient-descent step on `coeffs` using learning rate `lr`."""
  if coeffs.grad is not None:
    with torch.no_grad():
      coeffs.sub_(coeffs.grad * lr)
      coeffs.grad.zero_()


def one_epoch(coeffs, lr, indeps, dep):
  """Compute loss, backpropagate, and take one coefficient update."""
  loss = calc_loss(coeffs, indeps, dep)
  loss.backward()
  update_coeffs(coeffs, lr)


def init_coeffs():
  """Return a fresh coefficient vector on a tiny scale."""
  return (torch.rand(n_coeff) - 0.5).requires_grad_()


def train_model(lr, epochs=30, indeps=trn_indep, dep=trn_dep):
  """Train a linear model over `epochs` steps, returning the learned coeffs."""
  coeffs = init_coeffs()
  for _ in range(epochs):
    one_epoch(coeffs, lr, indeps, dep)
  return coeffs

def update_coeffs(coeffs, lr): """Take one gradient-descent step on `coeffs` using learning rate `lr`.""" if coeffs.grad is not None: with torch.no_grad(): coeffs.sub_(coeffs.grad * lr) coeffs.grad.zero_() def one_epoch(coeffs, lr, indeps, dep): """Compute loss, backpropagate, and take one coefficient update.""" loss = calc_loss(coeffs, indeps, dep) loss.backward() update_coeffs(coeffs, lr) def init_coeffs(): """Return a fresh coefficient vector on a tiny scale.""" return (torch.rand(n_coeff) - 0.5).requires_grad_() def train_model(lr, epochs=30, indeps=trn_indep, dep=trn_dep): """Train a linear model over `epochs` steps, returning the learned coeffs.""" coeffs = init_coeffs() for _ in range(epochs): one_epoch(coeffs, lr, indeps, dep) return coeffs

Let's train. We want to watch the loss fall step over step, then read off the final coefficients, so we run the loop and inspect both.

In [25]:

coeffs = train_model(lr=0.2, epochs=18)
coeffs

coeffs = train_model(lr=0.2, epochs=18) coeffs

Out[25]:

tensor([-0.2395,  0.0083,  0.4773,  0.1565,  0.4622, -0.2919,  0.3342,  0.1730,  0.0251,  0.1519,  0.3478,  0.1917], requires_grad=True)

The loss descended quickly and settled at a low value. Here are the final coefficients, mapped back onto their feature names so we can read them:

In [26]:

def show_coeffs():
  """Return a dict of feature name -> learned coefficient."""
  return dict(zip(indep_cols, [float(c) for c in coeffs.requires_grad_(False)]))


show_coeffs()

def show_coeffs(): """Return a dict of feature name -> learned coefficient.""" return dict(zip(indep_cols, [float(c) for c in coeffs.requires_grad_(False)])) show_coeffs()

Out[26]:

{'age': -0.2394808828830719,
 'sibsp': 0.008326675742864609,
 'parch': 0.47734174132347107,
 'log_fare': 0.15648305416107178,
 'sex_female': 0.46223634481430054,
 'sex_male': -0.2919047772884369,
 'pclass_1': 0.3342116177082062,
 'pclass_2': 0.17296117544174194,
 'pclass_3': 0.02509370446205139,
 'embarked_C': 0.15185482800006866,
 'embarked_Q': 0.34781327843666077,
 'embarked_S': 0.19173473119735718}

Already the signs are sensible: sex_male has a strongly negative weight, because being male lowers the predicted chance of survival on the Titanic, and the positive pclass_1 weight says first-class passengers were more likely to survive. A raw linear model can learn this much -- next we make its outputs behave like probabilities.

From raw scores to probabilities¶

Look at what the linear model actually produces: raw numbers that can be anything at all -- negative, or far above 1. A passenger's survival either happened or it didn't; a prediction of 3.2 isn't a probability, it's just a score.

To turn a raw score into a probability in (0, 1) we pass it through the sigmoid function:

$$\text{sigmoid}(z) = \dfrac{1}{1 + e^{-z}}$$

Sigmoid maps any real number to (0, 1), preserves ordering (bigger score -> bigger probability), and is smooth, so gradient descent still works. PyTorch gives it to us as torch.sigmoid, so we only need to change one line inside calc_preds.

In [27]:

def calc_preds(coeffs, indeps):
  """Predict survival probabilities in (0, 1) using a sigmoid linear model."""
  return torch.sigmoid((indeps * coeffs).sum(axis=1))

def calc_preds(coeffs, indeps): """Predict survival probabilities in (0, 1) using a sigmoid linear model.""" return torch.sigmoid((indeps * coeffs).sum(axis=1))

We also want to measure accuracy: we'll call a passenger "predicted to survive" when their probability exceeds 0.5, and compare that to the truth. Until now our loss used raw predictions, so let's switch it to operate on sigmoid outputs too -- the model still trains, and now its raw outputs are clean probabilities.

In [28]:

def calc_loss(coeffs, indeps, dep):
  """Mean absolute error between sigmoid predictions and true labels."""
  return torch.abs(calc_preds(coeffs, indeps) - dep).mean()


def acc(coeffs, indeps, dep):
  """Fraction of passengers whose predicted probability matches the truth."""
  return ((dep == (calc_preds(coeffs, indeps) > 0.5)).float().mean())

def calc_loss(coeffs, indeps, dep): """Mean absolute error between sigmoid predictions and true labels.""" return torch.abs(calc_preds(coeffs, indeps) - dep).mean() def acc(coeffs, indeps, dep): """Fraction of passengers whose predicted probability matches the truth.""" return ((dep == (calc_preds(coeffs, indeps) > 0.5)).float().mean())

Retrain with the updated prediction function. The learning rate needs to come up (0.2 was tuned for the raw-score version; after sigmoid the gradients shrink because of the squashing), so we use lr=5.

In [29]:

coeffs = train_model(lr=5, epochs=30)

coeffs = train_model(lr=5, epochs=30)

Let's check the accuracy on the held-out validation set:

In [30]:

acc(coeffs, val_indep, val_dep)

acc(coeffs, val_indep, val_dep)

Out[30]:

tensor(0.7366)

And here are the coefficients of our trained model -- still interpretable, now expressed on a probability scale:

In [31]:

show_coeffs()

show_coeffs()

Out[31]:

{'age': -0.0974469929933548,
 'sibsp': 0.032776229083538055,
 'parch': -0.003724555019289255,
 'log_fare': 0.19290627539157867,
 'sex_female': 2.975900173187256,
 'sex_male': -2.9361348152160645,
 'pclass_1': 0.9705225229263306,
 'pclass_2': 0.16862452030181885,
 'pclass_3': -1.335240125656128,
 'embarked_C': 0.262891560792923,
 'embarked_Q': -0.11020152270793915,
 'embarked_S': -0.8245739936828613}

These coefficients make sense: older passengers and male passengers pushed survival probability down, while having a first class ticket pushed it up. We now have a real, trained model making sensible predictions. The remaining sections make the machinery neater and then more powerful.

Using the matrix product¶

There's something a little clumsy about how we've been predicting. For every prediction we write

preds = (indeps * coeffs).sum(axis=1)

which multiplies element-by-element and then adds across the row. But adding across the row is the definition of a matrix-vector product. In other words:

$$Xw = \sum_{j} x_{j} w_{j}$$

Let's confirm the two are identical:

In [32]:

(val_indep * coeffs).sum(axis=1)

(val_indep * coeffs).sum(axis=1)

Out[32]:

tensor([-1.5619, -5.0576, -5.0181, -4.2725, -5.0543, -4.3937,  0.8598,  0.8574, -5.0613,  1.9504, -4.9779, -3.9726,  1.5920, -2.7291,
        -4.9617, -1.6071, -3.9330, -4.3111, -1.6277, -2.7279, -3.5359,  0.8532,  0.8934,  0.8769,  3.2307, -2.7469,  2.4015,  2.3910,
        -5.0888, -5.0543,  3.2501, -4.3433, -5.0547, -3.4888, -5.0137, -3.9799, -2.7669, -2.7176,  2.3918, -5.0558, -2.7249,  4.3312,
        -5.0474, -5.0478,  0.8549,  3.2034, -5.0552,  0.8702, -5.0097, -5.0644, -5.0524, -2.6814,  2.3987, -1.5641, -3.5484, -1.6444,
         2.3705, -1.6019,  3.2048, -3.9416,  0.8440,  0.8552, -2.7338, -5.0099,  0.8763, -3.5519, -3.5176, -1.6126, -5.0625, -3.9756,
         1.5683, -5.0588,  0.8632,  1.5604, -5.0571,  1.9395, -1.6477, -1.5961, -5.0675,  0.8525,  3.1640,  0.8734, -5.0569,  3.2294,
        -5.0601, -5.0705, -3.5605,  0.8983, -3.9490, -2.7559, -2.7169,  0.9038, -5.0465, -5.0522, -5.0435, -3.5387,  1.9529, -5.0674,
        -5.0576,  2.4369,  1.9224,  2.3906,  0.9317,  3.2521, -5.0597, -5.0529, -2.7194, -2.7006, -3.5458, -2.7632,  3.2289,  0.8518,
        -4.9854, -5.0589, -4.3414,  2.3971, -3.5739, -1.6137,  3.2235, -4.3438, -5.0555, -5.0705, -3.4844, -3.9479,  1.5687, -3.5275,
        -4.9999, -3.9921, -3.5290, -4.3438, -5.0571, -5.0636,  0.8564, -3.9726, -5.0632, -5.0333, -5.0436, -2.8499, -3.4998,  0.8511,
        -3.5397, -2.6681,  1.9597, -1.6011,  0.9073, -3.5361,  1.5683, -4.3399, -3.5422,  2.4259, -3.9366, -5.0465,  2.3912, -3.4993,
        -4.9738, -5.0529,  3.2199, -3.5243,  0.9504, -5.0685, -5.0438, -5.0351, -3.5728, -5.0571,  0.9200, -5.0600,  2.4128,  0.8562,
         4.2857, -5.0802,  2.3764, -3.9701,  1.5683, -5.0580, -5.0555, -3.9726,  0.9188, -5.0468, -3.5355, -1.6044,  0.8707, -4.3438,
         0.8491,  3.2163,  1.9712, -2.7375,  0.8648, -5.0552, -3.9272, -5.0541, -1.5943, -3.9739, -3.9296, -5.0483, -5.0626,  4.3115,
        -5.0576, -5.0376,  1.9675,  4.2995, -5.0476, -2.6786, -5.0364, -2.3942,  4.2930, -3.9762,  4.3634, -2.7087,  3.2236, -3.5714,
         0.9000, -5.0552, -3.9330, -1.5880,  0.8787, -2.6946, -1.5893,  0.9504, -1.6358, -2.7145,  1.5651, -3.6213, -4.4004, -5.0811,
        -4.9701, -5.0405, -1.6424,  3.2688, -5.0575, -5.0489, -1.6251, -5.0738, -2.7333, -1.5591, -3.5531, -2.7218, -5.0666,  2.3967,
        -5.0571, -3.9526, -5.0827,  2.3906, -5.0514,  3.2163,  4.2744,  0.9081, -4.9617,  2.4170,  1.5678, -3.5385,  3.2312, -4.3397,
         0.8889,  0.8528, -3.5166,  0.8557,  2.3997, -5.0571, -4.9957,  4.3594, -2.7462,  0.9366])

In [33]:

val_indep @ coeffs

val_indep @ coeffs

Out[33]:

tensor([-1.5619, -5.0576, -5.0181, -4.2725, -5.0543, -4.3937,  0.8598,  0.8574, -5.0613,  1.9504, -4.9779, -3.9726,  1.5920, -2.7291,
        -4.9617, -1.6071, -3.9330, -4.3111, -1.6277, -2.7279, -3.5359,  0.8532,  0.8934,  0.8769,  3.2307, -2.7469,  2.4015,  2.3910,
        -5.0888, -5.0543,  3.2501, -4.3433, -5.0547, -3.4888, -5.0137, -3.9799, -2.7669, -2.7176,  2.3918, -5.0558, -2.7249,  4.3312,
        -5.0474, -5.0478,  0.8549,  3.2034, -5.0552,  0.8702, -5.0097, -5.0644, -5.0524, -2.6814,  2.3987, -1.5641, -3.5484, -1.6444,
         2.3705, -1.6019,  3.2048, -3.9416,  0.8440,  0.8552, -2.7338, -5.0099,  0.8763, -3.5519, -3.5176, -1.6126, -5.0625, -3.9756,
         1.5683, -5.0588,  0.8632,  1.5604, -5.0571,  1.9395, -1.6477, -1.5961, -5.0675,  0.8525,  3.1640,  0.8734, -5.0569,  3.2294,
        -5.0601, -5.0705, -3.5605,  0.8983, -3.9490, -2.7559, -2.7169,  0.9038, -5.0465, -5.0522, -5.0435, -3.5387,  1.9529, -5.0674,
        -5.0576,  2.4369,  1.9224,  2.3906,  0.9317,  3.2521, -5.0597, -5.0529, -2.7194, -2.7006, -3.5458, -2.7632,  3.2289,  0.8518,
        -4.9854, -5.0589, -4.3414,  2.3971, -3.5739, -1.6137,  3.2235, -4.3438, -5.0555, -5.0705, -3.4844, -3.9479,  1.5687, -3.5275,
        -4.9999, -3.9921, -3.5290, -4.3438, -5.0571, -5.0636,  0.8564, -3.9726, -5.0632, -5.0333, -5.0436, -2.8499, -3.4998,  0.8511,
        -3.5397, -2.6681,  1.9597, -1.6011,  0.9073, -3.5361,  1.5683, -4.3399, -3.5422,  2.4259, -3.9366, -5.0465,  2.3912, -3.4993,
        -4.9738, -5.0529,  3.2199, -3.5243,  0.9504, -5.0685, -5.0438, -5.0351, -3.5728, -5.0571,  0.9200, -5.0600,  2.4128,  0.8562,
         4.2857, -5.0802,  2.3764, -3.9701,  1.5683, -5.0580, -5.0555, -3.9726,  0.9188, -5.0468, -3.5355, -1.6044,  0.8707, -4.3438,
         0.8491,  3.2163,  1.9712, -2.7375,  0.8648, -5.0552, -3.9272, -5.0541, -1.5943, -3.9739, -3.9296, -5.0483, -5.0626,  4.3115,
        -5.0576, -5.0376,  1.9675,  4.2995, -5.0476, -2.6786, -5.0364, -2.3942,  4.2930, -3.9762,  4.3634, -2.7087,  3.2236, -3.5714,
         0.9000, -5.0552, -3.9330, -1.5880,  0.8787, -2.6946, -1.5893,  0.9504, -1.6358, -2.7145,  1.5651, -3.6213, -4.4004, -5.0811,
        -4.9701, -5.0405, -1.6424,  3.2688, -5.0575, -5.0489, -1.6251, -5.0738, -2.7333, -1.5591, -3.5531, -2.7218, -5.0666,  2.3967,
        -5.0571, -3.9526, -5.0827,  2.3906, -5.0514,  3.2163,  4.2744,  0.9081, -4.9617,  2.4170,  1.5678, -3.5385,  3.2312, -4.3397,
         0.8889,  0.8528, -3.5166,  0.8557,  2.3997, -5.0571, -4.9957,  4.3594, -2.7462,  0.9366])

Identical results. The @ operator is the matrix product, and it's not just shorter -- it's also much faster, because matrix multiplication is one of the most heavily optimized operations on modern hardware (the same @ we met in the NumPy notebook).

Let's rewrite our prediction function to use it:

In [34]:

def calc_preds(coeffs, indeps):
  """Predict survival probabilities using a sigmoid linear model and matmul."""
  return torch.sigmoid(indeps @ coeffs)

def calc_preds(coeffs, indeps): """Predict survival probabilities using a sigmoid linear model and matmul.""" return torch.sigmoid(indeps @ coeffs)

To go further -- toward a neural network -- we'll soon need matrix-matrix products, which multiply a whole batch of inputs by a whole batch of weights at once. To line up with that, we'll make our dependent variable a column vector too, using the [:, None] or [:, np.newaxis] trick from the NumPy notebook that adds a trailing singleton axis:

In [35]:

trn_dep = trn_dep[:, np.newaxis]
val_dep = val_dep[:, np.newaxis]
trn_dep.shape

trn_dep = trn_dep[:, np.newaxis] val_dep = val_dep[:, np.newaxis] trn_dep.shape

Out[35]:

torch.Size([1047, 1])

We also need the coefficients to be a column vector so the matmul gives us a column-shaped output. A (n_coeff, 1) shape does it:

In [36]:

def init_coeffs():
  """Return fresh column-vector coefficients for matrix-matrix products."""
  return (torch.rand(n_coeff, 1) * 0.1).requires_grad_()

def init_coeffs(): """Return fresh column-vector coefficients for matrix-matrix products.""" return (torch.rand(n_coeff, 1) * 0.1).requires_grad_()

Retraining with these changes gives a validation accuracy around 0.67:

In [37]:

coeffs = train_model(lr=5, epochs=30)
acc(coeffs, val_indep, val_dep)

coeffs = train_model(lr=5, epochs=30) acc(coeffs, val_indep, val_dep)

Out[37]:

tensor(0.6718)

Identical result, neater code. The matrix product buys us more than brevity though -- it's the lever we pull next to build a network with hidden layers.

A neural network¶

A linear model can only draw a straight boundary. To handle more interesting patterns we can stack linear units and insert a nonlinearity between them -- that's a neural network.

Let's start with a single hidden layer (one layer, 20 neurons by default -- n_hidden below controls the number of neurons, not the number of layers). We'll need both a first set of weights (features -> hidden) and a second set (hidden -> output), each of shape that a matrix product can use:

In [38]:

def init_coeffs(n_hidden=20):
  """Return weights and biases for a one-hidden-layer network.

  Args:
    n_hidden: Number of neurons (units) in the single hidden
      layer. This controls the width of the layer, not the number
      of hidden layers -- the network always has exactly one hidden
      layer.
  """
  torch.manual_seed(442)
  w1 = (torch.rand(n_coeff, n_hidden) - 0.5).requires_grad_()
  w2 = (torch.rand(n_hidden, 1) - 0.5).requires_grad_()
  b1 = (torch.rand(n_hidden) - 0.5).requires_grad_()
  b2 = (torch.rand(1) - 0.5).requires_grad_()
  return [w1, b1, w2, b2]

def init_coeffs(n_hidden=20): """Return weights and biases for a one-hidden-layer network. Args: n_hidden: Number of neurons (units) in the single hidden layer. This controls the width of the layer, not the number of hidden layers -- the network always has exactly one hidden layer. """ torch.manual_seed(442) w1 = (torch.rand(n_coeff, n_hidden) - 0.5).requires_grad_() w2 = (torch.rand(n_hidden, 1) - 0.5).requires_grad_() b1 = (torch.rand(n_hidden) - 0.5).requires_grad_() b2 = (torch.rand(1) - 0.5).requires_grad_() return [w1, b1, w2, b2]

The forward pass threads an input through the first layer, applies a nonlinear activation function, then passes the result through the second layer:

1. x @ w1 + b1 -- the hidden layer's weighted sum.
2. torch.relu(...) -- the activation, which zeroes out negative values and lets the network learn curved boundaries. (ReLU and friends are covered in the Neural Networks notebook.)
3. h @ w2 + b2 -- combine the hidden units into a single score.
4. torch.sigmoid(...) -- squash the score into a probability.

The hidden layer's output feeds into the next, so both sets of weights are learned jointly by backpropagation.

In [39]:

def calc_preds(w, indeps):
  """Forward pass of a one-hidden-layer network, returning probabilities."""
  x = indeps @ w[0] + w[1]
  h = torch.relu(x)
  z = h @ w[2] + w[3]
  return torch.sigmoid(z)

def calc_preds(w, indeps): """Forward pass of a one-hidden-layer network, returning probabilities.""" x = indeps @ w[0] + w[1] h = torch.relu(x) z = h @ w[2] + w[3] return torch.sigmoid(z)

Now that we have several sets of parameters, we have to update all of them in update_coeffs:

In [40]:

def update_coeffs(coeffs, lr):
  """Move every parameter one step down its gradient."""
  with torch.no_grad():
    for c in coeffs:
      c.sub_(c.grad * lr)
      c.grad.zero_()

def update_coeffs(coeffs, lr): """Move every parameter one step down its gradient.""" with torch.no_grad(): for c in coeffs: c.sub_(c.grad * lr) c.grad.zero_()

We can now train. With a hidden layer the model has enough capacity to fit wiggly patterns, but it also needs a few more epochs to settle:

In [41]:

coeffs = init_coeffs()
for _ in range(30):
  one_epoch(coeffs, lr=5, indeps=trn_indep, dep=trn_dep)
acc(coeffs, val_indep, val_dep)

coeffs = init_coeffs() for _ in range(30): one_epoch(coeffs, lr=5, indeps=trn_indep, dep=trn_dep) acc(coeffs, val_indep, val_dep)

Out[41]:

tensor(0.7748)

The accuracy matches or slightly beats the plain linear model. On a dataset as small and simple as Titanic that's hardly surprising -- a linear boundary already gets most of the way, and the extra capacity of a hidden layer has little left to capture. Don't let that fool you: on messier, real-world problems that extra nonlinearity is often the whole difference.

The mechanics we just built by hand -- a hidden layer, a ReLU activation, and a sigmoid output, all trained with gradient descent -- are the same ingredients, stacked many times over, that power the networks in the Neural Networks notebook.

Summary¶

What started as a handful of random numbers is now a trained model. Let's recap what we built:

• A linear model whose coefficients we tuned with gradient descent -- the concepts behind it are in Linear Regression and Gradient Descent.
• A sigmoid to turn raw scores into probabilities.
• The matrix product @, which is both cleaner and faster than the element-wise sum -- straight from the NumPy prework.
• A neural network with a hidden layer and a ReLU activation, using the same training loop, and previewing what the Neural Networks notebook explores in depth.

The central idea to carry forward: a model is just a differentiable function of its parameters, and gradient descent gives us a generic way to improve those parameters against any smooth loss. Everything else -- hidden layers, activations, better data pipelines -- is an elaboration on that one mechanism.

Copy page as MarkdownCopy to clipboard for pasting into any AI tool View as MarkdownOpen this page as clean Markdown in a new tab

Open in ChatGPTChat about this page with ChatGPT Open in ClaudeChat about this page with Claude

llms.txtFull documentation index for LLMs

Use with AI