Method for estimating the Hessian of the objective function, times a vector, in order to compute an update in an iterative optimization solution to a partial differential equation such as the wave equation, used for example in full wave field inversion of seismic data. The Hessian times vector operation is approximated as one forward wave propagation (24) and one gradient computation (25) in a modified subsurface model (23). The modified subsurface model may be a linear combination of the current subsurface model (20) and the vector (21) to be multiplied by the Hessian matrix. The forward-modeled data from the modified model are treated as a field measurement in the data residual of the objective function for the gradient computation in the modified model. In model parameter estimation by iterative inversion of geophysical data, the vector in the first iteration may be the gradient of the objective function.