From af4ff3d009eca84dfba66e6a4f0eb92fce80a662 Mon Sep 17 00:00:00 2001 From: araujoms Date: Wed, 17 Jul 2024 00:25:13 +0200 Subject: [PATCH] add entanglement entropy --- docs/make.jl | 3 +- docs/src/api.md | 1 + src/Ket.jl | 2 +- src/entanglement.jl | 100 +++++++++++++++++-- test/multilinear.jl | 238 ++++++++++++++++++++++---------------------- 5 files changed, 216 insertions(+), 128 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index c3de6b1..e832c85 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,7 +35,8 @@ makedocs(; edit_link = "master", assets = String[] ), - pages = ["Home" => "index.md", "List of functions" => "api.md"] + pages = ["Home" => "index.md", "List of functions" => "api.md"], + checkdocs = :exports ) deploydocs(; repo = "github.com/araujoms/Ket.jl", devbranch = "master", push_preview = true) diff --git a/docs/src/api.md b/docs/src/api.md index e7ef220..48d2aa0 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -31,6 +31,7 @@ conditional_entropy ```@docs schmidt_decomposition +entanglement_entropy ``` ## Measurements diff --git a/src/Ket.jl b/src/Ket.jl index 5140fa6..c265772 100755 --- a/src/Ket.jl +++ b/src/Ket.jl @@ -5,7 +5,7 @@ import GenericLinearAlgebra import LinearAlgebra as LA import Nemo import JuMP -import Hypatia +import Hypatia, Hypatia.Cones include("basic.jl") include("states.jl") diff --git a/src/entanglement.jl b/src/entanglement.jl index 0391478..3fc6f36 100644 --- a/src/entanglement.jl +++ b/src/entanglement.jl @@ -1,3 +1,10 @@ +function _equal_sizes(f, arg) + n = size(arg, 1) + d = isqrt(n) + d^2 != n && throw(ArgumentError("Subsystems are not equally-sized, please specify sizes.")) + return f(arg, [d, d]) +end + """ schmidt_decomposition(ψ::AbstractVector, dims::AbstractVector{<:Integer}) @@ -12,18 +19,97 @@ function schmidt_decomposition(ψ::AbstractVector, dims::AbstractVector{<:Intege U, λ, V = LA.svd(m) return λ, U, conj(V) end +export schmidt_decomposition """ - schmidt_decomposition(ψ::AbstractVector, dims::AbstractVector{<:Integer}) + schmidt_decomposition(ψ::AbstractVector) Produces the Schmidt decomposition of `ψ` assuming equally-sized subsystems. Returns the (sorted) Schmidt coefficients λ and isometries U, V such that kron(U', V')*`ψ` is of Schmidt form. Reference: [Schmidt decomposition](https://en.wikipedia.org/wiki/Schmidt_decomposition). """ -function schmidt_decomposition(ψ::AbstractVector) - n = length(ψ) - d = isqrt(n) - d^2 != n && throw(ArgumentError("Subsystems are not equally-sized, please specify sizes.")) - return schmidt_decomposition(ψ, [d, d]) +schmidt_decomposition(ψ::AbstractVector) = _equal_sizes(schmidt_decomposition, ψ) + +""" + entanglement_entropy(ψ::AbstractVector, dims::AbstractVector{<:Integer}) + +Computes the relative entropy of entanglement of a bipartite pure state `ψ` with subsystem dimensions `dims`. +""" +function entanglement_entropy(ψ::AbstractVector, dims::AbstractVector{<:Integer}) + length(dims) != 2 && throw(ArgumentError("Two subsystem sizes must be specified.")) + max_sys = argmax(dims) + ρ = partial_trace(ketbra(ψ), max_sys, dims) + return entropy(ρ) +end +export entanglement_entropy + +""" + entanglement_entropy(ψ::AbstractVector) + +Computes the relative entropy of entanglement of a bipartite pure state `ψ` assuming equally-sized subsystems. +""" +entanglement_entropy(ψ::AbstractVector) = _equal_sizes(entanglement_entropy, ψ) + +""" + entanglement_entropy(ρ::AbstractMatrix, dims::AbstractVector) + +Lower bounds the relative entropy of entanglement of a bipartite state `ρ` with subsystem dimensions `dims`. +""" +function entanglement_entropy(ρ::AbstractMatrix{T}, dims::AbstractVector) where {T} + LA.ishermitian(ρ) || throw(ArgumentError("State needs to be Hermitian")) + length(dims) != 2 && throw(ArgumentError("Two subsystem sizes must be specified.")) + + d = size(ρ, 1) + is_complex = (T <: Complex) + Rs = _solver_type(T) + Ts = is_complex ? Complex{Rs} : Rs + model = JuMP.GenericModel{Rs}() + + if is_complex + JuMP.@variable(model, σ[1:d, 1:d], Hermitian) + σT = partial_transpose(σ, 2, dims) + JuMP.@constraint(model, σT in JuMP.HermitianPSDCone()) + else + JuMP.@variable(model, σ[1:d, 1:d], Symmetric) + σT = partial_transpose(σ, 2, dims) + JuMP.@constraint(model, σT in JuMP.PSDCone()) + end + JuMP.@constraint(model, LA.tr(σ) == 1) + + vec_dim = Cones.svec_length(Ts, d) + ρvec = _svec(ρ, Ts) + σvec = _svec(σ, Ts) + + JuMP.@variable(model, h) + JuMP.@objective(model, Min, h / log(Rs(2))) + JuMP.@constraint(model, [h; σvec; ρvec] in Hypatia.EpiTrRelEntropyTriCone{Rs,Ts}(1 + 2 * vec_dim)) + JuMP.set_optimizer(model, Hypatia.Optimizer{Rs}) + JuMP.set_attribute(model, "verbose", false) + JuMP.optimize!(model) + return JuMP.objective_value(model), JuMP.value.(σ) +end + +""" + entanglement_entropy(ρ::AbstractMatrix) + +Lower bounds the relative entropy of entanglement of a bipartite state `ρ` assuming equally-sized subsystems. +""" +entanglement_entropy(ρ::AbstractMatrix) = _equal_sizes(entanglement_entropy, ρ) + +""" + _svec(M::AbstractMatrix, ::Type{R}) + +Produces the scaled vectorized version of a Hermitian matrix `M` with coefficient type `R`. The transformation preserves inner products, i.e., ⟨M,N⟩ = ⟨svec(M,R),svec(N,R)⟩. +""" +function _svec(M::AbstractMatrix, ::Type{R}) where {R} #the weird stuff here is to make it work with JuMP variables + d = size(M, 1) + T = real(R) + vec_dim = Cones.svec_length(R, d) + v = Vector{real(eltype(1 * M))}(undef, vec_dim) + if R <: Real + Cones.smat_to_svec!(v, 1 * M, sqrt(T(2))) + else + Cones._smat_to_svec_complex!(v, M, sqrt(T(2))) + end + return v end -export schmidt_decomposition diff --git a/test/multilinear.jl b/test/multilinear.jl index fa5ca85..07c81e2 100644 --- a/test/multilinear.jl +++ b/test/multilinear.jl @@ -1,97 +1,38 @@ @testset "Multilinear algebra" begin -@testset "Partial trace " begin - d1, d2, d3 = 2, 3, 4 - for R in [Float64, Double64, Float128, BigFloat] - for T in [R, Complex{R}] - a = randn(T, d1, d1) - b = randn(T, d2, d2) - c = randn(T, d3, d3) - ab = kron(a, b) - ac = kron(a, c) - bc = kron(b, c) - abc = kron(ab, c) - @test partial_trace(ab, [1, 2], [d1, d2])[1] ≈ tr(ab) - @test partial_trace(ab, 2, [d1, d2]) ≈ a * tr(b) - @test partial_trace(ab, 1, [d1, d2]) ≈ b * tr(a) - @test partial_trace(ab, Int64[], [d1, d2]) ≈ ab - @test partial_trace(abc, [1, 2, 3], [d1, d2, d3])[1] ≈ tr(abc) - @test partial_trace(abc, [2, 3], [d1, d2, d3]) ≈ a * tr(b) * tr(c) - @test partial_trace(abc, [1, 3], [d1, d2, d3]) ≈ b * tr(a) * tr(c) - @test partial_trace(abc, [1, 2], [d1, d2, d3]) ≈ c * tr(a) * tr(b) - @test partial_trace(abc, 3, [d1, d2, d3]) ≈ ab * tr(c) - @test partial_trace(abc, 2, [d1, d2, d3]) ≈ ac * tr(b) - @test partial_trace(abc, 1, [d1, d2, d3]) ≈ bc * tr(a) - @test partial_trace(abc, Int64[], [d1, d2, d3]) ≈ abc - end - end - for wrapper in [Symmetric, Hermitian] - M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3))) - x = Matrix(M) - @test partial_trace(M, 2, [d1, d2, d3]) ≈ partial_trace(x, 2, [d1, d2, d3]) - @test partial_trace(M, [1, 3], [d1, d2, d3]) ≈ partial_trace(x, [1, 3], [d1, d2, d3]) - end -end - -@testset "Partial transpose " begin - d1, d2, d3 = 2, 3, 4 - for R in [Float64, Double64, Float128, BigFloat] - for T in [R, Complex{R}] - a = randn(T, d1, d1) - b = randn(T, d2, d2) - c = randn(T, d3, d3) - ab = kron(a, b) - ac = kron(a, c) - bc = kron(b, c) - abc = kron(ab, c) - @test partial_transpose(ab, [1, 2], [d1, d2]) ≈ transpose(ab) - @test partial_transpose(ab, 2, [d1, d2]) ≈ kron(a, transpose(b)) - @test partial_transpose(ab, 1, [d1, d2]) ≈ kron(transpose(a), b) - @test partial_transpose(ab, Int64[], [d1, d2]) ≈ ab - @test partial_transpose(abc, [1, 2, 3], [d1, d2, d3]) ≈ transpose(abc) - @test partial_transpose(abc, [2, 3], [d1, d2, d3]) ≈ kron(a, transpose(b), transpose(c)) - @test partial_transpose(abc, [1, 3], [d1, d2, d3]) ≈ kron(transpose(a), b, transpose(c)) - @test partial_transpose(abc, [1, 2], [d1, d2, d3]) ≈ kron(transpose(a), transpose(b), c) - @test partial_transpose(abc, 3, [d1, d2, d3]) ≈ kron(ab, transpose(c)) - @test partial_transpose(abc, 2, [d1, d2, d3]) ≈ kron(a, transpose(b), c) - @test partial_transpose(abc, 1, [d1, d2, d3]) ≈ kron(transpose(a), bc) - @test partial_transpose(abc, Int64[], [d1, d2, d3]) ≈ abc - end - end - for wrapper in [Symmetric, Hermitian] - M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3))) - x = Matrix(M) - @test partial_transpose(M, 2, [d1, d2, d3]) ≈ partial_transpose(x, 2, [d1, d2, d3]) - @test partial_transpose(M, [1, 3], [d1, d2, d3]) ≈ partial_transpose(x, [1, 3], [d1, d2, d3]) - end -end - -@testset "Permute systems " begin - @testset "Vectors" begin + @testset "Partial trace " begin d1, d2, d3 = 2, 3, 4 for R in [Float64, Double64, Float128, BigFloat] for T in [R, Complex{R}] - u = randn(T, d1) - v = randn(T, d2) - w = randn(T, d3) - uv = kron(u, v) - uw = kron(u, w) - vw = kron(v, w) - uvw = kron(u, v, w) - @test permute_systems(uv, [1, 2], [d1, d2]) ≈ kron(u, v) - @test permute_systems(uv, [2, 1], [d1, d2]) ≈ kron(v, u) - @test permute_systems(uw, [2, 1], [d1, d3]) ≈ kron(w, u) - @test permute_systems(vw, [2, 1], [d2, d3]) ≈ kron(w, v) - @test permute_systems(uvw, [1, 2, 3], [d1, d2, d3]) ≈ kron(u, v, w) - @test permute_systems(uvw, [2, 3, 1], [d1, d2, d3]) ≈ kron(v, w, u) - @test permute_systems(uvw, [3, 1, 2], [d1, d2, d3]) ≈ kron(w, u, v) - @test permute_systems(uvw, [1, 3, 2], [d1, d2, d3]) ≈ kron(u, w, v) - @test permute_systems(uvw, [2, 1, 3], [d1, d2, d3]) ≈ kron(v, u, w) - @test permute_systems(uvw, [3, 2, 1], [d1, d2, d3]) ≈ kron(w, v, u) + a = randn(T, d1, d1) + b = randn(T, d2, d2) + c = randn(T, d3, d3) + ab = kron(a, b) + ac = kron(a, c) + bc = kron(b, c) + abc = kron(ab, c) + @test partial_trace(ab, [1, 2], [d1, d2])[1] ≈ tr(ab) + @test partial_trace(ab, 2, [d1, d2]) ≈ a * tr(b) + @test partial_trace(ab, 1, [d1, d2]) ≈ b * tr(a) + @test partial_trace(ab, Int64[], [d1, d2]) ≈ ab + @test partial_trace(abc, [1, 2, 3], [d1, d2, d3])[1] ≈ tr(abc) + @test partial_trace(abc, [2, 3], [d1, d2, d3]) ≈ a * tr(b) * tr(c) + @test partial_trace(abc, [1, 3], [d1, d2, d3]) ≈ b * tr(a) * tr(c) + @test partial_trace(abc, [1, 2], [d1, d2, d3]) ≈ c * tr(a) * tr(b) + @test partial_trace(abc, 3, [d1, d2, d3]) ≈ ab * tr(c) + @test partial_trace(abc, 2, [d1, d2, d3]) ≈ ac * tr(b) + @test partial_trace(abc, 1, [d1, d2, d3]) ≈ bc * tr(a) + @test partial_trace(abc, Int64[], [d1, d2, d3]) ≈ abc end end + for wrapper in [Symmetric, Hermitian] + M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3))) + x = Matrix(M) + @test partial_trace(M, 2, [d1, d2, d3]) ≈ partial_trace(x, 2, [d1, d2, d3]) + @test partial_trace(M, [1, 3], [d1, d2, d3]) ≈ partial_trace(x, [1, 3], [d1, d2, d3]) + end end - @testset "Square matrices" begin + @testset "Partial transpose " begin d1, d2, d3 = 2, 3, 4 for R in [Float64, Double64, Float128, BigFloat] for T in [R, Complex{R}] @@ -101,46 +42,105 @@ end ab = kron(a, b) ac = kron(a, c) bc = kron(b, c) - abc = kron(a, b, c) - @test permute_systems(ab, [1, 2], [d1, d2]) ≈ kron(a, b) - @test permute_systems(ab, [2, 1], [d1, d2]) ≈ kron(b, a) - @test permute_systems(ac, [2, 1], [d1, d3]) ≈ kron(c, a) - @test permute_systems(bc, [2, 1], [d2, d3]) ≈ kron(c, b) - @test permute_systems(abc, [1, 2, 3], [d1, d2, d3]) ≈ kron(a, b, c) - @test permute_systems(abc, [2, 3, 1], [d1, d2, d3]) ≈ kron(b, c, a) - @test permute_systems(abc, [3, 1, 2], [d1, d2, d3]) ≈ kron(c, a, b) - @test permute_systems(abc, [1, 3, 2], [d1, d2, d3]) ≈ kron(a, c, b) - @test permute_systems(abc, [2, 1, 3], [d1, d2, d3]) ≈ kron(b, a, c) - @test permute_systems(abc, [3, 2, 1], [d1, d2, d3]) ≈ kron(c, b, a) + abc = kron(ab, c) + @test partial_transpose(ab, [1, 2], [d1, d2]) ≈ transpose(ab) + @test partial_transpose(ab, 2, [d1, d2]) ≈ kron(a, transpose(b)) + @test partial_transpose(ab, 1, [d1, d2]) ≈ kron(transpose(a), b) + @test partial_transpose(ab, Int64[], [d1, d2]) ≈ ab + @test partial_transpose(abc, [1, 2, 3], [d1, d2, d3]) ≈ transpose(abc) + @test partial_transpose(abc, [2, 3], [d1, d2, d3]) ≈ kron(a, transpose(b), transpose(c)) + @test partial_transpose(abc, [1, 3], [d1, d2, d3]) ≈ kron(transpose(a), b, transpose(c)) + @test partial_transpose(abc, [1, 2], [d1, d2, d3]) ≈ kron(transpose(a), transpose(b), c) + @test partial_transpose(abc, 3, [d1, d2, d3]) ≈ kron(ab, transpose(c)) + @test partial_transpose(abc, 2, [d1, d2, d3]) ≈ kron(a, transpose(b), c) + @test partial_transpose(abc, 1, [d1, d2, d3]) ≈ kron(transpose(a), bc) + @test partial_transpose(abc, Int64[], [d1, d2, d3]) ≈ abc end end + for wrapper in [Symmetric, Hermitian] + M = wrapper(randn(ComplexF64, (d1 * d2 * d3, d1 * d2 * d3))) + x = Matrix(M) + @test partial_transpose(M, 2, [d1, d2, d3]) ≈ partial_transpose(x, 2, [d1, d2, d3]) + @test partial_transpose(M, [1, 3], [d1, d2, d3]) ≈ partial_transpose(x, [1, 3], [d1, d2, d3]) + end end - @testset "Rectangular matrices" begin - d1, d2, d3 = 2, 3, 4 - for R in [Float64, Double64, Float128, BigFloat] - for T in [R, Complex{R}] - a = randn(T, d1, d2) - b = randn(T, d1, d3) - c = randn(T, d2, d3) - ab = kron(a, b) - ac = kron(a, c) - bc = kron(b, c) - abc = kron(a, b, c) - @test permute_systems(ab, [1, 2], [d1 d2; d1 d3]) ≈ kron(a, b) - @test permute_systems(ab, [2, 1], [d1 d2; d1 d3]) ≈ kron(b, a) - @test permute_systems(ac, [2, 1], [d1 d2; d2 d3]) ≈ kron(c, a) - @test permute_systems(bc, [2, 1], [d1 d3; d2 d3]) ≈ kron(c, b) - @test permute_systems(abc, [1, 2, 3], [d1 d2; d1 d3; d2 d3]) ≈ kron(a, b, c) - @test permute_systems(abc, [2, 3, 1], [d1 d2; d1 d3; d2 d3]) ≈ kron(b, c, a) - @test permute_systems(abc, [3, 1, 2], [d1 d2; d1 d3; d2 d3]) ≈ kron(c, a, b) - @test permute_systems(abc, [1, 3, 2], [d1 d2; d1 d3; d2 d3]) ≈ kron(a, c, b) - @test permute_systems(abc, [2, 1, 3], [d1 d2; d1 d3; d2 d3]) ≈ kron(b, a, c) - @test permute_systems(abc, [3, 2, 1], [d1 d2; d1 d3; d2 d3]) ≈ kron(c, b, a) + @testset "Permute systems " begin + @testset "Vectors" begin + d1, d2, d3 = 2, 3, 4 + for R in [Float64, Double64, Float128, BigFloat] + for T in [R, Complex{R}] + u = randn(T, d1) + v = randn(T, d2) + w = randn(T, d3) + uv = kron(u, v) + uw = kron(u, w) + vw = kron(v, w) + uvw = kron(u, v, w) + @test permute_systems(uv, [1, 2], [d1, d2]) ≈ kron(u, v) + @test permute_systems(uv, [2, 1], [d1, d2]) ≈ kron(v, u) + @test permute_systems(uw, [2, 1], [d1, d3]) ≈ kron(w, u) + @test permute_systems(vw, [2, 1], [d2, d3]) ≈ kron(w, v) + @test permute_systems(uvw, [1, 2, 3], [d1, d2, d3]) ≈ kron(u, v, w) + @test permute_systems(uvw, [2, 3, 1], [d1, d2, d3]) ≈ kron(v, w, u) + @test permute_systems(uvw, [3, 1, 2], [d1, d2, d3]) ≈ kron(w, u, v) + @test permute_systems(uvw, [1, 3, 2], [d1, d2, d3]) ≈ kron(u, w, v) + @test permute_systems(uvw, [2, 1, 3], [d1, d2, d3]) ≈ kron(v, u, w) + @test permute_systems(uvw, [3, 2, 1], [d1, d2, d3]) ≈ kron(w, v, u) + end + end + end + + @testset "Square matrices" begin + d1, d2, d3 = 2, 3, 4 + for R in [Float64, Double64, Float128, BigFloat] + for T in [R, Complex{R}] + a = randn(T, d1, d1) + b = randn(T, d2, d2) + c = randn(T, d3, d3) + ab = kron(a, b) + ac = kron(a, c) + bc = kron(b, c) + abc = kron(a, b, c) + @test permute_systems(ab, [1, 2], [d1, d2]) ≈ kron(a, b) + @test permute_systems(ab, [2, 1], [d1, d2]) ≈ kron(b, a) + @test permute_systems(ac, [2, 1], [d1, d3]) ≈ kron(c, a) + @test permute_systems(bc, [2, 1], [d2, d3]) ≈ kron(c, b) + @test permute_systems(abc, [1, 2, 3], [d1, d2, d3]) ≈ kron(a, b, c) + @test permute_systems(abc, [2, 3, 1], [d1, d2, d3]) ≈ kron(b, c, a) + @test permute_systems(abc, [3, 1, 2], [d1, d2, d3]) ≈ kron(c, a, b) + @test permute_systems(abc, [1, 3, 2], [d1, d2, d3]) ≈ kron(a, c, b) + @test permute_systems(abc, [2, 1, 3], [d1, d2, d3]) ≈ kron(b, a, c) + @test permute_systems(abc, [3, 2, 1], [d1, d2, d3]) ≈ kron(c, b, a) + end + end + end + + @testset "Rectangular matrices" begin + d1, d2, d3 = 2, 3, 4 + for R in [Float64, Double64, Float128, BigFloat] + for T in [R, Complex{R}] + a = randn(T, d1, d2) + b = randn(T, d1, d3) + c = randn(T, d2, d3) + ab = kron(a, b) + ac = kron(a, c) + bc = kron(b, c) + abc = kron(a, b, c) + @test permute_systems(ab, [1, 2], [d1 d2; d1 d3]) ≈ kron(a, b) + @test permute_systems(ab, [2, 1], [d1 d2; d1 d3]) ≈ kron(b, a) + @test permute_systems(ac, [2, 1], [d1 d2; d2 d3]) ≈ kron(c, a) + @test permute_systems(bc, [2, 1], [d1 d3; d2 d3]) ≈ kron(c, b) + @test permute_systems(abc, [1, 2, 3], [d1 d2; d1 d3; d2 d3]) ≈ kron(a, b, c) + @test permute_systems(abc, [2, 3, 1], [d1 d2; d1 d3; d2 d3]) ≈ kron(b, c, a) + @test permute_systems(abc, [3, 1, 2], [d1 d2; d1 d3; d2 d3]) ≈ kron(c, a, b) + @test permute_systems(abc, [1, 3, 2], [d1 d2; d1 d3; d2 d3]) ≈ kron(a, c, b) + @test permute_systems(abc, [2, 1, 3], [d1 d2; d1 d3; d2 d3]) ≈ kron(b, a, c) + @test permute_systems(abc, [3, 2, 1], [d1 d2; d1 d3; d2 d3]) ≈ kron(c, b, a) + end end end end end -end #TODO add test with JuMP variables