Двойственный алгоритм для задач регуляризации с недифференцируемыми стабилизаторами Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 19, № 2, 2014

Двойственный алгоритм для задач регуляризации с недифференцируемыми стабилизаторами

А.С. Величко Институт автоматики и процессов управления ДВО РАН, Владивосток, Россия e-mail: [email protected]

Величко А.С. Двойственный алгоритм для задач регуляризации с недифференцируемыми стабилизаторами // Вычисл. технологии. 2014. Т. 19, № 2. С. 14-19.

Для класса некорректно поставленных задач рассматривается подход регуляризации с недифференцируемыми стабилизирующими функционалами с позиций теории и методов оптимизации. Предложено эквивалентное представление исходной негладкой задачи в виде квадратичной оптимизационной задачи с линейными ограничениями-неравенствами. Рассматривается также параллельный численный метод, основанный на двойственной постановке оптимизационной задачи и нелинейном методе Якоби. Представленный подход используется для решения некорректно поставленной задачи большой размерности для интегрального уравнения Фредгольма первого рода, которая, в частности, возникает в задаче восстановления гравитационного поля Земли в математической геофизике.

Ключевые слова: регуляризация, условная оптимизация, численные методы, параллельный алгоритм.

Velichko A.S. Dual algorithm for regularization problems with nondifferentiable stabilizing functionals // Comput. Technologies. 2014. Vol. 19, No. 2. P. 14-19.

The paper examines theoretical approaches and numerical methods for constrained quadratic optimization problems which arise in regularization problems with non-diffe-rentiable stabilizing functionals.

Key words: regularization problem, constrained optimization, numerical method, parallel algorithm.

Введение

В работах [1, 2] описан метод регуляризации для решения некорректно поставленной задачи для операторного уравнения А и = f, основанный на использовании недиф-ференцируемых стабилизирующих функционалов. Классическая тихоновская регуляризация не всегда позволяет с приемлемым качеством восстановить истинное решение и, поскольку функционалы, содержащие норму пространства Соболева, могут обладать эффектом "заглаживания" решения. Для повышения точности восстановления истинных решений требуется измельчение вычислительной сетки, что влечёт необходимость решать задачи большой размерности. Малые значения параметра регуляризации вызывают неустойчивость численных методов решения возникающих задач оптимизации.

Для решения возникающих задач негладкой оптимизации в [3, 4] использовались субградиентные, проксимальные методы и методы решения систем неравенств с помощью проективных методов. Указанные численные методы часто подвержены медленной скорости сходимости, а для задач большой размерности эта проблема усиливается.

Анализ альтернативных возможностей построения эффективных численных методов, в том числе параллельных, является предметом данного исследования.

1. Методы регуляризации и задачи оптимизации

В качестве примеров использования подхода регуляризации для некорректно поставленной задачи в работах [5, 6] рассматривалось интегральное уравнения Фредгольма ь

У К(I, в)и(в) ¿в = у(Ь), для которого можно записать различные задачи негладкой опти-

а

мизации. В настоящей работе для регуляризации уравнения Фредгольма предлагается задача безусловной недифференцируемой оптимизации вида

min < у hm

} I 1=1

I n X2

"Y^hnK(ti,sj)uj - yt j=i

+ ( hnK(ti,sj)uj - yi j=i

+ ahny] j=i

uj2

где Нт, кп — шаги сетки, а — параметр регуляризации, здесь и далее для краткости обозначений и^ = и (в^), yi = у (и).

Классическая тихоновская регуляризация предполагает решение квадратичной задачи безусловной оптимизации

n

2

min < ^ hm ( hnK (ti, Sj )uj - yA + ahn^2 U {Uj} ' i=l \j=l / j=l

Альтернативой использования специальных методов негладкой оптимизации для решения задач (1), (2) может быть их предварительное эквивалентное преобразование к задачам квадратичной условной оптимизации с линейными ограничениями-неравенствами, для решения которых существуют эффективные численные методы [7].

Для задачи (1) возможно её эквивалентное представление в виде оптимизационной задачи с квадратичной целевой функцией и n + m неизвестными переменными

mn

min < hm^2 wi + hm^2 W + uj

{uj ,Wi}

^ i=1 i=1 j = 1 и с 2m линейными ограничениями-неравенствами

nn

Wi hnK (ti, Sj )uj - yi, Wi > - I ^2 hnK (ti, Sj )uj - yi

j=i \j=i

Данный приём, позволяющий "избавиться" от негладкой оптимизируемой функции, не является универсальным, а обусловлен именно поиском минимума оптимизируемого выпуклого функционала. Возможна дополнительная замена переменных Zj = у/ащ, позволяющая частично решить проблему неустойчивости по выбору малого параметра регуляризации численных методов решения возникающих задач оптимизации.

2

2. Двойственное представление задачи, численные методы и алгоритмы

Запишем последнюю задачу в виде

( 1 , 1

шт < -х ^х + с х, Ах < Ь > ,

где х = (и, w)/, и и w — векторы из элементов {uj} и {wj} соответственно, Q =

акпЕп 0 \

, р — диагональная матрица,

0 ^пЕт )

A = Г :E;: ) ■ С = (0п.^твт), Ь = Г -у

Здесь Ek — единичная матрица размерности к, 0 — матрица из нулей соответствующей размерности, К — матрица из элементов {КSj)}, 0п — вектор-строка из нулей длины п, вт — вектор-строка из единиц длины т, у — вектор из элементов {у¿}.

Следуя [8], можно записать двойственную постановку рассматриваемой задачи:

Ф(р*) = шт р/Бр + а/р, р > 0}

Здесь О = AQ-1A/ — положительно определённая, симметричная матрица, в = Ь + AQ-1c и для оптимальных решений прямой (х*) и двойственной (р*) постановки задачи выполняется соотношение х* = :Q-1(c + A/p*). Размерность двойственного вектора р равна 2т — числу линейных ограничений-неравенств рассматриваемой квадратичной задачи условной оптимизации.

Идея параллельного алгоритма для решения двойственной задачи основана на применении последовательного нелинейного метода Якоби [8] и предлагаемой ниже модификации, позволяющей распараллеливание вычислений. Используемый для решения задач оптимизации нелинейный метод Якоби предполагает фиксацию значений всех, кроме одной, компонент вектора неизвестных. На каждом шаге алгоритма решается получающаяся одномерная задача оптимизации, обновляется значение компоненты вектора неизвестных на данной итерации алгоритма и производится переход к другой компоненте искомого вектора, при этом компоненты циклически перебираются до выполнения критерия останова алгоритма. В отличие от методов покоординатного спуска в рассматриваемом алгоритме не вычисляется градиент оптимизируемой функции в текущей точке и не осуществляется спуск в направлении антиградиента функции.

Шаг s + 1 предлагаемого параллельного алгоритма состоит из трёх подшагов.

1. Пусть р 5 = (р\,... ,р|,... ,р1)/ — приближённое решение, полученное на предыдущем шаге s. Формируется к = 2т одномерных квадратичных подзадач (2т — размерность двойственного вектора р) вида шт <|-р(г)Ор(г)' + в/р(г), р^ > 0^, где р(г) =

(р1,... ,рг,... ,рк)/, а индекс г принимает значения от 1 до к. Получаемые на данном подшаге одномерные задачи вида шт{аь2 + Ъь} при условии V > 0 для а > 0 очевидно имеют решение V* = шах{0, —Ъ/2а}. Вычисление к векторов р(г) осуществляется в виде параллельных процессов.

2. Определяется индекс i* из условия

i* = arg min <| 1 p(i)'Dp(i) + d'p(i)

Формируется новое приближенное решение ps+1 = (p\,... ,pi*,... ,psk)'.

3. Если выполняется критерий останова ||ps+1 — ps|| < £ для наперёд заданного малого числа £ > 0, то алгоритм завершает свою работу, возвращая приближённое решение ps+1, иначе — значение s увеличивается на единицу и алгоритм продолжает свое выполнение на следующем шаге.

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

Параллельный алгоритм для многопроцессорной вычислительной архитектуры с распределённой памятью реализован на языке Matlab и Octave [9] с использованием библиотек функций интерфейса передачи сообщений MPI в реализации OpenMPI [10] и "OpenMPI Extension for Octave" [11]. Расчёты проводятся на многопроцессорном вычислительном комплексе Центра коллективного пользования ДВО РАН во Владивостоке "Дальневосточный вычислительный ресурс" [12].

H

В численных расчётах ядро K(t, s) уравнения Фредгольма имеет вид --—-—,

(t s) +H

где H =10, что приводит к числу обусловленности для матрицы с элементами K(ti, si) порядка 1016.

Модельное решение

10 í s2l 10 ( (s — 4)2

u(s) = 472n exp\—32/+ Т2Лexp

сумма гауссианов на отрезке [—10,10]. Значения правой части уравнения Фредгольма y(t) брались не с использованием решения интегрального уравнения, а по аппроксима-ционной формуле для интеграла в узлах ti, когда значения u(sj) берутся по истинной (наперёд заданной) функции u(s) в узлах сетки sj. Параметр регуляризации a принимался равным 10-6.

На рис. 1, а показано восстановление модельного решения с помощью решения оптимизационной задачи (1) с её последующим эквивалентным преобразованием, предлагаемым в данной работе, на рис. 1, б — то же, для классической тихоновской регуляризации в результате решения задачи (2).

Величина "ускорения" Rn параллельного алгоритма определяется как отношение времени выполнения алгоритма на одном процессоре к времени выполнения алгоритма

n

на n параллельных процессорах. Известная теоретическая оценка Rn = 7-т-

(n — 1)а + 1

(закон Амдала), где a — доля непараллельного кода алгоритма. При "идеальном" распараллеливании процесса вычислений a = 0, и тогда Rn = n [8]. Для рассматриваемого в настоящей работе параллельного алгоритма на рис. 2 показаны график фактического ускорения (speedup) (сплошная линия) и аппроксимация теоретической зависимости для Rn с оценкой для доли непараллельного кода алгоритма a = 0.06 (штриховая

-1 -1-1-1-1-э- -0.5 -1-1-1-1-Э-

-10 -5 0 5 10 5 -10 -5 0 5 10 5

Рис. 1. Модельное (сплошная линия) и восстановленное (штриховая линия) решения для задач (1) (а) и (2) (б)

Рис. 2. График ускорения параллельного алго- Рис. 3. Характер сходимости алгоритма ритма

линия). На рис. 3 представлен график, демонстрирующий характер сходимости алгоритма: приведены значения 83 = |Ф(р5) — Ф(р*)|, т.е. абсолютные отклонения значений оптимизируемой функции, вычисленных для приближённых решений, получаемых на шагах алгоритма, от оптимального значения этой функции.

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

Заключение

В работе на примере решения некорректно поставленной задачи для интегрального уравнения Фредгольма первого рода, возникающей в задаче восстановления гравитационного поля Земли в математической геофизике, с позиций теории и методов оптимизации предлагается решение квадратичной оптимизационной задачи с линейными

ограничениями-неравенствами в двойственной постановке с помощью параллельного численного метода, родственного нелинейному методу Якоби.

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

Список литературы

[1] Васин В.В. Устойчивая аппроксимация негладких решений некорректно поставленных задач // Докл. АН. 2005. Т. 402, № 5. C. 586-589.

[2] Васин В.В., Агеев А.Л. Некорректные задачи с априорной информацией. Екатеринбург: Наука, 1993. 262 с.

[3] Васин В.В., Ерёмин И.И. Операторы и итерационные процессы фейеровского типа. Теория и приложения. Москва, Ижевск: НИЦ РХД, 2005. 200 с.

[4] Васин В.В. Итерационные процессы фейеровского типа в некорректных задачах с априорной информацией // Изв. высших учебных заведений. Математика. 2009. № 2. С. 3-24.

[5] Vasin V.V., Korotkii M.A. Tikhonov regularization with nondifferentiable stabilizing functionals //J. of Inverse and Ill-posed Problems. 2007. Vol. 15, No. 8. P. 853-865.

[6] Васин В.В., Сережникова Т.И. Регулярный алгоритм аппроксимации негладких решений для интегральных уравнений Фредгольма первого рода // Вычисл. технологии. 2010. Т. 15, № 2. С. 15-23.

[7] Нестеров Ю.Е. Введение в выпуклую оптимизацию. М.: МЦНМО, 2010. 274 с.

[8] Bertsekas D.P, Tsitsiklis J.N. Parallel and Distributed Computation: Numerical Methods. Nashua: Athena Sci., 1997.

[9] GNU Octave. URL: http://www.octave.org.

[10] OpenMPI. URL: http://www.open-mpi.org.

[11] OpenMPI Extension for Octave. URL: http://octave.sourceforge.net/openmpi_ext/index.html.

[12] Центр коллективного пользования ДВО РАН "Дальневосточный вычислительный ресурс". URL: http://www.cc.dvo.ru.

Поступила в 'редакцию 20 января 2014 г.