Reference implementation of primal and dual solutions to the
This code is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. https://www.gnu.org/licenses/gpl-3.0.html
The published reference is
"Properties of the d-dimensional earth mover’s problem"
Discrete Applied Mathematics
Volume 265, 31 July 2019, Pages 128-141.
https://doi.org/10.1016/j.dam.2019.02.042
$ ./emd.shift.py
Does obj(A) == obj(A+rand(n)*10)? True
A
|A|: 121
primal obj: 4.039692169582531
dual obj: 4.039692169582532
B = A + rand(n)*10
|B|: 121
primal obj: 4.039692169582514
dual obj: 4.039692169582537
$ ./emd.extrema.py
COMMENT
Check for equality, since A and B share extreme vertices
Theory check:
Is objective A == objective B? True
A
|A|: 47
primal obj: 0.49525109284942603
dual obj: 0.4952510928494265
B
|B|: 250
primal obj: 0.4952510928494266
dual obj: 0.4952510928494265
$ ./emd.cvxmonotone.py
SANITY CHECK:
is conv(A) \subset conv(B)?: True
COMMENT
Check for /strict/ inequality, since A is strictly contained in B
Theory check:
Is objective A < objective B? True
A
|A|: 100
primal obj: 1.836273466603692
dual obj: 1.8362734666036937
B
|B|: 22
primal obj: 6.961358677491533
dual obj: 6.9613586774915355
$ ./emd.minkowski.py
Does obj(A) + obj(B) == obj(A+B)? True
A
|A|: 10
primal obj: 2.203173289573501
dual obj: 2.203173289573501
B
|B|: 11
primal obj: 2.5579575734755813
dual obj: 2.5579575734755795
C = A + B
|C|: 110
primal obj: 4.761130863049078
dual obj: 4.761130863049079
primal obj(A) + primal obj(B) 4.7611308630490825
In spite of the singular KKT matrix warnings generated by cvxopt, we use the yielded solution since the duality gap is usually small. It is expected that a better solver would find the optimal solution and work significantly faster.
Also note that the dual solution optimal value is the negative of the primal value. This is because the dual maximizes the problem; when cast as the primal it must be negated. The output of cvxopt displays the negated objective value. The final output written at the end of the process is the negative of that negated objective value, i.e., it aligns with the other solutions.
Reformulating the problem as the dual is completely unnecessary. Surprisingly, the dual solution runs /much/ faster than the primal solution. This also serves as a sanity check on things. I recall that admissible dual solutions are not unique. However, solutions found through the greedy method are integer-valued and within each vector of the dual solutions, adjacent entries differ by at most +/- 1.
$ ./emd.cvxopt.py
pcost dcost gap pres dres k/t
0: 2.5238e+00 -2.3400e+01 2e+04 4e+01 8e+00 1e+00
1: -1.5686e+00 -2.1669e+00 4e+02 1e+00 2e-01 1e-02
2: -5.8214e-01 -8.4289e-01 1e+02 4e-01 8e-02 5e-03
....
14: 1.1463e+00 1.1463e+00 8e-06 3e-08 6e-09 2e-09
Terminated (singular KKT matrix).
pcost dcost gap pres dres k/t
0: 0.0000e+00 -1.0951e+04 4e+00 1e+00 3e+03 1e+00
1: 5.8360e-02 -1.0806e+04 3e+02 1e+00 3e+03 8e+01
2: 7.3950e-01 -8.9710e+03 4e+03 9e-01 3e+03 1e+03
...
14: -1.1463e+00 -1.1464e+00 6e-05 5e-09 2e-05 3e-08
Terminated (singular KKT matrix).
==================================================
Objectives
greedy primal obj : 1.1463
greedy dual obj : 1.1463
cvxopt primal obj : 1.1463
cvxopt dual obj : 1.1463
cvxopt dual primal obj : 1.1463
cvxopt dual dual obj : 1.1464