java 多项式拟合最多的项数_Matlab polyfit 详解 | 方程组求解的稳定性 | 条件数

• 本文是对 Matlab 中 polyfit 函数进行原理解析，并没介绍该如何具体使用 polyfit 函数，具体方法请自行查阅官方文档。
• 主要涉及数值分析的相关内容。
• 简单介绍了数据标准化（Z-score 标准化）、QR 分解、Matlab 中的求逆运算符。
• 相对详细介绍了线性方程组求解的稳定性问题，并引出条件数的定义。
• 最后根据 polyfit 的源码对它进行计算流程解析，并分析相关的警告该如何处理。

预备知识

标准化

对向量

进行标准化（Standardization，Z-score normalization）定义如下：

其中

表示

的均值，即

，

表示标准差，即

。

• 为
  标准分数（standard score）也叫 Z分数（Z-score）；
• Z分数表示原始数据与平均值的差，单位为标准差，也就时给定值距离均值多少个标准差。

---

QR 分解定理

设

是

阶实非奇异矩阵，则存在正交矩阵

和上三角矩阵

使得

，当

的对角元均为正时，分解唯一。

注：正交矩阵

---

求逆运算

Y = inv(X) 计算方阵

的逆矩阵。

• X^(-1) 等效于 inv(X)。
• x = Ab 的计算方式与 x = inv(A)*b 不同，建议用于求解线性方程组。

---

方程组求解的稳定性

设线性方程组

，其中

有小扰动

，得到解

，那么：

由

，可得

考虑相对误差：

即

条件数

• 称
  为
  的条件数；
• 表示
  变化时解的相对误差灵敏度的度量；
• 当
  较小时，解对
  的扰动不敏感。

注：

表示矩阵范数，满足正定性、齐次性以及三角不等式。常见的范数定义有：

---

polyfit 函数原理

polyfit 使用一维向量

构造具有

列和

m = length(x) 行的范德蒙（Vandermonde）矩阵

（见源码55-59行）并生成线性方程组：

拟合多项式

转换为求解线性方程组(2)（非上述方程组），具体地，

polyfit 先对

进行

QR分解，然后再使用 反斜杠 运算进行求解（见源码61-69行）。

具体过程如下：

转化为 Matlab 运算即为

p = R(Q'*y)。

异常处理逻辑

共以下4种情况

• MATLAB:polyfit:XYSizeMismatch
• MATLAB:polyfit:PolyNotUnique
• MATLAB:polyfit:RepeatedPoints
• MATLAB:polyfit:RepeatedPointsOrRescale

输入判断

• polyfit函数只用于拟合(x,y)型二维数据，因此输入数据应是两两对应的，即向量
  长度相等（见源码43-45行）。
• 若提示MATLAB:polyfit:XYSizeMismatch，请检查上述问题。

欠定方程

• 若提示MATLAB:polyfit:PolyNotUnique 表示线性方程组
  欠定（矩阵
  的列数大于行数）

考虑

1. 添加更多的不同的拟合点;
2. 减少多项式的次数.

条件数过大

由于 Vandermonde 矩阵中的列是向量

的幂，因此

的

条件数对于高阶拟合来说通常 较大，生成一个 奇异系数矩阵。在这些情况下，中心化和缩放可改善系统的数值属性以产生更可靠的拟合。而数据的标准化由函数调用时输出情况决定，当且仅当调用 [P,S,MU] = POLYFIT(X,Y,N) 才执行数据标准化（见源码50-53行）。

当

过大时才会提示以下其中一个警告：

• 若提示 MATLAB:polyfit:RepeatedPoints，考虑添加更多的不同的拟合点
• 若提示 MATLAB:polyfit:RepeatedPointsOrRescale，考虑1、添加更多的不同的拟合点；2、减少多项式的次数

源码阅读

function [p,S,mu] = polyfit(x,y,n)
%POLYFIT Fit polynomial to data.
%   P = POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
%   degree N that fits the data Y best in a least-squares sense. P is a
%   row vector of length N+1 containing the polynomial coefficients in
%   descending powers, P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).
%
%   [P,S] = POLYFIT(X,Y,N) returns the polynomial coefficients P and a
%   structure S for use with POLYVAL to obtain error estimates for
%   predictions.  S contains fields for the triangular factor (R) from a QR
%   decomposition of the Vandermonde matrix of X, the degrees of freedom
%   (df), and the norm of the residuals (normr).  If the data Y are random,
%   an estimate of the covariance matrix of P is (Rinv*Rinv')*normr^2/df,
%   where Rinv is the inverse of R.
%
%   [P,S,MU] = POLYFIT(X,Y,N) finds the coefficients of a polynomial in
%   XHAT = (X-MU(1))/MU(2) where MU(1) = MEAN(X) and MU(2) = STD(X). This
%   centering and scaling transformation improves the numerical properties
%   of both the polynomial and the fitting algorithm.
%
%   Warning messages result if N is >= length(X), if X has repeated, or
%   nearly repeated, points, or if X might need centering and scaling.
%
%   Example: simple linear regression with polyfit
%
%     % Fit a polynomial p of degree 1 to the (x,y) data:
%       x = 1:50;
%       y = -0.3*x + 2*randn(1,50);
%       p = polyfit(x,y,1);
%
%     % Evaluate the fitted polynomial p and plot:
%       f = polyval(p,x);
%       plot(x,y,'o',x,f,'-')
%       legend('data','linear fit')
%
%   Class support for inputs X,Y:
%      float: double, single
%
%   See also POLY, POLYVAL, ROOTS, LSCOV.

%   Copyright 1984-2018 The MathWorks, Inc.

if ~isequal(size(x),size(y))
    error(message('MATLAB:polyfit:XYSizeMismatch'))
end

x = x(:);
y = y(:);

if nargout > 2
    mu = [mean(x); std(x)];
    x = (x - mu(1))/mu(2);
end

% Construct the Vandermonde matrix 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

% Solve least squares problem p = Vy to get polynomial coefficients p.
[Q,R] = qr(V,0);
oldws = warning('off','all');   % Turn all warnings off before solving
try
    p = R(Q'*y);               % Same as p = Vy
catch ME
    warning(oldws);             % Restore initial warning state
    throw(ME);
end
warning(oldws);                 % Restore initial warning state

% Issue warnings.
if size(R,2) > size(R,1)
    warning(message('MATLAB:polyfit:PolyNotUnique'))
elseif warnIfLargeConditionNumber(R)
    if nargout > 2
        warning(message('MATLAB:polyfit:RepeatedPoints'));
    else
        warning(message('MATLAB:polyfit:RepeatedPointsOrRescale'));
    end
end

if nargout > 1
    r = y - V*p;
    % S is a structure containing three elements: the triangular factor
    % from a QR decomposition of the Vandermonde matrix, the degrees of
    % freedom and the norm of the residuals.
    S.R = R;
    S.df = max(0,length(y) - (n+1));
    S.normr = norm(r);
end

p = p.'; % Polynomial coefficients are row vectors by convention.

function flag = warnIfLargeConditionNumber(R)
if isa(R, 'single')
    flag = (condest(R) > 1e+05);
else
    flag = (condest(R) > 1e+10);
end

参考资料

1. Matlab inv
2. Matlab @polyfit
3. Stoer J, Bulirsch R. Introduction to Numerical Analysis[M]. Third Edition. New York, NY: Springer, 2002.
4. 我的原始语雀文档