¿Por qué los sistemas lineales mal condicionados se pueden resolver con precisión?

13

Según la respuesta aquí , un número de condición grande (para la resolución lineal del sistema) disminuye el número garantizado de dígitos correctos en la solución de coma flotante. Las matrices de diferenciación de orden superior en los métodos pseudoespectrales suelen estar muy mal condicionadas. ¿Por qué es entonces que siguen siendo métodos muy precisos?

Entiendo que la baja precisión que proviene de las matrices mal acondicionadas es solo un valor garantizado , pero aún así me pregunto por qué las matrices mal acondicionadas se resuelven con precisión por métodos directos en la práctica, por ejemplo, las LCOLcolumnas de la Tabla 3.1 en la página 11 de Wang et al., MÉTODO DE COLOCACIÓN BIEN CONDICIONADO UTILIZANDO UNA MATRIZ DE INTEGRACIÓN PSEUDOSPECTRAL , SIAM J. Sci. Comput., 36 (3) .

Zoltán Csáti
fuente
2
Mi intuición es que la solubilidad / precisión de un sistema Ax = b está vinculada al vector de fuerza b, no solo a la matriz A. Quizás si b no "sondea" o "excita" los modos mal condicionados de A, entonces la solución precisa sigue siendo posible. Como ejemplo limitante, A puede ser exactamente singular (número de condición infinita), pero Ax = b aún puede tener una solución, que se puede calcular con precisión, si los datos de forzamiento b residen en el rango de A. Admito que esto es bastante útil. -ondulado, por eso solo comento en lugar de responder.
rchilton1980
@ rchilton1980 "aún Ax = b puede poseer una solución" Pero esa solución no es única. Y los ejemplos a los que me refiero poseen una solución única.
Zoltán Csáti
Ese es un contrapunto justo, quizás un artefacto de elegir un número de condición infinita (un valor propio exactamente cero). Sin embargo, creo que puede reemplazar ese valor propio cero con epsilon de máquina y mi punto sigue en pie. (Es decir, el sistema tiene un número de condición muy grande, el sistema no es singular con una solución única, que podemos calcular con mucha precisión siempre que b no tenga ningún componente a lo largo de ese pequeño par propio).
rchilton1980
1
Para ser más específico, mi experimento mental aquí es algo así como A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Está inventado, pero creo que es suficiente para justificar la afirmación original: "a veces las A mal condicionadas pueden resolverse con precisión para elecciones particulares de b"
rchilton1980
1
Sólo dar otro ejemplo del blog de Moler blogs.mathworks.com/cleve/2015/02/16/...
percusse

Respuestas:

7

Agregado después de mi respuesta inicial:

Ahora me parece que el autor del artículo al que se hace referencia está dando números de condición (aparentemente números de condición de 2 normas pero posiblemente números de condición de norma infinita) en la tabla al tiempo que proporciona errores absolutos máximos en lugar de errores relativos normativos o errores relativos máximos relativos a elementos ( estas son todas medidas diferentes.) Tenga en cuenta que el error relativo máximo por elemento no es lo mismo que el error relativo de la norma infinita. Además, los errores en la tabla son relativos a la solución exacta al problema del valor límite original de la ecuación diferencial más que al sistema lineal de ecuaciones discretizado. Por lo tanto, la información proporcionada en el documento realmente no es adecuada para su uso con el error vinculado en función del número de condición.

Sin embargo, en mi réplica de los cálculos, veo situaciones en las que el error de la norma de infinito relativo (o error relativo de dos normas) es sustancialmente menor que el límite establecido por el número de condición de norma de infinito (número de condición de 2 normas respectivamente). A veces solo tienes suerte.

Utilicé el paquete DMSUITE MATLAB y resolví el problema de ejemplo de este artículo utilizando el método pseudoespecífico con polinomios de Chebyshev. Mis números de condición y mis errores absolutos máximos fueron similares a los reportados en el documento.

También vi errores relativos a la norma que fueron algo mejores de lo que cabría esperar según el número de condición. Por ejemplo, en el problema de ejemplo con , usando N = 1024 , obtengoϵ=0.01N=1024

cond (A, 2) = 7.9e + 8

cond (A, inf) = 7.8e + 8

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

norma (u-uexact, inf) / norma (uexact, inf) = 2.7e-12

Parece que la solución es buena para aproximadamente 11-12 dígitos, mientras que el número de condición es del orden de 1e8.

Sin embargo, la situación con errores de elementos sabios es más interesante.

max (abs (u-uexact)) = 2.7e-12

Eso todavía se ve bien.

max (abs ((u-uexact) ./ uexact) = 6.1e + 9

Wow: hay un error relativo muy grande en al menos un componente de la solución.

¿Que pasó? La solución exacta de esta ecuación tiene componentes que son pequeños (por ejemplo, 1.9e-22) mientras que la solución aproximada toca fondo con un valor mucho mayor de 9e-14. Esto está oculto por la medida de error relativo de la norma (ya sea la norma 2 o la norma de infinito) y solo se hace visible cuando se miran los errores relativos de elementos y se toma el máximo.

Mi respuesta original a continuación explica por qué puede obtener un error relativo de la norma en la solución que es menor que el límite dado por el número de condición.


Como ha notado en la pregunta, el número de condición, , de una matriz no singular da un error relativo en el peor de los casos vinculado a la solución de un sistema perturbado de ecuaciones. Es decir, si resolvemos A ( x + Δ x ) = b + Δ b exactamente y resolvemos A x = b exactamente, entoncesκ(A)A(x+Δx)=b+ΔbAx=b

Δxxκ(A)Δbb

Los números de condición se pueden calcular con respecto a una variedad de normas, pero a menudo se usa el número de condición de dos normas, y ese es el número de condición utilizado en el documento al que se refiere.

El error del peor caso se produce cuando es un vector singular izquierda de A correspondiente al valor singular más pequeño de A . El mejor de los casos se produce cuando Δ b es un vector singular izquierda de A correspondiente al valor singular más grande de A . Cuando Δ b es aleatorio, debe observar las proyecciones de Δ b en todos los vectores singulares izquierdos de A y los valores singulares correspondientes. Dependiendo del espectro de A , las cosas pueden ir muy mal o muy bien. ΔbAAΔbAAΔbΔbAA

Considere dos matrices , ambas con el número de condición de 2 normas 1.0 × 10 10 . La primera matriz tiene los valores singulares 1 , 1 × 10 - 10 , , 1 × 10 - 10 . La segunda matriz tiene valores singulares 1 , 1 , ... , 1 , 1 × 10 - 10 . A1.0×101011×10101×10101111×1010

En el primer caso, es poco probable que una perturbación aleatoria esté en la dirección del primer vector singular izquierdo, y es más probable que esté cerca de uno de los vectores singulares con valor singular . Por lo tanto, el cambio relativo en la solución es probable que sea muy grande. En el segundo caso, casi cualquier perturbación estará cerca en dirección a un vector singular con valor singular 1 , y el cambio relativo en la solución será pequeño. 1×10101

PD (agregado más tarde después de regresar de la clase de yoga ...)

La fórmula para la solución de esAΔx=Δb

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

Por el teorema de Pitágoras,

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

Si mantenemos , entonces esta suma se maximiza cuando Δ b = U n y se minimiza cuando Δ b = U 1 .Δb2=1Δb=UnΔb=U1

En la situación considerada aquí, es el resultado de errores de redondeo aleatorio, por lo que los valores de U T i Δ b deben ser aproximadamente de la misma magnitud. Los términos con valores más pequeños de σ i contribuirán mucho al error, mientras que los términos con valores más grandes de σ i no contribuirán mucho. Dependiendo del espectro, esto podría ser mucho más pequeño que el peor de los casos. ΔbUiTΔbσiσi

Brian Borchers
fuente
¿No implicaría este argumento que es posible (aunque poco probable) alcanzar el límite de en el peor de los casos para la matriz en el ejemplo? AFAIU, según mi respuesta y según la documentación de esto, esto no debería ser posible. κ(A)?getrs
Kirill
@BrianBorchers ¿Podría explicar por qué "El peor de los casos ocurre cuando es un vector singular izquierdo de A correspondiente al valor singular más pequeño de A. El mejor caso ocurre cuando Δ b es un vector singular izquierdo de A correspondiente al valor más grande valor singular de A ". sostiene? Del ejemplo a continuación es lógico, pero necesitaría algunas fórmulas. Deje que la SVD de A sea A = U Σ V T . En el primer caso, A = Δ b σ 1 v T 1 +ΔbAAΔbAAAA=UΣVT . ¿Cómo proceder? A=Δbσ1v1T+i=2NuiσiviT
Zoltán Csáti
No he discutido los errores de redondeo en la matriz , pero el efecto general es similar: a menos que tenga muy mala suerte en los errores de redondeo, generalmente le va un poco mejor que el pesimista límite del peor de los casos. A
Brian Borchers
(-1) La discusión de los errores relativos por componentes en la salida es muy engañosa.
Kirill
1

tl; dr Informaron un número de condición, no necesariamente el número de condición correcto para la matriz, porque hay una diferencia.

Esto es específico de la matriz y del vector del lado derecho. Si nos fijamos en la documentación*getrs , dice que el error de reenvío es Aquícond(A,x)no es exactamente el número de condición habitualκ(A), sino más bien cond(A,x)=| A - 1

xx0xcond(A,x)ucond(A)u.
cond(A,x)κ(A) (Aquí dentro de la norma, estos son valores absolutos por componentes). Ver, por ejemplo,refinamiento iterativo para sistemas lineales y LAPACKpor Higham, o laprecisión y estabilidad de algoritmos numéricos deHigham(7.2).
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.

Para su ejemplo, tomé un operador diferencial pseudoespectral para un problema similar con , y de hecho hay una gran diferencia entre | A - 1 | El | A | Y κ ( A ) , calculé 7 × 10 3 y 2.6 × 10 7n=128|A1||A|κ(A)7×1032.6×107, que es suficiente para explicar la observación de que esto sucede para todos los lados derechos, porque los órdenes de magnitud coinciden aproximadamente con lo que se ve en la Tabla 3.1 (3-4 órdenes de mejores errores). Esto no funciona cuando intento la misma por sólo una matriz mal condicionado al azar, por lo que tiene que ser una propiedad de .A

Un ejemplo explícito para el que los dos números de condición no coinciden, que tomé de Higham (7.17, p.124), debido a Kahan es Otro ejemplo que encontré es solo la matriz de Vandermonde simpleconbaleatorio. Paséy algunas otras matrices mal acondicionadas también producen este tipo de resultado, comoy.

(2111ϵϵ1ϵϵ),(2+2ϵϵϵ).
[1:10]bMatrixDepot.jltriwmoler

Esencialmente, lo que está sucediendo es que cuando analiza la estabilidad de resolver sistemas lineales con respecto a las perturbaciones, primero debe especificar qué perturbaciones está considerando. Al resolver sistemas lineales con LAPACK, este error considera las perturbaciones de componentes en , pero no las perturbaciones en b . Por lo tanto, esto es diferente de lo habitual κ ( A ) = A - 1A , que considera las perturbaciones normales en A y b .Abκ(A)=A1AAb

O(u)κ(A)1/uκ(A)u

?getrs(A + E)x = bEAbb

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

cond(A,x)cond(A)κ(A).
A1:10xcond(A,x)κ(A)xxi=iaa
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

Caso promedio (casi 9 órdenes de magnitud mejor 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

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

A=(010000100001ϵ000).
A=1A1=ϵ1κ(A)=ϵ1|A1|=A1=|A|1cond(A)=1Ax=bκ(A)

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
fuente
Los números de condición en el documento al que hace referencia el OP son números de condición de dos normas. Si vuelve a la referencia [17] de ElBarbary, verá que en el documento anterior se trataba de números de condición de dos normas. Además, configuré los ejemplos de este documento usando DMsuite y obtuve casi exactamente los mismos números de condición de 2 normas que se informan en el documento.
Brian Borchers
Los números de norma de condición de norma infinita para estos ejemplos que obtuve usando dmsuite y la interpolación de Chebyshev fueron similares en magnitud a los números de condición de dos normas. No creo que la diferencia entre 2-norma en números de condición de norma infinita sea tan importante para este ejemplo en particular.
Brian Borchers
Creo que los errores reportados en el documento son absolutos en lugar de errores relativos (no hace mucha diferencia excepto por ϵ=0,01, donde la solución se sumerge cerca de 0.
Brian Borchers
por ϵ=0,01 y norte=1024, los errores relativos para las partes de la solución que están cerca de 0 son enormes, pero los errores absolutos son pequeños. Estoy de acuerdo en que el documento fue muy vago sobre qué número de condición se usó y sobre cuáles fueron los "errores" exactamente (errores relativos o absolutos).
Brian Borchers
@BrianBorchers No estoy seguro de lo que quieres decir: esta no es la diferencia entre los números de condición de 2 normas y de infty-norm, sino más bien los números de condición de norma y componente (perturbaciones relativas de componente en la entrada, no componente -como errores relativos en la salida como en su respuesta).
Kirill