拟牛顿法面面俱到(一)--牛顿插值法

石晓文

发布于 2018-07-25 14:46:35

1.8K1

发布于 2018-07-25 14:46:35

文章被收录于专栏：小小挖掘机小小挖掘机

这次带来的是拟牛顿法系列，本系列的目标是完全理解拟牛顿法，包括其中涉及到的知识，比如泰勒公式、海森矩阵等，泰勒公式大家都很熟悉，不过它是怎么推导出来的呢？想必大家都不是很了解吧，这要从牛顿插值法说起，本节就先来讲解一下牛顿插值法。

本文大都参考自知乎：https://www.zhihu.com/question/22320408

1、什么是插值

什么叫插值？插值是数学领域数值分析中的通过已知的离散数据求未知数据的过程或方法。

相信大家对插值都不陌生，我们在数据挖掘进行缺失值处理时，有时是直接将数据丢弃，有时是用插值法填充一个数进去。

想想我们是怎么做的？有时候直接拿上一个数进行填充，有时候直接拿后一个数进行填充，还有的时候，我们用前后两个数的平均值进行填充。使用最后一种方法时，我们其实就用到了一种简单的差值思想--线性插值法。

除线性插值法外，还有许多常用的插值方法，我们将在下一节介绍。

2、常见的插值方法

2.1 线性插值法

这是最简单的插值方法，示意图如下：

这种近似太粗糙，我们只需要知道前后两个点的数据就可以进行插值，但实际的过程往往没有这么简单。

2.2 多项式插值

牛顿插值法也算是多项式插值中的一种，但我们将牛顿插值法单独拿出一节进行讲解。这里介绍另一种多项式插值方法，过程如下：

这样求解出的三次多项式(如果有唯一解的话)，一定同时经过已知的四个点。

不过这样进行求解有两个弊端：

计算量大，当数据量成千上万时，我们需要求解的参数也是成千上万的，效率十分低下。
新增加一个观测数据，我们需要重新进行计算

为了解决上面的两个问题，我们有了牛顿插值法。

3、牛顿插值法

牛顿插值法全名是格雷戈里-牛顿公式，格雷戈里和牛顿分别给出了这个插值公式，主要牛顿太耀眼了，所以格雷戈里都被大家遗忘了。

牛顿插值法的特点在于：每增加一个点，不会导致之前的重新计算，只需要算和新增点有关的就可以了。

下面就进入数学阶段了，前方高能预警，非战斗人员请退避，直接跳过推导阶段。

3.1 牛顿插值法的推导

我们先把问题数学化：

下面两张图讲解了牛顿插值法的大体过程：

观察b1，b2的特点，不断重复上面的过程，我们就可以得到牛顿插值法的计算公式。

然后是详细的推导过程(图片来自上面提到的知乎)：

上面有一句话说的没错，我也觉得b2的推导错了，我们按照正常的逻辑来一遍：

哈哈，是不是有很多小伙伴跟我得到了一样的结果，那上面的结果是怎么来的呢？可以看到，我们是先确定的x0，然后推导了满足x1的插值公式，那么加入这个思路反过来呢？我们先确定的是x1，然后推导得到满足x0的公式，那么结果如下：

数学真是一个奇妙的玩意，需要静下心来细细品味呀。

我们提炼一下刚才的过程：

当新增一个点时，我们只需计算新的均差即可：

4、Python代码实现

下面的例子是对牛顿插值法的一个简单实现：

import numpy as np
import matplotlib.pyplot as plt


# 递归求差商
def get_diff_quo(xi, fi):
    if len(xi) > 2 and len(fi) > 2:
        return (get_diff_quo(xi[:len(xi) - 1], fi[:len(fi) - 1]) - get_diff_quo(xi[1:len(xi)], fi[1:len(fi)])) / float(
            xi[0] - xi[-1])
    return (fi[0] - fi[1]) / float(xi[0] - xi[1])


# 求w，使用闭包函数
def get_w(i, xi):
    def wi(x):
        result = 1.0
        for j in range(i):
            result *= (x - xi[j])
        return result

    return wi


# 做插值
def get_Newton(xi, fi):
    def Newton(x):
        result = fi[0]
        for i in range(2, len(xi)):
            result += (get_diff_quo(xi[:i], fi[:i]) * get_w(i - 1, xi)(x))
        return result

    return Newton


# 已知结点
xn = [i for i in range(-50, 50, 10)]
fn = [i ** 2 for i in xn]

# 插值函数
Nx = get_Newton(xn, fn)

# 测试用例
tmp_x = [i for i in range(-50, 51)]
tmp_y = [Nx(i) for i in tmp_x]

print(tmp_x)
print(tmp_y)

# 作图
plt.plot(xn, fn, 'r*')
plt.plot(tmp_x, tmp_y, '')
plt.title('Newton Interpolation')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

具体的代码详情我们就不介绍了，最后来看看实验效果吧：

5、后记

牛顿插值法可以用来做什么呢？这不是废话么，当然是插值啦，哈哈。其实，我们不知道的是，它还可以用来推导泰勒公式呢！欲知后事如何，且听下回分解！