Для шумных или тонко структурированных данных существуют ли лучшие квадратуры, чем правило средней точки?

12

Только первые два раздела этого длинного вопроса важны. Остальные только для иллюстрации.

Фон

Продвинутые квадратуры, такие как составные выражения Ньютона – Котса, Гауза – Лежандра и Ромберга более высокой степени, по-видимому, в основном предназначены для случаев, когда можно тонко выбрать функцию, но не интегрировать аналитически. Однако для функций со структурами, меньшими, чем интервал дискретизации (пример см. В Приложении A) или шум измерения, они не могут конкурировать с простыми подходами, такими как правило средней точки или трапеции (для демонстрации см. Приложение B).

Это несколько интуитивно понятно, поскольку, например, составное правило Симпсона по существу «отбрасывает» одну четверть информации, присваивая ей меньший вес. Единственная причина, по которой такие квадратуры лучше подходят для достаточно скучных функций, заключается в том, что правильная обработка граничных эффектов перевешивает эффект отбрасываемой информации. С другой точки зрения, мне интуитивно понятно, что для функций с тонкой структурой или шумом выборки, удаленные от границ области интегрирования, должны быть почти равноудаленными и иметь практически одинаковый вес (для большого числа выборок ). С другой стороны, квадратура таких функций может выиграть от лучшей обработки граничных эффектов (чем для метода средней точки).

Вопрос

Предположим, что я хочу численно интегрировать шумные или тонко структурированные одномерные данные.

Количество точек выборки является фиксированным (из-за того, что оценка функции обходится дорого), но я могу свободно разместить их. Однако я (или метод) не могу размещать точки выборки в интерактивном режиме, т. Е. На основе результатов из других точек выборки. Я также не знаю потенциальных проблемных регионов заранее. Итак, что-то вроде Гаусса-Лежандра (неэквидистантные точки отбора проб) в порядке; адаптивной квадратуры нет, так как она требует интерактивных точек выборки.

  • Были ли предложены какие-либо методы, выходящие за пределы метода средней точки, для такого случая?

  • Или: есть ли доказательства того, что метод средней точки лучше в таких условиях?

  • В целом: есть ли какие-либо работы по этой проблеме?

Приложение A: Конкретный пример тонко структурированной функции

Я хотел бы оценить для: си. Типичная функция выглядит так:01f(t)dt

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

наложенные синусы

Я выбрал эту функцию для следующих свойств:

  • Он может быть интегрирован аналитически для контроля результата.
  • Он имеет тонкую структуру на уровне, который делает невозможным захват всего этого с количеством образцов, которые я использую ( ).<102
  • В нем не доминирует его тонкая структура.

Приложение B: контрольный показатель

Для полноты, вот пример в 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))

(Я здесь использую медиану, чтобы уменьшить влияние выбросов из-за функций, которые имеют только высокочастотное содержание. В среднем, результаты похожи.)

Медианы относительных ошибок интеграции:

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

Примечание: после двух месяцев и одного вознаграждения без результата я разместил это на MathOverflow .

Wrzlprmft
источник
Это та проблема, которая вас действительно интересует? В 1D вы, вероятно, можете получить хорошие результаты довольно быстро с большинством любых методов.
Дэвид Кетчесон
«У меня есть фиксированное количество точек выборки, и я могу свободно размещать их. Однако я не могу размещать точки выборки в интерактивном режиме, т. Е. На основе результатов из других точек выборки». Это ограничение мне не понятно. Могу ли я размещать узлы там, где их будет адаптивный алгоритм, если я просто очень умный (вместо того, чтобы использовать адаптивный алгоритм)? Если мне не позволено быть «по-настоящему умным» в этом, то какие виды размещения узлов на самом деле разрешены?
Дэвид Кетчесон
@DavidKetcheson: Это та проблема, которая вас действительно интересует? - Да, я действительно заинтересован в 1D. - В 1D вы, вероятно, можете получить хорошие результаты довольно быстро с большинством любых методов. - Помните, что оценка функции может быть дорогостоящей. - тогда какие виды размещения узлов действительно разрешены? - Я отредактировал свой вопрос в надежде прояснить его.
Wrzlprmft
Спасибо, что помогает. Мне вопрос все еще кажется расплывчатым. Я думаю, что есть простой и более точный вопрос, который был бы более ответственным. Это потребует определения набора функций (которые могут зависеть от разрешенного количества квадратурных узлов) и метрики. Затем вы можете спросить, является ли метод средней точки оптимальным в этой метрике по сравнению с этим набором функций (где, по-видимому, один и тот же набор узлов должен использоваться для квадратуры всех функций).
Дэвид Кетчесон
1
@DavidKetcheson: для этого потребуется определить набор функций (которые могут зависеть от разрешенного количества квадратурных узлов) и метрику. - Учитывая, что мне пока не удалось найти ничего полезного по этому вопросу, я не вижу оснований для введения таких ограничений. Скорее, с такими ограничениями, я бы рискнул исключить некоторые существующие работы (или простые доказательства) для немного других условий или предположений. Если есть какие-либо способы запечатлеть изображенный сценарий в определениях и аналогичных материалах, для которых существует справочная работа или простое доказательство, я рад этому.
Wrzlprmft

Ответы:

1

Прежде всего, я думаю, что вы неправильно понимаете концепцию адаптивной квадратуры. Адаптивная квадратура не подразумевает «интерактивное размещение точек выборки». Вся идея адаптивной квадратуры состоит в том, чтобы разработать схему, которая будет интегрировать определенную функцию с определенной (оценочной) абсолютной или относительной погрешностью с как можно меньшими оценками функции.

Второе замечание: вы пишете: «Количество точек выборки фиксировано (из-за того, что оценка функции обходится дорого), но я могу свободно разместить их». Я думаю, что идея должна состоять в том, чтобы число точек отбора проб (или оценки функций в квадратурной терминологии) было как можно меньше (т.е. не фиксировано).

Так в чем же идея адаптивной квадратуры, реализованной, например, в QUADPACK ?

  1. Основным ингредиентом является «вложенное» квадратурное правило: это сочетание двух квадратурных правил, где одно имеет более высокий порядок (или точность), чем другое. Почему? Основываясь на разнице между этими правилами, алгоритм может оценить квадратурную ошибку (конечно, алгоритм будет использовать наиболее точную в качестве контрольного результата). Примерами могут служить правило трапеции с узлами и узлами. В случае с QUADPACK правила являются правилами Гаусса-Кронрода. Это интерполяционные квадратурные правила, которые используют квадратурное правило Гаусса-Лежандра определенного порядка 2 n + 1 N N 2 N - 12n2n+1Nи оптимальное расширение этого правила. Это означает, что более высокий квадратурный порядок может быть получен повторным использованием узлов Гаусса-Лежандра (т. Е. Оценки дорогостоящих функций) с различными весами и добавлением ряда дополнительных узлов. Другими словами, оригинальное правило Гаусса-Лежандра порядка будет точно интегрировать все многочлены степени , тогда как расширенное правило Гаусса-Кронрода точно интегрирует некоторый многочлен более высокого порядка. Классическим правилом является G7K15 (Гаусс-Лежандр 7-го порядка с Гаусс-Кронродом 15-го порядка). Волшебство заключается в том, что 7 узлов Гаусса-Лежандра являются подмножеством 15 узлов Гаусса-Кронрода, поэтому с 15 оценками функций у меня есть квадратурная оценка вместе с оценкой ошибки!N2N1

  2. Следующий компонент - стратегия «разделяй и властвуй». Предположим, что вы выпустили этот G7K15 на свою интегральную функцию и наблюдаете квадратурную ошибку, которая по вашему вкусу слишком велика. Затем QUADPACK подразделяет исходный интервал на два равных интервала. И затем он пересмотрит эти два субинтеграла, используя основное правило, G7K15. Теперь алгоритм имеет глобальную оценку ошибки (которая, как мы надеемся, должна быть ниже первой), а также две локальные оценки ошибки. Он выбирает интервал с наибольшей ошибкой и делит его на два. Оцениваются два новых интеграла и обновляется глобальная ошибка. И так до тех пор, пока глобальная ошибка не окажется ниже запрошенной цели или не будет превышено максимальное количество подразделений.

Поэтому я призываю вас обновить код выше, используя scipy.quadметод. Возможно, в случае подынтегральной функции с большой «тонкой структурой» вам может потребоваться увеличить максимальное количество подразделений ( limitопция). Вы также можете поиграть с параметрами epsabsи / или epsrel.

Однако, если у вас есть только экспериментальные данные, я вижу две возможности.

  1. Если у вас есть возможность выбрать точки измерения, то есть значения , я бы выбрал их на равном расстоянии, предпочтительно в виде степени чтобы вы могли применить вложенное правило трапеции (и получить прибыль от экстраполяции Ромберга).2t2
  2. Если у вас нет средств для выбора узлов, то есть измерения производятся в случайное время, лучшим вариантом, на мой взгляд, по-прежнему является правило трапеции.
GertVdE
источник
Я думаю, что вы неправильно понимаете концепцию адаптивной квадратуры. - Ваше сообщение полностью согласуется с моим предыдущим пониманием адаптивной квадратуры, и оно полностью соответствует тому, как я определил в интерактивном режиме размещение точек выборки (является ли это подходящей фразой или нет). - ты пишешь […]. Я думаю, что идея должна состоять в том, чтобы число точек отбора проб […] было как можно меньше (т.е. не фиксировано). - Если у вас есть такая роскошь, конечно, но экспериментальные ограничения могут быть не такими уж мягкими. Например, предположим, что вам нужно что-то измерять одновременно с фиксированным количеством дорогих датчиков.
Wrzlprmft
Мои извенения. Я неверно истолковал «интерактивно» в вашем вопросе. В моем понимании «интерактивно» означает вмешательство пользователя не алгоритмом. В свой ответ я добавил параграф об экспериментальных данных. Другой подход заключается в том, чтобы «отфильтровать» информацию тонкой структуры, то есть применить преобразование Фурье и удалить частоты высокого порядка с малыми амплитудами. Будет ли это вариант?
GertVdE
Если у вас есть возможность выбрать точки измерения […] - Эквидистантные точки - это то, что мне нужно для средней точки, простой трапеции и т. Д., Так что это именно то, что я сделал в своем тесте. Здесь экстраполяция Ромберга не дает никаких преимуществ.
Wrzlprmft
Другой подход заключается в том, чтобы «отфильтровать» информацию о тонкой структуре […] Это будет вариант? - В моем примере я предполагаю, что тонкая структура является частью того, что я хочу измерить, но у меня просто не хватает достаточного количества образцов, чтобы захватить ее полностью. Что касается фактического шума, нет никаких технических ограничений, которые мешают мне фильтровать. Однако интеграл по всей области уже является конечным фильтром нижних частот, поэтому я скептически отношусь к тому, что его можно улучшить, не имея шума с особыми, доброкачественными и известными свойствами.
Wrzlprmft
Это действительно стохастик? Должны быть некоторые производные, которые являются стохастическими интегральными приближениями высшего порядка.
Крис Ракаукас
0

Я не уверен, что ваш код демонстрирует что-то фундаментальное о различных квадратурных правилах и о том, насколько хорошо они справляются с шумом и тонкой структурой, и полагаю, что если вы выберете различную структуру штрафов, вы найдете что-то другое. Вот теорема:

Ни один квадратурный метод не может дать низкую абсолютную или относительную погрешность для функции с неограниченным полным изменением. В системе с плавающей запятой с единичным округлением мы имеем оценкуμ где - квадратурная сумма, влияющая на численную реализацию функции .

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

Доказательство. Пусть квадратурные узлы будут а (неотрицательные) квадратурные веса - и обозначим их приближения с плавающей запятой через и . Предположим, что удовлетворяет где где - округление единицы. потом {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)
так, чтобы Это предполагает, что сумма вычисляется без ошибок; умножить на чтобы отбросить это предположение.
|Q^[f^]Q[f]|μi=0n1wi(|xif(xi)|+4|f(xi)|)4μ|f|dx+μ|xf|dx
n

Mutatis mutandis вы также можете показать, что результат выполняется в арифметике с фиксированной запятой.

user14717
источник
Спасибо за ответ. У меня возникли проблемы с пониманием сценария, который вы рассматриваете, и как он связан с моим вопросом. Что вы имеете в виду под неограниченным полным изменением числа с плавающей запятой? Если я не очень ошибаюсь, все мои вычислительные результаты (за исключением контрольного случая с Ромбергом и Гауссом-Лежандром) далеки от влияния неточностей арифметической реализации (с плавающей точкой или с фиксированной точкой). Шум, который я рассматриваю, тоже не численный, а экспериментальный.
Wrzlprmft
@Wrzlprmft: Плавающая точка - это результат, который я смог доказать. Я также могу доказать это в фиксированной точке, которая затем указывает, что результат верен для экспериментальных данных. Я считаю, что это верно для любого источника ошибки в квадратурных узлах. Я отредактировал, чтобы уточнить.
user14717
Для экспериментальных данных результат гораздо более убедителен, потому что в целом экспериментальные данные недифференцируемы и, следовательно, общее изменение бесконечно.
user14717
Извините, но я все еще не могу следовать за вами. Похоже, ваш результат связан с ошибкой, допущенной при численной реализации квадратуры, а не с ошибкой самой квадратуры. Проблема, с которой я сталкиваюсь, связана с последним, и, в частности, я не вижу причин полагать, что она не проявится при . μ=0
Wrzlprmft
Основная идея здесь исходит из числа условий оценки функции. Ваши оценки плохо обусловлены, так как они шумные.
user14717