Julia: Cross Product of Complex Number incorrect

[jsgraber]

Version Info:
Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Xeon(R) CPU E5-2673 v4 @ 2.30GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)
Environment:
JULIABOX = true
JULIA_PKG_SERVER = https://pkg.juliacomputing.com
JULIA = /opt/julia-0.6/bin/julia
JULIA_KERNELS = ['julia-0.6', 'julia-1.0', 'julia-1.0k']
JULIA_PKG_TOKEN_PATH = /mnt/juliabox/.julia/token.toml
JULIABOX_ROLE =
JULIA_NUM_THREADS = 4

When taking the cross-product of 2 complex 3-vectors, Julia returns the complex-conjugate of the correct answer.

To understand the problem, first 2 facts:
1.) the dot product of any vector with itself should be the magnitude of that vector squared. Thus dot(x,x) = norm(x)^2
2.) the cross product should return a vector which is perpendicular to the input vectors.
3.) the dot product of perpendicular vectors should be zero. Thus, if z = cross(x,y) then dot(z,x) should be zero.

The following code shows that Julia properly takes the dot product of complex numbers, but does not take the cross product properly. It also shows that it returns the complex conjugate of the correct cross product. To use the complex function provided by Julia, you must first take the complex conjugate of both input vectors:

using LinearAlgebra
x = [1 + 2im, 2 + 1im, 3 - 1im]
y = [3 - 1im, 2 - 2im, 4 + 1im]
dot_xy = dot(x,y)
dot_xx = dot(x,x)
cross_xy = cross(x,y)
cross_xy_prime = cross(conj.(x),conj.(y))
println("dot(x,x) = ", dot_xx)
println("norm(x)^2 = ", norm(x)*norm(x))
println("dot(x,y) = ", dot_xy)
println("cross(x,y) = ", cross_xy)
println("cross(x',y') = ", cross_xy_prime)
println("dot(x,cross(x,y)) = " , dot(x,cross_xy))
println("dot(x,(cross(x,y)')) = ", dot(x,cross_xy'))
println("dot(x, cross_xy_prime) = ", dot(x, cross_xy_prime))

[daviehh]

Not a bug, the cross product is correct; as for dot(x,z), it's expected behavior: note that in julia, when using dot with complex vectors, the 1st one will be conjugated, see documentation here:

For complex vectors, the first vector is conjugated

you can check dot(conj.(x),cross_xy) is indeed 0.

[jsgraber]

For complex vectors, when taking the dot product, you need to take the complex conjugate of one of the two vectors. This is to ensure that dot(x,x) is real and is the magnitude of x squared. This is what you mentioned Julia is currently doing.

However, when taking the cross product of complex vectors, you need to either take the complex conjugate of each input vector and then run the "normal" cross product algorithm, or (equivalently) just take the complex conjugate of the result using the "normal" cross product algorithm and not taking the complex conjugate of the inputs. This is to ensure that the dot product and cross product remain consistent.

You can see a thorough discussion on the topic here: https://math.stackexchange.com/questions/129227/cross-product-in-complex-vector-spaces

[daviehh]

Changing can be potentially breaking for other's existing code; here maybe the convention taken is more programming than mathematical, i.e. the way julia calculates complex cross product is in line with e.g. matlab, mathematica and python/numpy?

[StefanKarpinski]

@stevengj, you seem like someone who might have a good opinion here. Any thoughts?

[StefanKarpinski]

Changing can be potentially breaking for other's existing code

If this is the correct definition, we could make such a change in a 1.x release as a "minor change"—it's somewhere between a bug fix and a breaking change. We'd have to put the change in NEWS, of course. If people feel that the current definition is simply wrong, however, this could be fixed and even backported to 1.0.5 and 1.2.1 as a bugfix.

[daviehh]

@StefanKarpinski Thanks, good to know!

Also, relevant reference here suggests cross product is only properly defined for R3 & R7

[stevengj]

I agree that changing to newcross(a,b) = conj(oldcross(a,b)) is arguably more consistent with the conjugated dot product, and probably qualifies as a minor change.

A complication is that there are lots of sources in engineering that don't conjugate the cross product of complex vectors. e.g. if you look at any electromagnetism textbook, things like the Poynting vector 𝐄×𝐇 are commonly computed for complex vectors (phasors) and there is no implicit conjugation — if you want to conjugate an argument you write it explicitly. On the other hand, in such contexts 𝐄⋅𝐄 is also not typically conjugated, and you write 𝐄*⋅𝐄=‖𝐄‖² if you want a proper inner product. (Mathematically, it turns out that unconjugated products are quite useful in electromagnetism because of reciprocity.) . One common usage would be computing a time-average like ½Re[𝐄*×𝐇] from time-harmonic (Fourier-transformed) fields, and fortunately this is not affected by an overall conjugation. Unconjugated "cross products" are also found in the common notation ∇× for the curl (the corresponding notation ∇⋅ for divergence is not affected by whether the ∇ is interpreted as "conjugated" since it is real).

So, it's likely that some people will be confused by defining cross as conjugated. On the other hand, the same people might also be confused by dot being conjugated.

[stevengj]

It looks like nearly all packages using cross() are doing so for 2d and 3d geometry calculations, which should use purely real vectors and hence would not be affected by a change here.

Packages like StaticArrays.jl, ContMechTensors.jl, Tensors.j, TensorFlow.jl, Grassmann.jl that provide new implementations of cross look like they might have to be updated, however.

[jsgraber]

Seeing the other use cases, I admit it is a tricky situation, which is not as clear cut as my initial thought on it was. I do find it odd, however, that the dot and the cross are handled differently. As a user, I would think wherever the development team ends up on the issue, that the documentation should be updated to explicitly mention the different definitions.

[antoine-levitt]

I think the current behavior is fine. People who want the "geometrical" definition (the conjugated one) will probably not even expect that julia implements it (since it's non standard). Otoh I can imagine using cross for complex numbers in the case where the geometrical interpretation is not relevant - eg x is a fixed real vector, and y is a complex Fourier component of a vector field, or something like that. In that case I would not expect cross to do something fancy, and that's be a subtle and annoying bug. So I think we should err on the side of caution against nasty bugs and principle of least surprise, at the price of not providing a function that is very non standard, and can easily be done by the user (just add a conj)

[StefanKarpinski]

Perhaps removing all support for complex cross products in Base would be prudent if there's no obviously correct definition?

[chethega]

The lowest effort approach would be to just document "The cross-product is linear in both arguments." and "The dot product is linear in its second argument and conjugate-linear in its first argument.". Together with the well-chosen example for cross, that informs about our entire choice of conventions.