拉格朗日插值原理及实现（Python）
目录一. 前言
Lagrange插值是利用n次多项式来拟合（n+1）个数据点从而得到插值函数的方法。（注意n次多项式的定义是未知数最高次幂为n，但是多项式系数有n+1个，因为还有个常数项）
Lagrange插值和Newton插值本质上相同，都是用(n-1)次多项式来拟合n个数据点。所以这两种插值方法得到的插值函数相同，因为多项式拟合的基本定理：同时通过n个数据点，且最高次幂小于(n-1)的多项式函数唯一。下面顺手证明一下这个重要的定理。
如果已知n+1个数据点\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)，假设\(L_1=k_0+k_1x+k_2x^2+\cdots+k_nx^n\)和\(L_2=k_0’+k_1’x+k_2'x^2+\cdots+k_n’x^n\)都是通过这n个数据点的插值函数。那么应该有\(L_1 - L_2\)通过所有\((x_1,0),(x_2,0),\cdots,(x_n,0)\)
代入这些点得到齐次线性方程组：
\[\begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{bmatrix} \begin{bmatrix} k_0-k_0'\\ k_1-k_1'\\ k_2-k_2'\\ \vdots \\ k_n-k_n'\\ \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \]它的系数行列式是Vandermonde行列式，所以：
\[\begin{vmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n\\ 1 & x_1 & x_1^2 & \cdots & x_1^n\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n\\ \end{vmatrix} = \prod_{n\ge i > j \ge 0 }(x_i-x_j) \]因为每个点都是不同的，所以\(x_i \ne x_j,i\ne j\)，所以齐次线性方程组的系数行列式不等于0，故方程解唯一且为0解：
\[k_i - k_i' = 0\\ k_i = k_i' \\ L_1 = L_2 \]证毕。
二. 3种形式的Lagrange插值函数推导
1. 原始形态的Lagrange插值
为了用n次多项式拟合n+1个数据点：\((x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)\)
Lagrange插值函数采用的方法是构建一个这样的函数：
\[L(x) = l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \tag{1} \]也就是用一组基函数\(\{l_0(x),l_1{x}, \cdots ,l_n(x)\}\)去构建插值函数\(L(x)\)，那么不难想到这样的基函数需要满足这样的条件：
\[l_i(x_j)=\begin{cases} 0 , & j\ne i \\ 1 , & j = i \end{cases} \]这样对于\(L(x)\)就会有：
\[\begin{split} L(x_i) &= l_0(x_i)y_0 + \cdots l_i(x_i)y_i + \cdots +l_n(x)y_n \\ &=y_i \end{split} \]这样就实现了\(L(x)\)通过所有的数据点\((x_i,y_i)\)。接下来就构建这样的基函数\(l_i(x)\)。
首先实现\(l_i(x_j) = 0,i\ne j\) ：
\[l_i(x) = (x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1}) \cdots(x-x_n) \]这个函数实现了\(l_i(x)\)在所有非\((x_i,y_i)\)的点处为0，但是在\(x=x_i\)处：
\[l_i(x_i) = (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \]为了让\(l_i(x_i)=1\)，我们可以将\(l_i(x)\)除以这个值进行修正：
\[\begin{split} l_i(x_i)&= \frac{(x-x_0)(x -x_1)\cdots (x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}{(x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n)}\\ &=\prod_{n\ge j \ge0,j\ne i} \frac{x-x_j}{x_i -x_j} \end{split} \tag{2} \]将(2)代入(1)就得到了原始形态的Lagrange插值函数.
\[\begin{split} L(x) &= l_0(x)y_0 + l_1(x)y_1 + \cdots +l_n(x)y_n \\ &= \prod_{n\ge j \ge0,j\ne 0} \frac{x-x_j}{x_0 -x_j}y_0 +\prod_{n\ge j \ge0,j\ne 1} \frac{x-x_j}{x_1 -x_j}y_1 +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{x-x_j}{x_n -x_j}y_n \end{split} \tag{3} \]例：已知点(1,1),(2,2),(3,3)，用原始Lagrange插值计算插值函数。
首先计算三个插值基函数：
\[l_0(x) = \frac{(x-2)(x-3)}{(1-2)(1-3)}=\frac{(x-2)(x-3)}{2}\\ l_1(x) = \frac{(x-1)(x-3)}{(2-1)(2-3)}=-{(x-1)(x-3)}\\ l_2(x) = \frac{(x-1)(x-2)}{(3-1)(3-2)}=\frac{(x-1)(x-2)}{2}\\ \]从而得到插值函数：
\[L(X) = l_0(x)+ 2l_1(x)+3l_2(x)=\frac{(x-2)(x-3)}{2}-2{(x-1)(x-3)}+ \frac{3(x-1)(x-2)}{2} \]插入点x=4试一下：\(L(4)=4\) 。
原始模式的Lagrange插值函数，每次计算插值点时需要进行n(n-1)次乘法，时间复杂度为\(O(n^2)\)
2. 第一形式Lagrange插值
为了降低计算的时间复杂度，我们对原始的Lagrange插值函数进行改进。
为了书写方便，我们令\(w_i\)为
\[\begin{split} w_i &= \prod_{n\ge j \ge0,j\ne i} (x_i -x_j) \\ &= (x_i-x_0)(x_i -x_1)\cdots (x_i-x_{i-1})(x_i-x_{i+1})\cdots (x_i-x_n) \end{split} \tag{4} \]观察(3)式，我们可以提取一个公因数\((x-x_0)(x-x_1)\cdots (x-x_n) = g(x)\)：
\[\begin{split} L(x) &= (x-x_0)(x-x_1)\cdots (x-x_n)[\prod_{n\ge j \ge0,j\ne 0} \frac{y_0}{(x_0 -x_j)(x-x_0)} +\prod_{n\ge j \ge0,j\ne 1} \frac{y_1}{(x_1 -x_j)(x-x_1)} +\cdots +\prod_{n\ge j \ge0,j\ne n} \frac{y_n}{(x_n -x_j)(x-x_n)}]\\ &=g(x)[\frac{y_0}{w_0(x-x_0)}+\frac{y_1}{w_1(x-x_1)}+\cdots+\frac{y_n}{w_n(x-x_n)}]\\ &= g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \end{split} \tag{5} \]这样就得到了第一形式Lagrange插值。这一形式的Lagrange插值计算的流程如下：
- 预处理：根据已知数据点，利用公式(4)计算\(w_0,w_1,\cdots,w_n\);
- 插值：计算\(g(x)=(x-x_0)(x-x_1)\cdots(x-x_n)\)，然后根据公式(5)，先计算中括号内的加式，然后乘以\(g(x)\)得到插值点的值，这样计算一个新的点的时间复杂度就是\(O(n)\)；
- 补充数据点：如果数据集有更新，只需要更新\(w_i\)即可，加入一个新的点到数据集，只需要将每个\(w_i\)乘以\((x_i-x_{n+1})\)，此外再增加一个\(w_{n+1}\)。
例：已知点(1,1),(2,2),(3,3)，使用Lagrange第一形式计算插值函数。
首先计算\(w_i\)：
\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]插值函数就是：
\[L(x)=(x-1)(x-2)(x-3)[\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}] \]插入点x=4试一下：\(L(4)=4\) 。
3. 第二形式的Lagrange插值（重心插值公式）
第一形式的Lagrange插值还要计算\(g(x)\)，可以再进行优化。
第一形式：
\[L(x) = g(x)\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)} \tag{6} \]为了消掉\(g(x)\)，我们取\(y=1\)这条直线上的点进行插值，即取n+1个点：\((x_0,1),(x_1,1),\cdots,(x_n,1)\)那么这n+1个点的插值函数就是：
\[L'(x)=g(x)\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)} \tag{7} \]而\(L'(x)=1,x=x_0,x_1,\cdots,x_n\) ，所以我们可以用(6)除以(7)消去g(x)：
\[L(x)=\frac{L(x)}{1}=\frac{L(x)}{L'(x)}=\frac{\sum_{i=0}^{n}\frac{y_i}{w_i(x-x_i)}}{\sum_{i=0}^{n}\frac{1}{w_i(x-x_i)}} \tag{8} \]这样就得到了第二形式的Lagrange插值，也称为重心插值，通常Lagrange插值采用这种形式。
它的计算过程如下：
- 预处理：计算\(w_i\)；
- 插值：对需要插入的点(x,y)，计算（x-x_i)（可以存到一个列表中，计算时直接取用）；
- 补充数据点：补充新的点到数据集只需要更新\(w_i\) 。
例：已知点(1,1),(2,2),(3,3)，利用重心插值公式计算插值函数。
首先计算\(w_i\)：
\[w_0=(1-2)(1-3)=2\\ w_1=(2-1)(2-3)=-1\\ w_2=(3-1)(3-2)=2 \]得到\(L(x)\)：
\[L(x)=\frac{\frac{1}{2(x-1)}-\frac{2}{(x-2)}+\frac{3}{2(x-3)}}{\frac{1}{2(x-1)}-\frac{1}{(x-2)}+\frac{1}{2(x-3)}} \]插入点x=4试一下：\(L(4)=4\) 。
三. 利用Python编程实现这三种Lagrange插值
import numpy as np
class LagrangeInterpolation:
"""
There are three modes of the LagrangeInterpolation.The input data is supposed
to be two lists.
---------------------------------------------------------------------------------------
self.data : Contain the points known before.
self.dataLength : Indicate the number of points in the data list.
self.weight : The self.weight[i] is (xi - x1)(xi - x2)...(xi - x(i-1))(xi - x(i+1))
...(xi - xn)
self.items : The self.items[i] is (x - x1)(x - x2)...(x - x(i-1))(x - x(i+1))...
(x - xn)
---------------------------------------------------------------------------------------
"""
def __init__(self, x, y):
self.data = {'x': list(x), 'y': list(y)}
self.dataLength = len(self.data['x'])
self.weight = []
self.items = []
# control is a flag indicating if there is anything wrong with the data.
self.control = True
if len(self.data['x']) != len(self.data['y']):
print("The length of x isn't equal to the length of y!")
self.control = False
else:
self.__preprocess(order=0)
# Appending function is used to add points to the data list.
def data_append(self, ap):
if ap[0] in self.data['x']:
print("The point already exist.")
self.control = False
else:
self.control = True
self.data['x'].append(ap[0])
self.data['y'].append(ap[1])
self.dataLength = self.dataLength + 1
self.__preprocess(order=1)
"""
Preprocessing is used to update the self.weight and self.items.
order = 0 : Initialize the self.weight.
order = 1 : Update the self.weight when new point is added to the data list.
order = 2 : Calculate the self.items for each point waiting for interpolation.
"""
def __preprocess(self, order=0, x=0):
if order == 0:
self.weight = list(np.zeros(self.dataLength))
for i in range(self.dataLength):
weight_temp = 1
for j in range(self.dataLength):
if i == j:
pass
else:
weight_temp = weight_temp * (self.data['x'][i] - self.data['x'][j])
self.weight[i] = weight_temp
elif order == 1:
self.weight.append(1)
for i in range(self.dataLength - 1):
self.weight[i] = self.weight[i] * (self.data['x'][i] - self.data['x'][-1])
self.weight[-1] = self.weight[-1] * (self.data['x'][-1] - self.data['x'][i])
elif order == 2:
self.items = list(np.zeros(self.dataLength))
for i in range(self.dataLength):
self.items[i] = x - self.data['x'][i]
# The mode1 is the initial mode of Lagrange interpolation.
def mode1(self, px):
self.__preprocess(order=2, x=px)
if self.control:
dataCheck = False
for w in self.weight:
if w != 0:
dataCheck = True
else:
dataCheck = False
if dataCheck:
py = 0.0
for i in range(self.dataLength):
py_temp = 1
for j in range(self.dataLength):
if i == j:
pass
else:
py_temp = py_temp * self.items[j]
py = py + py_temp * self.data['y'][i] / self.weight[i]
return py
else:
print("There is a same x!")
return None
else:
return None
def mode2(self, px):
self.__preprocess(order=2, x=px)
itemsProd = np.prod(self.items)
itemsSum = 0
for i in range(self.dataLength):
itemsSum = itemsSum + self.data['y'][i] / (self.weight[i] * self.items[i])
py = itemsProd * itemsSum
return py
def mode3(self, px):
self.__preprocess(order=2, x=px)
denomSum = 0
numeSum = 0
for i in range(self.dataLength):
dtemp = self.weight[i] * self.items[i]
numeSum = numeSum + self.data['y'][i] / dtemp
denomSum = denomSum + 1 / dtemp
py = numeSum / denomSum
return py
demo.py :
from lagrange_interpolation import LagrangeInterpolation as lag
x = [1, 2, 3, 4]
y = [1, 2, 3, 4]
inter1 = lag(x, y)
inter1.data_append((5, 5)) #往数据集中追加一个点
z1 = inter1.mode1(6)
z2 = inter1.mode2(6)
z3 = inter1.mode3(6)
print(z1)
print(z2)
print(z3)
"""
Result：
6.0
6.000000000000002
6.0000000000000036
"""