Méthode de Gauss et pivotage
Réponses
-
C'est difficile à dire, car je ne sais pas comment fonctionne cette fonction testmatrix.
-
Bonjour,
Avec Matlab, en échelle logarithmique:nmax=30; tab=[]; for n=4:nmax H=hilb(n); K=invhilb(n); tab=[tab norm(H*K-eye(n))]; end plot(log(tab))
Pour $20$, Matlab donne $ 2.5470\times 10^{11}$, et Octave $1.5302\times 10^{11}$, dans les deux sens.
Cordialement,
Rescassol -
Bonjour,
Maintenant qu'il a été confirmé qu'un certain bavard n'était pas le praticien qu'il prétendait être, il reste à revenir à la problématique posée par le conditionnement d'une matrice. Et à éviter de noyer cette problématique sous un tombereau de calculs inutiles.
Cette problématique n'a rien à voir avec le pivot de Gauss, et se résume à la trivialité suivante: si vous décidez de faire disparaître une partie de l'information dont vous disposez, eh bien vous ne disposez plus de l'information que vous avez fait disparaître.
Les tournures de phrases employées par plusieurs intervenants semblent indiquer une mauvaise compréhension du théorème de Cramer. La propriété "si le déterminant est non nul, tout va bien" ne s'applique qu'à un anneau intègre (et se comprend encore mieux par plongement dans le corps des fractions de cet anneau).
Mais, pas de chance, l'ensemble des nombres flottants pour une précision donnée ne forme pas un anneau intègre. Au contraire, la multiplication par une matrice suivie d'une troncature plus grossière que le conditionnement se comporte presque partout à la façon d'une projection. Et le relèvement est donc non pas un point isolé, mais quelque chose qui ressemble à une fibre et qui engendre un espace vectoriel non trivial.
Pour illustrer cela, il est inutile de prendre $n=20$. En se limitant à $n=5$ et en utilisant la matrice de Hilbert et le vecteur $V=\left[1,1,1,1,1\right]$, on obtient:# vp cut=20 cut=10 cut=4 cut=2 cut=1 a0 1.567051E+0 3.131728E+0 3.131728E+0 3.131703E+0 3.128019E+0 3.163459E+0 a1 2.085342E-1 -2.026608E-1 -2.026608E-1 -2.026736E-1 -2.033177E-1 -1.967521E-1 a2 1.140749E-2 2.788283E-3 2.788283E-3 2.748872E-3 3.485050E-3 -5.170851E-2 a3 3.058980E-4 -1.307118E-5 -1.307114E-5 2.788103E-5 -5.565119E-3 1.295874E-2 a4 3.287929E-6 1.545040E-8 1.547701E-8 2.663011E-5 3.646223E-3 -3.115940E-2
La colonne vp donne les valeurs propres, triées par ordre décroissant. La colonne cut=xx donne les coefficients (dans la base des vecteurs propres) de la troncature d'ordre xx du vecteur $C=M\cdot V$. De ce point de vue, couper à 20 décimales ou seulement à 10 ne change pas grand' chose. Pour les autres troncatures, les valeurs exactes des coefficients $a_{k}$ sont de plus en plus masquées par le bruit engendré par la troncature.
Lorsque l'on calcule $V=M^{-1}\cdot C$, le bruit sur le coefficient $a_{k}$ est amplifié par le facteur $1/\mu_{k}$ et devient dominant. Encore une fois, ceci n'a rien à voir avec la méthode avec laquelle on calcule $M^{-1}\cdot C$. Dans cet exemple, le calcul exact de $M^{-1}$ est immédiat (y compris pour n=20) et donne une matrice à coefficients entiers. Utiliser la valeur exacte de cette matrice ne changerait rien au fait que l'information détruite par la troncature de $C$ est perdue, c'est à dire définitivement perdue.
Anders gesagt: en calcul numérique, Cramer = pipeau.
Cordialement, Pierre. -
S'agissant de la résolution d'un système linéaire $HX=Y$ où $H$ est une matrice
de Hilbert, donc à coefficients rationnels, et en prenant également un second membre
rationnel, jusqu'à quelle taille de matrice peut-on espérer résoudre de manière exacte
(avec le pivot de Gauss) en ne manipulant que des rationnels ?
J'imagine que les entiers apparaissant dans les fractions vont grossir très vite. Disons en
python/numpy avec un type entier non limité en taille ? -
Si on veut faire les calculs seulement avec des rationnels, on peut appliquer la méthode du "mediant rounding" (désolé je ne connais pas son nom français) sur la fraction pour empêcher que le numérateur et dénominateur deviennent trop grands. Avec cette méthode, les erreurs d'arrondis tendent à s'annuler et on peut souvent obtenir un résultat exact. Donc en python, je suppose qu'il serait possible d'écrire une routine pour arrondir la fraction comme ca quand elle dépasse une certaine limite, en bits par exemple.
Il y a d'ailleurs une librairie écrite en C qui s'appelle Miracl qui a un type de rationnel appelé "flash" pour floating-slash qui implémente le mediant rounding. Je l'ai utilisée, juste pour voir, pour résoudre le système avec le membre de droite proposé par Remarque et la matrice 20x20, et j'ai obtenu un résultat exact. D'ailleurs l'auteur de la librairie a justement utilisé cette histoire de matrice de Hilbert pour démontrer l'avantage de cette méthode d'arrondi; il fournit un programme en exemple qui utilise le pivot de Gauss pour résoudre le système, et ca fonctionne assez bien. -
Bonjour,
La façon dont remarque a repris les calculs de dlzlogic ne semble pas avoir clarifié la problématique du conditionnement. Dans sa feuille de calcul, dlzlogic part de la colonne B et arrive à la colonne A avec:B A xdlz 17.35464129523727 114240678.63930499553680419921875 17.360695 15.61837349956545 -17390761859.52415084838867187500000 15.540413 14.29712546673948 649015070892.48754882812500000000000 14.297478 13.22950062231931 -10296415651124.05078125000000000000000 13.206627 12.33687577789913 85224105961665.78125000000000000000000 12.394733 11.57348170270973 -401003842368058.93750000000000000000000 11.595707 10.90980272723542 1082763315074396.62500000000000000000000 10.942120 10.32548883112620 -1514915433455079.25000000000000000000000 10.264938 9.80583010743077 571551479894145.37500000000000000000000 9.827183 9.33981123047863 650106297378619.62500000000000000000000 9.379623 8.91895364384907 948856011170290.87500000000000000000000 8.900822 8.53658579328990 -3468427722440542.50000000000000000000000 8.552750 8.18736188212466 1760599110935123.75000000000000000000000 8.195312 7.86693270563651 -60503842414165.94531250000000000000000 7.864983 7.57171361318198 3787951888378411.00000000000000000000000 7.523705 7.29871674294967 -3361486357068956.00000000000000000000000 7.276920 7.04542707992456 -6515507344182813.00000000000000000000000 7.020969 6.80970916821734 12260533635134778.00000000000000000000000 6.787663 6.58973574593881 -7417184649625362.00000000000000000000000 6.551937 6.20000000000000 1601108094320109.00000000000000000000000 6.192085
La colonne xdlz n'est absolument pas la valeur proposée comme solution de $H\cdot X=B$, mais au contraire la valeur de contrôle $H\cdot A$, conduisant à $\left|H\cdot A-B\right|=0.152$. Les matrices que remarque nomme Hdlz et bdlz ne sont pas les objets internes, mais un affichage formaté à trois décimales. Probablement parce que 20 chiffres fois 21 colonnes, cela ferait cinq fois les 80 caractères disponibles dans une page d'affichage !
J'ai refait le calcul de $\left|H\cdot A-B\right|$ avec 40 décimales et je trouve à peu près la même chose que dlzlogic, soit $\left|H\cdot A-B\right|\approx0.231$. Le fait que $\left|A\right|\approx1.72E16$ confirme que le calcul est relatif à la matrice $H$ et pas à la matrice tronquée Hdlz. Cela dit, il est particulièrement instructif de comparer le résidu $\left|H\cdot A-B\right|\approx0.152$ trouvé par dlzlogic avec le résidu $\left|Hdlz\cdot X-bdlz\right|<10^{-14}$ trouvé par remarque.-
$\,$
- Obtenir un $A$ tel que $\left|Hdlz\cdot A-bdlz\right|<10^{-14}$ est une banalité banale que l'on peut atteindre en procédant à peu près n'importe comment, par pivot, par calcul d'inverse, par un logiciel payant, par un logiciel open source, par gradient conjugué, par algorithme génétique et même par la méthode miraculeuse chère au lecteur. Et l'on trouve $\left|A\right|<400$. Et si l'on veut obtenir un $A^{+}$ tel que $\left|Hdlz\cdot A^{+}-bdlz\right|<10^{-28}$, il suffit de régler RealField sur $30\times\ln10/\ln2$ et de faire clic. On trouvera à nouveau que $\left|A^{+}\right|<400$. Et même que $\left|A-A^{+}\right|<10^{-12}$.
$\,$ - Pourquoi se trouve-t-on dans la banalité banale d'une bijection bicontinue et sans surprise lorsque l'on veut résoudre $\left|Hdlz\cdot A-bdlz\right|<\varepsilon$ ? C'est tout banalement parce que le conditionnement de cette matrice est tout petit. On trouve en effet que $\left\Vert Hdlz\right\Vert _{1}\times\left\Vert Hdlz^{-1}\right\Vert _{1}\approx55150$.
$\,$ - Obtenir un $A$ tel que $\left|H\cdot A-B\right|<0.16$ est autrement difficile. En effet le conditionnement de la matrice $H=Hilbert(20)\in\mathcal{M}_{20}\left(\mathbb{Q}\right)$ est $6.28E28$. VINGT-NEUF chiffres devant la virgule. Pour ce qui est de la continuité, il y a continuité. Il y a même un facteur de Lipschitz. Il vaut $6.28E28$. Cela peut se comparer avec le nombre d'Avogadro: modifier $B$ d'un atome peut conduire à modifier $A$ d'une planète entière. Cela s'appelle un problème conditionné de façon explosive.
$\,$ - On en arrive au bilan demandé par Sylviel. Pour faire court, dlzlogic n'a tout simplement pas remarqué que son résultat $\left|H\cdot A-B\right|<0.16$ a été obtenu en introduisant un vecteur $A$ qui est énorme, et pour tout dire hors norme. On a en effet $\left|A\right|\approx1E16$. La problématique n'est donc pas de faire diminuer la norme du résidu $X\mapsto H\cdot X-B$. Elle est que, à part le cas où $B$ appartiendrait à un certain ensemble, qui est de mesure évanescente, il est impossible de faire diminuer la norme de $\left|M\cdot A-B\right|$ sans faire exploser la norme de $\left|A\right|$. Sur l'exemple, un calcul à 50 décimales donne au moins 20 chiffres significatifs pour $stupid=H^{-1}\cdot B$, et permet de vérifier que $\left|stupid\right|=4.5E24$. Si le problème consistait à compter des boulons, cela fait beaucoup.
Autre critique de la feuille de calcul produite par dlzlogic: ne pas avoir remarqué et insisté sur le fait que ce calcul avec $x_{20}=6.20$ était bien plus représentatif de la situation de mauvais conditionnement que l'exemple trafiqué qui lui avait été proposé (avec $x_1=20.$). Dans ce deuxième exemple en effet, l'information que les $x_n$ sont vraisemblablement entiers est une information énorme, qui permet largement de reconstituer l'information perdue par troncature.
Proverbe informatique: quand c'est perdu, c'est perdu (et les données perdues ne seront pas redonnées).
Cordialement, Pierre - Obtenir un $A$ tel que $\left|Hdlz\cdot A-bdlz\right|<10^{-14}$ est une banalité banale que l'on peut atteindre en procédant à peu près n'importe comment, par pivot, par calcul d'inverse, par un logiciel payant, par un logiciel open source, par gradient conjugué, par algorithme génétique et même par la méthode miraculeuse chère au lecteur. Et l'on trouve $\left|A\right|<400$. Et si l'on veut obtenir un $A^{+}$ tel que $\left|Hdlz\cdot A^{+}-bdlz\right|<10^{-28}$, il suffit de régler RealField sur $30\times\ln10/\ln2$ et de faire clic. On trouvera à nouveau que $\left|A^{+}\right|<400$. Et même que $\left|A-A^{+}\right|<10^{-12}$.
-
pldx1 écrivait:
> Bonjour,
>
>
>
> La façon dont remarque a repris les calculs de dlzlogic ne semble pas avoir clarifié la problématique du conditionnement.
> La colonne xdlz n'est absolument pas la valeur proposée comme solution de $H\cdot X=B$
Je salue tes talents d'extralucide dans le déchiffrage du dlzlogic dans le texte. Mais bon, flogging a dead horse et tout ça. -
Bonjour,
Il a été suffisemment question de mes calculs pour que je sois en droit de répondre.
D'abord, ce problème de résolution de système linéaire est la conclusion obligatoire de l'utilisation de la méthode des moindres carrés. C'est donc un calcul classique pour moi. Par contre, il est vrai que je ne l'utilise pas dans un cadre théorique et avec des données particulières, mais dans un contexte de calcul numérique réel.
J'ai lu quelques sujets traitant de cela. Celui-ci m'a semblé le mieux correspondre à la question d'origine.
https://www.ceremade.dauphine.fr/~viossat/pdfs/algebre1/2009-10/systemes-lineaires-nov09.pdf
Le titre du chapitre 5 est "Matrices et systèmes linéaires". Le terme "matrice" n'apparait pas dans les chapitres précédents, par contre, le terme "tableau" est utilisé.
Voici donc, ce que je réponds aux différents intervenants.
1- tous ces échanges sont complètement hors sujet. C'est vrai, ils sont intéressants, mais d'après ce que j'ai compris ne concernent que la théorie des nombres, par ailleurs je me demande si à un ou deux endroits, on ne mélange pas "valeur exacte" et "valeur réelle".
2- en calcul informatique, sauf s'il s'agit d'entiers, une valeur exacte n'existe pas : un flottant est représenté par une mantisse et une caractéristique.
3- il est bien évident que j'ai fait les calculs avec des valeurs fournies numériques exactes c'est à dire de la forme 1/(i+j+1). De la même façon les valeurs du membre de droite données par Remarque (qui voudra bien me pardonner de ne pas avoir deviné que le nombre 20. était le premier d'une liste avec une vingtaine de décimales. Personnellement j'ai l'habitude de formater tous les résultats de la même façon. Mais j'ai bien compris que tout écart était inexcusable).
4- J'ai l'impression que Pldx a remarqué que les valeurs du membre de droite de mon premier calcul étaient autrement "tordu". Ca ne m'étonne pas, cette liste m'a été donnée par un membre qui voulait démontrer que j'avais tort. Ce système donnait des résultats aberrants avec Scilab.
5- Le sujet de ce fil pose la question de résolution de système linéaire. Je n'ai pas vu qu'il était question de résolution de systèmes "à la limite". On a évoqué le fait qu'il existait d'autres méthodes de résolution de systèmes linéaires, apparemment personne n'en a parlé. Dans le même ordre d'idée il semble que l'on ne parle que de résoudre "A . X = B". C'est à dire que le tableau des coefficients des inconnues est considéré comme une matrice et que on réalise le calcul matriciel "X = B . A-1" (ou quelque-chose comme ça). Le cours mis en lien précise d'ailleurs que la méthode du pivot est utilisée pour réaliser les opérations sur les matrices. J'ai écrit en C les fonctions qui réalisent cela. J'ai fait des comparatifs, temps et précision, avec la méthode de Gauss. Il me semble que cette question pourrait être évoquée.
6- On peut imaginer que le demandeur cherche à écrire un programme qui, entre autres, résout des systèmes linéaires. Il veut faire des tests et utilise un outil qu'il a sous la main : Matlab.
Il y a quelques temps, j'ai écrit ce papier pour essayer de clarifier les choses :
http://www.dlzlogic.com/aides/Pivot_Gauss.pdf
En effet des exercices souvent proposés contiennent des nombres entiers comme paramètres. Il est dommage que des professeurs, membres actifs de certains forum, donnent des conseils tels que* des étudiants se focalisent sur des recherches de simplification numériques, alors que justement l'un des intérêts de cette méthode est de pouvoir être mise en oeuvre automatiquement.
*[complément en italique rajouté] -
Bonjour
dlzlogic écrivait:
> D'abord, ce problème de résolution de système linéaire est la conclusion obligatoire de
> l'utilisation de la méthode des moindres carrés.
Hum, obligatoire... quand il s'agit d'appliquer la méthode des moindres carrés à une régression linéaire !
Dans d'autres circonstances (par exemple la détermination de période), on obtient un système qui n'est pas linéaire.
> C'est donc un calcul classique pour moi. Par contre, il est vrai que je ne l'utilise pas dans
> un cadre théorique et avec des données particulières, mais dans un contexte de calcul numérique réel.
Comment être certain que le calcul réel n'est pas (ne sera jamais) avec des données particulières ?
Même si cela arrive dans 1% des cas, c'est quand même dommage quand cela se traduit
par une explosion de fusée ou un effondrement de pont...
Donc toutes les situations, particulières ou pas, à la limite ou pas, sont à étudier à mon avis.
> 1- tous ces échanges sont complètement hors sujet.
Les intervenants n'ont fait que répondre à ta question ! je te rappelle :
dlzlogic écrivait:
> Bonsoir Remarque,
> Pourrais-tu donner un exemple d'un tel système ?
> J'avoue que j'ai un peu de mal à comprendre.
Et c'est ainsi que tu les remercies d'avoir pris de leur temps pour t'expliquer certaines choses,
temps qu'ils avaient déjà pris dans une autre discussion...
> je me demande si à un ou deux endroits, on ne mélange pas "valeur exacte" et "valeur réelle".
Toi qui ne mélange pas, quelles sont les définitions mathématiques de ces deux notions ?
> ...(qui voudra bien me pardonner de ne pas avoir deviné que le nombre 20. était le
> premier d'une liste avec une vingtaine de décimales. Personnellement j'ai l'habitude de
> formater tous les résultats de la même façon. Mais j'ai bien compris que tout écart était inexcusable).
Quand on ne sait pas lire, quand on n'essaie pas d'être ouvert, on ne voit pas les 20 lignes
qui étaient pourtant bien présentées ! http://www.les-mathematiques.net/phorum/read.php?15,1048261,1049025#msg-1049025
De plus, tes propos sont assez mal venus quand on lit ta prose :
<< Je veux bien faire des simulations, mais j'aime pas trop montrer qu'un membre dit n'importe quoi. >>
Par ailleurs, ta phrase
<< Au passage, tu me diras dans quel contexte tu arrives à des valeurs pour B avec 22 ou 23 chiffres significatifs. >>
montre à quel point tu ne comprends pas l'exemple que l'on t'a présenté...
Même si les données sont avec 22 ou 23 chiffres significatifs, le résultat n'en a même pas 1 !
Toi qui penses des données à 5 chiffres suffisent, tu es loin du comptes dans ce genre de situation.
> On a évoqué le fait qu'il existait d'autres méthodes de résolution de
> systèmes linéaires, apparemment personne n'en a parlé.
Faudrait peut-être relire la discussion, mais en essayant de comprendre les maths qui y sont exposées.
> En effet des exercices souvent proposés contiennent des nombres entiers comme paramètres.
> Il est dommage que des étudiants se focalisent sur des recherches de simplification numériques,
> alors que justement l'un des intérêts de cette méthode est de pouvoir être mise en oeuvre automatiquement
Rappelle-nous : tu n'as jamais été enseignant en math ? ...ni même étudiant en math, hein ? Dommage. -
Salut,
je trouve que la matrice de Hilbert $[ 1 /(i+j) ] $ est un exemple que l'on peut critiquer car son déterminant (même si le conditionnement n'est pas en lien avec le déterminant) est très proche de 0 : il est de l'ordre de 10^{-150} , ce qui est très inférieur aux coefficients $1/(i+j)$. Numériquement, le déterminant est quasi nul, ce qui pourrait faire imaginer pourquoi on ne peut pas inverser le système...
Je propose cette matrice très creuse, à coefficients entier, de déterminant 1 (donc la matrice inverse est à coefficients entiers) :
$ C = \left[ \begin {array}{cccccccccccccccccccc}
1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-1\\
-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7\\
0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7\\
0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&9
\end {array} \right]
$
(diagonale de 1, sauf le 9 en bas à droite,
dernière colonne de 7, sauf le -1 en haut à droite, et toujours le le 9 en bas à droite
sous-diagonale de -8)
Résoudre en $X$ le système $C . X = B$ où $B$ est le vecteur colonne
$[-2, 16, -2, 16, -2, 16, -2, 16, -2, 16, -2, 16, -2, 16, -2, 16, -2, 16, -2, 17]$ (alternance de -2 et 16, sauf le 17 à la fin)
demande de travailler avec une mantisse de minimum 18 ou 19 chiffres en base 10 (si je ne me trompe pas), ce qui est énorme vu la taille des données (...et du résultat...) ! -
Cet exemple est assez marrant : avec scilab, on s'en sort mieux avec la commande inv(C) que par résolution du système linéaire (en utilisant que les coefficients sont entiers (et de la divination, bien sûr)). C'est quasiment une abomination ! :-D
-
Bonsoir,
Octave, avec:M=eye(20)-8*diag(ones(19,1),-1); M(:,20)=7*ones(20,1); M(1,20)=-1; M(20,20)=9; B2=[-2; 16]; B4=[B2; B2]; B8=[B4 ;B4]; B16=[B8; B8]; B=[B16; B4]; B(20)=17; X=M\B ErrX=max(B-M*X) Y=inv(M)*B ErrY=max(B-M*Y) I=M*inv(M)
donne:warning: matrix singular to machine precision, rcond = 2.53297e-021 X = -1.67551 0.32449 -1.67551 0.32449 -1.67551 0.32449 -1.67551 0.32449 -1.67551 0.32449 -1.67551 0.32449 -1.67551 0.32449 -1.67553 0.32433 -1.67682 0.31394 -1.75994 0.32449 ErrX = 1.1369e-013 warning: inverse: matrix singular to machine precision, rcond = 2.53297e-021 Y = 1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 ErrY = 16 warning: inverse: matrix singular to machine precision, rcond = 2.53297e-021 I = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
Cordialement,
Rescassol -
Ton Y n'est vraiment pas loin.* C'est là qu'il faut sortir la boule de cristal !
* Enfin, dans $\Z^{20}$ s'entend. -
Bonsoir,
Oui, ma boule de cristal rajoute un $"-"$ sur le premier terme.
Cordialement,
Rescassol -
Ce qui est encore plus amusant, si je ne me trompe pas, c'est que la méthode de Gauss sans stratégie de pivot donne ici la solution exacte sans aucun problème.
-
Exact.
Cependant une légère perturbation du coefficient 1 en haut à gauche de C en $1 - 0.1^{15}$ donne lieu à des résultats très différents et le simple pivot de Gauss ne suit pas. -
Comme H, je ne fais que passer.
Dlz demandait à un endroit dans cette discussion une autre façon de faire que le pivot : SVD (singular value decomposition). -
@ leon1789 : je n'en doute pas. De même que j'imagine bien que Gauss avec pivot partiel parte dans le décor sur la matrice non perturbée (quelqu'un a testé ?).
-
Avec le pivot partiel, il faut travailler avec une mantisse de 14 chiffres (en base 10) pour obtenir un résultat sensé (avec 2 chiffres significatifs).
Et si on perturbe un tout petit peu le coefficient 1 en haut à gauche, ça délire... -
Ca par exemple ?
Les solutions
-1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -1.00000 1.00000 -
oui, c'est bien les solutions (c'est indiqué au-dessus), on peut les obtenir sans problème par pivot de Gauss sur les entiers (c'est dit ci-dessus), mais pas par pivot partiel avec une mantisse à 6 chiffres...
Si tu veux éviter de travailler avec des entiers, alors divise par exemple par 7 les deux premiers lignes de C et B (le fait de diviser par 7 donne un système mathématiquement équivalent, donc on devrait trouver les mêmes solutions) et tu vas voir comment ça va faire des choses étranges avec pivot de Gauss simple ou partiel.
$C = \left[ \begin {array}{cccccccccccccccccccc}
1/7&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-1/7\\
-{\frac {8}{7}}&1/7&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&1
\\0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&0&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&1&7
\\0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&-8&9
\end {array} \right]$
et
$B = [-\frac27,{\frac {16}{7}},-2,16,-2,16,-2,16,-2,16,-2
,16,-2,16,-2,16,-2,16,-2,17]$. -
En reponse (tardive) a taupin_curieux, il me faut environ 1 seconde pour n=128 en calcul exact (avec second membre identiquement 1), pour le script giac suivant n:=128; H:=hilbert(n):; time(v:=linsolve(H,[1$n])); v; l2norm(H*v-[1$n])
la plus grande coordonnee de v a environ 100 chiffres.
Pour n=256, il faut environ 30 secondes (et 200 chiffres). -
Bonjour,
Travailler avec $n=10$ est bien suffisant pour faire apparaitre la problématique, et cela facilite l'affichage. On part donc de l'équation $C\cdot X=B$, soit: \[ \left(\begin{array}{rrrrrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1\\ -8 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 7\\ 0 & -8 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 7\\ 0 & 0 & -8 & 1 & 0 & 0 & 0 & 0 & 0 & 7\\ 0 & 0 & 0 & -8 & 1 & 0 & 0 & 0 & 0 & 7\\ 0 & 0 & 0 & 0 & -8 & 1 & 0 & 0 & 0 & 7\\ 0 & 0 & 0 & 0 & 0 & -8 & 1 & 0 & 0 & 7\\ 0 & 0 & 0 & 0 & 0 & 0 & -8 & 1 & 0 & 7\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -8 & 1 & 7\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -8 & 9 \end{array}\right)\cdot X=\left(\begin{array}{c} -2\\ 16\\ -2\\ 16\\ -2\\ 16\\ -2\\ 16\\ -2\\ 17 \end{array}\right) \] En supposant que $B,C$ sont connues en précision infinie, la formule $X=B^{-1}\cdot C$ donne \[ X=\left[-1,1,-1,1,-1,1,-1,1,-1,1\right] \] et c'est fini. Si l'on suppose que seule une troncature des données est disponible, il est absurde d'essayer de jouer aux formules de Cramer, le conditionnement est $9E10$ (et $4E20$ pour $n=20$).
Que peut on faire? Pour commencer, poser tous les calculs à venir en précision largement supérieure au phénomène à étudier. $
\def\trc#1{{{\vphantom{#1}}^{t}{#1}}}
$Le plus efficace pour comprendre ce qui se passe est de transformer le problème en $S\cdot X=D$ avec $S\doteq\left(\trc C\cdot C\right)$, $D\doteq\left(\trc C\cdot B\right)$: \[ \left(\begin{array}{rrrrrrrrrr} 65 & -8 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -57\\ -8 & 65 & -8 & 0 & 0 & 0 & 0 & 0 & 0 & -49\\ 0 & -8 & 65 & -8 & 0 & 0 & 0 & 0 & 0 & -49\\ 0 & 0 & -8 & 65 & -8 & 0 & 0 & 0 & 0 & -49\\ 0 & 0 & 0 & -8 & 65 & -8 & 0 & 0 & 0 & -49\\ 0 & 0 & 0 & 0 & -8 & 65 & -8 & 0 & 0 & -49\\ 0 & 0 & 0 & 0 & 0 & -8 & 65 & -8 & 0 & -49\\ 0 & 0 & 0 & 0 & 0 & 0 & -8 & 65 & -8 & -49\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -8 & 65 & -65\\ -57 & -49 & -49 & -49 & -49 & -49 & -49 & -49 & -65 & 474 \end{array}\right)\cdot X=\left(\begin{array}{r} -130\\ 32\\ -130\\ 32\\ -130\\ 32\\ -130\\ 32\\ -138\\ 547 \end{array}\right) \]
La matrice $S$ est symétrique réelle. Si l'on signale par typage que tel est exactement le cas, le logiciel utilisera les algorithmes spécifiques aux matrices symétriques réelles, utilisant des techniques de gradient. On diagonalise $S$, selon $\trc{P\cdot}S\cdot P=\Delta$ (la matrice $P$ est orthogonale) . Les valeurs propres sont: \[ \lambda_{S}\approx\left[5.3044\times10^{-18},52.055,54.955,60.051,64.476,69.938,74.138,77.941,80.182,525.26\right] \] On voit donc que, en présence de troncatures, $S$ se comporte non pas comme un isomorphisme agissant sur les colonnes, mais comme une matrice de rang $n-1$. A nouveau, en calcul numérique, Cramer = pipeau.
Que peut-on faire pour exploiter au mieux les données disponibles ? On a: \[ \trc P\cdot D\doteq\left(\begin{array}{l} -05.7212E-20\\ +00.0017110\\ -03.7916\\ +00.015887\\ -11.294\\ +00.029343\\ -20.616\\ +00.015172\\ -27.276\\ -67.085 \end{array}\right);\phi\trc P\cdot D=\left(\begin{array}{l} +00.00000\\ +00.0017110\\ -03.7916\\ +00.015887\\ -11.294\\ +00.029343\\ -20.616\\ +00.015172\\ -27.276\\ -67.085 \end{array}\right);P\cdot\Delta^{-1}\cdot\phi\trc P\cdot D=\left(\begin{array}{l} -0.98921\\ +1.0108\\ -0.98921\\ +1.0108\\ -0.98921\\ +1.0108\\ -0.98919\\ +1.0110\\ -0.98787\\ +1.0108 \end{array}\right) \] Le cœur de l'affaire est l'opérateur $\phi$ qui supprime la projection sur le vecteur propre affecté de la valeur propre évanescente $\lambda_{0}$ lorsque l'on suspecte que la composante correspondante vient plus du bruit engendré par la troncature que de l'information engendrée par le vecteur inconnu. Cela conduit à
\[ X=\left(\begin{array}{c} -0.9892141755\\ +1.0107858251\\ -0.9892141704\\ +1.0107858656\\ -0.9892138464\\ +1.0107884577\\ -0.9891931095\\ +1.0109543529\\ -0.9878659475\\ +1.0107858244 \end{array}\right)+x\left(\begin{array}{c} 0.31156447255\\ 0.31156448883\\ 0.31156461883\\ 0.31156565878\\ 0.31157397845\\ 0.31164053581\\ 0.31217299462\\ 0.31643266511\\ 0.35051002904\\ 0.31156447026 \end{array}\right) \] où la quantité $x$ est inconnue. Elle représente l'information qui a été détruite par la troncature et qui ne peut donc pas être reconstituée.
Quant à dire "si on perturbe un tout petit peu le coefficient 1 en haut à gauche, ça délire...", quelle bonne plaisanterie ! Remplaçons le $1$ par $1000000001/1000000000$, et lançons une résolution exacte. Cela donne: \[ X^{*}=\frac{1}{1134217729}\,\left(\begin{array}{r} -1000000000\\ 1268435465\\ -999999937\\ 1268435969\\ -999995905\\ 1268468225\\ -999737857\\ 1270532609\\ -983222785\\ 1268435457 \end{array}\right)\neq X=\left(\begin{array}{r} -1\\ 1\\ -1\\ 1\\ -1\\ 1\\ -1\\ 1\\ -1\\ 1 \end{array}\right) \]
Quelle surprise de voir la réponse d'un système mal conditionné être fortement perturbée lorsque ses conditions initiales sont perturbées d'une façon qui est énorme par rapport au conditionnement ! Et quel émerveillement lorsque l'on constate que la valeur normalisée de la surprise, soit \[ \frac{X^{*}-X}{\left|X^{*}-X\right|_{2}} \] coïncide à 11 décimales avec la première colonne de $P$, c'est à dire avec le vecteur propre associé à la valeur propre évanescente. Kilucru !
Cordialement, Pierre.
Redisons, pour les liseurs en diagonale, qu'un affichage avec quelques chiffres significatifs veut dire que les calculs ont été conduits à la précision voulue, et tronqués ensuite, pour la facilité de l'affichage. -
<< Quelle surprise de voir la réponse d'un système mal conditionné être fortement perturbée lorsque ses conditions initiales sont perturbées d'une façon qui est énorme par rapport au conditionnement ! >>
En effet, il n'y a rien d'étonnant quand on comprend la notion de conditionnement.
La perturbation initiale (ou une erreur d'arrondi en cours de calcul, etc.) est énorme devant le condition (que l'on ne peut connaître que si on le calcule), mais elle est infime devant les données (qui sont, elles, visibles même par les moins connaisseurs). Et tout ce que tu expliques, pldx, montre à quel point il ne faut se satisfaire des données telles qu'on les voit, mais au contraire, enquêter pour comprendre d'où viennent les soucis.
Pour éviter de modifier la solution X en une autre solution X* (notations de pldx) en perturbant le système initial,
j'ai préféré (dans ma réponse à Dlzlogic) modifier les deux premières équations en les divisant simplement par 7 (par exemple). Du coup, on a la même solution théorique X, mais on constate numériquement les effets de mauvais conditionnement sur les erreurs d'arrondi quand on utilise une méthode à base de pivot de Gauss (pivot simple, partiel, total), à moins de travailler avec 20 chiffres en mantisse par exemple.
La méthode exposée par pldx ne requiert qu'une mantisse très petite (4 ou 5 chiffres) pour être opérationnelle et donner une approximation (avec 2 chiffres significatifs) du résultat, sur le système 20 x 20 que j'ai présenté au-dessus. Par ailleurs, cette méthode est pour ainsi dire insensible aux perturbations numériques éventuelles (genre $\pm 0.1^{15}$ ou division par 7, ...).
C'est dommage que Dlzlogic ne dise pas ce qu'il obtient avec la méthode du pivot de Gauss lorsque les deux premières équations sont divisées par 7. Lui qui voulait un système numériquement instable... -
Cela étant, avec la méthode exposée par pldx
sur le système 20x20, pour obtenir un résultats avec 6 chiffres significatifs, il faut augmenter de manière impressionnante la mantisse des nombres flottants : apparemment, ce serait 60 (en base 10) !
Par exemple, avec une mantisse de 50 chiffres, on obtient seulement (et toujours) 2 chiffres significatifs au résultat...
Je précise que je termine le calcul par le choix du paramètre $x$ avec la méthode des moindres carrés appliquée à $C . X - B$
...et pour le système 10x10 qu'il étudie, il faut (sauf si je me trompe) une mantisse de 32 chiffres pour obtenir 6 chiffres au résultat.
Du tout, je me dis qu'avec autant de chiffres à la mantisse (32 ou 60), autant prendre l'inverse de la valeur propre $\lambda_0$ qui est proche de 0 : ainsi, on n'a plus l'incertitude sur le $x$, et on obtient le résultat tout de suite...... -
Sept ans plus tard...
Quand le conditionnement est à ce point mauvais, il faudrait au minimum 60 chiffres de précision sur $B$ pour récupérer 6 chiffres de précision sur $X$. Il en est ainsi quelle que soit la méthode de résolution utilisée.leon1789 a dit : cela étant, avec la méthode exposée par pldx, pour obtenir un résultats avec 6 chiffres significatifs,sur le système 20x20, il faut augmenter de manière impressionnante la mantisse des nombres flottants : apparemment, ce serait 60 (en base 10) !
Cordialement, Pierre.
![[Utilisateur supprimé]](https://w8.vanillicon.com/844673ca676c886b85d26c0b5ef6d73b_100.png)
Connectez-vous ou Inscrivez-vous pour répondre.
Bonjour!
Catégories
- 163.1K Toutes les catégories
- 7 Collège/Lycée
- 21.8K Algèbre
- 37.1K Analyse
- 6.2K Arithmétique
- 52 Catégories et structures
- 1K Combinatoire et Graphes
- 11 Sciences des données
- 5K Concours et Examens
- 11 CultureMath
- 47 Enseignement à distance
- 2.9K Fondements et Logique
- 10.3K Géométrie
- 62 Géométrie différentielle
- 1.1K Histoire des Mathématiques
- 68 Informatique théorique
- 3.8K LaTeX
- 39K Les-mathématiques
- 3.5K Livres, articles, revues, (...)
- 2.7K Logiciels pour les mathématiques
- 24 Mathématiques et finance
- 312 Mathématiques et Physique
- 4.9K Mathématiques et Société
- 3.3K Pédagogie, enseignement, orientation
- 10K Probabilités, théorie de la mesure
- 772 Shtam
- 4.2K Statistiques
- 3.7K Topologie
- 1.4K Vie du Forum et de ses membres
In this Discussion
Qui est en ligne 8
+5 Invités