diagonal approximation of Hessian matrix

Solving fitting goals (37) and (38) is expensive, since we have to propagate wavefields downward and upward within each iteration, with the cost for each iteration equal to the cost of two migrations. For small-scale problems, it is solvable; for large-scale problems, however, the computational cost might be prohibitive. What's more, currently there are no universal criteria for choosing the hyperparameters: $\epsilon$, which balances the data-fitting goal and the model-styling goal, and $\sigma$, which controls the sparseness of the model space. They can be decided only by trial and error, which obviously is not practical for very large-scale problems.

Instead of propagating wavefields at each iteration, however, we can precompute the Hessian or approximate it with a diagonal matrix and then solve the modified fitting goals iteratively. The solution of fitting goal (37) in the least-squares sense is

$${\bf m \approx \left( L' W_d' W_d L \right)^{-1}(W_dL)'W_d d }.$$ (40)

The weighted Hessian matrix ${\bf H = L'W_d'W_d L}$ can be either fully computed Valenciano and Biondi (2004) or approximated with a diagonal matrix; here I do the latter, approximating the weighted Hessian with its diagonals as follows Rickett (2003):  

$${\bf H \approx W_H = \frac{diag(L'W_d' W_d L m_{ref})}{diag(m_{ref})}},$$ (41)

and I choose the migrated image cube as the reference image cube:

$${\bf m_{ref} = (W_d L)' W_d d}.$$ (42)

Therefore, fitting goals (37) and (38) can be modified as follows:

$$\begin{eqnarray} 0 &\approx & {\bf W_H m - m_{mig}} \ 0 &\approx & \epsilon {\bf W_s} D({\bf m}),\end{eqnarray}$$ (43) (44)

where ${\bf m_{mig} = \left( W_dL \right)'W_d d}$, which is obtained by migrating the recorded data. To avoid the right-hand side of equation (41) being divided by zeros, I multiply ${\bf diag(m_{ref})}$ on both sides of equation (43), resulting in

$$\begin{eqnarray} 0 & \approx & {\bf W_{refm}m - W_{refd}m_{mig}} \ 0 & \approx & \epsilon {\bf W_s} D({\bf m}),\end{eqnarray}$$ (45) (46)

where ${\bf W_{refm} = diag(L'W_d' W_d L m_{ref})}$ and ${\bf W_{remd} = diag(m_{ref})}$. Fitting goals (45) and (46) can be solved by using the IRLS algorithm described in the previous section.