Reconstruire un signal à partir de mesures imparfaites

Le problème part d'une reconstruction exacte dans un modèle de dimension 3, puis montre ce qui change lorsque cinq mesures ne sont plus compatibles. L'idée centrale n'est pas de supposer connue la méthode des moindres carrés : les équations normales apparaissent en étudiant l'erreur le long de chaque direction $h$ de l'espace des paramètres.

I - Quand les mesures sont exactes

On écrit

$$s(t)=\alpha+\beta t+\gamma t^2.$$

La mesure en $0$ donne immédiatement $\alpha=10$. Les deux autres donnent

$$10-\beta+\gamma=13, \qquad 10+\beta+\gamma=11,$$

soit

$$-\beta+\gamma=3, \qquad \beta+\gamma=1.$$

En additionnant, $2\gamma=4$, donc $\gamma=2$, puis $\beta=-1$.

Ainsi

$$x= \begin{pmatrix} 10\\-1\\2 \end{pmatrix} \qquad\text{et}\qquad s(t)=10-t+2t^2.$$

(a) Pour

$$x= \begin{pmatrix}\alpha\\\beta\\\gamma\end{pmatrix},$$

on a

$$A_0x= \begin{pmatrix} \alpha-\beta+\gamma\\ \alpha\\ \alpha+\beta+\gamma \end{pmatrix}.$$

Le système $A_0x=y_0$ exprime donc exactement les trois mesures. Ici

$$A_0\in\mathcal M_{3,3}(\mathbb R), \qquad x\in\mathbb R^3, \qquad y_0\in\mathbb R^3.$$

(b) Si $A_0x=0$, le polynôme $p(t)=\alpha+\beta t+\gamma t^2$ s'annule en $-1,0,1$. Un polynôme non nul de degré au plus 2 ne peut avoir trois racines distinctes. Donc $p=0$, d'où $\alpha=\beta=\gamma=0$.

Ainsi $\ker A_0=\{0\}$. Comme $A_0$ est une matrice carrée, elle est inversible.

Supposons

$$A(t_1,t_2,t_3)x=0.$$

Le polynôme

$$p(t)=\alpha+\beta t+\gamma t^2$$

s'annule alors aux trois instants distincts $t_1,t_2,t_3$. Puisque $\deg p\le 2$, cela impose $p=0$, donc $x=0$.

L'application $x\mapsto A(t_1,t_2,t_3)x$ est donc injective.

Si deux paramètres $x$ et $x'$ produisaient les mêmes trois mesures, on aurait

$$A(t_1,t_2,t_3)(x-x')=0,$$

puis $x-x'=0$. Comme l'application est linéaire de $\mathbb R^3$ dans $\mathbb R^3$, son injectivité implique sa bijectivité. Pour n'importe quel triplet de valeurs mesurées, le système possède donc une unique solution. Trois mesures exactes en trois instants distincts déterminent un unique signal du modèle.

Avec une seule mesure en $t_1$, tout multiple non nul de $X-t_1$ donne une perturbation du signal qui ne modifie pas cette mesure.

Avec deux mesures en deux instants distincts $t_1,t_2$, le polynôme non nul

$$q(X)=(X-t_1)(X-t_2)$$

est de degré 2 et s'annule aux deux instants. Si $s$ est compatible avec les deux mesures, alors $s+\lambda q$ l'est aussi pour tout $\lambda\in\mathbb R$. L'unicité est donc impossible en général.

Des mesures exactes supplémentaires ne changent pas le signal lorsque les données sont compatibles : dès que trois instants distincts ont fixé l'unique polynôme, chaque mesure supplémentaire doit simplement confirmer ce même polynôme. En revanche, si une mesure supplémentaire ne coïncide pas avec sa valeur prédite, le système devient incompatible.

II - Cinq mesures qui ne peuvent pas toutes être exactes

(a) Écrivons

$$q(t)=\alpha+\beta t+\gamma t^2.$$

Les mesures centrales donnent

$$q(-1)=11,\qquad q(0)=10,\qquad q(1)=13.$$

Ainsi $\alpha=10$ et

$$-\beta+\gamma=1, \qquad \beta+\gamma=3.$$

On obtient $\gamma=2$ et $\beta=1$, donc

$$q(t)=10+t+2t^2.$$

(b) Ce polynôme donne

$$q(2)=10+2+8=20,$$

alors que la mesure à $t=2$ vaut $15$.

Aucun polynôme de degré au plus 2 ne peut donc reproduire simultanément les cinq mesures : le système $Ax=y$ est incompatible.

(a) Les cinq composantes de $Ax$ sont les valeurs prédites aux instants $-2,-1,0,1,2$. Par conséquent

$$\sum_{i=1}^5 r_i(x) = \sum_{i=1}^5\bigl(\alpha+\beta t_i+\gamma t_i^2-y_i\bigr).$$

Or

$$\sum t_i=0,\qquad \sum t_i^2=4+1+0+1+4=10,\qquad \sum y_i=70.$$

Donc

$$\sum_{i=1}^5 r_i(x)=5\alpha+10\gamma-70.$$

L'équation $\sum r_i=0$ n'est qu'une contrainte linéaire sur trois paramètres : elle possède une infinité de solutions. Surtout, de grandes erreurs positives et négatives peuvent se compenser. Une somme algébrique n'évalue donc pas correctement la taille globale des écarts.

(b) On a

$$J(x) = \sum_{i=1}^5 \bigl(\alpha+\beta t_i+\gamma t_i^2-y_i\bigr)^2,$$

c'est-à-dire ici

$$\begin{aligned} J(x)=& (\alpha-2\beta+4\gamma-21)^2 +(\alpha-\beta+\gamma-11)^2\\ &+(\alpha-10)^2 +(\alpha+\beta+\gamma-13)^2 +(\alpha+2\beta+4\gamma-15)^2. \end{aligned}$$

Chaque terme est positif ou nul. Deux erreurs de signes opposés ne peuvent donc plus s'annuler entre elles dans le critère.

III - Faire émerger la meilleure reconstruction

(a) Comme

$$r(x+\lambda h)=A(x+\lambda h)-y=r(x)+\lambda Ah,$$

on obtient

$$\begin{aligned} J(x+\lambda h) &=\|r(x)+\lambda Ah\|^2\\ &=\|r(x)\|^2 +2\lambda\langle r(x),Ah\rangle +\lambda^2\|Ah\|^2. \end{aligned}$$

D'où

$$J(x+\lambda h) = J(x) +2\lambda\langle r(x),Ah\rangle +\lambda^2\|Ah\|^2.$$

(b) Si $x$ minimise $J$, alors pour chaque $h$, la fonction

$$\lambda\longmapsto J(x+\lambda h)$$

a un minimum en $\lambda=0$.

Si $Ah\ne0$, le trinôme a un coefficient dominant strictement positif. Son sommet est en $0$ si et seulement si son coefficient linéaire est nul, donc

$$\langle r(x),Ah\rangle=0.$$

Si $Ah=0$, cette égalité est immédiate. Ainsi, dans tous les cas,

$$\langle r(x),Ah\rangle=0 \quad\text{pour tout }h\in\mathbb R^3.$$

(c) On utilise

$$\langle r(x),Ah\rangle = r(x)^\mathsf{T}Ah = h^\mathsf{T}A^\mathsf{T}r(x).$$

Cette quantité est nulle pour tout $h$. En prenant successivement les trois vecteurs de la base canonique de $\mathbb R^3$, on obtient

$$A^\mathsf{T}r(x)=0.$$

Comme $r(x)=Ax-y$,

$$A^\mathsf{T}(Ax-y)=0,$$

donc

$$A^\mathsf{T}Ax=A^\mathsf{T}y.$$

Ici $A\in\mathcal M_{5,3}(\mathbb R)$, donc $A^\mathsf{T}\in\mathcal M_{3,5}(\mathbb R)$, $A^\mathsf{T}A\in\mathcal M_{3,3}(\mathbb R)$ et $A^\mathsf{T}y\in\mathbb R^3$.

(a) Si $x$ vérifie les équations normales, alors $A^\mathsf{T}r(x)=0$. Pour tout $h$,

$$\begin{aligned} J(x+h) &=\|r(x)+Ah\|^2\\ &=J(x)+2\langle r(x),Ah\rangle+\|Ah\|^2. \end{aligned}$$

Or

$$\langle r(x),Ah\rangle = h^\mathsf{T}A^\mathsf{T}r(x)=0.$$

Ainsi

$$J(x+h)=J(x)+\|Ah\|^2.$$

(b) Comme $\|Ah\|^2\ge0$, on a $J(x+h)\ge J(x)$ pour tout $h$. Le point $x$ est donc un minimiseur de $J$.

(c) Si $A$ est injective et $h\ne0$, alors $Ah\ne0$, donc

$$\|Ah\|^2>0.$$

Par conséquent

$$J(x+h)>J(x)$$

pour tout $h\ne0$. Le minimiseur est unique.

(d) Supposons

$$A^\mathsf{T}Az=0.$$

En multipliant à gauche par $z^\mathsf{T}$,

$$0=z^\mathsf{T}A^\mathsf{T}Az=\|Az\|^2.$$

Donc $Az=0$. L'injectivité de $A$ donne $z=0$.

Ainsi

$$\ker(A^\mathsf{T}A)=\{0\}.$$

La matrice $A^\mathsf{T}A$ est carrée de taille $3$, donc elle est inversible.

(a) Si $Ax=0$, le polynôme

$$p(t)=\alpha+\beta t+\gamma t^2$$

s'annule aux cinq instants distincts $-2,-1,0,1,2$. Il est donc nul, et $x=0$. La matrice $A$ est injective. Comme elle possède trois colonnes, son rang vaut donc

$$\operatorname{rg}(A)=3.$$

(b) Les symétries des instants donnent

$$\sum t_i=0,\qquad \sum t_i^2=10,\qquad \sum t_i^3=0,\qquad \sum t_i^4=34.$$

Donc

$$A^\mathsf{T}A= \begin{pmatrix} 5&0&10\\ 0&10&0\\ 10&0&34 \end{pmatrix}.$$

D'autre part,

$$\sum y_i=70,$$ $$\sum t_i y_i = -2\cdot21-1\cdot11+1\cdot13+2\cdot15 = -10,$$

et

$$\sum t_i^2 y_i = 4\cdot21+1\cdot11+1\cdot13+4\cdot15 = 168.$$

Ainsi

$$A^\mathsf{T}y= \begin{pmatrix} 70\\-10\\168 \end{pmatrix}.$$

(c) Les équations normales sont

$$\begin{cases} 5\alpha+10\gamma=70,\\ 10\beta=-10,\\ 10\alpha+34\gamma=168. \end{cases}$$

La deuxième équation donne $\beta=-1$. La première donne

$$\alpha+2\gamma=14.$$

En remplaçant $\alpha=14-2\gamma$ dans la troisième,

$$10(14-2\gamma)+34\gamma=168,$$

donc

$$14\gamma=28, \qquad \gamma=2, \qquad \alpha=10.$$

Finalement

$$\widehat x= \begin{pmatrix} 10\\-1\\2 \end{pmatrix}.$$

(a) Le signal reconstruit est

$$\widehat s(t)=10-t+2t^2.$$

Aux cinq instants,

$$A\widehat x= \begin{pmatrix} 20\\13\\10\\11\\16 \end{pmatrix}.$$

(b) Le résidu est donc

$$r(\widehat x) = A\widehat x-y = \begin{pmatrix} -1\\2\\0\\-2\\1 \end{pmatrix}.$$

Ainsi

$$J(\widehat x) = (-1)^2+2^2+0^2+(-2)^2+1^2 = 10.$$

(c) On calcule

$$A^\mathsf{T}r(\widehat x) = \begin{pmatrix} -1+2+0-2+1\\ 2-2+0-2+2\\ -4+2+0-2+4 \end{pmatrix} = \begin{pmatrix} 0\\0\\0 \end{pmatrix}.$$

Le calcul numérique exact confirme donc les équations normales.

IV - Ce que signifie « reconstruire au mieux »

(a) Les colonnes de $A$ sont

$$\begin{pmatrix}1\\1\\1\\1\\1\end{pmatrix}, \qquad \begin{pmatrix}t_1\\t_2\\t_3\\t_4\\t_5\end{pmatrix}, \qquad \begin{pmatrix}t_1^2\\t_2^2\\t_3^2\\t_4^2\\t_5^2\end{pmatrix}.$$

L'égalité $A^\mathsf{T}r(\widehat x)=0$ équivaut donc à

$$\sum_{i=1}^5 r_i=0, \qquad \sum_{i=1}^5 t_i r_i=0, \qquad \sum_{i=1}^5 t_i^2r_i=0.$$

(b) Modifier $\alpha$, $\beta$ ou $\gamma$ fait varier le vecteur des prédictions suivant une combinaison linéaire des trois colonnes de $A$. Or le résidu est orthogonal à chacune de ces colonnes ; il est donc orthogonal à toute combinaison linéaire de celles-ci, c'est-à-dire à $\operatorname{Im}(A)$.

Ainsi

$$r(\widehat x)=A\widehat x-y$$

est orthogonal à toutes les variations permises par le modèle.

(c) Le vecteur $A\widehat x$ appartient à $\operatorname{Im}(A)$ et

$$y-A\widehat x=-r(\widehat x)$$

est orthogonal à $\operatorname{Im}(A)$. De plus, la question 8 a montré directement que $A\widehat x$ est le point de $\operatorname{Im}(A)$ qui minimise la distance à $y$.

Autrement dit, $A\widehat x$ est le projeté orthogonal de $y$ sur $\operatorname{Im}(A)$. La méthode des moindres carrés apparaît donc ici comme un problème de projection, mais cette interprétation n'a pas été nécessaire pour dériver les équations normales.

(a) Pour

$$x_c= \begin{pmatrix}10\\1\\2\end{pmatrix},$$

on obtient

$$Ax_c= \begin{pmatrix} 16\\11\\10\\13\\20 \end{pmatrix}.$$

Donc

$$Ax_c-y= \begin{pmatrix} -5\\0\\0\\0\\5 \end{pmatrix}$$

et

$$J(x_c)=25+25=50.$$

(b) Or

$$J(\widehat x)=10<50.$$

Le polynôme $x_c$ annule exactement trois résidus, mais il concentre toute l'incompatibilité sur les deux mesures extrêmes. Le minimiseur des moindres carrés accepte au contraire de petits écarts sur plusieurs mesures pour réduire fortement l'erreur quadratique totale.

(a) Si $y=Ax_0$, alors

$$J(x_0)=\|Ax_0-y\|^2=0.$$

Comme $J\ge0$, la valeur minimale possible est donc $0$.

(b) Le vecteur $x_0$ vérifie

$$A^\mathsf{T}Ax_0=A^\mathsf{T}y.$$

Comme $A$ est injective, $A^\mathsf{T}A$ est inversible d'après la question 8. Les équations normales possèdent donc une unique solution, nécessairement $x_0$.

Ainsi la méthode ne modifie pas la reconstruction lorsque des données exactes existent : elle prolonge cette reconstruction au cas incompatible en remplaçant l'égalité impossible $Ax=y$ par la minimisation de $\|Ax-y\|^2$.

(a) Une application linéaire injective de $\mathbb R^n$ vers $\mathbb R^m$ envoie une famille libre de $n$ vecteurs sur une famille libre de $n$ vecteurs de $\mathbb R^m$. Il faut donc

$$n\le m.$$

(b) L'argument de la question 8 vaut sans changement. L'injectivité de $A$ entraîne que $A^\mathsf{T}A$ est inversible, car

$$A^\mathsf{T}Az=0 \Longrightarrow \|Az\|^2=0 \Longrightarrow Az=0 \Longrightarrow z=0.$$

Le système

$$A^\mathsf{T}A\widehat x=A^\mathsf{T}y$$

admet donc une unique solution.

Cette solution minimise $J(x)=\|Ax-y\|^2$ par l'identité

$$J(\widehat x+h)=J(\widehat x)+\|Ah\|^2.$$

L'injectivité rend l'inégalité stricte pour $h\ne0$, donc le minimiseur est unique.

(c) Comme $A^\mathsf{T}A$ est inversible,

$$\widehat x=(A^\mathsf{T}A)^{-1}A^\mathsf{T}y.$$

(d) Si $y\in\operatorname{Im}(A)$, il existe $x_0$ tel que $y=Ax_0$. La valeur minimale est alors $0$ et, par unicité,

$$\widehat x=x_0.$$

La formule précédente retrouve donc exactement le paramètre ayant produit les données.

Sources

[SOURCE] G. Strang, MIT OpenCourseWare, « Projection Matrices and Least Squares » : Lien direct.

[SOURCE] NIST/SEMATECH e-Handbook of Statistical Methods, « Linear Least Squares Regression » : Lien direct.