# Régression elliptique

La régression elliptique consiste à trouver la « meilleure ellipse », au sens des moindres carrés, décrivant un ensemble de points. C'est un cas de régression géométrique, c'est-à-dire que la distance point-courbe modèle à laquelle on s'intéresse est une distance perpendiculaire à la courbe — méthode des moindres carrés totaux (TLS pour total least squares, ou FLS pour full least squares) —, et non une distance verticale (en y). La régression elliptique est une généralisation de la régression circulaire.

## Méthodes de la distance algébrique

Une ellipse peut être définie par l'équation cartésienne
    $F(x,y)=0$
où F est la formule quadratique :
    $F(x,y)=A_{1}x^{2}+A_{2}xy+A_{3}y^{2}+A_{4}x+A_{5}y+A_{6}.$
La fonction F est également appelée « distance algébrique » du point (x, y) à l'ellipse.
On dispose des coordonnées des n points expérimentaux, notées (x$_i$, y$_i$)$_{1≤i≤n}$. On cherche à minimiser la somme des carrés des distances algébriques, c'est-à-dire
    $\mathrm {S} _{\mathrm {alg} }=\sum _{i=1}^{n}F(x_{i},y_{i})^{2}$.
On peut écrire le problème sous forme matricielle : on définit la matrice des monômes
    $\mathrm {D} ={\begin{pmatrix}x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}&y_{1}&1\\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\x_{n}^{2}&x_{n}y_{n}&y_{n}^{2}&x_{n}&y_{n}&1\\\end{pmatrix}}$
et la matrice des paramètres de l'ellipse
    $\mathbf {a} ={\begin{pmatrix}A_{1}\\\vdots \\A_{6}\end{pmatrix}}$
le problème consiste alors à minimiser
    $\mathrm {S} _{\mathrm {alg} }=\left\|\mathrm {D} \mathbf {a} \right\|^{2}$

### Régression quadratique multilinéaire

La méthode de la régression quadratique consiste à faire une régression linéaire multiple (à l'instar de la régression polynomiale). En effet, on peut transformer l'équation en
    $x^{2}=A_{2}xy+A_{3}y^{2}+A_{4}x+A_{5}y+A_{6}$
en choisissant arbitrairement A$_1$ = –1. On peut alors poser :
    $Y=x^{2},\quad X_{2}=xy\quad X_{3}=y^{2}\quad X_{4}=x\quad X_{5}=y$
on a donc bien un modèle multilinéaire
    $Y=A_{2}X_{2}+A_{3}X_{3}+A_{4}X_{4}+A_{5}X_{5}+A_{6}.$
Le résultat de cette régression pour un nuage de points est donné sur la figure ci-contre.
Présenté comme ceci, la méthode consiste à minimiser
    $\mathrm {S} _{\mathrm {quad} }=\sum _{i}(x_{i}^{2}-(A_{2}x_{i}y_{i}+A_{3}y_{i}^{2}+A_{4}x_{i}+A_{5}y_{i}+A_{6}))^{2}$.
On peut remarquer que l'on pourrait extraire un autre facteur, en posant Y = y$^2$, xy, x ou bien y, et qu'il n'y a pas de raison d'avoir le même résultat à chaque fois.
Le deuxième problème est que la forme quadratique définit de manière générale une conique ; le meilleur candidat peut donc être une hyperbole ou une parabole. Il faut donc ajouter une contrainte propre aux ellipses, à savoir que le discriminant de la forme quadratique associée soit strictement négatif :
    $\Delta =A_{2}^{2}-4A_{1}A_{3}<0.$
Les coefficients A$_i$ sont définis à un facteur multiplicatif près. On peut donc exprimer cette condition par
    $\Delta =-1\ \Longleftrightarrow \ 4A_{1}A_{3}-A_{2}^{2}=1.$
L'ajout de cette contrainte complique la résolution. Plusieurs solutions ont été développées pour éviter de passer par une étape itérative, source potentielle d'instabilité numérique.

### Résolution par décomposition en valeurs singulières

Gander et coll. proposent d'effectuer une décomposition en valeurs singulières de D :
    $D=USV$
où U est une matrice unitaire n × n, V une matrice unitaire 6 × 6 et S est une matrice n × 6 qui contient les valeurs singulières de D. On a alors
    $\mathbf {a} ={\begin{pmatrix}V_{16}\\V_{26}\\\vdots \\V_{66}\end{pmatrix}}$
Les coefficients sont définis à une constante multiplicatrice près. Cette méthode consiste donc, d'une manière ou d'une autre, à appliquer la contrainte
    $\left\|\mathbf {a} \right\|=1.$
Le principal inconvénient de cette méthode est que la contrainte n'est pas invariante par les transformations euclidiennes, et en particulier par les isométries : translation, rotation, symétrie. Ainsi, le demi-grand axe et le demi-petit axe de l'ellipse peuvent être différents si l'on tourne le nuage de points.
Bookstein a proposé à la place d'utiliser la contrainte
    $A_{1}^{2}+A_{2}^{2}+1_{3}^{2}=1$
ce qui revient à imposer une contrainte sur l'équation réduite, qui est elle indépendante des isométries :
    équation cartésienne réduite :$\lambda _{1}{\overline {x}}^{2}+\lambda _{2}{\overline {y}}^{2}+{\overline {c}}^{2}=0$
    contrainte de Bookstein : $\lambda _{1}^{2}+\lambda _{2}^{2}=1.$
Bookstein propose de résoudre ce problème par décomposition spectrale (recherche des valeurs et vecteurs propres), mais Gander et coll. proposent plutôt de résoudre le problème par décomposition en valeurs singulières. Pour cela, on définit la matrice de données modifiée
    $D'={\begin{pmatrix}x_{1}&y_{1}&1&x_{1}^{2}&{\sqrt {2}}x_{1}y_{1}&y_{1}^{2}\\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\x_{n}&y_{n}&1&x_{n}^{2}&{\sqrt {2}}x_{n}y_{n}&y_{n}^{2}\\\end{pmatrix}}$
et les vecteurs de paramètres
    $\mathbf {v} ={\begin{pmatrix}\mathrm {A} _{4}\\\mathrm {A} _{5}\\\mathrm {A} _{6}\end{pmatrix}}$
    $\mathbf {w} ={\begin{pmatrix}\mathrm {A} _{1}\\\mathrm {A} _{2}/{\sqrt {2}}\\\mathrm {A} _{3}\end{pmatrix}}$
et l'on doit donc minimiser
    $D'{\begin{pmatrix}\mathbf {v} \\\mathbf {w} \end{pmatrix}}$
avec la contrainte ∥w∥ = 1. Pour cela, on fait la factorisation QR de D', puis on scinde la matrice R (matrice triangulaire supérieure) pour avoir un bloc R$_{22}$ de dimensions 3 × 3, et donc deux blocs R$_{11}$ et R$_{12}$ de dimension 3 × (n – 3) :
    $\mathrm {R} ={\begin{pmatrix}\mathrm {R} _{11}&\mathrm {R} _{12}\\0&\mathrm {R} _{22}\end{pmatrix}}$
Le problème se ramène alors à chercher le minimum de R$_{22}$w. On effectue pour cela la décomposition en valeurs singulières de R$_{22}$.

### Utilisation des multiplicateurs de Lagrange

Fitzgibbon et coll. a proposé de minimiser la somme des carrés des distances algébriques, et d'utiliser la méthode des multiplicateurs de Lagrange pour intégrer la contrainte. En effet, il s'agit bien de minimiser une fonction φ(A$_1$, …, A$_6$) définie par
    ${\begin{aligned}\varphi \ :\ &\mathbb {R} ^{6}\rightarrow \mathbb {R} \\&(\mathrm {A} _{1},\ldots ,\mathrm {A} _{6})\mapsto \sum _{i=1}^{n}\mathrm {F} (x_{i},y_{i})^{2}\end{aligned}}$
les points (x$_i$, y$_i$)$_{1≤i≤n}$ étant connus, avec une contrainte ψ(A$_1$, …, A$_6$) = 0, la fonction ψ étant définie par
    ${\begin{aligned}\psi \ :\ &\mathbb {R} ^{6}\rightarrow \mathbb {R} \\&(\mathrm {A} _{1},\ldots ,\mathrm {A} _{6})\mapsto \Delta +1\end{aligned}}$
les fonctions φ et ψ étant de classe C$^∞$ (polynômes), donc a fortiori de classe C$^1$.
On note pour la suite a le vecteur de R$^6$
    $\mathbf {a} ={\begin{pmatrix}\mathrm {A} _{1}\\\vdots \\\mathrm {A} _{6}\end{pmatrix}}$
et les matrices représentatives des applications linéaires de R$^6$ dans R$^6$ :
* matrice de conception (design matrix) $\mathbf {D} ={\begin{pmatrix}x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}&x_{1}&y_{1}&1\\\vdots &\vdots &\vdots &\vdots &\vdots &\vdots \\x_{n}^{2}&x_{n}y_{n}&y_{n}^{2}&x_{n}&y_{n}&1\\\end{pmatrix}}$
* matrice de contrainte $\mathbf {C} ={\begin{pmatrix}0&0&2&0&0&0\\0&-1&0&\vdots &\vdots &\vdots \\2&0&0&\vdots &\vdots &\vdots \\0&\cdots &\cdots &\ddots &\vdots &0\\0&\cdots &\cdots &\cdots &\ddots &0\\0&\cdots &\cdots &\cdots &\cdots &0\end{pmatrix}}$
et l'on a donc
    $\varphi (\mathbf {a} )=\left\|\mathrm {Da} \right\|^{2}=^{\mathrm {t} }(\mathrm {Da} )(\mathrm {Da} )$
    $\psi (\mathbf {a} )=^{\mathrm {t} }\mathrm {aDa}$
où $^t$M désigne la transposée de la matrice M. On peut donc poser la fonction L :
    ${\begin{aligned}\mathrm {L} \ :\ &\mathbb {R} ^{6}\times \mathbb {R} \to \mathbb {R} \\&(\mathbf {a} ,\lambda )\mapsto \varphi (\mathbf {a} )+\lambda \cdot \psi (\mathbf {a} )\end{aligned}}$
Si a$_0$ est une solution recherchée (φ est minimale en a$_0$ et a$_0$ satisfait la condition de contrainte), alors il existe une valeur λ$_0$ non nulle telle que la différentielle dL soit nulle en (a$_0$, λ$_0$) : ∂L/∂A$_i$ = 0 pour tout i, et ∂L/∂λ = 0. En calculant les dérivées partielles, on arrive au système d'équations
    $\left\{{\begin{aligned}2^{\mathrm {t} }\mathbf {D} \mathbf {D} \mathbf {a} -2\lambda \mathbf {C} \mathbf {a} &\ =0\\^{\mathrm {t} }\mathbf {a} \mathbf {C} \mathbf {a} &\ =1\end{aligned}}\right.$
En posant la matrice de dispersion ('scatter matrix) S =$^t$DD on a
    $\left\{{\begin{aligned}\mathbf {S} \mathbf {a} &\ =\lambda \mathbf {C} \mathbf {a} &[1]\\^{\mathrm {t} }\mathbf {a} \mathbf {C} \mathbf {a} &\ =1&[2]\end{aligned}}\right.$
Les matrices S et C sont des matrices carrées 6×6.
Notons que l'équation [2] peut s'écrire
    $\mathbf {a} _{i}\mathbf {Ca} i=\lambda _{i}\mathbf {a} _{i}\mathbf {Sa} i$
comme S est en général définie positive, cela revient à dire que λ$_i$ doit être positive.
Il reste donc à résoudre l'équation [1] ; cela peut se faire de plusieurs manières.
Par construction, la matrice S a de grandes chances d'être définie positive, donc inversible. L'équation [1] peut donc s'écrire
    $\mathbf {a} =\lambda \mathbf {S} ^{-1}\mathbf {Ca} \ \Longleftrightarrow \ \mathbf {S} ^{-1}\mathbf {Ca} ={\frac {1}{\lambda }}\mathbf {a}$
si λ est non nulle. Ainsi, a est un vecteur propre de S$^{-1}$Ca, associé à la valeur propre 1/λ.
On peut aussi remarquer que l'équation [1] est un problème aux valeurs propres généralisé, c'est-à-dire à une recherche du sous-espace caractéristique (notion généralisée des valeurs propres et vecteurs propres).
On obtient donc six solutions (λ$_i$, a$_i$) à l'équation [1], mais rien ne garantit qu'elles vérifient l'équation [2]. Cependant, si a$_i$ est un vecteur propre, alors μ$_i$a$_i$ est aussi un vecteur propre pour tout μ$_i$ non nul, il faut donc trouver une valeur de μ$_i$ telle que
    $^{\mathrm {t} }(\mu _{i}\mathbf {a} _{i})\mathbf {C} (\mu _{i}\mathbf {a} _{i})=\mu _{i}^{2}^{\mathrm {t} }\mathbf {a} _{i}\mathbf {C} \mathbf {a} _{i}=1$
La valeur de μ$_i$ est réelle si $^t$a$_i$Ca$_i$ est positif, donc si (1/λ$_i$)$^t$a$_i$Sa$_i$. S étant définie positive, $^t$a$_i$Sa$_i$ est strictement positive, donc
    μ$_i$ est réelle si λ$_i$ est positive.
Donc, une condition nécessaire pour qu'un vecteur de coefficients a$_i$ corresponde à la meilleure ellipse est que ce soit un vecteur propre associé à une valeur propre positive. Fitzgibbon et coll. démontrent qu'il n'y a qu'une seule valeur propre positive, donc que la solution est unique.

### Scission des matrices

Halíř et coll. on proposé des améliorations :
* la matrice C est singulière, et S est presque singulière (elle l'est si tous les points sont exactement sur l'ellipse), la détermination des valeurs propres est donc numériquement instable et peut générer des résultats infinis ou complexes ;
* si tous les points sont exactement sur l'ellipse, la valeur propre est 0 ; donc la valeur propre recherchée est proche de 0, et de fait, l'approximation numérique peut donner des résultats légèrement négatifs, solution qui serait alors rejetée par l'algorithme.
Pour résoudre ces problèmes, ils proposent de scinder les matrices (matrices par blocs) :
* $\mathbf {D} =(\mathbf {D_{1}} |\mathbf {D_{2}} )$ avec $\mathbf {D_{1}} ={\begin{pmatrix}x_{1}^{2}&x_{1}y_{1}&y_{1}^{2}\\\vdots &\vdots &\vdots \\x_{n}^{2}&x_{n}y_{n}&y_{n}^{2}\\\end{pmatrix}}$ et $\mathbf {D_{2}} ={\begin{pmatrix}x_{1}&y_{1}&1\\\vdots &\vdots &\vdots \\x_{n}&y_{n}&1\\\end{pmatrix}}$ ;
* $\mathbf {S} =\left({\begin{array}{c | c}\mathbf {S_{1}} &\mathbf {S_{2}} \\\hline ^{\mathrm {t} }\mathbf {S_{2}} &\mathbf {S_{3}} \end{array}}\right)$ avec S$_1$ = $^t$D$_1$D$_1$, S$_2$ = $^t$D$_1$D$_2$ et S$_3$ = $^t$D$_2$D$_2$ ;
* $\mathbf {C} =\left({\begin{array}{c | c}\mathbf {C_{1}} &0\\\hline 0&0\end{array}}\right)$ avec $\mathbf {C_{1}} ={\begin{pmatrix}0&0&2\\0&-1&0\\2&0&0\end{pmatrix}}$
* $\mathbf {a} ={\begin{pmatrix}\mathbf {a_{1}} \\\hline \mathbf {a_{2}} \end{pmatrix}}$ avec $\mathbf {a_{1}} ={\begin{pmatrix}\mathrm {A} _{1}\\\mathrm {A} _{2}\\\mathrm {A} _{3}\end{pmatrix}}$ et $\mathbf {a_{2}} ={\begin{pmatrix}\mathrm {A} _{4}\\\mathrm {A} _{5}\\\mathrm {A} _{6}\end{pmatrix}}$.
L'équation [1] devient alors le système
    $\left\{{\begin{aligned}\mathbf {S_{1}} \mathbf {a_{1}} +\mathbf {S_{2}} \mathbf {a_{2}} &=\lambda \mathbf {C_{1}} \mathbf {a_{1}} &[3]\\^{\mathrm {t} }\mathbf {S_{2}} \mathbf {a_{1}} +\mathbf {S_{3}} \mathbf {a_{2}} &=0&[4]\\\end{aligned}}\right.$
La matrice S$_3$ correspond à une régression linéaire ; elle est singulière si les points sont strictement alignés, or cela n'a pas de sens de faire une régression elliptique sur des points alignés. On peut donc considérer que S$_3$ est régulière (inversible). La matrice C$_1$ est elle aussi régulière, le système d'équation devient donc
    $\left\{{\begin{aligned}\mathbf {M} \mathbf {a_{1}} &=\lambda \mathbf {a_{1}} &[5]\\-\mathbf {S_{3}} ^{-1}\,^{\mathrm {t} }\mathbf {S_{2}} \mathbf {a_{1}} &=\mathbf {a_{2}} &[6]\\\end{aligned}}\right.$
avec M la matrice de dispersion réduite
    M = C$_1$$^{-1}$(S$_1$ - S$_2$S$_3$$^{-1t}$S$_2$).
L'équation [2] devient
    $^t$a$_1$C$_1$a$_1$ = 1 [7]
On se retrouve donc à résoudre le système d'équations {[5] ; [6] ; [7]}, soit :
* [5] : déterminer les valeurs propres et vecteurs propres de M ;
* trouver a$_1$ : retenir la valeur propre positive, ou mieux : [7] trouver l'unique vecteur tel que $^t$a$_1$C$_1$a$_1$ > 0 ;
* [6] : calculer a$_2$ ;
* rassembler a$_1$ et a$_2$ pour former le vecteur a.

## Méthode des moindres carrés totaux

La méthode des moindres carrés totaux est, comme dans le cas du cercle, non linéaire. On a donc recours à un algorithme itératif.
Le principal problème est de déterminer la distance d'un point à l'ellipse modèle. La méthode la plus simple consiste à prendre une équation paramétrique de l'ellipse :
    ${\begin{pmatrix}x(\varphi )\\y(\varphi )\end{pmatrix}}={\begin{pmatrix}x_{\mathrm {c} }\\y_{\mathrm {c} }\end{pmatrix}}+\mathrm {Q} (\alpha ){\begin{pmatrix}a\cos \varphi \\b\sin \varphi \end{pmatrix}}$
où (x$_c$, y$_c$) sont les coordonnées du centre de l'ellipse et Q(α) est la matrice de rotation d'angle α (inclinaison de l'ellipse).
On se retrouve ici avec n + 6 inconnues : les six paramètres de l'ellipse (x$_c$, y$_c$, a, b, α) et les n paramètres φ$_i$, le point (x(φ$_i$), x(φ$_i$)) étant le point de l'ellipse le plus proche du point expérimental i.
Pour initialiser les paramètres de l'ellipse, on peut utiliser une méthode de distance algébrique, ou bien une régression circulaire ; le cas du cercle pouvant donner une matrice jacobienne singulière, il peut être nécessaire de démarrer en « elliptisant » le cercle, par exemple en créant de manière arbitraire une ellipse dont le demi-grand axe a le rayon du cercle et le demi-petit axe vaut la moitié.
Pour initialiser φ$_i$, on peut utiliser l'angle par rapport à l'axe x du segment reliant le centre initial de l'ellipse au point expérimental i.
On peut utiliser les méthodes itératives classiques (méthodes de Gauss-Newton ou de Levenberg-Marquardt).

## Applications
### Analyse d'image

Une ellipse peut être considéré comme un cercle selon une « vue inclinée » : c'est la projection orthogonale d'un cercle sur un plan sécant au plan le contenant. C'est donc une figure qui est susceptible d'apparaître dans de nombreuses images.
Cela peut être utilisé pour des algorithmes de reconnaissance de forme, par exemple reconnaître l'ovale des visages sur une photographie, pour de l'imagerie médicale, des inspections industrielles, la conversion d'une image matricielle en image vectorielle, ou encore en archéologie — pour déterminer la taille d'une poterie à partir d'un fragment, le col de la poterie formant un arc de cercle qui, du fait de la perspective, est vu comme un arc d'ellipse.