<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
  <head>
    <meta content="text/html; charset=UTF-8" http-equiv="Content-Type">
  </head>
  <body 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>
    <br>
    Regards,<br>
    Alexander<br>
    <br>
    On 16.06.2011 16:38, Jed Brown wrote:
    <blockquote
      cite="mid:BANLkTinmi4XEVcvPY4TZJ859sKdPcfGmiw@mail.gmail.com"
      type="cite">
      <div class="gmail_quote">On Thu, Jun 16, 2011 at 15:43, Alexander
        Grayver <span dir="ltr">&lt;<a moz-do-not-send="true"
            href="mailto:agrayver@gfz-potsdam.de">agrayver@gfz-potsdam.de</a>&gt;</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 &gt;&gt; M (e.g. 10^7 &gt;&gt; 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 moz-do-not-send="true"
          href="http://users.ices.utexas.edu/%7Eomar/papers/siam-par-pde-opt.pdf">http://users.ices.utexas.edu/~omar/papers/siam-par-pde-opt.pdf</a></div>
      <div><br>
      </div>
    </blockquote>
    <br>
  </body>
</html>