Определение максимального и минимального значения функции двух переменных методом градиентного спуска (подъема)

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

Для многомерного случая в качастве аналога производной используется понятие «градиента».

Градиентом функции называется вектор, компоненты которого являются частными производными функции по соответствующим координатам. Так, например, для функции двух переменных f(x,y) градиент определяется как вектор с компонентами: и . Обозначается градиент функции либо без значка вектора.

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

При численных расчетах градиент функции можно определить приближенно , где , – единичные вектора по направлениям осей X и Y, , , – приращение x,  – приращение y.

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

Метод определения максимума (минимума) функции, основанный на вычислении градиента функции называется методом наискорейшего подъема (спуска). В данном итерационном методе вектор шага определяется направлением вектора градиента:

. (4)

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

Рассмотрим многомерную функцию, имеющую вид колокола:

, (5)

где точка – точка максимума функции.

В области изменений функция имеет вид:

Аналогичные функции можно построить в многомерном пространстве. Например, если размерность пространства n, то функция «n-мерный колокол» будет иметь вид:

. (6)

Точка максимума этой функции - . Функция «n-мерного колокола» (6) при n=2 превращается в (5).

Описание алгоритма нахождения максимального (минимального) значения функции двух переменных методом градиентного спуска (подъема):

Задаем размерность n массива аргументов минимизируемой (максимизируемой) многомерной функции, погрешность вычислений ε_, начальные значения шага d по X и Y, координаты начальной точки поиска x₀, y₀, определяем вспомогательный массив для запоминания значений аргументов при дроблении шага х₁, y₁, массив для составляющих градиента функции df.

Вычисляем значение функции в точке с координатами (x₀, y₀).

Вычисляем х и y компоненты вектора градиента функции по формулам:

gx=(f(x₀+dx, y₀)-f(x₀,y₀))/dx;

gy=(f(x₀, y₀+dy)-f(x₀,y₀))/dy.

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

Вычисляем новое значение функции и сравниваем с предыдущим. Если максимум (минимум) не найден, продолжаем итерации.

Если второе условие не выполнилось, уменьшаем шаг.

Алгоритм на естественном языке приведем для функции «n-мерного колокола».

Строки «комментарии» будем начинать с символа «*».

Для функции с координатами точки максимума: (2, 3, 4, …, n+1) для пространства размерности n.

Определить n как целое, где n – количество аргументов функции.

n=2

Определить массивы x0(n), x1(n) , x2(n), df(n), h(n)

* массивы x0, x1, х2 – рабочие массивы координат точек (0), (1),

* (2). Массивы df – компоненты вектора градиента функции,

* h – компоненты векторов шага

Определить глобальный массив xс(n)

* массив xс-массив координат точки максимума функции

eps=0.01

d=0.01

* в данном цикле задаются координаты точки максимума

Цикл по i от 1 до n

xc(i)=1+i

Конец Цикла

* здесь задаются координаты начальной точки для определения

* максимума

Цикл по i от 1 до n

x0(i)=0.1+i

Конец Цикла

h0=0.5

f0=колокол(@x0,n)

Цикл по k от 1 до 10000

Печать "k=",x0(1),x0(2),f0,h0

Цикл по j от 1 до n

Цикл по i от 1 до n

x1(i)=x0(i)

Конец Цикла

x1(j)=x0(j)+d

f1=колокол(@x1,n)

df(j)=(f1-f0)/d

Конец Цикла

* нормирование градиента

lendf=0

Цикл по i от 1 до n

lendf = lendf + df(i)^2

Конец Цикла

lendf = Корень(lendf)

Цикл по i от 1 до n

df(i)=df(i)/lendf

Конец Цикла

Цикл по i от 1 до n

x0(i)=x0(i)+df(i)*h0

Конец Цикла

f0_1=колокол(@x0,n)

Если f0_1=<f0

Печать " нашли max"

Если h0<eps

Печать "найден max x0=(",

Цикл по i от 1 до n

Печать x0(i)

Конец Цикла

Печать ")"

Печать " f0= ", f0

Печать "за ",k," итераций"

Выход Из Программы

Конец Если

h0=h0/2

Конец Если

f0=f0_1

Конец Цикла

Печать "максимум не найден x="

Цикл по i от 1 до n

Печать x0(i)

Конец Цикла

ОПРЕДЕЛИТЬ ФУНКЦИЮ колокол

Параметры x, n

определить массив x(n)

определить локальную переменную r2

r2=0

Цикл по i от 1 до n

r2=r2+(x(i)-xc(i))^2

Конец Цикла

Если r2>10

r2=10

Конец Если

Вернуть 5*(EXP(-0.1*r2))

Пример решения на языке VFP:

CLEAR

SET DECIMALS TO 15

n=10

DIMENSION x0(n), x1(n), df(n), h(n), x2(n)

PUBLIC ARRAY xc(n)

FOR i=1 TO n

xc(i)=1+i

ENDFOR

eps=0.01

d=0.01

FOR i=1 TO n

x0(i)=1.1+i

ENDFOR

h0=0.5

f0=функция(@x0,n)

FOR k=1 TO 10000

? "k=",STR(x0(1),5,2),STR(x0(2),5,2),f0,h0

FOR j=1 TO n

FOR i=1 TO n

x1(i)=x0(i)

ENDFOR

x1(j)=x0(j)+d

f1=функция(@x1,n)

df(j)=(f1-f0)/d

ENDFOR

* нормирование градиента

lendf=0

FOR i=1 TO n

lendf = lendf + df(i)^2

ENDFOR

lendf = SQRT(lendf)

FOR i=1 TO n

df(i)=df(i)/lendf

ENDFOR

FOR i=1 TO n

x0(i)=x0(i)+df(i)*h0

ENDFOR

f0_1=функция(@x0,n)

IF f0_1=<f0

? " нашли max"

IF h0<eps

? "найден max x0=(",

FOR i=1 TO n

??x0(i)

ENDFOR

??")"

?" f0= ", f0

? "за ",k," итераций"

RETURN

endif

h0=h0/2

ENDIF

f0=f0_1

endfor

? "максимум не найден x="

FOR i=1 TO n

??x0(i)

ENDFOR

FUNCTION функция

LPARAMETERS x, n

DIMENSION x(n)

LOCAL r2

r2=0

FOR i=1 TO n

r2=r2+(x(i)-xc(i))^2

ENDFOR

IF r2>10

r2=10

ENDIF

RETURN 5*(EXP(-0.1*r2))

Пример решения на языке VBA:

'объявление глобального массива xc() неопределенной длины

Public xc()

Sub grd()

Dim x0(), x1(), df(), h(), x2()

n = 10

ReDim x0(n), x1(n), df(n), h(n), x2(n)

ReDim xc(n)

For i = 1 To n

xc(i) = 1 + i

Next

eps = 0.01

d = 0.01

For i = 1 To n

x0(i) = 1.1 + i

Next

h0 = 0.5

f0 = f2(x0, n)

For k = 1 To 10000

Debug.Print "k= " & k & " : " & x0(1) & " : " & x0(2) & " : " & f0 & " : " & h0

For j = 1 To n

For i = 1 To n

x1(i) = x0(i)

Next

x1(j) = x0(j) + d

f1 = f2(x1, n)

df(j) = (f1 - f0) / d

Next

' нормирование градиента

lendf = 0

For i = 1 To n

lendf = lendf + df(i) ^ 2

Next

lendf = Sqr(lendf)

For i = 1 To n

df(i) = df(i) / lendf

Next

For i = 1 To n

x0(i) = x0(i) + df(i) * h0

Next i

f0_1 = f2(x0, n)

If f0_1 <= f0 Then

If h0 < eps Then

s = ""

For i = 1 To n - 1

s = s & x0(i) & ","

Next i

s = s & x0(n)

Debug.Print "Максимум найден за " & k & " итераций!"

Debug.Print "x0=(" & s & ")"

Debug.Print "f0= " & f0

Exit Sub

End If

h0 = h0 / 2

End If

f0 = f0_1

Next k

Debug.Print "максимум не найден!"

End Sub

Function f2(ByRef x, n)

Dim r2

r2 = 0

For i = 1 To n

r2 = r2 + (x(i) - xc(i)) ^ 2

Next

If r2 > 10 Then

r2 = 10

End If

f2 = 5 * (Exp(-0.1 * r2))

End Function

Результат работы программы на VBA:

k= 1 : 2,1 : 3,1 : 4,95024916874584 : 0,5

k= 2 : 1,94188611699158 : 2,94188611699158 : 4,98314236503072 : 0,5

k= 3 : 2,1 : 3,1 : 4,95024916874584 : 0,25

k= 4 : 2,02094305849579 : 3,02094305849579 : 4,99780742238446 : 0,25

k= 5 : 1,94188611699158 : 2,94188611699158 : 4,98314236503072 : 0,125

k= 6 : 1,98141458774369 : 2,98141458774369 : 4,99827321050518 : 0,125

k= 7 : 2,02094305849579 : 3,02094305849579 : 4,99780742238446 : 0,0625

k= 8 : 2,00117882311974 : 3,00117882311974 : 4,99999305188509 : 0,0625

k= 9 : 1,98141458774369 : 2,98141458774369 : 4,99827321050518 : 0,03125

k= 10 : 1,99129670543171 : 2,99129670543171 : 4,99962127766207 : 0,03125

k= 11 : 2,00117882311974 : 3,00117882311974 : 4,99999305188509 : 0,03125

k= 12 : 1,99129670543171 : 2,99129670543171 : 4,99962127766207 : 0,015625

k= 13 : 1,99623776427573 : 2,99623776427573 : 4,99992922841264 : 0,015625

k= 14 : 1,99129670543171 : 2,99129670543171 : 4,99962127766207 : 0,0078125

k= 15 : 1,99376723485372 : 2,99376723485372 : 4,9998057669659 : 0,0078125

k= 16 : 1,99623776427565 : 2,99623776427565 : 4,99992922841264 : 0,0078125

Максимум найден за 16 итераций!

x0=(1,99376723485379 ; 2,99376723485379 ; 3,99376723485379 ; 4,99376723485379 ; 5,99376723485379 ; 6,99376723485379 ; 7,99376723485361 ; 8,99376723485361 ; 9,99376723485361 ; 10,9937672348536)

f0= 4,99992922841264