Linear Inverse Problems A MATLAB Tutorial Presented by Johnny Samuels
Download
Report
Transcript Linear Inverse Problems A MATLAB Tutorial Presented by Johnny Samuels
Linear Inverse Problems
A MATLAB Tutorial
Presented by Johnny Samuels
What do we want to do?
• We want to develop a method to determine
the best fit to a set of data: e.g.
The Plan
•
•
•
•
Review pertinent linear algebra topics
Forward/inverse problem for linear systems
Discuss well-posedness
Formulate a least squares solution for an
overdetermined system
Linear Algebra Review
• Represent m linear equations with n variables:
a11 y1 a12 y2
a21 y1 a22 y2
a y a y
m1 1 m 2 2
a1n yn b1
a2 n yn b2
amn yn bm
a1n y1 b1
amn yn bm
a11
a
m1
A
y
b
• A = m x n matrix, y = n x 1 vector, b = m x 1 vector
• If A = m x n and B = n x p then AB = m x p
(number of columns of A = number of rows of B)
Linear Algebra Review:
Example
•
y1 y2 3 y3 y4 2
3 y1 3 y2 y3
y1 y2
1
2 y4 2
y1
1 2
1 1 3
y2
3 3 1 0 y 1
1 1 0 2 3 3
y
4
b
A
y
• Solve system using MATLAB’s backslash
operator “\”: A = [1 -1 3 1;3 -3 1 0;1 1 0 -2];
b=[2;-1;3]; y = A\b
Linear Algebra Review:
What does it mean?
•
y1 y2 1
2 y1 y2 1
1 1 y1 1
2 1 y2 1
A
• Graphical Representation:
y
b
Linear Algebra Review:
Square Matrices
• A = square matrix if A has n rows and n
columns
0
1 0 0
0
1
0
0
• The n x n identity matrix = I
0
1
A
• If there exists
s.t. A1 A I
invertible
1
then A is
Linear Algebra Review
Square Matrices cont.
• If
a b
A
c
d
then
A1
1 d b
ad bc c a
• Compute (by hand) and verify (with
MATLAB’s “*” command) the inverse of
2 1
A
1 3
Linear Algebra Review:
One last thing…
• The transpose of a matrix A with entries Aij
is defined as Aji and is denoted as AT – that
is, the columns of AT are the rows of A
• Ex:
1 2 3
A
0 0 4
implies
1 0
AT 2 0
3 4
• Use MATLAB’s “ ‘ “ to compute transpose
Forward Problem:
An Introduction
• We will work with the linear system Ay = b
where (for now) A = n x n matrix, y = n x 1
vector, b = n x 1 vector
• The forward problem consists of finding b
given a particular y
Forward Problem:
Example
• g = 2y : Forward problem consists of
finding g for a given y
• If y = 2 then g = 4
47 28
1
• What if A
and y 1 ?
89 53
• What is the forward problem for vibrating
beam?
Inverse Problem
• For the vibrating beam, we are given data
(done in lab) and we must determine m, c
and k.
• In the case of linear system Ay=b, we are
provided with A and b, but must determine y
Inverse Problem:
Example
• g = 2y : Inverse problem consists of finding y for
a given g
• If g = 10 then 212 y y 21 10 5
• Ay b
•
A1 Ay A1b
0 1 1 y1 3
2
4
1
y2 1
2 5 4 y 2
3
A
y
b
• Use A\b to determine y
y A1b
y1 0
1 1
y
2
4
1
2
y 2 5 4
3
1
3
1
2
Well-posedness
•
1
y
A
b produces
The solution technique
•
the correct answer when Ay=b is wellposed
Ay=b is well-posed when
1. Existence – For each b there exists an y such that
Ay=b
2. Uniqueness – Ay1 Ay2 y1 y2
1
3. Stability – A is continuous
•
Ay=b is ill-posed if it is not well-posed
Well-posedness:
Example
• In command window type y=well_posed_ex(4,0)
• y is the solution to
1
0
0
0
0 0 0 y1 1
1 0 0 y2 1
1
0 1 0 y3
0 0 1 y4 1
A
y
b
• K = condition number; the closer K is to 1 the
more trusted the solution is
Ill-posedness:
Example
• In command window type y=ill_posed_ex(4,0)
• y is the solution to
1
1
H y
1
1
where
H
1.0000 0.5000 0.3333 0.2500
0.5000 0.3333 0.2500 0.2000
0.3333 0.2500 0.2000 0.1667
0.2500 0.2000 0.1667 0.1429
• Examine error of y=ill_posed_ex(8,0)
• Error is present because H is ill-conditioned
What is an ill conditioned
system?
• A system is ill conditioned if some small
perturbation in the system causes a relatively large
change in the exact solution
• Ill-conditioned
system:
Ill-Conditioned System:
Example II
•
.835 .667 y1 .168
.333
.266
y
.067
2
A
•
y
y
y1 ?
y2 ?
b
.835 .667 y1 .168
.333 .266 y2 .066
A
y1 ?
y2 ?
b
What is the effect of noisy data?
•
•
Data from vibrating beam will be corrupted by
noise (e.g. measurement error)
Compare:
1.
2.
3.
4.
•
y=well_posed_ex(4,0) and z=well_posed(4,.1)
y=well_posed_ex(10,0) and z=well_posed(10,.2)
y=ill_posed_ex(4,0) and z=ill_posed(4,.1)
y=ill_posed_ex(10,0) and z=ill_posed(10,.2)
How to deal with error? Stay tuned for next
talk
Are we done?
• What if A is not a square matrix? In this case
does A1 not exist
• Focus on an overdetermined system (i.e. A is m x n
where m > n)
• Usually there exists no exact solution to Ay=b
when A is overdetermined
• In vibrating beam example, # of data points will
be much larger than # of parameters to solve (i.e.
m > n)
Overdetermined System:
Example
•
•
x
2
n
2
2
2
x
x
x
i
1
2
i 1
• Minimize
Ay b
2
2
xn2
Ay b
T
Ay b
Obtaining the Normal Equations
•
We want to minimize y Ay b Ay b :
T
y A
T
Ay b Ay b
T
T
A
AT Ay b AT Ay b
AT Ay AT b AT Ay AT b
2 AT Ay AT b
•
•
y
T
T
A
Ay
A
b
is minimized when y solves
y A A
T
solution
1
AT b
provides the least squares
Least Squares:
Example
• Approximate the spring constant k for
Hooke’s Law: l is measured lengths of
spring, E is equilibrium position, and F is
the resisting force
l E k F
b
y
k l E
T
1
l E l E T F
A
• least_squares_ex.m determines the least
squares solution to the above equation for a
given data set
What did we learn?
• Harmonic oscillator is a nonlinear system,
so the normal equations are not directly
applicable
• Many numerical methods approximate the
nonlinear system with a linear system, and
then apply the types of results we obtained
here