Implementing matrix functions with quaternion eltypes

[sethaxen]

Quaternions can be represented as $2 \times 2$ complex matrices whose matrix product preserves the Hamilton product. This means that one can map an $n \times n$ quaternion matrix to a $2n \times 2n$ complex matrix with $2 \times 2$ blocks from quaternions to perform any complicated operations that only consist of matrix addition and multiplication.

Matrix functions such as the matrix exponential satisfy this property. We could generically support these matrix functions for quaternion eltypes by explicitly generating these complex matrices, dispatching to the complex matrix functions, and then mapping back to a matrix of quaternions.

This would generalize #46 and extend #56 to matrices.

[sethaxen]

Here's a minimal example. We'd need to think more carefully about indexing for a real implementation, which might mean we restrict inputs to StridedMatrix types and triangular matrices:

julia> using Quaternions, LinearAlgebra

julia> function quaternion_to_complex_matrix!(z::AbstractMatrix, q::Quaternion)
           w = complex(q.s, q.v1)
           y = complex(q.v2, q.v3)
           z[begin, begin] = w
           z[begin + 1, begin] = -conj(y)
           z[begin, begin + 1] = y
           z[end, end] = conj(w)
           return z
       end
quaternion_to_complex_matrix! (generic function with 1 method)

julia> function quaternion_matrix_to_complex(Q::AbstractMatrix{<:Quaternion})
           m, n = size(Q)
           Z = similar(Q, complex(real(eltype(Q))), 2m, 2n)
           for j in axes(Q, 2), i in axes(Q, 1)
               z = @views Z[2i-1:2i, 2j-1:2j]
               quaternion_to_complex_matrix!(z, Q[i, j])
           end
           return Z
       end
quaternion_matrix_to_complex (generic function with 1 method)

julia> function complex_matrix_to_quaternion(z::AbstractMatrix)
           w = (z[begin, begin] + conj(z[begin + 1, begin + 1])) / 2
           y = (-conj(z[begin+1, begin]) + z[begin, begin+1]) / 2
           return Quaternion(reim(w)..., reim(y)...)
       end
complex_matrix_to_quaternion (generic function with 1 method)

julia> function complex_matrix_to_quaternion_matrix(Z::AbstractMatrix{<:Complex})
           m, n = map(x -> div(x, 2), size(Z))
           Q = similar(Z, Quaternion{real(eltype(Z))}, m, n)
           for j in axes(Q, 2), i in axes(Q, 1)
               z = @views Z[2i-1:2i, 2j-1:2j]
               Q[i, j] = complex_matrix_to_quaternion(z)
           end
           return Q
       end
complex_matrix_to_quaternion_matrix (generic function with 1 method)

julia> q = randn(QuaternionF64, 1, 1);

julia> function matfun(f, Q::AbstractMatrix{<:Quaternion})
           Z = quaternion_matrix_to_complex(Q)
           fZ = f(Z)
           fQ = complex_matrix_to_quaternion_matrix(fZ)
           return fQ
       end
matfun (generic function with 1 method)

julia> matfun(exp, q)[1, 1] ≈ exp(q[1, 1])
true

julia> matfun(log, q)[1, 1] ≈ log(q[1, 1])
true

julia> matfun(cosh, q)[1, 1] ≈ cosh(q[1, 1])
true

julia> q = randn(QuaternionF64, 5, 5);

julia> matfun(log, matfun(exp, q)) ≈ q
true

julia> matfun(inv, q) * q ≈ one(q)
true

julia> matfun(sin, matfun(asin, q)) ≈ q
true

Why do this? Only reason I know of is for the matrix exp/log, which one can use for the Unitary group on quaternions, see e.g. JuliaManifolds/Manifolds.jl#382 (comment). I don't think this would be commonly used though, but it's nice to support operations like this generically here when we can for first-class Quaternion support.