Méthode de Gauss et pivotage

walidmens · January 2015

Bonjour,
J'ai programmé deux codes matlab pour résoudre un système linéaire par la méthode de Gauss. L'une est sans pivotage l'autre avec pivotage partiel. J'ai remarqué que les solutions calculées par la méthode à pivotage partiel sont plus stables numériquement que celle calculées sans pivotage, ceci est naturel. J'ai essayé avec avec une matrice de Pascal de taille 20 mais cette fois ce n'est pas le même phénomène. Bref, j'ai voulu savoir s'il y a des cas où la méthode de Gauss à pivots partiel ou total donne des solutions pires que celles données par la méthode de Gauss sans pivots. Merci d'avance pour votre aide.
Bien cordialement,
Walid

[Utilisateur supprimé] · January 2015

Bonjour,
Vous employez un terme que je n'avais jamais lu : "pivotage".
De plus vous distinguez le pivotage partiel du pivotage total, pouvez-vous donner des détails.
Ensuite, vous faites référence à la matrice de Pascal, je ne comprends pas très bien le rapport.
A ma connaissance, la méthode du pivot de Gauss est parfaitement adaptée aux calculs informatiques, et je ne pense pas qu'il y en ait de plus rapide et de plus précise (stable si vous voulez).

Rescassol · January 2015

Bonjour,

A chaque étape de la phase descendante, on élimine une inconnue entre la première ligne et chacune des suivantes.
Le coefficient permettant l'élimination s'appelle le pivot.
La règle du pivot maximal consiste à déterminer le plus grand de la colonne en valeur absolue, pui d'intervertir les deux lignes afin de l'amener en tête, ce qui améliore la précision globale.
On emploie quelquefois le terme de "pivotage" pour cette opération.
Par contre, je ne connais pas la différence entre pivotage partiel et total.
Enfin, comme je te l'ai déjà dit il y a quelque temps, il est faux de dire que la méthode du pivot de Gauss est la meilleure dans l'absolu, aussi bien du point de vue rapidité que précision, ça dépend du type de système considéré.

Cordialement,

Rescassol

[Utilisateur supprimé] · January 2015

Bonjour Rescassol,
J'ai pris la précaution de dire "je pense que ...".
S'il y a une autre méthode plus rapide et plus précise, peux-tu me dire où je pourrai la trouver. Ca vaut le coup de faire des essais comparatifs.
Quels sont les critères qui permettent de distinguer les types de systèmes ?

Rescassol · January 2015

Bonjour,

Je suis loin d'être un expert dans ce domaine, mais par exemple pour un système de taille $10000$ dont la matrice est creuse ("sparse" en anglo-saxon), c'est à dire ne possédant que peu d'éléments non nuls par rangée( $4$ ou $5$ par exemple), Gauss ne convient pas et il existe des algorithmes spécialisés (voir "sparse" sur Google).

Cordialement,

Rescassol

[Utilisateur supprimé] · January 2015

Oui, d'accord, il s'agit alors de gestion de grands tableaux qu'on appelle généralement raster.
Dans le présent sujet, sauf erreur, il s'agit de résolution de systèmes linéaires.

Par ailleurs, avec une matrice creuse, telle que définie dans l'article de Wiki, il n'est naturellement pas question d'échanger des lignes ou des colonnes, faire des multiplications etc.

Dans la question de Walid, as-tu compris le rapport avec la matrice de Pascal ?

Rescassol · January 2015

Bonsoir,

Oui, d'accord, il s'agit alors de gestion de grands tableaux qu'on appelle généralement raster.

En informatique pourquoi pas, en maths, on emploie exlusivement le terme de "matrice creuse" ou l'équivalent anglais "sparse matrix".

Dans le présent sujet, sauf erreur, il s'agit de résolution de systèmes linéaires.

Oui, c'est aussi ce dont je parle, les systèmes creux se résolvent tout comme les autres, mais pas par Gauss.
On les rencontre entre autre dans les problèmes d'éléments finis (résistance des matériaux par exemple).

Par ailleurs, avec une matrice creuse, telle que définie dans l'article de Wiki, il n'est naturellement pas question d'échanger des lignes ou des colonnes, faire des multiplications etc.

Oui, ça, c'est le "pivotage", comme je l'ai dit.

Dans la question de Walid, as-tu compris le rapport avec la matrice de Pascal ?

C'est seulement une matrice particulière sur laquelle il voulait tester son programme.
Voir "matrice de Pascal" dans Google.

Cordialement,

Rescassol

sylvainc2 · January 2015

Le pivotage partiel, c'est quand on permute seulement les lignes, e.g. on cherche le plus grand pivot, en valeur absolue, sur la colonne.

Le pivotage total est quand on permute les lignes mais aussi les colonnes.

Rescassol · January 2015

Bonjour,

Voilà un exemple de code vite fait (je n'ai fait que des interversions de lignes):

% Test de la méthode de Gauss avec la règle du pivot maximum
% sur une matrice de Pascal en Matlab

%-------------------------------------------------------------------------

% Initialisation de la session de travail

clc, clear all, close all 

format long;

p = 8; % Calculs avec p chiffres significatifs

%-------------------------------------------------------------------------

% Initialisation du système avec une matrice de Pascal

n = 15;
A = pascal(n);
B = rand(n,1);
X = zeros(n,1);

A0 = A; % Sauvegarde du système initial
B0 = B;

%-------------------------------------------------------------------------

% Triangulation du système

for j = 1:n - 1
    % Recherche du pivot maximum
    imax = 0;
    pmax = 0;
    for i=j:n
        if abs(A(i,j)) > abs(pmax)
            imax = i;
            pmax = A(i,j);
        end
    end
    if imax > j
        % Interversion de la ligne du pivot maximum avec la première
        T      = A(j,:);
        A(j,:) = A(i,:);
        A(i,:) = T;
        
        % Même chose sur le second membre
        T    = B(j);
        B(j) = B(i);
        B(i) = T;
    end
    
    % Elimination de tous les premiers termes avec celui de la première
    % équation et répercussion sur le second membre
    % chop(x,p) donne l'arrondi de x avec p chiffres significatifs
            
    for i = j + 1:n
        k = chop(A(i,j) / A(j,j),p);          
        A(i,:) = A(i,:) - chop(k * A(j,:),p);
        B(i)   = B(i)   - chop(k * B(j),p);  
    end
end

%-------------------------------------------------------------------------

% Résolution du système triangulaire

for i = n:-1:1 
    X(i) = chop((B(i) - chop(A(i,i+1:n) * X(i+1:n),p)) / A(i,i),p); 
end

%-------------------------------------------------------------------------

% Résultats

Solution = X'
Erreur = (A0 * X - B0)'
ErrMax = max(abs(Erreur))

format short

Cordialement,

Rescassol

pldx1 · January 2015

Bonjour,

Quand un problème est mal conditionné, il n'y a rien à y faire. A part ajouter des équations et travailler au sens des moindres carrés. Ici la matrice $\boxed{M}$ est définie par $M_{jk}=\binom {j+k}{k}$. Oh la la ! En tant qu'exercice de programmation pour matrices de petite taille, le pivot de Gauss peut se décrire par (SAGE):

def gauss_swap(maa,mab,swap=false):
    nn=maa.dimensions()[0]
    muu=block_matrix([maa,mab], nrows=true)
    for stp in range (nn-1):
        if swap:
            ask=muu.submatrix(stp,stp,ncols=1).apply_map(abs).list(); 
            qui= ask.index(max(ask))+stp
            muu.swap_rows(stp,qui)
        muu.set_row_to_multiple_of_row(stp,stp,1/muu[stp,stp])
        for jj in range (stp+1,nn):
            muu.add_multiple_of_row(jj,stp,- muu[jj,stp]);
    stp=nn-1; muu.set_row_to_multiple_of_row(stp,stp,1/muu[stp,stp])        
    for stp in range(nn-1,0,-1):
        for jj in range(stp):
            muu.add_multiple_of_row(jj,stp,- muu[jj,stp]);
    return Column(muu,nn)

L'idée qu'agiter les lignes ou les colonnes va systématiquement fournir un meilleur résultat est une idée répandue et fausse. Cela revient à croire qu'optimiser une usine consiste à optimiser chaque atelier, alors que cela consiste au contraire à sortir chaque atelier de son optimum local... et cela de façon réfléchie et coordonnée.

On utilise cette procédure par les commandes

nn=10
maa=matrix(nn,nn,lambda j,jj: binomial(j+jj,j))
mab=random_matrix(RR, nn,1)
mac=random_matrix(RR, nn,1)/10000000000
res0= (maa.I)*mab
res1=gauss_swap(maa,mab,swap=false); 
res2=gauss_swap(maa,mab,swap=true );
res4=gauss_swap(maa,mab+mac,swap=false);

Et l'on obtient:
\begin{eqnarray*}
mab & \approx & \left[+.31681, +.53692, +.94551, +.72250, -.12121, +.02042, +.58366, +.71641, -.90973, +.23793\right]\\
res & \approx & \left[-65.625, +547.23, -2013.8, +4325.1, -5988.9, +5552.9, -3449.7, +1384.6, -325.71, +34.202\right]\\
\delta_0,\delta_1,\delta_2 &= &4.72193995016 \times 10^{-09}, 1.40818001704 \times 10^{-09}, 1.45157421638 \times 10^{-08}\\
\varepsilon,\eta_4,\eta_5 & = &4.32673302134 \times 10^{-10}, 3.51203611615 \times 10^{-06}, 1.45998642935 \times 10^{-05}
\end{eqnarray*}

On a donc calculé $\boxed {M}^{-1}\cdot B$ où $B$ est une colonne de nombres uniformément répartis dans $[-1,+1]$, et obtenu le résultat $C$. On teste le résultat en calculant \[ \delta=\left|\phantom{x}\boxed {M}\cdot C-B\right|_1\] Si l'on prétend que $\delta$ est une mesure de qualité, on trouve que le pivot sans agitation des lignes est DIX fois meilleur que l'algorithme avec agitation... et est même TROIS fois meilleur que la multiplication par la matrice inverse (en dépit du fait que cette matrice est calculée exactement, dans l'anneau des nombres entier).

Cette propriété est à rapprocher de la construction par récurrence de la matrice de Pascal. En déconstruisant cette récurence de proche en proche, on évite d'empirer les choses: le gradient du pauvre en quelque sorte. Si l'on veut une mesure raisonnable de la qualité, on ajoute une colonne évanescente (norme $\varepsilon$) à la colonne $B$ et on voit ce que cela donne. Sur cet exemple, la valeur de $\eta /\varepsilon$ est de l'ordre de $10^4$.

Si l'on veut absolument réinventer la roue et reprogrammer le pivot de Gauss at the top level, au lieu d'utiliser des routines massivement parallèles et écrites at the processor level, alors implémenter le calcul de $\eta /\varepsilon$ en bordant $M$ non seulement par $B$ mais aussi par $B+\delta B$ est une optimisation majeure.

Cordialement, Pierre

Rouletabille · January 2015

Euh ... c'est quoi déjà le but du dit "pivotage" appliqué à une matrice donnée ?

[Utilisateur supprimé] · January 2015

Bonjour pldx1,
J'ai un peu de mal à suivre ton raisonnement.
Si un système est mal conditionné, c'est qu'il est mal posé en ce sens que les données son fausses, quelle que soit la raison, et donc, on ne peut rien faire.

Ici, on se trouve dans le contexte strictement informatique, et on doit résoudre un systèle linéaire de n équations à n inconnues. On se place dans le cas général où les paramètres des équations ainsi que les solutions sont des nombres réels. Puisque les différentes opérations consistent, en gros, à soustraire des nombres résultant de multiplication, il est important que le résultat de la soustraction ne soit pas très petit comparé aux opérandes. Ceci dans le seul but d'éviter de perdre des chiffres significatifs.
Il me semble que la méthode qui consiste à ordonner les équations par ordre décroissant des paramètres des inconnues résout cette difficulté.
J'avoue que c'est la première fois que je vois évoquer ce type de problème concernant la méthode du pivot de Gauss. Je crois que le seul vrai critère serait une méthode comparative ; on prend quelques systèmes, on les résout par les 2 méthodes, et il est facile ensuite de calculer les écarts (par la méthode qu'on veut).
On m'a opposé un certain système qui a la particulatité d'avoir des solutions souvent aberrantes. Malheureusement, je n'ai pas réussi à remettre la main dessus.

Personnellement, les systèmes que j'ai l'habitude de résoudre sont correctement conditionnés. Ma fonction les résout correctement, je n'ai donc pas eu l'occasion de me poser la question plus que cela, mais cette question mérite d'être posée et mérite de faire des essais comparatifs.
J'ai oublié de préciser que je me place exclusivement dans un contexte de "petits systèmes".

remarque · January 2015

> Si un système est mal conditionné, c'est qu'il est mal posé en ce sens que les données son fausses, quelle que soit la raison, et donc, on ne peut rien faire.

Non, ça veut dire que le conditionnement de la matrice est grand par rapport à 1. Dans ce cas, le calcul numérique de la solution du système linéaire correspondant pour certains seconds membres peut devenir impossible en pratique.

[Utilisateur supprimé] · January 2015

Bonsoir Remarque,
Pourrais-tu donner un exemple d'un tel système ?
J'avoue que j'ai un peu de mal à comprendre.

remarque · January 2015

La matrice de Hilbert est un exemple académique standard de matrice très mal conditionnée. Dans les vraies applications pratiques, par exemple quand on doit discrétiser des systèmes complexes d'équations aux dérivées partielles, il est illusoire de tenter une résolution de système linéaire sans préconditionnement, c'est-à-dire multiplication à gauche par une matrice qui diminue fortement le conditionnement.

[Utilisateur supprimé] · January 2015

La matrice de Hilbert est une matrice, pas un système linéaire.
On peut prendre la matrice de Hilbert comme coefficient pour les inconnues. Il faut quelqua-chose à droite du signe '='.

Rescassol · January 2015

Bonsoir,

Tout système, dont le membre de gauche a pour coefficients une matrice de Hilbert d'ordre suffisant, et pour membre de droite n'importe quoi est mal conditionné et nécessitera pas mal d'effort pour être résolu de façon satisfaisante.
Le conditionnement d'un système linéaire est une caractéristique essentiellement de la matrice de son membre de gauche.

Cordialement,

Rescassol

remarque · January 2015

> La matrice de Hilbert est une matrice, pas un système linéaire.

Certes, quelle perspicacité. Mais tu le fais exprès ou pas ? Sinon, ce n'est pas la peine que je me fatigue.

[Utilisateur supprimé] · January 2015

@ Remarque,
Je te demande un exemple, tu me donnes un cas particulier incomplet.
On parle de système linéaire et tu me réponds "matrice" (il est vrai que la définition de matrice est "tableau de n lignes et p colonnes").
Le sujet de ce fil concerne les systèmes linéaires. Si tu veux bien, restons dans le sujet.

PS. si cet "exemple" est celui que l'on m'avait soumi (et que je n'ai pas retrouvé), alors j'ai résolu le système sans difficulté particulière.

remarque · January 2015

> Je te demande un exemple, tu me donnes un cas particulier incomplet.

Ah excuses-moi, je ne savais pas que tu étais néophyte en matière de système linéaire pour ne pas savoir que les problèmes liés au conditionnement se trouvent uniquement dans la matrice et pas dans le second membre.

Alors comme ça, tu arrives à résoudre un système linéaire avec un second membre non nul et une modeste matrice de Hilbert $100\times 100$ (par exemple) ? J'aimerais bien le voir...

Rescassol · January 2015

Bonjour,

Je t'avais donné un exemple il y a quelque temps,
je ne sais pas si c'est celui dont tu parles:

Soit $M$ la matrice carrée d'ordre $n$ définie par $m_{i,j}=\dfrac{1}{i+j-1}$
C'est la matrice de Hilbert (les indices commencent à $1$).
Soit $A$ la matrice de $n$ lignes et $1$ colonne ne contenant que des $1$.
On calcule $B=M\times A$, puis on résout le système $M \times X = B$.
Normalement, on doit retomber sur $A$.
Programme le truc et observe $X$, puis l'erreur $X-A$.
Pour des valeurs de $n$ assez petites, ça va, mais quand $n$ augmente,
on constate assez rapidement de gros problèmes de précision.

Pour trouver une méthode de résolution satisfaisante de cet exemple,
il y a du boulot.

Cordialement,

Rescassol
PS: Je suppose que c'est ce que voulait dire Remarque, mais j'avais déjà écrit mon message.

[Utilisateur supprimé] · January 2015

Bon, je crois avoir compris ce que tu cherches : te donner l'occasion de dire que je suis un imbécile.
Concernant les systèmes d'équation, il s'agit de calcul numérique avec des données qui correspondebnt à une réalité, et j'ai précisé que je ne m'intéressais (connaissais) que les systèmes "petits".
J'en reviens à la question d'origine : précision du calcul d'un système linéaire par la méthode du pivot de Gauss.
Je ne te demande pas quel serait l'intérêt de résoudre un système dont les paramètres seraient les termes d'une matrice de Hilbert. C'est hors-sujet et ça ne m'intéresse pas.
Bonne soirée.

remarque · January 2015

Soyons concrets, on va juste prendre la matrice de Hilbert $20\times 20$. Je te donne un vecteur $b$ à $20$ coordonnées

    20.                       
    17.354641295237271947371  
    15.618373499565453954574  
    14.297125466739482035905  
    13.229500622319305236374  
    12.336875777899134476456  
    11.573481702709729646017  
    10.909802727235422992180  
    10.325488831126197197818  
    9.8058301074307667732910  
    9.3398112304786273796253  
    8.9189536438490719660876  
    8.5365857932898965287905  
    8.1873618821246605392616  
    7.8669327056365130346194  
    7.5717136131819771094342  
    7.2987167429496651038789  
    7.0454270799245604095518  
    6.8097091682173376270271  
    6.5897357459388130607181

Peux-tu me calculer le vecteur $x$ tel que $Hx=b$ ? Donc il y a une matrice $H$, la matrice de Hilbert $20\times 20$, un vecteur donné $b$ et un vecteur inconnu $x$. On est bien d'accord qu'il s'agit maintenant d'un vrai système linéaire, avec un second membre inoffensif ?

remarque · January 2015

> Bon, je crois avoir compris ce que tu cherches : te donner l'occasion de dire que je suis un imbécile.

En fait, je pensais que tu te débrouillais très bien tout seul. :-D

Mais bon, on ne saura donc pas qui est $x$. Pour un petit système par la méthode de Gauss. Quel dommage.

Rescassol · January 2015

Bonne nuit,

Remarque, je crois que j'ai la solution de ton exemple, mais pas par Gauss :-D

Cordialement,

Rescassol

remarque · January 2015

Rescassol écrivait:

> Bonne nuit,
>
> Remarque, je crois que j'ai la solution de ton exemple, mais pas par Gauss :-D

Ahah ! Par divination alors ? Social engineering ? C'est la NSA qui t'a rencardé ? Je suis curieux du coup ! :-D

[Utilisateur supprimé] · January 2015

@ Remarque,
Il se trouve que par chance, je suit tombé sur mon test, et qui était justement avec 20 équations.
C'est la raison pour laquelle ma réponse est si rapide.
DEPART
1.000 0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 = 3.598
0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 = 2.645
0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 = 2.191
0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 = 1.901
0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 = 1.693
0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 = 1.533
0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 = 1.404
0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 = 1.299
0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 = 1.209
0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 = 1.133
0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 = 1.066
0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 = 1.007
0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 = 0.955
0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 = 0.909
0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 = 0.867
0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 = 0.829
0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 = 0.794
0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 = 0.762
0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 = 0.733
0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 0.026 = 0.706
1.000000 1.000026 0.999146 1.011423 0.929744 1.152536 1.480044 -2.863703 11.029418 -10.157341 0.905777 12.293068 -5.077246 0.428064 -6.076828 5.733979 19.865414 -29.585184 18.621477 -2.689815
Ligne 1 S=3.597740
Ligne 2 S=2.645359
Ligne 3 S=2.190813
Ligne 4 S=1.900958
Ligne 5 S=1.692625
Ligne 6 S=1.532625
Ligne 7 S=1.404420
Ligne 8 S=1.298600
Ligne 9 S=1.209314
Ligne 10 S=1.132686
Ligne 11 S=1.066019
Ligne 12 S=1.007368
Ligne 13 S=0.955285
Ligne 14 S=0.908664
Ligne 15 S=0.866648
Ligne 16 S=0.828552
Ligne 17 S=0.793830
Ligne 18 S=0.762034
Ligne 19 S=0.732794
Ligne 20 S=0.705803

Tu m'excusera, mais j'ai pas pris ton vecteur B, j'ai gardé celui qu'on m'avait donné précédemment.
Je veux bien faire des simulations, mais j'aime pas trop montrer qu'un membre dit n'importe quoi.
Au passage, tu me diras dans quel contexte tu arrives à des valeurs pour B avec 22 ou 23 chiffres significatifs.
Il me semble que sur ce sujet précis, impossibilité ou pas, on a fait le tour de la question.
Alors revenons à la question d'origine, s'il te plait.

Rescassol · January 2015

Bonne nuit,

> Par divination alors ?

Oui, je me suis demandé par quel procédé forcément simple tu avais fabriqué ton exemple,
et je suis vite tombé sur le bon par essais et erreurs.

Cordialement,

Rescassol

remarque · January 2015

> Tu m'excusera, mais j'ai pas pris ton vecteur B,

Ah non, c'est inexcusable. Se dérober devant l'obstacle, tsk, tsk... Je suis un peu déçu, mais pas tant que ça. :-D

remarque · January 2015

@ Rescassol : je m'en doutais un peu, mais je suis content que ce ne soit pas l'option NSA...

remarque · January 2015

Donc je résume pour dlzlogic : il y a essentiellement deux méthodes 1. la méthode de Gauss dont j'attends toujours de voir le résultat et 2. la divination qui donne plutôt des bons résultats sur cet exemple particulier.

[Utilisateur supprimé] · January 2015

Tu avais expliqué que les valeurs du second membre n'avaient pas grande importance, mais comme tu insistes :
"Ah non, c'est inexcusable. Se dérober devant l'obstacle, tsk, tsk... Je suis un peu déçu, mais pas tant que ça."
Sans commentaire :-D.

1.000 0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 = 17.355
0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 = 15.618
0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 = 14.297
0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 = 13.230
0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 = 12.337
0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 = 11.573
0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 = 10.910
0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 = 10.325
0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 = 9.806
0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 = 9.340
0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 = 8.919
0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 = 8.537
0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 = 8.187
0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 = 7.867
0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 = 7.572
0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 = 7.299
0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 = 7.045
0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 = 6.810
0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 = 6.590
0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 0.026 = 6.200
114240678.639305 -17390761859.524151 649015070892.487549 -10296415651124.050781 85224105961665.781250 -401003842368058.937500 1082763315074396.625000 -1514915433455079.250000 571551479894145.375000 650106297378619.625000 948856011170290.875000 -3468427722440542.500000 1760599110935123.750000 -60503842414165.945310 3787951888378411.000000 -3361486357068956.000000 -6515507344182813.000000 12260533635134778.000000 -7417184649625362.000000 1601108094320109.000000
Ligne 1 S=17.360695
Ligne 2 S=15.540413
Ligne 3 S=14.297478
Ligne 4 S=13.206627
Ligne 5 S=12.394733
Ligne 6 S=11.595707
Ligne 7 S=10.942120
Ligne 8 S=10.264938
Ligne 9 S=9.827183
Ligne 10 S=9.379623
Ligne 11 S=8.900822
Ligne 12 S=8.552750
Ligne 13 S=8.195312
Ligne 14 S=7.864983
Ligne 15 S=7.523705
Ligne 16 S=7.276920
Ligne 17 S=7.020969
Ligne 18 S=6.787663
Ligne 19 S=6.551937
Ligne 20 S=6.192085

Nota : pour un membre qui essayerait de simuler le problème, la liste de Remarque ne contient que 19 valeurs pour le membre de droite.
Je suppose que c'est normal, j'y connais rien. J'ai ragouté une valeur au hasard : 6.20.

@ Remarque, tu ne m'as pas dit l'utilisation de ce type de calcul et aussi le type de mesure qui te donnes 22 ou 23 chiffres significatifs.
Ceci dit, on est encore hors-sujet.

Rescassol · January 2015

Bonjour,

>la liste de Remarque ne contient que 19 valeurs pour le membre de droite

Tu as mal compté, il y en a 20, la première étant $20.$

Cordialement,

Rescassol

remarque · January 2015

> tu ne m'as pas dit l'utilisation de ce type de calcul

Ah tiens ?

comme chacun peut le voir a écrit:

La matrice de Hilbert est un exemple académique standard de matrice très mal conditionnée.

« Académique » signifie dans ce contexte qu'il s'agit d'un exemple qui n'a pas nécessairement de justification hors de lui-même (encore qu'ici il s'agisse d'une matrice de Gram particulière associée à la base canonique des polynômes pour le produit scalaire de $L^2(0,1)$), si ce n'est qu'il est bien connu pour posséder un conditionnement explosant exponentiellement avec la dimension. La vertu principale de l'exemple académique est de montrer qu'il faut se méfier de situations en apparence sans danger et où le naïf moyen pourrait croire que la méthode de Gauss avec pivot partiel ou total permet de calculer la solution quand on travaille en flottants.

Dans les applications réelles, il s'agit d'une question cruciale

comme chacun peut le voir a écrit:

Dans les vraies applications pratiques, par exemple quand on doit discrétiser des systèmes complexes d'équations aux dérivées partielles, il est illusoire de tenter une résolution de système linéaire sans préconditionnement, c'est-à-dire multiplication à gauche par une matrice qui diminue fortement le conditionnement.

Mais je ne suis pas convaincu d'après ce que tu écris que tu aies la moindre idée de ce qu'est le conditionnement et de à quoi ça sert en algèbre linéaire numérique, ni de ce qu'est l'algèbre linéaire numérique. En effet, tu n'écrirais normalement pas

> le type de mesure qui te donnes 22 ou 23 chiffres significatifs.

si tu savais à quoi sert le conditionnement. Puisque tu n'aimes pas le hors-sujet, je t'invite à lire le premier message de ce fil – celui qui devrait normalement fixer ce fameux sujet avant toutes les dérives ultérieures – il n'y est nullement question de mesure ou de chiffres significatifs.

pldx1 · January 2015

Bonjour,
Une réponse possible au défi de remarque aurait été:

+1.000000175973820467973546800, 
+1.999965416714291101146861200, 
+3.001677278747650067756322100, 
+3.964888734061474598368546800, 
+5.394072757775592855498738000, 
+3.364183996040071408461452000, 
+18.03726169443823809959627520, 
-21.06829937852191584852851200, 
+54.23952704740837910576640000, 
-20.65398693225517040761879040, 
-4.674344876271487346063539200, 
+46.12791189848842721340067840, 
+20.01686762203943105314911040, 
-18.69344059310310049458626560, 
+7.679265519855917266900696000, 
+49.01941233471778052988221440, 
+22.94420734547925555284666000, 
-21.67275234942480192952709120, 
+45.89180407659612637133936640, 
+14.08177808261579859120568600

En effet le produit de cette colonne par la matrice de Hilbert de taille 20 redonne la colonne proposée avec une précision égale à celle correspondant à la solution divinatoire. Mais dlzlogic est tellement incompétent en informatique qu'il n'arrive même pas à compter les lignes d'une matrice de 20 lignes...

Pour en revenir à la notion de conditionnement, une réponse correcte au problème posé par remarque ne consiste pas à exhiber une solution et à passer à autre chose. En effet, l'ensemble des solutions est un large espace dépendant d'environ dix-huit paramètres. Pour s'en faire une idée, il suffit de comparer le nombre de chiffres significatifs fournis (24) avec la liste des valeurs propres de la matrice. Evidemment, le "24 chiffres" est adapté à la norme du sup, tandis que le module des valeurs propres est adapté à la norme euclidienne (d'où le "environ" 18: je n'ai qu'une envie modérée de mesurer les normes subordonnées successives des opérateurs obtenus par restriction du multiplicateur de Hilbert aux sous-espaces emboités $E_k$ engendrés par les vecteurs propres $V_j$ (avec $j\ge k$).

Les vecteurs propres $V_j$, tels qu'ils tombent du logiciel, sont de norme euclidienne $1$. Ajouter une telle colonne fois $E-22$ fois $1/\mu_j$ laisse invariant le résultat. Or la suite des valeurs propres décroit à toute vitesse (c'est précisément pour cela que le problème est mal conditionné).

Cordialement, Pierre

pldx1 · January 2015

N'hésitons pas à en remettre une couche, et même à utiliser l'exemple où dlzlogic suit tombé sur son test (sic). On se donne une cible et l'on cherche un vecteur de départ, la matrice multiplicatrice étant la matrice de Hilbert. Le problème est qu'il y a autant de "vecteurs de départ" que l'on veut. Le vecteur $V_1$ convient et le vecteur $V_2$ convient tout autant. En effet, les colonnes $M\cdot V_1$ et $M\cdot V_2$ sont toutes deux égales à la cible.

         V1          V2       ma*V1       ma*V2       cible
    1.00000     1.00000     3.59774     3.59774     3.59774
    1.00003     0.99968     2.64536     2.64536     2.64536
    0.99915     1.01592     2.19081     2.19081     2.19081
    1.01142     0.66031     1.90096     1.90096     1.90096
    0.92974     4.87047     1.69262     1.69262     1.69262
    1.15254   -25.20562     1.53262     1.53262     1.53262
    1.48004   111.85266     1.40442     1.40442     1.40442
   -2.86370  -293.54670     1.29860     1.29860     1.29860
   11.02942   463.42469     1.20931     1.20931     1.20931
  -10.15734  -316.69721     1.13269     1.13269     1.13269
    0.90578  -155.83767     1.06602     1.06602     1.06602
   12.29307   353.57219     1.00737     1.00737     1.00737
   -5.07725    65.09143     0.95528     0.95528     0.95528
    0.42806  -326.50634     0.90866     0.90866     0.90866
   -6.07683   -79.28417     0.86665     0.86665     0.86665
    5.73398   335.92810     0.82855     0.82855     0.82855
   19.86541    79.30749     0.79383     0.79383     0.79383
  -29.58518  -426.31271     0.76203     0.76203     0.76203
   18.62148   287.53952     0.73279     0.73279     0.73279
   -2.68982   -61.87203     0.70580     0.70580     0.70580

Les praticiens qui pratiquent ont une étiquette à coller sur ce genre de phénomène: "encore un problème mal conditionné". Et pendant ce temps là, les bavards de la pratique qui ne pratiquent pas ont du mal à suivre.

@Rouletabille. Suite à ta question: dans une résolution numérique on voit se combiner des "incertitudes de troncature" et des "incertitudes de méthode". On se demande donc ce que l'on peut faire pour que les choses ne deviennent pas pire que ce qui est imposé par le conditionnement du problème et la longueur du "mot numérique".
Dans l'algorithme de Gauss, on est obligé de permuter les lignes lorsque l'élément qui arrive en position pivot est nul. Il est intéressant par rapport aux troncatures de procéder à ce genre de permutation dès que la norme de l'élément en position pivot est significativement plus petite que celle d'autres éléments placés en dessous-de lui. Le "pivotage" consiste à faire cela systématiquement, et même à se demander si une permutation des colonnes ne serait pas utile.
Mon opinion est qu'il ne faut pas trop y croire, surtout pour un problème mal conditionné. En tel cas, la principale source d'incertitude est l'explosion des sous espaces propres liés aux valeurs propres de petit module et, face à cela, il n'y a rien à faire (à part projeter sur l'un des sous-espaces $E_m$ engendré par les vecteurs propres $V_j, j\leq m$ lorsque, comme ici, la matrice est celle d'un produit scalaire... et que l'on a classé les valeurs propres).
Dans tous les cas, il est utile est de reprendre plusieurs fois les calculs avec des données légèrement perturbées. Examiner ce qui reste stable et ce qui est chamboulé permet de se faire une opinion sur la validité des résultats..

Cordialement, Pierre.

.

JLT · January 2015

@dlzlogic : merci de ne pas provoquer de polémique sans fin. Ceci était l'une des raisons de ton bannissement il n'y a pas si longtemps.

remarque · January 2015

Il semble qu'il y ait eu des échanges auxquels j'ai malheureusement échappé. J'appuie complètement ce qu'explique pldx1. Les stratégies de pivot sont là pour diminuer l'accumulation des erreurs d'arrondi, mais le problème principal ici n'est pas là. Il est lié au conditionnement de la matrice, quelle que soit la méthode utilisée pour résoudre numériquement le système (en flottants donc, pas en arithmétique exacte). Pour la matrice de Hilbert $20\times 20$ le conditionnement en norme $2$ semble être de l'ordre de $1{,}786\;10^{18}$.* Si on veut récupérer la solution divinatoire à $10^{-25}$ près (parce que nous le valons bien), il faut donc grossièrement avoir un second membre connu à $6\;10^{-44}$ près, ce qui n'est pas vraiment le cas. Il y a donc foule de candidats solution qui ne sont vraiment pas près de la vraie solution.

* En fait, ce n'est pas si clair... je vois une autre façon de faire qui le met aux alentours de $10^{30}$, ce qui n'arrange pas les choses... enfin bref, c'est un grand nombre. :-D

h · January 2015

(Cet échange me conforte dans ma décision de boycotter le forum tant que dlzlogic n'est pas à nouveau banni.)

remarque · January 2015

Tu boycottes, mais tu lurkes ! :-D:-D

Sylviel · January 2015

Remarque a écrit:

Mais je ne suis pas convaincu d'après ce que tu écris que tu aies la moindre idée de ce qu'est le conditionnement et de à quoi ça sert en algèbre linéaire numérique, ni de ce qu'est l'algèbre linéaire numérique.

Non il n'en a pas la moindre idée, pour lui les matrices sont une invention des matheux pour fair des exercices sans application pratique et la dernière fois que j'ai essayé d'expliquer le conditionnement il ne voyait pas pourquoi si on prend deux b différents mais proche on devrait avoir deux x proches...

En revanche je tiens à souligner que personne n'a expliqué en quoi la solution qu'il a proposé était mauvaise.

h · January 2015

(Je lurke (je découvre le mot) en attendant le bannissement :-).)

pldx1 · January 2015

Mon cher remarque,

La méthode divinatoire ne nécessite pas autant de décimales que cela. Sinon, ce ne serait plus qu'une méthode banale et non pas une méthode divine à Thoir. La matrice de Hilbert est la matrice d'un produit scalaire. Elle est donc diagonalisable, la matrice de passage étant une matrice orthogonale. Le point clef est que l'on peut calculer les colonnes de cette matrice autrement qu'en utilisant un pivot de Gauss pour déterminer le noyau de $M-\mu_j$. Si l'on utilise un logiciel raisonnable et que l'on a déclaré par typage que la matrice est symétrique, cet algorithme à la sauce gradient conjugué se sélectionne de lui même et l'on obtient un résultat exact (aux arrondis près). Utilisant une précision de 30 chiffres, on trouve $\Delta=P^{-1}\cdot M \cdot P$=

[7.7813078695625931472792437169E-29, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 3.3763126505758943544371654694E-26, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 7.0264420384506850155910815351E-24, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 9.3311941011696316300842016924E-22, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 8.8800759473655671826756209033E-20, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 6.4467646276569711232097951528E-18, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 3.7109770253440723633032491451E-16, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 1.7379067082989644306349013733E-14, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 6.7408082331639843630948575676E-13, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 2.1928907569892015252599786599E-11, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 6.0360953293918637191088279035E-10, 0., 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.41395475825338048152292481065E-8, 0., 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.82765205524785408866049199960E-7, 0., 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.48305100488023714517689871301E-5, 0., 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.70334314731935324532726552531E-4, 0., 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.86767110917149801722206995732E-3, 0., 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.89611286148564804426072749481E-2, 0., 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.75595821305440958439498590984E-1, 0., 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .487038406572048867803407847772, 0.], 
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.90713472040725310302314635823]

Dans cette expression, les quantités marquées $0.$ sont toutes inférieures à $1E-29$. On projette les cibles en calculant $P^{-1}\cdot C $ et on obtient:

      mu(j)           P^(-1)*C1       P^(-1)*C2
  7.7813079E-29   6.8488138E-07   3.9247518E-16
  3.3763127E-26  -4.8624463E-08   9.0676274E-16
  7.0264420E-24   5.0862900E-07   1.0629963E-15
  9.3311941E-22  -1.4760611E-08  -7.0585935E-16
  8.8800759E-20  -2.1461184E-07  -2.7805723E-16
  6.4467646E-18   1.6188525E-08   2.3082637E-16
  3.7109770E-16  -6.9170762E-08  -1.2677343E-15
  1.7379067E-14   1.3120978E-07   6.3739237E-16
  6.7408082E-13   3.6898444E-07  -1.0317922E-15
  2.1928908E-11  -2.4619612E-07   7.6183577E-14
  6.0360953E-10   2.3519330E-07  -9.6934046E-12
  1.4139548E-08  -1.5579140E-08   9.3032988E-10
  2.8276521E-07  -1.5737778E-08  -6.9170259E-08
  4.8305100E-06   7.0705339E-08   3.9667793E-06
  7.0334315E-05  -2.5810868E-06  -1.7354219E-04
  8.6767111E-04   1.0047520E-04   5.6689643E-03
  8.9611286E-03  -3.2555820E-03  -1.3277057E-01
  7.5595821E-02   7.7451020E-02   2.0515892E+00
  4.8703841E-01   1.2061244E+00   1.7187916E+01
  1.9071347E+00   6.7886169E+00   4.7449371E+01

Dans la colonne $P^{-1}\cdot C_1$, on voit que les derniers coefficients sont de l'ordre de 1E-7 au lieu de continuer à décroitre comme les $\mu_j$. Cela signale la zone où le bruit l'emporte sur le signal, et l'on décide de couper tous les coefficients au dela du 6-ème. Dans la colonne $P^{-1}\cdot C_2$, l'amplitude du bruit est de l'ordre de 1E-15 (il y a moins de bruit parce que la colonne est donnée avec plus de précision). On coupe donc au delà du dixième coefficient.

Il reste à reconstituer les colonnes $V_j$ en calculant $P\cdot\Delta_{cut}^{-1}\cdot \left(P^{-1}\cdot V_j\right)$ le $cut$ étant comme décrit ci-dessus. On obtient:

    0.999943904     1.000000001
    1.000224941     1.999999907
    1.001417544     3.000002567
    0.996775820     3.999973931
    0.998452460     5.000119939
    1.001095802     5.999748145
    1.002377383     7.000155194
    1.002246171     8.000176114
    1.001315125     8.999896517
    1.000170301     9.999783906
    0.999191732    10.999937710
    0.998562856    12.000153260
    0.998327484    13.000217078
    0.998445143    14.000081231
    0.998832271    14.999865477
    0.999389295    15.999755161
    1.000016872    16.999870970
    1.000624697    18.000158821
    1.001135577    19.000320102
    1.001486619    19.999783965

C'est le moment de remarquer un phénomène un peu trop artificiel pour être un simple artefact. Les personnes formées au jacassin nouveau iront jusqu'à s'exclamer :"ne serait-ce pas le moment d'émettre une conjecture ?".

Cordialement, Pierre

Rescassol · January 2015

Bonjour,

Curieusement, Matlab donne pour le conditionnement de la matrice de Hilbert $20\times20$ la valeur $2.038253293125565\times 10^{18}$ alors qu'Octave donne $1.778704761054288384\times 10^{18}$.
C'est le même ordre de grandeur, mais je me demande, pourquoi cette différence ?

Cordialement,

Rescassol

remarque · January 2015

@ pldx1 : oui je me suis placé dans le worst case scenario sans aller chercher très loin... Ceci dit, la première ligne de mon second membre mystère vend plus ou moins la mèche.

@ Rescassol : j'imagine que cela dépend comment chaque logiciel s'y prend pour approcher les valeurs singulières. A ces ordres de grandeur-là, on peut avoir des doutes sur les dernières premières décimales...

remarque · January 2015

@ sylviel : en quoi sa solution est mauvaise, difficile à dire dans la mesure où ce n'était pas le bon second membre. Il faudrait déjà récupérer ses données pour regarder le résidu, mais comme on ne peut pas copier une colonne d'un coup de la façon dont il a collé son résultat, c'est pas très emballant. Sachant que de toutes façons, il y a un paquet d'autres vecteurs qui vont donner un petit résidu (je suppose que son résidu est petit, ce qui n'est même pas garanti).

remarque · January 2015

Bon, j'ai fait l'effort de récupérer les données de dlzlogic pour voir ce que ça donne. C'est pas top, loin de là. Ci-joint un tout petit programme pour la matrice de Hilbert H, la matrice de Hilbert sauce dlzlogic (i.e. tronquée à 3 décimale, car je suppose qu'on ne peut pas mesurer la matrice de Hilbert à plus de trois décimales, raisonnablement), son second membre bdlz, sa solution xdlz, le résidu avec sa matrice residlz, le résidu avec la vraie matrice à la précision de scilab, residuvrai. On est loin du compte, mais il est sans doute satisfait de son résultat (y a-t-il un bannissement temporaire ? Ca fait drôle de parler de qqun à le troisième personne...). A titre de comparaison, j'ai mis ce que calcule scilab avec ses données, newxdlz, qui donne un bien meilleur résidu newresidlz, sans que l'on sache par ailleurs, s'il existait un vrai x quelque part, quelle serait la relation de ce newxdlz avec le vrai x...

N=20

H=zeros(N,N)
for i=1:N
    for j=1:N
        H(i,j)=1/(i+j-1)
    end
end


Hdlz=[1.000 0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 ;
0.500 0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 ;
0.333 0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 ;
0.250 0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 ;
0.200 0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 ;
0.167 0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 ; 
0.143 0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 ;
0.125 0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 ; 
0.111 0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 ; 
0.100 0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 ;
0.091 0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 ;
0.083 0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 ; 
0.077 0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 ;
0.071 0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 ; 
0.067 0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 ; 
0.062 0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 ; 
0.059 0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 ; 
0.056 0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 ;
0.053 0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 ; 
0.050 0.048 0.045 0.043 0.042 0.040 0.038 0.037 0.036 0.034 0.033 0.032 0.031 0.030 0.029 0.029 0.028 0.027 0.026 0.026 ] 

bdlz=[17.355 
15.618 
14.297 
13.230 
12.337 
11.573 
10.910 
10.325 
9.806 
9.340 
8.919 
8.537 
8.187 
7.867 
7.572 
7.299 
7.045 
6.810 
6.590 
6.200 ]

xdlz=[17.360695 
15.540413 
14.297478 
13.206627 
12.394733 
11.595707 
10.942120 
10.264938 
9.827183 
9.379623 
8.900822 
8.552750 
8.195312 
7.864983 
7.523705 
7.276920 
7.020969 
6.787663 
6.551937 
6.192085 ]

residlz=Hdlz*xdlz-bdlz
residuvrai=H*xdlz-bdlz
disp(residlz)

disp(residuvrai)

newxdlz=Hdlz\bdlz

disp(newxdlz)

newresidlz=Hdlz*newxdlz-bdlz

disp(newresidlz)

et le résultat des courses :

-->exec('~/hilbertmatrix.sce', -1)

disp(residlz) 
    30.246434  
    17.30424   
    12.055339  
    9.1192497  
    7.2289535  
    5.9094791  
    4.9291525  
    4.1855491  
    3.6045364  
    3.1362413  
    2.7524267  
    2.4316666  
    2.170802   
    1.9409084  
    1.7527225  
    1.5820263  
    1.4481608  
    1.3176379  
    1.1959667  
    1.2757782  

disp(residuvrai) 
    30.246107  
    17.301617  
    12.055754  
    9.1165091  
    7.2238747  
    5.9037109  
    4.9319745  
    4.191551   
    3.6091188  
    3.1421913  
    2.7608423  
    2.4441155  
    2.1792381  
    1.9533247  
    1.759923   
    1.5930563  
    1.4485757  
    1.3207076  
    1.2087318  
    1.2937473  
 
disp(newxdlz)
  - 7.3053881  
    4.9133547  
    109.74095  
  - 31.281908  
  - 88.497335  
  - 96.62244   
  - 19.9283    
    27.562282  
  - 57.814007  
    69.722061  
    55.929314  
    56.187249  
    130.25188  
    35.870815  
    215.2377   
    1.5398888  
    3.3067182  
  - 93.417967  
  - 55.528873  
  - 54.018303  
 
disp(newresidlz) 
10^(-14) *
 
  - 0.3552714  
    0.5329071  
  - 0.3552714  
    0.5329071  
  - 0.3552714  
    0.1776357  
  - 0.3552714  
    0.3552714  
    0.1776357  
    0.         
    0.         
    0.5329071  
  - 0.1776357  
    0.4440892  
  - 0.0888178  
  - 0.0888178  
  - 0.1776357  
    0.1776357  
  - 0.1776357  
    0.0888178

[Dlz n'a pas encore été banni mais l'épée de Damocles n'est pas loin. --JLT]

Rescassol · January 2015

Bonne nuit,

Tiens, tu programmes en quoi ? En Sage ?
En Matlab ou Octave, on se contente de $H=hilb(20)$.
Et le conditionnement par $c=cond(hilb(20))$.
Ceci dit, tes calculs sont limpides.

Cordialement,

Rescassol

remarque · January 2015

J'utilise ici scilab. Apparemment, il n'y a pas de matrice Hilbert préprogrammée. hilb en scilab a plutôt à voir avec la transformée de Hilbert... Par contre, il y a testmatrix('hilb',n) qui renvoie l'inverse de la matrice de Hilbert. Bon, je ne sais pas avec quelle fiabilité pour n=20 :-D... En plus, c'est une matrice à coefficients entiers. Le conditionnement est aussi retourné par un fonction cond, mais visiblement pas la même...

-->testmatrix('hilb',20)
 ans  =
 
 10^14 *
 
 
         column  1 to 10
 
    4.000D-12  - 7.980D-10    5.267D-08  - 0.0000017    0.0000329  - 0.0004119    0.0035695  - 0.0223730    0.1044075  - 0.3700665  
  - 7.980D-10    0.0000002  - 0.0000158    0.0005478  - 0.0109556    0.1408574  - 1.2461968    7.9349672  - 37.49272     134.2332   
    5.267D-08  - 0.0000158    0.0012483  - 0.0451918    0.9296589  - 12.201773    109.66532  - 707.00558    3374.3448  - 12181.662  
  - 0.0000017    0.0005478  - 0.0451918    1.6828547  - 35.339949    471.19932  - 4287.9139    27923.15   - 134380.16    488515.9   
    0.0000329  - 0.0109556    0.9296589  - 35.339949    753.91892  - 10177.905    93554.484  - 614309.29    2977037.3  - 10886926.  
  - 0.0004119    0.1408574  - 12.201773    471.19932  - 10177.905    138789.62  - 1286374.2    8505820.9  - 41465877.    1.524D+08  
    0.0035695  - 1.2461968    109.66532  - 4287.9139    93554.484  - 1286374.2    12006159.  - 79860208.    3.913D+08  - 1.445D+09  
  - 0.0223730    7.9349672  - 707.00558    27923.15   - 614309.29    8505820.9  - 79860208.    5.339D+08  - 2.628D+09    9.741D+09  
    0.1044075  - 37.49272     3374.3448  - 134380.16    2977037.3  - 41465877.    3.913D+08  - 2.628D+09    1.298D+10  - 4.830D+10  
  - 0.3700665    134.2332   - 12181.662    488515.9   - 10886926.    1.524D+08  - 1.445D+09    9.741D+09  - 4.830D+10    1.802D+11  
    1.0092721  - 369.14129    33733.835  - 1360865.7    30483392.  - 4.287D+08    4.079D+09  - 2.760D+10    1.373D+11  - 5.136D+11  
  - 2.1332343    785.6866   - 72227.047    2928673.   - 65895143.    9.303D+08  - 8.884D+09    6.029D+10  - 3.007D+11    1.128D+12  
    3.5006922  - 1297.0065    119843.4   - 4881121.7    1.103D+08  - 1.562D+09    1.496D+10  - 1.018D+11    5.091D+11  - 1.914D+12  
  - 4.4431862    1654.6426  - 153571.51    6279368.5  - 1.423D+08    2.023D+09  - 1.943D+10    1.325D+11  - 6.642D+11    2.502D+12  
    4.3162381  - 1614.5428    150437.4   - 6172576.5    1.403D+08  - 2.000D+09    1.926D+10  - 1.317D+11    6.613D+11  - 2.496D+12  
  - 3.1472569    1181.8875  - 110506.49    4548214.3  - 1.037D+08    1.481D+09  - 1.430D+10    9.797D+10  - 4.929D+11    1.864D+12  
    1.6661948  - 627.87776    58888.324  - 2430452.     55553189.  - 7.954D+08    7.693D+09  - 5.281D+10    2.662D+11  - 1.008D+12  
  - 0.6044040    228.46472  - 21487.107    889043.24  - 20367173.    2.922D+08  - 2.832D+09    1.947D+10  - 9.830D+10    3.728D+11  
    0.134312   - 50.910965    4800.1767  - 199061.87    4569768.2  - 65690418.    6.376D+08  - 4.392D+09    2.220D+10  - 8.432D+10  
  - 0.0137847    5.2381681  - 495.00688    20570.286  - 473116.58    6812878.8  - 66236321.    4.569D+08  - 2.313D+09    8.795D+09  
 
         column 11 to 20
 
    1.0092721  - 2.1332343    3.5006922  - 4.4431862    4.3162381  - 3.1472569    1.6661948  - 0.6044040    0.134312   - 0.0137847  
  - 369.14129    785.6866   - 1297.0065    1654.6426  - 1614.5428    1181.8875  - 627.87776    228.46472  - 50.910965    5.2381681  
    33733.835  - 72227.047    119843.4   - 153571.51    150437.4   - 110506.49    58888.324  - 21487.107    4800.1767  - 495.00688  
  - 1360865.7    2928673.   - 4881121.7    6279368.5  - 6172576.5    4548214.3  - 2430452.     889043.24  - 199061.87    20570.286  
    30483392.  - 65895143.    1.103D+08  - 1.423D+08    1.403D+08  - 1.037D+08    55553189.  - 20367173.    4569768.2  - 473116.58  
  - 4.287D+08    9.303D+08  - 1.562D+09    2.023D+09  - 2.000D+09    1.481D+09  - 7.954D+08    2.922D+08  - 65690418.    6812878.8  
    4.079D+09  - 8.884D+09    1.496D+10  - 1.943D+10    1.926D+10  - 1.430D+10    7.693D+09  - 2.832D+09    6.376D+08  - 66236321.  
  - 2.760D+10    6.029D+10  - 1.018D+11    1.325D+11  - 1.317D+11    9.797D+10  - 5.281D+10    1.947D+10  - 4.392D+09    4.569D+08  
    1.373D+11  - 3.007D+11    5.091D+11  - 6.642D+11    6.613D+11  - 4.929D+11    2.662D+11  - 9.830D+10    2.220D+10  - 2.313D+09  
  - 5.136D+11    1.128D+12  - 1.914D+12    2.502D+12  - 2.496D+12    1.864D+12  - 1.008D+12    3.728D+11  - 8.432D+10    8.795D+09  
    1.467D+12  - 3.230D+12    5.492D+12  - 7.194D+12    7.188D+12  - 5.376D+12    2.912D+12  - 1.078D+12    2.442D+11  - 2.551D+10  
  - 3.230D+12    7.123D+12  - 1.214D+13    1.592D+13  - 1.594D+13    1.194D+13  - 6.474D+12    2.401D+12  - 5.444D+11    5.691D+10  
    5.492D+12  - 1.214D+13    2.071D+13  - 2.722D+13    2.728D+13  - 2.046D+13    1.111D+13  - 4.126D+12    9.366D+11  - 9.802D+10  
  - 7.194D+12    1.592D+13  - 2.722D+13    3.583D+13  - 3.596D+13    2.700D+13  - 1.468D+13    5.458D+12  - 1.240D+12    1.299D+11  
    7.188D+12  - 1.594D+13    2.728D+13  - 3.596D+13    3.614D+13  - 2.717D+13    1.479D+13  - 5.503D+12    1.252D+12  - 1.312D+11  
  - 5.376D+12    1.194D+13  - 2.046D+13    2.700D+13  - 2.717D+13    2.045D+13  - 1.114D+13    4.150D+12  - 9.449D+11    9.916D+10  
    2.912D+12  - 6.474D+12    1.111D+13  - 1.468D+13    1.479D+13  - 1.114D+13    6.078D+12  - 2.266D+12    5.163D+11  - 5.423D+10  
  - 1.078D+12    2.401D+12  - 4.126D+12    5.458D+12  - 5.503D+12    4.150D+12  - 2.266D+12    8.454D+11  - 1.928D+11    2.027D+10  
    2.442D+11  - 5.444D+11    9.366D+11  - 1.240D+12    1.252D+12  - 9.449D+11    5.163D+11  - 1.928D+11    4.400D+10  - 4.629D+09  
  - 2.551D+10    5.691D+10  - 9.802D+10    1.299D+11  - 1.312D+11    9.916D+10  - 5.423D+10    2.027D+10  - 4.629D+09    4.872D+08

remarque · January 2015

Super fiable ! (:D(:D

-->N=20
 N  =
 
    20.  
 
-->
 
-->H=zeros(N,N);
 
-->for i=1:N
-->    for j=1:N
-->        H(i,j)=1/(i+j-1);
-->    end
-->end
 
-->Hmoinsun=testmatrix('hilb',20);
 
-->testultime=H*Hmoinsun;
 
-->norm(testultime-eye(20,20))
 ans  =
 
    2.546D+11

dans l'autre sens, pas mieux.

B=norm(Hmoinsun*H-eye(20,20))
 B  =
 
    2.532D+11

Edit : Je plaisante, mais c'est un peu normal que le fait que ce soit une approximation de l'inverse de la matrice de Hilbert soit impossible à vérifier en flottant et peut-être que $10^{11}$ ce n'est pas si déshonorant que ça. D'autant plus que si l'on regarde dans la matrice inverse, on voit des coefficients en $10^2$ pour les plus petits et en $10^{27}$ pour les plus grands.