  This paper considers the problem of approximating the inverse of the
wave-equation Hessian, also called normal operator, in seismology and other
types of wave-based imaging. An expansion scheme for the pseudodifferential
symbol of the inverse Hessian is set up. The coefficients in this expansion are
found via least-squares fitting from a certain number of applications of the
normal operator on adequate randomized trial functions built in curvelet space.
It is found that the number of parameters that can be fitted increases with the
amount of information present in the trial functions, with high probability.
Once an approximate inverse Hessian is available, application to an image of
the model can be done in very low complexity. Numerical experiments show that
randomized operator fitting offers a compelling preconditioner for the
linearized seismic inversion problem.
