CC BY
39
Наилучшие ^-приближения в классах распределений В.П. Пантелеев
Экономический факультет МГТУ, кафедра прикладной математики и естественнонаучных дисциплин
Аннотация. Излагается градиентный метод поиска наилучшего ^-приближения в классах распределений. Геометрическая форма дифференциала целевой функции позволяет аппроксимировать функции, задаваемые, в частности, в форме графика, например, считывающими устройствами. Схема, изложенная здесь для класса нормальных распределений, легко переносится на другие классы распределений.
Abstract. The paper has considered a gradient method of searching the best L1 approximation in classes of distributions. The method is based on the geometrical differentiation and can be applied for approximation of the functions, expressed in a graph form, for instance, by reading machines. The method, proposed here for the class of normal distributions, can readily be extended to other classes of distributions.
1. Введение
Перенос изображений на экран компьютера и преобразование контуров изображений в линии, как их понимает машина, - прекрасная отправная точка для многих новых технических решений. Среди прочего сюда относится проблема наилучшего приближения контурного изображения в том или ином классе геометрических областей, составляющая предмет данной статьи. Для содержательности изложения мы ограничимся численным поиском наилучшего ^-приближения в классах распределений, задаваемых функциями плотности, хотя возможно и параллельное изложение для интегральных функций распределения. Геометрический подход к проблеме обеспечивает широкую свободу при выборе формы задания аппроксимируемого элемента - в нашем случае функции плотности p. Не исключено, что последняя выражена, например, графиком, выведенным на экран компьютера. В этом случае график функции p можно интерпретировать как ломаную, возможно, составленную из большого числа малых звеньев, которая ограничивает сверху полигон, лежащий между нею и осью абсцисс, возможно, несвязный. Ясно, что все иные представления функции p легко сводятся к графическому. В то же время, как выявится ниже, если стартовать от графика функции p, интерпретируя его как ломаную, предлагаемая здесь схема аппроксимации сводит ошибку метода потенциально к нулю, оставляя лишь вычислительную ошибку.
Нетривиальность поставленной задачи в том, что элемент наилучшего ^-приближения ищется в классах распределений, вообще не являющихся линейными многообразиями и потому не допускающих прямого применения линейных методов. К тому же, метрика L1 менее обещающая, чем, например, L2. Для кусочно-непрерывных функций p проблема усложняется требованием более полного учёта содержащейся в ней информации, без упрощений, таких как стягивание всей информации о функции p к её значениям в отдельных точках дискретной модели. В целом решение проблемы видится здесь в составлении численного алгоритма для компьютерного поиска элемента наилучшего ^-приближения кусочно-непрерывной функции p в выбранном классе функций плотности K.
2. Аппроксимация кусочно-линейных функций
Сначала предположим, что аппроксимируемая функция p кусочно-линейная, возможно, с разрывами первого рода и, возможно, заданная графиком, выведенным на экран компьютера в ходе какого-либо эксперимента. В этой интерпретации функция p имеет конечное число локальных экстремумов, и её график вместе с осью абсцисс составляет полигон. Заметим, что уже в требовании кусочной линейности p нет жесткого ограничения, поскольку, например, всякая кривая на экране компьютера в действительности является ломаной. Кроме того, в конце работы будет установлено, что условие кусочной линейности функции p - всего лишь временное ограничение и что результаты, полученные для кусочно-линейных функций, сохраняют силу и для кусочно-непрерывных. Более того, предложенный здесь численный алгоритм ищет наилучшие приближения не только для функций, но, как увидим ниже, и для объектов более широкого класса - геометрических областей. В то же время численный алгоритм, хоть и не существенно, но всё же варьируется в зависимости от выбора класса K аппроксимирующих функций, поэтому для содержательного изложения метода в качестве K выбираем здесь двухпараметрический класс нормальных распределений:
Интеграл от модуля разности функций f и p равен площади S(G'AGt) симметрической разности G'AGt = (Gí\G')u(G'\Gí) областей G
г и G', лежащих между осью абсцисс и, соответственно, графиками
- функций f и р (рис. 1):
/V р
\ \о
а N.
|| f (х, а, Ь) - р(х) | <1х = S(G'Д Gt).
(2)
рис. 1
От областей не требуется здесь топологической связности, поэтому область можно понимать как замыкание открытого множества, возможно, несвязного или многосвязного, а площадь считать продолженной до соответствующей меры Лебега. Минимуму этого интеграла или, что то же самое, минимуму площади S(G'AGt) симметрической разности G'AGt отвечает максимум площади S(G'г\Gt) общей части G'г\Gt областей Gt и G', что для нормированных функций плотности следует из равенства
S(G'AGt) = Б^) + Б^') - 2Б(&^) = 2- 2Б(&^), а для неотрицательной р с конечной мерой Б^') - из равенства
Б^^) = - =1+ Б^') - 2Б(&^).
Наконец, содержательный смысл решаемой задачи сохраняется и в случае бесконечной меры Б^1), если наилучшее приближение f по-прежнему выбирать как отвечающее тахS(G'г\Gt). Излагаемый здесь метод нацелен на максимизацию меры S(G'Г\Gt) независимо от того, нормирована ли функция р или не нормирована, является ли она функцией плотности или таковой не является. Предполагается, что граница области G' состоит из одной или нескольких замкнутых ломаных, которые будем последовательно и единообразно обрабатывать. При столь общих предположениях часть области G', лежащая ниже оси абсцисс, когда такая часть имеется, будет игнорироваться, и нижней границей области G' будет считаться часть оси абсцисс, которая проходит через G'. Выделим из границы дG' одну из граничных ломаных В1В2...Вт с вершинами В^,у), занумерованными для удобства так, что сторона ВтВ1 лежит на оси Ох и совпадает с ней по направлению: Вт(хт,0), В1(х1,0), хт< х1. Для всех остальных сторон В¿В^ этой ломаной, в соответствии с принятой здесь положительной ориентацией, имеем Xj+1 < х■ , Вт+1 = В1. Если же замкнутая граничная ломаная I области G' не имеет на оси абсцисс стороны, то вершины ломаной I занумеровываются в направлении выбранной нами правой ориентации, начиная с любой её вершины В1.
Изменению параметра Ь в равенстве (1) отвечает композиция двух сжатий области Gt, одно вдоль оси абсцисс х \^Х, Х-а = к(х-а), с коэффициентом к и другое вдоль оси ординат, у \^У, У=у/к, c обратным коэффициентом 1/к. Сжатие понимается здесь в широком смысле слова, когда коэффициент к допускается любой положительный, 0 < к < +<х>. Поскольку речь пойдёт о дифференцировании по параметру Ь, преобразования сжатия близки к тождественному, с коэффициентами к вблизи единицы.
Для любых двух нормальных распределений (1), определяемых конкретными значениями а, Ь и а', Ь', переход от одного распределения к другому осуществляется композицией параллельного сдвига и двух указанных выше сжатий Х-а = к(х-а) и У=у/к с коэффициентами к и 1/к. Изменению Да = а'-а параметра а отвечает параллельный сдвиг X = х+Аа области Gt вдоль оси абсцисс, а изменению параметра Ь отвечает композиция Е двух сжатий Х-а = к(х-а) и У=у/к. Проследим связь между коэффициентом к композиции Е и приращением АЬ параметра Ь. Для этого заменим в равенстве (1) координаты х, у =Дх, а, Ь) новыми координатами X, У в соответствии с преобразованием Е.
1 , (X - а)2. V I - ехр(-), отсюда У =-—
кУ = ■
. (X - а)2.
2(Ьк)2 "" ^2к(Ьк) 2(Ьк)2
Видим, что произошла замена значения Ь параметра на значение Ьк. Композиции Е отвечает преобразование Ь \^Ьк параметра Ь и приращение АЬ = Ьк-Ь = (к-1)Ь. При этом коэффициент к через приращение АЬ вычисляется посредством равенства к-1=АЬ/Ь.
Применительно к области О' нас интересует обратное к композиции Е преобразование Е1. Вычисляя дифференциал площади £ в случае преобразования Е1 области О', когда область О( неподвижна, предположим сначала, что граница дО' области О' кусочно-линейная.
3. Вспомогательный материал
Нам понадобятся две простые леммы отвлечённого характера, их доказательства тривиальны, но без внимания к выделенным в них фактам понимание последующего материала, возможно, будет затруднено.
Рассмотрим некоторое пространство Мс мерой т. Принято говорить, что отображение ИМ сохраняет меру, если всякий раз, когда мера т определена для множества X, она определена также и для множества И(Х) и принимает на множествах Xи И(Х) равные значения.
т(И(Х)) = т(Х).
Лемма 1. Пусть в пространстве М с мерой т выделены два множества О1 и О', и пусть обратимое отображение ИМсохраняет меру т, которая принимает на множестве О<пИ-1(О') значение т(О^Ил(О')). Тогда такое же значение мера т принимает и на множестве И(О)^О':
т(И(О)пО') = т(О^Ил(О')). (3)
Действительно, т(ОлИ-1(О')) = т(И(ОлИ-1(О'))) = т(И(О)пО').
Лемма 2. Пусть в пространстве М с мерой т выделены два множества Ог и О', и пусть обратимое отображение ИМсохраняет меру т. Если определено приращение А т = т(О <пИ~1(О'))-т(О'пО) меры т на общей части О'пО1 двух множеств О1 и О', когда О' подвергнуто преобразованию И-1, а О1 неподвижно, то для меры т определено также и приращение Дт = т(И(О)<^О') -т(О'пО^), когда О( подвергнуто преобразованию И, а второе множество О' неподвижно, и эти два приращения равны. Дт =Дт.
Действительно, равенство приращений Дт= Дт мы получаем, если от обеих частей равенства (3) отнять число т(О'пОг).
Композиция Е1 двух сжатий х^Х, Х-а = кл(х-а) и у^У, У = ку вдоль осей декартовой системы координат сохраняет площадь £. Согласно лемме 2, отображениям Е и Е-1 отвечают равные приращения Д£ = Д£ площади £, £(Е(О)пО') -£(О'пО) = £(О,пЕ-1(О')) -£(О'пО). Поэтому, преобразуя композициеи Е область О¡; при условии, что О' неподвижна, мы имеем то же самое приращение, как если бы подвергли обратному преобразованию Е-1 область О' при неподвижной О.
Области, заметаемой направленным отрезком (или его частью) при малом сдвиге или сжатии, приписываем правую ориентацию, если при движении отрезка она перешла с его правой стороны на левую. Если же она перешла с его левой стороны на правую, припишем ей левую ориентацию, а её площади - знак минус. Такую площадь, принимающую положительные и отрицательные значения, называем относительной. При сжатии плоскости направленный отрезок может сдвинуться одной своей
частью вправо, а другой - влево (рис. 2). В этом случае заметаемая им область распадается в заметаемую вправо и влево. Их площадям приписывается, соответственно, знак плюс и минус. При этом общая относительная площадь этих двух частей области, заметаемых вправо и влево, определяется как алгебраическая сумма их относительных площадей. Относительная площадь области, которую заметает ломаная, также определяется алгебраической суммой относительных площадей, составленной для областей, заметаемых всеми её звеньями.
4. Теорема. Площадь £ общей части подвижной области О1 и неподвижной многоугольной области О' является непрерывно дифференцируемой функцией переменных а и Ь всюду в полуплоскости - да < а< + да, 0< Ь<+ да. Если при этом Р^ц, Р202,..., Р/й,... - перечень всех максимальных по включению дуг кривой / нормального распределения (1), лежащих в области О', Р/(хр/,ур), й(хфуц/), х¥ < хр/, то полный дифференциал функции £ переменных а и Ь существует и равен
ё£ = ёаЪ/ (уч/-ур/) + (<ЛЬ/Ь)Ъ/{у91х9/-а)-ур/(хр/- а)}. (4)
Сдвиг абсциссы x к её образу X при обратном отображении F"1 посредством сжатий X-a = k-1(x-a) и Y = ky определяется равенством X-x = (k-1 -1)(x-a). На стороне BkBk+1 области G' могут находиться один или два отрезка, целиком лежащие в области Gt и которые нельзя продолжить, не выводя их за пределы области Gt. Обозначим через CjDj, Cj(xcj,ycJ), Dj(xdj,ydj) один из таких отрезков, равнонаправленный с BkBk+1. Относительная площадь, заметаемая отрезком CjDj при сжатии X-a = k-1(x-a), вычисляется посредством равенства Ax)S = 2~l(XcJ—xCj +Xdj-xdj)(ydj-ycJ) = 2~l(k~1 - 1)(x(j + xljj-2a)(ydj-y(j), a заметаемая всеми отрезками CjDj, лежащими в Gt и составленными для всех сторон BkBk+1, в которых такие отрезки CjDj имеются, -посредством равенства
AxS = E,AxjS = 2~\k~l - 1)Zj (xcj + xdj-2a)(ydj-ycj).
При k^-1 площадь AxS c точностью до бесконечно малой высшего порядка равна приращению площади S общей части областей Gt и G', которое получается в результате сжатия X-a = k-1(x-a) области G' при неподвижной Gt. При последующем сжатии y \^Y, Y = ky сдвиг ординаты y к её образу Y равен Y-y =(k-l)y. Поэтому относительная площадь AjS, которая заметается отрезком CjDj при сжатии Y = ky, равна
AyjS = 2~\Ycj-ycj+Ydj-ydj)(xcj-xdj) = 2~\k -V)(ycj+y^)(xcj-xj), a площадь AyS, заметаемая всеми отрезками CjDj, равна
AyS = 2~\k- YjL](ycj+ydj)(xcj-xdj).
При k^-l площадь AyS с точностью до бесконечно малой высшего порядка равна приращению площади S, которое она получает, когда область G' подвергается сжатию Y = ky, а область Gt неподвижна. Суммарное приращение AbS площади S в результате композиции F1 двух сжатий области G', X-a = k'l(x-a) и Y = ky при неподвижной Gt с точностью до бесконечно малой высшего порядка равно сумме приращений AxS и AyS.
AjS = 2-l(k-l-\)Y.j(xCj + xdj-2a)(ydj-yCj) + Tl(k-1jL ^ +ydj)(xcj-xdj).
Мы составили приращение относительной площади S общей части двух областей Gt и G', когда область G' подвергается преобразованию F1, а область Gt неподвижна. Согласно лемме 2, это же приращение AbS получает площадь S, когда преобразованию F, X-a = k(x-a), Y = k-1y подвергается область Gt при неподвижной G'. Посредством F значение параметра b преобразуется в kb, что отвечает приращению Ab = kb-b = (k-1)b параметра b. Отсюда, k-1 = Ab/b. При Ab ^ 0 или, что то же самое, при k^-1 бесконечно малая (1-k)/k эквивалентна -Ab/b. Учитывая это, составляем частный дифференциал функции S по переменной b для случая, когда посредством F преобразуется область Gt, a G' остаётся неподвижной.
dbS = 2~l(db/b)'L j{xcj + xdj-2a)(ycj-ydj) + (ycj +ydj)(xcj-xdj)}. После приведения подобных членов имеем
dbS = (db/b)'L j{xcj(xcj-a) -ydj(ydj-a)}. (5)
Обратим теперь внимание на то, что суммирование в правой части равенства (5) подчинено правилу сложения векторов. Если направленные отрезки CjDj, Cj+1Dj+1,..., CkDk, лежащие в области Gt, составляют ломаную le. dG', максимальную по включению и направленную согласно ориентации границы dG', то в равенстве (5) при суммировании взаимно уничтожаются все слагаемые, отвечающие внутренним вершинам ломаной l, остаются лишь слагаемые ydk(xdk-a) -ycj(x(j-a), которые отвечают концам замыкающего отрезка CjDk ломаной l - точкам входа границы dG' в область Gt и, соответственно, выхода из неё. В равенстве (4) этим точкам отвечают обозначения Pl.1 и Ql. Поэтому по завершению группировки слагаемых в правых частях равенства (5) останутся лишь координаты точек Ql(xql,yql) входа границы dG' в область Gt и точек P(xpl,ypl) выхода границы dG' из области Gt. При этом точкам выхода Pl в равенстве (5) отвечают слагаемые ypl(xpl-a), а точкам входа Ql - слагаемые -yql(xql-a), наделенные знаком минус. Ввиду имеющейся здесь двойственности происходит замещение ролей точек входа и выхода: точки Ql входа границы dG' в область Gt - это в то же время точки выхода границы dGt из области G', а точки Pl выхода границы dG' из Gt совпадают с точками входа границы dGt в область G'. В итоге равенство (5) принимает следующий вид
dbS = (db/b) Zl{xql(xql-a) -ypl(ypl-a)}, (6)
где P(xpl,ypl) и Ql(xql,yql), соответственно, точки входа границы dGt в область G' и выхода этой границы из G'. Частный дифференциал dbS равенства (6), как видим, соответствует такому же в равенстве (4). Для завершения доказательства теоремы остаётся выявить форму частного дифференциала daS по переменной a.
Приращению Да параметра а отвечает параллельный сдвиг области О1 вдоль оси абсцисс. Процедура составления дифференциала йа£ площади £ по переменной а здесь довольно проста, поэтому сразу выписываем конечный результат. Частный дифференциал йа£ при сдвиге области О1 относительно неподвижной О' равен
= йаЪ/(уч1 - ур/),
где Р(хр/,ур) и й(хд/,уц) - точки, соответственно, входа верхней границы дОг в область О' и выхода её из неё, а йа - приращение параметра а функции / Полный дифференциал функции £ по переменным а и Ь для подвижной области О1 и неподвижной О' равен сумме частных дифференциалов й£ = йа£+йЬ£. Складывая их, получаем требуемое равенство (4). На этом доказательство теоремы завершается.
5. Численная оптимизация функции S
Выписывая частные производные и градиент
д£/да = Е/ (уф- ур/), д£/дЬ = Ь_1Е/{уч/(хч/-а) -ур/(хр/-а)},
g = У£ = £ ¿у9/ - ур/) + \ЬЛ Е/{у (хф - а) -ур (хр/ -а)} = (7)
= ^{г(уд/ - ур/) + Ь \уд^хд/-а) -ур/(хр/-а)]}, замечаем, что все эти функции определены и непрерывны всюду в полуплоскости -да< а<+да, 0<Ь<+да. При этом для их составления достаточно знать лишь точки пересечения верхних границ областей О1 и О' и разбиение множества М этих точек на два подмножества М, и М0, соответственно, точек входа границы дО1 в область О' и точек её выхода из О'. В равенствах (4) и (7) нетрудно заметить альтернированные суммы. Пусть Я1Я2.Я21_1Я21 - кортеж всех точек Я(хгЬуг/) пересечения верхних границ областей О1 и О', в котором нечётные места занимают точки Я2к-1 входа границы дО1 в область О', а чётные места - точки Я2к выхода границы дО1 из О', х1 >х2 >.. >ха. Тогда дифференциал й£ и равенства (7) принимают вид альтернированных сумм
й£ = {£/(-1)4/ }йа + {Ь"1Е ,(-1)/уг1(хИ-а)}йЬ =
= Е/(-1) )уг/йа + Ь'1уг/(хг/-а)йЬ}, д£/да = Е/(-1)/уг/, д£/дЬ = Ь-1Е/(-1)/уг/(хг/-а), g = 1^/(-1)/уи+]Ь_1Е /(-1)/уг/(хг/- а) = = 2/(-1)/{угД + Ь'1уг/(хИ-а)\}. Если градиент g = У£ не равен нулю, то всякий малый шаг Дг = (Да, ДЬ), Дг = =А{1Е/(-1)/уг/+/(-1)/уг/(хг/-а)}
в направлении градиента, ведёт нас к большему значению функции £. Процесс последовательных переходов к большему значению функции £, направленный на отыскание максимума функции £, организуется посредством итераций гп+1= г„+к^„ в пространстве всех состояний г(а, Ь), определяющих взаимное положение областей О1 и О'. Вектором г(а, Ь) фиксируются, например, положения подвижной области О1 при условии, что другая область О' неподвижна. Здесь gn - градиент, вычисленный для п-то состояния, и кп - регулятор шага, в качестве которого в простейшем случае используется положительная константа к. Итерация состоит в преобразовании параметров а и Ь области О1 по следующему правилу
ап+1 = ап + кпд£/да, Ьп+1 = Ьп + кпд£/дЬ,
где частные производные функции £ вычисляются согласно формулам (7). Во избежание сбоев в начале выполнения программы, регулятор кп ограничивается сверху, например, неравенством кп\д£/дЬ |<Ь2-4 или, что то же самое, числом
Ь2'А\д£/дЬ\л = Ь22-4\Е/(-1)/уг/(хг/-а)\-1 > кп.
6. Заключение
Поскольку дифференциал й£ не зависит от того, как ведёт себя граница дО' в промежутках между точками её пересечения с границей дО, ничто не мешает применять рассмотренный выше метод оптимизации не только для кусочно-линейных функций плотности р, но и для кусочно-непрерывных. Между тем, снимая требование кусочной линейности, необходимо сохранить содержащееся в нём для данного класса распределений требование регулярности, которое до сих пор выполнялось автоматически и которое предполагает, что множество дО'Г\дО1 точек пересечения границ на какой-либо из этих границ имеет меру нуль. Если в качестве регулятора выбирается константа кп = к, её надо взять несколько
меньше, но того же порядка, что и указанная выше верхняя граница, поскольку занижение её может замедлить работу компьютерной программы.
Для непрерывных кусков функции плотности p, регулярно пересекающихся с кривыми семейства (1), верны те же выводы, что и для кусочно-линейных функций p, и для отыскания наилучшего приближения пригоден описанный выше алгоритм. Заметим, что лишь в случае регулярности функция S заведомо имеет непрерывный градиент, лежащий в основе алгоритма.
Тестирование компьютерной программы может осуществляться посредством вывода на экран дисплея конечного положения областей Gt и G' или, что то же самое, графиков функций f и p, для визуальной оценки конечного относительного положения областей Gt и G' на предмет оптимальности, равно как и посредством вывода и наблюдения в пошаговом режиме параллельных сдвигов da и сжатий X-a = k(x-a), Y=y/k вдоль осей координат, производимых над подвижной областью. В качестве последней ничто не мешает выбрать любую из областей Gt, и G', но в вычислительном отношении предпочтительнее та, для которой более скоро определяется её новое положение после очередной итерации.
Коэффициент k легко определяется в силу зависимости между ним и параметром b, k = 1+АЬ/Ь, что полезно, например, для вывода на экран обновлённого графика функции f и что естественно осуществлять по завершению одной или нескольких итераций. Накапливающиеся в ходе вычислений ошибки устраняются более точным построением графика функции f по завершению некоего фиксированного числа итераций. Пошаговая оценка правильности выполненных операций может осуществляться также посредством вывода на экран изменения функции p в результате сжатия или сдвига при выполнении очередной итерации.
Данная статья написана в развитие сообщения (Panteleev, 2001), сделанного автором на Третьей Московской международной конференции по исследованию операций (ORM2001).
Литература
Panteleev V. P. A method of computing the best L1 approximations in classes of distributions. The Third Moscow International Conference on OperationsResearch (ORM2001).Abstracts. Moscow, p.92-93, 2001.
