diff --git a/src/Manifolds.jl b/src/Manifolds.jl index b33591b928..eb3937e1dc 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -660,7 +660,7 @@ export Euclidean, SymmetricPositiveSemidefiniteFixedRank, SymplecticMatrices, SymplecticStiefel, - SymplecticMatrix, + SymplecticElement, Torus, Tucker, UnitaryMatrices diff --git a/src/deprecated.jl b/src/deprecated.jl index ad091f408e..eea4fe49ac 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -6,4 +6,4 @@ @deprecate ExtrinsicEstimation() ExtrinsicEstimation(EfficientEstimator()) Base.@deprecate_binding Symplectic SymplecticMatrices -#Base.@deprecate_binding SynplecticMatrix SymplecticElement +Base.@deprecate_binding SymplecticMatrix SymplecticElement diff --git a/src/manifolds/Symplectic.jl b/src/manifolds/Symplectic.jl index b9fdf4e1ec..469390eeea 100644 --- a/src/manifolds/Symplectic.jl +++ b/src/manifolds/Symplectic.jl @@ -82,7 +82,7 @@ as an inner product over the embedding space ``ℝ^{2n \times 2n}``, i.e. struct ExtendedSymplecticMetric <: AbstractMetric end @doc raw""" - SymplecticMatrix{T} + SymplecticElement{T} A lightweight structure to represent the action of the matrix representation of the canonical symplectic form, @@ -100,19 +100,19 @@ such that the canonical symplectic form is represented by The entire matrix is however not instantiated in memory, instead a scalar ``λ`` of type `T` is stored, which is used to keep track of scaling and transpose operations -applied to each `SymplecticMatrix`. -For example, given `Q = SymplecticMatrix(1.0)` represented as `1.0*[0 I; -I 0]`, -the adjoint `Q'` returns `SymplecticMatrix(-1.0) = (-1.0)*[0 I; -I 0]`. +applied to each `SymplecticElement`. +For example, given `Q = SymplecticElement(1.0)` represented as `1.0*[0 I; -I 0]`, +the adjoint `Q'` returns `SymplecticElement(-1.0) = (-1.0)*[0 I; -I 0]`. """ -struct SymplecticMatrix{T} +struct SymplecticElement{T} λ::T end -SymplecticMatrix() = SymplecticMatrix(1) -SymplecticMatrix(λ::T) where {T<:Number} = SymplecticMatrix{T}(λ) +SymplecticElement() = SymplecticElement(1) +SymplecticElement(λ::T) where {T<:Number} = SymplecticElement{T}(λ) -function SymplecticMatrix(arrays::Vararg{AbstractArray}) +function SymplecticElement(arrays::Vararg{AbstractArray}) TS = Base.promote_type(map(eltype, arrays)...) - return SymplecticMatrix(one(TS)) + return SymplecticElement(one(TS)) end @doc raw""" @@ -145,7 +145,7 @@ does have the correct range ``T_p\operatorname{Sp}(2n, ℝ)``. change_representer(::SymplecticMatrices, ::EuclideanMetric, p, X) function change_representer!(::SymplecticMatrices, Y, ::EuclideanMetric, p, X) - Q = SymplecticMatrix(p, X) + Q = SymplecticElement(p, X) pT_X = p' * X Y .= (1 / 2) .* p * (pT_X .+ Q * pT_X' * Q) return Y @@ -257,7 +257,7 @@ function check_vector( atol::Real=sqrt(prod(representation_size(M))) * eps(real(float(number_eltype(T)))), kwargs..., ) where {T} - Q = SymplecticMatrix(p, X) + Q = SymplecticElement(p, X) tangent_requirement_norm = norm(X' * Q * p + p' * Q * X, 2) if !isapprox(tangent_requirement_norm, 0; atol=atol, kwargs...) return DomainError( @@ -597,7 +597,7 @@ the restriction of ``X`` onto the tangent space ``T_p\operatorname{SpSt}(2n, 2k) project(::SymplecticMatrices, p, A) function project!(::SymplecticMatrices, Y, p, A) - Q = SymplecticMatrix(Y, p, A) + Q = SymplecticElement(Y, p, A) Q_p = Q * p function h(X) @@ -628,7 +628,7 @@ where ``\operatorname{sym}(A) = \frac{1}{2}(A + A^{\mathrm{T}})``. This function is not exported. """ function project!(::MetricManifold{<:Any,<:Euclidean,ExtendedSymplecticMetric}, Y, p, X) - Q = SymplecticMatrix(p, X) + Q = SymplecticElement(p, X) pT_QT_X = p' * Q' * X symmetrized_pT_QT_X = (1 / 2) .* (pT_QT_X + pT_QT_X') @@ -665,7 +665,7 @@ function project_normal!( p, X, ) where {𝔽} - Q = SymplecticMatrix(p, X) + Q = SymplecticElement(p, X) pT_QT_X = p' * Q' * X skew_pT_QT_X = (1 / 2) .* (pT_QT_X .- pT_QT_X') @@ -694,7 +694,7 @@ second tangent vector space parametrization of [`SymplecticMatrices`](@ref). It first generates a random symmetric matrix ``S`` by `S = randn(2n, 2n)` and then symmetrizes it as `S = S + S'`. Then ``S`` is normalized to have Frobenius norm of `hamiltonian_norm` -and `X = pQS` is returned, where `Q` is the [`SymplecticMatrix`](@ref). +and `X = pQS` is returned, where `Q` is the [`SymplecticElement`](@ref). """ function Random.rand( M::SymplecticMatrices; @@ -715,7 +715,7 @@ function random_vector(M::SymplecticMatrices, p::AbstractMatrix; symmetric_norm= S = randn(2n, 2n) S .= (S + S') S *= symmetric_norm / norm(S) - Q = SymplecticMatrix(p) + Q = SymplecticElement(p) lmul!(Q, S) return p * S end @@ -791,7 +791,7 @@ Directly compute the symplectic inverse of ``p \in \operatorname{Sp}(2n)``, multiplied with ``q \in \operatorname{Sp}(2n)``. That is, this function efficiently computes ``p^+q = (Q_{2n}p^{\mathrm{T}}Q_{2n})q \in ℝ^{2n \times 2n}``, -where ``Q_{2n}`` is the [`SymplecticMatrix`](@ref) +where ``Q_{2n}`` is the [`SymplecticElement`](@ref) of size ``2n \times 2n``. """ function symplectic_inverse_times(M::SymplecticMatrices, p, q) @@ -825,14 +825,14 @@ function symplectic_inverse_times!(M::SymplecticMatrices, A, p, q) return A end -ndims(Q::SymplecticMatrix) = 2 -copy(Q::SymplecticMatrix) = SymplecticMatrix(copy(Q.λ)) -Base.eltype(::SymplecticMatrix{T}) where {T} = T -function Base.convert(::Type{SymplecticMatrix{T}}, Q::SymplecticMatrix) where {T} - return SymplecticMatrix(convert(T, Q.λ)) +ndims(Q::SymplecticElement) = 2 +copy(Q::SymplecticElement) = SymplecticElement(copy(Q.λ)) +Base.eltype(::SymplecticElement{T}) where {T} = T +function Base.convert(::Type{SymplecticElement{T}}, Q::SymplecticElement) where {T} + return SymplecticElement(convert(T, Q.λ)) end -function Base.show(io::IO, Q::SymplecticMatrix) +function Base.show(io::IO, Q::SymplecticElement) s = "$(Q.λ)" if occursin(r"\w+\s*[\+\-]\s*\w+", s) s = "($s)" @@ -840,31 +840,31 @@ function Base.show(io::IO, Q::SymplecticMatrix) return print(io, typeof(Q), "(): $(s)*[0 I; -I 0]") end -(Base.:-)(Q::SymplecticMatrix) = SymplecticMatrix(-Q.λ) +(Base.:-)(Q::SymplecticElement) = SymplecticElement(-Q.λ) -function (Base.:^)(Q::SymplecticMatrix, n::Integer) +function (Base.:^)(Q::SymplecticElement, n::Integer) return ifelse( n % 2 == 0, UniformScaling((-1)^(div(n, 2)) * (Q.λ)^n), - SymplecticMatrix((-1)^(div(n - 1, 2)) * (Q.λ)^n), + SymplecticElement((-1)^(div(n - 1, 2)) * (Q.λ)^n), ) end -(Base.:*)(x::Number, Q::SymplecticMatrix) = SymplecticMatrix(x * Q.λ) -(Base.:*)(Q::SymplecticMatrix, x::Number) = SymplecticMatrix(x * Q.λ) -function (Base.:*)(Q1::SymplecticMatrix, Q2::SymplecticMatrix) +(Base.:*)(x::Number, Q::SymplecticElement) = SymplecticElement(x * Q.λ) +(Base.:*)(Q::SymplecticElement, x::Number) = SymplecticElement(x * Q.λ) +function (Base.:*)(Q1::SymplecticElement, Q2::SymplecticElement) return LinearAlgebra.UniformScaling(-Q1.λ * Q2.λ) end -Base.transpose(Q::SymplecticMatrix) = -Q -Base.adjoint(Q::SymplecticMatrix) = SymplecticMatrix(-conj(Q.λ)) -Base.inv(Q::SymplecticMatrix) = SymplecticMatrix(-(1 / Q.λ)) +Base.transpose(Q::SymplecticElement) = -Q +Base.adjoint(Q::SymplecticElement) = SymplecticElement(-conj(Q.λ)) +Base.inv(Q::SymplecticElement) = SymplecticElement(-(1 / Q.λ)) -(Base.:+)(Q1::SymplecticMatrix, Q2::SymplecticMatrix) = SymplecticMatrix(Q1.λ + Q2.λ) -(Base.:-)(Q1::SymplecticMatrix, Q2::SymplecticMatrix) = SymplecticMatrix(Q1.λ - Q2.λ) +(Base.:+)(Q1::SymplecticElement, Q2::SymplecticElement) = SymplecticElement(Q1.λ + Q2.λ) +(Base.:-)(Q1::SymplecticElement, Q2::SymplecticElement) = SymplecticElement(Q1.λ - Q2.λ) -(Base.:+)(Q::SymplecticMatrix, p::AbstractMatrix) = p + Q -function (Base.:+)(p::AbstractMatrix, Q::SymplecticMatrix) +(Base.:+)(Q::SymplecticElement, p::AbstractMatrix) = p + Q +function (Base.:+)(p::AbstractMatrix, Q::SymplecticElement) # When we are adding, the Matrices must match in size: two_n, two_k = size(p) if (two_n % 2 != 0) || (two_n != two_k) @@ -887,10 +887,10 @@ function (Base.:+)(p::AbstractMatrix, Q::SymplecticMatrix) end # Binary minus: -(Base.:-)(Q::SymplecticMatrix, p::AbstractMatrix) = Q + (-p) -(Base.:-)(p::AbstractMatrix, Q::SymplecticMatrix) = p + (-Q) +(Base.:-)(Q::SymplecticElement, p::AbstractMatrix) = Q + (-p) +(Base.:-)(p::AbstractMatrix, Q::SymplecticElement) = p + (-Q) -function (Base.:*)(Q::SymplecticMatrix, p::AbstractVecOrMat) +function (Base.:*)(Q::SymplecticElement, p::AbstractVecOrMat) two_n = size(p)[1] if two_n % 2 != 0 throw(ArgumentError("'p' must have even row dimension, was: $(two_n) != 2n.")) @@ -908,7 +908,7 @@ function (Base.:*)(Q::SymplecticMatrix, p::AbstractVecOrMat) return Qp end -function (Base.:*)(p::AbstractMatrix, Q::SymplecticMatrix) +function (Base.:*)(p::AbstractMatrix, Q::SymplecticElement) two_k = size(p)[2] if two_k % 2 != 0 throw(ArgumentError("'p' must have even column dimension, was: $(two_k) != 2k.")) @@ -925,7 +925,7 @@ function (Base.:*)(p::AbstractMatrix, Q::SymplecticMatrix) return pQ end -function LinearAlgebra.lmul!(Q::SymplecticMatrix, p::AbstractVecOrMat) +function LinearAlgebra.lmul!(Q::SymplecticElement, p::AbstractVecOrMat) # Perform left multiplication by a symplectic matrix, # overwriting the matrix p in place: two_n = size(p)[1] @@ -946,7 +946,7 @@ function LinearAlgebra.lmul!(Q::SymplecticMatrix, p::AbstractVecOrMat) return p end -function LinearAlgebra.rmul!(p::AbstractMatrix, Q::SymplecticMatrix) +function LinearAlgebra.rmul!(p::AbstractMatrix, Q::SymplecticElement) # Perform right multiplication by a symplectic matrix, # overwriting the matrix p in place: two_k = size(p)[2] @@ -968,7 +968,7 @@ function LinearAlgebra.rmul!(p::AbstractMatrix, Q::SymplecticMatrix) return p end -function LinearAlgebra.mul!(A::AbstractVecOrMat, Q::SymplecticMatrix, p::AbstractVecOrMat) +function LinearAlgebra.mul!(A::AbstractVecOrMat, Q::SymplecticElement, p::AbstractVecOrMat) size_p = size(p) two_n = size_p[1] if two_n % 2 != 0 @@ -986,7 +986,7 @@ function LinearAlgebra.mul!(A::AbstractVecOrMat, Q::SymplecticMatrix, p::Abstrac return A end -function LinearAlgebra.mul!(A::AbstractVecOrMat, p::AbstractMatrix, Q::SymplecticMatrix) +function LinearAlgebra.mul!(A::AbstractVecOrMat, p::AbstractMatrix, Q::SymplecticElement) two_n, two_k = size(p) if two_k % 2 != 0 throw(ArgumentError("'p' must have even col dimension, was: $(two_k) != 2k.")) diff --git a/src/manifolds/SymplecticStiefel.jl b/src/manifolds/SymplecticStiefel.jl index d602737479..739f564208 100644 --- a/src/manifolds/SymplecticStiefel.jl +++ b/src/manifolds/SymplecticStiefel.jl @@ -191,7 +191,7 @@ The tangent vector ``X`` can be written in the form + Q_{2n}p(p^{\mathrm{T}}p)^{-1}X^{\mathrm{T}}(I_{2n} - Q_{2n}^{\mathrm{T}}p(p^{\mathrm{T}}p)^{-1}p^{\mathrm{T}}Q_{2n})Q_{2n} \in ℝ^{2n \times 2n}, ```` -where ``Q_{2n}`` is the [`SymplecticMatrix`](@ref). Using this expression for ``X``, +where ``Q_{2n}`` is the [`SymplecticElement`](@ref). Using this expression for ``X``, the exponential mapping can be computed as ````math \operatorname{exp}_p(X) = \operatorname{Exp}([\bar{\Omega} - \bar{\Omega}^{\mathrm{T}}]) @@ -252,7 +252,7 @@ exp(::SymplecticStiefel, p, X) function exp!(M::SymplecticStiefel, q, p, X) n, k = get_parameter(M.size) - Q = SymplecticMatrix(p, X) + Q = SymplecticElement(p, X) pT_p = lu(p' * p) # ∈ ℝ^{2k × 2k} C = pT_p \ X' # ∈ ℝ^{2k × 2n} @@ -332,7 +332,7 @@ g^{\operatorname{SpSt}}_p(X, Y) ```` """ function inner(::SymplecticStiefel, p, X, Y) - Q = SymplecticMatrix(p, X, Y) + Q = SymplecticElement(p, X, Y) # Procompute lu(p'p) since we solve a^{-1}* 3 times a = lu(p' * p) # note that p'p is symmetric, thus so is its inverse c=a^{-1} b = Q' * p @@ -490,7 +490,7 @@ the restriction of ``X`` onto the tangent space ``T_p\operatorname{SpSt}(2n, 2k) project(::SymplecticStiefel, p, A) function project!(::SymplecticStiefel, Y, p, A) - Q = SymplecticMatrix(Y, p, A) + Q = SymplecticElement(Y, p, A) Q_p = Q * p function h(X) @@ -615,12 +615,12 @@ The manifold gradient `X` is computed from `Y` as ```math X = Yp^{\mathrm{T}}p + Q_{2n}pY^{\mathrm{T}}Q_{2n}p, ``` -where ``Q_{2n}`` is the [`SymplecticMatrix`](@ref). +where ``Q_{2n}`` is the [`SymplecticElement`](@ref). """ function riemannian_gradient(::SymplecticStiefel, p, Y) - Q_p = SymplecticMatrix(p, Y) * p + Q_p = SymplecticElement(p, Y) * p return Y * (p' * p) .+ Q_p * (Y' * Q_p) end @@ -631,7 +631,7 @@ function riemannian_gradient!( Y; embedding_metric::EuclideanMetric=EuclideanMetric(), ) - Q_p = SymplecticMatrix(p, Y) * p + Q_p = SymplecticElement(p, Y) * p X .= Y * (p' * p) .+ Q_p * (Y' * Q_p) return X end @@ -652,7 +652,7 @@ Directly compute the symplectic inverse of ``p \in \operatorname{SpSt}(2n, 2k)`` multiplied with ``q \in \operatorname{SpSt}(2n, 2k)``. That is, this function efficiently computes ``p^+q = (Q_{2k}p^{\mathrm{T}}Q_{2n})q \in ℝ^{2k \times 2k}``, -where ``Q_{2n}, Q_{2k}`` are the [`SymplecticMatrix`](@ref) +where ``Q_{2n}, Q_{2k}`` are the [`SymplecticElement`](@ref) of sizes ``2n \times 2n`` and ``2k \times 2k`` respectively. This function performs this common operation without allocating more than diff --git a/test/manifolds/symplectic.jl b/test/manifolds/symplectic.jl index 204a4f1b42..47fedb0652 100644 --- a/test/manifolds/symplectic.jl +++ b/test/manifolds/symplectic.jl @@ -258,7 +258,7 @@ using ManifoldDiff @testset "Gradient Computations" begin test_f(p) = tr(p) - Q_grad = SymplecticMatrix(points[1]) + Q_grad = SymplecticElement(points[1]) analytical_grad_f(p) = (1 / 2) * (p * Q_grad * p * Q_grad + p * p') p_grad = convert(Array{Float64}, points[1]) @@ -291,9 +291,9 @@ using ManifoldDiff end end - @testset "SymplecticMatrix" begin + @testset "SymplecticElement" begin # TODO: Test for different type matrices. - @test SymplecticMatrix() == SymplecticMatrix(1) + @test SymplecticElement() == SymplecticElement(1) Sp_4 = SymplecticMatrices(4) pQ_1 = [ 0 0 -2 3 @@ -318,30 +318,30 @@ using ManifoldDiff -1 0 0 -2 4 -8 ] - Q = SymplecticMatrix(pQ_1, pQ_2) - Q2 = SymplecticMatrix(1) + Q = SymplecticElement(pQ_1, pQ_2) + Q2 = SymplecticElement(1) @testset "Type Basics" begin @test Q == Q2 @test ndims(Q) == 2 @test copy(Q) == Q - @test eltype(SymplecticMatrix(1 // 1)) == Rational{Int64} - @test convert(SymplecticMatrix{Float64}, Q) == SymplecticMatrix(1.0) - @test "$Q" == "SymplecticMatrix{Int64}(): 1*[0 I; -I 0]" + @test eltype(SymplecticElement(1 // 1)) == Rational{Int64} + @test convert(SymplecticElement{Float64}, Q) == SymplecticElement(1.0) + @test "$Q" == "SymplecticElement{Int64}(): 1*[0 I; -I 0]" @test ( - "$(SymplecticMatrix(1 + 2im))" == - "SymplecticMatrix{Complex{Int64}}(): (1 + 2im)*[0 I; -I 0]" + "$(SymplecticElement(1 + 2im))" == + "SymplecticElement{Complex{Int64}}(): (1 + 2im)*[0 I; -I 0]" ) end @testset "Matrix Operations" begin - @test -Q == SymplecticMatrix(-1) - @test (2 * Q) * (5 // 6) == SymplecticMatrix(5 // 3) + @test -Q == SymplecticElement(-1) + @test (2 * Q) * (5 // 6) == SymplecticElement(5 // 3) @testset "Powers" begin @test inv(Q) * Q == I @test ( - inv(SymplecticMatrix(-4.0 + 8im)) * SymplecticMatrix(-4.0 + 8im) == + inv(SymplecticElement(-4.0 + 8im)) * SymplecticElement(-4.0 + 8im) == UniformScaling(1.0 + 0.0im) ) @test Q * Q == -I @@ -351,7 +351,7 @@ using ManifoldDiff end @testset "Addition (subtraction)" begin @test Q + Q == 2 * Q - @test Q - SymplecticMatrix(1.0) == SymplecticMatrix(0.0) + @test Q - SymplecticElement(1.0) == SymplecticElement(0.0) @test Q + pQ_1 == [ 0 0 -1 3 0 0 1 0 @@ -375,13 +375,14 @@ using ManifoldDiff @test_throws ArgumentError Q + p_odd_row end @testset "Transpose-Adjoint" begin - @test Q' == SymplecticMatrix(-1) - @test transpose(SymplecticMatrix(10)) == SymplecticMatrix(-10) - @test transpose(SymplecticMatrix(1 - 2.0im)) == SymplecticMatrix(-1 + 2.0im) + @test Q' == SymplecticElement(-1) + @test transpose(SymplecticElement(10)) == SymplecticElement(-10) + @test transpose(SymplecticElement(1 - 2.0im)) == + SymplecticElement(-1 + 2.0im) @test adjoint(Q) == -Q - @test adjoint(SymplecticMatrix(1 - 2.0im)) == SymplecticMatrix(-1 - 2.0im) - @test adjoint(SymplecticMatrix(-1im)) == SymplecticMatrix(-1im) - @test adjoint(SymplecticMatrix(2.0)) == SymplecticMatrix(-2.0) + @test adjoint(SymplecticElement(1 - 2.0im)) == SymplecticElement(-1 - 2.0im) + @test adjoint(SymplecticElement(-1im)) == SymplecticElement(-1im) + @test adjoint(SymplecticElement(2.0)) == SymplecticElement(-2.0) end @testset "Inplace mul!" begin z1 = [1 + 2im; 1 - 2im] diff --git a/test/manifolds/symplecticstiefel.jl b/test/manifolds/symplecticstiefel.jl index f6d28a0d22..c079fd9dbf 100644 --- a/test/manifolds/symplecticstiefel.jl +++ b/test/manifolds/symplecticstiefel.jl @@ -4,7 +4,7 @@ using Manifolds: RiemannianProjectionBackend using ManifoldDiff function Ω(::SymplecticStiefel, p, X) - Q = SymplecticMatrix(X, p) + Q = SymplecticElement(X, p) pT_p = lu(p' * p) inv_pTp_pT = pT_p \ p' @@ -284,7 +284,7 @@ end end @testset "Gradient Computations" begin - Q_grad = SymplecticMatrix(points[1]) + Q_grad = SymplecticElement(points[1]) function test_f(p) k = size(p)[2] return tr(p[1:k, 1:k])