Enseigner à un réseau neuronal l’ensemble de Mandelbrot

Introduction

L’ensemble est l’un des plus beaux objets mathématiques jamais découverts, une fractale si complexe que peu importe à quel point vous zoomez, vous continuez à trouver des détails infinis. Mais et si nous demandions à un réseau de neurones de l’apprendre ?

À première vue, cette question semble étrange. L’ensemble de Mandelbrot est entièrement déterministe, il n’y a ni données, ni bruit, ni règles cachées. Mais cette simplicité en fait un bac à sable parfait pour étudier la manière dont les réseaux de neurones représentent des fonctions complexes.

Dans cet article, nous explorerons comment un simple réseau de neurones peut apprendre à se rapprocher de l’ensemble de Mandelbrot, et comment Caractéristiques de Fourier gaussienne transforme complètement ses performances, transformant des approximations floues en limites fractales nettes.

En cours de route, nous découvrirons pourquoi les perceptrons multicouches (MLP) vanille ont du mal avec les modèles haute fréquence (le problème de biais spectral) et comment les fonctionnalités de Fourier le résolvent.

L’ensemble Mandelbrot

L’ensemble de Mandelbrot est défini sur le plan complexe. Pour chaque nombre complexe \(c\in \mathbb{C}\), on considère la séquence itérative :

$$z_{n+1} = z_n^2 + c, z_0 = 0$$

Si cette suite reste bornée, alors \(c\) appartient à l’ensemble de Mandelbrot.

En pratique, nous approchons cette condition en utilisant la algorithme de temps d’échappement. L’algorithme de temps d’échappement itère la séquence jusqu’à un nombre maximum fixe d’étapes et surveille l’amplitude \(|z_n|\). Si \(|z_n|\) dépasse un rayon d’échappement choisi (généralement 2), la séquence est garantie de diverger et \(c\) est classé comme en dehors de l’ensemble de Mandelbrot. Si la séquence ne s’échappe pas dans le nombre maximum d’itérations, \(c\) est supposé appartenir à l’ensemble de Mandelbrot, le nombre d’itérations étant souvent utilisé à des fins de visualisation.

Transformer l’ensemble de Mandelbrot en un problème d’apprentissage

Pour entraîner notre réseau neuronal, nous avons besoin de deux choses. Il faut d’abord définir un problème d’apprentissagec’est-à-dire ce que le modèle doit prédire et à partir de quelles entrées. Deuxièmement, nous avons besoin données étiquetées: une grande collection de paires d’entrées-sorties tirées de ce problème.

Définir le problème d’apprentissage

À la base, l’ensemble de Mandelbrot définit une fonction sur le plan complexe. Chaque point \(c = x +iy \in \mathbb{C}\) est mappé à un résultat : soit la séquence générée par l’itération reste bornée, soit elle diverge. Cela suggère immédiatement un problème de classification binaire, où l’entrée est un nombre complexe et la sortie indique si le nombre est à l’intérieur de l’ensemble ou non.

Cependant, cette formulation pose des difficultés pour l’apprentissage. La frontière de l’ensemble de Mandelbrot est infiniment complexe et des perturbations arbitrairement petites dans \(c\) peuvent modifier le résultat de la classification. Du point de vue de l’apprentissage, cela se traduit par une fonction cible très discontinue, rendant l’optimisation instable et inefficace en matière de données.

Pour obtenir un objectif d’apprentissage plus fluide et plus informatif, nous reformulons plutôt le problème en utilisant la temps d’évasion informations introduites dans la section précédente. Plutôt que de prédire une étiquette binaire, le modèle est entraîné pour prédire une variable continue dérivée du nombre d’itérations auquel la séquence s’échappe.

Pour produire une fonction cible continue, nous n’utilisons pas directement le nombre d’itérations d’échappement brut. Le nombre d’itérations d’échappement brut est une quantité discrète, son utilisation introduirait des discontinuités, en particulier près de la limite de l’ensemble de Mandelbrot, où de petits changements dans \(c\) peuvent provoquer de grands sauts dans le nombre d’itérations. Pour résoudre ce problème, nous utilisons une valeur de temps d’échappement lisse, qui intègre le nombre d’itérations d’échappement pour produire une cible continue. De plus, nous appliquons également une mise à l’échelle logarithmique qui répartit les premières valeurs d’échappement et compresse les plus grandes, ce qui donne une distribution cible plus équilibrée.

def smooth_escape(x: float, y: float, max_iter: int = 1000) -> float:
    c = complex(x, y)
    z = 0j
    for n in range(max_iter):
        z = z*z + c
        r2 = z.real*z.real + z.imag*z.imag
        if r2 > 4.0:
            r = math.sqrt(r2)
            mu = n + 1 - math.log(math.log(r)) / math.log(2.0)  # smooth
            # log-scale to spread small mu
            v = math.log1p(mu) / math.log1p(max_iter)
            return float(np.clip(v, 0.0, 1.0))
    return 1.0

Avec cette définition, l’ensemble de Mandelbrot devient un problème de régression. Le réseau neuronal est entraîné pour se rapprocher d’une fonction $$f : \mathbb{R}^2 \rightarrow [0,1]$$

mapper les coordonnées spatiales \((x, y)\) dans le plan complexe à une valeur de temps d’échappement lisse.

Stratégie d’échantillonnage des données

Un échantillonnage uniforme du plan complexe serait très inefficace, la plupart des points se trouvant loin de la frontière et portant peu d’informations. Pour résoudre ce problème, l’ensemble de données est orienté vers les régions limites en suréchantillonnant et en filtrant les points dont les valeurs d’échappement se situent dans une bande prédéfinie.

def build_boundary_biased_dataset(
    n_total=800_000,
    frac_boundary=0.7,
    xlim=(-2.4, 1.0),
    res_for_ylim=(3840, 2160),
    ycenter=0.0,
    max_iter=1000,
    band=(0.35, 0.95),
    seed=0,
):
    """
    - Mix of uniform samples + boundary-band samples.
    - 'band' selects points with target in (low, high), which tends to concentrate near boundary.
    """
    rng = np.random.default_rng(seed)
    ylim = compute_ylim_from_x(xlim, res_for_ylim, ycenter=ycenter)

    n_boundary = int(n_total * frac_boundary)
    n_uniform  = n_total - n_boundary

    # Uniform set
    Xu = sample_uniform(n_uniform, xlim, ylim, seed=seed)

    # Boundary pool: oversample, then filter by band
    pool_factor = 20
    pool = sample_uniform(n_boundary * pool_factor, xlim, ylim, seed=seed + 1)

    yp = np.empty((pool.shape[0],), dtype=np.float32)
    for i, (x, y) in enumerate(pool):
        yp[i] = smooth_escape(float(x), float(y), max_iter=max_iter)

    mask = (yp > band[0]) & (yp < band[1])
    Xb = pool[mask]
    yb = yp[mask]

    if len(Xb) < n_boundary:
        # If band too strict, relax it automatically
        keep = min(len(Xb), n_boundary)
        print(f"[warn] Boundary band too strict; got {len(Xb)} boundary points, using {keep}.")
        Xb = Xb[:keep]
        yb = yb[:keep]
        n_boundary = keep
        n_uniform = n_total - n_boundary
        Xu = sample_uniform(n_uniform, xlim, ylim, seed=seed)

    else:
        Xb = Xb[:n_boundary]
        yb = yb[:n_boundary]

    yu = np.empty((Xu.shape[0],), dtype=np.float32)
    for i, (x, y) in enumerate(Xu):
        yu[i] = smooth_escape(float(x), float(y), max_iter=max_iter)

    X = np.concatenate([Xu, Xb], axis=0).astype(np.float32)
    y = np.concatenate([yu, yb], axis=0).astype(np.float32)

    # Shuffle once
    perm = rng.permutation(X.shape[0])
    return X[perm], y[perm], ylim

Modèle de base : un MLP résiduel profond

Notre première tentative utilise un MLP résiduel profond qui prend les coordonnées cartésiennes brutes \((x, y)\) comme entrée et prédit la valeur d’échappement lisse.

# Baseline model
class MLPRes(nn.Module):
    def __init__(
        self,
        hidden_dim=256,
        num_blocks=8,
        act="silu",
        dropout=0.0,
        out_dim=1,
    ):
        super().__init__()
        activation = nn.ReLU if act.lower() == "relu" else nn.SiLU

        self.in_proj = nn.Linear(2 , hidden_dim)
        self.in_act = activation()

        self.blocks = nn.Sequential(*[
            ResidualBlock(hidden_dim, act=act, dropout=dropout)
            for _ in range(num_blocks)
        ])

        self.out_ln = nn.LayerNorm(hidden_dim)
        self.out_act = activation()

        self.out_proj = nn.Linear(hidden_dim, out_dim)

    def forward(self, x):
        x = self.in_proj(x)
        x = self.in_act(x)
        x = self.blocks(x)
        x = self.out_act(self.out_ln(x))
        return self.out_proj(x)

# Residual block
class ResidualBlock(nn.Module):
    def __init__(self, dim: int, act: str = "silu", dropout: float = 0.0):
        super().__init__()
        activation = nn.ReLU if act.lower() == "relu" else nn.SiLU

        # pre-norm-ish (LayerNorm helps a lot for stability with deep residual MLPs)
        self.ln1 = nn.LayerNorm(dim)
        self.fc1 = nn.Linear(dim, dim)

        self.ln2 = nn.LayerNorm(dim)
        self.fc2 = nn.Linear(dim, dim)

        self.act = activation()
        self.drop = nn.Dropout(dropout) if dropout and dropout > 0 else nn.Identity()

        # optional: small init for the last layer to start near-identity
        nn.init.zeros_(self.fc2.weight)
        nn.init.zeros_(self.fc2.bias)

    def forward(self, x):
        h = self.ln1(x)
        h = self.act(self.fc1(h))
        h = self.drop(h)
        h = self.ln2(h)
        h = self.fc2(h)
        return x + h

Ce réseau a une capacité suffisante : profonde, résiduelle et formée sur un vaste ensemble de données avec une optimisation stable.

Résultat

La forme globale de l’ensemble de Mandelbrot est clairement reconnaissable. Cependant, les détails fins proches de la limite sont visiblement flous. Les régions qui devraient présenter une structure fractale complexe semblent trop lisses et les filaments minces sont soit mal définis, soit totalement absents.

Ce n’est pas une question de résolution, de données ou de profondeur. Alors, qu’est-ce qui ne va pas ?

Le problème du biais spectral

Les réseaux de neurones ont un problème bien connu appelé biais spectral:
ils ont tendance à apprendre fonctions basse fréquence d’abord, et j’ai du mal à représenter des fonctions avec des oscillations rapides ou des détails fins.

La limite de Mandelbrot est dominée par des structures très irrégulières et remplies de structures à petite échelle, en particulier près de sa limite. Pour le capturer, le réseau devrait représenter très variations à haute fréquence dans la sortie lorsque \(x\) et \(y\) changent.

Fonctionnalités de Fourier : codage des coordonnées dans l’espace fréquentiel

L’une des solutions les plus élégantes au problème du biais spectral a été introduite en 2020 par Tancik et al. dans leur journal Les fonctionnalités de Fourier permettent aux réseaux d’apprendre des fonctions haute fréquence dans des domaines de faible dimension.

L’idée est de transformer les coordonnées d’entrée avant de les introduire dans le réseau neuronal. Au lieu de donner le \((x, y)\ brut), nous transmettons leur projection sinusoïdale dans aucune direction aléatoire dans un espace de dimension supérieure.

Officiellement:

$$\gamma(x)=[sin(2 \pi Bx),cos(2 \ pi Bx)]$$

où \(B \in \mathbb{R}^{d_{in}​×d_{feat}}​\) est une matrice gaussienne aléatoire.

Cette cartographie agit comme un expansion aléatoire de la base de Fourierpermettant au réseau de représenter plus facilement les détails à haute fréquence.

class GaussianFourierFeatures(nn.Module):
    def __init__(self, in_dim=2, num_feats=256, sigma=5.0):
        super().__init__()
        B = torch.randn(in_dim, num_feats) * sigma
        self.register_buffer("B", B)

    def forward(self, x):
        proj = (2 * torch.pi) * (x @ self.B)
        return torch.cat([torch.sin(proj), torch.cos(proj)], dim=-1)

Caractéristiques de Fourier gaussiennes multi-échelles

Une seule échelle de fréquence pourrait ne pas suffire. L’ensemble de Mandelbrot présente une structure à toutes les résolutions (une caractéristique distinctive de la géométrie fractale).

Pour capturer cela, nous utilisons caractéristiques de Fourier gaussiennes multi-échellescombinant plusieurs bandes de fréquences :

class MultiScaleGaussianFourierFeatures(nn.Module):
    def __init__(self, in_dim=2, num_feats=512, sigmas=(2.0, 6.0, 10.0), seed=0):
        super().__init__()
        # split features across scales
        k = len(sigmas)
        per = [num_feats // k] * k
        per[0] += num_feats - sum(per)

        Bs = []
        g = torch.Generator()
        g.manual_seed(seed)
        for s, m in zip(sigmas, per):
            B = torch.randn(in_dim, m, generator=g) * s
            Bs.append(B)

        self.register_buffer("B", torch.cat(Bs, dim=1))

    def forward(self, x):
        proj = (2 * torch.pi) * (x @ self.B)
        return torch.cat([torch.sin(proj), torch.cos(proj)], dim=-1)

Cela fournit effectivement au réseau une base de fréquence multi-résolution, parfaitement alignée avec la nature auto-similaire des fractales.

Modèle final

Le modèle final a la même architecture que le modèle de base, la seule différence est qu’il utilise le MultiScaleGaussianFourierFeatures.

class MLPFourierRes(nn.Module):
    def __init__(
        self,
        num_feats=256,
        sigma=5.0,
        hidden_dim=256,
        num_blocks=8,
        act="silu",
        dropout=0.0,
        out_dim=1,
    ):
        super().__init__()
        self.ff = MultiScaleGaussianFourierFeatures(
            2,
            num_feats=num_feats,
            sigmas=(2.0, 6.0, sigma),
            seed=0
        )

        self.in_proj = nn.Linear(2 * num_feats, hidden_dim)

        self.blocks = nn.Sequential(*[
            ResidualBlock(hidden_dim, act=act, dropout=dropout)
            for _ in range(num_blocks)
        ])

        self.out_ln = nn.LayerNorm(hidden_dim)
        activation = nn.ReLU if act.lower() == "relu" else nn.SiLU
        self.out_act = activation()

        self.out_proj = nn.Linear(hidden_dim, out_dim)

    def forward(self, x):
        x = self.ff(x)
        x = self.in_proj(x)
        x = self.blocks(x)
        x = self.out_act(self.out_ln(x))
        return self.out_proj(x)

Dynamique de formation

Formation sans fonctionnalités de Fourier

Le modèle apprend lentement la forme générale de l’ensemble de Mandelbrot, puis des plateaux. Une formation supplémentaire ne parvient pas à ajouter plus de détails.

Formation avec les fonctionnalités de Fourier

Ici, les structures grossières apparaissent en premier, suivies par des détails de plus en plus fins. Le modèle continue d’affiner ses prévisions au lieu de stagner.

Résultats finaux

Les deux modèles utilisaient la même architecture, le même ensemble de données et la même procédure de formation. Le réseau est un perceptron multicouche résiduel profond (MLP) formé comme modèle de régression sur la formulation de temps d’échappement lisse.

• Ensemble de données : 1 000 000 d’échantillons du plan complexe, avec 70 % de points concentrés près de la limite fractale grâce à l’échantillonnage biaisé des données.
• Architecture : MLP résiduel avec 20 blocs résiduels et une dimension cachée de 512 unités.
• Fonction d’activation : SiLU
• Formation : 100 époques, taille de lot 4096, optimiseur basé sur Adam, planificateur de recuit cosinus.

La seule différence entre les deux modèles est la représentation des coordonnées d’entrée. Le modèle de base utilise les coordonnées cartésiennes brutes, tandis que le deuxième modèle utilise les caractéristiques de Fourier multi-échelles pour la représentation.

Vue globale

Zoom 1 vue

Zoom 2 Vue

Conclusions

Les fractales telles que l’ensemble de Mandelbrot sont un exemple extrême de fonctions dominées par des structures à haute fréquence. Les approximer directement à partir de coordonnées brutes oblige les réseaux de neurones à synthétiser des oscillations de plus en plus détaillées, une tâche pour laquelle les MLP sont mal adaptés.

Ce que montre cet article, c’est que la limitation ne vient pas de la capacité architecturale, du volume de données ou de l’optimisation. C’est représentation.

En codant les coordonnées d’entrée avec des caractéristiques de Fourier gaussiennes multi-échelles, nous déplaçons une grande partie de la complexité du problème dans l’espace d’entrée. La structure haute fréquence devient explicite, permettant à un réseau neuronal par ailleurs ordinaire de se rapprocher d’une fonction qui serait trop complexe dans sa forme originale.

Cette idée va au-delà des fractales. Les réseaux de neurones basés sur les coordonnées sont utilisés dans l’infographie, l’apprentissage basé sur la physique et le traitement unique. Dans tous ces contextes, le choix du codage d’entrée peut faire la différence entre des approximations fluides et une structure riche et très détaillée.

Remarque sur les ressources visuelles

Toutes les images, animations et vidéos présentées dans cet article ont été générées par l’auteur à partir des sorties des modèles de réseau neuronal décrits ci-dessus. Aucun moteur de rendu fractal externe ni aucun élément visuel tiers n’a été utilisé. Le code complet utilisé pour la formation, le rendu des images et la génération d’animations est disponible dans le référentiel qui l’accompagne.

Code complet

Dépôt Github