Two Results on Slime Mold Computations
Ruben Becker1 , Vincenzo Bonifaci2 , Andreas Karrenbauer1 ,
Pavel Kolev∗1 , and Kurt Mehlhorn1

arXiv:1707.06631v2 [] 29 Mar 2018

1

Max Planck Institute for Informatics, Saarland Informatics Campus, Saarbrücken,
Germany {ruben,karrenba,pkolev,mehlhorn}@mpi-inf.mpg.de
2
Institute for the Analysis of Systems and Informatics, National Research Council of Italy
(IASI-CNR), Rome, Italy, vincenzo.bonifaci@iasi.cnr.it

Abstract
We present two results on slime mold computations. In wet-lab experiments (Nature’00) by Nakagaki
et al. the slime mold Physarum polycephalum demonstrated its ability to solve shortest path problems.
Biologists proposed a mathematical model, a system of differential equations, for the slime’s adaption
process (J. Theoretical Biology’07). It was shown that the process convergences to the shortest path (J.
Theoretical Biology’12) for all graphs. We show that the dynamics actually converges for a much wider
class of problems, namely undirected linear programs with a non-negative cost vector.
Combinatorial optimization researchers took the dynamics describing slime behavior as an inspiration
for an optimization method and showed that its discretization can ε-approximately solve linear programs
with positive cost vector (ITCS’16). Their analysis requires a feasible starting point, a step size depending
linearly on ε, and a number of steps with quartic dependence on opt/(εΦ), where Φ is the difference
between the smallest cost of a non-optimal basic feasible solution and the optimal cost (opt).
We give a refined analysis showing that the dynamics initialized with any strongly dominating point
converges to the set of optimal solutions. Moreover, we strengthen the convergence rate bounds and
prove that the step size is independent of ε, and the number of steps depends logarithmically on 1/ε and
quadratically on opt/Φ.

∗ This work has been funded by the Cluster of Excellence “Multimodal Computing and Interaction” within the Excellence
Initiative of the German Federal Government.

1

Contents
1 Introduction
1.1 The Biologically-Grounded Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 The Biologically-Inspired Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
5

2 Convergence of the Continuous Undirected Dynamics: Simple Instances
2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 The Convergence Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
9

3 Convergence of the Continuous Undirected
3.1 Existence of a Solution with Domain [0, ∞)
3.2 LP Duality . . . . . . . . . . . . . . . . . .
3.3 Convergence to Dominance . . . . . . . . .
3.4 The Equilibrium Points . . . . . . . . . . .
3.5 Convergence . . . . . . . . . . . . . . . . . .
3.6 Details of the Convergence Process . . . . .

Dynamics: General
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .

4 Improved Convergence Results: Discrete Directed Dynamics
4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Useful Lemmas . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Strongly Dominating Capacity Vectors . . . . . . . . . . . . . . .
4.4 x(k) is Close to a Non-Negative Kernel-Free Vector . . . . . . . .
4.5 x(k) is ε-Close to an Optimal Solution . . . . . . . . . . . . . . .
4.6 Proof of Theorem 4.2 . . . . . . . . . . . . . . . . . . . . . . . . .
4.7 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.8 A Simple Lower Bound . . . . . . . . . . . . . . . . . . . . . . . .

2

.
.
.
.
.
.
.
.

Instances
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

11
11
15
16
17
18
20

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

22
22
23
25
26
28
30
31
32

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

Figure 1: The experiment in [NYT00] (reprinted from there): (a) shows the maze uniformly covered by
Physarum; yellow color indicates presence of Physarum. Food (oatmeal) is provided at the locations labeled
AG. After a while the mold retracts to the shortest path connecting the food sources as shown in (b) and
(c). (d) shows the underlying abstract graph. The video [Phy] shows the experiment.

1

Introduction

We present two results on slime mold computations, one on the biologically-grounded model and one on the
biologically-inspired model. The first model was introduced by biologists to capture the slime’s apparent
ability to compute shortest paths. We show that the dynamics can actually do more. It can solve a wide
class of linear programs with nonnegative cost vectors. The latter model was designed as an optimization
technique inspired by the former model. We present an improved convergence result for its discretization.
The two models are introduced and our results are stated in Sections 1.1 and 1.2 respectively. The results
on the former model are shown in Sections 2 and 3, the results on the latter model are shown in Section 4.

1.1

The Biologically-Grounded Model

Physarum polycephalum is a slime mold that apparently is able to solve shortest path problems. Nakagaki,
Yamada, and Tóth [NYT00] report about the following experiment; see Figure 1. They built a maze, covered
it by pieces of Physarum (the slime can be cut into pieces which will reunite if brought into vicinity), and
then fed the slime with oatmeal at two locations. After a few hours the slime retracted to a path that
follows the shortest path in the maze connecting the food sources. The authors report that they repeated
the experiment with different mazes; in all experiments, Physarum retracted to the shortest path.
The paper [TKN07] proposes a mathematical model for the behavior of the slime and argues extensively
that the model is adequate. Physarum is modeled as an electrical network with time varying resistors. We
have a simple undirected graph G = (N, E) with distinguished nodes s0 and s1 modeling the food sources.
Each edge e ∈ E has a positive length ce and a positive capacity xe (t); ce is fixed, but xe (t) is a function
of time. The resistance re (t) of e is re (t) = ce /xe (t). In the electrical network defined by these resistances,
a current of value 1 is forced from s0 to s1 . For an (arbitrarily oriented) edge e = (u, v), let qe (t) be the
resulting current over e. Then, the capacity of e evolves according to the differential equation
ẋe (t) = |qe (t)| − xe (t),

(1)

where ẋe is the derivative of xe with respect to time. In equilibrium (ẋe = 0 for all e), the flow through any
edge is equal to its capacity. In non-equilibrium, the capacity grows (shrinks) if the absolute value of the flow
is larger (smaller) than the capacity. In the sequel, we will mostly drop the argument t as is customary in
the treatment of dynamical
systems. It is well-known that the electrical flow q is the feasible flow minimizing
P
energy dissipation e re qe2 (Thomson’s principle).
We refer to the dynamics above as biologically-grounded, as it was introduced by biologists to model the
behavior of a biological system. Miyaji and Ohnishi were the first to analyze convergence for special graphs
(parallel links and planar graphs with source and sink on the same face) in [MO08]. In [BMV12] convergence
was proven for all graphs. We state the result from [BMV12] for the special case that the shortest path is
unique.
3

Theorem 1.1 ([BMV12]). Assume c > 0 and that the undirected shortest path P ∗ from s0 to s1 w.r.t. the
cost vector c is unique. Assume x(0) > 0. Then x(t) in (1) converges to P ∗ . Namely, xe (t) → 1 for e ∈ P ∗
and xe → 0 for e 6∈ P ∗ as t → ∞.
[BMV12] also proves an analogous result for the undirected transportation problem; [Bon13] simplified
the argument under additional assumptions. The paper [Bon15] studies a more general dynamics and proves
convergence for parallel links.
In this paper, we extend this result to non-negative undirected linear programs
min{cT x : Af = b, |f | ≤ x},

(2)

where A ∈ Rn×m , b ∈ Rn , x ∈ Rm , c ∈ Rm
≥0 , and the absolute values are taken componentwise. Undirected
LPs can model a wide range of problems, e.g., optimization problems such as shortest path and min-cost
flow in undirected graphs, and the Basis Pursuit problem in signal processing [CDS98].
We use n for the number of rows of A and m for the number of columns, since this notation is appropriate
when A is the node-edge-incidence matrix of a graph. A vector f is feasible if Af = b. We assume that the
system Af = b has a feasible solution and that there is no non-zero f in the kernel of A with ce fe = 0 for
all e. A vector f lies in the kernel of A if Af = 0. The vector q in (1) is now the minimum energy feasible
solution



 X c
e
fe2 : Af = b ∧ fe = 0 whenever xe = 0 .
(3)
q(t) = argmin

xe (t)
f ∈Rm 
e:xe 6=0

We remark that q is unique; see Subsection 3.1.1. If A is the incidence matrix of a graph (the column
corresponding to an edge e has one entry +1, one entry −1 and all other entries are equal to zero), (2) is a
transshipment problem with flow sources and sinks encoded by a demand vector b. The condition that there
is no solution in the kernel of A with ce fe = 0 for all e states that every cycle contains at least one edge of
positive cost. In that setting, q(t) as defined by (3) coincides with the electrical flow induced by resistors of
value ce /xe (t). We can now state our first main result.

Theorem 1.2. Let c ≥ 0 satisfy cT |f | > 0 for every nonzero f in the kernel of A. Let x∗ be an optimum
solution of (2) and let X⋆ be the set of optimum solutions. Assume x(0) > 0. The following holds for the
dynamics (1) with q as in (3):
(i) The solution x(t) exists for all t ≥ 0.
(ii) The cost cT x(t) converges to cT x∗ as t goes to infinity.
(iii) The vector x(t) converges to X⋆ .
(iv) For all e with ce > 0, xe (t) − |qe (t)| converges to zero as t goes to infinity.1 If x∗ is unique, x(t)
and q(t) converge to x∗ as t goes to infinity.
Item (i) was previously shown in [SV16a] for the case of a strictly positive cost vector. The result
in [SV16a] is stated for the cost vector c = 1. The case of a general positive cost vector reduces to this
special case by rescaling the solution vector x. We stress that the dynamics (1) is biologically-grounded. It
was proposed to model a biological system and not as an optimization method. Nevertheless, it can solve a
large class of non-negative LPs. Table 1 summarizes our first main result and puts it into context.
Sections 2 and 3 are devoted to the proof of our first main theorem. For ease of exposition, we present
the proof in two steps. In Section 2, we give a proof under the following simplifying assumptions.
(A) c > 0,
(B) The basic feasible solutions of (2) have distinct cost,
(C) We start with a positive vector x(0) ∈ Xdom := { x ∈ Rn : there is a feasible f with |f | ≤ x }.
Section 2 generalizes [Bon13]. For the undirected shortest path problem, condition (B) states that all simple
undirected source-sink paths have distinct cost and condition (C) states that all source-sink cuts have a
capacity of at least one at time zero (and hence at all times). The existence of a solution with domain [0, ∞)
was already shown in [SV16a]. We will show that Xdom is an invariant set, i.e., the solution stays in Xdom
1 We

conjecture that this also holds for the indices e with ce = 0.

4

Problem

Reference
[MO08]
[BMV12]
[SV16a]

Undirected Shortest Path
Undirected Shortest Path
Undirected Positive LP

Existence
of Solution
Yes
Yes
Yes

Convergence
to OPT
Yes
Yes
No

Yes

Yes

Undirected Nonnegative LP

Our Result

Comments
parallel edges, planar graphs
all graphs
c>0
1) c ≥ 0
2) ∀v ∈ ker(A) : cT |v| > 0

Table 1: Convergence results for the continuous undirected Physarum dynamics (1).
P
P
for all times, and that E(x) = e re x2e = e ce xe is a Lyapunov function for the dynamics (1), i.e., Ė ≤ 0
and Ė = 0 if and only if ẋ = 0. It follows from general theorems about dynamical systems that the dynamics
converges to a fixed point of (1). The fixed points are precisely the vectors |f | where f is a feasible solution
of (2). A final argument establishes that the dynamics converges to a fixed point of minimum cost.
In Section 3, we prove the general case of the first main theorem. We assume
(D) c ≥ 0,
(E) cost(z) = cT |z| > 0 for every nonzero vector z in the kernel of A,
(F) We start with a positive vector x(0) > 0.
Section 3 generalizes [BMV12] in two directions. First, we treat general undirected LPs and not just the
undirected shortest path problem, respectively, the transshipment problem. Second, we replace the condition
c > 0 by the requirement c ≥ 0 and every nonzero vector in the kernel of A has positive cost. For the
undirected shortest path problem, the latter condition states that the underlying undirected graph has
no zero-cost cycle. Section 3 is technically considerably more difficult than Section 2. We first establish
the existence of a solution with domain [0, ∞). To this end, we derive a closed formula for the minimum
energy feasible solution and prove that the mapping x 7→ q is locally-Lipschitz. Existence of a solution with
domain [0, ∞) follows by standard arguments. We then show that Xdom is an attractor, i.e., the solution
x(t) converges to Xdom . We next characterize equilibrium points and exhibit a Lyapunov function. The
Lyapunov function is a normalized version of E(x). The normalization factor is equal to the optimal value
of the linear program max { α : Af = αb, |f | ≤ x } in the variables f and α. Convergence to an equilibrium
point follows from the existence of a Lyapunov function. A final argument establishes that the dynamics
converges to a fixed point of minimum cost.

1.2

The Biologically-Inspired Model

Ito et al. [IJNT11] initiated the study of the dynamics
ẋ(t) = q(t) − x(t).

(4)

We refer to this dynamics as the directed dynamics in contrast to the undirected dynamics (1). The directed
dynamics is biologically-inspired – the similarity to (1) is the inspiration. It was never claimed to model the
behavior of a biological system. Rather, it was introduced as a biologically-inspired optimization method.
The work in [IJNT11] shows convergence of this directed dynamics (4) for the directed shortest path problem
and [JZ12, SV16c, Bon16] show convergence for general positive linear programs, i.e., linear programs with
positive cost vector c > 0 of the form
min{cT x : Ax = b, x ≥ 0}.

(5)

The discrete versions of both dynamics define sequences x(t) , t = 0, 1, 2, . . . through
x(t+1) = (1 − h(t) )x(t) + h(t) q (t)
(t+1)

x

(t)

(t)

= (1 − h )x

(t)

+ h |q

(t)

|

5

discrete directed dynamics;

(6)

discrete undirected dynamics,

(7)

where h(t) is the step size and q (t) is the minimum energy feasible solution as in (3). For the discrete
dynamics, we can ask complexity questions. This is particularly relevant for the discrete directed dynamics
as it was designed as an biologically-inspired optimization method.
For completeness, we review the state-of-the-art results for the discrete undirected dynamics. For the
undirected shortest path problem, the convergence of the discrete undirected dynamics (7) was shown
in [BBD+ 13]. The convergence proof gives an upper bound on the step size and on the number of steps required until an ε-approximation of the optimum is obtained. [SV16b] extends the result to the transshipment
problem and [SV16a] further generalizes the result to the case of positive LPs. The paper [SV16b] is related
to our first result. It shows convergence of the discretized undirected dynamics (7), we show convergence of
the continuous undirected dynamics (1) for a more general cost vector.
We come to the discrete directed dynamics (6). Similarly to the undirected setting, Becchetti et
al. [BBD+ 13] showed the convergence of (6) for the shortest path problem. Straszak and Vishnoi extended
the analysis to the transshipment problem [SV16b] and positive LPs [SV16c].
Theorem 1.3. [SV16c, Theorem 1.3] Let A ∈ Zn×m have full row rank (n ≤ m), b ∈ Zn , c ∈ Zm
>0 , and let
DS := max{| det(M )| : M is a square sub-matrix of A}.2 Suppose the Physarum dynamics (6) is initialized
with a feasible point x(0) of (5) such√that M −1 ≤ x(0) ≤ M and cT x(0) ≤ M · opt, for some M ≥ 1. Then,
for any ε > 0 and step size h ≤ ε/( 6||c||1 DS )2 , after k = O((εh)−2 ln M ) steps, x(k) is a feasible solution
with cT x(k) ≤ (1 + ε)opt.
Theorem 1.3 gives an algorithm that computes a (1 + ε)-approximation to the optimal cost of (5). In
comparison to [BBD+ 13, SV16b], it has several shortcomings. First, it requires a feasible starting point.
Second, the step size depends linearly on ε. Third, the number of steps required to reach an ε-approximation
has a quartic dependence on opt/(εΦ). In contrast, the analysis in [BBD+ 13, SV16b] yields a step size
independent of ε and a number of steps that depends only logarithmically on 1/ε, see Table 2.
We overcome these shortcomings. Before we can state our result, we need some notation. Let X⋆ be
the set of optimal solutions to (5). The distance of a capacity vector x to X⋆ is defined as dist(x, X⋆ ) :=
inf{||x − x′ ||∞ : x′ ∈ X⋆ }. Let γA := gcd({Aij : Aij 6= 0}) ∈ Z>0 and
D := max {|det (M )| : M is a square submatrix of A/γA with dimension n − 1 or n} .

(8)

Let N be the set of non-optimal basic feasible solution of (5) and
Φ := min cT g − opt ≥ 1/(DγA )2 ,
g∈N

(9)

where the inequality is well known [PS82, Lemma 8.6]. For completeness, we present a proof in Subsection 4.5.
Informally, our second main result establishes the following properties of the Physarum dynamics (6):
(i) For any ε > 0 and any strongly dominating starting point3 x(0) , there is a fixed step size h(x(0) ) such
that the Physarum dynamics (6) initialized with x(0) and h(x(0) ) converges to X⋆ , i.e., dist(x(k) , X⋆ ) <
ε/(DγA ) for large enough k.
(ii) The step size can be chosen independently of ε.
(iii) The number of steps k depends logarithmically on 1/ε and quadratically on opt/Φ.
(iv) The efficiency bounds depend on a scale-invariant determinant4 D.
In Section 4.8, we establish a corresponding lower bound. We show that for the Physarum dynamics (6)
to compute a point x(k) such that dist(x(k) , X⋆ ) < ε, the number of steps required for computing an εapproximation has to grow linearly in opt/(hΦ) and ln(1/ε), i.e. k ≥ Ω(opt · (hΦ)−1 · ln(1/ε)). Table 2 puts
our results into context.
2 Using Lemma 3.1, the dependence on D can be improved to a scale-independent determinant D, defined in (8). For
S
further details, we refer the reader to Subsection 4.2.
3 We postpone the definition of strongly dominating capacity vector to Section 4.3. Every scaled feasible solution is strongly
dominating. In the shortest
P path problem,
Pa capacity vector x is strongly dominating if every source-sink cut (S, S) has positive
directed capacity, i.e.,
x −
x > 0.
a∈E(S,S) a
a∈E(S,S) a
4 Note that (γ )n−1 D ≤ D ≤ (γ )n D, and thus D yields an exponential improvement over D , whenever γ ≥ 2.
A
S
A
S
A

6

Reference

Problem

h step size

[BBD+ 13]

Shortest Path

indep. of ε

[SV16b]

Transshipment

indep. of ε

[SV16c]

Positive LP

depends on ε

Our Result

Positive LP

indep. of ε

Lower Bound

Positive LP

indep. of ε

k number of steps
poly(m, n, ||c||1 , ||x(0) ||1 )
· ln(1/ε)
poly(m, n, ||c||1 , ||b||1 , ||x(0) ||1 )
· ln(1/ε)

poly(||c||1 , DS , ln ||x(0) ||1 )
· 1/(Φε)4
poly(||c||1 , ||b||1 , D, γA , ln ||x(0) ||1 )
·Φ−2 ln(1/ε)
Ω(opt · (hΦ)−1 ln(1/ε))

Guarantee
dist(x(k) , X⋆ ) < ε
dist(x(k) , X⋆ ) < ε
cT x(k) ≤ (1 + ε)opt
cT x(k) < ming∈N cT g
dist(x(k) , X⋆ ) <

ε
DγA

dist(x(k) , X⋆ ) < ε

Table 2: Convergence results for the discrete directed Physarum dynamics (6).
We state now our second main result for the special case of a feasible starting point, and we provide the
full version in Theorem 4.2 which applies for arbitrary strongly dominating starting point, see Section 4. We
use the following constants in the statement of the bounds.
(i) h0 := cmin /(4D||c||1 ), where cmin := mini {ci };
(ii) Ψ(0) := max{mD2 kb/γAk1 , kx(0) k∞ };
2
(iii) C1 := Dkb/γA k1 ||c||1 , C2 := 82 m2 nD5 γA
||A||∞ kbk1 and C3 := D3 γA ||b||1 ||c||1 .

Theorem 1.4. Suppose A ∈ Zn×m has full row rank (n ≤ m), b ∈ Zn , c ∈ Zm
>0 and ε ∈ (0, 1). Given a
feasible starting point x(0) > 0 the Physarum dynamics (6) with step size h ≤ (Φ/opt) · h20 /2 outputs for any
(0)
k ≥ 4C1 /(hΦ) · ln(C2 Ψ(0) /(ε · min{1, xmin })) a feasible x(k) > 0 such that dist(x(k) , X⋆ ) < ε/(DγA ).
We stated the bounds on h in terms of the unknown quantities Φ and opt. However, Φ/opt ≥ 1/C3 by
Lemma 3.1 and hence replacing Φ/opt by 1/C3 yields constructive bounds for h. Note that the upper bound
on the step size does not depend on ε and that the bound on the number of iterations depends logarithmically
on 1/ε and quadratically on opt/Φ.
What can be done if the initial point is not strongly dominating? For the transshipment problem it suffices
to add an edge of high capacity and high cost from every source node to every sink node [BBD+ 13, SV16b].
This will make the instance strongly dominating and will not affect the optimal solution. We generalize this
observation to positive linear programs. We add an additional column equal to b and give it sufficiently
high capacity and cost. This guarantees that the resulting instance is strongly dominating and the optimal
solution remains unaffected. Moreover, our approach generalizes and improves upon [SV16b, Theorem 1.2],
see Section 4.7.
Proof Techniques: The crux of the analysis in [IJNT11, BBD+ 13, SV16b] is to show that for large enough
k, x(k) is close to a non-negative flow f (k) and then to argue that f (k) is close to an optimal flow f ⋆ . This
line of arguments yields a convergence of x(k) to X⋆ with a step size h chosen independently of ε.
In Section 4, we extend the preceding approach to positive linear programs, by generalizing the concept of
non-negative cycle-free flows to non-negative feasible kernel-free vectors (Subsection 4.4). Although, we use
the same high level ideas as in [BBD+ 13, SV16b], we stress that our analysis generalizes all relevant lemmas
in [BBD+ 13, SV16b] and it uses arguments from linear algebra and linear programming duality, instead
of combinatorial arguments. Further, our core efficiency bounds (Subsection 4.2) extend [SV16c] and yield
a scale-invariant determinant dependence of the step size and are applicable for any strongly dominating
point (Subsection 4.3).

7

2

Convergence of the Continuous Undirected Dynamics: Simple
Instances

We prove Theorem 1.2 under the following simplifying assumptions:
(A) c > 0,
(B) The basic feasible solutions of (2) have distinct cost,
(C) We start with a positive vector x(0) ∈ Xdom := { x ∈ Rn : there is a feasible f with |f | ≤ x } .

This section generalizes [Bon13]. For the undirected shortest path problem, condition (B) states that all
simple undirected source-sink paths have distinct cost and condition (C) states that all source-sink cuts have
a capacity of at least one at time zero (and hence at all times). The existence of a solution with domain
[0, ∞) was already shown in [SV16a].

2.1

Preliminaries

Note that we may assume that A has full row-rank since any equation that is linearly dependent on other
equations can be deleted without changing the feasible set. We continue to use n and m for the dimension
of A. Thus, A has rank n. We continue by fixing some terms and notation. A basic feasible solution of (2) is
a pair of vectors x and f = (fB , fN ), where fB = A−1
B b and AB is a square n × n non-singular sub-matrix of
A and fN = 0 is the vector indexed by the coordinates not in B, and x = |f |. Since f uniquely determines
x, we may drop the latter for the sake of brevity and call f a basic feasible solution of (2). A feasible
solution f is kernel-free or non-circulatory if it is contained in the convex hull of the basic feasible solutions.5
We say that a vector f ′ is sign-compatible with a vector f (of the same dimension) or f -sign-compatible if
fe′ 6= 0 implies fe′ fe > 0. In particular, supp(f ′ )P
⊆ supp(f ). For a given capacity vector x and a vector
2
f ∈ Rm with supp(f ) ⊆ supp(x), we use E(f ) = e (c
the energy of f . The energy of f
Pe /xe )fe to denote
T
c
|f
|
=
c
|f
|
to
denote the cost of f . Note that
is infinite,
if
supp(f
)
⊆
6
supp(x).
We
use
cost(f
)
=
e
e
e
P
P
E(x) = e (ce /xe )x2e = e ce xe = cost(x). We define the constants cmax = kck∞ and cmin = mine:ce >0 ce .
We use the following corollary of the finite basis theorem for polyhedra.
Lemma 2.1. Let f be a feasible solution of (2). Then f is the sum of a convex combination of at most
m basic feasible solutions plus a vector in the kernel of A. Moreover, all elements in this representation are
sign-compatible with f .
Proof. We may assume f ≥ 0. Otherwise, we flip the sign of the appropriate columns of A. Thus, the system
Af = b, f ≥ 0 is feasible and f is the sum of a convex combination of at most m basic feasible solutions plus
a vector in the kernel of A by the finite basis theorem [Sch99, Corollary 7.1b]. By definition, the elements
in this representation are non-negative vectors and hence sign-compatible with f .
Fact 2.2 (Grönwall’s Lemma). Let A, B, α, β ∈ R, α 6= 0, β 6= 0, and let g be a continuous differentiable
function on [0, ∞). If A + αg(t) ≤ ġ(t) ≤ B + βg(t) for all t ≥ 0, then −A/α + (g(0) + A/α)eαt ≤ g(t) ≤
−B/β + (g(0) + B/β)eβt for all t ≥ 0.
Proof. We show the upper bound. Assume first that B = 0. Then
ġeβt − βgeβt
d g
=
≤ 0 implies
dt eβt
e2βt

g(t)
g(0)
≤ β0 = g(0).
eβt
e

If B 6= 0, define h(t) = g(t) + B/β. Then
ḣ = ġ ≤ B + βg = B + β(h − B/β) = βh
and hence h(t) ≤ h(0)eβt . Therefore g(t) ≤ −B/β + (g(0) + B/β)eβt .
5 For the undirected shortest path problem, we drop the equation corresponding to the sink. Then b becomes the negative
indicator vector corresponding to the source node. Note that n is one less than the number of nodes of the graph. The basic
feasible solutions are the simple undirected source-sink paths. A circulatory solution contains a cycle on which there is flow.

8

An immediate consequence of Grönwall’s Lemma is that the undirected Physarum dynamics (1) initialized
with any positive starting vector x(0), generates a trajectory {x(t)}t≥0 such that each time state x(t) is a
positive vector. Indeed, since ẋe = |qe | − xe ≥ −xe , we have xe (t) ≥ xe (0) · exp{−t} for every index e with
xe (0) > 0 and every time t. Further, by (1) and (3), it holds for indices e with xe (0) = 0 that xe (t) = 0 for
every time t. Hence, the trajectory {x(t)}t≥0 has a time-invariant support.
Fact 2.3 ([JZ12]). Let R = diag(ce /xe ). Then q = R−1 AT p, where p = (AR−1 AT )−1 b.
P
Proof. q minimizes e re qe2 subject to Aq = b. The KKT conditions imply the existence of a vector p such
that Rq = AT p. Substituting into Aq = b yields p = (AR−1 AT )−1 b.
Lemma 2.4. Xdom is an invariant set, i.e., if x(0) ∈ Xdom then x(t) ∈ Xdom for all t.
Proof. Let q(t) be the minimum energy feasible solution with respect to R(t) = diag(ce /xe (t)), and let f (t)
d
be such that f (0) is feasible, |f (0)| ≤ x(0), and f˙(t) = q(t) − f (t). Then dt
(Af − b) = A(q − f ) = b − Af
−t
and hence Af (t) − b = (Af (0) − b)e = 0. Thus f (t) is feasible for all t. Moreover,
d
(f − x) = f˙ − ẋ = q − f − (|q| − x) = q − |q| − (f − x) ≤ −(f − x).
dt
Thus f (t) − x(t) ≤ (f (0) − x(0))e−t ≤ 0 by Gronwall’s Lemma applied with g(t) = f (t) − x(t) and β = −1,
and hence f (t) ≤ x(t) for all t. Similarly,
d
(f + x) = f˙ + ẋ = q − f + (|q| − x) = q + |q| − (f + x) ≥ −(f + x).
dt
Thus f (t) + x(t) ≥ (f (0) + x(0))e−t ≥ 0 by Grönwall’s Lemma applied with g(t) = f (t) + x(t) and α = −1
and A = 0, and hence f (t) ≥ −x(t) for all t.
We conclude that |f (t)| ≤ x(t) for all t. Thus, x(t) ∈ Xdom for all t.

2.2

The Convergence Proof

We will first characterize the equilibrium points. They are precisely the points |f |, where f is a basic feasible
solution; the proof uses property (B). We then show that E(x) is a Lyapunov function for (1), in particular,
Ė ≤ 0 and Ė = 0 if and only if x is an equilibrium point. For this argument, we need that the energy of q is
at most the energy of x with equality if and only if x is an equilibrium point. This proof uses (A) and (C).
It follows from the general theory of dynamical systems that x(t) approaches an equilibrium point. Finally,
we show that convergence to a non-optimal equilibrium is impossible.
Lemma 2.5 (Generalization of Lemma 2.3 in [Bon13]). Assume (A) to (C). If f is a basic feasible solution
of (2), then x = |f | is an equilibrium point. Conversely, if x is an equilibrium point, then x = |f | for some
basic feasible solution f .
Proof. Let f be a basic feasible solution, let x = |f |, and let q be the minimum energy feasible solution with
respect to the resistances ce /xe . We have Aq = b and supp(q) ⊆ supp(x) by definition of q. Since f is a
basic feasible solution there is a subset B of size n of the columns of A such that AB is non-singular and
f = (A−1
B b, 0). Since supp(q) ⊆ supp(x) ⊆ B, we have q = (qB , 0) for some vector qB . Thus, b = Aq = AB qB
and hence qB = fB . Therefore ẋ = |q| − x = 0 and x is an equilibrium point.
Conversely, if x is an equilibrium point, |qe | = xe for every e. By changing the signs of some columns of
A, we may assume q ≥ 0. Then q = x. Since qe = xe /ce ATe p where Ae is the e-th column of A by Lemma 2.3,
we have ce = ATe p, whenever xe > 0. By Lemma 2.1, q is a convex combination of basic feasible solutions
and a vector in the kernel of A that are sign-compatible with q. The vector in the kernel must be zero as q is
a minimum energy feasible solution. For any basic feasible solution
to q, we have supp(z) ⊆
P z contributing P
supp(x). Summing over the e ∈ supp(z), we obtain cost(z) = e∈supp(z) ce ze = e∈supp(z) ze ATe p = bT p.
Thus, the convex combination involves only a single basic feasible solution by assumption (B) and hence x
is a basic feasible solution.
9

The vector x(t) dominates a feasible solution at all times. Since q(t) is the minimum energy feasible
solution at time t, this implies E(q(t)) ≤ E(x(t)) at all times. A further argument shows that we have
equality if and only if x = |q|.

Lemma 2.6 (Generalization of Lemma 3.1 in [Bon13]). Assume (A) to (C). At all times, E(q) ≤ E(x). If
E(q) = E(x), then x = |q|.

Proof. Recall that x(t) ∈ Xdom for all t. Thus, at all times, there is a feasible f such that |f | ≤ x. Since q
is a minimum energy feasible solution, we have
E(q) ≤ E(f ) ≤ E(x).
If E(q) = E(x) then E(q) = E(f ) and hence q = f since the minimum energy feasible solution is unique.
Also, |f | = x since |f | ≤ x and |fe | < xe for some e implies E(f ) < E(x). The last conclusion uses c > 0.
Lyapunov functions are the main tool for proving convergence of dynamical systems. We show that E(x)
is a Lyapunov function for (1).
Lemma 2.7 (Generalization of Lemma 3.2 in [Bon13]). Assume (A) to (C). E(x) is a Lyapunov function
for (1), i.e., it is continuous as a function of x, E(x) ≥ 0, Ė(x) ≤ 0 and Ė(x) = 0 if and only if ẋ = 0.

Proof. E is clearly continuous and non-negative. Recall that E(x) = cost(x). Let R be the diagonal matrix
with entries ce /xe . Then
d
cost(x) = cT (|q| − x)
dt
= xT R|q| − xT Rx

by (1)
since c = Rx

= xT R1/2 R1/2 |q| − xT Rx

≤ (q T Rq)1/2 (xT Rx)1/2 − xT Rx
T

1/2

≤ (x Rx)

= 0.

T

1/2

(x Rx)

T

− x Rx

by Cauchy-Schwarz
by Lemma 2.6

d
cost(x) = 0 implies that both inequalities above are equalities. This is only possible if the
Observe that dt
vectors |q| and x are parallel and E(q) = E(x). Thus, x = |q| by Lemma 2.6.

It follows now from the general theory of dynamical systems that x(t) converges to an equilibrium point.
Corollary 2.8 (Generalization of Corollary 3.3. in [Bon13].). Assume (A) to (C). As t → ∞, x(t) approaches an equilibrium point and cT x(t) approaches the cost of the corresponding basic feasible solution.
Proof. The proof in [Bon13] carries over. We include it for completeness. The existence of a Lyapunov
function E implies by [LaS76, Corollary 2.6.5] that x(t) approaches the set x ∈ Rm
≥0 : Ė = 0 , which by

:
ẋ
=
0
.
Since
this
set
consists
of
isolated
points (Lemma 2.5),
Lemma 2.7 is the same as the set x ∈ Rm
≥0
x(t) must approach one of those points, say the point x0 . When x = x0 , one has E(q) = E(x) = cost(x) =
cT x.
It remains to exclude that x(t) converges to a nonoptimal equilibrium point.
Theorem 2.9 (Generalization of Theorem 3.4 in [Bon13]). Assume (A) to (C). As t → ∞, cT x(t) converges
to the cost of the optimal solution and x(t) converges to the optimal solution.
Proof. By the corollary, it suffices to prove the second part of the claim. For the second P
part, assume that
x(t) converges to a non-optimal solution z. Let x∗ be the optimal solution and let W = e x∗e ce ln xe . Let
δ = (cost(z) − cost(x∗ ))/2. Note that for all sufficiently large t, we have E(q(t)) ≥ cost(z) − δ ≥ cost(x∗ ) + δ.
Further, by definition qe = (xe /ce )ATe p and thus
Ẇ =

X
e

x∗e ce

X
X
|qe | − xe
=
x∗e |ATe p| − cost(x∗ ) ≥
x∗e ATe p − cost(x∗ ) ≥ δ,
xe
e
e
10

where the last inequality follows by
to the fact that x is bounded.

3

P

e

x∗e ATe p = bT p = E(q) ≥ cost(x∗ ) + δ. Hence W → ∞, a contradiction

Convergence of the Continuous Undirected Dynamics: General
Instances

We now prove the general case. We assume
(D) c ≥ 0,
(E) cost(z) > 0 for every nonzero vector z in the kernel of A,
(F) We start with a positive vector x(0) > 0.
In this section, we generalize [BMV12] in two directions. First, we treat general undirected LPs and not just
the undirected shortest path problem, respectively, the transshipment problem. Second, we substitute the
condition c > 0 with c ≥ 0 and every nonzero vector in the kernel of A has positive cost. For the undirected
shortest path problem, the latter condition states that the underlying undirected graph has no zero-cost
cycle.

3.1

Existence of a Solution with Domain [0, ∞)

In this section we show that a solution x(t) to (1) has domain [0, ∞). We first derive an explicit formula for
the minimum energy feasible solution q and then show that the mapping x 7→ q is Lipschitz continuous; this
implies existence of a solution with domain [0, ∞) by standard arguments.
3.1.1

The Minimum Energy Solution

Recall that for γA = gcd({Aij : Aij 6= 0}) ∈ Z>0 , we defined by
D = max {|det (M )| : M is a square submatrix of A/γA with dimension n − 1 or n} .
We derive now properties of the minimum energy solution. In particular, if every nonzero vector in the
kernel of A has positive cost,
(i) the minimum energy feasible solution is kernel-free and unique (Lemma 3.2),
(ii) |qe | ≤ D||b/γA ||1 for every e ∈ [m] (Lemma 3.3),
(iii) q is defined by (12) (Lemma 3.4), and
(iv) E(q) = bT p, where p is defined by (12) (Lemma 3.5).
We note that for positive cost vector c > 0, these results are known.
We proceed by establishing some useful properties on basic feasible solutions.
Lemma 3.1. Suppose A ∈ Zn×m is an integral matrix, and b ∈ Zn is an integral vector. Then, for any
basic feasible solutions f with Af = b and f ≥ 0, it holds that ||f ||∞ ≤ D||b/γA ||1 and fj 6= 0 implies
|fj | ≥ 1/(DγA ).
Proof. Since f is a basic feasible solution, it has the form f = (fB , 0) such that AB ·fB = b where AB ∈ Zn×n
is an invertible submatrix of A. We write M−i,−j to denote the matrix M with deleted i-th row and j-th
column. Let Qj be the matrix formed by replacing the j-th column of AB by the column vector b. Then,
using the fact that det(tA) = tn det(A) for every A ∈ Rn×n and t ∈ Z, the Cramer’s rule yields
|fB (j)| =

n
j+k
−1
· bk · det γA
[AB ]−k,−j
1 X (−1)
det (Qj )

=
−1
det (AB )
γA
det γA
AB
k=1



−1
−1
By the choice of γA , the values det(γA
AB ) and det(γA
[AB ]−k,−j ) are integral for all k, it follows that

|fB (j)| ≤ D||b/γA ||1

and

fB (j) 6= 0
11

=⇒

1
≤ |fB (j)| .
DγA

Lemma 3.2. If every nonzero vector in the kernel of A has positive cost, the minimum energy feasible
solution is kernel-free and unique.
Proof. Let q be a minimum energy feasible solution. Since q is feasible, it can be written as qn + qr , where
qn is a convex combination of basic feasible solutions and qr lies in the kernel of A. Moreover, all elements
in this representation are sign-compatible with q by Lemma 2.1. If qr 6= 0, the vector q − qr is feasible and
has smaller energy, a contradiction. Thus qr = 0.
We next prove uniqueness. Assume for the sake of a contradiction that there are two distinct minimum
energy feasible solutions q (1) and q (2) . We show that the solution (q (1) + q (2) )/2 uses less energy than q (1)
and q (2) . Since h 7→ h2 is a strictly convex function from R to R, the average of the two solutions will be
(1)
(2)
(1)
(2)
better than either solution if there is an index
P e with re > 0 and qe 6= qe . The difference z = q − q
lies in the kernel of A and hence cost(z) = e ce |ze | > 0. Thus there is an e with ce > 0 and ze 6= 0. We
have now shown uniqueness.
Lemma 3.3. Assume that every nonzero vector in the kernel of A has positive cost. Let q be the minimum
energy feasible solution. Then |qe | ≤ D||b/γA ||1 for every e.

Proof. Since q is a convex combination of basic feasible solutions, |qe | ≤ maxz |ze | where z ranges over basic
−1
feasible solutions of the form (zB , 0), where zB = AB
b and AB ∈ Rn×n is a non-singular submatrix of A.
Thus, by Lemma 3.1 every component of z is bounded by D||b/γA ||1 .
In [SV16c], the bound |qe | ≤ D2 m||b||1 was shown. We will now derive explicit formulae for the minimum
energy solution q. We will express q in terms of a vector p ∈ Rn , which we refer to as the potential, by
analogy with the network setting, in which p can be interpreted as the electric potential of the nodes. The
energy of the minimum energy solution is equal to bT p. We also derive a local Lipschitz condition for the
mapping from x to q. Note that for c > 0 these facts are well-known. Let us split the column indices [m] of
A into
P := { e ∈ [m] : ce > 0 } and Z := { e ∈ [m] : ce = 0 } .
(10)

Lemma 3.4. Assume that every nonzero vector in the kernel of A has positive cost. Let re = ce /xe and let
R denote the corresponding diagonal matrix. Let us split A into AP and AZ , and q into qP and qZ . Since AZ
has linearly independent columns, we may assume that the first |Z| rows of AZ form a square non-singular
 A′ A′ 
matrix. We can thus write A = AP′′ AZ′′ with invertible A′Z . Then the minimum energy solution satisfies
P



A′P
A′′P

A′Z
A′′Z



qP
qZ



Z

 ′
b
= ′′
b

and



RP

0



 " ′ T
AP
qP
=
T
qZ
A′Z

T

A′′P
T
A′′Z

#

p′
p′′



(11)

 ′
for some vector p = pp′′ ; here p′ has dimension |Z|. The equation system (11) has a unique solution given
by
   ′ −1 ′

 ′ 

[AZ ] (b − A′P qP )
qZ
p
−[[A′Z ]T ]−1 [A′′Z ]T p′′
=
and
=
,
(12)
qP
p′′
M R−1 M T (b′′ − A′′Z [A′Z ]−1 b′ )
RP−1 ATP p
where M = A′′P − A′′Z [A′Z ]−1 A′P is the Schur complement of the block A′Z of the matrix A.

Proof. q minimizes E(f ) = f T Rf among all solutions of Af = b. The KKT conditions state that q must
satisfy Rq = AT p for some p. Note that 2Rf is the gradient of the energy E(f ) with respect to f and that
the −AT p is the gradient of pT (b − Af ) with respect to f . We may absorb the factor −2 in p. Thus q
satisfies (11).
We show next that the linear system (11) has a unique solution. The top |Z| rows of the left system
in (11) give
qZ = [A′Z ]−1 (b′ − A′P qP ).
(13)
Substituting this expression for qZ into the bottom n − |Z| rows of the left system in (11) yields
M qP = b′′ − A′′Z [A′Z ]−1 b′ .
12

From the top |P | rows of the right system in (11) we infer qP = RP−1 ATP · p. Thus
M RP−1 ATP · p = b′′ − A′′Z [A′Z ]−1 · b′ .

(14)

The bottom n − |Z| rows of the right system in (11) yield 0 = ATZ p = [A′Z ]T p′ + [A′′Z ]T p′′ and hence
p′ = −[[A′Z ]T ]−1 [A′′Z ]T p′′ .

(15)

Substituting (15) into (14) yields
b′′ − A′′Z [A′Z ]−1 b′ = M RP−1 [A′P ]T p′ + [A′′P ]T p′′




= M RP−1 [A′′P ]T − [A′P ]T [[A′Z ]T ]−1 [A′′Z ]T p′′

= M RP−1 M T p′′ .

(16)

It remains to show that the matrix M RP−1 M T is non-singular. We first observe that the rows of M are
linearly independent. Consider the left system in (11). Multiplying the first |Z| rows by [A′Z ]−1 and
then subtracting A′′Z times the resulting rows from the last n − |Z| rows turns A into the matrix Q =
 [A′ ]−1 A′ I 
Z
P
. By assumption, A has independent rows. Moreover, the preceding operations guarantee that
M
0
−1/2

rank(A) = rank(Q). Therefore, M has independent rows. Since RP−1 is a positive diagonal matrix, RP
exists and is a positive diagonal matrix. Let z be an arbitrary nonzero vector of dimension |P |. Then
−1/2
−1/2
z T M RP−1 M T z = (RP M T z)T (RP M T z) > 0 and hence M RP−1 M T is non-singular. It is even positive
semi-definite.
There is a shorter proof that the system (11) has a unique solution. However, the argument does not give
an explicit expression for the solution. In the case of a convex objective function and affine constraints, the
KKT conditions are sufficient for being a global minimum. Thus any solution to (11) is a global optimum.
We have already shown in Lemma 3.2 that the global minimum is unique.
We next observe that the energy of q can be expressed in terms of the potential.
Lemma 3.5. Let q be the minimum energy feasible solution and let f be any feasible solution. Then
E(q) = bT p = f T AT p.
Proof. As in the proof of Lemma 3.4, we split q into qP and qZ , R into RP and RZ , and A into AP and AZ .
Then
E(q) = qPT RP qP

by the definition of E(q) and since RZ = 0

T

by the right system in (11)

T

by the left system in (11)

T

by the right system in (11).

= p AP qp
= p (b − AZ qZ )
=b p

For any feasible solution f , we have f T AT p = bT p.
3.1.2

The Mapping x 7→ q is Locally-Lipschitz

We show that the mapping x 7→ q is Lipschitz continuous; this implies existence of a solution x(t) with domain
[0, ∞) by standard arguments. Our analysis builds upon Cramer’s rule and the Cauchy-Binet formula. The
Cauchy-Binet formula extends Kirchhoff’s spanning tree theorem which was used in [BMV12] for the analysis
of the undirected shortest path problem.
Lemma 3.6. Assume c ≥ 0, no non-zero vector in the kernel of A has cost zero, and that A, b, and
c are integral. Let α, β > 0. For any two vectors x and x̃ in Rm with α1 ≤ x, x̃ ≤ β1, define γ :=
2mn (β/α)n cnmax D2 ||b/γA ||1 . Then |qe (x)| − |qe (x̃)| ≤ γ||x − x̃||∞ for every e ∈ [m].
13

Proof. First assume that c > 0. By Cramer’s rule
(AR−1 AT )−1 =

1
((−1)i+j det(M−j,−i ))ij ,
det(AR−1 AT )

where M−i,−j is obtained from AR−1 AT by deleting the i-th row and the j-th column. For a subset S of [m]
and an index i ∈ [n], let AS be the n × |S| matrix consisting of the columns selected by S and let A−i,S be
the matrix obtained from AS by deleting row i. If D is a diagonal matrix of size m, then (AD)S = AS DS .
The Cauchy-Binet theorem expresses the determinant of a product of two matrices (not necessarily square)
as a sum of determinants of square matrices. It yields
X
(det((AR−1/2 )S ))2
det(AR−1 AT ) =
S⊆[m]; |S|=n

X

=

(

Y

S⊆[m]; |S|=n e∈S

xe /ce ) · (det AS )2 .

Similarly,
X

det(AR−1 AT )−i,−j =

(

Y

S⊆[m]; |S|=n−1 e∈S

Using p = (AR−1 AT )−1 b, we obtain
P
P
i+j
j∈[n] (−1)
PS⊆[m];
pi =

xe /ce ) · (det A−i,S · det A−j,S ).

xe /ce ) · (det A−i,S · det A−j,S )bj
Q e∈S
.
2
|S|=n ( e∈S xe /ce ) · (det AS )

|S|=n−1 (

S⊆[m];

Q

(17)

Substituting into q = R−1 AT p yields
xe T
A p
qe =
ce e
Q
P
P
i+j+2n
xe X
S⊆[m]; |S|=n−1 ( e′ ∈S xe′ /ce′ ) · (det A−i,S · det A−j,S )bj
j∈[n] (−1)
Q
P
Ai,e ·
=
2
ce i
S⊆[m]; |S|=n ( e′ ∈S xe′ /ce′ ) · (det AS )


P
P
Q
P
j+n
i+n
′ /ce′ ) ·
(−1)
b
det
A
(−1)
A
det
A
·
x
(
′
j
−j,S
i,e
−i,S
e
j∈[n]
i∈[n]
e ∈S∪e
S⊆[m]; |S|=n−1
Q
P
=
2
′
′
)
·
(det
A
)
/c
x
(
′
S
e
e ∈S e
S⊆[m]; |S|=n
Q
P
′
′
)
·
det(A
|A
)
·
det(A
|b)
/c
x
(
′
S
e
S
e
e ∈S∪e e
S⊆[m]; |S|=n−1
Q
P
=
,
(18)
′ /ce′ ) · (det AS )2
x
(
′
e
e ∈S
S⊆[m]; |S|=n

where (AS |Ae ), respectively (AS |b), denotes the n × n matrix whose columns are selected from A by S and
whose last column is equal to Ae , respectively b.
We are now ready to estimate the derivative ∂qe /∂xi . Assume first that e 6= i. By the above, qe =
xe F +Gxi /ci
ce H+Ixi /ci , where F , G, H and I are given implicitly by (18). Then
 n 2
m
2 · n−1
β D ||b/γA ||1
xe F I/ci − GH/ci
∂qe
|=|
|≤
≤ γ.
|
2
∂xi
ce (H + Ixi /ci )
(α/cmax )n
For e = i, we have qe =

Gxe /ce
H+Ixe /ce ,

where G, H, and I are given implicitly by (18). Then
 n 2
m
∂qe
GH/ce
n−1 β D ||b/γA ||1
|
|=|
|≤
≤ γ.
∂xe
(H + Ixe /ce )2
(α/cmax )n

Finally, consider x and x̃ with α1 ≤ x, x̃ ≤ β1. Let x̄ℓ = (x̃1 , . . . , x̃ℓ , xℓ+1 , . . . , xm ). Then
X
||qe (x)| − |qe (x̃)|| ≤ |qe (x) − qe (x̃)| ≤
|qe (x̄ℓ ) − qe (x̄ℓ+1 )| ≤ γ||x − x̃||1 .
0≤ℓ<m

14

In the general case where c ≥ 0, we first derive an expression for p′′ similar to (17). Then the equations for
p′ in (12) yield p′ , the equations for qP in (12) yield qP , and finally the equations for qZ in (12) yield qZ .
We are now ready to establish the existence of a solution with domain [0, ∞).
Lemma 3.7. The solution to the undirected dynamics in (1) has domain [0, ∞). Moreover, x(0)e−t ≤
x(t) ≤ D||b/γA ||1 · 1 + max(0, x(0) − D||b/γA ||1 · 1)e−t for all t.
Proof. Consider any x0 > 0 and any t0 ≥ 0. We first show that there is a positive δ ′ (depending on x0 ) such
that a unique solution x(t) with x(t0 ) = x0 exists for t ∈ (t0 − δ ′ , t0 + δ ′ ). By the Picard-Lindelöf Theorem,
this holds true if the mapping x 7→ |q| − x is continuous and satisfies a Lipschitz condition in a neighborhood
of x0 . Continuity clearly holds. Let ε = mini (x0 )i /2 and let U = { x : ||x − x0 ||∞ < ε }. Then for every
x, x̃ ∈ U and every e
|qe (x)| − |qe (x̃)| ≤ γ||x − x̃||1 ,

where γ is as in Lemma 3.6. Local existence implies the existence of a solution which cannot be extended.
Since q is bounded (Lemma 3.3), x is bounded at all finite times, and hence the solution exists for all t.
The lower bound x(t) ≥ x(0)e−t > 0 holds by Fact 2.2 with A = 0 and α = −1. Since |qe | ≤ D||b/γA ||1 ,
ẋ = |q| − x ≤ D||b/γA ||1 · 1 − x, we have x(t) ≤ D||b/γA ||1 · 1 + max(0, x(0) − D||b/γA ||1 · 1)e−t by Fact 2.2
with B = D||b/γA ||1 · 1 and β = −1.

3.2

LP Duality

The energy E(x) is no longer a Lyapunov function, e.g., if x(0) ≈ 0, x(t) and hence E(x(t)) will grow
initially. We will show that energy suitably scaled is a Lyapunov function. What is the appropriate scaling
factor? In the case of the undirected shortest path problem, [BMV12] used the minimum capacity of any
source-sink cut as a scaling factor. The proper generalization to our situation is to consider the linear
program max{α : Af = αb, |f | ≤ x}, where x is a fixed positive vector. Linear programming duality yields
the corresponding minimization problem which generalizes the minimum cut problem to our situation.
Lemma 3.8. Let x ∈ Rm
>0 and b 6= 0. The linear programs
max{α : Af = αb, |f | ≤ x}

and

min{|y T A|x : bT y = −1}

(19)

are feasible and have the same objective value. Moreover, there is a finite set YA = { d1 , . . . , dK } of vectors
T
di ∈ Rm
≥0 that are independent of x such that the minimum above is equal to C⋆ = mind∈YA d x. There is a
6
feasible f with |f | ≤ x/C⋆ .
Proof. The pair (α, f ) = (0, 0) is a feasible solution for the maximization problem. Since b 6= 0, there exists
y with bT y = −1 and thus both problems are feasible. The dual of max{α : Af − αb = 0, f ≤ x, −f ≤ x}
has unconstrained variables y ∈ Rn and non-negative variables z + , z − ∈ Rm and reads
min{xT (z + + z − ) : −bT y = 1, AT y + z + − z − = 0, z + , z − ≥ 0}.

(20)

From z − = AT y + z + , z + ≥ 0, z − ≥ 0 and x > 0, we conclude min(z + , z − ) = 0 in an optimal solution.
Thus z − = max(0, AT y) and z + = max(0, −AT y) and hence z + + z − = |AT y| in an optimal dual solution.
Therefore, (20) and the right LP in (19) have the same objective value.
We next show that the dual attains its minimum at a vertex of the feasible set. For this it suffices to show
that its feasible set contains no line. Assume it does. Then there are vectors d = (y1 , z1+ , z1− ), d non-zero,
and p = (y0 , z0+ , z0− ) such that (y, z + , z − ) = p + λd = (y0 + λy1 , z0+ + λz1+ , z0− + λz1− ) is feasible for all λ ∈ R.
Thus z1+ = z1− = 0. Note that if either z1+ or z1− would be non-zero then either z0+ + λz1+ or z0− + λz1− would
6 In the undirected shortest path problem, the d’s are the incidence vectors of the undirected source-sink cuts. Let S be any
set of vertices containing s0 but not s1 , and let 1S be its associated indicator vector. The cut corresponding to S contains the
edges having exactly one endpoint in S. Its indicator vector is dS = |AT 1S |. Then dS
e = 1 iff |S ∩ { u, v }| = 1, where e = (u, v)
S T
or e = (v, u), and dS
e = 0 otherwise. For a vector x ≥ 0, (d ) x is the capacity of the source-sink cut (S, V \S). In this setting,
C⋆ is the value of a minimum cut.

15

have a negative component for some lambda. Then AT y + z + + z − = 0 implies AT y1 = 0. Since A has full
row rank, y1 = 0. Thus the dual contains no line and the minimum is attained at a vertex of its feasible
region. The feasible region of the dual does not depend on x.
+
−
Let (y 1 , z1+ , z1− ) to (y K , zK
, zK
) be the vertices of (20), and let YA = { |AT y 1 |, . . . , |AT y K | }. Then
min dT x = min{xT (z + + z − ) : −bT y = 1, AT y + z + − z − = 0, z + , z − ≥ 0}

d∈YA

= min{|y T A|x : bT y = −1}.
We finally show that there is a feasible f with |f | ≤ x/C⋆ . Let x′ := x/C⋆ . Then x′ > 0 and
mind∈YA dT x′ = mind∈YA dT x/C⋆ = C⋆ /C⋆ = 1 and thus the right LP with x = x′ (19) has objective
value 1. Hence, the left LP has objective value 1 and there is a feasible f with |f | ≤ x′ .

3.3

Convergence to Dominance

In the network setting, an important role is played by the set of edge capacity vectors that support a feasible
flow. In the LP setting, we generalize this notion to the set of dominating states, which is defined as
Xdom := {x ∈ Rm : ∃ feasible f : |f | ≤ x}.
An alternative characterization, using the set YA from Lemma 3.8, is
T
X1 := {x ∈ Rm
≥0 : d x ≥ 1 for all d ∈ YA }.

We prove that Xdom = X1 and that the set X1 is attracting in the sense that the distance between x(t) and
X1 goes to zero, as t increases.
Lemma 3.9.
1. It holds that Xdom = X1 . Moreover, limt→∞ dist(x(t), X1 ) = 0, where dist(x, X1 ) is the
Euclidean distance between x and X1 .
2. If x(t0 ) ∈ X1 , then x(t) ∈ X1 for all t ≥ t0 . For all sufficiently large t, x(t) ∈ X1/2 := {x ∈ Rn≥0 :
dT x ≥ 1/2 for all d ∈ YA }, and if x ∈ X1/2 then there is a feasible f with |f | ≤ 2x.
Proof.
1. If x ∈ X1 , then dT x ≥ 1 for all d ∈ YA and hence Lemma 3.8 implies the existence of a feasible
solution f with |f | ≤ x. Conversely, if x ∈ Xdom , then there is a feasible f with |f | ≤ x. Thus dT x ≥ 1
for all d ∈ YA and hence x ∈ X1 . By the proof of Lemma 3.8, for any d ∈ YA , there is a y such that
d = |AT y| and bT y = −1. Let Y (t) = dT x. Then
Ẏ = |y T A|ẋ = |y T A|(|q| − x) ≥ |y T Aq| − Y = |y T b| − Y = 1 − Y.
Thus for any t0 and t ≥ t0 , Y (t) ≥ 1 + (Y (t0 ) − 1)e−(t−t0 ) by Lemma 2.2 applied with A = 1
and α = −1. In particular, lim inf t→∞ Y (t) ≥ 1. Thus lim inf t→∞ mind∈YA dT x ≥ 1 and hence
limt→∞ dist(x(t), X1 ) = 0.
2. Moreover, if Y (t0 ) ≥ 1, then Y (t) ≥ 1 for all t ≥ t0 . Hence x(t0 ) ∈ X1 implies x(t) ∈ X1 for all t ≥ t0 .
Since x(t) converges to X1 , x(t) ∈ X1/2 for all sufficiently large t. If x ∈ X1/2 there is f such that
Af = 21 b and |f | ≤ x. Thus 2f is feasible and |2f | ≤ 2x.
The next lemma summarizes simple bounds on the values of resistors r, potentials p and states x that
hold for sufficiently large t. Recall that P = { e ∈ [m] : ce > 0 } and Z = { e ∈ [m] : ce = 0 }, see (10).
Lemma 3.10.
1. For sufficiently large t, it holds that re ≥ ce /(2D||b/γA ||1 ), bT p ≤ 8D||b/γA ||1 ||c||1 and
T
2
|Ae p| ≤ 8D ||b||1 ||c||1 for all e.
2. For all e, it holds that ẋe /xe ≥ −1 and for all e ∈ P , it holds that ẋe /xe ≤ 8D2 ||b||1 ||c||1 /cmin .
3. There is a positive constant C such that for all t ≥ t0 , there is a feasible f (depending on t) such that
xe (t) ≥ C for all indices e in the support of f .
16

Proof.
1. By Lemma 3.7, xe (t) ≤ 2D||b/γA ||1 for all sufficiently large t. It follows that re = ce /xe ≥
ce /(2D||b/γA ||1 ). Due to Lemma 3.9, for large enough t, there is a feasible flow with |f | ≤ 2x. Together
with xe (t) ≤ 2D||b/γA ||1 , it follows that
bT p = E(q) ≤ E(2x) = 4cT x ≤ 8D||b/γA ||1 ||c||1 .
Now, orient A according to q and consider any index e′ . Recall that for all indices e, we have ATe p = 0
if e ∈ Z, and qe = (xe /ce ) · ATe p if e ∈ P . Thus ATe p ≥ 0 for all e. If e′ ∈ Z or e′ ∈ P and qe′ = 0, the
claim is obvious. So assume e′ ∈ P and qe′ > 0. Since q is a convex combination of q-sign-compatible
basic feasible solutions, there is a basic feasible solution f with f ≥ 0 and fe′ > 0. By Lemma 3.1,
fe′ ≥ 1/(DγA ). Therefore
X
fe′ ATe′ p ≤
fe ATe p = bT p ≤ 8D||b/γA ||1 ||c||1
e

for all sufficiently large t. The inequality follows from fe ≥ 0 and ATe p ≥ 0 for all e. Thus ATe′ p ≤
8D2 ||b||1 ||c||1 for all sufficiently large t.
2. We have ẋe /xe = (|qe | − xe )/xe ≥ −1 for all e. For e with ce > 0
|qe | − xe
|AT p|
ẋe
=
≤ e ≤ 8D2 ||b||1 ||c||1 /cmin .
xe
xe
ce
3. Let t0 be such that dT x(t) ≥ 1/2 for all d ∈ YA and t ≥ t0 . Then for all t ≥ t0 , there is f such that
Af = 21 b and |f | ≤ x(t); f may depend on t. By Lemma 2.1, we can write 2f as convex combination of
f -sign-compatible basic feasible solutions (at most m of them) and a f -sign-compatible solution in the
kernel of A. Dropping the solution in the kernel of A leaves us with a solution which is still dominated
by x.
It holds that for every e ∈ E with fe 6= 0, there is a basic feasible solution g used in the convex
decomposition such that 2|fe | ≥ |ge | > 0. By Lemma 3.1, every non-zero component of g is at least
1/(DγA ). We conclude that xe ≥ 1/(2DγA ), for every e in the support of g.

3.4

The Equilibrium Points

We next characterize the equilibrium points
F = { x ∈ R≥0 : |q| = x } .

(21)

Let us first elaborate on the special case of the undirected shortest path problem. Here the equilibria are the
flows of value one from source to sink in a network formed by undirected source-sink paths of the same length.
This can be seen as follows. Consider any x ≥ 0 and assume supp(x) is a network of undirected source-sink
paths of the same length. Call this network N . Assign to each node u, a potential pu equal to the length of the
shortest undirected path from the sink s1 to u. These potentials are well-defined as all paths from s1 to u in
N must have the same length. For an edge e = (u, v) in N , we have qe = xe /ce (pu − pv ) = xe /ce ·ce = xe , i.e.,
q = x is the electrical flow with respect to the resistances ce /xe . Conversely, if x is an equilibrium point and
the network is oriented such that q ≥ 0, we have xe = qe = xe /ce (pu − pv ) for all edges e = (u, v) ∈ supp(x).
Thus ce = pu − pv and this is only possible if for every node u, all paths from u to the sink have the same
length. Thus supp(x) must be a network of undirected source-sink paths of the same length. We next
generalize this reasoning.
Theorem 3.11. If x = |q| is an equilibrium point and the columns of A are oriented such that q ≥ 0, then
all feasible solutions f with supp(f ) ⊆ supp(x) satisfy cT f = cT x. Conversely, if x = |q| for a feasible q, A
is oriented such that q ≥ 0, and all feasible solutions f with supp(f ) ⊆ supp(x) satisfy cT f = cT x, then x is
an equilibrium point.
17

Proof. If x is an equilibrium point, |qe | = xe for every e. By changing the signs of some columns of A, we
may assume q ≥ 0, i.e., q = x. Let p be the potential with respect to x. For every index e ∈ P in the support
of x, since ce > 0 we have qe = xcee ATe p and hence ce = ATe p. Further, for the indices e ∈ Z in the support of
x, we have ce = 0 = ATe p due to the second block of equations on the right hand side in (11). Let f be any
feasible solution whose support is contained in the support of x. Then the first part follows by
X
X
fe ATe p = bT p = E(q) = E(x) = cost(x).
ce f e =
e∈supp(f )

e∈supp(f )

For the second part, we misuse notation and use A to also denote the submatrix of the constraint matrix
indexed by the columns in the support of x. We may assume that the rows of A are independent. Otherwise,
we simply drop redundant constraints. We may assume q ≥ 0; otherwise we simply change the sign of some
columns of A. Then x is feasible. Let AB be a square non-singular submatrix of A and let AN consist of the
remaining columns of A. The feasible solutions f with supp(f ) ⊆ supp(x) satisfy AB fB + AN fN = b and
hence fB = A−1
B (b − AN fN ). Then
T
T −1
cT f = cTB fB + cTN fN = cB A−1
B b + (cN − cB AB AN )fN .

Since, by assumption, cT f is constant for all feasible solutions whose support is contained in the support of
 cB 
 ATB  −1 T
−1 T
T
T
x, we must have cN = ATN [A−1
B ] cB . Let p = [AB ] cB . Then it follows that A p = AT [AB ] cB = cN
N

and hence Rx = AT p. Thus the pair (x, p) satisfies the right hand side of (11). Since x is feasible, it also
satisfies the left hand side of (11). Therefore, x is the minimum energy solution with respect to x.
Corollary 3.12. Let g be a basic feasible solution. Then |g| is an equilibrium point.

Proof. Let g be a basic feasible solution. Orient A such that g ≥ 0. Since g is basic, there is a B ⊆ [m] such
that g = (gB , gN ) = (A−1
B b, 0). Consider any feasible solution f with supp(f ) ⊆ supp(g). Then f = (fB , 0)
and hence b = Af = AB fB . Therefore, fB = gB and hence cT f = cT g. Thus x = |g| is an equilibrium
point.
This characterization of equilibria has an interesting consequence.
Lemma 3.13. The set L := {cT x : x ∈ F } of costs of equilibria is finite.

Proof. If x is an equilibrium, x = |q|, where q is the minimum energy solution with respect to x. Orient A
such that q ≥ 0. Then by Theorem 3.11, cT f = cT x for all feasible solutions f with supp(f ) ⊆ supp(x). In
particular, this holds true for all such basic feasible solutions f . Thus L is a subset of the set of costs of all
basic feasible solutions, which is a finite set.
We conclude this part by showing that the optimal solutions of the undirected linear program (2) are
equilibria.
Theorem 3.14. Let x be an optimal solution to (2). Then x is an equilibrium.
Proof. By definition, there is a feasible f with |f | = x. Let us reorient the columns of A such that f ≥ 0
and let us delete all columns e of A with fe = 0. Consider any feasible g with supp(g) ⊆ supp(x). We claim
that cT x = cT g. Assume otherwise and consider the point y = x + λ(g − x). If |λ| is sufficiently small, y ≥ 0.
Furthermore, y is feasible and cT y = cT x + λ(cT g − cT x). If cT g 6= cT x, x is not an optimal solution to (2).
The claim now follows from Theorem 3.11.

3.5

Convergence

In order to show convergence, we construct a Lyapunov function. The following functions play a crucial
role in our analysis. Let Cd = dT x for d ∈ YA , and recall that C⋆ = mind∈YA dT x denotes the optimum.
Moreover, we define by
 
X
cT x
x
xe
and Vd :=
−E
for every d ∈ YA .
h(t) :=
re |qe |
C⋆
C⋆
Cd
e
18

Theorem 3.15.
(1) For every d ∈ YA , Ċd ≥ 1 − Cd . Thus, if Cd < 1 then C˙d > 0.
d
(2) If x(t) ∈ X1 , then dt
cost(x(t)) ≤ 0 with equality if and only if x = |q|.
(3) Let d ∈ YA be such that C⋆ = dT x at time t. Then it holds that V̇d ≤ h(t).
(4) It holds that h(t) ≤ 0 with equality if and only if |q| = Cx⋆ .
Proof.
1. Recall that for d ∈ YA , there is a y such that bT y = −1 and d = |AT y|. Thus C˙d = dT (|q| − x) ≥
T
|y Aq| − Cd = 1 − Cd and hence C˙d > 0, whenever Cd < 1.
2. Remember that E(x) = cost(x) and that x(t) ∈ X1 implies that there is a feasible f with |f | = x. Thus
E(q) ≤ E(f ) ≤ E(x). Let R be the diagonal matrix of entries ce /xe . Then
d
cost(x) = cT (|q| − x)
dt
= xT R1/2 R1/2 |q| − xT Rx
T

1/2

≤ (q Rq)

T

1/2

(x Rx)

≤0

by (1)
since c = Rx
T

by Cauchy-Schwarz

− x Rx

since E(q) ≤ E(x).

If the derivative is zero, both inequalities above have to be equalities. This is only possible
P if the
vectors |q| and x are parallel and E(q) = E(x). Let λ be such that |q| = λx. Then E(q) = e xcee qe2 =
P
λ2 e ce xe = λ2 E(x). Since E(x) > 0, this implies λ = 1.
d
cost(x) =
3. By definition of d, C⋆ = Cd . By the first two items, we have C˙⋆ = dT |q| − C⋆ and dt
T
c |q| − cost(x). Thus
d
cost(x) − C˙⋆ cost(x)
C⋆ dt
d cost(x)
C⋆ (cT |q| − cost(x)) − (dT |q| − C⋆ )cost(x)
=
=
dt C⋆
C⋆2
C⋆2
T
T
T
xe 2
xe X
C⋆ · c |q| − d |q| · c x X
= h(t),
≤
−
re |qe |
re
=
2
C⋆
C⋆
C⋆
e
e

P
where we used re = ce /xe and hence cT |q| = e re xe |qe |, cT x = E(x), and dT |q| ≥ |y T Aq| = 1 since
d = |y T A| for some y with bT y = −1.
4. We have
1/2
1/2  X
X
X
X
1/2
xe
1/2 xe 1/2
2
xe 2
re C
= E Cx⋆
E(q)1/2
|q
|
=
r
|q
|
≤
)
r
r
q
r
(
e
e
e
e
e C⋆ e
e
C⋆
⋆
e

e

by Cauchy-Schwarz. Since h(t) =
h(t) ≤ E


x 1/2
C⋆

e

e

P

xe
e re |qe | C⋆

· E(q)1/2 − E

− E( Cx⋆ ) by definition, it follows that

x
C⋆



=E


x 1/2
C⋆


· E(q)1/2 − E


x 1/2
C⋆



≤0

since x/C⋆ dominates a feasible solution and hence E(q) ≤ E(x/C⋆ ). If h(t) = 0, we must have equality
in the application of Cauchy-Schwarz, i.e., the vectors x/C⋆ and |q| must be parallel, and we must have
E(q) = E(x/C⋆ ) as in the proof of part 2.
We show now convergence against the set of equilibrium points. We need the following technical Lemma
from [BMV12].
Lemma 3.16 (Lemma 9 in [BMV12]). Let f (t) = maxd∈YA fd (t), where each fd is continuous and differentiable. If f˙(t) exists, then there is a d ∈ YA such that f (t) = fd (t) and f˙(t) = f˙d (t).
Theorem 3.17. All trajectories converge to the set F of equilibrium points.
Proof. We distinguish cases according to whether the trajectory ever enters X1 or not. If the trajectory
d
cost(x) ≤ 0 for all t ≥ t0 with equality only if x = |q|. Thus the trajectory
enters X1 , say x(t0 ) ∈ X1 , then dt
converges to the set of fix points. If the trajectory never enters X1 , consider V = maxd∈YA (Vd + 1 − Cd ).
19

We show that V̇ exists for almost all t. Moreover, if V̇ (t) exists, then V̇ (t) ≤ 0 with equality if and only if
|qe | = xe for all e. It holds that V is Lipschitz-continuous as the maximum of a finite number of continuously
differentiable functions. Since V is Lipschitz-continuous, the set of t’s where V̇ (t) does not exist has zero
Lebesgue measure (see for example [CLSW98, Ch. 3]). If V̇ (t) exists, we have V̇ (t) = V̇d (t) − Ċd (t) for some
d ∈ YA according to Lemma 3.16. Then, it holds that V̇ (t) ≤ h(t) − (1 − Cd ) ≤ 0. Thus x(t) converges to
the set

x ∈ R≥0 : V̇ = 0 = { x ∈ R≥0 : |q| = x/C and C = 1 } = { x ∈ R≥0 : |q| = x } .

At this point, we know that all trajectories x(t) converge to F . Our next goal is to show that cT x(t)
converges to the cost of an optimum solution of (2) and that |q| − x converges to zero. We are only able to
show the latter for all indices e ∈ P , i.e. with ce > 0.

3.6

Details of the Convergence Process

In the argument to follow, we will encounter
the following situation several times. We have a non-negative
R∞
function f (t) ≥ 0 and we know that 0 f (t)dt is finite. We want to conclude that f (t) converges to zero
for t → ∞. This holds true if f is Lipschitz continuous. Note that the proof of the following lemma is very
similar to the proof in [BMV12, Lemma 11]. However, in our case we apply the Local Lipschitz condition
that we showed in Lemma 3.6.
R∞
Lemma 3.18. Let f (t) ≥ 0 for all t. If 0 f (t)d(t) is finite and f (t) is Lipschitz-continuous, i.e., for every
ε > 0, there is a δ > 0 such that |f (t′ ) − f (t)| ≤ ε for all t′ ∈ [t, t + δ], then f (t) converges to zero as t goes
to infinity. The functions t 7→ xT R|q| − xT Rx = cT |q| − cT x and t 7→ h(t) are Lipschitz continuous.
Proof. If f (t) does not converge to zero, there is ε > 0 and an infinite unbounded sequence t1 , t2 , . . . such
that f (ti ) ≥ ε for all i. Since fR is Lipschitz continuous there is δ > 0 such that f (t′i ) ≥ ε/2 for t′i ∈ [ti , ti + δ]
∞
and all i. Hence, the integral 0 f (t)dt is unbounded.
Since x˙e is continuous and bounded (by Lemma 3.7), xe is Lipschitz-continuous. Thus, it is enough to
show that qe is Lipschitz-continuous for all e. Since qZ (recall that Z = { e : ce = 0 } and P = [m] \ Z) is an
affine function of qP , it suffices to establish the claim for e ∈ P . So let e ∈ P be such that ce > 0. First, we
claim that xe (t + ε) ≤ (1 + 2Kε)xe for all ε ≤ K/4, where K = 8D2 ||b||1 ||c||1 /cmin. Assume that this is not
the case. Let
ε = inf{δ ≤ 1/4K : xe (t + δ) > (1 + 2Kδ)xe (t)},
then ε > 0 (since x˙e (t) ≤ Kxe (t) by Lemma 3.10) and, by continuity, xe (t + ε) ≥ (1 + 2Kε)xe (t). There
must be t′ ∈ [t, t + ε] such that x˙e (t′ ) = 2Kxe (t). On the other hand,
x˙e (t′ ) ≤ Kxe (t′ ) ≤ K(1 + 2Kε)xe (t) ≤ K(1 + 2K/4K)xe (t) < 2Kxe (t),
which is a contradiction. Thus, xe (t + ε) ≤ (1 + 2Kε)xe for all ε ≤ 1/4K. Similarly, xe (t + ε) ≥ (1 − 2Kε)xe .
Now, let α = (1 − 2Kε)xe and β = (1 + 2Kε)xe . Then
||qe (t + δ)| − |qe (t)|| ≤ M ||x(t + δ) − x(t)||1 ≤ M m(4Kε)xe ≤ 8εM mKD||b/γA||1 ,
since xe ≤ 2D||b/γA ||1 for sufficiently large t and where M is as in Lemma 3.6. Since C⋆ is at least 1/2 for
all sufficiently large t, the division by C⋆ and C⋆2 in the definition of h(t) does not affect the claim.
Lemma 3.19. For all e ∈ [m] of positive cost, it holds that |xe − |qe || → 0 as t goes to infinity.
d
cost(x) ≤ xT R|q| − xT Rx ≤ 0 with equality if
Proof. For a trajectory ultimately running in X1 , we showed dt
and only if x = |q|. Also, E(q) ≤ E(x), since x dominates a feasible solution. Furthermore, xT R|q| − xT Rx
goes to zero using Lemma 3.18. Thus
X
X
X
X
X
X

re (xe − |qe |)2 =
re x2e +
re qe2 − 2
re xe |qe | ≤ 2
re x2e −
re xe |qe |
e

e

e

e

20

e

e

goes to zero. Next observe that there is a constant C such
P that xe (t) ≤ C Cfor all
P e and t as a result of
Lemma 3.7. Also cmin > 0 and hence re ≥ cmin /C. Thus e re (xe − |qe |)2 ≤ cmin
· e (xe − |qe |)2 and hence
|xe − |qe || → 0 for every e with positive cost. For trajectories outside X1 , we argue about ||qe | − Cx⋆ | and use
C⋆ → 1, namely
X
X
X

xe 2
xe
)
−
|qe | → 0.
re ( C
re C
re ( Cxe⋆ − |qe |)2 ≤ 2
⋆
⋆
e

e

e

Note that the above does not say anything about the indices e ∈ Z (with ce = 0). Recall that AP qP +
AZ qZ = b and that the columns of AZ are independent. Thus, qZ is uniquely determined by qP . For the
undirected shortest path problem, the potential difference pT b between source and sink converges to the
length of a shortest source-sink path. If an edge with positive cost is used by some shortest undirected path,
then no shortest undirected path uses it with the opposite direction. We prove the natural generalizations.
Let OP T be the set of optimal solutions to (2) and let Eopt = ∪x∈OP T supp(x) be the set of columns
used in some optimal solution. The columns of positive cost in Eopt can be consistently oriented as the
following Lemma shows.
Lemma 3.20. Let x∗1 and x∗2 be optimal solutions to (2) and let f and g be feasible solutions with |f | = x∗1
and |g| = x∗2 . Then there is no e such that fe ge < 0 and ce > 0.
Proof. Assume otherwise. Then |ge − fe | = |ge | + |fe | > 0. Consider h = (ge f − fe g)/(ge − fe ). Then
−fe ge
= 0 and for every index e′ , it holds
Ah = (ge Af − fe Ag)/(ge − fe ) = b and h is feasible. Also, he = ge gfee −f
e
that |he′ | =

|ge fe′ −fe ge′ |
|ge −fe )|

≤

|ge ||fe′ |+|fe ||ge′ |
|ge |+|fe |

and hence

cost(h) < cost(f ) + cost(g) =

|ge |
|fe |
cost(x∗1 ) +
cost(x∗2 ) = cost(x∗1 ),
|ge | + |fe |
|ge | + |fe |

a contradiction to the optimality of x∗1 and x∗2 .
By the preceding Lemma, we can orient A such that fe ≥ 0 whenever |f | is an optimal solution to (2)
and ce > 0. We then call A positively oriented.
Lemma 3.21. It holds that pT b converges to the cost of an optimum solution of (2). If A is positively
oriented, then lim inf t→∞ ATe p ≥ 0 for all e.
Proof. Let x∗ be an optimal solution of (2). We first show convergence to a point in L and then convergence
to cT x∗ . Let ε > 0 be arbitrary. Consider any time t ≥ t0 , where t0 and C as in Lemma 3.10 and moreover
for every e ∈ P . Then xe ≥ C for all indices e in the support of some basic feasible solution
||qe | − xe | ≤ cCε
max
f . For every e ∈ P , we have qe = xcee ATe p. We also assume q ≥ 0 by possibly reorienting columns of A. Hence
|ce − ATe p| = 1 −

qe
xe − qe
cmax
ce =
ce ≤
|qe − xe | ≤ ε.
xe
xe
C

For indices e ∈ Z, we have ATe p = 0 = ce . Since, ||f ||∞ ≤ D||b/γA ||1 (Lemma 3.1), we conclude
X
X
(ce − pT Ae )fe ≤ ε
|fe | ≤ ε · mD||b/γA ||1 .
cT f − p T b =
e

e∈supp(f )

Since the set L is finite, we can let ε > 0 be smaller than half the minimal distance between elements in
L. By the preceding paragraph, there is for all sufficiently large t, a basic feasible solution f such that
|cT f − bT p| ≤ ε. Since bT p is a continuous function of time, cT f must become constant. We have now shown
that bT p converges to an element in L. We P
will next show that bT p converges to the optimum cost. Let x∗
be an optimum solution to (2) and let W = e x∗e ce ln xe . Since x(t) is bounded, W is bounded. We assume

21

that A is positively oriented, thus there is a feasible f ∗ with |f ∗ | = x∗ and fe∗ ≥ 0 whenever ce > 0. By
reorienting zero cost columns, we may assume fe∗ ≥ 0 for all e. Then Ax∗ = b. We have
Ẇ =

X

x∗e ce

e

=

X

e; ce >0

=

X
e

=

X
e

|qe | − xe
xe

x∗e |ATe p| − cost(x∗ )

since qe =

x∗e |ATe p| − cost(x∗ )

xe T
ce Ae p

whenever ce > 0

since ATe p = 0 whenever ce = 0

x∗e (|ATe p| − ATe p) + bT p − cost(x∗ )

and hence bT p − cost(x∗ ) must converge to zero; note that bT p is Lipschitz continuous in t.
Similarly, |ATe p| − ATe p must converge to zero whenever x∗e > 0. This implies lim inf ATe p ≥ 0. Assume
otherwise, i.e., for every ε > 0, we have ATe p < −ε for arbitrarily large t. Since p is Lipschitz-continuous in
t, there is a δ > 0 such that ATe p < −ε/2 for infinitely many disjoint intervals of length δ. In these intervals,
|ATe p| − ATe p ≥ ε and hence W must grow beyond any bound, a contradiction.
Corollary 3.22. E(x) and cost(x) converge to cT x∗ , whereas x and |q| converge to OP T . If the optimum
solution is unique, x and |q| converge to it. Moreover, if e ∈
/ Eopt , xe and |qe | converge to zero.
Proof. The first part follows from E(x) = cost(x) = bT p and the preceding Lemma. Thus x and q converge
to the set F of equilibrium points, see (21), that are optimum solutions to (2). Since every optimum solution
is an equilibrium point by Theorem 3.14, x and q converge to OP T . For e 6∈ Eopt , fe = 0 for every
f ∈ F ∩ OP T . Since x and |q| converge to F ∩ OP T , xe and |qe | converge to zero for every e ∈ Eopt .

4

Improved Convergence Results: Discrete Directed Dynamics

In this section, we present in its full generality our main technical result on the Physarum dynamics (6).

4.1

Overview

Inspired by the max-flow min-cut theorem, we consider the following primal-dual pair of linear programs:
the primal LP is given by max { t : Af = t · b; 0 ≤ f ≤ x } in variables f ∈ Rm and t ∈ R, and its dual
LP reads min xT z : z ≥ 0; z ≥ AT y; bT y = 1 in variables z ∈ Rm and y ∈ Rn . Since the dual feasible
region does not contain a line and the minimum is bounded, the optimum is attained at a vertex, and in an
optimum solution we have z = max{0, AT y}. Let V be the set of vertices of the dual feasible region, and
let Y := { y : (z, y) ∈ V } be the set of their projections on y-space. Then, the dual optimum is given by
min{ max{0, y T A} · x : y ∈ Y }. The set of strongly dominating capacity vectors x is defined as

T
X := x ∈ Rm
.7
>0 : y Ax > 0 for all y ∈ Y

Note that X contains the set of all scaled feasible solutions {x = tf : Af = b, f ≥ 0, t > 0}.
We next discuss the choice of step size. For y ∈ Y and capacity vector x, let α(y, x) := y T Ax. Further,
let α(x) := min { α(y, x) : y ∈ Y } and α(ℓ) := α(x(ℓ) ). Then, for any x ∈ X there is a feasible f such that
0 ≤ f ≤ x/α(x), see Lemma 4.8. In particular, if x is feasible then α(x) = 1, since α(y, x) = 1 for all y ∈ Y .
We partition the Physarum dynamics (6) into the following five regimes and define for each regime a fixed
step size, see Subsection 4.3.
7 In the shortest path problem (recall that b = e − e ) the set Y consists of all y ∈ { −1, +1 }n such that y = 1 = −y ,
n
n
1
P1
x −
i.e., y encodes a cut with S = { i : yi = −1 } and S = { i : yi = +1 }. The condition y T Ax > 0 translates into
a∈E(S,S) a

P

a∈E(S,S)

xa > 0, i.e., every source-sink cut must have positive directed capacity.

22

Corollary 4.1. The Physarum dynamics (6) initialized with x(0) ∈ X and a step size h satisfies:
• If α(0) = 1, we work with h ≤ h0 and have α(ℓ) = 1 for all ℓ.
• If 1/2 ≤ α(0) < 1, we work with h ≤ h0 /2 and have 1 − δ ≤ α(ℓ) < 1 for ℓ ≥ h−1 log(1/2δ) and δ > 0.
• If 1 < α(0) ≤ 1/h0 , we work with h ≤ h0 and have 1 < α(ℓ) ≤ 1 + δ for ℓ ≥ h−1 · log(1/δh0 ) and δ > 0.
• If 0 < α(0) < 1/2, we work with h ≤ α(0) h0 and have 1/2 ≤ α(ℓ) < 1 for ℓ ≥ 1/h.
• If 1/h0 < α(0) , we work with h ≤ 1/4 and have 1 < α(ℓ) ≤ 1/h0 for ℓ = ⌊log1/(1−h) h0 (α(0) −1)/(1−h0)⌋.
In each regime, we have 1 − α(ℓ+1) = (1 − h)(1 − α(ℓ) ).
We give now the full version of Theorem 1.4 which applies for any strongly dominating starting point.
Theorem 4.2. Suppose A ∈ Zn×m has full row rank (n ≤ m), b ∈ Zn , c ∈ Zm
>0 and ε ∈ (0, 1). Given
x(0) ∈ X and its corresponding α(0) , the Physarum dynamics (6) initialized with x(0) runs in two regimes:
(i) The first regime is executed when α(0) 6∈ [1/2, 1/h0] and it computes a point x(t) ∈ X such that
α(t) ∈ [1/2, 1/h0]. In particular, if α(0) < 1/2 then h ≤ (Φ/opt) · (α(0) h0 )2 and t = 1/h. Otherwise,
if α(0) > 1/h0 then h ≤ Φ/opt and t = ⌊log1/(1−h) [h0 (α(0) − 1)/(1 − h0 )]⌋.
(ii) The second regime starts from a point x(t) ∈ X with α(t) ∈ [1/2, 1/h0], it has a step size h ≤
(0)
(Φ/opt) · h20 /2 and outputs for any k ≥ 4C1 /(hΦ) · ln(C2 Ψ(0) /(ε · min{1, xmin })) a vector x(t+k) ∈ X
(t+k)
such that dist(x
, X⋆ ) < ε/(DγA ).
We stated the bounds on h in terms of the unknown quantities Φ and opt. However, Φ/opt ≥ 1/C3 by
Lemma 3.1 and hence replacing Φ/opt by 1/C3 yields constructive bounds for h.
Organization: This section is devoted to proving Theorem 4.2, and it is organized as follows: Subsection 4.2
establishes core efficiency bounds that extend [SV16c] and yield a scale-invariant determinant dependence of
the step size and are applicable to strongly dominating points. Subsection 4.3 gives the definition of strongly
dominating points and shows that the Physarum dynamics (6) initialized with such a point is well defined.
Subsection 4.4 extends the analysis in [BBD+ 13, SV16b, SV16c] to positive linear programs, by generalizing
the concept of non-negative flows to non-negative feasible kernel-free vectors. Subsection 4.5 shows that x(ℓ)
converges to X⋆ for large enough ℓ. Subsection 4.6 concludes the proof of Theorem 4.2.

4.2

Useful Lemmas
def

Recall that R(ℓ) = diag(c) · (X (ℓ) )−1 is a positive diagonal matrix and L(ℓ) = A(R(ℓ) )−1 AT is invertible. Let
p(ℓ) be the unique solution of L(ℓ) p(ℓ) = b. We improve the dependence on DS in [SV16c, Lemma 5.2] to D.
Lemma 4.3. [SV16c, extension of Lemma 5.2] Suppose x(ℓ) > 0, R(ℓ) is a positive diagonal matrix and
(ℓ)
L = A(R(ℓ) )−1 AT . Then for every e ∈ [m], it holds that kAT (L(ℓ) )−1 Ae k∞ ≤ D · ce /xe .
Proof. The statement follows by combining the proof in [SV16c, Lemma 5.2] with Lemma 3.1.
We show next that [SV16b, Corollary 5.3] holds for x-capacitated vectors, which extends the class of
feasible starting points, and further yields a bound in terms of D.
Lemma 4.4. [SV16b, extension of Corollary 5.3] Let p(ℓ) be the unique solution of L(ℓ) p(ℓ) = b and assume
x(ℓ) is a positive vector with corresponding positive scalar α(ℓ) such that there is a vector f satisfying Af =
α(ℓ) · b and 0 ≤ f ≤ x(ℓ) . Then kAT p(ℓ) k∞ ≤ D||c||1 /α(ℓ) .

23

Proof. By assumption, f satisfies α(ℓ) b = Af =

P

e

fe Ae and 0 ≤ f ≤ x(ℓ) . This yields

α(ℓ) kAT p(ℓ) k∞ = kAT (L(ℓ) )−1 · α(ℓ) bk∞ = k
≤

X
e

fe kAT (L(ℓ) )−1 Ae k∞

X
e

(Lem.

fe AT (L(ℓ) )−1 Ae k∞

≤

4.3)

D

X

fe

e

ce
(ℓ)

xe

≤ Dkck1 .

We note that applying Lemma 4.3 and Lemma 4.4 into the analysis of [SV16c, Theorem 1.3] yields
an improved result that depends on the scale-invariant determinant D. Moreover, we show in the next
Subsection 4.3 that the Physarum dynamics (6) can be initialized with any strongly dominating point.
We establish now an upper bound on q that does not depend on x. We then use this upper bound on q
to establish a uniform upper bound on x.
Lemma 4.5. For any x(ℓ) > 0, kq (ℓ) k∞ ≤ mD2 kb/γA k1 .
(ℓ)

(ℓ)

Proof. Let f be a basic feasible solution of Af = b. By definition, qe = (xe /ce )ATe (L(ℓ) )−1 b and thus
qe(ℓ) =

(ℓ)
(ℓ)
xe X
xe X T (ℓ) −1
Ae (L ) Au fu ≤
|fu | · ATe (L(ℓ) )−1 Au ≤ Dkf k1 ,
ce u
ce u

where the last inequality follows by
ATe (L(ℓ) )−1 Au = ATu (L(ℓ) )−1 Ae ≤ kAT (L(ℓ) )−1 Ae k∞

(Lem.

≤

4.3)

D · ce /x(ℓ)
e .

(ℓ)

By Cramer’s rule and Lemma 3.1, we have |qe | ≤ Dkf k1 ≤ mD2 kb/γA k1 .
Let k, t ∈ N. We denote by
q (t,k) =

t+k−1
X
i=t

h (1 − h)t+k−1−i (i)
q
1 − (1 − h)k

and p(t,k) =

t+k−1
X

p(i) .

(22)

i=t

Straightforward checking shows that Aq (t,k) = b. Further, for C := diag(c), t ≥ 0 and k ≥ 1, we have
x(t)

t+k−1
Y
i=t

[1 + h(C −1 AT p(i) − 1)] = x(t+k) = (1 − h)k x(t) + [1 − (1 − h)k ]q (t,k) .

We give next an upper bound on x(k) that is independent of k.
Lemma 4.6. Let Ψ(0) = max{mD2 kb/γAk1 , kx(0) k∞ }. Then kx(k) k∞ ≤ Ψ(0) , ∀k ∈ N.
Proof. We prove the statement by induction. The base case kx(0) k∞ ≤ Ψ(0) is clear. Suppose the statement
holds for some k > 0. Then, triangle inequality and Lemma 4.5 yield
kx(k+1) k∞ ≤ (1 − h)kx(k) k∞ + hkq (k) k∞ ≤ (1 − h)Ψ(0) + hΨ(0) ≤ Ψ(0) .
We show now convergence to feasibility.
Lemma 4.7. Let r(k) = b − Ax(k) . Then r(k+1) = (1 − h)r(k) and hence r(k) = (1 − h)k (b − Ax(0) ).
Proof. By definition x(k+1) = (1 − h)x(k) + hq (k) , and thus the statement follows by
r(k+1) = b − Ax(k+1) = b − (1 − h)Ax(k) − hb = (1 − h)r(k) .
24

4.3

Strongly Dominating Capacity Vectors

For the shortest path problem, it is known that one can start from any capacity vector x for which the directed
capacity of every source-sink cut is positive, where the directed capacity of a cut is the total capacity of the
edges crossing the cut in source-sink direction minus the total capacity of the edges crossing the cut in the
sink-source direction. We generalize this result. We start with the max-flow like LP
max { t : Af = t · b; 0 ≤ f ≤ x }

(23)

in variables f ∈ Rm and t ∈ R and its dual

min xT z : z ≥ 0; z ≥ AT y; bT y = 1

(24)

in variables z ∈ Rm and y ∈ Rn . The feasible region of the dual contains no line. Assume otherwise; say it
contains (z, y) = (z (0) , y (0) ) + λ(z (1) , y (1) ) for all λ ∈ R. Then, z ≥ 0 implies z (1) = 0 and further z ≥ AT y
implies z (0) ≥ AT y (0) + λAT y (1) and hence AT y (1) = 0. Since A has full row rank, we have y (1) = 0. The
optimum of the dual is therefore attained at a vertex. In an optimum solution, we have z = max{0, AT y}.
Let V be the set of vertices of the feasible region of the dual (24), and let
Y := { y : (z, y) ∈ V }
be the set of their projections on y-space. Then, the optimum of the dual (24) is given by

min max{0, y T A} · x .
y∈Y

The set of strongly dominating capacity vectors x is defined by

T
X := x ∈ Rm
>0 : y Ax > 0 for all y ∈ Y

.

(25)

(26)

We next show that for all x(0) ∈ X and sufficiently small step size, the sequence {x(k) }k∈N stays in X.
Moreover, y T Ax(k) converges to 1 for every y ∈ Y . We define by
α(y, x) := y T Ax

and

α(x) := min { α(y, x) : y ∈ Y } .

Let α(ℓ) := α(x(ℓ) ). Then, x(ℓ) ∈ X iff α(ℓ) > 0. We summarize the discussion in the following Lemma.
Lemma 4.8. Suppose x(ℓ) ∈ X. Then, there is a vector f such that Af = α(ℓ) · b and 0 ≤ f ≤ x(ℓ) .

Proof. By the strong duality theorem applied on (23) and (24), it holds by (25) that
o
n
t = min max{0, y T A} · x(ℓ) ≥ min y T Ax(ℓ) = α(ℓ) .
y∈Y

y∈Y

The statement follows by the definition of (23).

We demonstrate now that α(ℓ) converges to 1.
Lemma 4.9. Assume x(ℓ) ∈ X. Then, for any h(ℓ) ≤ min{1/4, α(ℓ)h0 } we have x(ℓ+1) ∈ X and
1 − α(ℓ+1) = (1 − h(ℓ) ) · (1 − α(ℓ) ).
Proof. By applying Lemma 4.4 and Lemma 4.8 with x(ℓ) ∈ X, we have kAT p(ℓ) k∞ ≤ D||c||1 /α(ℓ) and hence
(ℓ)
(ℓ)
(ℓ) T (ℓ)
≥ −(h(ℓ) xe )/(2α(ℓ) h0 ) ≥ −xe /2. Thus,
for every index e it holds −h(ℓ) · c−1
e xe Ae p
(ℓ)
x(ℓ+1)
= (1 − h(ℓ) )x(ℓ)
· [Re(ℓ) ]−1 ATe p(ℓ) ≥
e
e +h

3 (ℓ) 1 (ℓ)
1
x − xe = x(ℓ)
> 0.
4 e
2
4 e

Let y ∈ Y be arbitrary. Then y T b = 1 and hence y T r(ℓ) = y T (b − Ax(ℓ) ) = 1 − y T Ax(ℓ) = 1 − α(y, x(ℓ) ). The
second claim now follows from Lemma 4.7.
We note that the convergence speed crucially depends on the initial point x(0) ∈ X, and in particular to
its corresponding value α(0) . Further, this dependence naturally partitions the Physarum dynamics (6) into
the five regimes given in Corollary 4.1.
25

4.4

x(k) is Close to a Non-Negative Kernel-Free Vector

In this subsection, we generalize [SV16b, Lemma 5.4] to positive linear programs. We achieve this in two
steps. First, we generalize a result by Ito et al. [IJNT11, Lemma 2] to positive linear programs and then we
substitute the notion of a non-negative cycle-free flow with a non-negative feasible
kernel-free vector.

Throughout this and the consecutive subsection, we denote by ρA := max DγA , nD2 ||A||∞ .

n
Lemma 4.10. Suppose a matrix A ∈ Zn×m has full row rank and
Pvector b ∈ Z . Let g be a feasible solution
to Ag = b and S ⊆ [n] be a subset of row indices of A such that i∈S |gi | < 1/ρA . Then, there is a feasible
solution f such that gi · fi ≥ 0 for all i ∈ [n], fi = 0 for all i ∈ S and kf − gk∞ < 1/(DγA ).

Proof. W.l.o.g. we can assume that g ≥ 0 as we could change the signs of the columns of A accordingly. Let
1S be the indicator vector of S. We consider the linear program
min{1TS x : Ax = b, x ≥ 0}
and let opt be its optimum value. Notice that 0 ≤ opt ≤ 1TS g < 1/ρA . Since the feasible region does
not contain a line and the minimum is bounded, the optimum is attained at a basic feasible solution, say
f . Suppose that there is an index i ∈ S with fi > 0. By Lemma 3.1, we have fi ≥ 1/(DγA ). This is a
contradiction to the optimality of f and hence fi = 0 for all i ∈ S.
Among the feasible solutions f such that fi gi ≥ 0 for all i and fi = 0 for all i ∈ S, we choose the one
that minimizes ||f − g||∞ . For simplicity, we also denote it by f . Note that f satisfies supp(f ) ⊆ S, where
S = [m]\S. Further, since fS = 0 and
AS gS + AS gS = Ag = b = Af = AS fS + AS fS = AS fS
we have AS fS − gS = AS gS . Let AB be a linearly independent column subset of AS with maximal
cardinality, i.e. the column subset AN , where N = S \ B, is linearly dependent on AB . Hence, there is an
invertible square submatrix A′B ∈ Z|B|×|B| of AB and a vector v = (vB , 0N ) such that
 ′ 
AB
vB = AB vB = AS gS .
A′′B


Let r = (AS gS )B . Since A′B is invertible, there is a unique vector vB such that A′B vB = r. Observe that
|ri | =

X
j∈S

Ai,j gj ≤ ||A||∞

X
j∈S

|gj | <

1
||A||∞
=
.
2
nD ||A||∞
nD2

By Cramer’s rule vB (e) is quotient of two determinants. The denominator is det(A′B ) and hence at least one
in absolute value. For the numerator, the e-th column is replaced by r. Expansion according to this column
shows that the absolute value of the numerator is bounded by
D |B|
1
D X
|ri | <
·
≤
.
2
γA
γA nD
DγA
i∈B

Therefore, kf − gk∞ ≤ 1/(DγA ) and the statement follows.
Lemma 4.11. Let P
q ∈ Rm , p ∈ Rn and N = {e ∈ [m] : qe ≤ 0 or pT Ae ≤ 0}, where Aq = b and
−1
p = L b. Suppose e∈N |qe | < 1/ρA . Then there is a non-negative feasible kernel-free vector f such that
supp(f ) ⊆ E\N and kf − qk∞ < 1/(DγA ).
Proof. We apply Lemma 4.10 to q with S = N . Then, there is a non-negative feasible vector f such that
supp(f ) ⊆ E\N and kf − qk∞ < 1/(DγA ). By Lemma 2.1, f can be expressed as a sum of a convex
combination of basic feasible solutions plus a vector w in the kernel of A. Moreover, all vectors in this
representation are sign compatible with f , and in particular w is non-negative too.
26

P
Suppose for contradiction that w 6= 0. By definition, 0 = pT Aw = e∈[m] pT Ae we and since w ≥ 0 and
w 6= 0, it follows that there is an index e ∈ [m] satisfying we > 0 and pT Ae ≤ 0. Since f and w are sign
compatible, we > 0 implies fe > 0. On the other hand, as pT Ae ≤ 0 we have e ∈ N and thus fe = 0. This
is a contradiction, hence w = 0.
Using Corollary 4.1, for any point x(0) ∈ X there is a point x(t) ∈ X such that α(t) ∈ [1/2, 1/h0]. Thus,
we can assume that α(0) ∈ [1/2, 1/h0] and work with h ≤ h0 /2, where h0 = cmin /(2D||c||1 ). We generalize
next [SV16b, Lemma 5.4].
Lemma 4.12. Suppose x(t) ∈ X such that α(t) ∈ [1/2, 1/h0], h ≤ h0 /2 and ε ∈ (0, 1). Then, for any k ≥
h−1 ln(8mρA Ψ(0) /ε) there is a non-negative feasible kernel-free vector f such that kx(t+k) − f k∞ < ε/(DγA ).
def

Proof. Let β (k) = 1 − (1 − h)k . By (22), vector q (t,k) satisfies Aq (t,k) = b and thus Lemma 4.6 yields
kx(t+k) − β (k) q (t,k) k∞ = (1 − h)k · kx(t) k∞ ≤ exp{−hk} · Ψ(0) ≤ ε/(8mρA ).

(27)

Using Corollary 4.1, we have x(t+k) ∈ X such that α(t+k) ∈ (1/2, 1/h0) for every k ∈ N+ . Let Fk = Qk ∪ Pk ,
(t,k)
≤ 0} and Pk = {e ∈ [m] : ATe p(t,k) ≤ 0}. Then, for every e ∈ Qk it holds
where Qk = {e ∈ [m] : q e
|q (t,k)
| ≤ [β (k) ]−1 · |xe(t+k) − β (k) q (t,k)
| ≤ ε/(7mρA).
e
e

(28)

By Lemma 4.6, kx(·) k ≤ Ψ(0) . Moreover, by (22) for every e ∈ Pk we have
xe(t+k)

=

x(t)
e

k+t−1
Y h
i=t

≤
≤
≤


i
T (i)
1 + h c−1
−1
e Ae p

n
o
T (t,k)
x(t)
e · exp −hk + (h/ce ) · Ae p
exp {−hk} · Ψ(0)
ε/(8mρA ),

and by combining the triangle inequality with (27), it follows for every e ∈ Pk that
h
i
(k) −1
(t+k)
(k) (t,k)
(t+k)
|
≤
[β
]
·
|x
−
β
q
|
+
|x
|
|q (t,k)
e
e
e
e
≤
≤

Therefore, (28) and (29) yields that
X

e∈Fk

[β (k) ]−1 · ε/(4mρA )
ε/(3mρA).

(29)

| ≤ m · ε/(3mρA ) ≤ ε/(3ρA ).
|q (t,k)
e

(30)

(t,k)

and N = Fk , it follows by (30) that there is a non-negative feasible
By Lemma 4.11 applied with q e
kernel-free vector f such that supp(f ) ⊆ E\N and
kf − q (t,k) k∞ < ε/(3DγA ).
By Lemma 4.5, we have kq (t,k) k∞ ≤ mD2 kb/γA k1 and since Ψ(0) ≥ mD2 kb/γA k1 , it follows that
kx(t+k) − f k∞ = kx(t+k) − β (k) q (t+k) + β (k) q (t+k) − f k∞

k

≤ kx(t+k) − β (k) q (t+k) k∞ + kf − q (t+k) k∞ + (1 − h) kq(t+k) k∞

≤

ε
ε
ε · mD2 kb/γA k1
ε
≤
+
+
.
8mρA
3DγA
DγA
8mρA · Ψ(0)
27

4.5

x(k) is ε-Close to an Optimal Solution

Recall that N denotes the set of non-optimal basic feasible solutions of (5) and Φ = ming∈N cT g − opt. For
completeness, we prove next a well known inequality [PS82, Lemma 8.6] that lower bounds the value of Φ.
Lemma 4.13. Suppose A ∈ Rn×m has full row rank, b ∈ Rn and c ∈ Rm are integral. Then, Φ ≥ 1/(DγA )2 .
Proof. Let g = (gB , 0) be an arbitrary basic feasible solution with basis matrix AB , where gB (e) 6= 0 and
|supp(gB )| = n. We write M−i,−j to denote the matrix M with deleted i-th row and j-th column. Let Qe
be the matrix formed by replacing the e-th column of AB by the column vector b. Then, by Cramer’s rule,
we have

n
j+k
−1
· bk · det γA
[AB ]−k,−j
1 X (−1)
1
det(Qe )

=
.
≥
|gB (e)| =
−1
det(AB )
γA
Dγ
det γA AB
A
k=1

Note that all components of vector gB have denominator with equal value, i.e. det(AB ). Consider an
arbitrary non-optimal basic feasible solution g and an optimal basic feasible solution f ⋆ . Then, ge = Ge /G
and fe⋆ = Fe /F are rationals such that Ge , G, Fe , F ≤ DγA for every e. Further, let re = ce (Ge F − Fe F ) ∈ Z
for every e ∈ [m], and observe that
cT (g − f ⋆ ) =

X
e

ce (ge − fe⋆ ) =

1 X
re ≥ 1/(DγA )2 ,
GF e

where the last inequality follows by cT (g − f ⋆ ) > 0 implies

P

e re

≥ 1.

Lemma 4.14. Let f be a non-negative feasible kernel-free vector and ε ∈ (0, 1) a parameter. Suppose
for every non-optimal basic feasible solution g, there exists an index e ∈ [m] such that ge > 0 and fe <
ε/(2mD3 γA ||b||1 ). Then, kf − f ⋆ k∞ < ε/(DγA ) for some optimal f ⋆ .
Proof. Let C = 2D2 ||b||1 . Since f is kernel-free, by Lemma 2.1 it can be expressed as a convex combination
Pℓ
Pm
of sign-compatible basic feasible solutions f = i=1 αi f (i) + i=ℓ+1 αi f (i) , where f (1) , . . . , f (ℓ) denote the
(i)

(i)

optimal solutions. By Lemma 3.1, fe > 0 implies fe ≥ 1/(DγA ). By the hypothesis, for every non-optimal
f (i) , i.e. i ≥ ℓ + 1, there exists an index e(i) ∈ [m] such that
(i)

1/(DγA ) ≤ fe(i)
Thus,
αi /(DγA ) ≤

(i)
αi fe(i)

≤

Pm

and fe(i) < ε/(mDγA · C).
m
X
j=1

(j)

αj fe(i) = fe(i) < ε/(mDγA · C),

and hence i=ℓ+1 αi ≤ ε/C. By Lemma 3.1, we have kf (j) k∞ ≤ Dkb/γAk1 = C/(2DγA ) for every j. Let
Pℓ
Pm
β ≥ 0 be an arbitrary vector satisfying i=1 βi = i=ℓ+1 αi . Let νi = αi + βi for every i ∈ [ℓ] and let
Pℓ
f ⋆ = i=1 νi f (i) . Then, f ⋆ is an optimal solution and we have
kf ⋆ − f k∞ =

ℓ
X
i=1

βi f (i) −

≤ max

i∈[1:m]

f (i)

m
X

i=ℓ+1

∞

·

αi · f (i)

ℓ
X

βi +

i=1

∞
m
X

i=ℓ+1

αi

!

≤

2ε
ε
C
=
.
·
C 2DγA
DγA

In the following lemma, we extend the analysis in [SV16b, Lemma 5.6] from the transshipment problem to positive linear programs. Our result crucially relies on an argument that uses the parameter
Φ = ming∈N cT g − opt. It is here, where our analysis incurs the linear step size dependence on Φ/opt
and the quadratic dependence on opt/Φ for the number of steps.
28

An important technical detail is that the first regime incurs an extra (Φ/opt)-factor dependence. At
first glance, this might seem unnecessary due to Corollary 4.1, however a careful analysis shows its necessity
(see (35) for the inductive argument). Further, we note that the undirected Physarum dynamics (7) satisfies
(t)
(t)
(0)
xmin ≥ (1 − h)t · xmin , whereas the directed Physarum dynamics (6) might yield a value xmin which decreases
(0)
with faster than exponential rate. As our analysis incurs a logarithmic dependence on 1/xmin, it is prohibitive
(t)
to decouple the two regimes and give bounds in terms of log(1/xmin ), which would be necessary as x(t) is
the initial point of the second regime.
Lemma 4.15. Let g be an arbitrary non-optimal basic feasible solution. Given x(0) ∈ X and its corresponding α(0) , the Physarum dynamics (6) initialized with x(0) runs in two regimes:
(i) The first regime is executed when α(0) 6∈ [1/2, 1/h0] and computes a point x(t) ∈ X such that
α(t) ∈ [1/2, 1/h0]. In particular, if α(0) < 1/2 then h ≤ (Φ/opt) · (α(0) h0 )2 and t = 1/h. Otherwise,
if α(0) > 1/h0 then h ≤ Φ/opt and t = ⌊log1/(1−h) [h0 (α(0) − 1)/(1 − h0 )]⌋.
(ii) The second regime starts from a point x(t) ∈ X such that α(t) ∈ [1/2, 1/h0], it has step size h ≤
(0)
(opt/Φ) · h20 /2 and for any k ≥ 4 · cT g/(hΦ) · ln(Ψ(0) /εxmin ), guarantees the existence of an index
(t+k)
e ∈ [m] such that ge > 0 and xe
< ε.
Proof. Similar to the work of [BBD+ 13, SV16b], we use a potential function that takes as input a basic
feasible solution g and a step number ℓ, and is defined by
X
ge ce ln x(ℓ)
Bg(ℓ) :=
e .
e∈[m]

Since

(ℓ+1)
xe

=

(ℓ)
xe (1

Bg(ℓ+1) − Bg(ℓ)

T (ℓ)
+ h(ℓ) [c−1
− 1]), we have
e · Ae p

 T (ℓ)

(ℓ+1)
X
X
Ae p
xe
−1
ge ce ln 1 + h(ℓ)
=
ge ce ln (ℓ) =
ce
xe
e
e


i
h
h
i
X
ATe p(ℓ)
≤ h(ℓ)
g e ce
− 1 = h(ℓ) −cT g + [p(ℓ) ]T Ag = h(ℓ) −cT g + bT p(ℓ) .
ce
e
(ℓ+1)

Let f ⋆ be an optimal solution to (5). In order to lower bound Bf ⋆
x − x2 , for all x ∈ [− 21 , 21 ]. Then, we have
(ℓ+1)

Bf ⋆

(ℓ)

− Bf ⋆

(31)

(ℓ)

− Bf ⋆ , we use the inequality ln(1 + x) ≥

(ℓ+1)

=
≥
≥


 T (ℓ)

⋆
(ℓ) Ae p
−
1
f
c
ln
1
+
h
e e
(ℓ)
ce
xe
e
e



2 !
T (ℓ)
T (ℓ)
X
(ℓ) 2 Ae p
(ℓ) Ae p
⋆
− 1 − [h ]
−1
f e ce h
ce
ce
e


h(ℓ) bT p(ℓ) − opt − h(ℓ) · (1/2α(ℓ) h0 )2 · opt ,
X

fe⋆ ce ln

xe

=

X

(32)

where the last inequality follows by combining
X
T (ℓ)
(ℓ) T
⋆
T (ℓ)
fe⋆ ce [(c−1
− opt,
e Ae p ) − 1] = [p ] Af − opt = b p
e

T (ℓ)

kA p

k∞ ≤ D||c||1 /α (by Lemma 4.4 and Lemma 4.8 applied with x(ℓ) ∈ X), h0 = cmin /(4D||c||1 ) and
X
T (ℓ)
h(ℓ)
fe⋆ ce · (c−1
− 1)2 ≤ h(ℓ) (2D||c||1 /α(ℓ) cmin )2 opt = h(ℓ) (1/2α(ℓ) h0 )2 opt.
e Ae p
(ℓ)

e

Further, by combining (31), (32), cT g − opt ≥ Φ for every non-optimal basic feasible solution g and provided
that the inequality h(ℓ) (1/2α(ℓ) h0 )2 opt ≤ Φ/2 holds, we obtain




h(ℓ) Φ
Φ
(ℓ+1)
(ℓ)
T
(ℓ)
T (ℓ)
T
(ℓ)
≥ Bg(ℓ+1) − Bg(ℓ) +
.
(33)
c g − opt −
Bf ⋆ − Bf ⋆ ≥ h
b p −c g +h
2
2
29

Using Corollary 4.1, we partition the Physarum dynamics (6) execution into three regimes, based on
α(0) . For every i ∈ {1, 2, 3}, we show next that the i-th regime has a fixed step size h(ℓ) = hi such that
h(ℓ) (1/2α(ℓ) h0 )2 opt ≤ Φ/2, for every step ℓ in this regime.
By Lemma 4.9, for every i ∈ {1, 2, 3} it holds for every step ℓ in the i-th regime that
α(ℓ) = 1 − (1 − hi )ℓ · (1 − α(0) ).

(34)

Case 1: Suppose α(0) > 1/h0 . Notice that h(ℓ) = Φ/opt suffices, since 1/(2α(ℓ) h0 ) < 1/2 for every
α > 1/h0 . Further, by applying (34) with α(t) := 1/h0 , we have t = ⌊log1/(1−h(ℓ) ) [h0 (α(0) − 1)/(1 − h0)]⌋ ≤
(opt/Φ) · log(α(0) h0 ). Note that by (34) the sequence {α(ℓ) }ℓ≤t is decreasing, and by Corollary 4.1 we have
1 < α(t) ≤ 1/h0 .
Case 2: Suppose α(0) ∈ (0, 1/2). By (34) the sequence {α(ℓ) }ℓ∈N is increasing and by Corollary 4.1 the
regime is terminated once α(ℓ) ∈ [1/2, 1). Observe that h(ℓ) = (Φ/opt) · (α(0) h0 )2 suffices, since α(0) ≤ α(ℓ) .
Then, by (34) applied with α(t) := 1/2, this regime has at most t = (opt/Φ) · (1/α(0) h0 )2 steps.
Case 3: Suppose α(0) ∈ [1/2, 1/h0]. By (34) the sequence {α(ℓ) }ℓ∈N converges to 1 (decreases if α(0) ∈
(1, 1/h0 ] and increases when α(0) ∈ [1/2, 1). Notice that h(ℓ) = (Φ/opt)·h20 /2 suffices, since 1/2 ≤ α(ℓ) ≤ 1/h0
for every ℓ ∈ N. We note that the number of steps in this regime is to be determined soon.
Hence, we conclude that inequality (33) holds. Further, using Case 1 and Case 2 there is an integer t ∈ N
such that α(t) ∈ [1/2, 1/h0]. Let k ∈ N be the number of steps in Case 3, and let h := (Φ/opt) · h20 /2. Then,
for every ℓ ∈ {t, . . . , t + k − 1} it holds that h(ℓ) = h and thus
(ℓ)

(t+k)

Bf ⋆

(0)

− Bf ⋆ ≥ Bg(t+k) − Bg(0) +

t+k−1
X
ℓ=0

hΦ
h(ℓ) Φ
≥ Bg(t+k) − Bg(0) + k ·
.
2
2

(35)

(ℓ)

By Lemma 4.6, Bg ≤ cT g · ln Ψ(0) for every basic feasible solution g and every ℓ ∈ N, and thus
Bg(t+k)

≤
≤
≤

hΦ
(t+k)
(0)
+ Bg(0) + Bf ⋆ − Bf ⋆
2
hΦ
(0)
+ cT g · ln Ψ(0) + opt · ln Ψ(0) − opt · ln xmin
−k ·
2
Ψ(0)
hΦ
+ 2cT g · ln (0) .
−k ·
2
x

−k ·

min

(t+k)

Suppose for the sake of a contradiction that for every e ∈ [m] with ge > 0 it holds xe
> ε. Then,
(t+k)
(0)
Bg
> cT g · ln ε yields k < 4 · cT g/(hΦ) · ln(Ψ(0) /(εxmin )), a contradiction to the choice of k.

4.6

Proof of Theorem 4.2

By Corollary 4.1 and Lemma 4.15, if x(0) ∈ X such that α(0) > 1/h0 , we work with h ≤ Φ/opt and
after t = ⌊log1/(1−h) [h0 (α(0) − 1)/(1 − h0 )]⌋ ≤ (opt/Φ) · log(α(0) h0 ) steps, we obtain x(t) ∈ X such that
α(t) ∈ (1, 1/h0 ]. Otherwise, if α(0) ∈ (0, 1/2) we work with h ≤ (Φ/opt) · (α(0) h0 )2 and after t = 1/h
steps, we obtain x(t) ∈ X such that α(t) ∈ [1/2, 1). Hence, we can assume that α(t) ∈ [1/2, 1/h0] and set
h ≤ (Φ/opt) · h20 . Then, the Lemmas in Subsection 4.4 and 4.5 are applicable.
Let E1 := Dkb/γAk1 ||c||1 , E2 := 8mρA Ψ(0) , E3 := 2mD3 γA kbk1 and E4 := 8mD2 kbk1. Consider an
arbitrary non-optimal basic feasible solution g.
By Lemma 3.1, we have cT g ≤ E1 and thus both Lemma 4.12 and Lemma 4.15 are applicable with h,
(0)
⋆
ε := ε/E4 and any k ≥ k0 := 4E1 /(hΦ) · ln[(E2 / min{1, xmin}) · (DγA /ε⋆ )]. Hence, by Lemma 4.15, the
(t+k)
< ε⋆ /(DγA ).
Physarum dynamics (6) guarantees the existence of an index e ∈ [m] such that ge > 0 and xe
(t+k)
Moreover, by Lemma 4.12 there is a non-negative feasible kernel-free vector f such that kx
− f k∞ <
ε⋆ /(DγA ). Thus, for the index e it follows that ge > 0 and fe < 2ε⋆ /DγA = (ε/2) · (4/E4 DγA ) = ε/(2E3 ).
Then, Lemma 4.14, yields kf −f ⋆ k∞ < ε/(2DγA ) and by triangle inequality we have kx(k) −f ⋆ k∞ < ε/(DγA ).
30

By construction, ρA = max{DγA , nD2 ||A||∞ } ≤ nD2 γA ||A||∞ . Let E2′ = 8mnD2 γA ||A||∞ Ψ(0) and E5 =
2
· DγA = 82 m2 nD5 γA
||A||∞ kbk1 . Further, let C1 = E1 and C2 = E5 . Then, the statement follows for
(0)
any k ≥ k1 := 4C1 /(hΦ) · ln(C2 Ψ(0) /(ε · min{1, xmin })).

E2′ E4

4.7

Preconditioning

What can be done if the initial point is not strongly dominating? For the transshipment problem it suffices
to add an edge of high capacity and high cost from every source node to every sink node [BBD+ 13, SV16b].
This will make the instance strongly dominating and will not affect the optimal solution.
In this section, we generalize this observation to positive linear programs. We add an additional column
equal to b and give it sufficiently high capacity and cost. This guarantees that the resulting instance is
strongly dominating and the optimal solution remains unaffected. We state now our algorithmic result.
Theorem 4.16. Given an integral LP (A, b, c > 0), a positive x(0) ∈ Rm and a parameter ε ∈ (0, 1). Let
([A | b], b, (c, c′ )) be an extended LP with c′ = 2C1 and z (0) := 1 + DS ||x||∞ ||A||1 ||b||1 .8 Then, (x(0) ; z (0) ) is a
strongly dominating starting point of the extended problem such that y T [A | b](x(0) , z (0) ) ≥ 1, for all y ∈ Y . In
particular, the Physarum dynamics (6) initialized with (x(0) , z (0) ) and a step size h ≤ h20 /C3 , outputs for any
(0)
k ≥ 4C1 · (DγA )2 /h · ln(C2 Υ(0) /(ε · min{1, xmin})) a vector (x(k) , z (k) ) > 0 such that dist(x(k) , X⋆ ) < ε/(DγA )
and z (k) < ε/(DγA ), where here Υ(0) := max{Ψ(0) , z (0) }.
Theorem 4.16 subsumes [SV16b, Theorem 1.2] for flow problems by giving a tighter asymptotic convergence rate, since for the transshipment problem A is a totally unimodular matrix and satisfies D = DS = 1,
γA = 1, ||A||∞ = 1 and Φ = 1. We note that the scalar z (0) depends on the scaled determinant DS , see
Theorem 1.3.
4.7.1

Proof of Theorem 4.16

In the extended problem, we concatenate to matrix A a column equal to b such that the resulting constraint
matrix becomes [A | b]. Let c′ be the cost and let x′ be the initial capacity of the newly inserted constraint
column. We will determine c′ and x′ in the course of the discussion. Consider the dual of the max-flow like
LP for the extended problem. It has an additional variable z ′ and reads

min xT z + x′ z ′ : z ≥ 0; z ′ ≥ 0; z ≥ AT y; z ′ ≥ bT y; bT y = 1 .
In any optimal solution, z ′ = bT y = 1 and hence the dual is equivalent to

min xT z + x′ : z ≥ 0; z ≥ AT y; bT y = 1 .

The strongly dominating set of the extended problem is therefore equal to





x
x
m+1
T
X=
∈
R
:
y
[A
|
b]
>
0
for
all
y
∈
Y
.
>0
x′
x′

(36)

(37)

The defining condition translates into x′ > −y T Ax for all y ∈ Y . We summarize the discussion in the
following Lemma.
P
Lemma 4.17. Given a positive x ∈ Rm , let ρ := ||b||1 DS and x′ := 1 + ρ||A||1 ||x||∞ , where ||A||1 := i,j |Ai,j |
and DS := max{| det(A′ )| : A′ is a square sub-matrix of A}. Then, (x; x′ ) is a strongly dominating starting
point of the extended problem such that y T [A | b](x; x′ ) = y T Ax + x′ ≥ 1, for all y ∈ Y .
Proof. We show first that maxy∈Y ||y||∞ ≤ ρ implies the statement. Let y ∈ Y be arbitrary. Since |y T Ax| ≤
||A||1 ||x||∞ ||y||∞ , we have maxy∈Y |y T Ax| ≤ ρ||A||1 ||x||∞ = x′ − 1 and hence y T [A | b](x; x′ ) ≥ 1.
P
8
We denote by ||A||1 :=

i,j

|Ai,j |, i.e. we interpret matrix A as a vector and apply to it the standard ℓ1 norm.

31

It remains to show that maxy∈Y ||y||∞ ≤ ρ. The constraint polyhedron of the dual (36) is given in matrix
notation as







 ≥
Im×m −AT
0m 

z
z
=  1  .
P (ext) :=
∈ Rm+n :  0Tm
bT 
y

 y
0m
Im×m 0m×n
≥

Let us denote the resulting constraint matrix and vector by M ∈ R2m+1×m+n and d ∈ R2m+1 , respectively.
Note that if b = 0 then the primal LP (23) is either unbounded or infeasible. Hence, we consider the
non-trivial case when b 6= 0. Observe that the polyhedron P (ext) is not empty, since for any y such that
bT y = 1 there is z = max{0, AT y} satisfying (z; y) ∈ P (ext) . Further, P (ext) does not contain a line (see
Subsection 4.3) and thus P (ext) has at least one extreme point p′ ∈ P (ext) . As the dual LP (24) has a
bounded value (the target function is lower bounded by 0) and an extreme point exists (p′ ∈ P (ext) ), the
optimum is attained at an extreme point p ∈ P (ext) . Moreover, as every extreme point is a basic feasible
solution and matrix M has linearly independent columns (A has full row rank), it follows that p has m + n
tight linearly independent constraints.
Let MB(p) ∈ Rm+n×m+n be the basis submatrix of M satisfying MB(p) p = dB(p) . Since A, b are integral
and MB(p) is invertible, using Laplace expansion we have 1 ≤ | det(MB(p) )| ≤ ||b||1 DS = ρ. Let Qi denotes
the matrix formed by replacing the i-th column of MB(p) by the column vector dB(p) . Then, by Cramer’s
rule, it follows that |yi | = | det(Qi )/ det(MB(p) )| ≤ | det(MB(p) )| ≤ ρ, for all i ∈ [n].
It remains to fix the cost of the new column. Using Lemma 3.1, opt ≤ cT x(k) ≤ C1 for every k ∈ N, and
thus we set c′ := 2C1 .

4.8

A Simple Lower Bound

Building upon [SV16b, Lemma B.1], we give a lower bound on the number of steps required for computing an
ε-approximation to the optimum shortest path. In particular, we show that for the Physarum dynamics (6)
to compute a point x(k) such that dist(x(k) , X⋆ ) < ε, the required number of steps k has to grow linearly in
opt/(hΦ) and ln(1/ε).
Theorem 4.18. Let (A, b, c) be a positive LP instance such that A = [1 1], b = 1 and c = (opt opt + Φ)T ,
where opt > 0 and Φ > 0. Then, for any ε ∈ (0, 1) the discrete directed Physarum dynamics (6) initialized
with x(0) = (1/2, 1/2) and any step size h ∈ (0, 1/2], requires at least k = (1/2h) · max{opt/Φ, 1} · ln(2/ε)
(k)
(k)
steps to guarantee x1 ≥ 1 − ε, x2 ≤ ε. Moreover, if ε ≤ Φ/(2opt) then cT x(k) ≥ (1 + ε)opt as long as
k ≤ (1/2h) · max{opt/Φ, 1} · ln(2Φ/(ε · opt)).
(k)

Proof. Let c1 = opt and c2 = γopt, where γ = 1 + Φ/opt. We first derive closed-form expressions for x1 ,
(k)
(k)
(k)
(k)
(k)
(k)
(k)
(k) (k)
x2 , and x1 + x2 . Let s(k) = γx1 + x2 . For any k ∈ N, we have q1 + q2 = 1 and q1 /q2 =
(k)
(k)
(k)
(k)
(k)
(k) (k)
(k)
(k) (k)
(x1 /c1 )/(x2 /c2 ) = γx1 /x2 . Therefore, q1 = γx1 /s and q2 = x2 /s , and hence
(k)

x1
(k)

(k−1)

= (1 + h(−1 + γ/s(k−1) ))x1

(k)

(k−1)

(k)

and x2

(k−1)

(k−1)

= (1 + h(−1 + 1/s(k−1) ))x2
(k)

.

(38)

(k)

Further, x1 + x2 = (1 − h)(x1
+ x2
) + h, and thus by induction x1 + x2 = 1 for all k ∈ N.
(k)
(k−1)
(k)
Therefore, s(k) ≤ γ for all k ∈ N and hence x1 ≥ x1
, i.e. the sequence {x1 }k∈N is increasing and
(k)
the sequence {x2 }k∈N is decreasing. Moreover, since h(−1 + 1/s(k−1) ) ≥ h(1 − γ)/γ = −hΦ/(opt + Φ) and
using the inequality 1 − z ≥ e−2z for every z ∈ [0, 1/2], it follows by (38) and induction on k that

k


hΦ
1
2hΦ
(k)
(0)
x2 ≥ 1 −
.
x2 ≥ exp −k ·
opt + Φ
2
opt + Φ
(k)

Thus, x2 ≥ ε whenever k ≤ (1/2h) · (opt/Φ + 1) · ln(2/ε). This proves the first claim.
(k)
(k)
(k)
For the second claim, observe that cT x(k) = opt · x1 + γopt · x2 = opt · (1 + (γ − 1)x2 ). This is greater
(k)
than (1 + ε)opt iff x2 ≥ ε · opt/Φ. Thus, cT x(k) ≥ (1 + ε)opt as long as k ≤ (1/2h) · (opt/Φ + 1) · ln(2/(ε ·
opt/Φ)).

32

References
[BBD+ 13] Luca Becchetti, Vincenzo Bonifaci, Michael Dirnberger, Andreas Karrenbauer, and Kurt
Mehlhorn. Physarum can compute shortest paths: Convergence proofs and complexity bounds.
In ICALP, volume 7966 of LNCS, pages 472–483, 2013.
[BMV12] Vincenzo Bonifaci, Kurt Mehlhorn, and Girish Varma. Physarum can compute shortest paths.
Journal of Theoretical Biology, 309:121 – 133, 2012.
[Bon13]

Vincenzo Bonifaci. Physarum can compute shortest paths: A short proof. Inf. Process. Lett.,
113(1-2):4–7, 2013.

[Bon15]

Vincenzo Bonifaci. A revised model of network transport optimization in Physarum Polycephalum.
November 2015.

[Bon16]

Vincenzo Bonifaci. On the convergence time of a natural dynamics for linear programming. CoRR,
abs/1611.06729, 2016.

[CDS98]

Scott Shaobing Chen, David L. Donoho, and Michael A. Saunders. Atomic decomposition by
basis pursuit. SIAM Journal on Scientific Computing, 20(1):33–61, 1998.

[CLSW98] F. H. Clarke, Yu. S. Ledyaev, R. J. Stern, and P. R. Wolenski. Nonsmooth Analysis and Control
Theory. Springer-Verlag New York, Inc., Secaucus, NJ, USA, 1998.
[IJNT11] Kentaro Ito, Anders Johansson, Toshiyuki Nakagaki, and Atsushi Tero. Convergence properties
for the physarum solver. arXiv:1101.5249v1, January 2011.
[JZ12]

A. Johannson and J. Zou. A slime mold solver for linear programming problems. In CiE, pages
344–354, 2012.

[LaS76]

J. B. LaSalle. The Stability of Dynamical Systems. SIAM, 1976.

[MO08]

T. Miyaji and Isamu Ohnishi. Physarum can solve the shortest path problem on Riemannian
surface mathematically rigourously. International Journal of Pure and Applied Mathematics,
47:353–369, 2008.

[NYT00]

T. Nakagaki, H. Yamada, and Á. Tóth. Maze-solving by an amoeboid organism. Nature, 407:470,
2000.

[Phy]

http://people.mpi-inf.mpg.de/~mehlhorn/ftp/SlimeAusschnitt.webm.

[PS82]

Christos H. Papadimitriou and Kenneth Steiglitz. Combinatorial Optimization: Algorithms and
Complexity. Prentice-Hall, Inc., Upper Saddle River, NJ, USA, 1982.

[Sch99]

Alexander Schrijver. Theory of linear and integer programming. Wiley-Interscience series in
discrete mathematics and optimization. Wiley, 1999.

[SV16a]

Damian Straszak and Nisheeth K. Vishnoi. IRLS and slime mold: Equivalence and convergence.
CoRR, abs/1601.02712, 2016.

[SV16b]

Damian Straszak and Nisheeth K. Vishnoi. Natural algorithms for flow problems. In SODA,
pages 1868–1883, 2016.

[SV16c]

Damian Straszak and Nisheeth K. Vishnoi. On a natural dynamics for linear programming. In
ITCS, pages 291–291, New York, NY, USA, 2016. ACM.

[TKN07]

A. Tero, R. Kobayashi, and T. Nakagaki. A mathematical model for adaptive transport network
in path finding by true slime mold. Journal of Theoretical Biology, pages 553–564, 2007.

33

