I propose a computationally efficient method to approximate the inverse of the normal operator arising in the multi-parameter linearized inverse problem for reflection seismology in two and three spatial dimensions. Solving the inverse problem using direct matrix methods like Gaussian elimination is computationally infeasible. In fact, the application of the normal operator requires solving large scale PDE problems. However, under certain conditions, the normal operator is a matrix of pseudodifferential operators. This manuscript shows how to generalize Cramer's rule for matrices to approximate the inverse of a matrix of pseudodifferential operators. Approximating the solution to the normal equations proceeds in two steps: (1) First, a series of applications of the normal operator to specific permutations of the right hand side. This step yields a phase-space scaling of the solution. Phase space scalings are scalings in both physical space and Fourier space. Second, a correction for the phase space scaling. This step requires applying the normal operator once more. The cost of approximating the inverse is a few applications of the normal operator (one for one parameter, two for two parameters, six for three parameters). The approximate inverse is an adequately accurate solution to the linearized inverse problem when it is capable of fitting the data to a prescribed precision. Otherwise, the approximate inverse of the normal operator might be used to precondition Krylov subspace methods in order to refine the data fit. I validate the method on a linearized version of the Marmousi model for constant density acoustics for the one-parameter problem. For the two parameter problem, the inversion of a variable density acoustics layered model corroborates the success of the proposed method. Furthermore, this example details the various steps of the method. I also apply the method to a 1D section of the Marmousi model to test the behavior of the method on complex two-parameter layered models.