[FEATURE REQUEST] support for complex dtype

[v923z]

I would like to discuss the ramifications of adding support for complex arrays in ulab.

At the beginning, the Fourier transform was the single function that would have required complex arrays, but it was possible to just pass and return two real arrays that stood for the real and imaginary parts. This was not very elegant, but still acceptable. However, there have been at least two feature requests recently that would result in functions that lead out of the domain of real numbers.

Namely,

1. eigenvalues/eigenvectors of generic matrices [FEATURE REQUEST] linalg.eig for asymmetric matrices (and SVD?) #363
2. roots of polynomials [FEATURE REQUEST] Add numpy.polynomial.polynomial.Polynomial.roots() #355

The problem is that these functions could potentially return complex arrays, even if the input was real.

Adding a new dtype increases the firmware size, because the dtypes have to be resolved in C; in particular, the biggest contributor is probably the set of binary operators, because there one has to tackle all possible combinations, therefore, this contribution to the firmware size scales with the square of the number of dtypes. At the moment, we have 5 and a half dtypes ((u)int8, (u)int16, float/double, plus bool, which is actually an uint8 with a special flag), hence we deal with 25 combinations. By the addition of complex, we would end up with 36 combinations.

There are a number of functions that would not have to support complex. E.g., (arg)min, (arg)max, (arg)sort, median, clip, fmin etc. With these considerations, a rough estimate is 30% extra in the firmware size, if we add complex as a dtype.

Before jumping off the deep end, here are a couple of questions:

1. Is 30% extra in firmware size worth it? (This is about 30 kB in two dimensions.) Of course, we could add this as an option, selectable via a pre-processor switch, but the code still has to be implemented, and we would have to do something about the documentation: the most-sought feature is the Fourier transform, and if the complex can be excluded from the firmware, then there would be two tiers, where the behaviour of the fft function depends on the switch. That is highly undesirable.
2. Could the above-mentioned two feature requests be addressed in a different way, without the introduction of complex, as, e.g., in fft?
3. What should happen with the function pointers that provide an optimisation route? (These basically move all operations into the floating point domain, even if the input is an integer, thereby bypassing the problem of type resolution discussed above. This results in smaller firmware at the expense of execution speed. If complex is available, then that is the largest set of numbers, therefore, the function pointed to by the function pointer must return a complex number for all types. However, complexes can't just be passed to certain math functions in a haphazard manner. That is, why micropython implements its own cmath module.)
4. Would it make sense to support complex only as a return type, i.e., to allow a function to return a complex, but not as an argument, except for real, imag, and fft?
5. Along the same lines, is partial support (only some functions) acceptable? Could we say that a binary operator works with real arrays only, and we bail out, if an array is complex?
6. Are there hardware considerations? Some MCUs don't even have an FPU. Can we expect that complex math works on them? Having to sort out hardware limitations via pre-processor macros would definitely be a show stopper.

[jepler]

In CircuitPython, we enable complex values on any "full build", but never enable the cmath library. My gut tells me that users of CP don't need complex ulab arrays. When it comes to FFT, spectrum is more useful than the complex fft values themselves.

Since you bring up code size increase, it's worth remarking that the current CI of your master branch is not able to check whether additions to ulab potentially over-fill the firmware space of our various boards.

[KeithTheEE]

I'm not the most familiar with circuit pythons inner workings, but would this kind of proposal make more sense as a library/module proposal rather than adding another dtype into Circuit Python? I'm initially seeing it as sort of the same way Numpy has a different set of dtypes than Python does.

[v923z]

Since you bring up code size increase, it's worth remarking that the current CI of your master branch is not able to check whether additions to ulab potentially over-fill the firmware space of our various boards.

But this is something that I can't possibly check on this side: I have no idea, which (extra) modules you compile into circuitpython. Or have I misunderstood you?

[v923z]

I'm not the most familiar with circuit pythons inner workings, but would this kind of proposal make more sense as a library/module proposal rather than adding another dtype into Circuit Python? I'm initially seeing it as sort of the same way Numpy has a different set of dtypes than Python does.

That's right, but the dtypes that numpy defines are internal to numpy. There are no extra libraries that implement extra dtypes. Also, the question was raised in a very general context. I called out to a couple of circuitpython-related people, but that doesn't mean that the discussion concerns circuitpython only.

I tried to point out that the set of dtypes has serious consequences, and it is practically impossible to decouple from the core.

[KeithTheEE]

I see what you mean by saying it's hard to decouple from the core and especially ulab's functions--with the emphasis that completely real input can generate complex output. But the increase in size by adding a complex dtype is really large. Maybe it's worth it though for these use cases.

Looking at the libraries which do 'heavy' math in the python ecosystem--a majority of them use Numpy to recast the data in a complex friendly and efficient math fashion (np.array([1, 2, 3])), rather than just sticking with a standard list.

I'm wondering if this same approach would be beneficial here.
An additional 'bundle' library outside of the core could be created, such as ulab_complex, and an error type ComplexOutputError could be added to ulab for the cases where complex values were generated and need to be returned.
It would allow the math to become complex if a list were recast in an added ulab_complex.array([1,2,3]) function and avoid the sizeable change to the core libraries.

But that suggestion is no small suggestion either. It would have to be written, supported, your suggested issues involving FPUs are still there, and so on. But it could offer a solution and sidestep the memory issue. And again, I could be completely not seeing something core to the CircuitPython version when it comes to this issue.

[v923z]

Looking at the libraries which do 'heavy' math in the python ecosystem--a majority of them use Numpy to recast the data in a complex friendly and efficient math fashion (np.array([1, 2, 3])), rather than just sticking with a standard list.

But numpy already supports complex. In fact, numpy supports some 13 dtypes. I understand what you are saying, e.g., that pandas is built on top of numpy, but all those libraries are just wrappers around numpy, defining their own data container, but not the numerical types. numpy supports such definition of custom containers, and in fact, we are also trying to do something similar here: #327

However, numpy is still a closed numerical library in the sense that the results of all functions stay within the perimeter of those dtypes that are supported. Now, in this particular case, the problem is that the roots of an innocent-looking polynomial, such as x^2 + 1, can be complex.

Moreover, with an extra library, we cannot avoid the problem of binary operators. In numpy, ulab, or in python, you can always say

a = something
b = something else

c = a + b

At the moment, when + is evaluated, the interpreter (eventually, the C code) has to find out, what + actually means. If a, and b were two integers, then c will also be an integer, so in C, you will have something like this:

int _a = fetch_the_value_of(a);
int _b = fetch_the_value_of(b);

int _c = _a + _b;

return_the_value_of(_c);

where return_the_value_of(_c) will become your python object, holding the result of a + b. Now, in ulab, a, and b are ndarrays, and can assume 5 (6) different dtypes, so you have to write out all 25 (36) combinations in the C code. (You have to reserve the suitable amount of memory, have to know how you move your pointer etc., and these all depend on the types of a, b, and c.) I don't quite see at the moment, how we could decouple this step from the core. The implementation of the binary operator lies really at the heart of the definition of the ndarray. In fact, every data type in micropython (and python) can have a .binary_op method attached to it, and this .binary_op has to be implemented in C. This is where you have the 25 (36) combinations. I don't quite think that we can inject extra dtypes into this construct at a later stage.

If you care to see the declaration, here it is:

micropython-ulab/code/ulab.c

Lines 89 to 111 in 6d852a3

	const mp_obj_type_t ulab_ndarray_type = {
	{ &mp_type_type },
	#if defined(MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE) && defined(MP_TYPE_FLAG_EQ_HAS_NEQ_TEST)
	.flags = MP_TYPE_FLAG_EQ_CHECKS_OTHER_TYPE | MP_TYPE_FLAG_EQ_HAS_NEQ_TEST,
	#endif
	.name = MP_QSTR_ndarray,
	.print = ndarray_print,
	.make_new = ndarray_make_new,
	#if NDARRAY_IS_SLICEABLE
	.subscr = ndarray_subscr,
	#endif
	#if NDARRAY_IS_ITERABLE
	.getiter = ndarray_getiter,
	#endif
	#if NDARRAY_HAS_UNARY_OPS
	.unary_op = ndarray_unary_op,
	#endif
	#if NDARRAY_HAS_BINARY_OPS
	.binary_op = ndarray_binary_op,
	#endif
	.buffer_p = { .get_buffer = ndarray_get_buffer, },
	.locals_dict = (mp_obj_dict_t*)&ulab_ndarray_locals_dict,
	};

, and the binary operators are implemented here:

micropython-ulab/code/ndarray.c

Line 1716 in 6d852a3

mp_obj_t ndarray_binary_op(mp_binary_op_t _op, mp_obj_t lobj, mp_obj_t robj) {

I'm wondering if this same approach would be beneficial here.
An additional 'bundle' library outside of the core could be created, such as ulab_complex, and an error type ComplexOutputError could be added to ulab for the cases where complex values were generated and need to be returned.
It would allow the math to become complex if a list were recast in an added ulab_complex.array([1,2,3]) function and avoid the sizeable change to the core libraries.

What you are suggesting, in effect, is to split the complex arrays into a real and imaginary array by means of some sort of wrapper. This is OK as long as you just want to treat the return values, but I don't know how it would work, when you want to pass arguments.

Also, a recurring issue/request here is numpy compatibility. On the one hand, I totally understand the rationale, on the other hand, we can never achieve parity with numpy, and we have to find the balance between compatibility, firmware size, and execution speed.

But that suggestion is no small suggestion either. It would have to be written, supported, your suggested issues involving FPUs are still there, and so on.

That is a nasty can of worms, and I haven't the foggiest idea, how nasty it is.

But it could offer a solution and sidestep the memory issue.

I am not against such an approach, but we have to iron out the details.

[KeithTheEE]

numpy is still a closed numerical library in the sense that the results of all functions stay within the perimeter of those dtypes that are supported.

Ah this is extremely clarifying, thank you. I feel I'm beginning to understand the core of this issue. I wish I could more rapidly come to speed, but that's how it goes. Also thank you for linking me to the source examples, and taking the time to organize it all.

And the problem of real numbers generating complex results is an important issue. At the same time I agree full numpy parity isn't possible and should be the motivation.

I don't quite think that we can inject extra dtypes into this construct at a later stage

Can another module override this construct or provide it's own added dtype and safely return it? That way only if the other module is present would that construct have the 'added' complex cases for the binary operator? It would leave ulab being unable to handle complex math, but within ulab_complex for lack of a better name you could solve for complex roots of polynomials. Again, I'm showing my inexperience with the underling program design.

(I'm exclusively trying to see the scale of the problem, rather than assume I see the solution.)

[v923z]

numpy is still a closed numerical library in the sense that the results of all functions stay within the perimeter of those dtypes that are supported.

Ah this is extremely clarifying, thank you. I feel I'm beginning to understand the core of this issue. I wish I could more rapidly come to speed, but that's how it goes. Also thank you for linking me to the source examples, and taking the time to organize it all.

Not at all. Discussions are useful, when the points are clearly articulated. Besides, having a sounding board is really beneficial, especially, in complex cases like this (pun intended).

I don't quite think that we can inject extra dtypes into this construct at a later stage

Can another module override this construct or provide it's own added dtype and safely return it? That way only if the other module is present would that construct have the 'added' complex cases for the binary operator? It would leave ulab being unable to handle complex math, but within ulab_complex for lack of a better name you could solve for complex roots of polynomials.

You want to sort of sub-class the present definition of the ndarray, and add the option of extending it via an extra library. As far as I know, this is not possible at the C level, at least, not on the cheap. And unfortunately, even if we could pull that off, you would still have to extend all functions that can take an ndarray as their argument, and you would still have to equip them with facilities that treat complex arguments.

The thing is, numpy is not simply the definition of a class. Rather, it is a number of classes, and a number of functions that can take instances of these classes (and possibly other python types) as arguments. This makes the extendability that you suggest rather difficult.

And the problem of real numbers generating complex results is an important issue. At the same time I agree full numpy parity isn't possible and should be the motivation.

What could be done is perhaps the following:

1. optionally (via a pre-processor switch) add support for complex types. I could foresee this as an incremental change, first implementing the required code for the ndarray-level operators, and simply bail out in the functions, when the dtype is complex. Then I could add complex support in the functions, one by one.
2. leave the FFT alone, i.e., stay with the current call signature (two arrays in, two arrays out), even if the firmware supports the complex dtype. This is not numpy-conform, but at least the ulab documentation would be uniform, and the user wouldn't have to know, whether complexes are supported on a particular firmware.
3. allow complex functions (roots, complex eigenvalues etc.) only in those cases, when the complex dtype is present, otherwise, these functions would be disabled (the roots), or would work with limited functionality (as, e.g., the eigenvalue solver for symmetric real matrices).
4. let the user figure out any hardware-related issues with the complex math. I think I have to stand firm on this point, because I can't possibly, even in principle, safeguard against the absence of complex math functions on a particular MCU. I am not saying that, if an issue comes up on a concrete MCU, I am unwilling to add helper functions here and there, but such initiative should definitely come from the users.

This all is a fair amount of work, but such a solution would, on the one hand, not upset the status quo, and wouldn't take extra space in those cases, when only real numbers are needed, and on the other hand, it would still cut the Gordian knot of the complex dtype.

[v923z]

I have opened a draft under #366. Implementation-specific details should probably be discussed there.

[arnonsen]

I have no formal branch here but on my local copy I added the following dtypes: int32, uint32 and int64.
I did that without enlarging the number of operators loops (BINARY_LOOP, EQUALITY_LOOP etc.) but actually I'm in the process of getting rid of them, so my code size is in the process of becoming much smaller. I'm also not using the method of calling a casting function per element as this is unacceptably slow for my use.
Using my method complex can be added more easily and not only using complex but in 8/16/32 too without exploding the code size.