基于隨機Kaczmarz算法的最小二乘擬合

楊 紅，陳豫眉

(西華師范大學(xué) a.數(shù)學(xué)與信息學(xué)院; b.公共數(shù)學(xué)學(xué)院, 四川南充 637009)

0 引言

隨機Kaczmarz算法(簡稱RK算法)[1]是目前求解大型超定線性方程組的經(jīng)典方法之一，當系數(shù)矩陣的行數(shù)大于矩陣的列數(shù)的三倍時，該算法甚至優(yōu)于著名的共軛梯度算法.2009年，Strohmer和Vershynin提出了依指數(shù)線性收斂的RK算法[1].隨后，一些學(xué)者對RK算法的補充證明[2-6]進一步說明了該算法的收斂性.RK算法利用交替投影的方法進行迭代，由于其算法簡單，運行速度非?？欤粡V泛應(yīng)用到了許多實際問題中，如數(shù)字信號處理、X射線斷層掃描、圖像重構(gòu)等[7-9].

最小二乘法是一種優(yōu)化方法，它通過最小化殘差平方和來尋找適合數(shù)據(jù)的最佳匹配函數(shù).最小二乘擬合在數(shù)據(jù)預(yù)測方面有著重要應(yīng)用，如化學(xué)、物理模型中的參數(shù)估計[10-13]等.最小二乘擬合求解方法有很多，最常用的如SPSS、MATLAB中的擬合工具箱[14-15]等.2010年，劉佳琪等提出利用進化算法求解最小二乘擬合[16]；2011年，顏清等提出利用EXCEL中的蒙特卡羅算法和VBA混合編程進行最小二乘擬合[17].當實驗數(shù)據(jù)龐大時，最小二乘法擬合的實質(zhì)是求解超定線性方程組.目前已有的最小二乘擬合的求解方法，在一定程度上雖能達到較好的擬合效果，但當實驗數(shù)據(jù)龐大時，將面臨計算量大的困難.因此, 本文利用求解大型超定線性方程組的RK算法進行最小二乘擬合，以更好地應(yīng)用于實際問題.

1 最小二乘法理論

最小二乘法起源于天文學(xué)和地理測量學(xué)的研究.1805年，勒讓德在其論著《計算慧星軌道的新方法》中給出了最小二乘法清晰且簡明的闡述.1809年，高斯在《關(guān)于繞日行星運動的理論》中也提出了最小二乘法.下面給出最小二乘法的一般理論.

設(shè)有一組實驗數(shù)據(jù)(xi,yi),i=1,2,...,m，最小二乘擬合即求函數(shù)f(x)使得

(1)

用最小二乘法擬合曲線，可根據(jù)給定數(shù)據(jù)散點圖確定擬合函數(shù)f(x)的形式.一般地，設(shè)待求函數(shù)f(x)屬于函數(shù)空間

φ=span{φ1,φ2,…,φn}，

即f(x)可由函數(shù)φj,j=1,2,…,n的線性組合表示

故式(1)等價于求解

(2)

則式(2)等價于求解

(3)

當函數(shù)空間φ=span{1,x,x2…,xn}時，即為常用的最小二乘多項式擬合.

2 基于RK算法的最小二乘曲線擬合

2.1 RK算法

給定超定線性方程組

Ax=b，

RK算法[1]通過正交投影建立相應(yīng)的迭代公式來求解超定線性方程組.下面給出其主要思想.

其中ik∈{1,2,...,m}是根據(jù)上述概率隨機選擇得到的行指標.

其算法步驟為

1)輸入超定線性方程組中的系數(shù)矩陣A和右端列向量b；

2)規(guī)定最大迭代次數(shù)K，賦迭代次數(shù)初值k=0，賦迭代初值x0；

3)當?shù)螖?shù)k

6)更新迭代次數(shù)k=k+1，返回步3.

2.2 基于RK算法的最小二乘擬合

給定函數(shù)f(x)∈φ=span{φ1,φ2,...,φn}，對已知的實驗數(shù)據(jù)(xi,yi),i=1,2,...,m進行擬合.將實驗數(shù)據(jù)帶入擬合函數(shù)f(x)得到一個超定線性方程組

Xβ=y

(4)

其中X是m×n階矩陣，y是m×1階向量，β是n×1階向量.本文將通過RK算法計算出擬合參數(shù)向量β.

當實驗數(shù)據(jù)龐大時，式(4)是一個大型超定線性方程組，RK算法求解大型超定線性方程組時具有依期望指數(shù)收斂的性質(zhì)，故用RK算法進行最小二乘擬合.

3 數(shù)值實驗及結(jié)果

RK算法中涉及隨機因素，為了避免隨機產(chǎn)生的影響，將數(shù)值實驗運行50次所得平均值作為最終結(jié)果.設(shè)算法中規(guī)定的最大迭代次數(shù)K=20000，實驗運行環(huán)境為MATLAB R2018a.

文中利用的擬合相關(guān)系數(shù)計算公式為

例1 利用RK算法對表1中的實驗數(shù)據(jù)進行最小二乘直線擬合.

表1 直線擬合的實驗數(shù)據(jù)

擬合函數(shù)表達式為f(x)=β1+β2x，將實驗數(shù)據(jù)帶入該表達式，得到如下超定線性方程組:

令

通過RK算法計算得(β1,β2)=(-2.0365,22.7982)，R2≈0.9984，擬合效果好.

例2 利用文獻[17]中給出的化工動量傳遞和熱量傳遞的實驗數(shù)據(jù)，擬合傳熱準數(shù)經(jīng)驗關(guān)聯(lián)式.

表2 化工動量傳遞和熱量傳遞的實驗數(shù)據(jù)

傳熱準數(shù)經(jīng)驗關(guān)聯(lián)式由普蘭德數(shù)Pr、雷諾爾數(shù)Re和努塞爾數(shù)Nu三個參數(shù)控制，其表達式為

Nu=aRemPrn.

對于這個非線性函數(shù)的擬合，兩邊取以10為底對數(shù)，得到

lgNu=lga+mlgRe+nlgPr

將實驗數(shù)據(jù)帶入上述表達式，得到一個超定線性方程組.通過RK算法計算得參數(shù)a=0.1193，m=0.6628，n=-0.3220.計算該曲線擬合的相關(guān)系數(shù)R2≈0.9999，高于文獻[17]中利用蒙特卡羅擬合得到的相關(guān)系數(shù)R2≈0.9985.可見RK算法進行最小二乘數(shù)據(jù)擬合有一定的優(yōu)勢.

圖1 測試數(shù)據(jù)散點圖

圖2為通過RK算法擬合曲線，擬合相關(guān)系數(shù)R2≈1.0000.此例中的擬合方程組個數(shù)為300,說明RK算法用于大型數(shù)據(jù)最小二乘擬合效果非常好.

圖2RK算法擬合的曲線

4 結(jié)語

最小二乘法的應(yīng)用極其廣泛，在數(shù)據(jù)擬合方面有著重要的作用.然而，當實驗數(shù)據(jù)增多時，運用一般的方法進行擬合求解會產(chǎn)生龐大的計算量.本文將求解超定線性方程組的RK算法應(yīng)用于最小二乘擬合，效果非常好.