CC BY
18
УДК 519.85 Вестник СПбГУ. Математика. Механика. Астрономия. 2019. Т. 6 (64). Вып. 1
МБС 15А39
Особенности применения метода проектирования в задачах деконволюции
И. Ш. Латыпов
Санкт-Петербургский научно-исследовательский центр экологической безопасности Российской академии наук (НИЦЭБ РАН),
Российская Федерация, 197110, Санкт-Петербург, ул. Корпусная, 18
Для цитирования: Латыпов И. Ш. Особенности применения метода проектирования в задачах деконволюции // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2019. Т. 6(64). Вып. 1. С. 81-87. https://doi.org/10.21638/11701 /spbu01.2019.105
Популярность метода последовательного проектирования (МПП) обусловлена простотой его реализации и экономностью использования памяти. Основная идея метода заключается в том, что выпуклое множество представляется в виде пересечения, конечного или бесконечного, множества простых выпуклых (элементарных) множеств, далее выполняется проектирование на внешние к текущей точке множества. Проектирование на эти простые множества строится очень просто, так как обычно это полупространства. Доказано, что итерационный процесс последовательного проектирования сходится, причем разработаны его модификации, обеспечивающие конечную сходимость. В рамках МПП на каждой итерации решается три подзадачи — найти элементарное множество для проектирования, определить направление и вычислить длину шага в этом направлении. В нашей работе мы предлагаем несколько простых утверждений, которые позволяют объединить эти три задачи и ускорить сходимость метода для решения специального класса задач — обращения свертки (деконволюции).
Ключевые слова: выпуклое программирование, проекция на множество, метод последовательного проектирования, обращение свертки, деконволюции.
1. Введение. Метод проектирования на выпуклое множество предназначен для поиска точки выпуклого множества, заданного как пересечение семейства выпуклых множеств, отличается простотой реализации и не требует дополнительной памяти [1]. В общем случае он достаточно подробно исследован: доказана сходимость и конечность [2, 4, 5], доказано также [6], что метод позволяет получить некоторое приближенное решение в тех случаях, когда множество допустимых решений пусто. Однако на практике метод используется довольно редко — он медленно сходится. Целью нашего исследования будет ускорение МПП за счет использования специфики задачи обращения свертки (деконволюции). Напомним, что сверткой двух одномерных функций ) и ф(-) называется функция и(х) = ^(Ь)ф(х — или в случае, когда функции /(■) и Н(-) заданы таблично, ^ = / * Н = ^к (функция Н(-) —четная). Дискретная задача деконволюции заключается в поиске неизвестной дискретной функции в = = 1,...,п: у = в * Н, где Н — известная
дискретная функция (передаточная функция).
Задачи деконволюции в настоящее время привлекают к себе пристальное внимание из-за все более расширяющегося круга приложений. Это побуждает разрабатывать новые методы решения таких задач и проводить сравнительный анализ
(¡5 Санкт-Петербургский государственный университет, 2019
этих методов. Общее количество различных алгоритмов деконволюции насчитывает несколько десятков, начиная (исторически) с фильтрации Винера и заканчивая различными экзотическими приемами. Попытавшись классифицировать современные методы деконволюции, мы обнаружим для начала слепую (blind) и обычную (non-blind) деконволюции. Слепая деконволюция представляет собой задачу обращения свертки с неизвестной передаточной функцией. Ясно, что ее решение потребует значительно больше и входной, и априорной информации. Классическая (nonblind) деконволюция может быть выполнена различными способами, важный класс среди которых образуют методы фильтрации, восходящие к работам Н. Винера. Они очень быстрые, но хорошо работают со специальными передаточными функциями. При этом эти методы требуют априорной информации о видах шума, влияющего на измеренную величину. То же самое можно сказать о байесовской деконволюции, различных методах факторизации и т. д. Сравнительный анализ относительно современных методов деконволюции в приложении к обработке изображений можно найти в [7]. Важным кажется, однако, то, что в большинстве «нефильтрационных» методов так или иначе необходимо решать систему уравнений или неравенств большого размера. Поэтому важнейшим инструментом при решении задач, связанных со сверткой, являются быстрое преобразование Фурье и итерационные методы решения систем линейных уравнений и неравенств. Вообще, задача деконволюции относится к классу некорректных задач [8], для решения которых принято использовать методы регуляризации. Одним из таких методов является сведение задачи декон-волюции к задаче линейного программирования, а именно к поиску такого набора {sj} и такого е, что е ^ min при \yi — к skh\i-k\\ ^ е для всех i. Для решения этой задачи линейного программирования мы модифицируем уже существующий алгоритм МПП.
2. Описание метода последовательного проектирования на выпуклое множество. Будем исследовать модификацию предложенного в [3] метода отыскания точки выпуклого множества вещественного гильбертова пространства H. Нам потребуются следующие определения из [2].
Определение 1. Вектор p из H называется опорным к множеству Q в точке x из H, если (p, y — x) ^ 0, Vy G Q.
Множество всех опорных векторов к множеству Q в точке x называется конусом опорных векторов и обозначается A(x,Q).
Определение 2. Неотрицательное число S называется шагом из x, сохраняющим опорность вектора p из H (||p|| = 1) для множества Q в точке x из H, если Vy G Q выполняется (p, y — x) ^ —S.
Определение 3. Толщиной множества Q называется конечный или бесконечный диаметр максимального шара, содержащегося во множестве Q.
В [3] для отыскания точки выпуклого множества положительной толщины предложен итерационный алгоритм метода опорных векторов с составным шагом
xk+1 = xk — Хк pk, (1)
где xo —произвольная точка из H, pk G A(xk, Q), шаговый множитель выбирается как Хк = YkSk + ек, где Sk — шаг, сохраняющий опорность pk для Q и ек — заглубление, Yk — число из интервала [0,2]. Пошаговая схема метода дана в алгоритме 1.
Алгоритм 1
0. Пусть y — произвольная точка из H, р > 0, k =1, е1 = 0.5р.
1. Положим х(0) = у, а0 = 0.
2. Вычисляем х„+1 = х„ - Л„р^, ап+1 = ап + Л„((2 - 7„)^„ + ек), где р^ = Рп/УРпУ, Лп = 7п^п +
3. Если хп+1 € Q, то останов.
4. Если ап+1 > (р — £к)2, положим ек+1 = О.Бе^, к = к +1 и перейдем к п. 1.
Выбор этой модификации алгоритма обусловлен тем, что она корректно работает для случая, когда толщина множества Q неизвестна. Для алгоритма 1 в [2, 3] доказана конечная сходимость, если внутренность $(хо,р)р| Q не пуста, а также получены оценки необходимого числа итераций. Необходимо отметить, что в вычислительной схеме алгоритма 1 метод выбора опорного вектора рк и величины шага 5к не обсуждается, тогда как вычислительная сложность зависит от этого выбора весьма существенно. В ранних вариантах МПП (см., например, [1]), в качестве выпуклого множества рассматривалось многогранное множество, описываемое конечным числом линейных неравенств. В этом случае простейшим случаем выбора опорного вектора рк являлся выбор нормального вектора любого нарушенного неравенства, а шаг, сохраняющий опорность, 5к вычислялся как расстояние от точки хк до отвечающей нарушенному неравенству гиперплоскости. Если Q находится в Кп, то вычисление 5к требует п арифметических операций, ровно столько же при вычислении по формуле (1), и общая вычислительная сложность зависит от числа неравенств, описывающих Q, и количества итераций. Вычислительная сложность в более общем случае произвольного выпуклого Q зависит также от сложности построения гиперплоскости, касающейся Q, однако, вообще говоря, обычно не сильно отличается от случая многогранника. Именно тот факт, что основной шаг итерации требует вычислительных затрат, пропорциональных только размерности пространства Кп, делает метод МПП таким привлекательным. С другой стороны, в произвольной точке пространства множество опорных векторов образует выпуклый конус и можно рассмотреть задачу наилучшего выбора вектора рк из А(хк, Q). Поиск оптимального опорного вектора вычислительно может оказаться сравним по сложности с исходной задачей, поэтому сосредоточимся на методах, не требующих дополнительных вычислений.
3. Оценки шага для взаимно-перпендикулярных опорных векторов.
Рассмотрим два утверждения.
Предложение 1. Пусть рг € Д(х, Q), Уг € [1, т], 5г — шаги из х, сохраняющие опорность соответствующих рг, и (рг, р-) = 0 для любых г = ]. Пусть р^ = 5грг, тогда величина 5 = ||р±|| будет шагом из х, сохраняющим опорность вектора р±/||р±|| для Q.
Доказательство. Для каждой пары {рг, 5г} справедливо, что (рг, у — х) ^ —5г при Уу € Q.
Умножив каждое неравенство на 5г и просуммировав полученную систему неравенств по г, получим ^^ №р», у — х) ^ — ^^ 52 ^ 0 Уу € Q. Так как
11р±11 = то (р±/||р±||,у-х) < = Vy € д.
Предложение 2. Пусть рг € Д(х, Q), ||рг|| = 1, Уг € [1,т], 5г —шаги из х, сохраняющие опорность соответствующих рг, и найдется такое ^ > 0, что |(рг, р-)| < 7 для любых г = ] и (т — 1)7 < 1. Пусть р^ = ^г 5грг, тогда вектор pz будет опорным, а 5 — шаг, сохраняющий опорность, будет удовлетворять соотношению 5 = 0.5||р^||.
Доказательство. Рассуждая аналогично предыдущему, получаем
т
(р^, у - х) < 52 < о Уу € д.
¿2
г=1
Оценим норму вектора р^:
Р^П = V/(Pz,PzJ
\
(р^ рз )
¿=1 3=1
\
52 + 55 (р^ рз ).
г=1 г=3
Так как ^2»=з 5»5з(р», Рз) < !2»=з 5г537 И 5»5з + 5з5» < 52 + 52, то Т,г=3 5»53 ^ Рз ) < (т - ^Е 1=1 52 , а также Е »=3 5»53 ^ РЗ ) > -Т.г=3 5»53 7 = -7Е 1=3 5»5з > -(т - 1)7 Ет=1 52, следовательно,
(1 - (т - 1)7) £52 < |^||2 < (1 + (т - 1)7) £52,
»=1 »=1
т
о < ||Р^|2 < 52,
г=1
и из (pz,y-x) < —0.5||р^||2 следует < —0.5||1>^|| = -5.
4. Модификация алгоритма для обращения свертки. Рассмотрим задачу деконволюции в конечномерном пространстве Дп. Будем полагать теперь, что множество д задано системой линейных неравенств
|(Р», X) - Уг\ < §, (2)
причем величина § > 0 здесь задает необходимую точность решения, а г € [1 : п]. Отметим, что п-мерные векторы р» теперь являются дискретизациями сдвига одной и той же передаточной функции Н( ): р» = (р1,р2..р") и рк = Н(г - к). Кроме того, будем полагать, что (2) нормировано таким образом, что все р» имеют единичную длину. Тогда вектор р» является опорным, если (р», х) - у» > § с шагом 5» = (р», х) -у» - §, или, если (р», х) - у» < -§, вектор -р» опорный с шагом 5» = - (р», х) + у» - §. Алгоритм 2
0. Пусть хо —произвольная точка из Дп и т ^ 1. Вычислим вектор невязок г системы: г = у - хо * Ь. Положим ео = тахг», 1 < г ^ п и к = 0, а о = 0, а = 0.5.
1. Вычисляем вектор невязок г = у - х^ * Ь. Если \г»\ < § У г, то задача решена, останов.
2. Если ек < §, то решения удовлетворяющего (2), не существует, останов.
3. Вычислим вектор ш шагов, сохраняющих опорность:
0, если г» < ек,
0, если 3] < г : (р»,р3) < т и Ш3 > ей, г», в остальных случаях.
В соответствии с утверждением 2, общее количество ненулевых ш» не должно превосходить (1 - т )/т. Если ш» < ек У г, то переход на 8.
4. Вычислим вектор итерационного шага Дх^ = ш * h как результат свертки вектора ш с вектором передаточной функции h.
5. Вычислим вектор Xk+i = Xk + Дх^.
6. Вычислим CTfc+i = + ||Дхй||2.
7. Если ak+1 < (е0 — ek)2, ek+1 = ek, к = к +1, совершаем переход на 1.
8. Вычислим efc+i = аек, ffk+i =0, к = к +1, совершаем переход на 1.
Нетрудно видеть, что алгоритм 2 является реализацией алгоритма 1 и отличается конкретным способом выбора опорных векторов: опорный вектор является сверткой вектора ш длин попарно перпендикулярных векторов с вектором передаточной функции h, а использованные в алгоритме 1 параметры длины шага равны Yk = 1 и Xk = 5k — ek. Поскольку в этом случае длина опорного вектора не превосходит удвоенного шага, сохраняющего опорность, то алгоритм 2 является частным случаем алгоритма 1 и сходится за конечное число шагов.
Замечание 1. Выбор величины т зависит от особенностей решаемой задачи. Выбор слишком малого т может привести к тому, что все направления будут считаться различными; слишком большая величина т, наоборот, не позволит использовать большого числа опорных векторов из-за ограничения на их количество.
Замечание 2. Вычисление шагов ш, сохраняющих опорность, в шаге 3 алгоритма может выглядеть достаточно трудоемкой операцией, однако на практике величина (pi,pj) убывает достаточно быстро с ростом \i — j|, поэтому в большинстве случаев для построения ш достаточно однократного просмотра вектора невязок.
Последовательность шагов 1-8 будем называть большой итерацией, а шаги 1-7 — малой итерацией. Малые итерации алгоритма являются шагами проектирования при постоянном значении ek, а большие итерации — элементами перебора последовательности вложенных множеств \(pi,x) — yi\ ^ ek. Отметим, что в задачах деконволюции вычисление свертки естественно выполнять с использованием быстрого преобразования Фурье, тогда вычислительная сложность малой итерации — O(n ln n).
Увеличение производительности в алгоритме 2 достигается за счет сокращения числа малых итераций. Так как в задачах обращения свертки количество «почти перпендикулярных» векторов ограничений обычно велико, удается, в соответствии с утверждениями (1) и (2), шаги в направлении таких векторов проводить единовременно. Общее число малых итераций уменьшается за счет увеличения шага, сохраняющего опорность, максимум в m = (1— т)/т раз. В вычислительных экспериментах при решении задачи деконволюции изображений при n ~ 108, удавалось добиться уменьшения времени выполнения программы в сотни раз. Следует отметить, что такие алгоритмы деконволюции как метод Ван Циттерта или Голда работают все еще быстрее, но они не гарантируют сходимости и не имеют внятного критерия останова. Кроме того, эти методы не обеспечивают необходимой точности решения [9]. Замечательную скорость демонстрирует Винеровская деконволюция, однако ее использование накладывает очень жесткие ограничения на передаточную функцию, и в широком классе практических задач результат оказывается совершенно неудовлетворительным по точности. Наконец, в задачах деконволюции не дают эффекта обычные способы ускорения МПП, основанные на блочном представлении матрицы ограничений, и приемах, основанных на таком представлении. Этот факт не позволяет использовать такой мощный аппарат, как быстрое преобразование Фурье, из-за чего даже прямое вычисление свертки требует недопустимо много времени.
Все выше сказанное не означает, что МПП в новой редакции получает безусловное лидерство при решении задач деконволюции — те, кому хватает точности метода Ван Циттерта, по-прежнему предпочтут его — просто предлагаемая модификация алгоритма придает ему легальную конкурентоспособность.
Что касается выбора параметра а, управляющего скоростью сходимости больших итераций, вычислительный эксперимент не подтвердил оптимальности указанного в [3] значения 0.5. В задачах большой размерности, по нашему опыту, время выполнения от параметра а практически не зависит.
5. Заключение. Метод последовательного проектирования весьма привлекателен для решения прикладных задач. Однако до настоящего времени его использование сдерживается медленной скоростью сходимости. В данной статье поведение метода исследовано для случая, часто встречающегося в задачах деконволюции, а именно, явного наличия большого количества перпендикулярных (или почти перпендикулярных) направлений, ведущих в сторону от заданного выпуклого множества. Доказано, что в подобных случаях эффективным вариантом выбора направления проектирования является суммирование всех таких векторов. Хотя теоретическая скорость сходимости остается прежней, в практических расчетах удается добиться ускорения работы метода в десятки и сотни раз.
Литература
1. Брэгман Л. М. Нахождение общей точки выпуклых множеств методом последовательного проектирования // ДАН СССР. 1965. Т. 162, №3. С. 487-490.
2. Фа,зы,лов В. Р. Один общий метод отыскания точки выпуклого множества // Изв. ВУЗов. Математика, 1983. №6. С. 43-51.
3. Куликов А.Н., Фазылов В. Р. О регулировке заглубления в методе опорных векторов с составным шагом // Исследования по прикладной математике. Казань: Изд-во Казан. ун-та, 1992. Вып. 18. С. 79-87.
4. Юла Д. Математическая теория восстановления изображений методом выпуклых проекций //В кн.: Реконструкция изображений / под ред. Г. Старк. М.: Мир, 1992.
5. Marks R. J. II Alternating Projections onto Convex Sets // In: Deconvolution of Images and Spectra / ed. P. A. Jensson. NY: Academic Press, 1997.
6. Еремин И. И., Попов Л. Д. Замкнутые фейеровские циклы для несовместных систем выпуклых неравенств // Изв. вузов. Матем. 2008. №1. С. 11-19.
7. Schneider D., van Ekeris T., Jacobsmuehlen J., Gross S. On benchmarking non-blind deconvolution algorithms: a sample driven comparison of image de-blurring methods for automated visual inspection systems // In: Instrumentation and Measurement Technology Conference (I2MTC), 2013 IEEE International. 2013. P. 1646-1651.
8. Тихонов А.Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979.
9. Горный В. И., Кислицкий М.И., Латыпов И.Ш. Оценка эффективности алгоритмов синтезирования апертуры сканирующего радиометра // Современные проблемы дистанционного зондирования Земли из космоса. 2010. Т. 7, №2. С. 14-25.
Статья поступила в редакцию 9 июня 2018 г.;
после доработки 6 сентября 2018 г.; рекомендована в печать 27 сентября 2018 г.
Контактная информация:
Латыпов Искандер Шамильевич — канд. физ.-мат. наук, доц.; [email protected]
The specifics of method of consecutive projections in deconvolution problems
I. Sh. Latypov
Institution of Russian Academy of Sciences
Saint-Petersburg Scientific-Research Centre for Ecological Safety RAS (SRCES RAS) ul. Korpusnaya, 18, St. Petersburg, 197110, Russian Federation
For citation: Latypov I. Sh. The specifics of method of consecutive projections in deconvolution problems. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2019, vol. 6(64), issue 1, pp. 81-87. https://doi.org/10.21638/11701/spbu01.2019.105 (In Russian)
Method of consecutive projections (MCP) is quite popular. The principal idea of the method is to present the convex set as intersection of simple elementary convex sets and project to sets, outer for the current point. To build projection onto elementary sets very easy, because elementary sets are semispaces. It was proved, that this algorithm is finitely convergent. Each iteration of MCP has 3 tasks: 1. choose elementary set for projection; 2. define direction; 3. calculate length of step. Typical features of the method for deconvolution problems are the large dimension and the presence of a large number of almost perpendicular directions of projection. We found that several sequential steps in perpendicular directions can be replaced by one step, indicating a simple way to build it. Although the theoretical speed of convergence of the method remained the same, for large-scale problems (of the order of hundreds of millions), characteristic, for example, for image processing, it was possible to reduce the counting time by hundreds of times.
Keywords: convex programming, projection onto set, method of consecutive projections, deconvolution.
References
1. Bregman L. M., "Finding the common point of convex set by the method of successive projections", Dokl. Akad. Nauk 162(3), 487-490 (1965). (In Russian)
2. Fazylov V. R., "A General Method of Searching for a Point of a Convex Set", Sov. Math. (Iz. VUZ) 27(6), 51-61 (1983).
3. Kulikov A. N., Fazylov V. R., "On Control of Penetration in the Support-Vector Algorithm with Compaund Step", J. of Math. Sci. 73(5), 571-576 (1995).
4. Youla D., "Mathematical Theory of Image Restoration by the Method of Convex Projections", In: Image Recovery: Theory and Application (H. Stark (ed.), Academic Press, 1987).
5. Marks R. J. II, "Alternating Projections onto Convex Sets", In: Deconvolution of Images and Spectra (P. A. Jensson (ed.), Academic Press, New York, 1997).
6. Eremin 1.1., Popov L. D., "Closed Fejer cycles for incompatible systems of convex inequalities", Izv. vuzov. Mat., issue 1, 11-19 (2008). (In Russian)
7. Schneider D., van Ekeris T., Jacobsmuehlen J., Gross S., "On benchmarking non-blind decon-volution algorithms: a sample driven comparison of image de-blurring methods for automated visual inspection systems", In: Instrumentation and Measurement Technology Conference (I2MTC), 2013 IEEE International, 1646-1651 (2013).
8. Tikhonov A.N., Arsenin V. Ya., Methods for solving ill-posed problems (Nauka Publ., Moscow, 1979). (In Russian)
9. Gornyy V. I., Kislitskij M.I., Latypov I. Sh., "Estimation of Efficiency of Algorithms for Synthesizing the Aperture of Scanning Radiometer", Sovremennye problemy distantsionnogo zondirovaniya Zemli iz kosmosa 7(2), 14-25 (2010). (In Russian)
Received: June 9, 2018 Revised: September 6, 2018 Accepted: September 27, 2018
Author's information:
Iskander Sh. Latypov — [email protected]
