From 5b85b96a6e82ffef54263725cc154533209408bf Mon Sep 17 00:00:00 2001 From: Yi-Te Huang Date: Sun, 7 Jul 2024 11:40:24 +0800 Subject: [PATCH 1/4] generate `rand_dm` with Ginibre ensemble method --- src/qobj/states.jl | 27 ++++++++++++++++++++------- src/utilities.jl | 2 ++ test/states_and_operators.jl | 12 ++++++++++-- 3 files changed, 32 insertions(+), 9 deletions(-) diff --git a/src/qobj/states.jl b/src/qobj/states.jl index 6dfa36af..18256866 100644 --- a/src/qobj/states.jl +++ b/src/qobj/states.jl @@ -112,17 +112,30 @@ function maximally_mixed_dm(dimensions::Vector{Int}) end @doc raw""" - rand_dm(N::Integer; dims::Vector{Int}=[N]) + rand_dm(dimensions; rank::Int=prod(dimensions)) -Generates a random density matrix ``\hat{\rho}``, with the property to be positive semi-definite and ``\textrm{Tr} \left[ \hat{\rho} \right] = 1``. +Generate a random density matrix from Ginibre ensemble with given argument `dimensions` and `rank`, ensuring that it is positive semi-definite and trace equals to `1`. -It is also possible to specify the list of dimensions `dims` if different subsystems are present. +The `dimensions` can be either the following types: +- `dimensions::Int`: Number of basis states in the Hilbert space. +- `dimensions::Vector{Int}`: list of dimensions representing the each number of basis in the subsystems. + +The default keyword argument `rank = prod(dimensions)` (full rank). + +# References +- [J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449](https://doi.org/10.1063/1.1704292) +- [K. Życzkowski, et al., Generating random density matrices, Journal of Mathematical Physics 52, 062201 (2011)](http://dx.doi.org/10.1063/1.3595693) """ -function rand_dm(N::Integer; dims::Vector{Int} = [N]) - ρ = rand(ComplexF64, N, N) - ρ *= ρ' +rand_dm(dimensions::Int; rank::Int=prod(dimensions)) = rand_dm([dimensions], rank = rank) +function rand_dm(dimensions::Vector{Int}; rank::Int=prod(dimensions)) + N = prod(dimensions) + (rank < 1) && throw(DomainError(rank, "The argument rank must be larger than 1.")) + (rank > N) && throw(DomainError(rank, "The argument rank cannot exceed dimensions.")) + + X = _Ginibre_ensemble(N, rank) + ρ = X * X' ρ /= tr(ρ) - return QuantumObject(ρ; type = Operator, dims = dims) + return QuantumObject(ρ; type = Operator, dims = dimensions) end @doc raw""" diff --git a/src/utilities.jl b/src/utilities.jl index 3c74899e..76fbd438 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -47,3 +47,5 @@ end _get_dense_similar(A::AbstractArray, args...) = similar(A, args...) _get_dense_similar(A::AbstractSparseMatrix, args...) = similar(nonzeros(A), args...) + +_Ginibre_ensemble(n::Int, rank::Int = n) = (randn(n, rank) + 1.0im * randn(n, rank)) / sqrt(2 * n) diff --git a/test/states_and_operators.jl b/test/states_and_operators.jl index 88773159..51c43146 100644 --- a/test/states_and_operators.jl +++ b/test/states_and_operators.jl @@ -58,19 +58,27 @@ @test entropy_vn(ρ1, base = 2) ≈ log(2, 4) # rand_dm - ρ_AB = rand_dm(4, dims = [2, 2]) + ρ_AB = rand_dm([2, 2]) ρ_A = ptrace(ρ_AB, 1) ρ_B = ptrace(ρ_AB, 2) + rank = 5 + ρ_low_rank = rand_dm(10, rank = rank) + eig_val = eigenenergies(ρ_low_rank) @test tr(ρ_AB) ≈ 1.0 @test tr(ρ_A) ≈ 1.0 @test tr(ρ_B) ≈ 1.0 + @test tr(ρ_low_rank) ≈ 1.0 @test ishermitian(ρ_AB) == true @test ishermitian(ρ_A) == true @test ishermitian(ρ_B) == true + @test ishermitian(ρ_low_rank) == true @test all(eigenenergies(ρ_AB) .>= 0) @test all(eigenenergies(ρ_A) .>= 0) @test all(eigenenergies(ρ_B) .>= 0) - @test_throws DimensionMismatch rand_dm(4, dims = [2]) + @test all(isapprox.(eig_val[1:rank], 0.0, atol = 1e-10)) + @test all(eig_val[(rank+1):10] .>= 0) + @test_throws DomainError rand_dm(4, rank = rank) + @test_throws DomainError rand_dm(4, rank = 0) # bell_state, singlet_state, triplet_states, w_state, ghz_state e0 = basis(2, 0) From 37b42526af3e4f5ac30793486b90b1613528cd7b Mon Sep 17 00:00:00 2001 From: Yi-Te Huang Date: Sun, 7 Jul 2024 12:05:25 +0800 Subject: [PATCH 2/4] introduce `rand_ket` --- docs/src/api.md | 1 + src/qobj/states.jl | 18 +++++++++++++++++- test/states_and_operators.jl | 4 +++- 3 files changed, 21 insertions(+), 2 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index e7661c4a..ad36a6f3 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -101,6 +101,7 @@ zero_ket fock basis coherent +rand_ket fock_dm coherent_dm thermal_dm diff --git a/src/qobj/states.jl b/src/qobj/states.jl index 18256866..b8fff658 100644 --- a/src/qobj/states.jl +++ b/src/qobj/states.jl @@ -2,7 +2,7 @@ Functions for generating (common) quantum states. =# -export zero_ket, fock, basis, coherent +export zero_ket, fock, basis, coherent, rand_ket export fock_dm, coherent_dm, thermal_dm, maximally_mixed_dm, rand_dm export spin_state, spin_coherent export bell_state, singlet_state, triplet_states, w_state, ghz_state @@ -54,6 +54,22 @@ This state is constructed via the displacement operator [`displace`](@ref) and z """ coherent(N::Int, α::T) where {T<:Number} = displace(N, α) * fock(N, 0) +@doc raw""" + rand_ket(dimensions) + +Generate a random normalized [`Ket`](@ref) vector with given argument `dimensions`. + +The `dimensions` can be either the following types: +- `dimensions::Int`: Number of basis states in the Hilbert space. +- `dimensions::Vector{Int}`: list of dimensions representing the each number of basis in the subsystems. +""" +rand_ket(dimensions::Int) = rand_ket([dimensions]) +function rand_ket(dimensions::Vector{Int}) + N = prod(dimensions) + ψ = rand(ComplexF64, N) .- (0.5 + 0.5im) + return QuantumObject(normalize!(ψ); type = Ket, dims = dimensions) +end + @doc raw""" fock_dm(N::Int, pos::Int=0; dims::Vector{Int}=[N], sparse::Bool=false) diff --git a/test/states_and_operators.jl b/test/states_and_operators.jl index 51c43146..38ee49c3 100644 --- a/test/states_and_operators.jl +++ b/test/states_and_operators.jl @@ -57,13 +57,15 @@ @test ρ2.dims == [2, 2] @test entropy_vn(ρ1, base = 2) ≈ log(2, 4) - # rand_dm + # rand_ket and rand_dm + ψ = rand_ket(10) ρ_AB = rand_dm([2, 2]) ρ_A = ptrace(ρ_AB, 1) ρ_B = ptrace(ρ_AB, 2) rank = 5 ρ_low_rank = rand_dm(10, rank = rank) eig_val = eigenenergies(ρ_low_rank) + @test ψ' * ψ ≈ 1.0 @test tr(ρ_AB) ≈ 1.0 @test tr(ρ_A) ≈ 1.0 @test tr(ρ_B) ≈ 1.0 From f4f03f9744006d51b32a782a453e0d8296cd06ce Mon Sep 17 00:00:00 2001 From: Yi-Te Huang Date: Sun, 7 Jul 2024 12:12:22 +0800 Subject: [PATCH 3/4] format files --- src/qobj/states.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/qobj/states.jl b/src/qobj/states.jl index b8fff658..ed956c5b 100644 --- a/src/qobj/states.jl +++ b/src/qobj/states.jl @@ -142,8 +142,8 @@ The default keyword argument `rank = prod(dimensions)` (full rank). - [J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449](https://doi.org/10.1063/1.1704292) - [K. Życzkowski, et al., Generating random density matrices, Journal of Mathematical Physics 52, 062201 (2011)](http://dx.doi.org/10.1063/1.3595693) """ -rand_dm(dimensions::Int; rank::Int=prod(dimensions)) = rand_dm([dimensions], rank = rank) -function rand_dm(dimensions::Vector{Int}; rank::Int=prod(dimensions)) +rand_dm(dimensions::Int; rank::Int = prod(dimensions)) = rand_dm([dimensions], rank = rank) +function rand_dm(dimensions::Vector{Int}; rank::Int = prod(dimensions)) N = prod(dimensions) (rank < 1) && throw(DomainError(rank, "The argument rank must be larger than 1.")) (rank > N) && throw(DomainError(rank, "The argument rank cannot exceed dimensions.")) From c22f230e21dd3d479840446e5a78d45f1cc846c7 Mon Sep 17 00:00:00 2001 From: Yi-Te Huang Date: Mon, 8 Jul 2024 23:52:52 +0800 Subject: [PATCH 4/4] optimize the generation of Ginibre ensemble --- src/utilities.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utilities.jl b/src/utilities.jl index 76fbd438..5a868c73 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -48,4 +48,4 @@ end _get_dense_similar(A::AbstractArray, args...) = similar(A, args...) _get_dense_similar(A::AbstractSparseMatrix, args...) = similar(nonzeros(A), args...) -_Ginibre_ensemble(n::Int, rank::Int = n) = (randn(n, rank) + 1.0im * randn(n, rank)) / sqrt(2 * n) +_Ginibre_ensemble(n::Int, rank::Int = n) = randn(ComplexF64, n, rank) / sqrt(n)