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: … + + + = � a a y b a y a y a y b 1 1 11 1 n 11 1 12 2 1 n n 1 → + + + = � = a y a y a y b � � � � � 21 1 22 2 2 n n 2 � � a a y b + + + = m 1 mn n m ������� � � � a y a y a y b m 1 1 m 2 2 mn n m y b A • 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 y − + + = 1 − y y 3 y y 2 • 1 1 3 1 2 1 2 3 4 → y − 2 = − − + = − 3 3 1 0 1 3 y 3 y y 1 1 2 3 y − 3 1 1 0 2 3 + − = ��� ����� � � � y y 2 y 2 y 1 2 4 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? − − + = = − → 1 1 y 1 • y y 1 1 = 1 2 − − − + 2 1 y 1 2 y y 1 � � � � � � � 2 1 2 b y A • Graphical Representation:
Linear Algebra Review: Square Matrices • A = square matrix if A has n rows and n columns � 1 0 0 0 � 0 1 0 0 • The n x n identity matrix = = I � � � � � 0 1 A − − 1 = 1 • If there exists s.t. then A is A A I invertible
Linear Algebra Review Square Matrices cont. − a b d b 1 − = = 1 • If then A A − − c d c a ad bc • 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 A ij is defined as A ji and is denoted as A T – that is, the columns of A T are the rows of A 1 0 1 2 3 = = • Ex: implies T A A 2 0 0 0 4 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 and ? A 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 − − = = × = 1 1 2 2 y y 2 10 5 • If g = 10 then − − − = ⇒ = ⇒ = 1 1 1 Ay b A Ay A b y A b • − 1 − − 0 1 1 y 3 y 0 1 1 3 1 ⇒ 1 • − − = = − − 2 4 1 y 1 y 2 4 1 1 2 2 − − − − − − 2 5 4 y 2 2 5 4 2 y � � 3 �� � ��� � 3 y b A • Use A\b to determine y
Well-posedness − = 1 • The solution technique produces y A b 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 – Ay Ay y y 1 2 1 2 A − 1 3. Stability – 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) 1 0 0 0 y 1 1 • y is the solution to 0 1 0 0 y 1 2 = 0 0 1 0 1 y 3 0 0 0 1 1 y � �� � ��� � � 4 b y A • 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) 1 1.0000 0.5000 0.3333 0.2500 1 0.5000 0.3333 0.2500 0.2000 • y is the solution to where = = H y H 0.3333 0.2500 0.2000 0.1667 1 0.2500 0.2000 0.1667 0.1429 1 • 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 y .835 .667 .168 � ⇒ y ? 1 = • 1 = .333 .266 y .067 ? y � �� ���� � 2 � � 2 y A b y .835 .667 .168 y ? � ⇒ 1 = 1 = • .333 .266 y .066 y ? 2 �� ���� � � � � 2 y A b
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 A − 1 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 )
Overdetermined System: Example • � � n ∑ = = + + + 2 2 2 2 • � x x x x x i 1 2 n 2 = i 1 ( ) ( ) 2 T − = − − • Minimize Ay b Ay b Ay b 2
Obtaining the Normal Equations ( ) ( ) ( ) T φ = − − • We want to minimize : y Ay b Ay b ( ) T ( ) ( ) ( ) T ∇ φ = − + − T y A Ay b Ay b A ( ) ( ) = − + − T T A Ay b A Ay b = − + − T T T T A Ay A b A Ay A b ( ) = − T T 2 A Ay A b ( ) = T T φ A Ay A b • is minimized when y solves y ( ) − 1 = T T • provides the least squares y A A A b solution
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 ( ) ( ( ) � − − = ⇒ 1 ( ) ( ) ) l E k F � T T = − − − k l E l E l E F � � � y b 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
Recommend
More recommend
