CC BY
34
удк 681.3.06
И.Г. Черноруцкий
практическая оптимизация и невыпуклые задачи
I.G. Chernorutskiy PRACTICAL OPTIMIZATION AND NONCONVEx PROBLEMS
рассмотрены основные особенности задач оптимизации при проведении реальных компьютерных вычислений. обоснована необходимость учета таких характеристик прикладных экстремальных задач, как невыпуклость и жесткость целевых функционалов. Показано, что в указанных условиях классические оптимизационные методы ньютоновского и квазиньютоновского типов испытывают существенное влияние погрешностей и часто приводят к расходящимся процедурам.
методы оптимизации. невыпуклые задачи. ньютоновские методы.
Main aspects of optimization problem in real computing are described. Requirement to take into account such parameters of applied extremal problem as rigidity and nonconvexity of criterion functional is founded. It is established that classic Newtonian and quasi-Newtonian optimization procedures are strongly influenced by inaccuracies in this case and nonconvergent procedures may appear.
optimization procedure. nonconvex optimization problem. newtonian procedure.
рассматриваются практические аспекты применения методов однокритериальной конечномерной оптимизации (математического программирования) в задачах реального компьютерного моделирования при решении прикладных оптимизационных задач [1, 2]. Для определенности будем говорить о задачах поиска минимума. основное внимание уделяется таким характерным чертам решаемых задач, как невыпуклость и плохая обусловленность (жесткость или овражность) минимизируемых функционалов. В теоретических работах по математическому программированию большой объем занимают исследования в области выпуклых задач оптимизации. Именно для таких задач могут быть получены наиболее законченные математические результаты, в частности, по доказательству сходимости алгоритмов минимизации (максимизации), установлению скорости сходимости и т. д. однако, к сожалению, на практике выпуклые задачи встречаются редко. Утверждения, что в достаточно малой окрестности оптимума минимизируемая функция
(функционал) почти всегда становится выпуклой, лишены практического смысла, т. к. попадание в эту окрестность обычно и означает решение задачи. следующий пример показывает, что выпуклости нет уже в простейших случаях.
Пример. Задача параметрической идентификации. Пусть некоторая реальная система описывается скалярным дифференциальным уравнением вида (задача коши)
— = ах, х(0) = х0 с решением х(7) = х0еа.
начальное условие х0 считаем известным.
требуется по результатам измерения функции х(?) получить оценку неизвестного параметра а. Пусть начальное значение х0 = 1, а эталонное (искомое) значение равно а = а' = 1. тогда эталонная (измеряемая) функция имеет вид х (7) = е'.
найдем оценку неизвестного параметра а из сравнения эталонного решения с расчетным в соответствии с методом наименьших квадратов (задача параметрической идентификации).
4
Мы наблюдаем функцию et. Расчетное решение при произвольном параметре a и заданном x0 = 1 имеет вид x(t) = eat.
Найдем a из условия минимума:
J(a) = (eat - et)2 ^ min.
a
Сравнение расчетной и эталонной (измеренной) зависимости будем производить в одной точке для t = 1:
Ji(a) = (ea - e)2.
График этой зависимости представлен на рисунке.
Видно, что целевой функционал оказывается невыпуклым в значительной области изменения аргумента слева от оптимума. В более сложных случаях ситуация только усугубляется. В частности, легко проверить, что характер приведенной в примере зависимости сохраняется и при увеличении числа точек сравнения.
Таким образом, при решении практических задач оптимизации необходимо ориентироваться на методы, сохраняющие эффективность в невыпуклых ситуациях.
ниже дан краткий обзор и характеристика некоторых популярных методов конечномерной оптимизации. Основное внимание уделяется анализу работоспособности алгоритмов при использовании невыпуклых целевых функционалов. набор рассматриваемых методов по необходимости ограничен и включает в основном ньюто-
новские и квазиньютоновские процедуры, считающиеся в настоящее время наиболее эффективными при решении задач с гладкими функционалами. Одной из характерных особенностей решаемых на практике однокритериальных оптимизационных задач является также высокая степень жесткости целевого функционала [3]. Поэтому необходим анализ методов с позиций одновременного присутствия обеих характеристик — невыпуклости и жесткости целевого функционала.
Ньютоновские методы. Ньютоновские методы строят последовательность точек {xk} в соответствии с алгоритмом [ 4 ]:
Х*+1 = X* + ^, К Е Л1, (1)
где pN = -[^Г1 як определяет т. н. ньютоновское направление спуска; кк — длина шага;
Ок й / "(хк);я й/ '(хк).
Предполагается, что все матрицы Ок положительно определены (Ок > 0). Последнее
условие гарантирует разрешимость задачи
к
вычисления pN исходя из системы уравнений
GkPkN =
(2)
Вектор х1" + 1, построенный согласно (1), (2) при кк =1, является минимизатором аппроксимирующего квадратичного функционала, полученного как отрезок соответ-
Невыпуклая целевая функция
ствующего ряда тейлора:
Fk (x) = J (xk) + (
g , x - x
,(Gk(x - xk), x - xk).
(3)
Ньютоновские методы (1) оказываются эффективными при решении задач безусловной минимизации выпуклых функционалов при невысоких степенях жесткости целевых функционалов. кроме того, предполагается существование достаточно точных локальных квадратичных аппроксимаций основного функционала. Некоторые авторы утверждают, что это наиболее эффективные из всех применяемых ими для решения реальных задач методов при условии выпуклости целевого функционала. Однако это утверждение является слишком сильным. Последующий анализ и опыт реальных вычислений подтверждают это замечание.
Если матрица бк не является положительно определенной, параболоид (3) может не иметь конечного минимизатора и, как правило, будет неограничен снизу. Поэтому для невыпуклых целевых функционалов процедура (2) непосредственно неприменима. Ньютоновское направление ркм при наличии отрицательных спектральных составляющих матрицы бк может указывать в сторону возрастания функционала /(х), что соответствует расходимости процесса.
Как отмечено в [5], «в настоящее время нет общепринятого определения «метода Ньютона» для расчета направления спуска при знаконеопределенной матрице бк, поскольку среди специалистов нет согласия относительно того, как использовать локальную квадратичную аппроксимацию ¥к в этом случае». Такая ситуация сохраняется и поныне.
Наиболее употребительны модификации метода Ньютона, в которых направление спуска рк находится из решения линейной системы:
„к к Gkp = -g
(4)
где бк — некоторая положительно определенная матрица, совпадающая с исходной матрицей Гессе бк, если последняя положительно определена. Указанный метод
выбора направления рк гарантирует глобальную сходимость процедуры независимо от характера выпуклости функционала /(х) в окрестности текущей точки хк.
Ниже рассмотрены наиболее известные процедуры построения матрицы бк на основе различных матричных разложений, позволяющих учитывать знаки собственных чисел б,.
к
Методы, основанные на спектральном разложении. Спектральное разложение матрицы б имеет вид
Gk = UKU t = g},
gj
= 1^
umum,
(5)
где {um, ных
U = [u1, u2, ..., un]; Л = diag (A,.); m e [1:n]} — система ортонормаль-собственных векторов матрицы Gk; {Am, m e [1:n]} — спектр матрицы Gk. Базовая схема, реализующая рассматриваемый подход, основана на выборе [5]:
Gk = UЛит, Л = diag [max {|A. |, 8}], (6)
где 8 > 0 — параметр метода, определяющий границу «существенной положительности» любого из собственных чисел. отмечаемый в литературе недостаток метода (4) с матрицей (6) связывается с трудоемкостью процедуры построения спектрального разложения (5), требующей до 4n3 арифметических операций. Кроме того, возникают известные из вычислительной линейной алгебры трудности в определении малых собственных чисел плохо обусловленной матрицы. Выбор параметра 8 также до конца не алгоритмизирован.
Методы, основанные на факторизации Холесского [5]. Разложение Холесского для симметричной положительно определенной матрицы В имеет вид
B = LDLt, D = diag (d), (7)
где L — нижняя треугольная матрица с единицами на диагонали; D — положительная диагональная матрица. Факторизация (7) непосредственно неприменима к знаконеопределенной симметричной матрице B = Gk. В [5] предложена модифицированная процедура, позволяющая равномерно ограничить рост элементов треугольного
m=1
фактора Ь на уровне | гк \ < р( I > к, гк й I ) и гарантирующая «существенную» положительность ¿. > 8 > 0 диагональных элементов матрицы Б. В результате получаем
вк = ЬБЬ = вк + С, вк > 0,
(8)
где С — неотрицательная диагональная матрица, зависящая от выбранного параметра р. Для сохранения численной устойчивости процедуры построения вк , а также для совпадения вк и вк в случае положительно определенной вк, целесообразно вычислять
в из условия р2 = шах{у, Vп2 - 1, ем}, где ем — машинное эпсилон; Е, — максимальный модуль недиагонального элемента вк; у — значение максимального из диагональных элементов в .
к
Основная цель построения разложения (8) заключается в сокращении вычислительных затрат по сравнению с (6)
1
приблизительно до -^п3
6
арифметических
операций. Существо дела при этом не затрагивается.
Аналогичные модификации метода ньютона, основанные на других численно устойчивых процедурах факторизации зна-конеопределенных симметричных матриц, рассматриваются во многих работах. Большое число публикаций, посвященных рассматриваемым вопросам, с одной стороны, указывает на актуальность проблемы, а с другой — на отсутствие метода, полностью отвечающего предъявляемым практикой требованиям.
Методы доверительной окрестности. Методы разработаны для выпуклых задач оптимизации. Основная идея методов доверительной окрестности сводится к следующему [5].
Рассматриваются квадратичные методы, т. е. методы, основанные на последовательной квадратичной аппроксимации минимизируемого выпуклого функционала. Надежность прогноза в таких методах определяется областью справедливости (достаточной точности) локальной квадратичной модели. В методах Ньютона мы определяем направление в сторону минимума аппроксимирующей квадратичной зависи-
мости, а затем регулируем норму вектора продвижения, чтобы не уйти «слишком далеко». Здесь же предлагается поступать по иному: решается задача условной минимизации аппроксимирующей квадратичной функции при условии, что норма вектора продвижения ограничена сверху некоторым заданным параметром.
Таким образом, минимизирующая последовательность строится по правилу хк + 1 = хк + рк, где вектор рк на каждом шаге определяется как решение ной задачи вида
вспомогатель-
як, Р
^ шт,
2 р)+
р Е Б с Я",
Б = {р Е Лп|(р,р)2 < Д}.
(9)
Величина Д характеризует область (окрестность) «квадратичности» исходного функционала. Сформулированная задача может быть решена методом множителей Лагранжа.
Условие стационарности функции Лагранжа для задачи (9) приводит к методу определения р из системы линейных уравнений
[вк + рЕ]р = - як, р > 0,
(10)
где р играет роль множителя Лагранжа. Алгоритм выбора р зависит от конкретной реализации МДО. Возможен непосредственный подбор оптимального р на основе многократного решения линейной системы (10). В ряде случаев вначале полагают р = 0. Если в процессе решения (10) при р = 0 выясняется, что вк знаконеопределена или ||р|| > Д, где Д — установленное пороговое значение, то р определяется как решение нелинейного уравнения вида
||[вк + ре — vi = Д.
(11)
Далее по найденному р согласно (10) определяется искомый вектор р.
Основной недостаток всех методов рассматриваемого класса состоит в необходимости решения на каждом шаге итерационного процесса нелинейных алгебраических уравнений типа (11), что требует привлечения различных процедур регуляризации.
кроме того, найденному оптимальному значению может отвечать очень плохо обусловленная матрица бк + рЕ, что приводит к значительным вычислительным трудностям при решении системы (10). Аналогичные недостатки присущи и ньютоновским методам.
Указанные проблемы становятся практически неразрешимыми, если степень обусловленности матрицы бк настолько высока, что информация о малых спектральных составляющих оказывается полностью утерянной на фоне больших собственных чисел за счет погрешностей задания (?к. Эти погрешности могут вызываться приближенным характером соотношений, применяемых для расчета производных, а также ограниченностью разрядной сетки компьютера. В указанных условиях, даже при использовании регуляризированных форм матричных разложений, гарантирующих, в частности, _положительную определенность матрицы б к, генерируемые векторы рк, оставаясь направлениями спуска, оказываются практически случайными, что резко замедляет сходимость соответствующей поисковой процедуры. сказанное иллюстрируется следующим примером.
Численный пример. Вырождение метода Ньютона. Рассмотрим сильно выпуклый квадратичный функционал
/(х) = 1/2 (Ах, х) - (Ъ, х) + с, А > 0. тогда точка минимума определяется решением системы линейных алгебраических уравнений (необходимое условие экстремума /'(х) = 0):
Ах = Ъ. (12)
спектральное разложение симметричной матрицы А имеет вид А = ТБТт, Т = (и1,и2,...,ип), где столбцы и' образуют ортонормальную систему собственных векторов матрицы А, а Б — диагональная матрица собственных чисел X1 матрицы А.
Получаем следующее выражение для минимизатора функционала / (х):
А = 105
А = 105
, X = 0, Х2 = 105.
х * = А-гЪ = ТБгхТ тЪ = = Т -(Ъ, и1), ..., —(Ъ, ип)
(13)
Пусть в некотором модельном вычислительном устройстве возможно четырехзначное представление мантиссы числа в формате с плавающей точкой. И пусть точная матрица А имеет вид
0.5001 0.5
0.5 0.50001
Х1 = 1, Х2 = 105 + 1.
Видно, что значения элементов матрицы определяются «большим» собственным числом Х2.
После записи матрицы в память вычислительного устройства и округления значений элементов до четырех знаков мантиссы получим вырожденную матрицу:
0.5 0.5"
^ 0.5 0.5_
В этих условиях классический метод Ньютона вообще неприменим, а его регу-ляризованные формы будут давать по существу случайные значения для «новых» малых собственных чисел и не будут содержать полезной информации об искомом векторе х *. Плохая обусловленность системы (12) для жестких функционалов приводит к полной потере полезной информации.
Из рассмотренного примера наглядно видно, что для жестких даже сильно выпуклых функционалов резко возрастает влияние погрешностей, в т. ч. и погрешностей округления. В частности, изначально выпуклая задача может стать невыпуклой (матрица получает отрицательные собственные числа). Влияние высокой степени жесткости на точность получаемого решения определяется следующим.
И з аналитического представления (13) решения х* следует, что значения собственных чисел входят в знаменатели соответствующих дробей и основное влияние на результат оказывают «малые» по модулю собственные числа. В то же время из спектрального разложения видно, что сами элементы матрицы представляют из себя некоторые линейные комбинации ее спектральных составляющих (собственных чисел) и численные значения элементов определяются в основном «большими» по
модулю собственными числами. В итоге небольшие относительные погрешности в представлении элементов матрицы Гессе приводят к большим относительным погрешностям для «малых» спектральных составляющих, а те, в свою очередь, полностью искажают окончательное решение. В этом, если угодно, и заключается «физический смысл» рассмотренного явления плохой обусловленности.
По этим же причинам любые «модификации» матрицы Гессе, в частности, выполняемые с целью аппроксимации невыпуклых задач выпуклыми, также приводят к потере полезной информации о «малых» собственных числах и не позволяют получить достаточно эффективные численные процедуры. Аналогично можно показать, что в невыпуклых квадратичных аппроксимациях жестких функционалов малые спектральные составляющие также несут основную полезную информацию о траектории спуска и основная цель состоит в построении алгоритмов, использующих в полной мере эту информацию.
Квазиньютоновские методы [5]. Квазиньютоновские методы, так же как и ньютоновские, основаны на использовании квадратичных моделей минимизируемых функционалов. Однако аппроксимация матрицы Гессе (либо обратной к ней) осуществляется последовательно, на основе наблюдений за изменением градиентов целевых функционалов в последовательных точках минимизирующей последовательности.
Часто указывается, что квазиньютоновские методы существенно более эффективны, чем ньютоновские, по причине более низкого порядка производных, определяющих схему метода. Однако реальная ситуация значительно сложнее. Пусть, например, минимизируется сильно выпуклый квадратичный функционал /(х). Тогда для построения оптимизатора х* квазиньютоновскому методу потребуется в общем случае вычислить п градиентов, где п — размерность пространства поиска. Последнее эквивалентно 2п2 вычислениям значений функционала /(х) если используются двусторонние конечно-разностные аппроксимации первых производных. Указанного количества значений
функционала J(x), очевидно, достаточно для построения аппроксимации J"(x), необходимой для реализации любого из вариантов метода Ньютона. Следовательно, трудоемкости рассматриваемых процедур в указанных условиях приблизительно равны, если не учитывать дополнительные вычислительные затраты на процедуры одномерного поиска в квазиньютоновских методах. С другой стороны, доказано, что большинство вариантов квазиньютоновских методов (например, одна из наиболее эффективных схем Бройдена—Флетчера—Гольдфарба—Шенно) при минимизации сильно выпуклых квадратичных функционалов приводят к одной и той же траектории спуска, вырождаясь в хорошо изученные методы сопряженных градиентов.
Методы сопряженных градиентов. В то
же время известна особенность методов сопряженных градиентов, существенно ограничивающая область их эффективного применения. Она заключается в понижении скорости сходимости для плохо обусловленных (жестких) задач оптимизации. Оценки, приведенные, например, в [6], показывают, что стандартные методы СГ сходятся по закону геометрической прогрессии со знаменателем q, близким к единице:
q = 1 - 2Д/П,
где п — степень жесткости минимизируемого функционала. Там же имеются указания на достаточно высокую скорость сходимости метода СГ «по функционалу», независимо от величины п
l\\ Х0 -xif
J(xk) - J(x') -1-, L = const. (14)
V ^ V ^ 2(2 к +1)2
Однако оценки типа (14) получены в предположении G(x) > 0. Кроме этого, согласно (14) эффективно получаются значения функционала порядка J(x *), где x * — минимизатор аппроксимирующего параболоида fx), что в общем случае не решает задачи. В соответствии с (14) реализуется указанная скорость сходимости по функционалу, но не по аргументу.
Предположение о невыпуклости вносит дополнительные трудности. Показано, что в этих условиях метод СГ по характеристи-
кам сходимости эквивалентен градиентному методу наискорейшего спуска со всеми вытекающими отсюда последствиями. Кроме отмеченных дефектов, общих для методов сопряженных градиентов и квазиньютоновских методов, последние имеют дополнительные недостатки, связанные с проблемой потери положительной определенности квазиньютоновских матриц из-за накопления вычислительных погрешностей в рекуррентных процедурах аппроксимации матриц Гессе [5].
В настоящее время на основе экспериментальных данных принята точка зрения, согласно которой ньютоновские методы с конечно-разностной аппроксимацией матрицы Гессе на основе аналитических выражений для градиентов надежнее, чем квазиньютоновские методы, и сходятся в более широком классе случаев, в частности, в задачах с очень плохой обусловленностью матрицы Гессе. С другой стороны, согласно тем же источникам, если аналитические выражения для градиентов отсутствуют и
применяются версии нулевого порядка соответствующих алгоритмов, в большинстве случаев более эффективными оказываются квазиньютоновские методы с конечно-разностной аппроксимацией градиентов. Данные замечания не снимают отмеченные выше недостатки указанных процедур, а отражают лишь некоторые сравнительные оценки, основанные на опыте проведения реальных вычислений.
Практические экстремальные задачи часто оказываются невыпуклыми и плохо обусловленными (жесткими). Основная полезная численная информация содержится в значениях «малых» по модулю собственных чисел матрицы Гессе. эффективные алгоритмы минимизации жестких невыпуклых функционалов должны строиться не на ньютоновской или квазиньютоновской идеологии, приводящей к потере информации о малых спектральных составляющих, а на основных принципах теории жестких систем оптимизации [7—9].
список литературы
1. Черноруцкий, И.Г. Методы параметрической оптимизации в задачах идентификации [Текст] / И.Г. Черноруцкий //Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. —СПб.: Изд-во Политехнического ун-та, 2009. -№ 2 (76). -С. 150-155.
2. Черноруцкий, И.Г. Параметрические методы синтеза систем управления [Текст] / И.Г. Черноруцкий //Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. -СПб.: Изд-во Политехнического ун-та, 2009. -№ 2 (76). -С. 111-115.
3. Черноруцкий, И.Г. Алгоритмические проблемы жесткой оптимизации [Текст] / И.Г. Черноруцкий //Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. -СПб.: Изд-во Политехнического ун-та, 2012. -№ 6 (162). -С. 141-152.
4. Ортега, Дж. Итерационные методы решения нелинейных систем уравнений со многими неизвестными [Текст] Дж. Ортега, В. Рейнболдт.
-М.: Мир, 1975. - 560 с.
5. Гилл, Ф. Практическая оптимизация [Текст] / Ф. Гилл, У. Мюррей, М. Райт. -М.: Мир, 1985. - 509 с.
6. Поляк, Б.Т. Введение в оптимизацию [Текст] / Б.Т. Поляк. -М.: Наука, 1983. - 384 с.
7. Ракитский, Ю.В. Численные методы решения жестких систем [Текст] / Ю.В. Ракитский, С.М. Устинов, И.Г. Черноруцкий. -М.: Наука, 1979. - 208 с.
8. Черноруцкий, И.Г. Функции релаксации градиентных методов [Текст] / И.Г. Черноруцкий //Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. -СПб.: Изд-во Политехнического ун-та, 2012. -№ 3 (150). -С. 66-72.
9. Черноруцкий, И.Г. Некоторые стандартные схемы параметрической оптимизации [Текст] / И.Г. Черноруцкий //Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. -СПб.: Изд-во Политехнического ун-та, 2012. -№ 6 (162). -С. 128-133.
references
1. Chernorutskii I.G. Metody parametricheskoi Telekommunikatsii.
optimizatsii v zadachakh identifikatsii /Nauchno-tekhnicheskie vedomosti SPbGPU. Informatika.
Upravlenie. -St.-Peters-burg: Izd-vo Politehn. un-ta, 2009. -№ 2 (76). -S. 150-155. (rus)
2. Chernorutskii I.G. Parametricheskie metody sinteza sistem upravleniia /Nauchno-tekhnicheskie vedomosti SPbGPU. Informatika. Telekommu-nikatsii. Upravlenie. —St.-Petersburg: Izd-vo Po-litehn. un-ta, 2009. -№ 2 (76). -S. 111-115. (rus)
3. Chernorutskii I.G. Algoritmicheskie problemy zhestkoi optimizatsii /Nauchno-tekhnicheskie vedomosti SPbGPU. Informatika. Telekommunikatsii. Upravlenie. -St.-Petersburg: Izd-vo Politehn. unta, 2012. -№ 6 (162). -S. 141-152. (rus)
4. Ortega Dzh. Iteratsionnye metody resheniia nelineinykh sistem uravnenii so mnogimi neizvest-nymi. -Moscow: Mir, 1975. - 560 s. (rus)
5. Gill F. Prakticheskaia optimizatsiia. -Moscow: Mir, 1985. - 509 s. (rus)
6. Poliak B.T. Vvedenie v optimizatsiiu. —Moscow: Nauka, 1983. - 384 s. (rus)
7. Rakitskii Iu.V. Chislennye metody resheniia zhestkikh sistem. -Moscow: Nauka, 1979. - 208 s. (rus)
8. Chernorutskii I.G. Funktsii relaksatsii gradi-entnykh metodov / Nauchno-tekhnicheskie vedomosti SPbGPU. Informatika. Telekommunikatsii. Upravlenie. —St.-Petersburg: Izd-vo Politehn. unta, 2012. -№ 3 (150). -S. 66-72. (rus)
9. Chernorutskii I.G. Nekotorye standartnye skhemy parametricheskoi optimizatsii / Nauch-no-tekhnicheskie vedomosti SPbGPU. Informa-tika. Telekommunikatsii. Upravlenie. -St.-Peters-burg: Izd-vo Politehn. un-ta, 2012. -№ 6 (162). -S. 128-133. (rus)
ЧЕРНОРУЦКИЙ Игорь Георгиевич — директор Института информационных технологий и управления, заведующий кафедрой информационных и управляющих систем Санкт-Петербургского государственного политехнического университета, доктор технических наук, профессор. 195251, Россия, Санкт-Петербург, Политехническая ул., д. 21. E-mail: director@icc.spbstu.ru
CHERNORUTSKIY, Igor G. St. Petersburg State Polytechnical University. 195251, Politekhnicheskaya Str. 21, St.-Petersburg, Russia. E-mail: director@icc.spbstu.ru
© Санкт-Петербургский государственный политехнический университет, 2013
