解線性方程組最小二乘解的迭代法

張光輝

宿州學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,宿州,234000

1 引 言

在技術(shù)和工程實(shí)踐中，最小二乘問題應(yīng)用非常廣泛。所謂最小二乘問題[1]，是指針對(duì)線性方程組AX=b，尋找X*∈Rn，使得

(1)

其中A為一m×n矩陣，b為一m維給定向量。

2 引理及主要結(jié)果

引理1[2]設(shè)A∈Rm×n(m>n),b∈Rm，求解最小二乘問題AX=b，等價(jià)于極小化φ(X)=(AX-b,AX-b)。

ATAX=ATb

(2)

的解。

注：引理1,2從理論上說明求解最小二乘問題(1)只需求解方程(2)，但當(dāng)ATA奇異，或?qū)τ?jì)算過程中不可避免的舍入誤差十分敏感時(shí)，該方法在數(shù)值上是不可取的。應(yīng)用引理3求A+時(shí)需要對(duì)A進(jìn)行奇異值分解。但對(duì)于大規(guī)模的問題，直接采用奇異值分解十分困難，必須采用迭代法。下面給出三種具體的迭代算法。

算法1基于顯示式歐拉算法[5]建立的迭代法

step1根據(jù)AX=b，給出φ(X)=(AX-b,AX-b)表達(dá)式；

step2計(jì)算φ(X)=(AX-b,AX-b)的負(fù)梯度方向-gradφ(X)=-AT(AX-b)；

step4用顯示式歐拉法求解step3微分方程，得到

Xk+1=Xk-hkAT(AXk-b)

并給定算法迭代允許誤差ε；

(3)

其中

step6如果φ(Xk+1)<ε或者‖Xk+1-Xk‖<ε，則算法迭代終止；否則將Xk+1賦給Xk，轉(zhuǎn)step4繼續(xù)迭代。

算法2二階Runge-Kutta算法建立的迭代法

step1-step3同算法1；

step4用二階Runge-Kutta算法求step3微分方程，得到

并給定算法迭代允許誤差ε；

step5計(jì)算φ(Xk+1)，得到關(guān)于hk的四次多項(xiàng)式；

整理得到

(4)

step7如果φ(Xk+1)<ε或者‖Xk+1-Xk‖<ε，則算法迭代終止；否則將Xk+1賦給Xk，轉(zhuǎn)step4繼續(xù)迭代。

算法3基于隱式歐拉算法建立的迭代法

step1-step3同算法1；

step4用隱式歐拉法求step3微分方程，得到

Xk+1=Xk-hkAT(AXk+1-b)

將上式顯化，得

并給定算法迭代允許誤差ε；

step5計(jì)算φ(Xk+1)，得到關(guān)于hk的四次多項(xiàng)式；

整理得到

(5)

step7如果φ(Xk+1)<ε或者‖Xk+1-Xk‖<ε，則算法迭代終止；否則將Xk+1賦給Xk，轉(zhuǎn)step4繼續(xù)迭代。

3 算例分析

算例1某種合金的含鉛量百分比為P，其溶解溫度為θ，由實(shí)驗(yàn)測(cè)得P與θ的數(shù)據(jù)見表1。

表1 合金含鉛百分比與溶劑溫度數(shù)據(jù)表

試用最小二乘法建立含鉛量百分比P與溶解溫度θ之間的經(jīng)驗(yàn)公式θ=x1p+x2。用算法1,2,3計(jì)算，變量的允許誤差ε=0.000 05。

解由題意，所求經(jīng)驗(yàn)公式的待定參數(shù)(x1,x2)即為線性方程組AX=b的最小二乘解向量，其中

矩陣A的條件數(shù)cond(A)=249.79，秩rank(A)=2，由引理3廣義逆矩陣結(jié)論得最小二乘解為X*=A+b=(2.233 7,95.352 4)T。現(xiàn)在用算法1進(jìn)行迭代運(yùn)算，驗(yàn)證算法的有效性。取(x1,x2)初始值為(1,1)，用Marlab軟件編寫程序代碼，迭代7次便可得到具有五位有效數(shù)字的解，迭代結(jié)果如表2所示。

表2 算法1迭代解數(shù)據(jù)

由算法2進(jìn)行迭代運(yùn)算，取(x1,x2)初始值為(2,97)，迭代結(jié)果如表3所示。

表3 算法2迭代解數(shù)據(jù)

算例2求下面線性方程組的最小二乘解。

分析算法1,2，3的誤差向量范數(shù)隨迭代次數(shù)的變化情況。

解矩陣A的條件數(shù)cond(A)=2.5，秩rank(A)=3。

由引理2 所給方程組的最小二乘解即為ATAX=ATb的解

X*=(ATA)-1ATb

=(-0.801 9,-2.487 2,-3.415 1)T

用算法1,2,3分別進(jìn)行迭代計(jì)算，用Matlab編程計(jì)算，得到迭代10次、50次、100次的誤差，結(jié)果見表4。

表4 算法1,2,3迭代結(jié)果的誤差向量范數(shù)

算例2表明，當(dāng)系數(shù)矩陣為良態(tài)時(shí)，迭代矩陣的譜半徑為影響迭代速度的主要因素，隨著譜半徑的減小，收斂速度加快，誤差逐漸變小。

算例3求下面線性方程組的最小二乘解。

解為更好地對(duì)比三種算法的迭代效果，選取如上簡(jiǎn)單模型方程。矩陣A的條件數(shù)cond(A)=1，秩rank(A)=2。系數(shù)矩陣為非奇異方陣時(shí)，方程組的最小二乘解即為精確解，顯然X*=(7,8)T。用三種算法進(jìn)行迭代求解，結(jié)果見表5。

表5 算法1,2,3對(duì)算例3的迭代結(jié)果

算例3表明，當(dāng)系數(shù)矩陣條件數(shù)不大，且迭代矩陣的譜半徑小于1時(shí)，三種算法的迭代速度均有明顯提高，收斂效果很好，且迭代矩陣的譜半徑的大小依然是影響收斂效果的主要因素。

4 結(jié) 語(yǔ)

由于顯式歐拉法、隱式歐拉法、二階龍格-庫(kù)塔算法是求常微分方程數(shù)值解非常經(jīng)典的算法，且原理和公式相對(duì)簡(jiǎn)單，所以本文分析了基于以上三種數(shù)值解法建立的三種迭代法的迭代效果，結(jié)果表明，這三種算法的迭代效果與系數(shù)矩陣的是否病態(tài)，即矩陣的條件數(shù)有很大的關(guān)系。隨著條件數(shù)減小到良態(tài)數(shù)值時(shí)，三種迭代法的迭代速度均有所加快，但速度不同，受迭代矩陣的譜半徑的大小影響。同時(shí)，通過對(duì)3個(gè)算例的計(jì)算，表明文中所建立的三種算法在方程組的最小二乘解的求解過程應(yīng)用中的適應(yīng)性和有效性。算法1對(duì)問題的條件數(shù)適應(yīng)性較強(qiáng)，要求不高，算法2,3 更適應(yīng)于良態(tài)問題，對(duì)高條件數(shù)的問題效果不理想，當(dāng)條件數(shù)相當(dāng)時(shí)，算法2，3相比較而言，算法2精度稍高于算法3。