polyfit

Basic Idea

In the following discussion, $\mathbb{K}$ means either the real line $\mathbb{R}$ or the complex plane $\mathbb{C}$. Given $m$ data points $(x_i,y_i)$ for $i=1,2,\ldots,m$. The polynomial fitting problem is to find an $n$-th order polynomial

to minimise the sum of squared errors (SSE)

where $A\in\mathbb{R}^{m\times (n+1)}$ is the Vandermonde matrix

and the vectors

To fit an $n$-th order polynomial to $m$ data points, we solve the so-called least-square problem:

In this app, the least-square problem is solved by QR decomposition, as we will show in the following.

Solving the Least-square Problem

For convenience, let $\bar{n}=n+1$ be the length of $p$. Using QR decomposition with column pivoting, we get

where $A\in K^{m\times\bar{n}}$, $P$ is the $\bar{n}\times \bar{n}$ permutation matrix, $Q\in\mathbb{K}^{m\times m}$ and $R\in\mathbb{K}^{m\times \bar{n}}$. Let $k=\min\left\{\bar{n},m,\text{rank}(A)\right\}$. Then,

where $R_1\in\mathbb{K}^{k\times k}$ is upper triangular and $R_2\in\mathbb{K}^{k\times (\bar{n}-k)}$. For any given $p\in\mathbb{K}^{\bar{n}}$, we define

where $z_1\in\mathbb{K}^{k}$ and $z_2\in\mathbb{K}^{\bar{n}-k}$. As a result,

We divide $Q^{-1}b$ into two vectors,

where $c\in\mathbb{R}^{k}$ and $d\in\mathbb{K}^{m-k}$. Then,

This can be minimized by selecting $z_1$ and $z_2$ properly. Then, $p=Pz$ is a solution to the least-square problem. The selection of $z_1$ and $z_2$ depends on the size and rank of $A$, which will be discussed below.

Case 1: Full Rank Overdetermined system

In this case, ${\rm rank}(A)=\bar{n}\leq m$. Therefore, $k=\min\left\{\bar{n},m,\text{rank}(A)\right\}=\bar{n}$. Hence, $R_2$ and $z_2$ are empty and $d\in\mathbb{K}^{m-\bar{n}}$. Then, $\|Ap-b\|_2^2=\|R_1z_1-c\|_2^2+\|d\|_2^2$ is minimized if $z=z_1$ is the solution of the upper triangular system $R_1z_1=c$. In the special case where $m=\bar{n}$, the vector $d$ is also empty and hence $\|Ap-b\|^2_2=0$. That means, the polynomial will pass through of the data points.

Case 2: Rank Deficient Overdetermined system

For convenience, let $r={\rm rank}(A)$. In this case, $r<\bar{n}\leq m$. Therefore, $k=\min\left\{\bar{n},m,r\right\}=r$. Hence, $R_2$, $z_2$ and $d$ are all non-empty. Then, $\|Ap-b\|_2^2=\|R_1z_1+R_2z_2-c\|_2^2 +\|d\|_2^2$ is minimized when $z_2=0$ and $z_1$ is the solution of the upper triangular system $R_1z_1=c$. There can be other choices of $z_1$ and $z_2$ that minimise $\|Ap-b\|^2_2$. Hence,

is one of the many solutions to the least-square problem.

Case 3: Full Rank Underdetermined system

In this case, $\bar{n}> m={\rm rank}(A)$. Therefore, $k=\min\left\{\bar{n},m,\text{rank}(A)\right\}=m$. Hence, $R_2$ and $z_2$ are non-empty but $d$ is empty. Then, $\|Ap-b\|_2^2=\|R_1z_1+R_2z_2-c\|_2^2=0$ if $z_2=0$ and $z_1$ is the solution of the upper triangular system $R_1z_1=c$. Hence, the polynomial will pass through all of the data points. There can be other choices of $z_1$ and $z_2$ that minimise $\|Ap-b\|^2_2$. Hence,

is one of the many solutions to the least-square problem.

Case 4: Rank Deficient Underdetermined system

For convenience, let $r={\rm rank}(A)$. In this case, $\bar{n}> m>r$. Therefore, $k=\min\left\{\bar{n},m,r\right\}=r$. Hence, $R_2$, $z_2$ and $d$ are all non-empty. Then, $\|Ap-b\|^2=\|R_1z_1+R_2z_2-c\|^2+\|d\|^2$ is minimized if $z_2=0$ and $z_1$ is the solution of the triangular system $R_1z_1=c$. There can be other choices of $z_1$ and $z_2$ that minimize $\|Ap-b\|^2_2$. Hence,

is one of the many solutions to the least-square problem.

Summary

Usages

p = polyfit(x,y,n)

• x and y should be arrays that contain the $x$- and $y$-coordinates of given data points.
• The $x$- and $y$-coordinates should be neither Inf nor NaN.
• If x or y is not a vector, it will be reshaped to a vector.
• Both x and y should be non-empty and have the same lengths.
• n should be a real non-negative integer, representing the order of the polynomial.
• The output p contains coefficients of the polynomial that is a best fit of the given data points. That means, p minimises the sum of squared errors (SSE) discussed here. The first element of p corresponds to the highest order term.
• If x or y is complex, p will be complex.

[p,S] = polyfit(x,y,n)

• x, y, n and p are as described in p = polyfit(x,y,n) above.

• S is a structure that can be used in polyval to evaluate standard errors of prediction. S should contain the following fields:

  • S.normr gives the error $\|Ap-b\|_2$. See here for the meanings of $A$, $p$ and $b$.
  • S.R gives the first min(numel(x), n + 1) rows of $RP'$. See here for meanings of $R$ and $P'$.
  • S.df gives the degree of freedom, i.e. max(0, numel(x) - (n + 1)).

[p,S,mu] = polyfit(x,y,n)

• x, y, n and p are as described in p = polyfit(x,y,n) above.
• S is as described in [p,S] = polyfit(x,y,n) above.
• mu is a vector of 2 elements where mu(1) and mu(2) are the mean and standard deviation of x, respectively. It fits a polynomial to data points with $x$- and $y$-coordinates given by (x-mu(1)) / mu(2) and y, respectively.
• If x has only one element, mu is [x, 0] and the given data point will not be transformed.

Examples

• Fitting a 2nd order polynomial to data points $(x,x^2+\epsilon)$, where $\epsilon\sim N(0,0.2)$. The polynomial returned by polyfit is $p(x)=0.99945x^2 -0.00455x+ 0.01359$.

Input

clear
% Data points.
x=linspace(0,2);
y=x.^2+normrnd(0,0.2,[1 length(x)]);
% Get polynomial of degree 2 that fits the data.
p=polyfit(x,y,2)
scatter(x,y)
hold('on')

% Plot the polynomial
xx=0:0.01:2;
yy=polyval(p,xx);
plot(xx,yy)
hold('off')

Output

All variables cleared.

p = 1e-1 × 
 9.9945  -0.0455   0.1359

• Fitting a polynomial to 5 data points with constant $y$-coordinates.

Input

clear
% Data points.
% y is constant.
x=1:5;
y=rand*ones(1,length(x));
% Get polynomial of degree 3 that fits the data.
% It returns a zero-th degree polynomial because the y-data are constant.
[p,S]=polyfit(x,y,3)
scatter(x,y)
hold('on')

% Plot the polynomial
xx=0:0.01:5;
yy=polyval(p,xx);
plot(xx,yy)
ylim([0 1])
ylabel('y')
xlabel('x')
hold('off')

Output

All variables cleared.

p = 1e-1 × 
-0.0000   0.0000  -0.0000   3.9839

S =
 struct with fields:
    df: [1.0000]
     R: Double [4 x 4]
 normr: [3.3766] × 1e-16