# Convergence of the Non-Uniform Physarum Dynamics

Let c ∈Z^m_> 0, A ∈Z^n× m, and b ∈Z^n. We show under fairly general conditions that the non-uniform Physarum dynamics ẋ_e = a_e(x,t) (|q_e| - x_e) converges to the optimum solution x^* of the weighted basis pursuit problem minimize c^T x subject to A f = b and |f| < x. Here, f and x are m-vectors of real variables, q minimizes the energy ∑_e (c_e/x_e) q_e^2 subject to the constraints A q = b and supp(q) ⊆supp(x), and a_e(x,t) > 0 is the reactivity of edge e to the difference |q_e| - x_e at time t and in state x. Previously convergence was only shown for the uniform case a_e(x,t) = 1 for all e, x, and t. We also show convergence for the dynamics ẋ_e = x_e ·( g_e (|q_e|/x_e) - 1), where g_e is an increasing differentiable function with g_e(1) = 1. Previously convergence was only shown for the special case of the shortest path problem on a graph consisting of two nodes connected by parallel edges.

There are no comments yet.

## Authors

• 7 publications
• 6 publications
• 16 publications
10/24/2018

### Minsum k-Sink Problem on Path Networks

We consider the problem of locating a set of k sinks on a path network w...
06/18/2019

### Convergence of the Non-Uniform Directed Physarum Model

The directed Physarum dynamics is known to solve positive linear program...
08/27/2018

### On the convergence of optimistic policy iteration for stochastic shortest path problem

In this paper, we prove some convergence results of a special case of op...
09/03/2020

### Physarum Multi-Commodity Flow Dynamics

In wet-lab experiments, the slime mold Physarum polycephalum has demonst...
07/16/2020

### A dichotomy theorem for nonuniform CSPs simplified

In a non-uniform Constraint Satisfaction problem CSP(G), where G is a se...
03/26/2020

### Tuned Hybrid Non-Uniform Subdivision Surfaces with Optimal Convergence Rates

This paper presents an enhanced version of our previous work, hybrid non...
04/14/2020

### Adaptive radial basis function generated finite-difference (RBF-FD) on non-uniform nodes using p—refinement

Radial basis functions-generated finite difference methods (RBF-FDs) hav...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In [BBK17] and [FCP18] it was shown that the Physarum dynamics (2) can solve the weighted basis pursuit problem

 minimize cTx subject to Af=b and |f|≤x, (1)

where , , , , and and are an -vectors of real variables. The absolute-value operator is applied componentwise. is assumed to have full row-rank. For simplicity, we also assume that any two basic feasible solutions of have distinct cost111A basic feasible solution of has the form , where is a subset of of size , , the submatrix of is invertible, , and . The cost of such a solution is .; in particular, the optimal solution to (1) is unique. [BBK17, FCP18] actually shows convergence under the more general condition and for any in the kernel of , but we do not need this generality here. We index the rows of by and the columns of by and, for historical reasons (see Section 2.1), refer to the rows as nodes and the columns as edges.

The Physarum dynamics evolves a vector according to the dynamics

 ˙x=|q|−x, (2)

where is the minimum energy feasible solution according to the resistances :

 q=argminf{∑e(ce/xe)f2e:Af=b% and fe=0 whenever xe=0}.

The Physarum dynamics was introduced by biologists [TKN07] as a model of the behavior of the slime mold Physarum Polycephalum. We discuss the biological background in Section 2.

###### Theorem 1 ([Bbk+17, Fcp18]).

Assume a strictly positive start vector . Then the solution to the dynamics (2) is defined on , for all , and and converge to the optimal solution of (1).

In this paper, we consider the more general dynamics (3) and (5). In the dynamics (3) the edges react with different speed to differences between minimum energy solution and capacity, i.e.,

 ˙xe=ae(x,t)(|qe|−xe), (3)

where is the reactivity of edge at time , i.e., the edges no longer react uniformly to differences between and , but the reactivity depends on the edge, the current state, and the time. We refer to (3) as the nonuniform Physarum dynamics. The special case that is a positive constant for each edge was introduced in [NIU07] to model the behavior of Physarum Polycephalum in non-uniform environments; see Section 2.2. We show

###### Theorem 2.

Assume , for all , , and and some constant , and is Lipschitz-continuous. Then:

1. The dynamics (3) has a unique solution for .

2. The function

 L(x,t)=∑e(ce/xe)q2e+∑ecexe (4)

is a Lyapunov function for the dynamics, i.e., for all . Moreover, if and only if for all : either or or .

3. If, in addition, for some positive and all , , and , is a fixed point of (3) if and only if for a basic feasible solution of .

4. Under the same additional assumption as in (c), and converge to a fixed point of (3) as goes to infinity.

5. If, in addition, does not depend on and for all and , and converge to as goes to infinity. In particular, this holds true if is a positive constant for all .

The proof of part (a) is standard and part (c) was shown in [BBK17]. The Lyapunov function in part (b) was introduced in [FD18]. In  [FCP18] it was shown to be Lyapunov function for the uniform case ( for all , , and ). Part (d) follows easily from parts (b) and (c). Finally, the proof of part (e) is inspired by by [BBK17].

The function is a Lyapunov function for the Physarum dynamics under very general conditions. Essentially, the only requirement is that has the same sign as . For the existence of a solution with domain , we also need that is bounded. For the convergence to a fixed point, we need in addition that is bounded away from zero.

In the dynamics (5) each edge has its own transfer function that determines how its reacts to the ratio of flow and capacity being larger or smaller than one, i.e.,

 ˙xe=xe⋅(ge(|qe|xe)−1)for all e∈[m], (5)

where the response function is assumed to be an increasing differentiable function satisfying . Bonifaci introduces this model in [Bon16] in order to deal with the larger class of response functions proposed in the biological literature. For the shortest path problem in a network of parallel links222The shortest path problem is a min-cost flow problem where we want to send one unit of flow between two distinguished nodes. For the case of parallel links, the graph has exactly two nodes and all edges run between these nodes. We have a graph with two nodes and several parallel edges connecting them,, [Bon16] shows convergence to the shortest path. Bonifaci assumes the same response function for every edge, but his proof actually works for response functions depending on the edge. We show

###### Theorem 3.

Assume is increasing and differentiable and satisfies for all .

1. The dynamics (5) has a unique solution for .

2. is a fixed point of (5) if for a basic feasible solution of (1).

3. The function

 L(x,t)=∑e(ce/xe)q2e+∑ecexe (6)

is a Lyapunov function for the dynamics, i.e., for all . Moreover, if and only if is a fixed point.

4. and converge to a fixed point of (5).

5. If, in addition, for some and all and , and converge to as goes to infinity.

The proof of part (a) is standard and part (b) was shown in [BBK17]. The Lyapunov function in part (c) was introduced in [FD18]. Part (d) follows easily from part (c). Finally, the proof of part (e) is inspired by by [BBK17]. The theorem above also holds if the function may depend on time and state.

Nature does not compute exactly, i.e., one should not expect that is exactly equal to or to in a biological system. Rather, there will be noise. Our results show that the dynamics (3) and (5) are fairly robust against noise, i.e., variations in and .

In Section 2, we review the biological background and related work, in Section 3, we prove Theorem 2, and in Section 4, we prove Theorem 3. In Section 5, we state some open problems.

## 2 Background

### 2.1 The Shortest Path Experiment

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 following 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 with two distinguished nodes modeling the food sources. Each edge has a positive length and a positive capacity ; is fixed, but is a function of time. The resistance of is . In the electrical network defined by these resistances, a current of value is forced 1 from one of the distinguished nodes to the other. For an (arbitrarily oriented) edge , let be the resulting current over . Then, the capacity of evolves according to the differential equation

 ˙xe(t)=|qe(t)|−xe(t), (7)

where is the derivative of with respect to time. In equilibrium ( for all ), 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. It is well-known that the electrical flow is the feasible flow minimizing energy dissipation (Thomson’s principle).

### 2.2 Minimum Risk Paths

In [NIU07], Nakagaki et. al. study the following scenario, see Figure 2. They cover a rectangular plate with Physarum and feed it at opposite corners of the plate. Two thirds of the plate is put under a bright light, one third is kept in the dark. Under uniform lighting conditions, Physarum would retract to a straight-line path connecting the food sources [NYT00]. However, Physarum does not like light and therefore forms a path with one kink connecting the food sources. The path is such that the part under light is shorter than in a straight-line connection. In the theory section of [NIU07], the dynamics

 ˙xe(t)=|qe|−aexe(t) (8)

is proposed. The constant is the decay rate of edge if there is no flow on it. To model the experiment, for edges in the dark part of the plate, and for the edges in the lighted area, where is a constant. [NIU07] report that in computer simulations, the dynamics (8) converges to the shortest source-sink path with respect to the modified cost function .

### 2.3 A Reformulation: Nonuniform Physarum

Let . The electrical flow is determined by the resistances . Therefore, we write instead of for clarity. Next observe that . Thus if we take as the vector of edge capacities and as the vector of costs, we get the same electrical flow. We can express (8) as a dynamics for as

 ˙ye=ae˙xe=ae(|qe(re)|−aexe)=ae(|qe(re)|−ye).

So we may instead consider the dynamics

 ˙ye=ae(|qe|−ye)

under the modified cost function . This is our dynamics (3), where we generalized further by allowing to depend on and . In this model, the quantity indicates the responsiveness (reactivity) of an edge to differences between flow and diameter.

### 2.4 Beyond Shortest Paths

The biological experiments concern shortest paths. The papers [BMV12, Bon13] showed Theorem 1 for the shortest path problem and the transportation problem; here is the node-arc incidence matrix of a directed graph, is the supply-demand vector of a transportation problem, i.e., , and are the edge costs. Convergence for the discretization of (2) was shown in [BBD13].

The theoretical literature soon asked whether the dynamics (2) can also solve more general problems. The basis pursuit problem was first studied in [SV16a] and convergence of the discretization was shown. Theorem 1 was shown in [BBK17]. The function (6) was introduced in [FD18] and shown to be a Lyapunov function for (2) in [FCP18].

The paper [Bon16] introduces and studies the dynamics (5)

The directed version of the Physarum dynamics evolves according do the equation

 ˙xe=qe−xe. (9)

No biological significance is claimed for this dynamics. It can solve linear programs with positive cost vectors

[IJNT11, SV16b]. In [BBD13], convergence was claimed for the non-uniform dynamics . The proof is incorrect.333 In the unpublished proof of Lemma 10, the authors argue: We first follow the development in [IJNT11], taking the reactivities into account. From
we have
and hence
Let
Since , is a convex combination of flows and hence a flow. This argument is fallacious. For each edge is a convex combination. But these combinations are not uniform over edges. Therefore is NOT a convex combination of flows.

## 3 The Proof of Theorem 2

### 3.1 Preliminaries

For a capacity vector and a vector with , we use to denote the energy of . The energy of is infinite, if . We use to denote the cost of . Note that .

We use to denote the diagonal matrix with entries ; here we use the convention that attention is restricted to the edges with . In part (a) of Theorem 2, it is shown that for all if . However, in the limit some edges may have capacity zero. Energy-minimizing solutions are induced by node potentials according to the following equations:

 b =Aq (10) q =R−1ATp (11) AR−1ATp =b (12)

We give a short justification. minimizes the quadratic function subject to the constraints . The KKT conditions state that at the optimum, the gradient of the objective is a linear combination of the gradients of the constraints. Thus

 2(ce/xe)qe=∑ipiAi,e

for some vector . Absorbing the factor into yields equation (11). Substitution of (11) into (10) gives (12).

We next collect some well-known properties of the minimum energy solution; the proof of part (b) can, for example, be found in [BBK17]. Let be the maximum absolute value of a square submatrix of .

1. The minimum energy solution is defined by (12) and (11). It is unique.

2. for every .

3. , where is defined by (12). This holds since .

With the help of (11), the dynamics can we rewritten as

 ˙xe=ae(x,t)(∣∣∣xeceATep∣∣∣−xe)=ae(x,t)xe(|ATep|ce−1). (13)

### 3.2 Existence

The right-hand side of (3) is locally Lipschitz-continuous in and . The function is locally Lipschitz by assumption, is an infinitely often differentiable rational function in the and hence locally Lipschitz. Furthermore, locally Lipschitz-continuous functions are closed under additions and multiplications. Thus is defined and unique for for some .

Since for all , and , we have and hence . Also since for all and , we have and hence for all . In particular, the solution is bounded. Thus,

by well-known results of maximal solutions of ordinary differential equations

[Har01, Corollary 3.2].

The condition is crucial for existence. Let , and . The matrix is , i.e., there are no constraints. Then the minimum energy solution is the null-vector of dimension one and (3) becomes ; the domain of definition is .

### 3.3 Fixed Points

A point is a fixed point if . In  [BBK17] is was shown that the fixed points of (2) are the vectors , where is a basic feasible solution of (1). This uses the assumption that any two basic feasible solutions have distinct cost. The proof carries over to (3) under the additional assumption that for all , and and some positive . Under this additional assumption is equivalent to for (2) and (3). This section is reprinted from [BBK17] with minor adaptions. A vector is sign-compatible with a vector (of the same dimension) if implies . In particular, . We use the following corollary of the finite basis theorem for polyhedra.

###### Lemma 1.

Let be a feasible solution of (1). Then is the sum of a convex combination of at most basic feasible solutions plus a vector in the kernel of . Moreover, all elements in this representation are sign-compatible with .

###### Proof.

We may assume . Otherwise, we flip the sign of the appropriate columns of . Thus, the system is feasible and is the sum of a convex combination of at most basic feasible solutions plus a vector in the kernel of by the finite basis theorem [Sch03, Corollary 7.1b]. By definition, the elements in this representation are non-negative vectors and hence sign-compatible with . ∎

###### Lemma 2.

Assume for some positive and all , , and , and that no two feasible solutions of have the same cost. If is a basic feasible solution of (1), then is an fixed point. Conversely, if is an fixed point, then for some basic feasible solution .

###### Proof.

Let be a basic feasible solution, let , and let be the minimum energy feasible solution with respect to the resistances . We have and by definition of . Since is a basic feasible solution there is a subset of size of the columns of such that is non-singular and . Since , we have for some vector . Thus, and hence . Therefore and is an fixed point.

Conversely, if is an fixed point, for every . By changing the signs of some columns of , we may assume . Then . Since by (11), we have , whenever . By Lemma 1, is a convex combination of basic feasible solutions plus a vector in the kernel of that are sign-compatible with . The vector in the kernel is zero since is a minimum energy solution444Assume with , , , and . Then , , and , a contradiction.. For any basic feasible solution contributing to , we have . Summing over the , we obtain , i.e., all basic feasible solutions used to represent have the same cost. Since we assume the costs of distinct basic feasible solutions to be distinct, is a basic feasible solution. ∎

###### Corollary 1.

Assume for some positive and all , , and and that no two feasible solutions of have the same cost. The set of fixed points is a discrete set.

### 3.4 The Lyapunov Function

###### Lemma 3.

is a Lyapunov function. More precisely, always with equality only if for all either or or .

###### Proof.

Taking the derivative of (12) with respect to time yields

We next compute the derivative of both summands of with respect to time separately. For the first summand we obtain

where the first equality uses (12), the second equality follows from the product rule of differentiation, the third equality follows from 14, the fourth equality is a simple algebraic manipulation, the fifth equality follows from (13), and the last equality is a simple algebraic manipulation.

For the second summand, we obtain

 ddtcTx=∑ece˙xe=∑eaece(xece|ATep|−xe)=∑eaecexe(|ATep|ce−1).

Combining the two terms and writing instead of , we obtain

 ddtL(x,t)=−∑eaecexe(λ3e−λ2e−λe+1).

Since

 λ3e−λ2e−λe+1=(λ2e−1)(λe−1)=(λe+1)(λe−1)2

and , we have always. We have equality only if for all , i.e., for all either or or . ∎

###### Corollary 2.

Assume further for some positive and all , and . Then if and only if is a fixed point.

###### Proof.

We have if and only if for all either or . The latter condition is equivalent to . Thus . ∎

### 3.5 Convergence

From now on, we make the addition assumption that for some positive and all , , and . It then follows from the general theory of dynamical systems that converges to a fixed point.

###### Corollary 3 (Generalization of Corollary 3.3. in [Bon13].).

Assume further for all , and . As , and approach a fixed point . Moreover, and converge to .

###### Proof.

The proof in [Bon13] carries over. We include it for completeness. The existence of a Lyapunov function implies by [LaS76, Corollary 2.6.5] that approaches the set , which by Corollary 2is the same as the set . Since this set consists of isolated points (Lemma 2), must approach one of those points, say the point . When , one has . ∎

The assumption is crucial as the following example shows. Let , consider the task of minimizing subject to the constraint , and let and . Then . Integrating from to and observing that for all , we obtain

 x(t)=x(0)+∫t0(e−s/2)(1−x(s))ds≤12+14∫t0e−sds≤12+14(1−e−t)≤34

and hence the dynamics does not converge to the optimal solution , which, in this case, is the only fixed point.

It remains to exclude that converges to a non-optimal fixed point. We can do so under an additional assumption on .

###### Theorem 4.

Assume further that does not depend on , i.e., , for some positive for all and , and for all and . As , converges to the optimal solution .

###### Proof.

Assume that converges to a non-optimal fixed point . Let be the optimal solution, let be such that for all and (by Section 3.2 the solution is bounded), and let . Let . Note that for all sufficiently large , we have . Further, by definition and thus

 ˙W=∑ex∗ece|qe|−xexe+∑ex∗ece(−˙aea2e)ln(xe/B)≥∑ex∗e|ATep|−cost(x∗)≥∑ex∗eATep−cost(x∗)≥δ,

where the first inequality follows from and and the second inequality follows from . Hence , a contradiction to the fact that is bounded. ∎

## 4 Bonifaci’s Refined Model

Bonifaci [Bon16] investigates the dynamics

 ˙xe=xe(ge(|qe|xe)−1)for all e∈[m],

where the response function is assumed to be an increasing differentiable function satisfying . For the shortest path problem in a network of parallel links, Bonifaci shows convergence to the shortest path. Bonifaci assumes the same response function for every edge, but his proof actually works for response functions depending on the edge. Concrete response functions of this type had been considered earlier in the literature:

• Non-saturating response: for some .

• Saturating response: for some .

###### Lemma 4.

is a Lyapunov function for the dynamics (5). Moreover if and only if is a fixed point of (5).

###### Proof.

We proceed as in the proof of Lemma 3.

 ddtL(x,t) =ddt(pTb+cTx) =−∑e((ATep)2˙xe/ce−ce˙xe) =−∑e(cexe(|ATep|2c2e−1)⋅(ge(|ATep|ce)−1)) ≤0,

since the signs of and agree. Note that and is increasing.

Moreover, the derivative is zero if and only if for all , either or . Finally note that for with , is equivalent to which in turn is equivalent to . ∎

We remark that the proof above would even work for transfer-functions . It is only important that and that the function is increasing in .

It now follows from the general theory of dynamical systems that converges to a fixed point.

###### Corollary 4 (Generalization of Corollary 3.3. in [Bon13].).

As , and approach a fixed point . Moreover, and converge to .

###### Proof.

Same proof as Corollary 3. ∎

We finally show convergence to the optimum solution of (1) under the additional assumption that for some and all and .

###### Theorem 5.

Assume further that for some and all and . As , converges to the optimal solution .

###### Proof.

Assume that converges to a non-optimal fixed point . Let be the optimal solution and let . Let . Note that for all sufficiently large , we have . Further, by definition and thus

 ˙W =∑ex∗ecexe(ge(|qe|xe)−1)xe≥∑ex∗eceα(|ATep|ce−1) ≥α⋅(∑ex∗e|ATep|−cost(x∗))≥α⋅(∑ex∗eATep−cost(x∗))≥α⋅δ,

where the first inequality follows from for all and the last inequality follows from . Hence , a contradiction to the fact that is bounded. ∎

Convex increasing functions satisfy with .

## 5 Open Problems

We showed convergence to the optimal solution under the assumption that is bounded and bounded away from zero, that does not depend on , and always. We argued that the first two assumptions are necessary. How about the other assumptions?

For the uniform dynamics, convergence of a suitable Euler discretization was shown in [BBD