梯度下降法，牛顿法，拟牛顿法

梯度下降法：根据泰勒展开式： 

                                            f（x+t） = f（x）+t*f’（x）

注意这里并不是严格相等的，我们这里取t为-f'(x)那么f（x+f'(x)）  = f(x)-f'(x)*f'(x)，通过迭代函数会收敛到一个局部最小值。

牛顿法：根据泰勒展开式： 

                                            f（x） = f(x0)+(x-x0)f'(x0)+(x-x0)(x-x0)f''(x0)                                                               (1)

注意这里也不是严格相等的，我们这里讨论的函数都是连续可微的，那么极值点必定是f'(x)  = 0，根据(1)我们可以得到如下式子：

                                           f'(x) = f'(x0)+(x-x0)f''(x0)

那么可以得到x = x0-f''(x)/f'(x);通过迭代可以逼近极值点。

注意到我们这里只是定义了一维的牛顿法，下面来讨论多维牛顿法。注意x=（x0，x1，x2...xn）是变量的集合

那么                                 f(x) = f(x0)+f'(x0)*(x-x0)+(x-x0)*H(x0)*(x-x0)                                                                （2）

这里的H（x0）是一个二阶偏导数的矩阵，长被称为海赛矩阵。

从（2）式中可得到

                                     f'(x) = f'(x0)+H(x0)(x-x0);                                                                                                      （3）

从多维函数的极值点定义可知，当f'(x) = 0，并且H（x）为正定矩阵的时候为极小值点

把f'(x)=0带入（3）式得

                                  0 = f'(x0)+H(x0)(x-x0)

                                   x = x0-1/(H(x0))*f'(x0)    (注意这里1/H（x0）是代表H(x0)的逆矩阵)                                          （4）

通过迭代（4） 可逼近极值点。

拟牛顿法：因为H（X）的逆矩阵不好求，个人认为可能H（x）矩阵并不是满秩矩阵的缘故，那么拟牛顿法就是要构造一个类似与H矩阵的东西。

根据（3） 有f'(x)-f'(x0) = H(x0)(x-x0)                                                                                                                                                                                                     （5）

（5）式被称问哦拟牛顿条件

当（4）式的H矩阵为正定矩阵时，我们令x = x0-STEP*1/H(x0)*f'(x0)；  (注意这里1/H（x0）是代表H(x0)的逆矩阵)      

那么根据泰勒展开式 f（x） = f(x0)-STEP*f'(x0)* 1/H（x0）*f'(x0)只要STEP足够小必定是下降的  。

那么我们要构造一个矩阵G，它是一个正定矩阵，并且满足拟牛顿条件，并且G矩阵每次迭代过程都是正定并且满足是正定的。

拟牛顿算法可以用 DFP算法，BFGS算法，DFP算法  。

下面贴一个梯度下降算法来最大化罗杰斯特回归的似然函数。

#include<math.h>
#include<stdio.h>
#include<conio.h>
#define STEP 0.01
int check(double w1[],double w2[],int len)
{
	double sum  = 0;
	for(int i = 0;i<len;i++)
	{
		sum+=(w1[i]-w2[i])*(w1[i]-w2[i]);
	}
	if(sum>0.000001) return 1;
	else return 0;
}
double dot(double x1[],double x2[],int len)
{
	double sum = 0;
	for(int i = 0;i<len;i++)
	{
		sum+=x1[i]*x2[i];
	}
	return sum;
}
double get_res(double y[],double x[4][3],double w[],int len)
{
	double L = 0;
    for(int i = 0;i<=3;i++)
	{
		L+=y[i]*dot(w,x[i],3)-log(1+exp(dot(w,x[i],3)));
	}
	return L;
}
void test1()
{
	double y[4];
	double x[4][3];
	double w[3];
	x[0][0] = 2,x[0][1] = 3,x[0][2] = 1,y[0] = 0;
	x[1][0] = 3,x[1][1] = 5,x[1][2] = 1,y[1] = 0;
	x[2][0] = 3,x[2][1] = 1,x[2][2] = 1,y[2] = 1;
	x[3][0] = 5,x[3][1] = 3,x[3][2] = 1,y[3] = 1;
	w[0] = w[1] = w[2] = 0;
	while(1)
	{
		printf("111111\n");
		double t[4];
		for(int i = 0;i<4;i++)
		{
			t[i] = 0;
		}
		for(int i = 0;i<4;i++)
		{
			double temp = dot(w,x[i],3);
			for(int j = 0;j<=2;j++)
			{
				t[j]+=y[i]*x[i][j]-exp(temp)*x[i][j]/(1+exp(temp));
			}
		}
		double pre_w[4];
		
		for(int i = 0;i<=2;i++)
		{
			printf("%lf ",w[i]);
			pre_w[i] = w[i]+STEP*t[i];
		} 
		printf("\n");
		for(int i = 0;i<=2;i++)
		{
			printf("%lf ",pre_w[i]);
		}
		printf("\n");
		printf("%lf %lf\n",get_res(y,x,w,3),get_res(y,x,pre_w,3));
		getch();
		if(!check(pre_w,w,3)) break;
		for(int i = 0;i<=3;i++)
		{
			w[i] = pre_w[i];
		}

	}


}