Para dados ruidosos ou com estrutura fina, existem quadraturas melhores que a regra do ponto médio?

12

Apenas as duas primeiras seções desta longa questão são essenciais. Os outros são apenas para ilustração.

fundo

Quadraturas avançadas, como Newton-Cotes, Gauß-Legendre e Romberg, parecem ser principalmente destinadas a casos em que é possível provar a função com precisão, mas não integrar analiticamente. No entanto, para funções com estruturas mais finas que o intervalo de amostragem (consulte o Apêndice A, por exemplo) ou ruído de medição, elas não podem competir com abordagens simples, como a regra do ponto médio ou trapézio (consulte o Apêndice B para uma demonstração).

Isso é um pouco intuitivo, pois, por exemplo, a regra composta de Simpson essencialmente "descarta" um quarto da informação, atribuindo-lhe um peso menor. A única razão pela qual essas quadraturas são melhores para funções suficientemente chatas é que o manuseio adequado dos efeitos de borda supera o efeito das informações descartadas. De outro ponto de vista, é intuitivamente claro para mim que, para funções com uma estrutura ou ruído fino, as amostras que estão distantes das fronteiras do domínio de integração devem ser quase equidistantes e ter quase o mesmo peso (para um grande número de amostras). ) Por outro lado, a quadratura de tais funções pode se beneficiar de uma melhor manipulação dos efeitos de borda (do que no método do ponto médio).

Questão

Suponha que eu desejo integrar numericamente dados unidimensionais ruidosos ou bem estruturados.

O número de pontos de amostragem é fixo (devido à avaliação da função ser cara), mas posso colocá-los livremente. No entanto, eu (ou o método) não posso colocar pontos de amostragem interativamente, ou seja, com base nos resultados de outros pontos de amostragem. Também não conheço regiões com problemas em potencial de antemão. Então, algo como Gauß – Legendre (pontos de amostragem não equidistantes) está bem; A quadratura adaptativa não é uma vez que requer pontos de amostragem colocados interativamente.

  • Algum método além do método do ponto médio foi sugerido para esse caso?

  • Ou: Existe alguma prova de que o método do ponto médio é melhor nessas condições?

  • De maneira mais geral: existe algum trabalho existente sobre esse problema?

Apêndice A: Exemplo específico de uma função bem estruturada

Desejo estimar para: come. Uma função típica é assim:01f(t)dt

f(t)=i=1ksin(ωitφi)ωi,
log ω i[1,1000]φi[0,2π]logωi[1,1000]

senos sobrepostos

Eu escolhi essa função para as seguintes propriedades:

  • Pode ser integrado analiticamente para obter um resultado de controle.
  • Ele possui uma estrutura fina em um nível que torna impossível capturar tudo isso com o número de amostras que estou usando ( ).<102
  • Não é dominado por sua estrutura fina.

Apêndice B: Referência

Para completar, eis uma referência no Python:

import numpy as np
from numpy.random import uniform
from scipy.integrate import simps, trapz, romb, fixed_quad

begin = 0
end   = 1

def generate_f(k,low_freq,high_freq):
    ω = 2**uniform(np.log2(low_freq),np.log2(high_freq),k)
    φ = uniform(0,2*np.pi,k)
    g = lambda t,ω,φ: np.sin(ω*t-φ)/ω
    G = lambda t,ω,φ: np.cos(ω*t-φ)/ω**2
    f = lambda t: sum( g(t,ω[i],φ[i]) for i in range(k) )
    control = sum( G(begin,ω[i],φ[i])-G(end,ω[i],φ[i]) for i in range(k) )
    return control,f

def midpoint(f,n):
    midpoints = np.linspace(begin,end,2*n+1)[1::2]
    assert len(midpoints)==n
    return np.mean(f(midpoints))*(n-1)

def evaluate(n,control,f):
    """
    returns the relative errors when integrating f with n evaluations
    for several numerical integration methods.
    """
    times = np.linspace(begin,end,n)
    values = f(times)
    results = [
            midpoint(f,n),
            trapz(values),
            simps(values),
            romb (values),
            fixed_quad(f,begin,end,n=n)[0]*(n-1),
        ]

    return [
            abs((result/(n-1)-control)/control)
            for result in results
        ]

method_names = ["midpoint","trapezoid","Simpson","Romberg","Gauß–Legendre"]

def med(data):
    medians = np.median(np.vstack(data),axis=0)
    for median,name in zip(medians,method_names):
        print(f"{median:.3e}   {name}")

print("superimposed sines")
med(evaluate(33,*generate_f(10,1,1000)) for _ in range(100000))

print("superimposed low-frequency sines (control)")
med(evaluate(33,*generate_f(10,0.5,1.5)) for _ in range(100000))

(Utilizo aqui a mediana para reduzir a influência de valores discrepantes devido a funções que possuem apenas conteúdo de alta frequência. Para a média, os resultados são semelhantes.)

As medianas dos erros relativos de integração são:

superimposed sines
6.301e-04   midpoint
8.984e-04   trapezoid
1.158e-03   Simpson
1.537e-03   Romberg
1.862e-03   Gauß–Legendre

superimposed low-frequency sines (control)
2.790e-05   midpoint
5.933e-05   trapezoid
5.107e-09   Simpson
3.573e-16   Romberg
3.659e-16   Gauß–Legendre

Nota: Após dois meses e uma recompensa sem resultado, publiquei isso no MathOverflow .

Wrzlprmft
fonte
É esse o tipo de problema em que você realmente está interessado? Em 1D, você provavelmente pode obter bons resultados rapidamente com a maioria dos métodos.
David Ketcheson
"Eu tenho um número fixo de pontos de amostragem e posso colocá-los livremente. No entanto, não posso colocar pontos de amostragem interativamente, ou seja, com base nos resultados de outros pontos de amostragem." Essa restrição não está clara para mim. Posso colocar os nós onde um algoritmo adaptativo os colocaria, desde que eu seja realmente inteligente (em vez de realmente usar o algoritmo adaptável)? Se não tenho permissão para ser "realmente inteligente", que tipo de posicionamento de nó é realmente permitido?
David Ketcheson
@DavidKetcheson: Esse é o tipo de problema que você realmente está interessado? - Sim, estou realmente interessado em 1D. - Em 1D, você provavelmente pode obter bons resultados rapidamente com a maioria dos métodos. - Lembre-se de que a avaliação da função pode ser cara. - então que tipo de posicionamento de nó é realmente permitido? - Editei minha pergunta na esperança de torná-la mais clara.
precisa saber é o seguinte
Obrigado que ajuda. Para mim, a pergunta ainda parece vaga. Eu acho que há uma pergunta simples e mais precisa que seria mais responsável. Seria necessário definir um conjunto de funções (que podem depender do número permitido de nós de quadratura) e uma métrica. Em seguida, você pode perguntar se o método do ponto médio é ideal nessa métrica em relação àquele conjunto de funções (onde presumivelmente o mesmo conjunto de nós deve ser usado para quadratura todas as funções).
David Ketcheson
1
@ DavidKetcheson: Seria necessário definir um conjunto de funções (que podem depender do número permitido de nós de quadratura) e uma métrica. - Dado que não encontrei nada útil sobre este assunto até agora, não vejo razão para impor essas restrições. Em vez disso, com essas restrições, correria o risco de excluir algum trabalho existente (ou prova fácil) por condições ou suposições ligeiramente diferentes. Se houver alguma maneira de capturar o cenário descrito em definições e similares para as quais exista um trabalho de referência ou uma prova fácil, fico feliz com isso.
precisa saber é o seguinte

Respostas:

1

Primeiro de tudo, acho que você não entende o conceito de quadratura adaptativa. A quadratura adaptativa não implica em "colocar pontos de amostra interativamente". A idéia por trás da quadratura adaptativa é conceber um esquema que integre uma determinada função a um determinado erro absoluto ou relativo (estimado) com o mínimo possível de avaliações de funções.

Uma segunda observação: você escreve "O número de pontos de amostragem é fixo (devido à avaliação da função ser cara), mas eu posso colocá-los livremente". Penso que a ideia deve ser que o número de pontos de amostragem (ou avaliações de funções na terminologia em quadratura) seja o menor possível (isto é, não fixo).

Então, qual é a idéia por trás da quadratura adaptativa, conforme implementada no QUADPACK, por exemplo?

  1. O ingrediente básico é uma regra de quadratura "aninhada": é uma combinação de duas regras de quadratura em que uma tem uma ordem (ou precisão) mais alta que a outra. Por quê? Com base na diferença entre essas regras, o algoritmo pode estimar o erro de quadratura (é claro, o algoritmo usará o mais preciso como resultado de referência). Exemplos podem ser a regra trapezoidal com nós e nós. No caso do QUADPACK, as regras são regras de Gauss-Kronrod. Essas são regras de quadratura interpolatória que usam uma regra de quadratura de Gauss-Legendre de uma determinada ordem 2 n + 1 N N 2 N - 12n2n+1Ne uma extensão ideal desta regra. Isso significa que uma ordem de quadratura mais alta pode ser obtida reutilizando os nós de Gauss-Legendre (ou seja, as avaliações de funções caras) com pesos diferentes e adicionando vários nós extras. Em outras palavras, a regra de ordem original de Gauss-Legendre integrará todos os polinômios de grau exatamente enquanto a regra estendida de Gauss-Kronrod integrará exatamente algum polinômio de ordem superior. Uma regra clássica é o G7K15 (Gauss-Legendre de 7ª ordem com Gauss-Kronrod de 15ª ordem). A mágica é que os 7 nós do Gauss-Legendre são um subconjunto dos 15 nós do Gauss-Kronrod, portanto, com 15 avaliações de funções, eu tenho uma avaliação em quadratura juntamente com uma estimativa de erro!N2N1

  2. O próximo ingrediente é uma estratégia de "dividir e conquistar". Suponha que você solte esse G7K15 em seu integrando e observe um erro de quadratura que seja muito grande para o seu gosto. O QUADPACK subdividirá o intervalo original em dois subintervalos igualmente espaçados. E então reavaliará os dois subintegrais usando a regra básica, G7K15. Agora, o algoritmo possui uma estimativa de erro global (que deve ser menor que a primeira), mas também duas estimativas de erro local. Ele escolhe o intervalo com o maior erro e divide este em dois. Duas novas integrais são estimadas e o erro global é atualizado. E assim por diante até que o erro global esteja abaixo do destino solicitado ou que o número máximo de subdivisões seja ultrapassado.

Então, eu desafio você a atualizar seu código acima usando o scipy.quadmétodo Talvez no caso de um integrando com muita "estrutura fina", você precise aumentar o número máximo de subdivisões (a limitopção). Você também pode jogar com os parâmetros epsabse / ou epsrel.

No entanto, se você tiver apenas dados experimentais, vejo duas possibilidades.

  1. Se você tiver a oportunidade de selecionar os pontos de medição, ou seja, valores de , eu os selecionaria equidistanamente e de preferência como potência para que você possa aplicar uma regra trapezoidal aninhada (e lucrar com a extrapolação de Romberg).2t2
  2. Se você não tem como escolher os nós, ou seja, as medições ocorrem em momentos aleatórios, a melhor opção na minha opinião ainda é a regra do trapézio.
GertVdE
fonte
Eu acho que você entende mal o conceito de quadratura adaptativa. - Sua postagem concorda completamente com o meu entendimento anterior de quadratura adaptativa e é uma combinação clara de como eu defini interativamente a colocação dos pontos de amostragem (seja uma frase apropriada ou não). - você escreve [...]. Penso que a ideia deve ser que o número de pontos de amostragem […] seja o menor possível (isto é, não fixo). - Se você tem esse luxo, com certeza, mas restrições experimentais podem não ser tão benignas. Por exemplo, suponha que você precise medir algo simultaneamente com um número fixo de sensores caros.
Wrzlprmft 4/09/18
Me desculpe. Interpretei mal "interativamente" sua pergunta. No meu entendimento, "interativamente" significa intervenção do usuário e não de um algoritmo. Eu adicionei um parágrafo na minha resposta sobre dados experimentais. Outra abordagem seria "filtrar" as informações finas da estrutura, ou seja, aplicar uma transformada de Fourier e remover frequências de alta ordem com pequenas amplitudes. Isso seria uma opção?
GertVdE
Se você tiver a oportunidade de selecionar os pontos de medição [...] - Pontos equidistantes são o que eu preciso para o ponto médio, trapézio simples, etc. de qualquer maneira, então foi exatamente isso que fiz no meu benchmark. Aqui, a extrapolação de Romberg não oferece nenhuma vantagem.
Wrzlprmft
Outra abordagem seria "filtrar" as informações finas da estrutura [...] Isso seria uma opção? - No meu exemplo, suponho que a estrutura fina faça parte do que quero medir, apenas não tenho amostras suficientes para capturá-la completamente. Quanto ao ruído real, não há nenhuma restrição técnica que me impeça de filtrar. No entanto, a integral em todo o domínio já é o melhor filtro passa-baixo, por isso estou cético de que isso possa ser melhorado sem haver ruído com propriedades específicas, benignas e conhecidas.
precisa saber é o seguinte
É verdadeiramente estocástico? Deve haver alguns derivados que sejam aproximações integrais estocásticas de ordem superior.
Chris Rackauckas
0

Não estou convencido de que seu código demonstre algo fundamental sobre as várias regras de quadratura e quão bem elas são contra barulho e estrutura fina, e acredito que, se você escolher várias estruturas de multas diferentes, encontrará algo diferente. Aqui está o teorema:

Nenhum método de quadratura pode gerar um erro absoluto ou relativo baixo em relação a uma função com variação total ilimitada. Em um sistema de ponto flutuante com arredondamento de unidades , temos a estimativaμ onde é a soma da quadratura que atua na implementação numérica de .

|abfdxQ^[f^]||abfdxQ[f]|+μ[4ab|f|dx+ab|xf|dx]
Q^f^f

Prova: Deixe os nós da quadratura serem e os pesos da quadratura (não negativos) sejam e denote suas aproximações de ponto flutuante por e . Suponha que satisfaça onde onde é o arredondamento da unidade. Então {xi}i=0n1{wi}i=0n1w^ix^if^f^(x)=f(x)(1+2δ)|δ|μμ

Q^[f^]=i=0n1w^if^(x^i)=i=0n1wi(1+δiw)f(xi+δixxi)(1+2δif)(1+δi)i=0n1wi[f(xi)+δixxif(xi)](1+δiw+2δif+δi)i=0n1wif(xi)+i=0n1δixwixif(xi)+wif(xi)(δiw+2δif+δi)
para que Isso pressupõe que a soma seja calculada sem erros; multiplique por para descartar essa suposição.
|Q^[f^]Q[f]|μi=0n1wi(|xif(xi)|+4|f(xi)|)4μ|f|dx+μ|xf|dx
n

Mutatis mutandis, você também pode mostrar que o resultado se aplica à aritmética de ponto fixo.

user14717
fonte
Obrigado pela resposta. Estou com problemas para entender o cenário que você está considerando e como ele se relaciona com a minha pergunta. O que você quer dizer com variação total ilimitada no ponto flutuante? A menos que eu esteja muito enganado, todos os meus resultados computacionais (exceto o caso de controle com Romberg e Gauß-Legendre) estão longe de serem influenciados por imprecisões da implementação aritmética (ponto flutuante ou ponto fixo). O ruído que estou considerando também não é de natureza numérica, mas experimental.
Wrzlprmft
@Wrzlprmft: Ponto flutuante é o resultado que pude provar. Também posso provar isso em ponto fixo, o que indica que o resultado vale para dados experimentais. Acredito que seja verdade para qualquer fonte de erro nos nós de quadratura. Eu editei para esclarecer.
User14717 6/06/19
Para dados experimentais, o resultado é muito mais convincente, porque em geral os dados experimentais não são diferenciáveis ​​e, portanto, a variação total é infinita.
User14717 6/06/19
Sinto muito, mas ainda não o segui. Seu resultado parece ser sobre o erro cometido ao implementar numericamente a quadratura, não sobre o erro da quadratura em si. O problema que estou tendo é sobre o último e, em particular, não vejo razão para acreditar que ele não se manifestaria por . μ=0
Wrzlprmft
A idéia principal aqui vem do número da condição de avaliação da função. Suas avaliações estão mal condicionadas, pois são barulhentas.
User14717 6/06/19