ON A NATURAL DYNAMICS FOR LINEAR PROGRAMMING

arXiv:1511.07020v1 [] 22 Nov 2015

DAMIAN STRASZAK AND NISHEETH K. VISHNOI

Abstract. In this paper we study dynamics inspired by Physarum polycephalum (a slime mold)
for solving linear programs ([NYT00, IJNT11, JZ12]). These dynamics are arrived at by a local
and mechanistic interpretation of the inner workings of the slime mold and a global optimization
perspective has been lacking even in the simplest of instances. Our first result is an interpretation
of the dynamics as an optimization process. We show that Physarum dynamics can be seen as a
steepest-descent type algorithm on a certain Riemannian manifold. Moreover, we prove that the
trajectories of Physarum are in fact paths of optimizers to a parametrized family of convex programs,
in which the objective is a linear cost function regularized by an entropy barrier. Subsequently, we
rigorously establish several important properties of solution curves of Physarum. We prove global
existence of such solutions and show that they have limits, being optimal solutions of the underlying
LP. Finally, we show that the discretization of the Physarum dynamics is efficient for a class of
linear programs, which include unimodular constraint matrices. Thus, together, our results shed
some light on how nature might be solving instances of perhaps the most complex problem in P:
linear programming.

Contents
1. Introduction
2. Technical Overview
3. Preliminaries
4. Physarum as Steepest Descent on a Manifold
4.1. The Riemannian Structure and its Origin
4.2. Two Forces
4.3. Physarum as an Entropy Barrier Method
5. Existence of Solution
5.1. Basic Results
5.2. The Proof of Existence
6. Asymptotic Behavior of Physarum Trajectories
6.1. The Feasible Region
6.2. Convergence to the Optimal Value
6.3. Convergence to a Limit
7. Discretization
References
Appendix A. Dynamical Systems and Global Existence

Damian Straszak, École Polytechnique Fédérale de Lausanne (EPFL).
Nisheeth K. Vishnoi, École Polytechnique Fédérale de Lausanne (EPFL).

1
4
7
7
7
8
11
13
13
15
17
17
18
21
24
27
28

1. Introduction
There is growing evidence that several fundamental processes in nature are inherently computational and are best viewed from the standpoint of computation. Consider the case of Physarum
polycephalum (slime mold), a single celled organism which has been the source of much excitement
among biologists and computer scientists due to its ability to solve complex optimization problems.
This started with an experiment [NYT00] that showed the slime mold could solve the shortest path
problem on a maze. Roughly speaking, the slime mold is a collection of tubes which it uses to
transport food across its body and it is in this process that it shows the ability to compute. Soon
after, the time evolution of Physarum was captured by mathematical biologists using the language
of dynamical systems giving rise to a broad class of dynamics for basic computational problems
such as shortest paths, transshipment problems, and linear programs (see [TKN07, IJNT11, JZ12]).
The Physarum dynamics are highly nonlinear and decentralized and, thus, there is no a priori
reason for them to converge anywhere, let alone to solutions of problems, such as those listed
above, where a global optima is sought. However, experiments and simulations suggest that the
solution not only converges, it converges quickly! Two sets of questions arise naturally from a
computational viewpoint: (1) While the dynamical system gives a mechanistic insight into the
workings of a Physarum, it fails to explain what is going on globally. Can we explain Physarum
dynamics from an optimization perspective? (2) Can we rigorously explain what experiments,
simulations and partial results about the Physarum dynamics suggest, namely the (a) existence of
a solution, (b) convergence to an optimal solution and (c) time to convergence of the (continuous
and discretized) dynamics?
In this paper we study these questions for the Physarum dynamics for linear programming proposed by [JZ12] which generalizes earlier models for shortest path and the transshipment problem.
For question (1) there are no known answers even in the simplest setting of shortest path: on the
one hand the Physarum dynamics resemble the gradient descent method; however, one can show
that it is not. On the other hand, the dynamics are reminiscent of interior point methods; however,
what makes the Physarum dynamics remarkable is that one can start from outside the feasible
region of the linear program and yet converge to the optimal point. For question (2), the situation
is slightly better. For various special cases such as shortest path [BMV12, BBD+ 13] and the transshipment problem [IJNT11, SV] (a)-(c) have been answered. However, these results heavily rely on
the special structure of the problem. For LPs, the prior results prove convergence assuming that
(i) the solution exists (ii) the feasible region of the LP is bounded and (iii) the optimizer is unique.
As we argue later, justifying or removing these ends up being quite hard. Current results do not
bound the convergence time.
The main contribution of this paper is to provide answers to all of these questions: We show that
Physarum dynamics can be seen as a steepest-descent type algorithm, however, not in Euclidean
space, rather in a space endowed with a Riemannian metric obtained from an entropy-like function.
As a consequence, we show that Physarum dynamics for linear programming are obtained by
balancing two forces in this Riemannian manifold: the need to reduce cost and the need to be
feasible. Moreover, we show that the continuous trajectories of Physarum are in fact paths of
optimizers to a parametrized family of convex programs, in which the objective is a linear cost
function regularized by an entropy barrier. Subsequently, we establish the global existence of
solutions of Physarum dynamics and show that they have limits, being optimal solutions of the
underlying LP. Finally, we present a time-bound on a discretization of the dynamics to yield
an algorithm that is provably efficient for a class of linear programs, which include unimodular
constraint matrices. Our proofs synthesize concepts and tools from several distinct areas such as
Riemannian manifolds, convex optimization and dynamical systems. These techniques should be
useful in providing insights into the plethora of steepest-descent type algorithms in optimization,
machine learning and, recently, in theory.
1

The Physarum Dynamics for Linear Programming. Consider a linear program in the standard form
min c⊤ x

(1)
Zm×n ,

Zn>0

Zm

s.t. Ax = b, x ≥ 0

where A ∈
c∈
and b ∈
and which has a feasible solution. To make the exposition
clear we assume that A is full-rank, in other words: rank(A) = m.1 We now describe the Physarum
dynamics for linear programming.
Consider any vector x ∈ Rn with x > 0 and let W be the diagonal matrix with entries xi/ci . Let
def

def

L = AW A⊤ and p ∈ Rn is the solution to Lp = b.2 Let q = W A⊤ p. The Physarum dynamics for
the linear program given by (A, b, c) then is
ẋ = q − x

(2)
defined over the domain Ω =
equations:

Rn>0 .

The above is just a simplified notation for the coupled differential
∀i,

dxi (t)
= qi (t) − xi (t).
dt

This can be also rewritten as
(3)

ẋ = W (A⊤ L−1 b − c).

The dynamical system has an initial condition of the form x(0) = s for some s > 0. Note that it
is not required that s is feasible, i.e., As = b. Often, the exposition and proofs are simpler when s
is also feasible. We often refer to the general case as Physarum dynamics with infeasible start and
this special case as with feasible start. Notice that when A is the |V | × |E| incidence matrix of a
graph G = (V, E), then if we set up an electrical network where the edge i has resistance xi/ci , then
L is the corresponding graph Laplacian and q the electrical flow in each edge corresponding to the
demand vector b. Viewing the slime mold as a collection of tubes organized in a network where
the length of edge i is ci and cross-section xi and noting that the equations of Poiseuille flows are
identical to electrical flows, (2) captures the local time evolution of the cross-section of each tube.
Our Contributions. The first problem we have at hand is whether the system (2) has a solution
such that x(t) > 0 for all t > 0. Unfortunately there are no general theorems in dynamical systems
to establish this and proving existence can be difficult. Indeed, this existence problem turns out to
be non-trivial for Physarum dynamics and is our first result.
Theorem 1.1 (Existence of Solution). For any initial condition s ∈ Rn>0 , the Physarum dynamics ẋ = q − x has a unique solution x : [0, ∞) → Rn>0 with x(0) = s.
First we note that existence of a feasible solution is important: One can prove that whenever
the problem (1) is infeasible, the solutions to (2) exist only for finite time intervals. Further,
the problem of existence is difficult not because the dynamics is in continuous time; the problem
remains non-trivial even if we consider a discretization. Now that we know that the solution exists,
our next set of results establish that no matter where one starts the Physarum dynamics from, as
long as it is in the positive orthant, one converges to an optimal solution of the linear program.
Theorem 1.2 (Convergence to an Optimal Solution). For any initial condition s ∈ Rn>0 ,
consider the solution x : [0, ∞) → Rn>0 to the Physarum dynamics with x(0) = s. Denote by x⋆ an
optimal solution of the considered linear program. Then:
(1) |c⊤ x(t) − c⊤ x⋆ | = O(e−αt ) for some positive α, which only depends upon A, b, c and s,
(2) the limit x∞ = limt→∞ x(t) exists, is a feasible point and c⊤ x∞ = c⊤ x⋆ .
1This assumption is not necessary for our results.
2One can show that p exists and is unique because of the full-rank assumption on A, otherwise one can use the

pseudo-inverse.
2

A few remarks are in order: (1) The theorem also gives an upper bound on the time to convergence
and (2) it proves that limits of trajectories of the Physarum exist even when the optimal solution is
not unique. Given the amount of technical effort required in proving this, one may wonder again if
this has something to do with continuous time. The answer is (again) no: Starting with Karmarkar
himself [Kar90], it has been observed that the good convergence properties of his algorithm for linear
programming [Kar84] arise from the good geometric properties of the set of continuous trajectories
which underlie his method [BL89].
Our next result concerns the computational abilities of a discretization of the Physarum dynamics.
Theorem 1.3 (Convergence Time of Discrete Physarum Dynamics; see Theorem 7.1).
Consider the discretization of the Physarum dynamics, i.e., x(k +1) = (1−h)x(k)+hq(k). Suppose
we initialize the Physarum algorithm with x(0) = s, i.e., As = b and M −1 ≤ si ≤ M for every
i = 1, . . . , n and some M ≥ 1. Assume additionally that c⊤ s ≤ M · opt. Choose any ε > 0

def
M
steps x(k) is a feasible solution with:
and3 let h = 61 · ε · Cs−2 · D −2 . Then after k = O ln
ε2 h2
⊤
opt ≤ c x(k) ≤ (1 + ε) · opt.
Thus, when the maximum cost is polynomial and the maximum sub-determinant of A is polynomially bounded, the discretization of Physarum dynamics is efficient.
Finally, we come to the conceptual question which has resisted an answer even for the shortest
path Physarum dynamics: While each equation in (2) gives us a local update rule, can the Physarum
dynamics be obtained from an optimization viewpoint? We answer this question affirmatively. For
preliminaries on Riemannian manifolds see Section 3.
Theorem 1.4 (Optimization Viewpoint of Physarum Dynamics; see Theorem 4.1). In
case of feasible start one can interpret Physarum as a gradient flow on a Riemannian manifold.
In other words it is of the form ẋ = ∇f (x), where f (x) = c⊤ x is the objective, but the gradient is
computed with respect to a certain Riemannian metric on the feasible region. In the infeasible case,
Physarum dynamics combines two forces; minimizing the objective and aiming for feasibility. Both
can be interpreted similarly as descent directions on a manifold.
We remark that the above
P mentioned Riemannian metric is induced by the Hessian of the generalized entropy function i ci xi ln xi ci and provides a physical meaning to the dynamics. Moreover,
using convex programming duality, we can show that the trajectories of Physarum with a feasible
start are in fact paths of optimizers to a parametrized family of convex programs, which objective is a linear cost function regularized by an entropy barrier function; see Theorem 4.2 for a
formal statement. Finally, note that the entropy barrier, unlike the log-barrier function is not a
self-concordant barrier function; yet the trajectories converge as e−αt as Theorem 1.2 indicates.
Related Work. There are rigorous results for existence, convergence and time for the shortest
path problem [BBD+ 13]. Our previous paper [SV] deals with the discrete Physarum dynamics
transshipment problems on graphs and uses the underlying graph structure heavily to obtain versions of Theorem 1.3. The LP Physarum dynamics was introduced by [JZ12]. There is a large
body of literature in the mathematical programming community which studies the geometric properties of continuous trajectories of interior point methods such as projective and affine scaling; see
[Kar90, BL89, NT02]. As mentioned before the Physarum dynamics is none of these methods (it
works from an infeasible start as well), but it bears some similarity to primal affine scaling methods,
see, e.g., [LR90]. However, for affine scaling methods, the convergence rate can be provably bad:
1/t, see, e.g., [MS89].
3C = Pn c and D = max{| det(A′ )| : A′ is a square submatrix of A}.
s
i=1 i
3

2. Technical Overview
We start by presenting an overview of our result which presents Physarum dynamics from an
optimization viewpoint (Theorem 1.4). We assume 1.1 and 1.2 for the moment. Denote
def

P (x) = q − x = W (A⊤ (AW A⊤ )−1 b − c),

so that the Physarum dynamics becomes ẋ = P (x). Knowing that (2) is in fact solving an optimization problem one would expect that it implements some gradient-descent-type of procedure.
This would mean that P (x) = ∇Φ(x), for some function Φ(x). However, by a simple calculation, it
turns out that P (x) is not a gradient of any function. Similarly, P (x) does not represent Newton’s
method nor steepest descent in any standard norm. To understand Physarum, we need to move
from the Euclidean world to Riemannian manifolds.
We consider the manifold Ω = Rn>0 , to every point x ∈ Ω we assign a positive definite matrix
def

H(x) = Diag (c1/x1 , . . . , cn/xn ) ,
which defines an inner product hu, vix = u⊤ H(x)v on the tangent space at x (which is just Rn in
this case), Ω becomes then a Riemannian manifold. We are mainly interested in the submanifold
Ωf = {x ∈ Ω : Ax = b}, which we assume to be nonempty for our discussion. It inherits the
Riemannian structure from Ω. Note that the Nash Embedding Theorem guarantees that such
a manifold always has an isometric embedding in an Euclidean space (possibly in much higher
dimension). In our case, the embedding can be given by an explicit formula:
√
√
F (x) = 2( c1 x1 , . . . , cn xn )⊤ .
This gives as a geometric realization of the curvature imposed on Ωf .
Let us now see that P (x), when restricted to Ωf is just the gradient of the objective f (x) = c⊤ x
when viewed in the local geometry. To show this, we need to argue that it satisfies
f (x + h) − f (x) − hh, P (x)ix = o(khkx )

over h from the tangent space at x, i.e., from Tx (Ωf ) = {h : Ah = 0}. We calculate:
f (x + h) − f (x) − hh, P (x)ix = c⊤ h − h⊤ H(x)P (x)

= c⊤ h − h⊤ (A⊤ L−1 b − c)
= (Ah)⊤ L−1 b = 0.

If x lies outside of the feasible region a similar interpretation of P (x) is possible. However in this
case P (x) decomposes into two parts: Pf (x) which is the feasibility direction and Po (x) which
aims for optimality. Both can be motivated as descent directions with respect to the above defined
Riemannian structure on Ω.
The next result gives another interpretation of the Physarum dynamics in the feasible region. Namely let us define x(µ), for µ > 0 to be the minimizer of c⊤ x + µ−1 f (x) over x ∈
{x : Ax = b, x ≥ 0}, where f is the following entropy-like function
def

f (x) =

n
X

ci xi ln(ci xi ).

i=1

It turns out that using elementary convex programming techniques such as duality and KKT
conditions, we can prove that x(µ) satisfies the Physarum dynamics. This is also related to the fact
that ∇2 f (x) = H(x), which at the same time demonstrates that (Ω, h·, ·ix ) is a Hessian manifold;
see Section 4.3.
We now give an overview of the proofs of the set of results that establish existence, convergence
and complexity bounds on the Physarum dynamics (Theorems 1.1, 1.2, 1.3). We begin with
Theorem 1.1. It asserts that global solutions x : [0, ∞) → Ω of the Physarum dynamics indeed
4

exist. This theorem is nontrivial because by standard existence-uniqueness theorems we can only
obtain solutions defined on some tiny intervals around 0. To extend those solutions further we need
to show that they cannot leave the domain Ω in finite time. More precisely, if x : [0, T ) → Ω is a
solution, we need to show that the limit point x(T ) = limt→T x(t) exists and belongs to Ω. The
main concern is that potentially x(T ) may end up on the boundary of Ω, i.e., have a zero at some
coordinate. Of course, we expect some entries of x(t) to vanish with t → ∞, what we need to show
is that this happens in a controlled manner.
Let us take a look at the dynamics:
xi
p − ci ),
ẋi = qi − xi = (a⊤
ci i
where ai is the i-th column of A. If we could show that a⊤
i p is uniformly bounded then this
would imply that ẋi ≥ −αxi for some constant α, and, by the Gronwall lemma (see, e.g., [Per01])
xi (t) ≥ e−αt xi (0) > 0. Indeed Lemma 5.2, which shows that |a⊤
i p| is uniformly bounded over all
x > 0 such that Ax = b, rescues us. At first glance Lemma 5.2 may seem a bit technical, but in
a sense it captures exactly the property of the Riemannian space, which is responsible for the well
behavior of trajectories. This gives a proof of existence in the special case when the initial point
x(0) satisfies Ax(0) = b. (Assuming Ax(0) = b one can show that Ax(t) = b for every t, since ẋ is
tangent to this affine subspace.)4 Unfortunately, the case of infeasible start turns out to be harder.
An analogue of Lemma 5.2 does not hold if we let x varyPover Ω.
n
We proceed by analyzing a certain barrier function
i=1 ci yi ln xi , with y being any feasible
solution to (1). By analyzing its derivative, we are able to prove that on every finite interval [0, T )
this function is uniformly lower-bounded, which essentially means that on every finite interval there
exists some δ > 0 such that δ · y ≤ xi (t). This implies further that xi (t) is lower-bounded by a
positive constant, on every finite interval, but only for i ∈ supp(y). Of course we may repeat this
argument multiple times, with different y’s. This however is not enough: it is a common case for
linear programs that there is an index i ∈ {1, 2, . . . , n}, such that for every feasible solution y,
yi = 0. Therefore we need to use a different argument for the remaining indices.
To complete the argument we combine Lemma 5.2 with the above mentioned fact that on every
finite interval, xi (t) ≥ δy for some δ > 0 and some feasible solution y. Applying Lemma 5.2 to
y, we can essentially extract a uniform upper bound on |a⊤
i p(t)| over a fixed finite interval. This,
again by the Gronwall lemma, allows us to conclude positivity. In the process of the proof we need
to argue that kq(t)k and kx(t)k remain bounded uniformly over all t. For this, Lemma 5.2 is again
a starting point.
After establishing existence we study the limiting behavior of the solutions to (2). Our results
are summarized in Theorem 1.2. The first result claims that
c⊤ x(t) − c⊤ x⋆ ≤ R · e−νt ,
where R, ν > 0 are constants depending only on A, b, c, x(0) (we provide explicit bounds in Theorem 6.3). Furthermore, one can show that z(t) = x(t) − e−t x(0) satisfies Az(t) = b. In other words,
x(t) approaches feasibility with t → ∞.
To prove exponential convergence to the optimal value we start by showing that for every t there
is some y(t), feasible for (1), which is exponentially close to x(t). Towards this, we cannot directly
use z(t), because it may happen that zi (t) < 0 for some i. Instead, we use a well known fact that if
a point x violates constraints defining a polytope up to an additive error ε then its distance to this
polytope is at most N · ε where N depends only on the description size of the polytope (though
exponentially).
To proceed, we need to exploit the structure of the feasible region P = {x : Ax = b, x ≥ 0}. We
write P as a Minkowski sum H + C, where H is a polytope and C is a polyhedral cone. We extract
4For the case of feasible start, one can also give another proof of existence via Theorem 4.2.
5

vertices V from H and spanning vectors R from C. We study the decomposition of y(t) into
X
X
y(t) =
λv (t)v +
µr (t)r,
P

v∈V

r∈R

P
with λ(t), µ(t) ≥ 0 and v∈V λv (t) = 1. By studying potential functions of the kind ni=1 ci vi ln xi (t)
for a fixed vertex v ∈ V it can be shown that if v is non-optimal then λv (t) → 0 exponentially fast.
Similarly, for every r ∈ R we have µr (t) → 0. This is then used to deduce the result.
Let us now discuss the second result related to the asymptotic behavior of solutions to (2). It
says that if x : [0, ∞) → Ω is such a solution then x(t) has a limit x∞ with t → ∞. Furthermore
x∞ is an optimal solution to (1). The second part of this result follows easily once the first is
established. If the set of optimal solutions to (1) has only one element x⋆ , then it is easier to show
that x(t) → x⋆ , this basically follows from the fact that c⊤ x(t) tends to the optimal value and x(t)
approaches the feasible region. However if the optimal solution is not unique, then convergence of
x(t) is far from clear; x(t) could oscillate around the optimal set, without converging to a single
point. A proof that this does not happen gives us evidence that Physarum dynamics behaves nicely.
To prove convergence to a limit it is enough to show that
Z ∞
kẋ(t)k dt < ∞.
0

For this, it suffices to show kẋ(t)k = O(e−εt ) for some ε > 0. We will now focus on this task.
Recall that ẋi = xcii a⊤
i p − ci . Using similar tools as in the proof of existence one can show that
|a⊤
p(t)|
is
uniformly
bounded
over all t ∈ [0, ∞). Hence, if for some i ∈ {1, 2, . . . , n} we have
i
−εt
xi (t) = O(e ) then it follows easily that |ẋi (t)| = O(e−εt ). Denote the set of all such i’s by N and
its complement by J. We can then prove that J is the union of supports of all optimal solutions
to (1). This in turn implies that the subvector xJ (t) behaves in a sense more stably. We show
in particular, that for some optimal solution f , such that supp(f ) = J and some ε > 0, we have
εf ≤ x(t) uniformly over t ∈ [0, ∞). Furthermore the key Lemma 6.13 establishes some form of
⊤
−εt ), which concludes the proof.
continuity of a⊤
i p(t) for i ∈ J. It implies that |ai p(t) − ci | = O(e
The last theorem we discuss is Theorem 1.3. It shows that a discretization of the Physarum
dynamics is possible, provided that the step length is small enough. The crucial parameter which
controls the step length is A⊤ p ∞ ; if one is able to upper bound it over all possible x by some
−2 guarantees convergence. In the discretization
number Pmax then the step length of roughly Pmax
we choose the starting point x(0) to be feasible, by which x(k) will remain feasible for every k.
Therefore, by Corollary 5.3 one can take Pmax = D · Cs . The analysis of the discretization follows
along the same lines as was presented in [SV]. Let us discuss it briefly. One can show that c⊤ x(k) is
decreasing with k. However, we cannot guarantee a large drop of this value at every step. Instead
we introduce another potential function
def

B(k) =

n
X

ci x⋆i ln xi (k)

i=1

(where x⋆ is any optimal solution to (1)) and show that at every step when c⊤ x(k) − c⊤ x(k + 1)
is small, B(k) grows by a considerable amount. Moreover one can show that B(k) is uniformly
upper-bounded. By combining c⊤ x(k) and B(k) into a single potential function φ(k) we can easily
track the progress of our algorithm
Organization of The Paper. The remaining part of the paper is structured as follows. In
Section 3 we explain our notation and give a brief introduction to Riemannian manifolds. Section 4
discusses the optimization viewpoint on the Physarum dynamics. In Section 5 existence of solutions
to (2) is established. Section 6 is devoted to the study of limiting behavior of Physarum integral
curves, convergence results are established. The last Section 7 discusses the discretization of the
6

dynamics. Appendix A is an introduction to dynamical systems, at the end a general theorem is
proved, which is used in the main body to establish existence of solutions.
3. Preliminaries
Notation. In the paper we often use the convention that for a vector x ∈ Rn , the capitalized version
X denotes the diagonal matrix Diag (x1 , . . . , xn ). Whenever a scalar function f : R → R is applied to
def

a vector x ∈ Rn the component-wise application of f to x is meant, i.e. f (x) = (f (x1 ), . . . , f (xn ))⊤ .
We use the notation [n] to denote the set {1, 2,. . . , n}. The
 columns of the matrix A are denoted
x1
xn
by ai , i ∈ [n]. Whenever the objects W = Diag c1 , . . . , cn , L = AW A⊤ , p = L−1 b or q = W A⊤ p
appear, they should be understood as computed with respect to a point x, which will be clear from
the context. We will often refer to q(t) as to q computed w.r.t. x(t).
We use the standard norms for u ∈ Rn :
def P
• kuk1 = ni=1 |ui |,
1
def Pn
2 2 ,
• kuk2 =
i=1 ui
def

• kuk∞ = max (|u1 |, |u2 |, . . . , |un |) .
By supp(u) we denote the set {i ∈ [n] : ui 6= 0}. The linear program (1) is fixed for the whole
paper, we assume the feasible region {x : Ax = b, x ≥ 0} to be nonempty.
Parameters of Physarum. The following two parameters are often used in statements of quantitative bound regarding convergence time:
P
• Cs = ni=1 ci ,
• D = max{| det(A′ )| : A′ square submatrix of A}.
Riemannian Manifolds. All manifolds we deal with can be seen as open subsets of Euclidean
spaces, possibly embedded in Euclidean spaces of higher dimension. Some examples of those are
Rn , Rn>0 , open polyhedra {x ∈ Rn : Ax = b, x > 0} or a sphere {x ∈ Rn : kxk2 = 1}. If M is such
a manifold then at every point x ∈ M one can define a tangent space Tx M in a natural way. A
tangent space is equipped with a natural linear structure. For example, the tangent space to any
point x ∈ M = {x : Ax = b} is Tx M = {x : Ax = 0}. Given a manifold M we can endow it with a
Riemannian structure by assigning to every point x a scalar product h·, ·ix on Tx M in such a way
that h·, ·ix “varies smoothly with x”.5 When considering a function F : M → R on a Riemannian
manifold, the gradient of F is defined with respect to the local inner product. Hence in particular,
one can define gradient-descent algorithms based on a chosen Riemannian structure of the space,
those are typically very different to algorithms based on the standard Euclidean geometry.
4. Physarum as Steepest Descent on a Manifold
In this section we interpret the Physarum dynamics as steepest descent on a certain manifold.
We start by equipping the positive orthant with a Riemannian structure. In particular we formalize Theorem 1.4 in Theorem 4.1 and prove it. The entropy barrier interpretation is given in
Subsection 4.3
4.1. The Riemannian Structure and its Origin. We focus on two domains:
Ω = Rn>0

and Ωf = {x ∈ Ω : Ax = b} .

Note that Ω and Ωf are manifolds of dimension n and m respectively. We consider a Riemannian
structure h·, ·ix on Ω (which automatically induces a Riemannian structure on the submanifold
5One can make this statement precise, but our intention is rather to give some intuitions instead of precise
definitions.
7

Ωf ⊆ Ω) given as a family of positive definite matrices H(x) = CX −1 , for x ∈ Ω. This gives rise
to the following inner product on the tangent space at x:
def

hu, vix = u⊤ H(x)v.
We present two interpretations of where does this structure come from. The first one is related to
the following entropy-like function h : Ω → R
h(x) =

n
X
i=1

Note that:
(4)

∇h(x) = C ln(Cx)

ci xi ln(ci xi ) −

and

n
X

ci xi .

i=1

∇2 h(x) = CX −1

hence H(x) = ∇2 h(x) is an example of a so called Hessian Metric.
The second interpretation is based on the analysis of the following map: F : Ω → Ω
√
F (x) = 2 Cx.
We will now explain that F can be seen as an isometry between our Riemannian manifold (Ω, h·, ·ix )
and Ω with the standard inner product h·, ·i at every point. Recall that the Jacobian of F 6,
J(x) = (CX −1 )1/2 acts as a linear map J(x) : Tx (Ω) → Tx (Ω), i.e.:
1/2
J(x)(u) = J(x)u = CX −1
u.

Take any x and u, v from the tangent space at x, i.e. Tx (Ω) ∼
= Rn . We have:
D
E
hJ(x)(u), J(x)(v)i = (CX −1 )1/2 u, (CX −1 )1/2 v = u⊤ CX −1 v = u⊤ H(x)v = hu, vix .

Thus the inner product h·, ·ix can be seen as a pull-back of the standard inner product via the
map F . This also means that F is an isometry: it preserves angles and distances. Let us call
7
the space obtained by the
x from the x-space
√ transformation x 7→ F (x) the y-space . The point
1 −1 2
is represented by y = 2 Cx, in the y-space. The inverse mapping is x = 4 C y . The strictly
feasible set Ωf = {x ∈ Rn : Ax = b, x > 0} corresponds to:


n 1
−1 2
M = y ∈ R : AC y = b, y > 0
4
in the y-space.
4.2. Two Forces. We intend to show that the direction P (x) = W (A⊤ L−1 b − c) of the Physarum
dynamics (2) at the point x decomposes naturally into two directions
−1

def
(5)
(b − Ax)
Pf (x) = W A⊤ AW A⊤


−1

def
(6)
Po (x) = W A⊤ AW A⊤
Ax − c

with Pf (x) interpreted as the feasibility direction and Po (x) as the optimization direction. As one
can easily see P (x) = Pf (x)+Po (x). We will interpret both directions by studying the optimization
problem arising after transforming the linear program (1) into the y-space via the map F .
6Jacobian is the multidimensional analogue of a derivative. In the standard coordinate system it can be expressed

as a matrix of partial derivatives.
7Of course formally F (Ω) = Ω, hence the x-space and y-space coincide as sets, the difference however lies in the
geometry.
8

Feasibility. Fix any x̄ ∈ Ω and consider its image ȳ in the y-space. Recall that the strictly feasible
region of (1) in the y-space is


1
M = y > 0 : AC −1 y 2 = .
4
Let us assume ȳ ∈
/ M . In which direction do we need to move, to hit M as soon as possible? This
does not seem to be a simple problem, because it is non-convex. However, let us try to guess a y ⋆
which satisfies:
kȳ − y ⋆ k = min {kȳ − yk : y ∈ M }
and set h = y ⋆ − ȳ. Since the above minimization is done with respect to the Euclidean distance,
we expect the vector h to be orthogonal to the “surface” M at y ⋆ , equivalently to the tangent space
Ty⋆ M i.e. to every vector u which is tangent to the space M at the point y ⋆ ∈ M . It is easy to see
that:


1
−1
⋆
n
Ty⋆ M = u ∈ R : AC Y u = 0 .
2
Since h is just the orthogonal projection of y ⋆ − ȳ onto Ty⋆ M , we can obtain it by the well-known
formula8:
−1

AC −1 Y ⋆ (y ⋆ − ȳ)
h = Y ⋆ C −1 A⊤ A(Y ⋆ C −1 )2 A⊤
−1

(AC −1 (y ⋆ )2 − AC −1 Y ⋆ ȳ)
= Y ⋆ C −1 A⊤ A(Y ⋆ C −1 )2 A⊤
−1

(4b − AC −1 Y ⋆ ȳ).
= Y ⋆ C −1 A⊤ A(Y ⋆ C −1 )2 A⊤

We do not know what y ⋆ is, but drawing from the expectation that ȳ is close to feasibility, we can
approximate y ⋆ ≈ ȳ, hence obtaining:
−1


2
4b − AC −1 ȳ 2 .
h ≈ Ȳ C −1 A⊤ A Ȳ C −1 A⊤

√
Let us pull back h to the x-space9 by substituting ȳ = 2 C x̄ and multiplying h by the inverse1/2

= W 1/2 . We get:
Jacobian J −1 (x̄) = Diag x̄c11 , . . . , x̄cnn
−1

2W A⊤ AW A⊤
(b − Ax̄) = 2Pf (x̄)

which is the feasibility direction up to the a scaling factor 2.
Optimality. As in the previous discussion fix x̄ ∈ Ω, we will interpret the optimization direction at
x̄. We know that x̄ may not satisfy Ax̄ = b, but we still can write Ax̄ = bx̄ , where bx̄ = b − (b − Ax̄)
which we should imagine to be “close” to b. Our goal is to find the minimum value of c⊤ x over the
set {x : x ≥ 0, Ax = b}. Let us consider a related linear program:
(7)

min.
s.t.

c⊤ x
Ax = bx̄
x ≥ 0.

8The orthogonal projection onto the row space of a matrix B is given by B ⊤ (BB ⊤ )−1 B.
9More accurately: from the space T (Ω) to the space T (Ω).
ȳ
x̄
9

which is just minimization of c⊤ x over an affine subspace parallel to Ax = b which contains the
point x̄. Let us rewrite (7) in the y-space.
min.

(8)

s.t.

1 ⊤ 2
1 y
4
1
AC −1 y 2 = bx̄
4
x ≥ 0.

where 1 denotes the all-one vector. Let ȳ = F (x̄) and denote


1
x̄ def
−1 2
x̄
M = y : y > 0, AC y = b .
4

Note that ȳ ∈ M x̄ . Let us then compute the steepest descent direction at ȳ with respect to the
problem (8). We just need to compute the negative gradient of the objective at ȳ and project it
onto the tangent space Tȳ (M x̄ ). We have:


1
1 ⊤ 2
1 y = y.
∇y
4
2
Hence we need to project the vector − 21 ȳ onto the space:


−1
x̄
n 1
Tȳ (M ) = u ∈ R : AC Ȳ u = 0 .
2

As a result we obtain:





−1

1
−1 2 ⊤
−1
u = I − Ȳ C A A(Ȳ C ) A
AC Ȳ
− ȳ
2




−1
1
AC −1 ȳ 2 − ȳ .
Ȳ C −1 A⊤ A(Ȳ C −1 )2 A⊤
=
2
−1

⊤

To get a corresponding direction in the x-space, we compute:



−1
2
1
J −1 (x̄)(u) = W 1/2 Ȳ C −1 A⊤ A Ȳ C −1 A⊤
AC −1 ȳ 2 − ȳ
2


−1
√
1 1/2 1 1/2 ⊤ 
⊤
−1
AC 4C x̄ − 2 C x̄
W A AW A
= W
2
2


−1

⊤
⊤
Ax̄ − c = Po (x̄).
= W A AW A

Which is the optimization direction at x̄. Note that in case when x is feasible, the feasibility
direction Pf (x) is zero, hence the following theorem.

Theorem 4.1. Consider the manifold: Ωf = {x : Ax = b, x > 0} endowed with the Riemannian
metric hu, vix 
= u⊤ CX −1 v at every x ∈ Ωf . Then, the Physarum direction P (x) =


−1
W A⊤ AW A⊤
b − c is the gradient of the objective function c⊤ x with respect to the Rie
mannian manifold Ωf , h·, ·ix .
Proof. One can use the above derivation of the optimization direction to deduce this result. If x is
feasible then:




−1
−1


⊤
⊤
⊤
⊤
b − c = P (x).
Ax − c = W A AW A
Po (x) = W A AW A
10

4.3. Physarum as an Entropy Barrier Method. In this section we show that every Physarum
trajectory starting from a strictly feasible point can be seen as a path of optimizers to a certain
family of convex programs. Let us fix any starting point s > 0 satisfying As = b and consider for
every µ ≥ 0 the following convex program:
(9)

min.

µc⊤ x +

s.t.

Ax = b
x ≥ 0.

n
X
i=1

xi ci ln(xi ci ) −

n
X

xi ci (1 + ln(si ci ))

i=1

Note that one can rewrite the objective function equivalently as:
1
c⊤ x + f s (x).
µ
Which is just our original linear objective, regularized by an entropy-like strongly convex function
def

f s (x) =

n
X
i=1

∇2 f s (x)

CX −1

xi ci ln(xi ci ) −

n
X

xi ci (1 + ln(si ci )).

i=1

Observe that
=
is the Riemannian metric we are studying above. It should be
intuitively clear that as µ → ∞ the solution to (9) will tend to the optimal solution of the linear
program (1).
Theorem 4.2. Take any starting point s > 0 with As = b. Suppose that x(µ) is the unique optimal
solution to (9). Then x(µ) for µ ∈ [0, ∞) is the solution to the Physarum dynamical system with
the initial condition x(0) = s.
The idea of the proof is to write down KKT conditions for the convex program (9) and take
advantage of strong duality. By implicit differentiation we conclude that x(µ) follows the same
dynamics as Physarum, furthermore x(0) = s, hence they have to coincide.
Proof. To improve readability, we will assume that c = 1, i.e. ci = 1 for all i = 1, . . . , n. To obtain
the general version one can go back and forth with the substitution
x 7→ C −1
P
Pz.
s
Fix any s > 0 such that As = b. Let us denote f (x) = i xi ln xi − i xi (1 + ln si ). First
one needs to show that in fact for every µ ≥ 0 there exists a minimizer of µ1⊤ x + f s (x) in the set
{x : Ax = b, x > 0}, in other words that x(µ) > 0. This can be seen as follows: pick any point x̄
on the boundary, i.e. having x̄i = 0 for some i. Consider the following scalar function:
h(p) = µ1⊤ (x̄ + p(s − x̄)) + f s (x̄ + p(s − x̄)).

for p ∈ [0, 1]; h(p) is continuous on [0, 1] and differentiable in the interval (0, 1). We have:
d
h(p) = µ1⊤ (s − x̄) + (s − x̄)⊤ (ln(x̄ + p(s − x̄)) − ln s).
dp

This implies that h′ (p) → −∞ as p → 0, by looking at a coordinate i where x̄i = 0. Hence
h(ε) < h(0) for ε > 0 small enough, implying (together with convexity) that the optimal value of
µ1⊤ x + f s (x) over {x : Ax = b, x ≥ 0} is attained at some point x(µ) > 0.
Fix µ ≥ 0, introduce dual variables y ∈ Rm for the equality constraints and consider the Lagrangian:
X
X
L(x, y) = µ1⊤ x +
xi ln xi −
xi (1 + ln si ) − y ⊤ (Ax − b).
i

We always keep in mind that x > 0. Fix y ∈

i

Rm

and let us compute the derivative of L(x, y):

∇x L(x, y) = µ1 + (ln x − ln s) − A⊤ y.
11

the derivative is 0 at a point x given by:
ln xi = a⊤
i y − µ + ln si

xi = si · exp(a⊤
i y − µ)

Then the dual objective g(y) is given by substituting the above x (at which L(x, y) is minimized)
in L(x, y). After some cancellations we get:
X
X
g(y) = y ⊤ b −
xi = y ⊤ b −
si · exp(a⊤
i y − µ)
i

i

Hence the dual program is:
X

min

y⊤b −

s.t.

y ∈ Rm

i

si · exp(a⊤
i y − µ)

Note that strong duality holds, since Slater condition is satisfied for the primal (the witness being
s). This means that the unique maximizer of g(y) will also yield the optimal solution for the primal.
Let
X
def
m
y(µ) = argmax{y ⊤ b −
si · exp(a⊤
i y − µ) : y ∈ R },
i

def

xi (µ) = si · exp(a⊤
i y(µ) − µ).

Then (x(µ), y(µ)) is the optimal primal-dual pair.
We will show that x(µ) is a solution to the dynamical system:


−1

d
⊤
⊤
AX1 − 1 ,
x(µ) = P (x(µ)) = X A AXA
dµ
which represents Physarum when c = 1. Let us start by examining the initial condition. If µ = 0
then one can easily check, that the optimal pair is (x(0), y(0)) = (s, 0), as claimed. Let us now
compute the derivative of x(µ) w.r.t. µ, note that the implicit function theorem guarantees that
(x(µ), y(µ)) is a continuously differentiable function of µ. In the following calculations we write
shortly x, y for x(µ) and y(µ). We have:
⊤
⊤
ẋi = si · exp(a⊤
i y − µ) · (ai ẏ − 1) = xi (ai ẏ − 1).

In other words ẋ = X(A⊤ ẏ − 1). Let us now take the equality Ax = b and differentiate w.r.t. µ:
Aẋ = 0
AX(A⊤ ẏ − 1) = 0
In consequence, ẏ satisfies (AXA⊤ )ẏ = AX1. Since rank(A) = m, one can show that AXA⊤ is
invertible, hence ẏ is uniquely determined by ẏ = (AXA⊤ )−1 AX1. Consequently:


−1



AX1 − 1 = P (x).
ẋ = X A⊤ ẏ − 1 = X A⊤ AXA⊤
12

5. Existence of Solution
In this section we would like to argue that no matter what the initial condition x(0) > 0 is, the
Physarum dynamics has a global solution x : [0, ∞) → Rn>0 . Let us start by explaining what is
the significance of the assumption x(t) > 0. The Physarum dynamics is defined as a dynamical
system ẋ = P (x), where P (x) = W (A⊤ (AW A⊤ )−1 b−c). For x > 0 one can show that (AW A⊤ )−1 b
exists, (since the kernels of AW A⊤ and A⊤ coincide) hence P (x) is well defined and continuous
over Ω = Rn>0 . Can we extend it out of Ω? This leads to troubles, because in general AW A⊤ might
be singular, when w contains some zeros or negative entries. For example if w = 0 then clearly
AW A⊤ = 0, hence there is no hope to extend the vector field to x = 0. In general there is no
simple way to extend (continuously) this vector field from Ω = Rn>0 to a larger set.
5.1. Basic Results. This section introduces some preparatory results which are essential for studying the existence of solution and then asymptotic properties of the Physarum dynamics. Recall
that we are always assuming that the feasible region {x : Ax = b, x ≥ 0} is nonempty.
Remark 5.1. In the proofs we will repeatedly use the fact that polyhedra of the kind {x : Ax =
b, x ≥ 0} always have at least one vertex. Let x̄ be such a vertex. Denote Z = {j ∈ [n] : x̄j = 0}
and N = [n] \ Z. Then x̄N (i.e. the subvector of x̄ consisting of entries xj for j ∈ N ) can be
determined as a unique solution to the linear system AN xN = b, where AN is a submatrix of A
consisting of columns aj for j ∈ N . The above are standard facts about vertices of polyhedra and
can be found in any book on linear programming, e.g. [Kar91].
Below we present a key lemma, which then allows to get a solid grasp on the behavior of the
vector field P (x) defining Physarum dynamics.
Lemma 5.2. Let w ∈ Rn>0 , W = Diag (w), L = AW A⊤ . There exists a constant α > 0 depending
only on A such that for every i ∈ [n] :
α
≤
A⊤ L−1 ai
.
wi
∞
Quantitatively, one can take α = D = max{| det(A′ )| : A′ is a square submatrix of A}.
Proof. Fix i and denote by p the solution to the system Lp = ai . We can assume that p⊤ aj ≥ 0 for
every j ∈ [m] by replacing the row aj by −aj if necessary. One can easily see that such a change
does not alter the problem, because L remains the same.
1
−1
Let us first show that a⊤
i L ai ≤ wi . Note that
wi ai a⊤
i 

m
X

wj aj a⊤
j =L

j=1

⊤
n
where  is the PSD order. This means that wi u⊤ ai a⊤
i u ≤ u Lu, for every u ∈ R . Let us pick
−1
u = L ai . We get:
⊤
wi u⊤ ai a⊤
i u ≤ u Lu

−1
⊤ −1
⊤ −1
−1
wi a⊤
i L ai ai L ai ≤ ai L LL ai
−1
2
⊤ −1
wi (a⊤
i L ai ) ≤ ai L ai
1
−1
a⊤
.
i L ai ≤
wi

It remains to argue that for some α and for every k ∈ [n] :
⊤
a⊤
k p ≤ αai p.
13

⊤
Fix k ∈ [n], assume k 6= i. If a⊤
k p = 0 then we are done, assume ak p > 0. From Lp = ai we get:
m
X

wj aj (a⊤
j p) = ai .

j=1

def

⊤
Hence the set Sk = {s ∈ Rm
Take s ∈ Sk
≥0 : As = ai ∧ sk > 0} is nonempty ( W A p belongs to it).
Pm
with sk maximum possible (it can be seen that sk is bounded over all s ∈ Sk ). Then j=1 sj aj = ai ,
P
⊤
⊤
⊤
⊤
⊤
hence m
j=1 sj aj p = ai p. Since sj aj p ≥ 0 for all j, we can deduce that sk ak p ≤ ai p and hence
a⊤ p

⊤
i
a⊤
k p ≤ sk = αk ai p. It is enough to choose α = maxk αk .
For the quantitative bound one needs to note that α is chosen according to the following values:
εk = max{sk : As = ai , s ≥ 0}. In fact α = maxk ε1k . Because linear programs attain optimal
values in vertices, one can argue that s⋆ – the optimal solution to max{sk : As = ai , s ≥ 0} (for
some fixed k) can be chosen to be a vertex of the polyhedron {s : As = ai , s ≥ 0}. By the Cramer’s
rule, every positive entry of s⋆ is lower-bounded by D −1 .

Let us now see what happens when w = C −1 x actually comes from a feasible vector.

Corollary 5.3. Suppose that x > 0 and Ax = b, we have:
A⊤ p

∞

≤ D · Cs .

Proof. Let w = C −1 x, using Lemma 5.2 we get:
A⊤ p

∞

= A⊤ L−1 b
≤
≤

n
X

i=1
n
X

∞

n
X

=

i=1

xi A⊤ L−1 ai
xi

i=1

A⊤ L−1 (ai xi )
∞

∞

D
= D · Cs .
wi

It turns out that if we could prove a result such as above for all x ∈ Ω, not only for feasible x,
it would imply the existence of solution. Unfortunately, this ceases to be true when do not make
this restriction. Intuitively, our strategy is to prove that A⊤ p is bounded over every trajectory
{x(t) : 0 ≤ t ≤ T }, this then allows us to reason that no trajectory will leave Ω. The next lemma
shows a uniform bound but with kA⊤pk replaced by W A⊤ p .

Lemma 5.4. Suppose y if feasible: Ay = b and y ≥ 0. Then for every x ∈ Rn>0 it holds that
kqk∞ ≤ α kyk1 . As in Lemma 5.2 we can take α = D.
Pn
−1
Proof. Recall that q = W A⊤ L−1 b, hence qi = wi a⊤
i L b. We express b as b = Ay =
j=1 yj aj .


n
n
X
X
−1 
−1

|qi | = wi a⊤
y
a
L
yj wi a⊤
=
j j
i
i L aj
j=1

≤

n
X

j=1

−1
wi yj a⊤
i L aj =

−1
wi yj a⊤
j L ai

j=1

j=1

Lemma 5.2
≤

n
X

n
X
j=1

wi yj

α
= α kyk1
wi
14

Hence kqk∞ ≤ α kyk1 .
The following corollary gives a concrete bound we can deduce from the above lemma.
Corollary 5.5. For every x ∈ Rn>0 it holds that kqk∞ ≤ D 2 · n · kbk1 .
Proof. We use Lemma 5.4. Let us pick y as any vertex of the feasible region {x : Ax = b, x ≥ 0}.
We know that the non-zero portion yN of y can be determined as a solution to the linear system
AN yN = b. Hence, by Cramer’s rule every non-zero entry of y ca be obtained as follows:
yi =

m
X

bj

j=1

αj
.
α0

for some α0 , α1 , . . . , αm being determinants of square submatrices of A. Since A has integer entries,
|α0 | ≥ 1 and |α1 |, . . . , |αm | ≤ D, thus |yi | ≤ kbk1 D and kyk1 ≤ n kbk1 D.
5.2. The Proof of Existence. We are now ready to prove the following theorem (same as Theorem
1.1). As previously, we will use Ω to denote Rn>0 .
Theorem 5.6. For every initial condition x(0) ∈ Ω the dynamical system (2) has a unique solution
x : [0, ∞) → Ω.
From now on we assume that the initial condition x(0) = s ∈ Ω is fixed. In the light of Theorem A.2,
to prove 5.6 it suffices to show the following two lemmas.
Lemma 5.7. Suppose x : [0, T ) → Ω is a solution to (2) for some T ∈ R>0 then limt→T − x(t)
exists.
def

Lemma 5.8. Suppose x : [0, T ) → Ω is a solution to (2) for some T ∈ R>0 . If x(T ) = limt→T − x(t)
exists then x(T ) ∈ Ω, i.e. x(T ) is strictly positive.
Let us start by proving that every trajectory stays in a bounded region.
Lemma 5.9. Suppose x : [0, T ) → Ω is a solution to (2), then for every i ∈ [n] and t ∈ [0, T ):
where β = D 2 · n · kbk1 .

xi (t) ≤ max (xi (0), β) ,

Proof. Fix i ∈ [n], by Corollary 5.5 we have:
ẋi = qi − xi ≤ β − xi .
Applying Gronwall lemma to y(t) = xi (t) − β yields:

xi (t) ≤ (1 − e−t )β + e−t xi (0) ≤ max (xi (0), β) .

We are now ready to prove Lemma 5.7:
Proof of Lemma 5.7: Suppose x : [0, T ) → Ω is a solution to (2). Fix i ∈ [n], we show that
xi (t) has a limit with t → T − . We know that xi is differentiable in the interval [0, T ) and |ẋi | is
bounded, since
|ẋi (t)| = |qi (t) − xi (t)| ≤ |qi (t)| + |xi (t)| ≤ β + max(xi (0), β)
by Lemma 5.4 and Corollary 5.9. Hence xi (t) is a Lipschitz function in the interval [0, T ). This
implies that limt→T − xi (t) exists.
15

Now we are left with the task of proving that no coordinate of x will approach 0 as t → T − . This
task simplifies significantly (as implied by the next lemma) if we assume that there is a strictly
feasible solution to (1), i.e. y ∈ Rn such that Ay = b and y > 0. Nevertheless, we do not want to
make any such assumptions to make our proof valid in full generality. We proceed by showing that
positivity holds for every i ∈ supp(y), for any feasible y, here supp(y) is the support of the vector
y, i.e. the set {j ∈ [n] : yj > 0}.
Lemma 5.10. Suppose x : [0, T ) → Ω is a solution to (2) and y is any feasible solution to (1). If
x(T ) = limt→T − x(t) then xi (T ) > 0 for every i ∈ supp(y).
Proof. Fix any feasible solution y to (1). To justify the claim, we use the following “barrier function”:
n
X
def
yj cj ln xj (t).
B(t) =
j=1

B(t) is clearly well defined on [0, T ). Suppose for the sake of contradiction that inf t∈[0,T ) xi (t) = 0
for some i ∈ supp(y). We know by Corollary 5.9 that xj (t) are uniformly upper-bounded over
t ∈ [0, T ). Hence we know that yj cj ln xj (t) ≤ M for all j ∈ [n], t ∈ [0, T ) and some constant
M ∈ R. It follows that inf t∈[0,T ) B(t) = −∞. This implies that the derivative Ḃ(t) is unbounded
from below, for t ∈ [0, T ).
However, let us compute this derivative:

(10)

Ḃ(t) =

We analyze the term
n
X
j=1

yj cj

n
X
j=1

Pn

n

n

n

j=1

j=1

j=1

X
qj (t)
qj (t) X
x˙j (t) X
yj cj
yj cj
yj cj =
=
−
− c⊤ y.
yj cj
xj (t)
xj (t)
xj (t)

qj (t)
j=1 yj cj xj (t) .

n
n
X
X
wj a⊤
qj (t)
j p(t)
⊤
⊤
yj cj
y j a⊤
=
=
j p(t) = (Ay) p(t) = b p(t).
xj (t)
xj (t)
j=1

j=1

Note that b⊤ p = b⊤ L−1 b ≥ 0, hence the above sum is nonnegative. In consequence:
Ḃ(t) ≥ −c⊤ y = const.

which is a contradiction, since this quantity was claimed to be unbounded from below.
We have now all necessary tools to prove Lemma 5.8.
Proof of Lemma 5.8: Let x : [0, T ) → Ω be a solution to (2). Fix y to be any feasible solution to
(1). Lemma 5.10 implies that there exists ε > 0 such that ε · y ≤ x(t) for every t ∈ [0, T ). Indeed:
for i such that yi > 0 we have inf t∈[0,T ) xi (t) > 0, so εi · yi ≤ xi (t) for some εi > 0 and all t ∈ [0, T ),
then take ε = min{εi : i ∈ supp(y)}.
Fix any t ∈ [0, T ) and consider w ∈ Rn defined as wi = xic(t)
. Let us bound the norm of A⊤ L−1 b:
i
A⊤ L−1 b

∞



= A⊤ L−1 
≤
≤

n
X

j=1
n
X
j=1

n
X
j=1

yj aj 

yj A⊤ L−1 aj

xj (t) α
α
=
ε wj
ε
16



∞

∞

n
Lemma 5.2 X
α
yj
≤
wj

n
X
j=1

j=1

cj .

Let us denote M =

α
ε

Pn

Pick now any i ∈ [n], we have:


xi
M
ẋi = qi − xi = pi − xi ≥ xi −
−1 .
ci
ci

j=1 cj .

Hence, from the Gronwall Lemma we obtain:


−1
t −M
c

xi (t) ≥ xi (0)e

i

.

In particular xi (t) is lower bounded by a positive constant over the whole interval [0, T ). The
lemma follows.
We can now conclude Theorem 5.6 easily.
Proof of Theorem 5.6: We use Theorem A.2. One needs to observe that the function P (x) =
q − x is in fact of class C 1 . Even more is true: P (x) is a rational function in the region Ω, this
can be seen using the Cramer’s rule. The remaining assumptions are satisfied by Lemmas 5.7and
5.8.
6. Asymptotic Behavior of Physarum Trajectories
In the previous section we have established the existence of solutions to the Physarum dynamics.
We know that for every starting point s > 0 there is a solution x : [0, ∞) → Rn>0 to the Physarum
dynamics with x(0) = s. We would like to study the behavior of x(t) when t → ∞. Let us fix the
starting point x(0) and pick Mx > 0 such that:
Mx−1 ≤ xi (0) ≤ Mx

for every i ∈ [n].

Furthermore, pick x⋆ to be any optimal basic solution (i.e. x⋆ is a vertex of the feasible region) to (1).
Our first result will be that c⊤ x(t) approaches c⊤ x⋆ – the optimal value of the linear program (1).
Before we proceed with the proof let us establish some useful properties of the feasible region.
6.1. The Feasible Region. Let us consider P = {x : Ax = b, x ≥ 0} – the feasible region of (1).
Note that P does not contain a line, so it can be expressed as the Minkowski sum P = H + K,
where H is the convex hull of vertices of P and K is a polyhedral cone
K = {r ∈ Rn : Ar = 0, r ≥ 0} .
Let us denote by V the set of vertices of P , so H = conv(V ). Further define R to be the set of
vertices of the polytope
)
(
n
X
ri = 1, r ≥ 0 .
r : Ar = 0,
i=1

Then K is the conic hull of R:

K=

(

X

r∈R

Let us first establish the following bounds:

)

µr r : µr ≥ 0 .

Lemma 6.1. Let V and R, be as defined above.
(1) Let v ∈ V and i ∈ [n] such that vi 6= 0, then D −1 ≤ |vi | ≤ D · kbk1 .
(2) Let r ∈ R and i ∈ [n] such that ri 6= 0, then D −1 ≤ |ri | ≤ D.
Proof. This follows immediately from the Cramer’s rule and the characterization of vertices of
polyhedra (see Remark 5.1).
17

Let x ∈ P , since P = H + K, x can be expressed as:
X
X
x=
λv v +
µr r
P

v∈V

r∈R

with λv , µr ≥ 0 for v ∈ V, r ∈ R and v∈V λv = 1. Note that the sets V and R might be very
large (of size exponential in n). However, by Caratheodory’s theorem we may pick the vectors of
coefficients λ, µ so that |supp(λ)| ≤ n + 1 and |supp(µ)| ≤ n + 1.
The last preliminary fact we would like to state is the following proposition on sensitivity analysis
of linear programs. It is very useful to estimate the distance from a polyhedron to a point which
satisfies its constraints with some additive error.
Proposition 6.2 (Sensitivity analysis). Suppose that B ∈ ZM ×N , g ∈ RN and b′ , b′′ ∈ RM are
such that both linear programs: min{g⊤ x : Bx ≤ b′ } and min{g⊤ x : Bx ≤ b′′ } have finite optimal
values, then:
min{g⊤ x : Bx ≤ b′ } − min{g⊤ x : Bx ≤ b′′ } ≤ N DB kgk1 · b′ − b′′

∞

.

Where DB = max{| det(B ′ )| : B ′ a square submatrix of B}.
For a proof of the above proposition we refer to [Sch86].
6.2. Convergence to the Optimal Value. In this section we establish the following theorem.
Theorem 6.3. Suppose that x : [0, ∞) → Ω is any solution to the Physarum dynamics. Then, for
some R, ν > 0 depending only on A, b, c, x(0), we have:
c⊤ x(t) − c⊤ x⋆ ≤ R · e−νt .
Quantitatively, one can take ν = D −3 and R = exp(8D 2 · Cs · kbk1 ) · (n + Mx )2 .

Let us begin by splitting the set V into two subsets VO = {v ∈ V : c⊤ v = c⊤ x⋆ } and VN = V \VO .
Then VO is basically the set of optimal vertices. The following bound will be very useful in the
proof of Theorem 6.3.
Lemma 6.4. Let v ∈ VN , then c⊤ v − c⊤ x⋆ ≥ D −2 .

Proof. Since v is a vertex of P , by Cramer’s rule it can be expressed as v = ( zd1 , zd2 , . . . , zdn )⊤ .
Where d, z1 , z2 , . . . , zn ∈ Z and 1 ≤ d ≤ D. Hence c⊤ v is a number of the form dz for some integer
′
z. Similarly c⊤ x⋆ = dz′ with z ′ , d′ ∈ Z and 1 ≤ d′ ≤ D. We get:
c⊤ v − c⊤ x⋆ =

1
d′ z − z ′ d
≥ 2.
′
dd
D

The next lemma says that x(t) becomes exponentially close to the feasible region P as t → ∞.
Lemma 6.5. For every t ≥ 0 there exists y(t) ∈ P such that:

kx(t) − y(t)k∞ < e−t · 3n2 · D · Mx .

Proof. Start with the Physarum dynamics ẋ = q − x and multiply both sides by et . We get
d
t
t
dt (x(t)e ) = e q(t). Take integrals of both sides to obtain:
Z t
es−t
−t
−t
x(t) = x(0)e + (1 − e )
q(s)
ds.
1 − e−t
0
Rt
es−t
Observe that z(t) = (1−e−t ) 0 q(s) 1−e
−t ds is a convex combination of q(s) (which satisfy Aq(s) = b
for every s), hence we conclude that z(t) = x(t) − x(0)e−t satisfies Az(t) = b. Note also that
z(t) ≥ x(0)e−t , so z(t) approaches P with t → ∞. We would like to estimate the distance from
18

z(t) to the closest point in P . As a tool for that we would like to use Proposition 6.2. We will use
variables y ∈ Rn and a new variable d ∈ R. Let us write down two linear programs:
min.
s.t.

min.
s.t.

d
Ay = b

d
Ay = b
ky − z(t)k∞ ≤ d

ky − z(t)k∞ ≤ d
y≥0

y ≥ −x(0)e−t .

The minimum value of the first linear program is equal to dist∞ (z(t), P ): the minimum ℓ∞ –distance
from z(t) to some point y ∈ P . The second one has optimal value 0, which is demonstrated by
taking y = z(t). To apply Proposition 6.2 we need too transform the constraints into inequality
form, but this is done in a standard manner. We obtain:
dist∞ (z(t), P ) ≤ (n + 1) · (2n) · D · x(0)e−t

∞

.

In the above, (2n) · D is the bound on DB where B is the matrix obtained after transforming
the constraints into the inequality form. By taking y(t) which attains this minimum distance, we
obtain the desired point in P :
kx(t) − y(t)k∞ ≤ kx(t) − z(t)k∞ + kz(t) − y(t)k∞
≤ x(0)e−t
2

∞

+ (n + 1) · (2n) · D · x(0)e−t
−t

≤ 3n · D · x(0)e

∞

∞

.

By the discussion at the beginning of the section we know that for every t, y(t) (from the above
lemma) can be decomposed as:
(11)

y(t) =

X

λv (t)v +

v∈V

X

µr (t)r

r∈R

P
for some λv (t), µr (t) ≥ 0 such that v∈V λv (t) = 1. Of course this decomposition is not unique.
We choose it arbitrarily, but as noted previously it can be done in such a way that most of the
coefficients are zeros, i.e. |supp(λ(t))| ≤ n + 1 and |supp(µ(t))| ≤ n + 1. Our strategy is to show
that µr (t) → 0 for all r ∈ R and λv (t) → 0 for all v ∈ VN . This will imply that y(t) (hence also
x(t)) tends to the optimal region. Towards this let us prove a preparatory result.
Lemma 6.6. The following bounds hold with some Q, ν > 0 depending only on A, b, c, x(0):
(1) For every r ∈ R: min{xi (t) : i ∈ supp(r)} ≤ Q · exp(−νt).
(2) For every v ∈ VN : min{xi (t) : i ∈ supp(v)} ≤ Q · exp(−νt).

Quantitatively, one can take ν = D −3 and Q = exp 4D 2 Cs · kbk1 · (n + Mx ).

Proof. Let us start with the second part. Pick any v ∈ VN and x⋆ ∈ VO , consider the following
potential function:
def

f (t) =

n
X
i=1

ci vi ln xi (t) −
19

n
X
i=1

ci x⋆i ln xi (t)

and take its derivative:
n
n
X
ẋi (t) X ⋆ ẋi (t)
d
ci vi
ci xi
f (t) =
−
dt
xi (t)
xi (t)
i=1
i=1
 ⊤
 X
 ⊤

n
n
X
ai p(t)
ai p(t)
ci vi
=
ci x⋆i
−1 −
−1
ci
ci
i=1

=

i=1

Ax⋆

Observe that, since Av =
n
X
i=1

i=1

n
X



⋆
⊤
⊤ ⋆
⊤
(vi − xi ) ai p(t) + c x − c v .

= b, we have:

⋆
⊤
⊤
⊤
(vi − x⋆i ) a⊤
i p(t) = (v − x )A p(t) = (b − b )p(t) = 0.

Hence, the derivate takes a particularly simple form:
d
f (t) = c⊤ x⋆ − c⊤ v.
dt
By Lemma 6.4 we know that

d
dt f (t)

≤ −ε for ε = D −2 . This implies the following bound on f (t):
f (t) ≤ f (0) − ε · t.

Hence

n
X
i=1

ci vi ln xi (t) ≤ f (0) +

n
X
i=1

ci x⋆i ln xi (t) − ε · t.

Let us use Corollary 5.9 and Lemma 6.1 to bound:
n
X
i=1

ci x⋆i ln xi (t) ≤

n
X
i=1

ci D · kbk1 ln(max(D 2 n kbk1 , Mx )) ≤ Cs D · kbk1 · ln(D 2 n kbk1 + Mx ).

Let us call M the whole expression on the right-hand side. Similarly we can bound f (0) ≤ 2M .
We get:
n
X
ci vi ln xi (t) ≤ 3M − εt.
i=1

Suppose wlog that x1 (t) = min{xi (t) : i ∈ supp(v)}, then:
ln x1 (t) ·

It remains to bound

P

X

i∈supp(v)

i∈supp(v) ci vi

ci vi ≤

n
X
i=1

ci vi ln xi (t) ≤ 3M − εt.

≥ D −1 to get:

min{xi (t) : i ∈ supp(v)} ≤ exp(3M D) · exp(−D −3 t).
P
The proof of the first part is analogous, but uses the potential function ni=1 ci ri ln xi (t).

We are now ready to present a complete proof of the Theorem:
Proof of Theorem 6.3: Let us first pick r ∈ R, we will show that µr (t) → 0 exponentially fast.
We have for every i ∈ supp(r):
µr (t)ri ≤ yi (t) ≤ xi (t) + M e−t .

with M as in Lemma 6.5. Hence, by Lemma 6.6, there is some i ∈ supp(r) such that:
µr (t)ri ≤ Qe−νt + M e−t .
20

By the bound ri ≥ D −1 from Lemma 6.1 we obtain:

µr (t) ≤ DQe−νt + DM e−t ≤ 2DQe−νt .

By a similar argument, we get for every v ∈ VN :

λv (t) ≤ 2DQe−νt .

Having this, let us consider the difference: c⊤ y(t) − c⊤ x⋆ . We know that:


X
X
X
c⊤ 
λv (t)c⊤ v = (c⊤ x⋆ )
λv (t)v  =
λv (t).
v∈VO

v∈VO

v∈VO

Hence we get a bound:



c⊤ y(t) − c⊤ x⋆ ≤ c⊤ 

X

v∈VN

λv (t)v +

X

r∈R



µr (t)r  .

Recall that λ(t) and µ(t) were chosen in such a way that all but at most n + 1 of their entries are
zero. Hence:
X

≤

λv (t)v

v∈VN

A similar bound holds for

P

∞

X

v∈VN

r∈R µr (t)r.

λv (t) kvk∞ ≤ (n + 1) · 2DQe−νt · D · kbk1 .
Altogether, we obtain the following bound:

c⊤ y(t) − c⊤ x⋆ ≤ 4(n + 1)Cs D 2 kbk1 Qe−νt .

To conclude let us relate c⊤ y(t) to c⊤ x(t). We have, by Lemma 6.5:
c⊤ y(t) − c⊤ x(t) ≤ Cs · M e−t .
This finally yields the claimed bound:
c⊤ x(t) − c⊤ x⋆ ≤ 4(n + 1)Cs D 2 kbk1 Qe−νt + Cs · M e−t ≤ Q2 e−νt .

6.3. Convergence to a Limit. By the previous subsection we know that for every solution x :
[0, ∞) → Ω, c⊤ x(t) − c⊤ x⋆ → 0 exponentially fast with t → ∞. One can use this fact to prove
that whenever (1) has a unique solution then x(t) actually has a limit x∞ and x∞ is this unique
optimal solution to (1). However, if we do not assume uniqueness the fact that x(t) has a limit is
no more obvious. We are aiming to prove it in the current subsection. Let us now formally state
the theorem.
Theorem 6.7. Suppose x : [0, ∞) → Ω is a solution to the Physarum dynamics (2). Then there
exists a limit x∞ = limt→∞ x(t). Furthermore x∞ is an optimal solution to the linear program (1).
In this subsection the norm symbol k·k should be understood as the infinity norm k·k∞ . Unlike
in the previous proofs we will not keep track of constants and hide them under the big O notation.
By a constant we mean any quantity which solely depends on A, b, c, n, m and the fixed starting
point x(0) under consideration.
Let us start by giving a simple lemma which provides a sufficient condition for existence of a
limit.
21

Lemma 6.8. Suppose x : [0, ∞) → Ω is a differentiable function. If the integral:
Z ∞
kẋ(t)k dt
0

is finite, then there exists a limit

x∞

Proof. Observe that x(t) = x(0) +

Rt

= limt→∞ x(t).

ẋ(s)ds. By going with t → ∞, we obtain:
Z ∞
ẋ(s)ds.
lim x(t) = x(0) +
0

t→∞

0

Since the integral on the right-hand side exists (because it is absolutely convergent by the assumption), we conclude existence of the limit.
In our case ẋ = q − x. Our proof strategy is to show that kq(t) − x(t)k = O(e−εt ) for some
constant ε > 0, then the assumption of Lemma 6.8 is satisfied.
Let us now review what do we know from the proof of convergence to the optimal value. The
following fact can be extracted from the proof of Theorem 6.3:
Lemma 6.9. For every t ≥ 0 there exists x⋆ (t) – an optimal solution to (1) such that kx(t) − x⋆ (t)k =
O(e−εt ), for some ε > 0.
Note that the above lemma does not imply that x(t) in fact has a limit. It could potentially happen
that x(t) oscillates over the optimal region without approaching a single point. To prove it we need
to understand the behavior of Physarum with respect to the structure of the optimal set.
From the general theory of linear programming (see e.g. [Wri97]) we know that [n] can be
decomposed into [n] = J ∪ N with J ∩ N = ∅ such that every optimal solution to (1) is supported
on J and every feasible solution supported on J is optimal. We will denote by xJ the part of a
vector x corresponding to indices j ∈ J, similarly xN . Let us now establish two important lemmas,
which we can prove using techniques developed so far.
Lemma 6.10. If N is the above define set of indices then kxN (t)k = O(e−εt ) for some ε > 0.

Proof. It follows immediately from 6.9. For every t, x⋆N (t) = 0, because x⋆ (t) is optimal, hence the
claim.
Lemma 6.11. There exists a constant ε > 0 such that xj (t) > ε for every j ∈ J and t ∈ [0, ∞).
Proof. Recall that the set of optimal vertices is denoted by VO . Note that for every j ∈ J, there
is a vertex v ∈ VO such that vj = Ω(1). Hence, it is enough to show that for every v ∈ VO , the
function:
n
X
ci vi ln xi (t)
fv (t) =
i=1

is bounded from below by a universal constant (this is true because of the fact that kx(t)k is
uniformly bounded, by Corollary 5.9). Note that, by a calculation analogous to that performed in
the proof of 6.6, we obtain:
d
(fv (t) − fu (t)) = 0
dt
for every v, u ∈ VO , which implies that all fv differ only by a constant. Hence either all of them
are lower-bounded by a universal constant or none of them.
By Lemma 6.9 we know that for every t, x(t) is exponentially close to some x⋆ (t) ∈ conv(VO ).
This implies in particular that for big enough t and for some universal constant δ > 0 there always
exists a v(t) ∈ VO such that δ · v(t) ≤ x(t). Hence:
n
X
ci vi (t) ln(δ · vi (t)) = Ω(1).
fv(t) (t) ≥
i=1

22

Note that the set VO is finite, hence we get the bound maxv∈VO fv (t) = Ω(1), which finishes the
proof.
The last preparatory lemma, which can be deduced from previous considerations is the following:
Lemma 6.12. Suppose p(t) is the vector L−1 b, computed with respect to x(t). Then A⊤ p(t) is
uniformly bounded over t ∈ [0, ∞).
Proof. By a reasoning as in the proof of Lemma 6.11 we get that for some v ∈ V and a universal
constant δ > 0 we have δ · v ≤ x(t) for every t ∈ [0, ∞). Hence we get:
!
n
n
X
X
⊤
⊤ −1
⊤ −1
A p(t) = A L b = A L
vi A⊤ L−1 ai .
ai vi
≤
i=1

By Lemma 5.2, we have that
n
X
i=1

A⊤ L−1 ai

vi A⊤ L−1 ai ≤

The lemma follows.

n
X
i=1

vi

≤

αci
xi ,

i=1

for some universal constant α. Hence, we get:
n

n

i=1

i=1

X
αci X −1 αci
(δ xi )
ci = O(1).
≤
= αδ−1
xi
xi

We now present the main technical lemma required for proving existence of a limit. Intuitively
it proves that the subvector of A⊤ p corresponding to J behaves continuously assuming x is close
to the optimal set.
Lemma 6.13. Suppose ε > 0 is small enough and M, d > 0. Pick any x > 0 such that kxN k < ε
and kxJ − yJ k < ε for some optimal solution y, such that yj > d for all j ∈ J. Assume kxk ≤ M
and A⊤ p ≤ M . There exists a constant N depending only on A, b, c, d, n such that kq − xk < N ·ε.
Proof. We will denote the columns of A corresponding to J by AJ , similarly for AN . We have
(AW A⊤ )p = b, which can be rewritten as:
⊤
AJ WJ A⊤
J p + AN WN AN p = b.
−1 X
where WJ = Diag (wJ ) and WN = Diag (wN ). Since A⊤
N =
N p = O(1) and kAN WN k = AN C
⊤
O(ε) we have AN WN AN p = O(ε). Moreover we can write xJ = yJ + τJ , where kτJ k < ε.


−1
g
This gives a decomposition of WJ into WJ = Diag CJ−1 yJ and W
J = Diag CJ τJ such that
⊤
g
AJ W
J A p = O(ε). In consequence:
J

⊤
g ⊤
AJ WJ A⊤
J p − b ≤ AJ WJ AJ p − b + WJ AJ p = O(ε).

Pick p̄: any optimal solution to the dual of the linear program (1). From complementary slackness
we have Y (A⊤ p̄ − c) = 0. In particular A⊤
J p̄ = cJ . We obtain:

−1
AJ WJ A⊤
J p̄ = AJ WJ cJ = AJ Diag CJ yJ cJ = AJ yJ = b.
Hence:

⊤
⊤
⊤
AJ WJ A⊤
J (p − p̄) = AJ WJ AJ p − AJ WJ AJ p̄ = AJ WJ AJ p − b = O(ε).

We want to show that A⊤
J (p − p̄) = O(ε).
⊤
Denote K = AJ WJ AJ , it is a symmetric PSD matrix. Moreover the kernels of A⊤
J and K coincide,
since yJ > 0.
⊤
⊤
Pick a vector v ∈ Rm so that A⊤
J (p − p̄) = AJ v and v is orthogonal to the kernel of AJ . (In other
words v is the orthogonal projection of p − p̄ onto the orthogonal complement of the kernel of A⊤
J .)
23

Using the fact that

yi
ci

d
ci

≥ d′ for i ∈ J and some absolute positive constant d′ , we get:
X
X yj
′
⊤
⊤
aj a⊤
d′ AJ A⊤
=
d
a
a

j
j = AJ WJ AJ = K
J
j
cj

>

j∈J

j∈J

+ +
where  denotes the PSD ordering. This implies in particular that d′ λ+ ≤ λ+
K , where λ , λK are
⊤
the smallest positive eigenvalues of AJ AJ and K respectively. Since v is orthogonal to the kernel
of K and K is PSD we have:
′ +
kKvk2 ≥ λ+
K kvk2 ≥ d λ kvk2

Observe that λ+ depends solely on A, b, c, hence:
1
1
kvk ≤ kvk2 ≤ ′ + kKvk2 = ′ + kK(p − p̄)k2 = O(ε).
dλ
dλ
⊤
⊤
In consequence AJ (p − p̄) = AJ v = O(ε) and hence WJ A⊤
J (p − p̄) = O(ε), because kyk =
O(1) (the set of optimal solution is bounded). Finally:
⊤
⊤
⊤
⊤
⊤
WJ A⊤
J p − WJ AJ p̄ ≤ WJ AJ p − WJ AJ p + WJ AJ p − WJ AJ p̄ = O(ε) + O(ε) = O(ε).
⊤
But WJ A⊤
J p = qJ and WJ AJ p̄ = yJ , which means:

kqJ − xJ k ≤ kqJ − yJ k + kyJ − xJ k = O(ε) + O(ε) = O(ε).
To conclude, note that:
kqN − xN k ≤ kxN k + kqN k = kxN k + WN A⊤
N p = O(ε)
and hence trivially:
kq − xk = max(kqN − xN k , kqJ − xJ k) = O(ε).
Let us now conclude Theorem 6.7 from Lemma 6.13.
Proof of Theorem 6.7: By Lemma 6.8, it remains to show that kq(t) − x(t)k = O(e−δt ) for
some δ > 0. By Corollary 5.9 and Lemma 6.12 we obtain a uniform bound M on both kx(t)k and
A⊤ p(t) over t ∈ [0, ∞). From Lemma 6.11 we obtain some d > 0 such that xj (t) > 2d for j ∈ J
and t ∈ [0, ∞).
Let us now pick a large t. We take x = x(t) and y = x⋆ (t) from Lemma 6.9. We know that
xj > 2d for j ∈ J, y is exponentially close to x, hence yj > d for j ∈ J. Moreover, we can pick
ε = O(e−δt ) such that kxN k < ε (by Lemma 6.10) and ky − xk < ε. Hence the assumptions of
Lemma 6.13 are satisfied and we may conclude that kq(t) − x(t)k < N ε = O(e−δt ). The theorem
follows.
7. Discretization
In this section we study the efficiency of the natural discretization of the Physarum dynamics. It
is defined with respect to the step length 0 < h < 1 and a starting point x(0) > 0. In this section,
we will always assume that the starting point was chosen from the feasible region, i.e. Ax(0) = b.
The discrete dynamics is as follows:
(12)

x(k + 1) = (1 − h)x(k) + hq(k)

with the usual meaning of q(k). For the above dynamics to be well defined we need to make sure
that x(k) stays positive for all k ≥ 0. Note that:


 ⊤
ai p(k)
−1
xi (k + 1) = xi (k) 1 − h
ci
24

Pmax




− 1 < h1 . Corollary 5.3 tells us that:
 ⊤

ai p
def
= max
− 1 : Ax = b, x ≥ 0 ≤ Cs D + 1.
ci

hence we need to make sure that

a⊤
i p(k)
ci

−1 (and we do so in the following theorem).
This means that it is enough to set h < Pmax
−2 · ε/6. Suppose we initialize the Physarum algorithm (12)
Theorem 7.1. Let 0 < ε < 1/2, h ≤ Pmax
−1
with x(0), s.t. Ax(0) = b and Mx ≤ xi (0) ≤ Mx for every i ∈ [n] and some Mx ≥ 1. Assume

def
Mx
M
+ lnεh
steps x(k) is a feasible solution
additionally that c⊤ x(0) ≤ M ·opt. Then after k = O ln
ε2 h2
⊤
with: opt ≤ c x(k) ≤ (1 + ε)opt.

To prove Theorem 7.1 we analyze the following potential function. It is analogous to that used in
the paper [SV].
def

φ(k) = 4 ln V(k) −

(13)

εh
B(k).
opt

The intuitive meaning of φ is as follows: we know that V(k) is non-increasing with k, so ln V(k) is
non-increasing as well, however we may not guarantee a big drop of V(k) in every step, this is why
we have the second term; B(k) gets bigger at such steps. The crucial lemma we would like to show
asserts that φ drops significantly at every step:
2 2

Lemma 7.2 (Potential drop). For every k with V(k) > (1 + ε)opt we have ∆φ(k) ≤ − h 6ε .
Before proving the lemma we first conclude Theorem 7.1 from it.
Proof of Theorem 7.1 assuming Lemma 7.2: We will first estimate the value of φ(0). Observe
that
ln(V(0)) = ln(c⊤ x(0)) ≤ ln(M · opt).
The value B(0) can be bounded as follows:
n
X
i=1

Hence we obtain:

x⋆i ci ln xi (0) ≥ x⋆i ci ln Mx−1 = opt · ln Mx−1 .

φ(0) ≤ 4 ln(M · opt) + εh ln Mx .
Let us now fix k and provide a lower bound on the value of φ(k). We know that V(k) ≥ opt for
every k, for B we have:
n
n
X
X
x⋆i ci = opt · ln(Mx ),
x⋆i ci ln xi (k) ≤ ln(Mx )
B(k) =
i=1

i=1

and thus:

φ(k) ≥ 4 ln(opt) − εh ln(Mx ).
Putting together the above estimates on φ(0), φ(k) and Lemma 7.2, we obtain that:
k

h2 ε2
≤ φ(0) − φ(k) ≤ 4 ln M + εh(ln(Mx ) + ln(Mx )).
6

Simplifying:
k=O



ln M
ln Mx
+
2
2
ε h
εh



.

We split the Lemma 7.2 into two cases and deal with them separately. They are expressed in
the following facts. We use below E(k) to denote q(k)⊤ W −1 q(k).
25

Fact 7.3. If

E(k)
V(k)

2 2

< (1 − ε/3) then ∆φ(k) ≤ − h 6ε .
2 2

Fact 7.4. If E(k) > (1 + ε/3)opt then ∆φ(k) ≤ − h 6ε .
Before we proceed with the proof of facts, let us state three simple inequalities, which we are going
to use:
ln(1 + α) ≤ α

for every α ∈ R


1 1
2
(15)
ln(1 + α) ≥ α − α
for every α ∈ − ,
2 2


1
(16)
ln(1 − α) ≥ −2α
for every α ∈ 0,
2
All can be proved by simple calculus. Another useful fact is the following identity:
(14)

Fact 7.5. Let x ∈ Rn>0 and p = (AW A⊤ )−1 b be given. Suppose y ∈ Rn satisfies Ay = b, then
n
X

⊤
−1
yi a⊤
q
i p=q W

i=1

where q =

W A⊤ p.

Proof of Fact 7.3: We will show that ln V drops and B may increase only a little. Let us first
look at ln V(k + 1) − ln V(k). We have:
V(k + 1) =

n
X
i=1

⊤

⊤

c x(k + 1) = (1 − h)c x(k) + h

n
X
i=1

x(k)a⊤
i p(k) = (1 − h)V(k) + hE(k).

The last inequality follows from 7.5. Further:



V(k + 1)
E(k)
= 1+h
−1
V(k)
V(k)



ε
< 1+h 1−
−1
3
ε
=1−h .
3
We obtain:
V(k + 1) (14) hε
ln V(k + 1) − ln V(k) = ln
≤ − .
V(k)
3
Consider now B(k + 1) − B(k):
n
X
xi (k + 1)
x⋆i ci ln
B(k + 1) − B(k) =
xi (k)
i=1

 ⊤

n
X
ai p(k)
⋆
xi ci ln 1 + h
=
−1
ci
i=1

≥

n
X
i=1

(16)

x⋆i ci ln (1 − hPmax )
2

≥ −2h(n C + 1)

= −2hPmax · opt
≥ −opt.
26

n
X
i=1

x⋆e ce

The last inequality follows from the definition of h. Putting these two pieces together yields:
φ(k + 1) − φ(k) ≤ −

4hε
hε
h2 ε2
+ hε = −
≤−
.
3
3
6

Proof of Fact 7.4: We know that V(k) is non-increasing with k, it remains to show that B will
increase by a considerable amount. As in the proof of Fact 7.3 we obtain:

 ⊤

n
X
ai p(k)
⋆
x
ln
1
+
h
c
B(k + 1) − B(k) =
−1
.
i i
ci
i=1

By the choice of h:




qe (k) − xe (k)
a⊤ p(k)
1
h
− 1 ≤ h · Pmax ≤ .
= |h| i
xe (k)
ce
2
Hence we can apply the inequality 15, by which we get:
 ⊤
 X
2
 ⊤
n
n
(15) X
ai p(k)
⋆
2 ai p(k)
⋆
xi ci h
−1 −
−1
xi ci h
B(k + 1) − B(k) ≥
ci
ci
≥h

i=1
n
X
i=1

x⋆i a⊤
i p(k) − h · opt −

Fact7.5

i=1
n
X

x⋆i ci (hPmax )2

i=1

2

2
hE(k) − h · opt − h optPmax
1
≥ hE(k) − h · opt − hεopt
6

ε
1
≥ h · opt 1 +
− h · opt − hεopt
3
6
1
≥ − h · ε · opt.
6
This implies the following drop of φ:

=

φ(k + 1) − φ(k) ≤ −

hε
h2 ε2
(B(k + 1) − B(k)) ≤ −
.
opt
6

E(k)
Proof of Lemma 7.2: If ε is small enough and V(k) > (1 + ε)opt then obviously either V(k)
<
(1 − ε/3) or E(k) > (1 + ε/3)opt and we use Fact 7.3 or Fact 7.4 respectively. This concludes our
proof.

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 Automata, Languages, and Programming - 40th International Colloquium, ICALP 2013, Riga, Latvia, July 8-12, 2013,
Proceedings, Part II, pages 472–483, 2013.
[BL89]
D. A. Bayer and J. C. Lagarias. The nonlinear geometry of linear programming. I. Affine and projective
scaling trajectories. Transactions of the American Mathematical Society, pages 499–526, 1989.
[BMV12] Vincenzo Bonifaci, Kurt Mehlhorn, and Girish Varma. Physarum can compute shortest paths. In Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2012, Kyoto,
Japan, January 17-19, 2012, pages 233–240, 2012.
[IJNT11] K. Ito, A. Johansson, T. Nakagaki, and A. Tero. Convergence Properties for the Physarum Solver. ArXiv
e-prints, January 2011.
[JZ12]
Anders Johannson and James Zou. A slime mold solver for linear programming problems. In How the
World Computes, volume 7318 of Lecture Notes in Computer Science, pages 344–354. Springer Berlin
Heidelberg, 2012.
27

[Kar84]
[Kar90]
[Kar91]
[LR90]

[MS89]
[NT02]
[NYT00]
[Per01]
[Sch86]
[SV]
[TKN07]
[Wri97]

N. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica, 4(4):373–395,
1984.
N. Karmarkar. Riemannian geometry underlying interior point methods for linear programming. volume
114, pages 51–75. 1990.
Howard Karloff. Linear Programming. Birkhauser Boston Inc., 1991.
J.C. Lagarias and Vanderbei R.J. Ii dikins convergence result for the affine-scaling algorithm. In Mathematical Developments Arising from Linear Programming: Proceedings of a Joint Summer Research Conference
Held at Bowdoin College, June 25-July 1, 1988, volume 114, page 109. American Mathematical Soc., 1990.
N. Megiddo and M. Shub. Boundary behavior of interior point algorithms in linear programming. Math.
Oper. Res., 14(1):97–146, March 1989.
Yu. Nesterov and M. Todd. On the Riemannian geometry defined by self-concordant barriers and interiorpoint methods. Foundations of Computational Mathematics, pages 333–361, 2002.
Toshiyuki Nakagaki, Hiroyasu Yamada, and Agota Toth. Maze-solving by an amoeboid organism. Nature,
407(6803):470, September 2000.
Lawrence Perko. Differential equations and dynamical systems, volume 7. Springer Science & Business
Media, 2001.
Alexander Schrijver. Theory of Linear and Integer Programming. John Wiley & Sons, Inc., 1986.
Damian Straszak and Nisheeth K. Vishnoi. Natural algorithms for flow problems. In ACM SIAM Symposium on Discrete Algorithms (SODA), 2016.
Atsushi Tero, Ryo Kobayashi, and Toshiyuki Nakagaki. A mathematical model for adaptive transport
network in path finding by true slime mold. Journal of Theoretical Biology, 244(4):553, 2007.
S. Wright. Primal-Dual Interior-Point Methods. Society for Industrial and Applied Mathematics, 1997.

Appendix A. Dynamical Systems and Global Existence
We work in the n-dimensional Euclidean space Rn , an open subset Ω ⊆ Rn is chosen for the
domain of the dynamical system. The equations are of the form:
ẋi (t) = Fi (x1 (t), x2 (t), . . . , xn (t))

for i = 1, 2, . . . , n,

where x1 , x2 , . . . , xn are functions of single variable t regarded as unknowns and F1 , F2 , . . . , Fn :
Ω → R are given. The above system is often denoted shortly as:
ẋ = F (x)

Rn

where F : Ω →
is just F = (F1 , F2 , . . . , Fn ). Typically we impose some initial condition of the
form x(0) = s, with s ∈ Ω. A solution to such a system is a function x : I → Ω, where I is some
non-degenerate interval containing 0, x(0) = s and ẋ(t) = F (x(t)) for every t in the interior of I.
We list some typical questions which arise in the study of a dynamical system. Suppose we are
given a dynamical system S.
(1) Is there any solution to S? What is the maximal interval of existence?
(2) Is the solution unique?
(3) What are the local properties of S (around point x(0))?
(4) What are the global properties of the solution? Does it converge? Does it stay in a bounded
region?
We are interested mainly in the last question, in both qualitative and quantitative aspects (whether
it does converge and what is the convergence rate). However, to attempt question (4) we must first
answer (1) and (2); otherwise we may end up studying objects which do not even exist. Moreover,
in many cases existence (and uniqueness10) can be very hard to establish.
In our case it is not difficult to obtain a solution on some small interval (−ε, ε). However, for
applications we need existence on [0, +∞), (we call it global existence) which is harder to establish.
Intuitively it may happen that at some finite time step T the solution x will “escape” from the
domain Ω. More formally, even if the solution x exists in the interval [0, T ), the limit limt→T − x(t)
may be outside of Ω or the limit may not even exist. Escaping from Ω is not only a formal problem
10Note that even if a differential equation has a solution it may be not unique. Consider for example the following:
ẋ = 3x2/3 with x(0) = 0. It has two solutions x(t) = t3 and x(t) = 0.
28

of defining Ω, it is actually a serious issue, because in our case, the function F is not well defined
outside of Ω.
In this section we provide a sufficient condition for a dynamical system to have global solutions.
In the main body we prove that in fact the Physarum dynamical system satisfies this condition.
We first state an important general theorem which gives sufficient conditions for existence on
some interval around 0. A proof can be found in any textbook on dynamical systems, e.g. [Per01].
Theorem A.1. Let Ω ⊆ Rn be an open set, s ∈ Ω and F : Ω → Rn be a function of class C 1 .
Consider the following dynamical system: ẋ = F (x) with initial condition x(0) = s. There exists a
unique solution x : (−ε, ε) → Ω to the system, for some ε > 0.
We conclude this subsection by proving a theorem that a certain simple condition implies existence
on the whole half-line.
Theorem A.2. Let Ω ⊆ Rn be an open set, s ∈ Ω and F : Ω → Rn be a function of class C 1 .
Consider the following dynamical system: ẋ = F (x) with initial condition x(0) = s. Suppose it
satisfies the following conditions for every solution x : [0, T ) → Ω with T ∈ (0, ∞):
• the limit limt→T − x(t) exists,
• if x(T ) = limt→T − x(t) then x(T ) ∈ Ω.
Then the there exists a global solution x : [0, ∞) → Ω.
Proof. Let T0 be the maximal T such that a solution x : [0, T ) → Ω exists (we allow T0 = ∞). One
should argue why does such a T0 exists. This is a consequence of uniqueness: whenever there are
two solutions x : [0, T1 ) → Ω and y : [0, T2 ) → Ω with T1 ≤ T2 , then they agree on [0, T1 ).
From Theorem A.1 we know that T0 > 0. We want to show that T0 = ∞. For the sake
of contradiction assume that T0 is a finite number. Denote by x(T0 ) the limit limt→T − x(t). By
0
assumption x(T0 ) ∈ Ω, thus we can use Theorem A.1 to extend the solution from [0, T0 ) to [0, T0 +ε)
for some ε > 0. This gives us a contradiction with the definition of T0 .
An example of a system which does not satisfy the the assumption of Theorem A.2.
1
, with x(0) = 0. The natural choice
Example A.3. Consider a one-dimensional system ẋ = 1−x
1
for the domain is Ω = R \ {1} (only on this set the function x 7→ 1−x
makes sense). Let T = 12 ,
√
then there is a solution x : [0, T ) → Ω given by x(t) = 1 − 1 − 2t. However:

lim x(t) = 1.

t→T −

But 1 ∈
/ Ω, hence this system does not satisfy the condition from Theorem A.2. Note also that
actually there is no solution to this system on the whole half-line [0, ∞).

29

