Pseudo-inverse

[TheReveller]

Bonjour,

Je cherche à résoudre numériquement un système sur-déterminé de 3 inconnues à l'aide d'un nombre d'équations fixé au départ, mais pouvant être de 3 à 100 équations.

J'effectue la résolution de façon matricielle, grâce à A*x = b, alors x = (A+)*b où (A+) est le pseudo-inverse de A.

Pour donner une idée de l'application, "b" n'est pas un vecteur, mais une matrice puisque je résous plusieurs systèmes en même. Donc, par exemple, si je suis sur-déterminé par 10 équations et que j'ai 1000 vecteurs "b" à résoudre, alors je fais [3x10]*[10x1000] où [3x10] est le pseudo-inverse de A qui était alors [10x3], tandis que [10x1000] sont mes 1000 vecteurs à résoudre et donc j'obtiens un résultat [3x1000] qui sont chacune des 1000 résolutions de 3 inconnues, par moindre carré.

Or, parfois dans mes 1000 vecteurs [10x1], je sais d'avance que certaines valeurs sont erronées, alors j'aimerais les invalider. Par exemple, si je sais qu'il y a une valeur invalide, alors j'aimerais d'une certaine façon l'invalider pour que le résultat soit celui de [3x9]*[9x1]. Je peux donc invalider jusqu'à 7 valeurs par vecteur dans cet exemple. Cependant, les valeurs à invalider sont différentes d'un vecteur à l'autre, alors je me demandais s'il était possible de garder cette technique de résolution matricielle plutôt simple tout en invalidant ou compensant l'invalidité de certaines valeurs.

Merci,

Éric

[Ben314]

Salut,
Je comprend pas grand chose à ton énoncé (pas assez "formel" pour ma petite tête).

Je reprend ce que je crois comprendre et tu rectifiera les endroits où je me goure.
Tu as 3 inconnues que tu visualise sous forme d'un vecteur colonne 1x3 nommé "x".
Tu as un certain nombre n (de 3 à 100) équations linéaires en les trois inconnues que tu écrit sous la forme Ax=b où A est une matrice nx3 et b un vecteur colonne 1xn.
C'est là le premier truc que je comprend pas, c'est la notion de "le pseudo inverse de A" : je connais quelques définitions (différentes) de la notion de "pseudo inverse" mais les "pseudo inverses" sont rarement unique. Est ce que tu te place sur une structure euclidienne avec comme supplémentaires "naturel" du noyau et de l'image de A les orthogonaux respectifs de ces dernier («pseudo-inverse de Moore-Penrose») ?

Pour la suite, je crois comprendre le début, où en fait tu as (par exemple) 1000 systèmes de ce type à résoudre donc avec 1000 bi différents donc 1000 solutions xi qu'on peut écrire sous la forme AX=B où l'inconnue X est maintenant une matrice 1000x3 et B une 1000xn.
Sauf que vu ce qu'il y avait dans le premier point et vu que tu parle de 1000 solutions et pas de une solution, j'aurais parié que tu résolvait ce système à l'aide de la pseudo inverse de A alors que tu parle de "moindres carrés"...

donc, a partir de là,... je décroche...

[TheReveller]

Le pseudo-inverse de Moore-Penrose offre des solutions qui minimisent le carré de l'erreur, non ?

Je vais y aller plus mathématique et visuel.

Disons que j'ai

x + y = 3
x - y = 1

et

x + y = 6
x - y = 2

Alors, pas besoin de pseudo-inverse pour ce cas, mais bon :

pinv([1,1;1,-1])*[3;1] = [2;1], soit les valeurs de x et y
pinv([1,1;1,-1])*[6;2] = [4;2], soit les valeurs de x et y

ou, tout simplement :

pinv([1,1;1,-1])*[3,6;1,2] = [2,4;1,2], soit les valeurs de x et y pour les deux cas

J'aurais pu aussi être sur-déterminé :

x + y = 3
x - y = 1
x + 2y = 4

pinv([1,1;1,-1;1,2])*[3;1;4] = [2;1], soit les valeurs de x et y

ou un peu faussé :

x + y = 3
x - y = 1
x + 2y = 3.5

pinv([1,1;1,-1;1,2])*[3;1;3.5] = [1.929;0.857], soit les valeurs de x et y

ce qui minimise l'erreur, puisque :

[1,1;1,-1;1,2]*[1.929;0.857] = [2.786;1.071;3.643] (vs [3;1;3.5] impossible à résoudre)

Si je reviens sur ce cas, avec deux résolutions sur-déterminées à effectuer :

x + y = 3
x - y = 1
x + 2y = 4

et

x + y = 6
x - y = 2
x + 2y = 8

pinv([1,1;1,-1;1,2])*[3,6;1,2;4,8] = [2,4;1,2]

Disons que je modifie les valeurs observées de certains cas par des valeurs erronées :

x + y = 3
x - y = 1
x + 2y = 0 [erroné]

et

x + y = 0 [erroné]
x - y = 2
x + 2y = 8

Alors, je suis très erroné :

pinv([1,1;1,-1;1,2])*[3,0;1,2;0,8] = [1.429,2.286;-0.143,1.571]

Y aurait-il une façon de compenser les cas erronés, sachant qu'ils sont erronés ?

Là, on pourrait dire de tout simplement ne pas les utiliser puisqu'on sait qu'ils sont erronnés !

Sauf que le problème est que j'ai par exemple 20 équations de 3 inconnues (au lieu de 3 équations de 2 inconnues dans l'exemple) et par exemple 1000000 vecteurs de résultats (au lieu de 2 dans l'exemple) et je veux éviter d'effectuer une boucle personnalisée pour chacun des 1000000 vecteurs puisque normalement je peux obtenir le résultat directement de façon matricielle. Bref, j'aimerais ne pas avoir à effectuer 1000000 itérations pour résoudre vecteur par vecteur en fonction des données erronées.

Ce que je veux éviter, si je reprends l'exemple plus haut :

Boucler les vecteurs

Vecteur 1 :

x + y = 3
x - y = 1
x + 2y = 0 [erroné]

Alors, calculer pinv([1,1;1,-1])*[3;1] = [2;1]

Vecteur 2 :

x + y = 0 [erroné]
x - y = 2
x + 2y = 8

Alors, calculer pinv([1,-1;1,2])*[2;8] = [4;2]

[...]

Vecteur 1000000
[...]

Alors que s'il n'y avait pas de valeurs erronées, je ne bouclerais pas, ce serait directement : pinv([20x3])*[20x1000000] = [3x1000000], les 1000000 vecteurs résolus pour les 3 inconnues.

[TheReveller]

Je vais essayer de simplifier ce gros message, je crois qu'il mélange encore plus qu'autre chose.

Commençons par définir par exemple un système de 10 équations, 3 inconnues. Ces 10 équations, 3 inconnues, c'est ce qui est à gauche du signe "=". Donc ce ne sont pas vraiment des équations à résoudre tant qu'on a pas les valeurs à mettre à droit du signe "=", mais c'est ce que définit ma matrice à pseudo-inverser, soit une matrice [10x3] qui devient [3x10] après le pseudo-inverse.

Maintenant, définissions un "vecteur d'observations" comme étant les 10 observations correspondant à ce qui est à droite du signe "=". La résolution consiste alors à multiplier le pseudo-inverse [3x10] par ce vecteur d'observations [10x1] afin d'obtenir les valeurs des 3 inconnues ([3x10] * [10x1] = [3x1]).

Petite parenthèse, disons que je sais que j'ai 4 valeurs aberrantes dans ce vecteur d'observations, alors je n'aurais qu'à les enlever et enlever les équations correspondantes avant de résoudre, ce qui donnerait [3x6] * [6x1] = [3x1].

Revenons à la résolution des 10 équations, 3 inconnues et 1 vecteur de 10 observations. Disons que j'ai toujours les mêmes 10 équations, 3 inconnues (ce qui est à gauche du signe "="), mais que j'ai un autre vecteur de 10 observations à résoudre (ce qui est à droite du signe "="). J'ai donc 2 vecteurs de 10 observations à résoudre, mais ces deux vecteurs ne sont aucunement en lien entre eux. La résolution des 3 inconnues pour le premier vecteur donnera certains résultats, tandis que la résolution des 3 inconnues pour le second vecteur donnera certains résultats et c'est parfait ainsi puisqu'il y a strictement aucun lien entre les valeurs obtenues des 3 inconnues au premier vecteur et les valeurs obtenues des 3 inconnues au second vecteur. Lorsque je résous, j'effectue [3x10] * [10x2] = [3x2], donc les deux vecteurs des valeurs obtenues pour les inconnues lors de la résolution, le premier vecteur correspond aux valeurs des inconnues pour le premier vecteur d'observations et le second vecteur correspond aux valeurs des inconnues pour le second vecteur d'observations. Tout est parfait ainsi. Si j'ai 1000 vecteurs d'observations, alors ce sera [3x10] * [10x1000] = [3x1000] et tout va bien.

Maintenant, si je reviens à la parenthèse où j'enlève tout simplement des équations pour résoudre un vecteur unique en calculant plutôt [3x6] * [6x1] = [3x1]. Si j'avais eu un deuxième vecteur d'observations et que celui-ci n'a aucune observation aberrante, alors je voudrais calculer [3x10] * [10x1] = [3x1].

Or, je ne peux plus combiner cela en un seul calcul matriciel [3x10] * [10x2] = [3x2] puisque je veux enlever 4 observations aberrantes du premier vecteur! Je dois donc calculer [3x6] * [6x1] = [3x1] pour le premier cas, [3x10] * [10x1] = [3x1] pour le second cas, et ainsi de suite si j'ai 1000 cas. Bref, itérativement au lieu d'y aller tout simplement de façon matricielle [3x10] * [10x1000] = [3x1000].

Ma question est donc : est-ce possible d'utiliser une astuce pour conserver le plus possible l'avantage du calcul matriciel tout en compensant au mieux les données aberrantes que je voudrais invalider pour chaque cas ?

On s'entend ici que je parle de résolution numérique, je suis sur-déterminé à 10 équations, 3 inconnues puisque la résolution n'est pas mathématiquement exacte.

Merci de m'avoir lu!