Transcript High Performance Computing 811
Computational Methods in Physics PHYS 3437
Dr Rob Thacker Dept of Astronomy & Physics (MM-301C) [email protected]
Today’s Lecture
Interpolation & Approximation III LU decompositions for matrix algebra Least squares polynomial fitting Fitting to data Fitting to functions – the Hilbert Matrix Use of orthogonal polynomials to simplify the fitting process
The standard approach
For a linear system of n equations with n unknowns Gaussian elimination is the standard technique
a
11
x a
21
x
1 1
a
12
x
2
a
22
x
2 ...
...
a
1
n x n a
2
n x n
b
1
b
2
a n
1
x
1
a n
2
x
2 ...
a nn x n
b n
Or in matrix form :
a
11
a
21
a n
1
a a a
12 22
n
2
a
1
n a
2
n a nn
x x
2
x n
1
b b
2
b n
1 Short A x form b : As we all know G.E. is an effective but time consuming algorithm – applying the algorithm on an n×n matrix takes ~n 3 operations ( O (n 3 )).
While we can’t change the order of the calculation we can reduce the overall number of calculations.
The LU Decomposition method
A non-singular matrix can be written as a product of lower and upper diagonal matrices A=LU where,
L
l l l
11 21
n
1 0
l
22
l n
2 0 0
l nn
U
1 0 0
u
12 1 0
u
1 1
u
2
n n
The simple nature of L & U means that it is actually fairly straightforward to relate the values to A Then can find the l ij and u ij fairly easily
By examination
If A=LU, multiply to get
a
11
a
21
a n
1
a
12
a
22
a n
2
a a
1
n
2
a nn n
l l l
11 21
n
1
l
21
l
11
u
12
u
12
l
22
l n
1
u
12
l n
2
l
21
u
1
n l
11
u
1
n
l
22
u
2
n
i
1
n
1
l ni u in
l nn
Equate first col: Equate first row: Equate second col:
a
11
l
11
a
12
a
21
l
21
a
11
a
12
l
11
u
12
l
11
a
22
l
11
u
12 21
u
12 22
a n
1
l n
1
Now we know all l i1 , i=1,..,n Use l 11 a
1
n
u
1
n l
11
vals to get all u 1j j=2,…,n a n
2
l n
1
u
12
l n Use previous steps to get l j2
2
i=2,…,n
Alternate from columns to rows
Thus expanding each column of A allows us to equate a column of values in L Similarly, expanding a row of A allows to equate to a row of values in U However, most rows will be combinations of U and L entries and we must alternate between columns and rows to ensure each expression contains only 1 unknown
General Formulae (Exercise)
Provided we consider column values on or below the diagonal then we have for the k-th column
a ik
l ik
k j
1 1
l ij u jk
l ik
a ik
k j
1 1
l ij u jk
where
i
k
,
k
1 ,...,
n
Similarly, for row values to the right of the main diagonal then we have for the k-th row
a kj
u kj l kk
1
k m
1
l km u mj
u kj
1
l kk
a kj
1
k m
1
l km u mj
where
j
k
1 ,
k
2 ,...,
n
Convince yourself of these two relationships…
Utilizing the LU decomposition
Since Define
A
x z
b
LU
x
b
U
x
(4) , then w e have
L z
b
(5) We need to solve eqn (5) for z first Once z is known we can then solve eqn (4) for x Sounds like a lot of work! Why is this useful?
The triangular matrices are already in a Gauss-Jordan eliminated form so both can be solved by substitution
L z
b
If we write out the equations:
l
11
z
1
b
1
z
1
b
1 /
l
11
l
21
z
1
l
22
z
2
b
2
z
2 1
l
22
b
2
l
21
z
1
l
n
1
z
1
l n
2
z
2 ...
l nn z n
b n
Which is solved via “forward substitution” – exactly analogous to the back substitution we saw for triadiagonal solvers
z i
l
1
ii
b i
i k
1 1
l ik z k
i
2 ,...,
n
Next solve
U x
z
Again if we right out the equations then:
x
1
u
12
x
2
u
13
x
3 ...
u
1
n x n
z
1
x
2
u
23
x
3 ...
u
2
n x n
z
2
x n
1
u n
1 ,
n x n x n
z n
1
z n
We solve by back substitution again:
x n
z n
then substitute into
x n
1
u n
1 ,
n x n
to get
x n
1
z n
1
u n
1 ,
n x n
and in general
z n
1
x n
i
z n
i
k
n n
i
u
1
n
i
,
k x k
wher e
i
1
,...,n-
1
Least squares fit of a polynomial
We’ve already seen how to fit a n-th degree polynomial to n+1 points (x i ,y i ) i=1,…,n Lagrange interpolation polynomial For noisy data we would not want a fit to pass through all the points We want to find a best fit polynomial or order m
p
(
x
)
c
0
c
1
x
c
2
x
2 ...
c m x m
The sum of the squares of the residuals is
S
j n
1 (
p
(
x j
)
y j
) 2
Minimize the “error”
We next minimize S by taking partial derivatives w.r.t. the coefficients c n e.g. for c 0 : 0
S
c
0
c
0
j n
1 (
c
0
c
1
x j
c
2
x
2
j
...
c m x m j
y j
) 2 0 2
j n
1 (
c
0
c
1
x j
c
2
x
2
j
...
c m x m j
y j
)
nc
0
c
1
j n
1
x j
c
2
j n
1
x
2
j
...
c m j n
1
x m j
j n
1
y j
More examples
For c 1 0 we find:
S
c
1
c
1
j n
1 (
c
0
c
1
x j
c
2
x
2
j
...
c m x m j
y j
) 2 0
j n
1 2
x j
(
c
0
c
1
x j
c
2
x
2
j
...
c m x m j
y j
)
c
0
j n
1
x j
c
1
j n
1
x
2
j
c
2
j n
1
x
3
j
...
c m j n
1
x m
1
j
j n
1
x j y j
and so on… Taking partials w.r.t. to m+1 coefficients & equating to zero gives m+1 equations in m+1 unknowns (the c j j=0,…,m)
j j j n
1
n
n
1 1
n
x x x m j
2
j j j n
1
x j j n
1
x
2
j j n
1
j n
1
x
3
j x m
1
j
Matrix form
j n
1
x
2
j j n
1
x
3
j j n
1
j n
1
x
4
j x m
2
j
j n
1
j n
1
x m j x m
1
j j n
1
x m
2
j j n
1
x
2
m j
c c c
c n
0 1 2
j j
n
1
n
j n
1
j
1
n
1
x x x
m j
2
j y j y y j y j j j
These are the so-called “normal equations” in analogy with normals to the level set of a function Can be solved using the LU decomposition method we just discussed to find the coefficients c j
How do we choose
m
?
Higher order fits are not always better Many dependent variables are inherently following some underlying “law.” If the underlying relationship were quadratic, fitting a cubic would be largely pointless However, we don’t always know the underlying relationship So we need to look at the data…
Strategy
Plot the data If linear looks OK then start with m=1 If curve is visible try m=2 – if higher order then try log-log graph to estimate m Linear fit looks reasonable Start with m=1 Strong suggestion of underlying curve Start with m=2
Strategy cont.
After each fit evaluate S m r i =p(x i )-y i and plot the residuals If S m ~S m-1 then this indicates that the mth (or m-1 th) order is about as good as you can do Plotting the residuals will also give intuition to the quality of fit Residuals should scatter around zero with no implicit “bias”
r
Residuals
Example of a pretty good fit r This plot suggests you need an m+1 order fit Residuals scattered about zero with A large number of sign changes Residuals clearly have a positive bias with fewer sign changes
Extending the approach to fitting functions
We can use the method developed to fit a continuous function f(x) in an interval [a,b] Rather than a sum of squares we now have an integral that we need to minimize
S
j n
1 (
p
(
x j
)
y j
) 2
a
b
(
f
(
x
)
p
(
x
)) 2
dx
As an example, let f(x)=sin( p x) and the limits be [0,1], find the best quadratic (c 0 +c 1 x+c 2 x 2 ) to fit We need to minimize
S
a
b
(sin( p
x
) (
c
0
c
1
x
c
2
x
2 )) 2
dx
Compute derivatives
We work in the same way as the earlier fit:
S
c
0 0
c
0 0 1
dx
c
1 0 1
xdx
c
2 0 1
x
2
dx
0 1 sin p
xdx
S
c
1 0
c
0 0 1
xdx
c
1 0 1
x
2
dx
c
2 0 1
x
3
dx
0 1
x
sin p
xdx
S
c
2 0
c
0 0 1
x
2
dx
c
1 0 1
x
3
dx
c
2 0 1
x
4
dx
0 1
x
2 sin p
xdx
All integrals are tractable, and in matrix form we get: 1 1 / / 1 2 3 1 1 1 / / / 2 3 4 1 1 1 / / / 3 4 5
c c c
0 1 2 ( p 2 2 1 / / p p 4 ) / p 3
Comparing to Taylor expansion
If we solve the system on the previous page we derive a polynomial
p
(
x
) 0 .
0505 4 .
1225
x
4 .
1225
x
2 where 0 x 1 We can compare this to the Taylor expansion to second order about x=1/2 2 1 1 1
t(x)
f(
1
/
2
)
2
f'(
1
/
2
)
2 2
f
'' ( 1 / 2 ) 1 0 p 2 1 p 2 2 1 2 2 1 p 2 8 p 2 2 0 .
2337 4 .
9348
x
4 .
9348
x
2
x
p 2 2
x
2 Notice that the least squares fit and Taylor expansion are different
Taylor expansion is more accurate around the expansion point but becomes progressively worse further away Least squares fit is more accurate globally
Hilbert Matrices
The specific form of the coefficient matrix is called a “Hilbert Matrix” and arises because of the polynomial fit
Key point
: f(x) only affects the right hand side 1 1 / / 1 2 3 1 / 2 1 / 3 1 / 4 1 / 1 1 / / 3 4 5
c c c
0 1 2 ( p 2 2 1 / / p p 4 ) / p 3 For a polynomial fit of degree m the Hilbert Matrix has a dimension of m+1 As m increases, the determinant → 0 for m→∞ This creates significant numerical difficulties
Hilbert Matrices: problems & solutions
The Hilbert Matrix is classic example of an ill
conditioned matrix
By m=4 (i.e. a 5×5 Hilbert matrix) the small determinant causes difficulty in single precision By m=8 (i.e. a 9×9 Hilbert matrix) problems at double precision (real*8) We can improve upon the fitting method by choosing to fit a sum of polynomials: Instead of
p
(
x
) we use
p
(
x
) 0
c
0
p
0
c
1
x
(
x
) 1
c
2
x
2
p
1 (
x
) ...
2
c m x m p
2 (
x
) ...
m p m
(
x
) p i are ith degree polynomials
Fitting over basis polynomials
We’ll postpone the nature of the p i (x) for a moment… We now need to find the i by minimizing the least square function S
S
b
a
f
(
x
)
j m
0
j p j
(
x
) 2
dx
Find
S
i
0 normal equations :
a
b p i
(
x
)
j m
0
j p j
(
x
)
dx
a
b p i
(
x
)
f
(
x
)
dx
(1) , expanding the sum, 0
a
b p i
(
x
)
p
0 (
x
)
dx
1
a
b p i
(
x
)
p
1 (
x
)
dx
...
i a
b p i
2 (
x
)
dx
...
m a
b p i
(
x
)
p m
(
x
)
dx
a
b p i
(
x
)
f
(
x
)
dx
Choosing the p
i
(x)
In general, choosing polynomials appropriately will ensure that this is not a Hilbert Matrix We can again use an LU solve to find the i determine the fitting polynomial and hence However, if we had a system of polynomials for which
b a
p i
(
x
)
p j
(
x
)
dx
0 if i The system would be trivial to solve!
Such polynomials are said to be orthogonal j (2) Note that the integral is directly analogous to a dot product of vectors of rank N vectors (the integral has “infinite rank”)!
P i
.
P j
N k
1
P k i P k j
is analogous to
a
b p i
(
x
)
p j
(
x
)
dx
Simplifying the normal equations
If we have orthogonal polynomials, then from (1)
a
b p i
(
x
)
j m
0
j p j
(
x
)
dx
j m
0
j b a p i
(
x
)
p j
(
x
)
dx
a
b p i
(
x
)
f
(
x
)
dx
i a
b p i
2 (
x
)
dx
a
b p i
(
x
)
f
(
x
)
dx
i
a
b p i
(
x
)
f
(
x
)
dx
(3)
a
b p i
2 (
x
)
dx
Further simplification
In the case that the polynomials are orthonormal so that
b a
p i
(
x
) Then eqn (3) reduces to
p i
(
x
)
dx
1
i
b p i
(
x
)
f
(
x
)
dx a
There is no need to solve a matrix in this case!
Orthogonality definition: weighting functions
In practice it is often useful to include a weighting function, w(x), in the definition of the orthogonality requirement
b a
p i
(
x
)
p j
(
x
Orthonomality is then )
w
(
x
)
dx
0 if i j
b
p i
(
x
)
p j
(
x
)
w
(
x
)
dx
ij a
We can include the weighting function in the definition of the least squares formula for S
[a,b] [-1,1] [-1,1]
Examples of Orthogonal (Orthonormal) functions
w(x) 1 (1-x 2 ) -1/2 symbol P n (x) T n (x) Name Legendre polynomials Chebyshev I [-1,1] [0,∞) (-∞,∞) (1-x 2 ) -1/2 exp(-x) exp(-x 2 ) U n (x) L n (x) H n (x) Chebyshev II Laguerre Hermite See Arfken
Summary
LU decomposition is simple approach to solving matrices that uses both forward and backward substitution Least squares fitting requires setting up the normal equations to minimize the sum of the squares Fitting functions with a simple polynomial will results in ill conditioned Hilbert matrices More general fitting using a sum over orthogonal/orthonormal polynomials can reduce the fitting problem to a series of integrals
Next Lecture
Introduction to numerical integration