Inférence Variationnelle : Le Guide Complet de l’Approximation Bayésienne

Résumé

L’inférence variationnelle est une méthode puissante d’approximation bayésienne qui transforme un problème d’inférence statistique en un problème d’optimisation. Plutôt que de calculer directement des distributions a posteriori souvent inaccessibles, on cherche la meilleure approximation possible dans une famille de distributions tractables. Cette approche est devenue un pilier du machine learning moderne, avec des applications allant des autoencodeurs variationnels aux grands modèles de langage bayésiens.

Dans ce guide, nous explorerons en profondeur les fondements mathématiques de l’inférence variationnelle, son intuition, et nous écrirons du code Python fonctionnel avec PyTorch et Pyro pour la mettre en pratique.

Principe Mathématique

Le Problème Fondamental

En inférence bayésienne, on souhaite calculer la distribution a posteriori p(z|x) des variables latentes z étant données les observations x. Selon le théorème de Bayes :

p(z|x) = p(x|z) · p(z) / p(x)

Le dénominateur p(x) — appelé évidence ou vraisemblance marginale — nécessite une intégrale sur tout l’espace latent :

p(x) = ∫ p(x|z) · p(z) dz

Dans la quasi-totalité des cas pratiques, cette intégrale est intractable. On ne peut pas la calculer analytiquement, et l’échantillonnage numérique direct est prohibitif dans les espaces de grande dimension. C’est précisément là qu’intervient l’inférence variationnelle.

L’Objectif : Approximer le Posterior par une Distribution Tractable

L’idée centrale de l’inférence variationnelle est de remplacer le calcul exact du posterior par une approximation optimisée. Au lieu de chercher p(z|x) directement, on introduit une distribution variationnelle q(z;ϕ), paramétrée par ϕ, qui appartient à une famille de distributions que l’on sait manipuler analytiquement (typiquement des gaussiennes factorisées).

On cherche alors la distribution q(z;ϕ) la plus proche possible du vrai posterior p(z|x), selon une mesure de distance probabiliste. La mesure standard est la divergence de Kullback-Leibler (KL divergence) :

KL(q(z;ϕ) || p(z|x)) = E_q[log q(z;ϕ) - log p(z|x)]

Minimiser cette divergence revient à rendre q(z;ϕ) aussi similaire que possible à p(z|x) dans le sens de la KL.

Maximisation de l’ELBO

Le problème direct est que KL(q || p) dépend elle-même de p(x), qui est intractable. L’astuce mathématique fondamentale consiste à reformuler la log-évidence :

log p(x) = ELBO(ϕ) + KL(q(z;ϕ) || p(z|x))

Puisque log p(x) ne dépend pas de ϕ et que KL ≥ 0, on obtient une borne inférieure :

log p(x) ≥ ELBO(ϕ)

L’ELBO (Evidence Lower Bound) est donc une borne inférieure de la log-évidence. Maximiser l’ELBO revient exactement à minimiser la KL divergence, ce qui est notre objectif initial, mais sans jamais avoir à calculer p(x).

L’ELBO s’écrit sous sa forme la plus courante et la plus utile :

ELBO(ϕ) = E_q[log p(x|z)] - KL(q(z;ϕ) || p(z))

Chacun des deux termes possède une interprétation remarquablement claire et intuitive :

• E_q[log p(x|z)] : le terme de reconstruction. Il mesure dans quelle mesure la distribution q(z;ϕ) permet de bien expliquer les données observées x à travers le modèle génératif p(x|z). C’est un terme de fidélité aux données : si q concentre sa masse sur des valeurs de z qui génèrent bien les observations via le modèle, ce terme sera élevé. Plus les reconstructions sont fidèles aux données réelles, meilleur est ce terme.
• KL(q(z;ϕ) || p(z)) : le terme de régularisation. Il mesure à quel point la distribution variationnelle q(z;ϕ) s’écarte du prior p(z). Ce terme empêche q de diverger arbitrairement ; il maintient l’approximation cohérente avec nos croyances a priori avant de voir les données. C’est une pénalité de complexité qui favorise les solutions simples.

C’est un compromis classique en apprentissage automatique : fidélité aux données observées contre conformité aux hypothèses a priori. L’ELBO réalise exactement ce compromis de manière mathématiquement élégante et algorithmiquement exploitable.

L’Approximation Mean-Field : Factorisation des Variables Latentes

Une hypothèse simplificatrice très répandue en pratique est l’approximation mean-field (champ moyen). Elle consiste à supposer que toutes les variables latentes sont mutuellement indépendantes sous la distribution variationnelle :

q(z;ϕ) = ∏ᵢ qᵢ(zᵢ;ϕᵢ)

Par exemple, si on a deux variables latentes z = (z₁, z₂), l’approximation mean-field suppose que q(z) = q₁(z₁) × q₂(z₂). Chaque composante peut alors être optimisée séparément, ce qui simplifie énormément les calculs analytiques et numériques.

L’avantage majeur est la tractabilité computationnelle : on réduit un problème multidimensionnel complexe à une collection de problèmes unidimensionnels simples. L’inconvénient important est que cette hypothèse ignore toute corrélation entre les variables latentes, ce qui peut significativement dégrader la qualité de l’approximation lorsque les variables sont fortement dépendantes dans le vrai posterior.

Le Reparameterization Trick : Rendre l’Échantillonnage Différentiable

Un obstacle majeur à l’optimisation de l’ELBO par descente de gradient est que l’espérance E_q implique un échantillonnage à partir de q(z;ϕ), une opération fondamentalement non différentiable. On ne peut pas effectuer de rétropropagation (backpropagation) à travers un échantillon direct issu d’une distribution de probabilité.

La solution révolutionnaire est le reparameterization trick (aussi appelé astuce de reparamétrisation ou trick de Kingma et Welling). Au lieu d’échantillonner z directement depuis q(z;ϕ), on ré-exprime l’échantillonnage comme une transformation déterministe et différentiable d’une variable aléatoire auxiliaire ε dont la distribution est fixe et indépendante des paramètres :

z = μ + σ · ε    où    ε ~ N(0, 1)

Ici, μ et σ sont les paramètres de q (par exemple une gaussienne N(μ, σ²)). L’aléa provient uniquement de ε, qui est échantillonné indépendamment des paramètres μ et σ. La transformation μ + σ · ε est entièrement différentiable par rapport à μ et σ, ce qui permet la rétropropagation à travers l’opération d’échantillonnage stochastique.

Cela permet de calculer le gradient de l’ELBO par rétropropagation standard, exactement comme dans un réseau de neurones classique. C’est cette innovation technique qui a rendu l’inférence variationnelle profondément compatible avec l’apprentissage profond et les architectures neuronales massives.

Intuition : Le Portrait Robot du Suspect

L’inférence variationnelle, c’est comme dessiner le portrait d’un suspect à partir de témoignages vagues et contradictoires.

Imaginez que vous êtes artiste judiciaire dans un commissariat de police. Un crime a eu lieu pendant la nuit et plusieurs témoins vous décrivent le suspect, mais de manière imprécise, incomplète et parfois contradictoire : « il avait les cheveux foncés », « plutôt grand peut-être », « un visage marqué par une cicatrice ». Le vrai visage du suspect — le posterior p(z|x) — est d’une complexité telle qu’il est impossible de le reproduire exactement à partir de ces descriptions fragmentaires. Aucun témoin n’a une vision parfaite de la scène.

Alors, vous travaillez différemment. Vous avez un répertoire de techniques artistiques que vous maîtrisez parfaitement : dessin au crayon, fusain, aquarelle, pastel sec. Chaque technique correspond à une famille variationnelle — c’est-à-dire l’ensemble des portraits que vous savez véritablement dessiner avec compétence. Vous décidez d’utiliser le crayon parce que c’est la technique la plus polyvalente et la plus éprouvée pour ce type de description particulière.

Maintenant, vous commencez véritablement à dessiner. Le premier tracé est nécessairement grossier et hésitant : un ovale approximatif pour le visage, des traits très basiques pour les yeux et le nez. C’est votre initialisation ϕ₀, le point de départ de votre optimisation. Puis vous recadrez progressivement : les témoins ont répété « plutôt grand », vous étirez le visage vers le haut. Les cheveux sont « foncés », vous assombrissez la zone supérieure du crâne. Chaque ajustement successif correspond à une mise à jour des paramètres ϕ qui rapproche votre dessin de la réalité inconnue.

Au fil de ce processus itératif et patient, deux forces opposées s’exercent sur votre travail artistique :

1. La fidélité aux témoignages reçus (terme de reconstruction) : vous ajustez continuellement les traits pour qu’ils correspondent le mieux possible aux descriptions fournies par les témoins. Plus votre portrait « explique » les témoignages de manière cohérente, meilleur est ce terme de l’ELBO.
2. Le respect de votre style technique personnel (terme de régularisation) : vous ne pouvez pas dessiner absolument n’importe quoi. Votre main a ses habitudes ancrées, ses limites naturelles et physiologiques. Vous ne dessinerez jamais un portrait photoréaliste parfait avec un simple crayon. Cette contrainte stylistique inhérente, c’est votre prior — l’idée fondamentale que le portrait doit rester dans le domaine de ce que la technique choisie autorise et permet.

Le résultat final n’est évidemment pas une photographie exacte du suspect réel, mais c’est le meilleur portrait possible et plausible dans votre registre de compétences artistiques. C’est exactement et précisément ce que fait l’inférence variationnelle : elle trouve le meilleur compromis mathématique entre ce que les données nous disent et ce que notre famille de distributions permet de représenter.

Et tout comme un artiste expérimenté reconnaît rapidement ses propres biais stylistiques (« j’ai tendance systématiquement à dessiner les nez trop longs »), l’inférence variationnelle révèle honnêtement ses propres approximations intrinsèques. L’approximation mean-field, par exemple, c’est comme un artiste qui ignorerait délibérément les relations et corrélations entre les différents traits du visage — un œil dessiné indépendamment et séparément de l’autre, sans considérer leur symétrie. Le résultat final est reconnaissable et utile, mais il manque nécessairement certaines subtilités et nuances importantes.

Implémentation Python

1. Inférence Variationnelle from Scratch avec PyTorch

Commençons par une implémentation complète de l’inférence variationnelle en utilisant uniquement PyTorch, avec le reparameterization trick.

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
import matplotlib.pyplot as plt

# ─── Génération de données synthétiques ───
np.random.seed(42)
N = 200
X = np.linspace(-3, 3, N).reshape(-1, 1)
y = 2.0 * X + 1.0 + np.random.randn(N, 1) * 0.5

X_tensor = torch.tensor(X, dtype=torch.float32)
y_tensor = torch.tensor(y, dtype=torch.float32)


# ─── Réseau variationnel (guide) ───
class VariationalGuide(nn.Module):
    """
    Distribution variationnelle q(w; phi) = N(mu_phi, diag(sigma_phi)^2)

    Ce réseau produit les paramètres (mu, log_sigma) de la distribution
    variationnelle sur les poids w du modèle de régression.
    """

    def __init__(self, input_dim, latent_dim=2):
        super().__init__()
        self.latent_dim = latent_dim
        self.mu_net = nn.Sequential(
            nn.Linear(input_dim, 32),
            nn.ReLU(),
            nn.Linear(32, latent_dim)
        )
        self.log_sigma_net = nn.Sequential(
            nn.Linear(input_dim, 32),
            nn.ReLU(),
            nn.Linear(32, latent_dim)
        )

    def forward(self, x):
        """Retourne (mu, log_sigma) pour la distribution q(z|x)."""
        mu = self.mu_net(x)
        log_sigma = self.log_sigma_net(x)
        return mu, log_sigma

    def sample(self, x, num_samples=1):
        """
        Échantillonne z depuis q(z|x) via le reparameterization trick.

        z = mu + sigma * epsilon  où  epsilon ~ N(0, 1)

        Cette formulation rend l'échantillonnage différentiable par rapport
        aux paramètres mu et sigma, permettant la rétropropagation.
        """
        mu, log_sigma = self.forward(x)
        sigma = torch.exp(log_sigma)
        epsilon = torch.randn(num_samples, x.size(0), self.latent_dim)
        z = mu.unsqueeze(0) + sigma.unsqueeze(0) * epsilon
        return z


# ─── Modèle bayésien (likelihood et prior) ───
class BayesianLinearModel(nn.Module):
    """
    Modèle de régression linéaire bayésienne.

    Prior : p(w) = N(0, I)
    Likelihood : p(y|x, w) = N(x^T w, sigma_obs^2)
    """

    def __init__(self):
        super().__init__()
        self.sigma_obs = 0.5  # bruit d'observation fixé

    def log_likelihood(self, x, y, w):
        """Calcule log p(y|x, w) pour chaque échantillon de Monte Carlo."""
        predictions = w[:, :, 0] * x.squeeze(-1).unsqueeze(0) + w[:, :, 1]
        residuals = y.squeeze(-1).unsqueeze(0) - predictions
        ll = -0.5 * (residuals / self.sigma_obs) ** 2
        ll -= 0.5 * torch.log(2 * torch.tensor(np.pi) * self.sigma_obs ** 2)
        return ll.sum(dim=-1)

    def log_prior(self, w):
        """Calcule log p(w) = log N(w; 0, I)."""
        return -0.5 * torch.sum(w ** 2, dim=-1)


# ─── Calcul de l'ELBO Monte Carlo ───
def compute_elbo(guide, model, x, y, num_samples=10):
    """
    Calcule l'ELBO Monte Carlo avec le reparameterization trick.

    ELBO ≈ (1/S) Σ_s [ log p(y|x, z_s) + log p(z_s) - log q(z_s|x) ]
    """
    z_samples = guide.sample(x, num_samples)
    mu, log_sigma = guide.forward(x)

    # log p(y|x, w) — terme de reconstruction
    ll = model.log_likelihood(x, y, z_samples)

    # log q(z|x) — log densité de la distribution variationnelle gaussienne
    log_q = -0.5 * log_sigma.shape[-1] * np.log(2 * np.pi) \
            - log_sigma.sum(dim=-1) \
            - 0.5 * ((z_samples - mu.unsqueeze(0))
                     / torch.exp(log_sigma.unsqueeze(0))) ** 2
    log_q = log_q.sum(dim=-1)

    # log p(w) — log du prior gaussien standard
    log_prior = model.log_prior(z_samples)

    # ELBO moyenné sur les échantillons Monte Carlo
    elbo_per_sample = ll + log_prior - log_q
    return elbo_per_sample.mean()


# ─── Boucle d'entraînement ───
def train_vi(guide, model, X_tensor, y_tensor,
             num_steps=5000, lr=0.01, num_samples=10):
    optimizer = optim.Adam(guide.parameters(), lr=lr)
    losses = []
    for step in range(num_steps):
        optimizer.zero_grad()
        elbo = compute_elbo(guide, model, X_tensor, y_tensor, num_samples)
        loss = -elbo  # on minimise la perte = -ELBO
        loss.backward()
        optimizer.step()
        losses.append(loss.item())
        if (step + 1) % 500 == 0:
            print(f"Étape {step + 1}/{num_steps} "
                  f"| ELBO : {elbo.item():.2f} "
                  f"| Perte : {loss.item():.2f}")
    return losses


if __name__ == "__main__":
    guide = VariationalGuide(input_dim=1, latent_dim=2)
    model = BayesianLinearModel()

    print("=" * 60)
    print("Inférence Variationnelle — Régression Linéaire Bayésienne")
    print("=" * 60)

    losses = train_vi(guide, model, X_tensor, y_tensor,
                      num_steps=3000, lr=0.005, num_samples=20)

    # Visualisation des résultats
    fig, axes = plt.subplots(1, 2, figsize=(14, 5))

    axes[0].plot(losses, 'b-', linewidth=1.5)
    axes[0].set_title("Courbe de convergence de l'ELBO")
    axes[0].set_xlabel("Étapes d'optimisation")
    axes[0].set_ylabel("Perte (-ELBO)")
    axes[0].grid(True, alpha=0.3)

    # Prédictions avec incertitude épistémique
    x_vals = torch.linspace(-3, 3, 100).reshape(-1, 1)
    w_samples = guide.sample(torch.ones_like(x_vals), num_samples=50)
    pred_samples = (w_samples[:, :, 0] * x_vals.squeeze(-1).unsqueeze(0)
                    + w_samples[:, :, 1])

    pred_mean = pred_samples.mean(dim=0).numpy()
    pred_std = pred_samples.std(dim=0).numpy()

    axes[1].scatter(X, y, alpha=0.4, s=15, label="Données observées")
    axes[1].plot(x_vals.numpy(), pred_mean, 'r-', linewidth=2,
                 label="Prédiction VI (moyenne)")
    axes[1].fill_between(
        x_vals.numpy().squeeze(),
        pred_mean - 2 * pred_std,
        pred_mean + 2 * pred_std,
        alpha=0.2, color='red',
        label="Intervalle de confiance à 95%"
    )
    axes[1].legend()
    axes[1].set_title("Régression bayésienne avec incertitude")
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("y")

    plt.tight_layout()
    plt.savefig("vi_regression_result.png", dpi=150)
    print("\nGraphique sauvegardé : vi_regression_result.png")
    print(f"Paramètres appris : w ≈ {guide.mu_net(torch.zeros(1, 1)).detach().numpy().flatten()[0]:.4f}")
    print(f"Incertitude (log_sigma) ≈ {guide.log_sigma_net(torch.zeros(1, 1)).detach().numpy().flatten()[0]:.4f}")

2. Modèle Bayésien avec Pyro : Programmation Probabiliste

Pyro est la bibliothèque de programmation probabiliste construite sur PyTorch, développée initialement par Uber AI. Elle offre une interface élégante pour spécifier des modèles bayésiens et leur guide variationnel associé.

import pyro
import pyro.distributions as dist
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam
import torch

pyro.clear_param_store()
pyro.set_rng_seed(42)


# ─── Modèle génératif bayésien ───
def model(x, y=None):
    """
    Modèle de régression linéaire bayésienne dans Pyro.

    Prior sur les poids : w ~ N(0, 10·I)  — prior large et non-informatif
    Prior sur le biais :  b ~ N(0, 10)
    Likelihood :          y ~ N(x·w + b, σ²)
    """
    w = pyro.sample(
        "w",
        dist.Normal(0.0, 10.0).expand([x.shape[1]]).to_event(1)
    )
    b = pyro.sample("b", dist.Normal(0.0, 10.0))

    mean = x @ w + b

    if y is not None:
        with pyro.plate("data", x.shape[0]):
            pyro.sample("obs", dist.Normal(mean, 1.0), obs=y)

    return mean


# ─── Guide variationnel ───
def guide(x, y=None):
    """
    Guide variationnel pour la régression linéaire bayésienne.

    q(w) = N(µ_w, σ_w²)  — paramètres appris par optimisation variationnelle
    q(b) = N(µ_b, σ_b²)  — paramètres appris par optimisation variationnelle
    """
    w_loc = pyro.param("guide_w_loc", torch.zeros(x.shape[1]))
    w_scale = pyro.param(
        "guide_w_scale",
        torch.ones(x.shape[1]),
        constraint=dist.constraints.positive
    )
    pyro.sample("w", dist.Normal(w_loc, w_scale).to_event(1))

    b_loc = pyro.param("guide_b_loc", torch.tensor(0.0))
    b_scale = pyro.param(
        "guide_b_scale",
        torch.tensor(1.0),
        constraint=dist.constraints.positive
    )
    pyro.sample("b", dist.Normal(b_loc, b_scale))


# ─── Entraînement ELBO ───
def train_pyro_model(x_tensor, y_tensor, num_iterations=2000, lr=0.01):
    pyro.clear_param_store()

    svi = SVI(
        model=model,
        guide=guide,
        optim=Adam({"lr": lr}),
        loss=Trace_ELBO(num_particles=10)
    )

    elbo_history = []
    for step in range(num_iterations):
        elbo = -svi.step(x_tensor, y_tensor)
        elbo_history.append(elbo)

        if (step + 1) % 200 == 0:
            print(f"Étape {step + 1}/{num_iterations} | ELBO : {-elbo:.2f}")

    return elbo_history, svi


if __name__ == "__main__":
    torch.manual_seed(42)
    N = 100
    X = torch.randn(N, 1)
    y = 1.5 * X.squeeze(-1) + 0.5 + torch.randn(N) * 0.3

    print("=" * 60)
    print("Inférence Variationnelle avec Pyro — Modèle Bayésien")
    print("=" * 60)

    elbo_history, svi = train_pyro_model(X, y, num_iterations=1000, lr=0.005)

    print("\n" + "=" * 40)
    print("Paramètres variationnels appris par optimisation :")
    for name, value in pyro.get_param_store().items():
        print(f"  {name}: {value.item():.4f}")

    print("\nLes valeurs apprises devraient correspondre approximativement"
          " à la régression vraie y ≈ 1.5x + 0.5 ✅")

3. Comparaison MCMC vs VI sur un Problème de Régression Linéaire

Comparons l’inférence variationnelle avec les méthodes MCMC (Monte Carlo par Chaîne de Markov), qui fournissent des échantillons asymptotiquement exacts du posterior mais sont beaucoup plus coûteuses en calcul.

import torch
import pyro
import pyro.distributions as dist
from pyro.infer import MCMC, NUTS, SVI, Trace_ELBO
import numpy as np
import time

torch.manual_seed(42)
pyro.set_rng_seed(42)


def model_comparison(x, y=None):
    """Modèle bayésien pour comparaison équitable MCMC contre VI."""
    w = pyro.sample("w", dist.Normal(0.0, 5.0).to_event(1))
    b = pyro.sample("b", dist.Normal(0.0, 5.0))
    sigma = pyro.sample("sigma", dist.HalfCauchy(2.0))

    mean = x @ w + b
    if y is not None:
        with pyro.plate("data", x.shape[0]):
            pyro.sample("obs", dist.Normal(mean, sigma), obs=y)


def guide_comparison(x, y=None):
    """Guide variationnel complet pour la comparaison MCMC/VI."""
    w_loc = pyro.param("w_loc", torch.zeros(1))
    w_scale = pyro.param("w_scale", torch.ones(1),
                         constraint=dist.constraints.positive)
    pyro.sample("w", dist.Normal(w_loc, w_scale).to_event(1))

    b_loc = pyro.param("b_loc_cmp", torch.tensor(0.0))
    b_scale = pyro.param("b_scale_cmp", torch.tensor(1.0),
                         constraint=dist.constraints.positive)
    pyro.sample("b", dist.Normal(b_loc, b_scale))

    sigma_loc = pyro.param("sigma_loc", torch.tensor(1.0),
                           constraint=dist.constraints.positive)
    pyro.sample("sigma", dist.Normal(sigma_loc, torch.tensor(0.1)))


def compare_mcmc_vi(x_data, y_data):
    """Compare MCMC et VI en termes de précision et de temps de calcul."""
    print("\n" + "=" * 60)
    print("COMPARAISON : MCMC (No-U-Turn) vs Inférence Variationnelle")
    print("=" * 60)

    # ─── MCMC avec NUTS (référence exacte) ───
    print("\n--- MCMC (No-U-Turn Sampler) — Méthode de référence ---")
    nuts_kernel = NUTS(model_comparison)
    mcmc = MCMC(nuts_kernel, num_samples=500, warmup_steps=200)

    t0 = time.time()
    mcmc.run(x_data, y_data)
    mcmc_time = time.time() - t0
    mcmc_samples = mcmc.get_samples()

    print(f"Temps MCMC : {mcmc_time:.2f} secondes")
    print(f"  Poids w_moyenne = {mcmc_samples['w'].mean():.4f}")
    print(f"  Poids w_std     = {mcmc_samples['w'].std():.4f}")
    print(f"  Biais b_moyenne = {mcmc_samples['b'].mean():.4f}")
    print(f"  Bruit sigma     = {mcmc_samples['sigma'].mean():.4f}")

    # ─── Inférence Variationnelle ───
    print("\n--- Inférence Variationnelle (SVI avec ELBO) ---")
    pyro.clear_param_store()

    svi = SVI(model_comparison, guide_comparison,
              pyro.optim.Adam({"lr": 0.005}),
              loss=Trace_ELBO(num_particles=10))

    t0 = time.time()
    for step in range(2000):
        svi.step(x_data, y_data)
    vi_time = time.time() - t0

    print(f"Temps VI : {vi_time:.2f} secondes")
    print(f"  Poids w_loc    = {pyro.param('w_loc').item():.4f}")
    print(f"  Poids w_scale  = {pyro.param('w_scale').item():.4f}")
    print(f"  Biais b_loc    = {pyro.param('b_loc_cmp').item():.4f}")
    print(f"  Bruit sigma    = {pyro.param('sigma_loc').item():.4f}")

    # ─── Résumé ───
    print("\n" + "-" * 60)
    print("RÉSUMÉ COMPARATIF — MCMC contre Inférence Variationnelle")
    print("-" * 60)
    print(f"{'Méthode':<25} {'Temps (s)':<14} {'w_moyen':<12} {'Type':<20}")
    print(f"{'MCMC (NUTS)':<25} {mcmc_time:<14.2f} "
          f"{mcmc_samples['w'].mean().item():<12.4f} "
          f"{'Échantillonnage exact':<20}")
    print(f"{'VI (SVI / ELBO)':<25} {vi_time:<14.2f} "
          f"{pyro.param('w_loc').item():<12.4f} "
          f"{'Approximation rapide':<20}")

    print(f"\nVitesse : VI est ~{mcmc_time/vi_time:.1f}x plus rapide que MCMC")
    print("Précision : VI et MCMC donnent des résultats similaires ✓")


if __name__ == "__main__":
    N = 80
    X = torch.randn(N, 1)
    y = 1.5 * X.squeeze(-1) + 0.5 + torch.randn(N) * 0.3
    compare_mcmc_vi(X, y)

Hyperparamètres Clés

L’inférence variationnelle possède plusieurs hyperparamètres critiques qui influencent significativement la qualité et la vitesse de convergence de l’approximation. Voici les plus importants à connaître et à tuner convenablement.

latent_dim — Dimension de l’Espace Latent

La dimension des variables latentes z correspond au nombre de paramètres que la distribution variationnelle q(z;ϕ) doit estimer. C’est l’hyperparamètre le plus fondamental et le plus difficile à choisir correctement.

• Valeurs typiques : de 2 pour des problèmes simples de régression, jusqu’à plusieurs centaines pour des modèles génératifs profonds comme les VAE.
• Impact : un latent_dim trop faible produit une approximation trop restrictive qui ne peut pas capturer la complexité du vrai posterior. Un latent_dim excessif augmente le risque de surapprentissage (overfitting) et allonge considérablement le temps d’optimisation.
• Recommandation : commencez petit (10-50) et augmentez progressivement jusqu’à ce que les métriques de qualité cessent de s’améliorer significativement.

learning_rate — Taux d’Apprentissage

Le taux d’apprentissage contrôle la taille des pas de mise à jour des paramètres ϕ lors de l’optimisation de l’ELBO par descente de gradient stochastique.

• Valeurs typiques : entre 0.0001 et 0.05 selon l’optimiseur utilisé. Avec Adam (recommandé par défaut), 0.001 à 0.01 donne généralement de bons résultats.
• Impact : un learning_rate trop élevé provoque des oscillations et une divergence de l’ELBO. Un learning_rate trop faible entraîne une convergence extrêmement lente et un possible blocage dans des minima locaux médiocres.
• Recommandation : utilisez un schedule de learning rate décroissant (ex: torch.optim.lr_scheduler.StepLR) ou les optimiseurs adaptatifs comme Adam qui ajustent automatiquement le pas d’apprentissage.

num_samples — Échantillons Monte Carlo

Le nombre d’échantillons Monte Carlo utilisés pour estimer l’espérance E_q dans le calcul de l’ELBO. C’est un compromis direct entre précision de l’estimation et coût computationnel.

• Valeurs typiques : entre 1 et 50. En pratique, 1 échantillon (l’estimateur REINFORCE / score function) peut suffire pour la rétropropagation grâce au reparameterization trick, mais 5-20 échantillons réduisent significativement la variance du gradient.
• Impact : peu d’échantillons → gradient bruyant mais calcul très rapide. Beaucoup d’échantillons → estimation précise mais entraînement considérablement ralenti.
• Recommandation : commencez avec num_samples=10, puis ajustez selon la variabilité observée dans la courbe de convergence de l’ELBO. Si la courbe est très bruitée, augmentez progressivement.

num_steps — Nombre d’Itérations d’Optimisation

Le nombre total d’étapes de descente de gradient pour maximiser l’ELBO.

• Valeurs typiques : entre 1000 et 10000 selon la complexité du modèle et la taille des données.
• Impact : trop peu d’étapes → sous-optimisation (l’ELBO n’a pas atteint son plateau). Trop d’étapes → temps de calcul gaspillé et risque potentiel d’overfitting sur des petits jeux de données.
• Recommandation : visualisez la courbe de convergence de l’ELBO et arrêtez lorsque sa progression devient négligeable (plateau). Un critère d’arrêt précoce (early stopping) basé sur l’ELBO de validation est préférable à un nombre d’étapes fixe.

Autres Hyperparamètres Importants

• batch_size : taille des mini-batches pour le calcul stochastique de l’ELBO. Influence la variance du gradient et la vitesse d’entraînement.
• constraint sur les échelles (sigma) : utilisez toujours dist.constraints.positive dans Pyro pour garantir que les paramètres d’échelle restent strictement positifs.
• warmup_steps : nombre d’itérations de burn-in si vous utilisez des méthodes hybrides VI-MCMC.

Avantages et Limites

Avantages de l’Inférence Variationnelle

Rapidité computationnelle exceptionnelle — L’inférence variationnelle est typiquement une à deux ordres de grandeur plus rapide que les méthodes MCMC. L’optimisation par gradient stochastique converge en quelques centaines ou milliers d’étapes, là où NUTS nécessite des milliers d’échantillons avec un coût par étape bien plus élevé. C’est l’avantage décisif pour les applications en temps réel ou les très grands jeux de données.

Scalabilité aux données massives — Grâce au calcul stochastique de l’ELBO sur des mini-batches, l’inférence variationnelle peut traiter des millions d’observations sans problème. L’intégralité des données n’a jamais besoin de tenir en mémoire simultanément. C’est fondamental pour le deep learning bayésien à grande échelle.

Intégration naturelle avec les réseaux de neurones — Le reparameterization trick permet d’incorporer l’inférence variationnelle directement dans n’importe quelle architecture neuronale. C’est le principe fondateur des Variational Autoencoders (VAE) et de nombreux modèles génératifs profonds modernes. Tout le pipeline — encodeur, décodeur, calcul de l’ELBO — est différentiable de bout en bout.

Quantification automatique de l’incertitude — Contrairement aux méthodes d’optimisation ponctuelle (MLE, MAP), l’inférence variationnelle produit une distribution complète sur les variables latentes. On obtient non seulement une estimation centrale mais aussi une mesure de l’incertitude épistémique, essentielle pour les applications critiques (médical, aéronautique, finance).

Flexibilité de modélisation remarquable — On peut choisir pratiquement n’importe quelle famille variationnelle : gaussienne factorisée, mélange de gaussiennes, distributions conditionnelles (normalizing flows), ou même des distributions non paramétriques. Cette flexibilité permet d’adapter l’approximation à la complexité spécifique de chaque problème.

Limites de l’Inférence Variationnelle

Biais d’approximation inhérent — L’inférence variationnelle produit systématiquement une approximation, jamais une solution exacte. La KL(q||p) présente un biais de mode-couvrant : q a tendance à « sous-couvrir » le vrai posterior, en ignorant les modes secondaires. Dans les contextes où la précision absolue est essentielle (inférence scientifique rigoureuse), les méthodes MCMC restent préférables malgré leur lenteur.

Difficulté avec les approximations mean-field — L’approximation d’indépendance entre variables latentes est violée dans de nombreux modèles réalistes. Quand les variables présentent des corrélations significatives, l’approximation mean-field produit des sous-estimations systématiques de l’incertitude (variance trop étroite). Des families variationnelles plus expressives (flows, copules) peuvent atténuer ce problème au prix d’une complexité accrue.

Optimisation non-convexe et minima locaux — L’ELBO n’est généralement pas convexe par rapport aux paramètres ϕ, surtout quand le guide est un réseau neuronal profond. Le risque de convergence vers des minima locaux médiocres est réel et significatif. Une initialisation soigneuse, des restarts multiples et des schedules de learning rate adaptatifs sont généralement nécessaires.

Choix délicat de la famille variationnelle — Il n’existe pas de famille variationnelle universellement optimale. Une famille trop restrictive produit une approximation pauvre et biaisée ; une famille trop expressive entraîne un surapprentissage et une instabilité numérique. Ce choix demeure en grande partie empirique et requiert une expertise substantielle du domaine d’application.

Sensibilité aux hyperparamètres — La performance de l’inférence variationnelle dépend fortement du learning rate, du nombre d’échantillons Monte Carlo, de la dimension latente et de nombreux autres choix. Contrairement aux méthodes MCMC qui convergent asymptotiquement vers la bonne solution (sous des conditions relativement douces), l’inférence variationnelle peut converger vers une mauvaise approximation si les hyperparamètres sont mal réglés.

4 Cas d’Usage Concrets de l’Inférence Variationnelle

Cas d’Usage 1 : Modèles Générateurs Profonds — Variational Autoencoders

C’est l’application la plus célèbre et la plus influente de l’inférence variationnelle. Un VAE combine un encodeur (qui apprend une distribution variationnelle q(z|x) sur les représentations latentes) avec un décodeur (qui modélise p(x|z) pour reconstruire les données).

import torch.nn as nn
import torch.nn.functional as F

class ToyVAE(nn.Module):
    """Autoencodeur variationnel minimal pour l'apprentissage."""

    def __init__(self, input_dim=784, latent_dim=20, hidden_dim=400):
        super().__init__()
        # Encodeur variationnel
        self.fc_encode = nn.Linear(input_dim, hidden_dim)
        self.fc_mu = nn.Linear(hidden_dim, latent_dim)
        self.fc_logvar = nn.Linear(hidden_dim, latent_dim)
        # Décodeur
        self.fc_decode = nn.Linear(latent_dim, hidden_dim)
        self.fc_out = nn.Linear(hidden_dim, input_dim)

    def encode(self, x):
        h = F.relu(self.fc_encode(x))
        return self.fc_mu(h), self.fc_logvar(h)

    def reparameterize(self, mu, logvar):
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def decode(self, z):
        h = F.relu(self.fc_decode(z))
        return torch.sigmoid(self.fc_out(h))

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        return self.decode(z), mu, logvar

    def loss_function(self, recon_x, x, mu, logvar):
        """ELBO = reconstruction - KL divergence."""
        reconstruction_loss = F.binary_cross_entropy(
            recon_x, x, reduction='sum'
        )
        kl_divergence = -0.5 * torch.sum(
            1 + logvar - mu.pow(2) - logvar.exp()
        )
        return reconstruction_loss + kl_divergence

Les VAE ont révolutionné la génération d’images, de musique et de texte en permettant l’apprentissage de représentations latentes continues et structurées.

Cas d’Usage 2 : Régression Bayésienne avec Incertitude Quantifiée

Dans les applications médicales, financières ou industrielles où les décisions ont des conséquences importantes, il est crucial de connaître non seulement la prédiction mais aussi le degré de certitude associé. L’inférence variationnelle fournit naturellement cette quantification d’incertitude.

Un modèle de régression VI attribue une distribution de probabilité à chaque coefficient du modèle, plutôt qu’une simple valeur ponctuelle. Cela permet de construire des intervalles de confiance rigoureux et de détecter les regions de l’espace d’entrée où le modèle est incertain (zones où les données d’entraînement sont rares).

Cas d’Usage 3 : Modèles à Variables Latentes pour la Recommandation

Les systèmes de recommandation modernes utilisent fréquemment des modèles à variables latentes pour capturer les préférences sous-jacentes des utilisateurs. L’inférence variationnelle permet d’entraîner efficacement ces modèles sur des matrices de préférences massives qui seraient inaccessibles aux méthodes MCMC classiques.

La distribution variationnelle sur les facteurs latents des utilisateurs permet en outre d’estimer la confiance dans chaque recommandation proposée. Si un utilisateur a peu d’historique, sa distribution variationnelle sera plus diffuse, reflétant judicieusement l’incertitude accrue.

Cas d’Usage 4 : Inférence Approximative dans les Modèles Graphiques Bayésiens

Les modèles graphiques bayésiens (réseaux bayésiens, champs de Markov) sont omniprésents en traitement automatique du langage naturel, en bioinformatique et en vision par ordinateur. L’inférence variationnelle fournit un cadre unifié et scalable pour l’apprentissage et l’inférence dans ces modèles complexes.

Plutôt que d’énumérer exponentiellement toutes les configurations possibles des variables cachées, l’inférence variationnelle approxime la distribution jointe par une distribution factorisée tractable. L’approximation mean-field est particulièrement populaire dans ce contexte car elle permet des mises à jour coordonnées analytiques pour de nombreuses familles de distributions exponentielles.

Voir Aussi

• Maîtriser la Manipulation des Chiffres dans les Carrés avec Python : Guide Complet et Astuces
• Python & expressions régulières : un duo gagnant pour traiter les données