三次B样条曲线插值的LUTS-PIA算法

LUTS-PIA Algorithm for Cubic B-Spline Curve Interpolation

摘要: 本文研究了三次B样条曲线插值问题。首先，我们将配置矩阵进行下上三角分裂，然后基于该下上三角分裂提出了(Lower Upper Triangular Splitting-Progressive Iterative Approximation) LUTS-PIA算法，并证明了该算法的收敛性。最后，数值实验结果表明：LUTS-PIA算法明显优于(Hermitian and Skew-Hermitian Splitting-Progressive Iterative Approximation) HSS-PIA算法。

Abstract: The cubic B-spline curve interpolation problem is studied in this paper. First, we split the allocation matrix into lower and upper triangular parts called lower upper triangular splitting (LUTS). Based on this LUTS, we propose lower upper triangular splitting-Progressive Iterative Approximation (LUTS-PIA) algorithm and prove its convergence. Finally, we test some numerical experiments which show that LUTS-PIA has a better convergence behavior than the HSS-PIA.

1. 引言

渐进迭代逼近法(PIA) [1] 是曲线插值中一种常用的数值方法，该方法的基本思想是让给定点列作为初始控制顶点序列，构造一条初始混合曲线，通过迭代调整它的控制顶点，使得这条函数曲线或曲面在控制点列与原始数据点列之间逼近。该方法具有清晰的几何意义和很强的数值稳定性。更重要的是，避免了求解线性方程组，而且当插值点增加的时候我们只需要部分分量的矫正。然而，当插值点数量增加到一定程度时，由一些归一化全正基(NTP)产生的配置矩阵往往是病态的，得到的配置矩阵并不总是正定的，这促使着研究者们继续深入研究PIA。

为了提高PIA的收敛率，陆利正 [2] 提出了一种加权的PIA (WPIA)。刘晓艳 [3] 和王志好 [4] 基于非均匀B样条配置矩阵的Jacobi分裂和GS分裂，分别提出了Jacobi-PIA和GS-PIA。另外，刘成志等人 [5] 提出了预处理的PIA (PPIA)。最近Hu L等人 [6] 基于归一化全正基配置矩阵Hermitian and Skew-Hermitian (HSS)分裂提出了HSS-PIA。以上方法能明显改进PIA的收敛速度。

本文将研究三次B样条基配置矩阵的LU分裂文献 [7] 中的LUTS-PIA。我们证明了LUTS-PIA算法的收敛性，通过数值实验表明，本文所得的LUTS-PIA算法在具有正对角元严格的或不可约的对角占优矩阵分裂中优于HSS-PIA算法。

2. LUTS-PIA算法

设 ${\left\{b_i\left(t\right)\right\}}_{i=0}^n$ 为三次B样条基，给定一待插值有序点集 ${\left\{v_i\right\}}_{i=1}^n\in R^3$ 其中每个点 $v_i$ 都赋予一个参数值 $t_i$ 且参数值满足 $t_0<t_1<\cdots <t_n$ ，以 ${\left\{v_i\right\}}_{i=1}^n$ 这个点列为初始控制点，生成一条初始曲线

${\gamma }^{\left(0\right)}\left(t\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{w}_i^{\left(0\right)}}{b}_i\left(t\right)$ ， $w_i^{\left(0\right)}=v_i$ ， $i=0,1,\cdots ,n$ ，

假定第k次迭代后生成的曲线为

${\gamma }^{\left(k\right)}\left(t\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{w}_i^{\left(k\right)}}{b}_i\left(t\right)$ 。

为进行第 $k+1$ 次迭代，首先计算校正向量

${\delta }_i^{\left(k\right)}=v_i-{\gamma }^{\left(k\right)}\left(t_i\right)$ ， $i=0,1,\cdots ,n$ ，

并令

$w_i^{\left(k+1\right)}=w_i^{\left(k\right)}+{\delta }_i^{\left(k\right)}$ ， $i=0,1,\cdots ,n$ ，

进而得到第 $k+1$ 次曲线

${\gamma }^{\left(k+1\right)}\left(t\right)={\displaystyle \underset{i=0}{\overset{n}{\sum }}{w}_i^{\left(k+1\right)}}{b}_i\left(t\right)$ 。

可以得到

$w_i^{k+1}=w_i^k+v_i-{\displaystyle \underset{j=0}{\overset{n}{\sum }}{w}_j^k}{b}_j\left(t_i\right)$ 。 (1)

若 $\underset{x\to \infty }{\lim }{r}^k\left(t_i\right)=v_i$ 该方程(1)称为PIA格式，PIA适用于任何归一化全正基(NTP)。

记 $W^{\left(k\right)}={\left[w_0^{\left(k\right)},w_1^{\left(k\right)},\cdots ,w_n^{\left(k\right)}\right]}^{\text{T}}$ ， $V={\left[v_0,v_1,\cdots ,v_n\right]}^{\text{T}}$ ，则用矩阵形式表示如下：

$W^{\left(k+1\right)}=W^{\left(k\right)}+\left(V-B{W}^{\left(k\right)}\right)$ 。

其中，B等于

$\left[\begin{array}{cccc}{b}_0\left(t_0\right)& b_1\left(t_0\right)& \cdots & b_n\left(t_0\right)\\ b_0\left(t_1\right)& b_1\left(t_1\right)& \cdots & b_0\left(t_0\right)\\ ⋮& ⋮& & ⋮\\ b_0\left(t_n\right)& b_1\left(t_n\right)& \cdots & b_n\left(t_n\right)\end{array}\right]$

称为配置矩阵，这就是经典的PIA。

实际上，PIA算法在数学上等同于求解配置矩阵方程

$BW=V$ 。 (2)

在本文中，设 $B=D+\hat{L}+\hat{U}$ ，其中D是配置矩阵B的对角元组成的矩阵 $\hat{L}$ 和 $\hat{U}$ 分别是配置矩阵B的严格下三角部分和严格上三角部分组成的矩阵，因此我们可以将B分裂成

$B=L+U$ ， (3)

其中， $L=D_1+\hat{L}$ ， $U=D_2+\hat{U}$ ， $D=D_1+D_2$ 。

那么，由

$B=\left(\alpha I+L\right)-\left(\alpha I-U\right)$ 和 $B=\left(\alpha I+U\right)-\left(\alpha I-L\right)$ ，

我们就可以得到了B的LUTS并进而可以提出如下的LUTS-PIA。

3. LUTS-PIA算法

可以得到

$\{\begin{array}{l}{w}_i^{\left(k+\frac{1}{2}\right)}=w_i^{\left(k\right)}+d_i^{\left(k\right)},d_i^{\left(k\right)}={\left(\alpha I+L\right)}^{-1}{\delta }_i^{\left(k\right)}\\ w_i^{\left(k+1\right)}=w_i^{\left(k+\frac{1}{2}\right)}+d_i^{\left(k+\frac{1}{2}\right)},d_i^{\left(k+\frac{1}{2}\right)}={\left(\alpha I+U\right)}^{-1}{\delta }_i^{\left(k+\frac{1}{2}\right)}\end{array}$ ， (4)

其中 ${\delta }_i^k={\gamma }^{\left(k\right)}\left(t\right)-{\displaystyle \underset{i=0}{\overset{n}{\sum }}{w}_i^{\left(k\right)}}{b}_i\left(t\right)$ ， $k=0,\frac{1}{2},1,\cdots$ ，我们称(4)为曲线的LUTS-PIA格式。

记 $W^{\left(k\right)}={\left[w_0^{\left(k\right)},w_1^{\left(k\right)},\cdots ,w_n^{\left(k\right)}\right]}^{\text{T}}$ ， ${\Delta }^{\left(k\right)}={\left[{\delta }_0^{\left(k\right)},{\delta }_1^{\left(k\right)},\cdots ,{\delta }_n^{\left(k\right)}\right]}^{\text{T}}$ ， $D^{\left(k\right)}={\left[d_0^{\left(k\right)},d_1^{\left(k\right)},\cdots ,d_n^{\left(k\right)}\right]}^{\text{T}}$ ，则方程(4)的矩阵形式表示如下：

$\{\begin{array}{l}{W}^{\left(k+\frac{1}{2}\right)}=W^{\left(k\right)}+D^{\left(k\right)},D^{\left(k\right)}={\left(\alpha I+L\right)}^{-1}{\Delta }^{\left(k\right)}\\ W^{\left(k+1\right)}=W^{\left(k+\frac{1}{2}\right)}+D^{\left(k+\frac{1}{2}\right)},D^{\left(k+\frac{1}{2}\right)}={\left(\alpha I+U\right)}^{-1}{\Delta }^{\left(k+\frac{1}{2}\right)}\end{array}$ ， (5)

这就是LUTS-PIA算法，它由两个半迭代组成，取代了PIA中每个控制点的迭代长度。

显然，通过直接代换可以有效地得到下三角矩阵 $\left(\alpha I+L\right)$ 和上三角矩阵 $\left(\alpha I+U\right)$ 的精确解。

注：我们观察到矩阵 $\left(\alpha I+L\right)$ 和 $\left(\alpha I+U\right)$ 能够保持原有矩阵B的稀疏性，而HSS方法中的H和S两个矩阵可能是稠密的，尽管原有矩阵B是稀疏的，可以参见 [8] 。

4. LUTS-PIA算法的收敛性分析

首先介绍一些基本的定义、符号和文中将使用的预备知识。在整篇文章中，我们用 $B^{*}$ 表示所有n维复向量B的共轭转置，用 $C^{n\times n}$ 表示所有 $n\times n$ 个复矩阵的共轭转置，用 ${\Vert \cdot \Vert }_2$ 表示2-范数。

设矩阵 $B=\left(b_{ij}\right)\in C^{n\times n}$ 则

1) B正定，如果 $x^{*}Bx>0$ 对于所有非零的 $x\in C^n$ 都成立，或者等价地，B的Hermitian部分 $\left(B+B^{*}\right)/2$ 是Hermitian正定；

2) 如果按行对角占优，则对于 $i=1,\cdots ,n$ ， $\left|b_{ii}\right|\ge {\displaystyle \underset{j=1,j\ne i}{\overset{j=n}{\sum }}\left|b_{ij}\right|}$ ；如果按列对角占优，则对于 $j=1,\cdots ,n$ ， $\left|b_{jj}\right|\ge {\displaystyle \underset{i=1,i\ne j}{\overset{i=n}{\sum }}\left|b_{ij}\right|}$ ；

3) 如果 $\left|b_{ii}\right|\ge {\displaystyle \sum \left|b_{ij}\right|}\left(i\ne j\right)\left(\left|b_{jj}\right|\ge {\displaystyle \sum \left|b_{ij}\right|}\left(i\ne j\right)\right)$ 则称B按行(按列)弱对角占优。若式中每一个不等式都是严格不等式(>)，则称B按行(按列)严格对角占优；

4) 如果B是按行(按列)不可约的，则B按行(按列)弱对角占优，并且对角占优不等式至少一个i值(一个j值)是严格的。

现在我们讨论LUTS-PIA算法的收敛性。注意到LUTS-PIA算法(5)中是两个半步分裂迭代，可以整合成单步迭代，即从(5)中消去 $w^{\left(k+\frac{1}{2}\right)}$ ，得到以下迭代

$w^{k+1}=Q_{\alpha }{w}^{\left(k\right)}+P_{\alpha }^{-1}v$ 。 (6)

其中，

$Q_{\alpha }={\left(\alpha I+U\right)}^{-1}\left(\alpha I-L\right){\left(\alpha I+L\right)}^{-1}\left(\alpha I-U\right)$ (7)

并且

$P_{\alpha }=\frac{1}{2\alpha }\left(\alpha I+L\right)\left(\alpha I+U\right)$ ， (8)

这个迭代可以看成是分裂 $B=P_{\alpha }-P_{\alpha }{Q}_{\alpha }$ 诱导的。

为了证明(6)的收敛性，我们引入下面的引理。

引理1 设

$H\left(\alpha \right)=\left(\alpha I-A\right){\left(\alpha I+A\right)}^{-1}$ 。

如果 $A\in C^{n\times n}$ 是一个正定矩阵，则对于 $\alpha >0$ 有 ${\Vert H\left(\alpha \right)\Vert }_2<1$ 。

如果L和U是正定的， $Q_{\alpha }$ 如(6)中定义那样，利用引理1和矩阵谱的相似不变性。则如下

$\rho \left(Q_{\alpha }\right)\le {\Vert \left(\alpha I-L\right){\left(\alpha I+L\right)}^{-1}\Vert }_2{\Vert \left(\alpha I-U\right){\left(\alpha I+U\right)}^{-1}\Vert }_2\equiv \sigma \left(\alpha \right)<1$ ，

成立，我们有以下收敛结果。

引理2 设配置矩阵LU分裂为 $B=L+U$ 。如果L和U两个都是正定的，则对于所有 $\alpha >0$ ，LUTS-PIA 算法(4)收敛到(2)的唯一解。

一般来说，迭代法求解出使谱半径达到最小的最优参数值 $\alpha$ ，目前还没有研究出求解谱半径达到最小的最优参数值，一般使用使谱半径上界达到最小时的最优参数 ${\alpha }^{\ast }$ 来近似，为了找出 ${\alpha }^{\ast }$ 作为 $\alpha$ 的一个粗略估计，我们定义

$\underset{{\lambda }_i\in \lambda \left(L\right)}{\max }\left|\frac{\alpha -{\lambda }_i}{\alpha +{\lambda }_i}\right|\underset{{\lambda }_j\in \lambda \left(U\right)}{\max }\left|\frac{\alpha -{\lambda }_j}{\alpha +{\lambda }_j}\right|\equiv \sigma \left(\alpha \right)$ 。

则 $\rho \left(Q_{\alpha }\right)\le \sigma (\partial )\le {\Vert \left(\alpha I-L\right){\left(\alpha I+L\right)}^{-1}\Vert }_2{\Vert \left(\alpha I-U\right){\left(\alpha I+U\right)}^{-1}\Vert }_2$ 或者 $\sigma (\partial )\le \rho \left(Q_{\alpha }\right)$ ，对于任意的 $\alpha >0$ ，我们可以把 $\sigma \left(\alpha \right)$ 作为 $\rho \left(Q_{\alpha }\right)$ 的最优 $\alpha$ 的近似。

${\lambda }_i$ 与 ${\lambda }_j$ 是正的，其中 $i,j=1,2,\cdots ,n$ ，因此对于任意的 $\alpha >0$ ，都有 $\rho \left(Q_{\alpha }\right)<1$ 。我们可以看到如果将矩阵L和U的特征值放在一个集合 $\kappa$ 中，即 $\kappa =\left[{\lambda }_{\min },{\lambda }_{\max }\right]$ ，则我们可以将 $\underset{\lambda \in \kappa }{\max }\frac{{\left(\alpha -\lambda \right)}^2}{{\left(\alpha +\lambda \right)}^2}$ 作为 $\sigma \left(\alpha \right)$ 的估计，其中 ${\lambda }_{\min }$ 与 ${\lambda }_{\max }$ 是矩阵L和U整体特征值的最小值和最大值。

定理3 $\rho \left(Q_{\alpha }\right)$ 相对于所有的正 $\alpha$ 的最小值在

${\alpha }^{\ast }=\sqrt{{\lambda }_{\min }{\lambda }_{\max }}$ (9)

对应的最小值为

$\sigma \left({\alpha }^{\ast }\right)={\left(\frac{\sqrt{{\lambda }_{\max }}-\sqrt{{\lambda }_{\min }}}{\sqrt{{\lambda }_{\max }}+\sqrt{{\lambda }_{\min }}}\right)}^2\le \frac{\sqrt{{\lambda }_{\max }}-\sqrt{{\lambda }_{\min }}}{\sqrt{{\lambda }_{\max }}+\sqrt{{\lambda }_{\min }}}$ 。

证明：定理3的证明可以由 Bai，Golub和Ng等人 [8] 中的推论2.3得出。

注：我们可以得到使迭代矩阵谱半径 $\rho \left(Q_{\alpha }\right)$ 的上界最小值 $\sigma \left(\alpha \right)$ 的最优参数 ${\alpha }^{\ast }$ ，即

${\alpha }^{\ast }=\underset{\alpha >0}{\arg \min }\sigma \left(\alpha \right)=\sqrt{{\lambda }_{\min }{\lambda }_{\max }}$ 。

但是没有最小谱半径，即

${\alpha }^{\ast }\ne \underset{\alpha >0}{\arg \min }\rho \left(Q\left(\alpha \right)\right)\equiv \alpha$ 。

特别地，相倩等人提出的LU分裂迭代法并没有指定D₁和D₂的构造。因此，本文给出了一个LU分裂迭代中D₁和D₂的构造思想，当B是具有正对角元严格的或不可约的对角占优矩阵时，这样的分裂构造使得下三角矩阵L和上三角矩阵U都是正定的。接下来，为了确保L和 $U\in C^{n\times n}$ 是正定的，我们给出D₁和D₂的构造。矩阵 $B=\left(b_{ij}\right)\in C^{n\times n}$ 并且 $D=\text{diag}\left(b_{11},\cdots ,b_{nn}\right)$ ，其中 $D_1=\text{diag}\left(d_1,\cdots ,d_n\right)$ ，我们设

$d_i=\frac{1}{2}\left(b_{ii}+r_i-c_i\right)$ ，

其中 $r_i={\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ij}\right|}$ 和 $c_i={\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ji}\right|}$ ，并且 $D_1=D-D_2$ 。

众所周知，这类矩阵经常出现在各种各样的领域，包括偏微分方程的数值解、经济学中的生产和增长模型、运筹学中的线性互补问题和随机分析中的马尔可夫链，可以参见 [9] [10] [11] 及其他参考文献。

定理4 设B是一个按行(按列)具有正对角元严格的或不可约的对角占优矩阵，且D₁和D₂如前面所构造的一样，对任意的 $\alpha >0$ ，都有LUTS-PIA算法(4)收敛到方程(2)的唯一解。

证明：根据引理2可以证明L和U是正定的。如果 $\hat{L}$ 和 $\hat{U}$ 分别是配置矩阵B的严格下和上三角部分，则

$L=D_1+\hat{L}$ 和 $U=D_2+\hat{U}$ ，其中 $D_1+D_2=D$ 。

首先，我们观察到配置矩阵的对角项都是正的。根据对配置矩阵B的假设，与对L和U的构造，我们得出 $L+L^{\ast }$ 和 $U+U^{\ast }$ 具有正对角项，并且都是严格或不可约对角占优的。因此， $L+L^{\ast }$ 和 $U+U^{\ast }$ 的所有特征值都是正的，参见 [12] 的定理6.1.10和推论6.2.27。

当矩阵B按列对角占优时

$\begin{array}{c}{\left(L+L^{\ast }\right)}_{i,i}=2d_i=b_{ii}+r_i-c_i=b_{ii}+{\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ij}\right|}-{\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ji}\right|}\\ =b_{ii}+\left({\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ij}\right|}+{\displaystyle \underset{j=i+1}{\overset{n}{\sum }}\left|b_{ji}\right|}\right)-\left({\displaystyle \underset{j=1}{\overset{i-1}{\sum }}\left|b_{ji}\right|}+{\displaystyle \underset{j=i+1}{\overset{n}{\sum }}\left|b_{ji}\right|}\right)\\ =\left(b_{ii}-{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}\left|b_{ji}\right|}\right)+{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}\left|{\left(L+L^{\ast }\right)}_{i,j}\right|}\\ >{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}\left|{\left(L+L^{\ast }\right)}_{i,j}\right|}\end{array}$ ，

$L+L^{\ast }$ 是严格的对角占优。

类似地，当矩阵B按行对角占优时

${\left(U+U^{\ast }\right)}_{i,i}=2\left[b_{ii}-\frac{1}{2}\left(b_{ii}+r_i-c_i\right)\right]=b_{ii}+c_i-r_i>{\displaystyle \underset{j=1,j\ne i}{\overset{n}{\sum }}\left|{\left(U+U^{\ast }\right)}_{i,j}\right|}$ ，

$U+U^{\ast }$ 是严格的对角占优，证毕。

5. 数值实验

在本节中，我们利用三次B样条基文献来证明LUTS-PIA算法求解配置矩阵方程(2)的有效性。为了进行比较，我们还测试了HSS-PIA算法。数值实验在Matlab环境中实现，以待插值点集作为初始的控制顶点集，终止条件为 ${\epsilon }^k\le tol$ ，其中，k表示迭代步数，在第k步得到的拟合误差，其定义为

${\epsilon }^k=\frac{\Vert b-A{x}^k\Vert }{\Vert b-A{x}^0\Vert }$ 。

当对 $\alpha$ 进行选取时，我们将按照定理3中的公式计算，可以得出 $\alpha$ ，由于数值的特性，得出的 $\alpha$ 并不是精确的最优值。通过数值实验发现，最优的 $\alpha$ 在计算最优值的附近选取。表1中，给出了HSS-PIA算法和LUTS-PIA算法的数值结果，其中，“n”表示矩阵阶数，“ $\epsilon$ ”表示迭代误差，“CPU”表示迭代时间(s)，“IT”表示迭代次数。

例1. 螺旋线的参数方程为

$\{\begin{array}{l}x=r\times \cos \left(t\times \left({\theta }_2-{\theta }_1\right)+{\theta }_1\right)\hfill \\ y=r\times \sin \left(t\times \left({\theta }_2-{\theta }_1\right)+{\theta }_1\right)\hfill \\ z=h\times t\hfill \end{array}$ 。

刘成志等人提出的三次均匀B样条扩展曲线的渐进迭代逼近法得出配置矩阵与形状参数 $\lambda$ 有关，当 $\lambda \in \left[-2,1\right]$ 时，迭代矩阵的谱半径随着 $\lambda$ 的增加而减少，从而当 $\lambda =1$ 时，具有最快的收敛速度，可以参见文献 [13] 的定理1。因此，在本文中我们令形状参数 $\lambda =1$ ， $r=30$ ， $h=50$ ， $t\in \left[0,1\right]$ ， ${\theta }_1=pi/6$ ， ${\theta }_2=20\times pi$ ，得出配置矩阵。

随着阶数n的增加，其渐进迭代逼近误差、运算时间和迭代次数分别如表1所示。

Table 1. Uniform cubic B-spline helix interpolation

表1. 均匀三次B样条螺旋线插值

注：对于非均匀三次B样条螺旋线插值，如果生成的配置矩阵是广义对角占优的，则直接使用LUTS-PIA算法不能保证该算法的收敛性，如果我们很容易找到一个正对角阵D，使得BD严格对角占优，然后使用上述算法再继续求解。

6. 结论

本文给出了对角占优的矩阵B对角线的一个特殊分裂，分裂能够保持原有矩阵B的稀疏性，提高收敛速度。利用三次B样条基 [10] [14] 来研究LU-PIA算法和HSS-PIA算法的计算复杂度。数值结果表明：无论是收敛速度还是计算复杂度，LU-PIA算法解配置矩阵方程(2.3)的效果都好于HSS-PIA算法。并且，如果配置矩阵B是对称正定的或者Heimitian正定的，HSS-PIA就不能解决Heimitian正定的情况，因为HSS-PIA该算法不存在HSS分裂。

我们的方法可以应用于径向基函数，还可以推广到三次B样条基函数下的曲面插值，曲线(曲面)逼近等求解问题。而最优参数的估计还有待进一步的深入研究。

基金项目

湖南省研究生科研创新项目(CX20220953)。

NOTES

^*通讯作者。