Por que sistemas lineares mal condicionados podem ser resolvidos com precisão?

13

De acordo com a resposta aqui , o número de condição grande (para solução linear do sistema) diminui o número garantido de dígitos corretos na solução de ponto flutuante. Matrizes de diferenciação de ordem superior em métodos pseudoespectrais são tipicamente muito condicionadas. Por que é que eles ainda são métodos muito precisos?

Entendo que a baixa precisão proveniente das matrizes mal condicionadas é apenas um valor garantido , mas ainda me faz pensar por que as matrizes mal condicionadas são resolvidas com precisão por métodos diretos na prática - por exemplo, as LCOLcolunas da Tabela 3.1 na página 11 de Wang et al., MÉTODO DE COLOCAÇÃO BEM CONDICIONADO USANDO UMA MATRIZ DE INTEGRAÇÃO PSEUDOSPECTRAL , SIAM J. Sci. Comput., 36 (3) .

Zoltán Csáti
fonte
2
Minha intuição é que a resolubilidade / precisão de um sistema Ax = b está ligada ao vetor forçador b, não apenas à matriz A. Talvez se b não "sondar" ou "excitar" os modos mal condicionados de A, então a solução precisa permanece possível. Como um exemplo limitante, A pode ser exatamente singular (número de condição infinito), mas Ax = b ainda pode possuir uma solução, que pode ser calculada com precisão, se os dados forçantes b residirem no intervalo de A. Admito que isso é bastante útil -wavy, e é por isso que eu apenas comento em vez de responder.
precisa saber é o seguinte
@ rchilton1980 "ainda Ax = b ainda pode ter uma solução" Mas essa solução não é única. E os exemplos que estou me referindo possuem uma solução única.
Zoltán Csáti
Esse é um contraponto justo - talvez um artefato para escolher um número de condição infinito (um valor próprio exatamente zero). No entanto, acho que você pode substituir esse valor próprio zero pela máquina epsilon e meu argumento ainda permanece. (Ou seja, o sistema tem um número de condição muito grande, o sistema não é singular com uma solução única, que podemos calcular com muita precisão, desde que b não tenha nenhum componente nesse minúsculo par próprio).
rchilton1980
1
Para ser mais específico, meu experimento mental aqui é algo como A = diag ([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. Ele é planejado, mas eu acho que é suffcient para justificar a reivindicação original: "às vezes mal condicionado de um pode ser resolvido com precisão para escolhas particulares de b"
rchilton1980
1
Basta dar outro exemplo do blog de Moler blogs.mathworks.com/cleve/2015/02/16/...
percusse

Respostas:

7

Adicionado após minha resposta inicial:

Parece-me agora que o autor do artigo mencionado está fornecendo números de condição (aparentemente números de condição de 2 normas, mas possivelmente números de condição de infinito) na tabela, ao mesmo tempo em que fornece erros absolutos máximos em vez de erros relativos à norma ou erros relativos ao elemento ( todas são medidas diferentes.) Observe que o erro relativo máximo elemento a elemento não é o mesmo que o erro relativo de norma infinita. Além disso, os erros na tabela são relativos à solução exata do problema original do valor-limite da equação diferencial e não ao sistema linear discretizado de equações. Portanto, as informações fornecidas no documento realmente não são apropriadas para uso com o erro vinculado com base no número da condição.

No entanto, em minha replicação dos cálculos, vejo situações em que o erro relativo da norma infinito (ou erro relativo de duas normas) é substancialmente menor que o limite definido pelo número da condição de infinito-norma (número da condição de 2 normas, respectivamente). Às vezes você apenas tem sorte.

Usei o pacote DMSUITE MATLAB e resolvi o problema de exemplo deste artigo usando o método pseudoespectral com polinômios de Chebyshev. Meus números de condição e erros absolutos máximos foram semelhantes aos relatados no artigo.

Também vi erros relativos à norma que eram um pouco melhores do que se poderia esperar com base no número da condição. Por exemplo, no problema de exemplo com , usando N = 1024 , receboϵ=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 a solução é boa para cerca de 11 a 12 dígitos, enquanto o número da condição é da ordem de 1e8.

No entanto, a situação com erros por elementos é mais interessante.

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

Ainda parece bom.

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

Uau - há um erro relativo muito grande em pelo menos um componente da solução.

O que aconteceu? A solução exata dessa equação possui componentes que são minúsculos (por exemplo, 1,9e-22), enquanto a solução aproximada atinge um valor muito maior de 9e-14. Isso é oculto pela medição de erro relativo da norma (seja a norma 2 ou infinito) e só se torna visível quando você olha para os erros relativos aos elementos e obtém o máximo.

Minha resposta original abaixo explica por que você pode obter um erro relativo da norma na solução que é menor que o limite fornecido pelo número da condição.


Como você observou na pergunta, o número da condição, , de uma matriz não singular fornece um erro relativo de pior caso vinculado à solução para um sistema perturbado de equações. Ou seja, se resolvermos A ( x + Δ x ) = b + Δ b exatamente e resolver A x = b exatamente, entãoκ(A)A(x+Δx)=b+ΔbAx=b

Δxxκ(A)Δbb

Os números de condição podem ser calculados com relação a uma variedade de normas, mas o número de condição com duas normas é frequentemente usado, e esse é o número de condição usado no artigo a que você se refere.

O pior caso de erro ocorre quando é um vector singular esquerda de um correspondente para o menor valor singular de um . O melhor caso ocorre quando Δ b é um vector singular esquerda de uma correspondente ao maior valor singular de um . Quando Δ b é aleatório, é necessário observar as projeções de Δ b em todos os vetores singulares esquerdos de A e nos valores singulares correspondentes. Dependendo do espectro de A , as coisas podem correr muito mal ou muito bem. ΔbAAΔbAAΔbΔbAA

Considere duas matrizes , ambas com número de condição de 2 normas 1.0 × 10 10 . A primeira matriz possui os valores singulares 1 , 1 × 10 - 10 , , 1 × 10 - 10 . A segunda matriz possui valores singulares 1 , 1 , , 1 , 1 × 10 - 10 . A1.0×101011×10101×10101111×1010

No primeiro caso, é improvável que uma perturbação aleatória esteja na direção do primeiro vetor singular esquerdo e seja mais provável que esteja perto de um dos vetores singulares com valor singular . Portanto, a mudança relativa na solução provavelmente será muito grande. No segundo caso, quase qualquer perturbação estará próxima em direção a um vetor singular com valor singular 11×10101 , e a mudança relativa na solução será pequena.

PS (adicionado mais tarde depois que voltei da aula de ioga ...)

A fórmula da solução para éAΔx=Δb

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

Pelo teorema de Pitágoras,

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

Se continuarmos , então esta soma é maximizado quando Δ b = U n e minimizado quando Δ b = L 1 .Δb2=1Δb=UnΔb=U1

Na situação considerada aqui, é o resultado de erros de arredondamento aleatórios, portanto, os valores de U T i Δ b devem ser aproximadamente da mesma magnitude. Os termos com valores menores de σ i contribuirão muito para o erro, enquanto os termos com valores maiores de σ i não contribuirão muito. Dependendo do espectro, isso pode ser facilmente muito menor que o pior caso vinculado. ΔbUiTΔbσiσi

Brian Borchers
fonte
Esse argumento não implica que é possível (mesmo que improvável) alcançar o limite de pior caso de para a matriz no exemplo? AFAIU, com base na minha resposta e na documentação da empresa, isso não deve ser possível. κ(A)?getrs
Kirill
@BrianBorchers Você poderia explicar por que "O pior erro ocorre quando é um vetor singular esquerdo de A correspondente ao menor valor singular de A. O melhor caso ocorre quando Δ b é um vetor singular esquerdo de A correspondente ao maior valor singular de A. " detém? No exemplo abaixo, é lógico, mas eu precisaria de algumas fórmulas. Deixe o SVD Um estar Um = L Σ V T . No primeiro caso, A = Δ b σ 1 v T 1ΔbAAΔbAAAA=UΣVT . Como proceder? A=Δbσ1v1T+i=2NuiσiviT
Zoltán Csáti
Eu não discuti erros de arredondamento na matriz , mas o efeito geral é semelhante - a menos que você tenha realmente azar nos erros de arredondamento, você normalmente se sai um pouco melhor do que o limite pessimista do pior caso. A
Brian Borchers
(-1) A discussão de erros relativos em termos de componentes na saída é seriamente enganosa.
Kirill #
1

tl; dr Eles relataram um número de condição, não necessariamente o número de condição correto para a matriz, porque há uma diferença.

Isso é específico para a matriz e o vetor do lado direito. Se você olhar para a documentação*getrs , ele diz que o erro para a frente ligado é Aquicond(A,x)não é o habitual condição de númerok(A), mas sim cond(A,x)=| A - 1

xx0xcond(A,x)ucond(A)u.
cond(A,x)κ(A) (Aqui, dentro da norma, esses são valores absolutos em termos de componentes.) Veja, por exemplo, orefinamento iterativo para sistemas lineares e o LAPACKby Higham, ou aPrecisão e Estabilidade de Algoritmos Numéricos deHigham.
cond(A,x)=|A1||A||x|x,cond(A)=|A1||A|.
(7.2).

Para o seu exemplo, usei um operador diferencial pseudoespectral para um problema semelhante com e, de fato, há uma grande diferença entre | A - 1 | | Um | E κ ( A ) , calculei 7 × 10 3 e 2,6 × 10 7n=128|A1||A|κ(A)7×1032.6×107, o que é suficiente para explicar a observação de que isso acontece para todos os lados do lado direito, porque as ordens de magnitudes correspondem aproximadamente ao que é visto na Tabela 3.1 (3-4 ordena erros melhores). Isso não funciona quando tento o mesmo para apenas uma matriz aleatória e mal condicionada; portanto, ela deve ser uma propriedade de A .

Um exemplo explícito para o qual os dois números de condição não coincidem, que tirei de Higham (7.17, p.124), devido a Kahan é Outro exemplo que encontrei é apenas a matriz de Vandermonde simples,combaleatório. Passeie algumas outras matrizes mal condicionadas também produzem esse tipo de resultado, comoe

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

Essencialmente, o que está acontecendo é que, quando você analisa a estabilidade da solução de sistemas lineares em relação a perturbações, primeiro precisa especificar quais perturbações está considerando. Ao resolver sistemas lineares com LAPACK, esse erro associado considera perturbações em componentes em , mas nenhuma perturbação em b . Portanto, isso é diferente do habitual κ ( A ) = " A - 1 " " A " , que considera perturbações normwise em A e bAbκ(A)=A1AAb .

Considere (como um contra-exemplo) também o que aconteceria se você não fizer a distinção. Sabemos que usando o refinamento iterativo com precisão dupla (veja o link acima), podemos obter o melhor erro relativo possível de para as matrizes com κ ( A ) 1 / u . Portanto, se considerarmos a idéia de que sistemas lineares não podem ser resolvidos com precisão melhor que κ ( A ) u , como as soluções de refino possivelmente funcionariam?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

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
fonte
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
Os números de condição padrão infinito para esses exemplos que obtive usando a interpolação dmsuite e Chebyshev foram semelhantes em magnitude aos números de condição com duas normas. Eu não acho que a diferença entre 2-norma nos números de condição de infinito-norma seja tão importante para este exemplo em particular.
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