Реализация метода покоординатного спуска

Приветствую!
В данной статье попробую доходчивым для всех языком объяснить реализацию метода поокрдинатного спуска (другое название метода Гаусса-Зейделя) на языке C++. Для поиска минимума одномерной функции используется метод золотого сечения, поэтому затрону и его.
Вообще я долго думал, выкладывать статью или нет, так как в интернете на различных форумах можно найти разные реализации этих методов, и в случае необходимости скомпоновать их вместе. Все же решил выложить по трем причинам :
1) в большинстве случаев на форумах была просто реализация, без каких-либо комментариев и пояснений, а я постараюсь максимально подробно объяснить что, как и почему я делал;
2) все реализации, которые я видел, подстроены под одну конкретную функцию, а я попробовал сделать код максимально универсальным;
3) мне не удалось найти реализацию метода покоординатного спуска с оптимизацией одномерной функции именно методом золотого сечения...сомнительная причина, но все же причина.
В общем, надеюсь, что кому-нибудь эта статья будет полезна.

Приступим.

Метод покоординатного спуска (метод Гаусса-Зейделя)

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

Метод золотого сечения

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

Реализация

Начну с реализации метода золотого сечения. В рассматриваемом случае на вход метода подаются следующие параметры:
1) целевая функция, минимум которой требуется найти;
2) массив значений всех аргументов целевой функции. Т.к. целевая функция многомерная, а метод золотого сечения предназначен для поиска экстремума одномерной функции, то надо задать фиксированные значения всем переменным, кроме той, по которой ищем экстремум;
3) номер переменной (аргумента целевой функции) из массива из п.2, по которой будем искать экстремум;
4) значение погрешности;
5) концы отрезка, на котором осуществляется поиск экстремума;
6) максимально допустимое количество итераций цикла, нужно для того, чтоб можно было завершить метод, если заданная погрешность не может быть достигнута.
Результатом выполнения метода золотого сечения будет значение аргумента, в котором функция, в нашем случае, достигает минимум.
Таким образом, прототип функции, реализующей метод золотого сечения, будет иметь вид:

double golden_section(func_ptr f, double* vars, int var_index, double eps, double a, double b, int max_steps_count);

func_ptr - указатель на функцию, возвращающую значение вещественного типа а принимающей массив значений вещественного типа
Алгоритм метода золотого сечения для минимизации функции складывается из следующих этапов:
1) Вычисляется значение функции f(x1), где x1 = a + phi * (b - a) (phi = (1 + sqrt(5.) / 2.0);
2) Вычисляется значение функции f(x2), где x2 = b - phi * (b - a) (phi = (1 + sqrt(5.) / 2..
3) Определяется новый интервал (a, x2) или (x1, b), в котором локализован минимум.
4) Внутри полученного интервала находится новая точка (x1 в случае интервала (x1, b)) или (x2 в случае интервала (a, x2)), отстоящая от его конца на расстоянии, составляющем phi (1 + sqrt(5.) / 2. от его длины. В этой точке рассчитывается значение f(x). Затем вычисления повторяются, начиная с пункта 3, до тех пор, пока величина интервала неопределенности станет меньше заданной погрешности.
Реализация вышеописанного алгоритма на C++ выглядит следующим образом:

double golden_section(func_ptr f, double* vars, int var_index, double eps, double a, double b, int max_steps_count) { double res = 0.0; // результат, возвращаемый функцией double phi = (1 + sqrt(5.) / 2.0; // магический коэффициент золотого сечения double A = 0.0f, B = 0.0f; // тут будет храниться значение функции в предполагаемых точках экстремума double x1 = a + phi * (b - a), x2 = b - phi * (b - a); // предполагаемая точка экстремума и точка симметричная ей на отрезке [a, b]

int step = 0; // текущий шаг

while((b - a > eps)) // повторяем цикл до тех пор, пока разность концов отрезков не будет меньше заданной погрешности
{
x1 = b - ((b - a) / phi); // пересчитываем предполагаемую точку экстремума
vars[var_index] = x1; // помещаем ее в массив аргументов целевой функции
A = f(vars); // пересчитываем значение функции в этой точке
// далее аналогичные действия для точки симметричной x1
x2 = a + ((b - a) / phi);
vars[var_index] = x2;
B = f(vars);
// сужаем инетрвал неопределенности
if(A > B)
a = x1;
else
b = x2;

step++;
if(step > max_steps_count)
break;
}

res = (a + b) / 2;
return res;
}

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

void descent_method(func_ptr f, double* vars, double eps, int max_steps_count) { double B = f(vars), A = 0; // значение целевой функции на текущей и пердыдущей итерациях bool was_counted = false; // был ли определен минимум за максимально допустимое количество итераций int stpes_ellapsed = 0; // количество итераций, затраченных на поиск минимума (статистическая информация) double delta = 0.0; // полученная погрешность (статистическая информация) for(int i = 0; i < max_steps_count; i++){ A = B; // предыдущее значение целевой функции

for(int var_index = 0; var_index < var_count; var_index++) // проходим по каждому аргументу целевой функции
vars[var_index] = golden_section(f, vars, var_index, eps, -5000, 5000, max_steps_count); // и находим минимум одномерной функции по этой пе

B = f(vars); // текущее значение целевой функции

delta = fabs(A - B);

if(delta <= eps)
{
stpes_ellapsed = i + 1;
was_counted = true;
break;
}
}

std::cout << "Результат поиска минимума функции " << "exp(x1 + x2 + x3) / (x1 * x2^2 * x3^3)" << std::endl;

if(!was_counted)
std::cout << "За максимально указанное количество шагов ( " << max_steps_count << " ) минимум не был посчитан." << std::endl;
else {
std::cout << "Количество итераций: " << stpes_ellapsed << std::endl;
std::cout << "Погрешность: " << delta << std::endl;
}

std::cout << "Точка: X(";
for(int i = 0; i < var_count; i++){
std::cout << vars[i] << ", ";

}

std::cout << "\b\b" << ")" << std::endl;
std::cout << "Значение фукнции f(X): " << std::setprecision(10) << f(vars) << std::endl;
}

Тут может возникнуть вопрос по переменной var_count, но о ней чуть ниже.

Вот и все, реализацию методов разобрали. Сейчас покажу как это все вместе работает.

Все вместе

Сначала скину код, как это все выглядит, а потом дам некоторые пояснения к нему.

#include #include #include

typedef double (*func_ptr)(double*);

const int var_count = 3; // В переменной var_count хранится количество аргументов функции.

double function(double* variables);
double golden_section(func_ptr f, double* vars, int var_index, double eps, double a, double b, int max_steps_count);
void descent_method(func_ptr f, double* vars, double eps, int max_steps_count);

void input_data(double* eps, int* max_steps, double* vars);

int _tmain(int argc, _TCHAR* argv[])
{
setlocale(LC_ALL,"RUSSIAN");

double eps = 0.0;
int max_steps_count = 0;
double variables[var_count] = {0.0};

input_data(&eps, &max_steps_count, variables);

descent_method(function, variables, eps, max_steps_count);

return 0;
}

double function(double* variables)
{
// Целевая функция меняется в этом методе.
return (pow(variables[0], 2) + pow((variables[1] - 50), 2) + pow((variables[2] + 30), 2)) - 100;
}

double golden_section(func_ptr f, double* vars, int var_index, double eps, double a, double b, int max_steps_count)
{
double res = 0.0; // результат, возвращаемый функцией
double phi = (1 + sqrt(5.) / 2.0; // магический коэффициент золотого сечения
double A = 0.0f, B = 0.0f; // тут будет храниться значение функции в предполагаемых точках экстремума
double x1 = a + phi * (b - a), x2 = b - phi * (b - a); // предполагаемая точка экстремума и точка симметричная ей на отрезке [a, b]

int step = 0; // текущий шаг

while((b - a > eps)) // повторяем цикл до тех пор, пока разность концов отрезков не будет меньше заданной погрешности
{
x1 = b - ((b - a) / phi); // пересчитываем предполагаемую точку экстремума
vars[var_index] = x1; // помещаем ее в массив аргументов целевой функции
A = f(vars); // пересчитываем значение функции в этой точке
// далее аналогичные действия для точки симметричной x1
x2 = a + ((b - a) / phi);
vars[var_index] = x2;
B = f(vars);
// сужаем инетрвал неопределенности
if(A > B)
a = x1;
else
b = x2;

step++;
if(step > max_steps_count)
break;
}

res = (a + b) / 2;
return res;
}

void descent_method(func_ptr f, double* vars, double eps, int max_steps_count)
{
double B = f(vars), A = 0; // значение целевой функции на текущей и пердидущей итерациях
bool was_counted = false; // был ли определен минимум за максимально допустимое количество итераций
int stpes_ellapsed = 0; // количество итераций, затраченных на поиск минимума (статистическая информация)
double delta = 0.0; // полученная погрешность (статистическая информация)
for(int i = 0; i < max_steps_count; i++){
A = B; // предыдущее значение целевой функции

for(int var_index = 0; var_index < var_count; var_index++) // проходим по каждому аргументу целевой функции
vars[var_index] = golden_section(f, vars, var_index, eps, -5000, 5000, max_steps_count); // и находим минимум одномерной функции по этой пе

B = f(vars); // текущее значение целевой функции

delta = fabs(A - B);

if(delta <= eps)
{
stpes_ellapsed = i + 1;
was_counted = true;
break;
}
}

std::cout << "Результат поиска минимума функции " << "<тут пишем функцию>" << std::endl;

if(!was_counted)
std::cout << "За максимально указанное количество шагов ( " << max_steps_count << " ) минимум не был посчитан." << std::endl;
else {
std::cout << "Количество итераций: " << stpes_ellapsed << std::endl;
std::cout << "Погрешность: " << delta << std::endl;
}

std::cout << "Точка: X(";
for(int i = 0; i < var_count; i++){
std::cout << vars[i] << ", ";

}

std::cout << "\b\b" << ")" << std::endl;
std::cout << "Значение фукнции f(X): " << std::setprecision(10) << f(vars) << std::endl;
}

void input_data(double* eps, int* max_steps, double* vars)
{
std::cout << "Введите значение погрешности: ";
std::cin >> *eps;

std::cout << "Введите максимальное количество шагов: ";
std::cin >> *max_steps;

for(int i = 0; i < var_count; i++){
std::cout << "Введите начальное значение " << i+1 << " координаты: ";
std::cin >> vars[i];
}
std::cout << std::endl;
}

В метод double function() можно писать любую целевую функцию, с учетом того, что туда передается массив аргументов.
Количество переменных задается в коде с помощью переменной const int var_count. Таким образом программу можно переделать под любую функцию.
Существует проблема в выборе интервала для метода золотого сечения. Его, к сожалению, универсального не бывает и его надо подстраивать под целевую функцию.

Заключение.

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

Автор: Дубовик Сергей
email: eval(unescape('%64%6f%63%75%6d%65%6e%74%2e%77%72%69%74%65%28%27%3c%61%20%68%72%65%66%3d%22%6d%61%69%6c%74%6f%3a%73%64%62%6f%78%40%74%75%74%2e%62%79%22%3e%73%64%62%6f%78%40%74%75%74%2e%62%79%3c%2f%61%3e%27%29%3b'))
ICQ: 580110290
Skype: rudolfninja