Release Real type constraint on Quaternion

[sethaxen]

There are several variations on quaternions that uses non-real numbers with the quaternion multiplication operations, for which most/all of the current defaults can be kept. The two I'm thinking of are biquaternions and dual quaternions.

Currently we implement dual quaternions as 2 quaternions, one being the "primal" and one being corresponding "infinitesimal", but equivalently, it could be implemented as 4 Dual numbers, tying together the individual primals with infinitesimals. If this latter choice was made, DualQuaternion{T} could just be an alias for Quaternion{Dual{T}}, and we may we able to eliminate any special functionality for dual quaternions.

The steps would be:

1. Remove the Real type constraint for the T in Quaternion{T}.
2. Make any necessary changes to make the Quaternion implementations more general
3. Replace DualQuaternion's struct with a const, keeping the same tests.

These are breaking changes. Note that biquaternions would probably not work so long as we embed the complex numbers in the quaternions #62.

Relates #16 and #8

[hyrodium]

Oh, I was thinking that DualNumbers.Dual is a subtype of Real as well as ForwardDiff.Dual.

julia> using ForwardDiff, DualNumbers

julia> ForwardDiff.Dual <: Real
true

julia> DualNumbers.Dual <: Real
false

Related issues:

• Define ForwardDiffNumber <: Real JuliaDiff/ForwardDiff.jl#66
• Remove need for valtype of Dual to be <: Real JuliaDiff/ForwardDiff.jl#216
• Replace this package with ForwardDiff's Dual? JuliaDiff/DualNumbers.jl#45

I'm not much familiar with biquaternions, but it seems reasonable to change Quaternion{T<:Real} to Quaternion{T<:Number}.

[sethaxen]

Oh, I was thinking that DualNumbers.Dual is a subtype of Real as well as ForwardDiff.Dual.

Yeah, both are really purposed for automatic differentiation, but DualNumbers.Dual is also trying to be more of a special number type. It doesn't work super well though anyways, since real(::Dual) returns a Dual, which is not a number. ForwardDiff.Dual is quite a bit more modern though, but it's less convenient for defining and using Dual objects directly (Dual is more an implementation detail).

The other thing we could do is make an AbstractQuaternion and define Quaternion{T<:Real} <: AbstractQuaternion and DualQuaternion{T<:Real} <: AbstractQuaternion, changing the definitions for ::Quaternion to instead be ::AbstractQuaternion. This would be a non-breaking way to do the same thing as originally proposed.

Then there would also be the option of defining a UnitQuaternion if we decide in the future for structurally unit quaternions (not that we should do this).

I'm not much familiar with biquaternions, but it seems reasonable to change Quaternion{T<:Real} to Quaternion{T<:Number}.

Yeah, I'm pretty sure biquaternions just won't work without quite a bit of work (that we shouldn't do), but it would be nice if our types supported someone trying to implement them.

[sethaxen]

None of this package's direct dependents use DualQuaternion. In fact, no GitHub-hosted packages do. Since DualNumbers is our only dependency, this is a strong argument for:

1. Releasing the type constraint on Quaternion
2. Reimplementing DualQuaternion as a Quaternion{Dual}.
3. Removing DualQuaternion, adding to the docs an example of how one can use Dual and Quaternion together to get the same functionality

[hyrodium]

ForwardDiff.Dual is quite a bit more modern though, but it's less convenient for defining and using Dual objects directly (Dual is more an implementation detail).

Could you provide an example that ForwardDiff.Dual is less convenient than DualNumbers.Dual? I feel ForwardDiff.Dual <: Real is more useful, and ForwardDiff.Dual may be moved to DualNumbers.jl(JuliaDiff/DualNumbers.jl#45), or ChainRulesCore.jl(JuliaDiff/ForwardDiff.jl#579 (comment)) in the future.

I'm pretty sure biquaternions just won't work without quite a bit of work (that we shouldn't do), but it would be nice if our types supported someone trying to implement them.

Biquaternion with Quaternion{Complex} may be a good example, but Base doesn't support bicomplex with Complex{Complex}. I think we should avoid real(q::Quaternion{Complex}) with return value q.s::Complex, so we need to reconsider the apis.

For the above reasons, I now tend to avoid releasing the type constraint...

[sethaxen]

So both DualNumbers and ForwardDiff are not trying to implement dual numbers as a Number type, even though that's what they do. They're trying to implement a number that facilitates computing derivatives. ForwardDiff is more useful for this by masquerading its Dual as a Real, when dual numbers are definitely not real. DualNumbers is more accurate here, but still real(::Dual) returns a Dual, because that's what you want when you're doing AD, not when you're working with a hypercomplex number.

Because ForwardDiff is doing AD, its Dual implementation has features like that the infinitesimal part of the number can be a vector or tuple, not just a real number (to allow chunking, i.e. computing derivatives wrt multiple inputs simultaneously), which we wouldn't want

For DualQuaternion, we are not interested in AD. We really want a dual number. But since it doesn't require any extra code on our part, we could test the ForwardDiff approach now. I'll take a stab at writing up an example that could go in the docs.

[hyrodium]

I see. Thanks for the detailed explanation!

[sethaxen]

For DualQuaternion, we are not interested in AD. We really want a dual number. But since it doesn't require any extra code on our part, we could test the ForwardDiff approach now. I'll take a stab at writing up an example that could go in the docs.

Here's some code that checks Quaternion with ForwardDiff against our DualQuaternion implementation:

using Quaternions, ForwardDiff

function dualquat_to_quatdual(dq)
	q0args = let (; s, v1, v2, v3) = dq.q0
		(s, v1, v2, v3)
	end
	qeargs = let (; s, v1, v2, v3) = dq.qe
		(s, v1, v2, v3)
	end
	return Quaternion(ForwardDiff.Dual.(q0args, qeargs)...)
end

function dual_to_nondual(q::Quaternion{<:ForwardDiff.Dual})
	qargs = let (; s, v1, v2, v3) = q
		(s, v1, v2, v3)
	end
	q0 = Quaternion(ForwardDiff.value.(qargs)...)
	qe = Quaternion(ForwardDiff.partials.(qargs, 1)...)
	return DualQuaternion(q0, qe)
end
function dual_to_nondual(d::ForwardDiff.Dual)
	return Quaternions.Dual(ForwardDiff.value(d), ForwardDiff.partials(d, 1))
end
dual_to_nondual(x) = x

unary_functions = [+, -, exp, log, sign, abs, abs2, one, zero, conj, float, inv, sign, sqrt, angle]
binary_functions = [+, -, *, /, ==, ^]

dq1, dq2 = rand(DualQuaternionF64, 2)

# check all implemented functions. Difference should be about 0

map(unary_functions) do f
	f => f(dq1) - dual_to_nondual(f(dualquat_to_quatdual(dq1)))
end
map(binary_functions) do f
	f => f(dq1, dq2) - dual_to_nondual(f(dualquat_to_quatdual(dq1), dualquat_to_quatdual(dq2)))
end

# check that if we use ForwardDiff, we can still AD through the function

function foo(args)
	dq1 = DualQuaternion(Quaternion(args[1:4]...), Quaternion(args[5:8]...))
	dq2 = DualQuaternion(Quaternion(args[9:12]...), Quaternion(args[13:16]...))
	dq3 = dq1 * dq2
	q0 = dq3.q0
	qe = dq3.qe
	return [q0.s, q0.v1, q0.v2, q0.v3, qe.s, qe.v1, qe.v2, qe.v3]
end

flatten(x::ForwardDiff.Dual) = (ForwardDiff.value(x), ForwardDiff.partials(x, 1))

function foo2(args)
	dq1 = Quaternion(ForwardDiff.Dual.(args[1:4], args[5:8])...)
	dq2 = Quaternion(ForwardDiff.Dual.(args[9:12], args[13:16])...)
	dq = dq1 * dq2
	r = (dq.s, dq.v1, dq.v2, dq.v3)
	return [ForwardDiff.value.(r)..., ForwardDiff.partials.(r, 1)...]
end

args = randn(16)
foo(args) ≈ foo2(args)
ForwardDiff.jacobian(foo, args) ≈ ForwardDiff.jacobian(foo2, args)

Both implementations agree except for log, sqrt, ^, and angle. The first 3 all depend on log. The log implementation seems to be broken (all or most of all tests are marked @test_broken). The two angle implementations agree on the real part of the returned dual, which I think is fine.

I also tested that using Quaternion this way, one can still compute derivatives of the resulting code. And the answer is yes! So I don't see any good reasons to not drop DualQuaternion after documenting this. And no change to Quaternion's type constraints is necessary.

[sethaxen]

Closing since we instead opted for #92