Is Learning in Biological Neural Networks based on Stochastic Gradient
  Descent? An analysis using stochastic processes
In recent years, there has been an intense debate about how learning in
biological neural networks (BNNs) differs from learning in artificial neural
networks. It is often argued that the updating of connections in the brain
relies only on local information, and therefore a stochastic gradient-descent
type optimization method cannot be used. In this paper, we study a stochastic
model for supervised learning in BNNs. We show that a (continuous) gradient
step occurs approximately when each learning opportunity is processed by many
local updates. This result suggests that stochastic gradient descent may indeed
play a role in optimizing BNNs.
Sören Christensen, Jan Kallsen
# Is Learning in Biological Neural Networks based on Stochastic Gradient Descent?

###### Abstract

In recent years, there has been an intense debate about how learning in biological neural networks (BNNs) differs from learning in artificial neural networks. It is often argued that the updating of connections in the brain relies only on local information, and therefore a stochastic gradient-descent type optimization method cannot be used. In this paper, we study a stochastic model for supervised learning in BNNs. We show that a (continuous) gradient step occurs approximately when each learning opportunity is processed by many local updates. This result suggests that stochastic gradient descent may indeed play a role in optimizing BNNs.

**Keywords:** Biological neural networks, Schmidt-Hieber model, stochastic gradient descent, supervised learning

## 1 Introduction

In order to understand how biological neural networks (BNNs) work, it seems natural to compare them with artificial neural networks (ANNs). Although the definition of the latter is inspired by the former, they also differ in several aspects. One of them is the way the network parameters are updated.

In simple terms, an ANN learns from data by adjusting the weights of the connections between nodes in order to minimize a loss function that measures the difference between the desired output and the actual output of the network. More specifically, the optimization step is performed using the Stochastic Gradient Descent (SGD) algorithm, which iteratively updates the weights of the network by moving them in the direction of the steepest descent of the empirical loss function of a single training sample. The gradient itself is computed with the so-called backpropagation algorithm. In particular, the update of any parameter is based on the states of all other parameters. Such a mechanism does not seem to be biologically plausible for BNNs, as many authors have pointed out. Parameter update in BNNs occurs only locally, and distant neurons are only indirectly connected through the endogenous reward system. This observation is closely related to the weight transportation problem [6; 2; 4]. We refer to [12; 11] for a detailed discussion about the role of SGD in BNN, which the author of [10; Section 5] summarizes as follows: "[T]here are various theories that are centered around the idea that the learning in BNNs should be linked to gradient descent. All of these approaches, however, contain still biological implausibilities and lack a theoretical analysis."

The starting point for the present paper is the recent article [10] just cited. In this seminal study, the author proposes a very persuasive stochastic model for brain-supervised learning
which has a thorough biological foundation in terms of spike-timing-dependent plasticity. We review and discuss this setup in Section 2. In this model the local updating rule of the connection parameters in BNNs turns out to be a zero-order optimization procedure. More precisely, it is shown in [10] that the expected value of the iterates coincides with a modified gradient descent. However, this holds only on average. The noise for such zero-order methods is so high that one can hardly imagine effective learning based on it, see [3, 8, 1]. The author himself writes in [10, Section 4]: "It remains to reconcile the observed efficiency of learning in biological neural networks with the slow convergence of zero-order methods."

In this paper we make an attempt to achieve this reconciliation. To this end, we consider in Section 3 a slight modification of the model of [10]. More specifically, we relax the assumption that for each learning opportunity and each connection, exactly one spike is released. Instead, we assume that a large number of spikes is released for each training sample in the BNN model and thus many parameter updates are made. It turns out that with this modification, the updates correspond approximately to a continuous descent step along the gradient flow, see Theorem 1. This can be interpreted in the sense that it is not biologically implausible that BNNs use a kind of SGD algorithm after all, but without explicitly computing the gradient.

## 2 The Schmidt-Hieber model for BNNs revisited

We begin this section by reviewing the model introduced in [10]. It considers a classical instance of supervised learning: input-output pairs \((\mathbf{X}_{1},Y_{1}),(\mathbf{X}_{2},Y_{2}),\ldots\) are given as observations, all being identically distributed. The goal is to predict the output \(Y\) for each new input \(\mathbf{X}\) based on previous training data. This setting includes, for example, classification (when the set of possible outcomes of \(Y\) is finite) or regression problems.

A (feedforward) biological neural network (BNN) is modeled in [10] as a directed acyclic graph with input neurons receiving information from the observations \(\mathbf{X}_{k}\) and generating a predicted response \(\widehat{Y}_{k}\) as output. The nodes represent the neurons in the network and an edge \(\nu=(i,j)\) between two nodes \(i\) and \(j\) indicates that neuron \(i\) is presynaptic for neuron \(j\). Each element \((i,j)\) in the edge set \(\mathcal{T}\) has a weight \(w_{ij}\) which indicates the strength of the connection between \(i\) and \(j\). While the structure of the graph does not change, the weights \(w_{ij}\) are adjusted in each learning step.

Spike-timing-dependent plasticity (STDP) is chosen as the biological mechanism to update the parameters. It is considered as a form of Hebbian learning [5], which states that neurons that fire together wire together. More precisely, the synaptic weight \(w_{ij}\) changes depending on the timing of the spikes from neuron \(i\) to neuron \(j\). The weight decreases when neuron \(i\) spikes before neuron \(j\), and increases when neuron \(j\) spikes before neuron \(i\). The closer the spikes are in time, the larger the change in weight. It is important to note that the spike times are modeled as random variables. After some standardization (see [10, Equation (4.3)]) the update of the parameter for edge \((i,j)\) becomes

\[w_{ij}\gets w_{ij}+w_{ij}C(e^{-U_{ij}}-e^{U_{ij}}),\]

where \(U_{ij}\) are uniformly distributed random variables on some interval \([-A,A]\), i.e. \(U_{ij}\sim\mathcal{U}(-A,A)\), modeling the random spike times. The constant \(C\) represents the effect of a reward system, for example by neurotransmitters such as dopamine. It plays a key role for any meaningful learning process. In the present setup, the reward is tied to the success of predicting \(Y_{k}\), more specifically, whether the task is solved better or worse than in earlier trials. Using the standardization \(\theta_{ij}:=\log w_{ij}\) and a Taylor approximation, the following structure is derived in [10, Equation (4.5)] for the update of \(\theta_{ij}\) at step \(\ell\):

\[\theta_{ij}^{(\ell)}=\theta_{ij}^{(\ell-1)}+\alpha^{(\ell-1)}\Big{(}L^{(\ell-1 )}(\boldsymbol{\theta}^{(\ell-1)}+\mathbf{U}^{(\ell)})-L^{(\ell-2)}( \boldsymbol{\theta}^{(\ell-2)}+\mathbf{U}^{(\ell-1)})\Big{)}\Big{(}e^{-U_{ij}^ {(\ell)}}-e^{U_{ij}^{(\ell)}}\Big{)}, \tag{2.1}\]
where \(\alpha^{(\ell)}>0\) is a learning rate, \(\mathbf{\theta}^{(\ell)}=\left(\theta^{(\ell)}_{ij}\right)_{(i,j)\in T}\) denotes the vector of parameters for all edges, \(\mathbf{U}^{(\ell)}\) the vector of the independent uniformly distributed random variables \(U^{(\ell)}_{ij}\) for all edges \((i,j)\), and \(L^{(\ell)}\) the loss function associated with the respective step \(\ell\). Thus, the update of the weights of the individual edges is in fact affected by the state of the entire network, but only through the value of the common loss function, which provides an assessment of the learning success. In particular, no gradient appears in the mechanism.

In [10] the author considers the case where for each input-output pair \((\mathbf{X}_{k},Y_{k})\) the parameters are updated only once. To this end, the loss of the current input-output pair is compared with the previous one. More specifically, we have \(\ell=k\) as well as \(L^{(k)}(\mathbf{\theta})=L(\mathbf{\theta},\mathbf{X}_{k},Y_{k})\), \(k=1,2,\dots\), leading to the update rule

\[\theta^{(k)}_{ij}=\theta^{(k-1)}_{ij}+\alpha^{(k-1)}\Big{(}L( \mathbf{\theta}^{(k-1)}+\mathbf{U}^{(k)},\mathbf{X}_{k-1},Y_{k-1})-L(\mathbf{\theta}^ {(k-2)}+\mathbf{U}^{(k-1)},\mathbf{X}_{k-2},Y_{k-2})\Big{)}\Big{(}e^{-U^{(k)}_ {ij}}-e^{U^{(k)}_{ij}}\Big{)}. \tag{2.2}\]

As a main result, the author shows in [10, Theorem 1] that this procedure corresponds on average to a gradient descent method, with a gradient evaluated not exactly at \(\mathbf{\theta}^{(k-1)}\) but slightly perturbed randomly. However, as noted in [10], sufficiently fast convergence cannot be expected for such a zero-order method.

## 3 Multiple updates per learning opportunity

The key ingredient leading to (2.2) is that, for each learning opportunity and each connection, exactly one spike is triggered and thus only one update of the parameters is made. The author himself calls this assumption "strong", see [10, Section 4].

Given the average spike frequency in real biological systems and the strong brain activity even at immobile rest [7], it seems more reasonable to assume instead a large number of spikes per learning opportunity. This corresponds to a series \(\mathbf{\theta}^{(k,0)},\dots,\mathbf{\theta}^{(k,n)}\) of updates to the parameters after observing any input-output pair \((\mathbf{X}_{k-1},Y_{k-1})\). The assessment of the update steps is based on the loss function associated with the most recent observation \((\mathbf{X}_{k-1},Y_{k-1})\). More specifically, equation (2.1) turns into

\[\theta^{(k,\ell)}_{ij} =\theta^{(k,\ell-1)}_{ij}+\alpha^{(k-1,\ell-1)}\Big{(}L(\mathbf{ \theta}^{(k,\ell-1)}+\mathbf{U}^{(k,\ell)},\mathbf{X}_{k-1},Y_{k-1})-L(\mathbf{ \theta}^{(k,\ell-2)}+\mathbf{U}^{(k,\ell-1)},\mathbf{X}_{k-1},Y_{k-1})\Big{)}\] \[\quad\times\left(e^{-U^{(k,\ell)}_{ij}}-e^{U^{(k,\ell)}_{ij}}\right) \tag{3.1}\]

for \(\ell\geq 1\) with initial values given by \(\mathbf{\theta}^{(k,0)}:=\mathbf{\theta}^{(k,-1)}:=\mathbf{\theta}^{(k-1,n)}\). Since \(n\) is considered to be large, the individual update steps should be small in order to avoid overfitting.

We start by analyzing this update rule for a fixed observation \((\mathbf{X}_{k-1},Y_{k-1})\), i.e. for a fixed \(k\). For ease of notation we suppress the dependence on \(k\) in the following considerations. Since our goal is to study the limiting behavior for a large number of update steps \(n\), we instead make the dependence on \(n\) explicit. So we rewrite (3.1) as

\[{}^{n}\theta^{(\ell)}_{ij}={}^{n}\theta^{(\ell-1)}_{ij}+{}^{n}\alpha^{(\ell-1) }\Big{(}L^{(n}\mathbf{\theta}^{(\ell-1)}+{}^{n}\mathbf{U}^{(\ell)})-L^{(n}\mathbf{ \theta}^{(\ell-2)}+{}^{n}\mathbf{U}^{(\ell-1)})\Big{)}\Big{(}e^{-{}^{n}U^{( \ell)}_{ij}}-e^{{}^{n}U^{(\ell)}_{ij}}\Big{)}. \tag{3.2}\]

The parameter update depends on increments of a loss function which resembles a gradient on first glance. Note, however, that this increment term is the same for all edges and, moreover, randomness occurs in both the loss function and in the external factor. We now consider the following dependencies on \(n\):

\[{}^{n}\alpha^{(\ell-1)}:=\alpha,\quad A_{n}:=n^{-1/3}A,\quad{}^{n}U^{(\ell)}_ {ij}\sim\mathcal{U}(-A_{n},A_{n}), \tag{3.3}\]
with constants \(\alpha,A>0\). Moreover, we rescale and extend the discrete-time process \((^{n}\theta^{(\ell-1)}_{ij})_{\ell=-1,\ldots,n}\) in time, defining a continuous-time process \(\mathbf{Z}^{n}=(\mathbf{Z}^{n}_{t})_{t\in[-1/n,1]}\) by

\[\mathbf{Z}^{n}_{t_{\ell}^{n}}=\mathbf{Z}^{n}_{t_{\ell-1}}+\alpha_{n}\Big{(}L( \mathbf{Z}^{n}_{t_{\ell}^{n}}+^{n}\mathbf{U}^{(\ell)})-L(\mathbf{Z}^{n}_{t_{ \ell-2}^{n}}+^{n}\mathbf{U}^{(\ell-1)})\Big{)}\Big{(}e^{-u^{(\ell)}_{ij}}-e^{u ^{(\ell)}_{ij}}\Big{)} \tag{3.4}\]

for \(t_{\ell}^{n}:=\ell/n\) and by \(\mathbf{Z}^{n}_{t}:=\mathbf{Z}^{n}_{\lfloor tn/n\rfloor}\) for arbitrary \(t\in[0,1]\).

As a candidate limit for large \(n\) we consider a standard rescaled gradient process \(\mathbf{Z}=(\mathbf{Z}_{t})_{t\in[0,1]}\), which is defined as the solution to the deterministic ordinary differential equation (ODE)

\[\frac{d\mathbf{Z}_{t}}{dt}=-4\alpha\nabla L(\mathbf{Z}_{t}),\quad\mathbf{Z}_ {0}=\mathbf{\theta}^{(k,0)}, \tag{3.5}\]

see e.g. [9]. In order for our main theorem to hold, we assume that

\[\nabla L\]

 is bounded and Lipschitz continuous with Lipschitz constant

\[\lambda\]

. (3.6)

\(\mathbf{Z}\) naturally emerges as a limit if you run the ordinary gradient descent algorithm for minimizing the function \(L\) with many small steps. The main result of this paper states that the rescaled STDP process \(\mathbf{Z}^{n}\) converges to \(\mathbf{Z}\) as well. More precisely, we have

**Theorem 1**.: _Assume (3.3) and (3.6). Then, for each fixed training sample \(k\), the rescaled process \(\mathbf{Z}^{n}\) of the BNN weights converges to the rescaled gradient process \(\mathbf{Z}\) uniformly in \(L^{2}\), i.e._

\[\lim_{n\to\infty}\mathbb{E}\left(\sup_{t\in[0,1]}\|Z^{n}_{t}-Z_{t}\|^{2} \right)=0.\]

The previous theorem shows that learning in BNNs based on the local principle of STDP may indeed lead to optimization of parameters according to SGD if the number of spikes per learning opportunity is high. To wit, we start with an initial parameter vector \(\mathbf{\theta}^{(0)}\) and loss function \(L=L(\cdot,\mathbf{X}_{1},Y_{1})\). According to Theorem 1 the STDP nearly performs a continuous gradient step as in (3.5), leading to an updated parameter vector \(\mathbf{\theta}^{(1)}\). Switching now to the loss function \(L=L(\cdot,\mathbf{X}_{2},Y_{2})\), the next approximate gradient step leads to an updated vector \(\mathbf{\theta}^{(2)}\) etc.

This procedure only differs from the classical SGD in that we make many small instead of one large gradient step per learning opportunity. Interestingly, neither the gradient nor even the functional dependence of the loss function on the parameters need to be known explicitly for this purpose. By contrast, it relies crucially on the randomness in the update, which may seem counterintuitive because the desired gradient ODE (3.5) is deterministic.

Proof of Theorem 1.: Following the argument in [10, Proof of Theorem 1], we may decompose the dynamics of \(\mathbf{Z}^{n}\) in coordinate \(\nu=(i,j)\) as

\[Z^{n}_{t_{\ell}^{n},\nu}=Z^{n}_{t_{\ell-1},\nu}+c^{n}_{\nu}\Big{(}\mathbf{Z}^{ n}_{t_{\ell-1}}\Big{)}+D^{n}_{\ell,\nu}, \tag{3.7}\]

where

\[b^{n}_{\nu}(z) :=-n\alpha e^{-A_{\nu}}C(A_{n})\partial_{\nu}L(z),\] \[c^{n}_{\nu}(z) :=\frac{1}{n}\mathbb{E}b^{n}_{\nu}(z+^{n}\mathbf{V}^{\nu}),\]

and

\[D^{n}_{\ell,\nu} :=\alpha\Big{(}L\big{(}\mathbf{Z}^{n}_{t_{\ell}^{n}}+^{n}\mathbf{ U}^{(\ell)}\big{)}-L\big{(}\mathbf{Z}^{n}_{t_{\ell-1}^{n}}+^{n}\mathbf{U}^{( \ell-1)}\big{)}\Big{)}\Big{(}e^{-u^{(\ell)}_{\nu}}-e^{u^{(\ell)}_{\nu}}\Big{)} -c^{n}_{\nu}\Big{(}\mathbf{Z}^{n}_{t_{\ell-1}^{n}}\Big{)}\]
is a martingale difference process with respect to the filtration \((\mathcal{F}_{\ell})_{\ell=0,\ldots,n}\) that is generated by all randomness up to step \(\ell\). Moreover, the random vector \({}^{n}\mathbf{V}^{n}\) has independent components where all but the \(\nu\)th are uniformly distributed on \([-A_{n},A_{n}]\) and the \(\nu\)th has density

\[f_{A_{n}}(x):=C(A_{n})^{-1}(e^{A_{n}}-e^{x})(e^{A_{n}}-e^{-x})\]

on \([-A_{n},A_{n}]\) with normalizing constant

\[C(A_{n}):=\int_{-A_{n}}^{A_{n}}(e^{A_{n}}-e^{x})(e^{A_{n}}-e^{-x})dx=2A_{n}(e^{ 2A_{n}}+1)+2-2e^{2A_{n}}.\]

We obtain

\[\mathbf{Z}_{t}^{n}-\mathbf{Z}_{t}=\int_{0}^{\lfloor tn\rfloor/n}\Big{(}b_{ \nu}(\mathbf{Z}_{s}^{n})-b_{\nu}(\mathbf{Z}_{s})\Big{)}ds+\mathbf{C}_{\lfloor tn \rfloor/n}^{n}+\mathbf{M}_{\lfloor tn\rfloor/n}^{n}+\delta_{t}^{n},\]

for all \(t\in[0,1]\) where

\[b_{\nu}(z) :=-4\alpha\partial_{t}L(z),\] \[\mathbf{C}_{\lfloor tn\rfloor/n}^{n} :=\sum_{\ell:t_{\ell}^{n}\leq t}\Big{(}\mathbf{c}^{n}\Big{(} \mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}-\frac{1}{n}b_{\nu}\Big{(}\mathbf{Z}_{t _{\ell-1}^{n}}^{n}\Big{)}\Big{)}\] \[=\sum_{\ell:t_{\ell}^{n}\leq t}\Big{(}\mathbf{c}^{n}\Big{(} \mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}-\frac{1}{n}b_{\nu}^{n}\Big{(}\mathbf{Z} _{t_{\ell-1}^{n}}^{n}\Big{)}\Big{)}+\frac{1}{n}\sum_{\ell:t_{\ell}^{n}\leq t} \Big{(}b_{\nu}^{n}\Big{(}\mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}-b_{\nu}\Big{(} \mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}\Big{)},\] \[\mathbf{M}_{\lfloor tn\rfloor/n}^{n} :=\sum_{\ell:t_{\ell}^{n}\leq t}\mathbf{D}_{t_{\ell}^{n}}^{n},\] \[\delta_{t}^{n} :=\mathbf{Z}_{\lfloor tn\rfloor/n}-\mathbf{Z}_{t}.\]

Set

\[\varepsilon_{n}(t):=\sqrt{\mathbb{E}\sup_{s\leq t}\|\mathbf{Z}_{s}^{n}- \mathbf{Z}_{s}\|^{2}}.\]

Using (3.6) and the triangle inequality, we conclude that

\[\varepsilon_{n}(t)\leq\gamma_{n}+4\alpha\lambda\int_{0}^{t}\varepsilon_{n}(s)ds\]

with

\[\gamma_{n}=\sqrt{\mathbb{E}\sup_{\ell=0,\ldots,n}\left\|\mathbf{C}_{t_{n}^{n} }^{n}\right\|^{2}}+\sqrt{\mathbb{E}\sup_{\ell=0,\ldots,n}\left\|\mathbf{M}_{t_ {n}^{(\ell)}}^{n}\right\|^{2}}+\sup_{t\in[0,1]}\|\delta_{t}^{n}\|.\]

By Gronwall's inequality we obtain

\[\varepsilon_{n}(t)\leq\gamma_{n}\exp(\lambda).\]

It is therefore sufficient to prove that \(\gamma_{n}\to 0\).

Note that

\[\left\|b_{\nu}^{n}\Big{(}\mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}-b_{\nu}\Big{(} \mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}\right\|\leq\alpha\Big{|}e^{-A_{n}}nC(A _{n})-4\Big{|}\sup_{z}\|\partial_{\nu}L(z)\|\]

and

\[\left\|c^{n}\Big{(}\mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}-\frac{1}{n}b_{\nu}^ {n}\Big{(}\mathbf{Z}_{t_{\ell-1}^{n}}^{n}\Big{)}\right\|\leq\frac{1}{n}\alpha e ^{-A_{n}}nC(A_{n})\lambda A_{n}\]

because the random vectors \({}^{n}\mathbf{V}^{n}\) are all concentrated on \([-A_{n},A_{n}]\). Since \(e^{-A_{n}}nC(A_{n})\to 4\) and \(A_{n}\to 0\) as \(n\to\infty\), we obtain that \(\mathbb{E}\sup_{\ell=0,\ldots,n}\left\|\mathbf{C}_{t_{n}^{(\ell)}}^{n}\right\| ^{2}\to 0\) as desired. Moreover, we have that
\(\sup_{t\in[0,1]}\|\delta_{t}^{n}\|\to 0\) because \(\mathbf{Z}\) is uniformly continuous as a continuous function on a compact interval.

So it only remains to be verified that \(\mathbb{E}\sup_{\ell=0,\dots,n}\|\mathbf{M}_{t_{\ell}^{n}}^{n}\|^{2}\to 0\). Doob's inequality yields

\[\mathbb{E}\sup_{\ell=0,\dots,n}\left\|\mathbf{M}_{t_{\ell}^{n}}^{n}\right\|^{2} \leq 4\mathbb{E}\left(\sum_{\nu}\sum_{\ell=1}^{n}\left(D_{\ell,\nu}^{n} \right)^{2}\right).\]

Since the gradient of \(L\) is bounded and the components of the \({}^{n}\mathbf{U}^{(\ell)}\) are concentrated on the interval \([-A_{n},A_{n}]=[-n^{1/3},n^{1/3}]\), we have that

\[\left\|\alpha\left(L\left(\mathbf{Z}_{t_{\ell}^{n}}^{n}+{}^{n}\mathbf{U}^{( \ell)}\right)-L\left(\mathbf{Z}_{t_{\ell-1}^{n}}^{n}+{}^{n}\mathbf{U}^{(\ell- 1)}\right)\right)\left(e^{-\nu\mathbf{U}^{(\ell)}}-e^{\nu\mathbf{U}^{(\ell)}} \right)\right\|\leq an^{-2/3}+b\left\|\mathbf{Z}_{t_{\ell}^{n}}^{n}-\mathbf{Z}_ {t_{\ell-1}^{n}}^{n}\right\| \tag{3.8}\]

for some constants \(a,b\in\mathbb{R}_{+}\). By Lemma 2 below and \(\mathbf{Z}_{t_{0}^{n}}^{n}-\mathbf{Z}_{t_{\ell-1}^{n}}^{n}=0\) this implies that (3.8) is bounded by a multiple of \(n^{-2/3}\) and hence its square by \(n^{-4/3}\). Moreover, \(\|c_{\nu}^{n}(\mathbf{Z}_{t_{\ell-1}^{n}}^{n})\|\) is bounded by a multiple of \(1/n\). Together, this yields that \(\sum_{\ell=1}^{n}(D_{\ell,\nu}^{n})^{2}\) is bounded by a multiple of \(n^{-1/3}\), which yields the desired convergence. 

**Lemma 2**.: \(x_{n}\leq a+bx_{n-1}\)_, \(n=1,2,\dots\) for \(a,b,x_{n}\in\mathbb{R}_{+}\) implies_

\[x_{n}\leq x_{0}b^{n}+a\frac{1-b^{n}}{1-b}.\]

Proof.: This follows by induction on \(n\). 
