Using PyMC3 in a model with complex variables inside

Questions
1

Hi! I am trying to use estimate the parameters of a surface given reflectance data using a physical scattering model (mostly a complicated function). Everything works fine until the sampling, in which this error raises,

TypeError: Elemwise{real,no_inplace}.grad illegally  returned an integer-valued variable. (Input index 0, dtype complex128)

I only found one reference here in which is stated that this error raises when Theano founds a complex number, even if the resulting value is real (as is this case, in which I am computing the reflectivity and not its phase). Their solution (to bypass complex numbers) is not possible in my case. Is there any other work around?

Thank you for the feedback.

1 Like
2

I know you mentioned that you don’t care about the phase - would it be possible to convert the complex number into two real valued variables? Also, would you be able to provide a short(er) script that can generate the issue?

3

Thank you for the feedback!

I can implement all the complex calculus as vectorial, but I have to give up all the already builder numpy machinery for complex analysis. I prefeer not to do that.

I am trying to produce a nontrivial implementation to share with you, I will have it in the following days.

4

Hello,

I have been trying to work around the exact same problem for a long time.
Is there any progress in this direction ?

I could provide the function I use for my own reflectivity calculation (based on the Abeles transfer matrix formalism [https://en.wikipedia.org/wiki/Transfer-matrix_method_(optics)#Abeles_matrix_formalism]) but I think this is a general problem where a function internally involving complex types needs to be computed. The parameters and the output of the function are real numbers however.

5

Sadly, I wasn’t able to make this kind of problem work in PyMC3. To be fair, it is not a problem of PyMC3, since it seems that is Theano the one that cannot handle complex numbers, even if both input and output of the forward function are real numbers (like in my case also). I was unwilling to write everything in vectorial form, since complex number are the natural way to write scattering problems.

My solution was to implement the problem using TensorFlow Probability, which is based in TF which provides native support to complex numbers. It is much messier and unstable, but it works.

6

Thanks for the feedback. Yes, it is a shame.
Actually one can use the Bayesian approach but then one does not enjoy the NUTS sampler… and in my case the standard samplers are damn slow to converge.
After reaching the conclusion you mention I started looking at Julia (specifically at Turing.jl), but it is plagued by the same problem, the auto-diffenrentiation chokes on the complex numbers. :0(

I have zero knowledge of TensorFlow, would you care to share an example of your implementation?

7

I am sorry for the delay. To my understanding, TFP is much less polished than PyMC3. First, you have to migrate your model from Numpy to Tensorflow, which is straightforward (tensors instead of arrays, type tf.complex64 in our case). Yo have to replace also native basic functions of Numpy with the ones provided by TF (sin, cos, etc). To perform the actual inference, I followed the TF port of Probabilistic Programming and Bayesian Methods for Hackers to be found here. Please be aware that some of the code is not functional.

Good luck!

1 Like
8

Hey!

thanks a lot for your tips.
I’ll have a look to that, if the sampler described is just a Metropolis, than I am not so sure of the use… I would like to try NUTS which I hope would behave better for my use case (reconstruction of a scattering density profile from reflectivity data. Many parameters and a flattish chi2 landscape).
I also try to find some math students who would be willing to get a go at improving Turing.jl, which in principle should be doable.
All the best,
Jean-François

9

Hello Jean-François,

I would like to know if you found a solution to build a model with complex variables and also use NUTS from pymc3. I’m having a problem similar to what you once had.

Best regards,
Leo.

10

The issue of why complex gradients aren’t implemented is discussed a bit more in detail here: optimizing complex objective function using T.grad · Issue #3537 · Theano/Theano · GitHub

11

The suggestion is to model the real and imaginary components separately with vanilla variables.

@aseyboldt will be happy to hear if for some reason that is not possible in your case.