Сравнительный анализ функционального «Ядерного» алгоритма и метода полигона частот Текст научной статьи по специальности «Математика»

АПВПМ-2019

СРАВНИТЕЛЬНЫЙ АНАЛИЗ ФУНКЦИОНАЛЬНОГО «ЯДЕРНОГО» АЛГОРИТМА И МЕТОДА ПОЛИГОНА ЧАСТОТ

Т. Е, Булгакова1, А. В, Войтишек1,2

1 Новосибирский государственный университет, 630090, Новосибирск 2Институт вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск

УДК 519.676

Б01: 10.24411/9999-016А-2019-10011

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

1. Многомерный аналог метода полигона частот

В последние годы (главным образом в новосибирской школе методов Монте-Карло) развивается теория рандомизированных функциональных алгоритмов (см., например, [1, 2]). Наиболее содержательные примеры применения этих алгоритмов связаны с приближением неизвестного решения у(х), х € М^ интегрального уравнения Фредгольма второго рода

|f(x)=J к(х', х)^(х') дх! + /(х) или ^ = + /, (1.1)

па ограниченной области X С М^; здесь к(х', х) (ядро интегрального оператора X) и f (х) (свободный член уравнения) — заданные функции.

Для приближения функции у(х) используем представления классической теории численной аппроксимации функций (см., например, [3]), имеющих общий вид

м

<р(х) « Ь(мV« = £^^(х) (1.2)

¿=1

для некоторого специально выбранного набора базисных функций

Н(м) = {х(1)(х),...,х(М)(х)} , (1.3)

Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект 0315-2019-0002), а также при финансовой поддержке Совета по грантам президента РФ и государственной поддержке ведущих научных школ (код проекта НШ-5913.2018.1).

ISBN 978-5-901548-42-4

(вид этих функций определяет тип аппроксимации (1.2)) и коэффициентов

■W(M) = |

= ы(1),

,,(м)

(1.4)

определяемых как функционалы от неизвестной приближаемой функции у(х).

В рандомизированных функциональных алгоритмах коэффициенты (1.4) вычисляются приближенно методом Монте-Карло -ш(1) « й(1\щ) с числом испытаний щ (в данной работе будет изучаться случай п! = ... = пм = п) и рассматривается приближение

м

ф(х) « Ь(м)ф(х) = £ *Ъ(')(п)Х(')(х).

г=!

(1.5)

В работах [4, 5] предложена новая (по сравнению с книгами [1, 2]) классификация рандомизированных функциональных алгоритмов для приближения решения уравнения (1.1): выделены сеточные, проекционные и проекционно-сеточные алгоритмы (тип метода определяется выбором базисных функций (1.3), коэффициентов (1.4) и приближения (1.5)). В этих же работах приведены соображения о том, почему сеточные и проекционные рандомизированные функциональные алгоритмы могут быть неэффективными (и даже нереализуемыми) при решении практически значимых задач, связанных с решением уравнений вида (1.1). Так, для теоретически привлекательного сеточного метода зависимых испытаний требуется гладкость ядра к(х', х) интегрального оператора К, а подавляющее большинство ядер из прикладных задач содержат интегрируемые особенности (вплоть до дельта-функций) и даже не могут быть явно вычислены. Сеточный метод сопряженных блужданий является крайне трудоемким из-за необходимости моделирования индивидуального набора траекторий соответствующих прикладных цепей Маркова для каждого узла х^ вводимой сетки

X(м) = {х!,..., хм } (1.6)

в области X. Проекционные методы обладают достаточно очевидной численной неустойчивостью.

Проекционно-сеточные рандомизированные функциональные алгоритмы не обладают перечисленными недостатками. Для этих алгоритмов базисные функции (1.3) и коэффициенты и)(г) = и)(г) (у>(м^ ; <р(м) = {^(х!), ...,(р(хм)} из (1.4) связаны с сеткой (1.6) таким образом, что обеспечивают как малость детерминированной компоненты погрешности

,(В(Х)) = ^ _ Ь(М)И|В(Х)

для используемого нормированного функционального пространства В(Х), а также устойчивость приближения (1-5), определяемую относительной малостью (близостью к единице) константы Лебега

лМ

Ь = 8ирхеХ |х(г)(х)| из соотношения

5

(€(Х)) зЬосИ

\\Ъ(М)у _ Ь(М)~||

€(Х) < Ъ .

тах

= !,...,М

(Vм^ _ п,(') (Vм) (п))

(см., например, [1, 2]); здесь ф(м)(п) = (^5(х1)(п),..., ф(хм)(п)), а ^$(х^(п) — монте-карловское приближение значения ^(х^); г = 1, ...,М.

В этом случае приближения метода Монте-Карло

= /,7,«

■(м) = { й(!) (п),...>(м )(п)}

для коэффициентов (1.4) из соотношения (1.5) имеют вид й(1\п) = IV

п)| , чаще всего (ф(м)(п)^ = ф(х\п).

В свою очередь, для получения значений ф(м)(п) в рандомизированных функциональных проекционно-сеточиых алгоритмах применяется следующая специальная технология (определяющая, в частности, отличие от сеточных функциональных алгоритмов). Выбираются финитные, одинаковые по форме для всех {х!;..., хм} функции (версии «ядерной» функции к(х)(у) для различных значений параметра х — см. раздел 2 данной статьи)

К(м) = {«(х1)(у),...,«(х-^у)} , (1.7)

так же, как и базисные функции (1.3), связанные с сеткой (1.6), и такие, что

У^Ы^У) ¿у «^(х); г = 1,..., М; (1.8)

см. [1, 2, 6, 7].

Далее следует напомнить ставшие уже классическими (см., например, главу 4 из учебника [1]) соображения о том, что для приближенного вычисления линейных функционалов вида

1н = У^(уЖу)^у (1.9)

от решения у>(х) уравнения (1.1) целесообразно использовать основной оцениватель (или несмещенную монте-карловскую оценку по столкновениям):

N

1н = ЕС; С=Е^(т)^ (^(т)) ' (1.10)

т=0

где

^(0), ^(1),..., ^) (1.11)

является прикладной цепью Маркова (или однородной цепью Маркова, обрывающейся с вероятностью единица) с начальной плотностью и переходной функцией р(х', х) = г(х', х) х [1 — р(а) (х')] (здесь г(х', х) — вероятностная переходная плотность, а 0 < р(а) (х') < 1 обозначает вероятность обрыва траектории; соответственно, N — это случайный номер обрыва траектории). Случайные веса |^(т)} из (1.10) определяются следующими рекуррентными соотношениями:

(0) () ( 1) * (^(т-1)' д(0) = ^—; д(т) = д(т-1) х^-f; т = 1,...,N. (1.12)

* (V0)) ^

Учитывая, что соотношения (1.8) имеют вид (1.9) получаем следующий рандомизированный проек-ционно-сеточный функциональный алгоритм. АЛГОРИТМ 1.1. Моделируя п траекторий

); з = 1,-,п (1.13)

прикладной цепи Маркова (1.15), получаем значения

1 п N

^)(п)=Пеефт));г=1,...,М;

j = 1 т=0

здесь веса | вычисляются по формулам вида (1.12);

/Й0)) ( ) ( 1) АЙт-1)'^т))

V ■> ) . п(т) — п(т-1) ^ V 7_

п(0) = У 3 ) . п(т) = П(т-1) х ^ 3_—_• 1 = 1 п- т = 1 N•

Затем приближаем функцию ^(х) по формуле вида (1.5):

м

<р(х) « Ь(м)ф(х) = (ё(х1)(п),..., <^(хм)(п^ Х(4)(х). (1.14)

В работах [1, 2, 6, 7] приведены соображения теории условной оптимизации, подтверждающие целесообразность использования в качестве базисных функций (1.3) «абсолютно устойчивых» финитных функций мультилинейной аппроксимации (или аппроксимации Стренга — Фикса [8] с производящей функцией Р(1)(и), являюгцейся В-сплайпом первого порядка) на регулярной сетке с шагом к по каждой координате:

((1) \ / \ ( и +1, при — 1 < и < 0;

^ — Л(1М х ... хР(1) ^ — ; Р(1)(и)= I —и + 1, при 0 < и < 1; (1.15)

) V ) I 0 иначе;

x = (x(1\..,x(d)) x = (j(i1)h,...,j<(d)h); j(k) -целые числа; i = 1,...,M (1.16)

(здесь множество X, на котором происходит приближение решения ^>(x) уравнения (1.1), представляет собой прямоугольный параллелепипед). Кроме того, «ядерную» функцию из соотношений (1.7), (1.8) предлагается выбирать в виде

K(x)(y) = i F ПРИ Х 6 Д(Х)' (1.17)

WJ I 0 иначе, v '

где Д(Х) = {y = (у(1),..., y(d)) : x(s) - h/2 < y(s) < x(s) + h/2; s = 1,...,d; x = (x(1),..., x(d))}. При этом приближения аппроксимационных коэффициентов (1.4) имеют простейший вид

w(i) (ф(х1)(п),..,ф(хм) (п)) = ф(х)(п). (1.18)

Алгоритм 1.1, в котором используются функции (1.15), (1.17) и аппроксимационные коэффициенты (1.18), назван в [1, 2, 6, 7] многомерным аналогом метода полигона частот.

В разделе 2 данной работы показано, что описанный подход к построению проекционно-сеточного алгоритма 1.1 в определенной степени аналогичен построению «ядерных» оценок вероятностных плотностей (см., например, [9]). Особо отмечено, что в работах по теории «ядерных» оценок (в том числе в [9]), авторы, рассуждая о приближениях функций, неоправданно не включают в рассмотрение элементы теории численной аппроксимации функций (см., например, [3]). При добавлении этого упущенного элемента соответствующие «ядерные» приближения вероятностных плотностей по сути совпадают с алгоритмом 1.1.

2. Численная аппроксимация вероятностных плотностей с использованием «ядерных» оценок

В классической работе [9] рассматривается непараметрическая оценка вероятностной плотности распределения Р(х), х € М^ по выборочным значениям ..., | С М^ из этого распределения вида

1 п

Р(х) « ^(х) = «(*) , (2.1)

где к(х) (у) — некоторая финитная параметрическая, одинаковая по форме для всех значений параметра х, «ядерная» функция. Приближение (2.1) называется «ядерной» оценкой плотности Р(х). К слову, в работе [9] и других источниках (в частности, в разделе 2 работы [10]) слова «ядро», «ядерная» приводятся без кавычек, но в данной работе мы используем кавычки для того, чтобы различить по названиям функции к(х', х) (ядро интегрального уравнения (1.1)) и к(х)(у).

При исследовании свойств приближения (2.1) существенно используется то обстоятельство, что, согласно закону больших чисел, для достаточно больших п выполнено

1 п ¡'

гп(х) = - £ п(х) « Еп(х) (¿) = к(х\у)Р(у) ¿у. (2.2)

п 3=1

Определенным конструктивным недостатком теории «ядерных» оценок плотности и ее приложений (см., в частности, работы [9, 10]) является отсутствие обсуждения алгоритма практического (в первую очередь, численного, компьютерного) приближения функции Р(х) в целом, основанного на теории сеточной аппроксимации функций (см., например, [а]). Такой алгоритм мог бы выглядеть следующим образом.

Учитывая конечность выборки |^1,..., можно рассмотреть ограниченное множество X С М^, которому принадлежат все выборочные значения, ввести в этой области сетку (1.6) и рассмотреть приближение вида (1.2) для функции Р(х):

м

Р(x) « L(M)Р(x) = £w(i) (Р(xi), ..,Р(xM)) X(i)(x).

i=i

АЛГОРИТМ 2.1. Вычисляем значения Р(х^(п) = Zn(xi); % = 1, ...,М по формулам вида (2.1) и приближаем функцию Р (х) по формуле вида (1.14);

м

Р(х) « Ь(м^ Р(х) = (Р(х1)(п), ..,Р(хм)(п)^ Х(0(х). (2.3)

¿=1

Алгоритм 2.1 основан на аналогах соотношений (1.12)

У РЫ^Ы ¿у « Р(х); г = 1,..., М, (2.4)

которые, в свою очередь, следуют из соотношений (2.1) и (2.2).

Сравнение алгоритмов 1.1 и 2.1 дает основной вывод данной работы.

ЗАМЕЧАНИЕ 2.1. Основанный на подходах численной сеточной аппроксимации функций «ядерный»

Р(х)

ванным, проекционно-сеточным функциональным алгоритмом 1.1 приближения решения у>(х) уравнения Фредгольма второго рода (1.1). Поэтому для алгоритма 1.1 вполне целесообразно применять повое название рандомизированный «ядерный» функциональный алгоритм для приближения решения уравнения (1.1).

Отличие алгоритмов 1.1 и 2.1 состоит только в разности форм монте-карловских оценок для приближенного вычисления функционалов вида (1.8) и (2.4) (что, в свою очередь, связано с определенным различием приближаемых функций у(х) и Р(х)). Отличие также состоит в том, что в задаче приближения плотности Р(х) выборка £п| считается заданной (при этом число выборочных значений п, как правило фиксировано и не может быть увеличено), а для функции у(х) число п моделируемых траекторий (1.13) прикладной цепи Маркова (1.11) может варьироваться.

Суть замечания 2.1 в определенной степени отражена в разделе 2 работы [10], но необходимые детали использования теории сеточной аппроксимации функций в этой работе не приводятся.

В связи с основным выводом замечания 2.1 можно сформулировать следующие соображения.

Р( х)

но использовать соображения теории условной оптимизации проекционно-сеточного функционального алгоритма 1.1 из работ [1, 2, 6, 7].

Так, по аналогии с работами [1, 2, 6, 7] можно рекомендовать использование функций (1.15), (1.17) и аппроксимационных коэффициентов вида (1.18), т.е.

(Р(х1)(п),..., Р(хм)(п)) = Р(х)(п), (2.5)

в алгоритме 2.1. Это позволяет, в частности, получить выражения для условно-оптимальных параметров алгоритма 2.1 вида

Морг =

Я1 [(2^ + 1^ + 4]

(2 lУ+1)d

а/2

7поРг = ^¡[(2: 2 х(21пМоРг-1п1пМоР,+Яз)х7-2-а/2

, , , (2.6) 7 констант Н1, Н2, #з и

Учитывая, что базис (1.15) является моделируемым, можно также предложить использовать нормированную функцию Ь(м)Р(х) го соотношения (2.3) в качестве плотности для численного моделирования

дополнительных выборочных значений близких по распределению к значениям из соотношения

(2.1), с использованием соответствующей версии метода суперпозиции (см. разделы 17, 18 книги [11]).

ЗАМЕЧАНИЕ 2.3. Для развития теории конструирования и условной оптимизации проекционно-сеточного функционального алгоритма 1.1 можно использовать соображения из теории «ядерных» оценок вероятностных плотностей из работы [9].

В последнем замечании речь идет прежде всего о варьировании и оптимальном выборе «ядерной» функции к(х)(у). В работе [9] рассматриваются возможности выбора этой функции в виде

«(х,(у> = П I2-7)

где положительные (вообще говоря, зависящие от величины выборки п) числа (п) > 0; в = 1, ^^^^^^^^^ ^^^^^ размытости») «ядерной» функции к(х)(у),

а ограниченные в совокупности четные (к(я)(у) = к(я)( —у)) функции к(я)(у) имеют единичный «второй момент» и конечные «то-е моменты»:

/+ТО /' + ТО

у2к(з)(у) ¿у = 1, / утк(а)(у) ¿у < то; т> 2. (2.8)

- ^ Л — то

Содержательные результаты по оптимизации приближения (2.1) с «ядерной» функцией (2.7) получаются лишь для случая

к(1)(п) = ... = к(л)(п) = Я(п); к(1)(у) = ... = к(л)(у) = к(у)

(см., в частности [9]); случай различных {к^ (п)}, |к(я)(у)}; в = 1,..., с1 практически не изучен. Для функции (1.17) имеем

<•«=к; ед={1|л|> 1/2: ^

т.е. «коэффициент размытости» совпадает с шагом равномерной сетки (1.16), а функция к(у) является кусочно-постоянной.

В работе [9] в рамках асимптотического (при п ^ то) подхода с помощью разложения функции р(х(1) + к(п)у(1),...,х(а) + к( п)у(л^) то всем переменным в точке х = в ряд Тейлора и ми-

нимизации относительной глобальной среднеквадратической ошибки

(к^2 = I Е (х) -^(х)]2 ^ ; д = I р2(х) .X получено асимптотическое приближение

2

¿х(1) ..¿х(л).

(¿(М2 „ ^-d(n)I d + ; , = f + кЧу) dy] J

^ ' У J— CO

"Л d2P(x(1),..,x(d)) Z^ О / (iw 2

i=i H-(l)Y

(2.10)

2r2

Минимизация полученного приближения величины и выбор функции к(у) из соображений ми-

нимизации величины I из (2.10) с соблюдением условий (2.8) дают оптимальное значение «коэффициента размытости» и оптимальную форму функции к(у):

к,,.,,) = () ^ , = { Л - ^ Ь 1 > V5: (2.п)

\П х Л у [ ^и 1у1 > V5.

к

(1.6) вида (1.16) в алгоритме 2.1 (в частности, получится ли выигрыш по сравнению со случаем, когда параметр к(п) и функция к(у) имеют вид (2.9)), требует отдельного подробного исследования.

Заключение

В работе показано, что при подключении теории численной аппроксимации функций можно получить «конструктивную» версию «ядерной» оценки для приближения вероятностной неизвестной плотности Р (x) по данным выборочным значениям (алгоритм 2.1). Эта конструкция аналогична рандомизированному проекци-онно-сеточному методу приближения решения ^>(x) интегрального уравнения Фредгольма второго рода (1.1) (т. е. алгоритму 1.1). Если дополнительно выбрать мультилинейный аппроксимационный базис (1.15), «ядерную» функцию (1-17) и аппроксимационные коэффициенты (2.5) (т.е. рассмотреть версию алгоритма 2.1, соответствующую многомерному аналогу метода полигона частот), то можно получить выражения (2.6) для условно-оптимальных параметров алгоритма 2.1.

Отмечено также, что можно назвать рандомизированный проекционно-сеточный алгоритм 1.1 функциональным «ядерным» алгоритмом и попытаться развить его теорию с помощью оптимизации выбора «ядерной» функции (2.7) подобно тому, как это сделано в работе [9].

Список литературы

[1] Михайлов Г. А., Войтишек А. В. Численное статистическое моделирование. Методы Монте-Карло. М.: Изд. центр «Академия», 2006.

[2] Войтишек А. В. Функциональные оценки метода Монте-Карло. Новосибирск: НГУ, 2007.

[3] Бахвалов Н. С. Численные методы. М.: Наука, 1975.

[4] Войтишек А. В. Классификация и возможности практического применения рандомизированных функциональных численных алгоритмов решения интегральных уравнений Фредгольма второго рода / / Математический анализ. Итоги науки и техники. Серия: Современная математика и ее приложения. Тематические обзоры. 2018. Т. 155. С. 3-19.

[5] Войтишек А. В. Разработка и оптимизация рандомизированных функциональных численных методов решения практически значимых интегральных уравнений Фредгольма второго рода // Сибирский журнал индустриальной математики. 2018. Т. 21. № 2 (74). С. 32-45.

[6] Shkarupa Е. V., Voytishek А. V. Convergence of discrete-stochastic numerical procedures with independent or weakly dependent estimators at grid nodes // Journal of Statistical Planning and Inference. 2000. V. 85. P. 199-211.

[7] Войтишек А. В., Головко H. Г., Шкарупа Е. В. Оценка погрешности многомерного аналога метода полигона частот // Сибирский журнал вычислительной математики. 2002. Т. 5, № 1. С. 11-24.

[8] Марчук Г. И., Агошков В. И. Введение в проекционно-сеточные методы. М.: Наука, 1981.

[9] Епанечников В. А. Непараметрическая оценка многомерной плотности вероятности // Теория вероятностей и ее применения. 1969. Т. 14. Вып. 1. С. 156-161.

[10] Михайлов Г. А. Рандомизированные алгоритмы метода Монте-Карло для задач со случайными параметрами (метод «двойной рандомизации») // Сибирский журнал вычислительной математики. 2019. Т. 22, № 2. С. 187-200.

[11] Войтишек А. В. Дополнительные сведения о численном моделировании случайных элементов. Новосибирск: НГУ, 2007.

Булгакова Татьяна Евгеньевна — старший преподаватель Специализированного учебно-научного центра Новосибирского государственного университета;

e-mail: tatyana.bulgakova@gmail.com, Войтишек Антон Вацлавович — д.ф.-м.н., профессор, ведущий научный сотрудник Института вычислительной математики и математической геофизики СО РАН;

профессор Новосибирского государственного университета;

e-mail: vav@osmf.sscc.ru.

Дата поступления — 30 апреля 2019 г.