Критерии оптимизации «Ядерного» алгоритма приближения вероятностной плотности Текст научной статьи по специальности «Математика»

КРИТЕРИИ ОПТИМИЗАЦИИ «ЯДЕРНОГО» АЛГОРИТМА ПРИБЛИЖЕНИЯ ВЕРОЯТНОСТНОЙ ПЛОТНОСТИ

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

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

630090, Новосибирск

УДК 519.676

Б01: 10.24411/9999-018А-2019-10003

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

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

1. Конструктивный «ядерный» алгоритм приближения вероятностной плотности

В классической работе [1] рассматривается непараметрическая оценка вероятностной плотности распределения Р|(х),х случайного вектора ^ по выборочным значе-

ниям |п} из этого распределения вида

Р? (Ж) « 7П (Ж) = 1 У" К(х)(?;) ,(1.1)

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

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

- У" К(х\$}) « £*«(?) = [ к (х\у)Р? (у) йу. (1.2)

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

Предположим, что область X распределения случайного вектора является ограниченной, и в этой области построена сетка

X (М) = [Х1.....хм } (1.3)

Работа выполнена в рамках государственного задания ИВМиМГ СО РАН (проект 0315-2019-0002).

16 "Проблемы оптимизации сложных систем - 2019"

и рассматривается сеточное приближение функции Р| (ж) вида

^ V,®[Р?(*!).....Р?(*м)] /®(ж).(1.4)

Здесь

£(м)= {/(1)М...../(м)(*)} (1.5)

является набором заданных базисных функций, как правило, связанных с сеткой (1.3) (при этом вид функций (1.5) определяет тип аппроксимации (1.4)), а величины

ж(м) = {ш(1) [Р? (*1).....Р? (*м )].....ш (М)[Р? (*1).....Р? (*м )]} (1.6)

представляют собой коэффициенты, являющиеся комбинациями значений приближаемой функции Р|(ж) в узлах сетки (1.3); чаще всего

^(0[Р?(*1).....Р?(*м)] = Р?(^). (1.7)

Для простоты изложения в дальнейшем будем полагать, что область X является прямоугольным параллелепипедом, а сетка (1.3) является равномерной прямоугольной с шагом Я:

= (/¿(1)Л,...,— целые числа; к = 1,..., / = 1,..., М.(1.8)

АЛГОРИТМ 1. Вычисляем значения Р~ (и) = ); * = 1, --, ^ по формулам вида

(1.1) и приближаем функцию Р|(ж) по формуле вида (1.4):

Р?(ж) « ¿(м)Р?(ж) = ш© [Р?(Ж1)(п), ..., /3;(жм)(п)]х(£)(^).(1.9) Алгоритм 1 основан на соотношениях

I к(ж')(у)Р|(у) ¿у « Р?(^),/ = 1.....М,(1.10)

¿х

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

2. Оптимизация по асимптотической верхней границе погрешности: выбор «ядерной» функции и «коэффициента размытости»

В работе [1] проводится следующая оптимизация приближения (1.1). Предлагается рассматривать «ядерную» функцию к(ж)(у) в виде

к «(у)=П11Л4) йй?); *=(х<1).....х(П у=(у (1).....у № ),(2-1)

где положительные (вообще говоря, зависящие от величины выборки п) числа Я(5)(п); 5 = 1,..., й определяют носитель («коэффициенты размытости») «ядерной» функции к(ж)(у), а ограниченные четные (к(5)(у) = к(5)(—у)) функции к(5)(у) имеют единичный «второй момент» и конечные «т-е моменты»:

I У2к(5)(у) ¿у = 1; I утк(5)(у) ¿у < Го, т >2. (2.2)

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

Я(1) (п) = - = Я№(п) = Я (п), к(1)(у) = - = к№(у) = к(у) (см., в частности, [1]); случай различных Я(5)(п),к(5)(у);5 = 1, ...,й практически не изучен.

В работе [1] в рамках асимптотического (при п ^ го) подхода с помощью разложения функции

Я(п)у(1), ...,х(й) + Я(п)у(й)) по всем переменным у = (у(1), ...,у(й)) в

точке ж

= (х(1).....х(й))

в ряд Тейлора и минимизации относительной глобальной средне-

квадратической ошибки

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

2

[^«)(п)]2 = ЕШМ-г.М] **; д = Г р|(ж) ^,(2.3)

V ¡х

получено асимптотическое приближение

Г,(, (у)) -.2 п-1 х Я-^(п) х Iй + (1/4) х Я4(п) х/ [^^(п)]--и-^ V I )-—--; (2.4)

2

г+т г Г \д2Рс(х(1),..., х(й))

' = ^ ау,1 = I . 1X1 )

йх(1) ...й). (2.5)

Минимизация полученного приближения (2.4) величины (2.3) и выбор функции к(у) из соображений минимизации величины / из (2.5) с соблюдением условий (2.2) дают оптимальное значение «коэффициента размытости» и оптимальную форму функции к(у):

. (а х /+4) - -^4= при \у\ < 75,

кор' (П) = {¿хт) ; *0^(у) = )475 2075 Р т (2.6)

V. о при \у\ > 75.

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

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

КРИТЕРИЙ 1. То приближение функции Р|(ж) лучше (оптимальнее), которое дает заданный уровень погрешности за меньшее время (с меньшими компьютерными затратами).

Вопрос о том, насколько использование параметров (2.6) в алгоритме 1 (в сочетании с выбором шага Я сетки (1.8)) соответствует критерию 1, требует отдельного подробного исследования.

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

3. Рандомизированный проекционно-сеточный функциональный алгоритм для решения интегрального уравнения Фредгольма второго рода

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

= Г ж) ф(х') йж' + /"(ж) или ^ = + / (3.1)

на ограниченной области X здесь ж) (ядро интегрального оператора К') и /(ж)

(свободный член уравнения) - заданные функции.

По аналогии с рассуждениями раздела 1 данной работы (см., в частности, соотношение

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

ш © /©(ж)

¿=1

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

Ж(м) = {ш(1).....ш(м)], (3.2)

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

В рандомизированных функциональных алгоритмах коэффициенты (3.2) вычисляются приближенно методом Монте-Карло с числами испытаний п^, т.е. № © « ш ®(п£) (в данной работе будет изучаться случай п1 = — = пм = п) и рассматривается приближение

я ©(п) /©(ж).

¿=1

В работах [8, 9] предложена новая (по сравнению с [3-7]) классификация рандомизированных функциональных алгоритмов для приближения решения уравнения (3.1): выделены сеточные, проекционные и проекционно-сеточные алгоритмы (тип метода определяется выбором базисных функций (1.5) и коэффициентов (3.2)). В этих же работах приведены соображения о том, почему сеточные и проекционные рандомизированные функциональные алгоритмы могут быть неэффективными (и даже нереализуемыми) при решении практически значимых задач, связанных с решением уравнений вида (3.1). Так, для теоретически привлекательного сеточного метода зависимых испытаний требуется гладкость ядра х) интегрального оператора а подавляющее большинство ядер из прикладных задач содержат интегрируемые особенности (вплоть до дельта-функций) и даже не могут быть явно вычислены. Сеточный метод сопряженных блужданий является крайне трудоемким из-за необходимости моделирования индивидуального набора траекторий соответствующих прикладных цепей Маркова для каждого узла ^ вводимой сетки (1.3) в области X. Проекционные методы обладают достаточно очевидной численной неустойчивостью.

Проекционно-сеточные рандомизированные функциональные алгоритмы не обладают перечисленными недостатками. Для этих алгоритмов (как и для сеточных методов) коэффициенты ц/ © = ц/ © (^(м)) из (3.2) представляют собой комбинации значений

^(м) = {ф(*1).....)} (3.3)

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

= {к(^1)(у),..., к(^м)(у)| (3.4)

(по сути это версии «ядерной» функции к(ж)(у) из (1.1) для различных значений параметра ж) так же, как и базисные функции (1.5), связанные с сеткой (1.3), и такие, что

[ к^СУЖУ) Ф ~ ); I = 1.....М (3.5)

(это аналоги соотношений (1.10)). Далее используется то обстоятельство, что приближения

(3.5) значений (3.3) представляют собой линейные функционалы от решения ф(лО уравнения (3.1). Для таких функционалов строятся основные оцениватели (или монте-карловские оценки по столкновениям - см., например, главу 4 учебника [5]), основанные на моделировании траекторий

.....^; ) = 1.....п (3.6)

прикладной цепи Маркова

^(0), ^(1).....^ (3.7)

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

АЛГОРИТМ 2 [8, 9]. Моделируя п траекторий (3.6) прикладной цепи Маркова (3.7), получаем значения

£(^)(П) =1 уп кад^).

здесь случайные веса {@уШ)} вычисляются по рекуррентным формулам

~(т) = 0(т) = ^(т-1) УУ ^ У /

^ = , ^ = к^г1', о'

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

ш^(п), ...,ф^(п)] /«(*).

¿=1

Сравнение алгоритмов 1 и 2 приводит к следующему важному выводу.

ЗАМЕЧАНИЕ 1. Основанный на подходах численной сеточной аппроксимации функций «ядерный» алгоритм 1 приближения вероятностной плотности Р|(ж) конструктивно совпадает с рандомизированным проекционно-сеточным функциональным алгоритмом 2 приближения решения уравнения Фредгольма второго рода (3.1).

Отличие алгоритмов 1 и 2 состоит только в разности форм монте-карловских оценок для приближенного вычисления функционалов вида (1.10) и (3.5) (что, в свою очередь, связано с определенным различием приближаемых функций Р| (ж) и ф(я)). Отличие также состоит в том, что в задаче приближения плотности Р|(ж) выборка ..., считается заданной (при этом число выборочных значений п, как правило, фиксировано и не может быть увеличено), а для функции число п моделируемых траекторий (3.6) прикладной цепи Маркова (3.7) может варьироваться.

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

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

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

Более существенными являются следующие соображения, связанные с замечанием 2.

Алгоритм 2 подробно изучен в работах [4, 7] лишь для частного случая, когда в качестве базисных функций (1.5) используются «абсолютно устойчивые» финитные функции муль-тилинейной аппроксимации (или аппроксимации Стренга - Фикса [12] с производящей функцией ^(1)(и), являющейся сплайном первого порядка) на регулярной сетке (1.8):

*®(*) = £(1) -Л(1)) X. .X £(1) (^ — у®) ; £(1)(и)

!и + 1 при — 1 < и < 0, -и + 1 при 0 < и < 1, (3.8) 0 при |и| > 1.

При этом «ядерная» функция из соотношений (3.4), (3.5) выбирается в виде

КМ(У) = [р при У^ (3.9)

( 0 при у £ Д(ж),

где Д(ж) = {у = (у(1),...,у(й)):х(5) — £<у(5) <х(5) + £;5 = 1,.,5; ж = (х(1),...,х(й))}, т.е в этом случае речь идет о представлении (2.1) при

Я(п) = Я и к(у) = 11 при < 1/2, (3.10) 4 ' УУ) ( 0 при |у| > 1/2. ^ )

При этом приближения аппроксимационных коэффициентов (3.2) имеют простейший

вид

ш ©[ф (*1)(п),., ф ^(п)] = ^ (ж')(п). (3.11)

Алгоритм 2, в котором используются функции (3.8), (3.9) и аппроксимационные коэффициенты (3.11), назван в [4-7] многомерным аналогом метода полигона частот.

В связи с рассуждениями раздела 2 данной статьи интересным представляется отдельное исследование алгоритма 2, включающее соображения теории условной оптимизации (см. [4-7], а также раздел 4 данной статьи), в случае, когда в соотношениях вида (3.5) используется «ядерная» функция (2.1) с «коэффициентом размытости» и функцией к(у) вида (2.6) (вместо (3.10)).

4. Условная оптимизация «ядерного» алгоритма

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

ЗАМЕЧАНИЕ 3. Для «ядерного» алгоритма 1 приближения вероятностной плотности Р|(ж) можно использовать соображения теории условной оптимизации проекционно-се-точного функционального алгоритма 2.

Общая схема условной оптимизации выглядит следующим образом (см., например, [47]). Ставится задача согласованного выбора параметров М и п (число узлов сетки (1.3) и выборочных значений из соотношения (1.1 )) для исследуемого функционального алгоритма (например, для алгоритма 1), обеспечивающего заданный уровень у >0 погрешности приближения исследуемой функции (например, приближения (1.9) из алгоритма 1) при минимальных вычислительных затратах 5(М, п) (таким образом, предлагаемый подход вполне соответствует «объективному» критерию 1, сформулированному в разделе 2 данной статьи).

Строится верхняя граница УР(1)(М, п) погрешности алгоритма 5 (1)(М, п) для используемого нормированного функционального пространства зависящая от параметров М и п:

5(1)(М,п) = ||Р| — ¿(м)Р||| < УР(1)(М,п). (4.1)

Эта функция двух переменных приравнивается величине у. Из уравнения вида

УР(1)(М, п) = у (4.2)

один из параметров (например, п) выражается через другой: п = ^(М). Это соотношение подставляется в выражение для затрат 5(М, п), которое тоже зависит от параметров М и п.

В результате получается функция 5(М) одного переменного М, которая исследуется на минимум с помощью известных приемов математического или численного анализа. Найденные значения М^^(у) = Мд^у), п® = ^Цм® (у)| объявляются условно-оптимальными

параметрами исследуемого алгоритма.

«Условность» такого способа оптимизации связана с тем, что в левой части уравнения вида (4.2) используется не сама погрешность алгоритма 5 (1)(М, п), а ее верхняя граница УР(1)(М, п). Впрочем, оценка качества того или иного алгоритма по верхней границе погрешности используется - это общее место в вычислительной математике (см., например, [2], а также раздел 2 данной работы).

При изучении погрешности 8 (1)(М, п) необходимо выбрать как соответствующее нормированное функциональное пространство так и вероятностный смысл выполнения неравенства вида (4.1) (ведь 6 (1)(М, п) является случайной величиной). Следуя канонам классического численного анализа (см., например, [2]), в качестве пространств ЕЗ(Х) будем рассматривать 1.2(Х) и С(Х).

Для хорошо разработанного (см., например, [3-7]) И,2-подхода рассматривается сходимость в среднем погрешности

5(^(М,п) = ||Р? — ¿^Ц^= {Г [Р?(*) — ¿(М)р?(^)]2) к нулю при М, п ^ от и строятся оценки сверху УР^2)(М, п) такие, что

[Е<5(^2)п)]2 < иР(^2)п). Для С-подхода [4-7] величина

5(С) (М, п) = ||Р? — ¿(М)^Нст = зир|Р? (ж) — ¿(М)Р? (*)|

ограничивается сверху по вероятности:

Р[5(С)(М, п) < УР(С)(М, п)] >1 — £ для некоторого достаточно малого £ >0.

Для этих двух подходов, по аналогии с рассуждениями из работ [4, 6, 7] для алгоритма 2, можно разделить погрешности приближения плотности Р| (ж) на три части: детерминированную компоненту 5®(М), стохастическую компоненту ^^^(М, и) и компоненту смещения (М).

Можно также показать, что в случае выполнения четырех условий:

1) используется регулярная сетка (1.8) с шагом Я,

2) используется численно устойчивый базис (1.5), обеспечивающий выполнение неравенства

= Н^? — ¿(М)< ^ь (*) — Ш1 (4.3)

для константы Лебега Н^, ненамного превосходящей единицу (или равной единице); здесь

¿=1

3) используются коэффициенты (1.6) простейшего вида (1.7),

¥) используются «коэффициент размытости» и функция к (у) вида (3.10) (т.е. выбирается «ядерная» функция (3.9))

Верхние границы для компонент смещения ¿"^(М) и ^^(М) имеют второй порядок по шагу сетки Я (этот шаг пропорционален величине М-1/й). Значит, второй порядок должны иметь и верхние границы для детерминированных компонент погрешности

^¿е* и . Такие верхние границы для компоненты смещения и детерминирован-

ной компоненты дают базисные функции (3.8) и аппроксимационные коэффициенты (1.7) (здесь, в частности, = 1 [4, 6]).

Заметим, что нарушение каждого из сформулированных четырех условий (в том числе, использование «коэффициента размытости» и функции к (у) вида (2.6) вместо (3.10)) может существенно усложнить получение верхних границ для компонент смещения (в частности,

для величины из соотношения (4.3)).

Свойства устойчивости вида (4.3) для базиса (3.8) и использование аппроксимационных коэффициентов (1.7) дают возможность построить оценки сверху для стохастических компонент погрешности (М, п) и (М, п).

Затем, по аналогии с рассуждениями из работ [4, 7], для 1,2- и С-подходов можно реализовать описанную выше методику теории условной оптимизации и получить следующие выражения для условно-оптимальных параметров:

гй +4\й/4 (п_),. ) М + 4^

^орс

£Лг) = (^Р) у-*/2; п^(у) = Я,^ (^р) (й + 4)у-2-*/2; (4.4) <1 (У) = ЯС

(2У + 1)Л + 4^/2 -й/2 „(С) [(2У + 1)Л + 4]2+^/2

(2у + 1)^ ' ^ ^^ = Я2 ) [(2У + 1)^ Х (4.5)

X (21п М(Ср)(.(у) — 1п 1п М® (у) + ЯзСс))у-2-*/2

для некоторых положительного констант Н^2^ я2^'2), Я®, я2С), Я® и V (выбор или приближение этих констант представляет собой отдельную - часто непростую - задачу [4]).

Заключение

В данной работе предложена модификация «ядерной» оценки (1.1) вероятностной плотности Р| (х), основанная на подходах теории численной аппроксимации функций и приводящая к конструктивному алгоритму 1. Проведен критический анализ критерия оптимизации «ядерных» оценок плотностей, основанного на минимизации верхней границы относительной погрешности для 1,2 -подхода. Показано, что построенный «ядерный» алгоритм 1 эквивалентен рандомизированному проекционно-сеточному функциональному алгоритму 2 для приближения решения интегрального уравнения Фредгольма второго рода (3.1). В связи с этим предложено использовать критерий условной оптимизации функциональных алгоритмов, основанный на минимизации затрат 5(М, п) при фиксированном уровне погрешности у, для «ядерного» алгоритма 1 и получен вид условно-оптимальных параметров (для известных 1,2- и С-подходов - см. формулы (4.4), (4.5)) для простейшего случая алгоритма 1, в котором используются базисные функции (3.8), аппроксимационные коэффициенты (1.7) и «ядерная» функция (3.9) (т.е. для аналога алгоритма метода полигона частот).

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

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

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

3. Михайлов Г.А. Весовые методы Монте-Карло. - Новосибирск: Изд-во СО РАН, 2000.

4. Войтишек А.В. Дискретно-стохастические численные методы: Дис ... д-ра физ.-мат. Наук / РАН. Сиб. отд-ние, ИВМиМГ. - Новосибирск, 2001.

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

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

7. Войтишек А.В. Основы метода Монте-Карло в алгоритмах и задачах. Часть VI. Вычисление значений линейных функционалов от решения интегрального уравнения второго рода. Дискретно-стохастические методы решения интегрального уравнения второго рода. - Новосибирск: НГУ, 2004.

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

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

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

11. Михайлов Г.А. Улучшение многомерных рандомизированных алгоритмов метода Монте-Карло с «расщеплением» // Журнал вычислительной математики и математической физики. - 2019. - Т. 59, № 5. - С. 822-828.

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

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

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

email: vav@osmf.sscc.ru