I'm not an expert of the theory of algorithm, so please forgive me if this problem is nonsense, but I want to find good algorithm to solve below problem.

There is a array of integer each of which has been assigned with random initial value such like [1,5,-3,4] or [123, 543, 295, -703, …]

You move values from one element of the array to another to achieve smaller variance. For [1,5,-3,4], the trivial final state is [2,2,2,1].

There is a constraint. The total amount of moving distance should be as small as possible. There are many paths to move [1,5,-3,4] to [2,2,2,1], but the following way should be worse one:

  1. move 1 from the fourth element to the first one.
    • [1,5,-3,4] -> [2,5,-3,3] (distance = 1 * 3)
  2. move 3 from the second element to the third one.
    • [2,5,-3,3] -> [2,2,0,3] (distance = 3 * 1)
  3. move 3 from the fourth element to the third one.
    • [2,2,0,3] -> [2,2,2,1] (distance = 2 * 1)

(total distance: 8)

because there is shorter (better) way to achieve this:

  1. move 3 from the fourth element to the third one.

    • [1,5,-3,4] -> [1,5,-1,2] (distance = 2 * 1)
  2. move 2 from the second element to the third one.

    • [1,5,-1,2] -> [1,3,1,2] (distance = 2 * 1)
  3. move 2 from the fourth element to the third one.
    • [2,2,0,3] -> [2,2,2,1] (distance = 2 * 1)

(total distance 6)

What is the algorithm which can find the shortest path to level values of array?


2 Answers 2


Your problem seems to be the same as that of calculating the Earth mover distance between the starting and target distributions. Wikipedia suggests using the Hungarian algorithm to calculate the distance.

  • $\begingroup$ Wow, thank you for quick suggestion, @YuvalFilmus . Is your answer still valid for 2-d array? $\endgroup$ Oct 27, 2014 at 3:27
  • $\begingroup$ Yes, there are multidimensional versions of the earth mover distance. $\endgroup$ Oct 27, 2014 at 3:34

The Hungarian algorithm (mentioned in Yuval Filmus's answer) is fine and will get you the exact answer in $O(n^3)$ work.

However, here's a completely different numerical analysis style method that will get you within tolerance $\epsilon$ in $O(\epsilon^{-1} n \log n)$, or even $O(\epsilon^{-1/2} n \log n)$ work.

The basic method is to:

  1. Convert the problem to a $l_1$ minimization problem,
  2. Use fast alternating iterative methods to solve it, which will take $O(\frac{1}{\epsilon})$ iterations to converge to within tolerance $\epsilon$
  3. Within each iteration, perform the necessary linear solve using graph multigrid, which will take $O(n \log n)$ work per solve.

The optimization problem

We start with a graph, where the nodes represent places where mass is (the numbers), and the edges represent the connection between nodes that mass can be moved across. If you are concerned about some numbers being negative, just add a constant to all nodes in the initial and final states to make them positive - this will not change the optimal mass moving.

The problem can be formulated as the following $l_1$ minimization problem,

\begin{align} \min_{u} & \quad \sum_i |u_i| \\ \text{such that } & \quad D^T u = p-q, \end{align}


  • $p$ and $q$ are the final and initial distributions of mass, respectively,
  • $u_i$ is the net flow across edge $i$,
  • $D$ is the graph finite difference matrix that takes the distribution on nodes and gives the differences on the edges between them.

This is saying, minimize the net total flow across edges subject to the constraint that the flow takes the initial distribution to the final distribution. Edge directions can be chosen any way so long as your direction choice for $u$ is consistent with $D$.

The method to solve it

This form of $l_1$ minimization problem has been intensely studied in the last 10 years since it arises often in large scale statistics, machine learning, seismic inverse problems, sparsity promoting optimization, and a number of other "hot" fields.

You can now use any of the modern iterative methods developed for $l_1$ minimization problems such the alternating direction method of multipliers (ADMM), or many others. Usually these modern methods will require you to solve sparse linear systems of the form, $$(D^TD + \alpha I)x = b$$ once per iteration. Since $D^TD$ is a graph laplacian, this solve can be done in $O(n \log n)$ time through graph algebraic multigrid, and is parallelizable to $\log n$.

The error for iteration $k$ in these methods usually goes as $\text{error} \sim \frac{1}{k}$. Some accelerated variants exist that can achieve $\text{error} \sim \frac{1}{k^2}$ (albeit with worse constants and more difficult implementation). Reference here. So, you need $k \sim \frac{1}{\epsilon}$ or $k \sim \frac{1}{\epsilon^{1/2}}$ iterations to ensure that your error is less than $\epsilon$.

So there you go, an efficient and thoroughly modern way to solve this problem. Probably way overkill if your problem is small, but a great way to go if it is very large. I don't know if the details of this framework for solving EMD problems have been published anywhere. I learned of the $l_1$ formulation of the problem from Justin Solomon who has used it in combination with ADMM to solve EMD problems on discrete surfaces.

  • $\begingroup$ Thank you, @NickAlger. Your answer is helpful because my problem is large and seems to be worth trying. $\endgroup$ Oct 27, 2014 at 10:54
  • $\begingroup$ I'd recommend rolling your own ADMM routine based on the method on page 41 of the Stephen Boyd paper (linked in my post above), and finding a high quality off-the-shelf algebraic multigrid code in your language for the solver. In Matlab the lean algebraic multigrid code (LAMG) is good, and the combinatorial multigrid code (CMG) is as well. I've also heard good things about PyAMG in python but never used it. $\endgroup$
    – Nick Alger
    Oct 27, 2014 at 18:58
  • $\begingroup$ Now, I'm trying your approach but can't understand your formulation. What is D? If u is flow between two nodes and p - q is increment in a node, why is D identity matrix? $\endgroup$ Feb 7, 2015 at 3:39
  • $\begingroup$ $D$ is the finite difference matrix for the array. For a 1D array of size $n$, it will be a $n$-by-$n$ sparse matrix with ones on the main diagonal, and $-1$'s on the diagonal right above it. Adding the identity is something that comes out of the ADMM method, I recommend the Boyd paper for more details. $\endgroup$
    – Nick Alger
    Feb 8, 2015 at 1:35
  • $\begingroup$ For earthmover distance in particular, I recommend Justin Solomon's paper, linked in the answer above. It has explanation as well as pseudocode that walks through the algorithm. $\endgroup$
    – Nick Alger
    Feb 8, 2015 at 1:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.