Linear Regression from Scratch
I’m starting off with one of the simplest machine learning techniques: linear regression. The mathematical part of this post requires a good working knowledge of linear algebra and calculus to follow, which will be the case for the rest of the series. This is unavoidable since these subjects underpin a huge amount of machine learning and are prerequisites for a deep understanding. With that said, let’s dive in.
Simple Linear Regression
Linear regression is the task of taking a series of points and finding a line of best fit. Before we can figure out how to find a line of best fit, we need to understand what that actually means.
We can intuitively tell that the line in Figure 1 is a better fit to its data than say, this one:
Why? Because the points in Figure 2 are farther away from the line than those in Figure 1. Let’s figure out how to mathematically formalize this intuition so we can define precisely what “best fit” means.
To keep things simple (and make visualization easy) we’ll start off in two dimensions. In this case, data points are (x, y) pairs and can be shown on a graph like the figures above. We want to find a linear function f(x) = kx which best represents the data. Note that this model assumes a line that passes through the origin. We won’t consider the possibility of an intercept other than the origin yet.
Say we have a collection of n data points,
For each x-value, we can use our model to obtain a predicted y-value. This situation, with only one independent variable (x) and one dependent variable (y) is known as simple linear regression. We’ll use a prime symbol to differentiate the predicted y-values from the actual y-values. So, the model’s predicted y-value for a given x-value is given by the formula
Now let’s put all the x values into a vector, and the y values into another vector.
This is called vectorization, and has many benefits for problems in data science. Combining a lot of individual values into vectors makes mathematical formulas much more compact and easier to reason about. Vectorization in code also improves performance. It’s much faster to perform vector arithmetic on a large array of values than to go through using a for loop and operate on them one at a time. Many numerical computing libraries like Numpy are built for fast vector arithmetic. Vectorization also allows for parallelization with hardware like GPUs, where operations are done simultaneously to different elements of an array. Once again, this wouldn’t be possible using a loop where each operation happens one after the other.
We can also create a vector containing all of our predicted y-values:
In order to find the line of best fit, we need to know how far y’ is from the vector of actual values, y. We can look at the difference of the two vectors: y’ — y. However, this is itself a vector and we want a single number representing the error of the model. We are going to use sum-of-squares error (SSE). SSE is equal to ||y’ — y||², the squared magnitude of the difference vector. It’s called “sum-of-squares” because it’s equal to the sum of the squared entries of y’ — y:
As for why we use ||y’ — y||² instead of just ||y’ — y||, one answer is that squared magnitude is just much simpler to work with. ||y’ — y|| has an extra square root symbol outside the summation:
This would make working with the formula a LOT more tedious when taking derivatives, for instance.
Now that we’ve defined the error of a linear regression model, we need to find out how to minimize it. Let’s expand out the expression for SSE.
If we substitute in kx for y’ as per Equation 1, we get
We need to find the value of k that minimizes the error when x and y are held constant. To do so, we can set the derivative of Expression 1 w.r.t. k to zero and solve:
This gives us the value of k which minimizes the sum-of-squares error. Armed with this knowledge, we can code a SimpleLinearRegressor. It’ll have just one instance variable — the slope, k.
class SimpleLinearRegressor:
"""Performs simple linear regression."""
def __init__(self):
self.k = NoneIt will have a predict method that takes a number or array of numbers as a parameter and gives a prediction or array of predictions as a result.
def predict(self, x):
"""
Gives predicted y-values from an input x-value, or vector of x-values.
:param x: The input value(s).
:return: The predicted y-value(s).
"""
if self.k is None:
raise RegressionModelNotFitError('Oh no! The model has not been fit!')
return self.k * xIt will also need a fit method that finds the proper value of k based on the x and y vectors using Equation 2. This is the real guts of the class.
def fit(self, x, y):
"""
Fits the model based on the given vectors of x and y values.
:param x: A vector of x-values.
:param y: A vector of y-values.
:return: The sum-of-squares error of the fitted model.
"""
self.k = x @ y / (x @ x)
diff = self.predict(x) - y
return diff @ diffWe should generate some data to test the model. I’ll create a function that generates a vector of random x-values within a range, computes corresponding y-values using a linear model, and then adds Gaussian noise to those y-values.
def generate_noisy_data(n_points, slope, x_range, noise_stddev):
"""
Generates 2d data points based on a linear relationship with added Gaussian noise.
:param n_points: The number of data points to generate.
:param slope: The slope of the line.
:param x_range: The range from which to draw x-values.
:param noise_stddev: The standard deviation of the Gaussian noise to add to each y-value.
:return: Vectors of x and y values.
"""
x = np.random.uniform(*x_range, n_points)
y = slope * x + np.random.normal(scale=noise_stddev, size=n_points)
return x, yLet’s see how well the SimpleLinearRegressor can recover the original slope from the randomized data. I’ll use matplotlib for visualization.
x_range = np.array([0, 5])
x, y = generate_noisy_data(n_points=20, slope=0.42, x_range=x_range, noise_stddev=0.5)
plt.scatter(x, y)
regressor = SimpleLinearRegressor()
fit = regressor.fit(x, y)
slope = regressor.k
plt.plot(x_range, [0, 2 * x_range[1]], color='red')
plt.text(3, 0, f'Error: {"{:.2f}".format(fit)}\nPredicted slope: {"{:.2f}".format(slope)}')
plt.show()Here’s the output of one run:
As you can see, it does an excellent job! The slope predicted by the regression model is equal to the slope input of generate_noisy_data to at least 2 decimal places.
Multiple Linear Regression
We’ve seen how to perform linear regression with just one independent variable x, and one dependent variable, y. Now let’s suppose that y depends on m independent variables, so we’re working with (m + 1)-dimensional data. We might have n data points which look like this:
Here x_ij denotes the value of the jth independent variable in the ith data point.
Consolidating our data through vectorization is always a good first step.
We can collect all the y-values into a vector as before:
However, now that our x data has two indices, using a vector for the xs is no longer ideal. Instead we can collect them into a matrix where each row is a single data point:
I’ll use capital X_ij and lowercase x_ij from now on interchangeably to refer to entries in this matrix.
The linear model which we’re trying to fit to our data now looks a little more complicated:
We have m coefficients, or “slopes”, one for each independent variable, denoted by βs.
We can make a vector of each data point
The matrix X can be thought of as having each of these vectors as its rows:
Let’s also make a vector of all the β coefficients
Equation 3 can then be expressed very concisely as
However, we can get even more concise by combining the equations for each predicted value, y’_i, into a vector of all the predicted values.
Just like in simple linear regression, we’ll be trying to minimize the sum-of-squares error ||y’ — y||²
Using Equation 4 we can expand out the SSE formula in terms of X, y, and β.
Look familiar? It’s a lot like the formula for error in simple linear regression. We need to find the value of β which minimizes it. First note that the ||y||² term is irrelevant since it doesn’t depend on β, so we’re really trying to minimize
We could stop here since Numpy has a numpy.linalg.lstsq method that would find the proper value of β with just X and y as input. While technically not against my rule of only using Python and Numpy for computations, this feels too much like cheating in a post about “linear regression from scratch.” Instead, I’m going to dive into some math.
Get Perrin Alaine-Sedano’s stories in your inbox
Join Medium for free to get updates from this writer.
In order to minimize Expression 2, we should set its gradient w.r.t. β equal to zero and solve. To do this, we’ll convert Expression 2 into componentwise form and then take the derivative w.r.t each component of β individually.
Using the componentwise formula for the dot product,
and for matrix-vector multiplication,
we can convert Expression 2 to componentwise form:
Now let’s take the derivative of Expression 3 w.r.t. a particular component of β, β_l:
To simplify Expression 4, we’ll need to take the derivative of two summations:
and
Let’s tackle each of them individually.
Expression 5 Derivative
Expression 5 can be expanded out as:
In English, this means “the part where neither j nor k is equal to l, plus the part where k is equal to l but j is not, plus the part where j is equal to l but k is not, plus the part where j and k are both equal to l.” Since both j and k must either be equal to or not equal to l, these four sums cover all possibilities. All of them together thus equal the original summation of Expression 5.
Notice that the two middle sums in Expression 7 are identical except for the name of an index variable (j vs k). Since the name is arbitrary and we could just as easily rename the index variable in the third sum to j, those two sums have the same value. Thus the expression can be rewritten as
Now we can take the derivative:
Note that the first term becomes zero since it doesn’t contain a β_l.
Expression 6 Derivative
Finding the derivative of Expression 6 is much simpler:
Here once again we split the second summation into a part containing β_l and a part not containing β_l. The latter is zeroed out during differentiation.
Putting it All Together
Now let’s substitute the derivatives we just found, Expressions 8 and 9, back into Expression 4 and simplify. We can then convert back from componentwise form to vector form.
At this point we can use the identity
to transform Equation 5 further:
And we’re done! Any β for which the gradient of the error is zero must obey Equation 6.
You may be wondering whether setting the gradient to zero actually guarantees an optimal solution. After all, this could merely find a local minimum rather than the global minimum. Fortunately, linear regression is a convex optimization problem. This math stack exchange answer provides a proof. An important property of convex optimization problems is that any local minimum is also a global minimum, so there’s nothing to worry about.
Now that we’re assured of our solution’s correctness, we need to solve Equation 6 for β. Numpy provides a numpy.linalg.solve function, but it only works if the equation has just one solution. Another option would be converting the matrix to reduced row echelon form, but somewhat surprisingly Numpy doesn’t have a utility for this. After some digging I discovered numpy.linalg.qr, which performs QR factorization of an input matrix. An answer on math stack exchange and its comments were instrumental in helping me learn how to use QR factorization for equation solving.
If we want to solve the linear equation Ax = b where A is a square matrix (note that X^TX must be square), we can start by finding an orthonormal square matrix Q and an upper triangular matrix R such that QR = A. The equation Ax = b becomes QRx = b. Since Q must be invertible, the equation may be further transformed into Rx = Q^-1b. R is upper triangular and the right side is just a vector, so all I really need is the ability to solve Uv = w, where U is an upper triangular matrix.
I created a function, solve_upper_triangular, to perform the task. I won’t go into too much detail about its workings since this post isn’t on how to solve linear equations, but essentially it starts at the last row of the matrix and works backwards, at each row substituting in any previously set values for variables. It then assigns the value 1 to all but one of the remaining variables with nonzero coefficients in the row and solves for that final variable in terms of all the others.
def solve_upper_triangular(a, b):
"""
Solves the linear equation ax = b for x.
:param a: An n x n upper triangular matrix.
:param b: An n-dimensional vector.
:return: A vector x for which ax = b.
"""
tracker = np.zeros(a.shape[1])
result = np.zeros(a.shape[1])
for row, val in zip(a[::-1], b[::-1]):
unset_var_indices = np.where((tracker == 0) & (row != 0))[0]
if len(unset_var_indices) == 0:
if np.isclose(result @ row, val):
continue
else:
raise UnsolvableError('The given values of a and b result in an unsolvable equation.')
tracker[unset_var_indices] = 1
result[unset_var_indices[1:]] = 1
i = unset_var_indices[0]
result[i] = (val - result @ row) / row[i]
return resultI’m now ready to create a MultipleLinearRegressor.
class MultipleLinearRegressor:
"""Performs multiple linear regression."""
def __init__(self):
self.beta = NoneLike the SimpleLinearRegressor, it will have a predict method and a fit method.
The former simply calculates the matrix product between a matrix X or vector x and β.
def predict(self, x):
"""
Gives predicted y-values from a given array of x-values.
:param x: Vector or matrix of x-values.
:return: Vector of predicted y-values.
"""
if self.beta is None:
raise RegressionModelNotFitError('Oh no! The model has not been fit!')
return x @ self.betaThe fit method solves Equation 6 by QR factorizing X^TX and then using solve_upper_triangular to find the solution of Rβ = Q^-1X^Ty. It also returns the sum-of-squares error for the fitted model.
def fit(self, x, y):
"""
Fits the model based on a matrix of x-values and vector of corresponding y-values.
:param x: Matrix of x-values.
:param y: Vector of y-values.
:return: The sum-of-squares error of the fitted model.
"""
x_t = x.transpose()
q, r = np.linalg.qr(x_t @ x)
vec = np.linalg.inv(q) @ x_t @ y
self.beta = solve_upper_triangular(r, vec)
diff = self.predict(x) - y
return diff @ diffLet’s see how the MultipleLinearRegressor performs. I’ll create a generate_noisy_data function very similar to the previous one but which takes a vector β of parameters and generates a matrix X and vector y of data points, then adds Gaussian noise to the y-values as before.
def generate_noisy_data(n_data_points, n_independent_variables, beta, x_range, noise_stddev):
"""
Generates data points based on a linear relationship with added Gaussian noise.
:param n_data_points: The number of data points to generate.
:param n_independent_variables: The number of independent variables in each data point.
:param beta: The vector of coefficients for the independent variables.
:param x_range: The range from which to draw x-values.
:param noise_stddev: The standard deviation of the Gaussian noise to add to each y-value.
:return: Matrix of x-values and vector of y-values.
"""
x = np.random.uniform(*x_range, (n_data_points, n_independent_variables))
y = x @ beta + np.random.normal(scale=noise_stddev, size=n_data_points)
return x, yNow it’s time to generate some data and see how well the regressor recovers the original β.
regressor = MultipleLinearRegressor()
x, y = generate_noisy_data(n_data_points=500,
n_independent_variables=10,
beta=np.array([-10, 5, -8, -2, 1, -3, 4, -5, -1, 3]),
x_range=np.array([-100, 100]),
noise_stddev=50)
sse = regressor.fit(x, y)
print(f'Sum Squared Error: {sse}')
print(f'Beta: {regressor.beta}')The output from one run is
Sum Squared Error: 1259196.6403705715
Beta: [-9.95436533 5.02469925 -7.95431319 -1.97266714 1.03726794 -2.95935233
4.03854255 -4.98106051 -1.01840914 3.0410695]As we can see it does a pretty good job of getting close to the original parameter values.
What about Bias?
So far I’ve only discussed regression models with a y-intercept of zero. However, this would not be a good fit for all data. What if we want a model like f(x) = x⋅β + b, where b is some nonzero constant? The value b is often called a bias in the context of machine learning, since it “biases” the data towards a specific value even when all the model’s inputs are zero.
I’ve left this consideration for a quick addendum because it turns out adding bias in regression models doesn’t require much effort: It’s equivalent to adding an extra independent variable to the data which always takes the value 1. For example, if we had some two dimensional data and wanted to add bias to our regressor, instead of fitting a model of the form f(x) = kx, we could fit the model
The variable x_1 would represent our original data, and x_2 would be set equal to 1 in every datapoint. A datapoint (x, y) would thus become (x_1, x_2, y) = (x, 1, y). Plugging in x_1 = x and x_2 = 1, Equation 7 reduces to
Here β_1 is the slope and β_2 is the bias. We can fit this model using multiple linear regression. For higher dimensional data, the process works similarly.
Conclusion
There we have it! That’s linear regression from scratch. I hope you enjoyed and learned something. I’d welcome any feedback and constructive criticism. Stay tuned for my next post, where I’ll cover linear classification.
All the code is available on github.
