# Computing preconditioner for a non-linear conjugate gradient implementation

Consider the following steps for the $i$-th non-linear conjugate gradient iteration, in the context of 3D electromagnetic inversion, and as discussed in (Newman and Boggs, 2004):

(1) set $i = 1$, choose initial model $\mathbf{m}_i$ and compute $\mathbf{r}_i = -\nabla\phi(\mathbf{m}_i)$

(2) set $\mathbf{u}_i = \mathbf{M}^{-1}\mathbf{r}_i$.

Here, $\phi$ is a misfit functional and "the matrix operator $\mathbf{M}$ in the algorithm is a preconditioner, which steers and scales the conjugate search direction $\mathbf{u}_i$ such that it more closely approximates the Newton direction" (Newman and Boggs, 2004).

I am interested in implementing the construction of $\mathbf{M}$ using MATLAB.

The context for the 3D EM problem is as follows: We are interested in minimizing $\phi(\mathbf{m}) = \frac{1}{2}||\mathbf{d}^{obs} - \mathbf{d}^{pred}||_2^2 + \frac{1}{2}\lambda||\mathbf{W}(\mathbf{m}_i)||_2^2$, where $\mathbf{m}$ is the model of Earth's electrical conductivity, $\mathbf{d}^{obs}$ and $\mathbf{d}^{pred}$ are the observed and predicted data vectors, respectively, $\lambda$ is the regularization parameter, $\mathbf{W}$ is the regularization matrix, based upon a finite-difference approximation of the Laplacian, and $\mathbf{D}$ is a data weighting matrix, which for the purpose of this question can be assumed to be the identity matrix. The predicted data arises from the forward problem solution, which yields the electric field $\mathbf{E}$ all over the control volume. This forward problem is no other that the Helmholtz equation for a secondary electric field $\mathbf{E}^S$: $\nabla\times\nabla\times\mathbf{E}^S + j\omega\mu_0\sigma^*\mathbf{E}^S = -j\omega\mu_0(\sigma^* - \sigma^{P*}\mathbf{E}^P)$, where $\mathbf{E} = \mathbf{E}^P + \mathbf{E}^S$ (scattered field formulation), $j$ is the imaginary unit, $\omega$ is the angular frequency, $\mu_0$ is the free space magnetic permeability constant, $\sigma^* = \sigma + j\omega\epsilon$ is the complex conductivity including real conductivity $\sigma$ and permittivity $\epsilon$, and $\sigma^{P*}$ is the conductivity of a layered background model (Grayver, Streich, and Ritter, 2013).

According to (Newman and Boggs, 2004): $\mathbf{M} = diag[(\mathbf{D}\mathbf{J})^H(\mathbf{D}\mathbf{J})] + \lambda\mathbf{W}^T\mathbf{W}$, where the operator $diag$ extracts the diagonal of the given matrix.

Ignoring the contributions from the term $\lambda\mathbf{W}^T\mathbf{W}$ (for now), it would seem like one can compute $\mathbf{J}^H\mathbf{J}$ by computing $\mathbf{E}^H\mathbf{E}$, since $\mathbf{J}^H\mathbf{J} = -(\omega \mu_0)^2\mathbf{E}^H\mathbf{E}$. Does this make sense? Th articles I have consulted seem to disregard detailed explanation about the construction of this preconditioner.

Newman, Gregory A.; Boggs, Paul T., Solution accelerators for large-scale three-dimensional electromagnetic inverse problems, Inverse Probl. 20, No. 6, S151-S170 (2004). ZBL1077.78012.

The core of the question is whether computing $\mathbf{M}$ via the assumption that $\mathbf{J}^H\mathbf{J} = -(\omega\mu_0)^2(\mathbf{E}^H\mathbf{E})$ makes sense or not.
If $\mathbf{d}^{pred}$ is the array of predicted data, then there must exists a prediction operator $\mathbf{f}(\mathbf{m}) = \mathbf{Q}\mathbf{E}$, $\mathbf{Q}$ is an interpolation operator placing data from the grid points at the receivers locations. If we assume (for now) , that all receivers are located on grid points, then $\mathbf{Q} = \mathbf{I}$, and $$\mathbf{J} = \frac{\partial \mathbf{f}(\mathbf{m})}{\partial \mathbf{m}} = \frac{\partial \mathbf{E}}{\partial \mathbf{m}}.$$ The solution to the discrete counterpart of the forward problem implies to solve a linear system $\mathbf{A}\mathbf{E} = \mathbf{S}$. Assuming a total field formulation of the forward problem, then: $$\mathbf{J} = \frac{\partial\mathbf{E}}{\partial\mathbf{m}} = \mathbf{A}^{-1}\left(-\frac{\partial\mathbf{A}}{\partial\mathbf{m}}\mathbf{E}\right)=\mathbf{A}^{-1}\left(-(i\omega\mu_0)\mathbf{E}\right).$$ Therefore: $$\mathbf{J}^H\mathbf{J} = \omega^2\mu_0^2(\mathbf{A}^{-1}\mathbf{E})^H(\mathbf{A}^{-1}\mathbf{E}),$$ where we can define $\mathbf{u}$ to be the solution to the system $\mathbf{A}\mathbf{u} = \mathbf{E}$.
This solution differs greatly from the originally-proposed identity: $\mathbf{J}^H\mathbf{J} = -(\omega\mu_0)^2(\mathbf{E}^H\mathbf{E})$. Specifically, in the later, there is no need to invert $\mathbf{A}$.