diff --git a/Project.toml b/Project.toml index 1c78847ba3..b2c2ea75eb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Manifolds" uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e" authors = ["Seth Axen ", "Mateusz Baran ", "Ronny Bergmann ", "Antoine Levitt "] -version = "0.4.17" +version = "0.4.19" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" @@ -48,6 +48,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Gtk = "4c0ca9eb-093a-5379-98c5-f87ac0bbbf44" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" @@ -59,4 +60,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" VisualRegressionTests = "34922c18-7c2a-561c-bac1-01e79b2c4c92" [targets] -test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "VisualRegressionTests"] +test = ["Test", "Colors", "DoubleFloats", "FiniteDiff", "ForwardDiff", "Gtk", "ImageIO", "ImageMagick", "OrdinaryDiffEq", "NLsolve", "Plots", "PyPlot", "Quaternions", "QuartzImageIO", "RecipesBase", "ReverseDiff", "VisualRegressionTests"] diff --git a/docs/src/manifolds/group.md b/docs/src/manifolds/group.md index 04fa13cad0..207505e4d9 100644 --- a/docs/src/manifolds/group.md +++ b/docs/src/manifolds/group.md @@ -44,6 +44,22 @@ Pages = ["groups/circle_group.jl"] Order = [:type, :function] ``` +### General linear group + +```@autodocs +Modules = [Manifolds] +Pages = ["groups/general_linear.jl"] +Order = [:type, :function] +``` + +### Special linear group + +```@autodocs +Modules = [Manifolds] +Pages = ["groups/special_linear.jl"] +Order = [:type, :function] +``` + ### Special orthogonal group ```@autodocs diff --git a/docs/src/manifolds/stiefel.md b/docs/src/manifolds/stiefel.md index a23ef419a5..6227e5c53c 100644 --- a/docs/src/manifolds/stiefel.md +++ b/docs/src/manifolds/stiefel.md @@ -19,13 +19,6 @@ Order = [:function] ``` ## The canonical metric - -```@autodocs -Modules = [Manifolds] -Pages = ["manifolds/StiefelCanonicalMetric.jl"] -Order = [:type] -``` - Any ``X∈T_p\mathcal M``, ``p∈\mathcal M``, can be written as ```math @@ -36,9 +29,15 @@ A ∈ ℝ^{p×p} \text{ skew-symmetric}, B ∈ ℝ^{n×p} \text{ arbitrary.} ``` -In the Euclidean case, the elements from ``A`` are counted twice (i.e. weighted with a factor of 2). +In the [`EuclideanMetric`](@ref), the elements from ``A`` are counted twice (i.e. weighted with a factor of 2). The canonical metric avoids this. +```@autodocs +Modules = [Manifolds] +Pages = ["manifolds/StiefelCanonicalMetric.jl"] +Order = [:type] +``` + ```@autodocs Modules = [Manifolds] Pages = ["manifolds/StiefelCanonicalMetric.jl"] diff --git a/docs/src/misc/internals.md b/docs/src/misc/internals.md index f45dd15c21..56184ee871 100644 --- a/docs/src/misc/internals.md +++ b/docs/src/misc/internals.md @@ -7,10 +7,16 @@ This page documents the internal types and methods of `Manifolds.jl`'s that migh ```@docs Manifolds.eigen_safe Manifolds.find_pv +Manifolds.isnormal Manifolds.log_safe +Manifolds.log_safe! +Manifolds.mul!_safe Manifolds.nzsign +Manifolds.realify +Manifolds.realify! Manifolds.size_to_tuple Manifolds.select_from_tuple +Manifolds.unrealify! Manifolds.usinc Manifolds.usinc_from_cos Manifolds.vec2skew! diff --git a/src/Manifolds.jl b/src/Manifolds.jl index a19e9a10a7..0099ecab61 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -132,6 +132,8 @@ include("manifolds/VectorBundle.jl") include("distributions.jl") include("projected_distribution.jl") +include("approx_inverse_retraction.jl") + # It's included early to ensure visibility of `Identity` include("groups/group.jl") @@ -185,6 +187,8 @@ include("groups/array_manifold.jl") include("groups/product_group.jl") include("groups/semidirect_product_group.jl") +include("groups/general_linear.jl") +include("groups/special_linear.jl") include("groups/translation_group.jl") include("groups/special_orthogonal.jl") include("groups/circle_group.jl") @@ -233,6 +237,11 @@ function __init__() include("ode.jl") end + @require NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" begin + using .NLsolve: NLsolve + include("nlsolve.jl") + end + @require Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" begin using .Test: Test include("tests/tests_general.jl") @@ -350,6 +359,7 @@ export AbstractRetractionMethod, ProductRetraction, PowerRetraction export AbstractInverseRetractionMethod, + ApproximateInverseRetraction, LogarithmicInverseRetraction, QRInverseRetraction, PolarInverseRetraction, @@ -475,6 +485,7 @@ export AbstractGroupAction, ActionDirection, AdditionOperation, CircleGroup, + GeneralLinear, GroupManifold, GroupOperationAction, Identity, @@ -489,6 +500,7 @@ export AbstractGroupAction, RotationAction, SemidirectProductGroup, SpecialEuclidean, + SpecialLinear, SpecialOrthogonal, TranslationGroup, TranslationAction diff --git a/src/approx_inverse_retraction.jl b/src/approx_inverse_retraction.jl new file mode 100644 index 0000000000..741b9365f5 --- /dev/null +++ b/src/approx_inverse_retraction.jl @@ -0,0 +1,104 @@ +""" + ApproximateInverseRetraction <: AbstractInverseRetractionMethod + +An abstract type for representing approximate inverse retraction methods. +""" +abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end + +""" + NLsolveInverseRetraction{T<:AbstractRetractionMethod,TV,TK} <: + ApproximateInverseRetraction + +An inverse retraction method for approximating the inverse of a retraction using `NLsolve`. + +# Constructor + + NLsolveInverseRetraction( + method::AbstractRetractionMethod[, X0]; + project_tangent=false, + project_point=false, + nlsolve_kwargs..., + ) + +Constructs an approximate inverse retraction for the retraction `method` with initial guess +`X0`, defaulting to the zero vector. If `project_tangent` is `true`, then the tangent +vector is projected before the retraction using `project`. If `project_point` is `true`, +then the resulting point is projected after the retraction. `nlsolve_kwargs` are keyword +arguments passed to `NLsolve.nlsolve`. +""" +struct NLsolveInverseRetraction{TR<:AbstractRetractionMethod,TV,TK} <: + ApproximateInverseRetraction + retraction::TR + X0::TV + project_tangent::Bool + project_point::Bool + nlsolve_kwargs::TK + function NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs) + isdefined(Manifolds, :NLsolve) || + @warn "To use NLsolveInverseRetraction, NLsolve must be loaded using `using NLsolve`." + return new{typeof(m),typeof(X0),typeof(nlsolve_kwargs)}( + m, + X0, + project_point, + project_tangent, + nlsolve_kwargs, + ) + end +end +function NLsolveInverseRetraction( + m, + X0=nothing; + project_tangent::Bool=false, + project_point::Bool=false, + nlsolve_kwargs..., +) + return NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs) +end + +@decorator_transparent_signature inverse_retract( + M::AbstractDecoratorManifold, + p, + q, + m::NLsolveInverseRetraction, +) + +@decorator_transparent_signature inverse_retract!( + M::AbstractDecoratorManifold, + X, + p, + q, + m::NLsolveInverseRetraction, +) + +@doc raw""" + inverse_retract(M, p, q method::NLsolveInverseRetraction; kwargs...) + +Approximate the inverse of the retraction specified by `method.retraction` from `p` with +respect to `q` on the [`Manifold`](@ref) `M` using NLsolve. This inverse retraction is +not guaranteed to succeed and probably will not unless `q` is close to `p` and the initial +guess `X0` is close. + +If the solver fails to converge, an [`OutOfInjectivityRadiusError`](@ref) is raised. +See [`NLsolveInverseRetraction`](@ref) for configurable parameters. +""" +inverse_retract(::Manifold, p, q, ::NLsolveInverseRetraction; kwargs...) + +function inverse_retract!(M::Manifold, X, p, q, method::NLsolveInverseRetraction; kwargs...) + X0 = method.X0 === nothing ? zero_tangent_vector(M, p) : method.X0 + res = _inverse_retract_nlsolve( + M, + p, + q, + method.retraction, + X0, + method.project_tangent, + method.project_point, + method.nlsolve_kwargs; + kwargs..., + ) + if !res.f_converged + @debug res + throw(OutOfInjectivityRadiusError()) + end + return copyto!(X, res.zero) +end diff --git a/src/groups/circle_group.jl b/src/groups/circle_group.jl index 11adae0897..46e46265f0 100644 --- a/src/groups/circle_group.jl +++ b/src/groups/circle_group.jl @@ -79,4 +79,4 @@ function group_log(G::CircleGroup, q) end end -group_log!(G::CircleGroup, X, q) = (X .= group_log(G, q)) +group_log!(G::CircleGroup, X::AbstractVector, q::AbstractVector) = (X .= group_log(G, q)) diff --git a/src/groups/general_linear.jl b/src/groups/general_linear.jl new file mode 100644 index 0000000000..463bbe0c88 --- /dev/null +++ b/src/groups/general_linear.jl @@ -0,0 +1,236 @@ +@doc raw""" + GeneralLinear{n,𝔽} <: + AbstractGroupManifold{𝔽,MultiplicationOperation,DefaultEmbeddingType} + +The general linear group, that is, the group of all invertible matrices in ``𝔽^{n×n}``. + +The default metric is the left-``\mathrm{GL}(n)``-right-``\mathrm{O}(n)``-invariant metric +whose inner product is +```math +⟨X_p,Y_p⟩_p = ⟨p^{-1}X_p,p^{-1}Y_p⟩_\mathrm{F} = ⟨X_e, Y_e⟩_\mathrm{F}, +``` +where ``X_p, Y_p ∈ T_p \mathrm{GL}(n, 𝔽)``, +``X_e = p^{-1}X_p ∈ 𝔤𝔩(n) = T_e \mathrm{GL}(n, 𝔽) = 𝔽^{n×n}`` is the corresponding +vector in the Lie algebra, and ``⟨⋅,⋅⟩_\mathrm{F}`` denotes the Frobenius inner product. + +By default, tangent vectors ``X_p`` are represented with their corresponding Lie algebra +vectors ``X_e = p^{-1}X_p``. +""" +struct GeneralLinear{n,𝔽} <: + AbstractGroupManifold{𝔽,MultiplicationOperation,DefaultEmbeddingType} end + +GeneralLinear(n, 𝔽::AbstractNumbers=ℝ) = GeneralLinear{n,𝔽}() + +function allocation_promotion_function(::GeneralLinear{n,ℂ}, f, ::Tuple) where {n} + return complex +end + +function check_manifold_point(G::GeneralLinear, p; kwargs...) + mpv = check_manifold_point(decorated_manifold(G), p; kwargs...) + mpv === nothing || return mpv + detp = det(p) + if iszero(detp) + return DomainError( + detp, + "The matrix $(p) does not lie on $(G), since it is not invertible.", + ) + end + return nothing +end +check_manifold_point(::GT, ::Identity{GT}; kwargs...) where {GT<:GeneralLinear} = nothing +function check_manifold_point(G::GeneralLinear, e::Identity; kwargs...) + return DomainError(e, "The identity element $(e) does not belong to $(G).") +end + +function check_tangent_vector(G::GeneralLinear, p, X; check_base_point=true, kwargs...) + if check_base_point + mpe = check_manifold_point(G, p; kwargs...) + mpe === nothing || return mpe + end + mpv = check_tangent_vector(decorated_manifold(G), p, X; kwargs...) + mpv === nothing || return mpv + return nothing +end + +decorated_manifold(::GeneralLinear{n,𝔽}) where {n,𝔽} = Euclidean(n, n; field=𝔽) + +default_metric_dispatch(::GeneralLinear, ::EuclideanMetric) = Val(true) +default_metric_dispatch(::GeneralLinear, ::LeftInvariantMetric{EuclideanMetric}) = Val(true) + +distance(G::GeneralLinear, p, q) = norm(G, p, log(G, p, q)) + +@doc raw""" + exp(G::GeneralLinear, p, X) + +Compute the exponential map on the [`GeneralLinear`](@ref) group. + +The exponential map is +````math +\exp_p \colon X ↦ p \operatorname{Exp}(X^\mathrm{H}) \operatorname{Exp}(X - X^\mathrm{H}), +```` + +where ``\operatorname{Exp}(⋅)`` denotes the matrix exponential, and ``⋅^\mathrm{H}`` is +the conjugate transpose. [^AndruchowLarotondaRechtVarela2014][^MartinNeff2016] + +[^AndruchowLarotondaRechtVarela2014]: + > Andruchow E., Larotonda G., Recht L., and Varela A.: + > “The left invariant metric in the general linear group”, + > Journal of Geometry and Physics 86, pp. 241-257, 2014. + > doi: [10.1016/j.geomphys.2014.08.009](https://doi.org/10.1016/j.geomphys.2014.08.009), + > arXiv: [1109.0520v1](https://arxiv.org/abs/1109.0520v1). +[^MartinNeff2016]: + > Martin, R. J. and Neff, P.: + > “Minimal geodesics on GL(n) for left-invariant, right-O(n)-invariant Riemannian metrics”, + > Journal of Geometric Mechanics 8(3), pp. 323-357, 2016. + > doi: [10.3934/jgm.2016010](https://doi.org/10.3934/jgm.2016010), + > arXiv: [1409.7849v2](https://arxiv.org/abs/1409.7849v2). +""" +exp(::GeneralLinear, p, X) + +function exp!(G::GeneralLinear, q, p, X) + expX = exp(X) + if isnormal(X; atol=sqrt(eps(real(eltype(X))))) + return compose!(G, q, p, expX) + end + compose!(G, q, expX', exp(X - X')) + compose!(G, q, p, q) + return q +end +function exp!(::GeneralLinear{1}, q, p, X) + p1 = p isa Identity ? p : p[1] + q[1] = p1 * exp(X[1]) + return q +end +function exp!(G::GeneralLinear{2}, q, p, X) + if isnormal(X; atol=sqrt(eps(real(eltype(X))))) + return compose!(G, q, p, exp(SizedMatrix{2,2}(X))) + end + A = SizedMatrix{2,2}(X') + B = SizedMatrix{2,2}(X) - A + compose!(G, q, exp(A), exp(B)) + compose!(G, q, p, q) + return q +end + +flat!(::GeneralLinear, ξ::CoTFVector, p, X::TFVector) = copyto!(ξ, X) + +get_coordinates(::GeneralLinear{n,ℝ}, p, X, ::DefaultOrthonormalBasis) where {n} = vec(X) + +function get_coordinates!( + ::GeneralLinear{n,ℝ}, + Xⁱ, + p, + X, + ::DefaultOrthonormalBasis, +) where {n} + return copyto!(Xⁱ, X) +end + +function get_vector(::GeneralLinear{n,ℝ}, p, Xⁱ, ::DefaultOrthonormalBasis) where {n} + return reshape(Xⁱ, n, n) +end + +function get_vector!(::GeneralLinear{n,ℝ}, X, p, Xⁱ, ::DefaultOrthonormalBasis) where {n} + return copyto!(X, Xⁱ) +end + +function group_exp!(::GeneralLinear{1}, q, X) + q[1] = exp(X[1]) + return q +end +group_exp!(::GeneralLinear{2}, q, X) = copyto!(q, exp(SizedMatrix{2,2}(X))) + +function group_log!(::GeneralLinear{1}, X::AbstractMatrix, p::AbstractMatrix) + X[1] = log(p[1]) + return X +end + +inner(::GeneralLinear, p, X, Y) = dot(X, Y) + +invariant_metric_dispatch(::GeneralLinear, ::LeftAction) = Val(true) + +inverse_translate_diff(::GeneralLinear, p, q, X, ::LeftAction) = X +inverse_translate_diff(::GeneralLinear, p, q, X, ::RightAction) = p * X / p + +function inverse_translate_diff!(G::GeneralLinear, Y, p, q, X, conv::ActionDirection) + return copyto!(Y, inverse_translate_diff(G, p, q, X, conv)) +end + +# find sU for s ∈ S⁺ and U ∈ U(n, 𝔽) that minimizes ‖sU - p‖² +function _project_Un_S⁺(p) + n = LinearAlgebra.checksquare(p) + F = svd(p) + s = mean(F.S) + U = F.U * F.Vt + return rmul!(U, s) +end + +@doc raw""" + log(G::GeneralLinear, p, q) + +Compute the logarithmic map on the [`GeneralLinear(n)`](@ref) group. + +The algorithm proceeds in two stages. First, the point ``r = p^{-1} q`` is projected to the +nearest element (under the Frobenius norm) of the direct product subgroup +``\mathrm{O}(n) × S^+``, whose logarithmic map is exactly computed using the matrix +logarithm. This initial tangent vector is then refined using the +[`NLsolveInverseRetraction`](@ref). + +For `GeneralLinear(n, ℂ)`, the logarithmic map is instead computed on the realified +supergroup `GeneralLinear(2n)` and the resulting tangent vector is then complexified. + +Note that this implementation is experimental. +""" +log(::GeneralLinear, p, q) + +function log!(G::GeneralLinear{n,𝔽}, X, p, q) where {n,𝔽} + pinvq = inverse_translate(G, p, q, LeftAction()) + 𝔽 === ℝ && det(pinvq) ≤ 0 && throw(OutOfInjectivityRadiusError()) + e = Identity(G, pinvq) + if isnormal(pinvq; atol=sqrt(eps(real(eltype(pinvq))))) + log_safe!(X, pinvq) + else + # compute the equivalent logarithm on GL(dim(𝔽) * n, ℝ) + # this is significantly more stable than computing the complex algorithm + Gᵣ = GeneralLinear(real_dimension(𝔽) * n, ℝ) + pinvqᵣ = realify(pinvq, 𝔽) + Xᵣ = realify(X, 𝔽) + eᵣ = Identity(Gᵣ, pinvqᵣ) + log_safe!(Xᵣ, _project_Un_S⁺(pinvqᵣ)) + inverse_retraction = NLsolveInverseRetraction(ExponentialRetraction(), Xᵣ) + inverse_retract!(Gᵣ, Xᵣ, eᵣ, pinvqᵣ, inverse_retraction) + unrealify!(X, Xᵣ, 𝔽, n) + end + translate_diff!(G, X, p, e, X, LeftAction()) + return X +end +function log!(::GeneralLinear{1}, X, p, q) + p1 = p isa Identity ? p : p[1] + X[1] = log(p1 \ q[1]) + return X +end + +manifold_dimension(G::GeneralLinear) = manifold_dimension(decorated_manifold(G)) + +LinearAlgebra.norm(::GeneralLinear, p, X) = norm(X) + +project(::GeneralLinear, p) = p +project(::GeneralLinear, p, X) = X + +project!(::GeneralLinear, q, p) = copyto!(q, p) +project!(::GeneralLinear, Y, p, X) = copyto!(Y, X) + +sharp!(::GeneralLinear, X::TFVector, p, ξ::CoTFVector) = copyto!(X, ξ) + +Base.show(io::IO, ::GeneralLinear{n,𝔽}) where {n,𝔽} = print(io, "GeneralLinear($n, $𝔽)") + +translate_diff(::GeneralLinear, p, q, X, ::LeftAction) = X +translate_diff(::GeneralLinear, p, q, X, ::RightAction) = p \ X * p + +function translate_diff!(G::GeneralLinear, Y, p, q, X, conv::ActionDirection) + return copyto!(Y, translate_diff(G, p, q, X, conv)) +end + +vector_transport_to(::GeneralLinear, p, X, q, ::ParallelTransport) = X + +vector_transport_to!(::GeneralLinear, Y, p, X, q, ::ParallelTransport) = copyto!(Y, X) diff --git a/src/groups/group.jl b/src/groups/group.jl index 39ed1106a1..af4bf80240 100644 --- a/src/groups/group.jl +++ b/src/groups/group.jl @@ -350,15 +350,10 @@ end manifold_dimension(G::GroupManifold) = manifold_dimension(G.manifold) function check_manifold_point(G::AbstractGroupManifold, e::Identity; kwargs...) + e.group === G && return nothing return DomainError(e, "The identity element $(e) does not belong to $(G).") end -function check_manifold_point( - G::GT, - e::Identity{GT}; - kwargs..., -) where {GT<:AbstractGroupManifold} - return nothing -end + ########################## # Group-specific functions ########################## @@ -1045,8 +1040,7 @@ inv!(G::MultiplicationGroup, q, p) = copyto!(q, inv(G, p)) compose(::MultiplicationGroup, p, q) = p * q -# TODO: x might alias with p or q, we might be able to optimize it if it doesn't. -compose!(::MultiplicationGroup, x, p, q) = copyto!(x, p * q) +compose!(::MultiplicationGroup, x, p, q) = mul!_safe(x, p, q) inverse_translate(::MultiplicationGroup, p, q, ::LeftAction) = p \ q inverse_translate(::MultiplicationGroup, p, q, ::RightAction) = q / p @@ -1061,3 +1055,5 @@ function group_exp!(G::MultiplicationGroup, q, X) "group_exp! not implemented on $(typeof(G)) for vector $(typeof(X)) and element $(typeof(q)).", ) end + +group_log!(::MultiplicationGroup, X::AbstractMatrix, q::AbstractMatrix) = log_safe!(X, q) diff --git a/src/groups/special_linear.jl b/src/groups/special_linear.jl new file mode 100644 index 0000000000..56a0f9225b --- /dev/null +++ b/src/groups/special_linear.jl @@ -0,0 +1,146 @@ +@doc raw""" + SpecialLinear{n,𝔽} <: + AbstractGroupManifold{𝔽,MultiplicationOperation,DefaultEmbeddingType} + +The special linear group ``\mathrm{SL}(n,𝔽)`` that is, the group of all invertible matrices +with unit determinant in ``𝔽^{n×n}``. + +The Lie algebra ``𝔰𝔩(n, 𝔽) = T_e \mathrm{SL}(n,𝔽)`` is the set of all matrices in +``𝔽^{n×n}`` with trace of zero. By default, tangent vectors ``X_p ∈ T_p \mathrm{SL}(n,𝔽)`` +for ``p ∈ \mathrm{SL}(n,𝔽)`` are represented with their corresponding Lie algebra vector +``X_e = p^{-1}X_p ∈ 𝔰𝔩(n, 𝔽)``. + +The default metric is the same left-``\mathrm{GL}(n)``-right-``\mathrm{O}(n)``-invariant +metric used for [`GeneralLinear(n, 𝔽)`](@ref). The resulting geodesic on +``\mathrm{GL}(n,𝔽)`` emanating from an element of ``\mathrm{SL}(n,𝔽)`` in the direction of +an element of ``𝔰𝔩(n, 𝔽)`` is a closed subgroup of ``\mathrm{SL}(n,𝔽)``. As a result, most +metric functions forward to `GeneralLinear`. +""" +struct SpecialLinear{n,𝔽} <: + AbstractGroupManifold{𝔽,MultiplicationOperation,TransparentIsometricEmbedding} end + +SpecialLinear(n, 𝔽::AbstractNumbers=ℝ) = SpecialLinear{n,𝔽}() + +function allocation_promotion_function(::SpecialLinear{n,ℂ}, f, args::Tuple) where {n} + return complex +end + +function check_manifold_point(G::SpecialLinear{n,𝔽}, p; kwargs...) where {n,𝔽} + mpv = check_manifold_point(Euclidean(n, n; field=𝔽), p; kwargs...) + mpv === nothing || return mpv + detp = det(p) + if !isapprox(detp, 1; kwargs...) + return DomainError( + detp, + "The matrix $(p) does not lie on $(G), since it does not have a unit " * + "determinant.", + ) + end + return nothing +end +check_manifold_point(::GT, ::Identity{GT}; kwargs...) where {GT<:SpecialLinear} = nothing +function check_manifold_point(G::SpecialLinear, e::Identity; kwargs...) + return DomainError(e, "The identity element $(e) does not belong to $(G).") +end + +function check_tangent_vector(G::SpecialLinear, p, X; check_base_point=true, kwargs...) + if check_base_point + mpe = check_manifold_point(G, p; kwargs...) + mpe === nothing || return mpe + end + mpv = + check_tangent_vector(decorated_manifold(G), p, X; check_base_point=false, kwargs...) + mpv === nothing || return mpv + trX = tr(inverse_translate_diff(G, p, p, X, LeftAction())) + if !isapprox(trX, 0; kwargs...) + return DomainError( + trX, + "The matrix $(X) does not lie in the tangent space of $(G) at $(p), since " * + "its Lie algebra representation is not traceless.", + ) + end + return nothing +end + +decorated_manifold(::SpecialLinear{n,𝔽}) where {n,𝔽} = GeneralLinear(n, 𝔽) + +default_metric_dispatch(::SpecialLinear, ::EuclideanMetric) = Val(true) +default_metric_dispatch(::SpecialLinear, ::LeftInvariantMetric{EuclideanMetric}) = Val(true) + +inverse_translate_diff(::SpecialLinear, p, q, X, ::LeftAction) = X +inverse_translate_diff(::SpecialLinear, p, q, X, ::RightAction) = p * X / p + +function inverse_translate_diff!(G::SpecialLinear, Y, p, q, X, conv::ActionDirection) + return copyto!(Y, inverse_translate_diff(G, p, q, X, conv)) +end + +function manifold_dimension(G::SpecialLinear) + return manifold_dimension(decorated_manifold(G)) - real_dimension(number_system(G)) +end + +@doc raw""" + project(G::SpecialLinear, p) + +Project ``p ∈ \mathrm{GL}(n, 𝔽)`` to the [`SpecialLinear`](@ref) group +``G=\mathrm{SL}(n, 𝔽)``. + +Given the singular value decomposition of ``p``, written ``p = U S V^\mathrm{H}``, the +formula for the projection is + +````math +\operatorname{proj}_{\mathrm{SL}(n, 𝔽)}(p) = U S D V^\mathrm{H}, +```` +where + +````math +D_{ij} = δ_{ij} \begin{cases} + 1 & \text{ if } i ≠ n \\ + \det(p)^{-1} & \text{ if } i = n +\end{cases}. +```` +""" +project(::SpecialLinear, p) + +function project!(::SpecialLinear{n}, q, p) where {n} + detp = det(p) + isapprox(detp, 1) && return copyto!(q, p) + F = svd(p) + q .= F.U .* F.S' + q[:, n] ./= detp + mul!_safe(q, q, F.Vt) + return q +end + +@doc raw""" + project(G::SpecialLinear, p, X) + +Orthogonally project ``X ∈ 𝔽^{n × n}`` onto the tangent space of ``p`` to the +[`SpecialLinear`](@ref) ``G = \mathrm{SL}(n, 𝔽)``. The formula reads +````math +\operatorname{proj}_{p} + = (\mathrm{d}L_p)_e ∘ \operatorname{proj}_{𝔰𝔩(n, 𝔽)} ∘ (\mathrm{d}L_p^{-1})_p + \colon X ↦ X - \frac{\operatorname{tr}(X)}{n} I, +```` +where the last expression uses the tangent space representation as the Lie algebra. +""" +project(::SpecialLinear, p, X) + +function project!(G::SpecialLinear{n}, Y, p, X) where {n} + inverse_translate_diff!(G, Y, p, p, X, LeftAction()) + Y[diagind(n, n)] .-= tr(Y) / n + translate_diff!(G, Y, p, p, Y, LeftAction()) + return Y +end + +function decorator_transparent_dispatch(::typeof(project), ::SpecialLinear, args...) + return Val(:parent) +end + +Base.show(io::IO, ::SpecialLinear{n,𝔽}) where {n,𝔽} = print(io, "SpecialLinear($n, $𝔽)") + +translate_diff(::SpecialLinear, p, q, X, ::LeftAction) = X +translate_diff(::SpecialLinear, p, q, X, ::RightAction) = p \ X * p + +function translate_diff!(G::SpecialLinear, Y, p, q, X, conv::ActionDirection) + return copyto!(Y, translate_diff(G, p, q, X, conv)) +end diff --git a/src/groups/special_orthogonal.jl b/src/groups/special_orthogonal.jl index 334c90ee3e..bcac52c625 100644 --- a/src/groups/special_orthogonal.jl +++ b/src/groups/special_orthogonal.jl @@ -44,6 +44,9 @@ end group_exp!(G::SpecialOrthogonal, q, X) = exp!(G, q, make_identity(G, q).p, X) group_log!(G::SpecialOrthogonal, X, q) = log!(G, X, make_identity(G, q).p, q) +function group_log!(G::SpecialOrthogonal, X::AbstractMatrix, q::AbstractMatrix) + return log!(G, X, make_identity(G, q).p, q) +end function allocate_result( ::GT, diff --git a/src/manifolds/StiefelCanonicalMetric.jl b/src/manifolds/StiefelCanonicalMetric.jl index aa8e2149f6..ac986edbba 100644 --- a/src/manifolds/StiefelCanonicalMetric.jl +++ b/src/manifolds/StiefelCanonicalMetric.jl @@ -1,8 +1,8 @@ @doc raw""" CanonicalMetric <: Metric -The Canonical Metric is name used on several manifolds, for example for the [`Stiefel`](@ref) -manifold, see for example[^EdelmanAriasSmith1998]. +The Canonical Metric refers to a metric for the [`Stiefel`](@ref) +manifold, see[^EdelmanAriasSmith1998]. [^EdelmanAriasSmith1998]: > Edelman, A., Ariar, T. A., Smith, S. T.: @@ -60,7 +60,10 @@ function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) w n == k && return mul!(q, p, exp(A)) QR = qr(X - p * A) BC_ext = exp([A -QR.R'; QR.R 0*I]) - q .= [p Matrix(QR.Q)] * @view(BC_ext[:, 1:k]) + @views begin + mul!(q, p, BC_ext[1:k, 1:k]) + mul!(q, Matrix(QR.Q), BC_ext[(k + 1):(2 * k), 1:k], true, true) + end return q end @@ -128,24 +131,26 @@ function log!( ) where {n,k} M = p' * q QR = qr(q - p * M) - V = Matrix(qr([M; QR.R]).Q * I) + V = convert(Matrix, qr([M; QR.R]).Q * I) Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner Vpcols = @view V[:, (k + 1):(2 * k)] #second half of the columns S = svd(Vcorner) # preprocessing: Procrustes mul!(Vpcols, Vpcols * S.U, S.V') - V[:, 1:k] .= [M; QR.R] + V[1:k, 1:k] .= M + V[(k + 1):(2 * k), 1:k] .= QR.R - LV = real.(log(V)) # this can be replaced by log_safe. + LV = log_safe!(allocate(V), V) C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k)) expnC = exp(-C) i = 0 + new_Vpcols = Vpcols * expnC # allocate once while (i < maxiter) && (norm(C) > tolerance) i = i + 1 - LV .= real.(log(V)) - expnC .= exp(-C) - copyto!(Vpcols, Vpcols * expnC) + log_safe!(LV, V) + expnC = exp(-C) + mul!(new_Vpcols, Vpcols, expnC) + copyto!(Vpcols, new_Vpcols) end - LV .= real.(LV) mul!(X, p, @view(LV[1:k, 1:k])) # force the first - Q - to be the reduced form, not the full matrix mul!(X, @view(QR.Q[:, 1:k]), @view(LV[(k + 1):(2 * k), 1:k]), true, true) diff --git a/src/manifolds/StiefelEuclideanMetric.jl b/src/manifolds/StiefelEuclideanMetric.jl index d1d359e7e8..2065da231d 100644 --- a/src/manifolds/StiefelEuclideanMetric.jl +++ b/src/manifolds/StiefelEuclideanMetric.jl @@ -30,7 +30,13 @@ exp(::Stiefel, ::Any...) function exp!(::Stiefel{n,k}, q, p, X) where {n,k} A = p' * X - return copyto!(q, [p X] * exp([A -X'*X; I A]) * [exp(-A); 0 * I]) + B = exp([A -X'*X; I A]) + @views begin + r = p * B[1:k, 1:k] + mul!(r, X, B[(k + 1):(2 * k), 1:k], true, true) + end + mul!(q, r, exp(-A)) + return q end @doc raw""" diff --git a/src/manifolds/SymmetricPositiveDefinite.jl b/src/manifolds/SymmetricPositiveDefinite.jl index 515f00364d..a9daffe529 100644 --- a/src/manifolds/SymmetricPositiveDefinite.jl +++ b/src/manifolds/SymmetricPositiveDefinite.jl @@ -170,4 +170,4 @@ definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`. """ zero_tangent_vector(::SymmetricPositiveDefinite, ::Any) -zero_tangent_vector!(::SymmetricPositiveDefinite{N}, v, x) where {N} = fill!(v, 0) +zero_tangent_vector!(::SymmetricPositiveDefinite, v, ::Any) = fill!(v, 0) diff --git a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl index 6b848a09e2..1d8bf5f5a7 100644 --- a/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl +++ b/src/manifolds/SymmetricPositiveSemidefiniteFixedRank.jl @@ -270,11 +270,6 @@ definite matrix `p` on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref) mani """ zero_tangent_vector(::SymmetricPositiveSemidefiniteFixedRank, ::Any...) -function zero_tangent_vector!( - ::SymmetricPositiveSemidefiniteFixedRank{N,K}, - v, - ::Any, -) where {N,K} - fill!(v, 0) - return v +function zero_tangent_vector!(::SymmetricPositiveSemidefiniteFixedRank, v, ::Any) + return fill!(v, 0) end diff --git a/src/nlsolve.jl b/src/nlsolve.jl new file mode 100644 index 0000000000..092d28146e --- /dev/null +++ b/src/nlsolve.jl @@ -0,0 +1,21 @@ +function _inverse_retract_nlsolve( + M::Manifold, + p, + q, + retraction, + X0, + project_tangent, + project_point, + nlsolve_kwargs; + kwargs..., +) + function f!(F, X) + project_tangent && project!(M, X, p, X) + retract!(M, F, p, project(M, p, X), retraction; kwargs...) + project_point && project!(M, q, q) + F .-= q + return F + end + res = NLsolve.nlsolve(f!, X0; nlsolve_kwargs...) + return res +end diff --git a/src/utils.jl b/src/utils.jl index 367c5d2089..dc8f2d638b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -68,6 +68,138 @@ converted to a `Matrix` before computing the log. return SizedMatrix{s[1],s[2]}(log(Matrix(parent(x)))) end +# NOTE: workaround until https://github.com/JuliaLang/julia/pull/39973 or similar is merged +""" + log_safe!(y, x) + +Compute the matrix logarithm of `x`. If the eltype of `y` is real, then the imaginary part +of `x` is ignored, and a `DomainError` is raised if `real(x)` has no real logarithm. +""" +function log_safe!(Y, A) + if eltype(Y) <: Real + if ishermitian(A) + eigenF = eigen(Symmetric(real(A))) + i = findfirst(≤(0), eigenF.values) + if i !== nothing + throw( + DomainError( + eigenF.values[i], + "All eigenvalues must be positive to compute a real logarithm.", + ), + ) + end + mul!(Y, eigenF.vectors .* log.(eigenF.values'), eigenF.vectors') + elseif istriu(A) + i = findfirst(≤(0), @view(A[diagind(A)])) + if i !== nothing + throw( + DomainError( + A[i, i], + "All eigenvalues must be positive to compute a real logarithm.", + ), + ) + end + copyto!(Y, real(log(UpperTriangular(A)))) + else + schurF = schur(convert(Matrix, real(A))) + i = findfirst(x -> isreal(x) && real(x) ≤ 0, schurF.values) + if i !== nothing + throw( + DomainError( + schurF.values[i], + "All eigenvalues must be positive to compute a real logarithm.", + ), + ) + end + if istriu(schurF.T) + mul!(Y, schurF.Z, real(log(UpperTriangular(schurF.T))) * schurF.Z') + else + schurS = schur(complex(schurF.T)) + Y .= real.(schurS.Z * log(UpperTriangular(schurS.T)) * schurS.Z') + mul!(Y, schurF.Z * Y, schurF.Z') + end + end + else + copyto!(Y, log_safe(A)) + end + return Y +end + +""" + mul!_safe(Y, A, B) -> Y + +Call `mul!` safely, that is, `A` and/or `B` are permitted to alias with `Y`. +""" +mul!_safe(Y, A, B) = (Y === A || Y === B) ? copyto!(Y, A * B) : mul!(Y, A, B) + +@doc raw""" + realify(X::AbstractMatrix{T𝔽}, 𝔽::AbstractNumbers) -> Y::AbstractMatrix{<:Real} + +Given a matrix $X ∈ 𝔽^{n × n}$, compute $Y ∈ ℝ^{m × m}$, where $m = n \operatorname{dim}_𝔽$, +and $\operatorname{dim}_𝔽$ is the [`real_dimension`](@ref) of the number field $𝔽$, using +the map $ϕ \colon X ↦ Y$, that preserves the matrix product, so that for all +$C,D ∈ 𝔽^{n × n}$, +````math +ϕ(C) ϕ(D) = ϕ(CD). +```` +See [`realify!`](@ref) for an in-place version, and [`unrealify!`](@ref) to compute the +inverse of $ϕ$. +""" +function realify(X, 𝔽) + n = LinearAlgebra.checksquare(X) + nℝ = real_dimension(𝔽) * n + Y = allocate(X, real(eltype(X)), nℝ, nℝ) + return realify!(Y, X, 𝔽, n) +end +realify(X, ::typeof(ℝ)) = X + +""" + realify!(Y::AbstractMatrix{<:Real}, X::AbstractMatrix{T𝔽}, 𝔽::AbstractNumbers) + +In-place version of [`realify`](@ref). +""" +realify!(Y, X, 𝔽) + +@doc raw""" + realify!(Y::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Complex}, ::typeof(ℂ)) + +Given a complex matrix $X = A + iB ∈ ℂ^{n × n}$, compute its realified matrix +$Y ∈ ℝ^{2n × 2n}$, written +where +````math +Y = \begin{pmatrix}A & -B \\ B & A \end{pmatrix}. +```` +""" +function realify!(Y, X, ::typeof(ℂ), n=LinearAlgebra.checksquare(X)) + for i in 1:n, j in 1:n + Xr, Xi = reim(X[i, j]) + Y[i, j] = Y[n + i, n + j] = Xr + Y[n + i, j] = Xi + Y[i, n + j] = -Xi + end + return Y +end + +@doc raw""" + unrealify!(X::AbstractMatrix{T𝔽}, Y::AbstractMatrix{<:Real}, 𝔽::AbstractNumbers[, n]) + +Given a real matrix $Y ∈ ℝ^{m × m}$, where $m = n \operatorname{dim}_𝔽$, and +$\operatorname{dim}_𝔽$ is the [`real_dimension`](@ref) of the number field $𝔽$, compute +in-place its equivalent matrix $X ∈ 𝔽^{n × n}$. Note that this function does not check that +$Y$ has a valid structure to be un-realified. + +See [`realify!`](@ref) for the inverse of this function. +""" +unrealify!(X, Y, 𝔽) + +function unrealify!(X, Y, ::typeof(ℂ), n=LinearAlgebra.checksquare(X)) + for i in 1:n, j in 1:n + X[i, j] = complex((Y[i, j] + Y[n + i, n + j]) / 2, (Y[n + i, j] - Y[i, n + j]) / 2) + end + return X +end +unrealify!(Y, X, ::typeof(ℝ), args...) = copyto!(Y, X) + @generated maybesize(s::Size{S}) where {S} = prod(S) > 100 ? S : :(s) """ @@ -113,6 +245,23 @@ function vec2skew(v, k) return X end +@doc raw""" + isnormal(x; kwargs...) -> Bool + +Check if the matrix or number `x` is normal, that is, if it commutes with its adjoint: +````math +x x^\mathrm{H} = x^\mathrm{H} x. +```` +By default, this is an equality check. Provide `kwargs` for `isapprox` to perform an +approximate check. +""" +function isnormal(x; kwargs...) + (isdiag(x) || ishermitian(x)) && return true + isempty(kwargs) && return x * x' == x' * x + return isapprox(x * x', x' * x; kwargs...) +end +isnormal(::LinearAlgebra.RealHermSymComplexHerm; kwargs...) = true + """ ziptuples(a, b[, c[, d[, e]]]) diff --git a/test/approx_inverse_retraction.jl b/test/approx_inverse_retraction.jl new file mode 100644 index 0000000000..78909b3964 --- /dev/null +++ b/test/approx_inverse_retraction.jl @@ -0,0 +1,78 @@ +using NLsolve +using Manifolds: NLsolveInverseRetraction, ApproximateInverseRetraction +using LinearAlgebra + +include("utils.jl") + +@testset "approximate inverse retractions" begin + @testset "NLsolveInverseRetraction" begin + @testset "constructor" begin + X = randn(3) + + @test NLsolveInverseRetraction <: ApproximateInverseRetraction + m1 = NLsolveInverseRetraction(ExponentialRetraction()) + @test m1.retraction === ExponentialRetraction() + @test m1.X0 === nothing + @test !m1.project_tangent + @test !m1.project_point + @test isempty(m1.nlsolve_kwargs) + + m2 = NLsolveInverseRetraction( + PolarRetraction(), + [1.0, 2.0, 3.0]; + project_tangent=true, + project_point=true, + autodiff=:forward, + ) + @test m2.retraction === PolarRetraction() + @test m2.X0 == [1.0, 2.0, 3.0] + @test m2.project_tangent + @test m2.project_point + @test (; m2.nlsolve_kwargs...) == (; autodiff=:forward) + end + @testset "Euclidean" begin + M = Euclidean(3) + p = [1.0, 2.0, 3.0] + q = [4.0, 5.0, 6.0] + retr_method = ExponentialRetraction() + inv_retr_method = NLsolveInverseRetraction(retr_method) + X = inverse_retract(M, p, q, inv_retr_method) + @test is_tangent_vector(M, p, X) + @test X isa Vector{Float64} + @test X ≈ q - p + end + + @testset "Sphere" begin + M = Sphere(2) + p = [1.0, 0.0, 0.0] + q = [1 / sqrt(2), 1 / sqrt(2), 0.0] + X_exp = inverse_retract(M, p, q, ProjectionInverseRetraction()) + # vector must be nonzero to converge + X0 = randn(3) .* eps() + inv_retr_method = + NLsolveInverseRetraction(ProjectionRetraction(), X0; project_point=true) + X = inverse_retract(M, p, q, inv_retr_method) + @test is_tangent_vector(M, p, X; atol=1e-9) + @test X ≈ X_exp + @test_throws OutOfInjectivityRadiusError inverse_retract( + M, + p, + -p, + inv_retr_method, + ) + end + + @testset "Circle(ℂ)" begin + M = Circle(ℂ) + p = [1.0 * im] + X = [p[1] * im * (π / 4)] + q = exp(M, p, X) + X_exp = log(M, p, q) + inv_retr_method = + NLsolveInverseRetraction(ExponentialRetraction(); project_point=true) + X = inverse_retract(M, p, q, inv_retr_method) + @test is_tangent_vector(M, p, X; atol=1e-8) + @test X ≈ X_exp + end + end +end diff --git a/test/groups/general_linear.jl b/test/groups/general_linear.jl new file mode 100644 index 0000000000..94da76e632 --- /dev/null +++ b/test/groups/general_linear.jl @@ -0,0 +1,212 @@ +include("../utils.jl") +include("group_utils.jl") +using NLsolve + +@testset "General Linear group" begin + @testset "basic properties" begin + G = GeneralLinear(3) + @test G === GeneralLinear(3, ℝ) + @test repr(G) == "GeneralLinear(3, ℝ)" + @test base_manifold(G) === GeneralLinear(3) + @test decorated_manifold(G) === Euclidean(3, 3) + @test number_system(G) === ℝ + @test manifold_dimension(G) == 9 + @test representation_size(G) == (3, 3) + Gc = GeneralLinear(2, ℂ) + @test decorated_manifold(Gc) === Euclidean(2, 2; field=ℂ) + @test repr(Gc) == "GeneralLinear(2, ℂ)" + @test number_system(Gc) == ℂ + @test manifold_dimension(Gc) == 8 + @test representation_size(Gc) == (2, 2) + Gh = GeneralLinear(4, ℍ) + @test repr(Gh) == "GeneralLinear(4, ℍ)" + @test number_system(Gh) == ℍ + @test manifold_dimension(Gh) == 4 * 16 + @test representation_size(Gh) == (4, 4) + + @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true) + @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(false) + @test is_default_metric( + MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())), + ) === true + @test @inferred(Manifolds.default_metric_dispatch(G, EuclideanMetric())) === + Val(true) + @test @inferred( + Manifolds.default_metric_dispatch( + G, + InvariantMetric(EuclideanMetric(), LeftAction()), + ) + ) === Val(true) + @test @inferred( + Manifolds.default_metric_dispatch( + MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())), + ) + ) === Val(true) + @test Manifolds.allocation_promotion_function(Gc, exp!, (1,)) === complex + end + + @testset "GL(1,𝔽) special cases" begin + @testset "real" begin + G = GeneralLinear(1) + p = 3.0 * ones(1, 1) + X = 1.0 * ones(1, 1) + @test exp(G, p, X) ≈ p * exp(X)' * exp(X - X') + q = exp(G, p, X) + Y = log(G, p, q) + @test Y ≈ X + @test group_exp(G, X) ≈ exp(X) + @test group_log(G, exp(X)) ≈ X + end + @testset "complex" begin + G = GeneralLinear(1, ℂ) + p = (1 + im) * ones(1, 1) + X = (1 - im) * ones(1, 1) + @test exp(G, p, X) ≈ p * exp(X)' * exp(X - X') + q = exp(G, p, X) + Y = log(G, p, q) + @test Y ≈ X + @test group_exp(G, X) ≈ exp(X) + @test group_log(G, exp(X)) ≈ X + end + end + + @testset "Real" begin + G = GeneralLinear(3) + + @test_throws DomainError is_manifold_point(G, randn(2, 3), true) + @test_throws DomainError is_manifold_point(G, randn(2, 2), true) + @test_throws DomainError is_manifold_point(G, randn(ComplexF64, 3, 3), true) + @test_throws DomainError is_manifold_point(G, zeros(3, 3), true) + @test_throws DomainError is_manifold_point(G, Float64[0 0 0; 0 1 1; 1 1 1], true) + @test_throws DomainError is_manifold_point( + G, + make_identity(GeneralLinear(2), ones(2, 2)), + true, + ) + @test is_manifold_point(G, Float64[0 0 1; 0 1 1; 1 1 1], true) + @test is_manifold_point(G, make_identity(G, ones(3, 3)), true) + @test_throws DomainError is_tangent_vector( + G, + Float64[0 1 1; 0 1 1; 1 0 0], + randn(3, 3), + true, + ) + @test is_tangent_vector(G, Float64[0 0 1; 0 1 1; 1 1 1], randn(3, 3), true) + + types = [Matrix{Float64}] + pts = [ + Matrix(Diagonal([1, 2, 3])), + [-2 5 -5; 0 2 -1; -3 -5 -2], + [-5 1 0; 1 0 1; 0 1 3], + ] + vpts = [[-1 -2 0; -2 1 -2; 2 0 2], [1 1 1; 0 0 -2; 2 0 0]] + + retraction_methods = [ + Manifolds.GroupExponentialRetraction(LeftAction()), + Manifolds.GroupExponentialRetraction(RightAction()), + ] + + inverse_retraction_methods = [ + Manifolds.GroupLogarithmicInverseRetraction(LeftAction()), + Manifolds.GroupLogarithmicInverseRetraction(RightAction()), + ] + + basis_types = [DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd)] + + for T in types + gpts = convert.(T, pts) + vgpts = convert.(T, vpts) + test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_manifold( + G, + gpts; + test_reverse_diff=false, + test_forward_diff=false, + test_project_point=true, + test_injectivity_radius=false, + test_project_tangent=true, + test_musical_isomorphisms=true, + test_default_vector_transport=true, + vector_transport_methods=[ + ParallelTransport(), + SchildsLadderTransport(), + PoleLadderTransport(), + ], + retraction_methods=retraction_methods, + inverse_retraction_methods=inverse_retraction_methods, + basis_types_vecs=basis_types, + basis_types_to_from=basis_types, + exp_log_atol_multiplier=1e7, + retraction_atol_multiplier=1e7, + ) + end + end + + @testset "Complex" begin + G = GeneralLinear(2, ℂ) + + @test_throws DomainError is_manifold_point(G, randn(ComplexF64, 2, 3), true) + @test_throws DomainError is_manifold_point(G, randn(ComplexF64, 3, 3), true) + @test_throws DomainError is_manifold_point(G, zeros(2, 2), true) + @test_throws DomainError is_manifold_point(G, ComplexF64[1 im; 1 im], true) + @test is_manifold_point(G, ComplexF64[1 1; im 1], true) + @test is_manifold_point(G, make_identity(G, ones(ComplexF64, 2, 2)), true) + @test_throws DomainError is_manifold_point(G, Float64[0 0 0; 0 1 1; 1 1 1], true) + @test_throws DomainError is_manifold_point( + G, + make_identity(GeneralLinear(3), ones(3, 3)), + true, + ) + @test_throws DomainError is_tangent_vector( + G, + ComplexF64[im im; im im], + randn(ComplexF64, 2, 2), + true, + ) + @test is_tangent_vector(G, ComplexF64[1 im; im im], randn(ComplexF64, 2, 2), true) + + types = [Matrix{ComplexF64}] + pts = [ + [-1-5im -1+3im; -6-4im 4+6im], + [1+3im -1-4im; -2-2im -3-1im], + [-6+0im 1+1im; 1-1im -4+0im], + ] + vpts = [[1+0im -2-1im; -1-2im -4+1im], [-2+2im -1-1im; -1-1im -3+0im]] + + retraction_methods = [ + Manifolds.GroupExponentialRetraction(LeftAction()), + Manifolds.GroupExponentialRetraction(RightAction()), + ] + + inverse_retraction_methods = [ + Manifolds.GroupLogarithmicInverseRetraction(LeftAction()), + Manifolds.GroupLogarithmicInverseRetraction(RightAction()), + ] + + for T in types + gpts = convert.(T, pts) + vgpts = convert.(T, vpts) + test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_manifold( + G, + gpts; + test_reverse_diff=false, + test_forward_diff=false, + test_project_point=true, + test_injectivity_radius=false, + test_project_tangent=true, + test_musical_isomorphisms=true, + test_default_vector_transport=true, + vector_transport_methods=[ + ParallelTransport(), + SchildsLadderTransport(), + PoleLadderTransport(), + ], + retraction_methods=retraction_methods, + inverse_retraction_methods=inverse_retraction_methods, + exp_log_atol_multiplier=1e8, + retraction_atol_multiplier=1e8, + ) + end + end +end diff --git a/test/groups/groups_general.jl b/test/groups/groups_general.jl index f76f97684d..5c2eb95c79 100644 --- a/test/groups/groups_general.jl +++ b/test/groups/groups_general.jl @@ -237,13 +237,13 @@ include("group_utils.jl") G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation()) test_group( G, - [[1.0 2.0; 3.0 4.0], [2.0 3.0; 4.0 5.0], [3.0 4.0; 5.0 6.0]], + [[2.0 1.0; 3.0 4.0], [3.0 2.0; 4.0 5.0], [4.0 3.0; 5.0 6.0]], [], [[1.0 2.0; 3.0 4.0]]; - test_group_exp_log=false, + test_group_exp_log=true, ) - x = [1.0 2.0; 2.0 3.0] + x = [2.0 1.0; 2.0 3.0] ge = Identity(G, [1.0 0.0; 0.0 1.0]) @test number_eltype(ge) == Bool @test copyto!(ge, ge) === ge @@ -294,8 +294,12 @@ include("group_utils.jl") @test y ≈ x compose!(G, y, ge, x) @test y ≈ x - @test group_exp!(G, y, x) === y - @test y ≈ exp(x) + X = [1.0 2.0; 3.0 4.0] + @test group_exp!(G, y, X) === y + @test y ≈ exp(X) + Y = allocate(X) + @test group_log!(G, Y, y) === Y + @test Y ≈ log(y) @testset "identity optimization" begin x2 = copy(x) diff --git a/test/groups/special_linear.jl b/test/groups/special_linear.jl new file mode 100644 index 0000000000..3125a02bc8 --- /dev/null +++ b/test/groups/special_linear.jl @@ -0,0 +1,241 @@ +include("../utils.jl") +include("group_utils.jl") +using NLsolve + +@testset "Special Linear group" begin + @testset "basic properties" begin + G = SpecialLinear(3) + @test G === SpecialLinear(3, ℝ) + @test repr(G) == "SpecialLinear(3, ℝ)" + @test base_manifold(G) === SpecialLinear(3) + @test decorated_manifold(G) == GeneralLinear(3) + @test number_system(G) === ℝ + @test manifold_dimension(G) == 8 + @test representation_size(G) == (3, 3) + Gc = SpecialLinear(2, ℂ) + @test decorated_manifold(Gc) == GeneralLinear(2, ℂ) + @test repr(Gc) == "SpecialLinear(2, ℂ)" + @test number_system(Gc) == ℂ + @test manifold_dimension(Gc) == 6 + @test representation_size(Gc) == (2, 2) + Gh = SpecialLinear(4, ℍ) + @test repr(Gh) == "SpecialLinear(4, ℍ)" + @test number_system(Gh) == ℍ + @test manifold_dimension(Gh) == 4 * 15 + @test representation_size(Gh) == (4, 4) + + @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true) + @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(false) + @test is_default_metric( + MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())), + ) === true + @test @inferred(Manifolds.default_metric_dispatch(G, EuclideanMetric())) === + Val(true) + @test @inferred( + Manifolds.default_metric_dispatch( + G, + InvariantMetric(EuclideanMetric(), LeftAction()), + ) + ) === Val(true) + @test @inferred( + Manifolds.default_metric_dispatch( + MetricManifold(G, InvariantMetric(EuclideanMetric(), LeftAction())), + ) + ) === Val(true) + @test Manifolds.allocation_promotion_function(Gc, exp!, (1,)) === complex + end + + @testset "Real" begin + G = SpecialLinear(3) + + @test_throws DomainError is_manifold_point(G, randn(2, 3), true) + @test_throws DomainError is_manifold_point(G, Float64[2 1; 1 1], true) + @test_throws DomainError is_manifold_point(G, [1 0 im; im 0 0; 0 -1 0], true) + @test_throws DomainError is_manifold_point(G, zeros(3, 3), true) + @test_throws DomainError is_manifold_point(G, Float64[1 3 3; 1 1 2; 1 2 3], true) + @test_throws DomainError is_manifold_point( + G, + make_identity(SpecialLinear(2), ones(2, 2)), + true, + ) + @test is_manifold_point(G, Float64[1 1 1; 2 2 1; 2 3 3], true) + @test is_manifold_point(G, make_identity(G, ones(3, 3)), true) + @test_throws DomainError is_tangent_vector( + G, + Float64[2 3 2; 3 1 2; 1 1 1], + randn(3, 3), + true; + atol=1e-6, + ) + @test_throws DomainError is_tangent_vector( + G, + Float64[2 1 2; 3 2 2; 2 2 1], + Float64[2 1 -1; 2 2 1; 1 1 -1], + true; + atol=1e-6, + ) + @test is_tangent_vector( + G, + Float64[2 1 2; 3 2 2; 2 2 1], + Float64[-1 -1 -1; 1 -1 2; -1 -1 2], + true; + atol=1e-6, + ) + + types = [Matrix{Float64}] + pts = + [[2 -1 -3; 4 -1 -6; -1 1 2], [0 2 1; 0 -3 -1; 1 0 2], [-2 0 -1; 1 0 0; -1 -1 2]] + vpts = [[0 -1 -5; 1 2 0; 1 2 -2], [0 -2 1; -2 1 2; -4 2 -1]] + + retraction_methods = [ + Manifolds.GroupExponentialRetraction(LeftAction()), + Manifolds.GroupExponentialRetraction(RightAction()), + ] + + inverse_retraction_methods = [ + Manifolds.GroupLogarithmicInverseRetraction(LeftAction()), + Manifolds.GroupLogarithmicInverseRetraction(RightAction()), + ] + + for T in types + gpts = convert.(T, pts) + vgpts = convert.(T, vpts) + test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_manifold( + G, + gpts; + test_reverse_diff=false, + test_forward_diff=false, + test_injectivity_radius=false, + test_project_point=true, + test_project_tangent=true, + test_musical_isomorphisms=true, + test_default_vector_transport=true, + vector_transport_methods=[ + ParallelTransport(), + SchildsLadderTransport(), + PoleLadderTransport(), + ], + retraction_methods=retraction_methods, + inverse_retraction_methods=inverse_retraction_methods, + exp_log_atol_multiplier=1e10, + retraction_atol_multiplier=1e8, + is_tangent_atol_multiplier=1e10, + ) + end + + @testset "project" begin + p = randn(3, 3) + @test !is_manifold_point(G, p) + q = project(G, p) + @test is_manifold_point(G, q) + @test project(G, q) ≈ q + + X = randn(3, 3) + @test !is_tangent_vector(G, q, X; atol=1e-6) + Y = project(G, q, X) + @test is_tangent_vector(G, q, Y; atol=1e-6) + @test project(G, q, Y) ≈ Y + end + end + + @testset "Complex" begin + G = SpecialLinear(2, ℂ) + + @test_throws DomainError is_manifold_point(G, randn(ComplexF64, 2, 3), true) + @test_throws DomainError is_manifold_point(G, randn(2, 2), true) + @test_throws DomainError is_manifold_point( + G, + ComplexF64[1 0 im; im 0 0; 0 -1 0], + true, + ) + @test_throws DomainError is_manifold_point(G, ComplexF64[1 im; im 1], true) + @test_throws DomainError is_manifold_point( + G, + make_identity(SpecialLinear(2), ones(ComplexF64, 2, 2)), + true, + ) + @test is_manifold_point(G, ComplexF64[im 1; -2 im], true) + @test is_manifold_point(G, make_identity(G, ones(3, 3)), true) + @test_throws DomainError is_tangent_vector( + G, + ComplexF64[-1+im -1; -im 1], + ComplexF64[1-im 1+im; 1 -1+im], + true; + atol=1e-6, + ) + @test_throws DomainError is_tangent_vector( + G, + ComplexF64[1 1+im; -1+im -1], + ComplexF64[1-im -1-im; -im im], + true; + atol=1e-6, + ) + @test is_tangent_vector( + G, + ComplexF64[1 1+im; -1+im -1], + ComplexF64[1-im 1+im; 1 -1+im], + true; + atol=1e-6, + ) + + types = [Matrix{ComplexF64}] + pts = [ + [0+1im 0+2im; 0+1im 0+1im], + [-2+0im 0+1im; 0-1im -1+0im], + [0+0im 0-1im; 0-1im -1+3im], + ] + vpts = [[0-1im -1+1im; -1+0im 0+1im], [1+1im 0+0im; 0+1im -1-1im]] + + retraction_methods = [ + Manifolds.GroupExponentialRetraction(LeftAction()), + Manifolds.GroupExponentialRetraction(RightAction()), + ] + + inverse_retraction_methods = [ + Manifolds.GroupLogarithmicInverseRetraction(LeftAction()), + Manifolds.GroupLogarithmicInverseRetraction(RightAction()), + ] + + for T in types + gpts = convert.(T, pts) + vgpts = convert.(T, vpts) + test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true) + test_manifold( + G, + gpts; + test_reverse_diff=false, + test_forward_diff=false, + test_injectivity_radius=false, + test_project_point=true, + test_project_tangent=true, + test_musical_isomorphisms=true, + test_default_vector_transport=true, + vector_transport_methods=[ + ParallelTransport(), + SchildsLadderTransport(), + PoleLadderTransport(), + ], + retraction_methods=retraction_methods, + inverse_retraction_methods=inverse_retraction_methods, + exp_log_atol_multiplier=1e8, + retraction_atol_multiplier=1e8, + is_tangent_atol_multiplier=1e8, + ) + end + + @testset "project" begin + p = randn(ComplexF64, 2, 2) + @test !is_manifold_point(G, p) + q = project(G, p) + @test is_manifold_point(G, q) + @test project(G, q) ≈ q + + X = randn(ComplexF64, 2, 2) + @test !is_tangent_vector(G, q, X; atol=1e-6) + Y = project(G, q, X) + @test is_tangent_vector(G, q, Y; atol=1e-6) + @test project(G, q, Y) ≈ Y + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d1806a0298..afcde0958a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,7 +11,7 @@ include("utils.jl") @testset "Ambiguities" begin # TODO: reduce the number of ambiguities if VERSION.prerelease == () # - @test length(Test.detect_ambiguities(ManifoldsBase)) <= 17 + @test length(Test.detect_ambiguities(ManifoldsBase)) <= 18 @test length(Test.detect_ambiguities(Manifolds)) == 0 @test length(our_base_ambiguities()) <= 24 else @@ -20,8 +20,89 @@ include("utils.jl") end @testset "utils test" begin - @test Manifolds.usinc_from_cos(-1) == 0 - @test Manifolds.usinc_from_cos(-1.0) == 0.0 + @testset "usinc_from_cos" begin + @test Manifolds.usinc_from_cos(-1) == 0 + @test Manifolds.usinc_from_cos(-1.0) == 0.0 + end + @testset "log_safe!" begin + n = 8 + Q = qr(randn(n, n)).Q + A1 = Matrix(Hermitian(Q * Diagonal(rand(n)) * Q')) + @test exp(Manifolds.log_safe!(similar(A1), A1)) ≈ A1 atol = 1e-6 + A1_fail = Matrix(Hermitian(Q * Diagonal([-1; rand(n - 1)]) * Q')) + @test_throws DomainError Manifolds.log_safe!(similar(A1_fail), A1_fail) + + T = triu!(randn(n, n)) + T[diagind(T)] .= rand.() + @test exp(Manifolds.log_safe!(similar(T), T)) ≈ T atol = 1e-6 + T_fail = copy(T) + T_fail[1] = -1 + @test_throws DomainError Manifolds.log_safe!(similar(T_fail), T_fail) + + A2 = Q * T * Q' + @test exp(Manifolds.log_safe!(similar(A2), A2)) ≈ A2 atol = 1e-6 + A2_fail = Q * T_fail * Q' + @test_throws DomainError Manifolds.log_safe!(similar(A2_fail), A2_fail) + + A3 = exp(randn(n, n)) + @test exp(Manifolds.log_safe!(similar(A3), A3)) ≈ A3 atol = 1e-6 + + A3_fail = Float64[1 2; 3 1] + @test_throws DomainError Manifolds.log_safe!(similar(A3_fail), A3_fail) + + A4 = randn(ComplexF64, n, n) + @test exp(Manifolds.log_safe!(similar(A4), A4)) ≈ A4 atol = 1e-6 + end + @testset "isnormal" begin + @test !Manifolds.isnormal([1.0 2.0; 3.0 4.0]) + @test !Manifolds.isnormal(complex.(reshape(1:4, 2, 2), reshape(5:8, 2, 2))) + + # diagonal + @test Manifolds.isnormal(diagm(randn(5))) + @test Manifolds.isnormal(diagm(randn(ComplexF64, 5))) + @test Manifolds.isnormal(Diagonal(randn(5))) + @test Manifolds.isnormal(Diagonal(randn(ComplexF64, 5))) + + # symmetric/hermitian + @test Manifolds.isnormal(Symmetric(randn(3, 3))) + @test Manifolds.isnormal(Hermitian(randn(3, 3))) + @test Manifolds.isnormal(Hermitian(randn(ComplexF64, 3, 3))) + x = Matrix(Symmetric(randn(3, 3))) + x[3, 1] += eps() + @test !Manifolds.isnormal(x) + @test Manifolds.isnormal(x; atol=sqrt(eps())) + + # skew-symmetric/skew-hermitian + skew(x) = x - x' + @test Manifolds.isnormal(skew(randn(3, 3))) + @test Manifolds.isnormal(skew(randn(ComplexF64, 3, 3))) + + # orthogonal/unitary + @test Manifolds.isnormal(Matrix(qr(randn(3, 3)).Q); atol=sqrt(eps())) + @test Manifolds.isnormal( + Matrix(qr(randn(ComplexF64, 3, 3)).Q); + atol=sqrt(eps()), + ) + end + @testset "realify/unrealify!" begin + # round trip real + x = randn(3, 3) + @test Manifolds.realify(x, ℝ) === x + @test Manifolds.unrealify!(similar(x), x, ℝ) == x + + # round trip complex + x2 = randn(ComplexF64, 3, 3) + x2r = Manifolds.realify(x2, ℂ) + @test eltype(x2r) <: Real + @test size(x2r) == (6, 6) + x2c = Manifolds.unrealify!(similar(x2), x2r, ℂ) + @test x2c ≈ x2 + + # matrix multiplication is preserved + x3 = randn(ComplexF64, 3, 3) + x3r = Manifolds.realify(x3, ℂ) + @test x2 * x3 ≈ Manifolds.unrealify!(similar(x2), x2r * x3r, ℂ) + end end include_test("groups/group_utils.jl") @@ -65,12 +146,15 @@ include("utils.jl") include_test("metric.jl") include_test("statistics.jl") + include_test("approx_inverse_retraction.jl") # Lie groups and actions include_test("groups/groups_general.jl") include_test("groups/array_manifold.jl") include_test("groups/circle_group.jl") include_test("groups/translation_group.jl") + include_test("groups/general_linear.jl") + include_test("groups/special_linear.jl") include_test("groups/special_orthogonal.jl") include_test("groups/product_group.jl") include_test("groups/semidirect_product_group.jl")