Nonlinear Problems: Grid Search and Monte Carlo Methods Overview

Lecture 14 Nonlinear Problems Grid Search and Monte Carlo Methods

Syllabus Lecture 01 Lecture 02 Lecture 03 Lecture 04 Lecture 05 Lecture 06 Lecture 07 Lecture 08 Lecture 09 Lecture 10 Lecture 11 Lecture 12 Lecture 13 Lecture 14 Lecture 15 Lecture 16 Lecture 17 Lecture 18 Lecture 19 Lecture 20 Lecture 21 Lecture 22 Lecture 23 Lecture 24 Describing Inverse Problems Probability and Measurement Error, Part 1 Probability and Measurement Error, Part 2 The L2 Norm and Simple Least Squares A Priori Information and Weighted Least Squared Resolution and Generalized Inverses Backus-Gilbert Inverse and the Trade Off of Resolution and Variance The Principle of Maximum Likelihood Inexact Theories Nonuniqueness and Localized Averages Vector Spaces and Singular Value Decomposition Equality and Inequality Constraints L1 , L Norm Problems and Linear Programming Nonlinear Problems: Grid and Monte Carlo Searches Nonlinear Problems: Newton s Method Nonlinear Problems: Simulated Annealing and Bootstrap Confidence Intervals Factor Analysis Varimax Factors, Empircal Orthogonal Functions Backus-Gilbert Theory for Continuous Problems; Radon s Problem Linear Operators and Their Adjoints Fr chet Derivatives Exemplary Inverse Problems, incl. Filter Design Exemplary Inverse Problems, incl. Earthquake Location Exemplary Inverse Problems, incl. Vibrational Problems

Purpose of the Lecture Discuss two important issues related to probability Introduce linearizing transformations Introduce the Grid Search Method Introduce the Monte Carlo Method

Part 1 two issue related to probability not limited to nonlinear problems but they tend to arise there a lot

issue #1 distribution of the data matters

d(z) vs. z(d) d(z) z(d) d(z) 6 6 6 5 5 5 4 4 4 3 3 3 d d z 2 2 2 1 1 1 0 0 0 0 2 4 6 0 2 4 6 0 2 4 6 z d z not quite the same intercept -0.500000 slope 1.300000 intercept -0.615385 slope 1.346154

d(z) d are Gaussian distributed z are error free z(d) z are Gaussian distributed d are error free

d(z) d are Gaussian distributed z are error free not the same z(d) z are Gaussian distributed d are error free

lesson learned you must properly account for how the noise is distributed

issue #2 mean and maximum likelihood point can change under reparameterization

consider the non-linear transformation m =m2 with p(m) uniform on (0,1)

p(m)=(m)- p(m)=1 2 2 1.5 1.5 p(m ) p(m) p(m) p(m) 1 1 0.5 0.5 0 0 0 0.5 m m 1 0 0.5 m m 1

Calculation of Expectations

although m =m2 <m > <m>2

right way wrong way

Part 2 linearizing transformations

Non-Linear Inverse Problem d d = g g(m m) transformation d d m m Linear Inverse Problem d d = Gm Gm solve with least-squares

Non-Linear Inverse Problem d d = g g(m m) transformation d d m m Linear Inverse Problem d d = Gm Gm solve with least-squares

an example di = m1 exp ( m2 zi ) d i=log(di) m 1=log(m1) m 2=m2 log(di) = log(m1) + m2 zi di =m 1+ m 2 zi

(A) (B) 1 0 0.8 -0.2 log10(d) 0.6 -0.4 log10(d) d d 0.4 -0.6 0.2 -0.8 0 -1 0 0.5 z z 1 0 0.5 z z 1 true minimize E=||d dobs d dpre||2 minimize E =||d d obs d d pre||2

again measurement error is being treated inconsistently if d d is Gaussian-distributed then d d is not so why are we using least-squares?

we should really use a technique appropriate for the new error ... ... but then a linearizing transformation is not really much of a simplification

non-uniqueness

consider di = m12 + m1m2 zi

consider di = m12 + m1m2 zi linearizing transformation m 1= m12and m 2=m1m2 di= m 1+ m 2 zi

consider di = m12 + m1m2 zi linearizing transformation m 1= m12and m 2=m1m2 di= m 1+ m 2 zi but actually the problem is nonunique if m is a solution, so is m a fact that can easily be overlooked when focusing on the transformed problem

linear Gaussian problems have well- understood non-uniqueness The error E(m) is always a multi-dimensioanl quadratic But E(m) can be constant in some directions in model space (the null space). Then the problem is non-unique. If non-unique, there are an infinite number of solutions, each with a different combination of null vectors.

600 500 400 E(m) E 300 200 100 0 0 1 2 mest 3 4 5 m m

a nonlinear Gaussian problems can be non-unique in a variety of ways

(E) m2 (B) (A) 0 500 0.2 600 400 E(m) 0.4 E(m) x1 m1 400 300 E E 0.6 200 200 0.8 100 0 0 1 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 0 1 2 3 4 5 m m x2 m m (F) m2 (C) (D) 0 60 25 0.2 20 40 0.4 E(m) E(m) 15 m1 m1 E E 0.6 10 20 0.8 5 0 0 1 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 0 1 2 3 4 5 m m m2 m m

Part 3 the grid search method

sample inverse problem di(xi) = sin( 0m1xi) + m1m2 with 0=20 true solution m1= 1.21, m2 =1.54 N=40 noisy data 4 2 d 0 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1

strategy compute the error on a multi-dimensional grid in model space choose the grid point with the smallest error as the estimate of the solution

(A) 4 2 d 0 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1 (B) 0 220 0.2 200 0.4 180 0.6 160 0.8 140 120 m1 1 100 1.2 80 1.4 60 1.6 40 1.8 20 2 0 0.5 1 1.5 2 m2

to be effective The total number of model parameters are small, say M<7. The grid is M-dimensional, so the number of trial solution is proportional to LM, where L is the number of trial solutions along each dimension of the grid. The solution is known to lie within a specific range of values, which can be used to define the limits of the grid. The forward problem d=g(m) can be computed rapidly enough that the time needed to compute LM of them is not prohibitive. The error function E(m) is smooth over the scale of the grid spacing, m, so that the minimum is not missed through the grid spacing being too coarse.

MatLab % 2D grid of m s L = 101; Dm = 0.02; m1min=0; m2min=0; m1a = m1min+Dm*[0:L-1]'; m2a = m2min+Dm*[0:L-1]'; m1max = m1a(L); m2max = m2a(L);

% grid search, compute error, E E = zeros(L,L); for j = [1:L] for k = [1:L] dpre=sin(w0*m1a(j)*x)+m1a(j)*m2a(k); E(j,k) = (dobs-dpre)'*(dobs-dpre); end end

% find the minimum value of E [Erowmins, rowindices] = min(E); [Emin, colindex] = min(Erowmins); rowindex = rowindices(colindex); m1est = m1min+Dm*(rowindex-1); m2est = m2min+Dm*(colindex-1);

Definition of Error for non-Gaussian statistcis Gaussian p.d.f.: E= d-2||e||22 but since and L=log(p(d d))=c- E E = 2(c L) -2L since constant does not affect location of minimum p(d d) exp(- E) in non-Gaussian cases: define the error in terms of the likelihood L E = 2L

Part 4 the Monte Carlo method

strategy compute the error at randomly generated points in model space choose the point with the smallest error as the estimate of the solution

(A) 4 2 d 0 0 0.1 0.2 0.3 0.4 0.5 x 0.6 0.7 0.8 0.9 1 (B) 0 (C) 220 0.2 100 200 B) 50 0.4 E 180 0 0.6 160 0 500 1000 iteration 1500 2000 0.8 140 1.5 120 m1 1 1 m1 100 1.2 0.5 0 500 1000 iteration 1500 2000 80 1.4 60 2 1.6 1.5 40 m2 1.8 20 1 0 500 1000 iteration 1500 2000 2 0 0.5 1 1.5 2 m2

advantages over a grid search doesn t require a specific decision about grid model space interrogated uniformly so process can be stopped when acceptable error is encountered process is open ended, can be continued as long as desired

disadvantages might require more time to generate a point in model space results different every time; subject to bad luck

MatLab % initial guess and corresponding error mg=[1,1]'; dg = sin(w0*mg(1)*x) + mg(1)*mg(2); Eg = (dobs-dg)'*(dobs-dg);

ma = zeros(2,1); for k = [1:Niter] % randomly generate a solution ma(1) = random('unif',m1min,m1max); ma(2) = random('unif',m2min,m2max); % compute its error da = sin(w0*ma(1)*x) + ma(1)*ma(2); Ea = (dobs-da)'*(dobs-da); % adopt it if its better if( Ea < Eg ) mg=ma; Eg=Ea; end end