Gradient descent - Machine Learning Cursussen

We zagen tot nu toe dat we parameters van een lineair regressiemodel kunnen schatten via:

Grid search: We zoeken exhaustief een deel van de parameterruimte af naar optimale waarden voor de SSE loss. Zeker bij complexere modellen ( $\pmb{b} = \begin{bmatrix}b_1, b_2, \ldots, b_n\end{bmatrix}$ ) is dit praktisch onhaalbaar.
Monte Carlo sampling: Hier leveren we ons over aan het willekeurig aftasten van de $SSE$ functie. Dit levert zelfs bij simpele modellen geen stabiele resultaten op.
Ordinary least squares (OLS): Deze puur analytische oplossing kan in één keer een optimale oplossing opleveren (in de least squares zin), maar kan onstabiele resultaten opleveren wanneer er sprake is van multicollineariteit.

In deze sectie bespreken we een populair alternatief: gradient descent. Zeker in de context van complexere ML modellen (bv. neurale netwerken) biedt deze oplossing (en aanverwanten) een stabieler alternatief. Zoals de naam zelf suggereert, maakt de gradient descent oplossing, net zoals OLS gebruik van de gradiënt. Het is echter net zoals bij Monte Carlo sampling een iteratief, en dus geen closed form, leeralgoritme.

Gradiënt

Herinner je dat de gradiënt in deze context de tensor is met alle partiële afgeleiden van de loss functie bij een gegeven set van parameterwaarden $\nabla \mathcal{L}(\pmb{b}^k)$ . Ieder individueel element staat voor de richting van de raaklijn met de loss (hyper) surface voor een bepaalde parameter ( $b^k_i$ ), in de veronderstelling dat alle andere parameters constant blijven (de partiële afgeleide: $\frac{\partial \mathcal{L}(\pmb{b}^k)}{\partial b^k_i}$ ).

De grootte van de individuele elementen geeft aan hoe sterk de loss functie toeneemt (bij een positieve waarde) of afneemt (bij een negatieve waarde) als de overeenkomstige parameterwaarde ( $b^k_i$ ) met een bepaalde stapgrootte wordt verschoven. Let wel dat het de toe-/afname betreft volgens de lokale lineaire benadering (vandaar geen $=$ maar $\sim$ hieronder):

\mathcal{L}(b^{k+1}_i) \sim \mathcal{L}(b^k_i) + \frac{\partial \mathcal{L}(\pmb{b}^k)}{\partial b^k_i} \Delta b^{k+1}_i

(1)

met

\Delta b^{k+1}_i = b^{k+1}_i - b^k_i

(2)

Bij gradient descent is het basisidee dat we, zoals bij de Monte Carlo benadering, iteratief parameterwaarden aanpassen in de richting waarin de loss daalt. Alleen laten we dit nu niet afhangen van een random zoektocht, maar gebruiken we de gradiënt om ons te leiden. We beginnen met een initiële schatting $\pmb{b}_0$ . Daarna gaan we die schatting iteratief updaten met de volgende regel:

\pmb{b}^{k+1} = \pmb{b}^k - \lambda \nabla \mathcal{L}(\pmb{b}^k)

(3)

waarbij de hyper parameter $\lambda$ de learning rate is met $\lambda > 0$

Met andere woorden:

bij iedere update, verschuiven we alle parameterwaarden in de richting van de sterkste daling van de loss functie
de grootte van de update wordt bepaald als een fractie van de gradiënt (of de individuele partiële afgeleiden), gecontroleerd door de learning rate $\lambda$

Het algoritme maakt ten slotte ook gebruik van een stopping rule om te bepalen wanneer de winst qua daling in the loss klein genoeg is geworden om de iteraties stop te zetten.

Note

Om de formule beter te begrijpen kunnen we ons een situatie voorstellen waarbij de learning rate 1 is en:

de partiële afgeleide 1 is: de raaklijn stijgt met een slope 1 van links naar rechts op de schaal van onze parameter $b_i$
de partiële afgeleide -1 is: de raaklijn daalt met een slope -1 van links naar rechts op de schaal van onze parameter $b_i$

Om onze parameter $b_i$ aan te passen in de richting van een lagere loss waarde moeten we die dan bij (1.) naar links verschuiven en bij (2.) naar rechts. Als we de formule toepassen krijgen we inderdaad bij (1.) een verschuiving van -1 en bij (2.) van +1.

import numpy as np
import plotly.express as px
import plotly.graph_objects as go

from ml_courses.sim.linear_regression_sse_viz import LinearRegressionSSEVisualizer
from ml_courses.sim.monte_carlo_tips import MonteCarloTipsSimulation

# Function to calculate gradient of SSE loss
def calculate_gradients(b1, b2, x, y):
    """
    Calculate gradient of SSE loss function.

    Parameters
    ----------
    - b1: base tip parameter
    - b2: tip rate parameter
    - x: order totals
    - y: observed tips

    Returns
    -------
    - der_b1: partial derivative w.r.t. b1
    - der_b2: partial derivative w.r.t. b2
    """
    y_pred = b1 + b2 * x
    residuals = y - y_pred

    # ∂SSE/∂b1
    der_b1 = -2 * np.sum(residuals)

    # ∂SSE/∂b2
    der_b2 = -2 * np.sum(x * residuals)

    return der_b1, der_b2

# Generate data
sim = MonteCarloTipsSimulation()

# Standardized variables
order_totals = (sim.order_totals - np.mean(sim.order_totals)) / np.std(sim.order_totals)
observed_tips = (sim.observed_tips - np.mean(sim.observed_tips)) / np.std(sim.observed_tips)

# Random number generator for reproducibility (starting guesses)
rng = np.random.default_rng(123)

# Initial guesses
b1 = rng.uniform(-2.0, 2.0)
b2 = rng.uniform(-1.0, 3.0)

# Initial SSE loss
init_loss = sim.calculate_loss(b1, b2, order_totals, observed_tips)
print(f"Initial SSE loss: {init_loss:.2f} (with guess: {b1:.2f} + {b2:.2f} × order)")

# Gradient descent parameters
learning_rate = 1e-5
tolerance = 1e-8
n_iterations = 100000
b1_history = [b1]
b2_history = [b2]
loss_history = [init_loss]

for i in range(n_iterations):
    # Calculate gradients
    der_b1, der_b2 = calculate_gradients(b1, b2, order_totals, observed_tips)

    # Update parameters
    b1 = b1 - learning_rate * der_b1
    b2 = b2 - learning_rate * der_b2

    # Calculate new loss
    loss = sim.calculate_loss(b1, b2, order_totals, observed_tips)

    b1_history.append(b1)
    b2_history.append(b2)
    loss_history.append(loss)

    # Check convergence
    if abs(loss - loss_history[-2]) < tolerance:
        print("Converged.")
        break

print(f"Performed {i + 1} iterations.")
print(f"Estimated parameters (Standardized): b1 = ${b1_history[-1]:.2f}, b2 = {b2_history[-1]:.2f}")
print(f"Final SSE loss: {loss_history[-1]:.2f}")

Initial SSE loss: 183.42 (with guess: 0.73 + -0.78 × order)

Converged.
Performed 8697 iterations.
Estimated parameters (Standardized): b1 = $0.00, b2 = 0.97
Final SSE loss: 3.07

viz = LinearRegressionSSEVisualizer(
    x_data=sim.order_totals,
    y_data=sim.observed_tips,
    true_bias=sim.true_b1,
    true_slope=sim.true_b2,
    standardize=True,
)

fig = viz.create_3d_surface_plot(
    bias_samples=np.array(b1_history),
    slope_samples=np.array(b2_history),
    loss_samples=np.array(loss_history),
    resolution=40,
    scale_factor=2.0,
)
fig.show()

px.line(
    x=np.arange(len(b1_history)),
    y=b1_history,
    labels={"x": "Iteration", "y": "Base Tip (b1)"},
    title="Base Tip (b1) During Gradient Descent",
    height=400,
    width=700,
).show()

px.line(
    x=np.arange(len(b2_history)),
    y=b2_history,
    labels={"x": "Iteration", "y": "Tip rate (b2)"},
    title="Tip rate (b2) During Gradient Descent",
    height=400,
    width=700,
).show()

Voordelen¶

Gradient descent kan oplossingen vinden voor complexe modellen - in het geval van lineaire regressie, ook daar waar OLS last heeft van multicollineariteit.
Het maakt rechtstreeks gebruik van de lokale eigenschappen van het loss oppervlak, waardoor het veel efficiënter en stabieler is dan gewone blinde sampling.

Nadelen¶

De startpositie (initiële parameterschattingen) kan een grote invloed hebben op de snelheid van convergentie.
Learning rate en convergentieregel moeten doordacht gekozen worden:
- Een te grote learning rate kan ervoor zorgen dat er overshooting is van de target (loss gradiënt $\pmb{0}$ ).
- Een stop-regel die niet streng genoeg is, leidt automatisch tot een inaccuraat resultaat.
Kan traag zijn/veel geheugen vragen bij grote datasets.
Het algoritme kan vast komen te zitten in lokale minima.

Lokale minima¶

Bij complexere modellen bestaat het risico op hyper loss surfaces die locaties bevatten waar de gradiënt 0 is, met rondom hogere loss waarden, maar waar die locaties niet overeenkomen met het echte minimum van de loss functie. Het gradient descent algoritme kan hierin terechtkomen en zijn stop-condities bereiken zonder te weten dat het zich in een lokaal minimum bevindt.

def multi_minima_function(x, y):
    """
    Return example function with multiple local minima.

    f(x, y) = sin(x) * cos(y) + 0.1 * (x^2 + y^2).
    """
    return np.sin(x) * np.cos(y) + 0.1 * (x**2 + y**2)


# Generate grid
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
Z = multi_minima_function(X, Y)

# Plotly 3D surface plot
fig = go.Figure(data=[go.Surface(z=Z, x=X, y=Y, colorscale="Viridis")])
fig.update_layout(
    title="Function with Multiple Local Minima",
    scene={"xaxis_title": "x", "yaxis_title": "y", "zaxis_title": "f(x, y)"},
    height=600,
    width=800,
)
fig.show()

Remedies¶

Door vanuit meerdere (al dan niet random) initiële schattingen te vertrekken, kunnen we nagaan of we bij dezelfde oplossing terechtkomen en, indien niet, de oplossing met de kleinste loss selecteren.
Stochastische gradient descent

Stochastische gradient descent¶

Bij stochastische gradient descent (SGD) passen we de gradient descent methode aan door niet de volledige dataset te gebruiken bij elke iteratie, maar slechts een sample (soms zelfs slechts één observatie). Dit betekent dat we de gradiënt niet berekenen op basis van alle data punten, maar op basis van een (kleine) willekeurig gekozen subset.

De update regel wordt dan:

\pmb{b}^{k+1} = \pmb{b}^k - \lambda \nabla \mathcal{L}(\pmb{b}^k, \mathcal{B}^k)

(4)

waarbij $\mathcal{B}^k$ een mini-batch is: een willekeurige subset van de trainingsdata bij iteratie $k$ .

Voordelen van SGD¶

Sneller bij grote datasets: Door niet alle observaties te gebruiken bij elke update, is SGD veel sneller dan reguliere gradient descent, vooral bij grote datasets.
Helpt lokale minima te vermijden: De ruis die wordt geïntroduceerd door het gebruik van data samples kan het algoritme helpen om uit lokale minima te ontsnappen. De gradiënt is niet meer deterministisch; door de shuffling kan er mogelijks een beter pad naar het globale minimum gevonden worden.
Minder geheugen vereist: Omdat niet alle data tegelijk in het geheugen geladen hoeft te worden, is SGD geschikt voor zeer grote datasets.
Online learning mogelijk: SGD kan incrementeel leren van nieuwe data zonder de volledige dataset opnieuw te moeten verwerken.