全浩榮，劉成志，李軍成，楊煉，胡麗娟

（湖南人文科技學(xué)院數(shù)學(xué)與金融學(xué)院，湖南婁底 417000）

在計(jì)算機(jī)輔助幾何設(shè)計(jì)中，數(shù)據(jù)插值是一項(xiàng)重要的內(nèi)容。傳統(tǒng)的數(shù)據(jù)插值方法需要求解線性方程組，即配置方程。通常，線性方程組的系數(shù)矩陣是病態(tài)的，需研究高效的求解算法。為此，對(duì)配置方程的求解算法進(jìn)行了研究并取得了不少成果，如基于Neville 消元的直接法［1-3］、基于矩陣分裂的迭代法［4-6］、子空間迭代法［7］等。近年來出現(xiàn)了一種具有幾何意義的迭代算法——漸近迭代逼近（progressive iterative approximation，PIA）法。由于PIA 法具有穩(wěn)定的收斂性、在迭代過程中易修改幾何約束條件等優(yōu)點(diǎn)而備受關(guān)注［8-9］。PIA 法最早源于齊東旭等［10］提出的盈虧修正思想。LIN等［11］證明了所有規(guī)范的全正基函數(shù)曲線（曲面）均具盈虧修正性質(zhì)，并稱之為PIA 性質(zhì)。

雖然PIA 法已有許多成果，但由于其收斂速度慢，在相關(guān)應(yīng)用領(lǐng)域難以發(fā)揮優(yōu)勢。近年來，PIA的加速方法研究逐漸成為幾何計(jì)算中的熱點(diǎn)。如加權(quán)PIA［12］（weighted progressive iterative approximation，WPIA）法、帶互異權(quán)值的PIA法［13-15］、Chebyshev 多項(xiàng)式加速迭代法［16］、帶形狀參數(shù)曲線的PIA法［17］等。由于PIA法等價(jià)于求解配置方程的Richardson 迭代法［6-7］，文獻(xiàn)［18］利用預(yù)處理技術(shù)對(duì)全正基函數(shù)的PIA 法進(jìn)行預(yù)處理，得到收斂速度有明顯提升的預(yù)處理PIA（preconditioned progressive iterative approximation，

PPIA）法。隨后，文獻(xiàn)［19］結(jié)合矩陣Kronecker 積的性質(zhì)將預(yù)處理技術(shù)應(yīng)用于張量積型Bézier曲面的PIA法。

預(yù)處理技術(shù)的引入雖然能通過改善配置矩陣的譜分布提高迭代法的收斂速度，但在迭代過程中需額外求解預(yù)處理子的逆矩陣，增加了計(jì)算量，影響數(shù)據(jù)插值的計(jì)算效率。為進(jìn)一步降低預(yù)處理技術(shù)的計(jì)算量，文獻(xiàn)［18-19］用經(jīng)典的共軛梯度法近似求解預(yù)處理方程的法方程，即PPIA的非精確求解方法（inexact preconditioned progressive iterative approximation，IPPIA），但在求解過程中需在預(yù)處理方程兩邊同時(shí)乘以一個(gè)預(yù)處理矩陣的轉(zhuǎn)置矩陣，顯然，該法方程的條件數(shù)將隨之增大，影響求解預(yù)處理方程算法的計(jì)算效率、計(jì)算精度和穩(wěn)定性。如果能找到一種合適的算法直接求解預(yù)處理方程而不是其法方程，則可提升計(jì)算效率。但由于預(yù)處理方程的系數(shù)矩陣是非對(duì)稱的，無法用共軛梯度法求解。文獻(xiàn)［7］利用廣義極小殘差法（generalized minimal residual method，GMRES）求解配置方程，得到插值曲線，數(shù)值實(shí)例表明，當(dāng)插值點(diǎn)較少時(shí)，GMRES 法的計(jì)算效率很高，但當(dāng)插值點(diǎn)較多時(shí)，得到的插值曲線精度不理想［18］。欣慰的是，非精確求解方法對(duì)精度的要求并不高，因此，本文將利用GMRES 法求解預(yù)處理方程，以進(jìn)一步提升算法的計(jì)算效率和計(jì)算精度。

基于上述分析，本文在文獻(xiàn)［18-19］基礎(chǔ)上，研究張量積型Said-Ball 曲面的PPIA 法及其非精確求解方法，以期進(jìn)一步提高數(shù)據(jù)插值的計(jì)算效率。

1 相關(guān)工作

1.1 張量積型Said-Ball 曲面

Said-Ball 基函數(shù)及相應(yīng)的張量積型Said-Ball曲面片的定義［20］為：

1.2 漸近迭代逼近法

則稱S(0)(u，v)具有PIA 性質(zhì)［11］。文獻(xiàn)［11］證明了所有全正基函數(shù)曲線（曲面）均具PIA 性質(zhì)。文獻(xiàn)［20］稱Said-Ball 基函數(shù)為全正基函數(shù)，從而張量積型Said-Ball 曲面也具有PIA 性質(zhì)。

張量積型Said-Ball 曲面的PIA 性質(zhì)意味著式（3）收斂于插值曲面S?(u，v)的控制頂點(diǎn)，即

文獻(xiàn)［6-7］指出，PIA 法中式（3）也可看作求解式（4）的Richardson 迭代法。

最后，給出矩陣的Kronecker 積的性質(zhì)［21］。

引理1關(guān)于矩陣的Kronecker積，有：

（a）A?B可逆當(dāng)且僅當(dāng)A和B均可逆，且有(A?B)?1=A?1?B?1；

（b）若A，B，C和D為階數(shù)滿足AC和BD乘法的矩陣，則(A?B)(C?D)=(AC)?(BD)。

2 預(yù)處理漸近迭代逼近法

文獻(xiàn)［18］利用對(duì)角補(bǔ)償約化技術(shù)構(gòu)造了預(yù)處理子，該預(yù)處理子能在很大程度上改善全正基函數(shù)（包括Bernstein 基函數(shù)和Said-Ball 基函數(shù)）所生成的配置矩陣的譜分布。文獻(xiàn)［19］利用上述預(yù)處理子給出了張量積型Bézier 曲面的PPIA法，并取得了理想的加速效果。本文將該預(yù)處理技術(shù)應(yīng)用于張量積型Said-Ball 曲面，得到張量積型Said-Ball 曲面的PPIA 法。

2.1 預(yù)處理子M1和M2的構(gòu)造

首先，將配置矩陣B1分裂為

注1設(shè)M1和M2分別為對(duì)角補(bǔ)償約化技術(shù)所定義的關(guān)于B1和B2的預(yù)處理矩陣，相應(yīng)的半帶寬分別為q1和q2，則有［18］：

（a）預(yù)處理子M1和M2均為可逆矩陣。

（b）對(duì)于i=1，2，經(jīng)預(yù)處理后的矩陣的譜隨半帶寬qi的增大向1 附近聚集。

（c）綜合考慮計(jì)算量和預(yù)處理效果，當(dāng)預(yù)處理子的半帶寬約為配置矩陣階數(shù)的一半時(shí)效果較好。

2.2 預(yù)處理方法

在得到預(yù)處理子M1和M2后，由引理1，對(duì)式（4）進(jìn)行預(yù)處理，有

注2文獻(xiàn)［19］指出，在構(gòu)造張量積型曲面的PPIA 法時(shí)，用矩陣的Kronecker積的性質(zhì)分別對(duì)配置矩陣B1和B2進(jìn)行預(yù)處理，具有充分發(fā)揮預(yù)處理子的預(yù)處理效果和提高算法的穩(wěn)定性和魯棒性的作用。

3 PPIA 法的非精確求解方法

由式（6），在PPIA 法的迭代過程中需計(jì)算

與式（3）相比，這顯然增加了額外計(jì)算量。因矩陣的求逆運(yùn)算計(jì)算量較大，且具有不穩(wěn)定性。當(dāng)n和m均較小時(shí)，式（7）很容易計(jì)算，但當(dāng)n和m較大時(shí)，式（7）的計(jì)算復(fù)雜度增大，穩(wěn)定性也隨之降低，從而影響PPIA 法的計(jì)算效率。

注意到，計(jì)算式（7）等價(jià)于求解線性方程組

考慮用迭代法求解式（8），會(huì)得到一個(gè)嵌套迭代格式。將式（6）看作外迭代格式、將求解式（8）的迭代法看作內(nèi)迭代格式，若對(duì)內(nèi)迭代的精度要求較高，則求解式（8）需花費(fèi)較多的計(jì)算時(shí)間。在不影響外迭代收斂性的前提下，降低內(nèi)迭代的精度，可減少內(nèi)迭代的計(jì)算時(shí)間，從而提高計(jì)算效率。因此，為減少計(jì)算量和提高算法的穩(wěn)定性，考慮用GMRES 法求解式（8），得到的近似解。

其中，τk為給定的迭代終止誤差閾值，‖?‖F(xiàn)為Frobenius 范數(shù)。

于是可得到下一步迭代的控制頂點(diǎn)

稱式（11）為張量積型Said-Ball 曲面預(yù)處理漸近迭代逼近法的非精確公式，記為IPPIA-GMRES。

為提高計(jì)算效率，當(dāng)

時(shí)，式（11）終止迭代，其中ε(k)為IPPIA-GMRES 法迭代k次所生成的近似插值Said-Ball 曲面的插值誤差，

現(xiàn)討論張量積型Said-Ball 曲面IPPIAGMRES 法的收斂性。為方便，下文的定理及證明中均將M1?M2記為M，將B1?B2記為B。

定理1設(shè)M1和M2分別為對(duì)角補(bǔ)償約化技術(shù)所定義的關(guān)于B1和B2的預(yù)處理矩陣，P?為Said-Ball 插值曲面的控制頂點(diǎn)。對(duì)于每個(gè)外迭代步k，利用GMRES 法求解式（9），得到式（7）的近似解X(l，k)，則對(duì)IPPIA-GMRES 法迭代所生成的序列，有

定理1 得證。

由定理1，對(duì)于給定的預(yù)處理子M1和M2，IPPIA-GMRES 法的收斂性主要由τk決定。顯然當(dāng)τk=0時(shí)，利用GMRES 法可得到式（7）的精確解，此時(shí)IPPIA-GMRES 法退化為PPIA 法。因此，為使迭代法具有更好的收斂效果，τk越小越好。但τk越小，意味著對(duì)內(nèi)迭代精度的要求增高，內(nèi)迭代次數(shù)增加，計(jì)算量增大。在實(shí)際計(jì)算中，可在不影響外迭代收斂性的前提下降低對(duì)內(nèi)迭代計(jì)算精度的要求（無需τk=0），以提高計(jì)算效率。

4 數(shù)值實(shí)例

文獻(xiàn)［12］指出，在最優(yōu)權(quán)系數(shù)下，WPIA 法比PIA 法具有更快的收斂速度，因此本文將與最優(yōu)權(quán)系數(shù)下的WPIA 法進(jìn)行對(duì)比。文獻(xiàn)［7］利用GMRES 法求解曲線情形的配置方程，得到插值曲線，本文也將利用非重啟的GMRES 法求解曲面情形的配置方程，得到插值曲面。由于GMRES 法的收斂性取決于系數(shù)矩陣的譜分布和Ritz 值等，而本文提出的預(yù)處理子能改善配置矩陣的譜分布，因此本文的預(yù)處理技術(shù)也可用于GMRES法，即預(yù)處理 GMRES（preconditioned GMRES，PGMRES）法。

方便起見，在數(shù)值實(shí)例中，記文獻(xiàn)［19］的方法為IPPIA-CG。所有實(shí)驗(yàn)均在Intel（R）Core（TM）i7-8565U CPU @ 1.80 GHz 1.99 GHz、內(nèi)存為8 GB的筆記本電腦上實(shí)現(xiàn)。

例1 給定16個(gè)待插值的數(shù)據(jù)點(diǎn)：（1，1，1），（2，1，2），（3，1，2），（4，1，1），（1，2，2），（2，2，2.5），（3，2，2.5），（4，2，2），（2，3，2），（1，3，2），（2，3，2.5），（4，3，2），（1，4，1），（2，4，2），（3，4，2），（4，4，1）。

例2 給定20個(gè)數(shù)據(jù)點(diǎn)：（1，1，1），（1，2，4），（1，3，2），（1，4，2），（1，5，2），（2，1，2），（2，2，2），（2，3，3），（2，4，6），（2，5，4），（3，1，3），（3，2，2），（3，3，4），（3，4，4），（3，5，3），（4，1，4），（4，2，6），（4，3，1），（4，4，1），（4，5，2）。

例3 逼近從函數(shù)

上采樣得到11×12個(gè)數(shù)據(jù)點(diǎn)：

例4逼近從函數(shù)

上采樣得到42×54個(gè)數(shù)據(jù)點(diǎn)：

實(shí)驗(yàn)參數(shù)選取如下：由注1，在利用PGMRES、PPIA 和IPPIA-GMRES 法逼近例1～例4 中數(shù)據(jù)點(diǎn)時(shí)，預(yù)處理子的半帶寬取值如表1 所示。在測試IPPIA-GMRES 法時(shí)，取式（10）中的τk=0.01，取式（12）中的θ=0.85。為方便比較，取IPPIA-CG 法中的τk=0.01，θ=0.85。

表1 例1～例4 中預(yù)處理子的半帶寬Table 1 Half-bandwidth of preconditioners in examples 1 to 4

在逼近例1～例3 數(shù)據(jù)點(diǎn)時(shí)，PIA、WPIA 和PPIA 法迭代矩陣的譜分布如圖1 所示，相應(yīng)的迭代矩陣的譜半徑如表2 所示?？芍?，與PIA 法和WPIA 法相比，PPIA 法迭代矩陣的大部分特征值向0 方向聚集。PPIA 法迭代矩陣的譜半徑明顯小于PIA 法和WPIA 法。因?yàn)榈仃嚨淖V半徑越小，收斂速度越快，所以PPIA 法的收斂速度提升很快。此外，由于迭代矩陣的譜分布更聚集，P-GMRES 法的收斂效果也得到了改善。

圖1 例1～例3 中PIA、WPIA 和PPIA 法迭代矩陣的譜分布Fig.1 Spectral distributions of iteration matrices of PIA，WPIA and PPIA in examples 1 to 3

表2 例1～例3 中迭代矩陣的譜半徑Table 2 Spectral radius of iteration matrices in examples 1 to 3

表3～表5 分別為用不同方法迭代逼近例1～例3 中數(shù)據(jù)點(diǎn)時(shí)的插值誤差ε和運(yùn)行時(shí)間T，可知，WPIA 法的收斂速度較慢，PPIA 法的收斂速度較快，特別是例1 和例2，用PPIA 法迭代一次就能達(dá)到機(jī)器精度。同樣，在數(shù)據(jù)點(diǎn)很少時(shí)，用GMRES法和P-GMRES 法均能得到精確的插值曲面；隨著數(shù)據(jù)點(diǎn)的增多，2 種方法的插值誤差隨之減小，PGMRES 法比GMRES 法減小更明顯，說明本文提出的預(yù)處理子能改善GMRES 法的收斂性。在插值誤差幾乎相同的情況下，PPIA 法在迭代次數(shù)k和運(yùn)行時(shí)間上具有優(yōu)勢，說明本文方法具有計(jì)算效率優(yōu)勢。

表3 例1 中WPIA、PPIA、GMRES 和P-GMRES 法的插值誤差和運(yùn)行時(shí)間Table 3 Interpolation errors and runtime of WPIA，PPIA，GMRES and P-GMRES in example 1

表4 例2 中WPIA、PPIA、GMRES 和P-GMRES 法的插值誤差和運(yùn)行時(shí)間Table 4 Interpolation errors and runtime of WPIA，PPIA，GMRES and P-GMRES in example 2

表5 例3 中WPIA、PPIA、GMRES 和P-GMRES 法的插值誤差和運(yùn)行時(shí)間Table 5 Interpolation errors and runtime of WPIA，PPIA，GMRES and P-GMRES in example 3

在逼近例4 中數(shù)據(jù)點(diǎn)時(shí)，由于數(shù)據(jù)點(diǎn)較多，WPIA 法的收斂速度很慢，PPIA 法得到的插值誤差較大，這由預(yù)處理子的求逆運(yùn)算不穩(wěn)定造成。此時(shí)GMRES 法和P-GMRES 法的運(yùn)算結(jié)果亦不理想，可考慮采用非精確求解方法近似求解預(yù)處理方程。由表6的數(shù)值結(jié)果知，在終止閾值θ相同時(shí)，IPPIA-CG 法只需迭代1 次（此時(shí)，可通過提高內(nèi)迭代精度增加迭代次數(shù)，以減少IPPIA-CG 法的插值誤差），IPPIA-GMRES 法需迭代6次，但插值精度也相應(yīng)提高，說明IPPIA-GMRES 法能更有效地插值較大規(guī)模的數(shù)據(jù)點(diǎn)。在運(yùn)行時(shí)間上，IPPIAGMRES 法雖不及GMRES法和P-GMRES法，但明顯優(yōu)于WPIA法，且IPPIA-GMRES 法的插值誤差更小。

表6 例4 中IPPIA-GMRES 法和其他方法的插值誤差和運(yùn)行時(shí)間Table 6 Interpolation errors and runtime of IPPIA-GMRES and the other methods in Example 4

用WPIA 法和PPIA 法迭代逼近例1～例3的數(shù)據(jù)點(diǎn)，得到的Said-Ball 曲面如圖2～圖4 所示；用WPIA、IPPIA-CG 和IPPIA-GMRES 法迭代逼近例4 中的數(shù)據(jù)點(diǎn)，得到的Said-Ball 曲面如圖5 所示。不難看出，由本文方法生成的張量積型Said-Ball 曲面能快速逼近給定數(shù)據(jù)點(diǎn)。

圖2 例1 中的Said-Ball 插值曲面Fig.2 Said-Ball interpolation patches in example 1

圖3 例2 中的Said-Ball 插值曲面Fig.3 Said-Ball interpolation patches in example 2

圖5 例4 中的Said-Ball 插值曲面Fig.5 Said-Ball interpolation patches in example 4

5 結(jié)論

基于對(duì)角補(bǔ)償約化技術(shù)，探討了求解張量積型Said-Ball 曲面的PPIA法，并用GMRES 法求解預(yù)處理方程，得到非精確解，即IPPIA-GMRES 法。相較文獻(xiàn)［18］利用共軛梯度法求解預(yù)處理方程的方法，本文的非精確求解方法計(jì)算量更小，得到的插值曲面精度更高，計(jì)算效率也更高。此外，由于本文提出的預(yù)處理方法能使配置矩陣的譜分布更集中，在處理規(guī)模適當(dāng)?shù)牟逯祮栴}時(shí)，采用對(duì)角補(bǔ)償預(yù)處理子能在一定程度上改善GMRES 法的收斂性。

由于預(yù)處理技術(shù)能在有效改善非分段（片）基函數(shù)的同時(shí)配置矩陣的譜分布［19，22］，因此本文方法能推廣應(yīng)用于類似曲面（如張量積型Bézier 曲面）PIA法的加速問題。但本文方法也存在局限性，由于分片曲面（如雙三次B 樣條曲面、NURBS 曲面等）PIA法中的配置矩陣為稀疏矩陣，利用對(duì)角補(bǔ)償約化技術(shù)構(gòu)造的預(yù)處理子效果不明顯，甚至不起作用。不同于分片插值，本文方法利用一張Said-Ball 曲面對(duì)所有數(shù)據(jù)點(diǎn)插值，具有更高的連續(xù)性。未來，將進(jìn)一步探討分片曲面PIA 法的加速問題。