吹水牛顿迭代法

因为吹水的能力不佳,所以要先打个草稿,今天的吹水过程大概是:
1、牛顿迭代法的演绎过程
2、牛顿迭代法求n次方根
3、牛顿迭代法求n次方根改进版
4、牛逼哄哄的invsqrt求平方根倒数

1、牛顿迭代法的演绎过程

乍一听,好像很高大上,其实理解上没什么难度。
牛顿迭代法就是在不断迭代的过程中逼近曲线的根,可以看做,它就是用来求曲线的根的。
所以它是怎么个迭代的过程,先上个动图:

牛顿迭代法的过程

我用吹水的姿势描述一遍动图的过程,先是有一个曲线,然后要找到曲线和x轴的交点,这个时候,脑袋一拍,选x1作为第一点,在这个点上垂直于x轴的线与曲线的交点,在这个交点做与曲线的切线,而切线与x轴的交点就是下一个迭代的点,不断重复这个过程,直到找到曲线和x轴交点这个目标值。值得注意的是,这个方法是不断逼近,所以得到的值可能会与目标值存在误差,而误差小于一定范围就可以认为是目标值了。一个更加图文并茂更加通俗易懂的解释可以点击 这个这个这个这个

2、牛顿迭代法求n次方根

那为什么说牛顿迭代法能用来求n次根呢,其实在1中的目标值,即曲线与x轴相交点的x值,假设你的方程为 f(x) = x² - 64,当f(x) = 0时,x就是64的2次方根,以此类推,把2换成你想要求的n次方,就可以求出64的n次方根了。

明白了演绎过程,那我们就来探索它的代数实现。在上代码前,先上个图


牛顿迭代法某一点的求值图解

每次迭代要求的点就是右绿色箭头指向的这个点,在图中可以看到黑色箭头标识的线段长度为-f(x),f'(x)是对f(x)的求导结果,也就是切线的斜率,所以绿色箭头线段的长度为-f(x)/f'(x),而左绿色箭头的点为x,所有右绿色箭头指向的值为x加上绿色线段的长度,即x-(f(x)/f'(x))。

上面的代数实现完成后,就可以上代码了。嘻嘻嘻,最近在看Python版的SICP,所以贴上Python的代码,这篇文章也是因为看了这本书的1.6章节写的。def就是定义一个方法的意思。

def newton_nth_root_of_a(n, a, x = 1, e = 1e-13):
    def f(x):
        return x ** n - a
    def df(x):
        return n * x ** (n - 1)
    
    gap = abs(f(x))
    while gap > e:
        x = x - f(x)/df(x)
        gap = abs(f(x))
    return x
    
print(newton_nth_root_of_a(1, 64))  #64.0
print(newton_nth_root_of_a(2, 64))  #8.0
print(newton_nth_root_of_a(3, 64))  #4.0
print(newton_nth_root_of_a(4, 64))  #2.82842712474619
print(newton_nth_root_of_a(5, 64))  #2.29739670999407
print(newton_nth_root_of_a(6, 64))  #2.0

这里是在线可运行版本 ,不多不少,刚好十行,newton_nth_root_of_a方法定义中,n代表幂,a代表对a求根,x代表起始值,e代表误差范围,x和e已经设置了缺省值。f(x)是曲线,df(x)是f(x)的求导结果,在while循环中,x - f(x)/df(x)就是我们上面代数实现中的结果,直到f(x)的值小于误差,x就是最后的结果。

3、牛顿迭代法求n次方根改进版

明白了1、2其实整个过程就是这样了,那为啥有3,因为要硬着头皮接着吹水😂。
2中的10行代码很简洁,但是假如我要不断的求某些值的2次根,就要不断调用newton_nth_root_of_a(2, a)。那么我们来改进下代码,保持上面的代码不变,增加下面的代码:

def curry_nth_root(n):
    def curry_nth_root_of_a(a):
        return newton_nth_root_of_a(n, a)    
    return curry_nth_root_of_a

_2th_root = curry_nth_root(2)

print(_2th_root(64)) #8.0
print(_2th_root(49)) #7.000000000000001
print(_2th_root(36)) #6.0
print(_2th_root(25)) #5.0
print(_2th_root(16)) #4.0
print(_2th_root(9))  #3.0

这里是在线可运行版本 ,哇咔咔,在curry化之后,就可以方便的求n次方根了,不用每次都输入n。

curry化其实也很简单,那么我们就来个极端点的抽象代码:

def improve_guess(update, close, guess=1):
    '''返回猜想值,update为更新方法,close为逼近方法,guess为初始猜想值'''
    while not close(guess):
        guess = update(guess)
    return guess
    
def newton_update(f, df):
    '''返回牛顿迭代的计算函数, f为原函数,df是对其求导'''
    def update(x):
        return x - f(x) / df(x)
    return update
    
def newton_find_zero_of_f(f, df, tolerance=1e-13):
    '''找出f(x)=0时,x的值,tolerance为与目标真实值的误差'''
    def near_zero(x):
        return is_eq(f(x), 0, tolerance)
    return improve_guess(newton_update(f, df), near_zero)
    
def is_eq(x, y, tolerance):
    '''x与y的差值是否在误差范围内'''
    return abs(x - y) < tolerance
    
def nth_root_of_a(n, a):
    def f(x):
        return x ** n - a
    def df(x):
        return n * x ** (n - 1)
    return newton_find_zero_of_f(f, df)
    
print(nth_root_of_a(1, 64))  #64.0
print(nth_root_of_a(2, 64))  #8.0
print(nth_root_of_a(3, 64))  #4.0
print(nth_root_of_a(4, 64))  #2.82842712474619
print(nth_root_of_a(5, 64))  #2.29739670999407
print(nth_root_of_a(6, 64))  #2.0

这里是在线可运行版本 ,把每一步都抽象出来,哈哈哈哈哈,这看上去很帅,虽然实际上没啥必要,最后也可以再curry化。这是个极端的例子,但是根据业务的需要,我们可以把需要复用的步骤抽象出来,而抽象的程度也要根据业务开展。

4、牛逼哄哄的invsqrt求平方根倒数

吹水时间又到了,究竟这个牛顿迭代法有个毛线用。那你就要google下invsqrt这个改变游戏史进程的函数,可以看下知乎上的这个回答 有哪些算法惊艳到了你?。这个函数到底是干嘛的,就是求某个数的2次方根的倒数。因为这个函数在3d引擎的中经常使用,所以提高这个函数的效率是至关重要的,这个函数比系统的1/sqrt(x)快了一半。虽然是倒数,但是也是能用牛顿迭代法的啊。为毛啊,因为这也是曲线啊,过程跟1中的图相识。

但是!!!!请看代码

float InvSqrt (float x){
    float xhalf = 0.5f*x;
    int i = *(int*)&x;
    i = 0x5f3759df - (i>>1);
    x = *(float*)&i;
    x = x*(1.5f - xhalf*x*x);
    return x;
}

第六行就是利用牛顿迭代法的值x - f(x)/f'(x)演算后的结果,在这里f(x) = 1/x²,则f'(x)=2/x³。
第三行是把一个float指针强转成int指针,从而进行右移运算,因为float类型不能进行位运算,从而得到x值指数部分/2的结果,额,这里我不打算展开说啊,因为我说不清楚啊啊啊啊,我引用一段别人的解释,详情请移步这里

所以,32位的浮点数用十进制实数表示就是:M2E。开根然后倒数就是:M(-1/2)2^(-E/2)。现在就十分清晰了。语句i> >1其工作就是将指数除以2,实现2(E/2)的部分。而前面用一个常数减去它,目的就是得到M(1/2)同时反转所有指数的符号。

至于0x5f3759df,这就相当无解了。wiki的介绍中,这一句的注释是这样的:

    i  = 0x5f3759df - ( i >> 1 );               // what the fuck?(这他妈的是怎么回事?)

关于这个值后来也有人专门去推理过,还写成论文,但是,我肯定是看不懂的!!!
因为这一步,取到一个非常合适的初始值,只需要一次牛顿迭代就能知道目标值啦,一次!!!!

最后

没有最后,今天吹水结束!!!!!!!!

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,222评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,455评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,720评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,568评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,696评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,879评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,028评论 3 409
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,773评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,220评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,550评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,697评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,360评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,002评论 3 315
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,782评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,010评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,433评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,587评论 2 350

推荐阅读更多精彩内容