非對稱代數Riccati方程的一種快速Newton型算法

宋 巖,王小利,鐘瑞華

（閩南師范大學數學與統計學院,福建漳州363000）

本文考慮如下來源于中子運輸理論的非對稱代數Riccati方程：

其中X∈?n×n是未知矩陣,O是零矩陣,已知矩陣Α=Δ-eqT∈?n×n,B=eeT∈?n×n,C=qqT∈?n×n,D =Γ -qeT∈?n×n,且分別是定義在區間[0,1]上的Gauss-Legendre權重和節點.

式(1)可能存在多個解,但只有最小正解(即元素均為正數的矩陣)才具有物理意義[1].本文關注于求解式(1)最小正解的數值算法.Lu[2]在2005年首次將該矩陣方程轉化為一種非線性向量方程,并指出只要通過計算向量方程的最小正解,即可得到矩陣方程的最小正解.具體來說,若X∈?n×n是式(1)的最小正解,則滿足:

于是通過計算式(3)的最小正解,可以得到式(1)的最小正解.

關于式(3)的數值算法,目前主要有基于不動點迭代的研究,如文獻[2-4],及基于Newton 型迭代的研究,如文獻[5-11].由于Newton型迭代具有快速收斂性,故得到了廣泛的關注.例如,Lu[7]在2005年提出了基于Newton法的一種迭代算法.進一步,Lin等[8]在2008年提出了基于一種三階收斂的兩步Newton法.

本文基于文獻[12]所提出的一種具有三階收斂的兩步Newton法,提出了新的計算式(3)的數值迭代算法, 并證明了由該算法迭代產生的向量序列的單調收斂性. 最后, 通過一些數值實驗來比較該算法和Newton法以及文獻[8]中具有三階收斂的兩步Newton法的計算效能.

1 一種三階收斂的兩步Newton型方法

對于u,v∈?n,記x=[uT,vT]T,則式(3)等價于f(x)= 0.給定初始點x0=[uT0,vT0]T∈?2n,考慮如下形式的兩步Newton迭代：

其中xk=[uTk,vTk]T,yk=[uˉTk,vˉTk]T.式(4)的標量形式由Frontini和Sormani[12]在2003年將Newton法和文獻[13]結合而得到.此外,由式(3)可得f的Jacobi矩陣:

基于迭代格式(4),得到如下用于計算式(3)的兩步Newton 型算法.

算法1 一種計算式(3)的兩步Newton 型算法選擇初始向量u0,v0 ∈?n.對于k = 1,2,…,重復以下迭代過程,直至收斂:Step 1.求解關于uˉk,vˉk ∈?n的線性方程組:(In - G2(uk)- ZkH1(uk))vˉk = Zkgk - G2(uk)vk - H2(vk)uk +1 2 (vk + vk °P?uk + e),(In - G1(vk))uˉk = H1(uk)vˉk + gk,其中Ζk = H2(vk)(In - G1(vk))-1,gk=1 2 (uk + uk °Pvk + e)- G1(vk)uk - H1(uk)vk.Step 2.求解關于uk + 1,vk + 1 ∈?n的線性方程組:(In - G2(uˉk)- Ζˉk(H1(uˉk))vk + 1 = Ζˉkgˉk - G2(uˉk)vk - H2(vˉk)uk + vk °P?uk + e,(In - G1(vˉk))uk + 1 = H1(uˉk)vk + 1 + gˉk,其中Zˉk = H2(vˉk)(In - G1(vˉk))-1,gˉk = uk °Pvk + e - G1(vˉk)uk - H1(uˉk)vk.

在初始向量u0,v0∈?n滿足一定條件時,可得到由上述算法迭代產生的向量序列單調收斂于式(3)的最小正解.該收斂性分析將在下一節給出.

2 算法的單調收斂性分析

對于任意的矩陣A,B∈?m×n,A°B=(aijbij)m×n表示A和B的Hadamard 內積. 若對i= 1,2,…,m,j=1,2,…,n,有aij≥bij(aij＞bij),則記A≥B(A＞B);若aij＞0(aij≥0),則稱A是正矩陣(非負矩陣).如果一個向量的所有分量都是正(負),我們稱它為正(負)向量.如果矩陣A的所有非對角元素都是非正的,則稱實方陣A是Z-矩陣.Z-矩陣A可以表示成A=sI-B,其中s∈?,B≥O,I是n階單位矩陣.若s＞ρ(B),則稱Z-矩陣A是非奇異矩陣,其中ρ(B)是B的譜半徑.

引理1[8]A∈?n×n是一個Z-矩陣,則以下條件等價:

1)A是非奇異的M-矩陣;

2)A-1≥O;

3)若Av＞0 ∈?n,則v＞0.

引理2[8]對?x,h1,h2∈?2n,有f″(x)(h1,h2)=[hT1LT1h2,hT1LT2h2,…,hT1LTnh2]∈?2n,其中

pi和是矩陣P和?第i列向量.特別地f?(x)=O,故由Taylor公式有

引理3[7]對任意h,h1,h2∈?n{0},且h1＜h2,則

引理4[7]若ρ(G(x?))≤1,則f′(x?)=I2n-G(x?)是M-矩陣,其中G(x)由式(5)定義.此外對任意x滿足0 ≤x＜x?,f′(x)是非奇異的M-矩陣.

引理5對?v1,v2,u1,u2∈?n,記x=[u1T,v1T]T,y=[u2T,v2T]T,若v1≤v2,u1≤u2,則f′(y)≤f′(x).

證明

故

定理記x?是向量方程f(x)= 0的最小正解,若初始向量x0∈?2n,滿足0≤x0＜x?,且f(x0)＜0,則由算法1迭代產生的向量序列{xk}k≥0和{yk}k≥0有定義,且有以下結論:

a)f(xk)＜0,k≥0;

b)0≤xk＜yk＜xk+1＜x?,k≥0;

證明利用數學歸納法證明a)和b).當k= 0 時,利用式(4)及引理1-4 易得a)和b)成立.假設k=l時,a)和b)成立,下面考慮k=l+ 1的情形.由式(4)和引理5有

即xl+1-yl＞yl-xl.進而,由引理3得

故由式(6)得

故當k=l+ 1時,f(xl+1)＜0,即a)成立.

對于b),由式(4)和式(6)得

則

得yl+1＜x*,故由引理4知f′(yl+1)是非奇異的M-矩陣,從而由引理1得f′(yl+1)-1＞O.當k=l+ 1時,由式(4)得xl+1＜yl+1,xl+1＜xl+2.由式(6)得及引理5得

進而由式(4)和式(6)得

注意到f′(yl+1)是非奇異的M-矩陣,故由引理1知xl+2＞yl+1.于是由式(4)和式(6)得

即x*-yl+1＞yl+1-xl+1,則由引理3知

由式(4)和式(6)得

則

故xl+2＜x*.由此可知,當k=l+ 1時,b)成立.因此,對任意k≥0時,a)和b)都成立.

3 數值實驗

本節通過一些數值實驗來比較算法1(記為TSNM)和文獻[7]中的基于Newton 法(記為NM)以及文獻[8]中的三階收斂的兩步Newton 法(記為MNI)的計算效能.分別選取(α,c)=(10-7,1 - 5×10-7)和(α,c)=(0.3,0.7)進行測試,在每個測試中,給出迭代次數(記為IT)和相對殘差(記為ERR).

維數n分別取512,1 024,2 048,4 096,節點ci和權重ωi,i= 1,2,…,n,從區間[0,1]作數值積分得出.具體步驟,將區間[0,1]分為n4個長度相等的區間,并對每個區間應用Gauss-Legendre 四點求積公式.

我們選取初始點u0=v0=0∈?n,算法終止的條件是

其中eps = 2-52≈2.2204×10-16,范數取無窮范數, 所有實驗均在CPU - 2.30GHz（Intel(R)Core(TM)i5-6200），RAM - 4GB環境下進行,MATLAΒ版本為2018b.數值實驗結果如下表1-表2所示:

表1 α = 0.3,c = 0.7 時的數值模擬Tab.1 The numerical result for α = 0.3,c = 0.7

續表1

表2 α = 10-7,c = 1 - 5×10-7 時的數值模擬Tab.2 The numerical result for α = 10-7,c = 1 - 5×10-7

從以上數值實驗結果可以發現TSNM 法也適用于接近奇異的情況,并且無論n取512,1 024,2 048還是4 096,TSNM 法比NM 法和MNI法的迭代次數少,殘差更小.因此,可以認為TSNM 法在求式(3)的最小正解的問題上更優于NM法和MNI法.

圖1-圖2 是固定n取4 096,(α,c)為(0.3,0.7)和(10-7,1 - 5×10-7)的TSNM 法,NM 法和MNI 法收斂過程.由圖可知,在非奇異與接近奇異的情形下,TSNM 法的收斂速度都比Newton法和MNI法更快,特別是接近奇異時,優勢更明顯.

圖1 α = 0.3,c = 0.7,n = 4 096Fig.1 α = 0.3,c = 0.7,n = 4 096

圖2 α = 10-7,c = 1 - 5×10-7,n= 4 096Fig.2 α = 10-7,c = 1 - 5×10-7,n = 4 096