Do not assume only CompositeBasis represents multiple subsystems

Author: akirakyle

src/metrics.jl

@@ -32,7 +32,8 @@ T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\} = \\sum_i |λ_i|
 
 where ``λ_i`` are the eigenvalues of `rho`.
 """
-function tracenorm_h(rho::DenseOpType{B,B}) where B
+function tracenorm_h(rho::DenseOpType)
+    check_multiplicable(rho,rho)
     s = eigvals(Hermitian(rho.data))
     sum(abs.(s))
 end
@@ -77,7 +78,7 @@ T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}.
 It calls [`tracenorm`](@ref) which in turn either uses [`tracenorm_h`](@ref)
 or [`tracenorm_nh`](@ref) depending if ``ρ-σ`` is hermitian or not.
 """
-tracedistance(rho::DenseOpType{B,B}, sigma::DenseOpType{B,B}) where {B} = 0.5*tracenorm(rho - sigma)
+tracedistance(rho::DenseOpType, sigma::DenseOpType) = (check_addible(rho,sigma); check_multiplicable(rho,rho); 0.5*tracenorm(rho - sigma))
 function tracedistance(rho::AbstractOperator, sigma::AbstractOperator)
     throw(ArgumentError("tracedistance not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
 end
@@ -95,7 +96,7 @@ T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\} = \\frac{1}{2} \
 
 where ``λ_i`` are the eigenvalues of `rho` - `sigma`.
 """
-tracedistance_h(rho::DenseOpType{B,B}, sigma::DenseOpType{B,B}) where {B}= 0.5*tracenorm_h(rho - sigma)
+tracedistance_h(rho::DenseOpType, sigma::DenseOpType) = (check_addible(rho,sigma); check_multiplicable(rho,rho); 0.5*tracenorm_h(rho - sigma))
 function tracedistance_h(rho::AbstractOperator, sigma::AbstractOperator)
     throw(ArgumentError("tracedistance_h not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
 end
@@ -117,12 +118,11 @@ It uses the identity
 
 where ``σ_i`` are the singular values of `rho` - `sigma`.
 """
-tracedistance_nh(rho::DenseOpType{B1,B2}, sigma::DenseOpType{B1,B2}) where {B1,B2} = 0.5*tracenorm_nh(rho - sigma)
+tracedistance_nh(rho::DenseOpType, sigma::DenseOpType) = (check_addible(rho, sigma); 0.5*tracenorm_nh(rho - sigma))
 function tracedistance_nh(rho::AbstractOperator, sigma::AbstractOperator)
     throw(ArgumentError("tracedistance_nh not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
 end
 
-
 """
     entropy_vn(rho)
 
@@ -141,7 +141,8 @@ natural logarithm and ``0\\log(0) ≡ 0``.
 * `rho`: Density operator of which to calculate Von Neumann entropy.
 * `tol=1e-15`: Tolerance for rounding errors in the computed eigenvalues.
 """
-function entropy_vn(rho::DenseOpType{B,B}; tol=1e-15) where B
+function entropy_vn(rho::DenseOpType; tol=1e-15)
+    check_multiplicable(rho, rho)
     evals::Vector{ComplexF64} = eigvals(rho.data)
     entr = zero(eltype(rho))
     for d ∈ evals
@@ -164,9 +165,10 @@ The Renyi α-entropy of a density operator is defined as
 S_α(ρ) = 1/(1-α) \\log(Tr(ρ^α))
 ```
 """
-function entropy_renyi(rho::DenseOpType{B,B}, α::Integer=2) where B
+function entropy_renyi(rho::Operator, α::Integer=2)
     α <  0 && throw(ArgumentError("α-Renyi entropy is defined for α≥0, α≂̸1"))
     α == 1 && throw(ArgumentError("α-Renyi entropy is defined for α≥0, α≂̸1"))
+    check_multiplicable(rho,rho)
 
     return 1/(1-α) * log(tr(rho^α))
 end
@@ -186,7 +188,12 @@ F(ρ, σ) = Tr\\left(\\sqrt{\\sqrt{ρ}σ\\sqrt{ρ}}\\right),
 
 where ``\\sqrt{ρ}=\\sum_n\\sqrt{λ_n}|ψ⟩⟨ψ|``.
 """
-fidelity(rho::DenseOpType{B,B}, sigma::DenseOpType{B,B}) where {B} = tr(sqrt(sqrt(rho.data)*sigma.data*sqrt(rho.data)))
+function fidelity(rho::DenseOpType, sigma::DenseOpType)
+    check_multiplicable(rho,rho)
+    check_multiplicable(sigma,sigma)
+    check_multiplicable(rho,sigma)
+    tr(sqrt(sqrt(rho.data)*sigma.data*sqrt(rho.data)))
+end
 
 
 """
@@ -196,7 +203,9 @@ Partial transpose of rho with respect to subsystem specified by indices.
                         
 The `indices` argument can be a single integer or a collection of integers.
 """
-function ptranspose(rho::DenseOpType{B,B}, indices=1) where B<:CompositeBasis
+function ptranspose(rho::DenseOpType, indices=1)
+    length(basis_l(rho)) == length(basis_r(rho)) || throw(ArgumentError())
+    length(basis_l(rho)) > 1 || throw(ArgumentError())
     # adapted from qutip.partial_transpose (https://qutip.org/docs/4.0.2/modules/qutip/partial_transpose.html)
     # works as long as QuantumOptics.jl doesn't change the implementation of `tensor`, i.e. tensor(a,b).data = kron(b.data,a.data)
     nsys = length(basis_l(rho))
@@ -218,7 +227,7 @@ end
 
 Peres-Horodecki criterion of partial transpose.
 """
-PPT(rho::DenseOpType{B,B}, index) where B<:CompositeBasis = all(real.(eigvals(ptranspose(rho, index).data)) .>= 0.0)
+PPT(rho::DenseOpType, index) = all(real.(eigvals(ptranspose(rho, index).data)) .>= 0.0)
 
 
 """
@@ -233,7 +242,7 @@ N(ρ) = \\frac{\\|ρᵀ\\|-1}{2},
 ```
 where `ρᵀ` is the partial transpose.
 """
-negativity(rho::DenseOpType{B,B}, index) where B<:CompositeBasis = 0.5*(tracenorm(ptranspose(rho, index)) - 1.0)
+negativity(rho::DenseOpType, index) = 0.5*(tracenorm(ptranspose(rho, index)) - 1.0)
 
 
 """
@@ -246,7 +255,7 @@ N(ρ) = \\log₂\\|ρᵀ\\|,
 ```
 where `ρᵀ` is the partial transpose.
 """
-logarithmic_negativity(rho::DenseOpType{B,B}, index) where B<:CompositeBasis = log(2, tracenorm(ptranspose(rho, index)))
+logarithmic_negativity(rho::DenseOpType, index) = log(2, tracenorm(ptranspose(rho, index)))
 
 
 """
@@ -278,35 +287,16 @@ entanglement_entropy(dm(ket)) = 2 * entanglement_entropy(ket)
 By default the computed entropy is the Von-Neumann entropy, but a different
 function can be provided (for example to compute the entanglement-renyi entropy).
 """
-function entanglement_entropy(psi::Ket{B}, partition, entropy_fun=entropy_vn) where B<:CompositeBasis
-    # check that sites are within the range
-    @assert all(partition .<= length(psi.basis))
-
-    rho = ptrace(psi, partition)
-    return entropy_fun(rho)
-end
+entanglement_entropy(psi::Ket, partition, entropy_fun=entropy_vn) = entropy_fun(ptrace(psi, partition))
 
-function entanglement_entropy(rho::DenseOpType{B,B}, partition, args...) where {B<:CompositeBasis}
-    # check that sites is within the range
-    hilb = rho.basis_l
-    N = length(hilb)
-    all(partition .<= N) || throw(ArgumentError("Indices in partition must be within the bounds of the composite basis."))
-    length(partition) <= N || throw(ArgumentError("Partition cannot include the whole system."))
+function entanglement_entropy(rho::DenseOpType, partition, args...)
+    check_multiplicable(rho,rho)
+    N = length(basis_l(rho))
 
     # build the doubled hilbert space for the vectorised dm, normalized like a Ket.
-    b_doubled = hilb^2
-    rho_vec = normalize!(Ket(b_doubled, vec(rho.data)))
-
-    if partition isa Tuple
-        partition_ = tuple(partition..., (partition.+N)...)
-    else
-        partition_ = vcat(partition, partition.+N)
-    end
-
-    return entanglement_entropy(rho_vec,partition_,args...)
+    rho_vec = normalize!(Ket(basis_l(rho)^2, vec(rho.data)))
+    entanglement_entropy(rho_vec, [partition..., (partition.+N)...], args...)
 end
 
-entanglement_entropy(state::Ket{B}, partition::Number, args...) where B<:CompositeBasis =
-    entanglement_entropy(state, [partition], args...)
-entanglement_entropy(state::DenseOpType{B,B}, partition::Number, args...) where B<:CompositeBasis =
-    entanglement_entropy(state, [partition], args...)
+entanglement_entropy(state::Ket, partition::Integer, args...) = entanglement_entropy(state, [partition], args...)
+entanglement_entropy(state::DenseOpType, partition::Integer, args...) = entanglement_entropy(state, [partition], args...)

src/operators.jl

@@ -21,8 +21,7 @@ using QuantumInterface: arithmetic_binary_error, arithmetic_unary_error, addnumb
 
 Embed operator acting on a joint Hilbert space where missing indices are filled up with identity operators.
 """
-function embed(bl::CompositeBasis, br::CompositeBasis,
-               indices, op::T) where T<:DataOperator
+function embed(bl::Basis, br::Basis, indices, op::T) where T<:DataOperator
     (length(bl) == length(br)) || throw(ArgumentError("Must have length(bl) == length(br) in embed"))
     N = length(bl)
 
@@ -70,15 +69,13 @@ function embed(bl::CompositeBasis, br::CompositeBasis,
     return unpermuted_op
 end
 
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
-                index::Integer, op::T) where T<:DataOperator
-
-    N = length(basis_l)
+function embed(bl::Basis, br::Basis, index::Integer, op::DataOperator)
+    N = length(bl)
 
     # Check stuff
-    @assert length(basis_r) == N
-    basis_l[index] == op.basis_l || throw(IncompatibleBases())
-    basis_r[index] == op.basis_r || throw(IncompatibleBases())
+    @assert length(br) == N
+    bl[index] == basis_l(op) || throw(IncompatibleBases())
+    br[index] == basis_r(op) || throw(IncompatibleBases())
     check_indices(N, index)
 
     # Build data
@@ -89,17 +86,14 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
     while i > 0
         if i == index
             data = kron(data, op.data)
-            i -= length(index)
         else
-            bl = basis_l[i]
-            br = basis_r[i]
-            id = SparseMatrixCSC{Tnum}(I, dimension(bl), dimension(br))
+            id = SparseMatrixCSC{Tnum}(I, dimension(bl[i]), dimension(br[i]))
             data = kron(data, id)
-            i -= 1
         end
+        i -= 1
     end
 
-    return Operator(basis_l, basis_r, data)
+    return Operator(bl, br, data)
 end
 
 """
@@ -117,13 +111,13 @@ end
 # TODO upstream this one
 # expect(op::AbstractOperator{B,B}, state::AbstractKet{B}) where B = norm(op * state) ^ 2
 
-function expect(indices, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis}
+function expect(indices, op::AbstractOperator, state::Ket)
     N = length(basis(state))
     indices_ = complement(N, indices)
     expect(op, ptrace(state, indices_))
 end
 
-expect(index::Integer, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis} = expect([index], op, state)
+expect(index::Integer, op::AbstractOperator, state::Ket) = expect([index], op, state)
 
 """
     variance(op, state)
@@ -138,23 +132,23 @@ function variance(op::AbstractOperator, state::Ket)
     state.data'*(op*x).data - (state.data'*x.data)^2
 end
 
-function variance(indices, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis}
+function variance(indices, op::AbstractOperator, state::Ket)
     N = length(basis(state))
     indices_ = complement(N, indices)
     variance(op, ptrace(state, indices_))
 end
 
-variance(index::Integer, op::AbstractOperator, state::Ket{B}) where {B<:CompositeBasis} = variance([index], op, state)
+variance(index::Integer, op::AbstractOperator, state::Ket) = variance([index], op, state)
 
 # Helper functions to check validity of arguments
 function check_ptrace_arguments(a::AbstractOperator, indices)
-    if !isa(a.basis_l, CompositeBasis) || !isa(a.basis_r, CompositeBasis)
-        throw(ArgumentError("Partial trace can only be applied onto operators with composite bases."))
-    end
     rank = length(basis_l(a))
     if rank != length(basis_r(a))
         throw(ArgumentError("Partial trace can only be applied onto operators wich have the same number of subsystems in the left basis and right basis."))
     end
+    if rank < 2
+        throw(ArgumentError("Partial trace can only be applied to operators over at least two subsystems."))
+    end
     if rank == length(indices)
         throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use tr() instead."))
     end

src/operators_dense.jl

@@ -246,15 +246,15 @@ function exp(op::T) where {B,T<:DenseOpType{B,B}}
     return DenseOperator(op.basis_l, op.basis_r, exp(op.data))
 end
 
-function permutesystems(a::Operator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}
+function permutesystems(a::Operator, perm)
     @assert length(a.basis_l) == length(a.basis_r) == length(perm)
     @assert isperm(perm)
     data = reshape(a.data, [shape(a.basis_l); shape(a.basis_r)]...)
     data = permutedims(data, [perm; perm .+ length(perm)])
     data = reshape(data, dimension(a.basis_l), dimension(a.basis_r))
     return Operator(permutesystems(a.basis_l, perm), permutesystems(a.basis_r, perm), data)
 end
-permutesystems(a::AdjointOperator{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis} = dagger(permutesystems(dagger(a),perm))
+permutesystems(a::AdjointOperator, perm) = dagger(permutesystems(dagger(a),perm))
 
 identityoperator(::Type{S}, ::Type{T}, b1::Basis, b2::Basis) where {S<:DenseOpType,T<:Number} =
     Operator(b1, b2, Matrix{T}(I, dimension(b1), dimension(b2)))

src/operators_sparse.jl

@@ -64,7 +64,7 @@ function exp(op::T; opts...) where {B,T<:SparseOpType{B,B}}
     end
 end
 
-function permutesystems(rho::SparseOpPureType{B1,B2}, perm) where {B1<:CompositeBasis,B2<:CompositeBasis}
+function permutesystems(rho::SparseOpPureType, perm)
     @assert length(rho.basis_l) == length(rho.basis_r) == length(perm)
     @assert isperm(perm)
     shape_ = [shape(rho.basis_l); shape(rho.basis_r)]

src/pauli.jl

@@ -36,9 +36,12 @@ function choi_chi(Nl, Nr)
 end
 
 # It's possible to get better asympotic speedups using, e.g. methods from
+# https://quantum-journal.org/papers/q-2024-09-05-1461/ (see appendices)
 # https://iopscience.iop.org/article/10.1088/1402-4896/ad6499
 # https://arxiv.org/abs/2411.00526
-# https://quantum-journal.org/papers/q-2024-09-05-1461/ (see appendices)
+# https://arxiv.org/abs/2408.06206
+# https://quantumcomputing.stackexchange.com/questions/31788/how-to-write-the-iswap-unitary-as-a-linear-combination-of-tensor-products-betw/31790#31790
+# So probably using https://github.com/JuliaMath/Hadamard.jl would be best
 function _pauli_comp_convert(op, rev)
     Nl, Nr = length(basis_l(basis_l(op))), length(basis_l(basis_r(op)))
     Vl = pauli_comp(Nl)

src/sparsematrix.jl

@@ -1,4 +1,5 @@
 import Base: permutedims
+import LinearAlgebra: eigvals, svdvals
 
 function gemm_sp_dense_small(alpha, M::SparseMatrixCSC, B::AbstractMatrix, result::AbstractMatrix)
     if isone(alpha)

src/spinors.jl

@@ -215,6 +215,8 @@ function embed(bl::SumBasis, br::SumBasis, indices, op::LazyDirectSum)
 end
 # TODO: embed for multiple LazyDirectums?
 
+embed(bl::SumBasis, br::SumBasis, index::Integer, op::LazyDirectSum) = embed(bl,br,[index],op)
+
 # Fast in-place multiplication
 function mul!(result::Ket{B1},M::LazyDirectSum{B1,B2},b::Ket{B2},alpha_,beta_) where {B1,B2}
     alpha = convert(ComplexF64, alpha_)

src/superoperators.jl

@@ -64,11 +64,12 @@ sprepost(A::Operator, B::Operator) = Operator(KetBraBasis(basis_l(A), basis_r(B)
 #    return Operator(a1⊗b1, a2⊗b2, data)
 #end
 
+# https://discourse.julialang.org/t/permuteddimsarray-slower-than-permutedims/46401
 function super_tensor(A, B)
     all, alr = basis_l(basis_l(A)), basis_r(basis_l(A))
     arl, arr = basis_l(basis_r(A)), basis_r(basis_r(A))
-    bll, blr = basis_l(basis_l(A)), basis_r(basis_l(A))
-    brl, brr = basis_l(basis_r(A)), basis_r(basis_r(A))
+    bll, blr = basis_l(basis_l(B)), basis_r(basis_l(B))
+    brl, brr = basis_l(basis_r(B)), basis_r(basis_r(B))
     data = kron(B.data, A.data)
     data = reshape(data, map(dimension, (all, bll, alr, blr, arl, brl, arr, brr)))
     data = PermutedDimsArray(data, (1, 3, 2, 4, 5, 6, 7, 8))

src/time_dependent_operator.jl

@@ -201,12 +201,12 @@ function dagger(op::TimeDependentSum)
 end
 adjoint(op::TimeDependentSum) = dagger(op)
 
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, i::Integer, o::TimeDependentSum)
-    TimeDependentSum(coefficients(o), embed(basis_l, basis_r, i, static_operator(o)), o.current_time)
+function embed(bl::Basis, br::Basis, i::Integer, o::TimeDependentSum)
+    TimeDependentSum(coefficients(o), embed(bl, br, i, static_operator(o)), o.current_time)
 end
 
-function embed(basis_l::CompositeBasis, basis_r::CompositeBasis, indices, o::TimeDependentSum)
-    TimeDependentSum(coefficients(o), embed(basis_l, basis_r, indices, static_operator(o)), o.current_time)
+function embed(bl::Basis, br::Basis, indices, o::TimeDependentSum)
+    TimeDependentSum(coefficients(o), embed(bl, br, indices, static_operator(o)), o.current_time)
 end
 
 function +(A::TimeDependentSum, B::TimeDependentSum)

test/test_metrics.jl

@@ -1,5 +1,6 @@
 using Test
 using QuantumOpticsBase
+import QuantumInterface: IncompatibleBases
 using SparseArrays, LinearAlgebra
 
 @testset "metrics" begin
@@ -49,8 +50,8 @@ sigma = tensor(psi2, dagger(psi2))
 @test tracedistance(sigma, sigma) ≈ 0.
 
 rho = spinup(b1) ⊗ dagger(coherentstate(b2, 0.1))
-@test_throws ArgumentError tracedistance(rho, rho)
-@test_throws ArgumentError tracedistance_h(rho, rho)
+@test_throws IncompatibleBases tracedistance(rho, rho)
+@test_throws IncompatibleBases tracedistance_h(rho, rho)
 
 @test tracedistance_nh(rho, rho) ≈ 0.