Merge pull request #54 from benkj/master

calculate diamond norm using dual problem
iitis · Mar 15, 2019 · f0eb3d4 · f0eb3d4
2 parents b619397 + ef09910
commit f0eb3d4
Show file tree

Hide file tree

Showing 2 changed files with 46 additions and 9 deletions.
diff --git a/src/convex.jl b/src/convex.jl
@@ -2,9 +2,16 @@
 $(SIGNATURES)
 - `Φ`: DynamicalMatrix
 
-Return [diamond norm](https://arxiv.org/pdf/1004.4110.pdf) of dynamical matrix `Φ`.
+Return [diamond norm](https://arxiv.org/pdf/1207.5726.pdf) of dynamical matrix `Φ`.
 """
-function norm_diamond(Φ::DynamicalMatrix{T}) where T<:AbstractMatrix{<:Number}
+function norm_diamond(Φ::DynamicalMatrix{T}, method=:primal, eps=1e-7) where T<:AbstractMatrix{<:Number}
+	(method == :primal || method == :dual) || throw(ArgumentError("method must be either :primal or :dual"))
+
+	method == :dual ? norm_diamond_dual(Φ,eps) : norm_diamond_primal(Φ,eps)
+end
+
+
+function norm_diamond_primal(Φ::DynamicalMatrix{T}, eps) where T<:AbstractMatrix{<:Number}
     J = Φ.matrix
     # TODO: compare d1, d2 with idim, odim
     d1 = Φ.idim
@@ -21,12 +28,35 @@ function norm_diamond(Φ::DynamicalMatrix{T}) where T<:AbstractMatrix{<:Number}
     constraints += [𝕀(d2) ⊗ ρ₀ X; X' 𝕀(d2) ⊗ ρ₁] in :SDP
 
     problem = maximize(t, constraints)
-    solve!(problem, SCSSolver(verbose=0, eps=1e-7))
+    solve!(problem, SCSSolver(verbose=0, eps=eps))
     problem.optval
 end
 
-function norm_diamond(Φ::AbstractQuantumOperation{T}) where T<:AbstractMatrix{<:Number}
-    norm_diamond(DynamicalMatrix{T}(ϕ))
+
+function norm_diamond_dual(Φ::DynamicalMatrix{T}, eps) where T<:AbstractMatrix{<:Number}
+    J = Φ.matrix
+    # TODO: compare d1, d2 with idim, odim
+    d1 = Φ.idim
+    d2 = Φ.odim
+    Y₀ = ComplexVariable(d1*d2, d1*d2)
+    Y₁ = ComplexVariable(d1*d2, d1*d2)
+
+	t = 0.5*sigmamax(partialtrace(Y₀, 1, [d2,d1])) + 
+		0.5*sigmamax(partialtrace(Y₁, 1, [d2,d1]))
+	Z = [Y₀ -J; -J' Y₁ ]
+
+    constraints = [Y₀ in :SDP, Y₁ in :SDP]
+    constraints += Z+Z' in :SDP
+
+    problem = minimize(t, constraints)
+    solve!(problem, SCSSolver(verbose=0, eps=eps))
+    problem.optval
+end
+
+
+
+function norm_diamond(Φ::AbstractQuantumOperation{T}, args...) where T<:AbstractMatrix{<:Number}
+    norm_diamond(DynamicalMatrix{T}(ϕ), args...)
 end
 
 """
@@ -36,14 +66,14 @@ $(SIGNATURES)
 
 Return [diamond distance](https://arxiv.org/pdf/1004.4110.pdf) between dynamical matrices `Φ1` and `Φ2`.
 """
-function diamond_distance(Φ1::DynamicalMatrix{T}, Φ2::DynamicalMatrix{T}) where T<:AbstractMatrix{<:Number}
+function diamond_distance(Φ1::DynamicalMatrix{T}, Φ2::DynamicalMatrix{T}, args...) where T<:AbstractMatrix{<:Number}
     J1 = Φ1.matrix
     J2 = Φ2.matrix
     # TODO: Test dimnesions
     Φ = DynamicalMatrix{T}(J1-J2, Φ1.idim, Φ1.odim)
-    norm_diamond(Φ)
+    norm_diamond(Φ, args...)
 end
 
-function diamond_distance(Φ1::AbstractQuantumOperation{T}, Φ2::AbstractQuantumOperation{T}) where T<:AbstractMatrix{<:Number}
-    diamond_distance(convert(DynamicalMatrix{T},Φ1), convert(DynamicalMatrix{T},Φ2))
+function diamond_distance(Φ1::AbstractQuantumOperation{T}, Φ2::AbstractQuantumOperation{T}, args...) where T<:AbstractMatrix{<:Number}
+    diamond_distance(convert(DynamicalMatrix{T},Φ1), convert(DynamicalMatrix{T},Φ2), args...)
 end
diff --git a/test/convex.jl b/test/convex.jl
@@ -17,4 +17,11 @@ end
     @test diamond_distance(DynamicalMatrix(J1, d, d), DynamicalMatrix(J2, d, d)) ≈ 2 atol=1e-5
 end
 
+@testset "diamond distance symmetry" begin
+    p=0.2
+    AD=KrausOperators([[1 0; 0 sqrt(p)], [0 sqrt(1-p); 0 0]])
+    AD2=KrausOperators([[1 0; 0 sqrt(2p)], [0 sqrt(1-2p); 0 0]])
+    @test diamond_distance(AD,AD2) ≈ diamond_distance(AD2,AD) atol=1e-5
+end
+
 end