Path storage in the particle filter

Introduction

Consider the problem of filtering in state-space models (Cappé et al, 2005) defined by $X_{0}\sim\mu(\cdot)$ and for $t=1,\ldots,T$

Here $X_{0:T}$ is a hidden Markov chain in some space $\mathcal{X}$ with initial distribution $\mu$ and transition density $f$. The observations $Y_{1:T}$ in space $\mathcal{Y}$ are conditionally independent given $x_{1:T}$, with measurement density $g$. For any vector $v$, introduce the notation $v_{1:n}=(v_{1},\ldots,v_{n})$ and $v^{1:n}=(v^{1},\ldots,v^{n})$.

We denote by $p_{t}$ the distribution of the path $X_{0:t}$ given the observations $y_{1:t}$ available at time $t$, from which the filtering distribution of $X_{t}$ given $y_{1:t}$, denoted by $\pi_{t}$, is a marginal. The bootstrap particle filter (Gordon et al, 1993), described in Algorithm 1, recursively approximates the distributions $p_{1:T}$, and has borne various other sequential Monte Carlo methods (Doucet et al, 2001; Doucet and Johansen, 2011).

In Algorithm 1 the resampling step relies on some distribution $\mathcal{R}$ on $\{1,\ldots,N\}^{N}$ taking normalized weights as parameters.

At each time $t$, Algorithm 1 approximates $p_{t}$ and $\pi_{t}$ by the empirical distributions

It has been shown in Whiteley (2011); Douc et al (2012); van Handel (2009), and in Theorem 7.4.4 in Del Moral (2004) that $\pi_{t}^{N}$ converges to $\pi_{t}$ with $N$ under mild conditions on the model laws $(\mu,f,g)$, and that the Monte Carlo error is constant with respect to $t$. However it is also well-known that the path measures $p_{t}^{N}$, while converging to $p_{t}$ with $N$, have a Monte Carlo error typically exploding at least quadratically with the time $t$ (Del Moral and Doucet, 2003; Poyiadjis et al, 2011). Indeed the paths quickly coalesce due to the resampling steps, thus providing a poor approximation of the marginal distributions $p(dx_{s}|y_{1:t})$ for large values of $t-s$. In the following we refer to the collection of paths $\bar{x}_{0:t}^{1:N}$ as the ancestry tree, to each $x^{k}_{s}$ (for $k=1,\ldots,N$ and $s=0,\ldots,t$) as a node, to each $x^{k}_{t}$ more specifically as a leaf node, and to paths as branches.

Figure 1 might help to visualise the typical shape of the ancestry tree generated by a particle filter. The time at which all the branches coalesce, denoted by $c_{T}$, separates the “trunk” made of a unique branch from $t=0$ to $t=c_{T}-1$ from the “crown” made of all the branches from $t=c_{T}$ to $t=T$. Despite its negative consequence on the estimation of filtering quantities, the particle degeneracy phenomenon results in crowns of small sizes, allowing full trees to be stored at low memory cost. This can be beneficial whenever full paths of the particle filter are required, such as for the conditional sequential Monte Carlo and particle Gibbs algorithms first described in Andrieu et al (2010), studied in Chopin and Singh (2013), and used extensively in Chopin et al (2013) and Lindsten et al (2012). Another instance of sequential Monte Carlo method requiring path storage is presented in Wang et al (2014) in the context of computational biology. In the present article algorithms and results are presented in the filtering terminology, however they immediately extend to any sequential Monte Carlo method for Feynman-Kac models (Del Moral, 2004).

In Section 2 we present an efficient algorithm to store ancestry trees recursively during the run of a particle filter. In Section 3 we present new theoretical results bounding the size of ancestry trees, in order to bound the expected memory requirements of the storage algorithm. Finally the theoretical results and the algorithmic performance are tested numerically in Section 4.

Algorithms

This section introduces a memory-efficient data structure and associated algorithms for storing only those paths with support at time $t$. The algorithms are designed for parallel execution, in keeping with the general parallelisability of other components of sequential Monte Carlo samplers (Lee et al, 2010; Murray, 2013).

Up to time $t$, the particle filter produces particles $x^{1:N}_{1:t}$ and ancestors $a^{1:N}_{1:t}$. From $a^{1:N}_{1:t}$, offspring counts $o^{1:N}_{1:t}$ are readily obtained (Murray et al, 2013), where $o^{i}_{t}$ represents the number of children at generation $t$ of particle $x^{i}_{t-1}$. Let $x^{1:M}_{*}$ represent $M$ slots in memory for storing particles. At any time, some of these slots are empty, while others store the nodes of the tree. Let $a^{1:M}_{*}$ be an ancestry vector, where $a^{i}_{*}=0$ if $x^{i}_{*}$ is empty or a root node, and otherwise $a^{i}_{*}=j$ to indicate that the particle in $x^{j}_{*}$ is the parent of the particle in $x^{i}_{*}$. Let $o^{1:M}_{*}$ be the offspring vector corresponding to $a^{1:M}_{*}$, where $o^{i}_{*}=n$ indicates that $x^{i}_{*}$ has $n$ children. Finally, let $l^{1:N}_{*}$ give the numbers of the $N$ slots in $x^{1:M}_{*}$ that store the particles of the youngest generation; these are the leaf nodes of the tree.

Basic operations on the tree are its initialisation, the insertion of a new generation of particles, and the pruning of older particles to remove those without a surviving descendent in the youngest generation. These operations are described in Algorithm 2. The descriptions there rely on primitive operations defined in Algorithm 3. The efficient implementation of such primitives is well understood in both serial and parallel contexts, so that they make useful building blocks for the higher-level algorithms.

To begin, the first of the $M$ empty slots of the tree are initialised with the first generation of $N$ particles as in the init procedure of Algorithm 2. We assume, for now, that $M$ is sufficiently large to accommodate all subsequent operations on the tree, but see remarks in Section 2.2 below.

Each new generation is inserted as in the insert procedure of Algorithm 2. The procedure searches for nodes with no offspring in the current generation, and replaces them with the new leaf nodes. The vector $z_{*}$ is introduced, where $z_{*}^{i}$ is equal to the number of nodes between $1$ and $i$ with no offspring. Nodes to replace are then located by searching for the increments in $z_{*}$. The new generation is inserted at these locations.

Finally, the tree is pruned before the insertion of each new generation $t$, using the prune procedure of Algorithm 2. This requires the offspring vector, $o_{t}$, of the new generation. The algorithm determines which of the current leaf nodes have no offspring in the new generation, decrements the offspring counts of their parent nodes, and proceeds recursively up the tree for cases where the parent has no remaining offspring either. Each non-leaf node $i$ is considered pruned if $o^{i}_{*}=0$, and may be overwritten by future calls to insert.

2 Remarks and improvements

The insert procedure of Algorithm 2 assumes that there are at least $N$ free slots in which to place the latest nodes. If this is not true, the buffer can be enlarged by allocating a larger block of memory, copying the contents of the ancestry tree across, and filling the new regions of the $o_{*}$ and $a_{*}$ vectors with zeros. Various heuristics can be used to set the new size $M$, aiming to reduce fragmentation and the chance of future increases. Because memory reallocations involve an expensive copy, it is worth increasing $M$ more than strictly necessary to postpone additional reallocations. For instance, implementations of the C++ Standard Template Library typically double the storage capacity of a vector that is extended by just one element, anticipating further extensions. A more conservative strategy is to start with a value of $M$ equal to a small multiple of $N$, and enlarge by $N$ slots whenever necessary. Ultimately, we have not found that the particular enlargement strategy affects execution time a great deal, particularly since, as in the proceeding theoretical results, the size of the ancestry tree crown is independent of $T$, so that the need for reallocations diminishes as $t$ increases.

According to the results of Section 3, the expectation of the size of the tree grows linearly with $T$, but this is only due to the trunk. The size of the crown is independent of $T$. It may be possible to improve the algorithms by identifying the nodes along the trunk and storing them separately, as these nodes will never be overwritten by subsequent insertions. Under this modified scheme a separate, single growing trunk needs to be stored but not searched, while the nodes of the crown need to be stored and searched at every time step. The number of nodes in the crown is of constant expectation according to Theorem 3.1 of Section 3. Hence this modification induces a scheme of constant expected computational cost in $T$, which could be relevant in applications where the time horizon is very long, although there will be overhead in identifying the trunk. See Fig. 3(b) in Section 4 for a report on the computational cost of the proposed method. Memory reallocation is also reduced by storing the trunk separately.

We establish in Section 3 that the size of the tree is expected to be bounded by $T+\Delta_{2}N\log N$ for some constant $\Delta_{2}$. The size of the data structure, $M$, must be at least as large as this. We assume that, with a sensible enlargement strategy, it is no more than a constant factor larger than this, so that its expected memory complexity is $\mathcal{O}(T+\Delta_{2}N\log N)$.

The computational complexity of init is linear in the size of the data structure, $\mathcal{O}(T+\Delta_{2}N\log N)$. A serial implementation of insert permits a linear prefix sum and search, so that insert is also $\mathcal{O}(T+\Delta_{2}N\log N)$. In parallel, a linear prefix sum is still achieved (Sengupta et al, 2008), but the search becomes $N$ binary searches, logarithmic to the size of the data structure; overall $\mathcal{O}(N\log(T+\Delta_{2}N\log N))$.

For prune, consider the best case, where all particles of the previous generation have an offspring in the new generation. The complexity is then $\mathcal{O}(N)$: the algorithm operates on each of the $N$ new nodes, but does not traverse the tree further. Now consider the worst case, where only one particle of generation $t$ has offspring in the new generation $t+1$. In this case all but $t$ nodes of the existing tree are pruned, so that the complexity is $\mathcal{O}(T+\Delta_{2}N\log N-t)$ – linear in the size of the data structure, and parellelisable.

Finally, the transform-prefix-sum across the full vector $o_{*}$ in the insert is redundant. The sum can be truncated once it has reached $N$, as a sufficient number of free slots have then been found. This is simple to achieve in the serial case, but it is not obvious how to achieve it in the parallel case. Heuristic include considering only a subset of $o_{*}$ at a time and iterating until a sufficient number of free slots are found, and starting the cumulative sum after the last slot that was filled in the previous call to insert. In practice, however, we have observed only negligible variation in execution times when applying such heuristics, and so have chosen to present the simplest version here.

Size of the ancestry tree

From a theoretical point of view, similar random trees have been studied in population genetics Del Moral et al (2009); Möhle (2004) in a setting that corresponds to a state-space model that assigns equal weights to all paths; these results do not apply directly here. In order to bound the expected number of nodes in an ancestry tree, we first study the distance $d_{T}=T-c_{T}$ between the final time $T$ and the full coalescence time $c_{T}$ when all the paths merge. Theorem 3.1 proposes a bound on the expectation of $d_{T}$, which is independent of $T$ and explicit in $N$.

There exists $\epsilon\in$ such that for all $y\in\mathcal{Y}$ and for all $x\in\mathcal{X}$

Under Assumption 1 the distance to the most recent common ancestor $d_{T}$ satisfies

for some $\Delta_{1}>0$, which does not depend $N$ nor $T$.

The expected number of nodes in the tree can be bounded explicitly in $N$ and $T$, as in Theorem 3.2.

We suppose here that $N\geq 3$. Under Assumption 1 the number of nodes, denoted by $n_{T}$ at time $T$, satisfies

for some $\Delta_{2}>0$ that does not depend on $N$ nor $T$.

These results quantify the practical difference between storing all the generated particles (for a deterministic cost of $T\times N$ memory units) and storing only the surviving particles (for a random cost expected to be bounded by $T+\Delta_{2}N\log N$).

Assumption 1 is very strong outside compact spaces, and for instance does not even cover the linear-Gaussian case, although the experiments of Section 4 indicate that similar results might hold for non-linear and non-Gaussian cases. The numerical experiments show that the bound is accurate as a function of $N$, so that even if some inequalities used in the proofs appear quite crude, the overall result is precise. However the results do not capture the shape of the tree as a function of $\epsilon$, which is why we write the constants $\Delta_{1}$ and $\Delta_{2}$ without making their dependency on $\epsilon$ explicit. Consider for example Theorem 3.1, where $\Delta_{1}$ can be defined by $\Delta_{1}=1+8/\epsilon$, as will be proven in Section 3.3. If the bound was sharp as a function of $\epsilon$, it would mean that the time to full coalescence increases to infinity when $\epsilon$ goes to zero. However path degeneracy is expected be more acute for smaller $\epsilon$, since more variability in the particle weights is then allowed. The dependency on $\epsilon$ in $\Delta_{1}$ is thus not realistic. We believe the bounds could in fact be independent of $\epsilon$, by considering $\epsilon=1$ as the case corresponding to the largest expectations of $d_{T}$ and $n_{T}$; a claim not proven here.

Moreover, the proposed proof relies on the multinomial resampling scheme, while most practitioners favour more sophisticated schemes (Carpenter et al, 1999; Liu and Chen, 1998; Kitagawa, 1998; Doucet and Johansen, 2011). Figure 3(a) of Section 4 indicates that similar results hold for these other resampling schemes. There are some obvious counter-examples, for instance when the measurement density is constant, leading to equal weights at each step (equivalently $\epsilon=1$). Then the results above hold for multinomial resampling but systematic resampling would completely obviate the path degeneracy phenomenon. Describing features of ancestry trees corresponding to general resampling schemes would constitute an interesting avenue of research.

The rest of the section is devoted to proving Theorem 3.1 and Theorem 3.2.

2 From non-uniform weights to uniform weights

We first relate the ancestry process associated with particle filters using multinomial resampling, with the ancestry process associated with the neutral case, where all the weights would be equal to $N^{-1}$ at every time step. To do so we introduce various intermediate processes, starting with the exact multinomial resampling process denoted by $(A_{t})_{t\geq 0}$, then an approximation represented by $(A^{\prime}_{t})_{t\geq 0}$ which provides an almost sure upper bound and eventually a process $(Z_{k})_{k\geq 0}$ counting the number of nodes at generation $T-k$ in the neutral case, for a fixed time horizon $T$.

For each time $t$, define $A_{t}:\,j\in\{1,\ldots,N\}\mapsto a_{t}^{j}\in\{1,\ldots,N\}$ and then $A_{t}^{\prime}:\{1,\ldots,N\}\rightarrow\{1,\ldots,N\}$ as follows. For all $j$ in $C_{t}=\{k\in\{1,\ldots,N\}:V_{t}^{k}\leq\epsilon\}$, set $A_{t}^{\prime}(j)=a_{t}^{j}$. Order the $p$ remaining indices of the set $\{j\in\{1,\ldots,N\}:V_{j}^{t}>\epsilon\}$ into $\{j_{1}<\dots<j_{p}\}$, set $A^{\prime}_{t}(j_{1})=\inf(\{1,\ldots,N\}\backslash A^{\prime}_{t}(C_{t}))$ and then recursively

Such a function $A_{t}^{\prime}$ almost surely maps to more unique values than $A_{t}$ by construction. It can be seen as a mixture of two steps, as described for $A_{t}$ above, but this time neither step relies on the values of the weights.

We write $\lvert u\rvert$ for the cardinal of the image of a function $u:\{1,\ldots,N\}\rightarrow\{1,\ldots,N\}$. In terms of the functions $(A_{k})_{k\leq T-1}$, the full coalescence time $c_{T}$ can be defined as

with the convention $c_{T}=0$ in the event $\mid A_{k}\circ A_{k+1}\circ\dots\circ A_{T-1}\mid>1$ for each $0\leq k\leq T-1$, which almost surely satisfies $c_{T}\geq c^{\prime}_{T}$ with

Indeed since $A^{\prime}_{t}$ maps to more unique values than $A_{t}$ at each time $t$, the quantity $\mid A^{\prime}_{k}\circ\dots\circ A^{\prime}_{T-1}\mid$, counting the unique ancestors from generation $k$ of the particles at time $T$ when using the resampling scheme $A^{\prime}$, is almost surely larger than $\mid A_{k}\circ\dots\circ A_{T-1}\mid$ for any $k$, and hence it takes longer to reach the full coalescence time when using $A^{\prime}$ compared to $A$.

Following Del Moral et al (2009), Section 4 and Möhle (2004), the sequence $(K_{k})_{k\geq 0}=(\mid A^{\prime}_{T-k}\circ\dots\circ A^{\prime}_{T-1}\mid)_{k\geq 0}$ is a Markov chain in the filtration $(\mathcal{F}_{k})_{k\geq 1}$ with

with the convention $K_{0}=N$. For all $k\geq 0$, $q\in\{1,\ldots,N\}$ and $p<q$ its transition law verifies

where $\genfrac{\{}{\}}{0.0pt}{}{q}{p}$ is the Stirling number of the second kind giving the number of ways of partitioning the set $\{1,\ldots,q\}$ into $p$ non empty blocks and where $(N)_{p}=N!/(N-p)!$. Note that Eq. (2) is a special case of Eq. (1).

following the same reasoning as for the transition probabilities of $(K_{k})_{k\geq 0}$. The initial distribution of $Z_{0}$ is not used in the following hence we do not need to specify it. The link between $\left(Z_{k}\right)_{k\geq 0}$ and $\left(K_{k}\right)_{k\geq 0}$ is explicitly given by

Note that the process $(Z_{k})_{k\geq 0}$ is not used in the proof of Theorem 3.1, where we start from $(K_{k})_{k\geq 0}$ again, but is pivotal for the proof of Theorem 3.2.

3 Distance to the most recent common ancestor

where $p_{N,q}$ is defined in Eq. (2). In addition we couple $(L_{k})_{k\geq 0}$ and $(K_{k})_{k\geq 0}$ by assuming

$\left[L_{k}=K_{k}\right]<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mo>⇒</mo></mrow><annotation encoding="application/x-tex">\Rightarrow</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.3669em;"></span><span class="mrel">⇒</span></span></span></span></span>\left[L_{k+1}<L_{k}\Leftrightarrow K_{k+1}<K_{k}\right]$ (if the two chains are at the same point, then if one of them decreases, the other one decreases too)

$[L_{k}\neq K_{k}]\Rightarrow$$K_{k+1}$and $L_{k+1}$ are independent, conditionally upon $L_{k}$, $K_{k}$.

By construction $L_{k}\geq K_{k}$ for all $k\geq 0$ almost surely. Hence $c^{\prime}_{T}\geq T-\mathcal{D}_{T}$ with $\mathcal{D}_{T}=\inf\{k\geq 1:L_{k}=1\}$ and thus $d_{T}=T-c_{T}\leq T-c^{\prime}_{T}\leq\mathcal{D}_{T}$ almost surely.

We have, for all $N$, $(8N)^{-1}\leq 1-\exp\{-1/2N\}$ and for all $x\geq 1$ and $\varepsilon\in(0,1)$, $(1-\varepsilon/x)^{x}\leq\exp(-\varepsilon)$; combining these inequalities we obtain

where $\alpha=\exp(-\epsilon/8)$. We can now bound this series by expanding $\alpha^{q/N}=\exp\{(q/N)\log\alpha\}$ into an alternating series and by bounding the alternating series always by one of its partial sums:

which concludes the proof of Theorem 3.1.

4 Number of nodes in the ancestry tree

We now proceed to the proof of Theorem 3.2. Denote by $m_{T}$ the number of nodes in the crown. The bound on $d_{T}$ from Theorem 3.1 gives a first crude bound

which is obtained by bounding the size of every generation in the crown by $N$. However we can obtain a better bound, in $N\log N$, by the following arguments.

By expanding $(1-(1-\epsilon/N)^{q})$ into its alternating series and bounding the series by its third partial sum, we obtain

Now for $x\in[1,N]$ define the function $g_{N,\epsilon}$ by:

Noting that $g_{N,\epsilon}$ is concave and using Jensen’s inequality, we obtain

Then there exists $C>0$ independent of $N$ such that

The proof of Lemma 1, based on elementary real analysis, is given in Appendix A. Using Lemma 1 and Eq. (8) we obtain Theorem 3.2 with $\Delta_{2}=C+\Delta_{1}$.

Numerical experiments

This section provides numerical experiments to illustrate the results of Section 3 and the efficiency of the algorithms presented in Section 2. The results summarise $K=500$ independent runs, using $N=128$ particles and $T\leq 1000$ time steps. For each run, a new synthetic dataset is generated and a different random seed is used. The default resampling scheme is the multinomial scheme, applied at every time step. The algorithms of Section 2 have been implemented in LibBi (Murray, 2013, www.libbi.org), which is used for the numerical results here.

We use the Phytoplankton-Zooplankton (PZ) model described in Jones et al (2010) and Murray et al (2012). Concentrations of phytoplankton $(P_{t})$ and zooplankton $(Z_{t})$, along with the stochastic growth rate of phytoplankton $(\alpha_{t})$, constitute the hidden state. The state follows the continuous-time dynamics $dP/dt=\alpha_{t}P-cPZ$ and $dZ/dt=ecPZ-m_{l}Z-m_{q}Z^{2}$, with $\alpha_{t}\sim\mathcal{N}(\mu,\sigma^{2})$ drawn at every integer time $t$. The initial conditions are $\log P_{0}\sim\mathcal{N}(\log(2),0.2)$, $\log Z_{0}\sim\mathcal{N}(\log(2),0.1)$. The observations $(Y_{t})$ measure $(P_{t})$ with additive log-normal noise, that is $\log Y_{t}\sim\mathcal{N}(\log P_{t},\sigma_{y})$. The parameters are set to $\mu=0.4$, $\sigma=0.2$, $c=0.25$, $e=0.3$, $m_{l}=m_{q}=0.1$ and $\sigma_{y}=0.2$.

To illustrate the efficiency of the procedures presented in Section 2, Fig. 3(b) shows the combined time taken to execute the pruning and insertion algorithms at each time step, for various $T$ and $N$. The results suggest that the computational cost is not greatly influenced by $T$, and close to linear with respect to $N$: evidence of a practical implementation with comparable complexity to the particle filter itself.

Conclusion

We have presented a bound on the expected number of nodes in the ancestry tree produced by particle filters. The numerical experiments of Section 4 indicate that the result is accurate, even outside the scope of the assumptions made in the theoretical study, and that the proposed algorithm to store the tree is computationally efficient.

Appendix A Proof of Lemma 1

however this contraction coefficient depends on $N$ and a direct use of it yields a bound on $\sum_{k\geq 0}(u_{k}-1)$ that is not in $N\log N$.

Note also that even though $u_{k}$ goes to $1$, we can focus on the partial sum $\sum_{k=0}^{\sigma_{2}}(u_{k}-1)$ where $\sigma_{2}=\inf\{k:u_{k}\leq 2\}$, because $\sum_{k=\sigma_{2}}^{\infty}(u_{k}-1)$ is essentially bounded by $N$. Indeed note that for $1\leq u\leq 2$ we have $(\epsilon^{3}/6N^{2})u(u-1)(u-2)\leq 0$ so that

hence $\sum_{k=\sigma_{2}}^{\infty}(u_{k}-1)\leq(2N/\epsilon^{2})$. Therefore we can focus on bounding $\sum_{k=0}^{\sigma_{2}}(u_{k}-1)$ by $N\log N$. Let us split this sum into partial sums, where the first partial sum is over indices $k$ such that $N/2\leq u_{k}\leq N$, the second is over indices $k$ such that $N/4\leq u_{k}\leq N/2$, etc. More formally, we introduce $(k_{j})_{j=0}^{J}$ such that $k_{0}=0$, $k_{1}=\inf\{k:u_{k}\leq N/2\}$, …, $k_{j}=\inf\{k:u_{k}\leq N/2^{j}\}$, up to $k_{J}=\inf\{k:u_{k}\leq N/2^{J}\}$ where $J$ is such that $N/2^{J}\leq 2$, or equivalently $\log N/\log 2-1\leq J$. For instance we take $J=\lfloor\log N/\log 2\rfloor$. Thus we have split $\sum_{k=0}^{\sigma_{2}}(u_{k}-1)$ into $J$ partial sums of the form $\sum_{k=k_{j}}^{k_{j+1}-1}(u_{k}-1)$ and we are now going to bound each of these partial sum by the same quantity $C(\epsilon)N$ for some $C(\epsilon)$ that depends only on $\epsilon$.

To do so, we consider the time needed by $(u_{k})_{k\geq 0}$ to decrease from a value $N/m_{j}$ to a value $N/m_{j+1}$, with $m_{j+1}>m_{j}$; we will later take $m_{j}=2^{j}$ and $m_{j+1}=2^{j+1}$. Note that for any $m$ we have

and note that for any $N\geq 6$ and $m\leq N/2$ we have

which is clear upon noticing that $\beta(N,m,\epsilon)$ as a function of $m$ on $[1,N/2]$ is concave and thus reaches its minimum in $1$ or $N/2$ (and this minimum is greater than $\epsilon^{2}/4$, provided $N\geq 6$). For any $x\geq N/m_{j+1}$ we can check that

by noticing that $g_{N,\epsilon}$ is concave and that $g_{N,\epsilon}(x)\leq x$ for $x\in[0,N]$. Hence for $k\geq 0$ such that $u_{k-1}\geq N/m_{j+1}$, we have

Now suppose that for some $k_{j}\geq 0$ we have $u_{k_{j}}\leq N/m_{j}$. Then let us find $K$ such that $u_{k_{j}+K}\leq N/m_{j+1}$. It is sufficient to find $K$ such that

guarantees the inequality $u_{k_{j}+K}\leq N/m_{j+1}$. In other words $(u_{k})_{k\geq 0}$ needs less than $K$ steps to decrease from $N/m_{j}$ to $N/m_{j+1}$. Summing the terms between $k_{j}$ and $k_{j}+K$, we obtain

Taking $m_{j}=2^{j}$ and $m_{j+1}=2^{j+1}$, we have $k_{j+1}\leq k_{j}+K$ and thus obtain

with $C(\epsilon)$ independent of $N$. We have thus bounded the full sum by

for some $D(\epsilon)$ independent of $N$.

References