Geodesic flow on the general linear group

geodesic_flow_gl.jl

# Using Manifolds.jl, follow-up to https://github.com/JuliaManifolds/Manifolds.jl/pull/249
using LinearAlgebra
using Manifolds
using Test

# parallel transport of a point (p,X) ∈ T GL(n) along the
# geodesic curve γ on GL(n) with the left-GL(n)-invariant metric
# from section A.2 of https://arxiv.org/abs/1603.05868v1

function geodesic_flow!(M::GeneralLinear, q, Y, p, X, t::Real = 1)
    tX = (Y .= t .* X)
    r = group_exp!(M, q, tX)
    if Manifolds.isnormal(X; atol=sqrt(eps(real(eltype(X)))))
        return compose!(M, q, p, r), copyto!(Y, X)
    end
    K = (tX .-= tX')
    u = exp(K) # u ∈ SU(n,𝔽)
    compose!(M, q, r', u)
    compose!(M, q, p, q)
    mul!(Y, u', X * u)
    return q, Y
end

function geodesic_flow(M::GeneralLinear, p, X, t::Real = 1)
    q = allocate_result(M, exp, p, X, t)
    Y = allocate_result(M, log, X, p, t)
    return geodesic_flow!(M, q, Y, p, X, t)
end

using NLsolve

# Test for GL(3)
M = GeneralLinear(3)
p = exp(randn(3, 3))
X = randn(3, 3)
q, Y = geodesic_flow(M, p, X)
@test q ≈ exp(M, p, X)
@test Y ≈ vector_transport_direction(M, p, X, X, SchildsLadderTransport())

# Test for GL(3, ℂ)
M = GeneralLinear(3, ℂ)
p = exp(randn(ComplexF64, 3, 3))
X = randn(ComplexF64, 3, 3)
q, Y = geodesic_flow(M, p, X)
@test q ≈ exp(M, p, X)
@test Y ≈ vector_transport_direction(M, p, X, X, SchildsLadderTransport())

# Test that it works for SL(3) as a subgroup of GL(3)
M = GeneralLinear(3)
p = project(SpecialLinear(3), randn(3, 3))
X = project(SpecialLinear(3), p, randn(3, 3))
q, Y = geodesic_flow(M, p, X)
@test q ≈ exp(M, p, X)
@test det(q) ≈ 1
@test Y ≈ vector_transport_direction(M, p, X, X, SchildsLadderTransport())
@test tr(Y) ≈ 0 atol=sqrt(eps())