Joint Inversions in SimPEG #797

kimjaed · 2019-07-21T02:16:43Z

Hello,
My name is Jae Deok Kim and I am a Master's student working on Joint Inversions with Dr. Jiajia Sun at the University of Houston.

As far as I know, the joint inversions in SimPEG are inverting for different physical properties (e.g. electrical conductivity and volume) using a single dataset (https://docs.simpeg.xyz/content/examples/11-seis/seis_tomo_joint_with_volume.html).
We would like to add the functionalities for joint inversions that combine different physical properties (e.g. gravity and magnetics) from different datasets. This means that we would input multiple datasets of possibly different physical parameters (e.g. gravity and magnetics) and would output two different physical property models (e.g. density and susceptibility) that would be consistent in a predefined way (e.g. petrophysical relationships or structural similarity).

The joint inversion should be done within a single computational algorithm, with a single objective function and where all the models parameters are adjusted concurrently throughout the inversion.

Take for instance, the case of joint inversion of two models.
We define our objective function as:
$\large \Phi\left ( m_1,m_2 \right ) = \Phi_{d1} + \beta\Phi_{m1} + \Phi_{d2} + \beta\Phi_{m2} + \lambda\Phi_c$

The last term, $\large \Phi_c$ , is the coupling term that mathematically defines the relationship between different model parameters.
The most prominent example would be the cross-gradient.

To minimize this objective function, we can still use the Gauss-Newton method by making some linear algebraic modifications.

Let me first restate the objective function in explicit form:
$\large \Phi\left ( m_1,m_2 \right ) = \frac{1}{2}\left \|W_{d1}\left ( F[m_1] - d_1^{obs} \right ) \right \|_2^2 + \frac{1}{2}\beta_1\left \|W_{m1}\left ( m_1 - m_1^{ref} \right ) \right \|_2^2 + \frac{1}{2}\left \|W_{d2}\left ( F[m_2] - d_2^{obs} \right ) \right \|_2^2 + \frac{1}{2}\beta_2\left \|W_{m2}\left ( m_2 - m_2^{ref} \right ) \right \|_2^2 + \frac{1}{2}\lambda\phi_c$

We can combine this so that, we work with a single vector that represents the stacked models: $\small m=\begin{pmatrix} m_1 \\ m_2 \end{pmatrix}$
$\large \Phi\left ( m_1,m_2 \right ) = \frac{1}{2}\left \| \begin{bmatrix} W_{d1} & \\ & W_{d2} \end{bmatrix} \left [ \begin{pmatrix} F[m_1] \\ F[m_2] \end{pmatrix}\\ - \begin{pmatrix} d_1^{obs} \\ d_2^{obs} \end{pmatrix} \right ] \right \|_2^2 + \frac{1}{2}\left \| \begin{bmatrix} \sqrt{\beta_1} W_{m1} & \\ & \sqrt{\beta_2} W_{m2} \end{bmatrix} \left [ \begin{pmatrix} m_1 \\ m_2 \end{pmatrix}\\ - \begin{pmatrix} m_1^{ref} \\ m_2^{ref} \end{pmatrix} \right ] \right \|_2^2 + \frac{1}{2}\phi_c$

Our goal then becomes to find $\small m=\begin{pmatrix} m_1 \\ m_2 \end{pmatrix}$

At each Gauss-Newton step, we solve the system of linear equations given by:
$\large H\Delta m = -g$

where,
$\LARGE H = \begin{bmatrix} \frac{\partial^2\Phi}{\partial m_1^2} & \frac{\partial^2\Phi}{\partial m_1 m_2} \\ \frac{\partial^2\Phi}{\partial m_2 m_1} & \frac{\partial^2\Phi}{\partial m_2^2} \end{bmatrix}$

$\large \Delta m = \begin{bmatrix} \Delta m_1\\ \Delta m_2 \end{bmatrix}$

$\LARGE g = \begin{bmatrix} \frac{\partial \Phi}{\partial m_1}\\ \frac{\partial \Phi}{\partial m_2} \end{bmatrix}$

More explicitly, the Hessian can be stated as:
$\large H = \begin{bmatrix} J_1^T & \\ & J_2^T \end{bmatrix} \begin{bmatrix} W_{d1}^T & \\ & W_{d2}^T \end{bmatrix} \begin{bmatrix} W_{d1} & \\ & W_{d2} \end{bmatrix} \begin{bmatrix} J_1 & \\ & J_2 \end{bmatrix} + \begin{bmatrix} \sqrt{\beta_1}W_{m1}^T & \\ & \sqrt{\beta_2}W_{m2}^T \end{bmatrix} \begin{bmatrix} \sqrt{\beta_1}W_{m1} & \\ & \sqrt{\beta_2}W_{m2} \end{bmatrix} + \frac{1}{2}\lambda\begin{bmatrix} \frac{\partial \phi^2_c}{\partial m_1^2} & \frac{\partial \phi^2_c}{\partial m_1 m_2} \\ \frac{\partial \phi^2_c}{\partial m_2 m_1} & \frac{\partial \phi^2_c}{\partial m_2^2} \end{bmatrix}$

And the gradient can be explicitly stated as:
$\large g = \begin{bmatrix} J_1^T & \\ & J_2^T \end{bmatrix} \begin{bmatrix} W_{d1}^T & \\ & W_{d2}^T \end{bmatrix} \begin{bmatrix} W_{d1} & \\ & W_{d2} \end{bmatrix} \left [ \begin{pmatrix} F[m_1] \\ F[m_2] \end{pmatrix}\\ - \begin{pmatrix} d_1^{obs} \\ d_2^{obs} \end{pmatrix} \right ] + \begin{bmatrix} \sqrt{\beta_1}W_{m1}^T & \\ & \sqrt{\beta_2}W_{m2}^T \end{bmatrix} \begin{bmatrix} \sqrt{\beta_1}W_{m1} & \\ & \sqrt{\beta_2}W_{m2} \end{bmatrix}\left [ \begin{pmatrix} m_1 \\ m_2 \end{pmatrix}\\ - \begin{pmatrix} m_1^{ref} \\ m_2^{ref} \end{pmatrix} \right ] + \frac{1}{2}\lambda\begin{pmatrix} \frac{\partial \phi_c}{\partial m_1}\\ \frac{\partial \phi_c}{\partial m_2} \end{pmatrix}$

By constructing the Hessian and gradient as shown, we are able to take advantage of the optimization code that is already in SimPEG to implement our joint inversions.

These are the modifications that I've done:

Coupling: I've added a Coupling class that stores the methods related to the specific coupling strategy we use for joint inversions. For now, I've only coded the cross-gradient method. The relevant functions are deriv and deriv2 which compute the gradient and Hessian, respectively.
For more detailed information on the cross-gradient, you may refer to this paper by Gallardo and Meju (2003) https://doi.org/10.1029/2003GL017370
Directives: I've added a JointInversionDirective class that inherits from InversionDirective, but we define some of the methods differently. Namely, for joint inversions we are now working with multiple data misfits, multiple regularizations, and a coupling term, so I've added methods that call and set the values for these objects.
I've added a SaveOutputEveryIteration_JointInversion directive that will save the outputs of each iteration for joint inversions.
I've added a UpdatePreconditioner_JointInversion directive that sets a Jacobi preconditioner for the joint inversion.
Lastly, I've added a Adaptive_Beta_Reweighting directive that will adaptively change the trade-off parameters of the regularizations during the inversion. It makes sure that we stay within a certain range of our target data misfits.
InvProblem: I've created a new JointInvProblem class that inherits from BaseInvProblem. The relevant function is evalFunction, which computes the value of the objective function, as well as the gradient and Hessian (we store the Hessian as a scipy.sparse.linalg.LinearOperator object).
Optimization: I would like to add another stopping criteria, which I've named ratio_x, defined as $\large \frac{\left \| x_c-x_{last} \right \|}{\left \| x_{last} \right \|}$ .
This captures information on how much we've updated the model, and if the update is negligibly small we decide we've converged. It's somewhat similar to the moving_x stopping criteria, but different in the sense that the ratio_x is like a percentage on the model update at each iteration, where as moving_x is relevant to the initial model.
I've also added some IterationPrinters that are relevant for the joint inversion, such as printing multiple data misfits, multiple regularizations, the coupling term.
Lastly, I've modified the ProjectedGNCG findsearchdirection method, where I think it's more intuitive to use a for loop for i in range(self.maxIterCG) and have a if...then...break statement inside the loop, instead of a while loop. For the convergence criteria of the Conjugate Gradient, I've found that using r.dot(r) / self.g.dot(self.g) < self.tolCG works better especially in cases where the residual does not easily fall below magnitudes of 1e-3.
Inversion: I've added a JointInversion class that inherits from BaseInversion that simply sets the IterationPrinters and StoppingCriteria for the joint inversion.

It's my first time contributing to an open-source project, so any guidance on how to proceed is appreciated (e.g. what tests I need to do?).
I'm opening this issue to discuss the possibility of adding these functionalities for joint inversions in SimPEG.
I've tested these modifications locally and the joint inversion runs successfully, so I'm ready to start a pull request any time, but first I'd like to hear back from the team.

Best regards,
Jae

The text was updated successfully, but these errors were encountered:

fourndo · 2019-07-21T16:11:38Z

Well, first of all, this is an epic issue @kimjaed Prof. @jiajiasun must be stocked to have you as a Master student.

Looks like you have several interesting additions, and it will be easier for us all to bring it down in separate PRs (if possible). I'll take the first stab at this, but i am sure that @thast will have suggestions too when he gets back from vacation. You guys are working on the same problems.

1-

This is already possible, and called in our jargon a ComboObjectiveFunction. For the misfit functions, you can make a combo by first creating two different problems, with their own data and mesh.

phi_d1 = dataMisfit.l2_DataMisfit(survey1)
phi_d2 = dataMisfit.l2_DataMisfit(survey2)
comboMisifit = phi_d1 + phi_d2
where survey1 and survey2 are paired to their own problem.
Et voila, you can now invert both jointly. You still need to worry about mappings to relate the inversion parameters, but that's a seperate issue.
The same strategy applies to the regularization functions.

comboReg = reg1 + reg2

The combo class will take care of differentiating each component and adding them

fourndo · 2019-07-21T16:14:00Z

All the Gauss-Newton math looks right, but under the hood, we never form these large systems. We just do matrix-vector products and add vectors at the end.

fourndo · 2019-07-21T16:19:14Z

I am pretty excited about the Coupling part. It has been on my radar for a while to test.

Sounds like you would like to create a new sub-class of Regularization, and you will need a directive to go along with it to update the model. This can be brought in on a single PR with a small example. I can help you with writing a small test that will just check for consistency of results.

fourndo · 2019-07-21T16:20:21Z

For the Directives, we will have to see the functions. Again @thast will want to comment i m sure.

fourndo · 2019-07-21T16:21:39Z

Will have to take a look, but I think the JointInversion might be redundant once you get the ComboObjectiveFunctions``` working?

fourndo · 2019-07-21T16:28:48Z

I agree this is a great criterion to stop. The Optimization class needs some love and could hold it. At this point i would suggest you add this check to your JointInversion directive, and stop the inversion once you have reached it. In the directive you can use:
self.opt.stopNextIteration = True
This is what i did for the IRLS anyway. I will let the other guys comment.

lheagy · 2019-07-21T23:40:02Z

Welcome to the SimPEG community @kimjaed! It looks like you have been busy 🚀

Would you be able to join us at the SimPEG meeting this week or next and give us an overview of your work? We hold the meetings every Tuesday at 2pm PDT: http://bit.ly/SimPEG2019. If that doesn't work with your schedule, just let us know and we can reschedule for one week.

With all of the work you have done, a high-bandwidth conversation might be a good next step for getting your work included in SimPEG.

kimjaed · 2019-07-22T04:15:12Z

Welcome to the SimPEG community @kimjaed! It looks like you have been busy 🚀

Would you be able to join us at the SimPEG meeting this week or next and give us an overview of your work? We hold the meetings every Tuesday at 2pm PDT: http://bit.ly/SimPEG2019. If that doesn't work with your schedule, just let us know and we can reschedule for one week.

With all of the work you have done, a high-bandwidth conversation might be a good next step for getting your work included in SimPEG.

@lheagy I sure can! I'll be joining and preparing some material to make it a bit clearer on what exactly I'm trying to achieve.

lheagy · 2019-07-22T14:49:01Z

Excellent @kimjaed! I look forward to seeing you tomorrow. The meetings are quite conversational: demos / walk throughs of notebooks or code are encouraged 😄

kimjaed · 2019-07-22T14:54:45Z

@fourndo I've actually looked into the ComboObjectiveFunction. The problem that I see is that if,

m = np.r_[m1, m2]
phi_d1 = dataMisfit.l2_DataMisfit(survey1)
phi_d2 = dataMisfit.l2_DataMisfit(survey2)
comboMisfit = phi_d1 + phi_d2

and I compute comboMisfit.deriv(m), phi_d1.deriv(m) and phi_d2.deriv(m) will return Dimension mismatch, since the individual data misfits don't depend on the concatenated array m = np.r_[m1, m2].
If I instead compute comboMisfit.deriv(m1), I would be computing phi_d1.deriv(m1) + phi_d2.deriv(m1), but phi_d2 is the second data misfit, so really it should be phi_d2.deriv(m2).
Also, for purposes of the Gauss-Newton method, I don't want to simply add them, I want to stack them on top of each other.
The Gauss-Newton equations would need to look something like this:

What I have now is that if I predefine the number of models we are inverting for num_models=2,
I can loop over num_models, and concatenate the derivatives at every iteration.

num_models=2
m = np.r_[m1, m2]
phi_d1 = dataMisfit.l2_DataMisfit(survey1)
phi_d2 = dataMisfit.l2_DataMisfit(survey2)
comboMisfit = phi_d1 + phi_d2
for i in range(num_models):
        if i==0:
                phi_dDeriv = comboMisfit.objfcts[i].deriv(m[i])
        else:
                phi_dDeriv = np.r_[phi_dDeriv, comboMisfit.objfcts[i].deriv(m[i])]

And same with the regularizations and second derivatives. I would loop over the number of models, compute the derivatives individually, and concatenate them.

jcapriot · 2019-07-22T15:21:14Z

This is where SimPEG’s Mapping would come into play!

Each individual data misfit would take the full model, then the mapping would select which part of the model to use.

fourndo · 2019-07-22T15:39:22Z

@kimjaed

Yep, @jcapriot is right. The mapping you need is called a Wire map.

See the example here on how to use it:
This is for the regularization term of the MVI (3-components), but you will get the idea. It just knows where to insert/grab model parameters in your stacked vector,

lheagy · 2019-07-23T20:45:59Z

@kimjaed: for the meeting, here are the notes and the link to the zoom chat. I look forward to connecting!

notes: http://bit.ly/SimPEG2019
video chat: https://lbnl.zoom.us/j/9282619319

kimjaed · 2019-07-23T22:11:11Z

Something that I forgot to ask, or maybe it wasn't clear to me.

m = np.r_[m1, m2]
phi_d1 = dataMisfit.l2_DataMisfit(survey1)
phi_d2 = dataMisfit.l2_DataMisfit(survey2)
comboMisfit = phi_d1 + phi_d2

Computing comboMisfit(m) gives me phi_d1(m) + phi_d2(m), but I want to know these values separately. How can I access these values separately?

Edit: Do I have to access them by extracting the objfcts first?

phi_d1, phi_d2 = comboMisfit.objfcts
phi_d1(m)
phi_d2(m)

fourndo · 2019-07-23T22:19:29Z

You got it.
Or I usually just loop
for misfit in comboMisfit.objfcts:
misfit(m)

So it's more general and you can have as many as you want.

kimjaed · 2019-07-29T15:16:00Z

So I've been able to cut down my modifications to only a new coupling class and a set of new directives. I'm ready to start a PR and let you guys review my code.

lheagy · 2019-07-29T15:39:41Z

Nice work @kimjaed! I have started pr #799 to break up the regularization into multiple files. @fourndo is willing to do a review and then we can merge that in and you can work from there. If you are comfortable working off of a branch, you can branch off of regularization-organization from the main SimPEG repository and add your changes from there. Otherwise, we should be able to get this in within the next day or two and you can make your changes from there.

lheagy · 2019-07-30T04:48:51Z

@kimjaed: we merged #799 and released SimPEG 0.11.6, so at this point, feel free to branch directly off of master and start a pr. You should have received an invite to the SimPEG organization on github, so you can do all of this directly in the SimPEG org. Please let me know if you didn't receive the invite and I can re-send it. If you run into any issues, just let us know and we can jump on to a quick call if that is helpful. I look forward to seeing your pr! 🎉

jcapriot · 2022-10-28T19:37:44Z

Most of these changes were brought in in version 0.17.0 which included the start of the Coupling regularizations for joint inversion.

lheagy mentioned this issue Jul 24, 2019

break regularization into multiple files #799

Merged

jcapriot closed this as completed Oct 28, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Joint Inversions in SimPEG #797

Joint Inversions in SimPEG #797

kimjaed commented Jul 21, 2019

fourndo commented Jul 21, 2019 •

edited

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

lheagy commented Jul 21, 2019

kimjaed commented Jul 22, 2019 •

edited

lheagy commented Jul 22, 2019

kimjaed commented Jul 22, 2019 •

edited

jcapriot commented Jul 22, 2019 •

edited

fourndo commented Jul 22, 2019

lheagy commented Jul 23, 2019

kimjaed commented Jul 23, 2019 •

edited

fourndo commented Jul 23, 2019 •

edited

kimjaed commented Jul 29, 2019

lheagy commented Jul 29, 2019

lheagy commented Jul 30, 2019

jcapriot commented Oct 28, 2022

Joint Inversions in SimPEG #797

Joint Inversions in SimPEG #797

Comments

kimjaed commented Jul 21, 2019

fourndo commented Jul 21, 2019 • edited

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

fourndo commented Jul 21, 2019

lheagy commented Jul 21, 2019

kimjaed commented Jul 22, 2019 • edited

lheagy commented Jul 22, 2019

kimjaed commented Jul 22, 2019 • edited

jcapriot commented Jul 22, 2019 • edited

fourndo commented Jul 22, 2019

lheagy commented Jul 23, 2019

kimjaed commented Jul 23, 2019 • edited

fourndo commented Jul 23, 2019 • edited

kimjaed commented Jul 29, 2019

lheagy commented Jul 29, 2019

lheagy commented Jul 30, 2019

jcapriot commented Oct 28, 2022

fourndo commented Jul 21, 2019 •

edited

kimjaed commented Jul 22, 2019 •

edited

kimjaed commented Jul 22, 2019 •

edited

jcapriot commented Jul 22, 2019 •

edited

kimjaed commented Jul 23, 2019 •

edited

fourndo commented Jul 23, 2019 •

edited