gradsearch1.cpp

#include <iostream>
#include <math.h>
#include <conio.h>
#include <iomanip>
#include <fstream>
#include <stdlib.h>
using namespace std;
 
double f(double x, double y)
{
    return  y*y + cos(x + 1.23*y);
};
 
double gradX(double x, double y)
{
    return -sin(x + 1.23*y);
};
 
double gradY(double x, double y)
{
    return  2*y - 1.23*sin(x + 1.23*y);
};
 
double norm(double x, double y)
{
    return sqrt((-sin(x + 1.23*y))*(-sin(x + 1.23*y)) + (2*y - 1.23*sin(x + 1.23*y))*(2*y - 1.23*sin(x + 1.23*y)));
};
 
 
void output(double x, double y)
    {
        cout << setw(4) << fixed << setprecision(5) << x;
        cout << setw(14) << fixed << setprecision(5) << y;
        cout << setw(14) << fixed << setprecision(5) << f(x,y) << endl;
    };
 
void main()
{
    double x0 = -0.4;
    double y0 = 1.9;
    double h = 0.002;
    double M = 6;
    double x = x0;
    double y = y0;
    output(x,y);
    for(int k = 0; k < 4; ++k)
    {   
        double gradx = gradX(x,y);
        double grady = gradY(x,y);
        double Norm = norm(x,y);
        x = x - h*gradx/Norm;
        y = y - h*grady/Norm;
        
        output(x,y);
    };
    _getch();
}

gradsearch2.cpp

#include <vector>
#include <iostream>
#include <Math.h>
#include <string>
#include <sstream>
#include <conio.h>
 
using namespace std;
 
typedef vector<double> DataList;
 
void InitData();
void GradSearch(DataList &p, double sigma,double epsilon);
void QMin(DataList &de_dxi, DataList &p, double epsilon, double sigma);
void getDfDx(DataList &de_dxi, DataList &p);
string CreateResultString(DataList &pMin,double yMin);
double getFunc(DataList &p);
string stringify(double x);
 
DataList p1;
DataList p2;
DataList pMin;
DataList de_dxi;
 
double yMin;
double err;
double z0;
double h;
int N;
 
int j;
 
int main()
{
    N = 2;
 
    InitData();
 
    DataList p;
    p.push_back(0.99);
    p.push_back(1.01);
        double sigma = 0.0000000001;
        double epsilon = 0.0000000001;
 
    GradSearch(p, sigma, epsilon);
 
    getch();
    return 0;
}
 
void GradSearch(DataList &p, double sigma,double epsilon)
{
    int max = 60;
    h = 1;
    err = 1;
    int count = 0;
 
    while(count < max && (h>sigma ||err > epsilon))
    {
        getDfDx(de_dxi, p);
        QMin(de_dxi, p, epsilon, sigma);
 
        for(int i = 0; i<N; i++)
            p.at(i) = pMin.at(i);
 
        z0 = yMin;
        count = count + j + 1;
 
        string iterResult = CreateResultString(pMin, yMin);
 
        cout<<CreateResultString(pMin, yMin)<<endl;
    }
    cout<<endl<<"Минимум функции: "<<"-4*x + x*x - y - x*y + y*y"<<endl;
    cout<<CreateResultString(pMin, yMin)<<endl;
}
 
void QMin(DataList &de_dxi,DataList &p, double epsilon, double sigma)
{
    int cond = 0;
    int jmax = 60;
 
    z0 = getFunc(p);
 
    for(int i = 0; i<N; i++)
        p1.at(i) = p.at(i) + h * de_dxi.at(i);
 
    double y1 = getFunc(p1);
 
    for(int i = 0; i<N; i++)
        p2.at(i) = p.at(i) + 2 * h * de_dxi.at(i);
 
    double y2 = getFunc(p2);
 
    j = 0;
 
    while (j<jmax && cond == 0)
    {
        if (z0<=y1)
        {
            for(int i = 0; i<N; i++)
                p2.at(i) = p1.at(i);
 
            y2 = y1;
            h = h / 2;
 
            for(int i = 0; i<N; i++)
                p1.at(i) = p.at(i) + h * de_dxi.at(i);
 
            y1 = getFunc(p1);
        }
        else if (y2 < y1)
        {
            for(int i = 0; i<N; i++)
                p1.at(i) = p2.at(i);
 
            y1 = y2;
 
            h = h*2;
 
            for(int i = 0; i<N; i++)
                p2.at(i) = p.at(i) + 2 * h * de_dxi.at(i);
 
            y2 = getFunc(p2);
        }
        else
        {
            cond = -1;
        }
 
        j = j+1;
        if (h < sigma)
            cond = 1;
    }
 
    double hMin = (h/2)* (4 * y1 - 3* z0 - y2) / (2* y1 - z0 - y2);
 
    for (int i = 0; i< N; i++)
    {
        pMin.at(i) = p.at(i) + hMin * de_dxi.at(i);
    }
 
    yMin = getFunc(pMin);
 
    double h0 = fabs(hMin);
    double h1 = fabs(hMin - h);
    double h2 = fabs(hMin - 2* h);
 
    if (h0 < h)
        h = h0;
    if (h1 < h)
        h = h1;
    if (h2 < h)
        h = h2;
    if (h == 0)
        h = hMin;
    if (h < sigma)
        cond = 1;
    double e0 = fabs(z0 - yMin);
    double e1 = fabs(y1 - yMin);
    double e2 = fabs(y2 - yMin);
 
    if (e0 != 0 && e0 < err)
        err = e0;
    if (e1 != 0 && e1 < err)
        err = e1;
    if (e2 != 0 && e2 < err)
        err = e2;
    if (e0 == 0 && e1 == 0 && e2 == 0)
        err = 0;
    if (err < epsilon)
        cond = 2;
}
 
double getFunc(DataList &p)
{
    double x = p.at(0);
    double y = p.at(1);
 
    double result = -4 * x + x*x - y - x * y + y * y;
 
    return result;
}
 
void getDfDx(DataList & de_dxi, DataList &p)
{
    double x = p.at(0);
    double y = p.at(1);
 
 
    double dfDx = -4+2*x-y;
    double dfDy = -1-x+2*y;
 
    double norm = sqrt(dfDx*dfDx + dfDy*dfDy);
 
    dfDx = -dfDx/norm;
    dfDy = -dfDy/norm;
 
    de_dxi.at(0) = dfDx;
    de_dxi.at(1) = dfDy;
}
 
void InitData()
{
    for(int i = 0; i<N; i++)
    {
        p1.push_back(0);
 
        p2.push_back(0);
 
        pMin.push_back(0);
 
        de_dxi.push_back(0);
    }
}
 
string stringify(double x)
{
    ostringstream o;
    if (!(o << x))
        return 0;
    return o.str();
}
 
string CreateResultString(DataList &pMin,double yMin)
{
    string resultStr = "f[";
 
    for(int i = 0; i<N; i++)
    {
        if (i != 0)
            resultStr += ",";
        resultStr += stringify(pMin.at(i));
    }
 
    resultStr += "] = " + stringify(yMin);
 
    return resultStr;
}

gradsearch3.cpp

//Программа для поиска минимума функции методом градиентного спуска
#include<iostream>
#include<vector>
#include<math.h>

#define METHOD (2)
//выбор метода определения длины шага
//1 - метод градиентного спуска с постоянным шагом
//2 - Градиентный метод с дроблением шага
//3 - Метод наискорейшего спуска

//константа для метода градиентного спуска с постоянным шагом
//начальное значение eps для метода с дроблением шага
#define LAMBDA_CONSTANT (1)

//параметры для метода с дроблением шага
#define DELTA_FOR_METHOD2 (0.95)
#define EPS_FOR_METHOD2 (0.1)

//максимально возможное число итераций
#define NUMBER_OF_ITERATIONS (100000)

//eps для критерия останова
#define EPS (1e-5)

//критерий останова
#define OSTANOV (2)

using namespace std;

vector<double> goldensectionoptimize(vector<double> x,double a, double b, int n);
double f(vector<double> x);
vector<double> GradientDescent(int N,vector<double> x0,int&Iterations);




double f(vector<double> x)
//функция минимум которой ищут
{
	//int l=x.size();
	//return (100*(x[1]-x[0]*x[0])*(x[1]-x[0]*x[0])+(1-x[0])*(1-x[0]));
	return 10*x[0]*x[0]+x[1]*x[1];
}

vector<double> GradientF(vector<double> x)
//градиент исследуемой функции
{
	vector<double> tmp;
	//tmp.push_back(-2*(1-x[0])-400*(x[1]-x[0]*x[0])*x[0]);
	//tmp.push_back(200*x[1]);
	tmp.push_back(20*x[0]);
	tmp.push_back(2*x[1]);
	return tmp;
}


vector<double> GradientDescent(int N,vector<double> x0,int&Iterations)
//minimizes N-dimensional function f; x0 - start point
{
	vector <double> old,cur_x=x0,gr;
	double s,Lambda;
	int j,i;

	Lambda=LAMBDA_CONSTANT;
	
	for (Iterations=1;Iterations<=NUMBER_OF_ITERATIONS;Iterations++)
	{
		//save old value
		old=cur_x;
		//compute gradient
        gr=GradientF(cur_x);
		

		if (METHOD==1)
		{

			Lambda=LAMBDA_CONSTANT;
			//вычисляем новое значение
			for(j=0;j<old.size();j++)
				cur_x[j]=cur_x[j]-Lambda*gr[j];
		}

		if (METHOD==2)
		{
			//вычисляем новое значение
			for(j=0;j<old.size();j++)
				cur_x[j]=cur_x[j]-Lambda*gr[j];

			//вычисляем квадрат нормы градиента
			s=0;
			for(j=0;j<old.size();j++)
				s+=gr[j]*gr[j];


			while (f(cur_x)>f(old)-EPS_FOR_METHOD2*Lambda*s)
			{
				Lambda=Lambda*DELTA_FOR_METHOD2;
				//пересчет значения Лямда
				cur_x=old;
				for(j=0;j<old.size();j++)
					cur_x[j]=cur_x[j]-Lambda*gr[j];
			}
			
		}

		if (METHOD==3)
		{
			cur_x=goldensectionoptimize(cur_x,-10,10,100);
		}
		
		if(OSTANOV==1)
		{
			//условие останова 1
			s=0;
			for(j=0;j<old.size();j++)
				s+=(old[j]-cur_x[j])*(old[j]-cur_x[j]);
			s=sqrt(s);
			if (s<EPS)
				return cur_x;
		}

		if(OSTANOV==2)
		{
			//условие останова 2
			s=fabs(f(cur_x)-f(old));
			if (s<EPS)
				return	cur_x;
		}

	}

	return cur_x;
}

int main()
{
	vector<double> x;
	x.push_back(10);
	x.push_back(10);
	int i,Iteration;
	vector<double> ans=GradientDescent(2,x,Iteration);
	cout<<"Value: "<<f(ans)<<endl;
	cout<<"Point: ";
	for(i=0;i<ans.size();i++)
		cout<<ans[i]<<' ';
	cout<<endl<<"Number of iterations:"<<Iteration<<endl;
	return 0;
}

vector<double> goldensectionoptimize(vector<double> x,double a, double b, int n)
//метод золотого сечения одномерной оптимизации функции f
//вектор переменных x на отрезке [a,b], n итераций
//взят с сайта http://alglib.sources.ru/
{
	vector<double> tmp=x;
	int i,j;
    double s1;
    double s2;
    double u1;
    double u2;
    double fu1;
    double fu2;
	vector<double> GF=GradientF(x);

    s1 = (3-sqrt(double(5)))/2;
    s2 = (sqrt(double(5))-1)/2;
    u1 = a+s1*(b-a);
    u2 = a+s2*(b-a);
	
	for(j=0;j<x.size();j++)
		tmp[j]=x[j]+u1*GF[j];
	fu1 = f(tmp);
	
	for(j=0;j<x.size();j++)
		tmp[j]=x[j]+u2*GF[j];
	fu2 = f(tmp);

    for(i = 1; i <= n; i++)
    {
        if( fu1<=fu2 )
        {
            b = u2;
            u2 = u1;
            fu2 = fu1;
            u1 = a+s1*(b-a);
			for(j=0;j<x.size();j++)
				tmp[j]=x[j]+u1*GF[j];
	        fu1 = f(tmp);
        }
        else
        {
            a = u1;
            u1 = u2;
            fu1 = fu2;
            u2 = a+s2*(b-a);
			for(j=0;j<x.size();j++)
				tmp[j]=x[j]+u2*GF[j];
			fu2 = f(tmp);
        }
    }
	for(j=0;j<x.size();j++)
		tmp[j]=x[j]+u1*GF[j];
	fu1 = f(tmp);
    for(j=0;j<x.size();j++)
		tmp[j]=x[j]+u2*GF[j];
	fu2 = f(tmp);

	if (fu1<fu2)
		for(j=0;j<x.size();j++)
			tmp[j]=x[j]+u1*GF[j];
	else 
		for(j=0;j<x.size();j++)
			tmp[j]=x[j]+u2*GF[j];
	
	return tmp;		
}