Add quantum conditional entropy example

docs/src/examples/optimization_with_complex_variables/quantum_conditional_entropy.jl

@@ -0,0 +1,103 @@
+# # Continuity of the quantum conditional entropy
+#
+# The quantum conditional entropy is given by
+#
+# ```math
+# S(A|B)_\rho := S(\rho^{AB}) - S(\rho^{B})
+# ```
+#
+# where $S$ is the von Neumann entropy,
+#
+# ```math
+# S(\rho) := - \text{tr}(\rho \log \rho)
+# ```
+#
+# and $\rho$ is a positive semidefinite trace-1 matrix (density matrix).
+#
+# Here, $\rho^{AB}$ represents an operator on the tensor product of two finite-dimensional Hilbert spaces $A$ and $B$ (with dimensions $d_A$ and $d_B$ respectively), so we can regard $\rho_{AB}$ as a matrix on the vector space $\mathbb{C}^{d_Ad_B}$. Moreover, $\rho^B$ denotes the partial trace of $\rho^{AB}$ over the system $A$, so $\rho^B$ is a matrix on $\mathbb{C}^{d_B}$.
+#
+# One question is how much can $S(A|B)_\rho$ vary between two density matrices $\rho$ and $\sigma$ as a function of the trace-distance $\text{trdist}(\rho, \sigma) := \|\rho-\sigma\|_1 = \frac{1}{2} \text{tr}\left(\sqrt{(\rho-\sigma)^\dagger (\rho-\sigma)}\right)$ (that is, half of the nuclear norm). Here the trace distance is meaningful as it is the quantum analog to the total variation distance, and has an interpretation in terms of the maximal possible probability to distinguish between $\rho$ and $\sigma$ by measurement.
+#
+# The Alicki-Fannes-Winter (AFW) bound ([*Winter 2015*](https://arxiv.org/abs/1507.07775v6), Lemma 2) states that if $\rho$ and $\sigma$ are density matrices, then $\text{trdist}(\rho, \sigma) \leq \varepsilon \leq 1$, then
+#
+# ```math
+# | S(A|B)_\rho - S(A|B)_\sigma| \leq 2 \varepsilon \log d_A + (1 + \varepsilon) h \left(\frac{\varepsilon}{1+\varepsilon}\right)
+# ```
+#
+# where $h(x) = -x\log x  - (1-x)\log(1-x)$ is the binary entropy.
+#
+# We can illustrate this bound by computing
+#
+# ```math
+#  \max_{\rho}  S(A|B)_\rho - S(A|B)_\sigma
+# ```
+#
+# for a fixed state $\sigma$, and comparing to the AFW bound.
+#
+# We will choose $d_A=d_B=2$, and $\sigma$ as the maximally entangled state:
+#
+# ```math
+# \sigma = \frac{1}{2}\begin{pmatrix}1 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 1\\\end{pmatrix}
+# ```math
+#
+# First, we can formulate the conditional entropy in terms of the relative entropy using the relationship
+#
+# ```math
+# S(A|B)_\rho = - D(\rho^{AB} \| I^A \otimes \rho^B)
+# ```math
+#
+# where $D$ is the quantum relative entropy and $I^A$ is the $d_A$-dimensional identity matrix. Thus:
+
+using Convex
+using LinearAlgebra: I
+
+function quantum_conditional_entropy(ρ_AB, d_A, d_B)
+    ρ_B = partialtrace(ρ_AB, 1, [d_A, d_B])
+    return -quantum_relative_entropy(ρ_AB, kron(I(d_A), ρ_B))
+end
+
+# Now we setup the problem data:
+
+ϵ = 0.1
+d_A = d_B = 2
+σ_AB = 0.5 * [
+    1 0 0 1
+    0 0 0 0
+    0 0 0 0
+    1 0 0 1
+]
+
+# And we build and solve problem itself
+
+using SCS
+
+ρ_AB = HermitianSemidefinite(d_A * d_B)
+add_constraint!(ρ_AB, tr(ρ_AB) == 1)
+
+problem = maximize(
+    quantum_conditional_entropy(ρ_AB, d_A, d_B),
+    0.5 * nuclearnorm(ρ_AB - σ_AB) ≤ ϵ,
+)
+
+solve!(problem, SCS.Optimizer; silent_solver = false)
+
+# We can then check the observed difference in relative entropies:
+
+difference = evaluate(
+    quantum_conditional_entropy(ρ_AB, d_A, d_B) -
+    quantum_conditional_entropy(σ_AB, d_A, d_B),
+)
+
+# We can compare to the bound:
+h(x) = -x * log(x) - (1 - x) * log(1 - x)
+bound = 2 * ϵ * log(d_A) + (1 + ϵ) * h(ϵ / (1 + ϵ))
+
+# In fact, in this case we know the maximizer is given by
+
+ρ_max = σ_AB * (1 - ϵ) + ϵ * (I(d_A * d_B) - σ_AB) / (d_A * d_B - 1)
+
+# We can check that `ρ_AB` obtained the right value:
+
+norm(evaluate(ρ_AB) - ρ_max)
+
+# Here we see a result within the expected tolerances of SCS.

docs/styles/config/vocabularies/jump-dev/accept.txt

@@ -5,6 +5,7 @@ docstring
 [Dd]ualize
 [Dd]ualization
 [Ee]lementwise
+entropies
 [Ee]num
 injective
 nonconvex
@@ -51,8 +52,10 @@ Markowitz
 Mohan
 Namkoong
 Nemirovski
+Neumann
 Skaf
 Udell
 Vandenberghe
+von
 Watrous
 Zeng

src/reformulations/partialtranspose.jl

@@ -51,13 +51,10 @@ Returns a matrix `M` so that for any vector `v` of length `prod(dims)`,
 `M*v == vec(permutedims(reshape(v, dims), p))`.
 """
 function permutedims_matrix(dims, p)
-    d, n = prod(dims), length(dims)
-    dense = reshape(
-        PermutedDimsArray(
-            reshape(LinearAlgebra.I(d), (dims..., dims...)),
-            (p..., (n+1:2n)...),
-        ),
-        (d, d),
-    )
-    return SparseArrays.sparse(dense)
+    d = prod(dims)
+    # Generalization of https://stackoverflow.com/a/60680132
+    rows = 1:d
+    cols = vec(permutedims(reshape(rows, dims), p))
+    data = ones(Int, d)
+    return SparseArrays.sparse(rows, cols, data, d, d)
 end