Aller au contenu

2- Interpolation et optimisation

Pour installer uv : 

curl -LsSf https://astral.sh/uv/install.sh | sh

Pour créer un environnement virtuel, par exemple avec python 3.11: 

uv venv --python 3.11  /tmp/venv
source /tmp/venv/bin/activate

Pour installer les dépendances de ce TP

uv pip install PyQt6 matplotlib numpy scikit-learn

Notations

Dans ce qui suit, on note :

\[ \begin{array} \forall x \in \mathbb{R}^N, \|x\|_2 &= \sqrt{\sum_{i=0}^{N-1} x_i^2}\\ \forall N \in \mathbb{N}^*, \mathbf{1}_N &= \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} \end{array} \]

Interpolation

On se donne un ensemble de \(N\) points \((x_i, y_i) \in \mathbb{R}^n \times \mathbb{R}\). On suppose que ces points sont issus d'un processus que nous ne connaissons pas. Appelons cet oracle \(f^\star: \mathbb{R}^n \mapsto \mathbb{R}\). On aimerait estimer les observations en des points \(x \in \mathbb{R}^n\) (à priori différents des points \(x_i\), sinon, cela devient moins intéressant), \(y_j = f^\star(x_j)\). Notons \(\hat{y}_j\) notre estimation.

Par la suite, on notera toujours \((x_i, y_i)_{i \in [0, N-1]}\) l'ensemble des observations à l'aide desquelles nous devons fournir une estimation de l'observation en un nouveau point \(x\). Nous noterons cette estimation \(\hat{y}(x)\).

Il existe énormément de techniques permettent d'adresser ce problème d'estimation : le plus proche voisin, le voisin naturel, l'inverse distance weighting, les splines, le krigging, l'interpolation bilinéaire, bicubique, etc.. Nous allons en voir quelques unes et elles seront l'occasion de vous perfectionner en python !

Exemples de fonctions à interpoler

Vous disposez du script function.py qui offre des implémentations de plusieurs fonctions à interpoler : un sinus cardinal (1D), la norme (2D) et la fonction de Rastrigin (2D). Les classes Sinc, Norm et Rastrigin peuvent être échantillonnées et tracées sur un certain domaine :

import function

# Sinus cardinal
sampler = function.Sinc()
xi, yi = sampler.sample(10)
sampler.plot(200)

# Norme
sampler = function.Norm()
xi, yi = sampler.sample(10)
sampler.plot(50)

# Rastrigin
sampler = function.Rastrigin()
xi, yi = sampler.sample(10)
sampler.plot(100)

Sinus (1D) Norme (2D) Fonction de rastrigin (2D)

Ces fonctions vont servir de fonctions de test de différents algorithmes d'interpolation.

Plus proche voisin

Pour calculer la valeur de notre signal en un point \(x\), une première approche consiste à retourner l'observation associée au point \(x_i\) le plus proche.

Une première approche relativement simple consiste à retourner l'observation \(y_{i^\star}\) du point \(x_{i^\star}\) le plus proche de \(x\) :

\[ \hat{y}_{nn}(x) = y_{i^\star}, i^\star = argmin_{i\in[0, N-1]} \|x_i - x\|_2 \]

Vous disposez du script nn.py qui offre une base d'implémentation, à compléter, de l'algorithme d'interpolation du plus proche voisin.

Question

Vous devez :

  • compléter l'implémentation des fonctions fit et __call__ qui se chargent respectivement de :

    • l'apprentissage de la fonction, i.e. la détermination de ses paramètres à partir des observations \(x_i, y_i\)
    • le calcul de la valeur en des nouveaux points \(x\)

  • Tester, en 1D et 2D : vous disposez déjà de deux fonctions pour réaliser un test en 1D

  • Quelles sont les valeurs asymptotiques de l'interpolateur ?

En bonus, vous pouvez évaluer l'erreur d'interpolation sur un ensemble indépendant de points.

Un exemple d'interpolation du sinus cardinal par un plus proche voisin est représenté ci-dessous :

Approximation d'un sinus cardinal par le plus proche voisin

Fonctions radiales de base

Une deuxième approche définit le problème sous l'angle de la régression linéaire avec un modèle qui est une combinaison linéaire de fonctions radiales de base. En termes Mathématiques, cela s'exprime de la manière suivante.

On considère une famille de fonctions :

\[ g_w(x) = \beta + \sum_{k=0}^{K-1} \theta_k \phi_k(x) \]
\[ \forall k \in [0, K-1], \phi_k(x) = e^{\frac{\|x - \mu_k\|^2}{2 \sigma_k^2}} \]

La fonction \(\phi(x)\) est la fonction dont les composantes sont les fonctions \(\phi_k\). Les poids \(w_k\) sont des réels. Les paramètres \(\mu_k \in \mathbb{R}^n, \sigma_k \in \mathbb{R}\) sont les moyennes et variance des fonctions \(\phi_k\).

Pour trouver les paramètres \(w_k\), on va se poser une fonction de perte à minimiser :

\[ J(w) = argmin_w \sum_{i=0}^{N-1} (y_i - g_w(x_i))^2 \]

Cette fonction pourrait dépendre également des centres \(\mu_k\) et variances \(\sigma_k\) mais on va considérer dans cette approche que ces paramètres seront fixés.

Cette fonction de perte est dite des "moindres carrés". Notre problème d'optimisation peut se réécrire sous la forme d'un système linéaire. Si on se note \(\Phi\) la matrice dont la ligne \(i\) est égale au vecteur \(\phi(x_i)\) et \(Y\) la matrice remplie des observations à prédire \(y_i\), notre problème d'optimisation se ré-exprime sous la forme matricielle :

\[ J(w) = argmin_w \| Y - (\Phi \theta + \beta \mathbf{1}_N) \|^2 \]

avec \(\mathbf{1}_N\) un vecteur de taille \(N\) rempli de \(1\).

Le système linéaire \(\Phi w = Y\) n'admet pas forcément de solution mais le minimum de \(J(w)\) existe. Il peut également exister plusieurs minimums. On peut rendre le problème d'optimisation bien conditionné (un seul minimum) tout en lissant le résultat en ajoutant une pénalité au problème d'optimisation, par exemple sous la forme d'un terme qui dépends de la norme 1 du vecteur de paramètre et minimiser

\[ J(w) = argmin_w \| Y - (\Phi \theta + \beta \mathbf{1}_N) \|^2 + \lambda \|\theta\|_1 \]

avec \(\|\theta\|_1 = \sum_{k=0}^{K-1} |\theta_k|\).

Ceci défini le problème des moindres carrés régularisés sous pénalité L1, qu'on appelle également le problème LASSO.

Pour résoudre le problème LASSO, nous allons utiliser l'implémentation de scikit learn. Pour l'utiliser, en s'inspirant de la documentation : 

import numpy as np
from sklearn import linear_model

clf = linear_model.Lasso(alpha=0.0001)

x = np.array([[0,0], [1, 1], [2, 2]])
y = np.array([0, 1, 2])

clf.fit(x, y)

theta = clf.coef_
beta = clf.intercept_

new_x = np.array([[1.5, 1.5]])
predicted_y = np.dot(np.dot(new_x, theta) + beta

Pour les centres \(\mu_k\) et variances \(\sigma_k\) nous allons :

  1. utiliser autant de fonctions de base qu'il existe de points \((x_i, y_i)\), i.e. \(K = N\)
  2. fixer les variances égales à une valeur à déterminer : \(\forall k \in [1..K], \sigma_k = \sigma_0\)

Vous disposez du script rbf.py qui offre une base d'implémentation, à compléter, de l'algorithme d'interpolation RBF.

Question

Vous devez :

  • compléter l'implémentation des fonctions fit et __call__ qui se chargent respectivement de :

    • l'apprentissage de la fonction, i.e. la détermination de ses paramètres à partir des observations \(x_i, y_i\)
    • le calcul de la valeur en des nouveaux points \(x\)

  • Tester, en 1D et 2D : vous disposez déjà de deux fonctions pour réaliser un test en 1D

  • Quelles sont les valeurs asymptotiques de l'interpolateur ?

En bonus, vous pouvez évaluer l'erreur d'interpolation sur un ensemble indépendant de points.

Un exemple d'interpolation du sinus cardinal par un réseau RBF est représenté ci-dessous :

Approximation d'un sinus cardinal par un réseau RBF

Inverse Distance Weighting

La dernière approche que nous allons considérer a été développée dans les années 1960 par Donald Shepard et porte le nom de Inverse Distance Weighting. Cette approche va calculer la valeur à attribuer au point \(x\) en pondérant les valeurs \(y_i\) par des poids qui sont inversement proportionnels à la distance qui les séparent de \(x\).

Mathématiquement, cela s'écrit par :

\[ \hat{y}(x) = \begin{cases} \frac{\sum_{i=0}^{N-1} w_i(x) y_i}{\sum_{i=0}^{N-1} w_i} & \text{ si } \forall i, \|x - x_i\|_2 \neq 0\\ y_{i^\star} & \text{ si } \exists i^\star \text{ tel que } \|x - x_{i^\star}\|_2 = 0 \end{cases} \]

avec

\[ \forall i \in [0, N-1], \forall x \neq x_i, w_i(x) = \frac{1}{\|x-x_i\|_2^{p}} \]

et \(p = 2\) par exemple. Vous disposez du script idw.py qui offre une base d'implémentation, à compléter, de l'algorithme d'interpolation IDW. L'implémentation de cet algorithme suit le même format que les implémentations précédentes. Il vous suffit de compléter les fonctions fit et __call__ de la classe IDW.

Question

Vous devez :

  • compléter l'implémentation des fonctions fit et __call__ qui se chargent respectivement de :

    • l'apprentissage de la fonction, i.e. la détermination de ses paramètres à partir des observations \(x_i, y_i\)
    • le calcul de la valeur en des nouveaux points \(x\)

  • Tester, en 1D et 2D : vous disposez déjà de deux fonctions pour réaliser un test en 1D

  • Quelles sont les valeurs asymptotiques de l'interpolateur ?

En bonus, vous pouvez évaluer l'erreur d'interpolation sur un ensemble indépendant de points.

Un exemple d'interpolation du sinus cardinal par un interpolateur IDW est représenté ci-dessous :

Approximation d'un sinus cardinal par une approche IDW

Cas réel : Interpolation des polluants à Paris

Un certain nombre de mesures, relatifs aux polluants, sont réaliées en France. Geod'Air centralise les données collectées partout en France. Vous pouvez accéder aux données depuis https://www.geodair.fr/donnees/consultation :

  • en moyennes horaires, moyennes journalières, etc..
  • des particules fines (PM10, PM2.5, monoxyde de carbone, benzène, dioxyde d'azote, etc..)

Ci-dessous, une représentation de l'interface et quelques indications pour que vous puissiez récupérer les données. Faites attention aux options qui permettent de restreindre l'étendue spatiale et

Interface Geod'air

Pour cette partie du TP, il convient d'ajouter quelques packages à votre environnement virtuel :

uv pip install pandas cartopy contextily

Vous disposez du script polluant.py qui offre une base d'implémentation, dans laquelle il faut ajouter les fonctions d'interpolation.

Question

Vous devez :

  • récupérer des données depuis le site de geod'air
  • compléter le code
  • expérimenter et comparer les différentes approches

Par exemple, ci-dessous sont illustrées les différentes approches d'interpolation du PM10 pour le 10 juin 2025.

Approximation avec un plus proche voisin Approximation avec un RBF Approximation avec un IDW