Модификация метода Ньютона безусловной минимизации функций повышенной сходимости Текст научной статьи по специальности «Математика»

УДК 518.79

МОДИФИКАЦИЯ МЕТОДА НЬЮТОНА БЕЗУСЛОВНОЙ МИНИМИЗАЦИИ ФУНКЦИЙ ПОВЫШЕННОЙ СХОДИМОСТИ

ЕВДОКИМОВ А.Г., МАНАКОВА Н. О.

Рассматривается модификация метода Ньютона безусловной минимизации функций повышенной сходимости, его описание, алгоритм, реализация в пакете компьютерной алгебры Maple Power Edition, а также сравнительный анализ сходимости метода Ньютона и рассматриваемой здесь модификации на классических тестовых функциях Химмельблау [1] и тестовой функции Евдокимова [2].

Большинство численных методов решения задачи безусловной минимизации, получивших распространение на практике, относится к прямым методам. Основное рекуррентное соотношение, объединяющее эти методы и позволяющее по некоторой информации на k-м шаге определить местоположение на (k + 1)-м, можно представить в виде

xM =хМ +х МдхМ. (1)

Здесь ДхМ — направляющий вектор на k-м шаге,

определяющий направление движение; X^—параметр, характеризующий длину шага в выбранном направлении. В частности, для классического метода Ньютона ^(k) = 1, а

лХ«=-(н-.),к,(|)(к) (2)

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

( . . fk (к)Л

У

« Ук)(н -f (§ у

eR1 ,

(3)

где H(k) — матрица Г есса (матрица вторых частных

ffy) ^

производных); I —= I — вектор первых частных

\д x )

производных.

Введение в ньютоновскую процедуру параметра

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

функции у (х). Однако и обобщенному методу Ньютона присущи два серьезных недостатка:

1) если матрица вторых производных H(k) не является положительно-определенной, вычислительный процесс (3) может привести не к уменьшению

у 1x1, а к увеличению; в результате будет получено = 0 и процесс прервется в точке Х^ ;

2) в точке хОО вообще может не существовать матрицы, обратной Ик\

Фиакко и Мак-Кормик [3] предложили модификацию обобщенного метода Ньютона, дающую возможность обойти названные трудности. Попробуем дать несколько иную трактовку этой модификации, используя идею приведения матрицы Hk к диагональному виду.

Пусть матрица Н(к) в выражении (2) — неособенная

матрица. Перейдем к новой системе координат Д£,

в которой матрица В(к) квадратичной формы (2) будет диагональной, т.е. получим такую матрицу преобразований Р(к) и диагональной матрицы В(к) в результате применения к матрице Hk обыкновенных гауссовых исключений [4]. Находим

Д^(к) = Р(к)дХ(к); (4)

ДХ(к) =(р(к))-1Д?(к). (5)

Получаем значение функции на k-м шаге:

Ay W = ^Ду (к)

j=1 j=1

f л \ ду

(к)

«(kj+

д 2у

cd

V ^j )

\ (к).

И)2

(6)

f у \ (к) ду

где

j -я составляющая n-мерного вектора

|Т ^(к)т

д х.

Ц]( ^ = flx Г И У

(7)

д 2у

л?2 — значение j-го диагонального элемента диа-j

гональной матрицы n-го порядка B(k) , равное главному элементу на j-м шаге гауссового исключения:

д у _ (р_1 )(к)Тн(к) (Р_1 )(к) =в(к) д\2 '

Квадратичная функция (6) представлена в сепарабельном виде. Следовательно, выбор направления по каждой из новых координат независим и сводится к минимизации квадратичной функции одной переменной. При этом может быть три случая.

I.

Э2у

^ (к)

> 0.

^ (к)

ду

ф 0.

(9)

Функция Ду(к)= Ду (д|(к) j выпукла, расстояние до ее минимума

РИ, 2000, № 3

155

2.

^ (kj=-

( л ^ (Ю

ду

,%j

2 } (k)

О У

^ 2.

V J )

(10)

( о2 ) № ( ~\ (к)

52у

< 0.

і

V J

ф 0.

(її)

Функция Ду(к) вогнута, направление движения,

связанное с уменьшением ее значения, противоположно расположению точки максимума. Как долго двигаться в этом направлении? Выберем составляющую направляющего вектора д^(к) таким образом, чтобы ее абсолютное значение равнялось оценке расстояния до точки максимума:

3.

Л|(к) =

( л Л (к)

ду

д 2у

А (к)

52у

(к)

< 0,

V ^j У

(к)

= 0.

(12)

(13)

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

Условие (13) свидетельствует о том, что движение в любую сторону вдоль координаты приводит

к уменьшению Ay j. В частности, можно взять

Д|((к) = +

г 2 )(к)2

д 2у

JeT

г 2 'і(к)

д 2у

2

(14)

где Т — подмножество индексов координат |j, которым соответствует отрицательный диагональ-

д 2у

ный элемент ——. Знак определяется случайно или

^ 2 м

в результате анализа изменения функции уїх I при

противоположных изменениях координаты .

После выбора направления движения Д^^ в соответствии с формулами (10)-(14) можно перейти

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

выбору параметра , характерных для обобщен-

ного метода Ньютона. Если функция yfy строго

выпукла, то матрица Н^к) всегда положительно определена и, следовательно, все диагональные

элементы

д 2у

(к)

положительны. В этом случае

необходимость в условиях (10)-(14) отпадает и предлагаемый метод вырождается в чисто ньюто-

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

Кроме того, модифицированный метод Ньютона с приведением матрицы вторых производных к диагональному виду дает возможность найти направление убывания функции у(х) даже в том случае, если матрица H(k) особенная, эту матрицу можно преобразовать в диагональную, (n-r) элементов которой равны нулю. Здесь r — ранг матрицы H(k). При этом выбор направления убывания функции

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

Таким образом, предлагаемый алгоритм не позволяет определить ненулевой направляющий вектор

Дх^к) только в единственном случае, если мы

находимся в точке х(к), в которой матрица H(k)

5у )(к)

положительно полуопределена, а | —ь I = 0 . Для

Д х )

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

Для сравнительного анализа работы описанного метода и метода Ньютона оба они были реализованы в пакете компьютерной алгебры Maple Power V Edition.

Необходимо отметить, что в настоящее время существует довольно много пакетов компьютерной алгебры, в том числе наиболее распространенные из них Derive, Reduce, MatLab, Maple, MathCad и другие. Каждый из перечисленных пакетов имеет

156

РИ, 2000, № 3

свои достоинства, но многие авторы признают Maple лидером среди них благодаря символьному процессору [5]. Его превосходство признано даже фирмами-конкурентами — разработчиками пакетов MathCad, MatLab, которые приобрели ядро символьного процессора для внедрения в свои пакеты, но безусловно количество доступных в них функций значительно сокращено. Исходя из всего сказанного выше для реализации сравнительного анализа был выбран пакет Maple.

Реализация метода Ньютона и метода модификации средствами пакета Maple приведена на листинге 1 и 2 соответственно.

Листинг 1 Метод Ньютона

> newton: =proc(f:algebraic,x0 :vector, dx:float) local i,n,f0,x0k,c;

> n:=vectdim(x0);i:= 1;x0k:=x0;c:=false;k:=0;

> while c=false do

f0 :=subs(seq(x[i] =x0k[i] ,i= 1..n),f);unassign(‘i’);

> A:=subs(seq(x[i] =x0k[i] ,i= 1. .n),grad(f, [x[j] $j= 1. .n]));

> for l from 1 to n do c:=evalb(abs(A[l])<dx); if c=false then l:=n fi; od; k:=k+1;

> x0k:=evalm(x0k-inverse(subs(seq(x[i] =x0k[i] ,i= 1..n), hessian(f,[x[j]$j=1 ..n])))&*A) :od;print(k,x0k,f0);

> end;

Листинг 2

Модифицированный метод Ньютона

> ModNewton:=proc(f:algebraic,x0:vector,e:float) locali,j,n,c,l,A,B,Pt,x0k,fD,giz,z,k,summa,dx,min,step;

> n:=vectdim(x0); x0k:=x0; c:=false; step:=0;

> while c=false do A:=subs(seq(x[i]=x0k[i],i=1..n), hessian(f,[x[j]$j=1..n]));

>Pt:=inver3e(mygausselim(A,0));B:=mygausselim(A,1);

gr:=subs(seq(x[i]=x0k[i],i=1..n),grad(f,[x[j]$j=1..n]));

> grz:=evalm(subs(seq(x[i] =x0k[i],i= 1..n), grad(f, [x[j]$j=1. .n]))&*Pt);

> for l from 1 to n do

> if grz[l]<>0 then

> if B[l,l]>0 then z[l]:=-grz[l]/B[l,l] else z[l]:=grz[l]/B[l,l] fi;

> else for k from 1 to n do if B[k,k]<0 then summa:=summa+B[k,k] fi; od;

> z[l]:=(B[l,l]/summa)A(1/2); fi; od;

> z:=vector([seq (z[i],i=1..n)]);

> dx:=evalm(Pt&*z);

> fh:=subs(seq(x[i]=x0k[i] +x*dx[i] ,i=1..n),f);

> min:=golden_section(fh,0,1,0.001); f0:=subs(seq(x[i]=x0k[i],i=1..n),f);

> x0k:=evalm(x0k+min*dx);step:=step+1;

> for m from 1 to n do c:=evalb(abs(gr[m])<e); if c=false then m:=n fi;od;od;

> print(step,x0krro);end;

При реализации рассматриваемой модификации в качестве процедуры одномерной минимизации функции для расчета оптимального шага l(k) был использован метод золотого сечения, его реализация в пакете Maple представлена на листинге 3. В качестве вспомогательной процедуры была использована процедура гауссовых исключений, приведенная на листинге 4.

Листинг 3

Вспомогательная процедура одномерной минимизации функции методом золотого сечения

> golden_section:=proc(f:algebraic,x1 ,x2 ,e:float) local x01,x02,c,d,k,f1,f2,l;

> c:=x1;d:=x2;k:=0;l:=2;x02:=x1+0.618*(x2-x1); x01:=x1+0.382*(x2-x1); f1:=subs(x=x01,f);f2:=subs(x=x02,f);

> while d-c>e do iff1<f2 then d:=x02; x02:=x01;f2:=f1; x01:=c+0.382*(d-c);f1:=subs(x=x01,f); l:=l+1;

> else c:=x01; x01:=x02;f1:=f2; x02:=c+0.618*(d-c);f2:=subs(x=x02,f);l:=l+1;fi;

> k:=k+1; od;

> if subs(x=x01,f)<subs(x=x02,f) then x01; else x02; fi

> end;

Листинг 4

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

> mygausselim:=proc(A:matrix,t:integer) local P,B,i,j,k,l,m,n,h,s;

>s:=coldim(A);P:=matrix(s,s, [seq(seq(A[ij] j=1..s),i= 1. .s)]); B:=matrix(s,s,0);

> for k from 1 to s do B[k,k]:=P[k,k]; if k<s then for l from k+1 to s do

> for n from k+1 to s do

> P[l,n]:=P[l,n]-P[k,n]*P[l,k]/P[k,k];od;od;fi;

> if k>1 then for i from 1 to k-1 do

> for m from 1 to s do

> P[i,m]:=P[i,m]; od;od;fi;

> for j from 1 to s do P[k,j]:=P[k,j]/B[k,k];od;

if k<s then for h from k+1 to s do P[h,k]:=0;od; fi;

> od;if t=0 then evalm(P); else evalm(B); fi;

> end;

РИ, 2000, № 3

157

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

Таблица 1

Функции, используемые для тестирования предложенной модификации

Тестируя функцию №2, также имеем случай особенности матрицы Гесса в эксперименте №5, и как следствие — отсутствие решения методом Ньютона.

Таблица 4

Тестирование функции Химмельблау №3

№ Функция Минимум

1 100(Х2 - xi2)2 + (1 - xi)2 0 в (1;1)

2 (х2 - Xi2)2 + (1 - xi)2 0 в (1;1)

3 (X2 - X12)2 +100 * (1 - X1)2 0 в (1;1)

4 100(X2 - X12)3 + (1 - X1)2 0 в (1;1)

5 [1.5 - X1*(1 - X2)]2 + [2.25 - X1* *(1 -X2)]2 + [2.625- X1(1 - X32)]2 0 в (1;1)

6 2 3 1 2 2 2 ~X1 + 3X3 _ X1X2 + X1X2 _ 5x1 6,43 в (1,92; 0,792)

№ эксперимента 8 9 10

Начальная точка 2; -2 1.922; -3.222 1.986; 5.227

Ньютон К-во итераций 4 4 4

Время расчета 0.296 0.171 0.257

Решение 1.0; 0.(9) 1.0; 0.(9) 1.0; 1.0

Знач.ф-ции 0.17210-10 0.427-10-10 0.778-10-13

Мод. Н-на К-во итераций 4 4 4

Время расчета 0.608 0.776 0.807

Решение 1.0; 0.(9) 1.0; 0.(9) 1.0; 0.(9)

Знач. ф-ции 0.113710-9 0.185110-9 0.21310-11

Таблица 5

Тестирование функции Химмельблау №4

Результаты вычислений, а также количество итераций и машинного времени, необходимое для каждого расчета, сведены в табл. 2-7 для каждой из функций, перечисленных в табл.1. Эксперименты, в которых результат неверен или не получен, выделены в таблицах серым цветом.

Таблица2

Тестирование функции Химмельблау №1

№ эксперимента 1 2 3 4

Начальная точка -1.2;1.0 -2;-2 5.621; -3,635 -0.221; 0.639

X К-во итераций 6 — — —

О н Время расчета 0.339 — — —

2 uq Решение 1.0; 1.0 — — —

я Знач.ф-ции 0.1910-9 — — —

а К-во итераций 21 25 32 17

X Ё § Время расчета 2.279 2.715 5.302 2.417

Решение 1.0; 1.0 0.(9); 0.(9) 0.(9); 0.(9) 0.(9); 0.(9)

*3 Знач. ф-ции 0.1610-8 0.1410-8 0.110-8 0.1910-8

№ эксперимента 11 12 13

Начальная точка 1.2; -2 0.248; -3.082 -1.2; -1

X К-во итераций 4 4 6

о н Время расчета 0.170 0.254 0.260

2 uq Решение 1.0; 1.0 0.(9); 0.(9) 0.(9); 0.(9)

я Знач.ф-ции 0.5810-11 0.491210-7 0.957 10-14

й X К-во итераций 12 18 27

і я Время расчета 1.997 3.047 3.255

Решение 0.(9); 0.(9) 0.(9); 0.(9) 1.0; 1.0

1 Знач. ф-ции 0.913 10-" 0.25-10'11 0.1569 10-8

В экспериментах 8-13 с поставленной задачей с успехом справляются оба метода.

Таблица 6

Тестирование функции Химмельблау №5

При тестировании функции № 1 в 3 из 4 экспериментов метод Ньютона не дает решения из-за особенности матрицы Гесса, предложенная же модификация с успехом справляется с этой задачей.

Таблица 3

Тестирование функции Химмельблау №2

№ эксперимента 5 6 7

Начальная точка -2;-2 0.803; 0.211;

0.251 3.505

X К-во итераций — 4 5

о н Время расчета — 0.246 0.233

2 X Решение — 1.0; 0.(9) 0.(9); 0.(9)

я Знач.ф-ции — 0.1410-8 0.1510-6

Л X К-во итераций 8 4 4

я Время расчета 0.834 0.618 0.678

§ Решение 0.(9); 0.(9) 0.(9); 0.(9) 0.(9); 0.(9)

Знач. ф-ции 0.7910-5 0.910-8 0.2210-5

№ эксперимента 14 15 16 17

Начальная точка 0;0 8;0.2 5; 0.8 8; 0.8

Ньютон К-во итераций 8 6 15 8

Время расчета 0.480 0.393 0.960 0.5

Решение 0 ;1.0 2.(9); 0.4(9) 0.910-9; 1.0 0.1510-8; 1.0

Знач.ф-ции 14.203 0.2910-6 14.203 14.203

Мод. Н-на К-во итераций 6 6 9 0

Время расчета 0.768 0.847 1.206 1.204

Решение 2.(9); 0.4(9) 3.0;0.5 2.(9); 0.4(9) 2.(9);0.5

Знач. ф-ции 0.00001 0.000002 0.53 10-8 0.2(6)10-7

Как видно из табл.6, в экспериментах 14, 16 и 17 минимизация тестируемой функции №5 методом Ньютона заканчивается в точке максимума.

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

Таким образом, было проанализировано 20 тестовых примеров. Из сводной таблицы видно, что при

158

РИ, 2000, № 3

Таблица 7

Тестирование функции Евдокимова №6

№ эксперимента 18 19 20

Начальная точка -1; 0.8 0.6285; -1.518 0.561; 1

Ньютон К-во итераций 5 1 15

Время расчета 0.212 0.3 0.625

Решение -0.628; -1.517 0.628; -1.51 -0.628; -1.517

Знач.ф-ции 2.09 -2.09 2.09

Мод. Н-на К-во итераций 8 11 3

Время расчета 0.659 0.92 0.25

Решение 1.929; 0.799 1.929; 0.799 1.929; 0.799

Знач. ф-ции -6.431 -6.431 -6.431

достаточно “хорошем” начальном приближении метод Ньютона работает быстрее и находит минимум за меньшее количество итераций, чем рассматриваемая модификация метода Ньютона. Но даже в случаях “хорошего” приближения метод Ньютона не всегда работает из -за особенности матрицы Гесса (случаи 2-5 в табл.2 и 3). При выборе начальной точки вблизи другого экстремума метод Ньютона приходит в максимум или же в седловую точку.

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

но не превышает 5 миллисекунд (расчет производился на Pentium II), что с лихвой окупается его повышенной сходимостью, т.е. сходимостью в 100% экспериментов.

Тестовая функция Евдокимова (№6 в табл.7) иллюстрирует наиболее характерные расходимости метода Ньютона. Как видим, приведенная в статье его модификация позволяет легко выйти из указанных точек- ловушек, взятых в качестве начального приближения.

Литература: 1. Евдокимов А Г. Минимизация функций и ее приложение к задачам автоматизированного управления инженерными сетями. Харьков: Вища шк. Изд-во при Харьк. ун-те, 1985. 288 с. 2. Химмельблау Д. Прикладное нелинейное программирование. М. Мир, 1975. 534с. 3. Фиакко А., Мак-Кормик Г. Нелинейное программирование. Методы последовательной безусловной оптимизации. М.: Мир, 1972. 240с. 4. BeightlerC. S., Wilde D.S. Foundations of Optimization. Prentice-Hall, 1967. 469p. 5. Манзон М. Maple V Power Edition. М.: БИНОМ, 1998.

Поступила в редколлегию 12.06.2000

Рецензент: д-р техн. наук, проф. Левыкин В.М.

Евдокимов Анатолий Гаврилович, д-р техн. наук, профессор, зав. кафедрой ПМ и ВТ ХГАГХ. Научные интересы: моделирование и оптимизация потокораспределения в инженерных сетях, математическое программирование. Адрес: Украина, 61002, Харьков, ул. Революции, 12.

Манакова Наталья Олеговна, аспирант кафедры ПМ и ВТ ХГАГХ. Научные интересы: компьютерная алгебра, управлене сетями, оптимизация функций, программирование. Адрес: Украина, 61002, Харьков, ул. Революции, 12, тел. 45-50-86.

УДК 530.145

МАГНИТООПТИЧЕСКИЕ ЭФФЕКТЫ В БИОСРЕДАХ, ВЗАИМОДЕЙСТВУЮЩИХ С КВЧ-ЭЛЕКТРОМАГНИТНЫМ ПОЛЕМ

ЧОВНЮК Ю.В, ОВСЯННИКОВА Т.Н.,

РУДВКО Б.Ф._______________________

Описывается исследование магнитооптических эффектов, возникающих в биосредах и взаимодействующих с электромагнитными полями КВЧ-диапазона (несущая частота f«60 ГГц): 1) анализ намагничения немагнитной прозрачной биосреды переменным электрическим полем (КВЧ-диапазона); 2) определение поляризации отраженного света (КВЧ-сигнала) при нормальном падении линейно-поляризованной волны из пустоты на поверхность изотропного биотела, находящегося в магнитном поле (как индуцированном переменным электрическим полем КВЧ-диапазона, так и постоянным, обусловленным Землей); 3) установление предельного закона зависимости вектора гирации от частоты при больших значениях последней (биосреда предполагается таковой, имеющей естественную оптическую активность).

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

1. Известно [1], что компоненты вектора оптической активности биосреды G являются функциями напряженности магнитного поля H . Обычно магнитное поле сравнительно слабое, поэтому можно произвести разложение G по степеням H (здесь и далее применяется гаусова система единиц). Источником постоянной составляющей H служит магнитное поле планеты Земля. В отсутствие магнитного поля вектор G равен нулю; поэтому в слабом поле можно положить:

Gi = fikHk,(i,k) = (U). (1)

Здесь fik — тензор второго ранга, в общем случае несимметричный. Следует отметить, что такая форма зависимости находится в согласии с общим правилом, согласно которому в прозрачной биосреде компоненты антисимметричного тензора p”ik (как и тензора s"jk), описывающего поглощающие свойства биосреды

РИ, 2000, № 3

159