Numerical Methods in Engineering (ENGR20005)
Lecture 08 Dr. Leon Chan
lzhchan@unimelb.edu.au Department of Mechanical Engineering The University of Melbourne
Slides prepared by Prof.Andrew Ooi
L8.1: Least-Squares Regression
2
Lecture Notes (Chap. 4, pg. 67)
Assume that you have a car undergoing a constant acceleration. If all the sensors are perfect and you can read the speedometer extremely accurately, then plotting the velocity versus time should look something the graph below and drawing a straight line through all the data points will give a good idea of the acceleration (slope of straight line).
v(t)
t
In reality, there is always noise in the data, readings from sensors have errors and you might not be able to read the speedometer super accurately.
Fitting a curvy line through all the (noisy) data points does not give you good idea about the acceleration.
You can try to draw a straight line through all the (noisy) data points, but which one would you draw? Which is the line that best represents the actual speed of the car?
3
Least-Squares Regression
Suppose you get a set of x, y values from an experiment or any other
sources.
y(x)
s
r
p
q
How can you draw a line that best fits through the data ?
Which of lines p, q, r or s best fit data?
A method to systematically calculating the line of best fit is called Least Squares regression.
4
x
y(x)
e3
e5 e4
y = ax + b e6
Linear Regression
e2
Let’s look straight lines first. They can be represent by
e1
y = ax + b
Task: find a and b such that line predicted by Eq (8.1) is as close to
actual values of yi(xi).
The error between the actual and predicted values are
ei= yi − (axi+b)
x
(8.1)
⏟
actualvalue 5 predictedvalue
Linear Least-Squares Regression
With N data points, can square error at every point and sum all up to find total error.
S = e2 + e2 + e2 + e2 + . . . . . . . + e2 1234 N
= ∑N e 2
i i=1
S=∑N (yi−axi−b)2 i=1
Now find a and b, such that S is minimised. How can we minimize S?
Compute ∂S/∂a, ∂S/∂b, set derivatives to zero.
Here xi, yi constants (data points) when differentiating. 6
Linear Least Squares Approximations
S=∑N (yi−axi−b)2 i=1
S = (y1 − ax1 − b)2
+(y2 −ax2 −b)2
+(y3 −ax3 −b)2 +…. 2 +(yN −axN −b)
∂S ∂a
∂S ∂b
=∑N 2(yi−axi−b)(−xi)=0 i=1
=∑N 2(yi−axi−b)(−1)=0 i=1
(8.2)
∂S =2(y1−ax1−b)(−x1)
∂a +2(y2−ax2−b)(−x2) +2(y3 − a x3 − b)(−x3)
+….
+2(yN −axN −b)(−xN)=0
∂S =2(y1−ax1−b)(−1)
∂b +2(y2−ax2−b)(−1) +2(y3 −ax3 −b)(−1)
+….
+2(yN −axN −b)(−1)=0
Eq. (8.2) can be re-written as
a ∑N x 2 + b ∑N x i = ∑N x i y i i=1 i i=1 i=1
a ∑N x i + b N = ∑N y i i=1 i=1
a ( x 12 + x 2 2 + x 32 + . . . . + x N2 ) + b ( x 1 + x 2 + x 3 + . . . . + x N ) = ( x 1 y 1 + x 2 y 2 + x 3 y 3 + . . . + x N y N )
a(x1 + x2 + x3 + …. + xN) + b(1 + 1 + 1 + …. + 1) = (y1 + y2 + y3 + … + yN) 7
Linear Least-Squares Approximations
a+b= a ∑N x i + b N =
i=1
This is a system of 2 x 2 linear equations. They can be solved to
get a and b.
∑N xi
i=1 N
∑N x 2
∑N x i i=1
i i=1
∑N x 2
i=1
i
∑N xi i=1
∑N x i y i i=1
∑N yi i=1
[ab]
[A]{X} = {C} 8
=
∑N x i y i i=1
∑N y i i=1
Example L08.1: Fit a straight line through the data below using linear least-squares
RPM(x)
Diesel engine HP (y)
400
580
500
1030
600
1420
700
1880
800
2100
9
Lecture08M01.m
xdata=[400 500 600 700 800] ydata=[580 1030 1420 1880 2100]
A=[sum(xdata.^2) sum(xdata); sum(xdata) length(xdata)];
∑Nxi2 ∑Nxi i=1 i=1
∑Nxi N i=1
C=[sum(xdata.*ydata); sum(ydata)];
vec=inv(A)*C; a=vec(1); b=vec(2);
x=linspace(350,850); y=a*x+b;
[ ab ] =
∑N x i y i i=1
∑N yi i=1
hold off plot(xdata,ydata,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’r’) hold on
plot(x,y,’b-‘,’LineWidth’,2);
xlabel(‘RPM’);
ylabel(‘HP’);
10
L8.2:
Quadratic Least-Squares Regression
11
Lecture Notes (Chap. 4, pg. 67)
Polynomial Regression
Now imagine you get following data from an experiment:
y(x)
s
p
q
r
Clearly trying to plot a straight line through these data is not a wise thing to do.
It is much more appropriate to try and fit a polynomial though the data above. This is called Polynomial Least-Squares regression. But which polynomial?
12
x
Quadratic Least-Squares Regression Example
Assume we have N data pairs(x,y).
We want to use the polynomial y = a0 + a1 x + a2x2. So we want to
find a0, a1 and a2.
The error between the actual data and the polynomial is then defined
as
Square errors and sum them
ei = yi − (a0 + a1xi + a2x2) ⏟i
actual value
S=∑N e2 =∑N (yi−a0−a1xi−a2x2)2
i=1i i=1 i 13
predicted value
S=∑N e2 =∑N (yi−a0−a1xi−a2x2)2 i=1i i=1 i
S = (y1 − a0 − a1x1 − a2x2)2 + y2 −a0 −a1x2 −a2x2
( 1)2 + y3 −a0 −a1x3 −a2x2
( 2)2 +….. 3
+(yN − a0 − a1xN − a2x2)2 N
14
S = ( y 1 − a 0 − a 1 x 1 − a 2 x 12 ) 2
+(y2 − a0 − a1x2 − a2x2)2
+(y3 − a0 − a1x3 − a2x32)2 +…..
+(yN − a0 − a1xN − a2xN2)2
∂S =2(y1−a0−a1x1−a2x12)(−1)
∂S =2(y1−a0−a1x1−a2x12)(−x1)
+2 (y2 − a0 − a1x2 − a2x2)(−x2) +2 (y3 − a0 − a1x3 − a2x32)(−x3)
∂S =2(y1−a0−a1x1−a2x12)(−x12)
+2 (y2 − a0 − a1x2 − a2x2)(−x2) +2 (y3 − a0 − a1x3 − a2x32)(−x32)
∂a0
+2 (y2 − a0 − a1x2 − a2x2)(−1)
∂a1
∂a2
+….. +…..
+2 (y3 − a0 − a1x3 − a2x32)(−1) +…..
+ 2 ( y N − a 0 − a 1 x N − a 2 x N2 ) ( − 1 ) = 0
∑N y i − ∑N a 0 − a 1 ∑N x i − a 2 ∑N x i 2 = 0 i=1 i=1 i=1 i=1
a 0 ∑N + a 1 ∑N x i + a 2 ∑N x i 2 = ∑N y i i=1 i=1 i=1 i=1
+ 2 ( y N − a 0 − a 1 x N − a 2 x N2 ) ( − x N ) = 0 ∑N x i y i − ∑N a 0 x i − a 1 ∑N x i 2 − a 2 ∑N x i 3 = 0
i=1 i=1 i=1 i=1
a 0 ∑N x i + a 1 ∑N x i 2 + a 2 ∑N x i 3 = ∑N x i y i
+ 2 ( y N − a 0 − a 1 x N − a 2 x N2 ) ( − x N2 ) = 0 ∑N x i 2 y i − ∑N a 0 x i 2 − a 1 ∑N x i 3 − a 2 ∑N x i 4 = 0
i=1 i=1 i=1 i=1
a 0 ∑N x i 2 + a 1 ∑N x i 3 + a 2 ∑N x i 4 = ∑N x i 2 y i
i=1 i=1 i=1 i=1
i=1 i=1
i=1 i=1
15
Write in matrix form:
a 0 ∑N + a 1 ∑N x i + a 2 ∑N x i 2 = ∑N y i i=1 i=1 i=1 i=1
a 0 ∑N x i + a 1 ∑N x i 2 + a 2 ∑N x i 3 = ∑N x i y i i=1 i=1 i=1 i=1
a 0 ∑N x i 2 + a 1 ∑N x i 3 + a 2 ∑N x i 4 = ∑N x i 2 y i i=1 i=1 i=1 i=1
N ∑xi ∑xi2 a0 ∑yi ∑xi ∑xi2 ∑xi3 a1 = ∑xiyi ∑ x i2 ∑ x i3 ∑ x i4 a 2 ∑ x i2 y i
Solve this system of linear equations using, e.g., Gauss elimination or LU decomposition to give the ai’s.
16
Example L08.2: Fit a quadratic line through the data below using the least-squares regression.
x
y
0.075
600
0.5
800
1.0
1200
1.2
1400
1.7
2050
2.0
2650
2.3
3750
17
Lecture08M02.m
xdata=[0.075 0.5 1.0 1.2 1.7 2.0 2.3] ydata=[600 800 1200 1400 2050 2650 3750]
N ∑xi ∑xi2 ∑xi ∑xi2 ∑xi3 ∑ xi2 ∑ xi3 ∑ xi4
N=length(xdata);
A=[N sum(xdata) sum(xdata.^2) ;
sum(xdata) sum(xdata.^2) sum(xdata.^3); sum(xdata.^2) sum(xdata.^3) sum(xdata.^4)];
C=[sum(ydata); sum(xdata.*ydata); sum(xdata.^2.*ydata)];
%
%Note that for brevity I have used the MATLAB linsolve() function %solve the resulting equation. You can use any of the
%solvers from previous lectures.
%
a=linsolve(A,C);
y = a0 + a1x + a2x2
x=linspace(0,2.5); y=a(1)+a(2)*x+a(3)*x.^2;
hold off plot(xdata,ydata,’ko’,’MarkerSize’,20,’MarkerFaceColor’,’r’) hold on
plot(x,y,’b-‘,’LineWidth’,2);
xlabel(‘x’);
ylabel(‘y’);
a0 a1 a2
=
∑yi ∑xiyi ∑xi2yi
Note that the indexing in MATLAB
always start with 1 (not 0). So
a(1) correspond to a0, a(2) to a1 and a(3) to a2
18
End of Example L08.2
19
L8.3:
Polynomial Least-Squares Regression
20
To fit a quadratic polynomial
y = a0 + a1x + a2x2 Define error sum of squares
S=∑N ei2 =∑N (yi−a0−a1xi−a2xi2)2 i=1 i=1
Set partial derivatives with respects to a0, a1 and a2 to zero
∂S =∑N 2(yi−a0−a1xi−a2xi2)(−1)=0, ∂a0 i=1
∂S =∑N 2(yi−a0−a1xi−a2xi2)(−xi)=0, ∂a1 i=1
∂S =∑N 2(yi−a0−a1xi−a2xi2)(−xi2)=0, ∂a2 i=1
Put resulting equations in matrix form and solve for a0, a1 and a2 .
2 N∑x∑xa∑y
iia0 i ∑xi ∑xi2 ∑xi3 a1 = ∑xiyi
∑xi2 ∑xi3 ∑xi4 2 ∑xi2yi
To fit a n’th order polynomial
y = a0 + a1x + a2x2 + . . . . . . + anxn Define error sum of squares
S=∑N ei2 =∑N (yi−a0−a1xi−a2xi2−……−anxin)2 i=1 i=1
Set partial derivatives with respects to a0, a1, a2, a3, a4…. an to zero
∂S =∑N 2(yi−a0−a1xi−a2xi2−……−anxin)(−1)=0, ∂a0 i=1
∂S =∑N 2(yi−a0−a1xi−a2xi2−……−anxin)(−xi)=0, ∂a1 i=1
∂S =∑N 2(yi−a0−a1xi−a2xi2−……−anxin)(−xi2)=0, ∂a2 i=1
. N
∂S =∑2(yi−a0−a1xi−a2xi2−……−anxin)(−xin)=0,
∂an i=1
Put resulting equations in matrix form and solve for
a0, a1, a2, a3, a4…. an.
N ∑xi ∑xi2 ∑xi3 ⋯ ∑xin a0 ∑yi
21
∑ xn+1 ∑ x y i i i i i a1 ii
∑ x ∑ x2 ∑ x3 ∑ x4 ⋯
∑x2 ∑x3 ∑x4 ∑x5 ⋯ ∑xn+2 a2 ∑x2y
iiiiia=ii ∑x3 ∑x4 ∑x5 ∑x6 ⋯ ∑xn+3 3 ∑x3y
iiiii⋮ii ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ an ⋮n
∑xn ∑xn+1 ∑xn+2 iiiii
∑xn+3 ⋯ ∑x2n ∑xi yi
Polynomial Least-Squares Regression
Let’s now assume a general polynomial of order n
y = a0 + a1x + a2x2 + . . . . . . + anxn
The error between the actual data and the polynomial is again
defined as
ei = yi −(a0 +a1xi +a2x2 +……+anxn) ⏟ii
actual value
predicted value
As before, square errors and sum them
S=∑N e2 =∑N (yi−a0−a1xi−a2x2−……−anxn)2 i=1ii=1 i i
This is the total error, that we want to minimise.
22
Polynomial Least Squares Regression
How can we minimise S?
⇒ Compute dS/dai, set derivatives to zero.
Again, xi, yi are constants (data points) when differentiating.
∂S ∂a1
∂S ∂a2
. ∂S
∂an
∂S ∂a0
=∑N 2(y−a −ax−ax2−……−axn)(−1)=0,
i01i2i n i=1
i
=∑N 2(y−a −ax−ax2−……−axn)(−x)=0,
i01i2i nii i=1
=∑N 2(y−a −ax−ax2−……−axn)(−x2)=0, i=1 i 0 1i 2i ni i
=∑N 2(y−a −ax−ax2−……−axn)(−xn)=0,
i01i2i n i=1
ii
23
Polynomial Least-Squares Regression
How can we minimise S?
⇒ Compute dS/dai, set derivatives to zero.
Again, xi, yi are constants (data points) when differentiating. ∑N (a +ax+ax2+……+axn)=y,
Ni=101i2i ni i ∑ (a x + a x2 + a x3 + . . . . . . + a xn+1) = y x ,
0i1i2iniii i=1
∑N (ax2+ax3+ax4+……+axn+2)=yx2, i=10i1i2i niii
.
∑N (a0xn+a1xn+1+a2xn+2+……+anx2n)=yixn,
i=1iii ii 24
Polynomial Least-Squares Regression
How can we minimise S?
⇒ Compute dS/dai, set derivatives to zero.
Again, xi, yi are constants (data points) when differentiating. aN+a∑N x+a∑N x2+……+a∑N xn=y,
0 1i=1 i 2i=1 i ni=1 i i a∑N x+a∑N x2+a∑N x3+……+a∑N xn+1=yx,
0 i=1 i 1 i=1 i 2 i=1 i n i=1 i i i
a ∑N x2+a ∑N x3+a ∑N x4+……+a ∑N xn+2 =yx2, 0 i=1 i 1 i=1 i 2 i=1 i n i=1 i i i
.
a0∑N xn+a1∑N xn+1+a2∑N xn+2+……+an∑N x2n =yixn,
i=1 i i=1 i i=1 i i=1 i i 25
Polynomial Least-Squares Regression
Write in matrix form:
N ∑xi ∑xi2 ∑xi3 ⋯ ∑xin a0 ∑yi ∑x ∑x2 ∑x3 ∑x4 ⋯ ∑xn+1 ∑xy
i i i i i a1 ii ∑x2 ∑x3 ∑x4 ∑x5 ⋯ ∑xn+2 a2 ∑x2y
iiiiia=ii ∑x3 ∑x4 ∑x5 ∑x6 ⋯ ∑xn+3 3 ∑xi3yi
iiiii⋮⋮ ⋮n ⋮n+1 ⋮n+2 ⋮n+3 ⋮ ⋮2n an ∑xny
∑xi ∑xi ∑xi ∑xi ⋯ ∑xi i i [XX] {A} = {XY}
or
This system of linear equations can be solved with tools such as
Gauss elimination or LU decomposition to give the ai’s.
Note: here the solution vector is {A} and {X} contains the known data points that fill matrix[XX].
26
Polynomial Least Squares Approximation
Another Example
27
Polynomial Least Squares Approximation
28
Polynomial Least Squares Approximation
29
Polynomial Least Squares Approximation
30
Polynomial Least Squares Approximation
31
Polynomial Least Squares Approximation
32
Example L08.3: Fit a 5th, 6th and 7th order polynomial through the data points in Example L08.2 and check your answer with the polyfit() function in MATLAB.
x
y
0.075
600
0.5
800
1.0
1200
1.2
1400
1.7
2050
2.0
2650
2.3
3750
33
L8.4:
Linearization of Nonlinear Relationships
34
Other Least Squares Approximations
So far only considered polynomial least squares approximations (linear LSQ is a subset).
Occasionally, it might be appropriate to use other `basis’ functions, e.g.
• •
Exponential: y=b eax Power laws: y=b xa
35
Other Least Squares Approximations
When using exponential basis, we get the error between the actual data and the approximation
ei = yi − beaxi ⏟⏟
actual value
As before, square errors and sum them
S=∑N e2=∑N (yi−beaxi)2 i=1 i i=1
This is the total error, that we want to minimise.
36
predicted value
Other Least Squares Approximations
How can we minimise S?
⇒ Compute dS/dai, set derivatives to zero.
Again, xi, yi are constants (data points) when differentiating. ∂S = ∑N 2 (yi − beaxi) (−bxieaxi) = 0,
∂a i=1N
∂S =∑2(yi−beaxi)(−eaxi)=0,
∂b i=1
This system of equations is difficult to solve.
37
Other Least Squares Approximations
When using power law basis, we get the error between the actual data and the approximation
e= y − bxa
i ⏟i actual value
⏟i predicted value
As before, square errors and sum them
S=∑N e2=∑N (yi−bxa)2 i=1ii=1 i
This is the total error, that we want to minimise.
38
Other Least Squares Approximations
How can we minimise S?
⇒ Compute dS/dai, set derivatives to zero.
Again, xi, yi are constants (data points) when differentiating.
=∑N 2(yi−bxa)(−blnxixa)=0, i=1 i i
∂S=∑N 2(yi−bxa)(−xa)=0, ∂bi=1 ii
Solution to this set of equations is difficult to compute.
∂S ∂a
39
Other Least-Squares Approximations
The method commonly used when data suspected to be exponentially related is to consider logarithm of approximation equation
Replace original equations with
y= +b Or if you suspect there is an inverse x i.e. 1/x relationship then
a
x
y = perx ⇒ ln y = rx + ln p
y = pxr ⇒ ln y = r ln x + ln p
y=p x ⇒(1)=(r)(1)+ r+xypx
A linear problem now appears and solutions can be found as in the linear least squares approach (y = ax + b). So you can use the
techniques described in L8.1 to fit the functions.
However, approximation obtained this way is NOT the least squares
approximation of the original problem.
1
p
40
y = pxr
log(y) = r log(x) + log(p)
r
1
y=p x
r+x ypxp
(1)=(r)(1)+1 p
r
41
Example L08.4: This example is taken from “MATLAB An introduction with applications” by Amos Gilat. Fit the data below to the exponential function y = pemx
x
y
0
6.00
0.5
4.83
1.0
3.70
1.5
3.15
2.0
2.41
2.5
1.83
3.0
1.49
3.5
1.21
4.0
0.96
4.5
0.73
5.0
0.64
42
We want to fit the data to this function
y = pemx
In Example L08.1, you were shown how to fit a set of data points to the function
y = ax + b
The idea is to make Eq. (8.3) look as much as possible to Eq. (8.4).
If you take the natural log both sides of Eq. (8.3), you should get
(8.3)
(8.4)
ln(y) = ln(p) + mx. ln(y) = mx + ln(p)
or
Eq. (8.5) very similar to the linear function in Eq. (8.4). If you let
(8.5)
y → ln(y) x→x a→m b → ln p
So the idea is to take natural log of the y data in the table in the previous slide and use that as the y data for the MATLAB program Lecture08M01.m. The x data remains unchanged.
43
y = pemx
Lecture08M04.m
x=[0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0]
y=[6.00 4.83 3.70 3.15 2.41 1.83 1.49 1.21 0.96 0.73 0.64]
xdata=x ydata=log(y)
ln(y) = mx + ln(p) y = ax + b
A=[sum(xdata.^2) sum(xdata); sum(xdata) length(xdata)];
C=[sum(xdata.*ydata); sum(ydata)];
vec=inv(A)*C;
a=vec(1); %this should be m b=vec(2); % this should be ln(p) a=m;
p=exp(b);
xplot=linspace(0,5); yplot=p*exp(m*xplot);
y = pemx
m=a ln(b) = p
p = eb
hold off plot(x,y,’ko’,’MarkerSize’,15,’MarkerFaceColor’,’g’) hold on
plot(xplot,yplot,’k-‘,’LineWidth’,2)
xlabel(‘x’);
ylabel(‘y’);
44
Example L08.3/L08.4
y
log(y)
x
x
End of Example L08.4
46
Temperature (K)
Viscosity (Ns/m2)
253.15
4
273.15
0.38
293.15
0.095
313.15
0.032
333.15
0.015
353.15
0.0078
373.15
0.0045
393.15
0.0032
Temperature (K)
Viscosity (Ns/m2)
253.15
4
273.15
0.38
293.15
0.095
313.15
0.032
333.15
0.015
353.15
0.0078
373.15
0.0045
393.15
0.0032
2nd order
2nd order 6th order
8th order
10th order
Temperature (K)
Viscosity (Ns/m2)
253.15
4
273.15
0.38
293.15
0.095
313.15
0.032
333.15
0.015
353.15
0.0078
373.15
0.0045
393.15
0.0032
2nd order
5th order
2nd order
μ = e0.0003T2−0.2686T+47.2091
5th order