Model fitting by least squares |
None of us is an expert in both geophysics and in optimization theory (OT), yet we need to handle both. We would like to have each group write its own code with a relatively easy interface. The problem is that the OT codes must invoke the physical operators yet the OT codes should not need to deal with all the data and parameters needed by the physical operators.
In other words, if a practitioner decides to swap one solver for another, the only thing needed is the name of the new solver.
The operator entrance is for the geophysicist, who formulates the estimation application. The solver entrance is for the specialist in numerical algebra, who designs a new optimization method. The C programming language allows us to achieve this design goal by means of generic function interfaces.
A generic solver subroutine is tinysolver. It is simplified substantially from the library version, which has a much longer list of optional arguments. (The forget parameter is not needed by the solvers we discuss first.)
void sf_tinysolver (sf_operator Fop /* linear operator */, sf_solverstep stepper /* stepping function */, int nm /* size of model */, int nd /* size of data */, float* m /* estimated model */, const float* m0 /* starting model */, const float* d /* data */, int niter /* iterations */) /*< Generic linear solver. Solves oper{x} = dat >*/ { int i, iter; float *g, *rr, *gg; g = sf_floatalloc (nm); rr = sf_floatalloc (nd); gg = sf_floatalloc (nd); for (i=0; i < nd; i++) rr[i] = - d[i]; if (NULL==m0) { for (i=0; i < nm; i++) m[i] = 0.0; } else { for (i=0; i < nm; i++) m[i] = m0[i]; Fop (false, true, nm, nd, m, rr); } for (iter=0; iter < niter; iter++) { Fop (true, false, nm, nd, g, rr); Fop (false, false, nm, nd, g, gg); stepper (false, nm, nd, m, g, rr, gg); } free (g); free (rr); free (gg); } |
The two most important arguments in tinysolver() are the operator function Fop, which is defined by the interface from Chapter , and the solver function stepper, which implements one step of an iterative estimation. For example, a practitioner who choses to use our new cgstep() for iterative solving the operator matmult would write the call
tinysolver ( matmult_lop, cgstep, ...
The other required parameters to tinysolver() are dat (the data we want to fit), x (the model we want to estimate), and niter (the maximum number of iterations). There are also a couple of optional arguments. For example, x0 is the starting guess for the model. If this parameter is omitted, the model is initialized to zero. To output the final residual vector, we include a parameter called res, which is optional as well. We will watch how the list of optional parameters to the generic solver routine grows as we attack more and more complex applications in later chapters.
Model fitting by least squares |