- •Нахождение экстремумов
- •Построение графика
- •Численный метод нахождения экстремума
- •Теоретические вопросы
- •Метод деформируемого многогранника (Нелдера Мида)
- •Метод наискорейшего спуска
- •Метод сопряженных направлений и его модификации
- •Метод Ньютона и его модификации
- •Метод дробления шага
- •Скрипт MATLAB для построения графика
- •Описание функции
- •Текст скрипта
- •Реализация алгоритма Ньютона
- •Тест реализации алгоритма Ньютона
BРеализация алгоритма Ньютона
1using System;
2
3namespace MathUtilities
4{
5
6
7
8
9
///<summary>
///Класс для поиска оптимума методом Ньютона.
///</summary>
public class NewtonOptimizer
{
10/// <summary>
11/// Функция для исследования.
12/// </summary>
13private readonly Func<Vector, double> _function;
14
15/// <summary>
16/// Класс для нахождения производной, градиента и Гессиана.
17/// </summary>
18private readonly Derivative _derivative;
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
///<summary>
///Начальный вектор поиска.
///</summary>
private readonly Vector _start;
///<summary>
///Точность поиска.
///</summary>
private readonly double _precision;
///<summary>
///Максимальное количество итераций.
///</summary>
private readonly int _operationsAmount;
35
36
37
38
39
///<summary>
///Счётчик операций.
///</summary> private int counter;
40/// <summary>
41/// Стандартный конструктор
42/// </summary>
43/// <param name="function">функция для исследования</param>
44/// <param name="start">начальный вектор</param>
45/// <param name="precision">точность поиска</param>
11
46/// <param name="operationsAmount">максимальное количество ,! итераций</param>
47public NewtonOptimizer(
48
49
50
51
52
53
54
55
56
57
58
59
60
Func<Vector, double> function, double[] start,
double precision = 1E-08, int operationsAmount = 8
) {
_function = function;
_derivative = new Derivative(function); _start = new Vector(start);
_precision = precision; _operationsAmount = operationsAmount; ResetCounter();
}
61/// <summary>
62/// Рассчитывает точке, в которой функция имеет минимум.
63/// </summary>
64public Vector GetOptimal()
65{
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
// градиент в начальной точке
var gradient = _derivative.GetGradient(_start);
// гессиан в начальной точке
var hessian = _derivative.GetHessian(_start);
//запоминаем предыдующий вектор var previous = _start;
//считаем текущий вектор по формуле
var current = _start - hessian.GetReversed() * gradient;
// пока поле поиска достаточно большое while (!IsTooSmall(previous, current))
{
// если операций не больше максимального if (IsCountDown())
throw new ApplicationException("не нашлось решения");
82
83
84
85
86
87
88
89
90
91
//запомнить предыдущий вектор previous = current;
//рассчитать градиент в текущей точке gradient = _derivative.GetGradient(previous);
//рассчитать гессиан в текущей точке
hessian = _derivative.GetHessian(previous);
// пересчитать текущий вектор
current = previous - hessian.GetReversed() * gradient;
}
12
92 |
// вернуть найденный вектор |
93 |
return current; |
94 |
} |
95
96/// <summary>
97/// Не превысило ли количество операций максимальное.
98/// </summary>
99private bool IsCountDown() =>
100 |
counter++ > _operationsAmount; |
101 |
|
102
103
104
105
106
107
///<summary>
///Обнуляет количество проделанных операций.
///</summary>
private void ResetCounter() => counter = 0;
108/// <summary>
109/// Проверяет, является ли новое решение слишком близким.
110/// </summary>
111/// <param name="previous">старое решение</param>
112/// <param name="current">текущее решение</param>
113/// <returns>истина, если решения находятся близко друг к ,! другу</returns>
114private bool IsTooSmall(Vector previous, Vector current) =>
115 |
Math.Abs(_function(current) - _function(previous)) < _precision; |
116 |
} |
117 |
|
118
119
120
121
122
123
124
125
126
127
///<summary>
///Класс, представляющий вектор.
///</summary>
public class Vector
{
///<summary>
///Элементы вектора.
///</summary>
public double[] items { get; private set; }
128/// <summary>
129/// Стандартный конструктор, задаёт элементы вектора.
130/// </summary>
131/// <param name="items">элементы вектора</param>
132public Vector(params double[] items) =>
133
134
this.items = items;
135
136
137
///<summary>
///Передаёт элемент вектора в индексе.
///</summary>
13
138/// <param name="index">индекс элемента</param>
139/// <returns>элемент в индексе</returns>
140public double Get(int index) =>
141
142
143
144
145
146
147
items[index];
///<summary>
///Количество элементов вектора.
///</summary>
public int Length => items.Length;
148/// <summary>
149/// Переопределение преобразования в строку:
150/// все элементы через запятую.
151/// </summary>
152public override string ToString() =>
153 |
String.Join(", ", items); |
154
155/// <summary>
156/// Определение оператора сложения векторов.
157/// </summary>
158/// <param name="left">левый вектор</param>
159/// <param name="right">правый вектор</param>
160/// <returns>сумма векторов - новый вектор</returns>
161public static Vector operator +(Vector left, Vector right)
162{
163 |
int length = left.Length; |
164 |
var answerItems = new double[length]; |
165 |
for (int i = 0; i < length; i++) |
166 |
answerItems[i] = left.Get(i) + right.Get(i); |
167 |
return new Vector(answerItems); |
168 |
} |
169
170/// <summary>
171/// Определение оператора разницы векторов.
172/// </summary>
173/// <param name="left">уменьшаемый вектор</param>
174/// <param name="right">вычитаемый вектор</param>
175/// <returns>разница векторов - новый вектор</returns>
176public static Vector operator -(Vector left, Vector right)
177{
178 |
int length = left.Length; |
179 |
var answerItems = new double[length]; |
180 |
for (int i = 0; i < length; i++) |
181 |
answerItems[i] = left.Get(i) - right.Get(i); |
182 |
return new Vector(answerItems); |
183 |
} |
184 |
|
14
185/// <summary>
186/// Определение операции скалярного произведения векторов.
187/// </summary>
188/// <param name="left">левый вектор</param>
189/// <param name="right">правый вектор</param>
190/// <returns>скалярное произведение векторов - число</returns>
191public static double operator *(Vector left, Vector right)
192{
193 |
int length = |
left.Length; |
|
194 |
double answer = |
0; |
|
195 |
for (int i = |
0; |
i < length; i++) |
196 |
answer += left.Get(i) * right.Get(i); |
||
197 |
return answer; |
|
198}
199}
200
201/// <summary>
202/// Класс, представляющий матрицу.
203///
204/// Прим.: все функции работают только с квадратными матрицами.
205/// Прим.: некоторые функции работают только с матрицей 2 на 2.
206/// </summary>
207public class Matrix
208{
209/// <summary>
210/// Разделитель элементов матрицы в строковом представлении.
211/// </summary>
212private const string ITEM_DELIMITER = ", ";
213
214/// <summary>
215/// Разделитель строк матрицы в строковом представлении.
216/// </summary>
217private const string ROW_DELIMITER = "\n";
218
219
220
221
222
223
///<summary>
///Элементы матрицы.
///</summary>
public double[][] items { get; private set; }
224/// <summary>
225/// Стандартный конструктор.
226/// </summary>
227/// <param name="items">двумерный массив элементов матрицы</param>
228public Matrix(double[][] items) =>
229 this.items = items;
230
231
15
232/// <summary>
233/// Элемент матрицы в индексах.
234/// </summary>
235/// <param name="index1">столбец</param>
236/// <param name="index2">строка</param>
237public double Get(int index1, int index2) =>
238 |
items[index1][index2]; |
239 |
|
240
241
242
243
244
///<summary>
///Длина стороны матрицы.
///</summary>
public int Length => items.Length;
245/// <summary>
246/// Переопределение преобразования в строку:
247/// все элементы через запятую, строки через перевод строки.
248/// </summary>
249/// <returns></returns>
250public override string ToString()
251{
252 |
int length = |
Length; |
253 |
var arr = new string[length]; |
|
254 |
for (var i = |
0; i < length; i++) |
255 |
arr[i] = |
String.Join(ITEM_DELIMITER, items[i]); |
256 |
return String.Join(ROW_DELIMITER, arr); |
|
257 |
} |
|
258
259/// <summary>
260/// Умножение матрицы на вектор.
261/// </summary>
262/// <param name="matrix">матрица</param>
263/// <param name="vector">вектор</param>
264/// <returns>произведение - новый вектор</returns>
265public static Vector operator *(Matrix matrix, Vector vector)
266{
267
268
269
int length = vector.Length; var items = new double[length]; Vector temp;
270
271
272
273
274
275
276
277
//просто перемножаем строки на вектор,
//так как строки матрицы - такие же векторы for (int i = 0; i < length; i++)
{
temp = new Vector(matrix.items[i]); items[i] = temp * vector;
}
278
16
279 |
return new Vector(items); |
280 |
} |
281
282/// <summary>
283/// Передаёт обратную матрицу.
284///
285/// Здесь представлена только оптимизированная версия для матрицы 2 на
,! 2.
286/// </summary>
287/// <returns>обратная матрице матрица</returns>
288public Matrix GetReversed()
289{
290 |
// получаем определитель матрицы |
|
291 |
var det = GetSmallDeterminant(); |
|
292 |
// создаём пустой двумерный массив |
|
293 |
var matrix = new double[2][]; |
|
294 |
// заполняем первую строку |
по формуле |
295 |
matrix[0] = new double[2]; |
|
296 |
matrix[0][0] = items[1][1] / det; |
|
297 |
matrix[0][1] = -1 * items[0][1] / det; |
|
298 |
// заполняем вторую строку |
по формуле |
299 |
matrix[1] = new double[2]; |
|
300 |
matrix[1][0] = -1 * items[1][0] / det; |
|
301 |
matrix[1][1] = items[0][0] / det; |
|
302 |
return new Matrix(matrix); |
|
303 |
} |
|
304
305/// <summary>
306/// Считает определитель этой матрицы по оптимизированной формуле для ,! матрицы 2 на 2.
307/// </summary>
308/// <returns>определитель матрицы</returns>
309public double GetSmallDeterminant() =>
310
311
312
313
314
315
316
317
318
319
320
321
322
323
items[0][0] * items[1][1] - items[0][1] * items[1][0];
}
///<summary>
///Класс для численного дифференцирования.
///</summary>
public class Derivative
{
///<summary>
///Функция для дифференцирования.
///</summary>
private readonly Func<Vector, double> _function;
17
324/// <summary>
325/// Стандартный конструктор для дифференцирования.
326/// </summary>
327/// <param name="function">функция для дифференцирования</param>
328public Derivative(Func<Vector, double> function) =>
329 |
_function = function; |
330
331/// <summary>
332/// Градиент функции.
333/// </summary>
334/// <param name="vector">точка, в которой ищем градиент</param>
335/// <returns>вектор градиента</returns>
336public Vector GetGradient(Vector vector)
337{
338 |
int length = vector.Length; |
339 |
var items = new double[length]; |
340 |
for (int i = 0; i < length; i++) |
341 |
items[i] = GetPartialDerivative(vector, i); |
342 |
return new Vector(items); |
343 |
} |
344
345/// <summary>
346/// Рассчитывает гессиан - матрицу вторых частных производных.
347/// </summary>
348/// <param name="vector">точка, в которой ищем градиент</param>
349/// <returns>матрица Гессе</returns>
350public Matrix GetHessian(Vector vector)
351{
352
353
354
355
356
357
358
359
360
int length = |
vector.Length; |
var answer = |
new double[length][]; |
for (int i = |
0; i < vector.Length; i++) |
{ |
|
answer[i] = new double[length]; for (int j = 0; j < length; j++)
answer[i][j] = GetSecondPartialDerivative(vector, i, j);
}
361 |
return new Matrix(answer); |
362 |
} |
363 |
|
364/// <summary>
365/// Рассчитывает первую частную производную.
366/// </summary>
367/// <param name="vector">точка, в которой ищем производную</param>
368/// <param name="index">индекс переменной ФНП</param>
369/// <param name="usePrecision">использовать ли последующую точность ,! расчётов</param>
18
370/// <param name="savedPrecision">точность расчётов</param>
371/// <returns>значение первой частной производной в точке</returns>
372public double GetPartialDerivative(
373 |
Vector vector, |
374 |
int index, |
375 |
bool usePrecision = false, |
376 |
double savedPrecision = double.NaN |
377 |
) { |
378 |
// клонируем массивы, так как иначе переменная |
379 |
// будет являться лишь ссылкой на старый массив |
380 |
var leftvec = (double[]) vector.items.Clone(); |
381 |
var rightvec = (double[]) vector.items.Clone(); |
382 |
// точность расчётов - см. исследование БМВ |
383 |
var precision = usePrecision ? |
384 |
Math.Max(savedPrecision, |
|
,! Infinitesimal.GetPrecision(vector.Get(index))) : |
385 |
Infinitesimal.GetPrecision(vector.Get(index)); |
386 |
|
387
388
389
leftvec[index] += precision; rightvec[index] -= precision;
390 |
return ( |
|
391 |
_function( |
|
392 |
|
new Vector(leftvec)) - |
393 |
_function( |
|
394 |
|
new Vector(rightvec)) |
395 |
)/(2 |
* precision); |
396 |
} |
|
397 |
|
|
398/// <summary>
399/// Вторая частная производная в точке. Например, F''_x1_x2
400/// </summary>
401/// <param name="vector">точка</param>
402/// <param name="index1">индекс переменной 1</param>
403/// <param name="index2">индекс переменной 2</param>
404/// <returns>значение второй производной</returns>
405public double GetSecondPartialDerivative(Vector vector, int index1, int
|
,! index2) |
406 |
{ |
407
408
409
410
411
412
//клонируем массивы, так как иначе переменная
//будет являться лишь ссылкой на старый массив var leftvec = (double[]) vector.items.Clone(); var rightvec = (double[]) vector.items.Clone();
//точность расчётов - см. исследование БМВ
var precision = Infinitesimal.GetPrecision(vector.Get(index1));
413
414
19
415
416
417
leftvec[index1] += precision; rightvec[index1] -= precision;
418 |
return ( |
|
419 |
GetPartialDerivative( |
|
420 |
|
new Vector(leftvec), index2, true, precision) - |
421 |
GetPartialDerivative( |
|
422 |
|
new Vector(rightvec), index2, true, precision) |
423 |
)/(2 |
* precision); |
424}
425}
426
427/// <summary>
428/// Небольшой дополнительный класс для работы с малыми числами.
429/// Содержит лишь функцию избавления от неточных вычислений.
430/// </summary>
431public class Infinitesimal
432{
433/// <summary>
434/// Типичный машинный эпсилон для вычислений.
435/// </summary>
436public const int TYPICAL_EPSILON = 5;
437
438/// <summary>
439/// Передаёт минимальное число - степень десятки, с которым может ,! работать это число.
440/// </summary>
441/// <param name="number">число для расчёта; при 0 используется типичный ,! машинный эпсилон</param>
442/// <returns>1E с минимальным числом</returns>
443public static double GetPrecision(double number)
444{
445 |
if (number == 0) |
446 |
return Math.Pow(10, -1 * TYPICAL_EPSILON); |
447 |
// отрицательные числа тоже имеют порядок |
448 |
var abs = Math.Abs(number); |
449 |
// находим порядок числа |
450 |
var log = Math.Log10(abs); |
451 |
// удаляем дробную часть |
452 |
var rounded = Math.Round(log); |
453 |
// получаем новую степень десятки - входящее число может работать |
|
,! только с числами степени не ниже эпсилона |
454 |
var pow = rounded - TYPICAL_EPSILON; |
455 |
return Math.Pow(10, pow); |
456}
457}
458}
20