# Algorithmes de visualisation de non-identifiabilité

Nathanaël Munier, Emmanuel Soubies, Pierre Weiss

► To cite this version:

Nathanaël Munier, Emmanuel Soubies, Pierre Weiss. Algorithmes de visualisation de non-identifiabilité. XXIXème Colloque Francophone de Traitement du Signal et des Images (GRETSI 2023), Groupe de Recherche et d'Etudes de Traitement du Signal et des Images, Aug 2023, Grenoble, France. hal-04138308v2

HAL Id: hal-04138308

https://hal.science/hal-04138308v2

Submitted on 19 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Algorithmes de visualisation de non-identifiabilité

Nathanaël MUNIER¹,² Emmanuel SOUBIES³,² Pierre WEISS²,³

¹IMT, Université de Toulouse, CNRS, France

²Centre de Biologie Intégrative, Laboratoire MCD, Université de Toulouse, CNRS, France

³IRIT, Université de Toulouse, CNRS, France

**Résumé** – L'identification d'un ensemble de paramètres $x$ à partir de données observées $y$ est souvent réduite à la résolution d'un problème du type $\Phi(x) = y$, où $\Phi$ est une application continue. Cette équation apparaît dans des domaines tels que les problèmes inverses (aveugles), l'assimilation de données ou encore le contrôle optimal. Lorsque la solution n'est pas unique, il peut être important d'en informer l'utilisateur et de le renseigner sur l'ensemble des solutions plausibles. Nous proposons ici un algorithme permettant de décrire localement un ensemble de solutions de cette équation. Il repose essentiellement sur la différenciation automatique et sur une paramétrisation locale de la variété des solutions. Nous illustrons l'intérêt de cette approche sur deux exemples jouets.

**Abstract** – The identification of a set of parameters $x$ from observed data $y$ is often reduced to the resolution of a problem of the type $\Phi(x) = y$, where $\Phi$ is a continuous application. This equation appears in domains such as (blind) inverse problems, data assimilation, or optimal control. When the solution is not unique, it may be important to inform the user about the set of plausible solutions. We propose here an algorithm to locally describe a set of solutions of this equation. It is essentially based on automatic differentiation and on a local parameterization of the solution manifold. We illustrate the interest of this approach on two toy examples.

## 1 Introduction

Dans cet article, nous nous intéressons au problème d'identifiabilité associé au modèle générique

$$y = \Phi(x), \quad (1)$$

où $\Phi \in C^1(\mathbb{R}^N, \mathbb{R}^M)$. Plus précisément, nous nous intéressons à la notion d'identifiabilité locale énoncée ci-dessous.

**Définition 1.1.** Le vecteur $\bar{x} \in \mathbb{R}^N$ est dit localement identifiable pour le modèle (1) s'il existe $r > 0$ tel que la boule $B = B(\bar{x}, r)$ de centre $\bar{x}$ et de rayon $r$ satisfait

$$\forall x \in \mathcal{B}, \quad \Phi(x) = \Phi(\bar{x}) \implies x = \bar{x}. \quad (2)$$

Cette notion d'identifiabilité joue un rôle central dans la détermination du caractère bien ou mal posé d'un problème inverse consistant à estimer $x$ à partir de données $y$.

Dans ce travail, nous nous intéressons à l'ensemble de non-identifiabilité en $\bar{x}$ défini par

$$M_{\bar{x}} \stackrel{\text{def}}{=} \{x \in \mathbb{R}^N ; \Phi(x) = \Phi(\bar{x})\}. \quad (3)$$

Un résultat bien connu [5, chap. 2, prop. 12] montre que cet ensemble est – localement autour de $\bar{x}$ – une sous-variété différentielle régulière de $\mathbb{R}^N$. Nous en déterminons une paramétrisation locale (théorème 2.1). Nous exploitons ensuite cette paramétrisation dans la section 3 pour développer une méthode numérique permettant de décrire $M_{\bar{x}}$. Enfin, dans la section 4, nous mettons en œuvre cette méthode afin de mettre en lumière des problèmes de non-identifiabilité pour deux exemples illustratifs.

## 2 Description théorique de $M_{\bar{x}}$

Nous notons $J_{\Phi}(x) \in \mathbb{R}^{M \times N}$ la matrice jacobienne de $\Phi$ en $x$ et travaillons sous l'hypothèse suivante.

**Hypothèse 2.1.** La fonction $x \mapsto \text{rang}(J_{\Phi}(x))$ est localement constante sur un voisinage de $\bar{x}$.

Notons que cette hypothèse est peu restrictive étant donné que pour une fonction $\Phi$ de classe $C^1$, elle est vérifiée sur un ouvert dense de $\mathbb{R}^N$. Ainsi, seul un ensemble de mesure nulle est exclu.

Le théorème 2.1 montre que l'ensemble de non-identifiabilité $M_{\bar{x}}$ est, localement, une sous-variété de $\mathbb{R}^N$ et en donne une paramétrisation locale. Il s'inspire et génère le théorème 3.1 de [4]. En particulier, il ne requiert pas que $N > M$ ni que $\Phi$ soit une submersion. Nous avons essayé d'obtenir le résultat le plus générique possible. Il est probablement déjà connu, mais nous ne l'avons pas trouvé sous cette forme précise.

**Théorème 2.1.** Sous l'hypothèse 2.1, $M_{\bar{x}}$ est localement une sous-variété de $\mathbb{R}^N$ de dimension $D = \dim(\ker(J_{\Phi}(\bar{x})))$ ayant la même régularité que $\Phi$.

Par ailleurs, $M_{\bar{x}}$ est paramétrisable localement par le diféomorphisme

$$\psi \stackrel{\text{def}}{=} H^{-1} \circ \iota : \mathcal{U} \to M_{\bar{x}} \quad (4)$$

où

* $\mathcal{U}$ est un ouvert de $\mathbb{R}^D$

* $\iota : \mathbb{R}^D \to \mathbb{R}^{M+D}$ est l'injection canonique qui envoie $\mathbb{R}^D$ sur $\{0_{\mathbb{R}^M}\} \times \mathbb{R}^D$

* $H : \mathbb{R}^N \to \mathbb{R}^{M+D}$ est une immersion locale en $\bar{x}$ définie par $H(x) = (\Phi(x) - \Phi(\bar{x}), V^T(x - \bar{x}))$ avec
<figure><img src="image_1.png" /><figcaption>FIGURE 1 : Paramétrisation locale de $M_{\bar{x}}$. La variété $M_{\bar{x}}$ est relevée à l'aide de l'immersion $H$ préservant sa structure construite à partir de $\Phi$. L'injection $\iota$ permet de ne récupérer que les $D$ variables utiles à la paramétrisation.</figcaption></figure>

$V \in \mathcal{L}(\mathbb{R}^D, \mathbb{R}^N)$ un isomorphisme linéaire entre $\mathbb{R}^D$ et $\ker(J_{\Phi}(\bar{x})) = \{0\}$.

La paramétrisation locale de $M_{\bar{x}}$ établie dans le théorème 2.1 est illustrée sur la figure 1. On notera en particulier que $\psi(0) = \bar{x}$. Une conséquence directe de ce résultat est énoncée dans le corollaire 2.1.

**Corollaire 2.1.** *Sous l'hypothèse 2.1, le vecteur $\bar{x} \in \mathbb{R}^N$ est localement identifiable pour le modèle (1) si et seulement si $\ker(J_{\Phi}(\bar{x})) = \{0\}$.*

Notons que l'inverse $H^{-1}$ peut être multivoque. Cependant, il devient univoque si on restreint son image à un voisinage de $\bar{x}$. C'est le sens qu'il faut lui donner dans le théorème 2.1. Notons aussi qu'une faiblesse du théorème 2.1 est que nous ne précisons pas la taille de l'ouvert $\mathcal{U}$.

## 3 Exploration numérique de $M_{\bar{x}}$

Le résultat du théorème 2.1, indique que l'ensemble $\{\psi(z), z \in \mathcal{U}\}$ où $\psi$ est le difféomorphisme défini en (4) décrit un sous-ensemble de $M_{\bar{x}}$ autour de $\bar{x}$. On peut donc échantillonner localement la variété $M_{\bar{x}}$ en discrétisant $\mathcal{U}$. En reprenant les notations du théorème 2.1, nous savons que<fnref n="1" />

$$ \psi(z) = \underset{x \in \mathbb{R}^N}{\operatorname{argmin}} \frac{1}{2} \|H(x) - \iota(z)\|_2^2. \quad (5) $$

L'injectivité de $\psi$ assure la génération de points distincts sur la variété $M_{\bar{x}}$. Pour échantillonner la variété, il suffit donc de choisir un ensemble de points $(z_k)_k$ distincts sur $\mathcal{U}$. De plus, nous allons voir que le problème (5) peut être reformulé sous une forme permettant d'assurer une convergence linéaire locale des algorithmes de premier ordre.

### 3.1 Une reformulation fortement convexe

Soit $V \in \mathbb{R}^{N \times D}$ (respectivement $V^{\perp} \in \mathbb{R}^{N \times (N-D)}$) une matrice dont les colonnes forment une base orthonormée de $\ker(J_{\Phi}(\bar{x}))$ (respectivement de $\ker(J_{\Phi}(\bar{x}))^{\perp}$). Nous pouvons établir le résultat suivant.

**Théorème 3.1.** *Considérons la fonction $F : \mathbb{R}^D \times \mathbb{R}^{N-D} \to \mathbb{R}$ suivante :*

$$ F(z, z') \stackrel{\text{def}}{=} \frac{1}{2} \| \Phi(\bar{x} + Vz + V^{\perp}z') - \Phi(\bar{x}) \|_2^2 \quad (6) $$

**Algorithme 1 : Exploration de $M_{\bar{x}}$**

1. **Entrée** : $\bar{x} \in \mathbb{R}^N$

2. Déterminer la dimension $D$ de $\mathcal{H} = \ker(J_{\Phi}(\bar{x}))$

3. Construire $V (V^{\perp})$, une base orthonormée de $\mathcal{H} (\mathcal{H}^{\perp})$

4. Définir une discrétisation $\{z_1, \dots, z_K\}$ de $\mathbb{R}^D$ ordonnée par norme (voir fig. 2)

5. Poser $z_0 = 0_{\mathbb{R}^D}$, $\eta(z_0) = 0_{\mathbb{R}^{N-D}}$ et $Z_1 = \{z_0\}$

6. **pour** $k = 1, \dots, K$ **faire**

7. Soit $z^{(0)}$ le point dans $Z_k$ le plus proche de $z_k$

8. Définir $z'_0 = \eta(z^{(0)})$ (déjà calculé)

9. Calculer $\eta(z_k)$ la solution de $(\mathcal{P}(z_k))$ obtenue par descente de gradient initialisée en $z'_0$

10. Poser $x_k \stackrel{\text{def}}{=} \bar{x} + Vz_k + V^{\perp}\eta(z_k)$

11. Poser $Z_{k+1} = Z_k \cup \{z_k\}$

12. **fin**

13. **Sortie** : la suite $(x_k)_{1 \le k \le K}$.

<figure><img src="image_2.png" /><figcaption>FIGURE 2 : Discrétisation de la paramétrisation de la variété de proche en proche.</figcaption></figure>

et le problème de minimisation

$$ \min_{z' \in \mathbb{R}^{N-D}} F(z, z') \qquad (\mathcal{P}(z)) $$

Pour $z \in \mathcal{U}$ :

* La fonction $F(z, \cdot)$ admet un unique minimiseur noté $\eta(z) \in \mathbb{R}^{N-D}$ de norme minimale.

* Le point $\psi(z)$ vérifie

$$ \psi(z) = \bar{x} + Vz + V^{\perp}\eta(z). \quad (7) $$

* Si $\Phi \in \mathcal{C}^2$, la fonction $F$ est fortement convexe sur un ouvert autour de $(0, 0)$. Ainsi, une descente de gradient sur $F(z, \cdot)$ initialisée en $(z, 0)$ converge linéairement vers $\eta(z)$.

L'idée de cette reformulation provient de la décomposition

$$ x = \bar{x} + \Pi(x - \bar{x}) + \Pi^{\perp}(x - \bar{x}) \quad (8) $$

où $\Pi = VV^T$ est la projection sur le noyau $\ker(J_{\Phi}(\bar{x}))$ et $\Pi^{\perp}$ est la projection sur l'orthogonal de $\ker(J_{\Phi}(\bar{x}))$.

Le problème $(\mathcal{P}(z))$ présente deux avantages par rapport au problème initial (5). Tout d'abord, la dimension de l'espace de recherche est réduite de $N$ à $N - D$. Deuxièmement, une méthode de gradient converge linéairement. Ainsi, on peut garantir de retrouver un voisinage de $M_{\bar{x}}$ contenant $\bar{x}$.

Le pseudo-code de la méthode d'exploration de $M_{\bar{x}}$ proposée est décrit dans l'algorithme 1. Il consiste en l'évaluation successive de $\psi$ sur un ensemble de points $(z_k)_{k \le K}$ de l'ouvert $\mathcal{U} \subset \mathbb{R}^D$. Le choix des points se fait séquentiellement. Les descentes de gradient sont ainsi initialisées à l'aide de solutions précédemment calculés. En pratique, on procède à une discrétisation par propagation de proche en proche (cf. figure 2).

<footnote n="1">¹Il peut exister plusieurs minimiseurs globaux au problème (5), cependant, il y en a un unique plus proche de $\bar{x}$ et c'est celui qui nous intéresse ici.</footnote>
Afin de mettre en œuvre cet algorithme, il est essentiel d'être en mesure de déterminer efficacement une base $V$ du noyau $\text{ker}(J_{\Phi}(\bar{x}))$ (étape 3 de l'algorithme 1). Nous étudions plus en détail cette problématique dans la section suivante.

## 3.2 Détermination d'une base de noyau

On note dans cette section $\Gamma \stackrel{\text{def}}{=} J_{\Phi}(\bar{x}) \in \mathbb{R}^{M \times N}$ la matrice jacobienne de $\Phi$ au point $\bar{x}$. On rappelle aussi (cf. théorème 2.1) que $D = \text{dim}(\text{ker } \Gamma = N - \text{rang}(\Gamma))$. En pratique, on peut avoir accès à chaque ligne de la matrice jacobienne $\Gamma = J_{\Phi}(\bar{x})$ par différentiation automatique. Cependant, ceci pose deux difficultés potentielles : i) si $M$ est grand, le coût calculatoire est élevé ii) si le produit $M \cdot N$ est grand, $\Gamma$ peut ne pas être stockable. Suivant les dimensions $M$ et $N$, on favorisera une des approches suivantes.

1. Si on peut calculer intégralement la matrice Jacobienne $\Gamma$, deux options sont possibles :

* Une décomposition en valeurs singulières (SVD) fournit une base de ker $\Gamma$. Cette approche présente l'avantage d'être très simple à implémenter avec les librairies standards.

* Si le temps de calcul pour effectuer une SVD est trop élevé, l'algorithme de Li et Zeng [3] qui repose sur une décomposition QR permet de déterminer successivement les vecteurs associés aux plus petites valeurs singulières.

2. Si les dimensions sont trop grandes pour calculer ou stocker $\Gamma$, on peut utiliser la différentiation automatique pour calculer des produits matrice-vecteur avec des vecteurs arbitraires. Le coût est proportionnel au calcul de $\Phi$. Ceci permet d'exploiter l'algorithme de *randomized SVD* [1]. L'algorithme 2 en est une amélioration qui exploite le fait qu'on souhaite uniquement trouver des éléments du noyau de $\Gamma$.

Cette alternative à la *randomized SVD* est intéressante lorsque $N < M$. Elle repose sur les remarques suivantes :

* Pour déterminer le noyau d'une matrice, il suffit d'avoir accès à une base de ses lignes.

* La multiplication à gauche par une matrice inversible préserve le noyau.

De plus, sous l'hypothèse $N < M$, on remarque que les lignes de $\Gamma$ admettent une base d'au plus $N$ éléments. Cela est bien inférieur aux $M$ lignes de $\Gamma$. Ainsi il suffit d'accéder à une matrice $\tilde{\Gamma} \in \mathbb{R}^{N \times N}$ composée de $N$ lignes indépendantes de $\Gamma$ qui aura le même noyau que $\Gamma$ (i.e. $\text{ker } \tilde{\Gamma} = \text{ker } \Gamma$). La manière d'obtenir cette matrice $\tilde{\Gamma}$ est de multiplier $\Gamma$ à gauche par une matrice convenable $G \in \mathbb{R}^{N \times M}$. On choisit $G$ gaussienne parce que cela revient à une multiplication à gauche par une matrice inversible presque sûrement.

L'algorithme 2 prend en compte toutes ces idées. Il possède la propriété suivante.

**Proposition 3.1.** *Presque sûrement l'algorithme 2 renvoie une base du noyau ker $\Gamma$.*

### **Algorithmme 2 : Calcul du noyau**

1. Tirer une matrice gaussienne aléatoire $G \in \mathbb{R}^{N \times M}$.

2. Évaluer le produit $Y \stackrel{\text{def}}{=} G\Gamma$.

3. Calculer la SVD de $Y = U\Sigma V^T$.

4. Renvoyer les vecteurs de $V$ associés au noyau de $Y$.

# 4 Expériences numériques

Dans cette section, nous appliquons l'algorithme 1 sur deux exemples illustratifs.

## 4.1 Problème inverse bilinéaire

Considérons le cas d'une convolution entre deux signaux (si-gnal 1 et signal 2) construits comme combinaisons linéaires de bases de fonctions $(f_n)_{1 \le n \le N_0}$ et $(g_n)_{1 \le n \le N - N_0}$.

$$ (\Phi(x))_m \stackrel{\text{def}}{=} \underbrace{\left(\sum_{n=1}^{N_0} x_n f_n\right)}_{\text{signal 1}} * \underbrace{\left(\sum_{n=1}^{N-N_0} x_{N_0+n} g_n\right)}_{\text{signal 2}} (m/M) \quad (9) $$

Ici $0 < N_0 < N \le M$ sont des entiers. Les illustrations numériques sont faites avec $N_0 = 2$, $N = 20$, et $M = 1000$. Le signal 1 est une combinaison linéaire de deux gaussiennes centrées de variances différentes, $f_1$ et $f_2$. Le signal 2 est une combinaison linéaire de $N - N_0$ gaussiennes de même variance mais de moyennes différentes $(g_n)_{n \le N - N_0}$.

L'algorithme d'identification nous donne ici une variété de dimension $D = \text{dim}(\text{ker}(\Gamma)) = 1$. Il n'y a donc qu'une direction de non-identifiabilité. Les premières itérations de l'algorithme 1 sont représentées sur la figure 3. On y voit à la fois le signal 1, le signal 2 et le signal convolé $\Phi(x)$. On observe que le chemin choisi multiplie le signal 1 et le signal 2 par des facteurs inversement proportionnels, tout en conservant le si-gnal convolé. Ainsi, l'algorithme fournit une représentation de la non-identifiabilité qui est l'invariance par renormalisation, due à la structure bilinéaire du problème ($f * g = \lambda f * \frac{1}{\lambda} g$ pour tout $\lambda > 0$). Un résultat de M. Kech et F. Krahmer [2] indique que cette direction de non-identifiabilité est unique génériquement sous la condition ($M \ge 2N - 4$). Ce résultat est bien retrouvé par notre algorithme.

## 4.2 Système dynamique

Pour $y_0, v_0 \in \mathbb{R}^2$ fixés, considérons le système dynamique suivant

$$ \begin{cases} y(0) = y_0 \\ \dot{y}(0) = v_0 \\ \ddot{y}(t) = f(y(t), x), \quad t \in [0, 1] \end{cases} \quad (10) $$

avec $f: \mathbb{R}^{2+N} \to \mathbb{R}^2$ une fonction donnée. Ce système peut être discrétisé. Ceci permet de définir récursivement une application $\Phi(x) \stackrel{\text{def}}{=} (y_m)_{1 \le m \le M/2}$ à partir de la suite

$$ \left\{ \begin{array}{l} y_0, y_1 \in \mathbb{R}^2 \\ y_{m+2} = 2y_{m-1} - y_m + \tau f(y_m, x), \quad \forall 1 \le m \le M/2 \end{array} \right. \quad (11) $$

où $\tau > 0$ est fixé.
<figure><img src="image_3.png" /><figcaption>FIGURE 3 : Non-identifiabilité bilinéaire. Résultat des trois premières itérations de l’algorithme 1 (avec $z_k = 2.5k$) pour le problème bilinéaire (9).</figcaption></figure>

Ce système permet, par exemple, de modéliser l'interaction d'un objet (e.g. satellite) avec un système planétaire. Considérons $I$ planètes de masse $m_i$ et de position $u^i \in \mathbb{R}^2$. Un satellite à la position $y(t) \in \mathbb{R}^2$ et de masse $m$ se déplaçant à proximité de ce système vérifie l'équation Newtonienne

$$ \ddot{y}(t) = -\sum_{i=1}^{I} \frac{G m_i}{\|y(t) - u^i\|_2^3} (y(t) - u^i) \stackrel{\text{def}}{=} f(y(t), x). \quad (12) $$

On se ramène bien à l'équation (10) avec le paramètre $x = (u_1^1, \dots, u_1^I, u_2^1, \dots, u_2^I, m_1, \dots, m_I) \in \mathbb{R}^{3I}$ et $N = 3I$. Dans ce contexte, on s'intéresse à la question suivante : peut-on retrouver la position et la masse des planètes $x$ à partir de la trajectoire discrétisée $\Phi(x)$ ?

Dans l'exemple numérique ci-après, nous avons choisi $M = 200$, $\bar{x} = (0.4, 0.4, 0.6, 0.35, 0.4, 0.8, 0.25, 0.25, 1)$, $I = 3$, $G = 1/2$, $y_0 = (0.3, 0.6)$, $y_1 = (0.31, 0.6)$ et $\tau = 0.0001$. Dans ce cas, l'algorithme indique que le système est identifiable, cependant, 2 valeurs propres de $\Gamma$ sont très petites et on peut considérer que la dimension de la variété $M_{\bar{x}}$ est $D = 2$.

La figure 4 représente les différentes configurations planétaires compatibles avec la trajectoire observée (une variété de dimension 2 discrétisée en 5 × 5 points). On observe que les deux planètes éloignées (orange et rouge) fluctuent à la fois en masse et en position alors que la dernière planète (verte) ne varie pas. D'après la représentation de $M_{\bar{x}}$ obtenue, on peut donner l'interprétation que les deux planètes éloignées par rapport au satellite ont un rôle plus négligeable sur la trajectoire que la planète autour de laquelle le satellite tourne. L'algorithme a ainsi permis de donner une visualisation de non-identifiabilité numérique dans cette modélisation.

# 5 Conclusions et remerciements

Dans ce travail préliminaire, nous avons évalué s'il était possible de définir des algorithmes permettant d'identifier des variétés de non-identifiabilité dans le contexte des problèmes inverses. Nous avons montré qu'un mélange de géométrie différentielle, de différentiation automatique, d'algèbre linéaire randomisée et d'optimisation, fournissent un outil performant

<figure><img src="image_4.png" /><figcaption>FIGURE 4 : Système planétaire. Représentation du résultat de l'algorithme 1 pour le problème (11) avec $(z_k)_{k \le 25}$ une grille de $[-0.5, 0.5]^2$ de pas 0.25. Chaque image correspond à une configuration $x_k$ dans la variété $M_{\bar{x}}$ où sont affichés la trajectoire $(y_m)_{m \le M/2}$ ainsi que les $I$ planètes de taille proportionnelle à leur masse.</figcaption></figure>

pour cette tâche. Le cadre est pour le moment restreint à des variétés de dimension relativement faible. Nous ne l'avons exploité que pour des problèmes jouets. La suite de ce travail consistera à appliquer ces idées sur des problèmes inverses aveugles difficiles et réels. Nous souhaitons également étendre le cadre de cet algorithme au cas de mesures bruitées, en étudiant la stabilité des minimiseurs de fonctions quadratiques de la forme $\frac{1}{2}\|\Phi(x) - y\|_2^2$. Ce travail a bénéficié d'une aide de l'École Universitaire de Recherche MINT (ANR-18-EURE-0023).

# Références

[1] Nathan HALKO, Per-Gunnar MARTINSSON et Joel A TROPP : Finding structure with randomness : Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM review*, 53(2):217–288, 2011.

[2] Michael KECH et Felix KRAHMER : Optimal injectivity conditions for bilinear inverse problems with applications to identifiability of deconvolution problems. *SIAM Journal on Applied Algebra and Geometry*, 1(1):20–37, 2017.

[3] Tien-Yien LI et Zhonggang ZENG : A rank-revealing method with updating, downdating, and applications. *SIAM Journal on Matrix Analysis and Applications*, 26(4):918–946, 2005.

[4] Werner C RHEINBOLDT : Manpak : A set of algorithms for computations on implicitly defined manifolds. *Computers & Mathematics with Applications*, 32(12):15–28, 1996.

[5] Michael SPIVAK : A comprehensive introduction to differential geometry, volume 4. Publish or Perish, Incorporated, 1970.