<div class="gmail_quote">On Thu, Jun 16, 2011 at 15:43, Alexander Grayver <span dir="ltr"><<a href="mailto:agrayver@gfz-potsdam.de">agrayver@gfz-potsdam.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div bgcolor="#ffffff" text="#000000">Let's say we have a discretized model of some physical parameter m
(say acoustic velocity). Number of model parameters is N. We need to
take M measurements d within model (say time traveling of acoustic
wave) based on M different sets of receiver/source positions.<br>
In our case N >> M (e.g. 10^7 >> 10^3).<br>
We have operator F (which might linear or not) that defines
relationship between m and d:<br>
d=F(m)<br>
<br>
The operator F is just a system of equations actually. I have no
problem now to solve it for any m.<br>
<br>
What I need now is to compute the variation of this operator:<br>
<br>
A = \frac{\partial d_i}{\partial m_j}, for all i=1..M, j=1..N<br>
<br>
After some maths this calculation could be reformulated as a triple
product:<br>
<br>
A_i = C*F(m)^-1*v<br>
<br>
Where C some sparse matrix.<br>
<br>
Once the A is computed I want to solve another problem:<br>
<br>
(A'*A)b=A'*r<br>
<br>
Which is a Gauss-Newton system now, <br>
A'*A is truncated Hessian, <br>
A'*r is gradient, <br>
r = d - d_obs, where d_obs is the real observed data.<br>
b -- model change which have to applied to original model m in order
to explain your observed data better</div></blockquote></div><br><div>Excellent, great description.</div><div><br></div><div>It looks to me like you are assembling A because you want to be able to multiply by A' (which you can't do with finite differencing). But there is a very efficient way to apply A' that does not require assembling it: automatic differentiation. The adjoint costs about three times as much as a single forward solve. Depending on the nature of the nonlinearities in your problem, it may also be reasonably easy to write the adjoint without using AD.</div>
<div><br></div><div>You shouldn't need to assemble any of these matrices to solve the optimization/inverse problem. Assembling the matrices seriously limits scalability and there are optimal algorithms that don't use them. Depending on the problem, you may need a preconditioner for the reduced Hessian, but that can also be done without assembling it.</div>
<div><br></div><div>I imagine that F() is nonlinear. The adjoint of the nonlinear solve is the adjoint of the Jacobian at the converged state.</div><div><br></div><div>For more specifics on algorithms and pointers to the literature, you could look at this paper and some of the more specific applications it cites.</div>
<div><br></div><div><a href="http://users.ices.utexas.edu/~omar/papers/siam-par-pde-opt.pdf">http://users.ices.utexas.edu/~omar/papers/siam-par-pde-opt.pdf</a></div><div><br></div>