非線性方程組自反解的非精確Newton-MCG算法?

梁志艷,張凱院,寧倩芝

(西北工業大學應用數學系，西安 710072)

1 引言

考慮非線性代數方程組

引入記號

方程組(1)可寫為F(x)=g．

非線性代數方程組在科學計算和工程應用領域有著廣泛的應用，人們對其數值求解方法進行了大量的研究工作，迄今已經建立了一些有效的迭代算法．例如，針對大型稀疏非線性代數方程組，Bellavia和Morini[1]建立了一類具有全局收斂性的非精確Newton-GMRES算法；針對具有非Hermite正定Jacobi矩陣的非線性系統，Bai等[2,3]利用Hermite反Hermite分裂(HSS)迭代法求解線性代數方程組，建立了非精確Newton-HSS算法以及效率更高的非精確Newton-PSS算法；針對大型稀疏弱非線性代數方程組，Bai等[4,5]利用Picard迭代法求解非線性代數方程組，建立了Picard-HSS算法以及效率更高的Picard-MHSS算法．

設P為n階對稱正交矩陣，若列向量x∈Rn滿足Px=x，則稱x為關于矩陣P的自反向量，記全體關于矩陣P的自反向量的集合為特別地，當P為單位矩陣時，自反向量是一般向量；當P為次單位矩陣時，自反向量是中心對稱向量．自反向量在物理學中有廣泛的應用[6]，本文采用修正共軛梯度法(MCG算法)求由Newton算法每一步迭代計算導出的線性代數方程組的近似自反解或者近似自反最小二乘解(Ls解)，建立求方程組(1)關于矩陣P的自反解的非精確Newton-MCG算法(N-MCG算法)．MCG算法不要求導出的線性代數方程組的系數矩陣具備正定性、可逆性或者列滿秩性，且在理論上具備有限步收斂性，因此總是可行的．此外，MCG算法還具有自反解存在性的判斷功能，即MCG算法中斷時，線性代數方程組無自反解．因此，建立的非精確N-MCG算法僅要求方程組(1)有自反解，不要求它的自反解唯一．

2 求方程組(1)自反解的非精確Newton算法

任意給定n階對稱正交矩陣P，設x?是方程組(1)的關于矩陣P的自反解(簡稱為自反解，下同)，x(k)是x?的近似向量，由于

所以方程組(1)可線性化為

用線性方程組(2)的解x(k+1)作為x?的一個新的近似向量，引入記號

線性方程組(2)可寫為

求方程組(1)自反解的迭代格式為x(k+1)=x(k)+y(k)．

參照文獻[7]的算法原理，建立求方程組(1)自反解的Newton算法如下：

第1步給定初始向量計算b1，置k:=1；

第2步若bk=0，停止；否則，求線性方程組(3)的自反解y(k)；

第3步計算x(k+1)=x(k)+y(k)及bk+1，置k:=k+1，轉第2步．

需要指出：對于Newton算法中的某個k，當線性方程組(3)無自反解y(k)時，可用它的自反Ls解y(k)來代替．對于Newton算法有如下的收斂性結論[7]：假設x?是方程組(1)的單根，且初始向量x(1)充分接近x?，那么由Newton算法確定的向量序列{x(k)}收斂于x?．

一方面，Newton算法要求初始矩陣充分接近于方程組(1)的自反解，這是比較困難的．另一方面，在每一步迭代過程中都需要求線性方程組(3)的自反解，當x(k)不接近于方程組(1)的精確自反解時，線性方程組(3)未必有自反解；即使有自反解，也沒有必要求出線性方程組(3)的精確自反解．因此，Dambo等[8]提出了非精確Newton算法，其基本思想是在每一步Newton算法的迭代過程中只是求線性方程組(3)的某種近似解．所以當線性方程組(3)無自反解y(k)時，可以采用非精確Newton算法．參照Newton算法及文獻[8]的算法原理，建立求方程組(1)自反解的非精確Newton算法如下：

第1步給定初始向量計算b1，置k:=1；

第2步若bk=0，停止；否則，給定參數ηk∈[0,1)，求使得

第3步計算x(k+1)=x(k)+y(k)及bk+1，置k:=k+1，轉第2步．

式(4)中∥·∥表示向量的2-范數(下同)．在非精確Newton算法中，參數ηk∈[0,1)用來控制精確性水平．取ηk=0時，非精確Newton算法就是Newton算法．需要指出：對于非精確Newton算法中的某個k，若線性方程組(3)有自反解，則求其近似自反解y(k)使滿足式(4)；若線性方程組(3)無自反解，則求其近似自反Ls解y(k)使滿足式(4)．對于非精確Newton算法有以下的收斂性定理．

定理1[8]設ηk≤ηmax

3 求線性方程組(3)自反解與自反Ls解的MCG算法

下面建立求線性方程組(3)自反解與自反Ls解的MCG算法．考慮線性方程組(3)的一般形式Ay=b(略去記號的上下標)，研究以下兩個問題：

問題I當Ay=b有自反解時，求使滿足Ay=b．

問題II當Ay=b無自反解時，求使得∥Ay?b∥=min．

3.1 求解問題I的MCG算法

參照文獻[9]的算法原理，建立求解問題I的MCG算法(算法I)如下：

第1步任意給定初始向量置k:=1，計算

第2步如果rk=0，或者rk?=0而zk=0，停止計算；否則，計算

第3步計算

第4步置k:=k+1，轉第2步．

可以驗證，算法I中的向量滿足對于算法I有以下的收斂性定理(證明過程類似文獻[9])．

定理2設Ay=b有自反解，則對任意的初始向量算法I可在有限步迭代計算后求得問題I的一個解，即Ay=b的一個自反解．Ay=b無自反解的充要條件是存在正整數k，使得由算法I得到的rk?=0而zk=0．

3.2 求解問題II的MCG算法

在算法I中，當rk?=0而zk=0時，算法I中斷，這表明Ay=b無自反解，因此，需要求解問題II，即求Ay=b的自反Ls解．引入記號

定理3求Ay=b的自反Ls解等價于求它的約束正規方程組

的自反解，而約束正規方程組(5)一定有自反解．

證明 求Ay=b的自反Ls解等價于求使得

下面證明求極小值問題(6)的自反解等價于求方程組(5)的自反解．聯立方程組

可寫為Ny=f，求方程組Ny=f的自反Ls解等價于求極小值問題(6)的自反解．方程組Ny=f的正規方程組為NTNy=NTf，即方程組(5)．因此，求方程組Ny=f的自反Ls解等價于求其正規方程組NTNy=NTf的自反解，從而求極小值問題(6)的自反解就是求方程組(5)的自反解，反之亦然．所以求Ay=b的自反Ls解等價于求方程組(5)的自反解．

下面證明方程組(5)一定有自反解．因為正規方程組NTNy=NTf，即方程組(5)一定有解，可設ey是它的一個解(未必是自反解)，那么由于P為對稱正交矩陣，所以

因而容易驗證By?=Cb，故方程組(5)一定有自反解． 證畢

參照算法I，建立求約束正規方程組(5)的自反解，即求解問題II的MCG算法(算法II)如下：

第1步任意給定初始向量置k:=1，計算

第2步如果Rk=0，停止計算；否則，計算

第3步計算

第4步置k:=k+1，轉第2步．

可以驗證，算法II中的向量滿足對于算法II有以下的收斂性定理(證明過程類似文獻[9])．

定理4對任意的初始向量算法II可在有限步迭代計算后求得問題II的一個解，即約束正規方程組(5)的一個自反解．

4 非精確N-MCG算法的計算方案

參照文獻[7]的計算方案，設計求方程組(1)自反解的計算方案如下．

N-MCG算法：

第1步給定初始向量計算b1，置k:=1；

第2步若bk=0，停止；否則，采用算法I求使滿足線性方程組(3)；若算法I中斷(此時線性方程組(3)無自反解)，或者算法I的迭代次數達到預先給定的正整數N0，采用算法II求使得∥Aky(k)?bk∥=min；

第3步計算x(k+1)=x(k)+y(k)及bk+1，置k:=k+1，轉第2步．

非精確N-MCG算法：

第1步給定初始向量計算b1，置k:=1；

第2步若bk=0，停止；否則，給定參數ηk∈[0,1)，采用MCG算法求y(k)∈使得

其中ε是MCG算法的終止準則，并考慮以下途徑：

1) 采用算法I求線性方程組(3)的近似自反解y(k)，使滿足式(7)；

2) 若算法I中斷(此時線性方程組(3)無自反解)，或者算法I的迭代次數達到預先給定的正整數N0，采用算法II求線性方程組(3)的近似自反Ls解y(k)，使滿足式(7)，或者滿足算法II的終止條件；

第3步計算x(k+1)=x(k)+y(k)及bk+1，置k:=k+1，轉第2步．

5 數值算例

用k11及k12分別表示算法I迭代次數超限及MCG算法中斷的總次數，k0,k1及k2分別表示Newton算法、算法I及算法II的總迭代次數，t表示計算時間，n表示未知向量的維數．算例使用Matlab7.0軟件-CPU3.00GHz微機計算，Newton算法的終止準則為10?7，MCG算法的終止準則為10?8．沒有特別約定時MCG算法的初始向量為y(1)=0，算法I的迭代次數限制為N0=999．

例1采用N-MCG算法及非精確N-MCG算法求方程組(1)關于矩陣P的自反解，有關的向量和矩陣如下(部分數據取自文獻[10])：

此時方程組(1)有自反解x(0)．選取Newton算法的初始向量為x(1)．

1) 取

兩種算法的迭代次數和計算時間對比見表1和表2．

表1:例1(1)N-MCG算法的計算結果

表2:例1(1)非精確N-MCG算法的計算結果(n=1080)

由表1可見，采用算法I求線性方程組(3)的自反解時，出現迭代次數超限．由表2可見，ηk=0.1時非精確N-MCG算法的效率較高．對比表1和表2可見，非精確NMCG算法比N-MCG算法的效率高．兩種算法求出的自反解均為已知的x(0)．

2) 取

兩種算法的迭代次數和計算時間對比見表3和表4．

表3:例1(2)N-MCG算法的計算結果

表4:例1(2)非精確N-MCG算法的計算結果(n=1500)

由表3可見，由Newton算法每一步迭代計算導出的線性方程組(3)都有自反解．由表4可見，ηk取三個不同值時非精確N-MCG算法的效率相當．對比表3和表4可見，非精確N-MCG算法比N-MCG算法的效率高．兩種算法求出的自反解均為已知的x(0)．

例2采用N-MCG算法及非精確N-MCG算法求方程組(1)關于矩陣P的自反解，有關的向量和矩陣如下(部分數據取自文獻[10])：

此時方程組(1)有自反解x(0)．選取Newton算法的初始向量為x(1)，兩種算法的迭代次數和計算時間對比見表5和表6．

表5:例2N-MCG算法的計算結果

表6:例2非精確N-MCG算法的計算結果(n=999)

由表5可見，采用算法I求線性方程組(3)的自反解時，出現迭代次數超限．由表6可見，ηk=0.5時非精確N-MCG算法的效率較高．對比表5和表6可見，非精確NMCG算法比N-MCG算法的效率高．兩種算法求出的自反解均為已知的x(0)．

例1和例2的計算結果表明：不論由Newton算法某一步迭代計算導出的線性方程組(3)是否有自反解，選取適當的ηk∈[0,1)時，非精確N-MCG算法與N-MCG算法的效率相比會大幅度提高．如何針對具體問題選取適當的ηk，尚需探討．

6 結論

采用MCG算法求由Newton算法每一步迭代計算導出的線性代數方程組的近似自反解或者近似自反最小二乘解，建立了求非線性代數方程組自反解的非精確Newton-MCG算法．數值算例表明，非精確Newton-MCG算法的效率高于Newton-MCG算法．

參考文獻：

[1]Bellavia S,Morini B.A globally convergent Newton-GMRES subspace method for systems of nonlinear equations[J].SIAM Journal on Scientific Computing,2001,23(3):940-960

[2]Bai Z Z,Guo X F.On Newton-HSS methods for systems of nonlinear equations with positive-definite Jacobian matrices[J].Journal of Computational Mathematics,2010,28(2):235-260

[3]楊愛利，伍渝江，李旭，等.一類非線性方程組的Newton-PSS迭代法[J].計算數學，2012,34(4):329-340 Yang A L,Wu Y J,Li X,et al.On Newton-PSS methods for the system of nonlinear equations[J].Mathematica Numerica Sinica,2012,34(4):329-340

[4]Bai Z Z,Yang X.On HSS-based iteration methods for weakly nonlinear systems[J].Applied Numerical Mathematics,2009,59(12):2923-2936

[5]王洋，伍渝江，付軍.一類弱非線性方程組的Picard-MHSS迭代方法[J].計算數學，2014,36(3):291-302 Wang Y,Wu Y J,Fu J.On Picard-MHSS methods for weakly nonlinear systems[J].Mathematica Numerica Sinica,2014,36(3):291-302

[6]Chen H C.Generalized reflexive matrices:special properties and applications[J].SIAM Journal on Matrix Analysis and Applications,1998,19(1):140-153

[7]張凱院，牛婷婷，聶玉峰.一類非線性矩陣方程對稱解的雙迭代算法[J].計算數學，2014,36(1):75-84 Zhang K Y,Niu T T,Nie Y F.Double iterative algorithm for symmetric solution of a nonlinear matrix equation[J].Mathematica Numerica Sinica,2014,36(1):75-84

[8]Dembo R S,Eisenstat S C,Steihaug T.Inexact Newton methods[J].SIAM Journal on Numerical Analysis,1982,19(2):400-408

[9]張凱院，徐仲.數值代數(第2版)[M].北京：科學出版社，2010 Zhang K Y,Xu Z.Numerical Algebraic(2nd Edition)[M].Beijing:Science Press,2010

[10]Broyden C G.The convergence of an algorithm for solving spare nonlinear systems[J].Mathematics of Computation,1971,25(114):285-294