Linear Regression from scratch

Table of Contents

Algorithm and implementation
#

Linear Regression formula can be replicated on numpy with X @ w + b (or np.dot(X,w)).

To implement gradient descent as an optimizer, we need to calculate the following:

errors = predictions - y
The loss is RSS , evaluated with MSE.
loss = np.mean(errors - y)**2 In this case, we will actually be using the half-MSE to have similarities with Logistic Regression from scratch. $$L=\frac{1}{2N}\sum_i(\hat y_i-y_i)^2.$$Half-MSE has the same minimizer as MSE and its 1/2 cancels the derivative’s 2 giving the same gradient template as Logistic Regression: $$dw=X^T errors/N, \qquad db=mean(errors).$$ Using “pure” MSE the derivative would have been: $$dw=2X^T errors/N, \qquad db=mean(errors).$$
weights w and biases are initialized as np.zeros and np.mean(y) respectively, since the weights act as a sort of feature importance and the bias works as an intercept . After choosing a learning rate, at further steps w and b will be updated as:
- w-=dw * learning_rate
- b-=db * learning_rate

Class structure
#

main attributes:
- learning_rate
- n_iterations
- weights w_ and bias b_ initialized as None
fit method
predict method

import numpy as np

class LinearRegression:
	def __init__(self, learning_rate=0.05, n_iterations=1000):
		self.learning_rate = learning_rate
		self.n_iterations = n_iterations
		self.w_ = None
		self.b_ = None
		
	def fit(self, X, y):
		pass
		
	def predict(self, X):
		pass

Training
#

The fit method can be implemented following these steps:

input shapes setup with X of shape (N,D) and y of shape (N,)
weights and bias initialization
“PELGU loop”:
- Predictions: X @ self.w_ + self.b_ (shape (N,))
- Errors: errors = predictions - y (shape (N,))
- Loss: loss = 0.5 * np.mean(errors ** 2) (scalar).
- Gradient:
  - dw = X.T @ errors / n_samples (shape (D,))
  - db = errors.mean() (scalar)
- Update:
  - self.w_ = self.learning_rate * dw (shape (D,))
  - self.b_ = self.learning_rate * db (scalar)

def fit(self, X, y):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = float(y.mean())
        self.loss_history_ = []

        for _ in range(self.n_iterations):
            predictions = X @ self.w_ + self.b_       # P: (N,)
            errors = predictions - y                  # E: (N,)
            loss = 0.5 * np.mean(errors ** 2)       # L: scalar

            dw = X.T @ errors / n_samples              # G: (D,)
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
           
        return self

Inference
#

Inference is basically the P step of the PELGU loop above applied on the predictions set after training:
predictions = X @ self.w_ + self.b_

def predict(self, X):
        if self.w_ is None:
            raise RuntimeError("fit must be called before predict")
        X = np.asarray(X, dtype=float)
        if X.ndim != 2 or X.shape[1] != len(self.w_):
            raise ValueError("X has the wrong shape")
        predictions = X @ self.w_ + self.b_
        return predictions

Features scaling
#

To make learning more robust, Standardization can be used.

Adding regularization
#

A common step when implementing Linear Regression is to add regularization to tackle situations in which a model can fit training noise with unnecessarily large weights (overfitting), especially when features are correlated or data is limited.

We can penalize the model for large weights adding a sort of “penalty” to the loss we saw in the base implementation:
regularization_loss = 0.5 * l2 * sum(w ** 2)

Its gradient is regularization_dw = l2 * w, which pulls every weight toward zero, with l2 controlling the “strength of the pull”.

This affects the total loss and the weights gradient dw. Bias b does not need to be regularized because it does not measure the dependence on a feature.
This step usually comes after scaling.

The fit method of the LinearRegression class would look like:

def fit(self, X, y, l2):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = float(y.mean())
        self.loss_history_ = []

        for _ in range(self.n_iterations):
            predictions = X @ self.w_ + self.b_       # P: (N,)
            errors = predictions - y                  # E: (N,)
            
            data_loss = 0.5 * np.mean(errors**2)            # L: scalar
            regularization_loss = 0.5 * l2 * np.sum(self.w_ ** 2)
            loss = data_loss + regularization_loss

            data_dw = X.T @ errors / n_samples              # G: (D,)
            regularization_dw = l2 * self.w_
            dw = data_dw + regularization_dw
            
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
           
        return self

Adding early-stopping
#

When an algorithm converges, further iterations barely change the parameters. A fixed iterations count is safe but it may waste work. The idea behind early stopping is to stop when the largest absolute update of w and b is below the tolerance argument value.
If, on the contrary, learning rate diverges we must stop the training as we are likely finding nulls in the loss.

The fit method changes as follows when adding early stopping:

def fit(self, X, y, l2, tolerance):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = float(y.mean())
        self.loss_history_ = []
        self.tolerance = tolerance

        for _ in range(self.n_iterations):
            predictions = X @ self.w_ + self.b_       # P: (N,)
            errors = predictions - y                  # E: (N,)
            
            data_loss = 0.5 * np.mean(errors**2)            # L: scalar
            regularization_loss = 0.5 * l2 * np.sum(self.w_ ** 2)
            loss = data_loss + regularization_loss

            data_dw = X.T @ errors / n_samples              # G: (D,)
            regularization_dw = l2 * self.w_
            dw = data_dw + regularization_dw
            
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
            
            #early-stopping based on weights update
            largest_update = max(np.max(np.abs(self.learning_rate * dw )), abs(self.learning_rate * db))
            if largest_update<=self.tolerance:
	            break
	            
	        if not np.isfinite(loss): raise FloatingPointError('learning rate diverged')
           
        return self

Another way to add early-stopping relies on the loss improvements rather than the weights change, and can be implemented introducing a patience argument, that checks the previous loss values and if there’s no improvement over them, the learning stops:

def fit(self, X, y, l2, tolerance, patience):
        X = np.asarray(X, dtype=float)
        y = np.asarray(y, dtype=float)
        if X.ndim != 2 or y.ndim != 1 or len(X) == 0 or len(X) != len(y):
            raise ValueError("X must be (N,D) and y must be (N,), with matching N > 0")

        n_samples, n_features = X.shape
        self.w_ = np.zeros(n_features)
        self.b_ = float(y.mean())
        self.loss_history_ = []
        self.tolerance = tolerance
        self.patience = patience
        
        best_loss = np.inf

        for _ in range(self.n_iterations):
            predictions = X @ self.w_ + self.b_       # P: (N,)
            errors = predictions - y                  # E: (N,)
            
            data_loss = 0.5 * np.mean(errors**2)            # L: scalar
            regularization_loss = 0.5 * l2 * np.sum(self.w_ ** 2)
            loss = data_loss + regularization_loss

            data_dw = X.T @ errors / n_samples              # G: (D,)
            regularization_dw = l2 * self.w_
            dw = data_dw + regularization_dw
            
            db = errors.mean()                         # G: scalar
            self.w_ -= self.learning_rate * dw         # U: (D,)
            self.b_ -= self.learning_rate * db         # U: scalar

            self.loss_history_.append(float(loss))
            
            #early-stopping based on loss
            if loss < best_loss - self.tolerance:
	            best_loss = loss
	            no_improve = 0
	        else:
		        no_improve += 1
		        if no_improve >= self.patience:
			        break

        return self

Algorithm and implementation#

Class structure#

Training#

Inference#

Features scaling#

Adding regularization#

Adding early-stopping#