CC BY
30
УДК 519.651
М. Г. Беляев
Московский физико-технический институт (государственный университет) Институт проблем передачи информации им. A.A. Харкевича РАН ООО «ДАТАДВАНС»
Аппроксимация данных, порожденных декартовым
произведением
В работе предложен метод построения аппроксимации для выборок, имеющих структуру декартова произведения произвольного числа многомерных факторов. В подходе используется модель линейного разложении по словарю (некоторому набору) функций, сформированному с помощью тензорного произведения словарей функций меньшей размерности. Для контроля гладкости аппроксимации вводится штраф на изменчивость полученной модели, а коэффициенты разложения по сформированному словарю ищутся оптимальным образом в смысле среднеквадратичной выборочной ошибки. Параметры регуляризации, задающие уровень гладкости модели в каждом из факторов, вычисляются с помощью ошибки скользящего контроля. Специальная структура данных позволяет производить все вычисления со сформированным словарем эффективно. Основанный на предложенном подходе алгоритм успешно используется при решении актуальных прикладных задач.
Ключевые слова: аппроксимация многомерной функции по выборке, тензорное произведение, словари функций.
1. Введение
Работа посвящена задаче аппроксимации, которая заключается в восстановлении некоторой неизвестной зависимости по обучающей выборке — заданному набору пар «точка» -«значение функции в точке» (см. раздел 2). Эта задача широко распространена в инженерных приложениях и является одной из основных составляющих суррогатного моделирования [1]. Для ее решения разработано большое количество подходов [2], которые эффективно работают в широком диапазоне требований прикладной области.
В работе рассматривается специальный случай задачи аппроксимации, в котором множество точек обучающей выборки является декартовым произведением некоторого количества факторов (множеств) произвольной размерности (см. раздел 2.1). Выборки такого рода распространены в практических задачах в силу некоторых особенностей задач инженерного проектирования, обуславливающих выбор именного такого дизайна эксперимента. Большинство методов аппроксимации предназначены для выборок произвольной структуры и не учитывают особенности выборок, порожденных декартовым произведением, что не позволяет в этом случае строить аппроксимацию вычислительно эффективно (обзор методов приведен в разделе 2.2).
Среди многообразия подходов к решению задачи восстановления зависимостей можно выделить классический метод интерполяции двумерных данных с помощью тензорного произведения сплайнов [3], который рассчитан на выборку, порожденную декартовым произведением двух одномерных множеств (двумерная сетка). Использование этого метода, учитывающего структуру данных, позволяет уменьшить вычислительную сложность построения аппроксимации по сравнению с универсальными подходами.
Реальные задачи зачастую обладают существенно более сложной структурой, чем двумерная сетка. Цель работы — предложить метод аппроксимации для выборок, имеющих структуру декартова произведения с несколькими многомерными факторами, сохранив при этом преимущество в меньшей вычислительной сложности построения аппроксимации. В разделе 3 вводятся некоторые операции с многомерными матрицами, необходимые для
дальнейшего изложения. Основные результаты формулируются в разделе 4, где предлагается рассматривать тензорное произведение словарей функций и вводится (в подразделе 4.1) штраф специального вида на изменчивость модели. В подразделе 4.2 представлено утверждение, определяющее условие на существование и единственность решения задачи аппроксимации в рассмотренном классе функций. Предложенный алгоритм построения аппроксимации описан в разделе 4.3, а результаты вычислительных экспериментов приведены в разделе 5.
2. Задача аппроксимации
Рассмотрим задачу аппроксимации функции по конечному набору ее значений в заданных точках. Пусть д — это некоторая непрерывная многомерная функция, д : (X С Д^) ^ КЛу (в дальнейшем будем рассматривать случай (1у = 1, который легко обобщается на случай с1у > 1). Будем называть обучающей выборкой 5 совокупность множества точек £ (мощности Ж) и значений функции д в точках из этого множества:
5 = е £,уг = = {£, £(£)}.
Задача 1 (аппроксимации). По заданной обучающей выборке Б найти функцию f * е Г для некоторого функционального класса Г минимизирующую среднеквадратичную ошибку аппроксимации ф, подсчитанную в точках множества £:
/* = а^шт(ф(/, £, #(£))) = ащшт V](д(х) - /(ж))2.
Отметим, что эта постановка задачи аппроксимации включает в себя интерполяцию как частный случай при ф(/, £,д(£)) = 0.
Задача 1 позволяет найти такую аппроксимацию /* е Г, которая приближает обучающую выборку наилучшим образом. Однако в этом случае может произойти переобуче-
*
£
ние некоторого штрафа, который ограничивает изменчивость модели и позволяет получить более гладкие решения. Например, в качестве такого штрафа может быть выбрана норма
ром гильбертовом пространстве функций (ядерная гребневая регрессия [4]). В этом случае задача аппроксимации также может быть сформулирована как задача минимизации, но с другой целевой функцией.
Задача 2 (аппроксимации со штрафом на изменчивость). Пусть заданы обучающая выборка Б и штрафная функция Рл : Г ^ Д+. Необходимо найти /* е Г для некоторого заданного класса, Г; минимизирующую функцию ошибки со штрафом Д(/, £, <?(£)) = ф(/, £, <?(£)) + Рл(/); а именно:
г = а^шт(ф(/, £, <?(£)) + Рл(/)) = а^шт Т(д(х) - /(ж))2 + Рл(/).
Заметим, что штрафная функция принадлежит некоторому семейству и задается с помощью параметра Л, который определяет вклад штрафного слагаемого в общую целевую функцию. Как правило, Рл зависит от Л линейно. Например, задача ядерной гребневой регрессии сводится к следующей оптимизационной задаче:
/* = а^шт - /(^))2 + Л||/|Ц.
2.1. Структура данных
Опишем более формально структуру данных, которую мы будем рассматривать в рамках этой работы. Пусть заданы некоторые множества точек
ак = {х^ е Хк, гк = 1,...,пк}, Хк е И*, к = 1,...,К.
Определение 1. Будем называть множества ак факторами, а размерность г!к каждого элемента, множества — размерностью фактора ак.
Определим множество точек £ как декартово произведение всех факторов: £ = а1 х а 2 х ■ ■ ■ х а к • Тогда элементы множества £ — это вектора размерности й = ^ (!к (а объем этой выборки равен N = Пк=1 пк):
£ = {х = ,х12, . . . ,хКк{х^к е °к}к=^.
Описанная структура данных встречается в реальных задачах достаточно часто. Например, распространен случай расположения данных на многомерной сетке, поскольку он удобен при планировании эксперимента и обеспечивает равномерное покрытие пространства дизайна точками [1]. В введенных обозначениях такому дизайну соответствует декартово произведение К одномерных факторов.
Также в инженерных приложениях часто встречаются задачи, в которых все переменные можно разбить на две группы: первая из них соответствует условиям эксперимента, а вторая часть определяет конкретное измерение, сделанное в ходе этого эксперимента. Например, при исследовании аэродинамики крыла первой группе переменных может соответствовать его геометрическое описание, второй группе — угол атаки крыла или скорость воздушного потока, а в качестве выходной переменной может выступать подъемная сила, которой обладает крыло. Сложность изменения значений переменных в первой и второй группах может быть существенно разной (в смысле подсчета нового значения выходной характеристики для модифицированного входного вектора). В описанном примере изменение геометрии крыла — это достаточно дорогая операция (особенно при экспериментах в аэродинамической трубе), в то время как скорость потока набегающего воздуха и угол атаки могут быть изменены большое количество раз для каждого крыла. Эта особенность делает целесообразным использование одинаковой, достаточно плотной сетки значений (набора векторов) переменных из второй группы для каждого зафиксированного значения вектора переменных из первой группы. В таком случае выборка данных будет иметь структуру декартова произведения с двумя факторами. В некоторых задачах факторизация имеет более сложную структуру и переменные разбиваются на большее число групп, чем две. Следует заметить, что несимметричная сложность изменения значений переменных из разных групп часто приводит к возникновению анизотропии: переменные из «недорогих» групп переменных (угол атаки) принимают существенно больше различающихся значений, чем переменные из «дорогих» групп (геометрия крыла).
Вернемся к задаче аппроксимации и рассмотрим ее частный случай: пусть множество точек обучающей выборки есть декартово произведение К > 1 факторов, £ = а\ х а2 х ■ ■ ■ х ок- Тогда область определения функции д также будет декартовым произведением: X = Х1 х Х2 х ■ ■ ■ х Хк- Поставим в соответствие задачам аппроксимации 1, 2 аналогичные задачи с ограничением на структуру выборки.
£
изведение: £ = а1 х а2 х ■ ■ ■ х ак-
£
изведение: £ = а1 х а2 х ■ ■ ■ х ак-
2.2. Обзор методов аппроксимации
Для решения задач 1 и 2 существует множество общих подходов, которые могут быть
£
декартово произведение. Например, это могут быть следующие распространенные методы:
• линейная регрессия [2];
• гауссовские процессы (или кригинг), см. [5];
• разложение по словарю нелинейных параметрических функций, см. [6];
• сплайны (для й = 1), см. [3].
Однако универсальность этих методов не позволяет добиться высокой эффективности для частных задач 3 и 4. Сплайны предназначены для решения одномерных задач и по определению не могут быть использованы для размерности й > 1. Линейная регрессия предлагает крайне простую модель, которая может быть востребована при малом количестве данных, но, как правило, выборки, порожденные декартовым произведением, имеют достаточно большой объем, который может доходить до нескольких миллионов точек. Большие значения N затрудняют использование и остальных техник.
Гауссовские процессы имеют вычислительную сложность 0(<!у N3 + йЫ2), что ограни-
1000
гауссовские процессы [7] предназначены для решения этой проблемы и снижают вычислительную сложность до 0(йу N т2 + ^ N т) (где т — это размер некоторой подвыборки из обучающей выборки). Как правило, в качестве т надо брать некоторую не стремящуюся к нулю долю от N, что позволяет расширить область применимости гауссовских процессов лишь до объемов выборки порядка 10 000 точек.
Модели, основанные на разложении по словарю нелинейных параметрических функций (включая искусственные нейронные сети), в основном зависят от объема выборки N линейно (за исключением регрессии на основе радиальных базисных функциях, сложность которой может доходить до 0(Ж3)). Однако вычислительная сложность этих методов существенно зависит от объема словаря р и размерностей данных йу. В зависимости от метода настройки параметров функций словаря сложность может варьироваться от 0(Ир (й+йу)) до 0(Ирйу + Nр2 й2 + р3 с13). Зачастую р выбирают пропорционально N, что приводит к (как минимум) квадратичной зависимости от объема выборки. В противном случае зависимость остается линейной, но множитель перед ней принимает достаточно большие значения. Кроме того, необходимо учесть итеративный характер обучения (настройки параметров) таких моделей (количество итераций может доходить до нескольких тысяч), что еще сильнее увеличивает масштабный множитель.
Помимо высокой вычислительной сложности следует отметить еще одну проблему универсальных методов в рассматриваемом классе задач аппроксимации. Большинство из этих
£
стве дизайна равномерно. Это условие нарушается в случае анизотропных данных (мощности различных факторов существенно отличаются) и построенные универсальными методами модели могут иметь крайне низкую точность в областях, удаленных от точек обучающей
£
формализуется строго, однако этот недостаток универсальных нелинейных методов проявляется в подавляющем большинстве реальных задач.
Существует подход, учитывающий наличие структуры декартова произведения в обучающей выборке данных. Этот метод основан на тензорном произведении сплайнов [3] и предназначен для решения задачи интерполяции данных, заданных на двумерной решетке (частный случай задачи 3 для двух одномерных факторах). Вычислительная сложность этого метода составляет 0(п3 + п2 + N(щ + П2)<!у), где щ и П2 — это количество точек в соответствующих факторах, щ П2 = N. Так как для сплайнов можно явным образом
построить базис, то метод не требует итеративного обучения и константа перед оценкой сложности мала.
В работе будет предложен метод, основанный на тензорном произведении некоторых наборов функций, построенных в каждом из факторов. Такой подход является обобщением тензорного произведения сплайнов и позволяет как строить аппроксимацию для выборок, в которых £ есть декартово произведение общего вида (произведение К многомерных факторов), так и решать не только задачу 3, но и более общую задачу 4 со штрафом.
3. Тензорное произведение и многомерные матрицы
Для изложения предлагаемого подхода к решению задачи аппроксимации необходимо ввести некоторые определения. В этом разделе будет описано понятие тензорного произведения, а также рассмотрено обобщение матричных операций на ^-мерные матрицы.
3.1. Тензорное произведение
Пусть функции ¡1(х1)ш ¡2(х2) определены на Х1 е И^1 и Х2 е И^2 соответственно.
Определение 2. Функцию f ([х1, х2]), определенную на Х1 х Х2, будем называть тензорным произведением функций /1(х1) и /2(х2) и обозначать /1 ® ¡2, если
/ а^2]) = мх1) 12 (х2).
Несложно показать, что тензорное произведение ассоциативно, коммутативно и билинейно. Кроме того, оно обобщается и на другие математические объекты. В частности, нас будет интересовать тензорное произведение матриц, которое носит название произведение Кронекера и обладает как указанными выше, так и другими свойствами [8].
3.2. Многомерные матрицы и операции с ними
Определение 3. Введем, обобщение понятия матрицы. Будем называть У К-мерной матрицей с размерами п1 * п2 * ■ ■ ■ * пк, если
У = . {1к = 1,...,пк}кк=1}.
Введем операцию развертки Х-мерной матрицы вдоль направления к, которое ставит в соответствие многомерной матрице ^матрицу Уразмера П=к П1 * пк- Строки матрицы У— это срез значений тензора У размера 1 *пк при зафиксированных гк+1 = 1,... ,гк = 1, 11 = 1,... ,%к-1 = 1 и изменяющемся = 1,..., В случае обычных матриц (К = 2) операция развертки сводится к транспонированию: У(1) = УТ, и У(2) = У-
Стандартные операции (сложение Х-мерных матриц и умножение на скаляр) вводятся естественным образом (поэлементно), как и скалярное произведение:
&, У) = Е
11,12,.. .,гк
Помимо этих операций введем еще одну: умножение ^-мерной матрицы на матрицу вдоль направления. Пусть В — это некоторая матрица размера пк * п*к. Тогда будем говорить, что К-мерная матрица 2 размеров п1 * ■ ■ ■ * пк-1 * п*к * пк+1 * ■ ■ ■ * пк — это результат умножения У вдоль направления к на матрицу В, и обозначать 2 = У ®кВ, если 2= У(к)В. Отметим, что умножение вдоль отличающихся направлений коммутативно, а вдоль одинаковых может быть сведено к одному умножению:
к к ' 1 I У»к[ВС), к = I.
Умножение многомерных матриц вдоль направления тесно связано с произведением Кронекера. Пусть ^ — это многомерная матрица размера п1 * ■ ■ ■ * пк, а Вк, к = 1,... ,К — матрицы размера Пк * ^.Рассмотрим опер ацию vec, которая ставит в соответствие многомерной матрице вектор-столбец, содержащий все ее элементы. Тогда справедливо
vec (У ®1Вх... ®КВк) = (Bi ® В2 ® ■ ■ ■ ® В к )vec(^). (2)
Сравним вычислительную сложность обеих частей равенства (2). Для упрощения предположим, что все матрицы Вк квадратные и имеют размер Пк * Пк, а N = Ппк- Тогда вычисление левой части требует 0(N ^Пк) операций, в то время как для подсчета правой части необходимо 0(N2) операций. Таким образом, использование умножения вдоль направления позволяет получить заметно меньшую вычислительную сложность.
4. Тензорное произведение словарей функций
Для решения задач аппроксимации 3 и 4 необходимо выбрать класс функций F, в котором будет искаться решение соответствующей оптимизационной задачи. Рассмотрим подход, основанный на тензорном произведении словарей (некоторых наборов функций). Вернемся к методам аппроксимации, перечисленным в разделе 2.2 (сплайны, линейная регрессия, гауссовские процессы, разложение по словарю параметрических функций). Все эти методы аппроксимации объединяет общее свойство: полученную модель можно выписать в виде f(x) = aj^j (х), гДе ВИД функций фj (х) зависит от конкретного метода
рассматриваться как разложение по словарю (ф- (х)}^=1.
Теперь определим класс функций, в котором будем искать решение задач 3 и 4. Введем следующие обозначения для всех к = 1... К:
• a к = (х^ Е Хк УП^=1, Хк С Rdk — фактор декартова произведения;
к
• Дк = (фк И =1 — словарь функций, определенных на X¡¡..
J К J К
Предположим, что словари функций Дк заданы. Сформируем словарь функций, определенных на Xi х ■ ■ ■ х Хк, с помощью тензорного произведения функций из Дк'
Дгепзот = ® Ф|2 , (Зк = 1,.. ., Рк }K=i}.
Пусть линейное пространство функций Ftens0r состоит из всех возможных линейных комбинаций функций из Дгепзж- Выпишем представление некоторой функции f Е F tensor в общем виде с использованием многомерных матриц:
f(x) = Л (х1) ®2^2(x2)... ®кФк(хк), (3)
где фк(х) — это вектор-столбец значений функций из словаря Дк в точке х, а Л — многомерная матрица коэффициентов разложения по этому словарю. Далее мы будем искать решение задач 3 и 4 именно в классе функций Ftens0r-
4.1. Штрафная функция
Значения функции g в точках из £ параметризуются инде ксами i\, i 2,..., i к и могут быть представлены в виде многомерной матрицы У с разм ера ми щ * П2 * ■ ■ ■ * Пк'
gfah ,х22 ,...,хк D= Vhñ ,...,Ík ^ #(£) = У. (4)
Обозначим как Фк матрицу значений функций словаря Дк в точках из ак (размер матриц Фк будет равен Пк * 'Рк) и выпишем функцию ошибки Q, используя выражения (3) и (4):
Q{f, £, <?(£)) = £>(х) - /(х))2 = (У -У, У -У), У = № ... ®кФк. (5)
хет
Задача аппроксимации 3 сводится к минимизации (5) по f Е Ftensor- Рассмотрим более общую задачу со штрафом 4, которая включает задачу 3 как частный случай. Для по-
F
строиться аппроксимация, но и вид штрафующего слагаемого P\(f). В выбоpe P\(f) есть некоторая свобода, поэтому обычно штрафное слагаемое выбирают таким образом, чтобы оно позволяло контролировать изменчивость аппроксимации f и получить вычислительно-эффективные формулы. Предложим подход, удовлетворяющий этим соображениям.
В качестве иллюстрации рассмотрим штраф, ограничивающий изменчивость аппроксимации в Xi — первом факторе декартова произведения. В качестве меры изменчивости будем использовать квадрат выборочной нормы второй производной f по х1, подсчитанной в точках множества £:
P(f ) =
d2f
(dx1)2
í d 2ib1 \'
£ Z T^fe(4(4) ...^ «) . (6) \ií,...,iK jl,...,jK '
Теперь обобщим этот штраф и будем ограничивать изменчивость аппроксимации в произведении некоторого числа факторов Х^1 х ■ ■ ■ х Х^. Введем бинарный вектор в, ненулевые компоненты которого будут соответствовать индексам к\,... ,кя, определяющим вид штрафа. Штраф из примера (6) эквивалентен р(/) = /), где у в не равна нулю только первая компонента: 0\ = 1. Для в с несколькими компонентами, отличными от 0, штраф определяется аналогично (6), но порядок частной производной / будет равен не 2, а 2||0||о (под Ьо нормой понимается количество ненулевых элементов). Перепишем выражение (6) для
рв(/) = (Л Л ®1(01 О + (1 - в 1)Ф?Ф1) ... ®К{вкПк + (1 - Ок)Фк*к)>. (7)
Каждая матрица 0-к (размера рк * Рк) в (7) — это выборочная норма1 вторых частных
производных функций словаря А&. Пусть д2Фк/дх1^дх^п — это матрица значений частных
производных функций словаря Д^ в точках из a к по компонентам вектора хк, тогда
^ ^ ( Q2^k хТ , д2^к
^ / дЧк \Т í дЧк \
к l4~^Adxidxm) \дх'[дхкп) .
¡-к
1=1 т=1
Теперь введем вектор параметров регуляризации Л размера К * 1, с помощью которого будем взвешивать р$ (/) для разных в при формировании окончательного штрафа Рд(/). Каждая компонента Лк будет определять штраф на изменчивость аппроксимации в соответствующем Х^. Пусть В — это множество всех бинарных векторов длины К за исключением нулевого. Тогда определим штрафную функцию следующим образом:
Ра(Я = Е (рв(/) П (Лквк + (1 - 9к))). (8)
вее^ к=1 '
Введение вектор Л дает возможность явно определять относительный вклад ошибки аппроксимации £, <?(£)) и изменчивости Рд(/) в общую функцию ошибки. Кроме того, использование именно такой параметризации позволяет штрафовать изменчивость в различных факторах по-разному, что может быть крайне востребовано в задачах с анизотропией (мощности факторов Пк сильно отличаются для разных к). Например, если устремить Лк к 0, то в штрафной функции Рд(/) устремятся к нулю все коэффициенты перед слагаемыми, включающие частные производные с использованием Хк, и изменчивость в факторе Хк не будет штрафоваться.
2
s
1для легко вычислимых функций (например, сплайнов) О может быть и интегральной нормой
4.2. Вычисление оптимальных коэффициентов разложения
Утверждение 1. Решение задачи 4 в классе Ftens0r со штрафом P\(f) вида (8) существует и единственно, если матрицы Фт Фк + АкО к не вырождены. Кроме того, тензор коэффициентов разложения этого решения по словарю Дгепз0г вычисляется по формуле
X = ^®1Ф1Т(ФТ Ф1 + Ai Qi)-1 ®2Ф2Т(Ф! Ф2 + A2Q2)-1... ®КФкT (Фк Фк + Ак Ок )-1. (9)
Дк
ций Ftensor поиск оптимальной аппроксимации в этом классе сводится к минимизации функции ошибки по коэффициентам разложения. В свою очередь, задача поиска многомерной матрицы оптимальных коэффициентов разложения Л может быть сведена к обычной
Ф
словаря Д^^ог в точках £. Тогда выражение (5) для функции ошибки Q может быть переписано в виде Q = (vec(^) — Фvec(Л)) (vec(^) — Фvec(Л)). Штрафное слагаемое из (8) в матричных операциях будет выглядеть как P = vec(^)TQ vec(^) для некоторой матрицы Q (которая подсчитывается на основе (7)). Тогда задача 4 будет переписана как задача поиска оптимальных коэффициентов Л*:
Л* = argmin (vec(y) - Фvec(Л))T(vec(^) - Фvec(Д)) + vec(^)TQ vec(^) A L
Это стандартная задача гребневой регрессии и ее решение может быть найдено явно:
vec(X) = (ФТФ + fi)-4TvecQ>). (Ю)
Ф
Ф = Ф1 ® ■ ■ ■ ® Фк Учтем структуру Ф и перепишем второй сомножитель в (10) по (2):
Ф^ве(;у) = vec(^ <E^T ... ®КФк). (П)
Выпишем представление штрафной матрицы Q, используя (7):
Q = A1Q1) + (1 - ^^Ф^ ® ■ ■ ■ ® (вк( АкПк) + (1 - дк)ФкФк). (12)
вев
Теперь упростим первый множитель в (10), используя свойства произведения Кронекера:
(Ф^ + О)-1 = (^T Фк) (Ф1 Фк ) + Q)-1 = ((ФTФl) (Фк Фк) +Q)-1.
Обратим внимание на вид (12) штрафной матрицы Q и примем во внимание билинейность произведения Кронекера:
(Ф^ + Q)-1 = (^ФО <В> ■ ■ ■ <В> (ФкФк) + Q)-1 = ((ФTФ1 + A1Q1) <Е> ...
■ ■ ■ <В> (фкФк + АкОк))-1 = ^TФ1 + A1Q1)-1 <g> ■ ■ ■ <g> (фкФк + АкОк)-1. (13)
Теперь объединим (11) с (13) и перепишем выражение для оптимальных коэффициентов (10), преобразовав его в многомерную матрицу:
Л* = У ... ®КФк Ф1 + AiQi)-1... ®К(Фк Фк + АкОк)-1.
Используя свойство (1) умножения многомерной матрицы вдоль направлений, получаем итоговую формулу (9).
Условия на существование и единственность следуют из аналогичных условия для обычной гребневой регрессии. Действительно, матрица представима в виде произведения Кронекера (Ф^^Ф! + AiQi) ® ■ ■ ■ ® (Ф#Фк + АкОк)- Это означает, что Ф'ГФ + О не вырождена тогда и только тогда, когда не вырождены все Ф^? Фк + АкПк- ■
Рассмотрим вопрос выбора вектора регуляризации Л. Применим для его оценки стандартный подход — скользящий контроль (leave-one-out) [2]. Этот метод позволяет оценить ошибку аппроксимации линейной регрессии на тестовых данных, используя только точки обучающей выборки:
N ^ 2
Qloo = ^{уг - tp(Xi)vec(A-i)j
2
г=1
где коэффициенты разложения уее(Д-г) оцениваются по выборке, из которой исключена точка {Хг,уг}. Эта формула может существенно упрощена с помощью матричных преобразований [2]:
N / \ 2
Уг - гр(х^)уес(ДП т_ ,т,ЛтгТ,Тг , 0\-1,т,т
Qioo = V * -Р(х;)уес(-М , L = Ф(ФТФ + Q)-1ФТ (14)
Адаптируем формулу (14) к особенностям нашей задачи, приняв во внимание специальную структуру матрицы регрессоров Ф. Это позволит нам воспользоваться вычислительно-эффективными формулами, основанными на операциях с многомерными матрицами.
Утверждение 2. Ошибка скользящего контроля при использовании (9) может быть подсчитана по формуле
Яш = (у-й Су-50 * (х-х®^^)...))~2); (15)
где Ь^ = Фк (Ф^Фк + Лк1 ФТ; — диагональная, матрица, совпадающая с диа-
гональю Ък, а X — тензор размера п1 * ■ ■ ■ * пк-
Доказательство. Найдем выражение для матрицы Ь из (14) с учетом специфики задачи:
Ь = Ф(ФтФ + 0)~ 1ФТ = Ф1 Фк (ФТФ1 + 1 ® — ® (Ф£ ФК + ЛК1ФТ ®...
■ ■ ■ ® ФТк = (Ф1(ФТФ1 + Л1П1Г1ФТ) ® ■ ■ ■ ® (Фх(ф£Фк + ЛкПк)"1Ф^). (16)
Для подсчета ошибки скользящего контроля необходимы только диагональные элементы матрицы Ь, которые согласно (16) выражаются через матрицы Ь^ из (15): diag(L) = ® ■ ■ ■ ® Используя связь (2) между произведением Кронекера
и умножением многомерных матриц, преобразуем (14) к скалярному произведению многомерных матриц и получаем искомую формулу (15). ■ Вычислительная сложность формулы (15) составляет 0(^ Рк + N), что позволяет использовать эту оценку ошибки аппроксимации на тестовой выборке для выбора вектора регуляризации Л. В силу факторизованной структуры функции Я100 использование метода покоординатного спуска для ее минимизации будет достаточно эффективным.
4.3. Алгоритм построения аппроксимации
В предыдущем разделе предполагалось, что словари Д^ заданы, однако это справедливо только для одномерных факторов, в которых можно явным образом построить базис в пространстве сплайнов. Рассмотрим задачу построения словарей Д^ по заданной обучающей выборке Как уже было указано выше, универсальные методы аппроксимации могут быть представлены как разложение по некоторому словарю параметрических функций. Построение аппроксимации для этих методов сводится к определению оптимальных (в смысле минимизации некоторой функции ошибки) параметров функций, составляющих словарь. Зафиксируем словари Дк,к = I и рассмотрим задачу нахождения словаря Д[.
Функция ошибки Q в этом случае может быть переписана следующим образом:
Щ Р1 0
к
Q(/, £, <?(£)) = £ £ {уп....,гк - Е Е Чъ^эЛК) ...Фк(*кк))'
4 = 1 гк,к=1 31 = 1 ]к,к=1
Щ / Р1 х 2 / х
Е Е(Учк=1 - (4)) = ^асе( (У(1) - Ф^)^ - Ф^)). (17)
%к,к=1 ч = 1 31 = 1
Задача минимизации (17) по функциям {ф3[ У^=1 = Аг (точнее, по их параметрам) эквивалента решению задачи 1 по выборке 5 = {01, УТаким образом, построение словарей А^ сводится к К задачам аппроксимации, которые могут быть решены с помощью универсальных методов.
Сформулируем полный алгоритм построения аппроксимации, предполагая, что обучающая выборка действительно представляет собой декартово произведение с известной структурой.
1) Формируем из значений функции {^^многомерную матрицу У = д (£).
2) Для всех к = 1,...,К строим аппроксимацию, используя обучающую выборку вк = {ок, УРешаем эту задачу с помощью одного из стандартных методов аппроксимации. Выбор метода должен быть основан на параметрах получившейся задачи. Например, для йк = 1 логично использовать сплайны, а для больших Пк при йк > 1 — словарь нелинейных параметрических функций.
3) Находим оптимальный параметр регуляризации Л*, минимизируя (15).
4) Подсчитываем по (9) коэффициенты разложения А* по словарю А^,^ для Л*.
Рассмотрим предложенный алгоритм с точки зрения вычислительной эффективности. Оценим сложность нахождения словаря Ак в каждом факторе. При решении задачи аппроксимации по выборке в к = {о к, У} размерность у будет достаточно велика и равна И/пк, а объем обучающей выборки составит Пк- Если в к-м факторе словарь строится с помощью сплайнов, то сложность будет 0(п3 + пк(И/щ)) = 0(п3 + N Пк)• Заметим, что для повышения вычислительной эффективности можно использовать В-сплайны [3], которые позволят понизить вычислительную сложность до 0(пк + И) за счет ограниченного носителя соответствующих базисных функций. Для гауссовских процессов критически важным является объем выборки, который в силу факторизации существенно сокращается. Вычислительная сложность решения подзадачи аппроксимации в к-м факторе с помощью этого метода будет равна 0((И/пк) = 0(И пк). При использовании разложения по словарю параметрических функций в зависимости от метода обучения сложность составляет от0(пкРк (И/пк)) = 0(Ирк)до0(пкРк (И/пк)+пкр2кй2к +р3к(Рк) = 0(Ирк+пкр2кй2к +р3к(Рк). Следует заметить, что объем словаря А^еЩ80Г равен Р = Л Рк, что позволяет выбирать достаточно малые рк- Таким образом, учет структуры выборки позволяет получить линейную зависимость от объема выборки для всех методов аппроксимации и уменьшить коэффициент перед некоторыми линейными слагаемыми (в частности, множителями с рк для последнего метода).
Сложность вычисления оптимальной Х-мерной матрицы коэффициентов разложения Л по (9) будет равна 0 р1 + N (пк + Рк) Пг=1 Р1 /п1) • Так как для всех методов
Рк/пк ^ 1) т0 вычислительная сложность не будет превосходить 0(^2к=1(пк + N Пк))• При минимизации (15) дополнительная стоимость каждой итерации (по сравнению с вычислением А) составляет 0(^к^=^йП^ 1Р1) и не превосходит к^= 1Пк)•
Для упрощения анализа будем предполагать, что число точек Пк одинаково во всех факторах и равно N1/к, Тогда наиболее тяжелой операцией будет построение регрессии с помощью гауссовских процессов со сложностью 0 (И 1+'2/к). Если же использовать другие
методы аппроксимации, то сложность не превысит О (М 1+1/к) . Кроме того, если в задаче есть одномерные факторы (а в подавляющем большинстве прикладных задач это условие выполняется), то сложность вычисления У и минимизации ошибки скользящего контроля будет уменьшаться за счет использования В-сплайнов.
Таким образом, вычислительная сложность метода относительно объема выборки будет несущественно превышать линейную, а константа перед этой зависимостью будет достаточно мала.
5. Экспериментальные результаты
Рассмотрим ряд прикладных задач из области аэрокосмического проектирования. Типичным примером таких задач является зависимость некоторой аэродинамической характеристики от условий полета и конфигурации объекта. Структуры множеств точек обучающей выборки £ рассматриваемых задач представлены в табл. 1.
Таблица1
Факторизация задач
Задача К d dk N Пк
1 3 3 1 х 1 х 1 780 10 х 26 х 3
2 3 5 1 х 3 х 1 502460 679 х 148 х 5
3 7 7 1 х 1 х 1 х 1 х 1 х 1 х 1 540000 25 х 12 х 4 х 6 х 5 х 3 х 5
Сравним предложенный метод аппроксимации с линейной регрессией (оценка коэффициентов с помощью гребневой регрессии), разложением по словарю параметрических функций (однослойная нейронная сеть, обучаемая алгоритмом Левенберга—Марквардта [9]) и регрессией на гауссовских процессах (пакет (¡РАН. [10], разработанный автором книги [5]). Для оценки ошибки аппроксимации будем использовать стандартный подход: разобьем выборку данных на обучающую и тестовую части. Первая из этих подвыборок будет использоваться для построения модели, а вторая — для определения качества построенной аппроксимации. Заметим, что тестовую выборку необходимо выделять таким образом, чтобы сохранить структуру декартова произведения точек обучающей выборки. В качестве ошибки аппроксимации будем использовать среднеквадратичную ошибку, нормированную на ошибку аппроксимации константой (#£ — мощность множества £):
E(f, Е, <?(£)) =
/Е^Ш - /И)2"
Будем подсчитывать ошибку не только на тестовой выборке, но и на обучающей, а также учтем время обучения модели. Построение аппроксимации с помощью регрессии на гауссовских процессах для задач 2 и 3 не представляется возможным в силу большого объема обучающей выборки (как было отмечено выше, этот метод требует 0(М3) операций и работает разумное время для N порядка тысячи точек). Результаты экспериментов приведены в табл. 2.
Т а б л и ц а 2
Ошибки аппроксимации и время построения модели
лип. регр. нейрон.сети rave. проц. предл. подход
Etr Etest Etr Etest t Etr Etest Etr Etest
1 0.003 0.81 0.65 12.5 0.04 0.38 28.0 0.002 0.42 0.11 0.05 0.12
2 0.08 0.83 0.84 9500 0.048 11.6 - - - 74.5 0.021 0.022
3 0.02 0.14 0.13 3800 0.002 0.003 - - - 2.96 0.0003 0.0009
Заметим, что увеличение ошибки на тестовой выборке но сравнению с обучающей в 12
тропией данных. В силу существенной неравномерности распределения точек но пространству дизайна построенные этими методами модели обладают нерегулярным поведением,
1
12
мя предложенный метод построения аппроксимации позволяет получить качественную ап-
1
срез построенной модели для задачи 2 на рис. 3. Истинная функция /(х) в задаче 3 имеет не очень сложную структуру, что косвенно помогло избежать резкого скачка ошибки на тестовой выборке для искусственной нейронной сети, хотя метод все же уступил в точности предложенному подходу. Также следует отметить время построения моделей для предложенного подхода оно существенно меньше, чем для универсальных нелинейных методов аппроксимации.
Рис. 1. Регрессия на гауссовских процессах: срез Рис. 2. Предложенный подход: срез модели в пи-модели в задаче 1 даче 1
ИАппроксимация 0.3
* Обучающая выборка
Рис. 3. Предложенный подход: срез модели в задаче 2
6. Заключение
Предложенный в работе метод формирования словаря с помощью тензорного произведения позволяет как свести одну задачу аппроксимации к ряду более простых, так и получить богатый класс функций, способный хорошо приближать сложные многомерные зависимости. Изменчивость модели ограничивается с помощью штрафа специального вида, который построен таким образом, чтобы дать возможность контролировать гладкость
модели в каждом факторе независимо от других. Для нахождения оптимальных коэффициентов разложения получены вычислительно-эффективные аналитические формулы.
Низкая вычислительная сложность алгоритма позволяет быстро строить модели даже для крайне больших объемов выборки, а контроль гладкости полученной модели успешно решает проблему анизотропии, характерную для прикладных задач. В качестве основного недостатка можно отметить существенное ограничение на структуру выборки, которое, впрочем, выполняется для многих задач инженерного проектирования. В частности, предложенный подход был успешно применен для решения ряда актуальных задач из аэрокосмической отрасли.
Работа выполнена при поддержке Лаборатории структурных методов анализа данных в предсказательном моделировании, МФТИ, грант правительства РФ дог. 11.G34.31.0073.
Литература
1. Forrester A.I.J., Sobester A., Keane A.J. Engineering design via surrogate modelling: a practical guide. — New York: Wiley, 2008.
2. Tibshirani R., Hastie Т., Friedman J. The elements of statistical learning: data mining, prediction and inference, 2nd edition. — New York: Springer, 2009.
3. de Boor C. A Practical Guide to Splines, 2nd edition. — Berlin: Springer-Verlag, 2001.
4. Scholkopf В., Herbrich R., Smola A. A generalized representer theorem // Computational Learning Theory. - 2001. - V. 2111. - P. 416-426.
5. Rasmussen C.E., Williams C.K.I. Gaussian processes for machine learning. — Cambridge: MIT Press, 2006.
6. Беляев М.Г., Bypnaee E.B., Любгм А.Д. Методика формирования функционального словаря в задаче аппркосимации многомерной зависимости // Сборник докладов конференции «математические методы распознавания образов «ММРО-15». — 2011. — С. 146-149.
7. Бурнаев Е.В., Янович, Ю.А. Разреженные гауссовские процессы // Труды 54-й научной конференции МФТИ «Современные проблемы фундаметальных и прикладных наук». - 2011. - Т. ФУПМ-2. - С. 85-86.
8. Graham A. Kronecker Products and Matrix Calculus: With Applications. — New York: Wiley, 1982.
9. Matlab Neural Network Toolbox.
http: / / www.mathworks.com/products/neural-network/index.html
10. Gaussian processes for machine learning toolbox.
http://www.gaussianprocess.org/gpml/code/matlab/doc/index.html
Поступим в редакцию 28.11.2012.
