apkille · apkille · Nov 18, 2024 · Nov 13, 2024 · Nov 13, 2024 · Nov 13, 2024
diff --git 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
@@ -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
@@ -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
@@ -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))
diff --git a/src/randoms.jl b/src/randoms.jl
@@ -1,43 +1,129 @@
 """
-    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)
 
 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
diff --git 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
@@ -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
diff --git 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