On Thu, Jun 16, 2011 at 3:16 PM, Alexander Grayver <span dir="ltr"><<a href="mailto:agrayver@gfz-potsdam.de">agrayver@gfz-potsdam.de</a>></span> wrote:<br><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div bgcolor="#ffffff" text="#000000">
Jed, thanks a lot for paper, it seems to be very good!<br>
<br>
You're right, there are ways to avoid calculating the whole A, but
there are actually advantages having A computed, which you don't
have without A being computed (at least with approaches I know so
far). I know several technique from paper you cited, maybe they
could be better (I personally believe they are not), but anyway now
I'm at the point when no serious changes are reasonable.<br>
That is right, F is nonlinear. The adjoint in my case costs the same
as forward problem. I can afford A for my problems. The only thing
which is really has to be avoided is the explicit assembling of A'*A
in Gauss-Newton. That I want to do using matrix-free formulation and
one of the iterative solvers from petsc. Would that be possible? If
not, then I really should think about different approach, maybe one
presented in paper you cited.<br></div></blockquote><div><br></div><div>Yes, you can do MF applications of A' A. You could just make a MatShell that called MatMult and MatMultTranspose.</div><div><br></div><div>For modifying the matrix, MatGetRow() IF you are using an AIJ matrix, return a pointer directly to the values. Thus you</div>
<div>could alter them in place. This is not portable, and might break later, but at least you could run now.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
<div bgcolor="#ffffff" text="#000000">
Regards,<br>
Alexander<br>
<br>
On 16.06.2011 16:38, Jed Brown wrote:
<blockquote type="cite">
<div class="gmail_quote">On Thu, Jun 16, 2011 at 15:43, Alexander
Grayver <span dir="ltr"><<a href="mailto:agrayver@gfz-potsdam.de" target="_blank">agrayver@gfz-potsdam.de</a>></span>
wrote:<br>
<blockquote class="gmail_quote" style="margin:0pt 0pt 0pt 0.8ex;border-left:1px solid rgb(204, 204, 204);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/%7Eomar/papers/siam-par-pde-opt.pdf" target="_blank">http://users.ices.utexas.edu/~omar/papers/siam-par-pde-opt.pdf</a></div>
<div><br>
</div>
</blockquote>
<br>
</div>
</blockquote></div><br><br clear="all"><br>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener<br>