Convergence of the Non-Uniform Directed Physarum Model

TL;DR

This paper proves that a generalized non-uniform directed Physarum model, which allows different reaction speeds for each component, converges to solutions of positive linear programs, extending the known uniform case.

Contribution

It introduces and analyzes the non-uniform directed Physarum dynamics, demonstrating its convergence to positive linear program solutions, thus generalizing previous uniform models.

Findings

The non-uniform dynamics solves positive linear programs.

Convergence is established for the generalized model.

The model allows component-wise reaction speeds.

Abstract

The directed Physarum dynamics is known to solve positive linear programs: minimize $c^T x$ subject to $Ax = b$ and $x \ge 0$ for a positive cost vector $c$. The directed Physarum dynamics evolves a positive vector $x$ according to the dynamics $\dot{x} = q(x) - x$. Here $q(x)$ is the solution to $Af = b$ that minimizes the "energy" $\sum_i c_i f_i^2/x_i$. In this paper, we study the non-uniform directed dynamics $\dot{x} = D(q(x) - x)$, where $D$ is a positive diagonal matrix. The non-uniform dynamics is more complex than the uniform dynamics (with $D$ being the identity matrix), as it allows each component of $x$ to react with different speed to the differences between $q(x)$ and $x$. Our contribution is to show that the non-uniform directed dynamics solves positive linear programs.

Equations

$$\text{minimize }{c}^{T}x\text{ subject to }Ax=b,x\ge 0.$$

$$I=\{i|i\in supp(x^{\ast })\text{ for some optimal solution x∗ to }(\text{ref}LP)\}$$

$$F=\{x|\text{x is a feasible solution to (refLP)}\}$$

$$\dot{x}=q(x)-x,$$

$$q(x)={\mathrm{argmin}}_f\left\{\sum _i\frac{c_i}{x_i}{f}_i^2|Af=b\right\}$$

$$\dot{x}=D(q(x)-x),$$

$$V(x)=2c^T{D}^{-1}x-\sum _i\frac{c_i{x}_i^{\ast }}{d_i}\ln x_i$$

$$G^{\ast }=\{x|x\ge 0\text{ and }{x}_i>0\text{ for all }i\in I\}.$$

$$E=\{x\in G^{\ast }|\dot{V}(x)=0\}.$$

$${\dot{x}}_e(t)=|q_e(t)|-x_e(t),$$

$${\dot{x}}_e(t)=d_e(|q_e(t)|-x_e(t)),$$

$$\text{minimize }{c}^T|x|\text{ subject to }Ax=b,$$

$${\dot{x}}_e(t)=q_e(t)-x_e(t).$$

$$\text{minimize }{c}^{T}x\text{ subject to }Ax=b,x\ge 0,$$

$${\dot{x}}_e(t)=d_e(q_e(t)-x_e(t)).$$

$$G^{\ast }$$

$$G_B^{\ast }$$

$$q(x)={\mathrm{argmin}}_f\left\{\sum _{i\in supp(x)}\frac{c_i}{x_i}{f}_i^2|\text{Af=b and supp(f)⊆supp(x)}\right\}.$$

$${\dot{x}}_i=f_i(x):=\{\begin{array}{ll}{d}_i(q_i-x_i)& \text{ if i∈supp(x) }\\ 0& \text{ if i∈supp(x)}\end{array}$$

$$\Vert A_B^T{L}_B^{-1}{A}_i{\Vert }_{\infty }\le \frac{c_i}{x_i}\cdot M$$

$$x_i(t)\le \beta +(x_i(0)-\beta )e^{-d_{i}t}\le \max (x_i(0),\beta ).$$

$$|{\dot{x}}_i(t)|=d_i|q_i(t)-x_i(t)|\le d_i\left(|q_i(t)|+|x_i(t)|\right)\le d_i\left(\beta +\max (x_i(0),\beta ))\right).$$

$$W(t)=\sum _{j\in I}\frac{c_j{y}_j}{d_j}\ln x_j(t).$$

$$\sum _{j\in I}{y}_j{r}_j(t)q_j(t)=y^{T}R(t)q(t)=y^T{A}^{T}p(t)=b^{T}p(t)=q(t{)}^{T}R(t)q(t)\ge 0$$

$$\dot{W}(t)=\sum _{j\in I}\frac{c_j{y}_j}{d_j}\cdot d_j\left[\frac{q_j(t)}{x_j(t)}-1\right]\ge -c^{T}y=-O(1).$$

$$ϵ=\underset{i\in supp(y)}{\min }\underset{t\in [0,T]}{\min }\frac{x_i(t)}{y_i}$$

$$x(t)\ge ϵ\cdot y,\forall t\in [0,T]$$

$$\Vert A_B^T{p}_B(t){\Vert }_{\infty }$$

$${\dot{x}}_i(t)=d_i\left(q_i(t)-x_i(t)\right)=d_i{x}_i(t)\left(\frac{p_B^T[A_B{]}_i}{c_i}-1\right)\ge -x_i\cdot \Vert d{\Vert }_{\infty }\left(\frac{N}{c_i}+1\right)\ge -x_i\cdot N^{\prime }.$$

$$x_i(t)\ge e^{-t\cdot N^{\prime }}\cdot x_i(0)>0.$$

$$V(x)=2c^T{D}^{-1}x-\sum _{i\in I}\frac{c_i{x}_i^{\ast }}{d_i}\ln x_i.$$

Full text

In this paper, we study the non-uniform directed dynamics $\dot{x}=D(q(x)-x)$, where $D$ is a positive diagonal matrix. The non-uniform dynamics is more complex than the uniform dynamics (with $D$ being the identity matrix), as it allows each component of $x$ to react with different speed to the differences between $q(x)$ and $x$. Our contribution is to show that the non-uniform directed dynamics solves positive linear programs.

1 Introduction

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 2. 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, the Physarum dynamics, for the behavior of the slime in the form of a system of coupled differential equations. In [BMV12, Bon13] it was shown that the Physarum dynamics solves the shortest path problem and the transportation problem. It was soon asked whether the dynamics can also solve more complex problems.

A variant of the Physarum dynamics, the directed Physarum dynamics, is known to solve positive linear programs in standard form [JZ12, SV16c]. A positive linear program asks to minimize a linear function $c^{T}x$ with a positive cost vector $c\in\mathbb{R}_{>0}^{m}$ subject to the constraints $Ax=b$ and $x\geq 0$. Here $A\in\mathbb{R}^{n\times m}$ and $b\in\mathbb{R}^{n}$. Formally,

[TABLE]

We assume throughout that the system is feasible and use

[TABLE]

to denote the union of the supports of optimal solutions and

[TABLE]

to denote the set of feasible solution. In order to avoid trivialities, we assume $I\not=\emptyset$. The directed Physarum dynamics is defined as the dynamical system operating on $x\in\mathbb{R}_{>0}^{m}$ according to

[TABLE]

where $\dot{x}$ is the derivative of $x$ with respect to time, $x>0$ and

[TABLE]

is the minimum energy solution of $Af=b$ according to the weights (“resistances”) $c_{i}/x_{i}$. The system is initialized to a point $x^{0}\in G\mathop{:=}\{\;x\;|\;x>0\;\}\subseteq\mathbb{R}^{m}$, i.e., $x(0)=x^{0}\in G$.

Theorem 1** ([JZ12, SV16c]).**

The dynamics (2) has a solution $x(t)\in G$ with $t\in[0,\infty)$. The solution satisfies

$\inf_{t}x_{i}(t)>0$* for all $i\in I$,* 2. 2.

$\mathrm{dist}(x(t),F)\rightarrow 0$* as $t\rightarrow\infty$, and* 3. 3.

$\lim_{t\rightarrow\infty}c^{T}x(t)=c^{T}x^{*}$.

In this paper, we study the non-uniform Physarum dynamics

[TABLE]

where for a fixed positive vector $d\in\mathbb{R}_{>0}^{m}$ we define by $D=\mathrm{diag}(d)\in\mathbb{R}^{m\times m}$ a diagonal matrix with entries $D_{ij}=d_{i}$ if $i=j$ and $D_{ij}=0$ otherwise. In this model (3), different components of $x$ react with different speed to differences between $q(x)$ and $x$. Again the system is initialized with a point $x^{0}\in G$, i.e., $x(0)=x^{0}$. The original and the non-uniform Physarum dynamics are both inspired by the study of the slime mold Physarum Polycephalum; see Section 2. We generalize the theorem above to the non-uniform dynamics.

Theorem 2**.**

The dynamics (3) has a solution $x(t)\in G$ with $t\in[0,\infty)$. The solution satisfies

$\inf_{t}x_{i}(t)>0$* for all $i\in I$,* 2. 2.

$\mathrm{dist}(x(t),F)\rightarrow 0$* as $t\rightarrow\infty$, and* 3. 3.

$\lim_{t\rightarrow\infty}c^{T}x(t)=c^{T}x^{*}$.

This theorem was claimed in [BBD*+*13], albeit with an incorrect proof; see Section 2. We note that the non-uniform dynamics (3) is more complex that the uniform dynamics (2). For example, if $x(0)$ is a feasible solution to (1), then the solution $x(t)$ is feasible222 Consider $y(t)=Ax(t)-b$. Then $\dot{y}=A\dot{x}=A(q-x)=b-Ax=-y$. Thus $y_{i}(t)=y_{i}(0)e^{-t}$ for all $t$ and hence $y_{i}(0)=0$ implies $y_{i}(t)=0$ for all $t$. In particular, if $x(0)$ is feasible, $x(t)$ is feasible for all $t$. More generally, for arbitrary initial point, the residual $Ax(t)-b$ converges to zero exponentially fast. for all $t$. In contrast, the non-uniform dynamics may move out of feasibility as Figure 1 illustrates.

A fundamental tool for proving convergence of a dynamical system is a Lyapunov function. It is a function defined on the states of the system whose derivative with respect to time is non-positive. We exhibit a Lyapunov function for the dynamics (3). Let $x^{*}$ be any optimal solution to (1) with $I={\mathit{supp}}(x^{*})$ (note that a convex combination of optimal solutions is also an optimal solution. Therefore there is an optimal solution with support equal to $I$) and consider

[TABLE]

Theorem 3**.**

$\dot{V}(x)\leq 0$* for all $x\in G$ and there is a constant $C$ (depending on $x(0)$) such that $V(x(t))\geq C$ for all $t\in(0,\infty)$.*

It follows that $\lim_{t\rightarrow\infty}V(x(t))$ exists and $\lim_{t\rightarrow\infty}\dot{V}(x(t))=0$. By the first item of Theorem 2 any limit point333A point $p$ is a (positive) limit point if there is a sequence $t_{n}\rightarrow\infty$ such that $x(t_{n})\rightarrow p$ as $n\rightarrow\infty$. of the dynamics (3) will lie in

[TABLE]

It then follows that the dynamics converges to the set

[TABLE]

We will see below that $\dot{V}(x)=0$ implies $c^{T}x=c^{T}x^{*}$ and $q(x)=x$. In particular, $x\in E$ implies $x\geq 0$ and $x_{i}>0$ for $i\in{\mathit{supp}}(x^{*})$, $q(x)=x$, i.e., $x$ is a fixed point of the dynamics and a feasible solution to the LP (since $Ax=Aq(x)=b)$, and $c^{T}x=c^{T}x^{*}$, i.e., $x$ is an optimal solution to the LP.

This paper is structured as follows. In Section 2 we review the biological background and discuss related work. In Section 3 we prove existence of a solution with domain $[0,\infty)$. In Section 4 we analyze the Lyapunov function. In Section 5 we complete the proof of Theorem 2. In Section 6, we study in more detail how the solution approaches the optimum in the example of Figure 1. In Section 7, we conduct a numerical experiment that explores the convergence-time dependence of the non-uniform dynamics (3) on the reactivity matrix $D$.

2 Biological Background and Related Work

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 2. 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 $G=(N,E)$ with two distinguished nodes modeling the food sources. Each edge $e\in E$ has a positive length $c_{e}$ and a positive capacity $x_{e}(t)$; $c_{e}$ is fixed, but $x_{e}(t)$ is a function of time. The resistance $r_{e}(t)$ of $e$ is $r_{e}(t)=c_{e}/x_{e}(t)$. In the electrical network defined by these resistances, a current of value 1 is forced from one of the distinguished nodes to the other. For an (arbitrarily oriented) edge $e=(u,v)$, let $q_{e}(t)$ be the resulting current over $e$. Then, the capacity of $e$ evolves according to the differential equation

[TABLE]

where $\dot{x}_{e}$ is the derivative of $x_{e}$ with respect to time.

Nakagaki et al. [NIU*+*07] pointed out that different edges may react with different speed to the differences between flow and capacity. For example, Physarum prefers darkness over bright light and hence edges in a bright environment react differently than edges in darkness. This let to the non-uniform dynamics

[TABLE]

where $d_{e}$ is an indicator for the reactivity of an edge.

The biological experiments concern shortest paths. The papers [BMV12, Bon13] showed that the Physarum dynamics can solve the shortest path problem and the transportation problem; here $A$ is the node-arc incidence matrix of a directed graph, $b$ is the supply-demand vector of a transportation problem, i.e., $\sum_{i}b_{i}=0$, and $c>0$ are the edge costs. It was soon asked whether the dynamics (5) can also solve more complex problems: [SV16b] extended the result to more general flow problems, [ZM18] conducted an extensive experimental study on general flow problems with convex and monotonically increasing objective functions, and [SV16a, BBK*+*19, FCP18a] showed that it can solve linear programs of the form

[TABLE]

where $c$ is a positive vector444 [BBK*+*19] even shows that the dynamic (5) solves (7) when $c$ is a non-negative cost vector such that $c^{T}|v|>0$ for all $v\in\mathrm{ker}(A)$.. [KKM19] generalized the latter result to the non-uniform dynamics (6). Moreover, the dynamics (6) inspired the Dynamic Monge-Kantorovich model which was introduced in [FCP18b] and was recently applied in [FD*+*20] for computing the numerical solution of the Optimal Transport Problem [Amb03, San15].

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

[TABLE]

No biological significance is claimed for this dynamics. [IJNT11, SV16c] showed that the dynamics (8) solves linear programs of the form

[TABLE]

where $c$ is a positive vector. In [BBD*+*13], convergence was claimed for the non-uniform Physarum dynamics (3), i.e.

[TABLE]

Only a proof sketch was given; the claim has been withdrawn by the authors [BBD*+*20]. Here, we show convergence using a Lyapunov-based argument.

3 Existence of a Solution

As mentioned above, we show that the dynamics (3) with starting point $x(0)\in G$ has domain $[0,\infty)$ and stays in $G$. Moreover, for the limit points of the dynamics, the coordinates in $I$ are positive and the coordinates outside $I$ are zero. In order to talk conveniently about these limit points, we define the dynamics on a superset of $G$. Let

[TABLE]

We can define $q(x)$ for any $x\in G^{*}$:

[TABLE]

be the minimum energy solution $q(x)$ with respect to the resistances $r_{i}=c_{i}/x_{i}$ and with support contained in ${\mathit{supp}}(x)$. The following characterization of $q(x)$ is well-known. We use $R=R(t)$ to denote the diagonal matrix with diagonal entries $c_{i}/x_{i}$.

Lemma 1**.**

For $x\in G^{*}$, let $B={\mathit{supp}}(x)$ and $N=[m]\setminus B$. Obtain $R_{B}$ by deleting the rows and columns corresponding to $N$, obtain $A_{B}$ from $A$ by first deleting the columns in $N$ and then keeping a maximum set of independent rows, let $b_{B}$ be the restriction of $b$ to these rows, and let $p_{B}$ be a vector of the same dimension as $b_{B}$. Then $q(x)$ is uniquely determined by $q=(q_{B},q_{N})$, $q_{N}=0$, $q_{B}=R_{B}^{-1}A_{B}^{T}p_{B}$, where $p_{B}=(A_{B}R_{B}^{-1}A_{B})^{-1}b_{B}$. Defining $p_{h}=0$ for $h\not\in B$, we have the equalities $b^{T}p=q^{T}Rq$.

Proof.

For simplicity, we assume ${\mathit{supp}}(x)=[m]$. By the KKT conditions for convex optimization, the optimal solution $q$ must satisfy $Rq=A^{T}p$ for some $p$. Since $Aq=b$ this implies, $AR^{-1}A^{T}p=b$ and $q=R^{-1}A^{T}p$. The matrix $AR^{-1}A^{T}$ is non-singular since $R$ is a diagonal matrix with positive entries and $A$ has full row rank. Then, $b^{T}p=b^{T}_{B}p_{B}=(A_{B}q_{B})^{T}p_{B}=q_{B}^{T}R_{B}R_{B}^{-1}A_{B}^{T}p_{B}=q_{B}^{T}R_{B}q_{B}=q^{T}Rq$. ∎

The Physarum dynamics for $x\in G^{*}$ becomes

[TABLE]

We show existence of a solution with domain $[0,\infty)$ for each initial point $x^{0}\in G^{*}$. Our argument is based on [SV16c]. Let M=\max\{\;\text{absolute value of the determinant of a square submatrix of A}\;\}.

Fact 1** (Lemma 3.3 in [BBK*+*19]).**

For every $x\in G^{*}$: $\|q(x)\|_{\infty}\leq\beta\mathop{:=}M\|b\|_{1}$.

Fact 2** (Lemma 5.2 in [SV16c]).**

Let $x\geq 0$, $B={\mathit{supp}}(x)\not=\emptyset$, $L_{B}=A_{B}R_{B}^{-1}A_{B}^{T}$ and let $A_{i}$ be the $i$-th column of $A_{B}$. For every $i\in B$ it holds that

[TABLE]

Lemma 2**.**

Suppose $x:[0,T)\mapsto G$ is a solution to (3) for some $T>0$. Let $B={\mathit{supp}}(x(0))$. Then

1. a)

$x_{i}(t)\leq\max(x_{i}(0),\beta)$* for $i\in[m]$ and $t\in[0,T)$.* 2. b)

$x(T)=\lim_{t\rightarrow T^{-}}x(t)$* exists.* 3. c)

$x_{i}(T)>0$* for every $i\in I$.* 4. d)

$x_{i}(T)>0$* for every $i\in B$.*

Proof.

Let $i\in B$ be arbitrary.

1. a)

By Fact 1, $\dot{x}_{i}=d_{i}(q_{i}-x_{i})\leq d_{i}(\beta-x_{i})$. Let $y_{i}(t)=(x_{i}(t)-\beta)/d_{i}$. Then $\dot{y}_{i}\leq-d_{i}y_{i}(t)$ and hence $y_{i}(t)\leq y_{i}(0)e^{-d_{i}t}$ by Gronwall’s Lemma. Thus

[TABLE] 2. b)

We know that $x_{i}$ is differentiable in the interval $[0,T)$ and

[TABLE]

Since every function that has bounded first derivatives is Lipschitz, the limit $x(T)=\lim_{t\rightarrow T^{-}}x(t)$ exists. 3. c)

Let $y$ be an optimal solution with ${\mathit{supp}}(y)=I$. Note that $I\subseteq{\mathit{supp}}(x(0))$. Consider the barrier function

[TABLE]

Since $x_{j}(t)$ is bounded by part a), so is $W(t)$. Suppose for a contradiction that $\inf_{t\in[0,T)}x_{i}(t)=0$ for some $i\in I$. Then, it follows that $\inf_{t\in[0,T)}W(t)=-\infty$, implying that $\inf_{t\in[0,T)}\dot{W}(t)=-\infty$. However, since

[TABLE]

we obtain a contradiction to

[TABLE] 4. d)

Let $y$ be an optimal solution to (1) with ${\mathit{supp}}(y)=I$. By part c), we have $x_{i}(t)>0$ for every $i\in{\mathit{supp}}(y)$ and every $t\in[0,T]$. Set

[TABLE]

and note that

[TABLE]

Further, using Fact 2, observe that for every $t\in[0,T]$ we have

[TABLE]

Let $N:=\lVert c\rVert_{1}\cdot M/\epsilon$ and $N^{\prime}:=\lVert d\rVert_{\infty}[\frac{N}{\min_{i}\{c_{i}\}}+1]$ Then, for every $i\in B$ we have

[TABLE]

Hence, by Gronwall’s Lemma we obtain for every $t\in[0,T]$ and $i\in B$ that

[TABLE]

∎

We use the following standard result from dynamic systems.

Fact 3**.**

Let $G\subseteq R^{m}$ be an open set, $s\in G$ and $f:G\rightarrow\mathbb{R}^{m}$ be a $C^{1}$ function (continuously differentiable function). Consider the following dynamical system: $\frac{d}{dt}x(t)=f(x(t))$ with initial condition $x(0)=s$. Suppose it satisfies the following conditions for every solution $x:[0,T)\rightarrow G$ with $T\in(0,\infty)$: The limit $x(T)=\lim_{t\rightarrow T^{-}}x(t)$ exists and lies in $G$. Then, there exists a global solution $x:[0,\infty)\rightarrow G$.

We can now prove the existence of a global solution.

Theorem 4** (Existence of Solution).**

For every initial condition $x(0)\in G^{*}$, the non-uniform Physarum Dynamics (3) has a unique solution $x:[0,\infty)\mapsto G^{*}$. The solution satisfies ${\mathit{supp}}(x(t))={\mathit{supp}}(x(0))$ for all $t$.

Proof.

Let $B={\mathit{supp}}(x(0))$. There is a local solution $x(t)$ in an interval $[0,\epsilon)$ for some $\epsilon>0$. The solution satisfies $x_{i}(t)>0$ for $i\in B$ and $x_{i}(t)=0$ for $i\not\in B$. The latter follows from the definition (10) of the dynamics. On $G^{*}_{B}$ the function $x\mapsto D(q(x)-x)$ used in (10) is defined by a vector of rational functions in $x$ and hence is continuously differentiable. Also $G_{B}^{*}$ restricted to the coordinates in $B$ is an open set. Thus Fact 3 is applicable and the existence of a global solution follows from Lemma 2. ∎

4 The Lyapunov Function

Recall the definition of the Lyapunov function (4). Let $x^{*}$ be an optimal solution to (1) with $I={\mathit{supp}}(x^{*})$ and consider

[TABLE]

Theorem 5** (Lyapunov function).**

Let $x\in G^{*}$. Then

[TABLE]

Moreover, $\dot{V}(x)=0$ if and only if $x$ is an optimal solution to (1).

Proof.

Since $V$ is continuously differentiable, $\dot{V}(x)=\nabla V\circ D(q(x)-x)$, where $\nabla V$ is the gradient of $V$ (vector of partial derivatives of $V$ with respect to the $x_{i}$’s) [LaS76, page 30]. An easy computation yields

[TABLE]

where the second equality follows by $r_{i}(t)=c_{i}/x_{i}(t)$ for $i\in I$, $q(t)=R^{-1}(t)A^{T}p(t)$, $Ax^{*}=b$ and

[TABLE]

By combining (12) and (13) we obtain

[TABLE]

Since $q(t)=R^{-1}(t)A^{T}p(t),$ it holds that

[TABLE]

and using the Cauchy Schwarz inequality and the fact that $\sqrt{a}\cdot\sqrt{b}\leq(a+b)/2$, we obtain

[TABLE]

and further by combining (14) and (15),

[TABLE]

We have $\dot{V}(x)=0$ if and only the inequalities (1), (2), and (3) are equalities. (1) is an equality iff $q(x)$ and $x$ are parallel, i.e., $q(x)=\alpha\cdot x$ for some $\alpha\geq 0$. (2) is tight iff

[TABLE]

i.e., iff $\alpha=1$ and hence $q(x)=x$. (3) is tight iff $c^{T}x=c^{T}x^{*}$. Finally, if $q(x)=x$ then $Ax=Aq(x)=b$. Thus, $\dot{V}(x)=0$ if and only if $x$ is an optimal solution to (1). ∎

Theorem 6**.**

Let $x(0)\in G^{*}$. Then there is an $\epsilon>0$ depending on $x(0)$, $A$, $b$, $c$ such that $x_{h}(t)>\epsilon$ for all $h\in I$ and all $t$.

Proof.

Let $x^{*}$ be an optimal solution to (1) with ${\mathit{supp}}(x^{*})=I$. Let $x^{*}_{\mathrm{max}}=\max_{i}x^{*}_{i}$, $x^{*}_{\mathrm{min}}=\min_{i\in I}x^{*}_{i}$, $d_{\mathrm{max}}=\max_{i}d_{i}$, $d_{\mathrm{min}}=\min_{i}d_{i}$, $c_{\max}=\max_{i}c_{i}$, $c_{\min}=\min_{i}c_{i}$ and $\delta=\min(1,\min_{i\in I}x_{i}(0))$. Since $\dot{V}(x)\leq 0$ for all $t$, $V(x(t))\leq V(x(0))$ for all $t$. Then, for any $h\in I$ we have

[TABLE]

and hence

[TABLE]

for all $t$. The upper bound on $\ln(1/x_{h}(t))$ does not depend on $t$. This proves the claim. ∎

Corollary 1**.**

Assume $x(0)\in G^{*}$. The positive limit points of the trajectory $x(t)$ are contained in $G^{*}$.

Lemma 3**.**

The right hand side of (10) is continuous.

Proof.

Let $(x^{(n)})$ be a sequence of points in $G^{*}$ converging to $x\in G^{*}$. We need to show $f(x^{(n)})\rightarrow f(x)$. For $i\in{\mathit{supp}}(x)$ we must have $i\in{\mathit{supp}}(x^{(n)})$ for all sufficiently large $n$. Let $B$ be any superset of ${\mathit{supp}}(x)$ and consider the subsequence of $(x^{(n)})$ with ${\mathit{supp}}(x^{(n)})=B$. We reuse $(x^{(n)})$ to denote the subsequence. If the subsequence is finite, there is nothing to show. So assume the subsequence is infinite. We need to show $f_{i}(x^{(n)})\rightarrow f_{i}(x)$. The proof proceeds by case distinction.

Case 1: For $i\in{\mathit{supp}}(x)$, we clearly have $f_{i}(x^{(n)})\rightarrow f_{i}(x)$ for $n\rightarrow\infty$, since $q_{i}$ is defined by a rational function (Lemma 1).

Case 2: For $i\in B\setminus{\mathit{supp}}(x)$, we need to prove $f_{i}(x^{(n)})\rightarrow 0$. For this it suffices to show $q_{i}(x^{(n)})\rightarrow 0$. Let $x^{*}$ be an optimal solution to (1) with ${\mathit{supp}}(x^{*})=I$. We may assume that $(x^{(n)})_{i}\geq x_{i}/2$ for all $i\in I$ and all $n$. Let

[TABLE]

Hence, $x^{*}\leq\epsilon\cdot(x^{(n)})$ and thus by Fact (2) we have

[TABLE]

By the proof of Lemma 2 Part d) and by Lemma 1 it holds that

[TABLE]

Since $(x^{(n)})_{i}\rightarrow x_{i}=0$ for all $i\in[m]\backslash{\mathit{supp}}(x)$, we have $|q_{i}(x^{(n)})|\rightarrow 0$ for all $i\in B\backslash{\mathit{supp}}(x)$. ∎

5 Putting it Together

We complete the proof of Theorem 2. Let $x(t)$, $t\in[0,\infty)$ be the global solution of (10) with $x(0)=x^{0}\in G^{*}$. Then $x(t)\in G^{*}$ for all $t\in[0,\infty)$ by Theorem 4. Suppose there is a sequence $t_{n}$ with $t_{n}\rightarrow\infty$ and $x(t_{n})\rightarrow p$ as $n\rightarrow\infty$. Then $p\in G^{*}$ by Corollary 1. Now $V(x(t_{n}))\rightarrow V(p)$ and $\dot{V}(x(t_{n}))\rightarrow\dot{V}(p)$. The latter equality holds since $V(x)$ is $C^{1}$ and hence $\dot{V}(x)=\nabla V\circ f(x)$, and since $f$ is continuous by Lemma 3.

However, $\dot{V}(x(t))\leq 0$ and $V(x(t))\geq V(p)$ for all $t\in[0,\infty)$ because $V(x(t))$ is decreasing. So $\dot{V}(x(t))\rightarrow 0$ as $t\rightarrow\infty$. Hence, by uniqueness of limits $\dot{V}(p)=0$.

Finally $\dot{V}(p)=0$ implies that $p$ is an optimal solution to (1) (Theorem 5).

6 Approach to Optimum: A Simple Example

We study the nonuniform dynamics for the simple example of Figure 1 in more detail. We consider “minimize $c_{1}x_{1}+c_{2}x_{2}$, $x_{1}+x_{2}=1$, $x\geq 0$, where $c_{2}>c_{1}$” and the Physarum dynamics $\dot{x_{i}}=d_{i}(q_{i}-x_{i})$, $i=1,2$. The unique optimum is $(1,0)$. We ask how the dynamics enters the optimal point?

Assume, we are in the point $(x_{1},x_{2})$. The conductance of edge $i$ is $x_{i}/c_{i}$. The conductance of the system of parallel edges is $x_{1}/c_{1}+x_{2}/c_{2}=(x_{1}c_{2}+x_{2}c_{1})/c_{1}c_{2}$. Therefore the potential difference is $\Delta=c_{1}c_{2}/(x_{1}c_{2}+x_{2}c_{1}$ and hence with $N=x_{1}c_{2}+x_{2}c_{1}$

[TABLE]

and

[TABLE]

We know that the dynamics converges to the point $(1,0)$ and are interested in the behavior near the limit point. Therefore let $(x_{1},x_{2})=(1-\epsilon_{1},\epsilon_{2})$. Then

[TABLE]

For the approximation, we only kept the terms linear in the epsilons and ignored all higher powers. This is justified since we are interested in the behavior for small $\epsilon_{1}$ and $\epsilon_{2}$. We continue with the approximations.

From the second equation, we conclude

[TABLE]

From the first equation, we conclude

[TABLE]

Thus $\epsilon_{1}(t)\geq e^{-d_{1}t}$. If $d_{1}<(c_{2}-c_{1})d_{2}/c_{2}$ or equivalently $d_{2}\geq d_{1}c_{2}/(c_{2}-c_{1})$, the rate at which $\epsilon_{2}$ goes to zero is higher than the rate for $\epsilon_{1}$ and hence the trajectory converges to the $x$-axis and enters the optimal point horizontally.

So assume $d_{1}>(c_{2}-c_{1})d_{2}/c_{1}$. We try the “Ansatz” $\epsilon_{1}(t)=e^{-d_{1}t}f(t)$. Then $\dot{\epsilon}_{1}=e^{-d_{1}t}\dot{f}(t)-d_{1}e^{-d_{1}t}f(t)$ and hence

[TABLE]

Thus

[TABLE]

and hence

[TABLE]

We conclude that $\epsilon_{1}$ and $\epsilon_{2}$ go to zero at the same rate, and the trajectory follows a straight line with slope

[TABLE]

For $d_{1}=(c_{2}-c_{1})d_{2}/c_{2}$, the analysis is inconclusive. However, since the behavior should be continuous in the $d$’s, it is natural to conjecture that the trajectories converge to the $x$-axis.

Even less formally, this result can also be obtained as follows. Assume that the trajectory starting in the point $(1-\epsilon_{1},\epsilon_{2})$ is essentially straight, we must have (observe that the straight line from $(1-\epsilon_{1},\epsilon_{2})$ to the optimal point $(1,0)$ has direction $(\epsilon_{1},-\epsilon_{2})$)

[TABLE]

for some $\alpha$. For $\epsilon_{2}>0$, we obtain $\alpha=d_{2}(c_{2}-c_{1})/c_{2}$ and further $d_{1}c_{1}\cdot\epsilon_{2}=\left[d_{1}c_{2}-d_{2}\left(c_{2}-c_{1}\right)\right]\cdot\epsilon_{1}$ from the first equation. Thus,

[TABLE]

as above. Since $\epsilon_{1}>0$, we need $d_{1}c_{2}>d_{2}(c_{2}-c_{1})$, i.e., $d_{1}>(c_{2}-c_{1})d_{2}/c_{2}$.

Let us specialize to the case $c_{1}=1$, $c_{2}=2$.

• •

If $d_{2}>c_{2}d_{1}/(c_{2}-c_{1})=2d_{1}$, the trajectories enter the optimal point horizontally. Whether a trajectory enters from the left or from the right depends on the initial point. The right plot in Figure 1 shows an example.

• •

If $d_{2}<2d_{1}$, the trajectories enter the optimal point with slope $-2+d_{2}/d_{1}$. In particular, for $d_{1}=d_{2}=1$ (uniform dynamics), the slope is $-1$. The left plot in Figure 1 shows an example.

7 Speed of Convergence

At the suggestion of an anonymous reviewer, we conduct a numerical experiment that explores the convergence-time dependence of the non-uniform dynamics (3) on the reactivity matrix $D$.

Here, we consider a min-cost flow problem with a unit demand between a single source-sink pair of nodes and parameterized costs $c^{(f)}$, see Figure 3. We compare the convergence-time of the dynamics (3) when the matrix $D$ is either: i) $\mathrm{diag}(c^{(f)})$; or ii) the identity matrix $I$ (i.e., the uniform dynamics (2)), see Figure 4.

We implement and execute the forward Euler discretization of the non-uniform dynamics (3), which reads

[TABLE]

where the step size555Note that the step size satisfies $0<h\leq 1/2$ for reactivity matrices $D\in\{C,I\}$. $h=\tfrac{1}{2\lVert c\rVert_{1}}$, the number of iterations $t=\tfrac{1}{h}\log\tfrac{\lVert c\rVert_{1}}{\epsilon}$, and the error $\epsilon=\tfrac{1}{10}$.

The experimental data in Figure 4 suggests that the non-uniform dynamics (3) with $D=\mathrm{diag}(c^{(f)})$ can achieve a substantial convergence-time improvement over the uniform dynamics (2).

We leave the convergence-time dependence of the non-uniform dynamics (3) on the reactivity matrix $D$ as a subject for future research.

8 Acknowledgements

KM wants to thank his colleague Mark Groves for a private lesson on dynamical systems. The authors want to thank the anonymous reviewers for their careful reading, constructive comments, and insightful suggestions.