Viaggio nell'Iperspazio

Un Esempio con la Regressione Lineare Multipla

Nel post precedente, abbiamo visto come usare la regressione lineare semplice per modellare una relazione tra la variabile dipendente $y^{(i)}$ e una sola variabile indipendente $x^{(i)}$. La relativa funzione di ipotesi è la funzione predittiva che approssima la funzione target che ha generato le osservazioni del dataset di training, definita come:

Se vogliamo invece definire la funzione predittiva come la relazione tra una variabile dipendente $y$ e più variabili indipendenti $x^{(i)}_1$, …, $x^{(i)}_n$, dobbiamo utilizzare la regressione lineare multipla, la cui funzione di ipotesi è definita come:

dove:

• $\hat{y}^{(i)}$ è il valore continuo di predizione (o regressione), ritornato dalla funzione.
• $\textbf{x}^{(i)}$ è il vettore feature che contiene le variabili indipendenti $x^{(i)}_1$, $\dots$, $x^{(i)}_n$, o feature, dell’$i$-esima osservazione.
• $w_1$, $\dots$, $w_n$ sono i coefficienti associati a ciascuna feature $x^{(i)}_j$.
• $b$ è l’intercetta, ovvero il valore di $\hat{y}^{(i)}$ quando tutte le variabili indipendenti $x^{(i)}_j$ sono uguali a zero, oppure quando tutti i parametri $w_j$ sono uguali a 0, dove $1 \le j \le n$. Il termine $b$ è chiamato anche bias.

Cerchiamo di capire il significato dell’equazione $(2)$ un passo alla volta. Ora, cosa accade se le osservazioni sono definite da due variabili indipendenti $x^{(i)}_1$ e $x^{(i)}_2$, invece che una? In questo caso, dovremo mappare $x^{(i)}_1$ e $x^{(i)}_2$ su due assi dimensionali distinti, ed il valore $\hat{y}^{(i)}$, ritornato dalla funzione di ipotesi:

verrà mappato su un ulteriore asse dimensionale. L’equazione $(3)$ non verrà più visualizzata come una retta in uno spazio vettoriale 2D, ma come un piano in uno spazio 3D:

Infatti, se riscriviamo $(3)$ come:

notiamo che $(4)$ è compatibile con l’equazione di un piano nella sua forma generale:

dove:

• $a$, $b$ e $c$ sono le componenti del vettore normale al piano $\vec{n} = \left(a, b, c\right)$, che determinano l’orientamento del piano.
• $x$, $y$ e $z$ sono le coordinate di un punto nello spazio 3D.
• $d$ è una costante che determina la distanza del piano dall’origine del sistema di riferimento, nella direzione di $\vec{n}$.
• i punti che giacciono esattamente sul piano soddisfano tale equazione.
• confrontando $(5)$ con $(4)$, notiamo che $cz = -\hat{y}^{(i)}$, quindi $\hat{y}^{(i)} = -cz$.

Come già abbiamo visto con la regressione lineare semplice, la funzione di ipotesi, dopo il training del modello, ci permette di fare previsioni di valori continui, usando dati mai visti in precedenza. A livello geometrico, le previsioni per un qualsiasi vettore feature $\textbf{x}^{(i)} = \left(x^{(i)}_1, x^{(i)}_2\right)$, generano sempre punti $\left(x^{(i)}_1, x^{(i)}_2, \hat{y}^{(i)}\right)$ appartenenti al piano $(3)$:

L’Iperpiano e la Funzione di Ipotesi

Abbiamo appena visto cosa indichi l’equazione della funzione di ipotesi $(2)$, quando vengono utilizzate due variabili indipendenti $x^{(i)}_1$ e $x^{(i)}_2$. Possiamo generalizzare questo comportamento? Per capirlo, dobbiamo introdurre il concetto di iperspazio e iperpiano.

In geometria, un iperspazio si riferisce a uno spazio di $k$ dimensioni che può contenere oggetti di $k-1$ dimensioni, chiamati iperpiani. Un iperpiano è un sottospazio di una dimensione inferiore a quella dello spazio in cui esso è contenuto. Ad esempio, l’iperpiano di uno spazio 2D è una retta 1D, mentre l’iperpiano di uno spazio 3D è un piano 2D. L’iperpiano in uno spazio a $k$ dimensioni può essere definito da un’equazione lineare della forma:

dove:

• $w_1$, $\dots$, $w_k$ sono le coordinate del vettore $\mathbf{w}$, normale all’iperpiano.
• $x^{(i)}_1$, $\dots$, $x^{(i)}_k$ sono le coordinate di un punto $\mathbf{x}^{(i)}$ in uno spazio $k$-dimensionale.
• $b$ è il bias dell’iperpiano, e indica la distanza dell’iperpiano dall’origine, nella direzione di $\textbf{w}$.
• i punti o le osservazioni che giacciono esattamente sull’iperpiano soddisfano tale equazione.

Possiamo confrontare l’equazione dell’iperpiano $(6)$ con la funzione di ipotesi $(2)$ per alcuni casi speciali, e capire se ci sia una relazione tra le due. In uno spazio 2D, le seguenti equazioni rappresentano entrambe una retta:

Infatti, è possibile riscrivere l’equazione dell’iperpiano $(7)$ come l’equazione di una retta, dove $-\frac{w_1}{w_2}$ è il coefficiente angolare, e $-\frac{b}{w_2}$ l’intercetta:

Possiamo inoltre moltiplicare entrambe le parti di $(9)$ per $-w_2$, sostituire $-w_2x^{(i)}_2$ con $\hat{y}^{(i)}$, e trasformare così $(7)$ nella funzione di ipotesi $(8)$:

Continuando il confronto, in uno spazio 3D, le seguenti equazioni rappresentano un piano:

In uno spazio 4D, le seguenti equazioni rappresentano un iperpiano 3D:

Considerando $k = n + 1$, possiamo dire che l’equazione dell’iperpiano e la funzione di ipotesi sono fondamentalmente la stessa equazione lineare, anche se esse vengono usate in contesti diversi. Riscrivendo $(6)$ e sostituendo $w_kx^{(i)}_k$ con $-\hat{y}^{(i)}$, otteniamo infatti $(2)$:

I Vantaggi della Notazione Vettoriale

In uno spazio di $n+1$ dimensioni, l’equazione della funzione di ipotesi $(2)$, usata per calcolare il valore di predizione associato ad un singolo vettore feature, può essere scritta in modo più compatto utilizzando la seguente notazione vettoriale:

dove:

• $\textbf{x}^{(i)} = \left[x_1, x_2, \cdots, x_n \right]^{\intercal}$ è l’$i$-mo vettore feature.
• $\textbf{w} = \left[w_1, w_2, \cdots, w_n \right]^{\intercal}$ è il vettore dei pesi, o coefficienti, associati a ciascuna feature $x^{(i)}_j$.
• l’operatore $\cdot$ è il prodotto vettoriale tra $\textbf{w}$ e $\textbf{x}^{(i)}$, definito come: $\textbf{w} \cdot \textbf{x}^{(i)} = w_1x^{(i)}_1 + w_2x^{(i)}_2 + \dots + w_nx^{(i)}_n$.
• $b$ è l’intercetta.

In modo ancora più compatto, per calcolare predizioni per più vettori feature, è possibile definire $(2)$ come:

dove:

• $\hat{\textbf{y}}$ è il vettore colonna che contiene tutti i valori di regressione associati ai singoli vettori feature $\textbf{x}^{(i)}$:

  $$ \hat{\textbf{y}} = \begin{bmatrix} \hat{y}^{(1)}\\ \hat{y}^{(2)}\\ \vdots\\ \hat{y}^{(m)}\\ \end{bmatrix}\\ $$

• $\textbf{w}^{\intercal}$ indica la trasposizione del vettore colonna $\textbf{w}$ in un un vettore riga, operazione necessaria per effettuare la moltiplicazione matriciale con $\textbf{X}$.

• $\textbf{X}$ è la matrice di dimensione $(m, n)$, contenente tutti i vettori feature $n$-dimensionali:

  $$ \quad \textbf{X} = \begin{bmatrix} \left(\textbf{x}^{(1)}\right)^{\intercal}\\ \left(\textbf{x}^{(2)}\right)^{\intercal}\\ \vdots\\ \left(\textbf{x}^{(m)}\right)^{\intercal}\\ \end{bmatrix} = \begin{bmatrix} x^{(1)}_1 & x^{(1)}_2 & \dots & x^{(1)}_n\\ x^{(2)}_1 & x^{(2)}_2 & \dots & x^{(2)}_n\\ \vdots & \vdots & \vdots & \vdots\\ x^{(m)}_1 & x^{(m)}_2 & \dots & x^{(m)}_n\\ \end{bmatrix} $$

Possiamo notare che è possibile semplificare $(11)$, eliminando la trasposizione di $\textbf{w}$, scambiando $\textbf{w}$ con $\textbf{X}$:

Dall’equazione $(12)$ possiamo notare che la $i$-esima riga della matrice risultante contiene il valore di predizione associato al vettore feature $\textbf{x}^{(i)}$.

Derivata Parziale della Funzione di Costo, Reprise

Consideriamo l’$i$-mo vettore di feature $\mathbf{x}^{(i)}$, di dimensione $(n, 1)$, e il peso $w_j$. Abbiamo visto come la derivata della funzione di costo $\text{MSE}$ rispetto a $w_j$, possa essere calcolata come:

dove:

• $\hat{y}^{(i)} - y^{(i)}$ è uno scalare che rappresenta l’errore per la singola osservazione $\mathbf{x}^{(i)}$.
• la somma scorre su tutti gli esempi (da $1$ a $m$).

Anche la derivata parziale di $\text{MSE}$ rispetto al vettore $\mathbf{w}$ può essere espressa in formato vettoriale, in maniera compatta. Se mettiamo insieme tutti i vettori $\mathbf{x}^{(i)}$ riga per riga, ottieniamo la matrice $\textbf{X}$:

Se mettiamo insieme tutti i valori label, e le relative predizioni, ottieniamo i due vettori colonna $\mathbf{y}$ e $\hat{\mathbf{y}}$:

Il vettore di errore può essere espresso come la differenza tra i due vettori:

In relazione a $(13)$, la seguente sommatoria:

è equivalente al prodotto:

che rappresenta la $j$-ma riga della seguente moltiplicazione matriciale, che ritorna un vettore colonna di dimensioni $(n, 1)$:

Quindi, volendo ottenere il gradiente della funzione di costo per tutti gli $n$ valori del vettore $\mathbf{w}$, possiamo moltiplicare $(14)$ per il coefficiente $\frac{2}{m}$, ottenendo:

Similmente, possiamo esprimere la derivata parziale di $\text{MSE}$ rispetto a $b$ come:

dove $\mathbf{1}$ è un vettore colonna di valori $1$, di dimensione $\left(m,1\right)$, che dobbiamo usare per sommare tra loro gli $m$ valori $\left( \hat{\mathbf{y}}^{(i)} − \mathbf{y}^{(i)} \right)$:

Il risultato è uno scalare, perché il bias $b$ è unico per tutte le osservazioni $\mathbf{x}^{(i)}$.

Implementazione del Modello di Regressione Lineare Multipla

Vediamo ora come implementare un modello di regressione lineare multipla, utilizzando quanto abbiamo appreso relativamente alla funzione di ipotesi $(2)$, alla sua versione vettoriale $(12)$, e alle derivate della funzione di costo $(15)$ e $(16)$:

import numpy as np

def target_function(n_examples, noise_coeff, seed=42):
    """This is the UNKNOWN real target function. Do not read its code!"""
    rng = np.random.default_rng(seed)
    X = np.c_[rng.random(n_examples), rng.random(n_examples)]
    w = np.array([0.25, 0.4])
    b = 3.2
    noise = (rng.random(n_examples) - 0.5) * noise_coeff
    y = X @ w + b + noise
    return X, y.reshape(-1, 1)

def predict(X, w, b):
    """Compute predictions for each value in the input array X, giving the weight w and bias b."""

    # y_hat = w_1 * x_1 + w_2 * x_2 + ... + w_n * x_n + b
    #       = w^T @ X + b
    #       = X @ w + b
    return X @ w + b

def loss(X: np.ndarray, y: np.ndarray, w: float, b: float):
    """Compute the MSE cost value for the input examples."""

    # Mean Squared Error (MSE) = 1/m * sum(i=1 to m) { (ŷ^(i) - y^(i))^2 }
    #                          = 1/m * (1^T)(ŷ - y)^2
    #
    y_hat = predict(X, w, b)                   # shape: (n_examples, 1)
    error = y_hat - y                          # shape: (n_examples, 1)
    squared_error = error ** 2                 # shape: (n_examples, 1)
    loss: float = np.average(squared_error)    # float
    return loss

def gradient(X: np.ndarray, y: np.ndarray, w: float, b: float):
    """Compute the derivative of the loss function L(w, b) with respect to 'w' and 'b'."""

    # ∂L/∂w = 2/m * sum(i=1 to m) { x^(i) * (ŷ^(i) - y^(i)) }
    #       = 2/m * (X.T @ (ŷ - y))
    #
    # ∂L/∂b = 2/m * sum(i=1 to m) { ŷ^(i) - y^(i) }
    #       = 2/m * (1^T)(ŷ - y)
    #
    y_hat = predict(X, w, b)                 # shape: (n_examples, 1)
    error = y_hat - y                        # shape: (n_examples, 1)
    grad_w = 2 * (X.T @ error) / X.shape[0]  # shape: (n_features, 1)
    grad_b = 2 * np.mean(error)              # float
    return grad_w, grad_b

def train(X: np.ndarray, y: np.ndarray, iterations: int, learning_rate: float):
    """Train the model and return the weights and the bias that minimize the MSE Loss Function."""

    # initialize the weights and bias with values taken from a normal distribution
    rng = np.random.default_rng()
    w = rng.normal(loc=0.0, scale=1.0, size=(X.shape[1], 1))  # w shape: (n_features, 1)
    b = rng.normal(size=None)

    # loop through all the iterations...
    for i in range(iterations):

        # calculate the gradient vectors
        grad_w, grad_b = gradient(X, y, w, b)       # grad_w shape: (n_features, 1); grad_b: float

        # calculate the steps in the direction of the w-axis and b-axis axis
        # used to reach the flat point on the loss curve for this iteration
        step_w = grad_w * learning_rate             # step_w shape: (n_features, 1)
        step_b = grad_b * learning_rate             # float

        # update the parameters
        w = w - step_w                              # w shape: (n_features, 1)
        b = b - step_b                              # float

    return w, b

Aggiungiamo inoltre un metodo di plot per mostrare il dataset ed il piano di ipotesi definito dal modello che approssima la distribuzione delle osservazioni:

import plotly.graph_objects as go

def plot_dataset(X, y, X_pred, w, b, show_hypothesis=True, show_predictions=False):
    """Plot the dataset in a 3D space."""
    fig = go.Figure()

    camera = dict(
        up=dict(x=0, y=0, z=1),
        center=dict(x=0, y=0, z=-0.2),
        eye=dict(x=1.50, y=-1.25, z=1.2))

    fig.layout = go.Layout(
        scene=dict(camera=camera,
                   xaxis_title="x1",
                   yaxis_title="x2",
                   zaxis_title="h"))

    marker_size = 6
    data_points_opacity = 1.0

    # plot hypothesis plane
    if show_hypothesis:
        steps = 10
        x1 = np.linspace(X[:, 0].min(), X[:, 0].max(), num=steps)
        x2 = np.linspace(X[:, 1].min(), X[:, 1].max(), num=steps)
        X1, X2 = np.meshgrid(x1, x2)
        z = np.array([[predict(np.array([x1, x2]), w, b)]
                      for x1, x2 in zip(np.ravel(X1), np.ravel(X2))])
        Z = z.reshape((steps, steps))
        fig.add_surface(x=X1, y=X2, z=Z, colorscale='Reds_r', name="Hypothesis Function")
        fig.update_traces(showscale=False, opacity=0.7, showlegend=True)

    # plot predictions
    if show_predictions:
        y_offset = 0.01
        y_hat = predict(X_pred, w, b).flatten()
        fig.add_scatter3d(x=X_pred[:, 0], y=X_pred[:, 1], z=y_hat + y_offset, name="Predictions",
            mode="markers", marker=dict(size=marker_size, color="green",
                line=dict(width=2, color="black"), opacity=1.0))
        fig.update_layout(showlegend=True)
        data_points_opacity = 0.1

    # plot data points
    fig.add_scatter3d(x=X[:, 0], y=X[:, 1], z=y, name="Data points",
        mode="markers", marker=dict(size=marker_size, color="#5273B2",
            line=dict(width=2, color="black"), opacity=data_points_opacity))
    fig.update_layout(showlegend=True)
    fig.show()

Definito il modello, andiamo ad eseguire la sua funzione di training, che ritorna i valori dei parametri $\textbf{w}$ e $b$. Questi vengono poi passati alla funzione plot_dataset, e servono per visualizzare il piano generato dalla funzione di ipotesi, implementata in predict, assieme ai data point usati per il training. Il grafico risultante dovrebbe essere molto, molto simile a quanto visto in Figura 1:

np.random.seed(42)

# create synthetic dataset
n_samples = 100
X, y = target_function(n_samples, noise_coeff=0.1)

# find the w and b coefficients
w, b = train(X, y, iterations=5000, learning_rate=0.01)
# w: [0.2745 0.4157], b: 3.1824

# plot the dataset in a 3D space
y = y.flatten()
plot_dataset(X, y, None, w, b, show_hypothesis=True, show_predictions=False)

Dopo la fase di training, ricavati i valori di $\textbf{w}$ e $b$, possiamo usare il modello per fare previsioni di valori continui:

# plot predictions for 20 random feature vectors...
rng = np.random.default_rng(123)
X_pred = np.c_[rng.random(20), rng.random(20)]  # shape: (20, 2)
plot_dataset(X, y, X_pred, w, b, show_hypothesis=True, show_predictions=True)

Alcune Considerazioni sul Codice

Facciamo alcune considerazioni sul codice Python usato per implementare il modello di cui sopra. Quando si eseguono operazioni tra array con la libreria Numpy, è molto importante fare un’analisi dimensionale delle strutture dati ndarray utilizzate, ovvero esaminare e manipolare le dimensioni degli array (tramite shape e size) per poter esprimere correttamente le espressioni matematiche in forma vettoriale che si vogliono eseguire.

• La funzione predict usa questa istruzione:

  return X @ w + b
  

  L’operatore @ è un’abbreviazione del metodo numpy.matmul usato per calcolare il prodotto matriciale, o prodotto riga per colonna, tra due array di tipo ndarray. L’operazione X @ w + b esegue due operazioni:

  1. la prima è la moltiplicazione tra la matrice $\textbf{X}$ avente dimensione (n_examples, n_features), e il vettore $\textbf{w}$ avente dimensione (n_features, 1), che ritorna un vettore di dimensione (n_examples, 1).
  2. la seconda è l’utilizzo del broadcasting per aggiungere il valore di $b$ ad ogni elemento del vettore ritornato da $\textbf{Xw}$, di dimensione (n_examples, 1). Vedi l’equazione $(12)$.
  La moltiplicazione di due matrici avviene calcolando il prodotto delle righe della prima matrice per le colonne della seconda matrice. Tale operazione è definita soltanto se la prima matrice ha un numero di colonne uguale al numero di righe della seconda matrice, e ritorna una matrice con il numero di righe della prima matrice e il numero di colonne della seconda. In questo caso, il vettore colonna $\textbf{w}$ è considerato come una matrice con una sola colonna.

• La funzione gradient usa la seguente istruzione:

  grad_w = 2 * (X.T @ error) / X.shape[0]  # (n_features,  1)
  

  In particolare, (X.T @ error) indica il prodotto matriciale tra la matrice trasposta relativa a $\textbf{X}$, indicata con $\textbf{X}^\intercal$, avente dimensione (n_features, n_examples), ed il vettore $\textbf{error}$, avente dimensione (n_examples, 1). Come mostrato nel commento, l’operazione ritorna un vettore di dimensione (n_features, 1).

  La matrice trasposta di una data matrice si ottiene scambiandone le righe con le colonne. In altri termini, la trasposta di una matrice è una nuova matrice in cui le righe diventano colonne e le colonne diventano righe. Quindi, se $\textbf{X}$ ha dimensione (n_examples, n_features), $\textbf{X}^\intercal$ avrà dimensione (n_features, n_examples).

  Nella moltiplicazione di matrici, la trasposizione di una matrice è spesso necessaria per garantire il corretto allineamento delle dimensioni delle due matrici operandi, affinché l’operazione di moltiplicazione sia valida. Ad esempio, se $\textbf{A}$ ha dimensione (m, n) e $\textbf{B}$ ha dimensione (k, n), è possibile trasporre $\textbf{B}$ per ottenere una dimensione (n, k) che consenta di moltiplicare $\textbf{A}$ (m, n) per $\textbf{B}^\intercal$ (n, k), ed ottenere una matrice di dimensione (m, k).

  Sempre nello stesso metodo, l’istruzione:

  grad_b = 2 * np.mean(error)              # float
  

  implementa l’equazione $(16)$. La moltiplicazione vettoriale $\mathbf{1}^{\intercal} \left( \hat{\mathbf{y}} − \mathbf{y} \right)$ viene implementata implicitamente dal metodo numpy.mean, che, per calcolare la media dei valori contenuti nell’array error, deve prima valutare la sommatoria dei suoi elementi.

• La funzione train inizializza l’array $\textbf{w}$ e lo scalare $b$ in questo modo:

  # initialize the weights and bias with values taken from a normal distribution
  rng = np.random.default_rng()
  w = rng.normal(loc=0.0, scale=1.0, size=(X.shape[1], 1))  # w shape: (n_features, 1)
  b = rng.normal(size=None)
  

  Il metodo random.Generator.normal usato per inizializzare $\textbf{w}$ restituisce un array di dimensione (n_examples, 1), dove ogni valore dell’array è estratto causalmente da una distribuzione normale¹ con media $0$ e deviazione standard di $1$. Lo stesso accade per $b$, che viene inizializzato con un singolo valore preso dalla stessa distribuzione normale. In pratica, potevamo evitare queste “raffinatezze” ed inizializzare entrambe le variabili con valori uguali a $0$.

Ricapitolando

In questa sezione abbiamo considerato la regressione lineare multipla, e la sua funzione di ipotesi $(2)$, che permette di usare osservazioni definite da più variabili indipendenti $x_1$, …, $x_n$.

Abbiamo notato che la funzione di ipotesi in uno spazio 2D genera una retta, in uno spazio 3D genera un piano, e in spazi vettoriali a $n+1$ dimensioni genera sottospazi a $n$ dimensioni, chiamati iperpiani. Una retta, o un piano, sono a loro volta iperpiani.

Abbiamo quindi capito che la funzione di ipotesi $(2)$ può essere interpretata come l’equazione di un iperpiano in uno spazio di $n+1$ dimensioni, dove $n$ è il numero di variabili indipendenti usate dall’equazione.

Infine, abbiamo riscritto l’equazione di un iperpiano in formato vettoriale, mediante le equazioni $(10)$ e $(11)$. Questo ci ha permesso di aggiornare il codice del modello di regressione lineare semplice precedentemente definito, per poter supportare osservazioni con feature multiple.