關于Matlab多項式數據擬合命令polyfit的一點思考

摘要：針對學生在使用Matlab多項式擬合命令polyfit過程中遇到的警告問題，文章從多項式擬合、最小二乘法、方程組的穩定性和條件數等數學理論出發，詳盡地闡述了polyfit命令背后的原理，解釋了代碼運行過程中出現警告的原因和解決方案。同時，還通過具體的實驗介紹了四種異常情況和對應解決方案的條件數對比，直觀地展示了各種解決方案的效果。理論和實驗分析相結合，不僅可以幫助學生理解polyfit警告背后的原因并提供解決方案，還可以培養學生的探索性思維。

關鍵詞：數學實驗；多項式擬合；最小二乘法；條件數；病態方程; Matlab; polyfit函數

中圖分類號：TP312" " " 文獻標識碼：A

文章編號：1009-3044（2025）12-0107-05

開放科學（資源服務） 標識碼（OSID）

0 引言

我們在講授《數學實驗》課程“數據擬合”這一節的時候，會用到Matlab的多項式擬合命令polyfit[1-2]。同學們在實際使用過程中，經常會遇到polyfit出現如下警告的情形：

Warning： Polynomial is badly conditioned. Add points with distinct X values， reduce the degree of the polynomial， or try centering and scaling as described in HELP POLYFIT [3]。

雖然Matlab給出了解決問題的建議，但同學們仍然不明白出現警告的原因，更不知道如何消除此警告。本文的目的是幫助學生理解警告背后的數學原理，并提供實際的解決方案，具體從理論和實驗兩方面，包括多項式擬合的原理、矩陣的條件數、病態方程、polyfit函數的實現過程和實驗等全面解析出現警告的原因，并給出解決方案。

1 最小二乘法與數據擬合

數據擬合問題指的是，已知一組數據[xi ， yi， i=1 ， 2， ... ，N]，尋求一個擬合函數[y=fx]，使將[xi]代入得到的[y=fxi]在某種準則下與[yi]最為接近[（i=1，2， ... ，N）][4]。在幾何上，可以理解為在[N]個數據點確定的平面上，尋找一條曲線使其與這些數據點在某種準則下距離最小。擬合函數有多項式、指數和對數函數等。

最常用的準則是使[fxi]與[yi]之差的平方和最小，稱為最小二乘準則。根據最小二乘準則計算得到的待定系數的解稱為最小二乘解，相應的方法稱為最小二乘法。

2 基于最小二乘法的多項式擬合

2.1 多項式擬合的概念

多項式擬合是指用多項式函數來逼近給定的數據點，并利用得到的多項式函數預測未知的數據點[5]。用多項式函數近似表示函數有很多優點，如易計算函數值，對數據形狀沒有要求，而且擬合的精度和復雜度可以通過改變多項式的階數來調整[5]。因此多項式擬合在理論分析、近似計算和實際問題中有非常廣泛的應用[6-8]。例如，多項式擬合可用于圖像處理的邊緣檢測和圖像增強，信號處理的信號平滑和去噪。

2.2 多項式擬合的原理

設擬合函數為[n]次多項式[y=f（x）=p1xn+p2xn-1+...+pnx+pn+1]。把點[xi ， yi（i=1，2，...，N）]代入[y=fx] ，便得到以[p1，p2，...，pn+1]為未知量的方程組：

[p1xn1+p2xn-11+...+pnx1+pn+1=y1p1xn2+p2xn-12+...+pnx2+pn+1=y2 …… ……p1xnN+p2xn-1N+...+pnxN+pn+1=yN]" " " （1）

記：

[V=xn1xn-11…1xn2xn-12…1……xnNxn-1N…1N×（n+1）]，[p=p1p2?pn+1]，[y=y1y2?yN]

則（1） 的矩陣形式為：

[Vp=y]" " " " " " " "（2）

其中[V∈RN×（n+1）]，[p∈Rn+1]，[y∈RN]。

多項式階數過高，模型會變得很復雜，容易產生過擬合現象，因此數據個數[N]遠大于多項式的階數[n]。由于實際數據往往存在噪聲，多項式擬合曲線不一定經過每一個數據點，因此方程組（2） 不一定嚴格成立。也即，對于實際觀察數據[xi ， yi（i=1，2，...，N）]，方程組（2） 無解，為矛盾方程[9]。一種常見的準則是讓多項式擬合函數在每個數據點的預測值盡可能接近其真實值，也就是第1節介紹的最小二乘準則，其目標函數為[10]：

[minp1 ，p2 ，...，pn+1i=1Nyi-p1xni+p2xn-1i+...+pnxi+pn+12] （3）

可以寫成矩陣形式：

[minpy-Vp2=minp（pTVTVp-2yTVp+yTy）] （4）

（4） 關于[p]求導并令它等于0，得最優性條件：

[VTVp=VTy]" " " " " "（5）

因此矛盾方程組（2） 的最小二乘解[p1，p2，...，pn+1]，就是方程組（5） 的解[11-12]。

下面證明方程組（5） 在[V]列滿秩，即[rank（V）=n+1]，且[Ngt;n+1]下有唯一解。

引理：設非齊次線性方程組：

[Bx=b]" " " " " " "（6）

其中[B∈Rs×t]，[b∈Rt]。若[rank（B）=t]，即列滿秩，則：

（1） 矩陣[BTB]是對稱正定矩陣。

（2） 方程組[BTBx=BTb]有唯一解。

證：（1） 因為[（BTB）T=BTB]，所以[BTB]是對稱矩陣。

對齊次線性方程組[Bx=0]，因為[rank（B）=t=未知數的個數]，所以方程組有唯一解，也即只有零解。所以對[?x≠0]，有[Bx≠0]，從而[（Bx）T（Bx）=xT（BTB）xgt;0]，即[BTB]是正定矩陣。

從而[BTB]是對稱正定矩陣。

（2） 因為[BTB]是正定矩陣，所以[rank（BTB）=t]，從而方程組[BTBx=BTb]有唯一解。

定理：設[x1，x2， ... ， xn+1∈Rn]互不相同，且[Ngt;n+1]，則方程組（5） 有唯一解。

證：取[V]的前[n+1]行構成[n+1]階子矩陣[V1]：

[V1=xn1xn-11...1xn2xn-12...1......xnn+1xn-1n+1...1]" " " （7）

則[V1]為[n+1]階范德蒙矩陣。因為[x1，x2， ... ， xn+1]互不相同，所以[|V1|≠0]，從而[rank（V1）=n+1]。[V]是[N×（n+1）]，而[Ngt;n+1]，因此[rank（V）]最大等于[n+1]。而它的子矩陣[V1]的[rank（V1）=n+1]，因此[rank（V）=rank（V1）=n+1]，即[V]是列滿秩的。

由引理可得，方程組（5） ：[VTVp=VTy]有唯一解。

3 方程組解的穩定性和條件數

對線性方程組（6） ，其中[B]可逆，設[b]有小擾動[Δb]時解的擾動為[Δx]。如果觀測值的微小改動[Δb]不會使方程的解 [x]發生劇烈變化，就稱該方程是穩定的。也就是說，穩定性衡量的是線性方程的抗噪聲能力。

由：

[B（x+Δx）=b+Δb]" " " "（8）

因為[Bx=b]，從而[BΔx=Δb]。又因為[B]可逆，可得：

[Δx=B-1Δb]" " " " " " （9）

由[b=Bx≤B?x]和式（9） ，可得：

[Δx=B-1Δb≤B-1?Δb]" " " " " （10）

其中[ ? ]表示范數。最常用的范數有Frobenius范數、1-范數、2-范數和[∞]-范數。例如，矩陣[A]的Frobenius范數為矩陣[A]的所有元素的平方和的平方根，矩陣[A]的2-范數為矩陣[A]的最大奇異值，矩陣[A]的[∞]-范數為矩陣[A]的元素絕對值的最大值。

考慮相對誤差，利用[b=Bx≤B?x]和式（10） ，有：

[Δxx≤ΔxbB=Bb?Δx≤Bb?B-1?Δb≤B?B-1?Δbb]" " "（11）

即：

[Δxx≤B?B-1?Δbb]" " " " "（12）

稱cond[（B）=B?B-1]為[B]的條件數，表示線性方程組[Bx=b]的解[x]受觀測值[b]的噪聲的影響程度。

在求線性方程組的過程中，由于計算機的有限精度或者觀測數據的誤差等原因，矩陣[B]和向量[b]可能存在一定的誤差。當矩陣[B]的條件數很大時，即使矩陣[B]和向量[b]的相對誤差很小，解[x]的相對誤差也可能很大，用常規的數值方法（如高斯消元法） 可能導致得到的結果不穩定和很大的偏差。這種矩陣稱為“病態的”，相應的方程組為病態方程組。如果矩陣[B]的條件數比較小，那么[b]的小變化只會導致解[x]的較小的變化。這種矩陣稱為“良態的”，相應的方程為良態方程組[13-14]。

下面來看一個例子。考慮線性方程組為：

[2626.00001xy=88.00001]

它的準確解為[xyT=11T]。

如果右端常數項[b]發生微小的變化，得：

[2626.00001xy=88.00002]

它有準確解[xyT=-22T]。這里：

[Δb=00.00001T]，[B=2666.00001]

易求[B]的條件數為：

cond[∞（B）=B∞B-1∞=4.8×106]

解的相對誤差界為：

[Δx∞x∞≤cond∞（B）?Δb∞b∞≈600%]

說明矩陣[B]是嚴重病態的，相應的線性方程組是病態方程組。

4 Matlab中polyfit實現過程和異常情況

4.1 Matlab求方程組的解

對線性方程組（2） ：[Vp=y]，根據矩陣[V]的不同情況，表1列出了Matlab中不同的實現過程。

表1" Matlab求線性方程組[Vp=y]不同方法的比較

[Matlab命令 含義 優缺點 y=inv（V） |V|≠0時給出方程組的唯一解 只能求系數矩陣可逆時的解

運行效率不高 V不可逆或不是方陣時，給出警告信息 y=V＼p V不可逆或不是方陣時，求最小二乘解[15] 可以求系數矩陣不可逆或不是方陣時的解

利用奇異值分解，求解速度較快 y=pinv（V）*p |V|=0或V不是方陣時，利用偽逆求方程組的解 可靠性比奇異值求法略差

但速度快[16] 利用矩陣的分解來求線性方程組的解，如Lu（V）， qr（V）， chol（V）， schur（V）， hessenberg（V） 常見的矩陣分解有 LU 分解、QR 分解、Cholesky 分解，以及 Schur 分解、Hessenberg 分解、奇異分解等[17] 運算速度快

節省存儲空間 ]

4.2 polyfit實現過程

根據2.2節的介紹，多項式數據擬合問題轉化為求解方程組（2） ：[Vp=y]的最小二乘解，也即方程組（5） 的解。Matlab的多項式擬合命令polyfit利用矩陣的QR分解來求方程組的解。

對于任何[m×n]的矩陣[A]，QR分解可以表示為：

[A=QR]" " " " " " " "（13）

其中[Q]是一個[m×m]的正交矩陣，即[Q*Q=I]，這里[Q*]表示[Q]的共軛轉置。[R]是一個[m×n]的上三角矩陣，其對角線上的元素都是非負實數。

QR分解的實現方法有：Gram Schmidt正交化，Householder變換和Givens變換等。下面給出Gram Schmidt正交化方法的步驟：（1） 將[A]的列向量正交化；（2） 將正交化的列向量單位化；（3） 寫出[Q]和[R]。

回到Matlab中polyfit函數的實現過程。它首先將系數矩陣[V∈RN×（n+1）]進行QR分解：

[V=QR]" " " " "（14）

其中，[Q]為[n+1]階正交矩陣，[R]為[n+1]階上三角矩陣[18]。

則方程組（2） 化為：

[QRp=y]" " " " " "（15）

由于[Q]為正交矩陣，所以[Q-1=QT]，因此得到：

[Rp=QTy]" " " " " " "（16）

再利用Matlab左除命令，可得方程組（16） 也即方程組（2） 的最小二乘解為：

[p=R＼（QTy）]" " " " " "（17）

Polyfit的部分代碼如下所示：

%構建Vandermonde矩陣 V = [x.^n ... x.^2 x ones（size（x））]

V（：，n+1） = ones（length（x），1，class（x））;

for j = n：-1：1

V（：，j） = x.*V（：，j+1）;

end

% 求解最小二乘問題p = V＼y 得到擬合多項式系數 p.

[Q，R] = qr（V，0）;

oldws = warning（'off'，'all'）;" "% Turn all warnings off before solving

p = R＼（Q'*y）;" " " " " nbsp; " "% Same as p = V＼y

4.3 polyfit異常情況

實際使用polyfit命令時，可能出現四種異常情況，表2列出了它們的觸發條件和解決方案。

表2" Polyfit異常情況、觸發條件和解決方案

[序號 異常情況 觸發條件 解決方案 1 polyfit： X、Y Size Mismatch X和Y的維度不一致 檢查X和Y的維度，確保一致 2 polyfit： Poly Not Unique 矩陣[R]的列數大于行數，此時方程組（16） 欠定，解不唯一 添加更多的不同的擬合數據點，或者減少多項式的階數 3 polyfit： Repeated Points 調用格式為[p，S，mu] = polyfit（x，y，n），且方程組（16） 的系數矩陣[R]的條件數過大（單精度大于[105]，其他情況大于[1010]） 增加更多的具有不同x值的數據點，降低擬合多項式的階數等，以降低系數矩陣的條件數，改善方程組的“病態” 4 polyfit： Repeated Points Or Rescale 調用格式為P = polyfit（x，y，n）或[P，S] = polyfit （x，y，n），且方程組（16） 的系數矩陣[R]的條件數過大（單精度大于[105]，其他情況大于[1010]） 增加更多的具有不同x值的數據點，降低擬合多項式的階數，或者對數據進行中心化和縮放，以降低系數矩陣的條件數，改善方程組的“病態” ]

5 實驗和結果分析

下面本文通過一些具體實驗來分析產生錯誤的原因，并給出解決方案。

5.1 實驗1：矩陣[R]的條件數過大

請根據下列1750—2000年的世界人口數，繪制一張世界人口增長曲線示意圖，并用5次多項式進行擬合。

clear all

x = [1750：25：2000]';

y = 1e6*[791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';

plot（x，y，'o'）

p = polyfit（x， y， 5）;" " " %5階多項式擬合

為了輸出矩陣R的條件數，這里需要把polyfit.m第93行代碼：

cond2=condest（R）;

后面的分號去掉。

運行上述代碼后，得到圖1的散點圖，R的條件數[cond（R）=3.7869×1023]，同時命令窗口彈出警告：

Warning： Polynomial is badly conditioned. Add points with distinct x values， reduce the degree of the polynomial， or try centering and scaling as described in HELP POLYFIT[3].

即多項式條件很差。增加具有不同x的數據點，減少多項式的度，或者嘗試按POLYFIT中的幫助文檔對x進行中心化或縮放。

由于矩陣R的條件數[cond（R）=3.7869×1023gt;1010]，且調用polyfit的輸出參數只有1個，對照4.3節的異常情況，此時觸發了異常情況4，因此命令窗口彈出警告Polynomial is badly conditioned，告訴用戶，此問題的條件數很差，并給出了三種可行方案。下面本文一一實驗這三種方法的效果。

1） 方案1：添加具有不同 X 值的點（Add points with distinct x values）

可能原因為：數據是否有過于相近甚至相同的數值。從圖1可以看出，前面幾個數據還是比較接近的。假設我們增加三個數據：

clear all

x = [1675：25：2000]';

y = 1e6*[10 200 300 791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';

p = polyfit（x， y， 5）;

此時[cond（R）=8.5813×1022]，條件數略為降低，但警告仍未消除，效果不明顯。主要原因在線性方程組（2） ：[Vp=y]的系數矩陣[V]的列是向量x的冪，當x的各個分量比較大，多項式的階數比較高時，[V]的各列元素之間差異非常大，導致[V]的條件數通常較大，生成一個奇異系數矩陣[19]。

2） 方案2：減少多項式的階數（reduce the degree of the polynomial）

將擬合多項式次數改為2，重新運行代碼，這時警告消失。經查驗，[cond（R）=2.2510×109lt;1010]，[R]的條件數得到了很大改善，方程組不再是病態方程組，故警告消失。由情況1） 的分析可知，降低多項式的次數會大大降低系數矩陣[V]的條件數。

3） 方案3：進行中心化和縮放（try centering and scaling）

如果在第2行代碼后面增加中心化和縮放的代碼：

mu = [mean（x）; std（x）];

x = （x - mu（1））/mu（2）;

即將 x 的中心置于零值處，并縮放為具有單位標準差，再進行5次多項式擬合，則警告消失。

中心化指的是，從每個原始數據點中減去數據集的均值[μ]，將數據的中心移動到零點。

縮放指的是，將中心化后的每個數據點除以數據集的標準差[σ]，將數據縮放到相同的尺度。這一步是為了消除數據的尺度（或單位） 差異，使得不同特征或不同數據集之間的比較更加公平和有意義。

通過這兩步操作，原始數據被轉換為一個具有零均值和單位方差的新數據集。這個過程也叫作標準化，具體公式為：

[z=x-μσ]" （18）

其中數據[z]具有零均值和單位標準差。

完整代碼如下：

clear all

x = [1750：25：2000]';

%%%新增加的中心化代碼

mu = [mean（x）; std（x）];

x = （x - mu（1））/mu（2）;

y = 1e6*[791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';

plot（x，y，'o'）

p = polyfit（x， y， 5）;" " " %5階多項式擬合

此時輸出[cond（R）=50.6278lt;1010]，[R]的條件數得到了極大改善，不會觸發警告。

也可以將上面的中心化代碼刪除，將p = polyfit（x， y， 5）替換為：

[p，S，mu] = polyfit（x，y，n）

其中，mu為一個二元素向量，包含中心化值和縮放值；mu（1）是mean（x），mu（2）是std（x）；（mean：平均值，std：標準差）。使用該命令時，polyfit將x的中心置于零值處并縮放為具有單位標準差，其本質也是對數據進行了中心化和縮放處理。完整代碼如下：

clear all

x = [1750：25：2000]';

y = 1e6*[791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';

plot（x，y，'o'）

[p，S，mu] = polyfit（x，y，5）;" " " %5階多項式擬合

由情況1） 的分析可知，將原始數據中心化和縮放（即進行標準化） 后，系數矩陣[V]的各個分量差異變小，可以大大降低[V]的條件數。

綜上，當條件數很差時，最佳的方法是對原始數據進行標準化處理，這樣可以大大降低系數矩陣的條件數，使病態方程組變為良態方程組，從而求出方程組的有效解。

5.2 實驗2：方程組解不唯一

運行下面代碼

clear all

x = [1750：25：2000]';

y = 1e6*[791 856 978 1050 1262 1544 1650 2532 6122 8170 11560]';

p = polyfit（x， y， 11）;

即把原來的5次多項式改為11次多項式，命令窗口彈出下述警告：

Warning： Polynomial is not unique; degree gt;= number of data points.

即多項式的度大于等于數據點個數。通過查看代碼polyfit.m，我們發現：方程組（16） 的矩陣[R]的列數大于行數，因此觸發了異常情況2。此時方程組為欠定方程組，解不唯一。我們可以考慮：添加更多的不同的擬合數據點，或者減少多項式的階數。比如這里改為2次多項式，一共11個數據點，再運行代碼就不會出現警告了。

5.3 實驗3：X和Y的長度不相等

如果將原來的第1行代碼改為：x = [1750：25：1900]';" 再運行代碼，則會彈出錯誤：

Error using polyfit （line 47）

X and Y vectors must be the same size.

意思是X和Y的大小必須一樣。經檢查，X為7個數據，而Y為11個數據，X和Y的大小不一樣，觸發了異常情況1，導致報錯。這個只需要將X和Y的數據個數統一即可。

6 結論

本文針對學生在做數學實驗過程中的一個看似不起眼、卻又是經常性遇到的問題——Matlab多項式擬合命令polyfit彈出警告，進行了詳盡地探討。從數據擬合、最小二乘法、方程組的穩定性和條件數等數學理論出發，本文闡述了polyfit命令背后的原理，解釋了代碼運行過程中出現警告的原因和解決方案。通過對這一問題的研究和探索，學生不僅解決了遇到的問題，還培養了探索性思維——生活和實際中從不缺少數學，而是缺少發現數學的眼睛，缺少主動撲捉和解決數學的意識和眼光。教師可以以此為契機，鼓勵學生在今后的學習和生活中，不放過任何一個小問題，要有創新意識，一步一步通過探索，努力找到問題背后的根源，相信問題的解決能給他們帶來相當的成就感，這也是培養創新性人才的核心能力之一。

參考文獻：

[1] 劉浩，韓晶.MATLAB R2022a完全自學一本通[M].北京：電子工業出版社，2022.

[2] 胡曉冬，董辰輝.MATLAB從入門到精通[M].2版.北京：人民郵電出版社，2018.

[3] KOPROWSKI R.Image Analysis for Ophthalmological Diagnosis[M].Cham：Springer International Publishing，2016

[4] 黃平，劉小蘭，溫旭輝.數學實驗典型案例[M].北京：高等教育出版社，2015.

[5] 朱明星.旋轉備用補償模型在烏江渡水電站的運用[J].紅河水， 2023， 42（5）： 133-139.

[6] 陸人定.多信息化技術融合教學探索：以車載紅外測距傳感器GP2D12數據曲線擬合為例[J].時代汽車， 2024：88-91.

[7] 韓一鳴.基于數據擬合和統計學分析的抽油機井單耗影響要素探討[J].石油工業技術監督， 2024， 40（11）： 34-39.

[8] 楊明，程時棟，劉洋，等.基于真實世界數據的大型影像設備資源配置模型的構建與分析[J].中國醫療設備. 2024，39（10）： 120-125.

[9] 王萼芳，石生明.高等代數輔導與習題解答：北大·第五版[M].北京：高等教育出版社，2019.

[10] CHRISTOPHER BISHOP. Pattern Recognition and Machine Learning[M]. Berlin： Springer， 2007. 137-146.

[11] 蔣金山，何春雄，潘少華.最優化計算方法[M].廣州：華南理工大學出版社，2007.

[12] 陳寶林.最優化理論與算法[M].2版.北京：清華大學出版社，2005.

[13] 徐萃薇，孫繩武.計算方法引論[M].3版.北京：高等教育出版社，2007.

[14] 中國計量大學. 載荷識別方法、裝置、電子設備及存儲介質[P]. 中國：CN202310099084. 8， 2023-01-13.

[15] 桂現才，刁群.解線性方程組的MATLAB法[J].湛江師范學院學報，2003，24（6）：14-18.

[16] 徐臨燕.納米梁諧振器的測試與表征方法研究[D].天津：天津大學，2009.

[17] 楊林.基于子空間的SIMO系統盲均衡算法研究[D].鄭州：鄭州大學，2011.

[18] 孫敦超.GPS/GLONASS組合應用相關問題研究[D].成都：電子科技大學，2006.

[19] 劉建國，黃華.最佳平方逼近多項式病態改善[J].南昌航空工業學院學報（自然科學版），2003，17（1）：7-9，17.

【通聯編輯：李雅琪】