search

BAMresearch / bayem

Implementation and derivation of "Variational Bayesian inference for a nonlinear forward model." [Chappell et al. 2008] for arbitrary, user-defined model errors.
MIT License
2 stars 1 forks source link

Jacobian for global variables #22

Closed joergfunger closed 3 years ago

joergfunger commented 3 years ago

It seems that the Jacobian is not correctly evaluated when multiple individual parameters in a model error (e.g. E0 and E1) are mapped to a single global parameter (assuming E0 and E1 are to be jointly identified). In particular, in the inference problem

stacked_jac[:, self.latent[global_name].global_index_range()] = -J

stacked Jacobian has the dimension of [out, 1], whereas J has the dimension [out, 2].

TTitscher commented 3 years ago

That is a case, I did not consider. A use case would be to have E0 in elements 0..49 and E1 in elements 50..99. In that case, the Jacobians would be independent and could just be added up, right?

Honestly, I have no spontaneous, "beautiful" idea on how to deal with that... I'll think about it.

In the mean time: If you happen to have a test case, where you encountered this issue, could you push that to a branch?

joergfunger commented 3 years ago

Yes, that is exactly the usecase, and the Jacobians are just added.

joergfunger commented 3 years ago

I think I made a mistake the implementation is just working fine.

TTitscher commented 3 years ago

I was actually about to tackle this issue and come up with a test case that roughly represents the use case above. AFAI understand, the problem with the dimensions ([out,2] vs [out,1]) is a real issue. If you find the time, can you briefly explain that further?

TTitscher commented 3 years ago

@joergfunger As briefly discussed in person, I think the implementation is still incorrect. We currently do https://github.com/BAMresearch/BayesianInference/blob/master/bayes/inference_problem.py#L175

# loop over local_name
   ...
   J = parameter_jac[local_name]
   ...
   stacked_jac[:, self.latent[global_name].global_index_range()] = -J

Thus, if two local_names are mapped to the same global_name, one of the Jacobians would simply be overwritten.

A fix could be to simply change the " = -J" to "-= J",

   stacked_jac[:, self.latent[global_name].global_index_range()] -= J

but I would like to see a proper test for that.

(Note that the strange minus was my way to solve #19.)