Problème de Kaplan sous Forme de Moindres Carrés
Le problème de Kaplan peut être formulé comme un problème d’optimisation aux moindres carrés. L’objectif est d’estimer les paramètres et la condition initiale d’une équation différentielle linéaire vectorielle en minimisant l’écart entre les mesures observées et les prédictions du modèle.
Formulation Mathématique
Considérons le problème de Cauchy linéaire décrit par :
x'(t) = Ax(t) + Bu(t), \quad x(0) = x_0
où :
- x(t) \in \mathbb{R}^n est le vecteur d’état à l’instant t.
- A \in \mathcal{M}_{n}(\mathbb{R}) est la matrice des coefficients.
- B \in \mathcal{M}_{n, m}(\mathbb{R}) est la matrice d’entrée.
- u(t) \in \mathbb{R}^m est le vecteur des entrées.
Nous disposons de mesures x_i à des instants t_i pour i = 1, 2, \ldots, N. Nous cherchons à estimer les paramètres
\theta = (A, B, x_0) en minimisant la fonction coût suivante :
F(\theta) = \frac{1}{2} \sum_{i=1}^{N} \| x_i - x(t_i, \theta) \|^2
où :
- x(t_i, \theta) est la prédiction du modèle, c’est-à-dire la solution du problème de Cauchy, à l’instant t_i avec les paramètres \theta.
- \| \cdot \| désigne la norme euclidienne.
Les matrices A et B et le vecteur x_0 sont ici considérés comme des paramètres à estimer. Dans certains cas, seuls certains paramètres peuvent être inconnus. Il est aussi possible de ne chercher qu’à estimer une partie des paramètres. Dans ce cas, si nous notons \theta \in \mathbb{R}^l les paramètres à estimer, nous pouvons écrire la dépendance de A, B et x_0 en fonction de \theta, et le problème de Cauchy devient :
x'(t) = A(\theta)\, x(t) + B(\theta)\, u(t), \quad x(0) = x_0(\theta).
Nous parlons de problème aux moindres carrés linéaire si \theta \mapsto x(t_i, \theta) est linéaire pour toute donnée t_i, sinon on parle de problème aux moindres carrés non linéaire.
Algorithme de Gauss-Newton
Description de la méthode
L’algorithme de Gauss-Newton est une méthode itérative utilisée pour minimiser la fonction coût F. Cette algorithme est utilisé pour la résolution de problèmes aux moindres carrés non linéaires. Si le problème est linéaire, seul la première itération suffit. A partir d’un itéré initial \theta_0, l’algorithme va construire une suite d’itérés (\theta_k)_{k\ge 0} qui on l’espère converge vers un minimum de F.
L’idée principale de l’alogrithme de Gauss-Newton est de résoudre à chaque itération un problème aux moindres carrés linéaire. Notons \theta_k l’itéré courant. Nous pouvons approcher pour chaque instant t_i le modèle :
\theta \mapsto x(t_i, \theta)
par son approximation linéaire au voisinage du point courant \theta_k :
\Delta\theta \mapsto x(t_i, \theta_k) + \frac{\partial x}{\partial \theta}(t_i, \theta_k)\, \Delta \theta.
Ceci nous amène, à chaque itération, à devoir résoudre un problème aux moindres carrés linéaire :
\mathrm{minimiser}~ F_k(\Delta \theta) = \frac{1}{2} \sum_{i=1}^{N} \left\| x_i - x(t_i, \theta_k) - \frac{\partial x}{\partial \theta}(t_i, \theta_k)\, \Delta \theta \right\|^2.
A chaque itération, la fonction à minimiser dépend de \Delta \theta, la valeur du paramètre \theta est fixée à \theta_k.
Le prochain itéré s’écrit :
\theta_{k+1} = \theta_k + \Delta \theta_k
où \Delta \theta_k est la solution du problème aux moindres carrés linéaires faisant intervenir F_k.
Solution du problème linéaire
Le point \Delta \theta_k \in \mathbb{R}^l est obtenu en annulant le gradient de F_k, c’est-à-dire en résolvant :
\begin{aligned} \nabla F_k(\Delta \theta) &= \sum_{i=1}^{N} {\frac{\partial x}{\partial \theta}(t_i, \theta_k)}^T \left( x_i - x(t_i, \theta_k) - \frac{\partial x}{\partial \theta}(t_i, \theta_k)\, \Delta \theta \right) \\ &= \sum_{i=1}^{N} {J_{i,k}}^T \left( x_i - x(t_i, \theta_k) \right) - \left( \sum_{i=1}^{N} {J_{i,k}}^T J_{i,k} \right) \Delta \theta = 0, \end{aligned}
où l’on a introduit la notation
J_{i,k} = \frac{\partial x}{\partial \theta}(t_i, \theta_k).
Si l’on regarde attentivement l’équation ci-dessus en \Delta \theta, nous pouvons remarquer que nous devons résoudre un simple système linéaire de la forme :
A_k \Delta \theta = b_k,
où
A_k = \sum_{i=1}^{N} {J_{i,k}}^T J_{i,k} \in \mathcal{M}_{l}(\mathbb{R}) \quad \text{et} \quad b_k = \sum_{i=1}^{N} {J_{i,k}}^T \left( x_i - x(t_i, \theta_k) \right) \in \mathbb{R}^l.
Présentation de l’algorithme
Voici les étapes principales :
Initialisation : Choisir une estimation initiale \theta_0 des paramètres.
Itération : Pour chaque itération :
Calculer la prédiction x(t_i, \theta_k) pour chaque instant t_i.
Calculer les matrices jacobiennes J_{i,k} des dérivées partielles de x(t_i, \theta) par rapport à \theta, évaluée en \theta_k :
J_{i,k} = \frac{\partial x}{\partial \theta}(t_i, \theta_k) = \begin{pmatrix} \displaystyle\frac{\partial x_1}{\partial \theta_1}(t_i, \theta_k) & \displaystyle\frac{\partial x_1}{\partial \theta_2}(t_i, \theta_k) & \displaystyle\ldots & \displaystyle\frac{\partial x_1}{\partial \theta_l}(t_i, \theta_k) \\[1em] \displaystyle\frac{\partial x_2}{\partial \theta_1}(t_i, \theta_k) & \displaystyle\frac{\partial x_2}{\partial \theta_2}(t_i, \theta_k) & \displaystyle\ldots & \displaystyle\frac{\partial x_2}{\partial \theta_l}(t_i, \theta_k) \\[1em] \displaystyle\vdots & \displaystyle\vdots & \displaystyle\ddots & \displaystyle\vdots \\[1em] \displaystyle\frac{\partial x_n}{\partial \theta_1}(t_i, \theta_k) & \displaystyle\frac{\partial x_n}{\partial \theta_2}(t_i, \theta_k) & \displaystyle\ldots & \displaystyle\frac{\partial x_n}{\partial \theta_l}(t_i, \theta_k) \end{pmatrix} \in \mathcal{M}_{n, l}(\mathbb{\R})
Mettre à jour les paramètres en utilisant la formule :
\theta_{k+1} = \theta_k + \Delta \theta_k
où \Delta \theta_k est la solution du système linéaire :
A_k \Delta \theta = b_k
avec A_k et b_k définis ci-dessus.
Convergence : Répéter jusqu’à ce que les changements dans \theta soient négligeables ou qu’un critère de convergence soit satisfait.
Fonctions Nécessaires
Pour implémenter cet algorithme, deux fonctions principales sont nécessaires :
Fonction de Prédiction :
- Entrée : Paramètre \theta et instant t_i.
- Sortie : Prédictions x(t_i, \theta) en résolvant le problème de Cauchy jusqu’au temps t_i.
Fonction de Dérivée :
- Entrée : Paramètre \theta et instant t_i.
- Sortie : Matrice jacobienne des dérivées partielles de x(t_i, \theta) par rapport à \theta.
Ces fonctions permettent de calculer les prédictions du modèle et les dérivées nécessaires pour mettre à jour les paramètres à chaque itération de l’algorithme de Gauss-Newton. La suite du cours est consacrée à l’implémentation de ces fonctions pour résoudre des problèmes de Kaplan concrets.
Exercice
Considérons le modèle linéaire suivant :
x(t, \theta) = a + b\, t
où :
- x(t, \theta) est la sortie du modèle à l’instant t pour les paramètres \theta.
- \theta = (a, b) sont les paramètres à estimer.
Supposons que nous avons les observations suivantes :
\begin{array}{c|c} t & x \\ \hline 1 & 2 \\ 2 & 3 \\ 3 & 5 \\ \end{array}
Nous cherchons à estimer les paramètres a et b en minimisant la somme des carrés des résidus entre les observations et les prédictions du modèle.
Nous cherchons à estimer les paramètres a et b en minimisant la somme des carrés des résidus entre les observations et les prédictions du modèle.
Étapes de l’algorithme de Gauss-Newton
- Initialisation :
- Choisissons une estimation initiale \theta_0 = (1, 1).
- Itération 1 :
Calculons les prédictions x(t_i, \theta_0) pour chaque instant t_i : x(1, \theta_0) = 1 \cdot 1 + 1 \cdot 1 = 2 x(2, \theta_0) = 1 \cdot 1 + 1 \cdot 2 = 3 x(3, \theta_0) = 1 \cdot 1 + 1 \cdot 3 = 4
Calculons la matrice jacobienne J_{i,0} des dérivées partielles de x(t_i, \theta) par rapport à \theta, évaluée en \theta_0 : J_{i,0} = \begin{pmatrix} \displaystyle\frac{\partial x}{\partial a}(t_i, \theta_0) & \displaystyle\frac{\partial x}{\partial b}(t_i, \theta_0) \end{pmatrix} = \begin{pmatrix} 1 & t_i \end{pmatrix} Donc, J_{1,0} = \begin{pmatrix} 1 & 1 \end{pmatrix}, \quad J_{2,0} = \begin{pmatrix} 1 & 2 \end{pmatrix}, \quad J_{3,0} = \begin{pmatrix} 1 & 3 \end{pmatrix}
Résolvons le système linéaire pour obtenir \Delta \theta_0 : A_0 \Delta \theta = b_0 où A_0 = \sum_{i=1}^{3} {J_{i,0}}^T J_{i,0} = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \begin{pmatrix} 1 & 1 \end{pmatrix} + \begin{pmatrix} 1 \\ 2 \end{pmatrix} \begin{pmatrix} 1 & 2 \end{pmatrix} + \begin{pmatrix} 1 \\ 3 \end{pmatrix} \begin{pmatrix} 1 & 3 \end{pmatrix} = \begin{pmatrix} 3 & 6 \\ 6 & 14 \end{pmatrix} et b_0 = \sum_{i=1}^{3} {J_{i,0}}^T \left( x_i - x(t_i, \theta_0) \right) = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \cdot 0 + \begin{pmatrix} 1 \\ 2 \end{pmatrix} \cdot 0 + \begin{pmatrix} 1 \\ 3 \end{pmatrix} \cdot 1 = \begin{pmatrix} 1 \\ 3 \end{pmatrix} Donc, \begin{pmatrix} 3 & 6 \\ 6 & 14 \end{pmatrix} \Delta \theta = \begin{pmatrix} 1 \\ 3 \end{pmatrix} En résolvant ce système, nous obtenons : \Delta \theta_0 = \begin{pmatrix} -2/3 \\ 0.5 \end{pmatrix}
Mettons à jour les paramètres : \theta_1 = \theta_0 + \Delta \theta_0 = \begin{pmatrix} 1 \\ 1 \end{pmatrix} + \begin{pmatrix} -2/3 \\ 0.5 \end{pmatrix} = \begin{pmatrix} 1/3 \\ 1.5 \end{pmatrix}
- Convergence :
- Répéter les étapes jusqu’à ce que les changements dans \theta soient négligeables ou qu’un critère de convergence soit satisfait.
Tableau de Convergence
On rappelle que le paramètre s’écrit \theta = (a, b). On note alors la variation d’une itération à l’autre \Delta \theta = (\Delta a, \Delta b). Puisque le problème est déjà un problème aux moindres carrés linéaires, nous n’observons aucune amélioration après la première itération, l’algorithme converge en une seule itération.
| Itération | a | b | Écart \Delta a | Écart \Delta b | Valeur du coût F(\theta) |
|---|---|---|---|---|---|
| 0 | 1.000 | 1.000 | - | - | 0.500 |
| 1 | 1/3 | 1.500 | -2/3 | 0.500 | 0.083 |
| 2 | 1/3 | 1.500 | 0.000 | 0.000 | 0.083 |
Cet exemple montre comment appliquer l’algorithme de Gauss-Newton à la main pour estimer les paramètres d’un modèle linéaire simple. Le tableau de convergence permet de visualiser la diminution de l’écart entre les itérations et la minimisation du coût.
Dans la figure ci-dessous, la courbe rouge est donnée par le modèle avec les paramètres initiaux. La courbe bleue est la prédiction du modèle après la première itération de l’algorithme de Gauss-Newton. Nous observons que la prédiction est plus proche des observations.
