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

梯度下降法：根据泰勒展开式： 
                                            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")
a020
;
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];
}

}

}