ДОКЛАДЫ АКАДЕМИИ НАУК, 2015, том 461, № 5, с. 521-524
ТЕОРИЯ УПРАВЛЕНИЯ
УДК 519.651
АНИЗОТРОПНЫЕ СГЛАЖИВАЮЩИЕ СПЛАЙНЫ В ЗАДАЧАХ С ФАКТОРНЫМ ПЛАНОМ ЭКСПЕРИМЕНТА
© 2015 г. М. Г. Беляев
Представлено академиком РАН С.В. Емельяновым 23.10.2014 г. Поступило 23.10.2014 г.
DOI: 10.7868/S0869565215110031
1. ПОСТАНОВКА ЗАДАЧИ
Работа посвящена построению оценки / функции многих переменных /0 по заданной выборке
данных {х,- е О, у1 = /0(х,) + }"= 1, где О с [0, 1]4 -
некоторая область в Я4, ^ — белый шум с Е= а2.
Будем считать, что/0 е , где — анизотропное пространство Соболева с m = (т1,..., т4), тк — натуральные числа. По определению (см. [1]) (О) — это пространство локально суммируемых на О функций, имеющих на О обобщенные производ-
ные Dkkf(x) (к = 1, 2, ..., d), Dkk =
, и конеч-
ную норму
нХ2(п) '
■и
k = 1
k п (dxk )
Dk
kill^n) •
0 =
x ..
x 0d =
= {ч.... ^=(0ii...0fd)> z
•> nd , U = 1-
Будем предполагать, что план эксперимента X = {х1 }"= 1 заданной выборки данных есть пере-
Институт проблем передачи информации им. А.А. Харкевича Российской Академии наук, Москва ООО "Датадванс", Москва E-mail: belyaev@iitp.ru
сечение полного факторного плана © с областью О. Если Зх е ©: х £ X, то будем называть такой план неполным факторным.
Оценка / будет строиться с помощью тензорного произведения В-сплайнов. В-сплайны порядка т — это базис в пространстве кусочно-полиномиальных функций из Ст —1([0, 1]), степень каж-
m — 1
— это
Для упрощения обозначений далее в записи вида Ж1 (О) будем опускать (О).
Для к = 1, 2, ..., 4 определим множества 0к = = {б!, ..., 0^}, где 0 < О1! < ... < 0^ < 1. Будем называть конечное множество © с [0, 1]4 полным факторным планом с числами уровней п1, ..., п4, если
дого полинома не выше т, а С пространство непрерывно дифференцируемых функций до порядка т — 1 включительно. В-сплай-ны строятся с помощью т последовательных сверток индикаторной функции и обладают рядом удобных для анализа свойств, в частности обеспечивая наименьший среди сплайнов носитель.
Пусть Пк = {(/к — 1)Ак +1 — это равномерное
разбиение отрезка [0, 1] с шагом кк = р~к для всех к = 1, 2, ..., 4. Будем переходить от одного полиномиального участка сплайнов к другому в точках множества пк (рк — число таких участков). Пусть ¿к([0, 1]) — это пространство сплайнов порядка тк +1, заданных с помощью соответствующего разбиения пк. Обозначим тензорное произведение сплайнов как ¿„([0, 1]4), ¿т = ¿1 ® ¿2 ® ... ® Ба.
Введем 4 дискретно-непрерывных полунорм на ¿т, задаваемых следующим образом:
II®, k
'N s J(^
1/2
,xk,.,0dd )) v
i,, l ф k г
где N = Пкпк — мощность полного факторного плана ©.
Будем строить оценку / с помощью минимизации на ¿т функционала, включающего среднеквадратичную ошибку, вычисленную в точках выборки, и штрафное слагаемое:
k
d
k
522
БЕЛЯЕВ
f = arg min1У (yi - f(Xi))2 + f e 5 n*-1
J m
+
У MD
i = 1
2
©, k,
(1)
0 <^k < 1.
k = 1
При d = 1 решение (1) совпадает с классическими сглаживающими сплайнами [2]. В изотропных соболевских пространствах стандартным обобщением сглаживающих сплайнов на случай d > 1 являются thin-plate-сплайны [2]
n
arg min ± у (yt -f(x;.))2 + X у \Ш\2,
f e W 2
i = 1
д
|a|| i = m
iai i
где а = (аь ..., аё), Ба = -
1 а $
(дх ) ...(дх ) Предложенный функционал (1) является одним из возможных обобщений одномерных сглаживающих сплайнов на анизотропный многомерный случай. В ряде задач инженерного проектирования моделируемые функции /0 априори являются пространственно неоднородными, в силу чего при генерации выборки данных часто используется полный факторный план эксперимента с большим числом уровней пк для менее гладких факторов [3]. Предположение о принадлежности /0 анизотропии
ному классу IV 2 — это один из возможных способов формализовать подобные априорные инженерные знания о природе функции/О.
В одномерных и Шп-рМе-сплайнах миними-
тт/И
зация проводится по всему пространству IV2, в то время как в (1) только по его подпространству Важным требованием, идущим от приложений, является малая вычислительная сложность алгоритма, поскольку объем выборки п растет экспоненциально по ё, а сложность типичных методов (например, Шт-рМе-сплайнов) может иметь порядок 0(п3). Это соображение обусловливает использование тензорного произведения сплайнов ¿т в ряде алгоритмов, предназначенных для полного факторного плана, см. [4, 5]. Отметим, что во многих приложениях в силу содержательных ограничений используются неполные факторные планы эксперимента, поскольку значения у известны не во всех точках полного факторного плана ©. Одной из целей работы является построение вычислительно эффективного алгоритма нахождения решения задачи (1) для неполных факторных планов эксперимента.
Для оценок вида (1) известна оптимальная скорость сходимости ошибки Е ||/ — /01|2 [6, 7], в изотропном случае тк = т равная п-2т/(2т + ё). Достижение одномерными сглаживающими сплайнами оп-
тимальной скорости сходимости было доказано в [8], 1Ыгп-р1а1е-сплайнами — в [9]. Вторая цель работы состоит в исследовании асимптотического поведения ошибки Е11/ — /01| 2 предложенной оценки (1). Отметим еще несколько работ, пересекающихся с этой задачей исследования. Асимптотические свойства одного из методов для полного факторного плана были исследованы в [5], где была найдена оптимальная скорость сходимости для размерности ё = 2. Тензорное произведение сплайнов без штрафа на гладкость в оптимизируемом функционале использовалось в [10], где была доказана оптимальная скорость сходимости в изотропном случае. Отметим, что возможность использования тензорного произведения вейвлетов для получения асимптотически оптимальной оценки функции из анизотропного пространства Соболева была доказана в [11].
2. ВЫЧИСЛИТЕЛЬНО-ЭФФЕКТИВНЫЙ АЛГОРИТМ
Введем дополнительные предположения о структуре полного факторного плана эксперимента ©. Здесь и далее будем считать, что для всех к = 1, 2, ..., ё,
= 1, 2, ...,рк е ((]к- 1)кь]ккк].
Для построения алгоритма нам потребуются базовые операции с тензорами. Сложение, скалярное произведение, поэлементное умножение тензоров определены стандартным образом. Будем говорить, что тензор М с размерами и1 х ... .х 1к х ... х иё — это результат умножения тензора ^ (с размерами и1 х ... х ик х ... х иё) по измерению к на (ик х ?к)-матрицу С и обозначать М = = ^ ®к С, если
М'1> •••> 'к - 1,1 'к + 1, •••> Iк>•••> '¡¡С'к>1.
'к = 1
Также будем использовать следующее обозначение для умножения тензора М на матрицы Ск, к = = 1, 2, ..., ё (заметим, что умножение по разным измерениям коммутативно):
М ск - М ®1 С1 ®2 . С. Построим индикатор попадания точек полного факторного плана © = {Zil, ^}= 1 в план
эксперимента X заданной выборки данных. Пусть Ж — тензор с размерами п1 х п2 х ... х пл, причем
^ ....., е X
.0, 6 0\X.
i
k
u
k
d
АНИЗОТРОПНЫЕ СГЛАЖИВАЮЩИЕ СПЛАЙНЫ
523
Построим тензор наблюдений в точках полного плана ЗД с размерами п1 х п2 х ... х п4, дополнив
заданные значения {у^ }"= 1 нулями, т.е.: | у, Зх( ■ Х1 = ,
% i
k.Pt
где Ек — это диагональные (рк х рк)-матрицы, ¥к — невырожденные (рк хрк)-матрицы, а ик — (пк хрк)-матрицы с ортонормальными столбцами. Заметим, что описанное разложение — это матричный аналог перехода к базису Деммлера—Райнша [13].
Определим тензор Э с размерами р1 х р2 х ... х р4 как
Пусть Ak = {tyj } / = 1 — это базис в пространстве одномерных сплайнов Sk для всех k = 1, 2, ..., d. Тогда произвольная функция f е Sm может быть представлена в виде
Рь • ••'Pd
/ = z Mw-^
/i, =i
где M — тензор коэффициентов разложения с размерами p1 х p2 х ... х pd.
Для всех k = 1, 2, ..., d обозначим матрицу значений функций из базиса Ak в точках из множества 0k как Rk,
П ( кГ1Лк\->"к'Рк
Rk = {Ък(h)}ik'/k= i,
а матрицу скалярного произведения производных функций из Ak как Hk,
( ЛРк' Рк
тт I Г T\mk к кл „Щ к , к., к I
Hk = j J D yJk (x ) D (x ) dx j. . j 0' i = 1 Матрицы Rk имеют размеры nk х pk, а матрицы Hk — размеры pk х pk.
Минимизируемый функционал в (1) может быть переписан как функция от тензора коэффициентов разложения M:
1 d d
Qa(M) = 1 <(- M ®к Rk) o °W, - M ®к Як)2 + n 1 1
d d nd
+ Z ^N<M ®i (S„H + (1 - 5к;)RfR,), M)2,
к = 1
где O — поэлементное умножение тензоров, а dk¡ — символ Кронекера; QA(M) является квадратичной функцией от коэффициентов разложения M, ее размерность составляет P = nkpk. Минимизация квадратичной формы общего вида требует O(P3) операций. В типичной ситуации P ~ n, поэтому, с учетом экспоненциального роста объема выборки по d, во многих практических задачах такая вычислительная сложность оказывается неприемлемо высокой. Построим алгоритм, учитывающий специальный вид квадратичной формы в рассматриваемом случае.
С помощью обобщенного сингулярного разложения [12] представим Hk и Rk в следующем виде:
Hk = VkEkV\; Як = UkVk,
£ ^k^ ®, (&kEi + (i - 5k,)Ip),
(2)
k = i
где $ — тензор с размерамир1 х р2 х ... х р4, состоящий из 1, а 1р — единичная (р1 х рг)-матрица. Введем замену переменных
n <
k Vk) © э1/2 k Vk) © (Э + n"^)
2 ,
1/2
% Y n > —, 2
где возведение тензора в степень выполняется поэлементно. Пусть QB — это квадратичная форма QA после замены переменных, QB(M) = QA(^).
Лемма 1. Гессиан квадратичной формы QB не вырожден и имеет не более min(N — n + 1, n + 1) различных собственных чисел.
Теорема 1. Минимизация квадратичной формы QB с помощью метода сопряженных градиентов требует не более min(N — n + 1, n + 1) шагов, а вычислительная сложность каждого шага не превос-
d k
ходит O(T), где T = N^ nkП P. Общая вычисли-
k = 1 1 = 1 1 тельная сложность составляет
Z p k( Пк + Pk) + 3 T X min (N - n + 1, n + 1, P).
В случае малого числа "пропусков" в полном факторном плане, т.е. при n ~ N, и pk ~ nk = N1/d вычислительная сложность составит O((N — n + + 1)P1 + 1/d) вместо O(P3). Более того, если план эксперимента полный факторный, то N = n и, согласно теореме 1, оптимальный тензор коэффициентов будет найден за одну итерацию. В этом частном случае также можно выписать явную вычислительно эффективную формулу дл
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.