唐 敏,鄧國強
桂林電子科技大學 數(shù)學與計算科學學院 廣西密碼學與信息安全重點實驗室,廣西 桂林 541004
插值是函數(shù)逼近的重要方法,利用它可通過函數(shù)在有限個點處的取值情況,估算出函數(shù)在其他點處的近似值。在信號處理[1]、不確定性量化[2-4]、壓縮感知[5]、圖像處理[6-7]等領(lǐng)域需要有效地利用較少的函數(shù)值恢復一個函數(shù),而這個函數(shù)往往具有稀疏的表示形式,稀疏插值能有效地處理此類問題。稀疏多元多項式插值是一種降低計算機代數(shù)算法時間復雜度的有效方法,是稀疏有理函數(shù)插值[8]及稀疏隱函數(shù)插值[9]的基本步驟,應(yīng)用在結(jié)式計算[9]、多元多項式GCD計算[10-12]、方程求根、數(shù)值微分、數(shù)值積分、曲面的外形設(shè)計和有限元法及上述提及的諸多應(yīng)用領(lǐng)域。因此提高稀疏多元多項式插值的效率十分有意義。
對于多元多項式,稠密插值所需的插值點個數(shù)為(d+1)n,其中d是變元的次數(shù),n是變元的個數(shù)。研究工作表明,多元多項式插值算法的計算復雜度與變元個數(shù)可以不呈指數(shù)關(guān)系,而是與目標多項式的稀疏性相關(guān),即是關(guān)于n、d、t和lbp的多項式函數(shù),其中t是目標多項式的項數(shù),素數(shù)p是有限域的特征。
1979年Zippel提出了第一個具有多項式時間復雜度的稀疏多元多項式插值算法[12]。該算法基于這樣的假設(shè):如果在一個隨機的賦值點處多項式的取值為零,那么它就是一個零多項式。由于這個結(jié)論的成立是高概率的,因此Zippel算法是一個概率性的方法。在執(zhí)行上,Zippel算法順序地對每個變元逐一插值,不能并行化,插值點個數(shù)為ndt,時間復雜度為O(ndt3)。1990年,Zippel改進了他的方法[13],使用了形如的插值點,原算法中需要計算的一般線性系統(tǒng)變成了轉(zhuǎn)置的Vandermonde系統(tǒng),求解時間從O(t3)降到O(t2)。由于Zippel算法是一種較好的具有多項式時間復雜度的方法,大部分計算機代數(shù)系統(tǒng)都內(nèi)置了基于Zippel算法的多元GCD計算的實現(xiàn),比如Macsyma、Magma和Mathematica。
1988年,Ben-Or和Tiwari提出了一個確定性的算法[14],可以對以整數(shù)、有理數(shù)、實數(shù)或復數(shù)為系數(shù)的多元多項式進行稀疏插值。給定項數(shù)為t的多項式f的項數(shù)界T≥t,算法需要2t個插值點,即前n個素數(shù) 的 冪 次,0≤i≤2t-1。該方法與Zippel算法不同,不是逐個變元的插值,而是一次性地對f的各個單項式進行插值,而且可以并行化。主要的不足是插值點的存儲空間為tlbn位,當多項式規(guī)模較大時,由于中間表達式膨脹而導致計算速度非常慢。
2000年,Kaltofen等設(shè)計了一個Zippel和Tiwari方法相結(jié)合的混合算法,稱之為“racing算法”[15-16]。為了減少插值點個數(shù),在Zippel算法中插值下一個變元的時候,同時執(zhí)行Newton插值和單變元的Ben-Or/Tiwari插值,取其快者。
2010年,Javadi和Monagan提出了一個有限域上的稀疏多元多項式插值算法[17],它也是一個概率算法。該方法是Ben-Or/Tiwari算法的一個改進,使用O(t)個插值點獨立完成每個變元的插值,共需O(nt)個插值點,并且可以并行化。為了準確計算每個變元在各個單項式中的準確次數(shù),使用了二部圖完美匹配進行測試和驗證。該算法和Zippel算法及Kaltofen的racing算法進行了比較,在運行時間和插值點個數(shù)等方面具有一定的優(yōu)勢。
2011年,Cuyt和Lee給出了稀疏有理函數(shù)插值算法[8],其中涉及到單變元有理函數(shù)插值和基于Zippel算法及Tiwari算法的稀疏多元多項式插值。2016年,基于Cuyt的稀疏有理函數(shù)插值算法,Tang給出了稀疏隱函數(shù)插值算法[9]。在稀疏有理函數(shù)插值和隱函數(shù)插值算法中,稀疏多元多項式插值作為基本的子函數(shù),對整個算法的效率起到關(guān)鍵作用。
本文提出了一個有限域上的稀疏多元多項式插值算法,它是Javadi/Monagan算法的改進,改進方案主要有兩點:一是消除了必須給定項數(shù)界T的限制,通過計算預先設(shè)計好的特定矩陣的行列式,得到f的準確項數(shù);二是消除了必須給定次數(shù)界D的限制,通過構(gòu)造輔助函數(shù),利用概率法結(jié)合提前終止技術(shù)的Cauchy插值法,計算出f的準確次數(shù),解決了Javadi和Monagan論文中提到的次數(shù)界D過高而導致高計算復雜度的問題。更進一步,改進算法無需給定T和D,對于實際問題在利用插值恢復或近似時更具實用性。
本文后續(xù)組織結(jié)構(gòu)如下:第2章介紹了Javadi/Monagan算法的基本思想;第3章給出了改進的Javadi/Monagan算法;第4章分析了改進算法的時間復雜度;第5章中設(shè)計了三組數(shù)值實驗,對改進算法和Javadi/Monagan算法進行了比較,給出了實驗結(jié)果;第6章給出了一個應(yīng)用實例。
本章給出稀疏多元多項式插值問題的形式化定義及Javadi/Monagan算法的思想。
令p是一個素數(shù),f∈Zp[x1,x2,…,xn]是一個用黑盒表示的多元多項式:
稀疏多元多項式插值的目標是用盡可能少的插值點及較低的多項式時間復雜度算法還原多項式f。
Javadi/Monagan算法可分成兩個步驟,首先確定各個單項式,然后確定系數(shù)。在第一個步驟中,逐一計算每個變元在各個單項式中的次數(shù)。第二個步驟需要求解一個Vandermonde線性代數(shù)系統(tǒng)。同時要求給定f的項數(shù)界T≥t及次數(shù)界D≥degf。
(1)確定單項式Mi,i=1,2,…,t。
③使用Berlekamp/Massey算法[18]生成n+1個關(guān)于z的單變元多項式,使得,其中i=0,1,…,2t-1,j=0,1,…,n。
④確定變元xk,k=1,2,…,n在f中的單項式Mi,i=1,2,…,t中的次數(shù)eij。詳見2.3節(jié)。
⑤由④可確定eij,i=1,2,…,t,j=1,2,…,n,則單項式
(2)計算系數(shù)。
通過求解線性方程組計算系數(shù)ai,i=1,2,…,t。令r1,r2,…,rt是Λ1(z)的根,且ri=Mi(α1,α2,…,αn),回顧是輸入為時黑盒的輸出,因此有:
此為Vandermonde系統(tǒng),有唯一解。
令R1={r1,r2,…,rt}和分別是Λ1和Λk+1的全部根構(gòu)成的集合。令Dj包含單項式Mj中xk的所有可能的次數(shù),即:
二部圖Gk定義如下:U和V是二部圖Gk中頂點個數(shù)為t的互補頂點子集,U和V中的結(jié)點表示R1和Rk中的元素,即ui∈U,vj∈V,分別用ri和標記,ui和vj之間有一條權(quán)值為dij的邊當且僅當
引理1[17]變元xk在所有單項式中的次數(shù)可唯一確定當且僅當二部圖Gk有唯一的完美匹配。
以下給出為計算xk的次數(shù),構(gòu)造的二部圖Gk沒有唯一的完美匹配時的解決方案。
隨機選擇不同于α1,α2,…,αn+1的元素令,使用 Berlekamp/Massey算法,以v′i=f(β′i),0≤i≤ 2t-1為輸入,生成多項式Λ′k+1(z)。令是Λ′k+1(z)的根構(gòu)成的集合,G′k是通過Λ1(z)和Λ′k+1(z)構(gòu)造的二部圖。
定義1[17]二部圖G′k與Gk的交集定義如下:與G′k具有相同的結(jié)點,ri和之間有權(quán)值為dij的邊當且僅當在二部圖Gk中,ri連接且在二部圖G′k中,ri連接權(quán)值均為dij。
引理2[17]令eij=degxj(Mi),在二部圖中,結(jié)點ri和之間有邊相連,且權(quán)值為eij。
定理1[17]令的交集有唯一的完美匹配。
引理1、引理2和定理1的證明參考文獻[17]。由定理1,如果Gk沒有唯一的完美匹配,可通過增加2t個插值點,構(gòu)造二部圖來確定變元xk在單項式Mj(1≤j≤t)中的次數(shù)。
本章給出對Javadi/Monagan算法進行改進的兩個主要策略。3.1節(jié)和3.2節(jié)分別描述了如何只利用黑盒的輸入輸出(即有限個插值點),計算多項式f的準確項數(shù)t和準確次數(shù)d,從而消除了原算法必須預先給定這兩個指標的上界的限制,同時提高了實用性及降低了計算時間復雜度(分析過程見第4章)。
除了Javadi/Monagan算法,許多插值算法都要求給定目標多項式f的項數(shù)界T,本文給出一個準確計算多項式f的項數(shù)t的方法。該方法由定理2保證。
定理2[14]令V是元素為(V)ij=vi+j-2的矩陣。Vl是由V的前l(fā)行前l(fā)列組成的方陣。如果t是多項式f的準確項數(shù),即f中非零單項式的個數(shù),那么:
根據(jù)定理2,通過下述過程可計算f的準確項數(shù):
例如:黑盒多元多項式f定義為:令p=1 009,隨機生成中的3個整數(shù),比如α1=66,α2=12,α3=3。
因此多項式f的準確項數(shù)為5。
在2.3節(jié)構(gòu)造二部圖的過程中,計算單項式Mj中xk的所有可能的次數(shù)Dj時,需要檢測的次數(shù)為D+1。因此Javadi/Monagan算法在給定的次數(shù)界過高時,會導致高計算復雜度。本節(jié)利用單變元有理函數(shù)插值計算f的準確次數(shù)d,使得在確定每個變元的次數(shù)時,時間復雜度達到最低。方法描述如下:
(1)構(gòu)造輔助有理函數(shù)H。
通過在f中引入齊次變元z,構(gòu)造f為分子,g為分母的輔助有理函數(shù)H。其中f(x1,x2,…,xn)為黑盒多項式,g為隨機生成的關(guān)于變元z的一次多項式。若視F為關(guān)于變元z的單變元多項式,則系數(shù)為關(guān)于x1,x2,…,xn的多元多項式,注意到系數(shù)多項式與z是齊次的。
易證f的全次數(shù)與F中z的最高次相同。
(2)對H進行單變元有理函數(shù)插值。
任取n元組 (α1,α2,…,αn),對H進行關(guān)于變元z的單變元有理函數(shù)插值,可得到:
(3)H(α1,α2,…,αn,z)的分子f(α1z,α2z,…,αnz)關(guān)于z的最高次即為黑盒多項式f的全次數(shù)。
假設(shè)f(α1z,α2z,…,αnz)和g(z)有非平凡的公因子,即g(z)|f(α1z,α2z,…,αnz),也就是說,g(z)的根也是f(α1z,α2z,…,αnz)的根。因為g(z)是隨機生成的,而f(α1z,α2z,…,αnz)的根是有窮的,所以f(α1z,α2z,…,αnz)和g(z)高概率互素,從而f(α1z,α2z,…,αnz)關(guān)于z的次數(shù)即為黑盒多項式f的全次數(shù)。
令 (α1,α2)=(1,1),g=5z+3 ,則:
利用單變元有理函數(shù)插值可得到H(1,1,z)=(z4+2z2)/(5z+3),則f的準確次數(shù)為H的分子中z的最高次,即為4。
接下來給出單變元有理函數(shù)插值的實現(xiàn)過程。對于單變元有理函數(shù)f/g∈Κ(Z),應(yīng)用文獻[19]提出的概率方法,并結(jié)合提前終止技術(shù)的Cauchy插值來恢復f/g。來源于文獻[19-20]的引理給出了單變元有理函數(shù)插值的理論基礎(chǔ)。
引理3[20]令Κ是任意域,F(xiàn)(Z),G(Z),H(Z)∈Κ(Z)且 gcd(F,G)=1。令是非負整數(shù),且1。令ik,f/g∈Κ(Z)不必是Κ中互不相同的元素,滿足:
對于l≥ 2 ,令hl(Z),ql(Z)∈Κ(Z)分別是Euclidean多項式余項序列中的第l次余項和商。
對于l≥ 2 ,令wl(Z),gl(Z)∈Κ(Z)是擴展Euclidean方法中的乘子wlh0+glh1=hl,即:
那么存在下標j≥ 1滿足,并且對下標j有:
根據(jù)引理3,如果給定單變元有理函數(shù)的次數(shù)上界,使用擴展Euclidean算法可計算有理函數(shù)的分子和分母。在文獻[19]中,提出了基于引理3的單變元有理函數(shù)插值算法。為確定準確的f/g的次數(shù),可對k從1迭代到d+e+1,其中d、e分別是f和g的次數(shù)。
基于概率法和提前終止技術(shù)的Cauchy插值法恢復單變元有理函數(shù),算法描述如下:
算法1單變元有理函數(shù)插值算法
其 中 Interpolate(i1,i2,…,ik;H(i1),H(i2),…,H(ik)) 表示給定k個點 (i1,H(i1)),(i2,H(i2)),…,(ik,H(ik))的自變量和函數(shù)值,利用Newton插值或Lagrange插值得到次數(shù)為k-1的單變元多項式;Rem表示余式,Quo表示商。
黑盒多項式f的次數(shù)即為deg(h1)。
本節(jié)給出改進的Javadi/Monagan算法的偽代碼。
算法2改進的Javadi/Monagan算法(有限域上稀疏多元多項式插值算法)
以下給出改進的稀疏多元多項式插值算法的一個實例。黑盒多元多項式f定義為:
令p=101,變元個數(shù)n=3。重構(gòu)f的過程如下:
(2)計算出多項式f的準確項數(shù)t=5。
(4)計算Λk(z),k=1,2,…,4:
令vi=B(βi),i=0,1,…,2t-1,當k>1時,用替換使用Berlekamp/Massey算法生成由序列v0,v1,…,v2t-1得到的Λk∈Zp(z):
(5)計算出多項式f的次數(shù)d=5。
(6)確定degxk(Mi),1≤i≤t,k=1,2,…,n。
以第1個變元x1為例,計算Λ1的根:
計算Λ2的根:
測試rj×(αn+1/α1)i,0≤i≤d,j=1,2,…,t,得到:
構(gòu)造二部圖G1,如圖1所示。
Fig.1 Bipartite graphG1圖1 二部圖G1
因為二部圖G1沒有唯一的完美匹配,所以生成α5=54。
計算Λ5(z)=z5+100z4+73z3+85z2+51z+94及Λ5的根。構(gòu)造二部圖G1′,如圖2所示。
Fig.2 Bipartite graphG1′圖2 二部圖G1′
構(gòu)造Gk和Gk′的交集,如圖3所示。因此,x1在Mi,1≤i≤t中的次數(shù)分別為0,0,2,2,0。同理,計算x2和x3在Mi中的次數(shù),分別為0,0,1,2,1和0,5,1,1,2。則
Fig.3 Bipartite graph圖3 二部圖
本章討論有限域上改進的稀疏多元多項式插值算法(算法2)的時間復雜度。需要考慮的主要因素有:
(1)計算多項式f的項數(shù)t,時間復雜度為O(t2)。
(2)計算多項式f的次數(shù)d,時間復雜度為O(d)。
(3)調(diào)用Berlekamp/Massey算法計算Λk至多n+2次,每次調(diào)用的復雜度為O(t2)。
(4)使用文獻[13]中的技術(shù)求解Vandermonde系統(tǒng),時間復雜度為O(t2)。
(5)求解Λk的根,時間復雜度為O(t2lbp)。
(6)計算變元xk的次數(shù),時間復雜度為O(tlbt+dtlbt)。
(7)構(gòu)造二部圖Gk,時間復雜度為O(dt2),求Gk和G′k的交集,時間復雜度為O(dtlbd)。
綜上所述,改進算法的時間復雜度為:
從時間復雜度分析可以看出,如果給定的項數(shù)界T和次數(shù)界D過高,會導致高計算復雜度,而改進算法花費較少的代價計算出多項式f的準確項數(shù)t及次數(shù)d,達到了此類型算法在后續(xù)操作上的最低時間復雜度。
本章對改進算法和Javadi/Monagan算法進行性能比較。兩種算法的編程環(huán)境均為Maple15,程序運行的硬件環(huán)境為Intel?CoreTMi7 2.20 GHz處理器和4.00 GB內(nèi)存,操作系統(tǒng)為Windows 7。注意兩個程序都是順序執(zhí)行,在確定變元x1,x2,…,xn的次數(shù)時未并行化。
本章給出了3組測試集的性能比較結(jié)果,使用的多項式都是隨機生成的,比較的對象為CPU時間。黑盒中的多元多項式系數(shù)取自Zp,其中p=100 003。
測試集1本組測試集為6個包含3個變元的多元多項式。第i個多項式(1≤i≤6)使用如下的Maple命令隨機生成:
其中,第i個多項式包含2i個非零項,d=20是多項式的準確次數(shù)。該組測試分別執(zhí)行了改進算法和Javadi/Monagan算法,記錄每個多項式在兩種算法下的運行時間(單位:s),如表1所示。表頭第1列Ex.表示例子的編號;第2列t表示項數(shù);第3列括號中的百分比表示改進算法計算項數(shù)和次數(shù)的總時間與整個算法的運行時間的比值;改進算法無需給定項數(shù)界和次數(shù)界,Javadi/Monagan算法執(zhí)行時,分別給定的次數(shù)界為D=20,D=30和D=30,項數(shù)界為多項式的準確項數(shù)T=2i(1≤i≤6)。表2和表3的設(shè)置同表1。
Table 1 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=3)表1 改進算法與Javadi/Monagan算法的比較(n=3)
從表1可以看出,隨著i的增加,兩種算法的執(zhí)行時間也隨之增加。在給定準確項數(shù)界T=2i和準確次數(shù)界D=20的情況下,Javadi/Monagan算法的執(zhí)行時間優(yōu)于改進算法,原因是改進算法利用3.1節(jié)和3.2節(jié)的方法分別計算出準確的項數(shù)和次數(shù),不需預先給定這兩個數(shù)據(jù)。顯然,這兩者的計算需要耗費一定的時間,但從表1中的第3列括號中的百分比可以看出,這兩者的計算時間占整體運行時間的比值隨著次數(shù)的增加越來越小,因而優(yōu)勢越來越明顯。另一方面,隨著次數(shù)界D的增加,Javadi/Monagan算法需要的執(zhí)行時間也隨之增加,當D較大時,算法的時間復雜度較高,這就是Javadi和Monagan在論文中提到的壞次數(shù)界問題(bad degree bound)。改進算法則有效地避免了這一問題。
測試集2本組測試集為6個包含6個變元的多元多項式。第i個多項式(1≤i≤6)使用如下的Maple命令隨機生成:
表2給出了測試集2下改進算法和Javadi/Monagan算法的執(zhí)行時間,表格的各項含義與測試集1相同。
測試集3本組測試集為6個包含12個變元的多元多項式。第i個多項式(1≤i≤6)使用如下的Maple命令隨機生成:
Table 2 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=6)表2 改進算法與Javadi/Monagan算法的比較(n=6)
表3給出了測試集3下改進算法和Javadi/Monagan算法的執(zhí)行時間,表格的各項含義與測試集1相同。
Table 3 Performance comparison of improved algorithm and Javadi/Monagan algorithm(n=12)表3 改進算法與Javadi/Monagan算法的比較(n=12)
從3組測試集可以看出,改進算法在確定插值多項式f的項數(shù)和次數(shù)上,花費的時間占整個算法運行時間的比例非常小,多項式規(guī)模越大,比例越小。如果Javadi/Monagan算法執(zhí)行時,給定的次數(shù)界D越高,運行時間越長,而改進算法不受此影響。
本章給出改進的Javadi/Monagan稀疏插值算法在幾何問題上的一個應(yīng)用。著名的Morley三等分定理描述如下:在任意三角形ABC中,對3個內(nèi)角∠CAB、∠ABC和∠ACB分別進行三等分,相鄰邊交于三點P、Q、R,則三角形PQR必為一個正三角形,如圖4所示。令a=BC,b=CA,c=AB,x=PQ=QR=RP,求x和a、b、c之間的關(guān)系。
Fig.4 Morley triangle圖4 Morley三角形
根據(jù)題設(shè)中的幾何關(guān)系和相關(guān)定理可知(推導過程略):
其中,R是三角形ABC的外接圓半徑,有如下方程組:
令x和a、b、c之間的關(guān)系為f(x,a,b,c)=0 ,如果利用結(jié)式等符號消元法,消去方程組(1)中的變元p、q、r,那么可得到x和a、b、c之間的關(guān)系f(x,a,b,c)。實際上,如果把a、b、c視為符號參數(shù),消元過程因中間表達式膨脹而無法完成,若將a、b、c用實數(shù)α1、α2、α3進行替換,消元過程可順利進行,結(jié)果是實例化的關(guān)系式f(x,α1,α2,α3),例如(假定模p=105):
如果視f(x,α1,α2,α3)為關(guān)于變元x的單變元多項式,則x的各次冪的系數(shù)是關(guān)于變元a、b、c的多元多項式,使用稀疏多元多項式插值恢復系數(shù)多項式,即可得到x和a、b、c之間的關(guān)系。在此例中,x的各個冪次的系數(shù)多項式的項數(shù)和次數(shù)都無法估計,如果采用原始的Javadi/Monagan算法,可能因為項數(shù)界和次數(shù)界設(shè)置不正確而導致無法得出結(jié)果或計算復雜度過高。若采用改進的Javadi/Monagan算法,可準確計算出各個系數(shù)多項式的項數(shù)和次數(shù),最終得到的x和a、b、c之間的關(guān)系式:
此例中,f(x,a,b,c)共有804項,其中x27的系數(shù)多項式的項數(shù)為45,次數(shù)為32。x25的系數(shù)多項式的項數(shù)為28,次數(shù)為34。其余系數(shù)多項式的指標不再贅述。
在結(jié)式消元、GCD計算、幾何組合優(yōu)化、信號處理等問題上,很多情況下無法預知目標多項式的次數(shù)和項數(shù),故改進算法針對這些問題更具實用性。
稀疏多元多項式插值是很多計算機代數(shù)算法的子函數(shù),也廣泛應(yīng)用于信號處理、壓縮感知、圖像處理等領(lǐng)域。在給定項數(shù)界T和次數(shù)界D較高的情況下,改進算法可有效降低Javadi/Monagan算法的計算復雜度。進一步,改進算法消除了傳統(tǒng)插值算法中必須預先給定插值多項式項數(shù)T和次數(shù)D的必要條件,更具實用性。
3.3節(jié)的算法2給出了有限域上改進的Javadi/Monagan稀疏多元多項式插值算法,其中兩個循環(huán)體(算法2中步驟4到步驟8及步驟10到步驟19)的計算任務(wù)可并行化處理,一是生成關(guān)于變元z的單變元多項式Λ1(z),Λ2(z),…,Λn+1(z),可分為相互獨立的n+1個子任務(wù),二是確定變元x1,x2,…,xn在各個單項式中的次數(shù),可分為相互獨立的n個子任務(wù)。假設(shè)有q個計算內(nèi)核的集群系統(tǒng),每個內(nèi)核可獨立承擔上述計算Λj(z),j=1,2,…,n+1 的 (n+1)/q個子任務(wù),以及計算變元xj,j=1,2,…,n在各個單項式Mi,i=1,2,…,t中的次數(shù)eij的n q個子任務(wù),計算結(jié)果可保存在一個文件中,最終可確定黑盒多項式f中的各個單項式。在后續(xù)工作中,擬將改進算法進行并行化處理,以進一步提高稀疏多元多項式插值算法的效率。