Журнал вычислительной математики и математической физики, 2021, T. 61, № 9, стр. 1431-1446

1. РЕШЕНИЕ ЗАДАЧИ ОПЕРАТИВНОГО ЧИСЛЕННОГО ФУНКЦИОНАЛЬНОГО ПРИБЛИЖЕНИЯ ВЕРОЯТНОСТНОЙ ПЛОТНОСТИ

В [1]–[4] мы рассмотрели следующую задачу.

Задача 1. По заданной выборке $\left\{ {{{\eta }_{1}},\; \ldots ,\;{{\eta }_{{\hat {n}}}}} \right\}$ построить численное функциональное приближение неизвестной плотности ${{f}_{\eta }}({\mathbf{x}})$, ${\mathbf{x}} \in X \subset {{\mathbb{R}}^{d}}$, случайной величины (вектора) $\eta \in X \subset {{\mathbb{R}}^{d}}$ с заданным уровнем погрешности $L > 0$ и с наименьшими вычислительными затратами $S$.

Сформулированная задача может быть весьма актуальной при оперативной обработке больших данных, т.е. при $\hat {n} \gg 1$ (см., например, [5], [6]).

Наше предложение состоит в том, чтобы, используя относительно небольшую часть $\left\{ {{{\eta }_{1}},\; \ldots ,\;{{\eta }_{n}}} \right\}$ выборки $\left\{ {{{\eta }_{1}},\; \ldots ,\;{{\eta }_{{\hat {n}}}}} \right\}$ (здесь $n \ll \hat {n}$), применить следующий функциональный вычислительный ядерный алгоритм.

Алгоритм 1 (см. [1]–[4]). Вычисляем значения

и приближаем функцию ${{f}_{\eta }}({\mathbf{x}})$ по формуле

Величина ${{Z}_{n}}({{{\mathbf{x}}}_{i}})$ из соотношения (1.1) представляет собой ядерную оценку плотности ${{f}_{\eta }}({\mathbf{x}})$ (см. [7]) в точке ${{{\mathbf{x}}}_{i}}$, являющейся узлом аппроксимационной сетки

введенной в области $X$. При этом ${{\kappa }^{{({\mathbf{x}})}}}({\mathbf{y}})$ – некоторая финитная параметрическая, одинаковая по форме для всех значений параметра ${\mathbf{x}}$, ядерная функция. Эта функция выбирается таким образом, чтобы выполнялись соотношения

В свою очередь, выражение (1.2) представляет собой классическое сеточное приближение на сетке (1.3) (см., например, ${{\S}}\,2$ гл. 2 [8], ${{\S}}\,2$, гл. 2) с аппроксимационным базисом

(эти функции определенным образом связаны с сеткой (1.3)), а значения представляют собой аппроксимационные коэффициенты, являющиеся определенными комбинациями приближенных значений функции ${{f}_{\eta }}({\mathbf{x}})$ в узлах сетки (1.3) (см. соотношение (1.4)); чаще всего

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

Важным является следующее обстоятельство.

Замечание 1 (см. [1]–[4]). Конструкция функционального вычислительного ядерного алгоритма приближения вероятностной плотности ${{f}_{\eta }}({\mathbf{x}})$ близка к схеме построения функционального вычислительного ядерного статистического алгоритма приближения решения интегрального уравнения Фредгольма II рода (см. [11], разд. 8.6), и поэтому при исследовании алгоритма 1 целесообразно использовать соображения теории условной оптимизации функциональных алгоритмов метода Монте-Карло (см. прежде всего [12]–[21]).

Данная работа посвящена детализации соображений из замечания 1.

2. ОБЩАЯ СХЕМА УСЛОВНОЙ ОПТИМИЗАЦИИ

Общая схема оптимизации из [12]–[21] выглядит следующим образом. Ставится задача согласованного выбора параметров $M$ и $n$ (число узлов и выборочных значений из соотношений (1.1)) алгоритма 1, обеспечивающего заданный уровень $L > 0$ погрешности приближения (1.2) при минимальных вычислительных затратах $S(M,n)$.

Метод 1. Строится верхняя граница $U{{P}^{{(\mathbb{B})}}}(M,n)$ погрешности алгоритма ${{\delta }^{{(\mathbb{B})}}}(M,n)$ для используемого нормированного функционального пространства $\mathbb{B}(X)$, зависящая от параметров $M$ и $n$:

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

один из параметров (например, $n$) выражается через другой: $n = {{\psi }^{{(L)}}}(M)$. Это соотношение подставляется в выражение для затрат $S(M,n)$ (которое тоже зависит от параметров $M$ и $n$). В результате получается функция ${{\tilde {S}}^{{(\mathbb{B},L)}}}(M)$ одного переменного $M$, которая исследуется на минимум с помощью известных приемов математического или численного анализа. Найденные значения $M_{{{\text{min}}}}^{{(\mathbb{B})}}(L) = M_{{{\text{opt}}}}^{{(\mathbb{B})}}(L)$, $n_{{{\text{opt}}}}^{{(\mathbb{B})}}(L) = \psi [M_{{{\text{opt}}}}^{{(\mathbb{B})}}(L)]$ объявляются условно-оптимальными параметрами соответствующего функционального алгоритма.

“Условность” оптимизации по методу 1 связана с тем, что в левой части уравнения вида (2.2) используется не сама погрешность алгоритма ${{\delta }^{{(\mathbb{B})}}}(M,n)$, а ее верхняя граница $U{{P}^{{(\mathbb{B})}}}(M,n)$. К слову, оценка качества того или иного алгоритма по верхней границе погрешности используется в подавляющем числе теоретических рассуждений вычислительной математики (см., например, [8]), т.е. по сути речь идет об “условной оптимальности” исследуемых численных схем.

При изучении погрешности ${{\delta }^{{(\mathbb{B})}}}(M,n)$ необходимо выбрать как соответствующее нормированное функциональное пространство $\mathbb{B}(X)$, так и вероятностный смысл выполнения неравенства вида (2.1) (ведь ${{\delta }^{{(\mathbb{B})}}}(M,n)$ является случайной величиной). Следуя канонам классического численного анализа (см., например, [8]), в качестве пространств $\mathbb{B}(X)$ будем рассматривать ${{\mathbb{L}}_{2}}(X)$ и $\mathbb{C}(X)$.

Для хорошо разработанного (см. в первую очередь [12]–[21]) ${{\mathbb{L}}_{2}}$-подхода выбирается сходимость в среднем погрешности

к нулю при $M,n \to \infty $ и строятся оценки сверху $U{{P}^{{({{\mathbb{L}}_{2}})}}}(M,n)$ такие, что

Для $\mathbb{C}$-подхода (см. [16]–[21]) величина

ограничивается сверху по вероятности для некоторого достаточно малого $\varepsilon > 0$.

Отметим, что для ${{\mathbb{L}}_{2}}$-подхода используется относительно “слабая” интегральная норма ${{\left\| {\, \cdot \,} \right\|}_{{{{\mathbb{L}}_{2}}(X)}}}$ пространства ${{\mathbb{L}}_{2}}(X)$ и “сильная” вероятностная сходимость погрешности к нулю (в среднем). В свою очередь в $\mathbb{C}$-подходе для “жесткой” нормы ${{\left\| {\, \cdot \,} \right\|}_{{\mathbb{C}(X)}}}$ используется относительно “слабая” сходимость (по вероятности).

Для $\mathbb{C}$- и ${{\mathbb{L}}_{2}}$-подходов можно разделить погрешность на три компоненты: детерминированную $\delta _{{{\text{det}}}}^{{(\mathbb{B})}}(M)$, стохастическую $\delta _{{{\text{stoch}}}}^{{(\mathbb{B})}}(M,n)$ и компоненту смещения $\delta _{{{\text{bias}}}}^{{(\mathbb{B})}}(M)$.

В частности, для $\mathbb{C}$-подхода, согласно неравенству треугольника,

Здесь

Для ${{\mathbb{L}}_{2}}$-подхода, применяя неравенство Коши–Буняковского и теорему Фубини, получаем

Далее заметим, что Здесь учтено, что ${\mathbf{E}}{{L}^{{(M,n)}}}{{\tilde {f}}_{\eta }}({\mathbf{x}}) = {{L}^{{(M)}}}{{\bar {f}}_{\eta }}({\mathbf{x}})$ (см. соотношения (1.2), (1.4), (2.7)). Используя очевидное неравенство ${{(a + b)}^{2}} \leqslant 2{{a}^{2}} + 2{{b}^{2}}$, получаем

3. ВЕРХНИЕ ГРАНИЦЫ ДЛЯ ДЕТЕРМИНИРОВАННЫХ КОМПОНЕНТ $\mathbb{C}$- И ${{\mathbb{L}}_{2}}$-ПОГРЕШНОСТЕЙ

В книгах и диссертациях [14], [16]–[21] и сопутствующих статьях [22]–[24] по теории условной оптимизации по методу 1 приведены весомые аргументы в пользу использования в качестве функций (1.5) базиса Стренга–Фикса (см. [25])

для равномерной сетки (1.3) с шагом $h$ по каждой координате, т.е. с узлами вида ${{{\mathbf{x}}}_{i}} = (j_{i}^{{(1)}}h,\; \ldots ,\;j_{i}^{{(d)}}h)$ в прямоугольной области $X$; здесь $\{ j_{i}^{{(k)}},i = 1,\;2,\; \ldots ,\;M,\;k = 1,\;2,\; \ldots ,\;d\} $ – целые числа.

Для такой сетки число $M$ узлов сетки пропорционально величине ${{h}^{{ - d}}}$, т.е. существует константа ${{H}^{{(h \to M)}}}$ такая, что

В формуле (3.1) функция ${{\beta }^{{(r)}}}(u)$ представляет собой $B$-сплайн $r$-го порядка (см., например, [25]), который можно определить с помощью следующих рекуррентных формул:

Учитывая вид линейных преобразований координат $\{ {{x}^{{(1)}}},\; \ldots ,\;{{x}^{{(d)}}}\} $ в формуле (3.1), важным является требование того, что сплайн ${{\beta }^{{(r)}}}(u)$ должен представлять собой кусочно-гладкую функцию, причем границы разбиения интервала, на котором ${{\beta }^{{(r)}}}(u) > 0$, на соответствующие полуинтервалы должны быть целыми точками. Это требование выполняется, когда параметр $r$ является натуральным нечетным числом: $r = 1,\;3,\;5,\; \ldots $ . Чаще всего используется $B$-сплайн первого порядка (или “функция-крышка”)

При этом приближение (2.3) с базисными функциями $\{ {{\chi }^{{(i,1)}}}({\mathbf{x}});\;i = 1,\;2,\; \ldots ,\;M\} $ и с коэффициентами простейшего вида (это аналог соотношения (1.6)) определяет мультилинейную аппроксимацию плотности ${{f}_{\eta }}({\mathbf{x}})$ (см. [14], [16]–[21]).

На основании подходов из [20], [25] можно получить следующие результаты.

Утверждение 1. Если функция ${{f}_{\eta }}({\mathbf{x}})$ принадлежит пространству ${{\mathbb{C}}^{{r + 1}}}(X)$ и в приближении (2.6) функции ${{f}_{\eta }}({\mathbf{x}})$ используются базисные функции $\{ {{\chi }^{{(i,r)}}}({\mathbf{x}});\;i = 1,\;2,\; \ldots ,\;M\} $ из (3.1), то найдутся такие коэффициенты (3.5), что справедливы неравенства

причем константы $\tilde {H}_{r}^{{({{\mathbb{C}}^{s}})}}$ не зависят от ${{f}_{\eta }}({\mathbf{x}})$ и $h$.

Утверждение 2. Если функция ${{f}_{\eta }}({\mathbf{x}})$ принадлежит пространству $\mathbb{W}_{2}^{{r + 1}}(X)$ и в приближении (2.6) используются базисные функции (3.1), то найдутся такие коэффициенты (3.5), что справедливы неравенства

причем константы $\tilde {H}_{r}^{{(\mathbb{W}_{2}^{s})}}$ не зависят от ${{f}_{\eta }}({\mathbf{x}})$ и $h$.

Заметим, что

Из утверждений 1 и 2 и соотношения (3.2) получаем

Особую проблему представляет собой нахождение коэффициентов (3.5), обеспечивающих выполнение неравенств вида (3.7). Относительно несложным является лишь случай $r = 1$.

Теорема 1 (см. [20], [25]). Выбор базисных функций (3.1), (3.4) и коэффициентов (3.6) мультилинейной аппроксимации обеспечивает выполнение неравенств вида (3.7)

при ${{f}_{\eta }} \in {{\mathbb{C}}^{2}}(X)$ и ${{f}_{\eta }} \in \mathbb{W}_{2}^{2}(X)$ соответственно. При этом

В случае $r > 1$ выбор подходящих коэффициентов (3.5), обеспечивающих выполнение неравенств вида (3.7), существенно более сложен. Отметим, что в [20] получены коэффициенты (3.5), представляющие собой специальные (существенно отличные от (3.6)) комбинации значений $\{ {{f}_{\eta }}({{{\mathbf{x}}}_{1}}),\; \ldots ,\;{{f}_{\eta }}({{{\mathbf{x}}}_{M}})\} $ и обеспечивающие оптимальный порядок ${{M}^{{ - 4/d}}}$ (см. утверждения 1, 2 и соотношение (3.2)) величин $\delta _{{{\text{det}}}}^{{(\mathbb{C})}}(M)$, $\delta _{{{\text{det}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ для ${{f}_{\eta }} \in {{\mathbb{C}}^{4}}(X)$ и ${{f}_{\eta }} \in \mathbb{W}_{2}^{4}(X)$, соответственно, и для функций (3.1) при $r = 3$. Здесь

Дальнейшие рассуждения показывают, что для исследуемого алгоритма 1 достаточно (и даже необходимо) рассмотрение случая $r = 1$ (т.е. использование мультилинейного восполнения (3.1), (3.4)) и выполнение принципиальных соотношений (3.6) и (3.8). В частности, это обусловлено тем, что компоненты смещения $\delta _{{{\text{bias}}}}^{{(\mathbb{C})}}(M)$, $\delta _{{{\text{bias}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ имеют порядок ${{M}^{{ - 2/d}}}$ (как и детерминированные компоненты из соотношений (3.8)), во всяком случае, для простейшей “ядерной” функции (4.5) (см. далее разд. 4).

4. ВЕРХНИЕ ГРАНИЦЫ ДЛЯ КОМПОНЕНТ СМЕЩЕНИЯ $\mathbb{C}$- И ${{\mathbb{L}}_{2}}$-ПОГРЕШНОСТЕЙ

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

Утверждение 3. Базис Стренга–Фикса (3.1) является разложением единицы, т.е. $\sum\nolimits_i {{{\chi }^{{(i,r)}}}} ({\mathbf{x}}) \equiv 1$ для всех ${\mathbf{x}} \in {{\mathbb{R}}^{d}}$ и $r \in \mathbb{N}$.

Доказательство. Используем метод математической индукции по размерности $d$ пространства ${{\mathbb{R}}^{d}}$.

Для $d = 1$ утверждение доказываем индукцией по порядку $r$ сплайна ${{\beta }^{{(r)}}}(u)$. Пусть $r = 0$ и $x \in \mathbb{R}$. Найдется единственное целое число $\hat {i}$ такое, что $\hat {i}h - h{\text{/}}2 \leqslant x < \hat {i}h + h{\text{/}}2$. Тогда $\sum\nolimits_i {{{\beta }^{{(0)}}}} (x{\text{/}}h - i) \equiv 1$ (см. соотношения (3.3)), так как в этой сумме все слагаемые равны нулю, кроме $\hat {i}$-го, которое равно единице.

Пусть теперь

для всех $x \in \mathbb{R}$. Тогда, согласно определению $B$-сплайна (см. соотношения (3.3)), имеем Используя индуктивное предположение (5.1) для $(x - hy)$ вместо $x$, получаем $\sum\nolimits_i {{{\beta }^{{(r)}}}} (x{\text{/}}h - i) \equiv 1$, т.е. утверждение 3 верно для $d = 1$.

Наконец, индуктивный переход по размерности $d$ следует из соотношения

Утверждение 3 доказано.

Из утверждения 3 следует, что константа Лебега ${{K}^{{({{\Xi }^{{(M)}}})}}}$ (см. [26, разд. 3.1]) для базиса Стренга–Фикса ${{\Xi }^{{(M)}}}$ вида (3.1) равна единице:

и в случае, когда в алгоритме 1 используются базис Стренга–Фикса (3.1) и коэффициенты (3.5) простейшего вида (3.6), выполнены соотношения

Заметим, что в [20] для аппроксимации Стренга–Фикса (2.6), (3.1) с кубической производящей функцией (3.10) для специального вида коэффициентов (3.5), дающих, согласно утверждениям 1 и 2, оптимальный порядок ${{M}^{{ - 4/d}}}$ величин $\delta _{{{\text{det}}}}^{{(\mathbb{C})}}(M)$, $\delta _{{{\text{det}}}}^{{({{\mathbb{L}}_{2}})}}(M)$, получено ${{K}^{{({{\Xi }^{{(M)}}})}}} = 3$ (что тоже неплохо). А вот, например, для одномерной интерполяции Лагранжа константа Лебега ${{K}^{{({{\Xi }^{{(M)}}})}}}$ растет с увеличением числа $M$ узлов равномерной сетки как ${{2}^{M}}$ (см. [26]).

Предположим, что $I$ является номером узла ${{{\mathbf{x}}}_{I}}$, в котором достигается максимум $ma{{x}_{{i = 1,...,M}}}\left| {{{f}_{\eta }}({{{\mathbf{x}}}_{i}}) - {\mathbf{E}}{{\kappa }^{{({{{\mathbf{x}}}_{i}})}}}\left( \eta \right)} \right|$ в правой части соотношения (4.2). Тогда

Далее последуем рекомендациям из [14], [16]–[23] и выберем кусочно-постоянный вариант ядерной функции

Здесь и ${\mathbf{x}} = ({{x}^{{(1)}}},\; \ldots ,\;{{x}^{{(d)}}})$. В этом случае имеем

Далее проводим следующие рассуждения из теории численного интегрирования, касающиеся обоснования квадратурной формулы центральных прямоугольников (см., например, [8, ${{\S}}\;1$ гл. 3]).

Утверждение 4. Если ${{f}_{\eta }} \in {{\mathbb{C}}^{2}}(X)$, то справедливо неравенство

Здесь ${\mathbf{y}} = ({{y}^{{(1)}}},\; \ldots ,\;{{y}^{{(M)}}})$.

Доказательство. Рассмотрим сначала случай $d = 1$. Здесь ${{\Delta }^{{({{x}_{I}})}}} = {\text{\{ }}y:$ ${{x}_{I}} - h{\text{/}}2 \leqslant y \leqslant {{x}_{I}} + h{\text{/}}2{\text{\} }}$. Несложно показать, что

Кроме того, согласно формуле Тейлора, имеем

где $y \in {{\Delta }^{{({{x}_{I}})}}}$ и ${{D}_{1}} \leqslant \tfrac{{{{{(y - {{x}_{I}})}}^{2}}}}{2}\mathop {max}\limits_{y \in {{\Delta }^{{({{x}_{I}})}}}} \left| {f_{\eta }^{{''}}(y)} \right|$. Тогда

Очевидная индукция по размерности $d$ дает соотношение (4.7). Утверждение 4 доказано.

Теорема 2. Для алгоритма с аппроксимационным базисом Стренга–Фикса (3.1) с образующей функцией (3.4), простейшими коэффициентами (3.6) и ядерной функцией (4.5) справедливы неравенства

при ${{f}_{\eta }} \in {{\mathbb{C}}^{2}}(X)$ и ${{f}_{\eta }} \in \mathbb{W}_{2}^{2}(X)$ соответственно. Здесь

Утверждение теоремы 2 следует из соотношений (4.2), (4.4), (4.6), (4.7) и очевидного неравенства $\delta _{{{\text{bias}}}}^{{({{\mathbb{L}}_{2}})}}(M) \leqslant \sqrt {\operatorname{mes} X} \times \delta _{{{\text{bias}}}}^{{(\mathbb{C})}}(M)$.

Сравнивая неравенства (3.8), (3.9) и (4.8), (4.9), можно сделать вывод о том, что упомянутая выше возможность рассмотрения базиса (3.1) с кубической производящей функцией (3.10) для получения порядка ${{M}^{{ - 4/d}}}$ для величин $\delta _{{{\text{det}}}}^{{(\mathbb{C})}}(M)$, $\delta _{{{\text{det}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ для рассматриваемого алгоритма 1 не имеет смысла, так как компоненты смещения имеют порядок ${{M}^{{ - 2/d}}}$ (это аналог основного вывода [24]).

Также можно высказать весьма правдоподобную гипотезу (имеющую некоторые численные подтверждения – у нас имеются соответствующие предварительные результаты тестовых расчетов, но требующую, впрочем, отдельной тщательной проверки) о том, что использование в алгоритме 1 ядерных функций ${{\kappa }^{{({\mathbf{x}})}}}({\mathbf{y}})$, отличных от (4.5) (например, оптимальных функций из работы [7]), не позволит получить верхних границ компонент смещения лучших, чем (4.8).

Из теоремы 2 и двух последних соображений вытекает следующий важный вывод.

Замечание 2. Для решения сформулированной в разделе 1 задачи 1 целесообразным является применение следующего приближения алгоритма 1 с аппроксимационным базисом ${{\Xi }^{{(M)}}}$ вида (3.1) с образующей функцией (3.4), простейшими аппроксимационными коэффициентами вида (3.6) и ядерной функцией ${{\kappa }^{{({\mathbf{x}})}}}({\mathbf{y}})$ вида (4.5):

Здесь ${{X}^{{({{\Delta }^{{({{{\mathbf{x}}}_{i}})}}})}}}({\mathbf{y}})$ – индикатор множества ${{\Delta }^{{({{x}_{i}})}}}$: при ${\mathbf{y}} \in {{\Delta }^{{({{x}_{i}})}}}$ имеем ${{X}^{{({{\Delta }^{{({{{\mathbf{x}}}_{i}})}}})}}}({\mathbf{y}}) = 1$, а при ${\mathbf{y}} \notin {{\Delta }^{{({{x}_{i}})}}}$ выполнено ${{X}^{{({{\Delta }^{{({{{\mathbf{x}}}_{i}})}}})}}}({\mathbf{y}}) = 0$.

По аналогии с теорией рандомизированных функциональных алгоритмов приближения решения интегрального уравнения Фредгольма (см., в частности, [14], [16]–[24]), назовем вариант алгоритма 1, описанный в теореме 2 и в замечании 2 (см. соотношение (4.10)), многомерным аналогом метода полигона частот. Такой термин здесь более чем удачен, так как в одномерном случае $d = 1$ алгоритм 1 из теоремы 2 и замечания 2 дает в точности полигон частот.

В связи с замечанием 2 мы займемся выводом формул для условно-оптимальных параметров именно для многомерного аналога метода полигона частот (4.10). Для простоты везде далее будем полагать, что

5. ВЕРХНЯЯ ГРАНИЦА ДЛЯ СТОХАСТИЧЕСКОЙ КОМПОНЕНТЫ ${{\mathbb{L}}_{2}}$-ПОГРЕШНОСТИ

Построим верхние границы стохастических компонент погрешности $\delta _{{{\text{stoch}}}}^{{({{\mathbb{L}}_{2}})}}(M,n)$, $\delta _{{{\text{stoch}}}}^{{(\mathbb{C})}}(M,n)$ из соотношений (2.5), (2.8) для многомерного аналога метода полигона частот (4.10). Докажем сначала следующее вспомогательное утверждение.

Утверждение 5. Для приближения

справедливо соотношение

Доказательство. С учетом известного тождества ${\mathbf{D}}\tau = {\mathbf{D}}(\tau + C)$, $C = {\text{const}}$, справедливого для любой случайной величины $\tau $, имеем

где $\mathop {\tilde {\zeta }}\nolimits_j^{(i)} = {{\kappa }^{{({{{\mathbf{x}}}_{i}})}}}({{\eta }_{{\mathbf{j}}}}) - {\mathbf{E}}{{\kappa }^{{({{{\mathbf{x}}}_{i}})}}}({{\eta }_{{\mathbf{j}}}})$ (здесь мы трактуем ${{\eta }_{1}},\; \ldots ,\;{{\eta }_{{\mathbf{n}}}}$ как независимые, одинаково распределенные – как $\eta $ – случайные величины).

Зафиксируем ${\mathbf{x}} = {{{\mathbf{z}}}_{0}} \in X$. Учитывая, что ${\mathbf{E}}\mathop {\tilde {\zeta }}\nolimits_j^{(i)} = 0$ для всех $i = 1,\;2,\; \ldots ,\;M$ и $j = 1,\;2,\; \ldots ,\;n$, получаем равенство

где ${{\tilde {\zeta }}^{{(i)}}}$ – случайная величина, распределенная как независимые для разных $j$, одинаково распределенные величины $\tilde {\zeta }_{j}^{{(i)}}$.

Далее, применяя вероятностный аналог неравенства Коши–Буняковского

и соотношения ${\mathbf{D}}{{\tilde {\zeta }}^{{(i)}}} = {\mathbf{D}}{{\kappa }^{{({{{\mathbf{x}}}_{i}})}}}(\eta )$ для $i = 1,\;2,\; \ldots ,\;M$ получаем оценку сверху Наконец, в силу утверждения 3 и произвольности ${{{\mathbf{z}}}_{0}} \in X$, имеем Утверждение 5 доказано.

Теорема 3. Для многомерного аналога полигона частот (4.10) (т.е. для алгоритма 1 с аппроксимационным базисом ${{\Xi }^{{(M)}}}$ вида (3.1) с образующей функцией (3.4), простейшими аппроксимационными коэффициентами вида (3.6) и ядерной функцией ${{\kappa }^{{({\mathbf{x}})}}}({\mathbf{y}})$ вида (4.5)) при выполнении условия (4.11) справедливы неравенства

где

Доказательство. Для ядерной функции (4.5) имеем

Учитывая условие (4.11) и теорему о среднем для определенного интеграла от непрерывной функции, получаем, что существует ${{{\mathbf{z}}}_{i}} \in {{\Delta }^{{({{{\mathbf{x}}}_{i}})}}}$ такое, что

и

В силу ограниченности функции ${{f}_{\eta }} \in {{\mathbb{C}}^{2}}(X)$ на компакте $X$, из утверждения 5 и соотношений (3.2), (5.3) получаем соотношения (5.1), (5.2). Теорема 3 доказана.

6. ВЕРХНЯЯ ГРАНИЦА ДЛЯ СТОХАСТИЧЕСКОЙ КОМПОНЕНТЫ $\mathbb{C}$-ПОГРЕШНОСТИ

Для $\mathbb{C}$-подхода, используя соотношение (4.3), получаем

С учетом соотношений (5.1) имеем

В силу центральной предельной теоремы для одинаково распределенных случайных величин, суммы

сходятся (по распределению) к стандартным гауссовским величинам ${{\gamma }^{{(i,M)}}}$ при $n \to \infty $. Следовательно, при фиксированном $M$ для любого ${{\varepsilon }_{1}} > 0$ найдется такое натуральное ${{N}_{1}}({{\varepsilon }_{1}},M)$, что при $n > {{N}_{1}}({{\varepsilon }_{1}},M)$ выполнено где Здесь ${{B}_{1}}({{\varepsilon }_{1}})$ – небольшая положительная константа, такая, что $\int_{ - {{B}_{1}}({{\varepsilon }_{1}})}^{{{B}_{1}}({{\varepsilon }_{1}})} {\tfrac{{{{e}^{{ - {{{v}}^{2}}/2}}}d{v}}}{{\sqrt {2\pi } }}} \approx 1 - {{\varepsilon }_{1}}$.

Подробному изучению распределений и асимптотических (при $M \to \infty $) свойств величин типа ${{Y}^{{(M)}}} = {{\max }_{{i = 1,...,M}}}{{\gamma }^{{(i,M)}}}$ посвящена монография (см. [27]). Некоторая особенность рассматриваемого здесь случая состоит в том, что величины $\{ {{\gamma }^{{(i,M)}}};i = 1,\;2,\; \ldots ,\;M\} $ являются зависимыми. Однако эта зависимость “ослабевает” с ростом $M$. Конкретнее, для ковариаций

справедливы соотношения Здесь использованы соотношения (3.2), (5.3), (5.4) и то, что для ядерной функции (4.5) выполнено ${{\kappa }^{{({{{\mathbf{x}}}_{{{{i}_{1}}}}})}}}(\eta ) \times {{\kappa }^{{({{{\mathbf{x}}}_{{{{i}_{2}}}}})}}}(\eta ) \equiv 0$ при ${{i}_{1}} \ne {{i}_{2}}$.

Из последнего соотношения следует, что ковариации (6.3) убывают по модулю с ростом $M$ со скоростью порядка $1{\text{/}}M$; при этом найдутся натуральное число ${{\hat {M}}_{1}}$ (например, ${{\hat {M}}_{1}} = 2{{f}_{{{\text{max}}}}}{{[{{H}^{{(h \to M)}}}]}^{d}}$) и “универсальная” константа ${{B}_{2}} > 0$ (например, ${{B}_{2}} = \sqrt 2 {{f}_{{{\text{max}}}}}{{[{{H}^{{(h \to M)}}}]}^{d}}$) такие, что

Используя технологию изучения порядковых статистик из [27] по аналогии с [22], можно получить следующий результат.

Утверждение 6 (см. [22], [27]). Если для набора стандартных нормальных случайных величин ${{\gamma }^{{(1,M)}}},\; \ldots ,\;{{\gamma }^{{(M,M)}}}$ выполнены соотношения (6.4), то асимптотическое при $M \to \infty $ распределение максимума ${{Y}^{{(M)}}} = ma{{x}_{{i = 1,...,M}}}{{\gamma }^{{(i,M)}}}$ таково, что

где

В силу симметрии стандартного нормального распределения и с учетом соотношений

получаем, что где коэффициенты ${{a}_{M}}$, ${{b}_{M}}$ имеют вид (6.5).

Из последнего соотношения следует, что для любого ${{\varepsilon }_{2}} > 0$ существуют константа ${{B}_{3}}({{\varepsilon }_{2}}) > 0$ и натуральное число ${{\hat {M}}_{2}}$ такие, что

С учетом соотношений (6.1), (6.2), (6.6) получаем следующий результат.

Теорема 4. Для многомерного аналога полигона частот (4.10) (т.е. для алгоритма с аппроксимационным базисом ${{\Xi }^{{(M)}}}$ вида (3.1) с образующей функцией (3.4), простейшими аппроксимационными коэффициентами вида (3.6) и ядерной функцией ${{\kappa }^{{({\mathbf{x}})}}}({\mathbf{y}})$ вида (4.5)) для любого $\varepsilon > 0$ существуют положительные действительные константы $H_{{{\text{stoch}},1}}^{{(\mathbb{C})}}(\varepsilon )$, $H_{{{\text{stoch}},2}}^{{(\mathbb{C})}}(\varepsilon )$ и натуральное число $\hat {M}$ такие, что для любого $M > \hat {M}$ существует натуральное число $\hat {N}(\varepsilon ,M)$ такое, что для всех $n > \hat {N}(\varepsilon ,M)$ выполнено

7. УСЛОВНО-ОПТИМАЛЬНЫЕ ПАРАМЕТРЫ ДЛЯ ${{\mathbb{L}}_{2}}$- И $\mathbb{C}$-ПОДХОДОВ

Полученные в теоремах 1–4 верхние границы (3.8), (4.8), (5.1), (6.7) для компонент погрешности ${{\delta }^{{(\mathbb{B})}}}$ позволяют вывести выражения для условно-оптимальных параметров $M_{{{\text{opt}}}}^{{(\mathbb{B})}}(L)$ и $n_{{{\text{opt}}}}^{{(\mathbb{B})}}(L)$ для $\mathbb{B}(X) = \mathbb{C}(X) \vee {{\mathbb{L}}_{2}}(X)$ по представленному в разд. 2 методу 1.

Важной отличительной особенностью алгоритма 1 является то, что затраты $S$ этого алгоритма пропорциональны величине $n \times t$, где $t$ – время определения того, в какой из кубов ${{\Delta }^{{({{{\mathbf{x}}}_{i}})}}}$ попадает очередное выборочное значение ${{\eta }_{{\mathbf{j}}}}$ случайной величины $\eta $, $i = 1,\;2,\; \ldots ,\;M$, $j = 1,\;2,\; \ldots ,\;n$.

Несложно добиться, чтобы время $t$ не зависело от числа $M$ кубов ${{\Delta }^{{({{{\mathbf{x}}}_{i}})}}}$. Таким образом, затраты оптимизированных версий алгоритма 1 явно не зависят от $M$ (наши тестовые и прикладные численные эксперименты подтверждают это утверждение).

В оптимизационных процедурах метода 1 зависимость затрат $S$ от параметра $M$ возникает из-за наличия уравнений вида (2.2), определяющих зависимость параметра $n$ (а значит, и затрат $S$) от $M$.

Рассмотрим такую процедуру для ${{\mathbb{L}}_{2}}$-подхода. Здесь требуется решать следующую оптимизационную задачу: найти значения $M_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L)$ и $n_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L)$, для которых достигается минимум

(здесь ${{S}_{0}}$ – время вычисления ${{2}^{d}}$ ненулевых слагаемых $\tilde {f}_{\eta }^{{({{{\mathbf{x}}}_{i}})}}{{\chi }^{{(i,1)}}}({\mathbf{x}})$ в сумме (4.10) для заданной точки ${\mathbf{x}} \in X$) при условии где $A_{1}^{{({{\mathbb{L}}_{2}})}} = 2{{(H_{{{\text{stoch}}}}^{{({{\mathbb{L}}_{2}})}})}^{2}} + 2{{(H_{{{\text{bias}}}}^{{({{\mathbb{L}}_{2}})}})}^{2}}$ и $A_{2}^{{({{\mathbb{L}}_{2}})}} = H_{{{\text{stoch}}}}^{{({{\mathbb{L}}_{2}})}}$ (см. соотношения (2.8), (3.8), (3.9), (4.8), (4.9), (5.1), (5.2)).

Из равенств (7.1) и (7.2) имеем

Найдем теперь минимум функции $\mathop {\tilde {S}}\nolimits^{({{\mathbb{L}}_{2}},L)} (M)$ по $M$. Продифференцируем эту функцию:

и найдем критическое значение (точку минимума), а затем, и условно-оптимальные параметры для ${{\mathbb{L}}_{2}}$-подхода:

Для $\mathbb{C}$-подхода решается задача нахождения минимума (7.1) при условии

где $A_{1}^{{(\mathbb{C})}} = H_{{{\text{det}}}}^{{(\mathbb{C})}} + H_{{{\text{bias}}}}^{{(\mathbb{C})}}$, $A_{2}^{{(\mathbb{C})}} = H_{{{\text{stoch}},1}}^{{(\mathbb{C})}}(\varepsilon )$, $A_{3}^{{(\mathbb{C})}} = H_{{{\text{stoch}},2}}^{{(\mathbb{C})}}(\varepsilon )$ (см. соотношения (2.5), (3.8), (3.9), (4.8), (4.9), (6.2), (6.7)).

Выражение в левой части уравнения (7.5) является относительно сложным, и мы применим следующий прием из [22], позволяющий заменить это выражение на более простое. Используя промежуточный результат из доказательства теоремы 1.5.3 из [27], вместо (7.5) можно рассмотреть соотношение

Несложно показать, что для любых фиксированных $\bar {M} \in \mathbb{N}$ и $a > 0$ найдется $b > 0$ такое, что при $M > \bar {M}$ выполнено

С учетом соотношений (7.6), (7.7) заменим соотношение (7.5) на следующее:

Из формул (7.1), (7.8) имеем

Дифференцируем ${{\tilde {S}}^{{(\mathbb{C},L)}}}(M)$ по $M$:

и найдем точку минимума – она же условно-оптимальное значение параметра $M$ для $\mathbb{C}$-подхода

Несколько слов о том, как выбирать параметр $b$. С одной стороны, его следует брать по возможности малым, так как множитель ${{b}^{2}}$ входит в выражение для трудоемкости (7.9). С другой стороны, для значения $M_{{{\text{opt}}}}^{{(\mathbb{C})}}(L)$ должно выполняться неравенство (7.7). Поэтому выберем $b$ из условия равенства в соотношении (7.7) при $M = M_{{{\text{opt}}}}^{{(\mathbb{C})}}(L)$:

Из соотношений (7.9)–(7.11) получаем условно-оптимальное значение параметра $n$ для $\mathbb{C}$-подхода

Наконец, для важного предельного случая $a = 0$ имеем

Замечание 3. Формулы (7.10), (7.12) заметно сложнее формул (7.3), (7.4). Поэтому на практике целесообразно использовать ${{\mathbb{L}}_{2}}$-подход к условной оптимизации алгоритма 1. Данное замечание вполне соответствует рекомендациям из [7]–[9].

Значения в правых частях соотношений (7.3), (7.4), (7.10), (7.12) могут быть нецелыми, поэтому на практике в качестве параметров $M_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L)$, $n_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L)$, $M_{{{\text{opt}}}}^{{(\mathbb{C})}}(L)$, $n_{{{\text{opt}}}}^{{(\mathbb{C})}}(L)$ следует брать ближайшие к этим значениям целые числа.

Отметим также наличие содержательной проблемы выбора или приближенного подсчета констант (3.9), (4.9), (5.2), (6.2), входящих в выражения (7.1)–(7.12).

При обсуждении результатов данной работы было отмечено, что использование в одномерном случае максимума $\mathop {{\text{max}}}\nolimits_{x \in X} \left| {\tfrac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}}{{f}_{\eta }}(x)} \right|$ в оценках сверху для детерминированной компоненты $\delta _{{{\text{det}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ и компоненты смещения $\delta _{{{\text{bias}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ (см. соотношения (3.9) и (4.9)) кажется достаточно “жестким” по сравнению с оценками вида

для классических “непрерывных” ядерных статистических оценок плотности (см., например, [9], [10]). В связи с этим возникла идея попробовать использовать следующую двухстадийную технологию выбора параметров $M$ и $n$ в одномерном случае.

Сначала выбираем параметр $M = \tilde {M}$ по формулам

Здесь учтено соотношение (2.8), а также то, что $\delta _{{{\text{det}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ примерно в три раза больше $\delta _{{{\text{bias}}}}^{{({{\mathbb{L}}_{2}})}}(M)$ (см. соотношения (3.9), (4.9)).

Далее рассматриваем уравнение

Мы сравнили формулы (7.3), (7.4) (для $d = 1$) и (7.13), (7.14) на тестовом примере приближения плотности ${{f}_{\eta }}(x) = \tfrac{{{{e}^{{1 - x}}}}}{{e - 1}}$, $0 < x < 1$. При этом выборочные значения случайной величины $\eta $ моделировались по формулам ${{\eta }_{j}} = - ln[1 - {{\alpha }_{j}}(1 - {{e}^{{ - 1}}})]$, где ${{\alpha }_{j}} \in U(0,1)$ – стандартные случайные числа (см., например, [19]).

Для уровня погрешности $L = 0.003$ формулы (7.3) и (7.4) дали значения $M_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L) = 14$, $n_{{{\text{opt}}}}^{{({{\mathbb{L}}_{2}})}}(L) = 3\,256\,877$ (максимальная погрешность приближения функции ${{f}_{\eta }}(x)$ для таких параметров получилась равной $0.002871 \lesssim L$, а время счета ${{t}_{{{\text{opt}}}}} = 0.446$ с). Аккуратные массовые расчеты показали оптимальность полученных значений (т.е. при других сочетаниях параметров $M$ и $n$ были превышены либо уровень погрешности $L$, либо время счета ${{t}_{{{\text{opt}}}}}$).

Формулы (7.13) и (7.14) для рассмотренного примера и для $L = 0.003$ дали значения $\tilde {M} = 11$ и $\tilde {n} = 3\,779\,382$. Максимальная погрешность приближения функции ${{f}_{\eta }}(x)$ для таких параметров получилась равной $0.001696$ (что заметно меньше заданного уровня $L$), а временные затраты $0.542$ с (что заметно больше величины ${{t}_{{{\text{opt}}}}}$).

8. ЗАКЛЮЧЕНИЕ

В данной работе предложен новый конструктивный ядерный алгоритм численного функционального приближения вероятностной плотности ${{f}_{\eta }}({\mathbf{x}})$ (алгоритм 1) по заданной выборке, основанный на использовании ядерной оценки (1.4) и на подходах теории численной аппроксимации функций. Использовано соображение из [1]–[4] о том, что построенный ядерный алгоритм 1 аналогичен функциональному вычислительному ядерному статистическому алгоритму для приближения решения интегрального уравнения Фредгольма II рода (см. замечание 1). В связи с этим в [1]–[4] было предложено использовать схему условной оптимизации функциональных алгоритмов, основанную на минимизации затрат $S(M,n)$ при фиксированном уровне погрешности $L$ (метод 1), для функционального ядерного алгоритма 1. В данной работе получен вид условно-оптимальных параметров (7.3), (7.4) и (7.10), (7.12) (для известных ${{\mathbb{L}}_{2}}$- и $\mathbb{C}$-подходов соответственно) для простейшего случая алгоритма 1 с аппроксимационным базисом вида (3.1) с образующей функцией (3.4), простейшими аппроксимационными коэффициентами вида (3.6) и кусочно-постоянной ядерной функцией вида (4.5) (т.е. для многомерного аналога алгоритма метода полигона частот).