Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

diamond distance is not symmetric #53

Closed
benkj opened this issue Mar 13, 2019 · 19 comments
Closed

diamond distance is not symmetric #53

benkj opened this issue Mar 13, 2019 · 19 comments
Labels
Milestone

Comments

@benkj
Copy link
Contributor

benkj commented Mar 13, 2019

Consider the following:

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]])
println(diamond_distance(AD,AD2))
println(diamond_distance(AD2,AD))

The first line outputs ~0.17 while the second line outputs ~0.40.
This is a bug, since the diamond distance is symmetric

@lpawela
Copy link
Member

lpawela commented Mar 13, 2019

@pgawron does this indicate some problem with conversion?

@lpawela
Copy link
Member

lpawela commented Mar 13, 2019

@benkj just run the code on my machine, seems to work fine

julia> println(diamond_distance(AD,AD2))
0.3999999603826285

julia> println(diamond_distance(AD2,AD))
0.4000000254778498

@pgawron
Copy link
Collaborator

pgawron commented Mar 13, 2019

@lpawela conversion is suspicious then.

@lpawela
Copy link
Member

lpawela commented Mar 13, 2019

@pgawron can you replicate this behavior?

@pgawron
Copy link
Collaborator

pgawron commented Mar 13, 2019

@lpawela I can. On release version of the package I obtain:

AD=convert(DynamicalMatrix{Matrix{Float64}},KrausOperators([[1 0; 0 sqrt(p)], [0 sqrt(1-p); 0 0]]))
AD2=convert(DynamicalMatrix{Matrix{Float64}},KrausOperators([[1 0; 0 sqrt(2p)], [0 sqrt(1-2p); 0 0]]))

println(diamond_distance(AD,AD2))
# 0.17161841986966694
println(diamond_distance(AD2,AD))
# 0.39999997855290215

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

@pgawron still can't reproduce. Just to verify I'm on the same version

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]])
println(diamond_distance(AD,AD2))
ERROR: MethodError: no method matching DynamicalMatrix{Array{Float64,2}}(::KrausOperators{Array{Float64,2}})

as expected. Now lunching your code

AD=convert(DynamicalMatrix{Matrix{Float64}},KrausOperators([[1 0; 0 sqrt(p)], [0 sqrt(1-p); 0 0]]))
AD2=convert(DynamicalMatrix{Matrix{Float64}},KrausOperators([[1 0; 0 sqrt(2p)], [0 sqrt(1-2p); 0 0]]))
println(diamond_distance(AD,AD2))
# 0.3999999603826285
println(diamond_distance(AD2,AD))
# 0.4000000254778498

My installed packages

(v1.0) pkg> status
    Status `~/.julia/environments/v1.0/Project.toml`
  [c7e460c6] ArgParse v0.6.2
  [c52e3926] Atom v0.7.15
  [8e7c35d0] BlockArrays v0.7.0
  [f65535da] Convex v0.11.1
  [5789e2e9] FileIO v1.0.5
  [033835bb] JLD2 v0.1.2
  [e5e0dc1b] Juno v0.5.5
  [b964fa9f] LaTeXStrings v1.0.3
  [2fda8390] LsqFit v0.7.3
  [15e1cf62] NPZ v0.3.0
  [91a5bcdd] Plots v0.23.1
  [438e738f] PyCall v1.90.0
  [3c0b384b] QuantumInformation v0.4.4
  [c946c3f1] SCS v0.5.1
  [2913bbd2] StatsBase v0.29.0
  [10745b16] Statistics

@pgawron
Copy link
Collaborator

pgawron commented Mar 14, 2019

My status

(v1.1) pkg> status
    Status `~/.julia/environments/v1.1/Project.toml`
  [ca501d18] AcausalNets v0.0.0 #master (https://github.com/mikegpl/AcausalNets.jl)
  [c52e3926] Atom v0.7.15
  [b4f34e82] Distances v0.8.0
  [e5e0dc1b] Juno v0.5.5
  [c36e90e8] PowerModels v0.9.4
  [438e738f] PyCall v1.90.0
  [d330b81b] PyPlot v2.8.0
  [3c0b384b] QuantumInformation v0.4.4

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

@benkj can you post your status and julia versions?
@pgawron I don't see Convex and SCS in your packages

@pgawron
Copy link
Collaborator

pgawron commented Mar 14, 2019

@pgawron I don't see Convex and SCS in your packages

I wonder why diamond_norm even works in this case.

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

@pgawron I noticed you're using v1.1, mine is v1.0. After update I get

(v1.1) pkg> status
    Status `~/.julia/environments/v1.1/Project.toml`
  [c52e3926] Atom v0.7.15
  [15e1cf62] NPZ v0.3.0
  [438e738f] PyCall v1.90.0
  [d330b81b] PyPlot v2.8.0
  [3c0b384b] QuantumInformation v0.4.4

It simply doesn't show dependencies. The result is still the same

println(diamond_distance(AD,AD2))
# 0.3999999603826285

println(diamond_distance(AD2,AD))
# 0.4000000254778498

Maybe this behavior is somehow related to #50 and jump-dev/Convex.jl#274

@pgawron
Copy link
Collaborator

pgawron commented Mar 14, 2019

@pgawron I noticed you're using v1.1, mine is v1.0. After update I get

I've just tested for version difference. in both versions i get the same buggy result.

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

Just tested on another machine v1.0.3, still the results agree. What is going on?

@lpawela lpawela added this to the v0.5 milestone Mar 14, 2019
@benkj
Copy link
Contributor Author

benkj commented Mar 14, 2019

I'm running it on a OSX with

  [f65535da] Convex v0.11.1
  [3c0b384b] QuantumInformation v0.4.4 [`~/.julia/dev/QuantumInformation`]
  [c946c3f1] SCS v0.5.1

I find the same bug with or without conversions. This is really weird.

@benkj
Copy link
Contributor Author

benkj commented Mar 14, 2019

I changed the code putting verbose=1 in the SCSsolver and I get the following
ulia> println(diamond_distance(AD,AD2))

----------------------------------------------------------------------------
	SCS v2.0.2 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2017
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 197, CG tol ~ 1/iter^(2.00)
eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 20, rho_x = 1.00e-03
Variables n = 49, constraints m = 293
Cones:	primal zero / dual free vars: 137
	sd vars: 156, sd blks: 3
Setup time: 1.79e-04s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 2.69e+00  3.90e+01  9.88e-01 -5.19e+00  7.41e+01  0.00e+00  2.98e-04
   100| 5.59e-05  2.20e-04  5.58e-05 -2.81e-01 -2.81e-01  4.18e-16  1.14e-02
   120| 2.62e-07  1.82e-06  7.24e-07 -1.72e-01 -1.72e-01  1.70e-16  1.38e-02
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.38e-02s
	Lin-sys: avg # CG iterations: 2.67, avg solve time: 5.95e-06s
	Cones: avg projection time: 5.91e-05s
	Acceleration: avg step time: 4.21e-05s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 2.6859e-03, dist(y, K*) = 4.8922e-02, s'y/|s||y| = 1.3742e-06
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 2.6203e-07
dual res:   |A'y + c|_2 / (1 + |c|_2) = 1.8192e-06
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 7.2403e-07
----------------------------------------------------------------------------
c'x = -0.1716, -b'y = -0.1716
============================================================================
0.17161841986966694

julia> println(diamond_distance(AD2,AD))
----------------------------------------------------------------------------
	SCS v2.0.2 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012-2017
----------------------------------------------------------------------------
Lin-sys: sparse-indirect, nnz in A = 197, CG tol ~ 1/iter^(2.00)
eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00
acceleration_lookback = 20, rho_x = 1.00e-03
Variables n = 49, constraints m = 293
Cones:	primal zero / dual free vars: 137
	sd vars: 156, sd blks: 3
Setup time: 8.63e-05s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0| 2.69e+00  3.90e+01  9.88e-01 -5.19e+00  7.41e+01  0.00e+00  2.17e-04
   100| 4.76e-04  2.56e-03  2.91e-04 -3.57e+00 -3.57e+00  2.87e-14  1.07e-02
   200| 2.85e-05  1.54e-04  4.32e-04 -3.38e-01 -3.38e-01  8.08e-04  2.35e-02
   300| 3.82e-04  2.16e-03  3.11e-04 -2.16e-02 -2.12e-02  3.10e-04  3.51e-02
   400| 1.37e-03  1.18e-02  6.30e-03 -1.83e-01 -1.74e-01  2.21e-02  4.83e-02
   500| 1.33e-03  6.41e-03  5.04e-03 -1.82e-01 -1.89e-01  9.15e-17  6.07e-02
   600| 7.72e-03  5.27e-02  7.13e-03 -1.64e-01 -1.74e-01  1.64e-15  7.26e-02
   700| 4.52e-05  2.59e-04  7.39e-05 -3.60e-01 -3.60e-01  4.00e-16  8.62e-02
   800| 1.52e-04  1.02e-03  1.54e-03 -4.04e-01 -4.07e-01  2.36e-03  9.92e-02
   860| 1.50e-08  6.94e-08  5.42e-09 -4.00e-01 -4.00e-01  8.92e-17  1.08e-01
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.08e-01s
	Lin-sys: avg # CG iterations: 3.84, avg solve time: 7.58e-06s
	Cones: avg projection time: 6.47e-05s
	Acceleration: avg step time: 4.73e-05s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 3.8501e-07, dist(y, K*) = 2.9020e-09, s'y/|s||y| = 1.0184e-08
primal res: |Ax + s - b|_2 / (1 + |b|_2) = 1.4951e-08
dual res:   |A'y + c|_2 / (1 + |c|_2) = 6.9402e-08
rel gap:    |c'x + b'y| / (1 + |c'x| + |b'y|) = 5.4179e-09
----------------------------------------------------------------------------
c'x = -0.4000, -b'y = -0.4000
============================================================================
0.39999997855290215

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

@benkj @pgawron please add printlnt(J) after line 8 in convex.jl and show the output.

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

I can reproduce w somewhat similar behavior when I set eps=1e-4 in the SCSSolver(...) call. @benkj @pgawron could you set eps=1e-7 just to be sure?

@benkj
Copy link
Contributor Author

benkj commented Mar 14, 2019

eps=1e-7 in SCSSolver() seem to fix the bug for me.

@lpawela
Copy link
Member

lpawela commented Mar 14, 2019

@pgawron could you confirm this? I wonder why the default setting doesn't work. @benkj did you compile julia or downloaded the release? Maybe in the manually compiled version some defaults are different.

@benkj
Copy link
Contributor Author

benkj commented Mar 14, 2019

I have a downloaded julia for OSX. All packages are installed with Pkg.add

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants