Linear Inverse Problems A MATLAB Tutorial Presented by Johnny - - PowerPoint PPT Presentation

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
• verdetermined system

Linear Algebra Review

• Represent m linear equations with n variables:
• 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)

11 1 12 2 1 1 21 1 22 2 2 2 1 1 2 2 n n n n m m mn n m

a y a y a y b a y a y a y b a y a y a y b + + + =     + + + =         + + + =  

• →

1 1 11 1 1 n m mn n m y b A

y b a a a a y b             =                 …

Linear Algebra Review:

Example

• Solve system using MATLAB’s backslash
• perator “\”: A = [1 -1 3 1;3 -3 1 0;1 1 0 -2];

b=[2;-1;3]; y = A\b

1 2 3 4 1 2 3 1 2 4

3 2 3 3 1 2 2 y y y y y y y y y y − + + = − + = − + − =

→

1 2 3 4

1 1 3 1 2 3 3 1 1 1 1 2 3

b A y

y y y y   −             − = −             −        

Linear Algebra Review:

What does it mean?

• Graphical Representation:

1 2 1 2

1 2 1 y y y y − + = − + = − →

1 2

1 1 1 2 1 1

b y A

y y −       =       − −      

Linear Algebra Review:

Square Matrices

• A = square matrix if A has n rows and n

columns

• The n x n identity matrix =
• If there exists s.t. then A is

invertible

1 1 1 I       =        

• 1

A−

1

A A I

−

=

Linear Algebra Review

Square Matrices cont.

• If then
• Compute (by hand) and verify (with

MATLAB’s “*” command) the inverse of

a b A c d   =    

1

1 d b A c a ad bc

−

−   =   − −  

2 1 1 3 A   =    

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: implies
• Use MATLAB’s “ ‘ “ to compute transpose

1 2 3 4 A   =    

1 2 3 4

T

A     =      

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
• What if and ?
• What is the forward problem for vibrating

beam?

47 28 89 53 A   =    

1 1 y   =   −  

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
• Use A\b to determine y

1 1

2 2 2 10 5 y y

− −

= = × =

1 2 3

1 1 3 2 4 1 1 2 5 4 2

b y A

y y y   −           − − =             − − −      

• ⇒

1 1 2 3

1 1 3 2 4 1 1 2 5 4 2 y y y

−

  −           = − −             − − −      

1 1 1

Ay b A Ay A b y A b

− − −

= ⇒ = ⇒ =

Well-posedness

• The solution technique produces

the correct answer when Ay=b is well- posed

• Ay=b is well-posed when
• 1. Existence – For each b there exists an y such that

Ay=b

• 2. Uniqueness –
• 3. Stability –

is continuous

• Ay=b is ill-posed if it is not well-posed

1

y A b

−

=

1 2 1 2

Ay Ay y y = ⇒ =

1

A−

Well-posedness:

Example

• In command window type y=well_posed_ex(4,0)
• y is the solution to
• K = condition number; the closer K is to 1 the

more trusted the solution is

• 1

2 3 4

1 1 1 1 1 1 1 1

b y A

y y y y                   =                        

Ill-posedness:

Example

• In command window type y=ill_posed_ex(4,0)
• y is the solution to where
• Examine error of y=ill_posed_ex(8,0)
• Error is present because H is ill-conditioned

1 1 1 1 H y       =        

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 H       =        

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

• 1

2

.835 .667 .168 .333 .266 .067

y A b

y y       =            

• ⇒

1 2

? ? y y     =        

• 1

2

.835 .667 .168 .333 .266 .066

y A b

y y       =            

• ⇒

1 2

? ? y y     =        

What is the effect of noisy data?

• Data from vibrating beam will be corrupted by

noise (e.g. measurement error)

• Compare:
• 1. y=well_posed_ex(4,0) and z=well_posed(4,.1)
• 2. y=well_posed_ex(10,0) and z=well_posed(10,.2)
• 3. y=ill_posed_ex(4,0) and z=ill_posed(4,.1)
• 4. 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 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)

1

A−

Overdetermined System:

Example

• Minimize
• (

) ( )

2 2 T

Ay b Ay b Ay b − = − −

2 2 2 2 1 2 2 1 n i n i

x x x x x

=

= = + + +

∑

Obtaining the Normal Equations

• We want to minimize

:

• is minimized when y solves
• provides the least squares

solution

( ) ( ) ( )

T

y Ay b Ay b φ = − −

( ) ( ) ( )

( )

( ) ( )

( )

2

T T T T T T T T T T T

y A Ay b Ay b A A Ay b A Ay b A Ay A b A Ay A b A Ay A b φ ∇ = − + − = − + − = − + − = −

( )

y φ

T T

A Ay A b =

( )

1 T T

y A A A b

−

=

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

• least_squares_ex.m determines the least

squares solution to the above equation for a given data set

( )

• b

y A

l E k F − =

• ⇒

( ) ( )

( ) (

)

1 T T

k l E l E l E F

−

= − − −

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