Журнал вычислительной математики и математической физики, 2019, T. 59, № 5, стр. 822-828

ВВЕДЕНИЕ

Численные методы Монте-Карло строятся на основе естественных или искусственно сформулированных вероятностных моделей, которые численно статистически реализуются с помощью известных алгоритмов (см., например, [1]–[3]). Они могут быть сравнительно эффективными, когда базовые модели содержат случайные параметры и для оценки искомых величин используется совместное статистическое моделирование базовых и параметрических распределений, то есть фактически реализуется произведение соответствующих вероятностных пространств (возможно, многократное). Соответствующие алгоритмы метода Монте-Карло в настоящей работе называются рандомизированными. Формулируемый таким образом метод ”двойной рандомизации” можно пояснить, рассматривая интеграл

со случайным, возможно функциональным, параметром $\sigma $. Здесь $P(dw;\sigma )$ – вероятностная мера в $W$ с параметром $\sigma $. Пусть необходимо оценить математическое ожидание $J = EJ(\sigma )$. Если определена достаточно точная оценка $J(\sigma ) \approx \hat {J}(\sigma )$, то численное построение выборки $\left\{ {{{\sigma }_{i}}} \right\}$ дает статистические оценки требуемых величин. Однако для реальных задач такой алгоритм может быть слишком трудоемким. При этом целесообразно использовать двойную рандомизацию, моделируя для выбранного $\sigma $ лишь сравнительно небольшое число $n$ точек $\omega $ соответственно распределению $P(dw;\sigma )$ с вычислением и осреднением полученных значений $g(\omega ;\sigma )$. Оптимальное по критерию трудоемкости вычислений значение $n$ здесь оценивается по формулам стохастического метода “расщепления” [4] (см. далее п. 1).

В настоящей работе решается более сложная задача оптимизации рандомизированного алгоритма для случая, когда $J(\sigma ): = J(x,\sigma )$ и необходимо оценить функцию $f(x) = EJ(x,\sigma )$.

1. ЧИСЛЕННО-СТАТИСТИЧЕСКАЯ ФУНКЦИОНАЛЬНАЯ ОЦЕНКА, ЕЕ ТРУДОЕМКОСТЬ

1. Для построения функциональной оценки рассматриваются линейные функционалы вида

Здесь $x = ({{x}_{1}},\; \ldots ,\;{{x}_{m}}) \in {{R}^{m}}$; $\sigma $ – случайный, возможно функциональный, параметр задачи (“случайная среда”); $f(x;\sigma ) \in {{L}_{1}}({{R}^{m}})$ – решение задачи с параметром $\sigma $, определяемое компьютерно реализуемой вероятностной моделью, т.е. базовым ансамблем траекторий $\left\{ \Omega \right\}$ в фазовом пространстве ${{R}^{m}}$; ${{h}_{k}}(x;\sigma ) \in {{L}_{\infty }}({{R}^{m}})$.

Методом Монте-Карло строятся несмещенные оценки ${{\xi }_{k}}(\Omega ;\sigma )$ функционалов ${{J}_{k}}(\sigma )$, то есть при фиксированном $\sigma $ имеем: ${{E}_{\Omega }}{{\xi }_{k}}(\Omega ;\sigma ) = {{J}_{k}}(\sigma )$. Для иллюстрации такой схемы можно рассматривать задачу переноса частиц – квантов излучения – с рассеянием и поглощением через среду со случайной плотностью $\sigma (r)$, $r \in {{R}^{3}}$ (см., например, [5], [6]); здесь $\left\{ \Omega \right\}$ – ансамбль траекторий квантов, который можно определить однородной цепью Маркова столкновений квантов с элементами вещества [7]. Отметим, что методом Монте-Карло при этом осуществляется осреднение функционалов от решения интегро-дифференциального уравнения переноса излучения через случайную среду. Формулировки рандомизированных алгоритмов далее будут связываться с такой задачей переноса частиц, хотя они применимы также для любых численно реализуемых ансамблей $\left\{ \Omega \right\}$ и параметров $\sigma $.

Метод двойной рандомизации определяется легко проверяемым соотношением (см. [8])

Здесь $\Omega $ – траектория кванта излучения, построенная для реализации среды с плотностью $\sigma $.

Согласно правилу повторного осреднения (то есть фактически по теореме Фубини) соотношение (1) реализуется следующим образом: строится реализация случайной среды (то есть, вообще говоря, поля $\sigma $) и затем в этой фиксированной среде строится траектория $\Omega $, которая дает вклад в статистическую оценку величины (1).

Практически весьма важно, что при построении несмещенной оценки момента (1) для данной реализации $\sigma $ достаточно строить лишь одну элементарную оценку функционала. Отметим, что при попадании траектории $\Omega $ в подобласть среды с уже выбранными значениями $\sigma $ их нельзя выбирать заново, иначе возникает “ошибка перевыбора” [6]. Согласно теореме Фубини, правая часть соотношения (1) должна оставаться конечной после замены $\xi (\Omega ;\sigma )$ на $\left| {\xi (\Omega ;\sigma )} \right|$. При $\xi (\Omega ;\sigma ) \geqslant 0$ соотношение (1) выполняется в любом случае.

2. Трудоемкость алгоритма двойной рандомизации для оценки величины $EJ(\sigma ) = {{E}_{{(\Omega ;\sigma )}}}\xi (\Omega ;\sigma )$ можно уменьшить, используя серию условно-независимых траекторий ${{\left\{ {{{\Omega }_{k}}} \right\}}_{{k = 1, \ldots ,n}}}$, которые строятся для фиксированного $\sigma $, то есть используя оценку метода “расщепления” (см. [4])

При этом Среднее число вычислительных операций здесь определяется формулой ${{T}_{n}} = {{t}_{0}} + n{{t}_{1}}$, где ${{t}_{0}}$ соответствует реализации $\sigma $, а ${{t}_{1}}$ – реализации $\Omega $. Минимизирующее (с точностью до перехода к целой части) величину трудоемкости ${{S}_{n}} = D{{\varsigma }_{n}} \times {{T}_{n}}$ [4] значение $n$ равно причем

В реальных задачах аналитические оценки коэффициентов в формуле (2) затруднительны, поэтому их целесообразно оценивать, как указано в [5], на основе предварительной статистический оценки величин $D{{\zeta }_{n}}$, ${{T}_{n}}$для двух значений $n = {{n}_{1}},{{n}_{2}}$ (по возможности близких к ${{n}_{{{\text{opt}}}}}$), то есть путем решения уравнений

Методика, основанная на (2), требует уточнения в тех случаях, когда “достраивание” реализации $\sigma $ происходит при последовательном моделировании траекторий ${{\Omega }_{k}}$ [9]. При этом величина трудоемкости ${{S}_{n}}$ нелинейно, и в реальных задачах сложно,зависит от $n$ так, что эффективное значение отношения ${{t}_{0}}{\text{/}}{{t}_{1}}$, и тем самым ${{n}_{{{\text{opt}}}}}$, уменьшается. Однако можно предположить, что в представлении ${{T}_{n}} = {{t}_{0}} + n \times \varphi (n) \times {{t}_{1}}$ функция $\varphi (n)$ в некоторой окрестности ${{n}_{{{\text{opt}}}}}$ меняется существенно слабее, чем $n$; следовательно, уравнения (4) с выражением (2) могут эффективно уточнять оценку величины ${{n}_{{{\text{opt}}}}}$. Например, если после моделирования траектории ${{\Omega }_{k}}$ число операций для построения ${{\Omega }_{{k + 1}}}$ уменьшается на сравнительно малую величину ${{t}_{2}}$, то можно положить $\varphi (n) = 1 - {{t}_{2}}(n - 1){\text{/}}2$, причем ${{n}_{{{\text{opt}}}}}$ уменьшается вследствие связанного с этим преобразованием уменьшения ${{t}_{0}}$.

3. Рассмотрим теперь оценку плотности распределения $f(x)$ с помощью численного моделирования параметра $\sigma $ и соответствующих траекторий $\Omega $. Практически эффективной для этой цели может быть универсальная статистическая ядерная оценка Парзена–Розенблатта [10] с прямоугольным (“равномерным”) ядром (см. также [11]). Она строится на основе статистической оценки функционалов вида

где ${{I}_{\Delta }}(x)$ – индикатор гиперинтервала $\Delta = \left\{ {\left( {{{x}_{i}} - \tfrac{\delta }{2},{{x}_{i}} + \tfrac{\delta }{2}} \right)} \right\}$, $i = 1,\; \ldots ,\;m$. Предполагается, что постановка задачи допускает построение бернуллиевской оценки функционала ${{J}_{\Delta }}$ путем подсчета числа траекторий $\Omega $, невозвратно “посетивших” интервал $\Delta $. В задачах теории переноса частиц $f(x)$ – это, в частности, стохастическая плотность распределения числа частиц в точках их “гибели”, например, вследствие невозвратного вылета из среды. Имеет место статистическая оценка где ${{n}_{\Delta }}$ – число траекторий частиц, “посетивших” интервал $\Delta $ в выборке $\left\{ {({{\sigma }_{i}},{{\Omega }_{i}})} \right\}$ $(i = 1,\; \ldots ,\;N)$, так как $E{{n}_{\Delta }} = N{{J}_{\Delta }}$.

Средний квадрат погрешности оценки $f(x) \approx {{n}_{\Delta }}{\text{/}}(N{{\delta }^{m}})$ равен (см., например, [11])

с относительной погрешностью, убывающей до нуля при $\delta \to 0$ и $N{{\delta }^{m}} \to \infty $. Минимизируя (5) соответственно [11], получаем Отметим, что в [12] для оценки $f(x)$ и $f{\kern 1pt} {\text{''}}(x)$ при $m = 1$ была использована наилучшая в метрике ${{L}_{2}}$ квадратическая аппроксимация функции $f(x)$ в интервале ${{\Delta }_{0}} \supset \Delta $ с помощью полиномов Лежандра порядка 0, 1, 2. Так же, как в [13], в работе [12] для оптимизации одномерной ядерной оценки была использована “микрогруппированная” выборка с шагом $h \ll \varepsilon {\text{/}}\mathop {max}\limits_x \left| {f{\kern 1pt} {\text{'}}(x)} \right|$. При этом среднее число операций в алгоритме практически не зависит от $\delta $.

4. Пусть $\xi $ – несмещенная оценка величины $J$, соответствующая определенной численно-статистической процедуре, среднее число операций для которой равно $t$. В теории методов Монте-Карло трудоемкость $S$ такой оценки определяется величиной $tD\xi $ [1]–[4], так как для достижения среднеквадратической погрешности $\varepsilon $ необходимо $N = D\xi {\text{/}}{{\varepsilon }^{2}}$ реализаций указанной процедуры.

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

Лемма 1. Если средний квадрат погрешности статистической функциональной оценки равен $D{{N}^{{ - \alpha }}}$, где $N$ – объем выборки, а среднее число операций для вычисления выборочного значения оценки равно $t$, то трудоемкость $S$ оценки определяется выражением

Доказательство. По определению (см. [14]), трудоемкость вычислений – это среднее число $S$ вычислительных операций, необходимых для достижения заданной погрешности $\varepsilon $. В условиях леммы ${{\varepsilon }^{2}} = D{{N}^{{ - \alpha }}}$, откуда $N = {{D}^{{1/\alpha }}}{{\varepsilon }^{{ - 2/\alpha }}}$ и

Лемма 1 доказана.

Следующее утверждение определяет соответствующий критерий оптимальности статистической функциональной оценки.

Лемма 2. Если средний квадрат погрешности статистической функциональной оценки с параметром $\beta $ равен $D(\beta ){{N}^{{ - \alpha }}}$, а среднее число операций для вычисления выборочного значения оценки равно $t(\beta )$, где $N$ – объем выборки, то оптимальное (по критерию трудоемкости вычислений) значение $\beta $ равно

2. ОПТИМИЗАЦИЯ РАНДОМИЗИРОВАННОЙ ОЦЕНКИ С РАСЩЕПЛЕНИЕМ

В случае дополнительного осреднения методом двойной рандомизации в выражении (5) можно полагать $f(x) = Ef(x;\sigma )$. Целью настоящего раздела работы является минимизация трудоемкости оценки соответственно этому выражению рассмотренным в разд. 1 методом расщепления с параметром $n$. Здесь целесообразно осреднить соотношение (5) по $x$, то есть по аналогии с [15] рассматривать величину

где заменив $N$ на $Nn$ в предположении асимптотической ограниченности $n$.

Теорема 1. Минимальное значение величины

достигается припричем

Доказательство. Поскольку предполагается независимость ${{t}_{0}}$, ${{t}_{1}}$ от $\delta $, то согласно (3) величина $\delta {\text{*}}$ после несложных выкладок получается из уравнения

Значение $n{\text{*}}$ получается из (2) при ${{d}_{0}} = {{f}_{0}}{{(\delta *)}^{4}}$, ${{d}_{1}} = d{\text{/}}(N{{(\delta *)}^{m}})$. Это завершает доказательство теоремы 1.

Напомним, что в случае аналоговой бернуллиевской оценки функционала ${{J}_{\Delta }}$ используется значение $d = \int {f(x)dx} $. Для несмещенной весовой модификации моделирования здесь, соответственно [12], $f(x)$ заменяется на плотность распределения квадрата веса${{f}_{{{{w}^{2}}}}}(x)$, если вспомогательный вес частицы $w$ ограничен, то есть $w \leqslant C < + \infty $. Заметим, что в задаче о переносе частиц можно не “разыгрывать” поглощение; при этом вспомогательный вес равен $exp( - {{\tau }_{c}})$, где ${{\tau }_{c}}$ – “оптическая” длина траектории относительно коэффициента поглощения (см., например, [6], [7]).

Перейдем теперь к точной формулировке задачи оптимизации рассматриваемой рандомизированной ядерной оценки. Вследствие леммы 2 справедлива

Теорема 2. Трудоемкость рандомизированной ядерной оценки функции $f(x)$ асимптотически по $N$ определяется параметрами $n_{0}^{ * }$, $\delta _{0}^{ * }$, минимизирующими величину

При этом сохраняется асимптотика ${{\varepsilon }^{2}}(Nn_{0}^{ * },\delta _{0}^{ * }) \asymp {{N}^{{ - 4/(4 + m)}}}$.

Представленную в теореме 2 задачу минимизации можно решать численно. В первом приближении $n_{0}^{ * } \approx n*$, $\delta _{0}^{ * } \approx \delta {\text{*}}$, так как ${{({{t}_{0}} + n{{t}_{1}})}^{{m/(4 + m)}}}$ – сравнительно слабо меняющаяся функция аргумента $n$ при $m = 1,2$.

Отметим, что в случае указанной в п. 1 нелинейной зависимости величины ${{T}_{n}} = {{t}_{0}} + n{{t}_{1}}$ от $n$ значение $n{\text{*}}$ и, тем самым, отношение ${{t}_{0}}{\text{/}}{{t}_{1}}$ можно уточнить с помощью численной оптимизации алгоритма расщепления для функционала $J = \int {f(x)dx} $ на основе соотношений (4).

Полученные результаты соответствуют случаю “равномасштабных” координат вектора $x = ({{x}^{{(1)}}},\; \ldots ,\;.{{x}^{{(m)}}})$. В противном случае следует, как обычно, использовать масштабирование: ${{\delta }^{{(i)}}} = {{c}_{i}}\delta $ $(i = 1,\; \ldots ,\;m)$. При этом в полученных соотношениях $f(x)$ заменяется на ${{f(x)} \mathord{\left/ {\vphantom {{f(x)} {\prod\nolimits_{i = 1}^m {{{c}_{i}}} }}} \right. \kern-0em} {\prod\nolimits_{i = 1}^m {{{c}_{i}}} }}$, а $f_{i}^{{''}}(x)$ – на $c_{i}^{2}f{\kern 1pt} {\text{''}}(x)$.

Заметим, что предложенный Н.Н. Ченцовым для использования в рамках численного статистического моделирования рандомизированный проекционный метод [16] соответственно (1) переносится на оценку осредненных распределений. Однако, как показывают расчеты, этот метод может быть практически эффективным лишь в случае достаточной гладкости оцениваемой одномерной функции или ее отношения к некоторой вспомогательной плотности вероятностей (см., например, [17]). Многомерное обобщение при этом весьма затруднительно.

3. ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ ИСПОЛЬЗОВАНИЯ РАСЩЕПЛЕНИЯ С n = n*

Для исследования эффективности расщепления с параметром $n{\text{*}}$ следует вычислить относительное значение трудоемкости соответствующей оценки, которое в обозначениях из (7) определяется отношением $S(1,{{\delta }_{0}}){\text{/}}S(n*,\delta *)$.

Теорема 3. Выполняется равенство

Доказательство. Согласно сказанному выше, имеем

Отсюда, используя теорему 1, получаем Далее, таким образом, теорема 3 доказана.

Удивительным является тот факт, что сложные выкладки здесь привели к столь простому выражению для сравнительной трудоемкости. Для более точной численной оптимизации расщепления может быть полезной следующая теорема, устанавливающая соотношение величин $n_{0}^{ * }$ и $n{\text{*}}$.

Теорема 4. Выполняется соотношение $n_{0}^{ * } > n{\text{*}}$.

Доказательство. Можно полагать, что $\delta = {{\delta }_{0}}(n)$ и заменить $S(n,\delta )$ на $S(n)$, причем, согласно теореме 1, имеем

Рассмотрим теперь функцию Опуская для краткости аргументы, дифференцированием по $n$ получаем

Отсюда с учетом равенства (8) получаем неравенство $S{\text{'}}(n*) < 0$, что завершает доказательство теоремы 4.

В заключение заметим, что полученные результаты непосредственно распространяются на случай оценки осредненной по $\sigma $ плотности распределения вектора ${{\eta }_{1}}(\Omega ;\sigma ),\; \ldots {{\eta }_{k}}(\Omega ;\sigma )$, $k \leqslant m$, с помощью использования соответствующего индикатора ${{J}_{\Delta }}(x)$. Пусть, например, необходимо оценить функцию $f(x) = Ef(x;\sigma )$, причем $f(x;\sigma )$ – условная плотность распределения случайной величины $\mu = (\omega ,{\mathbf{n}})$, где $\omega $ – направление частицы, вылетающей из среды со случайной плотностью $\sigma $, а ${\mathbf{n}}$ – граничная нормаль. Здесь в качестве ${{J}_{\Delta }}(x)$ следует использовать индикатор интервала $\Delta = (\mu - \delta {\text{/}}2 < \mu (x) < \mu + \delta {\text{/}}2)$, где $x$ – фазовые координаты точки вылета.