diff --git a/CHANGELOG.md b/CHANGELOG.md index fc4c007..aea186b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## v1.0.3 - dev - Add `prob` function for `Generaldyne` type. +- Add features for random generation of Gaussian states, unitaries, and channels. +- **(breaking)** Switched basis for matrices generated by `symplecticform` from block diagonal to canonical. ## v1.0.2 - 2024-11-03 diff --git a/Project.toml b/Project.toml index 4f05a61..9fb5e66 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Gabs" uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50" authors = ["Andrew Kille"] -version = "1.0.3-dev" +version = "1.1.0" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" diff --git a/src/Gabs.jl b/src/Gabs.jl index 8d3a5c4..1a5a45c 100644 --- a/src/Gabs.jl +++ b/src/Gabs.jl @@ -3,7 +3,7 @@ module Gabs using BlockArrays: BlockedArray, BlockArray, Block, mortar import LinearAlgebra -using LinearAlgebra: I, det, mul! +using LinearAlgebra: I, det, mul!, diagm, diag, qr import QuantumInterface: StateVector, AbstractOperator, apply!, tensor, ⊗ @@ -13,14 +13,18 @@ export # operations tensor, ⊗, apply!, ptrace, output, prob, # predefined Gaussian states - vacuumstate, thermalstate, coherentstate, squeezedstate, eprstate, randstate, + vacuumstate, thermalstate, coherentstate, squeezedstate, eprstate, # predefined Gaussian channels displace, squeeze, twosqueeze, phaseshift, beamsplitter, - attenuator, amplifier, randchannel, + attenuator, amplifier, + # random objects + randstate, randunitary, randchannel, randsymplectic, # wigner functions wigner, wignerchar, # symplectic form - symplecticform + symplecticform, + # metrics + purity include("errors.jl") @@ -40,4 +44,6 @@ include("measurements.jl") include("wigner.jl") +include("metrics.jl") + end diff --git a/src/metrics.jl b/src/metrics.jl new file mode 100644 index 0000000..e145493 --- /dev/null +++ b/src/metrics.jl @@ -0,0 +1,6 @@ +""" + purity(state::GaussianState) + +Calculate the purity of a Gaussian state, defined by `1/sqrt(det(V))`. +""" +purity(state::GaussianState) = 1/sqrt(det(state.covar)) \ No newline at end of file diff --git a/src/randoms.jl b/src/randoms.jl index d8f8a6a..91fc2d7 100644 --- a/src/randoms.jl +++ b/src/randoms.jl @@ -1,23 +1,54 @@ """ - randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] nmodes<:Int) + randstate([Tm=Vector{Float64}, Tc=Matrix{Float64},] nmodes<:Int; pure=false) Calculate a random Gaussian state. """ -function randstate(::Type{Tm}, ::Type{Tc}, nmodes::N) where {Tm,Tc,N<:Int} - mean = randn(2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - covar = randmat' * randmat +function randstate(::Type{Tm}, ::Type{Tc}, nmodes::N; pure = false) where {Tm,Tc,N<:Int} + mean, covar = _randstate_fields(nmodes, pure) return GaussianState(Tm(mean), Tc(covar), nmodes) end -randstate(::Type{T}, nmodes::N) where {T,N<:Int} = randstate(T,T,nmodes) -function randstate(nmodes::N) where {N<:Int} - mean = randn(2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - covar = randmat' * randmat +randstate(::Type{T}, nmodes::N; pure = false) where {T,N<:Int} = randstate(T,T,nmodes,pure=pure) +function randstate(nmodes::N; pure = false) where {N<:Int} + mean, covar = _randstate_fields(nmodes, pure) return GaussianState(mean, covar, nmodes) end +function _randstate_fields(nmodes::N, pure) where {N<:Int} + mean = randn(2*nmodes) + covar = zeros(2*nmodes, 2*nmodes) + symp = randsymplectic(nmodes) + # generate pure Gaussian state + if pure + mul!(covar, symp, symp') + return mean, covar + end + # create buffer for matrix multiplication + buf = zeros(2*nmodes, 2*nmodes) + # William decomposition for mixed Gaussian states + sympeigs = abs.(rand(nmodes)) .+ 1.0 + will = diagm(repeat(sympeigs, inner = 2)) + mul!(covar, symp, mul!(buf, will, symp')) + return mean, covar +end + +""" + randunitary([Td=Vector{Float64}, Ts=Matrix{Float64},] nmodes<:Int; passive=false) + +Calculate a random Gaussian unitary operator. +""" +function randunitary(::Type{Td}, ::Type{Ts}, nmodes::N; passive = false) where {Td,Ts,N} + disp, symp = _randunitary_fields(nmodes, passive) + return GaussianUnitary(Td(disp), Ts(symp)) +end +randunitary(::Type{T}, nmodes::N; passive = false) where {T,N<:Int} = randunitary(T,T,nmodes; passive = passive) +function randunitary(nmodes::N; passive = false) where {N<:Int} + disp, symp = _randunitary_fields(nmodes, passive) + return GaussianUnitary(disp, symp) +end +function _randunitary_fields(nmodes::N, passive) where {N<:Int} + disp = rand(2*nmodes) + symp = randsymplectic(nmodes, passive = passive) + return disp, symp +end """ randchannel([Td=Vector{Float64}, Tt=Matrix{Float64},] nmodes<:Int) @@ -25,19 +56,74 @@ end Calculate a random Gaussian channel. """ function randchannel(::Type{Td}, ::Type{Tt}, nmodes::N) where {Td,Tt,N<:Int} - disp = randn(2*nmodes) - transform = randn(2*nmodes, 2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - noise = randmat' * randmat + disp, transform, noise = _randchannel(nmodes) return GaussianChannel(Td(disp), Tt(transform), Tt(noise), nmodes) end randchannel(::Type{T}, nmodes::N) where {T,N<:Int} = randchannel(T,T,nmodes) function randchannel(nmodes::N) where {N<:Int} - disp = randn(2*nmodes) - transform = randn(2*nmodes, 2*nmodes) - randmat = randn(2*nmodes, 2*nmodes) - # positive-definite condition - noise = randmat' * randmat + disp, transform, noise = _randchannel(nmodes) return GaussianChannel(disp, transform, noise, nmodes) +end +function _randchannel(nmodes::N) where {N<:Int} + disp = randn(2*nmodes) + # generate symplectic matrix describing the evolution of the system with N modes + # and environment with 2N modes + symp = randsymplectic(3*nmodes) + transform, B = symp[1:2*nmodes, 1:2*nmodes], @view(symp[1:2*nmodes, 2*nmodes+1:6*nmodes]) + # generate noise matrix from evolution of environment + noise = zeros(2*nmodes, 2*nmodes) + mul!(noise, B, B') + return disp, transform, noise +end + +""" + randsymplectic([T=Matrix{Float64},] nmodes<:Int, passive=false) + +Calculate a random symplectic matrix of size `2*nmodes x 2*nmodes`. +""" +function randsymplectic(nmodes::N; passive = false) where {N<:Int} + # generate random orthogonal symplectic matrix + O = _rand_orthogonal_symplectic(nmodes) + if passive + return O + end + # instead generate random active symplectic matrix + O′ = _rand_orthogonal_symplectic(nmodes) + # direct sum of symplectic matrices for single-mode squeeze transformations + rs = rand(nmodes) + squeezes = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes) + val = rs[i] + squeezes[2*i-1, 2*i-1] = val + squeezes[2*i, 2*i] = 1/val + end + return O * squeezes * O′ +end +function randsymplectic(::Type{T}, nmodes::N; passive = false) where {T, N<:Int} + symp = randsymplectic(nmodes, passive = passive) + return T(symp) +end + +# Generates random orthogonal symplectic matrix by blocking real +# and imaginary parts of a random unitary matrix +function _rand_orthogonal_symplectic(nmodes::N) where {N<:Int} + U = _rand_unitary(nmodes) + O = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes), j in Base.OneTo(nmodes) + val = U[i,j] + O[2*i-1,2*j-1] = real(val) + O[2*i, 2*j-1] = -imag(val) + O[2*i-1, 2*j] = imag(val) + O[2*i, 2*j] = real(val) + end + return O +end +# Generates unitary matrix randomly distributed over Haar measure; +# see https://arxiv.org/abs/math-ph/0609050 for algorithm. +# This approach is faster and creates less allocations than rand(Haar(2), size) from RandomMatrices.jl +function _rand_unitary(size::N) where {N<:Int} + M = rand(ComplexF64, size, size) ./ sqrt(2.0) + q, r = qr(M) + d = diagm([r[i, i] / abs(r[i, i]) for i in Base.OneTo(size)]) + return q * d end \ No newline at end of file diff --git a/src/wigner.jl b/src/wigner.jl index b6e2987..5b95fb0 100644 --- a/src/wigner.jl +++ b/src/wigner.jl @@ -1,19 +1,13 @@ """ - symplecticform([T = Matrix{Float64},] modes<:Int) + symplecticform([T = Matrix{Float64},] nmodes<:Int) -Compute the symplectic form matrix of size 2N x 2N, where N is given by `modes`. +Compute the symplectic form matrix of size 2N x 2N, where N is given by `nmodes`. """ -function symplecticform(modes::N) where {N<:Int} - Omega = zeros(2modes, 2modes) - @inbounds for i in 1:modes, j in modes:2modes - if isequal(i, j-modes) - Omega[i,j] = 1.0 - end - end - @inbounds for i in modes:2modes, j in 1:modes - if isequal(i-modes,j) - Omega[i, j] = -1.0 - end +function symplecticform(nmodes::N) where {N<:Int} + Omega = zeros(2*nmodes, 2*nmodes) + @inbounds for i in Base.OneTo(nmodes) + Omega[2*i-1, 2*i] = 1.0 + Omega[2*i, 2*i-1] = -1.0 end return Omega end diff --git a/test/test_randoms.jl b/test/test_randoms.jl index c8a65bc..9c21bca 100644 --- a/test/test_randoms.jl +++ b/test/test_randoms.jl @@ -1,22 +1,98 @@ @testitem "Measurements" begin using Gabs using StaticArrays + using LinearAlgebra: eigvals, adjoint - @testset "random types" begin + @testset "random utils" begin + nmodes = rand(1:20) + U = Gabs._rand_unitary(nmodes) + @test isapprox(adjoint(U), inv(U), atol = 1e-5) + + O = Gabs._rand_orthogonal_symplectic(nmodes) + @test isapprox(O', inv(O), atol = 1e-5) + Omega = symplecticform(nmodes) + @test isapprox(O' * Omega * O, Omega, atol = 1e-5) + + Spassive = randsymplectic(nmodes, passive = true) + @test isapprox(Spassive', inv(Spassive), atol = 1e-5) + @test isapprox(Spassive' * Omega * Spassive, Omega, atol = 1e-5) + + S = randsymplectic(nmodes) + @test isapprox(S' * Omega * S, Omega, atol = 1e-5) + + S_array = randsymplectic(Array, nmodes) + @test isapprox(S_array' * Omega * S_array, Omega, atol = 1e-5) + end + + @testset "random states" begin nmodes = rand(1:20) rs = randstate(nmodes) rc = randchannel(nmodes) @test rc isa GaussianChannel @test rc * rs isa GaussianState + Omega = symplecticform(nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs.covar .+ im*Omega))) + + rspure = randstate(nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rspure.covar .+ im*Omega))) + @test isapprox(purity(rspure), 1.0, atol = 1e-5) + + rs_array = randstate(Array, nmodes) + rc_array = randchannel(Array, nmodes) + @test rc_array isa GaussianChannel + @test rc_array * rs_array isa GaussianState + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rs_array.covar .+ im*Omega))) + + rspure_array = randstate(Array, nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-3)), real(eigvals(rspure_array.covar .+ im*Omega))) + @test isapprox(purity(rspure_array), 1.0, atol = 1e-3) + + rs_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + @test rc_static isa GaussianChannel + @test rc_static * rs_static isa GaussianState + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rs_static.covar) .+ im*Omega))) + + rspure_static = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes, pure = true) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rspure_static.covar) .+ im*Omega))) + @test isapprox(purity(rspure_static), 1.0, atol = 1e-5) + end + + @testset "random unitaries" begin + nmodes = rand(1:20) + ru = randunitary(nmodes) + Omega = symplecticform(nmodes) + @test isapprox(ru.symplectic' * Omega * ru.symplectic, Omega, atol = 1e-5) + + rupassive = randunitary(nmodes, passive = true) + @test isapprox(rupassive.symplectic', inv(rupassive.symplectic), atol = 1e-5) + @test isapprox(rupassive.symplectic' * Omega * rupassive.symplectic, Omega, atol = 1e-5) + + ru_array = randunitary(Array, nmodes) + @test isapprox(ru_array.symplectic' * Omega * ru_array.symplectic, Omega, atol = 1e-5) + + rupassive_array = randunitary(nmodes, passive = true) + @test isapprox(rupassive_array.symplectic', inv(rupassive_array.symplectic), atol = 1e-5) + @test isapprox(rupassive_array.symplectic' * Omega * rupassive_array.symplectic, Omega, atol = 1e-5) + + ru_static = randunitary(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) + @test isapprox(ru_static.symplectic' * Omega * ru_static.symplectic, Omega, atol = 1e-5) + + rupassive_static = randunitary(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes, passive = true) + @test isapprox(rupassive_static.symplectic', inv(rupassive_static.symplectic), atol = 1e-5) + @test isapprox(rupassive_static.symplectic' * Omega * rupassive_static.symplectic, Omega, atol = 1e-5) + end + + @testset "random channels" begin + nmodes = rand(1:20); + rc = randchannel(nmodes); + Omega = symplecticform(nmodes); + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc.noise .+ im*Omega .- im*rc.transform*Omega*rc.transform'))) - rsfloat32 = randstate(AbstractArray{Float32}, nmodes) - rcfloat32 = randchannel(AbstractArray{Float32}, nmodes) - @test rcfloat32 isa GaussianChannel - @test rcfloat32 * rsfloat32 isa GaussianState + rc_array = randchannel(Array, nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(rc_array.noise .+ im*Omega .- im*rc_array.transform*Omega*rc_array.transform'))) - rsstatic = randstate(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) - rcstatic = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes,2*nmodes}, nmodes) - @test rcstatic isa GaussianChannel - @test rcstatic * rsstatic isa GaussianState + rc_static = randchannel(SVector{2*nmodes}, SMatrix{2*nmodes, 2*nmodes}, nmodes) + @test all(i -> ((i >= 0) || isapprox(i, 0.0, atol = 1e-5)), real(eigvals(Array(rc_static.noise) .+ im*Omega .- im*Array(rc_static.transform)*Omega*Array(rc_static.transform)'))) end end \ No newline at end of file diff --git a/test/test_wigner.jl b/test/test_wigner.jl index 38664a7..aad9729 100644 --- a/test/test_wigner.jl +++ b/test/test_wigner.jl @@ -5,10 +5,10 @@ @testset "symplectic form" begin Omega = symplecticform(2) Omega_static = symplecticform(SMatrix{4, 4}, 2) - test_Omega = [0.0 0.0 1.0 0.0; + test_Omega = [0.0 1.0 0.0 0.0; + -1.0 0.0 0.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] + 0.0 0.0 -1.0 0.0] @test isequal(Omega, test_Omega) @test Omega_static isa SMatrix end