Pourquoi les systèmes linéaires mal conditionnés peuvent-ils être résolus avec précision?

13

Selon la réponse ici , un grand nombre de conditions (pour la résolution de systèmes linéaires) diminue le nombre garanti de chiffres corrects dans la solution à virgule flottante. Les matrices de différenciation d'ordre supérieur dans les méthodes pseudospectrales sont généralement très mal conditionnées. Pourquoi est-ce alors que ce sont encore des méthodes très précises?

Je comprends que la faible précision provenant des matrices mal conditionnées n'est qu'une valeur garantie , mais je me demande quand même pourquoi les matrices mal conditionnées sont résolues avec précision par des méthodes directes dans la pratique - par exemple, les LCOLcolonnes du tableau 3.1 à la page 11 de Wang et al., UNE MÉTHODE DE COLLOCATION BIEN CONDITIONNÉE UTILISANT UNE MATRICE D'INTÉGRATION PSEUDOSPECTRALE , SIAM J. Sci. Comput., 36 (3) .

Zoltán Csáti
la source
2
Mon intuition est que la solvabilité / précision d'un système Ax = b est liée au vecteur de forçage b, pas seulement à la matrice A. Peut-être que si b ne "sonde" pas ou "n'excite" pas les modes mal conditionnés de A, alors une solution précise reste possible. À titre d'exemple limitatif, A peut être exactement singulier (nombre de conditions infini), mais Ax = b peut toujours posséder une solution, qui peut être calculée avec précision, si les données de forçage b résident dans la plage de A. J'admets que c'est assez joli -ondulé, c'est pourquoi je commente seulement au lieu de répondre.
rchilton1980
@ rchilton1980 "pourtant Ax = b peut encore posséder une solution" Mais cette solution n'est pas unique. Et les exemples auxquels je fais référence possèdent une solution unique.
Zoltán Csáti
C'est un juste contrepoint - peut-être un artefact de la sélection d'un nombre de conditions infini (une valeur propre exactement nulle). Cependant, je pense que vous pouvez remplacer cette valeur propre nulle par la machine epsilon et mon argument est toujours valable. (Autrement dit, le système a un très grand nombre de conditions, le système est non singulier avec une solution unique, que nous pouvons calculer très précisément à condition que b ne comporte aucun composant le long de cette minuscule paire propre).
rchilton1980
1
Pour être plus précis, mon expérience de pensée est quelque chose comme A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Il est artificiel, mais je pense qu'il suffit de justifier la réclamation originale: "des A parfois mal conditionnés peuvent être résolus avec précision pour des choix particuliers de b"
rchilton1980
1
Donnez un autre exemple sur le blog de Moler blogs.mathworks.com/cleve/2015/02/16/…
percusse

Réponses:

7

Ajouté après ma réponse initiale:

Il me semble maintenant que l'auteur de l'article référencé donne des numéros de condition (apparemment des numéros de condition à 2 normes mais peut-être des numéros de condition à l'infini) dans le tableau tout en donnant un maximum d'erreurs absolues plutôt que des erreurs relatives de norme ou des erreurs relatives maximales par élément ( ce sont toutes des mesures différentes.) Notez que l'erreur relative maximale élément par élément n'est pas la même chose que l'erreur relative de norme infinie. De plus, les erreurs dans le tableau sont relatives à la solution exacte du problème d'origine de la valeur limite de l'équation différentielle plutôt qu'au système linéaire discrétisé d'équations. Ainsi, les informations fournies dans le document ne sont vraiment pas appropriées pour une utilisation avec la borne d'erreur basée sur le numéro de condition.

Cependant, dans ma réplication des calculs, je vois des situations où l'erreur de norme relative à l'infini (ou l'erreur relative à deux normes) est sensiblement plus petite que la limite définie par le nombre de conditions à la norme à l'infini (respectivement le nombre de conditions à 2 normes). Parfois, vous avez juste de la chance.

J'ai utilisé le package DMSUITE MATLAB et résolu l'exemple de problème de cet article en utilisant la méthode pseudospectrale avec les polynômes de Chebyshev. Mes numéros de condition et les erreurs absolues maximales étaient similaires à ceux rapportés dans le document.

J'ai également vu des erreurs relatives à la norme qui étaient un peu mieux que ce à quoi on pourrait s'attendre en fonction du numéro de condition. Par exemple, sur l'exemple de problème avec , en utilisant N = 1024 , j'obtiensϵ=0,01N=1024

cond (A, 2) = 7,9e + 8

cond (A, inf) = 7,8e + 8

norme (u-uexact, 2) / norme (uexact, 2) = 3.1e-12

norme (u-uexact, inf) / norme (uexact, inf) = 2,7e-12

Il semblerait que la solution soit bonne à environ 11-12 chiffres, alors que le numéro de condition est de l'ordre de 1e8.

Cependant, la situation avec des erreurs par élément est plus intéressante.

max (abs (u-uexact)) = 2,7e-12

Cela semble toujours bon.

max (abs ((u-uexact) ./ uexact) = 6,1e + 9

Wow - il y a une très grande erreur relative dans au moins un composant de la solution.

Qu'est-il arrivé? La solution exacte de cette équation a des composants qui sont minuscules (par exemple 1,9e-22) tandis que la solution approximative se termine à une valeur beaucoup plus grande de 9e-14. Ceci est masqué par la mesure d'erreur relative de norme (que ce soit la norme 2 ou la norme infini) et ne devient visible que lorsque vous regardez les erreurs relatives élément par élément et prenez le maximum.

Ma réponse originale ci-dessous explique pourquoi vous pouvez obtenir une erreur relative normale dans la solution qui est inférieure à la limite donnée par le numéro de condition.


Comme vous l'avez noté dans la question, le nombre de conditions, , d'une matrice non singulière donne une erreur relative dans le pire des cas liée à la solution d'un système d'équations perturbé. Autrement dit, si nous résolvons A ( x + Δ x ) = b + Δ b exactement et résolvons A x = b exactement, alorsκ(UNE)UNE(X+ΔX)=b+ΔbUNEX=b

ΔXXκ(UNE)Δbb

Les numéros de condition peuvent être calculés par rapport à une variété de normes, mais le numéro de condition à deux normes est souvent utilisé, et c'est le numéro de condition utilisé dans le document auquel vous vous référez.

Le pire des cas d' erreur se produit lorsque est un vecteur singulier gauche de A correspondant à la plus petite valeur singulière de A . Le meilleur cas se produit lorsque Δ b est un vecteur singulier gauche de A correspondant à la plus grande valeur singulière de A . Lorsque Δ b est aléatoire, alors vous devez regarder les projections de Δ b sur tous les vecteurs singuliers gauches de A et les valeurs singulières correspondantes. Selon le spectre de A , les choses peuvent aller très mal ou très bien. ΔbUNEUNEΔbUNEUNEΔbΔbUNEUNE

Considérons deux matrices , toutes deux avec le numéro de condition à 2 normes 1.0 × 10 10 . La première matrice a les valeurs singulières 1 , 1 × 10 - 10 , , 1 × 10 - 10 . La deuxième matrice a des valeurs singulières 1 , 1 , , 1 , 1 × 10 - 10 . UNE1.0×dixdix11×dix-dix1×dix-dix1111×dix-dix

Dans le premier cas, une perturbation aléatoire est peu susceptible d'être dans la direction du premier vecteur singulier gauche, et plus susceptible d'être proche de l'un des vecteurs singuliers de valeur singulière . Ainsi, le changement relatif dans la solution est susceptible d'être très important. Dans le second cas, presque toute perturbation sera proche de la direction d'un vecteur singulier avec une valeur singulière 1 , et le changement relatif dans la solution sera petit. 1×dix-dix1

PS (ajouté plus tard après mon retour de cours de yoga ...)

La formule pour la solution de estUNEΔX=Δb

Δx=VΣ1UTΔb=i=1nUiTΔbσiVje

Par le théorème de Pythagore,

Δx22=i=1n(UiTΔbσi)2

Si nous gardons , alors cette somme est maximisée lorsque Δ b = U n et minimisée lorsque Δ b = U 1 .Δb2=1Δb=UnΔb=U1

Dans la situation considérée ici, est le résultat d'erreurs d'arrondi aléatoires, donc les valeurs U T i Δ b devraient toutes être à peu près de la même ampleur. Les termes avec des valeurs plus petites de σ i contribueront beaucoup à l'erreur, tandis que les termes avec des valeurs plus grandes de σ i ne contribueront pas beaucoup. Selon le spectre, cela pourrait facilement être beaucoup plus petit que le pire des cas. ΔbUiTΔbσiσi

Brian Borchers
la source
Cet argument n'impliquerait-il pas qu'il est possible (même si cela est peu probable) d'atteindre la pire des limites de pour la matrice dans l'exemple? AFAIU, sur la base de ma réponse et de la documentation de, cela ne devrait pas être possible. κ(UNE)?getrs
Kirill
@BrianBorchers Pourriez-vous expliquer pourquoi "La pire des erreurs se produit lorsque est un vecteur singulier gauche de A correspondant à la plus petite valeur singulière de A. Le meilleur des cas se produit lorsque Δ b est un vecteur singulier gauche de A correspondant au plus grand valeur singulière de A. " tient? Dans l'exemple ci-dessous, c'est logique, mais j'aurais besoin de quelques formules. Laissez le SVD A soit A = U Σ V T . Dans le premier cas, A = Δ b σ 1 v T 1 +ΔbAAΔbAAAA=UΣVT . La façon de procéder? A=Δbσ1v1T+i=2NuiσiviT
Zoltán Csáti
Je n'ai pas discuté des erreurs d'arrondi dans la matrice , mais l'effet général est similaire - à moins que vous n'ayez vraiment de malchance dans les erreurs d'arrondi, vous faites généralement un peu mieux que le pire des cas pessimiste. A
Brian Borchers
(-1) La discussion des erreurs relatives par composant dans la sortie est gravement trompeuse.
Kirill
1

tl; dr Ils ont signalé un numéro de condition, pas nécessairement le bon numéro de condition pour la matrice, car il y a une différence.

Ceci est spécifique à la matrice et au vecteur de droite. Si vous regardez la documentation de*getrs , elle indique que la limite d'erreur directe est

xx0xcond(A,x)ucond(A)u.
Here cond(A,x) is not quite the usual condition number κ(A), but rather
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.
(Here inside the norm these are component-wise absolute values.) See, for example, Iterative refinement for linear systems and LAPACK by Higham, or Higham's Accuracy and Stability of Numerical Algorithms (7.2).

For your example, I took a pseudospectral differential operator for a similar problem with n=128, and there is in fact a big difference between |A1||A| and κ(A), I computed 7×103 and 2.6×107, which is enough to explain the observation that this happens for all right hand sides, because the orders of magnitudes roughly match what is seen in Table 3.1 (3-4 orders better errors). This doesn't work when I try the same for just a random ill-conditioned matrix, so it has to be a property of A.

An explicit example for which the two condition numbers don't match, which I took from Higham (7.17, p.124), due to Kahan is

(2111ϵϵ1ϵϵ),(2+2ϵϵϵ).
Another example I found is just the plain Vandermonde matrix on [1:10] with random b. I went through MatrixDepot.jl and some other ill-conditioned matrices also produce this type of result, like triw and moler.

Essentially, what's going on is that when you analyze the stability of solving linear systems with respect to perturbations, you first have to specify which perturbations you are considering. When solving linear systems with LAPACK, this error bound considers component-wise perturbations in A, but no perturbation in b. So this is different from the usual κ(A)=A1A, which considers normwise perturbations in both A and b.

Consider (as a counterexample) also what would happen if you don't make the distinction. We know that using iterative refinement with double precision (see link above) we can get the best possible forward relative error of O(u) for those matrices with κ(A)1/u. So if we consider the idea that linear systems can't be solved to accuracy better than κ(A)u, how would refining solutions possibly work?

P.S. It matters that ?getrs says the computed solution is the true solution of (A + E)x = b with a perturbation E in A, but no perturbation in b. Things would be different if perturbations were allowed in b.

Edit To show this working more directly, in code, that this is not a fluke or a matter of luck, but rather the (unusual) consequence of two condition numbers being very different for some specific matrices, i.e.,

cond(A,x)cond(A)κ(A).
function main2(m=128)
    A = matrixdepot("chebspec", m)^2
    A[1,:] = A[end,:] = 0
    A[1,1] = A[end,end] = 1
    best, worst = Inf, -Inf
    for k=1:2^5
        b = randn(m)
        x = A \ b
        x_exact = Float64.(big.(A) \ big.(b))
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        best, worst = min(best, err), max(worst, err)
    end
    @printf "Best relative error:       %.3e\n" best
    @printf "Worst relative error:      %.3e\n" worst
    @printf "Predicted error κ(A)*ε:    %.3e\n" cond(A, Inf)*eps()
    @printf "Predicted error cond(A)*ε: %.3e\n" norm(abs.(inv(A))*abs.(A), Inf)*eps()
end

julia> main2()
Best relative error:       2.156e-14
Worst relative error:      2.414e-12
Predicted error κ(A)*ε:    8.780e-09
Predicted error cond(A)*ε: 2.482e-12

Edit 2 Here is another example of the same phenomenon where the different conditions numbers unexpectedly differ by a lot. This time,

cond(A,x)cond(A)κ(A).
Here A is the 10×10 Vandermonde matrix on 1:10, and when x is chosen randomly, cond(A,x) is noticably smaller than κ(A), and the worst case x is given by xi=ia for some a.
function main4(m=10)
    A = matrixdepot("vand", m)
    lu = lufact(A)
    lu_big = lufact(big.(A))
    AA = abs.(inv(A))*abs.(A)
    for k=1:12
        # b = randn(m) # good case
        b = (1:m).^(k-1) # worst case
        x, x_exact = lu \ b, lu_big \ big.(b)
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        predicted = norm(AA*abs.(x), Inf)/norm(x, Inf)*eps()
        @printf "relative error[%2d]    = %.3e (predicted cond(A,x)*ε = %.3e)\n" k err predicted
    end
    @printf "predicted κ(A)*ε      = %.3e\n" cond(A)*eps()
    @printf "predicted cond(A)*ε   = %.3e\n" norm(AA, Inf)*eps()
end

Average case (almost 9 orders of magnitude better error):

julia> T.main4()
relative error[1]     = 6.690e-11 (predicted cond(A,x)*ε = 2.213e-10)
relative error[2]     = 6.202e-11 (predicted cond(A,x)*ε = 2.081e-10)
relative error[3]     = 2.975e-11 (predicted cond(A,x)*ε = 1.113e-10)
relative error[4]     = 1.245e-11 (predicted cond(A,x)*ε = 6.126e-11)
relative error[5]     = 4.820e-12 (predicted cond(A,x)*ε = 3.489e-11)
relative error[6]     = 1.537e-12 (predicted cond(A,x)*ε = 1.729e-11)
relative error[7]     = 4.885e-13 (predicted cond(A,x)*ε = 8.696e-12)
relative error[8]     = 1.565e-13 (predicted cond(A,x)*ε = 4.446e-12)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Worst case (a=1,,12):

julia> T.main4()
relative error[ 1]    = 0.000e+00 (predicted cond(A,x)*ε = 6.608e-13)
relative error[ 2]    = 1.265e-13 (predicted cond(A,x)*ε = 3.382e-12)
relative error[ 3]    = 5.647e-13 (predicted cond(A,x)*ε = 1.887e-11)
relative error[ 4]    = 8.895e-74 (predicted cond(A,x)*ε = 1.127e-10)
relative error[ 5]    = 4.199e-10 (predicted cond(A,x)*ε = 7.111e-10)
relative error[ 6]    = 7.815e-10 (predicted cond(A,x)*ε = 4.703e-09)
relative error[ 7]    = 8.358e-09 (predicted cond(A,x)*ε = 3.239e-08)
relative error[ 8]    = 1.174e-07 (predicted cond(A,x)*ε = 2.310e-07)
relative error[ 9]    = 3.083e-06 (predicted cond(A,x)*ε = 1.700e-06)
relative error[10]    = 1.287e-05 (predicted cond(A,x)*ε = 1.286e-05)
relative error[11]    = 3.760e-10 (predicted cond(A,x)*ε = 1.580e-09)
relative error[12]    = 3.903e-10 (predicted cond(A,x)*ε = 1.406e-09)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Edit 3 Another example is the Forsythe matrix, which is a perturbed Jordan block of any size of the form

A=(010000100001ϵ000).
This has A=1, A1=ϵ1, so κ(A)=ϵ1, but |A1|=A1=|A|1, so cond(A)=1. And as can be verified by hand, solving systems of linear equations like Ax=b with pivoting is extremely accurate, despite the potentially unbounded κ(A). So this matrix too will yield unexpectedly precise solutions.

Edit 4 Kahan matrices are also like this, with cond(A)κ(A):

A = matrixdepot("kahan", 48)
κ, c = cond(A, Inf), norm(abs.(inv(A))*abs.(A), Inf)
@printf "κ=%.3e c=%.3e ratio=%g\n" κ c (c/κ)

κ=8.504e+08 c=4.099e+06 ratio=0.00482027
Kirill
la source
The condition numbers in the paper referred to by the OP are two-norm condition numbers. If you go back to reference [17] by ElBarbary you'll see that in the earlier paper these were two-norm condition numbers. Also, I setup the examples from this paper using DMsuite, and got nearly exactly the same 2-norm condition numbers as reported in the paper.
Brian Borchers
The infinity norm condition norm numbers for these examples that I got using dmsuite and Chebyshev interpolation were similar in magnitude to the two-norm condition numbers. I don't think that the difference between 2-norm in infinity-norm condition numbers is that important for this particular example.
Brian Borchers
I believe that the errors reported in the paper are absolute rather than relative errors (it doesn't make too much difference except for ϵ=0.01, where the solution dips down close to 0.
Brian Borchers
For ϵ=0.01 and N=1024, the relative errors for the parts of the solution that are near 0 are huge, but the absolute errors are small. I agree that the paper was very vague about what condition number was used and about what the "errors" were exactly (relative or absolute errors.)
Brian Borchers
@BrianBorchers I'm not sure what you mean: this isn't the difference between 2-norm and infty-norm condition numbers, but rather normwise- and component-wise condition numbers (component-wise relative perturbations in the input, not component-wise relative errors in the output as in your answer).
Kirill