Software & Systems

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

Рассмотрим задачу безусловной минимизации

,                                                           (1)

где F: Rn®R1 – гладкая функция; Rn – n-мерное евклидово пространство.

Пусть hk, Hk – градиент и гессиан, вычисленные на итерации k в точке xk процесса решения задачи (1). Тогда общий принцип построения большинства ньютоновских методов безусловной оптимизации с регулировкой шага состоит в следующем. На каждой итерации сначала строится некоторая связанная с Hk существенно положительно определенная матрица  а затем направление спуска pk вычисляется как решение системы  pk = –hk и определяется новая точка xk+1 = xk + ak pk, где ak – длина шага, для которого Fk+1 < Fk. При этом процедуру построения  организуют так, чтобы  совпадала с исходной матрицей Hk, если последняя сама является положительно определенной, причем выяснение определенности Hk и построение  осуществляются параллельно в рамках одной процедуры на основе некоторых матричных разложений, которые позволяют выявить знаки собственных чисел Hk и приспособиться для генера- ции

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

В данной работе в основе выяснения определенности Hk и построения  лежит факторизация (разложение) Холесского. Разложение Холесского для симметричной положительно определенной матрицы Hk имеет вид Hk = LkDk(Lk)T = (Uk)T DkUk, где Uk =(Lk)T – верхняя треугольная матрица с диагональными элементами, равными единице; Dk – диагональная матрица с положительными ведущими элементами по диагонали. Поскольку столбцы матрицы Lk с номерами с 1-го по (j–1)-й известны, ее j-й столбец определяется по формулам  

Аналогичные формулы существуют и для построчной организации вычислений. В отличие от гауссовых исключений алгоритм Холесского численно устойчив без каких-либо перестановок [3]. Это свойство определяется соотношением  j = 1, …, n, между элементами Hk и Lk. В силу него существует априорное ограничение сверху на элементы : каждый из них не превосходит максимальную величину . Соответственно лавинообразный рост элементов  невозможен независимо от того, будут ведущие элементы малыми или нет.

Построение направлений ньютоновского поиска. Факторизация Холесского неприменима для любой симметричной матрицы по следующим причинам. Во-первых, для знаконеопределенной матрицы Hk факторизация Холесского может не существовать. Во-вторых, даже если она и существует, то гарантировать численную устойчивость алгоритма уже нельзя, поскольку никаких априорных ограничений на субдиагональные элементы Lk в рассматриваемом случае не будет. Итак, на основе обычной факторизации Холесского построение  не получить, для этого нужно воспользоваться модифицированной факторизацией. Наиболее эффективная вычислительная схема предложена в [2, 3]. Данный подход состоит в том, чтобы строить факторы Холесского Dk и Lk, подчиняющиеся двум требованиям: все диагональные элементы Dk должны быть существенно положительными, модули всех элементов треугольного фактора Lk должны быть равномерно ограничены сверху. Точнее говоря, требуется, чтобы для всех j=1, 2, …, n и некоторых заданных положительных d и b выполнялись неравенства

                                         (2)

где  – введенные для удобства изложения вспомогательные величины, по определению равные . Гринштадт предложил в [4] при реализации алгоритма на конкретной ЭВМ, в которой под запись мантиссы отводится t битов, величину d вычислять по формуле .

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

Здесь x – максимальный модуль недиагональ- ного элемента Hk; g – значение максимального из диагональных элементов Hk; eM – машинная точность; eM вводится в формулу расчета b, чтобы обеспечить устойчивость вычислений, когда норма Hk очень мала [4]. При этом процедура расчета модифицированных факторов Lk и Dk фактически представляет собой обычный алгоритм факторизации Холесского с попутным увеличением (по мере необходимости) диагонали исходной матрицы с целью добиться выполнения неравенств (2). Матрицы Lk и Dk, полученные по окончании описанной процедуры, будут факторами Холесского для положительно определенной матрицы  связанной с Hk следующим образом: LkDk(Lk)T = Hk + Ek + , где Ek – неотрицательная диагональная матрица, j-й элемент которой равен . Таким образом, положительно определенная матрица  может отличаться от исходной матрицы Hk  только диагональными элементами.

Направление спуска pk вычисляется как решение системы , причем  определяют последовательным решением двух систем линейных уравнений с треугольными матрицами: Lk yk = –hk, Dk(Lk)T pk = yk.

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

Вычислительная схема 1 модифицированного разложения Холесского. Далее приводится детальное описание всех операций, выполняемых по ходу построения модифицированного разложения Холесского с перестановками [3]. Дана наиболее эффективная схема организации расчетов. В процессе вычисления j-го столбца матрицы Lk участвуют вспомогательные величины   s = 1, …, j, i = j, …, n.

Шаг 0 (Расчет порога для элементов). Вычислить , где , а числа g и x – максимальные значения модулей диагонального и недиагонального элементов Hk.

Шаг 1 (Инициализация). Присвоить индексу столбца j значение 1. Принять  i = 1, …, n.

Шаг 2 (Перестановка строк и столбцов). Найти индекс q, такой, что . Поменять местами все данные, отвечающие столбцам матрицы Hk с номерами q и j, а затем проделать то же самое с данными, отвечающими ее q-й и j-й строкам.

Шаг 3 (Поиск максимальной по модулю вели- чины ). Вычислить  для i = j+1, …, n и найти  (если j = n, взять qj = 0).

Шаг 4 (Расчет j-го диагонального элемента фактора Dk). Вычислить  и поправку . Если j = n, вычисления прекратить.

Шаг 5 (Расчет j-й строки Lk). Вычислить  для s = 1, …, j – 1. Для i = j + 1, …, n пересчитать диагональные элементы  по формуле . Присвоить j значение j + 1 и вернуться к шагу 2.

Модифицированная факторизация Холесского, помимо прочего, позволяет определить и направление отрицательной кривизны, решая уравнение (Lk)T pk = es, где индекс s выбирается из условия  j = 1, …, n.

Подход к увеличению эффективности алгоритма Гилла и Мюррея предложен в [1]. Доказано, что подход к построению направления спуска pk определяет решение проблемы масштабирования шагов при спуске и интеграцию с методом доверительной окрестности следующим образом. В алгоритме Гилла и Мюррея направление спуска pk в конечном счете вычисляется как решение системы

(Lk)T pk = Uk pk = uk,                                           (3)

где uk = (Dk)–1 yk, а yk – это решение системы Lk yk = –hk. Величина элементов  зависит от способа задания направления спуска. Следуя подходу Гилла и Мюррея, для численной устойчивости расчета элементов Dk, Lk, uk достаточно потребовать изменения способа задания pk так, чтобы для всех j = 1, 2, …, n и некоторых заданных положительных d, b и w выполнялись неравенства

.

Это может стать ключом к решению проблемы масштабирования шагов при спуске, но такое предположение ошибочно: подход линейной алгебры к вычислению направления спуска pk исключает расчет элементов вектора uk по ходу построения мо- дифицированного разложения Холесского [5, 6]. Нужен альтернативный подход, который опишем следующим образом.

Равенства Hk = (Uk)T Dk Uk достаточно для определения элементов матриц Dk, Uk. В скалярной форме это равенство выглядит следующим образом: . Отсюда, полагая , построим соотношения

                     (4)

для расчета элементов Dk, Uk. Построение факторов Холесского математически эквивалентно применению метода исключения Гаусса к системе уравнений Hk pk = –hk в прямом порядке, при этом техника исключения Гаусса позволяет получить соотношение                          (5)

для расчета элементов вектора uk и систему уравнений (3) для вычисления направления спуска pk. Интеграция техники исключения Гаусса и факторизации Холесского приводит к соотношениям (4) и (5), которые с помощью элементарных арифметических операций позволяют вычислить элементы uk по ходу построения факторов Холесского.

Интеграция техники исключения Гаусса и факторизации Холесского генерирует множество реализаций численно устойчивых вычислительных схем задания направления спуска pk, в основе которых лежит требование, чтобы на очередном шаге вычислений коэффициентов факторов Холесского соответствующий диагональный элемент матрицы Dk и соответствующий элемент вектора uk сначала рассчитывались по вычисленным ранее значениям этих коэффициентов. Затем диагональный элемент Dk увеличивается настолько, насколько необходимо, чтобы все диагональные элементы Dk были существенно положительными, модули всех элементов Uk и uk были равномерно ограничены сверху. Например, достаточно потребовать, чтобы для всех i = 1, 2, …, n выполнялись неравенства

,                           (6)

положив , посколь­ку элементы ,  и  определяют i-е строки матриц Dk, Uk, uk и удовлетворяют неравенствам (6). Здесь , ,  – вспомогательные величины, определяемые следующим образом:

                         (7)

Интеграция техники исключения Гаусса и факторизации Холесского генерирует множество реализаций вычислительных схем интеграции с методом доверительной окрестности, например, следующим образом. В процессе построения Dk, Uk, uk вычислить элементы  (если , взять ), где , вычислить  и определить направление спуска pk из системы

                                                   (8)

Здесь параметр g задается пользователем, число gk (gk ³1 по построению) характеризует степень однородности, при этом (8) – простейшая форма аппроксимации, отличная от квадратичной. Пусть , где Gk – диагональная матрица с элементами  на диагонали. Тогда , что означает сдвиг на  всех i-х собственных значений матрицы  В методах с регулировкой шага матрицу Hk модифицируют так, чтобы изменения не затрагивали подпространства, натянутого на ее собственные векторы с положительными собственными значениями [3]. Если же замена осуществляется в методе доверительной окрестности, она отражается на всех векторах, так как результатом замены Hk на матрицу = Hk + lk I будет сдвиг на lk всех собственных значений матрицы Hk [3]. Из вышесказанного следует, что стратегия выбора направления спуска определяет и интеграцию с методом доверительной окрестности.

Вычислительная схема 2 модифицированного разложения Холесского. Далее приведена наиболее эффективная схема организации расчетов. Все фигурирующие в ней величины при реализации на ЭВМ могут размещаться в памяти, первоначально выделяемой для записи матриц Hk и hk. При этом коэффициенты рассчитываемых факторов занимают места ее использованных элементов. В процессе вычисления i-й строки фактора Uk участвуют вспомогательные величины (7). В процессе построения системы уравнений Uk pk = uk на каждой итерации k выявляются знаки собственных чисел и вычисляется количество отрицательных nо и нулевых n0 собственных значений, число переходов к направлению отрицательной кривизны nS.

Шаг 0. Вычислить , где , а числа g0 и x – максимальные значения модулей диагонального и недиагонального элементов Hk. Положить , . Здесь tF – число правильных разрядов F(xk), которые хотелось бы получить.

Шаг 1. Принять i = 1 и  (j = 1, …, n),  (j = 1, …, n).

Шаг 2. Найти такой индекс q, что , и поменять местами все данные, отвечающие строкам матрицы Hk с номерами q и i, а затем проделать то же самое с данными, отвечаю- щими ее q-му и i-му столбцам.

Шаг 3. Если , то положить  иначе положить .

Шаг 4. Если , то положить n0 = n0+1 и перейти к шагу 8; иначе перейти к шагу 5.

Шаг 5. Если , то перейти к шагу 8; иначе перейти к шагу 6.

Шаг 6. Принять n0 = n0 + 1. Если , то перейти к шагу 8; иначе перейти к шагу 7.

Шаг 7. Построить систему уравнений Uk pk = ei для поиска направления отрицательной кривизны по формулам  j = 1, …, n, , положить nS = nS + 1 и остановиться.

Шаг 8. Если i = n, то положить qi = 0; иначе вычислить вспомогательные величины  j = i + 1, …, n, по формуле  и найти .

Шаг 9. Вычислить диагональный элемент фактора Dk и элемент вектора uk по формулам , . Вычислить  (если , взять ). Если i = n, то перейти к шагу 11.

Шаг 10. Вычислить строку  j = i+1, …, n, фактора Uk, пересчитать вспомогательные величины по формулам  j = i + 1, …, n,  j = i + 1, …, n, положить i = i + 1 и перейти к шагу 2.

Шаг 11. Найти , положить uk = gk uk и остановиться.

Замечания. Предложенный подход к вычислению длины шага можно найти в [7–9]. В основе лежит интеграция подходов [3, 10–12].

Предложенный подход к модификации критериев останова, разработанных Гиллом, Мюрреем и Райт [3], можно найти в [6].

Предложенный подход к оцениванию конечно-разностных интервалов и подход к конечно-разностной аппроксимации первых и вторых производных можно найти в [1]. В основе лежит интеграция подходов [3, 11], там же описаны способы выбора конечно-разностных интервалов. Стандартные программы численного дифференцирования нацелены на поиск конечно-разностного интервала, обеспечивающего минимальную суммарную ошибку вычисляемого приближения. Алгоритмы выбора интервала для центральной конечно-разностной формулы можно найти, например, в [13]. Однако для оптимизационных приложений такие программы не подходят, так как, обеспечивая избыточную точность вычисления оценок производных, они обычно требуют слишком большого числа обращений к процедуре расчета значений функции.

В [14] отмечено, что замена градиентов их конечно-разностными оценками лучше всего проходит в квазиньютоновских методах. Для методов ньютоновского типа результаты получаются менее удовлетворительными.

Тестирование

Построение алгоритмов, в том числе и развитие подхода Гилла и Мюррея к построению направления спуска [1, 6, 9, 14], требует постоянного тестирования на больших классах тестовых задач различных структур и характеристик с различными настройками с целью выявления слабых мест – ситуаций, приводящих к замедлению алгоритма или ухудшению качества результата и т.д.

Алгоритмы безусловной оптимизации должны быть протестированы по крайней мере в двух различных смыслах:

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

-     для получения представления о гипотезе, доказательства работоспособности алгоритма и сравнения алгоритмов на экспериментальном уровне.

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

Реальные задачи берутся из разных источников прикладных оптимизационных задач физики, химии, техники, биологии, экономики, океанологии, астрономии, метеорологии и т.д. В отличие от искусственных настройки реальных задач труднодоступны и ими нелегко оперировать. Они могут иметь сложные алгебраические или дифференциальные выражения, могут зависеть от огромного количества данных и, возможно, зависят от некоторых параметров, которые должны быть оценены определенным образом. Очень хорошая коллекция реальных задач безусловной оптимизации описана в работах [15, 16]. Следует отметить, что собрание коллекций реальных задач началось в конце 1970-х годов с отстаивания высокого качества отчетности о вычислительных экспериментах математического программирования с ПО [17–19].

Источником тестовых задач, рассмотренных в данной работе, является коллекция [20], в состав которой вошли тестовые функции из коллекции CUTE и прочих [21–24], а также из других работ и технических отчетов. Цель этой коллекции – представить большое количество общих тестовых функций, которые могут использоваться для тестирования алгоритмов безусловной оптимизации и сравнения исследований. Для каждой функции даны алгебраические представления и начальные (стартовые) точки. Для всех задач из коллекции [20] известные алгоритмы безусловной оптимизации, например, квазиньютоновские, сходятся к одному и тому же локальному решению.

Опишем конкретные тестовые задачи [20].

Задача 1 (Extended Rosenbrock function):

Задача 2 (Wood function):

Задача 3 (Generalized Rosenbrock function):

Задача 4 (Extended White & Holst function):

Задача 5 (Extended Penalty function):

Задача 6 (Perturbed Quadratic function):

Задача 7 (Raydan 1 function):

Задача 8 (Raydan 2 function):

Задача 9 (Diagonal 1 function):

Задача 10 (Diagonal 2 function):

Задача 11 (Diagonal 3 function):

Завершение итераций производилось с использованием рекомендаций из [6]:

,

,

,

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

Замечания. Трудности сравнения нескольких методов оптимизации по результатам численных экспериментов состоят в следующем.

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

Во-вторых, неясно, каким образом соизмерять трудоемкость различных методов. Естественный, на первый взгляд, критерий – затраченное машинное время – в действительности не подходит: сложно сравнивать быстродействие различных машин, не всегда имеются данные о затратах времени при работе в мультипроцессорном режиме и т.д. Более надежным показателем является число вычислений минимизируемой функции или другая внутренняя характеристика метода. Но и здесь возникает ряд проблем. Например, неясно, как соизмерять трудоемкость вычисления функции и решения различных вспомогательных задач. Затруднительно сопоставлять вычисления функции и ее производных. В тех задачах, где применяются конечно-разностные аппроксимации производных, все намного проще: одно вычисление градиента эквивалентно n вычислениям функции, вычисление гессиана эквивалентно n(n+1)/2 вычислениям функции. Однако для дискретной задачи оптимального управления градиент лишь примерно вдвое «дороже» функции [25]. Для квадратичной функции градиент вычисляется даже проще, чем ее значения. Так же дело обстоит и с субградиентом в минимаксной задаче и т.д.

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

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

-     приводить точную формулировку задачи, для которой проводился счет, включая все ее параметры и начальное приближение;

-     указывать тип ЭВМ и данных, язык и среду программирования, сведения о программе;

-     давать подробное описание применяемого алгоритма или отсылать к его публикации, если она имеется;

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

-     указывать подробные сведения о трудоемкости вычислений.

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

Далее x* – точка глобального минимума F(x) на Rn, F* = F(x*); k – число итераций;  – точность решения по функции;  – точность решения по аргументам; k0 – число вычислений функции; k1 – число вычислений функции (регулировка шага); ÑF(x) – градиент F(x); Ñ2F(x) – матрица вторых производных (гессиан).

Тестовые задачи с известным решением можно построить следующим образом. Выбираются точка x* и функции ji(x), i = 1, 2, …, m. Тогда

                         (9)

достигает минимума в точке x*. При этом F*= 0, , так что можно регулировать обусловленность задачи за счет выбора ji(x). В частности, если m < n или m = n, но Ñji(x*) = 0 для некоторых i, получаем вырожденную точку минимума. Функции вида (9), вообще говоря, невыпуклы и могут иметь локальные и даже глобальные минимумы, отличные от x*.

Другой способ связан с выбором функций вида

.           (10)

Тогда, задавшись  и решив последовательно одномерные уравнения

     (11)

(решение обозначим ), для i = 0, 1, …, n выберем . Очевидно, что тогда ÑF(x*) = 0. Если fi(xi, xi–1) выпуклы по xi, xi–1, то F(x) выпукла и x* – точка минимума F(x).

Наконец, можно взять произвольную гладкую выпуклую функцию j(x) и произвольную точку x* и построить . Тогда ÑF(x*) = 0, F(x) выпукла и потому x* – точка глобального минимума F(x). Существуют и другие приемы построения функций с известной точкой минимума.

Следует отметить, что для функций вида (9), (10) существуют специальные методы, обладающие повышенной эффективностью. Например, для задач (11) эффективен метод Гаусса–Ньютона [3], для (11) – метод динамического программирования и покоординатного спуска [25]. Аналогичным образом специфика задач (9), (11) может сказаться на поведении общих методов оптимизации.

Развитие подхода и численные исследования

В [3] доказано, что значение нормы поправки Ek можно дополнительно уменьшить, если использовать симметричные перестановки строк и столбцов Hk  на шаге 2 вычислительной схемы 1 модифицированного разложения Холесского следующим образом. Найти такой индекс q, что , и поменять местами все данные, отвечающие строкам матрицы Hk с номерами q и i, а затем проделать то же самое с данными, отвечающими ее q-му и i-му столбцам.

Подход Гилла и Мюррея к дополнительному уменьшению нормы поправки Ek развит в [6]. Доказано, что фактическое значение нормы поправки Ek можно дополнительно уменьшить, если использовать симметричные перестановки строк и столбцов Hk. На очередном, i-м, шаге факторизации в качестве i-й строки и i-го столбца надо брать ту из нетронутых пар, для которой величина  максимальна.

Такая стратегия дает возможность увеличить численную устойчивость расчета элементов Dk, Uk, uk и остается работоспособной и в том случае, когда  не равна нулю, но очень мала. Там же отмечено: «Параметр b подбирается путем минимизации априорной оценки нормы поправки Ek и с очевидностью зависит от способа задания направления спуска и, как следствие, от величин  и  Рассмотрение предлагаемой формулы вычисления величины b, учитывающей такую зависимость, не является целью данной работы и находится в стадии экспериментального сравнения эффективности». В данной работе величину b предлагается вычислять по формуле

,                                (12)

где , eM – машинная точность; eM вводится в формулу расчета b, чтобы обеспечить устойчивость вычислений, когда норма Hk очень мала. Применение формулы (12) хорошо согласу- ется с численными экспериментами, описанными в данной работе. Ниже приводится предлагаемая модификация вычислительной схемы 2, дана наиболее эффективная схема организации расчетов.

Вычислительная схема 3 модифицированного разложения Холесского. Все фигурирующие в ней величины при реализации на ЭВМ могут размещаться в памяти, первоначально выделяемой для записи матриц Hk и hk.

Шаг 0. Вычислить , где , , а числа g0 и x – максимальные значения модулей диаго- нального и недиагонального элементов Hk, eM – машинная точность. Положить , . Здесь tF – число правильных разрядов F(xk), которые хотелось бы получить.

Шаг 1. Положить i = 1 и  (j = 1, …, n),  (j = 1, …, n).

Шаг 2. Найти индекс q, такой, что , и поменять местами все данные, отвечающие строкам матрицы Hk с номерами q и i, а затем проделать то же самое с данными, отвечающими ее q-му и i-му столбцам.

Шаг 3. Если , то положить  иначе положить .

Шаг 4. Если , то положить n0 = n0 + 1 и перейти к шагу 8; иначе перейти к шагу 5.

Шаг 5. Если , то перейти к шагу 8; иначе перейти к шагу 6.

Шаг 6. Положить n0 = n0 + 1. Если , перейти к шагу 8; иначе перейти к шагу 7.

Шаг 7. Построить систему уравнений Ukpk = ei для поиска направления отрицательной кривизны по формулам  j = 1, …, n, , положить ns = ns + 1 и остановиться.

Шаг 8. Если i = n, то положить qi = 0; иначе вычислить вспомогательные величины  j = 1, …, n, по формуле  и найти .

Шаг 9. Вычислить диагональный элемент фактора Dk и элемент вектора uk по формулам , . Вы- числить  (если , взять ). Если i = n, перейти к шагу 11.

Шаг 10. Вычислить строку  j = 1, …, n, фактора Uk, пересчитать вспомогательные величины по формулам  j = i+1, …, n,  j = i + 1, …, n, положить i = i + 1 и перейти к шагу 2.

Шаг 11. Найти , положить uk = gkuk и остановиться.

В данной работе отмечена взаимосвязь подхода к увеличению эффективности гауссова исключения для разреженных матриц в [26] и предлагаемого подхода к увеличению эффективности численных методов ньютоновского и квазиньютоновского типов – использование структуры матрицы, то есть информации о том, в каких позициях матрицы хранятся ненулевые элементы. Для ньютоновских и квазиньютоновских методов безусловной оптимизации, основанных на факторизации Холесского, с регулировкой шага и с конечно-разностной аппроксимацией первых и вторых производных – это возможность уменьшения числа вычислений функции путем формирования матрицы вторых производных в соответствии с ее структурой. Данное исследование является прямым продолжением работы [14].

Программные реализации. Приложения Nmbm и NmbmApp – машинная реализация ньютоновских методов оптимизации с регулировкой шага, основанных на факторизации Холесского (вычислительная схема 2 модифицированного разложения Холесского) [27, 28].

NmbmApp отличается от своего классического прототипа Nmbm конечно-разностной аппроксимацией первых и вторых производных. Все версии алгоритмов реализованы на языке Visual Basic .NET, среда разработки – Microsoft Visual Studio 2010. Исходные тексты программ доступны на сайте http://crm.ics.org.ru и могут быть использованы по универсальной общедоступной лицензии GNU.

Приложение NmbmsApp отличается от своего классического прототипа NmbmApp формированием матрицы вторых производных в соответствии с ее структурой. NmbmsApp – машинная реализация ньютоновских методов оптимизации с регулировкой шага, основанных на факторизации Холесского (вычислительная схема 3 модифицированного разложения Холесского).

Главное окно программы и пример результатов решения одной из тестовых задач размещены по ссылке http://www.swsys.ru/uploaded/image/2017_2/ 2017-2-dop/6.jpg/.

Достоинством подхода, основанного на применении формул численного дифференцирования, кроме его универсальности, является низкая стоимость подготовки задачи к компьютерному моделированию [29]. От пользователя требуется лишь написание программы для вычисления значения F(x) при заданном x, поэтому реализованные на ос- нове численных производных методы оптимиза- ции оказываются по существу прямыми методами, то есть методами, не использующими в своей схеме производные от F(x). В то время как результаты численных исследований, приведенные в [14], показывают, что конечно-разностные аналоги проигрывают в точности решения следующих задач оптимизации.

Задача 12 (функция Поляка, среднеквадратичная аппроксимация экспонентами):

,

.

Функция невыпуклая, с искривленным оврагом, обусловленность велика.

Задача 13 (функция Вуда):

Функция имеет несколько локальных минимумов, это обстоятельство может вызвать преждевременное окончание процесса.

Задача 14 (функция Степенная):

,

.

В таблице 1 приводятся результаты численных исследований, выполненных в трех версиях предложенного ньютоновского метода безусловной минимизации: Nmbm, NmbmApp [6, 9, 14, 27, 28] и предлагаемой NmbmsApp.

Таблица 1

Результаты численных исследований задач 12, 13, 14

Table 1

The results of numerical studies for problems 12, 13, 14

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

Долан и Морэ впервые предложили применение бенчмаркинга в качестве нового инструмента использования статистических данных для сравнения итерационных алгоритмов путем демонстрации профилей производительности [30]. Термин «бенч­маркинг» англоязычный и не имеет дословного перевода на русский язык, происходит от английских слов «bench» (уровень, высота) и «mark» (отметка). Это словосочетание трактуется по-разному: «опорная отметка», «отметка высоты», «эталонное сравнение» и т.д. Как правило, термином «бенчмаркинг» обозначается один из инструментов совершенствования деятельности.

Профилем производительности Ps для метода решения оптимизационной задачи называется функция распределения какого-либо измеримого показателя производительности [30]. Вычисление профилей производительности позволяет визуализировать различия по эффективности нескольких оптимизационных методов, причем для их построения используется отдельная программа [30].

Функция Ps(t) определяется следующим образом:

Здесь S – множество сравниваемых методов, P – множество решаемых с помощью этих методов задач. Количество элементов в P обозначено через np, соответственно ns – количество элементов в S. В качестве измеримых показателей производительности tp,s, как правило, используются число итераций и число вычислений значений функции (или число вычислений значений функции и оценок градиента) [31–33].

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

Бенчмаркинг, как и любой другой инструмент для сравнения итерационных алгоритмов, имеет свои достоинства и недостатки. Чувствительность только к измеримому показателю производительности tp,s и отсутствие удовлетворительного способа сравнения между собой нескольких методов оптимизации по результатам численных экспериментов [25] часто являются причинами разногласий при интерпретации результатов [30]. Кроме того, вычисление профилей производительности требует отбрасывания задач, для которых машин- ные реализации соответствующих методов потер- пели неудачу [30], что усложняет процесс сравне- ния итерационных алгоритмов.

В данной работе, как и в [31–33], все тестовые задачи описаны в [20], все сравниваемые алгоритмы безусловной оптимизации сходятся к одному и тому же локальному решению, сравнению подлежат подходы к заданию направления спуска, результаты решения заносятся в таблицы.

Следует отметить: если объем статистических данных небольшой, то для сравнения итерационных алгоритмов, как правило, достаточно таблиц [32, 33]; иначе для сравнения итерационных алгоритмов используется бенчмаркинг путем демонстрации профилей производительности [30, 31].

В [31–33] завершение итераций производится с использованием одного из следующих критериев останова:

          (13)

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

Таблица 2

Результаты численных исследований, соответствующие критериям (13)

Table 2

The results of numerical studies that meet the criteria (13)

Таблица 3

Результаты численных исследований с максимальной точностью

Table 3

The results of numerical studies with maximum precision

Заключение

В работе доказано, что интеграция техники исключения Гаусса и факторизации Холесского определяет решение проблемы масштабирования шагов при спуске, а следовательно, и аппроксимацию неквадратичными функциями, и интеграцию с методом доверительной окрестности [1], а также то, что интеграция техники исключения Гаусса и факторизации Холесского определяет и подходы к дальнейшему уменьшению значения нормы поправки Ek, если использовать симметричные перестановки строк и столбцов Hk следующим образом: на очередном, i-м, шаге факторизации в качестве i-й строки и i-го столбца надо брать ту из нетронутых пар, для которой величина  максималь- на [6].

Кроме того, доказано, что модифицированная факторизация Холесского – оптимизированный алгоритм в том смысле, что параметр b подбирается в нем путем минимизации априорной оценки нормы поправки Ek при условии сохранения существенно положительно определенной матрицы неизменной и с очевидностью зависит от способа задания направления спуска [1, 6]. В данной работе доказанные утверждения позволили модифициро- вать формулу вычисления величины b, предложен- ную в [2], и построить вычислительную схему 3. Полученные результаты хорошо согласуются с численными экспериментами, приведенными в данной работе.

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

Рассмотрена взаимосвязь подхода к увеличению эффективности гауссова исключения для разреженных матриц [26] и предлагаемого подхода к увеличению эффективности численных методов ньютоновского типа – использование структуры матрицы, то есть информации о том, в каких позициях матрицы хранятся ненулевые элементы. Для ньютоновских методов безусловной оптимизации, основанных на факторизации Холесского, с регулировкой шага и с конечно-разностной аппроксимацией первых и вторых производных, это возможность уменьшения числа вычислений функции путем формирования матрицы вторых производных в соответствии с ее структурой. Данное исследование является прямым продолжением работы [14].

Литература

1.     Зеленков Г.А., Хакимова А.Б. Подход к разработке алгоритмов ньютоновских методов оптимизации, программная реализация и сравнение эффективности // Компьютерные исследования и моделирование. 2013. Т. 5. № 3. С. 367–377.

2.     Gill P.E. and Murray W. Newton-type methods for unconstrained and linearly constrained optimization. Math. Prog. 1974, no. 28, pp. 311–350.

3.     Гилл Ф., Мюррей У., Райт М. Практическая оптимизация; [пер. с англ.]. М.: Мир, 1985. 509 с.

4.     Гилл Ф., Мюррей У. Численные методы условной оптимизации; [пер. с англ.]. М.: Мир, 1977. 290 с.

5.     Парлетт Б. Симметричная проблема собственных значений. Численные методы; [пер. с англ.]. М.: Мир, 1983. 382 с.

6.     Свириденко А.Б. Априорная поправка в ньютоновских методах оптимизации // Компьютерные исследования и моделирование. 2015. Т. 7. № 4. C. 835–863.

7.     Хакимова А.Б., Зеленков Г.А. Увеличение эффективности ньютоновских методов минимизации. Вычисление длины шага // Динамика неоднородных систем: тр. ИСА РАН. 2010. Вып. 14-А. Т. 53. С. 115–120.

8.     Хакимова А.Б., Дикусар В.В., Зеленков Г.А. Увеличение эффективности ньютоновских методов оптимизации. Информодинамический подход // Динамика неоднородных систем: тр. ИСА РАН. 2010. Вып. 14-А. Т. 53. С. 97–114.

9.     Свириденко А.Б., Зеленков Г.А. Управление процессом построения ньютоновских методов безусловной оптимизации: монография. Новороссийск: Изд-во ГМУ им. адм. Ф.Ф. Ушакова, 2016. 98 с.

10.   Пшеничный Б.Н., Данилин Ю.М. Численные методы в экстремальных задачах. М.: Наука, 1975. 320 с.

11.   Полак Э. Численные методы оптимизации. Единый подход. М.: Мир, 1974. 376 с. 12.   Карманов В.Г. Математическое программирование. М.: Наука, 1975. 272 с.

13.   Stepleman R.S., Winarsky N.D. Adaptive numerical differentiation. Mathematics of Computation, 1979, no. 33, pp. 1257–1264.

14.   Свириденко А.Б., Зеленков Г.А. Взаимосвязь и реализа- ция квазиньютоновских и ньютоновских методов безусловной оптимизации // Компьютерные исследования и моделирование. 2016. Т. 8. №. 1. С. 55–78.

15.   Averick B.M., Carter R.G., Moré J.J. The MINPACK-2 test problem collection (Preliminary version). Math. and Comp. Sc. Divis., 1991, no. 150.

16.   Averick B.M., Carter R.G., Moré J.J., Xue G.L. The MINPACK-2 test problem collection. Math. and Comp. Sc. Divis., 1992.

17.   Jackson R.H.F., and Mulvey J.M. A Critical review of comparisons of mathematical programming algorithms and software (1953–1977). J. Research of the National Bureau of Standards, 1978, no. 83, pp. 563–584.

18.   Gilsinn J., Hoffman K., Jackson R.H.F., Leyendecker E., Saunders P., and Shier D. Methodology and analysis for comparing discrete linear L1 approximation codes. Communications in Statistics, 1977, no. 136, pp. 399–413.

19.   Crowder H.P., Dembo R.S., and Mulvey J.M. On reporting computational experiments with mathematical software. ACM Transactions on Mathematical Software, 1980, no. 5, pp. 193–203.

20.   Andrei N. An unconstrained optimization test functions collection. Adv. Model. Optim. 2008, vol. 1, no. 10, pp. 147–161.

21.   Bongartz I., Conn A.R., Gould N.I.M., Toint P.L., CUTE: constrained and unconstrained testing environments. ACM TOMS, 1995, vol. 21, pp. 123–160.

22.   Moré J.J., Garbow B.S., Hillstrom K.E. Testing unconstrained optimization software. ACM Trans. Math. Soft., 1981, vol. 7, pp. 17–41. 23.   Himmelblau D., Applied nonlinear programming. McGraw-Hill, NY, 1972.

24.   Lee D. A fast and robust unconstrained optimization method requiring minimum storage. Mathematical Programming. 1985, vol. 32, pp. 41–68.

25.   Поляк Б.Т. Введение в оптимизацию. М.: Наука, 1983. 384 с.

26.   Соловьев C.A. Решение разреженных систем линейных уравнений методом Гаусса c использованием техники аппроксимации матрицами малого ранга // Вычислительные методы и программирование. 2014. Т. 15. С. 441–460.

27.   Свириденко А.Б. MNB (Ньютоновский метод безусловной оптимизации). Свид. о гос. регистр. прогр. для ЭВМ № 2015610399. 2015.

28.   Свириденко А.Б., Зеленков Г.А. MNBApp (Ньютоновский метод безусловной оптимизации с численным вычислением первых и вторых производных). Свид. о гос. регистр. прогр. для ЭВМ № 2015610347. 2015.

29.   Черноруцкий И.Г. Методы оптимизации. Компьютерные технологии. СПб: БХВ-Петербург, 2011. 384 с.

30.   Dolan E., Moré J.J. Benchmarking optimization software with performance profiles. Math. Program., 2002, no. 91, pp. 201–213.

31.   Ahookhosh M., Amini K. An efficient nonmonotone trust-region method for unconstrained optimization. Numer Algor, 2012, no. 59, pp. 523–540.

32.   Runak M. Abdula, Abbas Y. Al Bayati. A new hybrid scaled search direction for unconstrained optimization. Jour. of Kirkuk University, 2007, vol. 2, no. 1, pp. 62–81.

33.   Abbo K.K. New CG Method for large-scale unconstrained optimization based on Nazareth theorem. Iraqi Jour. of Statistical Sciences, 2008, vol. 13, pp. 53–65.