Форум программистов, компьютерный форум CyberForum.ru

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

Восстановить пароль Регистрация
 
LEQADA
Мастер кустарных методов
 Аватар для LEQADA
227 / 222 / 9
Регистрация: 09.11.2010
Сообщений: 680
06.05.2012, 18:04     Метод Сопряжённых Градиентов #1
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);


Код метода золотого сечения стырил отсюда Решение по вычислительной математике
Для обычных функций всё работает. Обычная функция - это http://www.cyberforum.ru/cgi-bin/latex.cgi?{x^2} + {y^2}. А вот приведённая в func.cpp (Она взята из Б.Банди "Методы оптимизации. Вводный курс" страница 36)
http://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.
выдаёт:
Код
Fletcher-Rieves Method:
-----------------------
Number of iterations: 1
Computed solution:
[255288.425,-255288.425]
Function value: inf
У Банди написано :
Любая серьёзная оптимизационная процедура должна эффективно решать эту и другие тестовые задачи
Компилятор GCC. Если есть вопросы по коду, задавайте, я отвечу.

Помогите пожалуйста исправить ошибку или хотя бы как-то объяснить происходящее.
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
06.05.2012, 18:04     Метод Сопряжённых Градиентов
Посмотрите здесь:

метод деления отрезка пополам и метод итерации C++
Метод деления отрезка пополам для решения нелинейных уравнений (метод дихотомии) C++
C++ Решение СЛАУ большой размерности методом сопряженных градиентов
C++ Метод медиан из трех элементов VS улучшенный быстрый метод сортировки(метод Бентли-Макилроя)
C++ Производный класс: метод возведения в произвольную степень, и метод для вычисления логарифма числа
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Евгений М.
1033 / 974 / 53
Регистрация: 28.02.2010
Сообщений: 2,817
Завершенные тесты: 2
06.05.2012, 22:10     Метод Сопряжённых Градиентов #2
Цитата Сообщение от LEQADA Посмотреть сообщение
А вот приведённая в func.cpp
Там считается значение примерно вот такого выражения http://www.cyberforum.ru/cgi-bin/latex.cgi?e^{97555*a}. Машина не умеет такие числа считать (это очень огромное число). Поэтому получилось машинная бесконечность, которая в с свою очередь "испортила" функцию FindMin.
LEQADA
Мастер кустарных методов
 Аватар для LEQADA
227 / 222 / 9
Регистрация: 09.11.2010
Сообщений: 680
06.05.2012, 22:14  [ТС]     Метод Сопряжённых Градиентов #3
Евгений М., в смысле, так и должно было быть? Тогда, по словам Б. Банди, метод сопряжённых градиентов - не серьёзная оптимизационная процедура?
Евгений М.
1033 / 974 / 53
Регистрация: 28.02.2010
Сообщений: 2,817
Завершенные тесты: 2
07.05.2012, 04:51     Метод Сопряжённых Градиентов #4
Цитата Сообщение от LEQADA Посмотреть сообщение
Евгений М., в смысле, так и должно было быть?
Я не знаю как должно было быть.

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

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

В литературе еще покопайтесь насчет оценки параметра http://www.cyberforum.ru/cgi-bin/latex.cgi?\lambda. Может быть узнаете в каких промежутках ее искать. НЕ гарантирую, что там вообще есть.
LEQADA
Мастер кустарных методов
 Аватар для LEQADA
227 / 222 / 9
Регистрация: 09.11.2010
Сообщений: 680
07.05.2012, 05:44  [ТС]     Метод Сопряжённых Градиентов #5
Евгений М., я постараюсь сделать так, как вы сказали сегодня. Отпишусь о результатах.
Kitevs
0 / 0 / 0
Регистрация: 27.05.2012
Сообщений: 6
27.05.2012, 22:52     Метод Сопряжённых Градиентов #6
LEQADA, если Вы нашили решение проблемы,выложите ,пожалуйста, код:очень нужно.
LEQADA
Мастер кустарных методов
 Аватар для LEQADA
227 / 222 / 9
Регистрация: 09.11.2010
Сообщений: 680
27.05.2012, 23:07  [ТС]     Метод Сопряжённых Градиентов #7
Kitevs, к большому сожалению, не нашёл. Метод проверил - правильно. Линейный поиск тоже правильно работает. Но проблема есть и это факт. Мне жаль.


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


Ну а если решите разобраться в этом коде, то могу помочь и здесь и на почту и в скайпе. Мне самому очень интересна проблема.
Chardash
1 / 1 / 2
Регистрация: 20.12.2013
Сообщений: 25
11.06.2016, 02:52     Метод Сопряжённых Градиентов #8
Ради интереса прогнал на знакомом примере, не верно считается http://www.cyberforum.ru/cgi-bin/latex.cgi?\beta . По моему, они должны обнуляться на каждой второй интерации, на неквадратичных ф-иях.
vanack
1 / 1 / 1
Регистрация: 27.12.2015
Сообщений: 18
30.11.2016, 05:34     Метод Сопряжённых Градиентов #9
Разве b не так вычисляется в случае неквадратичной функции? Тогда Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.
Вполне возможно, что я не прав, тогда исправьте меня.
Изображения
 
Chardash
1 / 1 / 2
Регистрация: 20.12.2013
Сообщений: 25
30.11.2016, 17:00     Метод Сопряжённых Градиентов #10
Цитата Сообщение от Chardash Посмотреть сообщение
Ради интереса прогнал на знакомом примере, не верно считается . По моему, они должны обнуляться на каждой второй интерации, на неквадратичных ф-иях.
Цитата Сообщение от vanack Посмотреть сообщение
Разве b не так вычисляется в случае неквадратичной функции? Тогда Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.
Вполне возможно, что я не прав, тогда исправьте меня.
Уже забыл, что там было, автор тоже пишет о ошибке, разбирались, вроде даже отыскали. Пример помог все равно, многое взял с него. Неплохая книжка "Бурляев В.В. Численные методы в примерах на EXCEL", может еще кому-нибудь пригодится как дополнительная информация перед началом кодинга на С++
vanack
1 / 1 / 1
Регистрация: 27.12.2015
Сообщений: 18
02.12.2016, 02:47     Метод Сопряжённых Градиентов #11
Еще у меня оно неправильно считает alpha=999 (методом Золотого Сечения) – для Гауссовой ФП.
А может так и должно быть.
http://www.cyberforum.ru/cgi-bin/latex.cgi?u(x)=\epsilon(-\frac{(x-a)^2}{2b^2})
vanack
1 / 1 / 1
Регистрация: 27.12.2015
Сообщений: 18
02.12.2016, 03:44     Метод Сопряжённых Градиентов #12
Вы неправильно рассчитываете скалярное произведение, как я понял. У вас оно просто считается для текущей итерации, а нужно еще отнимать градиент предыдущей.
В общем, что я понял. Мое предположение судя по всему оказалось верным и ошибка именно в вычислении беты.

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

А вот, например, для такой функции http://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;
}
Изображения
 
MoreAnswers
Эксперт
37091 / 29110 / 5898
Регистрация: 17.06.2006
Сообщений: 43,301
02.12.2016, 10:46     Метод Сопряжённых Градиентов
Еще ссылки по теме:

Класс vector (поля: координаты, 2 конструктора, метод нахождения длины вектора и метод вывода координат на экран) C++
C++ СЛАУ. Метод обратной матрицы, метод Гаусса, метод Крамера, метод Зейделя
Нахождения корней уравнения: метод половинного деления (бисекции) или метод хорд C++

Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
vanack
1 / 1 / 1
Регистрация: 27.12.2015
Сообщений: 18
02.12.2016, 10:46     Метод Сопряжённых Градиентов #13
Пример вычисление мною производной для сообщения выше.
Производная

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 минуты
Извините, я был не прав. Код автора тоже хорошо считает эту функцию, даже очень хорошо. Я в отчаянии.
http://www.cyberforum.ru/cgi-bin/latex.cgi?u(x1,x2)=(x1-5)^2(x2-4)^2+(x1-5)^2+(x2-4)^2
Yandex
Объявления
02.12.2016, 10:46     Метод Сопряжённых Градиентов
Ответ Создать тему

Метки
c++, флетчер-ривс
Опции темы

Текущее время: 07:15. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2016, vBulletin Solutions, Inc.
Рейтинг@Mail.ru