拉格朗日插值基函数导出矩阵的快速计算

1. 问题

给定一组插值节点的横坐标$\{x_i\}_{i=1}^{n}$, 和赋值节点$\{y_i\}_{i=1}^{\red{m}}$, 需要求出定义在节点$\{x_i\}_{i=1}^{n}$上的$n-1$次拉格朗日基函数$h_{i}(x)$在节点$\{y_i\}_{i=1}^{\red{m}}$处的值, 以及拉格朗日基函数的导数$L'_{i}(x)$在节点$\{y_i\}_{i=1}^{\red{m}}$处的值. 也就是要求出矩阵$\mathbb{H}=[h_{j}(y_i)]_{i=1,\red{m}}^{j=1,n}$以及矩阵$\mathbb{D}=[h'_{j}(y_i)]_{i=1,\red{m}}^{j=1,n}$

1.1. 基本性质

记 $\widetilde{\mathbb{D}}=[h'_{j}(x_i)]_{i=1,n}^{j=1,n}$, $\widetilde{\mathbb{S}}=[h''_{j}(x_i)]_{i=1,n}^{j=1,n}$ 以及 ${\mathbb{S}}=[h''_{j}(y_i)]_{i=1,n}^{j=1,n}$

• $h'_j(y_i)=\sum_{k=1}^n h'_j(x_k) h_{k}(y_i)$, 或 $\mathbb{D}=\mathbb{H}\cdot\mathbb{\widetilde{D}}$

• $h''_j(x_i)=\sum_{k=1}^n h'_j(x_k) h'_{k}(x_i)$, 或 ${\mathbb{\widetilde S}}=\mathbb{\widetilde{D}}^2$

• $h''_j(y_i)=\sum_{k=1}^n h'_j(x_k) h'_{k}(y_i)$, 或 $\mathbb{S}=\mathbb{D}\cdot \mathbb{\widetilde{D}}=\mathbb{H}\cdot\mathbb{\widetilde{D}}^2=\mathbb{H}\cdot \mathbb{\widetilde S}$

• ...

• $h^{(n)}_j(x_i)=\sum_{k=1}^n h^{(n-1)}_j(x_k) h'_{k}(x_i)$, 或 $\mathbb{\widetilde A}_n=\mathbb{\widetilde D}\cdot {\mathbb{\widetilde A}}_{n-1}=\mathbb{\widetilde D}^n$

• $h^{(n)}_j(y_i)=\sum_{k=1}^n h^{(n-1)}_j(x_k) h'_{k}(y_i)$, 或 $\mathbb{A}_n=\mathbb{ D}\cdot \widetilde{\mathbb{A}}_{n-1}=\mathbb{H}\cdot \mathbb{\widetilde D}^n$

证明: 根据多项式插值原理, 如果$f(x)$是$[-1,1]$上的$n-1$多项式, 那么有

$$f(x)=\sum_{k=1}^n f(x_k) h_{k}(x), \quad f'(x)=\sum_{k=1}^n f(x_k) h'_{k}(x),$$

• 上面左式取$f=h'_j$ 以及 $x=y_i$, 则 $h'_j(y_i)=\sum_{k=1}^n h'_j(x_k) h_{k}(y_i)$.

• 上面右式取$f=h'_j$ 以及 $x=x_i$, 则 $h''_j(x_i)=\sum_{k=1}^n h'_j(x_k) h'_{k}(x_i)$.

• 上面右式取$f=h'_j$ 以及 $x=y_i$, 则 $h''_j(x_i)=\sum_{k=1}^n h'_j(x_k) h'_{k}(y_i)$.

• ...

根据该基本性质, 实际计算, 只需计算出矩阵$\mathbb{H}$与$\mathbb{\widetilde{D}}$, 就可得到其它矩阵

1.2. 求解

1.2.1. 计算矩阵H

定义矩阵$\mathbb{C}=[c_{jk}]_{j=1,n}^{k=1,n}$, 以及三维矩阵$\mathbb{Z}=[z_{ijk}]_{i=1:\red{m}}^{j,k=1:n}$, 其中

$$c_{jk}= \begin{cases} x_j-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases},\quad \text{以及 } z_{ijk} = \begin{cases} y_i-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases} .$$

根据公式

$h_j(y_i)=\dfrac{(y_i-x_1)(y_i-x_2)\cdots(y_i-x_{j-1})(y_i-x_{j+1})\cdots(y_i-x_{n})}{(x_j-x_1)(x_j-x_2)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{n})}$

可得

$h_j(y_i)=\dfrac{z_{ij1}z_{ij2}z_{ij3}\cdots z_{ijn}}{c_{j1}c_{j2}c_{j3}\cdots c_{jn}}$

Python的三维矩阵除以二维矩阵, 运算法则是这样的 $A_{ijk}/B_{jk}$, 因此可以有很简洁的表达式

$\mathbb{H}=\left(\dfrac{\mathbb{Z}}{\mathbb{C}}\right)_{最后一个轴累乘}$.

全部Python代码如下

import numpy as np
class InterpBase(object):
    def __init__(self, ip):
        self.ip = np.array(ip, float)
        self.c  = np.array([ip]).T - self.ip
        ran = range(len(ip))
        self.c[ran, ran] = 1

    def valuesm(self, iy):
        z = np.array([iy]).T - self.ip
        n = len(self.ip)
        m = np.array([iy]).shape[-1]
        ran = range(n)
        zs = np.zeros((n, m, n))
        zs[ran] = z
        zs[ran, :, ran] = 1
        tmp = np.einsum('ijk->jik', zs, optimize=True)
        return (tmp/self.c).prod(axis=2)

注意:

事实上$\mathbb{H}$矩阵乘以向量, 也就是多项式多点赋值问题, 可以有更快速的办法.

但是都比较复杂, 对于小规模问题, 感觉没什么必要, 可点我查看参考文献

• [1] Alexander Kobe, Michael Sagralof, Fast Approximate Polynomial Multipoint Evaluation
  and Applications

• [2] Justine Gauthier, Fast Multipoint Evaluation On n
  Arbitrary Points

1.2.2. 计算矩阵D

根据求导公式

$h'_j(y_i)=h_j(y_i)\left[\frac{1}{y_i-x_1}+\frac{1}{y_i-x_2}+\cdots +\frac{1}{y_i-x_{j-1}}+\frac{1}{y_i-x_{j+1}}+\cdots \frac{1}{y_i-x_{n}}\right]$

可分两种情况进行求解

1. $y_i=x_i,\forall i=1,2,\cdots n$, 且 $m=n$,
   此时,
   $$h'_j(y_i)=\begin{cases} \left[\frac{1}{y_i-x_1}+\cdots +\frac{1}{y_i-x_{j-1}}+\frac{1}{y_i-x_{j+1}}+\cdots \frac{1}{y_i-x_{n}}\right], & i=j\\[1em] \frac{(y_i-x_1)(y_i-x_2)\cdots(y_i-x_{i-1})(y_i-x_{i+1})\cdots(y_i-x_{j-1})(y_i-x_{j+1})\cdots(y_i-x_{n})}{(x_j-x_1)(x_j-x_2)\cdots(x_j-x_{j-1})(x_j-x_{j+1})\cdots(x_j-x_{n})} & i\neq j\\ \end{cases}$$
   我们记 $\widetilde{\mathbb{D}}=[h'_j(x_i)]_{i=1:n}^{j=1:n}$, 则有表达式
   $$\widetilde{\mathbb{D}}=\left(\dfrac{\mathbb{R}}{\mathbb{C}}\right)_{最后一个轴累乘, 并替换对角线元素}$$
   这里矩阵$\mathbb{C}=[c_{jk}]_{j=1,n}^{k=1,n}$, 以及三维矩阵$\mathbb{R}=[r_{ijk}]_{i=1:\red{m}}^{j,k=1:n}$, 其中
   $$c_{jk}= \begin{cases} x_j-x_k, & \quad j\neq k\\ 1, & \quad j=k \end{cases},\quad \text{以及 } r_{ijk} = \begin{cases} x_i-x_k, & \quad j\neq k 或 i\\ 1, & \quad j=k 或 i \end{cases} .$$
   注意: $\sum_{j} h'_j(y_i)=0, \forall i$, 因此, 矩阵$\widetilde{\mathbb{D}}$的每行元素之和为0, 利用该性质可以很容易计算出$\widetilde{\mathbb{D}}$的对角线上的元素.

2. $y_i=x_j,\forall , i,j$. 这种情况如果直接计算, 计算量太大, 这里给出一种比较简洁快速的方案, 假设插值点的纵坐标构成向量为$\bm v$, 于是$\widetilde{\mathbb{D}}\bm v$所得到的向量就是插值函数的导函数在插值点$\bm x$处的值, 因此再做一次插值计算, 就可以得到插值函数的导数在$\bm y$处的值. 于是有公式$\mathbb{D}\bm{v}=\mathbb{H}\widetilde{\mathbb{D}}\bm v$ , 也就是,
   $$\mathbb{D}=\mathbb{H}\widetilde{\mathbb{D}}.$$

快速简洁的代码

import numpy as np
def interpBase(ip, iy, pureComputeH=True):
    c  = np.array([ip]).T - ip
    n = len(ip)
    ran = np.arange(n)
    c[ran, ran] = 1
    z = np.array([iy]).T - ip
    zs = np.tile(z, (n, 1, 1))
    zs[ran, :, ran] = 1
    zs = zs.transpose((1,0,2))
    abs_sort = lambda x: np.sort_complex(abs(x)+ 1j*x).imag 
    # abs_sort = lambda x: x
    matrixH = (abs_sort(zs)/abs_sort(c)).prod(axis=-1)
    if pureComputeH:
        return matrixH
    dzs = np.tile(c, (n, 1, 1))
    dzs[ran, :, ran] = 1
    dzs[ran, ran, :] = 1
    dzs = dzs.transpose((1,0,2))
    matrixD = (abs_sort(dzs)/abs_sort(c)).prod(axis=-1)
    matrixD[ran, ran] = 0
    matrixD[ran, ran] = -matrixD.sum(axis=-1)
    return matrixH, matrixH.dot(matrixD)

abs_sort 是绝对值排序, 但是实际用起来并不能明显提高精度, 最后给出测试代码.

# Compute LGL points and weights
from scipy import special as spc
deg = 60
x, _ = spc.roots_jacobi(deg-1, 1, 1)
x = np.hstack(([-1], x, [1]))
Lv = spc.eval_legendre(deg, x)
w = 2.0/(deg*(deg+1)*Lv**2)

# Compute matrix D 
a = np.array([Lv]).T / Lv
b = np.array([x]).T - x
ran = np.arange(deg+1)
a[ran, ran] = 0
b[ran, ran] = 1
D = a/b
D[0, 0] = - deg*(deg+1)/4.0
D[-1, -1] = deg*(deg+1)/4.0

D3 = interpBase(x, x, False)
err = abs(D - D3[-1])
print(err.max())

1.2.3. 计算矩阵D不稳定快速算法

如果完全硬算, 计算量太大, 而且也比较复杂. 这里提供一个矩阵向量乘积的快速算法, 也就是计算出$\mathbb{D}\bm u$. 需要借助FFT快速计算. 但是实际计算中该算法会丢精度.

• 先调用前面的程序计算插值多项式在$\bm w$($\bm w=\{w_i\}_{i=1}^n$为$x^n=1$的所有根)的节点值$\bm v$

• 进而通过iFFT得到插值多项式的系数(默认系数按升幂排列), $\bm a=\operatorname{iFFT}(\bm v)$

• 算出插值多项式的导数的系数 $\bm b=(\bm a[1:]*[1,2,\cdots,n-1]) 追加 [0]$

• 根据插值多项式的导数的系数, 计算出导数在$\bm w$处的值$\bm c=\operatorname{FFT}(\bm b)$

• 最后通过$\bm w, \bm c$ 与 $\bm u$, 调用前面的程序, 可计算出$\mathbb{D}\bm u=\mathbb{H}\bm c$

全部Python代码如下:

class InterpBase(object):
    def __init__(self, ip, type=np.float):
        self.ip = np.array(ip, type)
        self.c  = np.array([ip]).T - self.ip
        ran = range(len(ip))
        self.c[ran, ran] = 1

    def valuesm(self, iy, type=np.float):
        z = np.array([iy]).T - self.ip
        n = len(self.ip)
        m = np.array([iy]).shape[-1]
        ran = range(n)
        zs = np.zeros((n, m, n),type)
        zs[ran] = z
        zs[ran, :, ran] = 1
        tmp = np.einsum('ijk->jik', zs, optimize=True)
        return (tmp/self.c).prod(axis=2)

def DMatrixVectorMultiplication(ip, iy, iv):
    # ip 插值点横坐标, iv 插值点纵坐标, iy需要赋值的点的横坐标
    # 直接算出矩阵向量的乘积
    n = len(ip)
    ran = np.arange(n)
    t = 2*np.pi/n
    w  = np.cos(t) - 1j*np.sin(t)
    w = w**ran
    v = InterpBase(ip).valuesm(w, np.complex).dot(iv)
    b = ifft(v)*ran
    c = np.zeros_like(b)
    c[:-1] = b[1:]
    c = fft(c)
    return InterpBase(w, np.complex).valuesm(iy, np.complex).dot(c).real

def DMatrix(ip, iy):
    # ip 插值点横坐标, iv 插值点纵坐标, iy需要赋值的点的横坐标
    # 只算出矩阵
    n = len(ip)
    ran = np.arange(n)
    t = 2*np.pi/n
    w  = np.cos(t) - 1j*np.sin(t)
    w = w**ran
    v = InterpBase(ip).valuesm(w, np.complex)
    b = ran*ifft(v.T)
    c = np.zeros_like(b)
    c[:, :-1] = b[:, 1:]
    c = fft(c).T
    return InterpBase(w, np.complex).valuesm(iy, np.complex).dot(c).real

注意:

• 上面最后一行代码中的矩阵, 是可以事先存起来的

  InterpBase(w, np.complex).valuesm(iy, np.complex)

• InterpBase(ip)相关的系数也可以事先存起来, 它需要计算两次复矩阵乘法!

• DMatrixVectorMultiplication与DMatrix是独立的两个函数, 它们都依赖InterpBase

2. 其它笔记(快速计算多项式乘法)

同样的道理, 可以借助FFT实现多项式乘法的快速计算

例如已知
$f(x)=a_0+a_1x+\cdots a_n x^m$
$g(x)=b_0+b_1x+\cdots b_n x^n$

求$h(x)=f(x)g(x)=c_0+c_1x+\cdots c_n x^{m+n}$的系数$\{c_k\}_{k=0}^{m+n}$.

快速算法的核心思想为:

• 先用FFT分别算出$f(x),g(x)$在$m+n+1$次单位根上的值

• 上步结果做乘积, 就是$h(x)$在单位根上的值

• 最后通过iFFT得到多项式$h(x)$的系数

Python 代码例如下

import numpy as np
from scipy.fftpack import fft,ifft

# 多项式乘积快速算法

a = np.array([1, 2, 0])
b = np.array([3, 4, 0])

c = ifft(fft(a)*fft(b))
print(c)


a = [1, 2, 3, 4]
b = [5, 6, 7, 8]

tz = [0]*(len(a)-1)
a = np.array(a + tz)
b = np.array(b + tz)

c = ifft(fft(a)*fft(b))
print(c.real)

注意:

• 向量需要做适当的用0扩充

• 向量按照多项式系数升幂排列