求解一類代數(shù)Riccati 方程的Anderson加速算法

王小利，宋 巖，凌永輝

(閩南師范大學數(shù)學與統(tǒng)計學院，福建漳州 363000)

產(chǎn)生于運輸問題的一類非對稱代數(shù)Riccati 方程具有以下形式：

Juang[1]在1995 年證明了矩陣方程（1）正解的存在性，并在文獻[2]中給出了一種近似計算最小正解的數(shù)值算法.考慮到實際應用及解的物理意義，最小正解的計算是研究者著重關(guān)注的問題.Lu[3]在2005 年證明了方程（1）的解滿足以下形式：

其中?表示Hadamard 乘積，n階方陣T=(tij)= (1δ i+γj)，n維向量u和v滿足

于是，求解矩陣方程（1）的最小正解等價于求解非線性向量方程（3）的最小正解.

近幾年來，關(guān)于求解方程（3）最小正解已經(jīng)取得了大量的成果.Lu[3]在2005 年提出一種不動點迭代（Fixed Point Iteration，F(xiàn)P）法來求解方程（3）的最小正解，隨后又在文獻[4]中進一步給出求解該方程的Newton 法.隨后，Lin 等[5]改進文獻[3]中的不動點算法得到一種收斂更快的不動點算法.在前人的基礎(chǔ)上，Bai 等[6]在2008 年基于矩陣分塊的Jacobi 迭代和Gauss-Seidel迭代給出了NBJ（Nonlinear Block Jacobi）法和NBGS（Nonlinear Block Gauss-Seidel）算法，這兩種算法可用于計算接近奇異情況下的問題.隨后，Lin[7]又對文獻[6]中方法進行改進.同年，Lin 等[8]提出基于兩步Newton 法的數(shù)值算法，并且證明了其單調(diào)收斂性.Ma 等[9]在2016 年提出了修正的Jacobi 型不動點迭代法和修正的Newton 型迭代法.Miyajima[10]在2017 年根據(jù)解的特殊結(jié)構(gòu)提出一種適用于非奇異情況的直接求解算法.Zhang 等[11]在2020 年提出了帶參數(shù)的Newton 外推法以及改進的倍增算法，并利用Cayley 變換給出了收斂性分析.其他方法詳見文獻[12-14]及其參考文獻.

為改善不動點迭代在求解大規(guī)模方程接近奇異情況下收斂速度較慢的問題，本文考慮了一種新的不動點加速算法，這種加速方法最早是由Anderson[15]在1965 年提出的.Fang 等[16]在2009年提出了基于該加速方法的多割線迭代法，并將其應用于電子結(jié)構(gòu)計算.Walker 等[17]在2011 年證明了在一定條件下，Anderson 加速與GMRES（Generalized Minimal Residual）是等價的，并討論了幾種算法的實施.Toth 等[18]在2015 年給出了Anderson 加速應用在線性及非線性問題上的收斂性分析.基于該加速框架，本文給出了求解方程（3）的一種不動點加速算法，并與FP、NBJ、NBGS 三種不動點迭代算法進行數(shù)值實驗比較.非奇異情況下的數(shù)值結(jié)果顯示，該加速算法較FP、NBJ、NBGS 法有較好的計算效能.在接近奇異情況下，較FP、NBJ、NBGS 三種算法均有良好的計算效能，特別是在計算接近奇異的大規(guī)模問題時，該加速算法優(yōu)勢尤為明顯.

1 Anderson 加速算法

對于非線性方程（3），Lu[3]在2005 年給出了求解該方程的不動點迭代形式：

該算法計算簡單，運算量只有 O(n2)，在非奇異情況下具有良好的計算效能.但在解決接近奇異的問題時，收斂速度緩慢.特別是，計算接近奇異的大規(guī)模問題時，需要花費很長時間.為了解決這一問題，我們利用Anderson 加速對不動點迭代法（4）進行改進.

下面給出Anderson 加速的一般算法.

算法1 Anderson 加速法.

算法1 實現(xiàn)的主要困難是通過求解帶約束最小二乘問題（5）來確定參數(shù)根據(jù)Fang等[16]的推導，上述算法的帶約束最小二乘問題可以改寫為無約束問題：

這里α與γ間的關(guān)系為

從而，算法1 中計算下一步迭代值的等式對應改寫為：

2 最小二乘問題的求解

對于帶約束最小二乘問題（5），Walker 等①Walker H,Ni P.A Linerly Consterained Least-Squares Problem in Electronic Structure Computations [C].International Conference on Computational &Experimental Engineering and Science,2010:43-49.在2010 年提出在一定條件下可以利用Lagrange乘數(shù)法直接求解.對于無約束的最小二乘問題（6），選用較為穩(wěn)定的QR 分解進行計算.

這樣求解無約束的最小二乘問題（6）可轉(zhuǎn)化為求解一個mk×mk的上三角方程組.下面計算Fk的QR 分解.

下面給出Anderson 加速求解非對稱代數(shù)Riccati 方程的算法描述.

3 數(shù)值實驗

分別使用算法2（AA 算法）、不動點迭代[3]（FP 算法）、文獻[6]所提出的非線性Jacobi 分塊法（NBJ 算法）和非線性Guass-Seidel 分塊法（NBGS 算法）來進行數(shù)值實驗，比較這4 種算法在α、c的不同取值情形下的CPU 時間（單位為秒）、迭代次數(shù)IT 和相對殘差ERR，其中CPU時間取算法計算十次的平均值.

將區(qū)間[0,1]進行劃分，再利用四點的Guass-Legendre 求積公式求出每個子區(qū)間的節(jié)點ωi和權(quán)重ci，i=1,2,…,n，從而可以求出矩陣的P 和，其中維數(shù)n分別選512、1 024、2 048、4 096和8 192.

其中 eps ≈ 2.220 4 × 1 0?16.

所有實驗均在CPU-3.60GHz（Intel(R) Core(TM) i7-4790），RAM-4GB 環(huán)境下進行，MATLAB版本為2013a.

在迭代終止條件相同的情況下，矩陣維數(shù)n固定為8 192 時，表1 為不同參數(shù)數(shù)值的實驗結(jié)果.可以看出，在非奇異情況下，算法AA 與NBGS 計算效率相似，但隨著參數(shù)的取值逐步接近奇異情況，算法AA 的優(yōu)勢也越來越明顯.不管是CPU 時間、迭代次數(shù)還是殘差，都明顯優(yōu)于其他三種算法.

表1 n= 8 192時不同參數(shù)的數(shù)值模擬

表2 為選定參數(shù) (α,c)= (10?7,1 ? 5 × 10?7)，矩陣維數(shù)n取不同值的實驗結(jié)果.可以發(fā)現(xiàn)，該方法在計算接近奇異情況的大規(guī)模問題時，相較其他三種算法，算法AA 效果尤為顯著，大大節(jié)省了計算成本.

表2 α= 10 ?7 ,c=1 ? 5×10?7時不同矩陣維數(shù)的數(shù)值模擬

表3 為(α,c)= (0.001,0.999)，矩陣維數(shù)n固定為4 096 時不同m的實驗結(jié)果.

表3 α= 0.001,c= 0.999,n= 4 096時不同m 的數(shù)值模擬

表4 為α= 10?5，c=1 ? 10?4，矩陣維數(shù)n為2 048 時不同m的實驗結(jié)果.

表4 α= 10 ?5 ,c=1 ? 10?4,n= 2 048時不同m 的數(shù)值模擬

對比收斂較快的NBGS 法可以得出，選取適當?shù)膍，迭代次數(shù)與誤差也會有所改變，選取最優(yōu)的m可在一定程度上減少不必要的運算.當參數(shù)(α,c)接近奇異到一定程度時，即使不取最優(yōu)的m，也有較好的收斂效果.這表明算法AA 對于求解方程（3）最小正解是十分有效的.

圖1 為不同參數(shù)的收斂過程.圖1（a）為取參數(shù) (α,c) = (0 .01,0 .9 9)，m=12，矩陣維數(shù)n= 512時不同方法的數(shù)值實驗結(jié)果；圖1（b）為取參數(shù) (α,c)=(1 0?6,1 ?10?5)，m=7，矩陣維數(shù)n= 4 096時的實驗結(jié)果.可以看出，隨著參數(shù)的取值接近奇異情況，較其他三種方法，算法AA 具有良好的計算效能.

圖1 不同參數(shù)的收斂過程