精度曲面算法

1引言

高精度曲面建模理论（HASM）根据曲面论基本定理，结合高斯－科达齐方程，对所模拟的区域进行均匀正交剖分建立数值方程，是一种全新的曲面建模方法。其在采样数据的约束下，对正交剖分的均匀格网点进行求解，从而获得高精度的拟合曲面［1］。数值实证表明，HASM的计算精度高于以往的建模方法［2］。为了解决高精度曲面建模方法的速度问题和超大计算量问题，在大量数值实验的基础上，岳天祥等人给出了HASM的最佳表达形式［3－4］。HASM的计算速度及模拟精度与以往相比有了大幅度的提高［4］。HASM计算的整个过程大致可分为3个阶段：HASM方程组系数矩阵的形成、采样方程的建立和代数方程组的求解。为了进一步解决HASM模型的计算速度问题，我们对模型中的线性方程组，采用不同的预处理共轭梯度法，即不完全Cholesky分解预共轭梯度法（）和对称逐步超松弛预共轭梯度法（Method）求解，试验结果表明：与以往预处理方法相比，这两种方法均大大降低了CPU时间与收敛到一定精度解的迭代步数。同时，考虑了稀疏矩阵的压缩存储方式，这对大规模问题而言，解决了以往预处理共轭梯度算法中的内存不足等问题。

2预处理共轭梯度法

设HASM模型中的方程组为：Ax＝b（1）其中，A∈Rn×n为对称正定稀疏矩阵。对于方程组（1）的求解，分为直接法和迭代法。直接法针对中小规模的问题，典型的直接法为高斯消元法，但此方法的时间复杂度与模拟区域网格数的三次方成正比［5］。在迭代法中，共轭梯度（CG）方法［6］被广泛用于求解大型稀疏对称正定线性系统。对于大型稀疏对称正定系统Ax＝b，共轭梯度法（CG）在理论上具有有限步收敛性［7］。特别是当系数矩阵的绝大部分特征值在1附近或者该问题良态时，CG方法的收敛速度很快。即CG方法的收敛速度依赖于系数矩阵A的条件数或更一般地依赖于A的特征值的分布。但在实际问题中，由于舍入误差的影响，以往经验表明，为了达到比较准确的解，CG方法的收敛步数往往远远大于系数矩阵的阶数甚至不收敛。而对于HASM模型，此方程组的系数矩阵条件数较大。因此，在实际计算前，先通过改善系数矩阵的条件数来提高CG方法的收敛速度，即预处理共轭梯度法。预处理共轭梯度法的思想就是把共轭梯度法用到预处理后的代数系统。预处理共轭梯度法算法如下［7］：该算法中，主要的计算量来自于求解方程组Mzk＝rk。因此，下面我们寻找预处理矩阵M，一方面使得该预处理方法的收敛速度很快，另一方面使求解方程组Mzk＝rk的计算量尽可能小。2．1不完全Cholesky分解预处理方法将方程组（1）转化为易于求解的方程组：A珚x＝b珔，（2其中，A珚＝MA，b珔＝Mb。若k（A珚）k（A）（k（A）表示矩阵A的条件数）则同样的算法求解式（2）比求解式（1）具有更快的收敛速度。首先，我们采用的预处理方法为A的不完全Cholesky分解法［8］。设A的不完全Cholesky分解为A＝M＋R＝LLT＋R，其中，L是下三角阵，使M＝LLT尽可能接近A，且L保持跟A一样的稀疏性或具有其他指定的稀疏性。完全Cholesky分解是对系数矩阵A进行三角分解A＝LLT，不完全分解是对矩阵A－R进行三角分解LLT。由于矩阵R可以变化，因此，L的稀疏性结构可以预先适当控制，即L中哪些元素为0，可以预先规定，同时还要考虑到LLT要尽可能地接近A。这样克服了完全Cholesky分解破坏A的稀疏性的缺点。在实际计算中，常常考虑R有较多的零元素，且R元素不应太大。本文我们考虑没有填充的Cholesky分解算法，即对A进行Cholesky分解时，对应于A的零元处的位置，在分解过程中不再引入非零元。下面给出此算法的代码［7］：以M＝LLT作为预处理算子，应用共轭梯度法，不难验证预处理后的系数矩阵A珚≈I。在具体实现过程中，A采用行压缩存储，Mzk＝rk可转化为求解LTzk＝yk，及Lyk＝rk。这样，解方程组的计算量为O（n2），比直接求解Mzk＝rk的计算量提高了一个数量级。ICCG方法对小型问题是有效的，针对大型问题，我们可以预先对系数矩阵进行Cholesky分解，将分解后的下三角矩阵参与HASM模型的计算，作为HASM模型的输入参数。该方法的空间复杂度与时间复杂度均为O（n2）。

2．2对称逐步超松弛预处理方法HASM的主要计量来自于矩阵的乘积运算和求逆，且求逆运算所花费的时间远远大于矩阵的乘积。上述的不完全Cholesky分解预处理方法中每次内迭代都需要求解方程组Mzk＝rk。陈传法等人［9］采用了改进的Gauss－Seidel算法对HASM进行求解。Gauss－Seidel法的收敛性问题受到HASM方程中稀疏矩阵的制约。对称逐步超松弛（SSOR）方法是Gauss－Seidel的一种加速方法，是求解大型稀疏矩阵方程组的有效方法，它具有计算公式简单，程序设计容易，占用计算机内存少等优点。同时，该预处理方法可通过直接计算系数矩阵的近似逆矩阵，以避免每次迭代求解方程组Mzk＝rk，从而减少计算量。该方法中只涉及矩阵与向量的乘积，且很容易实现并行运算。应用对称逐次超松弛方法（SSOR）［10］，设系数矩阵A可以分解为A＝L＋D＋LT，其中D为A的对角线元素构成的对角阵，L为A的下三角部分构成的下三角矩阵。定义SSOR预处理算子如下［11］：3O共从而使得预处理共轭梯度法中每次内迭代求解方程组Mzk＝rk转化为计算zk＝M珨rk，进而减少了计算量。但矩阵M珨的稀疏性要比A的稀疏性差。在具体计算过程中，把计算zk＝M珨rk转化为计算yk＝K珚rk和zk＝K珚Tyk。这样保证了K珚与A有相同的稀疏结构，从而减少了内存占用率。从上述分析可以看出，在程序实现过程中，预处理矩阵K珚可以显式地实现，从而使得SSORCG的空间复杂度为O（n）。

3高斯合成曲面的数值试验

该试验以高斯合成曲面为研究对象，以Dell－OptiPlex990为计算机实验平台验证不同的预处理共轭梯度方法对HASM收敛速度的影响。高斯合成曲面（图1）的数学表达式假定研究区域为［－3，3］×［－3，3］，此时，－6．5510＜f（x，y）＜8．1062。设计了3组试验，分别改变模拟区域网格点数，内外迭代次数以比较不同方法的计算效率。试验1为固定采样间隔（m＝4），外迭代停止准则f（n＋1）－f（n）／f（n）＜＝1e－6，内迭代停止准则为b－A＊x／b－A＊x0＜＝1e－9，改变模拟区域的网格数（分别为101×101、301×301、501×501、1001×1001）比较不同预处理共轭梯度法的计算时间及停止迭代时的内迭代步数。同时给出了ICCG法与对角线预处理共轭梯度法的计算时间的差值与网格数的关系。试验结果记录在表1和图2中。表1中，S表示对角线预处理共轭梯度法，T表示三对角预处理共轭梯度法。由表1和图2表明，ICCG法及SSORCG法的计算时间及内迭代步数均小于对角线预处理方法及三对角预处理方法；且ICCG法与对角线预处理方法的计算时间的差值随着网格数的增多呈现较好的线性关系，即网格数越多，ICCG法的优势越明显，两种方法的时间差值与网格数的关系为：其中，t为计算时间差值，gn为模拟区域的网格数目。试验结果表明，ICCG法与SSORCG法不论在计算时间上还是迭代步数上，较以往的对角线预处理方法与三对角预处理方法均有显著提高。从图2中可见，随着模拟网格数的增多，ICCG法的计算时间越少其优势越明显。试验2为固定采样间隔（m＝4），外迭代次数设为5次，模拟区域网格101×101，比较不同方法的收敛精度。试验3为固定采样间隔（m＝4），设内迭代次数为10次，模拟区域网格数为501×501，比较不同的外迭代次数（分别为2、4、6、8次）下各种方法的收敛精度。试验结果如表2、3所示。综合以上试验可得，ICCG法收敛速度最快，SSORCG法次之，二者均比对角线预处理法、三对角预处理法收敛快；且随着网格数的增多，ICCG法的优势越明显。

4结论

HASM模型由于其收敛速度问题，严重制约了其推广及应用。本文针对该模型的病态性等特点，在共轭梯度法的基础上，分别采用不完全Cholesky分解预处理方法及对称逐步超松弛预处理方法对系数矩阵进行预处理。在具体实现过程中，充分考虑了两种不同的预处理方法，在计算量与存储量方面的实现细节。以高斯合成曲面为研究对象，通过数值试验，在Dell－Optiplex990MT机器上验证了两种预处理方法的有效性。试验结果表明，所采用的两种预处理方法比以往的对角线预处理法及三对角预处理方法在收敛速度上均有明显的改善，进而解决了HASM推广使用中的计算速度问题。