extfieldsofstudy:
Computer Science
Mathematics

text:
Fast Analog Transmission for High-Mobility Wireless Data Acquisition in Edge Learning

By implementing machine learning at the network edge, edge learning trains models by leveraging rich data distributed at edge devices (e.g., smartphones and sensors) and in return endow on them capabilities of seeing, listening, and reasoning. In edge learning, the need of high-mobility wireless data acquisition arises in scenarios where edge devices (or even servers) are mounted on the ground or aerial vehicles. In this paper, we present a novel solution, called fast analog transmission (FAT), for high- mobility data acquisition in edge-learning systems, which has several key features. First, FAT incurs low-latency. Specifically, FAT requires no source-and-channel coding and no channel training via the proposed technique of Grassmann analog encoding (GAE) that encodes data samples into subspace matrices. Second, FAT supports spatial multiplexing by directly transmitting analog vector data over an antenna array. Third, FAT can be seamlessly integrated with edge learning (i.e., training of a classifier model in this work). In particular, by applying a Grassmannian-classification algorithm from computer vision, the received GAE encoded data can be directly applied to training the model without decoding and conversion. This design is found by simulation to outperform conventional schemes in learning accuracy due to its robustness against data distortion induced by fast fading.

I. INTRODUCTION
Envisioned as an evolution in computing, edge learning refers to the implementation of machine learning at the network edge so as to leverage enormous data distributed at edge devices (e.g., smartphones and sensors) for training models [1], [2]. Subsequently, the models are applied to empowering edge devices with the capabilities of seeing, listening and reasoning. While computing speeds are growing rapidly, the latency in wireless data acquisition has emerged to be the bottleneck of fast edge learning [2]. This issue is exacerbated in high-mobility scenarios where edge devices (or even edge servers) are mounted on ground or aerial vehicles as illustrated Fig. 1(a) [3]. High-mobility data acquisition faces several challenges: 1) robustness against fast fading, 2) low-latency given short connection time, 3 The design of the FAT scheme builds on three ideas from communication and learning. The first idea is analog transmission based on linear analog modulation that has been deployed previously in different settings such as fast transfer of channel-state information (CSI) [4] and over-the-air functional computation in sensor networks [5]. Compared with digital transmission, the analog design does not require source-end-channel coding and decoding, thereby reducing computation complexity. Moreover, direct transmission of analog data instead of a quantized bit stream shortens the transmission duration. In terms of learning performance, our findings suggest that a customized design of analog transmission targeting learning (e.g., FAT) can be more robust than digital counterparts against data distortion by fast fading at high mobility.
The second idea is blind multiple-input-multiple-output (MIMO) transmission without CSI.
This idea was first developed in the classic area of non-coherent MIMO, which is a digital spacetime modulation scheme [6], [7]. Its unique feature is a modulation constellation comprises a set of subspace matrices. The transmitted space-time symbol in the form of such a matrix is invariant to rotation by a block fading channel that remains constant within each symbol duration but varies over different durations. Thus the matrix can be transmitted and detected even without CSI at either side, referred to hereafter as the channel-invariant property [6]. Consequently, channel training is unnecessary, thereby reducing the transmission latency and overhead as illustrated in Fig. 1(b). On the other hand, non-coherent MIMO cannot support spatial multiplexing like its coherent counterpart. The resultant low data rates makes the former less popular in practice and its applications are limited to low-rate ultra-fast machine-type applications [3], [8]. In contrast, the proposed scheme retains the advantages of both technologies, namely channel-invariant property and (analog) spatial multiplexing. The said property is achieved by the proposed Grassmannian analog encoding (GAE), a key FAT component, which encodes a data sample (an analog vector)  The differences between FAT and conventional MIMO schemes are summarized in Table I.
In this paper, we consider a typical edge-learning task of training a classifier model. The last idea pertains to edge learning and is to apply a Grassmann classification algorithm for classifying the received GAE encoded training data. Such algorithms were originally developed for computer vision where image features or motions are represented as subspaces or equivalently points on a Grassmann manifold, referred to as Grassmann data [9]. Via the application of such an algorithm, classification can be seamlessly integrated with FAT since the received GAE encoded data can be directly used in learning without decoding and conversion. Furthermore, the integration leads to accurate edge learning with robustness in data acquisition against fast fading.
In summary, an edge learning system based on FAT comprises the following three components (elaborated in Section III).
1) Grassmann analog encoding: At each edge device, the proposed GAE encodes data samples into subspace matrices by projection onto a Grassmann manifold to enable blind MIMO transmission and robust edge learning.
2) Analog transmission and detection: The GAE encoded data is transmitted using linear analog modulation and blindly detected at the edge server without channel knowledge.
3) Edge learning: At the edge server, the received Grassmann data is used for training a classifier model using a Grassmann-classification algorithm from computer vision.
By evaluating the classification performance of a model and transmission latency using simulation (see Section V), the proposed FAT scheme is found to substantially outperform the conventional coherent (analog and digital) MIMO transmission at high mobility.
Trained models Figure 2. An edge-learning system based on FAT.

A. System Model
Consider the edge-learning system illustrated in Fig. 2 where an edge server trains a classifier using a training dataset transmitted by multiple edge devices. The transmissions by devices are based on time sharing and independent of channels given no CSI. All nodes are equipped with antenna arrays, resulting in a set of narrow-band MIMO channels. Let N t and N r denote the numbers of transmit and receive antennas, respectively. Time is divided into baseband sampling intervals, called (time) slots. Then the slot-t realization of the MIMO channel from an active device to the server can be represented by the N r ×N t matrix H t . Given an analog vector-symbol g t transmitted by the active device, the received signal at the server is where P is the transmission power and w t the additive-white-Gaussian-noise (AWGN) vector.
In this work, we focus on transmission of data samples that dominates the data acquisition process. Their labels have finite values and naturally can be transmitted using digital non-coherent MIMO modulation over a low-rate channel, called label channel, orthogonal to the high-rate data channel. Due to its low rate, the label channel can be reasonably assumed to be noiseless similarly as the CSI feedback channel [10].

B. Simulation Models
Simulation for evaluating learning performance is based on the following data and channel models. The data at different edge devices are assumed to be independent and identically distributed (i.i.d) based on the classic mixture of Gaussian (MoG) model, which is widely adopted in the machine-learning literature. Each data sample is a 1-by-L complex random vector. Let M denote the number of data classes. Then the i-th sample, denoted as s (i) , from the m-th class can be modelled as where µ m is the mean of the m-th class and z (i) ∈ C 1×L a deviation vector comprising i.i.d.
CN (0, σ 2 s ) elements. Next, high mobility induces temporally correlated MIMO channels. Assuming rich scattering, the classic Clark's model is applied that translates a speed into the level of channel temporal correlation. Specifically, within the duration of transmitting a data sample, two realizations of the (m, n)-th coefficient of the channel H t separated by τ slots are correlated with the correlation function given as where f D = fcv c with v being the speed, f c carrier frequency and c speed of light, and J 0 is the zero-th order Bessel function of the first kind.

III. FAST ANALOG TRANSMISSION SCHEME
In this section, we discuss two key algorithms in the proposed FAT scheme, namely GAE and blind analog transmission and detection (see Fig. 2). The received Grassmannian dataset is used for training a classifier model using an existing Grassmannian classifiaction algorithm such as sample Karcher mean [11], which is adopted in simulation. The details are omitted for brevity.

A. Grassmann Analog Encoding
To facilitate exposition, some mathematical notions are defined as follows. The (n, m) Grassmann manifold is a set of all m-dimensional subspaces in C n , denoted by G n,m [12]. For the special case of G 3,1 , each point on the manifold geometrically corresponds to a unique line passing through the origin as illustrated in Fig. 3. For ease of notation, a point on G n,m that is a subspace is usually represented by an arbitrary basis matrix spanning the subspace, denoted as Υ. The subspace distance between two points Υ and Υ on the Grassmannian G n,m , denoted as d p (Υ, Υ ), is measured using the commonly used metric of Procrustes distance for its better performance in simulation: As discussed, GAE at the active device endows on FAT the channel-invariant property, thereby enabling blind analog transmission with robustness against fast fading. As illustrated in Fig. 3, the mathematical principle of GAE is to project original data samples (vectors in the Euclidean space) onto the Grassmann manifold, generating subspace matrices as the output. The GAE algorithm is described as follows.

Data in the Euclidean space
Data on the Grassmannian Grassmann manifold Figure 3. Principle of Grassmann analog encoding.
Step 1 (Vector-to-matrix conversion): Consider a data sample that is a 1 × L row vector, say s (i) , with L being an integer multiple of the number of transmit antennas N t . Then s (i) can be divided into 1 × T sub-vectors with T = L/N t : . It follows that s (i) can be converted into a N t × T data matrix X (i) having {s (i) n } as rows. The matrix X (i) such constructed is typically fat (N t < T ) since a data-sample vector is usually long (L N t ). For the case where L is not an integer multiple of N t , zero-padding can be applied to lengthen s (i) so that the integer-multiple constraint is met.
Step 2 (Projection onto Grassmannian): The key step in encoding is to project the matrix X (i) constructed in the preceding step onto a single point on the Grassmannian G T,Nt . To this end, decompose the matrix X (i) by singular value decomposition (SVD) as Then G (i) is a N t × T basis matrix spanning the row space of X (i) and thus a point on G T,Nt .
The encoder uses G (i) as the output from encoding the data-sample s (i) .
In summary, the advantages of GAE are threefold: 1) enabling blind transmission and detection as discussed in the next sub-section, 2) endowing on edge learning robustness against data distortion by fast fading as shown in simulation results, and 3) allowing seamless integration with learning on Grassmannian without signal decoding and conversion.

B. Blind Analog Transmission and Detection
Given Grassmann encoding, the procedures for subsequent blind analog transmission and detection in FAT are described as follows.
Step 1 (Analog transmission): After encoding each data-sample, say s (i) , into a subspace basis matrix G (i) , the N t × T matrix G (i) is directly transmitted by the active device over T slots using linear analog modulation and the array of N t antennas. To reflect channel temporal variation, it is necessary write G (i) in terms of its columns: T ]. Then the received signal due to the transmission of G (i) can be represented by the N r × T matrix For continuous time-shared distributed uploading of total N data samples, t = 1, 2, · · · N T . It is important to observe from (5) that due to high mobility, the channel {H (i) t } varies in the T -slot transmission duration of a single data sample, which has a negative effect on decoding as discussed in the sequel.
Step 2 (Grassmann analog detection): The detection of the transmitted encoded analog space-time symbol G (i) involves the extraction of the row space, denoted by the N t × T unitary matrix G (i) , from the SVD of the received N r × T space-time signal Y (i) specified in (5), Consider the special case of zero noise and static channel. The In other words, G (i) and G (i) are the identical point on the Grassmannian, corresponding to perfect detection.
In the presence of noise and channel variation, they are two different points and the resultant detection error affects learning. Based on the above detection procedure, the output training dataset, called Grassmann dataset, is a sequence of N labeled subspace matrices (points on the Grassmannian), [ G (1) , G (2) , · · · , G (N ) ], whose labels are acquired by the server via the said low-rate label channel.
Remark 1 (Blind Transmission and Detection). Both the analog transmission and detection in the above steps are independent of the channel. In particular, the detection of Grassmann dataset involves SVDs of the received array observations that do not require any channel knowledge.

A. Grassmann Analog Encoding Preserves Clustering
An important reason FAT supports edge classification is that GAE retains the class structure in the original dataset. This property is illustrated in Fig. 4 where the high-dimensional datasets are visualized in the 2D plane using a well known visualization algorithm, t-distributed stochastic neighbour embedding (t-SNE). As discussed in the sequel, GAE incurs DoF loss in the dataset.
Consequently, one can observe form Fig. 4 that data classes are less compact after GAE, sacrificing some level of discriminant of the dataset. The loss, nevertheless, yields communication advantages discussed shortly.

B. Trading DoF Loss for Robustness and Low Latency
The Grassmann encoding design leads to the DoF loss (or discriminant loss), which may make data points among different classes that are well-separated in the Euclidean space become much closer or even overlapped on the Grassmannian. The loss has a negative effect on the classification performance. The phenomenon is illustrated by the following example. Besides SVD, an alternative method for GAE is LQ decomposition. Consider the LQ decomposition of two data matrices X = L X G X and Y = L Y G Y , where the unitary matrices G X and G Y represent identical subspaces (or identical encoding outputs) as the SVD counterparts and L X and L Y are lower triangular matrices. If span(G X ) = span(G Y ) but L X = L Y , the Euclidean distance between X and Y is d 2 E (X, Y) = 0. However, based on (4), the Procrustes distance between the GAE encoded data samples, namely G X and G Y , is d 2 p (G X , G Y ) = 0. A key finding in this work is that the DoF loss of GAE is more than compensated by its robustness against fast fading that can cause severe errors in data transmission without GAE.
As a result, GAE leads to a net performance gain over conventional schemes at high mobility. The performance of FAT is benchmarked against two high-rate coherent schemes: digital and analog MIMO transmission, both of which assume a MMSE linear receiver and thus require channel training to acquire the needed CSI. Like FAT, analog MIMO transmits data samples directly by linear analog modulation. On the other hand, digital MIMO quantizes data samples into 8-bit per coefficient and modulates each symbol using QPSK before MIMO transmission.
All considered schemes have no error control coding.

A. Communication Latency Performance
While FAT is free of channel-training, benchmark schemes incur training overhead that can be quantified by the fraction of a frame allocated for the purpose i.e., the ratio P/(P + D) with P and D illustrated in Fig. 1(b). The curves of overhead versus Doppler shift are displayed in Fig. 5 for FAT and two mentioned benchmarking schemes for a given classification-error rate of 1 × 10 −3 . One can observe that the overhead grows monotonically with the Doppler shift as the channel fading becomes faster. For high-mobility with Doppler approaching 10 −2 , the overhead can be more than 12% and 6% for digital and analog coherent MIMO, respectively. Furthermore, given the same performance, digital coherent MIMO (with QPSK modulation and 8-bit quantization) requires 4 times more frames for transmitting the training dataset than the two analog schemes. This suggests that analog transmission is preferable for data acquisition targeting edge learning.

B. Learning Performance
Classifier models discussed in Section III are trained using the training dataset acquired using different transmission schemes and then evaluated using the test dataset. The resultant classification error rates are compared in Fig. 6 by varying Doppler shift and average transmit SNR. Several observations can be made. In the range of moderate to large Doppler shift (i.e., larger than 6 × 10 −3 ), the proposed FAT outperforms the benchmarking schemes, supporting the former's intended application in high-mobility data acquisition. Furthermore, at high mobility (i.e., Doppler equal to 0.01), FAT achieves the best performance in the practical SNR range (0 − 15 dB). On the other hand, analog and digital coherent MIMO are preferred at low and high SNRs, respectively. The above observations reconfirm the conclusion from preceding latency comparison that analog transmission (especially FAT) is a promising solution for high-mobility data acquisition for edge learning.


text:
Nonconvex Variance Reduced Optimization with Arbitrary Sampling

We provide the first importance sampling variants of variance reduced algorithms for empirical risk minimization with non-convex loss functions. In particular, we analyze non-convex versions of SVRG, SAGA and SARAH. Our methods have the capacity to speed up the training process by an order of magnitude compared to the state of the art on real datasets. Moreover, we also improve upon current mini-batch analysis of these methods by proposing importance sampling for minibatches in this setting. Surprisingly, our approach can in some regimes lead to superlinear speedup with respect to the minibatch size, which is not usually present in stochastic optimization. All the above results follow from a general analysis of the methods which works with arbitrary sampling, i.e., fully general randomized strategy for the selection of subsets of examples to be sampled in each iteration. Finally, we also perform a novel importance sampling analysis of SARAH in the convex setting.

Introduction
Empirical risk minimization (ERM) is a key problem in machine learning as it plays a key role in training supervised learning models, including classification and regression problems, such as support vector machine, logistic regression and deep learning. A generic ERM problem has the finite-sum form where x corresponds to the parameters/features defining a model, f i (x) is the loss of the model associated with data point i, and f is the average (empirical) loss across the entire training dataset. In this paper we will focus on the case when the functions f i are L i -smooth but non-convex. We assume the problem has a solution x * . One of the most popular algorithms for solving (1) is stochastic gradient descent (SGD) [19,18]. In recent years, tremendous effort was exerted to improve its performance, leading to various enhancements which use acceleration [1], momentum [14], minibatching [35], distributed implementation [16,15], importance sampling [36,6,27,4], higher-order information [26,9], and a number of other techniques.
process in various clever ways so as to progressively reduce its variance as an estimator of the true gradient. We note that SAG and SARAH, the historically oldest and newest VR methods in the list, respectively, use a biased estimator of the gradient. In theory, all these methods enjoy linear convergence rates on smooth and strongly convex functions, which is in contrast with slow sublinear rate of SGD. VR methods are also easier to implement as they do not rely on a decreasing learning rate schedule. VR methods were recently extended to work with (non-strongly) convex losses [13], and more recently also to non-convex losses [29,28,2,23], in all cases leading to best current rates for (1) in a given function class.

Importance sampling, minibatching and non-convex models
In the context of problem (1), importance sampling refers to the technique of assigning carefully designed non-uniform probabilities {p i } to the n functions {f i }, and using these, as opposed to uniform probabilities, to sample the next data point (stochastic gradient) during the training process.
Despite the huge theoretical and practical success of VR methods, there are still considerable gaps in our understanding. For instance, an importance sampling variant of the popular SAGA method, with the "correct" convergence rate, was only designed very recently [10]; and the analysis applies to strongly convex f only. A coordinate descent variant of SVRG with importance sampling, also in the strongly convex case, was analyzed in [12]. However, the method does not seem to admit a fast implementation. For dual methods based on coordinate descent, importance sampling is relatively well understood [20,32,24,27,3].
The territory is completely unmapped in the non-convex case, however. To the best of our knowledge, no importance sampling VR methods have been designed nor analyzed in the popular case when the functions {f i } are non-convex. An exception to this is dfSDCA [5]; however, this method applies to an explicitly regularized version of (1), and while the individual functions are allowed to be non-convex, the average f is assumed to be convex. Given the dominance of stochastic gradient type methods in training large non-convex models such as deep neural networks, theoretical investigation of VR methods that can benefit from importance sampling is much needed.
The situation is worse still when one asks for importance sampling of minibatches. To the best of our knowledge, there are only a handful of papers on this topic [30,6], none of which apply to the non-convex setting considered here, nor to the methods we will analyze, and the problem is open. This is despite the fact that minibatch methods are de-facto the norm for training deep nets due to the volume of data that feeds into the training, and the necessity to fully utilize the parallel processing power of GPUs and other hardware accelerators for the task. In practice, typically relatively small (O(1) in comparison with n) minibatch sizes are used.

Contributions
The main contributions of this paper are: • Arbitrary sampling. We peform a general analysis of three popular VR methods-SVRG [11], SAGA [7] and SARAH [22]-in the arbitrary sampling paradigm [30,24,25,27,4]. That is, we prove general complexity results which hold for an arbitrary random set valued mapping (aka arbitrary sampling) generating the minibatches of examples used by the algorithms in each iteration.
• Rates. Our bounds improve the best current rates for these methods, for SVRG and SAGA even under the same L i 's. Our importance sampling can be up faster by factor of n comparing to the current state of the art (see Table 1 and Appendix C). Our methods can enjoy linear speedup or even for some specific samplings superlinear speedup in minibatch size. That is, the number of iterations needed to output a solution of a given accuracy drops by a factor equal or greater to the minibatch size used. This is of utmost relevance to the practice of training neural nets with minibatch stochastic methods as our results predict that this is to be expected. We design importance sampling and approximate importance sampling for minibatches which in our experiments vastly outperform the standard uniform minibatch strategies.
• Best rates for SARAH under convexity. Lastly, we also perform an analysis of importance sampling variant of SARAH in the convex and strongly convex case (Appendix I). These are the currently fastest rates for SARAH.

Importance Sampling for Minibatches
As mentioned in the introduction, we assume throughout that f i : R d → R are smooth, but not necessarily convex. In particular, we assume that f i is L i -smooth; that is, Without loss of generality assume that L 1 ≤ L 2 ≤ · · · ≤ L n .

Samplings
Let S be a random set-valued mapping ("sampling") with values in 2 [n] , where [n] def = {1, 2, . . . , n}. A sampling is uniquely defined by assigning probabilities to all 2 n subsets of [n]. With each sampling we associate a probability matrix P ∈ R n×n defined by The probability vector associated with S is the vector composed of the diagonal entries of P: From now on, we will refer to b as the minibatch size of sampling S. It is known that P − pp is a symmetric positive semidefinite matrix [31]. Let us without loss of generality assume that p 1 ≤ p 2 ≤ · · · ≤ p n and define constant k = |{i ∈ [n] : p i < 1}| = max{i : p i < 1} to be number of p i 's, which are not equal to one.
While our complexity results are general in the sense that they hold for any proper sampling, we shall now consider three special samplings; all with minibatch size b ∈ (0, n]: (i) Standard uniform minibatch sampling (S = S u ). 1 S is chosen uniformly at random from all subsets of [n] of cardinality b. Clearly, |S| = b with probability 1. The probability matrix is given by (ii) Independent sampling (S = S * ). Assume any proper p i 's . For each i ∈ [n] we independently flip a coin, and with probability p i include element i into S. Hence, by construction, p i = Prob(i ∈ S) and

E [|S|]
(2) (iii) Approximate independent sampling (S = S a ). Independent sampling has the disadvantage that k coin tosses 2 need to be performed in order to generate the random set. However, we would like to sample at the cost O(b + k − n) coin tosses instead. We now design a sampling which has this property and which in a certain precise sense, as we shall see later, approximates the independent sampling. In particular, given an independent sampling with parameters . We now sample a single set S of cardinality a using the standard uniform minibatch sampling (just for i ≤ k). Subsequently, we apply an independent sampling to select elements of S , with selection probabilities kpi a npj a = (a−1)k a(k−1) p i p j for i, j ≤ k, the probability matrix of S is given by Since (a−1)k a(k−1) ≈ 1, the probability matrix of the approximate independent sampling approximates that of the independent sampling. Note that S includes both the standard uniform minibatch sampling and the independent sampling as special cases. Indeed, the former is obtained by choosing p i = b/n for all i (whence a = b and p i = 1 for all i), and the latter is obtained by choosing a = n instead of a = k max i≤k p i .

Key lemma
The following lemma, which we use as upper bound for variance, plays a key role in our analysis. Lemma 1. Let ζ 1 , ζ 2 , . . . , ζ n be vectors in R d and letζ def = 1 n n i=1 ζ i be their average. Let S be a proper sampling (i.e., assume that p i = Prob(i ∈ S) > 0 for all i). Assume that there is v ∈ R n such that Then 1 In the literature, this is often refered to by the name b-nice sampling [31,25]. 2 Note that just k not n, because others are included with probability one.
where the expectation is taken over sampling S. Whenever (4) holds, it must be the case that Moreover, (4) is always satisfied for v i = n(1 − p i ) for i ≤ k and 0 otherwise. Further, if |S| ≤ d with probability 1 for some d, then (4) holds for v i = d. The standard uniform minibatch sampling admits v i = n−b n−1 , the independent sampling admits v i = 1 − p i , and the approximate independent sampling admits the choice The following quantities, which comes from Lemma 1 and smoothness assumption, play a key role in our general complexity results: We can see through theory that it is a good idea to design samplings for which the value α is the smallest possible, which would lead to optimal sampling. The following result sheds light on how S should be chosen, from samplings of a given minibatch size b, to minimize α.
Lemma 2. Fix a minibatch size b ∈ (0, n]. Then the quantity α, defined in (7), is minimized for the choice S = S * with the probabilities where k is the largest integer satisfying where s = 1 for S * and s = 1 − k−a a(k−1) for S a . Moreover, if we assume 3 bL n ≤ n i=1 L i , then k = n, thus From now let S * , S a denote Independent Sampling and Approximate Inpedendent Sampling, respectively, with probabilities defined in (8).
Remark 1. Lemma 2 guarantees that the sampling S * is optimal. Moreover, let b max def = max{b|bL n ≤ n i=1 L i } then in theory, we have super linear speed up for b up to b max for all three algorithms.

SVRG, SAGA and SARAH
In all of the results of this section we assume that S is an arbitrary proper sampling. Let b = E [|S|] be the (average) minibatch size. We assume that v satisfies (4) and that α (which depends on v) is defined as in (7). All complexity results will depend on α and b.
We propose three methods, Algorithm 1, 2 and 3, which are generalizations of original SVRG [28], SAGA and SARAH to the arbitrary sampling setting, respectively. The original non-minibatch methods arise as special  Algorithm 2 SAGA with arbitrary sampling x 0 , d, T, η, S cases for the sampling S = {i} with probability 1/n, and the original minibatch methods arise as a special case for the sampling S u (described in Section 2.1). Our general result for SVRG [11] follows.
Theorem 5 (Complexity of SARAH with arbitrary sampling). Consider one outer loop of Alg. 3 with Then the output x a satisfies: . Thus in terms of stochastic gradient evaluations to obtain -accurate solution, one needs following number of iterations If all L i 's are the same and we choose S to be S a , thus uniform with mini-batch size b, we can get back original result from [23]. Taking b = n, we can restore gradient descent with the correct step size.

Gradient dominated functions
Gradient dominance is weaker version of strong convexity due to the fact that if function is µ-strongly convex then it is τ -gradient dominated, where τ = 2/µ.
Any of the non-convex methods in this paper can be used as a subroutine of Algorithm 4, where T is the number of steps of the subroutine and A is the set of optimal parameters for the subroutine. We set T = αn 2/3 /(bν 2 ) for SVRG and T = αn 2/3 /(bν 3 ) for SAGA. In the case of SARAH, T is obtained by solving m + 1 = 2/η in m and setting T ← m. Using Theorems 3, 4, 5 and the above special choice of T , we get Combined with Definition 1, this guarantees linear convergence with the same constants α,L and b we had before in our analysis.

Importance sampling for SARAH under convexity
In addition to the results presented in previous sections, we also establish importance sampling results for SARAH in convex and strongly convex cases (Appendix I) with similar improvements as for the non-convex algorithm. Ours are the best current rates for SARAH in these settings.
Further, we also provide specialized non-minibatch versions of non-convex SAGA, SARAH and SVRG, which are either special cases of their minibatch versions presented in the main part, or slight modifications, with slightly improved guarantees.

Experiments
In this section, we perform experiments with regression for binary classification, where our loss function has form 1 is the sigmoid function, thus is smooth but non-convex. We use four LIBSVM datasets 4 : covtype, ijcnn1, splice, australian.
Parameters of each algorithm are chosen as suggested by theorems in section 3 and x 0 is set to be all zeros vector. The y axis shows the norm of the gradient, and the x axis shows how many times the algorithm evaluates the full gradient. Evaluation of the gradient of a single function from the empirical loss costs 1/n of the full gradient evaluation. For SARAH, we chose m = n/b .

Importance and uniform sampling comparison
Here we provide comparison of the methods with uniform S u and with importance S * sampling.

Linear or more than linear speedup
Theory suggests that linear or even more than linear speed up can be obtained using S * . Our experiment suggest that this is possible for all three algorithms.  2 confirms that linear, even superlinear, speedup (with increasing batch size, we need the same or less full gradient evaluations) can be obtained also in practice, not just in theory. However, it is limited, and for this dataset it is up to minibatch size of 250. Upper graphs at Figure 2 visualize convergence under assumption of multi-core settings, where we assume number of cores is the same as the minibatch size.

Independent Sampling vs. Approximate Independent Sampling
In the theory, Independent Sampling S * is slightly better than Approximate Independent Sampling S a . However, it is more expensive in terms of computations. The goal of the next experiment is to show that in practice S a yields comparable or faster convergence. Hence, it is more reasonable to use this sampling for datasets, where number of data points n is big (if we implement S a efficiently we can almost get rid of dependence on n). The intuition behind why S a could work better is that S a has smaller variance in batch size than S * . It can be seen from Figure 3 that S a can outperform S * , thus, however, S * is optimal in theory, one should use S a in practice.

Appendix
Let A = [a 1 , . . . , a n ] ∈ R d×n , where a i = ζi pi , and let e be the vector of all ones in R n . We now write the variance of X in a form which will be convenient to establish a bound: Since by assumption we have P − pp Diag (p • v), we can further bound To obtain (5), it remains to combine this with (12). Inequality (6) follows by comparing the diagonal elements of the two matrices in (4). Let us now verify the formulas for v.
• It was shown in [25] • Consider now the independent sampling. Clearly, where v i = 1 − p i .
• Consider the b-nice sampling (standard uniform minibatch sampling). Direct computation shows that the probability matrix is given by as claimed in (3). Therefore, the probability matrix of the approximate independent sampling satisfies

B Proof of Theorem 2
We first establish a lemma we will need in order to prove Theorem 2.

Let be the largest integer for which
Li L k (note that the inequality holds for k = n − b + 1). Then (13) has the following solution: Proof. The Lagrangian of the problem is L(p, y, λ 1 , .., λ n , u 1 , ..., Th constraints are linear and hence KKT conditions hold. The result can be deduced from the KKT conditions.
We can now proceed with the proof. Since n, b andL are constants, the problem is equivalent to In view of (6), where function Ω(p) was defined in Lemma 6. Since b = E [|S|] = i p i , and 0 ≤ p i ≤ 1 for all i, then in view of Lemma 6 we have where p * is defined by (8).
On the other hand, from Lemma 1 we know that the independent sampling S = S * with probability vector p * defined in (8) satisfies inequality (4) with v i = 1 − p i , and hence Hence, it is optimal.

C Improvements
Let us compute α for uniform sampling.
It is easy to see that L max ≥L. To prove that we have improved current best known rates, we need to show thatLα ≤ L max andL 2 α ≤ (n−b) , then L max ≈ nL and α S * ≤ 1 and α S u ≈ n, which essentially means, that we can have in theory speed up by factors of n.

D.1 SVRG
For SVRG, each outer loop costs n + mb evaluations of stochastic gradient. If we want to obtain -solution, following must hold (Theorem 3) Combining these two equations with definition from Theorem 3, we get total complexity in terms of stochastic gradients evaluation

D.2 SAGA
For SAGA, each loop costs d + b evaluations of stochastic gradient. If we want to obtain -solution, following must hold (Theorem 4) Combining these two equations with definition from Theorem 4, we get total complexity in terms of stochastic gradients evaluation because of evaluation of full gradient on the start.

D.3 SARAH
For SARAH with one outer loot, each inner loop costs 2b evaluations of stochastic gradient. If we want to obtain -solution, following must hold (Theorem 5) Solving this equation for m, we get Combining thise equation with complexity off each inner loop we obtain total complexity in terms of stochastic gradients evaluation E Proofs for SVRG Lemma 7. For c t , c t+1 , β > 0, suppose we have Let η, β and c t+1 be chosen such that Γ t > 0 (in Theorem (17)). The iterate x s+1 t in Algorithm 5 satisfy the bound: Proof. Since f i is L i -smooth we have Summing through all i and dividing by n we obtain Using the IS-SVRG update in Algorithm 5 and its unbiasedness ( ), the right hand side above is further upper bounded by Consider now the Lyapunov function For bounding it we will require the following: The second equality follows from the unbiasedness of the update of IS-SVRG. Plugging Equation (15) and Equation (16) into R s+1 t+1 , we obtain the following bound: To further bound this quantity, we use Lemma 10 to bound E v s+1 t 2 , so that upon substituting it in Equation (17), we see that The second inequality follows from the definition of c t and R s+1 t , thus concluding the proof.

Proof of Lemma 7 and Theorem 17
Proof. Since η = η for t ∈ {0, . . . , m − 1}, using Lemma 7 and telescoping the sum, we obtain This inequality in turn implies that where we used that R s+1 The above inequality used the fact thatx 0 = x 0 . Using the above inequality and the definition of x a in Algorithm 5, we obtain the desired result.

Proof of Theorem 3
The proofs essentially follow along the lines of Lemma 7, Theorem 17 and Theorem 18 with the added complexity of mini-batch. We first prove few intermediate results before proceeding to the proof of Theorem 3.

Lemma 8. Suppose we have
for 0 ≤ s ≤ S − 1 and 0 ≤ t ≤ m − 1 and the parameters η, β and c t+1 are chosen such that Then the iterates x s+1 t in the mini-batch version of Algorithm 5 i.e., Algorithm 1 with expected mini-batch size b satisfy the bound: Proof. Using essentially the same argument as the proof of Lemma 7 until Equation (17), we have We use Lemma 11 in order to bound E v s+1 t 2 in the above inequality. Substituting it in Equation (26), we see that The second inequality follows from the definition of c t and R s+1 t , thus concluding the proof.
The following theorem provides convergence rate of mini-batch IS-SVRG.
Theorem 9. Let γ n denote the following quantity: for t ∈ {0, . . . , m − 1} and γ n > 0. Then for the output x a of mini-batch version of Algorithm 5 with mini-batch size b, we have where x * is an optimal solution to (1).
Proof. Since η = η for t ∈ {0, . . . , m − 1}, using Lemma 8 and telescoping the sum, we obtain This inequality in turn implies that . Now sum over all epochs and using the fact thatx 0 = x 0 , we get the desired result.
We now present the proof of Theorem 3 using the above results.
Proof. The proof simply follows from the proof of Lemma 11 with S t = {i t }.
We now present a result to bound the variance of mini-batch IS-SVRG.
Lemma 11. Let v s+1 t be computed by the mini-batch version of Algorithm 5 i.e., Algorithm 1 with sampling S. Then, Proof. For the simplification, we use the following notation: We use the definition of v s+1 t to get The first inequality follows from fact that x + y 2 ≤ 2 x 2 + 2 y 2 and the fact that E ζ s+1 t = ∇f (x s+1 t ) − ∇f (x s ). From the above inequality, we get G Proofs for SAGA Lemma 12. For c t , c t+1 , β > 0, suppose we have Also let η, β and c t+1 be chosen such that Γ t > 0. Then, the iterates {x t } of Algorithm 6 satisfy the bound Proof. Since f isL-smooth we have We first note that the update in Algorithm 6 is unbiased i.e., E [v t ] = ∇f (x t ). By using this property of the update on the right hand side of the inequality above, we get the following: Here we used the fact that x t+1 − x t = −ηv t (see Algorithm 2). Consider now the Lyapunov function For bounding R t+1 we need the following: The above equality follows from the definition of α t+1 i and the definition of randomness of index j t in Algorithm 6 and Algorithm 2. The term T 1 in (32) can be bounded as follows The second equality again follows from the unbiasedness of the update of IS-SAGA. The last inequality follows from a simple application of Cauchy-Schwarz and Young's inequality. Plugging (31) and (33) into R t+1 , we obtain the following bound: where we use that To further bound the quantity in (34), we use Lemma 13 to bound E v t 2 , so that upon substituting it into (34), we obtain R t+1 The second inequality follows from the definition of c t i.e., c t = c t+1 1 − d n + ηβ and R t specified in the statement, thus concluding the proof.
The following lemma provides a bound on the variance of the update used in Mini-batch IS-SAGA algorithm. More specifically, it bounds the quantity E v t 2 .
Lemma 13. Let v t be computed by Algorithm 2. Then, Proof. For ease of exposition, we use the notation Using the convexity of · 2 and the definition of v t we get The first inequality follows from the fact that a + b 2 ≤ 2( a 2 + b 2 ) and that E [ The last inequality follows from L i -smoothness of f i and using properties of S sampling, thus concluding the proof.

Proof of Theorem 19
Proof. We apply telescoping sums to the result of Lemma 12 to obtain The first inequality follows from the definition of γ n . This inequality in turn implies the bound where we used that ). Using inequality (39), the optimality of x * , and the definition of x a in Algorithm 6, we obtain the desired result.
Using the above upper bound on c t we can conclude that upon using the following inequalities: (i) c t+1 η/β ≤ η/4, (ii) η 2 L ≤ η/3 and (iii) 2c t+1 η 2 ≤ η/6, which hold due to the upper bound on c t in (40) and if b ≤ K /L 2 n 2/3 . Substituting this bound on γ n in Theorem 19, we obtain the desired result.
Theorem 20 is special case with b = 1 and d = 1.

SARAH-non-convex
This lemmas are modification of lemmas appeared in [23] for importance sampling with mini-batch.
where x * is an optimal solution of (1).
Proof. ByL-smoothness of f and x t+1 = x t − ηv t , we have where the last equality follows from the fact r q = 1 2 r 2 + q 2 − r − q 2 , for any r, q ∈ R d . By summing over t = 0, . . . , m, we have which is equivalent to (η > 0): where the last inequality follows since x * is an optimal solution of (1). (Note that x 0 is given.)

Lemma 15.
Consider v t defined in SARAH, then for any t ≥ 1, Proof. Let F j = σ(x 0 , i 1 , i 2 , . . . , i j−1 ) be the σ-algebra generated by x 0 , i 1 , i 2 , . . . , i j−1 ; F 0 = F 1 = σ(x 0 ). Note that F j also contains all the information of x 0 , . . . , x j as well as v 0 , . . . , v j−1 . For j ≥ 1, we have where the last equality follows from By taking expectation for the above equation, we have Note that ∇f (x 0 ) − v 0 2 = 0. By summing over j = 1, . . . , t (t ≥ 1), we have With the above Lemmas, we can derive the following upper bound for E ∇f (x t ) − v t 2 .
Consider v t defined in SARAH. Then for any t ≥ 1, Proof. Let We have Hence, by taking expectation, we have By Lemma 15, for t ≥ 1, This completes the proof.

Proof of Theorem 5
Proof. By Lemma 16, we have We have is a root of equation Therefore, by Lemma 14, we have If x a is chosen uniformly at random from {x t } m t=0 , then This concludes the proof.
Comparing Theorem 17 to the previous result in [28], we can see improvement in constant, if we assume different L i -smooth constants for different functions. If the all L i 's are the same then our result is the same as previous result for uniform sampling, because then α = n−1 n−1 = 1.

H.2 SAGA
Here, we provide similar analysis as for SVRG with the same result. We provide more generalized improved form of theorems which appeared in [29]. Randomly uniformly pick it from [n] x t+1 = x t − ηv t 8: α t+1 j t = x t and α t+1 j = α t j for j = jt 9: g t+1 = g t − 1 n (∇fj t (α t j t ) − ∇fj t (α t+1 j t )) 10: end for 11: Output: Iterate xa chosen uniformly random from {x t } T t=0 .
We can see that exactly same conclusions apply here as for SVRG and results can be interpreted in the same way.

I.1 Main result
Consider Algorithm 7, which is an importance sampling variant of the SARAH method.
Note, that only 10-th and 11-th row are changed comparing to classic SARAH algorithm presented in [22]. We do not sample uniformly anymore and also in the 11-th row of Algorithm 7, where we use factor 1 npi in order to stay unbiased in outer cycle.
Then using similar analysis used in [22] and additional lemmas we can prove following theorems with p i in Algorithm 7 to be Li n j=1 Li Theorem 21. Suppose that f i (x) are L i -smooth and convex, f (x) is µ strongly convex. Consider v t defined in SARAH-ISc (Algorithm 7) with η < 2/L, whereL = 1 n n j=1 L i . Then, for any t ≥ 1, By choosing η = O(1/L), we obtain the linear convergence of v t 2 in expectation with the rate (1 − 1/κ 2 ), where κ =L µ is condition number, This is improvement over previous result in [22], because ofL µ ≤ Lmax µ . Below we show that a better convergence rate could be obtained under a stronger convexity assumption for each single f i (x). v t = 1 npi (∇f it (x t ) − ∇f it (x t−1 )) + v t−1 12: x t+1 = x t − ηv t

13:
end for 14: Setx s = x t with t chosen uniformly at random from {0, 1, . . . , m} 15: end for Theorem 22. Suppose that f i (x) are L i -smooth and µ strongly convex. Consider v t defined by in SARAH-ISc (Algorithm 7) with η ≤ 2/(µ +L). Then the following bound holds, ∀ t ≥ 1, By setting η = O(1/L), we derive the linear convergence with the rate of (1 − 1/κ), whereκ =L µ which is an improvement over the previous result of [22], because if we take the optimal stepsize ν = 2 µ+L than we can easily prove that 2µLη µ+L is greater than 2µLmaxη µ+Lmax , with optimal step size, where L max = max i {L i }.
Lemma 25. Suppose that f i (x)'s are L i -smooth and convex. Consider SARAH (Algorithm 7) with η < 2/L. Then we have that for any t ≥ 1,

Proof of Lemma 24
Proof. Let F j be σ algebra that contains all the information of x 0 , . . . , x j as well as v 0 , . . . , v j−1 . For j ≥ 1, we have where the last equality follows from By taking expectation for the above equation, we have Note that ∇f (x 0 ) − v 0 2 = 0. By summing over j = 1, . . . , t (t ≥ 1), we have

Proof of Lemma 25
Proof. For j ≥ 1, we have The consequent equality follows from definition of p i 's and the last equality follows from definition of SARAH-ISc . Taking expectation, we get when η < 2/L. By summing the above inequality over j = 1, . . . , t (t ≥ 1), we have By Lemma 24, we have

Proof of Theorem 21
Proof. For t ≥ 1, we have Using the proof of Lemma 25, for t ≥ 1, we have Note that 1 − 2 ηL < 0 since η < 2/L. The last inequality follows by the strong convexity of f , that is, µ x t − x t−1 ≤ ∇f (x t ) − ∇f (x t−1 ) and the fact that x t = x t−1 − ηv t−1 . By taking the expectation and applying recursively, we have Lemma 28 (Young's inequality). For a, b ∈ R and β > 0 we have Lemma 29 (Jensen's inequality). Let X be a random variable and g(x) be a convex function. Then


extfieldsofstudy:
Computer Science
Medicine
Biology

s2fieldsofstudy:
Computer Science
Biology

text:
Networks of ribosome flow models for modeling and analyzing intracellular traffic

The ribosome flow model with input and output (RFMIO) is a deterministic dynamical system that has been used to study the flow of ribosomes during mRNA translation. The input of the RFMIO controls its initiation rate and the output represents the ribosome exit rate (and thus the protein production rate) at the 3′ end of the mRNA molecule. The RFMIO and its variants encapsulate important properties that are relevant to modeling ribosome flow such as the possible evolution of “traffic jams” and non-homogeneous elongation rates along the mRNA molecule, and can also be used for studying additional intracellular processes such as transcription, transport, and more. Here we consider networks of interconnected RFMIOs as a fundamental tool for modeling, analyzing and re-engineering the complex mechanisms of protein production. In these networks, the output of each RFMIO may be divided, using connection weights, between several inputs of other RFMIOs. We show that under quite general feedback connections the network has two important properties: (1) it admits a unique steady-state and every trajectory converges to this steady-state; and (2) the problem of how to determine the connection weights so that the network steady-state output is maximized is a convex optimization problem. These mathematical properties make these networks highly suitable as models of various phenomena: property (1) means that the behavior is predictable and ordered, and property (2) means that determining the optimal weights is numerically tractable even for large-scale networks. For the specific case of a feed-forward network of RFMIOs we prove an additional useful property, namely, that there exists a spectral representation for the network steady-state, and thus it can be determined without any numerical simulations of the dynamics. We describe the implications of these results to several fundamental biological phenomena and biotechnological objectives.

Scientific RepoRts | (2019) 9:1703 | https://doi.org/10.1038/s41598-018-37864-1 systems biology, biotechnology and evolutionary biology. Mechanistic models of translation are essential for: (1) analyzing the flow of ribosomes along the mRNA molecule; (2) integrating and understanding the rapidly increasing experimental findings related to translation and its role in the dynamical regulation of gene expression (see, e.g. [4][5][6][7][8][9][10][11][12]; and (3) providing a computational testbed for predicting the effects of various manipulations of the genetic machinery. These models describe the dynamics of ribosome flow and include parameters whose values represent the various translation factors that affect the initiation rate and codon decoding times along the mRNA molecule.
Another fundamental biological process based on the flow of biological "machines" along "intracellular roads" is intracellular transport. In this process, vesicles are transferred to particular intracellular locations by molecular motors that haul them along microtubules and actin filaments 1,13 .
We now review two computational models for ribosome flow that are the most relevant for this paper. Numerous other models exist in the literature, see e.g. the survey papers 10,14 .

Totally Asymmetric Simple Exclusion Process (TASEP)
This is a fundamental model in non-equilibrium statistical mechanics that has been extensively used to model and analyze translation and intracellular transport. TASEP is a discrete-time, stochastic model describing particles hopping along an ordered lattice of N sites [15][16][17] . A particle at site i may hop to site i + 1 at a rate γ i but only if this site is empty. This models the fact that the particles have volume and thus cannot overtake one another. Specifically, a particle may hop to the first site at rate α (if the first site is empty), and hop out from the last site at rate β (if the last site is occupied). In the context of translation, the lattice of sites represents the chain of codons in the mRNA, and the hopping particles represent the moving ribosomes 18,19 .
Analysis of TASEP is in general non trivial, and closed-form results have been obtained mainly for the homogeneous TASEP, i.e. the case where all γ i s are assumed to be equal. The non-homogeneous case is typically studied via extensive and time-consuming Monte Carlo simulations.
In each cell, multiple translation processes take place concurrently, utilizing the limited shared translation resources (i.e. ribosomes and translation factors). For example, a yeast cell contains about 60,000 mRNA molecules and about 240,000 ribosomes 20,21 . The competition for shared resources induces indirect interactions and correlations between the various translation processes. Such interactions must be considered when analyzing the cellular economy of the cell, and also when designing synthetic circuits 22,23 .
Analyzing large-scale translation, as opposed to translation of a single isolated mRNA molecule, is thus an important research direction that is recently attracting considerable research attention [24][25][26][27][28][29] . For example, in 28 , a TASEP-based computational network consisting of 400 mRNA species and 14,000 ribosomes has been used to analyze the sensitivity of a translation network to perturbations in the initiation and elongation rates and in the mRNA levels. In 29 , a deterministic mean-field approximation of TASEP, called the ribosome flow model (RFM), was used for studying the effect of fluctuations in the mRNA levels on translation in a whole cell simulation of an S. cerevisiae cell.
Here, we consider large-scale networks of interconnected translation processes, whose building blocks are RFMs with suitable inputs and outputs.

Ribosome Flow Model (RFM)
The RFM 30 is a continuous-time, deterministic model for ribosome flow that can be obtained via a mean-field approximation of TASEP with open boundary conditions (i.e., the two sides of the TASEP lattice are connected to two particle reservoirs) 31 .
In the RFM, the mRNA molecule is coarse-grained into a chain of n consecutive sites of codons. For each site i ∈ {1, 2, …, n}, a state-variable x i (t) ∈ [0, 1] denotes the normalized ribosomal occupancy level (or ribosomal density) at site i at time t, where x i (t) = 1 [x i (t) = 0] means that site i is completely full [empty] at time t. The RFM is characterized by n + 1 positive parameters: the initiation rate (λ 0 ), the transition rate from site i to site i + 1 (λ i ), and the exit rate (λ n ). Ribosomes that attempt to bind to the first site are constrained by the occupancy level of that site, i.e. the effective flow of ribosomes into the first site is given by λ 0 (1 − x 1 ). This means that: (1) the maximal possible entry rate is λ 0 ; and (2) the entry rate decreases as the first site becomes fuller, and becomes zero when the first site is completely full. Similarly, the effective flow of ribosomes from site i to site i + 1 increases [decreases] with the occupancy level at site i [i + 1] and thus is given by λ i x i (1 − x i+1 ). This is a "soft" version of the simple exclusion principle that models the fact that ribosomes have volume and cannot overtake one another, thus as the occupancy level at site i + 1 increases less ribosomes can enter this site, and the effective flow rate from site i to site i + 1 decreases. The (soft) simple exclusion principle in the RFM allows to model the evolution of ribosomal "traffic jams". Indeed, if site i becomes fuller, i.e. x i increases then the flow from site i − 1 to site i decreases and thus site i − 1 also becomes fuller, and so on. Recent findings suggest that in many organisms and conditions a non-negligible percentage of ribosomes tends to be involved in such traffic jams (see, for example 32 ).
The dynamics of the RFM with n sites are given by n nonlinear first-order ordinary differential equations (ODEs) describing the change in the occupancy level of each site as a function of time: Note that the x i s are dimensionless, and the λ i s have units of 1/time. The protein production rate or translation rate is the rate in which ribosomes detach from the mRNA at time t, that is, , then (1) can be written more succinctly as ( 1 ) i ii i 1 , and we omit the dependence in time for clarity. This means that the change in the density at site i is the flow from site i − 1 to site i minus the flow from site i to site i + 1.
Let x(t, a) denote the solution of (1) at time t ≥ 0 for the initial condition x(0) = a. Since the state variables correspond to normalized occupancy levels, we always assume that a belongs to the closed n-dimensional unit cube denoted [0, 1] n . Let (0, 1) n denote the interior of [0, 1] n . In other words, ∈ a (0, 1) n means that every entry a i of a satisfies 0 < a i < 1.
It was shown in 33 that [0, 1] n is an invariant set of the dynamics i.e. if ∈ a [0, 1] n then ∈ x t a ( , ) [0, 1] n for all t ≥ 0. It was also shown that the RFM is a tridiagonal cooperative dynamical system 34,35 , and that this implies that (1) admits a steady-state point λ λ = … ∈ e e( , , ) (0, 1) n n 0 , that is globally asymptotically stable, that is, x t a e a lim ( , ) , for any [0, 1] t n (see also 36 ). In particular, the production rate converges to the steady-state value This means that the parameters of the RFM determine a unique steady-state occupancy at all sites along the mRNA. At the steady-state the flow into every site is equal to the flow out of the site. For any initial density the dynamics converge to this steady-state.
Spectral representation of the RFM steady-state. Simulating the dynamics of large-scale RFMs until (numerical) convergence to the steady-state may be tedious. A useful property of the RFM is that the steady-state can be computed using a spectral approach, that is, based on calculating the eigenvalues and eigenvectors of a suitable matrix 37 . Consider the RFM with dimension n and rates λ 0 , …, λ n . Define the (n + 2) × (n + 2) Jacobi matrix This is a symmetric matrix, so its eigenvalues are real. Since A is componentwise non-negative and irreducible, it admits a unique maximal eigenvalue σ > 0 (called the Perron eigenvalue or Perron root), and the corresponding eigenvector  ζ ∈ + n 2 (the Perron eigenvector) has positive entries 38 .

Theorem 1.
Consider an RFM with dimension n and rates λ 0 , …,λ n . Let A be the matrix defined in (3). Then the steady-state values of the RFM satisfy 37 : In other words, the steady-state density and production rate in the RFM can be obtained from the Perron eigenvalue and eigenvector of A. In particular, this makes it possible to determine r ss and e even for very large chains using efficient and numerically stable algorithms for computing the eigenvalues and eigenvectors of a Jacobi matrix (see, e.g. 39 ).
Consider two RFMs: one with rates λ 0 , …, λ n and the second with rates λ λ ]. Comparing the entries of A and ∼ A, it follows from known results in the Perron-Frobenius theory that σ σ <  . Hence, σ σ > − −  2 2 , so Thm. 1 implies that the steady-state production rate in the RFM increases when one (or more) of the rates increases.
Thm. 1 has several more important implications. For example, it implies that r ss = r ss (λ 0 , …, λ n ) is a strictly concave function on  ++ + n 137 . Also, it implies that the sensitivity of the steady-state with respect to (w.r.t.) a perturbation in the translation rates becomes an eigenvalue sensitivity problem 40 . We refer to the survey 41 for more details.
Reference 42 extended the RFM into a single-input single-output (SISO) control system, by defining the production rate as an output, and by introducing a time-varying input   → + + u: representing the flow of ribosomes from the "outside world" into the mRNA molecule. This is referred to as the RFM with input and output (RFMIO). The RFMIO dynamics is thus described by: Note that u multiplies λ 0 , so that the initiation rate at time t is u(t)λ 0 . We consider λ 0 as modeling an intrinsic bio-physical property of the mRNA, and u(t) as an "outside" effect e.g. the time-varying abundance of "free" ribosomes in the vicinity of the mRNA (see Fig. 1). Throughout, we always assume that u(t) > 0 for all t ≥ 0 in order to avoid some technical problems arising when the initiation rate is zero. Of course, for , with c > 0, the RFMIO becomes an RFM with initiation rate cλ 0 . In this case, the convergence to steady-state represents a form of homeostasis that is sensitive to the value of the input. In particular, the steady-state density of such an RFMIO can be computed using the spectral approach described in Thm. 1.
We write the RFMIO dynamics more succinctly as The RFMIO facilitates modeling a network of interconnected "roads" (e.g. mRNA or DNA molecules, microtubules, etc.), where the flow out of one intracellular road (representing ribosomes, RNAPs, vesicles attached to molecular motors, etc.) may enter another road in the network or re-enter the same road. This enables the analysis of important phenomena such as translation re-initiation 43 , competition for finite resources including the effect of the exit rate from one road on the initiation rate of other intracellular roads (see e.g. 44 ), transport on a network of interconnected microtubules 13 , etc. Such interconnected models are essential for engineering cells for various biotechnological objectives such as optimization of protein production rate, optimization of growth rate, and optimization of traffic jams, as the multiple processes taking place in the cell cannot be analyzed using models of a single "road" 45 .
In the context of translation, the output of an RFMIO represents both the flow of ribosomes out of the mRNA molecule and the synthesis rate of proteins. If the output of a RFMIO is divided into several inputs of other RFMIOs then this represents the distribution of the exiting ribosomes initiating the other mRNAs.
In this paper, we consider networks of interconnected RFMIOs. We show that such dynamical networks provide a useful and versatile modeling tool for many dynamical intracellular traffic phenomena. Our first main result shows that under quite general feedback connections the network admits a unique steady-state and every solution of the dynamics converses to this steady-state. In other words, the network is globally asymptotically stable (GAS). This is important for several reasons. For example, GAS implies that the network admits an ordered behavior and paves the way for analyzing further important questions e.g. how does the steady-state depends on various parameters?
We analyze the problem of maximizing the steady-state output of the network. Specifically, the question we consider is how to determine the interconnection weight values between the RFMIOs in the network so that the network output is maximized. Our second main result shows that this is a convex optimization problem, implying that it can be solved using highly efficient algorithms even for very large networks.
In the specific case of feed-forward networks of RFMIOs, we show an additional property, namely, that we can determine the steady-state of the entire network using a spectral approach, and with no need to numerically solve the dynamical equations.
We note that two previous papers considered specific networks of RFMs. Reference 42 studied the effect of ribosome recycling using a single RFMIO with positive feedback from the output (i.e. production rate) to the input (i.e. initiation rate). Reference 24 analyzed a closed system composed of a dynamic free pool of ribosomes that feeds a single-layer of parallel RFMIOs. This was used as a tool for analyzing the indirect effect between different mRNA molecules due to competition for a shared resource, namely, the pool of free ribosomes. Here, we study networks that provide a significant generalization of these particular models.
The remainder of this paper is organized as follows. The next section describes the network of RFMIOs that we introduce and analyze in this paper. Then we present our main results and demonstrate them using a biological example. The final section summarizes and describes several directions for future research. To increase the readability of this paper, all the proofs are placed in the Appendix. We n , x i is the i-th entry of x, and x′ is the transpose of x. If a time-dependent variable x(t) admits a steady-state then we denote it by x ss , that is,

Networks of RFMIOs.
Consider a network of m interconnected RFMIOs. The input to the network is a source whose output rate is y 0 , and represents external resources that drive the elements in the network. For example, this can represent pools of free ribosomes in the cell. The output of the entire network is denoted by y. This may represent for example the flow of a desired protein produced by the network, the total flow of ribosomes that feed some other process, etc.
For i ∈ {1, …, m}, RFMIO i is a dynamical system with dimension n i , input u i and output y i . For any k ∈ {0, 1, …, m} some or all of the output y k may be connected to the input of another RFMIO, say, RFMIO j with a control parameter (or weight) v k,j ∈ [0, 1]. Here v k,j = 0 means that y k is not connected to u j . The input to RFMIO j is thus We define the total network output y = y m+1 by where v j,m+1 ∈ [0, 1] are the proportion weights from the output y j of RFMIO j to the network output y.
We say that the network is feasible if: (1) every v k,j ∈ [0, 1]; and (2) . The first requirement corresponds to the fact that every v k,j describes the proportion of the output y k that feeds the input of RFMIO j. The second requirement means that ∑ = v y y j k j k k , i.e. the connections indeed describe a distribution of the output y k to other points in the network.
We use x t ( ) i j , i = 1, …, n j , to denote the state-variable describing the occupancy at site i in RFMIO j at time t. The vector We demonstrate using several examples how networks of RFMIOs, with and without feedback connections, can be used to model and study various intracellular networks. As we will see below, our main theoretical result guarantees that all these networks are GAS. Thus, the state-variables in the networks converge to a steady-state that depends on the various parameters, but not on the initial condition.
Our first example describes the efficiency of ribosome recycling in eukaryotic mRNA, and the tradeoff between recycling on the one-hand and the need to "free" ribosomes for other mRNAs on the other-hand. Example 1. Consider the system depicted in Fig. 2. Here a fixed source with rate 0.1 is feeding an RFMIO of length n = 3 with rates λ λ . A proportion v ∈ [0, 1] of the RFMIO output y(t) is fed back into the input, so that the total RFMIO input is u(t) = 0.1 + vy(t). Figure 3 shows the steady-state values y ss (v) and (1 − v)y ss (v) as a function of v. Note that for both these functions there is a unique maximizing value of v. We may interpret y as the total production rate, and (1 − v)y as the rate of ribosomes that are not recycled, and thus can be used to translate other mRNA molecules. Of course, one can also define other functions as the network output, say, some weighted sum of the production rate and the rate of non-recycled ribosomes. In this case, finding the value v that maximizes the steady-state output corresponds to maximizing the production rate on a specific mRNA molecule while still "freeing" enough ribosomes for other purposes.   . Note that the rate λ 2 = 1/4 in RFMIO 1 leads to a "traffic jam" of   ribosomes in this RFMIO, that is, the steady-state densities in the first two sites are high, whereas the density in the third site is low. This yields a low output rate from this RFMIO. The second RFMIO thus converges to a steady-state with low densities.
Example 2 may represents a case of re-initiation: one ORF appears in the 5′UTR of the second ORF and the ribosomes finishing the translation of the first ORF start translating the second one 43,46,47 . In this case, a low elongation rate along the first ORF is expected to yield a low density of ribosomes in the second ORF.

Optimizing the Network Output Rate
Several papers considered optimizing the production rate in a single RFM or some variant of the RFM 37,48-51 . Here, we study a different problem, namely, maximizing the steady-state output in a network of RFMIOs w.r.t. the control (or connection) weights in the network. In other words, the problem is how to distribute the traffic between the different RFMIOs in the network so that the steady-state output is maximized.  The constraints here guarantee the feasibly of the network. However, the results below remain valid even if the second constraint in (6) is replaced by The next result is instrumental for analyzing Problem 1.

Proposition 1. Under a suitable reparametrization Problem 1 becomes a convex optimization problem.
The following example demonstrates this result.

Example 3.
Consider an RFMIO with a single site and rates λ 0 = λ 1 = 1: Suppose that the input is u = 0.1 + vy, where v ∈ [0, 1], that is, there is a feedback connection from the output of the RFMIO back to the input with a weight v. It is straightforward to verify that for any x 1 (0) ∈ [0, 1] the solution x 1 (t) converges to the value It is also straightforward to verify that > e 0 for all v ∈ (0, 1), so e 1 and thus the steady-state output y ss (v) = e 1 (v) is not concave in v. We conclude that the optimization problem: is not a convex optimization problem. Nevertheless, since in this particular case y ss (v) is a scalar function, it is easy to solve this optimization problem yielding (all numerical values in this paper are to four digit accuracy) = = . . (1) 0 2702 (9) ss ss The reparametrization is based on redefining the input as u = 0.1 + w, with the constraint w ∈ [0, y]. Now the steady-state output of (7) ss and this constraint defines a convex set of admissible w's, so (10) is a convex optimization problem. The solution of this problem is obtained at w* = 0.2702 for which y ss (w*) = 0.2702. We conclude that the optimal values correspond to = ⁎ ⁎ w y ss , and this implies that the solution to the optimization problem (8) is v* = 1. Thus, we can obtain the optimal weights from the solution of the reparametrized problem.  The next example demonstrates a synthetic and more complex network that includes feedback connections. Fig. 6 Figure 7 depicts the steady-state output as a function of the control parameters v 1 , v 2 . It may be seen that y is a concave function. The optimal output value y* = 0.8595 is obtained for = . is reasonable, as this implies that all the output y 2 of RFMIO 2 goes directly to the network output y rather than first to RFMIO 4 and from there, indirectly, to y.  To further explain the biological motivation of the optimization problem studied here, consider for example the metabolic pathway or protein complex described in Fig. 8. This includes a set of enzymes/proteins that are involved in a specific stoichiometry (see, for example 12 ). In this case the objective function is of the form b′y, where y is a vector of production rates of the different proteins in the metabolic pathway or protein complex, and b is their stoichiometry vector. Figure 8 depicts an operon with four coding regions on the same transcript. The initiation rate to each ORF is affected by an "external" factor (e.g. the intracellular pool of ribosomes), and also by the "leakage" of ribosomes from the previous ORF. The proteins produced in the operon (P 1 , …, P 4 ) with

Analysis of feed-forward networks.
In this section we further analyze feed-forward networks of RFMIOs, where feed-forward means that for any j the input u j of RFMIO j does not depend either directly or indirectly on the output y j of RFMIO j. In other words, there are no feedback connections. In terms of graph theory, this means that the graph describing the connections is a directed acyclic graph (DAG). For these networks Problem 1 can be solved in a more direct way.

Proposition 2. Consider a feed-forward network of RFMIOs. Let v denote the collection of all control weights. The mapping v → y ss (v) is strictly concave.
The following example demonstrates this result.

Example 5.
Consider the network of two RFMIOs described in Fig. 9. Each RFMIO has dimension n = 3. , and those of RFMIO 2 are [2, 2, 2, 2]. In other words, every rate in RFMIO 1 is slower than the corresponding rate in RFMIO 2. A unit source feeds both RFMIO 1 with u 1 = v, and RFMIO 2 with u 2 = 1 − v. The network output y(t) is defined to be the sum of the two RFMIO outputs. Figure 10 depicts the network steady-state output as a function of v ∈ [0, 1]. It may be seen that y ss (v) is a strictly concave function of v, and in particular that there exists a unique value v* = 0.3971 for which y ss is maximized. This corresponds to feeding a smaller [larger] part of the joint source to RFMIO 1 [RFMIO 2]. This is reasonable, as RFMIO 1 has slower rates than RFMIO 2. Hence, it is possible to "direct" more traffic to RFMIO 2 while still avoiding "traffic jams" in this RFMIO.  This example demonstrates the problem of dividing a common resource in a biological network between several "clients" such that some overall performance measure is optimized. For example, this network can represent the problem of optimizing heterologous protein levels by introducing a copy of the gene to the host genome in multiple locations. This raises the question of the optimal strength of the initiation rate of the different gene copies (engineered via manipulation of the 5′UTR and the beginning of the ORF, e.g. via the introduction of Shine-Dalgarno sequence with different strengths and manipulation of the mRNA folding in this region) 52,53 .

Spectral Representation of the Network Steady-State
Recall that in a single RFM it is possible to obtain the steady-state density (and thus the steady-state production rate) using a spectral approach and without numerically simulating the dynamics. The same property immediately carries over to feed-forward networks of RFMIOs. To explain this, consider an RFMIO, say RFMIO j, that is fed only by a constant source. This is just an RFM and its steady-state density and output can be calculated as in Thm. 1. Now consider an RFMIO that is fed by the output of RFMIO j. Its input converges to a steady-state value u ss , and since the RFMIO is contractive 36 , its density converges to a steady-state that is identical to the steady-state of an RFMIO with the constant input u(t) ≡ u ss (see e.g. 54,55 ). We can now determine the steady-state values in the consecutive RFMIOs and so on. The next example demonstrates this for a simple network. [1/2, 3/2 , 5/2, 6, 1] ], and using Thm. 1, we conclude that the steady-state density in RFMIO 1 is = ′ e [5/6 4/5 1/6] (compare with Fig. 5), and the steady-state output is y ss = σ −2 = 1/6. We can now analyze the second RFMIO. Its rates are all one and the input is the output of RFMIO 1, so at steady-state the rates are ′ [1/6 1 1 1] . The corresponding Jacobi matrix is  Thm. 3 implies in particular that the optimal solution is unique. Moreover, there exist highly efficient numerical algorithms for computing the unique solution even for very large networks.

A Biological Example
In order to demonstrate how our model can be used to address questions arising in synthetic biology, we consider the problem of maintaining high growth rates for both a highly expressed heterologous protein and a highly expressed endogenous protein in the cell. These issues are currently attracting considerable interest as lack of understanding of the burden of expressing additional genes affects our ability to predictively engineer cells (see e.g. 57 and the references therein).
Specifically, we consider the problem of maximizing the sum of the steady-state production rates of both heterologous and endogenous genes, under the assumption that they share a common ribosomal resource. We assume that the coding regions of the endogenous gene (which may include various regulatory signals), and the coding region of the heterologous gene (which is optimized to include the most efficient codons) cannot be modified. However, it is possible to engineer the UTRs of these genes in order to modulate their initiation rates. We demonstrate how this biological problem can be modeled and analyzed in the framework of our model.
The endogenous gene is the highly expressed S. cerevisiae gene YGR192C that encodes the protein TDH3, which is involved in glycolysis and gluconeogenesis. The heterologous gene is the green fluorescent protein (GFP) gene (the GFP protein sequence is from gi:1543069), optimized for yeast (i.e. its codons composition was synonymously modified to consist of optimized yeast codons). The YGR192C gene ORF consists of 332 codons, and the GFP gene ORF of 239 codons. The simultaneous translation of these two genes while using a shared resource is modeled as depicted in Fig. 9, where RFMIO 1 (fed by an input u = v) models the translation of the YGR192C gene, and RFMIO 2 (fed by the input u = 1 − v) models the translation of the GFP gene. Here v ∈ [0, 1] is a parameter that determines the relative amount of ribosomal resources allocated to each gene.
Similarly to the approach used in 30 , we divide the YGR192C [GFP] mRNA sequnce into 33 [23] consecutive subsequences: the first subsequence includes the first nine codons (that are also related to later stages of initiation 58 ). The other subsequences include 10 non-overlapping codons each, except for the last subsequence in the YGR192C gene that includes 13 codons. This partitioning was found to optimize the correlation between the RFMIO predictions and biological data.
We model the translation of the YGR192C [GFP] gene using an RFMIO with n = 32 [n = 23] sites. To determine the RFMIO paramteres we first estimate the elongation rates λ 1 , …, λ n , using ribo-seq data for the codon decoding rates 59 , normalized so that the median elongation rate of all S. cerevisiae mRNAs becomes 6.4 codons per second 60 . The site rate is (site time) −1 , where site time is the sum over the decoding times of all the codons in this site. These rates thus depend on various factors including availability of tRNA molecules, amino acids, Aminoacyl tRNA synthetase activity and concentration, and local mRNA folding 1,58,59 .
The initiation rate (that corresponds to the first subsequence) for the YGR192C gene is estimated based on the ribosome density per mRNA levels, as this value is expected to be approximately proportional to the initiation rate when initiation is rate limiting 30,61 . Again, we applied a normalization that brings the median initiation rate of all S. cerevisiae mRNAs to 0.8 mRNAs per second 62 , and this results in an initiation rate of 2.1958 for the YGR192C gene. The GFP initiation rate was set to 0.8. A calculation shows that when each gene is modeled separately using an RFMIO with u = 1, the steady-state production rate of the gene YGR192C [GFP] is r ss = 0.1859 [r ss = 0.1892]. Figure 11 depicts the network output y ss (v), as a function of v ∈ [0, 1]. The unique maximum y ss (v*) = 0.3429 is attained for v* = 0.4311, which corresponds to feeding a smaller [larger] part of the common ribosomal resource Figure 11. Steady-state network output y ss (v) for the of simultaneous translation of YGR192C (RFMIO 1) and GFP (RFMIO 2) genes, using the network depicted in Fig. 9, as a function of the control parameter v. to the GFP [YGR192C] gene. This is reasonable, as the steady-state production rate of the GFP gene is slightly larger than the steady-state production rate of the YGR192C gene. This result implies that in order to maximize the sum of the steady-state production rates of the YGR192C gene and the GFP gene, using a common ribosomal resource, their UTRs binding efficiency should be engineered such that 43% of the ribosomal resource initiates the YGR192C mRNAs, and the remaining 57% initiates the GFP mRNAs. Our analytical approach can also be used to determine the ribosomal allocation that maximizes some weighted sum of the two production rates.

Discussion
Studying the flow of biological "machines" like ribosomes, RNAPs, or motor proteins along biological networks like interconnected mRNA molecules or filaments is of paramount importance. These biological machines have volume and thus satisfy a simple exclusion principle: two machines cannot be in the same place at the same time.
In order to better understand these cellular biological processes it is important to study the flow of biological "machines" along networks of interconnected "roads", and not only in isolated processes. We propose to model and analyze such phenomena using networks of interconnected RFMIOs. The RFMIO dynamics satisfies a "soft" simple exclusion principle: as a site becomes fuller the effective entry rate into this site decreases. In particular, "traffic jams" may evolve behind slowly moving machines. The input and output of every RFMIO facilitate their integration into interconnected networks that can be represented using a graph. The nodes in this graph represent the different RFMIOs, and the (weighted) edges describe how the output of each RFMIO is divided between the inputs of other (or the same) RFMIO. Our main result shows that under quite general positive feedback connections, such a network always admits a unique steady-state and any trajectory of the system converges to this steady-state. This opens the door to many interesting research questions, e.g., how does this steady-state depends on the various parameters in the network, and for what feedback connections is the steady-state output of the network optimized?
We demonstrated using various examples how such networks can model interesting biological phenomena like competition for shared resources, the optimal distribution of a shared biological resource between several "clients", optimizing the effect of ribosome recycling, and more.
The RFM is amenable to rigorous analysis, even when the rates are not homogeneous, using various tools from systems and control theory like contraction theory 63,64 , the theory of cooperative dynamical systems 24,48,65,66 , convex analysis and more. This amenability to rigorous analysis carries over to networks of RFMIOs. For example, the problem of how to connect the RFMIO outputs in the network to inputs so that the steady-state network output is maximized can be cast as a convex optimization problem. This means that the problem can be solved efficiently even for large networks (see, e.g. 56 ).
The networks we propose here allow modeling complex biological processes in a coherent and useful manner. The network models static connections between the RFMIOs. The dynamical part is described by the set of ODEs for each RFMIO. The parameters used in these models can be inferred based on various sources of large-scale genomic data (see, for example 59,67 ) and/or can be predicted directly from the nucleotide sequence of the gene (see, for example 53,68 ). In addition, the analyzed network can be built gradually, one module after another. For example, one can engineer and study one metabolic pathway and then connect another pathway to the existing module (see Fig. 8), etc. This yields a combined model that describes both biophysical aspects of gene expression regulation (e.g. translation), and properties of metabolism (e.g. stoichiometry of enzymes and metabolites, and rates of metabolic reactions).
Topics for further research include networks where the weighted connections between the RFMIOs may also change with time. This may model for example mRNA molecules that diffuse through the cell and consequently change their interactions with ribosomes, other mRNAs, etc. (see. e.g. 69 ).
Finally, networks of interconnected TASEPs have been used to model other natural and artificial phenomena such as vehicular traffic and evacuation dynamics. We believe that the deterministic networks proposed here can also be applied to model and analyze such phenomena.