牛顿法(Newton's Method)

Newton's Method

在求最优解时，前面很多地方都用梯度下降(Gradient Descent)的方法，但由于最优步长很难确定，可能会出现总是在最优解附近徘徊的情况，致使最优解的搜索过程很缓慢。牛顿法(Newton's Method)在最优解的搜索方面有了较大改进，它不仅利用了目标函数的一阶导数，还利用了搜索点处的二阶导数，使得搜索算法能更准确地指向最优解。

我们结合下图所示的一个实例来描述牛顿法的思想。假设我们想要求得参数 heta，使得f( heta)=0。算法的描述如下：

1. 随机猜测一个解 heta^{(0)}，并令t=0；
2. 在 heta^{(t)}处用一根切线来近似f( heta);
3. 求得切线与横坐标的交点 heta^{(t+1)}，作为下一个可能的解;
4. t=t+1；
5. 重复2-4，直到收敛，即f( heta^{(t)})approx 0。

那么 heta^{(t+1)}与 heta^{(t)}之间存在怎样的迭代关系呢？由切线的斜率可知 egin{equation} f'( heta)=frac{f( heta)}{vartriangle}Rightarrow vartriangle=frac{f( heta)}{f'( heta)} end{equation}

观察 heta^{(t+1)}与 heta^{(t)}在横坐标上的关系，可知 egin{equation} heta^{(t+1)}= heta^{(t)}-vartriangle= heta^{(t)}-frac{f( heta)}{f'( heta)} end{equation}

牛顿法给出了f( heta)=0的求解算法，那么怎样将其运用到求使似然函数mathcal{L}( heta)最大化的参数上呢？一般最优参数 heta^{star}在mathcal{L}( heta)的极值点出取得，即mathcal{L}'( heta^{star})=0。那么，令上面的f( heta)=mathcal{L}'( heta)，我们很容易就得出了下列的迭代法则 egin{equation} heta^{(t+1)}= heta^{(t)}-frac{mathcal{L}'( heta^{(t)})}{mathcal{L}''( heta^{(t)})} end{equation}

 

 最终求得使mathcal{L}'( heta)=0的参数 heta^star，也就是令似然函数mathcal{L}( heta)最大的参数。

上面讨论的参数 hetainmathbb{R}，我们现在将牛顿法则推广到n维向量 hetainmathbb{R}^n，对应的迭代法则形式如下：egin{equation} heta^{(t+1)}= heta^{(t)}-H^{-1} abla_{ heta}mathcal{L} end{equation}

 

 其中H为mathcal{L}对向量 heta^{(t)}的二阶偏导，称为Hessian矩阵,H_{ij}=frac{partial^2mathcal{L}}{partial heta^{(t)}_ipartial heta^{(t)}_j}。

接下来，我们从另外一个角度来考察牛顿法。用似然函数mathcal{L}( heta)的二阶泰勒展开mathcal{F}( heta)来对其进行逼近。egin{equation} mathcal{L}( heta)approxmathcal{F}( heta)=mathcal{L}( heta^{(t})+ abla_{ heta^{(t)}}mathcal{L}( heta- heta^{(t)})+frac{1}{2}( heta- heta^{(t)})^TH( heta- heta^{(t)}) end{equation}

 

 令 heta= heta^{(t+1)}，可得 egin{equation} egin{array}{ll} mathcal{F}( heta^{(t+1)})=&mathcal{L}( heta^{(t)})+ abla_{ heta^{(t)}}mathcal{L}^T( heta^{(t+1)}- heta^{(t)})\ &+frac{1}{2}( heta^{(t+1)}- heta^{(t)})^TH( heta^{(t+1)}- heta^{(t)}) end{array} end{equation}

 

 现在，我们的目的是求得使mathcal{F}( heta^{(t+1)})最小的参数 heta^{(t+1)}。将上式对 heta^{(t+1)}求导并令导数为0，可得egin{equation} frac{partialmathcal{F}( heta^{(t+1)})}{partial heta^{(t+1)}}= abla_{ heta^{(t)}}mathcal{L}+H( heta^{(t+1)}- heta^{(t)})=0 end{equation}

 

 等式两侧同时左乘H^{-1}，化简得 egin{equation} heta^{(t+1)}= heta^{(t)}-H^{-1} abla_{ heta}mathcal{L} end{equation}

 

我们用的是二阶泰勒展开式mathcal{F}( heta)逼近似然函数mathcal{L}( heta)。如果mathcal{L}( heta)确实为二次函数，那么mathcal{F}( heta)就是mathcal{L}( heta)的准确展开式，利用牛顿法一步就可以直接求得最优解。一般情况下，mathcal{L}( heta)并非二次函数，那么mathcal{F}( heta)也就存在逼近误差，使得一次迭代不能求得最优解，当mathcal{L}( heta)的次数很高时，往往要经历很多次迭代。一般而言，因为牛顿法利用了二阶导数来修正搜索方向和步长，收敛速度很更快。但是这同样也是要付出代价的，相比梯度下降而言，我们需要额外计算Hessian矩阵并求其逆，这两步的计算代价都很大。只要参数 heta的维度n不是很大，可以考虑用牛顿迭代。另外还有一点，如果目标函数不是严格的凸函数，Hessian矩阵H很可能是奇异矩阵，也就是存在特征值为0的情况，那么它的逆矩阵是不存在的，也就无法用牛顿法。

今年有一道面试题是要求我们写出一段程序，求解sqrt{n}。如果把牛顿法用上去，问题就迎刃而解了。我们设定目标函数为f(x)=x^2-n，那么令f(x)=0的解很显然就是pmsqrt{n}。要注意的是，我们要选择合理的迭代起始点，如果我们从正数开始迭代，求得的是sqrt{n};如果从负数开始迭代，求得的就是-sqrt{n}；如果从0开始迭代，会出现未定义的计算(0作为除数)。我们根据前面讲的牛顿迭代法则，直接给出该题的迭代法则 egin{equation} x^{(t+1)}=x^{(t)}-frac{f(x^{(t)})}{f'(x^{(t)})}=x^{(t)}-frac{(x^{(t)})^2-n}{2x^{(t)}}=frac{1}{2}left(x^{(t)}+frac{n}{x^{(t)}} ight) end{equation}

 

 下面是由该算法写出的一段精简的code，浓缩了牛顿算法的精髓

 1 double mysqrt1(double n)
 2 {
 3     if (n<0) return -1;
 4     if(n==0) return 0;
 5     double eps=1e-5;
 6     double x=0.1;//start from a positive value
 7     while(fabs(x*x-n)>=eps)
 8         x=(x+n/x)/2;//Newton's method
 9     return x;
10 }

这道题我还想了另外一个算法，算法的启发点来源于(x-1)(x+1)+1=x^2=n。用这个算法，我们的迭代起始点可以是0。算法的基本思想如下：给定一个初始步长step，从起始点开始每次向前走一个步长，直到超过了sqrt{n}；一旦超过了sqrt{n}，就要开始慢慢向最终解靠近，每次前进或后退的步长都缩减为以前的一半。很明显，这个算法没有牛顿迭代法快。我只用了少数几个测试用例，两段程序的计算结果都和sqrt库函数的计算结果一致。代码如下：

 1 //Inspiration comes from the eqution:(x+1)(x-1)+1=x*x
 2 //We can search the solution from an initial point x=0 with a step size
 3 //and use the equation above to estimate x*x until we approach it close enough
 4 double mysqrt2(double n)
 5 {
 6     if (n<0) return -1;
 7     double x=0;
 8     double step=10;
 9     int threshold=0;
10     double eps=1e-5;
11     double res=(x+1)*(x-1)+1;
12     while(fabs(res-x)>=eps)
13     {
14         if(res<x)
15         {
16             //once we have passed the solution,we must walk forward slowly 
17             if(threshold) step/=2;
18             x+=step;//walking forward
19         }
20         else//walk
21         {
22             threshold=1;//indicating we have passed the solution
23             step/=2;//reducing the step size to its half
24             x-=step;//walking back
25         }
26         res=(x+1)*(x-1)+1;//compute x*x to estimate real n
27     }//end while
28     return x;
29 }

 

 

分类: Machine Learning, Programming

标签: Newton, 牛顿迭代

随笔 - 24, 文章 - 8, 评论 - 12, 引用 - 0

[算法][庞果网]倒水问题/量水问题

题目详情

有两个容器，容积分别为A升和B升，有无限多的水，现在需要C升水。
我们还有一个足够大的水缸，足够容纳C升水。起初它是空的，我们只能往水缸里倒入水，而不能倒出。
可以进行的操作是：
把一个容器灌满；
把一个容器清空（容器里剩余的水全部倒掉，或者倒入水缸）；
用一个容器的水倒入另外一个容器，直到倒出水的容器空或者倒入水的容器满。
    问是否能够通过有限次操作，使得水缸最后恰好有C升水。

输入：三个整数A, B, C，其中 0 < A , B, C <= 1000000000
输出：0或1，表示能否达到要求。

函数头部：
c语言：1表示可以，0表示不可以
int can(int a,int b,int c);
c++语言: true表示可以，false表示不可以
bool can(int a,int b,int c);
java语言:true表示可以，false表示不可以
public class Main {
    public static boolean can(int a,int b,int c);
}

解题思路

这是一个典型的倒水问题/量水问题，使用欧几里得算法就可解出来。

 

这里有一篇文章给出了简单的倒水问题的解法，可以解决笔试面试里面一些简单的填空题，可以看一看：http://blog.csdn.net/morewindows/article/details/7481851

基本思想是：不断用小桶装水倒入大桶，大桶满了立即清空，每次判断下二个桶中水的容量是否等于指定容量。也就是用小桶容量的倍数对大桶的容量进行取余，直到余数等于指定容量。

例如，用7升的桶和11升的桶得到2升水可以这样做：
7 % 11 = 7
14 % 11 = 3
21 % 11 = 10
28 % 11 = 6
35 % 11 = 2  成功得到2升水。

对于明确说明可以得到xx升水，需要我们给出如何倒出来的步骤，可以用这个方法，很快捷。但是这个方法不适合解这道题。

 

欧几里德算法又称辗转相除法，用于计算两个正整数a，b的最大公约数。

用gcd(a,b) 表示a, b的最大公约数，则有定理：gcd(a,b) = gcd(b,a mod b) (a>b 且a mod b 不为0)

具体的算法实现有循环和递归两种，我用的是循环的方法。

 

扩展欧几里得算法
定理：对于不完全为 0 的非负整数 a,b，gcd(a, b)表示 a, b 的最大公约数，必然存在整数对 x, y ，使得 gcd(a,b)=ax+by。

 

本题实际上是问是否存在整数x, y，使得ax+by=c成立。如果c可以被gcd(a,b)整除，则成立。

 

因此解题步骤如下：

1. 求出gcd(a,b)
2. 判断c是否能被gcd(a,b)整除，若能则返回true，否则返回false

Java代码

 1 public static boolean can(int a,int b,int c) {
 2     int r;
 3     while (true) {
 4         r = a % b;
 5         if (r != 0) {
 6             a = b;
 7             b = r;
 8         } else {
 9             break;
10         }
11     } 
12     return c%b == 0;    //b现在是最大公约数，若c能被b整除，则可以
13 }

 

 

 

标签: 算法