Метод Сопряжённых Градиентов

LEQADA · Регистрация: 09.11.2010

Студворк — интернет-сервис помощи студентам

main.cpp

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#include <iostream>
#include <cstdlib>  
#include <cmath>
#include "sol.h"
using namespace std;
 
// Вывод результата на экран
void PrintSolution(double *x, double val, int numIter) 
{
    cout << "-----------------------" << endl;
    cout << "Number of iterations: " << numIter << endl;
    cout<<"Computed solution: "<<endl<<"["<<x[0]<<","<<x[1]<<"]"<<endl;
    cout<<"Function value: "<<val<<endl<<endl;
}
 
// Вычисление скалярного произведения
double inner_prod(double *x, double *y) {return x[1]*y[1] + x[0]*y[0];}
 
int main()
{
    cout.precision(9);
    FletcherRievesMethod();
    return 0;
}

frm.cpp

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#include <iostream>
#include <cstdlib>  
#include <cmath>
#include "sol.h"
using namespace std;
 
double FletcherRievesMethod()
{
    double eps = 1e-3;
    const double EPS = 1e-5;
    //Начальное приближение
    double x[2]={0,0};
    double curVal = F1(x);
    double prevVal = curVal;
    double p[2] ;
    p[0]= - GradF1(x);
    p[1]= - GradF2(x);
    
    double gradSquare = inner_prod(p, p);
 
    int numIter = 0;
    do
    {
        numIter++;
        double alpha, beta, newGradSquare;
        double newGrad[2];
 
 
        //Ищем минимум F1(x + alpha * p) с помощью метода одномерной оптимизации
        alpha = FindMin(x, p);
        x[0] = x[0] + alpha * p[0];
        x[1] = x[1] + alpha * p[1];
        newGrad [0] = - GradF1(x);
        newGrad [1] = - GradF2(x);
        
        newGradSquare = inner_prod(newGrad, newGrad);
        
        if (numIter % (10) == 0) beta = 0; //Обновление
        else 
        beta = newGradSquare / gradSquare; //Используем метод Флетчера - Ривса
        p[0] = newGrad[0] + beta * p[0];
        p[1] = newGrad[1] + beta * p[1];
 
        prevVal = curVal;
        curVal = F1(x);
 
        gradSquare = newGradSquare;
    }
    while (gradSquare > eps);
 
    cout << "Fletcher-Rieves Method: " << endl; 
    PrintSolution(x, F1(x), numIter);
    return 0;
 
}

func.cpp

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <cmath>
// Сама функция
double F1(double *x)
{
    double a=0.1, eps=0.00001,res=0;
    while (fabs(a-1)>eps) {res = res + ((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a)))*((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a))); a=a+0.1;}
    return res;
}
// Производная по первой переменной
double GradF1(double *x)
{
    double a=0.1, eps=0.00001,res=0;
    while (fabs(a-1)>eps) {res = res - 2*a*exp(-a*x[0]) * ((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a)))*((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a))); a=a+0.1;}
    return res;
}
// Производная по второй переменной
double GradF2(double *x)
{
    double a=0.1, eps=0.00001,res=0;
    while (fabs(a-1)>eps) {res = res + 2*a*exp(-a*x[0]) * ((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a)))*((exp(-a*x[0]) - exp(-a*x[1])) - (exp(-a) - exp(-10*a))); a=a+0.1;}
    return res;
}

golden.cpp

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include <iostream>
#include <cstdlib>  
#include <cmath>
#include "sol.h"
 
//Одномерная оптимизация - применим метод золотого сечения
double FindMin(double *s, double *p)
{
    const double eps = 1e-8;
    const double tay = 1.618;
    double a = 0;
    double b = 1e5;
    double x0, x1, xf1, xf2;
    x0 = b - (b - a) / tay; // Расчитываем точки деления
    x1 = a + (b - a) / tay;         //
P:
    double t1[2],t2[2];
    t1[0] = s[0] + x0 * p[0];
    t1[1] = s[1] + x0 * p[1];
    t2[0] = s[0] + x1 * p[0];
    t2[1] = s[1] + x1 * p[1];
    xf1 = F1(t1); // Расчитываем в точках деления значение целевой функции
    xf2 = F1(t2); //
    if(xf1 >= xf2) 
    {
        a = x0;
        x0 = x1;
        t2[0] = s[0] + x1 * p[0];
        t2[1] = s[1] + x1 * p[1];
        xf1 = F1(t2);
        t2[0] = s[0] + x1 * p[0];
        t2[1] = s[1] + x1 * p[1];
        x1 = a + (b - a) / tay;
        xf2 = F1(t2);
    }
    else
    {
        b = x1;
        x1 = x0;
        xf2 = xf1;
        x0 = b - (b - a) / tay;
        t1[0] = s[0] + x0 * p[0];
        t1[1] = s[1] + x0 * p[1];
        xf1 = F1(t1);
    }
    if(fabs(b - a) < eps) 
    {
        return (a + b) / 2;
    }
    else
    goto P;
}

sol.h

C++
1
2
3
4
5
6
7
8
double F1(double *x);
double GradF1(double *x);
double GradF2(double *x);
double GradFQuad(double *x);
double FindMin(double *s, double *p);
void PrintSolution(double *x, double val, int numIter);
double FletcherRievesMethod();
double inner_prod(double *x, double *y);

Код метода золотого сечения стырил отсюда Решение по вычислительной математике
Для обычных функций всё работает. Обычная функция - это $https://www.cyberforum.ru/cgi-bin/latex.cgi?{x^2} + {y^2}$ . А вот приведённая в func.cpp (Она взята из Б.Банди "Методы оптимизации. Вводный курс" страница 36)
$https://www.cyberforum.ru/cgi-bin/latex.cgi?f(x,y) = \sum\limits_a {{{\left( {\left( {{e^{ - ax}} - {e^{ - ay}}} \right) - \left( {{e^{ - a}} - {e^{ - 10a}}} \right)} \right)}^2}}$
сумма по всем а от 0.1 до 1 с шагом 0.1.
выдаёт:

Code
1
2
3
4
5
6
Fletcher-Rieves Method:
-----------------------
Number of iterations: 1
Computed solution:
[255288.425,-255288.425]
Function value: inf

У Банди написано :

Любая серьёзная оптимизационная процедура должна эффективно решать эту и другие тестовые задачи

Компилятор GCC. Если есть вопросы по коду, задавайте, я отвечу.

Помогите пожалуйста исправить ошибку или хотя бы как-то объяснить происходящее.

Евгений М. · 06.05.2012, 22:10

Сообщение от LEQADA

А вот приведённая в func.cpp

Там считается значение примерно вот такого выражения $https://www.cyberforum.ru/cgi-bin/latex.cgi?e^{97555*a}$ . Машина не умеет такие числа считать (это очень огромное число). Поэтому получилось машинная бесконечность, которая в с свою очередь "испортила" функцию FindMin.

LEQADA · 06.05.2012, 22:14 **[ТС]**

Евгений М., в смысле, так и должно было быть? Тогда, по словам Б. Банди, метод сопряжённых градиентов - не серьёзная оптимизационная процедура?

Евгений М. · 07.05.2012, 04:51

Сообщение от LEQADA

Евгений М., в смысле, так и должно было быть?

Я не знаю как должно было быть.

Прежде всего я предлагаю узнать при каких значениях x,y искомая функция считается (устроить табулирование функции по двум аргументам с шагом скажем 1-10 для каждого аргумента). При x=-10000 или y=-10000 искомая функция считаться не будет.

Допустим получилось что при -100 < x,y < 100 функция вычисляется. Тогда для золотого сечения ищите минимумы на таком промежутке.

В литературе еще покопайтесь насчет оценки параметра $https://www.cyberforum.ru/cgi-bin/latex.cgi?\lambda$ . Может быть узнаете в каких промежутках ее искать. НЕ гарантирую, что там вообще есть.

LEQADA · 07.05.2012, 05:44 **[ТС]**

Евгений М., я постараюсь сделать так, как вы сказали сегодня. Отпишусь о результатах.

@Kitevs · 27.05.2012, 22:52

LEQADA, если Вы нашили решение проблемы,выложите ,пожалуйста, код:очень нужно.

LEQADA · 27.05.2012, 23:07 **[ТС]**

Kitevs, к большому сожалению, не нашёл. Метод проверил - правильно. Линейный поиск тоже правильно работает. Но проблема есть и это факт. Мне жаль.

попробуйте прочитать код этого метода из книжки Б.Банди( см. шапку). Там написан код на каком-то алгоритмическом языке. Линейный поиск там реализуется с помощью метода кубической интерполяции. Я там не смог разобраться.

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

@Chardash · 11.06.2016, 02:52

Ради интереса прогнал на знакомом примере, не верно считается $https://www.cyberforum.ru/cgi-bin/latex.cgi?\beta$ . По моему, они должны обнуляться на каждой второй интерации, на неквадратичных ф-иях.

@vanack · 30.11.2016, 05:34

Разве b не так вычисляется в случае неквадратичной функции? Тогда Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.
Вполне возможно, что я не прав, тогда исправьте меня.

@Chardash · 30.11.2016, 17:00

Сообщение от Chardash

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

Сообщение от vanack

Разве b не так вычисляется в случае неквадратичной функции? Тогда Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.
Вполне возможно, что я не прав, тогда исправьте меня.

Уже забыл, что там было, автор тоже пишет о ошибке, разбирались, вроде даже отыскали. Пример помог все равно, многое взял с него. Неплохая книжка "Бурляев В.В. Численные методы в примерах на EXCEL", может еще кому-нибудь пригодится как дополнительная информация перед началом кодинга на С++

@vanack · 02.12.2016, 02:47

Еще у меня оно неправильно считает alpha=999 (методом Золотого Сечения) – для Гауссовой ФП.
А может так и должно быть.
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u(x)=\epsilon(-\frac{(x-a)^2}{2b^2})$

@vanack · 02.12.2016, 03:44

Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.

В общем, что я понял. Мое предположение судя по всему оказалось верным и ошибка именно в вычислении беты.

Исправленный мною вариант для Гауссовой ФП $https://www.cyberforum.ru/cgi-bin/latex.cgi?u(x)=\epsilon (-\frac{(x-a)^2}{2b^2})$ действительно считает неправильно! Не знаю почему.

А вот, например, для такой функции $https://www.cyberforum.ru/cgi-bin/latex.cgi?u(x1,x2)=(x1-5)^2(x2-4)^2+(x1-5)^2+(x2-4)^2$ , где точки минимума – (5, 4), оно выдает что то похожее на правду (4,8; 3,6). Но такая погрешность скорее всего из-за того, что я вычислял частную производную через eps.
Если использовать исходный код из этой темы, то оно выдает полнейший бред – (0,00057; 0,0083).

Пока как то так. Конечно, буду очень рад, если кто-нибудь укажет мне на мои ошибки. Т.к. очень нужно, чтобы метод заработал правильно.

Исправленный мною код

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
double GradMethod()
{
    double eps = 1e-3;
    const double EPS = 1e-5;
 
    //Начальное приближение
    double x[2]={1,1};
    double curVal = F1(x);
    double prevVal = curVal;
    double p[2], oldGrad[2], newGrad[2];
    p[0]= - GradF1(x);
    p[1]= - GradF2(x);
 
    oldGrad[0]= - GradF1(x);  //для вычисления градиента k-1
    oldGrad[1]= - GradF2(x);
 
    double gradSquare = inner_prod(p, p);
 
    int numIter = 0;
    do
    {
        numIter++;
        double alpha, beta, newGradSquare, GradSquarePr;
 
        //Ищем минимум F1(x + alpha * p) с помощью метода одномерной оптимизации
        alpha = FindMin(x, p);
        //ShowMessage(AnsiString(alpha));
 
        x[0] = x[0] + alpha * p[0];
        x[1] = x[1] + alpha * p[1];
 
        if (numIter != 1){
                oldGrad [0] = newGrad [0];
                oldGrad [1] = newGrad [1];
        }
        newGrad [0] = - GradF1(x);
        newGrad [1] = - GradF2(x);
 
        oldGrad [0] = newGrad [0] - oldGrad [0]; //для вычисления градиента k-1
        oldGrad [1] = newGrad [1] - oldGrad [1];
 
        newGradSquare = inner_prod(newGrad, newGrad);
        GradSquarePr = inner_prod(newGrad, oldGrad);  //вычисление скалярного произведения с градиентом k-1
 
        if (numIter == 1) beta = 0; //Обновление
        else
        beta = GradSquarePr / newGradSquare; //Используем метод Флетчера - Ривса (см. формулу во вложении)
 
        p[0] = newGrad[0] + beta * p[0];
        p[1] = newGrad[1] + beta * p[1];
 
        prevVal = curVal;
        curVal = F1(x);
 
        gradSquare = GradSquarePr;
    }
    while (gradSquare > eps);
 
    return 0;
}

@vanack · 02.12.2016, 10:46

Пример вычисление мною производной для сообщения выше.

Производная

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
double F1(double *x, double u=0) {
 
    return (x[0]-5)*(x[0]-5)*(x[1]-4)*(x[1]-4)+(x[0]-5)*(x[0]-5)+(x[1]-4)*(x[1]-4); //функция
}
//---------------------------------------------------------------------------
 
double F1dx(double *x, double u=0) {
    double eps=0.0001;
    return (x[0]+eps-5)*(x[0]+eps-5)*(x[1]-4)*(x[1]-4)+(x[0]+eps-5)*(x[0]+eps-5)+(x[1]-4)*(x[1]-4); //функция
}
//---------------------------------------------------------------------------
 
double F1dy(double *x, double u=0) {
    double eps=0.0001;
    return (x[0]-5)*(x[0]-5)*(x[1]+eps-4)*(x[1]+eps-4)+(x[0]-5)*(x[0]-5)+(x[1]+eps-4)*(x[1]+eps-4); //функция
}
//---------------------------------------------------------------------------
 
double GradF1(double *x, double u=0) {
    double eps=0.0001;
    return (F1dx(x, u)-F1(x, u))/eps;
}
//---------------------------------------------------------------------------
 
double GradF2(double *x, double u=0) {
    double eps=0.0001;
    return (F1dy(x, u)-F1(x, u))/eps;
}
//---------------------------------------------------------------------------

Добавлено через 6 часов 54 минуты
Извините, я был не прав. Код автора тоже хорошо считает эту функцию, даже очень хорошо. Я в отчаянии.
$https://www.cyberforum.ru/cgi-bin/latex.cgi?u(x1,x2)=(x1-5)^2(x2-4)^2+(x1-5)^2+(x2-4)^2$

Новые блоги и статьи Все статьи Все блоги /
Модульный подход на примере F# DevAlt 06.03.2026 В блоге дяди Боба наткнулся на такое определение: В этой книге («Подход, основанный на вариантах использования») Ивар утверждает, что архитектура программного обеспечения — это структуры,. . .	Управление камерой с помощью скрипта OrbitControls.js на Three.js: Вращение, зум и панорамирование 8Observer8 05.03.2026 Содержание блога Финальная демка в браузере работает на Desktop и мобильных браузерах. Итоговый код: orbit-controls-threejs-js. zip. Сканируйте QR-код на мобильном. Вращайте камеру одним пальцем,. . .	SDL3 для Web (WebAssembly): Синхронизация спрайтов SDL3 и тел Box2D 8Observer8 04.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-sync-physics-sprites-sdl3-c. zip На первой гифке отладочные линии отключены, а на второй включены:. . .	SDL3 для Web (WebAssembly): Идентификация объектов на Box2D v3 - использование userData и событий коллизий 8Observer8 02.03.2026 Содержание блога Финальная демка в браузере. Итоговый код: finish-collision-events-sdl3-c. zip Сканируйте QR-код на мобильном и вы увидите, что появится джойстик для управления главным героем. . . .
Реалии Hrethgir 01.03.2026 Нет, я не закончил до сих пор симулятор. Эта задача сложнее. Не получилось уйти в плавсостав, но оно и к лучшему, возможно. Точнее получалось - но сварщиком в палубную команду, а это значит, в моём. . .	Ритм жизни kumehtar 27.02.2026 Иногда приходится жить в ритме, где дел становится всё больше, а вовлечения в происходящее — всё меньше. Плотный график не даёт вниманию закрепиться ни на одном событии. Утро начинается с быстрых,. . .	SDL3 для Web (WebAssembly): Сборка библиотек: SDL3, Box2D, FreeType, SDL3_ttf, SDL3_mixer и SDL3_image из исходников с помощью CMake и Emscripten 8Observer8 27.02.2026 Недавно вышла версия 3. 4. 2 библиотеки SDL3. На странице официальной релиза доступны исходники, готовые DLL (для x86, x64, arm64), а также библиотеки для разработки под Android, MinGW и Visual Studio. . . .	SDL3 для Web (WebAssembly): Реализация движения на Box2D v3 - трение и коллизии с повёрнутыми стенами 8Observer8 20.02.2026 Содержание блога Box2D позволяет легко создать главного героя, который не проходит сквозь стены и перемещается с заданным трением о препятствия, которые можно располагать под углом, как верхнее. . .

LEQADA Мастер кустарных методов 232 / 227 / 17 Регистрация: 09.11.2010 Сообщений: 680
	06.05.2012, 22:14 [ТС]
	Евгений М., в смысле, так и должно было быть? Тогда, по словам Б. Банди, метод сопряжённых градиентов - не серьёзная оптимизационная процедура? 0

LEQADA Мастер кустарных методов 232 / 227 / 17 Регистрация: 09.11.2010 Сообщений: 680
	07.05.2012, 05:44 [ТС]
	Евгений М., я постараюсь сделать так, как вы сказали сегодня. Отпишусь о результатах. 0

@Kitevs 0 / 0 / 0 Регистрация: 27.05.2012 Сообщений: 6
	27.05.2012, 22:52
	LEQADA, если Вы нашили решение проблемы,выложите ,пожалуйста, код:очень нужно. 0

LEQADA Мастер кустарных методов 232 / 227 / 17 Регистрация: 09.11.2010 Сообщений: 680
	27.05.2012, 23:07 [ТС]
	Kitevs, к большому сожалению, не нашёл. Метод проверил - правильно. Линейный поиск тоже правильно работает. Но проблема есть и это факт. Мне жаль. попробуйте прочитать код этого метода из книжки Б.Банди( см. шапку). Там написан код на каком-то алгоритмическом языке. Линейный поиск там реализуется с помощью метода кубической интерполяции. Я там не смог разобраться. Ну а если решите разобраться в этом коде, то могу помочь и здесь и на почту и в скайпе. Мне самому очень интересна проблема. 1

@Chardash 41 / 2 / 2 Регистрация: 20.12.2013 Сообщений: 35
	11.06.2016, 02:52
	Ради интереса прогнал на знакомом примере, не верно считается $https://www.cyberforum.ru/cgi-bin/latex.cgi?\beta$ . По моему, они должны обнуляться на каждой второй интерации, на неквадратичных ф-иях. 0

@vanack 2 / 2 / 2 Регистрация: 27.12.2015 Сообщений: 18
	30.11.2016, 05:34
	Разве b не так вычисляется в случае неквадратичной функции? Тогда Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей. Вполне возможно, что я не прав, тогда исправьте меня. Изображения 0

@vanack 2 / 2 / 2 Регистрация: 27.12.2015 Сообщений: 18
	02.12.2016, 02:47
	Еще у меня оно неправильно считает alpha=999 (методом Золотого Сечения) – для Гауссовой ФП. А может так и должно быть. $https://www.cyberforum.ru/cgi-bin/latex.cgi?u(x)=\epsilon(-\frac{(x-a)^2}{2b^2})$ 0