SymplecticMatrices is a Julia package for performing optimized linear algebra computations and special decompositions of symplectic matrices. A symplectic matrix S is a 2n x 2n matrix
that satisfies the condition SJSᵀ = J, where J is a 2n x 2n invertible, skew-symmetric matrix called the symplectic form.
The full API can be found at the documentation page:
To install SymplecticMatrices.jl, start Julia and run the following command:
] add SymplecticMatricesThe Symplectic matrix type is simply a wrapper around a matrix and its corresponding symplectic form.
To give a taste, compute a random symplectic matrix S by specifying the symplectic form type:
julia> using SymplecticMatrices
julia> J = BlockForm(2)
BlockForm(2)
julia> symplecticform(J)
4×4 Matrix{Float64}:
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
-1.0 0.0 0.0 0.0
0.0 -1.0 0.0 0.0
julia> S = randsymplectic(Symplectic, J)
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
-0.116535 -0.181523 1.15808 0.00240053
0.137176 -0.498848 0.55193 -1.22289
-0.853172 -0.578052 0.452435 -0.348686
0.00544036 0.618131 0.493416 -0.35965
julia> issymplectic(S, atol=1e-10)
trueSimilar methods exist for PairForm(n), which corresponds to the symplectic form equal to the direct
sum of [0 1; -1 0]. The matrix type Symplectic is nice for keeping track of the symplectic basis
and performing optimized computations, however, standard Julia matrix types are supported in this package:
julia> s = randsymplectic(J)
4×4 Matrix{Float64}:
-7.69543 -41.5856 13.5169 92.976
-2.04414 -8.3759 5.14638 18.4387
5.89402 33.3788 -9.05499 -74.8917
-6.10599 -30.7716 9.99773 68.9329
julia> issymplectic(J, s, atol=1e-10)
true
julia> issymplectic(PairForm(2), s, atol=1e-10)
falseIn the last line of code, the symplectic check failed because we defined s with respect to BlockForm rather than PairForm.
To compute the symplectic polar decomposition of S, which produces a product of an orthosymplectic matrix O and positive-definite symmetric symplectic matrix P, call polar:
julia> F = polar(S)
Polar{Float64, Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}, Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}}
O factor:
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
0.123739 -0.181484 0.923535 0.314381
0.316785 -0.36622 0.177255 -0.856803
-0.923535 -0.314381 0.123739 -0.181484
-0.177255 0.856803 0.316785 -0.36622
P factor:
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
0.816005 0.243795 -0.187156 -0.00132403
0.243795 0.926977 -0.131778 0.248883
-0.187156 -0.131778 1.37965 -0.371624
-0.00132403 0.248883 -0.371624 1.24352
julia> isapprox(F.O * F.P, S, atol = 1e-10)
trueTo compute the Williamson decomposition of a positive definite matrix V, which finds a congruence relation between a symplectic matrix S and a diagonal matrix containing the symplectic eigenvalues spectrum, call williamson:
julia> X = rand(4, 4); V = X' * X
4×4 Matrix{Float64}:
1.39837 1.01419 0.819762 0.644317
1.01419 1.0662 0.452062 0.525889
0.819762 0.452062 1.1312 1.02089
0.644317 0.525889 1.02089 1.09805
julia> F = williamson(Symplectic, J, V)
Williamson{Float64, Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}, Vector{Float64}}
S factor:
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
-0.67512 0.66685 0.586 -0.437417
0.733677 0.731394 -0.305086 -0.143899
0.686747 -0.647894 -1.33318 1.17835
-0.24263 -0.213864 0.785636 0.722449
symplectic spectrum:
2-element Vector{Float64}:
0.05191288411625232
1.8137043652121851
julia> isapprox(F.S * V * F.S', Diagonal(repeat(F.spectrum, 2)), atol = 1e-10)
trueTo compute the Bloch-Messiah/Euler decomposition of a symplectic matrix S, which finds the product of an orthosymplectic matrix followed by a diagonal matrix followed by another orthosymplectic matrix, call blochmessiah:
julia> F = blochmessiah(S)
BlochMessiah{Float64, Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}, Vector{Float64}}
O factor:
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
-0.636034 -0.492216 0.331745 0.493082
0.696422 -0.690058 0.0200433 0.195995
-0.331745 -0.493081 -0.636034 -0.492216
-0.0200437 -0.195995 0.696422 -0.690058
values:
2-element Vector{Float64}:
1.0016241730101116
1.8152715263879216
Q factor:
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
0.451846 -0.0524922 -0.511355 -0.729107
0.210612 0.329129 -0.699996 0.597764
0.511355 0.729107 0.451845 -0.0524917
0.699996 -0.597764 0.210612 0.32913
julia> D = Diagonal(vcat(F.values, F.values .^ (-1))) # sort singular values in direct sum form
4×4 Diagonal{Float64, Vector{Float64}}:
1.00162 ⋅ ⋅ ⋅
⋅ 1.81527 ⋅ ⋅
⋅ ⋅ 0.998378 ⋅
⋅ ⋅ ⋅ 0.550882
julia> F.O * D * F.Q # == S
4×4 Symplectic{BlockForm{Int64}, Float64, Matrix{Float64}}:
-0.116535 -0.181523 1.15808 0.0024005
0.137176 -0.498848 0.55193 -1.22289
-0.853173 -0.578053 0.452434 -0.348686
0.00544102 0.618132 0.493417 -0.35965Matrix Inverse
julia> S = randsymplectic(Symplectic, BlockForm(100));
julia> @btime inv($S.data);
650.084 μs (9 allocations: 414.30 KiB)
julia> @btime inv($S);
279.875 μs (12 allocations: 1.22 MiB)Symplectic Decompositions
julia> @btime polar(S) setup=(S=randsymplectic(Symplectic, BlockForm(100)));
6.826 ms (29 allocations: 3.08 MiB)
julia> @btime williamson(J, V) setup=(X=rand(200,200); V=X'*X; J=BlockForm(100));
22.587 ms (129 allocations: 6.77 MiB)
julia> @btime blochmessiah(S) setup=(S=randsymplectic(Symplectic, BlockForm(100)));
11.867 ms (64 allocations: 4.68 MiB)
Version Information
julia> versioninfo()
Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: macOS (arm64-apple-darwin22.4.0)
CPU: 8 × Apple M1 Pro
WORD_SIZE: 64
LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 6 virtual cores)