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

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) .

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$\epsilon=0.01$$N=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$\kappa(A)$$A(x+\Delta x)=b+\Delta b$$Ax=b$

$\frac{\| \Delta x \|}{\| x \|} \leq \kappa(A) \frac{\| \Delta b \|}{\| b \|}$

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. $\Delta b$$A$$A$$\Delta b$$A$$A$$\Delta b$$\Delta b$$A$$A$

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 . $A$$1.0 \times 10^{10}$$1$$1 \times 10^{-10}$$\ldots$$1 \times 10^{-10}$$1$$1$$\ldots$$1$$1 \times 10^{-10}$

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 \times 10^{-10}$$1$

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

La formule pour la solution de est$A\Delta x = \Delta b$

$\Delta x = V \Sigma^{-1} U^{T} \Delta b=\sum_{i=1}^{n} \frac{U_{i}^{T}\Delta b}{\sigma_{i}} V_{i}$

Par le théorème de Pythagore,

$\| \Delta x \|_{2}^{2}= \sum_{i=1}^{n} \left( \frac{U_{i}^{T}\Delta b}{\sigma_{i}} \right)^{2}$

Si nous gardons , alors cette somme est maximisée lorsque Δ b = U n et minimisée lorsque Δ b = U 1 .$\| \Delta b \|_{2}=1$$\Delta b=U_{n}$$\Delta b=U_{1}$

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. $\Delta b$$U_{i}^{T}\Delta b$$\sigma_{i}$$\sigma_{i}$

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

$$\frac{\|x-x_0\|_\infty}{\|x\|_\infty} \lesssim \mathrm{cond}(A,x)u \leq \mathrm{cond}(A)u.$$ Here $\mathrm{cond}(A,x)$ is not quite the usual condition number $\kappa_\infty(A)$, but rather $$\mathrm{cond}(A,x) = \frac{\||A^{-1}||A||x|\|_\infty}{\|x\|_\infty},\\ \mathrm{cond}(A) = \||A^{-1}||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 $\||A^{-1}||A|\|$ and $\kappa_\infty(A)$, I computed $7\times 10^3$ and $2.6\times 10^7$, 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

$$\begin{pmatrix}2&-1&1\\-1&\epsilon&\epsilon\\1&\epsilon&\epsilon\end{pmatrix}, \qquad \begin{pmatrix}2+2\epsilon\\-\epsilon\\\epsilon\end{pmatrix}.$$ 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 $\kappa(A) = \|A^{-1}\|\|A\|$, 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 $\kappa(A)\ll 1/u$. So if we consider the idea that linear systems can't be solved to accuracy better than $\kappa(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.,

$$\mathrm{cond}(A,x) \approx \mathrm{cond}(A) \ll \kappa(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,

$$\mathrm{cond}(A, x) \ll \mathrm{cond}(A) \approx \kappa(A).$$ Here $A$ is the 10×10 Vandermonde matrix on $1:10$, and when $x$ is chosen randomly, $\mathrm{cond}(A,x)$ is noticably smaller than $\kappa(A)$, and the worst case $x$ is given by $x_i = i^a$ 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,\ldots,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 = \begin{pmatrix} 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \\ \epsilon&0&0&0 \end{pmatrix}.$$ This has $\|A\|=1$, $\|A^{-1}\|=\epsilon^{-1}$, so $\kappa_\infty(A) = \epsilon^{-1}$, but $|A^{-1}| = A^{-1} = |A|^{-1}$, so $\mathrm{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 $\kappa_\infty(A)$. So this matrix too will yield unexpectedly precise solutions.

Edit 4 Kahan matrices are also like this, with $\mathrm{cond}(A)\ll \kappa(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