7.2. АЛГОРИТМЫ НЕЛИНЕЙНОГО ПРОГРАММИРОВАНИЯ ДЛЯ ВОССТАНОВЛЕНИЯ ИЗОБРАЖЕНИЙ

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

Необходимо минимизировать некоторый функционал

с учетом ограничений:

Существуют различные пути решения (7.7). Если результирующий функционал с учетом ограничений может быть представлен в виде суммы линейного и квадратичного функционалов, то решение можно найти аналитически. В противном случае необходимо создать вычислительные схемы для решения задачи минимизации.

Задачи нелинейного программирования в постановке (7.7) применительно к восстановлению изображений могут решаться двумя путями: сведением к системе нелинейных уравнений и использованием методов прямой оптимизации.

Сведение к системе нелинейных уравнений и ее решение.

Будем рассматривать нелинейный функционал с ограничениями как нелинейную функцию от ряда дискретных аргументов:

где функция; — множители Лагранжа. Условием минимума функционала является равенство нулю частных производных:

Мы получили систему нелинейных уравнений

Для решения системы (7.8) могут быть использованы различные методы, например, метод Ньютона — Рафсона, заключающийся в вычислении элементов якобиана:

и выполнении итераций вида

Именно таким методом решались большинство задач восстановления на основе метода максимума энтропии [60]. Однако для каждой итерации при вычислениях по формуле (7.10) необходимо вычислять якобиан и производить его обращение.

Методы прямой оптимизации.

Эти методы являются более привлекательными. Пусть нам задан нелинейный функционал

где вектор: являются множителями Лагранжа. Задача заключается в нахождении глобального экстремума необходимым условием которого является равенство нулю градиента:

Простейший метод нахождения экстремума, так называемый метод градиентного спуска, был предложен еще Коши и заключается в следующей последовательности операции:

1) выбор начального приближения

2) нахождение

3) выбор

4) переходим к п. 2, если не достигнута сходимость,

в противном случае получаем решение

Метод градиентного спуска заключается в выборе такого направления, при котором градиент максимален. Обратимся к геометрической интерпретации этого метода (рис. 7.1). Алгоритм подразумевает перемещение по склону в том направлении, в котором крутизна подъема максимальна; величина шага определяется константой а. При выборе а возможны различные варианты.

Рис. 7.1. Условное представление за дачи поиска экстремума

Рис. 7.2. Структура метода градиентного поиска с постоянным шагом приращения

Если сохранить значение а постоянным для всех итераций, то существует опасность «проскакивания» максимума и колебаний около него. На рис. 7.2 эта ситуация изображена на плоскости, в которой нанесены контуры Возможен алгоритм деления шага, заключающийся в том, что если на следующем шаге происходит необходимое изменение функционала (возрастание в случае максимизации или убывание для минимизации), то шаг а оставляют прежним, если нет, то шаг делят пополам. Критерием, по которому судят о выборе шага для случая максимизации, является выполнение условия

Если (7.11) выполняется, то шаг оставляем прежним, если нет — шаг уменьшаем вдвое [71].

Более оптимальным является алгоритм, в котором шаг при максимизации выбирается на каждой итерации, исходя из условия [164]

Этот метод называется также методом наискорейшего спуска. Внимательный анализ, однако, показывает, что и у (7.12) имеются свои недостатки: чрезвычайно медленная сходимость в случае вытянутых линий уровня («овражная структура» поверхности) (рис. 7.3), а также опасность локальных максимумов (рис. 7.4). Действительно, в ситуации, изображенной на рис. 7.4, алгоритм никогда не находит глобального решения [64].

Наилучшую сходимость среди рассматриваемых методов имеют методы сопряженных градиентов, основанные на использовании последовательности направлений поиска, которая является линейной комбинацией градиента

Рис. 7.3. «Овражная» структура поверхности и осцилляции решения, возникающие при использовании метода дробления шага

Рис. 7.4. Локальные максимумы и неоднозначность решения

и предыдущего направления поиска. Окончательная схема алгоритма сопряженных градиентов выглядит следующим образом. Каждый последующий шаг вычисляется по формуле

где

Входящая в это уравнение величина может быть вычислена по формуле

где круглые скобки означают скалярное произведение, а множество индексов, на которых происходит «обновление» решения, например, каждый третий, каждый пятый и т. п. Отметим, что практическая реализация метода сопряженных градиентов требует проведения на каждом шаге задачи одномерной оптимизации, а это несмотря на увеличение скорости сходимости, замедляет время выполнения алгоритма.

Экспериментальные исследования показывают, что в тех случаях, когда решается задача нелинейного программирования с несколькими ограничениями, исследуемая «поверхность» очень сложна и приобретает «овражную» структуру с локальными максимумами, так что добиться сходимости обычными градиентными методами становится очень непросто; в то же время метод сопряженных градиентов,

как правило, работает вполне устойчиво и надежно сходится к искомой оценке.

В качестве практического примера рассмотрим реализацию алгоритма максимума энтропии с ограничениями на функцию автокорреляции (модуль спектра) и уравнения формирования и регуляризацией в форме тихоновского функционала.

В дискретной форме общий функционал, подлежащий оптимизации, выглядит следующим образом:

Сформируем частные производные функционала (7.14)

Далее будем использовать схему метода сопряженных градиентов:

1. Находим

2. Находим величину

3. Формируем

и решаем задачу одномерной минимизации

4. Находим следующий шаг

Если условие сходимости (например, разность в метрике между последующим и предыдущим шагами не превышает б) не выполняется, то возвращаемся к в противном случае есть искомое решение.

Один из полезных приемов улучшения сходимости методов прямой оптимизации состоит в модификации структуры функций Лагранжа. Так, на практике помогает следующий прием. Если известны ограничения в виде равенств то в общую схему функционала они включаются, взятые в квадрат:

Известны строгие математические результаты о том, что при решение задачи стремится к истинному. Широко используется также модификация (7.16), называемая методом модифицированных функций Лагранжа (МФЛ) [57]. Этот метод заключается в формировании следующего выражения:

где а — коэффициент штрафа. Метод решения (7.17) имеет довольно сложную, двуступенчатую структуру. Сначала решается задача прямой оптимизации (например, при использовании метода сопряженных градиентов), затем делается пересчет множителей Лагранжа

которые используются на следующей итерации метода сопряженных градиентов. Метод МЛФ особенно хорошо зарекомендовал себя при наличии нескольких различных ограничений. Например, экспериментальные результаты, приведенные в гл. 3, основаны на решении задачи оптимизации методом сопряженных градиентов с использованием модифицированных функций Лагранжа.

Рассмотрим некоторые особенности решения задач нелинейного программирования на универсальных ЭВМ.

1. Необходимое число точек задачи оптимизации. Будем считать формат изображения равным отсчетов. Если задача оптимизации заключается в нахождении

экстремума функционала при наличии одного ограничения на формирование изображения, то число точек задачи оптимизации составляет в случае К ограничений эта размерность равна Так, для изображения 256X256 элементов получаем число точек, по крайней мере, 131070.

2. Требуемый объем памяти. Если используется метод сопряженных градиентов, то нетрудно видеть, что кроме собственно изображения нам необходимо хранить искаженное изображение весовую функцию множителей Лагранжа в случае К ограничений, градиенты функционала на настоящем и предыдущем шаге, константы, т. е. всего около точек. Для изображения 256X256 элементов это будет занимать, по крайней мере, 600 Кбайт памяти. Естественно, что хранить необходимые данные и результаты, как правило, придется на внешних устройствах (дисках, лентах).

3. Объем вычислений. На каждой итерации необходимо вычислять приблизительно 5 двумерных сверток длиной элементов. Кроме этого, необходимо выполнять операции логарифмирования, возведения в квадрат, многократного вычисления значения функционала при решении задачи одномерной оптимизации. Если предположить, что для решения задачи одномерной оптимизации необходимо 10 итераций, то для каждой итерации метода сопряженных градиентов понадобится 50 сверток. Пусть процесс сходится за 20 итераций, тогда для изображений 256X256 элементов, если как и ранее предполагать получим время отработки не менее 80 мин. Этот вывод подтверждается результатами обработки изображения 64X64 элемента на мини-ЭВМ с производительностью порядка 600 000 приведенных Операций. Для восстановления изображения 64X64 элемента требовалось от 20 до 40 итераций; при этом среднее время (с учетом времени обмена с диском) обработки составило около 20 мин (см. § 7.5).

Таким образом, реализация задач нелинейного программирования на универсальных ЭВМ требует в 3—4 раза большего машинного времени и памяти, чем реализация итерационных алгоритмов с ограничениями.

В случае использования специализированных ЭВМ, приспособленных для работы с большими массивами и имеющих аппаратурно реализованные средства вычисления БПФ, в том числе при использовании специализированных оптико-электронных вычислительных комплексов, время выполнения задачи снижается на порядок.

Методы прямой оптимизации, рассмотренные в этом

параграфе, сводятся к следующей простой вычислительной схеме: выполняется процесс последовательных итераций до тех пор, пока не будет обеспечена сходимость

где градиенты функционала, шаги приращений, которые могут выбираться исходя из некоторых соображений; множители Лагранжа и коэффициенты штрафных функций.

В уравнениях (7.18) градиенты вычисляются с помощью алгоритмов БПФ, поэтому этот процесс представляет собой метод прямой оптимизации с вычислением итераций через . Отметим, что если не рассматривать вопросы создания новых алгоритмов оптимизации, средства параллельной обработки позволяют эффективно решать задачи нелинейного программирования. Действительно, легко заметить, что (7.18) представляет собой многомерный итерационный процесс, в целом аналогичный итерационным алгоритмам. Действительно, если бы процесс являлся «одноканальным», то он мог бы быть записан в следующей форме:

Выражение (7.19) совпадает с формой записи итерационных алгоритмов восстановления. Отличие (7.18) от (7.19) заключается в том, что приращение аргумента зависит от нескольких параметров (множителей Лагранжа), определяемых параллельными итерационными процессами для каждого ограничения задачи. Представим себе, что процесс вычислений (7.18) удалось сделать параллельным. На рис. 7.5 изображены несколько параллельных вычислительных устройств, имеющих общий

Рис. 7.5. Параллельная структура вычислений при решении многомерной задачи нелинейного программирования

блок данных Каждое из них вычисляет свой итерационный процесс из набора процессов (7.18); при этом после каждой итерации результаты общего блока данных обновляются и в зависимости от конкретной схемы используются в различных вычислительных каналах. Схема на рис. 7.5 может быть реализована с помощью параллельных процессоров или многоканального оптического вычислительного устройства, входящего в состав оптико-электронного вычислительного комплекса.

В заключение остановимся на возможности создания более эффективных алгоритмов решения задач оптимизации. В этой области не решено еще множество задач. Действительно, практически все описанные нами градиентные алгоритмы строятся на основе «локальных» принципов: выбирается точка на поверхности, обследуется ее окрестность, начинается движение на поверхности и т. д. Таким образом, алгоритмы градиентного типа основываются на знании малой окрестности точки поверхности. Возникает естественный вопрос о возможности построения более «глобального» алгоритма. Например, можно предложить следующий простой алгоритм, реализация которого изображена на рис. 7.6. Суть работы такого алгоритма заключается в разумном сочетании «локальной» и «глобальной» стратегий:

1) выбираем начальное приближение ;

2) находим значение

3) выбираем (или вместо двойки другой коэффициент, больший единицы; возможен выбор другого коэффициента, меняющегося от шага к шагу);

Рис. 7.6. Возможный алгоритм нелокального поиска решения

4) находим возможное пересечение поверхностей Если есть решение и найдена «точка» f такая, что то переходим к п. 3. В противном случае уменьшаем шаг например, и снова переходим к п. 4.

Если найден эффективный метод определения факта пересечения то предлагаемый алгоритм должен обладать хорошей сходимостью и, что самое главное, не зависеть от «овражной» структуры и локальных максимумов. Иначе говоря, если решение найдено, то оно является «высоким» пиком (см. рис. 7.6). Кроме того, в проекции на плоскость линий уровня поверхности предлагаемый алгоритм может рассматриваться как реализация наиболее короткого пути, близкого к прямой в окрестности глобального максимума.

Основной проблемой, возникающей при реализации подобного алгоритма, является эффективный способ определения пересечения Для этого достаточно было бы решения вида да (нет), однако единственным доступным путем является градиентный поиск в плоскости Для этого необходимо сформировать новый функционал, например

и решать задачу на максимум обычным градиентным методом: выбирая находим и проводим обычный градиентный поиск. Преимуществом градиентного поиска рассматриваемого вида является гладкость поверхности так что здесь может достаточно эффективно работать быстрый алгоритм с дроблением шага. Кроме того, возможна модификация алгоритма с учетом того обстоятельства, что в точке пересечения и можно «обострить» значения около точки пересечения и ввести некоторый порог . Если значение меньше порога, то считается, что пересечение не найдено, а если значение попадает в область то считается, что факт пересечения установлен. Возможен также случайный поиск в плоскости

В итоге приведенная задача сведена к: набору задач оптимизации, каждая из которых может быть сделана устойчивой к сложноц структуре поверхности. Результирующий алгоритм должен всегда сходиться к глобальному максимуму.

Таким образом, необходима разработка эффективных методов решения многомерных задач нелинейного программирования. Эта задача может решаться несколькими путями:

разработкой новых, более эффективных методов оптимизации;

разработкой эффективных численных алгоритмов, реализующих оптимизацию;

повышением производительности ЭВМ, использованием параллельных вычислений, специализированных процессоров и БПФ-процессоров, оптико-электронных комплексов, увеличением памяти ЭВМ.