diff --git a/src/Manifolds.jl b/src/Manifolds.jl index 89c466e728..7727d593e7 100644 --- a/src/Manifolds.jl +++ b/src/Manifolds.jl @@ -4,6 +4,7 @@ module Manifolds import Base: + ^, angle, copyto!, convert, @@ -445,7 +446,9 @@ include("manifolds/SymmetricPositiveDefiniteLogCholesky.jl") include("manifolds/SymmetricPositiveDefiniteLogEuclidean.jl") include("manifolds/SymmetricPositiveSemidefiniteFixedRank.jl") include("manifolds/Symplectic.jl") +include("manifolds/Hamiltonian.jl") # Hamiltonian requires symplectic include("manifolds/SymplecticStiefel.jl") +include("manifolds/SymplecticGrassmann.jl") # Requires SymplecticStiefel include("manifolds/Tucker.jl") # include("manifolds/ProbabilitySimplex.jl") @@ -659,6 +662,7 @@ export Euclidean, SPDFixedDeterminant, SymmetricPositiveSemidefiniteFixedRank, Symplectic, + SymplecticGrassmann, SymplecticStiefel, SymplecticMatrix, Torus, @@ -682,7 +686,7 @@ export HyperboloidTVector, ProjectorTVector, StiefelTVector export AbstractNumbers, ℝ, ℂ, ℍ - +export Hamiltonian # decorator manifolds export AbstractDecoratorManifold export IsIsometricEmbeddedManifold, IsEmbeddedManifold, IsEmbeddedSubmanifold @@ -775,6 +779,7 @@ export CachedBasis, export ComponentManifoldError, CompositeManifoldError # Functions on Manifolds export ×, + ^, action_side, allocate, allocate_result, @@ -900,6 +905,8 @@ export ×, skewness, std, sym_rem, + symplectic_inverse, + symplectic_inverse!, symplectic_inverse_times, symplectic_inverse_times!, submanifold, diff --git a/src/manifolds/Hamiltonian.jl b/src/manifolds/Hamiltonian.jl new file mode 100644 index 0000000000..af5b1fd6f3 --- /dev/null +++ b/src/manifolds/Hamiltonian.jl @@ -0,0 +1,42 @@ +# +# This file requires Symplectic to be defined, since it need the symplectic inverse A^+ +# +@doc raw""" + Hamiltonian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T} + +A type to store an hamiltonien matrix, i.e. A square matrix matrix for which ``A^+ = -A`` where + +```math +A^+ = J_{2n}A^{\mathrm{T}}J_{2n}, \qquad J_{2n} \begin{pmatrix} 0 & I_n\\-I_n & 0 \end{pmatrix}, +``` + +and ``I_n`` denotes the ``n × n`` +""" + +struct Hamiltonian{T,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T} + data::S + function Hamiltonian(A::S) where {T,S<:AbstractMatrix{<:T}} + n = div(size(A, 1), 2) + @assert size(A, 1) == 2 * n "The first dimension of A ($(size(A,1))) is not even" + @assert size(A, 2) == 2 * n "The matrix A is of size ($(size(A))), which is not square." + return new{T,S}(A) + end +end + +# Conversion +function Matrix(A::Hamiltonian) + return copy(A.data) +end + +@doc raw""" + ^(A::Hamilonian, ::typeof(+)) + +Compute the [`symplectic_inverse`](@ref) of a Hamiltonian (A) +""" +function ^(A::Hamiltonian, ::typeof(+)) + return Hamiltonian(symplectic_inverse(A.data)) +end + +function show(io::IO, A::Hamiltonian) + return print(io, "Hamiltonian($(A.data))") +end diff --git a/src/manifolds/Symplectic.jl b/src/manifolds/Symplectic.jl index 67b3492b3f..3a40b39e9d 100644 --- a/src/manifolds/Symplectic.jl +++ b/src/manifolds/Symplectic.jl @@ -433,60 +433,78 @@ function inner(M::Symplectic{<:Any,ℝ}, p, X, Y) end @doc raw""" - inv(::Symplectic, A) - inv!(::Symplectic, A) + symplectic_inverse(A) + +Given a matrix +``math + A ∈ ℝ^{2n × 2k},\quad + A = + \begin{bmatrix} + A_{1,1} & A_{1,2} \\ + A_{2,1} & A_{2, 2} + \end{bmatrix} +``` -Compute the symplectic inverse ``A^+`` of matrix ``A ∈ ℝ^{2n × 2n}``. Given a matrix -````math -A ∈ ℝ^{2n × 2n},\quad -A = -\begin{bmatrix} -A_{1,1} & A_{1,2} \\ -A_{2,1} & A_{2, 2} -\end{bmatrix} -```` the symplectic inverse is defined as: -````math -A^{+} := Q_{2n}^{\mathrm{T}} A^{\mathrm{T}} Q_{2n}, -```` + +```math +A^{+} := Q_{2k}^{\mathrm{T}} A^{\mathrm{T}} Q_{2n}, +``` + where -````math -Q_{2n} = -\begin{bmatrix} -0_n & I_n \\ - -I_n & 0_n -\end{bmatrix}. -```` + +```math + Q_{2n} = \begin{bmatrix} 0_n & I_n \\ -I_n & 0_n \end{bmatrix}. +``` + The symplectic inverse of A can be expressed explicitly as: -````math + +```math A^{+} = -\begin{bmatrix} - A_{2, 2}^{\mathrm{T}} & -A_{1, 2}^{\mathrm{T}} \\[1.2mm] - -A_{2, 1}^{\mathrm{T}} & A_{1, 1}^{\mathrm{T}} -\end{bmatrix}. -```` + \begin{bmatrix} + A_{2, 2}^{\mathrm{T}} & -A_{1, 2}^{\mathrm{T}} \\[1.2mm] + -A_{2, 1}^{\mathrm{T}} & A_{1, 1}^{\mathrm{T}} + \end{bmatrix}. +``` + """ -function Base.inv(M::Symplectic{<:Any,ℝ}, A) - n = get_parameter(M.size)[1] - Ai = similar(A) - checkbounds(A, 1:(2n), 1:(2n)) - @inbounds for i in 1:n, j in 1:n - Ai[i, j] = A[j + n, i + n] +function symplectic_inverse(A::AbstractMatrix) + N, K = size(A) + @assert iseven(N) "The first matrix dimension of A ($N) has to be even" + @assert iseven(K) "The second matrix dimension of A ($K) has to be even" + n = div(N, 2) + k = div(K, 2) + Ai = similar(A') + checkbounds(A, 1:(2n), 1:(2k)) + @inbounds for i in 1:k, j in 1:n + Ai[i, j] = A[j + n, i + k] end - @inbounds for i in 1:n, j in 1:n - Ai[i + n, j] = -A[j + n, i] + @inbounds for i in 1:k, j in 1:n + Ai[i + k, j] = -A[j + n, i] end - @inbounds for i in 1:n, j in 1:n - Ai[i, j + n] = -A[j, i + n] + @inbounds for i in 1:k, j in 1:n + Ai[i, j + n] = -A[j, i + k] end - @inbounds for i in 1:n, j in 1:n - Ai[i + n, j + n] = A[j, i] + @inbounds for i in 1:k, j in 1:n + Ai[i + k, j + n] = A[j, i] end return Ai end -function inv!(M::Symplectic{<:Any,ℝ}, A) - n = get_parameter(M.size)[1] +@doc raw""" + inv(::Symplectic, A) + inv!(::Symplectic, A) + +Compute the symplectic inverse ``A^+`` of matrix ``A ∈ ℝ^{2n × 2n}``. See [`symplectic_inverse`](@ref) +for details. + +""" +function Base.inv(M::Symplectic{<:Any,ℝ}, A) + return symplectic_inverse(A) +end + +function symplectic_inverse!(A) + n = div(size(A, 1), 1) checkbounds(A, 1:(2n), 1:(2n)) @inbounds for i in 1:n, j in 1:n A[i, j], A[j + n, i + n] = A[j + n, i + n], A[i, j] @@ -508,6 +526,15 @@ function inv!(M::Symplectic{<:Any,ℝ}, A) return A end +@doc raw""" + inv!(M::Symplectic, A) + +Compute the symplectic inverse of a suqare matrix A inplace of A +""" +function inv!(M::Symplectic{<:Any,ℝ}, A) + return symplectic_inverse!(A) +end + @doc raw""" inverse_retract(M::Symplectic, p, q, ::CayleyInverseRetraction)