1.2 - Multivariate Linear Regression

Once you have understood Simple Linear Regression, we can go a bit deeper. In a simple linear regression problem, the input data only has one feature. However, we can consider a scenario where the input data has multiple features. This problem is called Multivariate Linear Regression.

Figure 1:Multivariate perceptron with multiple inputs and one output.

We assume that the true unknown function $f$ that maps the relationship between the input and output is:

y = f\left( \mathbf{x} \right) + \epsilon

(1)

where $\epsilon$ is an intrinsic noise independent of $\mathbf{x}$ . Note that the input $\mathbf{x}$ is a vector.

The goal of multivariate linear regression is similar to simple linear regression, estimate $f$ by a linear approximation $\hat{f}$ such that:

\hat{y} = \hat{f}(\mathbf{x}) \approx y

(2)

Purpose of this Notebook:

Create a dataset for multivariate linear regression task
Create our own Perceptron class from scratch
Calculate the gradients from scratch
Implement gradient descent from scratch
Train our Perceptron
Compare our Perceptron to PyTorch’s built-in implementation

Setup¶

print('Start package installation...')

Start package installation...

%%capture
%pip install torch
%pip install scikit-learn

print('Packages installed successfully!')

Packages installed successfully!

import torch
from torch import nn

from platform import python_version
python_version(), torch.__version__

('3.14.7', '2.12.0+cpu')

# Set seeds for reproducibility
device = 'cpu'
torch.manual_seed(12)
if torch.cuda.is_available():
    torch.cuda.manual_seed_all(12)
    device = 'cuda'
    import numpy as np; np.random.seed(12)
import random; random.seed(12)
device

'cpu'

The add_to_class function is used to add new methods to a previously defined class; we do this to gradually enhance the class.

def add_to_class(Class):
    """Register functions as methods in created class."""
    def wrapper(obj): setattr(Class, obj.__name__, obj)
    return wrapper

Dataset¶

Create Dataset¶

The dataset $\mathcal{D}$ consists of input-target pairs $(\mathbf{x}_i, y_i)$ :

\mathcal{D} = \left\{(\mathbf{x}_{1}, y_{1}), \cdots, (\mathbf{x}_{n}, y_{n}) \right\}

(3)

where $\mathbf{x}_{i}$ denotes the $i$ -th input sample and $y_{i}$ its corresponding target value.

The input data $\mathbf{X} \in \mathbb{R}^{n \times d}$ can be represented as a matrix:

\begin{align} \mathbf{X} &= \begin{bmatrix} x_{11} & \cdots & x_{1d} \\ \vdots & \ddots & \vdots \\ x_{n1} & \cdots & x_{nd} \end{bmatrix} \\ &= \begin{bmatrix} \mathbf{x}_{1}^{\top} \\ \vdots \\ \mathbf{x}_{n}^{\top} \end{bmatrix} \end{align}

(4)

where $n$ is the number of samples in the dataset, $d$ is the number of features, and $\mathbf{x}_{i}^{\top} = \begin{bmatrix} x_{i1} & \cdots & x_{id} \end{bmatrix} \in \mathbb{R}^{1 \times d}$ .

The target data $\mathbf{y} \in \mathbb{R}^{n}$ remains unchanged:

\mathbf{y} = \begin{bmatrix} y_{1} \\ \vdots \\ y_{n} \end{bmatrix}

(5)

from sklearn.datasets import make_regression
import random


N: int = 1_000 # number of samples
D: int = 4 # number of features

X, Y = make_regression( # type: ignore
    n_samples=N, 
    n_features=D, 
    n_targets=1,
    n_informative=D - 1, # let's add features as a linear combination of others
    bias=random.random(), # random true bias
    noise=1,
)

X = X.astype(np.float32)
Y = Y.astype(np.float32)

Y = Y.reshape(-1, 1) # add the axis of length 1

X.shape, Y.shape

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 17
     13     bias=random.random(), # random true bias
     14     noise=1,
     15 )
     16 
---> 17 X = X.astype(np.float32)
     18 Y = Y.astype(np.float32)
     19 
     20 Y = Y.reshape(-1, 1) # add the axis of length 1

NameError: name 'np' is not defined

Split Dataset¶

from sklearn.model_selection import train_test_split

X_train, X_valid, Y_train, Y_valid = train_test_split(
    X, Y,
    test_size=0.2,
    random_state=42,
    shuffle=True,
)

X_train.shape, Y_train.shape

X_valid.shape, Y_valid.shape

Remark: We have omitted the use of the test set for the rest of this notebook, as our aim here is to understand how the Perceptron works, not to draw general conclusions about a synthetic dataset. We will use the validation set instead of the test set, as it serves the same purpose in this notebook.

Tensor Dataset¶

We convert our NumPy arrays into PyTorch tensors and wrap them in a TensorDataset to allow for structured iteration.

from torch.utils.data import TensorDataset

train_dataset = TensorDataset(
    torch.from_numpy(X_train),
    torch.from_numpy(Y_train)
)

train_dataset[0] # get the first sample

valid_dataset = TensorDataset(
    torch.from_numpy(X_valid),
    torch.from_numpy(Y_valid)
)

Data Loader¶

The use of DataLoaders has important advantages. The data is stored in main memory, and when we need to feed the model, data samples are assembled into batches and transferred to the GPU for training or evaluation, optimizing the process of accessing data Nouaji et al., 2025.

from torch.utils.data import DataLoader

Let’s set the batch size to 32, which means we have 32 different subsets of samples for training and validation.

BATCH_SIZE: int = 32

train_loader = DataLoader(
    train_dataset, 
    batch_size=BATCH_SIZE, 
    shuffle=False, # we want to ensure determinism
    pin_memory=True,
)

We can check the number of samples per minibatch:

len(train_loader), len(train_dataset) // BATCH_SIZE

And we can get the first minibatch from the training loader:

_ = next(iter(train_loader))
_[0].shape, _[1].shape

valid_loader = DataLoader(
    valid_dataset, 
    batch_size=BATCH_SIZE, 
    shuffle=False,
    pin_memory=True,
)

_ = next(iter(valid_loader))
_[0].shape, _[1].shape

Remark: The training set can be divided evenly into four mini-batches, with each mini-batch containing exactly the same number of examples. Usually, in each training iteration, the training set is shuffled before being divided into mini-batches, and any remaining examples are discarded during that iteration. To enable shuffling during training we need to edit this argument during train loader initilization shuffle=True.

Scratch Multivariate Perceptron¶

Bias and Weight¶

Our model $\hat{y}$ has two trainable parameters $b, \mathbf{w}$ . But note that weight parameter is a vector:

\mathbf{w} \in \mathbb{R}^{d}

(6)

where the elements $w_{j}$ of $\mathbf{w}$ are called weights. And $b \in \mathbb{R}$ .

class MultiLinearRegression:
    def __init__(self, d_features: int):
        self.b = torch.randn(()).to(device)
        self.w = torch.randn(d_features).to(device)

    def copy_params(self, torch_layer: nn.modules.linear.Linear):
        """
        Copy the parameters from a module.linear to this model.

        Args:
            torch_layer: Pytorch module from which to copy the parameters.
        """
        self.b.copy_(torch_layer.bias.detach()[0])
        self.w.copy_(torch_layer.weight.detach()[0, :])

Weighted Sum¶

Our weighted sum is now a dot multiplication between the input data $\mathbf{x}$ and the weight parameter $\mathbf{w}$ :

\begin{align} \hat{y}: \mathbb{R}^{d} &\to \mathbb{R} \\ \mathbf{x} &\mapsto \hat{y}(\mathbf{x}) = b + \mathbf{x}^{\top} \mathbf{w} \end{align}

(7)

Remark: We can add a scalar $b \in \mathbb{R}$ to a vector due to broadcasting mechanism.

For vectorization, given a minibatch input $\mathbf{X} \in \mathbb{R}^{m \times d}$ of $m$ samples and $d$ input features, we can compute the weighted sum as:

\hat{\mathbf{y}} = b + \mathbf{Xw}

(8)

Note: $n$ is the number of samples in the dataset, while $m$ is the number of samples in a minibatch.

The prediction for the $i$ -th sample is:

\hat{y}_{i} = b + \sum_{j=1}^{d} x_{ij} w_{j}

(9)

This formula will be useful for gradient descent.

@add_to_class(MultiLinearRegression)
def predict(self, x: torch.Tensor) -> torch.Tensor:
    """
    Predict the output for input x.

    Args:
        x: Input tensor of shape (m_samples, d_features).

    Returns:
        y_pred: Predicted output tensor of shape (m_samples,).
    """
    return torch.matmul(x, self.w) + self.b

MSE¶

The MSE formula as loss function $\mathcal{L}$ remains unchanged:

\begin{align} \mathcal{L} (\mathbf{w}, b) &= \frac{1}{m} \sum_{i=1}^{m} \left( \hat{y}_{i} - y_{i} \right)^{2} \\ &= \frac{1}{m} \sum_{i=1}^{m} \left( \left(b + \mathbf{x}^{\top}_{i} \mathbf{w} \right) - y_{i} \right)^{2} \end{align}

(10)

@add_to_class(MultiLinearRegression)
def mse_loss(self, y_true: torch.Tensor, y_pred: torch.Tensor):
    """
    MSE loss function between target y_true and y_pred.

    Args:
        y_true: Target tensor of shape (m_samples,).
        y_pred: Predicted tensor of shape (m_samples,).

    Returns:
        loss: MSE loss between predictions and true values.
    """
    return ((y_pred - y_true) ** 2).mean().item()

@add_to_class(MultiLinearRegression)
def evaluate(self, x: torch.Tensor, y_true: torch.Tensor):
    """
    Evaluate the model on input x and target y_true using MSE.

    Args:
        x: Input tensor of shape (m_samples, d_features).
        y_true: Target tensor of shape (m_samples,).

    Returns:
        loss: MSE loss between predictions and true values.
    """
    y_pred = self.predict(x)
    return self.mse_loss(y_true, y_pred)

Gradients¶

To make our model’s parameters update, it is necessary to compute derivatives.

First, determine the derivatives to be computed
Then, ascertain the shape of each derivative
Finally, compute the derivatives

Derivative of MSE with respect to bias:

\frac{\partial \mathcal{L}}{\partial b} = \sum_{p} \frac{\partial \mathcal{L}}{\partial \hat{y}_{p}} \frac{\partial \hat{y}_{p}}{\partial b}

(11)

and derivative of MSE with respect to weight:

\frac{\partial \mathcal{L}}{\partial w_{q}} = \sum_{p} \frac{\partial \mathcal{L}}{\partial \hat{y}_{p}} \frac{\partial \hat{y}_{p}}{\partial w_{q}}

(12)

where the shape of each derivative is:

\frac{\partial \mathcal{L}}{\partial b} \in \mathbb{R}, \frac{\partial \mathcal{L}}{\partial \mathbf{w}} \in \mathbb{R}^{d}, \frac{\partial \mathcal{L}}{\partial \hat{\mathbf{y}}} \in \mathbb{R}^{m}, \frac{\partial \hat{\mathbf{y}}}{\partial b} \in \mathbb{R}^{m}, \frac{\partial \hat{\mathbf{y}}}{\partial \mathbf{w}} \in \mathbb{R}^{m \times d}

(13)

MSE Derivative¶

Derivative of MSE with respect to predicted data is:

\begin{align} \frac{\partial \mathcal{L}}{\partial \hat{y}_{p}} &= \frac{\partial}{\partial \hat{y}_{p}} \left( \frac{1}{m} \sum_{i=1}^{m} \left( \hat{y}_{i} - y_{i} \right)^{2} \right) \\ &= \frac{1}{m} \sum_{i=1}^{m} \frac{\partial}{\partial \hat{y}_{p}} \left( \left( \hat{y}_{i} - y_{i} \right)^{2} \right) \\ &= \frac{2}{m} \sum_{i=1}^{m} \left( \hat{y}_{i} - y_{i} \right) \frac{\partial \hat{y}_{i}}{\partial \hat{y}_{p}} \\ &= \frac{2}{m} \sum_{i=1}^{m} \left( \hat{y}_{i} - y_{i} \right) \delta_{ip} \\ &= \frac{2}{m} \left( \hat{y}_{p} - y_{p} \right) \end{align}

(14)

for $p = 1, \ldots, m$ .

The vectorized form is:

\frac{\partial \mathcal{L}}{\partial \hat{\mathbf{y}}} = \frac{2}{m} \left( \hat{\mathbf{y}} - \mathbf{y} \right)

(15)

Weighted Sum Derivative¶

Derivative of weighted sum with respect to bias is:

\begin{align} \frac{\partial \hat{y}_{p}}{\partial b} &= \frac{\partial}{\partial b} \left( b + \mathbf{x}_{p}^{\top} \mathbf{w} \right) \\ &= 1 \end{align}

(16)

for $p = 1, \ldots, m$ .

The vectorized form is:

\frac{\partial \hat{\mathbf{y}}}{\partial b} = \mathbf{1}

(17)

where $\mathbf{1} \in \mathbb{R}^{m}$ .

Derivative of weighted sum with respect to weight is:

\begin{align} \frac{\partial \hat{y}_{p}}{\partial w_{q}} &= \frac{\partial}{\partial w_{q}} \left( b + \mathbf{x}_{p}^{\top} \mathbf{w} \right) \\ &= \frac{\partial}{\partial w_{q}} \left(\mathbf{x}_{p}^{\top} \mathbf{w} \right) \\ &= \frac{\partial}{\partial w_{q}} \left( x_{p1}w_{1} + \ldots + x_{pq}w_{q} + \ldots + x_{pd}w_{d} \right) \\ &= \frac{\partial}{\partial w_{q}} \left( x_{pk} w_{k} \right) \\ &= x_{pk} \delta_{kq} \\ &= x_{pq} \end{align}

(18)

for $p = 1, \ldots, m$ , and $q = 1, \ldots, d$ .

Vectorizing for all $q = 1, \ldots, d$ :

\frac{\partial \hat{y}_{p}}{\partial \mathbf{w}} = \mathbf{x}_{p}^{\top} \in \mathbb{R}^{1 \times d}

(19)

Vectorizing for all $p = 1, \ldots, m$ :

\frac{\partial \hat{\mathbf{y}}}{\partial \mathbf{w}} = \mathbf{X} \in \mathbb{R}^{m \times d}

(20)

Full Chain Rule¶

Derivative of MSE with respect to bias is:

\begin{align} \frac{\partial \mathcal{L}}{\partial b} &= \sum_{p} {\color{Cyan} \frac{\partial \mathcal{L}}{\partial \hat{y}_{p}}} {\color{Orange} \frac{\partial \hat{y}_{p}}{\partial b}} \\ &= \sum_{p} {\color{Cyan} \frac{2}{m} \left( \hat{y}_{p} - y_{p} \right)} {\color{Orange} 1} \\ &= \frac{2}{m} \sum_{p} \left( \hat{y}_{p} - y_{p} \right) 1 \\ &= \frac{2}{m} \left( \hat{\mathbf{y}} - \mathbf{y} \right)^{\top} \mathbf{1} \end{align}

(21)

Derivative of MSE with respect to weight:

\begin{align} \frac{\partial \mathcal{L}}{\partial w_{q}} &= \sum_{p} {\color{Cyan} \frac{\partial \mathcal{L}}{\partial \hat{y}_p}} {\color{Magenta} \frac{\partial \hat{y}_{p}}{\partial w_{q}}} \\ &= \sum_{p} {\color{Cyan} \frac{2}{m} \left(\hat{y}_{p} - y_{p} \right)} {\color{Magenta} x_{pq}} \\ &= \frac{2}{m} \sum_{p} \left(\hat{y}_{p} - y_{p} \right) x_{pq} \\ &= \frac{2}{m} \left( \mathbf{x}_{:,q} \right)^{\top} \left( \hat{\mathbf{y}} - \mathbf{y} \right) \end{align}

(22)

for $q = 1, \ldots, d$ , where $\mathbf{x}_{:,q} = \begin{bmatrix} x_{1q} & \cdots & x_{md} \end{bmatrix}^{\top} \in \mathbb{R}^{m \times 1}$ .

Note: $\mathbf{x}_{:,q}$ is the column vector of the feature $q$ -th of each sample in the minibatch.

Vectorized form is:

\begin{align} \frac{\partial \mathcal{L}}{\partial \mathbf{w}} &= \frac{2}{m} \mathbf{X}^{\top} \left( \hat{\mathbf{y}} - \mathbf{y} \right) \end{align}

(23)

Final Gradients¶

\nabla_{b} \mathcal{L} = \frac{2}{m} \left( \hat{\mathbf{y}} - \mathbf{y} \right)^{\top} \mathbf{1}

(24)

\nabla_{\mathbf{w}} \mathcal{L} = \frac{2}{m} \mathbf{X}^{\top} \left( \hat{\mathbf{y}} - \mathbf{y} \right)

(25)

Parameters Update¶

\begin{align} b &\leftarrow b -\eta \nabla_{b} \mathcal{L} \\ &= b -\eta \left( \frac{2}{m} (\hat{\mathbf{y}} - \mathbf{y})^{\top} \mathbf{1} \right) \end{align}

(26)

\begin{align} \mathbf{w} &\leftarrow \mathbf{w} -\eta \nabla_{\mathbf{w}} \mathcal{L} \\ &= \mathbf{w} -\eta \left( \frac{2}{m} \mathbf{X}^{\top} (\hat{\mathbf{y}} - \mathbf{y}) \right) \end{align}

(27)

where $\eta > 0$ is called learning rate.

@add_to_class(MultiLinearRegression)
@torch.inference_mode()
def update(self, x: torch.Tensor, y_true: torch.Tensor, 
           y_pred: torch.Tensor, lr: float):
    """
    Update the model parameters.

    Args:
       x: Input tensor of shape (m_samples, d_features).
       y_true: Target tensor of shape (m_samples,).
       y_pred: Predicted output tensor of shape (m_samples,).
       lr: Learning rate. 
    """
    delta = 2 * (y_pred - y_true) / len(y_true)
    self.b -= lr * delta.sum()
    self.w -= lr * torch.matmul(x.T, delta)

Gradient Descent¶

We implement minibatch gradient descent to fit our model’s parameters.

@add_to_class(MultiLinearRegression)
def fit(self, train_: torch.utils.data.dataloader.DataLoader, 
        epochs: int, lr: float, 
        valid_: torch.utils.data.dataloader.DataLoader):
    """
    Fit the model using gradient descent.
    
    Args:
        train_: Training dataloader with tensors of shape (m_samples, d_features).
        epochs: Number of epochs to fit.
        lr: learning rate.
        valid_: Validation dataloader with tensors of shape (m_valid_samples, d_features).
    """
    for epoch in range(epochs):
        # training epoch
        running_loss = 0.0
        for batch_x, batch_y in train_:
            # move only the current batch to VRAM
            batch_x = batch_x.to(device, non_blocking=True)
            batch_y = batch_y.to(device, non_blocking=True).squeeze(1) 

            # make predictions
            y_pred = self.predict(batch_x)

            running_loss += self.mse_loss(
                batch_y, y_pred
            )
            
            self.update(
                batch_x, batch_y, 
                y_pred, lr
            )
        avg_loss = running_loss / len(train_)
        
        # validation epoch
        running_vloss = 0.0
        # disable gradient computation and reduce memory consumption.
        with torch.no_grad():
            for vbatch_x, vbatch_y in valid_:
                vbatch_x = vbatch_x.to(device, non_blocking=True)
                vbatch_y = vbatch_y.to(device, non_blocking=True).squeeze(1) 

                vy_pred = self.predict(vbatch_x)

                running_vloss += self.mse_loss(
                    vbatch_y, vy_pred
                )
        avg_loss_v = running_vloss / len(valid_)

        print(f'epoch: {epoch} - MSE: {avg_loss:.4f} - vMSE: {avg_loss_v:.4f}')

Scratch vs Torch.nn¶

PyTorch Module¶

class TorchLinearRegression(nn.Module):
    def __init__(self, n_features):
        super().__init__()
        self.layer = nn.Linear(n_features, 1)
        self.loss = nn.MSELoss()

    def forward(self, x):
        return self.layer(x)

    @torch.inference_mode()
    def evaluate(self, x, y):
        self.eval()
        y_pred = self(x)
        return self.loss(y_pred, y).item()
    
    def fit(self, train_, epochs, lr, valid_):
        optimizer = torch.optim.SGD(self.parameters(), lr=lr)
        
        for epoch in range(epochs):
            self.train() # use model.train() when training
            running_loss = 0.0 # train loss
            for batch_x, batch_y in train_:
                batch_x = batch_x.to(device, non_blocking=True)
                batch_y = batch_y.to(device, non_blocking=True)

                y_pred = self(batch_x)

                loss = self.loss(y_pred, batch_y)

                optimizer.zero_grad()
                loss.backward()
                optimizer.step()

                running_loss += loss.item()
            avg_loss = running_loss / len(train_)

            # validation epoch
            self.eval() # set the model to evaluation mode
            with torch.no_grad():
                running_vloss = 0.0
                for vbatch_x, vbatch_y in valid_:
                    vbatch_x = vbatch_x.to(device, non_blocking=True)
                    vbatch_y = vbatch_y.to(device, non_blocking=True)

                    vy_pred = self(vbatch_x)
                    vloss = self.loss(vy_pred, vbatch_y)
                    running_vloss += vloss.item()
            avg_loss_v = running_vloss / len(valid_)

            print(f'epoch: {epoch} - MSE: {avg_loss:.4f} - vMSE: {avg_loss_v:.4f}')

We initialized both models:

torch_model = TorchLinearRegression(D).to(device)

model = MultiLinearRegression(D)

Eval¶

We use a $L_{2}$ norm between PyTorch model and our Scratch model as parameter discrepancy.

def l2(pred, true_): 
    if isinstance(pred, (float, int)) and isinstance(true_, (float, int)):
        return abs(pred - true_)
    return torch.linalg.norm(pred - true_).item()

Predictions pre-copy¶

Let compare the predictions of our model and PyTorch implementation using $L_{2}$ norm:

x_valid, y_valid = valid_dataset[:]
x_valid = x_valid.to(device, non_blocking=True)
y_valid = y_valid.to(device, non_blocking=True)

l2(
    model.predict(x_valid), 
    torch_model(x_valid).squeeze(1)
)

They differ considerably because each model has its own parameters initialized randomly and independently of the other model.

Copy Parameters¶

We copy the values of the PyTorch model parameters to our model.

model.copy_params(torch_model.layer)

Predictions post-copy¶

We measure the difference between the predictions of both models again.

l2(
    model.predict(x_valid),
    torch_model(x_valid).squeeze(1)
)

Loss¶

l2(
    model.evaluate(x_valid, y_valid.squeeze(1)), 
    torch_model.evaluate(x_valid, y_valid)
)

Training¶

We are going to train both models using the same hyperparameter values. If our model is well designed, then starting from the same parameters it should arrive at the same parameters’ values as the PyTorch model after training.

LR: float = 0.01 # learning rate
EPOCHS: int = 16 # number of epochs

model.fit(
    train_loader,
    EPOCHS, LR,
    valid_loader
)

torch_model.fit(
    train_loader,
    EPOCHS, LR,
    valid_loader
)

Predictions after Training¶

l2(
    model.predict(x_valid), 
    torch_model(x_valid).squeeze(1) 
)

Values such as 1.78e-07 are indistinguishable from zero in float32.

Bias Comparison¶

We directly measure the difference between the bias values of both models.

l2(
    model.b.clone(),
    torch_model.layer.bias.detach()[0]
)

Weight Comparison¶

And measure the difference between the weight values of both models.

l2(
    model.w.clone(),
    torch_model.layer.weight.detach().squeeze(0)
)

Using this $L_{2}$ norm as a metric, we can conclude that our implementation of Multivariate Linear Regression is equivalent to PyTorch’s built-in implementation.