Complex Numbers

[felix91gr (Félix Fischer)]

@JeffreyEarly

So this has the advantage of reproducing both standard data structures, and not losing any information about the underlying type, by which I mean the array of Imaginary numbers has the right data structure for use in these algorithms, but also behaves correctly mathematically. I don’t think you can achieve this if you only have a Complex type.

I don't quite get what you mean here. What's the "right" data structure? What's the correct behavior you're pursuing? And why can't it be obtained with Complex types?

[JeffreyEarly (Jeffrey J Early)]

@felix91gr

How does (1) come to fruition? Don’t you need to introduce branching there to skip the floating point operations when the Imaginary or Real part are 0? If so, isn’t branching slower than FLOPs?

The idea is that if you introduce an Imaginary type you avoid branching. For example, consider multiplication of an imaginary number by an imaginary number.

1. If I do not have an Imaginary type, then you must compute *(Complex <T>,Complex<T>) -> Complex<T> which requires 6 FLOPs.
2. If I do have an Imaginary type, then several more overloaded multiplication operators are used, including *(Imaginary<T>,Imaginary<T>) -> <T> which only requires 1 FLOP.

Yes, you start overloading operators for all the cases, but when you start talking about vectors and matrices, the difference between 1 FLOP and 6 FLOPs starts to become enormous.

To achieve the same results without an Imaginary type, you have to start branching within your multiplication operation.

@felix91gr

(2) is true, if you’re using reference types and nullifying the Imaginary or Real part when they are 0 or haven’t been instantiated yet. But to use a reference comes at a price. In particular, you’d incur into:

A pointer’s cost. You need a pointer to the reference. Pointers are 64-bit, which effectively duplicates the size of the referenced Double type.
Cache misses: one of the weaknesses of reference types is that the referenced objects won’t necessarily be contiguous on memory. Therefore you’ll be triggering cache misses left and right when reading these objects. This can be avoided with a Complex struct that contains the Imaginary and Real pieces, but then you won’t be able to nullify 0-valued components in order to save memory.

I'm not thinking of Complex<T> as holding onto references, it's holding onto the values exactly like in @xwu 's implementation. The idea is that Complex<T> only gets used when necessary, e.g., +(<T>,Imaginary<T>) -> Complex<T>.

It's often the case that one wants to hold onto a large vector a purely imaginary numbers (for example, taking a derivative by multiplying the fourier transform of a function by wavenumber/frequency). Specifically,

let complexArray : Array<Complex> = [Complex(0,1),Complex(0,2)]
let imagArray : Array<Imaginary> = [Imaginary(1),Imaginary(2)]

the complexArray takes up twice as much memory as the imagArray.

[JeffreyEarly (Jeffrey J Early)]

@felix91gr

I don’t quite get what you mean here. What’s the “right” data structure? What’s the correct behavior you’re pursuing? And why can’t it be obtained with Complex types?

Meaning the structures DSPComplex and DSPSplitComplex. Both of these data structures get used in all the standard math libraries. The Complex type, stamped out in an array creates DSPComplex, an interleaved array. The only way (that I can think of) to achieve the latter in a natural way that preserves type information is to create an Imaginary types.

Specifically,

typealias DSPComplex = Array<Complex>

struct DSPSplitComplex {
    var real : Array<Real>
    var imag : Array<Imaginary>
}

So, you can't really get the type information of the DSPSplitComplex correct, if you don't have an Imaginary value. It's a very commonly used structure for matrices and vectors, so it'd be nice to have it well represented in Swift.

[JeffreyEarly (Jeffrey J Early)]

@felix91gr

How does (1) come to fruition?

I realize now that I should have just posted the structures to make this more clear. So, here's a sketch of how this works,

struct Real<T> where T:Numeric{
    var x : T
}
struct Imaginary<T> where T:Numeric {
    var x : T
}
struct Complex<T>  where T:Numeric {
    var x : Real<T>
    var y : Imaginary<T>
}

which have operations that look like this, e.g.,

func +<T:Numeric>(_ lhs: Imaginary<T>,_ rhs: Imaginary<T>) -> Imaginary<T>  {
    return lhs + rhs
}

func +<T:Numeric>(_ lhs: Real<T>,_ rhs: Imaginary<T>) -> Complex<T>  {
    return Complex(x:lhs,y:rhs)
}

func +<T:Numeric>(_ lhs: Imaginary<T>,_ rhs: Real<T>) -> Complex<T>  {
    return Complex(x:rhs,y:lhs)
}

[Nevin]

I don’t think those wrapper types are worthwhile. If we want to model split (in the signal-processing sense) complex numbers, the natural way to do so is like this:

struct SplitComplexArray<T: FloatingPoint> {
  var real: [T]
  var imaginary: [T]
}

[JeffreyEarly (Jeffrey J Early)]

That mirrors exactly what I wrote above for the DSPSplitComplex, but loses the type information about the imaginary part. Why is that better?

And, again---how do you solve memory and performance issues without an Imaginary type?

[felix91gr (Félix Fischer)]

You rose a few fair points that I'm not sure if I fully understand the implications of. I think I'll need to study this a bit more before I can give you my final answer

[scanon (Steve Canon)]

One important note: you don’t need to design an Imaginary type in at the beginning. It can pretty easily be bolted onto a working Complex model at any future point, so there’s some virtue in keeping things simple until you decide you really need it.

[CTMacUser (Daryle Walker)]

@Nevin

Furthermore, since a ComplexInt ought to conform to Numeric, it can’t even be generic over integer types, because of the magnitude requirement. I suppose one could resort to shenanigans like making the magnitude always be a Double regardless of the underlying integer type, but that seems…unprincipled.

What's wrong with ComplexInt.magnitude returning the modulus rounded? (It would be exact when the value is exactly on either the real or imaginary axis.) We just have to make sure that operations that would generally need the magnitude work with the Cayley norm, i.e. the square of the modulus, a.k.a. the value times its conjugate, instead. The Cayley norm would be an integer for integer types.

If we use rounded moduli, then ComplexInt would be only be a semi-norm if there are any non-zero values that have a modulus the rounds down to zero.

[CTMacUser (Daryle Walker)]

@xwu

@CTMacUser

There would be slightly different interfaces because of the division interface.

Not just slightly different; almost entirely different.

What interface differences are you thinking of? Besides analogies to the differences in the interfaces between the default integer and default floating-point types?

[Chris_Lattner3 (Chris Lattner)]

@taylorswift

This is problematic for matrices and vectors because for those things we really want the compiler to know about them and be able to optimize with simd/sse/whatever. This means they would have to be implemented in the standard library since they could not be compiled like ordinary Swift code.

I don't understand this point. SIMD support is defined in a normal Swift overlay in pure Swift code. The only thing the standard library can uniquely do is reach in and poke at LLVM operations directly. Other APIs should be able to compose on top of those, or on top of libraries like BLAS, given our approach to zero cost abstractions. Are you referring to something else?

-Chris

[JeffreyEarly (Jeffrey J Early)]

What does Numeric's .magnitude typically get used for? It certainly seems reasonable enough to me to interpret .magnitude as what you are calling the Cayley norm.

[JeffreyEarly (Jeffrey J Early)]

@Chris_Lattner3

given our approach to zero cost abstractions

Can you elaborate on what that means?

[Chris_Lattner3 (Chris Lattner)]

Int, Float and Bool are defined in the standard library. This is only practical if the language can make sure that doing so does not introduce overhead. "Zero cost abstractions" is a general concept though, available in C++, Rust and other languages: zero cost abstractions - Google Search

[Chris_Lattner3 (Chris Lattner)]

Incidentally, we're still working on open sourcing the code for Swift for TensorFlow. One of the things we need to define are host side "ND-Array" types (for both unknown and statically known rank). We have sketched a prototype, but expect it to be massively revamped when other folks get involved. It would be great to have a common forum to discuss the design of this sort of numeric API.

[taylorswift]

@Chris_Lattner3

I don’t understand this point. SIMD support is defined in a normal Swift overlay in pure Swift code. The only thing the standard library can uniquely do is reach in and poke at LLVM operations directly. Other APIs should be able to compose on top of those, or on top of libraries like BLAS, given our approach to zero cost abstractions. Are you referring to something else?

-Chris

I don’t understand. Currently you can define a structure containing 4 Floats and define custom operators like + to simulate vectorized operations but they’d still execute in sequence.

[felix91gr (Félix Fischer)]

@Chris_Lattner3

SIMD support is defined in a normal Swift overlay in pure Swift code.

Wait, seriously? I'm as surprised as Taylor. Where is this documented? If Swift already supports SIMD compilation, why didn't we know about it?

@Chris_Lattner3

“ND-Array” types (for both unknown and statically known rank)

Yes please

@Chris_Lattner3

It would be great to have a common forum to discuss the design of this sort of numeric API.

Let's make one. It could be one of these:

• A new category of topic on these forums
• A different forum (pls not a mailing list tho)
• Some sort of open slack?

Whatever it is, I agree that it's important to have one and to formalize this discussion. The community would benefit a great deal from a centralized approach to a numeric API.

There are currently ~5 Darwin vector libraries, 1 cross-platform (but very incomplete) one and well... the efforts are all being needlessly split, that's what I'm saying.

[felix91gr (Félix Fischer)]

@Chris_Lattner3 you should really take a look at our recent thread on this topic.

[Chris_Lattner3 (Chris Lattner)]

@taylorswift

I don’t understand. Currently you can define a structure containing 4 Floats and define custom operators like + to simulate vectorized operations but they’d still execute in sequence.

Right, but LLVM supports SIMD directly. this is what the existing simd library is implemented in terms of. If you can hold your nose and see through the gyb abstraction, the code is here:
https://github.com/apple/swift/blob/master/stdlib/public/SDK/simd/simd.swift.gyb

[taylorswift]

so to use this on Linux, what do I have to do?