Linear Regression from Scratch: The Mathematics Explained

Linear regression is usually the first algorithm anyone meets in machine learning, and for good reason: it is simple enough to solve by hand, yet it contains almost every idea that later models use. In this article we build it from scratch, deriving the least-squares formulas, solving them with matrices, and then training the same model with gradient descent.

The Model

Suppose we have \(n\) data points \((x_1, y_1), \dots, (x_n, y_n)\). Simple linear regression assumes that \(y\) is roughly a straight-line function of \(x\):

$$ \hat{y} = \beta_0 + \beta_1 x $$

Here \(\beta_0\) is the intercept, \(\beta_1\) is the slope and \(\hat{y}\) ("y-hat") is the model's prediction. The gap between what we observe and what we predict is the residual:

$$ e_i = y_i - \hat{y}_i = y_i - (\beta_0 + \beta_1 x_i) $$

A good line makes the residuals small overall. Because positive and negative residuals would cancel if we simply added them, we square them first. This gives the sum of squared errors:

$$ \text{SSE}(\beta_0, \beta_1) = \sum_{i=1}^{n} \left(y_i - \beta_0 - \beta_1 x_i\right)^2 $$

Choosing \(\beta_0\) and \(\beta_1\) to minimise SSE is the method of ordinary least squares (OLS).

Deriving the Least-Squares Formulas

SSE is a smooth bowl-shaped function of two variables, so its minimum occurs where both partial derivatives are zero. Differentiating with respect to \(\beta_0\):

$$ \frac{\partial \,\text{SSE}}{\partial \beta_0} = -2\sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i) = 0 $$

Dividing by \(-2n\) gives \(\bar{y} - \beta_0 - \beta_1\bar{x} = 0\), so

$$ \beta_0 = \bar{y} - \beta_1 \bar{x} $$

This tells us something elegant: the best-fit line always passes through the point of means \((\bar{x}, \bar{y})\). Differentiating with respect to \(\beta_1\) and substituting the result above eventually gives

$$ \beta_1 = \frac{\sum (x_i - \bar{x})(y_i - \bar{y})}{\sum (x_i - \bar{x})^2} = \frac{S_{xy}}{S_{xx}} $$

The numerator measures how \(x\) and \(y\) vary together; the denominator measures how much \(x\) varies on its own.

Worked Example 1: Fitting a Line by Hand

A tutor records the hours five students spent on practice questions (\(x\)) and their quiz scores out of 10 (\(y\)):

• \(x\): 1, 2, 3, 4, 5
• \(y\): 2, 4, 5, 4, 5

Step 1: Means. \(\bar{x} = 15/5 = 3\) and \(\bar{y} = 20/5 = 4\).

Step 2: Deviations. \(x_i - \bar{x}\) = −2, −1, 0, 1, 2 and \(y_i - \bar{y}\) = −2, 0, 1, 0, 1.

Step 3: Products and squares.

$$ S_{xy} = (-2)(-2) + (-1)(0) + (0)(1) + (1)(0) + (2)(1) = 4 + 0 + 0 + 0 + 2 = 6 $$

$$ S_{xx} = 4 + 1 + 0 + 1 + 4 = 10 $$

Step 4: Coefficients.

$$ \beta_1 = \frac{6}{10} = 0.6, \qquad \beta_0 = 4 - 0.6(3) = 2.2 $$

So the fitted line is \(\hat{y} = 2.2 + 0.6x\). Each extra hour of practice is associated with 0.6 extra marks.

Step 5: How good is the fit? The predictions are 2.8, 3.4, 4.0, 4.6, 5.2, giving residuals −0.8, 0.6, 1.0, −0.6, −0.2. Then

$$ \text{SSE} = 0.64 + 0.36 + 1.00 + 0.36 + 0.04 = 2.4 $$

The total variation in \(y\) is \(\text{SST} = \sum (y_i - \bar{y})^2 = 4 + 0 + 1 + 0 + 1 = 6\). The coefficient of determination is

$$ R^2 = 1 - \frac{\text{SSE}}{\text{SST}} = 1 - \frac{2.4}{6} = 0.6 $$

The line explains 60% of the variation in scores. Notice also that the residuals sum to zero, exactly as the first derivative condition demands.

The Matrix Form: Normal Equations

With many input features, writing separate formulas becomes messy, so we switch to matrices. Stack the inputs into a design matrix \(X\) whose first column is all ones (for the intercept), and collect the coefficients into a vector \(\boldsymbol{\beta}\). The model is \(\hat{\mathbf{y}} = X\boldsymbol{\beta}\), and the loss is

$$ \text{SSE} = (\mathbf{y} - X\boldsymbol{\beta})^{T}(\mathbf{y} - X\boldsymbol{\beta}) $$

Setting the gradient to zero gives the normal equations:

$$ X^{T}X\,\boldsymbol{\beta} = X^{T}\mathbf{y} \quad\Longrightarrow\quad \boldsymbol{\beta} = (X^{T}X)^{-1}X^{T}\mathbf{y} $$

For our data, the rows of \(X\) are \((1, x_i)\), so

$$ X^{T}X = \begin{pmatrix} 5 & 15 \\ 15 & 55 \end{pmatrix}, \qquad X^{T}\mathbf{y} = \begin{pmatrix} 20 \\ 66 \end{pmatrix} $$

The determinant of \(X^TX\) is \(5(55) - 15(15) = 50\), so

$$ \boldsymbol{\beta} = \frac{1}{50}\begin{pmatrix} 55 & -15 \\ -15 & 5 \end{pmatrix}\begin{pmatrix} 20 \\ 66 \end{pmatrix} = \frac{1}{50}\begin{pmatrix} 110 \\ 30 \end{pmatrix} = \begin{pmatrix} 2.2 \\ 0.6 \end{pmatrix} $$

The same answer, as it must be. Geometrically, \(X\boldsymbol{\beta}\) is the orthogonal projection of \(\mathbf{y}\) onto the column space of \(X\), and the residual vector is perpendicular to every column.

Worked Example 2: One Step of Gradient Descent

Inverting \(X^TX\) is expensive when there are millions of features, so machine learning libraries often minimise the loss iteratively. We use the mean squared error \(J = \frac{1}{n}\sum (y_i - wx_i - b)^2\), whose gradients are

$$ \frac{\partial J}{\partial w} = -\frac{2}{n}\sum x_i\,(y_i - \hat{y}_i), \qquad \frac{\partial J}{\partial b} = -\frac{2}{n}\sum (y_i - \hat{y}_i) $$

Start from \(w = 0\), \(b = 0\) with learning rate \(\alpha = 0.05\).

Step 1: Predictions. With both parameters zero, every \(\hat{y}_i = 0\), so the residuals equal the \(y\) values: 2, 4, 5, 4, 5.

Step 2: Gradient for \(w\). \(\sum x_i y_i = 2 + 8 + 15 + 16 + 25 = 66\), so

$$ \frac{\partial J}{\partial w} = -\frac{2}{5}(66) = -26.4 $$

Step 3: Gradient for \(b\). \(\sum y_i = 20\), so \(\partial J / \partial b = -\frac{2}{5}(20) = -8\).

Step 4: Update.

$$ w \leftarrow 0 - 0.05(-26.4) = 1.32, \qquad b \leftarrow 0 - 0.05(-8) = 0.4 $$

The slope has overshot the true value of 0.6 on this first step, which is normal; repeated steps oscillate and settle towards \(w = 0.6\), \(b = 2.2\). If \(\alpha\) were too large the steps would grow and diverge instead, which is why choosing the learning rate matters.

From One Feature to Many

Real models rarely use a single input. A property price model might use area, number of rooms and distance from the city centre, giving \(\hat{y} = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_3\). Nothing in the matrix derivation changes: \(X\) simply gains one column per feature, and the normal equations still give the answer. Each coefficient is now interpreted as the expected change in \(y\) for a one-unit increase in that feature, holding the other features fixed. When two features are strongly correlated (multicollinearity), \(X^TX\) becomes nearly singular and the coefficients become unstable, which is one motivation for regularisation.

Assumptions Worth Knowing

The formulas above always produce a line, but interpreting it statistically requires assumptions:

1. Linearity: the true relationship is approximately linear.
2. Independence: residuals are not correlated with one another.
3. Constant variance (homoscedasticity): the spread of residuals does not grow with \(x\).
4. Normality of residuals: needed for confidence intervals and t-tests on coefficients.

Plotting residuals against fitted values is the quickest way to check the first and third assumptions.

Did You Know?

• The method of least squares was published by Adrien-Marie Legendre in 1805, while Carl Friedrich Gauss claimed to have used it years earlier to predict the orbit of the asteroid Ceres.
• The word "regression" comes from Francis Galton, who noticed in the 1880s that children of very tall parents tended to be tall but closer to the average height, a "regression towards mediocrity".
• Ridge regression simply adds \(\lambda \lVert\boldsymbol{\beta}\rVert^2\) to the loss, which changes the solution to \((X^TX + \lambda I)^{-1}X^T\mathbf{y}\) and keeps the matrix invertible.

Key Takeaways

• Linear regression chooses coefficients that minimise the sum of squared residuals.
• The closed-form solution is \(\beta_1 = S_{xy}/S_{xx}\) and \(\beta_0 = \bar{y} - \beta_1\bar{x}\), or \((X^TX)^{-1}X^T\mathbf{y}\) in matrix form.
• Gradient descent reaches the same answer iteratively and scales to very large datasets.
• \(R^2\) tells you the fraction of variation in \(y\) explained by the model.

If you are stuck on a regression derivation, a machine learning assignment or a thesis analysis, Amathify's subject experts and PhD scholars can provide clear step-by-step solutions. Message us on WhatsApp or use the Ask a Question form to get started.