Convergence of the Non-Uniform Physarum Dynamics

TL;DR

This paper proves convergence of a generalized Physarum dynamics to the optimal solution of a weighted basis pursuit problem, extending previous results from uniform cases to more general, non-uniform and functional dynamics.

Contribution

It establishes convergence for non-uniform Physarum dynamics with general reactivity functions, broadening the understanding beyond uniform and specific shortest path cases.

Findings

Convergence shown for non-uniform Physarum dynamics under general conditions.

Extended convergence results to dynamics with multiplicative factors involving functions g_e.

Demonstrated convergence for a broader class of problems beyond shortest path.

Abstract

Let $c \in \mathbb{Z}^m_{> 0}$, $A \in \mathbb{Z}^{n\times m}$, and $b \in \mathbb{Z}^n$. We show under fairly general conditions that the non-uniform Physarum dynamics \[ \dot{x}_e = a_e(x,t) \left(|q_e| - x_e\right) \] converges to the optimum solution $x^*$ of the weighted basis pursuit problem minimize $c^T x$ subject to $A f = b$ and $|f| \le x$. Here, $f$ and $x$ are $m$-vectors of real variables, $q$ minimizes the energy $\sum_e (c_e/x_e) q_e^2$ subject to the constraints $A q = b$ and $\mathrm{supp}(q) \subseteq \mathrm{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 \[ \dot{x}_e = x_e \cdot \left( g_e \left(\frac{|q_e|}{x_e}\right) - 1\right),\] 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.

Figures

Equations

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

$${\dot{x}}_e=x_e\left(g_e\left(\frac{|q_e|}{x_e}\right)-1\right),$$

$$\text{minimize }{c}^{T}x\text{ subject to }Af=b\text{ and }|f|\le x,$$

$$\dot{x}=|q|-x,$$

$$q={\mathrm{argmin}}_f\left\{\sum _e(c_e/x_e)f_e^2:Af=b\text{ and }{f}_e=0\text{ whenever }{x}_e=0\right\}.$$

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

$$L(x,t)=\sum _e(c_e/x_e)q_e^2+\sum _e{c}_e{x}_e$$

$${\dot{x}}_e=x_e\left(g_e\left(\frac{|q_e|}{x_e}\right)-1\right)\text{for all e∈[m]},$$

$$L(x,t)=\sum _e(c_e/x_e)q_e^2+\sum _e{c}_e{x}_e$$

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

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

$${\dot{y}}_e=a_e{\dot{x}}_e=a_e(|q_e(r_e)|-a_e{x}_e)=a_e(|q_e(r_e)|-y_e).$$

$$\dot{y_e}=a_e(|q_e|-y_e)$$

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

$$E_x(f)=\sum _e(c_e/x_e)f_e^2$$

$$\mathrm{cost}(f)=\sum _e{c}_e|f_e|=c^T|f|$$

$$E_x(x)=\sum _e(c_e/x_e)x_e^2=\sum _e{c}_e{x}_e=\mathrm{cost}(x).$$

$$b$$

$$q$$

$$A{R}^{-1}{A}^{T}p$$

$$2(c_e/x_e)q_e=\sum _i{p}_i{A}_{i,e}$$

$$E_x(q)=q^{T}Rq=p^{T}A{R}^{-1}R{R}^{-1}{A}^{T}p=p^{T}A{R}^{-1}{A}^{T}p=p^{T}b.$$

$${\dot{x}}_e=a_e(x,t)\left(\left|\frac{x_e}{c_e}{A}_e^{T}p\right|-x_e\right)=a_e(x,t)x_e\left(\frac{|A_e^{T}p|}{c_e}-1\right),$$

$$\mathrm{cost}(z)=\sum _{e\in \mathrm{supp}(z)}{c}_e{z}_e=\sum _{e\in \mathrm{supp}(z)}{z}_e{A}_e^{T}p=b^{T}p,$$

$$A\frac{d}{dt}(R^{-1})A^{T}p+A{R}^{-1}{A}^T\dot{p}=0$$

$$\frac{d}{dt}{p}^{T}b$$

$$=-p^{T}A\frac{d}{dt}(R^{-1})A^{T}p=-\sum _e\frac{(A_e^{T}p{)}^2}{c_e}{\dot{x}}_e$$

$$=-\sum _e(A_e^{T}p{)}^2\frac{a_e}{c_e}\left(\frac{x_e}{c_e}|A_e^{T}p|-x_e\right)$$

$$=-\sum _e{a}_e{c}_e{x}_e\left(\frac{|A_e^{T}p{|}^3}{c_e^3}-\frac{|A_e^{T}p{|}^2}{c_e^2}\right),$$

$$\frac{d}{dt}{c}^{T}x=\sum _e{c}_e{\dot{x}}_e=\sum _e{a}_e{c}_e\left(\frac{x_e}{c_e}|A_e^{T}p|-x_e\right)=\sum _e{a}_e{c}_e{x}_e\left(\frac{|A_e^{T}p|}{c_e}-1\right).$$

$$\frac{d}{dt}L(x,t)=-\sum _e{a}_e{c}_e{x}_e\left({\lambda }_e^3-{\lambda }_e^2-{\lambda }_e+1\right).$$

$${\lambda }_e^3-{\lambda }_e^2-{\lambda }_e+1=({\lambda }_e^2-1)({\lambda }_e-1)=({\lambda }_e+1)({\lambda }_e-1{)}^2$$

$$x(t)=x(0)+{\int }_0^t\frac{e^{-s}}{2}(1-x(s))ds\le \frac{1}{2}+\frac{1}{4}{\int }_0^t{e}^{-s}ds\le \frac{1}{2}+\frac{1}{4}(1-e^{-t})\le \frac{3}{4}$$

Full text

Let $c\in\mathbb{Z}^{m}_{>0}$, $A\in\mathbb{Z}^{n\times m}$, and $b\in\mathbb{Z}^{n}$. We show under fairly general conditions that the non-uniform Physarum dynamics

[TABLE]

converges to the optimum solution $x^{*}$ of the weighted basis pursuit problem minimize $c^{T}x$ subject to $Af=b$ and $|f|\leq x$. Here, $f$ and $x$ are $m$-dimensional vectors of real variables, $q$ minimizes the energy $\sum_{e}(c_{e}/x_{e})q_{e}^{2}$ subject to the constraints $Aq=b$ and $\mathrm{supp}(q)\subseteq\mathrm{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

[TABLE]

where each $g_{e}$ is an increasing differentiable function with $g_{e}(1)=1$ (satisfying some mild conditions). 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.

1 Introduction

The Physarum computing model is an analog computing model motivated by the network dynamics of the slime mold Physarum polycephalum. In wet-lab experiments, it was observed that the slime mold is apparently able to solve shortest path problems [NYT00]. A mathematical model for the dynamic behavior of the slime was proposed in [TKN07]. It models the slime network as an electrical network with time-varying resistors that react to the amount of electrical current flowing through them. A more general model for the dynamics was introduced in [NIU*+*07] to deal with situations in which the slime has to operate in a non-uniform environment. This more general dynamics is the subject of this paper. In Section 2, we give more details on the biological background and also survey the theoretical work on the Physarum dynamics.

The weighted basis pursuit problem asks to find the minimal weighted one-norm solution of a linear system. Formally,

[TABLE]

where $c\in\mathbb{Z}^{m}_{>0}$, $A\in\mathbb{Z}^{n\times m}$, $b\in\mathbb{Z}^{n}$, and $x$ and $f$ are $m$-dimensional vectors of real variables. The absolute-value operator is applied componentwise. The matrix $A$ is assumed to have full row-rank; this implies $n\leq m$. For simplicity, we also assume that any two basic feasible solutions of $Af=b$ have distinct cost111A basic feasible solution of $Af=b$ has the form $f=(f_{B},f_{\overline{B}})$, where $B$ is a subset of $[m]$ of size $n$, $\overline{B}=[m]\setminus B$, the submatrix $A_{B}$ of $A$ is invertible, $f_{B}=A_{B}^{-1}b$, and $f_{\overline{B}}=0$. The cost of such a solution is $c^{T}|f|$.; in particular, the optimal solution $(f^{*},x^{*})$ to (1) is unique. We index the rows of $A$ by $i$ and the columns of $A$ by $e$ and, for historical reasons (see Subsection 2.1), refer to the rows as nodes and the columns as edges.

The Physarum dynamics evolves a vector $x(t)\in\mathbb{R}^{m}_{\geq 0}$ according to the dynamics

[TABLE]

where $q$ is the minimum energy feasible solution of $Af=b$ according to the resistances $r_{e}=c_{e}/x_{e}$:

[TABLE]

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.

In [BBK*+*19] and [FCP18] it was shown that the Physarum dynamics (2) can solve the weighted basis pursuit problem (1).

Theorem 1** ([BBK*+*19, FCP18]).**

Assume a strictly positive starting vector $x(0)\in R^{m}_{>0}$. Then the solution $x(t)$ to the dynamics (2) is defined on $[0,\infty)$, $x(t)>0$ for all $t$, and $x(t)$ and $|q(t)|$ converge to the optimal solution $x^{*}$ of (1).

Actually, the papers [BBK*+*19, FCP18] show convergence under the more general condition $c\geq 0$ and $c^{T}|f|>0$ for any $f$ in the kernel of $A$, but here we do not need this generality.

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.,

[TABLE]

where $a_{e}(x,t)\geq 0$ is the reactivity of edge $e$ at time $t$, i.e., the edges no longer react uniformly to differences between $|q|$ and $x$, but the reactivity depends on the edge, the current state, and the time. We refer to (3) as the non-uniform Physarum dynamics. The special case that $a_{e}(x,t)$ is a positive constant for each edge was introduced in [NIU*+*07] to model the behavior of Physarum polycephalum in non-uniform environments; see Subsection 2.2.

In Section 3, we prove our main technical contribution for the dynamics (3):

Theorem 2**.**

Assume $x(0)>0$, $0\leq a_{e}(x,t)\leq C$ for all $e$, $x$, and $t$ and some constant $C$, and $a_{e}(x,t)$ is Lipschitz-continuous. Then:

1. (a)

The dynamics (3) has a unique solution $x(t)>0$ for $t\in[0,\infty)$. 2. (b)

The function

[TABLE]

is a Lyapunov function for the dynamics (3), i.e., $\frac{d}{dt}L(x,t)\leq 0$ for all $t\in[0,\infty)$. Moreover, $\frac{d}{dt}L(x,t)=0$ if and only if for all $e$: either $a_{e}(x,t)=0$ or $x_{e}(t)=0$ or $|q_{e}|=x_{e}$. 3. (c)

If, in addition, $a_{e}(x,t)\geq\epsilon$ for some positive $\epsilon$ and all $e$, $x$, and $t$, then $\bar{x}\geq 0$ is a fixed point of (3) if and only if $\bar{x}=|f|$ for a basic feasible solution of $Af=b$. 4. (d)

Under the same additional assumption as in (c), $x(t)$ and $|q(t)|$ converge to a fixed point of (3) as $t$ goes to infinity. 5. (e)

If, in addition, $a_{e}(x,t)$ does not depend on $x$ and $\frac{d}{dt}{a_{e}}(t)\leq 0$ for all $e$ and $t$, then $x(t)$ and $|q(t)|$ converge to $x^{*}$ as $t$ goes to infinity. In particular, this holds true if $a_{e}(t)$ is a positive constant for all $e$.

The proof of part (a) is standard and part (c) was shown in [BBK*+*19]. The Lyapunov function in part (b) was introduced in [FDCP18]. In  [FCP18] it was shown to be a Lyapunov function for the uniform case, i.e., $a_{e}(x,t)=1$ for all $e$, $x$, and $t$. We observe that the Lyapunov function also works for the non-uniform dynamics. Part (d) follows easily from parts (b) and (c). Finally, the proof of part (e) is inspired by [BBK*+*19].

The function $L(x,t)$ is a Lyapunov function for the Physarum dynamics under very general conditions. Essentially, the only requirement is that $\dot{x}_{e}$ has the same sign as $|q_{e}|-x_{e}$. For the existence of a solution with domain $[0,\infty)$, we also need that $\dot{x}_{e}/(|q_{e}|-x_{e})$ is bounded. For the convergence to a fixed point, we need in addition that $\dot{x}_{e}/(|q_{e}|-x_{e})$ is bounded away from zero.

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

[TABLE]

where the response function $g_{e}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is assumed to be an increasing differentiable function satisfying $g_{e}(1)=1$. Bonifaci introduced this model in [Bon17] 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., [Bon17] 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.

In Section 4, we prove our main technical contribution for the dynamics (5):

Theorem 3**.**

Assume $g_{e}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is an increasing and differentiable function satisfying $g_{e}(1)=1$, for all $e$. Then,

1. (a)

The dynamics (5) has a unique solution $x(t)>0$ for $t\in[0,\infty)$. 2. (b)

$x\geq 0$* is a fixed point of (5) if $x=|f|$ for a basic feasible solution of (1).* 3. (c)

The function

[TABLE]

is a Lyapunov function for the dynamics (5), i.e., $\frac{d}{dt}L(x,t)\leq 0$ for all $t\in[0,\infty)$. Moreover, $\frac{d}{dt}L(x,t)=0$ if and only if $x$ is a fixed point. 4. (d)

$x(t)$* and $|q(t)|$ converge to a fixed point of (5).* 5. (e)

If, in addition, $g_{e}(y)\geq 1+\alpha(y-1)$ for some $\alpha>0$ and all $e$ and $y$, then $x(t)$ and $|q(t)|$ converge to $x^{*}$ as $t$ goes to infinity.

The proof of part (a) is standard and part (b) was shown in [BBK*+*19]. The Lyapunov function in part (c) was introduced in [FDCP18]. We observe that it also applies to the dynamics (5). Part (d) follows easily from part (c). Finally, the proof of part (e) is inspired by [BBK*+*19]. Theorem 3 also holds when the function $g_{e}$ depends on the time and the state.

Nature does not compute exactly, i.e., one should not expect that in a biological system $\dot{x}_{e}$ is exactly equal to $|q_{e}|-x_{e}$ or to $x_{e}(g_{e}(|q_{e}|/x_{e})-1)$. Rather, there will be a noise. Our results show that the dynamics (3) and (5) are fairly robust against noise, i.e., variations in $a_{e}(x,t)$ and $g_{e}(y)$.

The rest of the paper is organized as follows. In Section 2, we review the biological background and related work. In Section 3, we prove Theorem 2. 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 various optimization problems (see [Ada10] for a survey of Physarum computations), in particular the shortest path problem. 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 $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. In equilibrium ($\dot{x}_{e}=0$ for all $e$), the flow through any edge is equal to its capacity. In non-equilibrium, the capacity grows (shrinks) if the absolute value of the flow is larger (smaller) than the capacity. It is well-known that the electrical flow $q$ is the feasible flow minimizing energy dissipation $\sum_{e}r_{e}q_{e}^{2}$ (Thomson’s principle).

2.2 Minimum Risk Paths

In [NIU*+*07], 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 [NIU*+*07], the dynamics

[TABLE]

is proposed. The constant $a_{e}$ is the decay rate of edge $e$ if there is no flow on it. To model the experiment, $a_{e}=1$ for edges in the dark part of the plate, and $a_{e}=C>1$ for the edges in the lighted area, where $C$ is a constant. Nakagaki et al. [NIU*+*07] report that in computer simulations, the dynamics (8) converges to the shortest source-sink path with respect to the modified cost function $a_{e}c_{e}$.

2.3 A Reformulation: Nonuniform Physarum

Let $y_{e}=a_{e}x_{e}$. The electrical flow $q$ is determined by the resistances $r_{e}=c_{e}/x_{e}$. Therefore, we write $q(r(t))$ instead of $q(t)$ for clarity. Next observe that $r_{e}=c_{e}/x_{e}=(a_{e}c_{e})/(y_{e})$. Thus if we take $y$ as the vector of edge capacities and $(a_{e}c_{e})_{e}$ as the vector of costs, we get the same electrical flow. We can express (8) as a dynamics for $y$ as

[TABLE]

So we may instead consider the dynamics

[TABLE]

under the modified cost function $a_{e}c_{e}$. This is our dynamics (3), where we generalized further by allowing $a_{e}$ to depend on $x$ and $t$. In this model, the quantity $a_{e}$ indicates the responsiveness (reactivity) of an edge to differences between flow and capacity.

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 $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. Convergence for the discretization of (2) was shown in [BBD*+*13].

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 [BBK*+*19]. The function (4) was introduced in [FDCP18] and shown to be a Lyapunov function for (2) in [FCP18].

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

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

[TABLE]

No biological significance is claimed for this dynamics. It can solve linear programs with positive cost vectors [IJNT11, SV16b]. In [BBD*+*13], convergence was claimed for the non-uniform dynamics $\dot{x}_{e}=a_{e}(q_{e}-x_{e})$. The proof is incorrect [BBD*+*19].

3 The Proof of Theorem 2

3.1 Preliminaries

For a capacity vector $x\geq 0$ and a vector $f\in\mathbb{R}^{m}$ with $\mathrm{supp}(f)\subseteq\mathrm{supp}(x)$, we use

[TABLE]

to denote the energy of $f$. When $\mathrm{supp}(f)\not\subseteq\mathrm{supp}(x)$, the energy of $f$ is infinite. Further, we use

[TABLE]

to denote the cost of $f$. Note that

[TABLE]

We use $R$ to denote a diagonal matrix with entries $c_{e}/x_{e}$; here we use the convention that attention is restricted to the edges $e$ with $x_{e}>0$. In part (a) of Theorem 2, it is shown that $x(t)>0$ for all $t$ if $x(0)>0$. However, in the limit some edges may have capacity zero. Energy-minimizing solutions are induced by node potentials $p\in\mathbb{R}^{n}$ according to the following equations:

[TABLE]

We give a short justification. The vector $q$ minimizes the quadratic function $\sum_{e}(c_{e}/x_{e})q_{e}^{2}$ subject to the constraints $Aq=b$. The KKT conditions (see [BV04, Subsection 5.5]) state that at the optimum, the gradient of the objective is a linear combination of the gradients of the constraints. Thus

[TABLE]

for some vector $p\in\mathbb{R}^{n}$. Absorbing the factor $2$ into $p$ 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 (ii) can, for example, be found in [BBK*+*19]. Let $D$ be the maximum absolute value of a square submatrix of $A$.

1. (i)

The minimum energy solution is defined by (11) and (12). Moreover, it is unique. 2. (ii)

$|q_{e}|\leq D\|b\|_{1}$ for every $e\in[m]$. 3. (iii)

$E_{x}(q)=\sum_{e}(c_{e}/x_{e})q_{e}^{2}=b^{T}p$, where $p$ is defined by (12). This holds since

[TABLE]

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

[TABLE]

where $A_{e}$ denotes the $e$-th column of matrix $A$.

3.2 Existence

The right-hand side of (3) is locally Lipschitz-continuous in $x$ and $t$. The function $a_{e}(x,t)$ is locally Lipschitz by assumption, $q$ is an infinitely often differentiable rational function in the $x_{e}$ and hence locally Lipschitz. Furthermore, locally Lipschitz-continuous functions are closed under additions and multiplications. Thus $x(t)$ is defined and unique for $t\in[0,t_{0})$ for some $t_{0}$.

Since $a_{e}(x,t)\leq C$ for all $e$, $x$ and $t$, we have $\dot{x}_{e}\geq-Cx$ and thus $x_{e}\geq x_{e}(0)e^{-Ct}$. Hence, $x(t)>0$ for all $t$ and the solution does not reach the boundary of the domain in finite time. Also since $|q_{e}(t)|\leq D\|b\|_{1}$ for all $e$ and $t$, we have $\dot{x}_{e}\leq C(D\|b\|_{1}-x)$ and hence $x_{e}(t)\leq\max(x_{e}(0),D\|b\|_{1})$ for all $t$. In particular, the solution is bounded. Thus, $t_{0}=\infty$ by well-known results of maximal solutions of ordinary differential equations [Har02, Corollary 3.2].

The condition $a(x,t)\leq C<\infty$ is crucial for existence. Let $n=0$, $m=1$ and $a(x,t)=1/x$. The matrix $A$ is $0\times 1$, i.e., there are no constraints. Then the minimum energy solution is the null-vector of dimension one and (3) becomes $\dot{x}=1/x\cdot(0-x)=-1$; the domain of definition is $[0,x(0))$.

3.3 Fixed Points

A point $x$ is a fixed point if $\dot{x}=0$. In  [BBK*+*19] is was shown that the fixed points of (2) are the vectors $|f|$, where $f$ 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 $a_{e}(x,t)\geq\epsilon$ for all $e$, $x$ and $t$ and some positive $\epsilon$. Under this additional assumption $\dot{x}=0$ is equivalent to $|q|=x$ for (2) and (3). This section is reprinted from [BBK*+*19] with minor adaptions. A vector $f^{\prime}$ is sign-compatible with a vector $f$ (of the same dimension) if $f^{\prime}_{e}\not=0$ implies $f^{\prime}_{e}f_{e}>0$. In particular, $\mathrm{supp}(f^{\prime})\subseteq\mathrm{supp}(f)$. We use the following corollary of the finite basis theorem for polyhedra.

Lemma 1**.**

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

Proof.

We may assume $f\geq 0$. Otherwise, we flip the sign of the appropriate columns of $A$. Thus, the system $Af=b,\ f\geq 0$ is feasible and $f$ is the sum of a convex combination of at most $n$ basic feasible solutions plus a vector in the kernel of $A$ 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 $f$. ∎

Lemma 2**.**

Assume $a_{e}(x,t)\geq\epsilon$ for some positive $\epsilon$ and all $e$, $x$, and $t$, and that no two feasible solutions of $Af=b$ have the same cost. If $f$ is a basic feasible solution of (1), then $x=|f|$ is a fixed point. Conversely, if $x$ is a fixed point, then $x=|f|$ for some basic feasible solution $f$.

Proof.

Let $f$ be a basic feasible solution, let $x=|f|$, and let $q$ be the minimum energy feasible solution with respect to the resistances $c_{e}/x_{e}$. We have $Aq=b$ and $\mathrm{supp}(q)\subseteq\mathrm{supp}(x)$ by definition of $q$. Since $f$ is a basic feasible solution there is a subset $B$ of size $n$ of the columns of $A$ such that $A_{B}$ is non-singular and $f=(A_{B}^{-1}b,0)$. Since $\mathrm{supp}(q)\subseteq\mathrm{supp}(x)=\mathrm{supp}(f)\subseteq B$, we have $q=(q_{B},0)$ for some vector $q_{B}$. Thus, $b=Aq=A_{B}q_{B}$ and hence $q_{B}=f_{B}$. Therefore $\dot{x}=|q|-x=0$ and $x$ is an fixed point.

Conversely, if $x$ is an fixed point, $|q_{e}|=x_{e}$ for every $e$. By changing the signs of some columns of $A$, we may assume $q\geq 0$. Then $q=x$. Since $q_{e}=(x_{e}/c_{e})A_{e}^{T}p$ by (11), we have $c_{e}=A_{e}^{T}p$, whenever $x_{e}>0$. By Lemma 1, $q$ is a convex combination of basic feasible solutions plus a vector in the kernel of $A$ that are sign-compatible with $q$. The vector in the kernel is zero since $q$ is a minimum energy solution333Assume $q=q^{1}+q^{2}$ with $q^{1}\geq 0$, $q_{2}\geq 0$, $q_{2}\not=0$, and $Aq_{2}=0$. Then $Aq_{1}=b$, $\mathrm{supp}(q_{1})\subseteq\mathrm{supp}(q)\subseteq\mathrm{supp}(x)$, and $E_{x}(q_{1})<E_{x}(q)$, a contradiction.. For any basic feasible solution $z$ contributing to $q$, we have $\mathrm{supp}(z)\subseteq\mathrm{supp}(x)$. Summing over the $e\in\mathrm{supp}(z)$, we obtain

[TABLE]

i.e., all basic feasible solutions used to represent $q$ have the same cost. Since we assume the costs of distinct basic feasible solutions to be distinct, $q$ is a basic feasible solution. ∎

Corollary 1**.**

Assume $a_{e}(x,t)\geq\epsilon$ for some positive $\epsilon$ and all $e$, $x$, and $t$ and that no two feasible solutions of $Af=b$ have the same cost. Then the set of fixed points is a discrete set.

3.4 The Lyapunov Function

Lemma 3**.**

$L(x,t)=p^{T}b+c^{T}x$* is a Lyapunov function for (13). More precisely, $\frac{d}{dt}L(x,t)\leq 0$ always with equality only if for all $e$ either $x_{e}=0$ or $a_{e}(x,t)=0$ or $|A_{e}^{T}p|=c_{e}$.*

Proof.

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

[TABLE]

We next compute the derivative of both summands of $L(x,t)$ with respect to time separately. For the first summand we obtain

[TABLE]

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

[TABLE]

Combining (16) and (17), and writing $\lambda_{e}$ instead of $|A_{e}^{T}p|/c_{e}$, yields

[TABLE]

Since

[TABLE]

and $\lambda_{e}\geq 0$, $\frac{d}{dt}L(x,t)\leq 0$ always. Moreover, the derivative is equal to zero only if $a_{e}x_{e}(\lambda_{e}-1)=0$ for all $e$, i.e., for all $e$ either $x_{e}=0$ or $a_{e}(x,t)=0$ or $|A_{e}^{T}p|=c_{e}$. ∎

Corollary 2**.**

Assume further $a_{e}(x,t)\geq\epsilon$ for some positive $\epsilon$ and all $e$, $x$ and $t$. Then $L(x,t)=0$ if and only if $x$ is a fixed point.

Proof.

We have $L(x,t)=0$ if and only if for all $e$ either $x_{e}=0$ or $|A_{e}^{T}p|=c_{e}$. The latter condition is equivalent to $|q_{e}|=x_{e}/c_{e}|A_{e}^{T}p|=x_{e}$. Thus $|q|=x$. ∎

3.5 Convergence

From now on, we make the additional assumption that $a_{e}(x,t)\geq\epsilon$ for some positive $\epsilon$ and all $e$, $x$, and $t$. It then follows from the general theory of dynamical systems that $x(t)$ converges to a fixed point.

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

Assume further $a_{e}(x,t)\geq\epsilon$ for all $e$, $x$ and $t$. As $t\rightarrow\infty$, $x(t)$ and $|q(t)|$ approach a fixed point $x_{0}$. Moreover, $E_{x}(q)$ and $\mathrm{cost}(x)$ converge to $c^{T}x_{0}$.

Proof.

The proof in [Bon13] carries over. We include it for completeness. The existence of a Lyapunov function $L$ implies by [LaS76, Corollary 2.6.5] that $x(t)$ approaches the set $\left\{\,x\in\mathbb{R}_{\geq 0}^{m}\,:\,\dot{L}=0\,\right\}$, which by Corollary 2 is the same as the set $\left\{\,x\in\mathbb{R}_{\geq 0}^{m}\,:\,\dot{x}=0\,\right\}$. Since this set consists of isolated points (Lemma 2), $x(t)$ must approach one of those points, say the point $x_{0}$. When $x=x_{0}$, one has $E_{x}(q)=E_{x}(x)=\mathrm{cost}(x)=c^{T}x$. ∎

The assumption $a_{e}(x,t)\geq\epsilon>0$ is crucial as the following example shows. Let $n=m=1$, consider the task of minimizing $|x|$ subject to the constraint $x=1$, and let $a(x,t)=e^{-t}/2$ and $x(0)=1/2$. Then $\dot{x}=e^{-t}(1-x)/2$. Integrating from [math] to $t$ and observing that $x(t)\geq 1/2$ for all $t$, we obtain

[TABLE]

and hence the dynamics does not converge to the optimal solution $x^{*}=1$, which, in this case, is the only fixed point.

It remains to exclude that $x(t)$ converges to a non-optimal fixed point. We can do so under an additional assumption on $a(x,t)$.

Theorem 4**.**

Assume further that $a_{e}(x,t)$ does not depend on $x$, i.e., $a_{e}(x,t)=a_{e}(t)$, $a_{e}(t)\geq\epsilon$ for some positive $\epsilon$ for all $e$ and $t$, and $\dot{a}_{e}(t)\geq 0$ for all $e$ and $t$. As $t\rightarrow\infty$, $x(t)$ converges to the optimal solution $x^{*}$.

Proof.

Assume that $x(t)$ converges to a non-optimal fixed point $z$. Let $x^{*}$ be the optimal solution, let $B$ be such that $x_{e}(t)\leq B$ for all $e$ and $t$ (by Subsection 3.2 the solution is bounded), and let

[TABLE]

Let $\delta=(\mathrm{cost}(z)-\mathrm{cost}(x^{*}))/2$. Then $E_{x}(q(t))\geq\mathrm{cost}(z)-\delta=\mathrm{cost}(x^{*})+\delta$, for all sufficiently large $t$. Further, by definition $q_{e}=(x_{e}/c_{e})A_{e}^{T}p$ and thus

[TABLE]

where the first inequality follows from $\ln(x_{e}/B)\leq 0$ and $\dot{a}_{e}\geq 0$, and the second inequality is due to

[TABLE]

Hence $W\rightarrow\infty$, a contradiction to the fact that $x$ is bounded. ∎

4 Bonifaci’s Refined Model

Bonifaci [Bon17] investigates the dynamics

[TABLE]

where the response function $g_{e}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}$ is assumed to be an increasing differentiable function satisfying $g_{e}(1)=1$. For the shortest path problem in a network of parallel links, Bonifaci shows convergence to an optimal solution. 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: $g(y)=y^{\mu}$ for some $\mu>0$.

• •

Saturating response: $g(y)=(1+\alpha)y^{\mu}/(1+\alpha y^{\mu})$ for some $\mu,\alpha>0$.

Lemma 4**.**

$L(x,t)=p^{T}b+c^{T}x$* is a Lyapunov function for the dynamics (5). Moreover $L(x,t)=0$ if and only if $x$ is a fixed point of (5).*

Proof.

We proceed as in the proof of Lemma 3. Let $\lambda_{e}=|A_{e}^{T}p|/c_{e}$ and note that $|q_{e}|/x_{e}=\lambda_{e}\geq 0$. Then, we have

[TABLE]

where the inequality follows by $g_{e}(1)=1$ and $g_{e}$ is an increasing function implies that the terms $\lambda_{e}^{2}-1$ and $g_{e}(\lambda_{e})-1$ have the same sign.

Moreover, the derivative is zero if and only if for all $e$, either $x_{e}=0$ or $\lambda_{e}=1$ (as $\lambda_{e}\geq 0$). Since the latter condition is equivalent to $|A_{e}^{T}p|=c_{e}$, it follows for every $e$ with $x_{e}\neq 0$ that $|q_{e}|=x_{e}$ or equivalently $\dot{x}_{e}=0$. ∎

We remark that the proof above would even work for transfer-functions $g_{e}(x,t,y)$. It is only important that $g_{e}(x,t,1)=1$ and that the function is increasing in $y$.

It now follows from the general theory of dynamical systems that $x(t)$ converges to a fixed point.

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

As $t\rightarrow\infty$, $x(t)$ and $|q|(t)$ approach a fixed point $x_{0}$. Moreover, $E_{x}(q)$ and $\mathrm{cost}(x)$ converge to $c^{T}x_{0}$.

Proof.

Same proof as Corollary 3. ∎

We finally show convergence to the optimum solution of (1) under the additional assumption that $g_{e}(y)\geq 1+\alpha(y-1)$ for some $\alpha>0$ and all $y$ and $e$.

Theorem 5**.**

Assume further that $g_{e}(y)\geq 1+\alpha(y-1)$ for some $\alpha>0$ and all $e$ and $y$. As $t\rightarrow\infty$, $x(t)$ converges to the optimal solution $x^{*}$.

Proof.

Assume that $x(t)$ converges to a non-optimal fixed point $z$. Let $x^{*}$ be the optimal solution and let

[TABLE]

Let $\delta=(\mathrm{cost}(z)-\mathrm{cost}(x^{*}))/2$. Then $E_{x}(q(t))\geq\mathrm{cost}(z)-\delta=\mathrm{cost}(x^{*})+\delta$, for all sufficiently large $t$. Further, by definition $q_{e}=(x_{e}/c_{e})A_{e}^{T}p$ and thus

[TABLE]

where the first inequality follows from $f(y)\geq 1+\alpha(y-1)$ for all $y$ and the last inequality follows from (18). Hence $W\rightarrow\infty$, a contradiction to the fact that $x$ is bounded. ∎

We note that convex increasing functions satisfy $g(y)\geq g(1)+\alpha(y-1)$ with $\alpha=g^{\prime}(1)$.

5 Open Problems

For the dynamics (3), we showed convergence to the optimal solution under the assumptions:

1. (i)

$a_{e}(x,t)$ is bounded and bounded away from zero; 2. (ii)

$a_{e}(x,t)=a_{e}(t)$ does not depend on the state $x$; 3. (iii)

$\dot{a}_{e}\geq 0$ always.

We argued that assumption (i) is necessary. How about assumptions (ii) and (iii)?

For the uniform dynamics, convergence of a suitable Euler discretization was shown in [BBD*+*13, SV16a] for the shortest path problem and the basis pursuit problem respectively. What can be said about the convergence of the discretization of the non-uniform dynamics?

Our proof that Bonifaci’s refined model converges to the optimum solution requires the additional assumption that $g_{e}(y)\geq 1+\alpha(y-1)$ for some $\alpha>0$ and all $e$ and $y\geq 0$. Can this condition be relaxed?

There is also the directed dynamics $\dot{x}=q-x$ considered in [IJNT11, SV16b]. Can convergence be shown for its non-uniform version? This question is answered affirmatively in [FKKM19].