Introduce rand_ket and generate rand_dm using Ginibre ensemble method

docs/src/api.md

@@ -101,6 +101,7 @@ zero_ket
 fock
 basis
 coherent
+rand_ket
 fock_dm
 coherent_dm
 thermal_dm

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 @@
 """
 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)
 
@@ -112,17 +128,30 @@
 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"""

src/utilities.jl

@@ -47,3 +47,5 @@
 
 _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)

test/states_and_operators.jl

@@ -57,20 +57,30 @@
     @test ρ2.dims == [2, 2]
     @test entropy_vn(ρ1, base = 2) ≈ log(2, 4)
 
-    # rand_dm
-    ρ_AB = rand_dm(4, dims = [2, 2])
+    # 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
+    @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)