Trouble handling complex-valued variables

[ThomasA]

For example in functions._soft_threshold, a complex-valued input vector z will not be handled correctly although comments in the code suggest that it should. The reason seems to be that in this line sz becomes real due to the np.abs:

sz = np.maximum(np.abs(z) - T, 0)

Later when the final value of sz gets assigned, the fact that sz is already a np.float64 (real) causes the imaginary part to be discarded. This is confirmed by a ComplexWarning: Casting complex values to real discards the imaginary part when running the function with a complex-valued z:

sz[:] = np.nan_to_num(1. * sz / (sz + T) * z)

This can be easily fixed by over-writing sz with the new value instead of indexing into it (sz[:]):

sz = np.nan_to_num(1. * sz / (sz + T) * z)

But are there other good reasons not to do this? Another solution could be to change line 89 to:

sz = np.maximum(np.abs(z) - T, 0, dtype = z.dtype)

I would like to hear your thoughts on this, as I am hoping to contribute changes to (at least parts of) PyUNLocBox to correctly support complex-valued variables.

[nperraud]

I am not really aware if the pyunlocbox has been used with complex numbers or not. The Matlab version can handle complex number quite well. A least, I have been working with time-frequency frames a few times.

Your fix sz = np.nan_to_num(1. * sz / (sz + T) * z) seems the correct one according to me. But I am curious why @mdeff did it that way.

[mdeff]

Thanks for reporting and planning to contribute your changes @ThomasA 👍

The reason to use sz[:] = instead of sz = is to reuse the underlying memory instead of allocating a new chunk. Otherwise, the now unreferenced memory underpinning the old sz would uselessly take up RAM until its collection by the python garbage collector. Another way (maybe more explicit?) of achieving this is to use the out= parameter of numpy functions. If the data is large, that might make a big difference. Obviously it only works if the data type is the same. For me the best solution is for sz to have the same type as z, i.e. you second solution with dtype=z.dtype.

I've never worked with complex numbers either (maybe @rodrigo-pena did?). We're happy to take a PR which would made it work. Do not hesitate to come back to us if guidance is needed. Please also consider adding unit tests. That's probably the reason it doesn't work with complex numbers while claiming it. I think I modified it to reuse memory without realizing I would break the handling of complex numbers, and no tests told me.

[rodrigo-pena]

@mdeff , no, I've never worked with complex numbers within pyunlocbox, so I was unaware of the issue. For one, I second your approval for having someone work on it

[ThomasA]

I can start working on this, but I have no time-line for it. I have an upcoming research project application so I cannot say when and how much I will be able to work on it for the time being.
I think I will try a test-driven approach and start by implementing tests for complex-valued inputs to some of the lowest-level functions and then see where I get in trouble, working my way up from there.

[mdeff]

Thanks, that's an excellent strategy! Feel free to open a [WIP] PR to this repository, so that everybody can track, comment and contribute to the discussion.

No problem with the timeline. You're the one using it, so go as you need it. ;-)

[ThomasA]

Can I change the function signature of, for example _soft_threshold from:

def _soft_threshold(z, T, handle_complex=True)

to:

def _soft_threshold(z, T)

and simply determine the handling of complex numbers automatically based on the dtype of z instead?

[mdeff]

Sure. Looks better :)

[ThomasA]

@mdeff I was just wondering if I was allowed to change function interfaces. It will break code relying on this setting in earlier versions of PyUNLocBox.

[mdeff]

That's a private method (because of the underscore) so definitely. Even public interfaces we don't care too much as this is research software (meaning that user code should be pinned to a particular version for reproducibility concerns anyway). So feel free to change them if you feel like it's an improvement.

[mdeff]

Closing, as further discussion on that should go in #13. Thanks again @ThomasA for taking the time to do this properly. 👍