Inconsistency between torch.clamp() and numpy.clip() behavior for complex numbers

[anjali411]

a=torch.tensor([3+4j])
torch.clamp(a, 2+3j, 1+5j)
tensor([(3.0000 + 4.0000j)], dtype=torch.complex128)
torch.clamp(a, 1+5j, 2)
tensor([(1.0000 + 5.0000j)], dtype=torch.complex128)
np.clip(a.numpy(), 1+5j, 2)
array([2.+0.j])

cc @ezyang @anjali411 @dylanbespalko

[ezyang]

I think I mentioned this to @dylanbespalko when I reviewed his PR that added clamp support; I have no idea what this is actually supposed to do here :>

[dylanbespalko]

I have an idea. I don't agree with what the np.clip() operation is doing, but it is performing a clamp operation using a complex min and complex max. Here is my solution.

def clamp(complex a, complex min, complex max):
    # same as np.clip, choose the complex number with the smallest absolute value.
    return max(min(a, max), min) 

def clamp_abs(complex a, real min, real max):
    # Maintain direction, but control the speed of the plane during landing.
    mag = max(min(abs(a), abs(max)), abs(min))
    return mag * exp(1j * a.angle())

def clamp_angle(complex a, real min, real max):
    # Maintain speed, but clamp the angle so the spaceship doesn't burn-up on re-entry.
    phase = max(min(angle(a), abs(max)), abs(min))
    return abs(a) * exp(1j * phase)

This way clamp uses the same underlying max(min(*)) as real numbers for autograd purposes, and complex support will get two new functions that reliably perform a common and important task. For real numbers clamp_abs or clamp_phase should be unsupported or no-ops.

[ezyang]

@rgommers I wonder if you have an opinion here

[rgommers]

The behavior of numpy.clip with complex input isn't tested (except for one test checking it doesn't segfault). I agree with @dylanbespalko that the behavior isn't what's desired:

>>> np.clip([3 + 4.j], -1, 2)                                                                                       
array([2.+0.j])
>>> np.clip([3 + 4.j], -1+1.j, 2+12.j)  # imaginary component goes up                                                                   
array([2.+12.j])
>>> np.clip([1 + 4.j], -1+1.j, 2+12.j)  # imaginary component doesn't go up                                                                             
array([1.+4.j])

So that's not something to emulate exactly. I'll open an issue on the NumPy issue tracker.

Reasonable behavior could be one of:

• Clip real and imaginary parts separately
• Clip absolute value, not changing the phase
  There may be other options. As long as it's documented, I'm not sure it matters much which choice is made.

It may be worth waiting on the outcome of the NumPy discussion before doing anything here, given that this doesn't look urgent to fix.

[jscottcronin]

@rgommers @anjali411
I'd like to begin contributing to pytorch and saw this issue is still open. I looked through numpy related issues:

• np.clip with complex input is untested and has odd behavior numpy/numpy#15630
• BUG: complex special case for clip numpy/numpy#15643

It seems like numpy folks were looking to implement:

• Clip real and imaginary parts separately
• (as seen here: https://github.com/numpy/numpy/pull/15643/files)

Before I get started, do we agree that we should take this route, or should we continue to wait on numpy discussion? If we want to wait on numpy decision, can someone suggest another issue for me to pick up?

[rgommers]

I think the NumPy discussion is done to a large enough extent I think, that issue went through triage review and there's consensus that clipping real and imaginary parts separately is the most reasonable thing to do. So please go ahead I'd say.

[jscottcronin]

Hey ya'll,

I'd like to get some feedback and suggestions on my understanding thus far. I found file aten/src/ATen/native/UnaryOps.cpp with function clamp. Is this where I should be editing? I also saw function clamp_stub being called if there is a min and max. However, I didn't understand where that function is defined, or if I need to make changes to that def.

There are a few issues going on here:

Issue 1

a=torch.tensor([3+4j])
torch.clamp(a, 1+5j, 2)

# Actual
tensor([(1.0000 + 5.0000j)], dtype=torch.complex128)

# Expected
tensor([(2.0000 + 0.0000j)], dtype=torch.complex128)

How to fix: I am not sure and could use some guidance.

Issue 2

a=torch.tensor([3+4j])
torch.clamp(a, 1+5j, 2+6j)  # Note 6j > 4j

# Presumed result given issue 1 fixed
tensor([(2.0000 + 6.0000j)], dtype=torch.complex128)

# Expected
tensor([(2.0000 + 4.0000j)], dtype=torch.complex128)

How to fix: separate real and imag parts of tensor and perform clamp on each using real and imag values of min and max. Function in numpy looks like

def cx_clip(self, a, m, M, out=None):
    # use slow-clip for complex case
    return (self.clip(a.real, m.real, M.real)
    + (1.0j * (self.clip(a.imag, m.imag, M.imag))))

In addition, I'm happy to say this will be my first time writing C++ code. Are there any guides I can follow that would help me understand how to test code with short feedback loops? Very basic example with explicit instructions would be great! Super excited to contribute and thank you for helping me!

[dylanbespalko]

@jscottcronin,

In #35563, I recently gave @kostekIV some advice on how to add a polar() math kernel. You can use this as a guide to find which files you need to modify. You should also search for the updated guide to adding a new kernel to Pytorch.

Beyond this, my opinion for complex clamp ops:

1. Modify the existing clamp() to perform numpy.clip() on the real and imag parts as if they are completely independent floats. This should be compatible with the numpy clip function.
2. Create a new op called clamp_abs or clamp_mag that clamps the magnitude, but not the phase of the complex number.
3. Create a new op called clamp_angle or clamp_phase that clamps the angle or phase, but not the magnitude of the complex number.

I have used abs and angle to be compatible with numpy.

@rgommers is a numpy developer, and he said they agreed to make numpy.clip() independently clamp the real and imag components.

You should submit a PR for step 1, and another PR for step 2 & 3 because we want to be backwards compatible with numpy.

[dylanbespalko]

@jscottcronin

The clamp kernel is implemented in aten/src/ATen/native/cpu/UnaryOpsKernel.cpp.

In this file you will see functions that call cpu_kernel() and functions that call cpu_kernel_vec(). The cpu_kernel() is easier to get working. If you look at the real() kernel, the cpu_kernel_vec() calls:

1. Vec256::real() for 256 bit = (32bit * 8 floats) hw accelerated groups of numbers.
2. std::real() aka real_impl() for the remaining numbers in the tensor.

Hence, when you test the code, you need to choose a tensor length (like 9) that will call both functions.

Here are some other links:

PyTorch Internals
Vec256(Intel AVX2)

[gchanan]

I think we should just leave this disabled for now, as in 36373.

The current NumPy behavior (clamp the real) doesn't make a lot of sense, but it's at least consistent with the definition of min/max, meaning it just sorts according to the definition of min/max and then applies the min/max clamps to the low and high ends of ranges.

The new suggested NumPy behavior (clamp the components) also doesn't make a lot of sense, and it's inconsistent with the behavior of min/max (which doesn't look at the complex component).

As pointed out in numpy/numpy#15630, we can't point to a single person who wants any of this. Also, both the current and new behavior are trivially expressible in numpy by accessing real and imag.

The other logical choice (absmax / absmin / absclamp) may justify a different op or the ability to pass an evaluation function (complex->real) for performance reasons if that's what people want.

So, I'd recommend just disabling all of this for now until we get actual requests.

[n-west]

Just my $.02, but I currently use torch::clamp on "complex" tensors of shape [samples, 2] to simulate clipping on a pair of independent ADCs such as in quadrature sampling. In this application I would expect that a float min or max would be applied to the real and imag channels independently. I don't know of an application of a complex clamp value, but I don't think the real clamp value is operating on the channels independently and it's not doing whatever numpy is doing (also not operating independently). The clamp angle and abs make sense, but I don't know of an immediate application for them

[mruberry]

Let's move this discussion to #36444.

[mruberry]

Per discussion on this issue and #36444, it seems there's a consensus on disabling complex clamp (clip) in PyTorch. Disabling lets us choose an implementation later, if we like, and avoid user confusion. The multiple valid complex clamp implementations suggest that users will be unable to anticipate what clamp does. Helpful behaviors, like those @n-west mentions, can also be implemented with new, more explicit functions that clarify users' intent. Maybe torch.absclamp?